diff --git a/CMakeLists.txt b/CMakeLists.txt
index e786a2b1bbb2a8bce1c0e55724b46eb0cb7cfe20..cf9556bdf30cf44a77cd0846b4bc006a60784ebf 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -192,6 +192,7 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/filter_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_test.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/Gauss_field_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_output_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/get_rfields.cpp
@@ -233,6 +234,7 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/joint_acc_vel_stats.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/filter_test.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/Gauss_field_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_output_test.hpp
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 6a728545d4b9bf53727404809150a56e2637c974..6c5e5a9fbd6c9e61da21034aa82ed2e15928a934 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -257,6 +257,9 @@ class TEST(_code):
dest = 'TEST_class',
help = 'type of simulation to run')
subparsers.required = True
+ parser_Gauss_field_test = subparsers.add_parser(
+ 'Gauss_field_test',
+ help = 'test incompressible Gaussian random fields')
parser_filter_test = subparsers.add_parser(
'filter_test',
help = 'plain filter test')
@@ -272,7 +275,8 @@ class TEST(_code):
parser_test_interpolation = subparsers.add_parser(
'test_interpolation',
help = 'test velocity gradient interpolation')
- for parser in ['parser_filter_test',
+ for parser in ['parser_Gauss_field_test',
+ 'parser_filter_test',
'parser_field_test',
'parser_symmetrize_test',
'parser_field_output_test',
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..521706f7b6809a712553f02f999c41b82016ab96
--- /dev/null
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -0,0 +1,157 @@
+/**********************************************************************
+* *
+* 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 = 0)
+{
+ TIMEZONE("generate gaussian random field");
+ thread_local std::mt19937_64 rgen;
+ std::normal_distribution<rnumber> rdist;
+ rgen.seed(rseed + output_field->clayout->myrank*omp_get_max_threads() + omp_get_thread_num());
+ output_field->real_space_representation = false;
+ *output_field = 0.0;
+ kk->CLOOP_K2([&](
+ ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2)
+ {
+ if (k2 > 0)
+ switch(fc)
+ {
+ case ONE:
+ output_field->cval(cindex,0) = rdist(rgen);
+ output_field->cval(cindex,1) = rdist(rgen);
+ break;
+ case THREE:
+ for (int cc = 0; cc<3; cc++)
+ {
+ output_field->cval(cindex,cc,0) = rdist(rgen);
+ output_field->cval(cindex,cc,1) = rdist(rgen);
+ }
+ break;
+ case THREExTHREE:
+ for (int cc = 0; cc<3; cc++)
+ for (int ccc = 0; ccc<3; ccc++)
+ {
+ output_field->cval(cindex,cc,ccc,0) = rdist(rgen);
+ output_field->cval(cindex,cc,ccc,1) = rdist(rgen);
+ }
+ 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);
+ DEBUG_MSG("L2norm = %g\n", this->vector_field->L2norm(this->kk));
+ this->kk->template rotate_divfree<rnumber>(this->vector_field->get_cdata());
+ 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);
+ 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>;
+
diff --git a/cpp/full_code/Gauss_field_test.hpp b/cpp/full_code/Gauss_field_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3854bb1f48f8b3a724a67d77fe370372901e299a
--- /dev/null
+++ b/cpp/full_code/Gauss_field_test.hpp
@@ -0,0 +1,70 @@
+/**********************************************************************
+* *
+* 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 *
+* *
+**********************************************************************/
+
+
+
+#ifndef GAUSS_FIELD_TEST_HPP
+#define GAUSS_FIELD_TEST_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "kspace.hpp"
+#include "field.hpp"
+#include "full_code/test.hpp"
+
+/** \brief A test of incompressible Gaussian random fields.
+ *
+ */
+
+template <typename rnumber>
+class Gauss_field_test: public test
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ double filter_length;
+
+ /* other stuff */
+ kspace<FFTW, SMOOTH> *kk;
+ field<rnumber, FFTW, ONE> *scalar_field;
+ field<rnumber, FFTW, THREE> *vector_field;
+
+ Gauss_field_test(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ test(
+ COMMUNICATOR,
+ simulation_name){}
+ ~Gauss_field_test(){}
+
+ int initialize(void);
+ int do_work(void);
+ int finalize(void);
+ int read_parameters(void);
+};
+
+#endif//GAUSS_FIELD_TEST_HPP
+
diff --git a/cpp/kspace.cpp b/cpp/kspace.cpp
index 86d692a25c099ab5d97dd0fd5b247533815e61b6..610df6de6d6997646b0403ffe46c33e2e3c0de1f 100644
--- a/cpp/kspace.cpp
+++ b/cpp/kspace.cpp
@@ -527,9 +527,9 @@ void kspace<be, dt>::dealias(typename fftw_interface<rnumber>::complex *__restri
template <field_backend be,
kspace_dealias_type dt>
template <typename rnumber>
-void kspace<be, dt>::force_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a)
+void kspace<be, dt>::project_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a)
{
- TIMEZONE("kspace::force_divfree");
+ TIMEZONE("kspace::project_divfree");
this->CLOOP_K2(
[&](ptrdiff_t cindex,
ptrdiff_t xindex,
@@ -558,6 +558,83 @@ void kspace<be, dt>::force_divfree(typename fftw_interface<rnumber>::complex *__
std::fill_n((rnumber*)(a), 6, 0.0);
}
+/** \brief Rotate vector modes perpendicular to wavenumber
+ * This is different from project because it maintains the energy of the field,
+ * I want it in order to be able to generate random fields with prescribed
+ * spectra.
+ */
+template <field_backend be,
+ kspace_dealias_type dt>
+template <typename rnumber>
+void kspace<be, dt>::rotate_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a)
+{
+ TIMEZONE("kspace::rotate_divfree");
+ this->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+ if (k2 > 0)
+ {
+ double cosTheta;
+ double usize;
+ // first compute real part of dot product u0.k
+ cosTheta = (this->kx[xindex]*((*(a + cindex*3 ))[0]) +
+ this->ky[yindex]*((*(a + cindex*3+1))[0]) +
+ this->kz[zindex]*((*(a + cindex*3+2))[0]) );
+ // now compute size of initial velocity vector
+ usize = sqrt(pow((*(a + cindex*3 ))[0], 2) +
+ pow((*(a + cindex*3+1))[0], 2) +
+ pow((*(a + cindex*3+2))[0], 2) +
+ pow((*(a + cindex*3 ))[1], 2) +
+ pow((*(a + cindex*3+1))[1], 2) +
+ pow((*(a + cindex*3+2))[1], 2));
+ // finalize computation of cos Theta
+ if (usize != 0){
+ cosTheta /= (usize * sqrt(k2));
+ // now compute cross product
+ // cross product for complex vectors is complex conjugate of regular cross product
+ double cp0[3], cp1[3];
+ cp0[0] = (*(a + cindex*3+1))[0]*this->kz[zindex] - (*(a + cindex*3+2))[0]*this->ky[yindex];
+ cp0[1] = (*(a + cindex*3+2))[0]*this->kx[xindex] - (*(a + cindex*3+0))[0]*this->kz[zindex];
+ cp0[2] = (*(a + cindex*3+0))[0]*this->ky[yindex] - (*(a + cindex*3+1))[0]*this->kx[xindex];
+ cp1[0] = -((*(a + cindex*3+1))[1]*this->kz[zindex] - (*(a + cindex*3+2))[1]*this->ky[yindex]);
+ cp1[1] = -((*(a + cindex*3+2))[1]*this->kx[xindex] - (*(a + cindex*3+0))[1]*this->kz[zindex]);
+ cp1[2] = -((*(a + cindex*3+0))[1]*this->ky[yindex] - (*(a + cindex*3+1))[1]*this->kx[xindex]);
+ double cpsize = sqrt(cp0[0]*cp0[0] + cp0[1]*cp0[1] + cp0[2]*cp0[2] +
+ cp1[0]*cp1[0] + cp1[1]*cp1[1] + cp1[2]*cp1[2]);
+ cp0[0] /= cpsize;
+ cp0[1] /= cpsize;
+ cp0[2] /= cpsize;
+ cp1[0] /= cpsize;
+ cp1[1] /= cpsize;
+ cp1[2] /= cpsize;
+ double sinTheta = sqrt(1 - cosTheta*cosTheta);
+ // store initial velocity
+ double u0[3], u1[3];
+ u0[0] = (*(a + cindex*3+0))[0];
+ u1[0] = (*(a + cindex*3+0))[1];
+ u0[1] = (*(a + cindex*3+1))[0];
+ u1[1] = (*(a + cindex*3+1))[1];
+ u0[2] = (*(a + cindex*3+2))[0];
+ u1[2] = (*(a + cindex*3+2))[1];
+ // store final result
+ // cross product is between two complex vectors, that's why it looks more complicated
+ (*(a + cindex*3+0))[0] = u0[0]*sinTheta + ((cp0[1]*u0[2] - cp0[2]*u0[1]) - (cp1[1]*u1[2] - cp1[2]*u1[1]))*cosTheta;
+ (*(a + cindex*3+1))[0] = u0[1]*sinTheta + ((cp0[2]*u0[0] - cp0[0]*u0[2]) - (cp1[2]*u1[0] - cp1[0]*u1[2]))*cosTheta;
+ (*(a + cindex*3+2))[0] = u0[2]*sinTheta + ((cp0[0]*u0[1] - cp0[1]*u0[0]) - (cp1[0]*u1[1] - cp1[1]*u1[0]))*cosTheta;
+ (*(a + cindex*3+0))[1] = u1[0]*sinTheta - ((cp0[1]*u1[2] - cp0[2]*u1[1]) + (cp1[1]*u0[2] - cp1[2]*u0[1]))*cosTheta;
+ (*(a + cindex*3+1))[1] = u1[1]*sinTheta - ((cp0[2]*u1[0] - cp0[0]*u1[2]) + (cp1[2]*u0[0] - cp1[0]*u0[2]))*cosTheta;
+ (*(a + cindex*3+2))[1] = u1[2]*sinTheta - ((cp0[0]*u1[1] - cp0[1]*u1[0]) + (cp1[0]*u0[1] - cp1[1]*u0[0]))*cosTheta;
+ }
+ }
+ }
+ );
+ if (this->layout->myrank == this->layout->rank[0][0])
+ std::fill_n((rnumber*)(a), 6, 0.0);
+}
+
template <field_backend be,
kspace_dealias_type dt>
template <typename rnumber,
@@ -1047,8 +1124,18 @@ template double kspace<FFTW, SMOOTH>::L2norm<double, THREE>(
template double kspace<FFTW, SMOOTH>::L2norm<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a);
+template void kspace<FFTW, SMOOTH>::project_divfree<float>(
+ typename fftw_interface<float>::complex *__restrict__ a);
+template void kspace<FFTW, SMOOTH>::project_divfree<double>(
+ typename fftw_interface<double>::complex *__restrict__ a);
+
template void kspace<FFTW, SMOOTH>::force_divfree<float>(
typename fftw_interface<float>::complex *__restrict__ a);
template void kspace<FFTW, SMOOTH>::force_divfree<double>(
typename fftw_interface<double>::complex *__restrict__ a);
+template void kspace<FFTW, SMOOTH>::rotate_divfree<float>(
+ typename fftw_interface<float>::complex *__restrict__ a);
+template void kspace<FFTW, SMOOTH>::rotate_divfree<double>(
+ typename fftw_interface<double>::complex *__restrict__ a);
+
diff --git a/cpp/kspace.hpp b/cpp/kspace.hpp
index ea82d54be7b9f915f5ec475c3c19ba4b10161acc..1377acebc1d897dee02fb1fb3de59a3480e2455b 100644
--- a/cpp/kspace.hpp
+++ b/cpp/kspace.hpp
@@ -220,7 +220,14 @@ class kspace
}
}
template <typename rnumber>
- void force_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a);
+ void project_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a);
+ // TODO: can the following be done in a cleaner way?
+ template <typename rnumber>
+ void force_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a){
+ this->template project_divfree<rnumber>(a);
+ }
+ template <typename rnumber>
+ void rotate_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a);
};
#endif//KSPACE_HPP