Skip to content
Snippets Groups Projects
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>;