-
Lukas Bentkamp authoredLukas Bentkamp authored
Gauss_field_test.cpp 7.49 KiB
/**********************************************************************
* *
* Copyright 2019 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
* This file is part of TurTLE. *
* *
* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
**********************************************************************/
#include <string>
#include <cmath>
#include <random>
#include "Gauss_field_test.hpp"
#include "scope_timer.hpp"
template <typename rnumber,
field_backend be,
field_components fc,
kspace_dealias_type dt>
int make_gaussian_random_field(
kspace<be, dt> *kk,
field<rnumber, be, fc> *output_field,
const int rseed,
const double slope,
const double k_cutoff,
const double coefficient)
{
TIMEZONE("make_gaussian_random_field");
// initialize a separate random number generator for each thread
std::vector<std::mt19937_64> rgen;
std::uniform_real_distribution<rnumber> rdist;
rgen.resize(omp_get_max_threads());
// seed random number generators such that no seed is ever repeated if we change the value of rseed.
// basically use a multi-dimensional array indexing technique to come up with actual seed.
for (int thread_id=0; thread_id < omp_get_max_threads(); thread_id++)
{
int current_seed = (
rseed*omp_get_max_threads()*output_field->clayout->nprocs +
output_field->clayout->myrank*omp_get_max_threads() +
thread_id);
DEBUG_MSG("in make_gaussian_random_field, thread_id = %d, current_seed = %d\n", thread_id, current_seed);
rgen[thread_id].seed(current_seed);
}
output_field->real_space_representation = false;
*output_field = 0.0;
// inside loop use only thread-local random number generator
kk->CLOOP_K2([&](
ptrdiff_t cindex,
ptrdiff_t xindex,
ptrdiff_t yindex,
ptrdiff_t zindex, double k2)
{
if (k2 > 0)
switch(fc)
{
case ONE:
{
double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]);
output_field->cval(cindex,0) = coefficient * cos(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
output_field->cval(cindex,1) = coefficient * sin(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
break;
}
case THREE:
for (int cc = 0; cc<3; cc++)
{
double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]);
output_field->cval(cindex,cc,0) = coefficient * cos(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
output_field->cval(cindex,cc,1) = coefficient * sin(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
}
break;
case THREExTHREE:
for (int cc = 0; cc<3; cc++)
for (int ccc = 0; ccc<3; ccc++)
{
double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]);
output_field->cval(cindex,cc,ccc,0) = coefficient * cos(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
output_field->cval(cindex,cc,ccc,1) = coefficient * sin(temp_phase) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
}
break;
}
});
output_field->symmetrize();
return EXIT_SUCCESS;
}
template <typename rnumber>
int Gauss_field_test<rnumber>::initialize(void)
{
TIMEZONE("Gauss_field_test::initialize");
this->read_parameters();
this->scalar_field = new field<rnumber, FFTW, ONE>(
nx, ny, nz,
this->comm,
FFTW_ESTIMATE);
this->vector_field = new field<rnumber, FFTW, THREE>(
nx, ny, nz,
this->comm,
FFTW_ESTIMATE);
this->kk = new kspace<FFTW, SMOOTH>(
this->scalar_field->clayout, this->dkx, this->dky, this->dkz);
if (this->myrank == 0)
{
hid_t stat_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDWR,
H5P_DEFAULT);
this->kk->store(stat_file);
H5Fclose(stat_file);
}
return EXIT_SUCCESS;
}
template <typename rnumber>
int Gauss_field_test<rnumber>::finalize(void)
{
TIMEZONE("Gauss_field_test::finalize");
delete this->vector_field;
delete this->scalar_field; delete this->kk;
return EXIT_SUCCESS;
}
template <typename rnumber>
int Gauss_field_test<rnumber>::read_parameters()
{
TIMEZONE("Gauss_field_test::read_parameters");
this->test::read_parameters();
hid_t parameter_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
H5P_DEFAULT);
this->filter_length = hdf5_tools::read_value<double>(parameter_file, "/parameters/filter_length");
H5Fclose(parameter_file);
return EXIT_SUCCESS;
}
template <typename rnumber>
int Gauss_field_test<rnumber>::do_work(void)
{
TIMEZONE("Gauss_field_test::do_work");
make_gaussian_random_field(this->kk, this->scalar_field);
make_gaussian_random_field(this->kk, this->vector_field, 3);
DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata(), true);
DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
make_gaussian_random_field(this->kk, this->vector_field, 3);
DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
return EXIT_SUCCESS;
}
template class Gauss_field_test<float>;
template class Gauss_field_test<double>;