add symmetrize test files

bfps/cpp/full_code/symmetrize_test.cpp

+#include <string>
+#include <cmath>
+#include <random>
+#include "symmetrize_test.hpp"
+#include "scope_timer.hpp"
+
+
+template <typename rnumber>
+int symmetrize_test<rnumber>::initialize(void)
+{
+    this->read_parameters();
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int symmetrize_test<rnumber>::finalize(void)
+{
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int symmetrize_test<rnumber>::read_parameters()
+{
+    this->test::read_parameters();
+    hid_t parameter_file = H5Fopen(
+            (this->simname + std::string(".h5")).c_str(),
+            H5F_ACC_RDONLY,
+            H5P_DEFAULT);
+    this->random_seed = hdf5_tools::read_value<int>(
+            parameter_file, "/parameters/random_seed");
+    H5Fclose(parameter_file);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int symmetrize_test<rnumber>::do_work(void)
+{
+    // allocate
+    field<rnumber, FFTW, THREE> *scal_field0 = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            DEFAULT_FFTW_FLAG);
+    field<rnumber, FFTW, THREE> *scal_field1 = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            DEFAULT_FFTW_FLAG);
+    std::default_random_engine rgen;
+    std::normal_distribution<rnumber> rdist;
+    rgen.seed(1);
+    kspace<FFTW,SMOOTH> *kk = new kspace<FFTW, SMOOTH>(
+            scal_field0->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_field0
+    *scal_field0 = 0.0;
+    scal_field0->real_space_representation = false;
+    kk->CLOOP_K2(
+            [&](ptrdiff_t cindex,
+                ptrdiff_t xindex,
+                ptrdiff_t yindex,
+                ptrdiff_t zindex,
+                double k2){
+            //if (k2 < kk->kM2)
+            {
+                scal_field0->cval(cindex, 0, 0) = rdist(rgen);
+                scal_field0->cval(cindex, 0, 1) = rdist(rgen);
+                scal_field0->cval(cindex, 1, 0) = rdist(rgen);
+                scal_field0->cval(cindex, 1, 1) = rdist(rgen);
+                scal_field0->cval(cindex, 2, 0) = rdist(rgen);
+                scal_field0->cval(cindex, 2, 1) = rdist(rgen);
+            }
+            if (k2 > 0)
+            {
+                scal_field0->cval(cindex, 0, 0) /= sqrt(k2);
+                scal_field0->cval(cindex, 0, 1) /= sqrt(k2);
+                scal_field0->cval(cindex, 1, 0) /= sqrt(k2);
+                scal_field0->cval(cindex, 1, 1) /= sqrt(k2);
+                scal_field0->cval(cindex, 2, 0) /= sqrt(k2);
+                scal_field0->cval(cindex, 2, 1) /= sqrt(k2);
+            }
+            else
+            {
+                scal_field0->cval(cindex, 0, 0) = 0;
+                scal_field0->cval(cindex, 0, 1) = 0;
+                scal_field0->cval(cindex, 1, 0) = 0;
+                scal_field0->cval(cindex, 1, 1) = 0;
+                scal_field0->cval(cindex, 2, 0) = 0;
+                scal_field0->cval(cindex, 2, 1) = 0;
+            }
+            });
+    // dealias (?!)
+    //kk->template low_pass<rnumber, THREE>(scal_field0->get_cdata(), kk->kM);
+    kk->template dealias<rnumber, THREE>(scal_field0->get_cdata());
+    // make the field divergence free
+    kk->template force_divfree<rnumber>(scal_field0->get_cdata());
+    // apply symmetrize to scal_field0
+    scal_field0->symmetrize();
+
+
+    // make copy in scal_field1
+    // this MUST be made after symmetrizing scal_field0
+    // (alternatively, we may symmetrize scal_field1 as well before the ift-dft cycle
+    scal_field1->real_space_representation = false;
+    *scal_field1 = scal_field0->get_cdata();
+
+    // go back and forth with scal_field1, to enforce symmetry
+    scal_field1->ift();
+    scal_field1->dft();
+    scal_field1->normalize();
+
+    // now compare the two fields
+    double max_diff = 0;
+    ptrdiff_t ix, iy, iz;
+    double k_at_max_diff = 0;
+    double a0, a1;
+
+    kk->CLOOP_K2(
+            [&](ptrdiff_t cindex,
+                ptrdiff_t xindex,
+                ptrdiff_t yindex,
+                ptrdiff_t zindex,
+                double k2){
+            double diff_re0 = scal_field0->cval(cindex, 0, 0) - scal_field1->cval(cindex, 0, 0);
+            double diff_re1 = scal_field0->cval(cindex, 1, 0) - scal_field1->cval(cindex, 1, 0);
+            double diff_re2 = scal_field0->cval(cindex, 2, 0) - scal_field1->cval(cindex, 2, 0);
+            double diff_im0 = scal_field0->cval(cindex, 0, 1) - scal_field1->cval(cindex, 0, 1);
+            double diff_im1 = scal_field0->cval(cindex, 1, 1) - scal_field1->cval(cindex, 1, 1);
+            double diff_im2 = scal_field0->cval(cindex, 2, 1) - scal_field1->cval(cindex, 2, 1);
+            double diff = sqrt(diff_re0*diff_re0 + diff_re1*diff_re1 + diff_re2*diff_re2 +
+                               diff_im0*diff_im0 + diff_im1*diff_im1 + diff_im2*diff_im2);
+            double amplitude0 = (scal_field0->cval(cindex, 0, 0)*scal_field0->cval(cindex, 0, 0) +
+                                 scal_field0->cval(cindex, 1, 0)*scal_field0->cval(cindex, 1, 0) +
+                                 scal_field0->cval(cindex, 2, 0)*scal_field0->cval(cindex, 2, 0) +
+                                 scal_field0->cval(cindex, 0, 1)*scal_field0->cval(cindex, 0, 1) +
+                                 scal_field0->cval(cindex, 1, 1)*scal_field0->cval(cindex, 1, 1) +
+                                 scal_field0->cval(cindex, 2, 1)*scal_field0->cval(cindex, 2, 1));
+            double amplitude1 = (scal_field1->cval(cindex, 0, 0)*scal_field1->cval(cindex, 0, 0) +
+                                 scal_field1->cval(cindex, 1, 0)*scal_field1->cval(cindex, 1, 0) +
+                                 scal_field1->cval(cindex, 2, 0)*scal_field1->cval(cindex, 2, 0) +
+                                 scal_field1->cval(cindex, 0, 1)*scal_field1->cval(cindex, 0, 1) +
+                                 scal_field1->cval(cindex, 1, 1)*scal_field1->cval(cindex, 1, 1) +
+                                 scal_field1->cval(cindex, 2, 1)*scal_field1->cval(cindex, 2, 1));
+            double amplitude = sqrt((amplitude0 + amplitude1)/2);
+            if (amplitude > 0)
+            if (diff/amplitude > max_diff)
+            {
+                max_diff = diff / amplitude;
+                ix = xindex;
+                iy = yindex + scal_field0->clayout->starts[0];
+                iz = zindex;
+                k_at_max_diff = sqrt(k2);
+                a0 = sqrt(amplitude0);
+                a1 = sqrt(amplitude1);
+            }
+            });
+    DEBUG_MSG("found maximum relative difference %g at ix = %ld, iy = %ld, iz = %ld, wavenumber = %g, amplitudes %g %g\n",
+            max_diff, ix, iy, iz, k_at_max_diff, a0, a1);
+
+    // deallocate
+    delete kk;
+    delete scal_field1;
+    delete scal_field0;
+    return EXIT_SUCCESS;
+}
+
+template class symmetrize_test<float>;
+template class symmetrize_test<double>;
+

bfps/cpp/full_code/symmetrize_test.hpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2018 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                                 *
+*                                                                     *
+**********************************************************************/
+
+
+
+#ifndef SYMMETRIZE_TEST_HPP
+#define SYMMETRIZE_TEST_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "kspace.hpp"
+#include "field.hpp"
+#include "full_code/test.hpp"
+
+/** \brief A class for testing basic field class functionality.
+ */
+
+template <typename rnumber>
+class symmetrize_test: public test
+{
+    public:
+        int random_seed;
+
+        symmetrize_test(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            test(
+                    COMMUNICATOR,
+                    simulation_name){}
+        ~symmetrize_test(){}
+
+        int initialize(void);
+        int do_work(void);
+        int finalize(void);
+        int read_parameters(void);
+};
+
+#endif//SYMMETRIZE_TEST_HPP
+