-
Cristian Lalescu authoredCristian Lalescu authored
field_test.cpp 5.10 KiB
/******************************************************************************
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
* This file is part of bfps. *
* *
* bfps 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. *
* *
* bfps 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 bfps. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
******************************************************************************/
#include <string>
#include <cmath>
#include <random>
#include "field_test.hpp"
#include "scope_timer.hpp"
template <typename rnumber>
int field_test<rnumber>::initialize(void)
{
TIMEZONE("field_test::initialize");
this->read_parameters();
return EXIT_SUCCESS;
}
template <typename rnumber>
int field_test<rnumber>::finalize(void)
{
TIMEZONE("field_test::finalize");
this->read_parameters();
return EXIT_SUCCESS;
}
template <typename rnumber>
int field_test<rnumber>::read_parameters()
{
TIMEZONE("field_test::read_parameters");
this->test::read_parameters();
// in case any parameters are needed, this is where they should be read
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 field_test<rnumber>::do_work(void)
{
TIMEZONE("field_test::do_work");
// allocate
field<rnumber, FFTW, ONE> *scal_field = new field<rnumber, FFTW, ONE>(
this->nx, this->ny, this->nz,
this->comm,
FFTW_ESTIMATE);
field<rnumber, FFTW, ONE> *scal_field_alt = new field<rnumber, FFTW, ONE>(
this->nx, this->ny, this->nz,
this->comm,
FFTW_ESTIMATE);
std::default_random_engine rgen;
std::normal_distribution<rnumber> rdist;
rgen.seed(2);
//auto gaussian = std::bind(rgen, rdist);
kspace<FFTW,SMOOTH> *kk = new kspace<FFTW, SMOOTH>(
scal_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);
kk->store(stat_file);
H5Fclose(stat_file);
}
// fill up scal_field
scal_field->real_space_representation = true;
scal_field->RLOOP(
[&](ptrdiff_t rindex,
ptrdiff_t xindex,
ptrdiff_t yindex,
ptrdiff_t zindex){
scal_field->rval(rindex) = rdist(rgen);
});
*scal_field_alt = scal_field->get_rdata();
double L2r = scal_field->L2norm(kk);
scal_field->dft();
double L2c = scal_field->L2norm(kk);
scal_field->ift();
scal_field->normalize();
DEBUG_MSG("L2r = %g, L2c = %g\n",
L2r, L2c / scal_field->npoints);
double max_error = 0;
scal_field->RLOOP(
[&](ptrdiff_t rindex,
ptrdiff_t xindex,
ptrdiff_t yindex,
ptrdiff_t zindex){
double tval = fabs(scal_field->rval(rindex) - scal_field_alt->rval(rindex));
if (max_error < tval)
max_error = tval;
});
DEBUG_MSG("maximum error is %g\n", max_error);
scal_field->dft();
kk->template dealias<rnumber, ONE>(scal_field->get_cdata());
scal_field->symmetrize();
scal_field->normalize();
L2c = scal_field->L2norm(kk);
scal_field->ift();
L2r = scal_field->L2norm(kk);
DEBUG_MSG("L2r = %g, L2c = %g\n",
L2r, L2c);
// deallocate
delete kk;
delete scal_field;
delete scal_field_alt;
return EXIT_SUCCESS;
}
template class field_test<float>;
template class field_test<double>;