diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
index 82f5751ebb168d5160bc54c7857fddfe24a9fc7c..5fdc3506f87f5fc7e964bb779007d303cbfd0ef7 100644
--- a/TurTLE/test/test_Gaussian_field.py
+++ b/TurTLE/test/test_Gaussian_field.py
@@ -22,9 +22,9 @@ def main():
n = 256
slope = -5./3.
k_cutoff = 30.
- func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
+ func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
total_energy = quad(func, 1, k_cutoff*4)[0]
- coeff = 3./2./total_energy
+ coeff = 1./total_energy
bin_no = 100
rseed = int(time.time())
@@ -44,10 +44,10 @@ def main():
'--spectrum_coefficient', str(coeff),
'--field_random_seed', str(rseed)] +
sys.argv[1:])
+ plot_stuff(c.simname, total_energy = total_energy)
+ return None
- plot_stuff(c.simname)
-
-def plot_stuff(simname):
+def plot_stuff(simname, total_energy = 1.):
df = h5py.File(simname + '.h5', 'r')
for kk in ['spectrum_slope',
'spectrum_k_cutoff',
@@ -59,7 +59,7 @@ def plot_stuff(simname):
k_cutoff = df['parameters/spectrum_k_cutoff'][()]
slope = df['parameters/spectrum_slope'][()]
bin_no = df['parameters/histogram_bins'][()]
-
+
f = plt.figure()
# test spectrum
a = f.add_subplot(121)
@@ -67,7 +67,7 @@ def plot_stuff(simname):
print('dk: {}'.format(df['kspace/dk'][()]))
phi_ij = df['statistics/spectra/velocity_velocity'][0] / df['kspace/dk'][()]
energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
-
+
a.scatter(kk[1:-2], energy[1:-2])
a.plot(kk[1:-2], coeff*kk[1:-2]**slope*np.exp(-kk[1:-2]/k_cutoff), ls='--', c='C0')
a.set_xscale('log')
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 657cc0cb3d7d9c630ac7e0a949c9cde6dae4ce96..fdd8c2b4eb748021de91c016a4b48957f9432cc6 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -29,6 +29,7 @@
#include <cstdlib>
#include <algorithm>
#include <cassert>
+#include <random>
#include "field.hpp"
#include "scope_timer.hpp"
#include "shared_array.hpp"
@@ -2032,6 +2033,76 @@ field<rnumber, be, fc> &field<rnumber, be, fc>::operator=(
return *this;
}
+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::normal_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;
+ DEBUG_MSG("slope: %g\n", slope);
+ // 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:
+ {
+ output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ break;
+ }
+ case THREE:
+ for (int cc = 0; cc<3; cc++)
+ {
+ output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ }
+ 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[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ }
+ break;
+ }
+ });
+ output_field->symmetrize();
+
+ return EXIT_SUCCESS;
+}
+
template class field<float, FFTW, ONE>;
template class field<float, FFTW, THREE>;
template class field<float, FFTW, THREExTHREE>;
@@ -2189,3 +2260,49 @@ template int joint_rspace_3PDF<double, FFTW>(
const std::vector<double>,
const std::vector<double>);
+
+template int make_gaussian_random_field<
+ float,
+ FFTW,
+ ONE,
+ SMOOTH>(
+ kspace<FFTW, SMOOTH> *kk,
+ field<float, FFTW, ONE> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient);
+template int make_gaussian_random_field<
+ double,
+ FFTW,
+ ONE,
+ SMOOTH>(
+ kspace<FFTW, SMOOTH> *kk,
+ field<double, FFTW, ONE> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient);
+template int make_gaussian_random_field<
+ float,
+ FFTW,
+ THREE,
+ SMOOTH>(
+ kspace<FFTW, SMOOTH> *kk,
+ field<float, FFTW, THREE> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient);
+template int make_gaussian_random_field<
+ double,
+ FFTW,
+ THREE,
+ SMOOTH>(
+ kspace<FFTW, SMOOTH> *kk,
+ field<double, FFTW, THREE> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient);
+
diff --git a/cpp/field.hpp b/cpp/field.hpp
index af2136b193c4b22d65972fac287c31ac7d20ec38..0f74d66c0b548e9ef496bb60bae1b21171db937c 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -374,5 +374,20 @@ int joint_rspace_3PDF(
const std::vector<double> max_f2_estimate,
const std::vector<double> max_f3_estimate);
+/** \brief Generate a Gaussian random field
+ *
+ * */
+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,
+ const double slope = -5./3.,
+ const double k2_cutoff = 10.,
+ const double coefficient = 1.);
+
#endif//FIELD_HPP
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
index 7c2a2f746685be066df8c4e79a3a8c3bd9f5b696..6251d558efec714385c134ebde681d940f7512eb 100644
--- a/cpp/full_code/Gauss_field_test.cpp
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -32,75 +32,7 @@
-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::normal_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;
- DEBUG_MSG("slope: %g\n", slope);
- // 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:
- {
- output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
- output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
- break;
- }
- case THREE:
- for (int cc = 0; cc<3; cc++)
- {
- output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
- output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
- }
- 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[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
- output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
- }
- break;
- }
- });
- output_field->symmetrize();
- return EXIT_SUCCESS;
-}
template <typename rnumber>
int Gauss_field_test<rnumber>::initialize(void)
@@ -173,7 +105,7 @@ int Gauss_field_test<rnumber>::do_work(void)
this->spectrum_coefficient);
/// impose divergence free condition while maintaining the energy of the field
- this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
+ this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata(), true);
/// initialize statistics file.
hid_t stat_group, stat_file;
diff --git a/cpp/full_code/Gauss_field_test.hpp b/cpp/full_code/Gauss_field_test.hpp
index 68b9891809cf0cd016707047ca492827432a5a2c..855b19279c27da7a058b5c010ef2197308d317ef 100644
--- a/cpp/full_code/Gauss_field_test.hpp
+++ b/cpp/full_code/Gauss_field_test.hpp
@@ -39,20 +39,6 @@
*
*/
-
-
-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,
- const double slope = -5./3.,
- const double k2_cutoff = 10.,
- const double coefficient = 1.);
-
template <typename rnumber>
class Gauss_field_test: public test
{
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index fb65329b776d117a036cd71fa133a2d2e3463cfa..00ba67f10a24fd738fd1f4cce64814749b280fec 100644
--- a/cpp/full_code/kraichnan_field.cpp
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -28,7 +28,6 @@
#include "kraichnan_field.hpp"
#include "scope_timer.hpp"
#include "fftw_tools.hpp"
-#include "Gauss_field_test.hpp"
template <typename rnumber>
@@ -65,12 +64,12 @@ int kraichnan_field<rnumber>::initialize(void)
this->velocity = new field<rnumber, FFTW, THREE>(
nx, ny, nz,
this->comm,
- fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->velocity->real_space_representation = false;
+ // the velocity layout should be the same as the vorticity one
this->kk = new kspace<FFTW, SMOOTH>(
this->velocity->clayout, this->dkx, this->dky, this->dkz);
- // IS THE VELOCITY LAYOUT THE SAME AS THE VORTICITY ONE?
// initialize particles
this->ps = particles_system_builder(
@@ -128,13 +127,13 @@ int kraichnan_field<rnumber>::step(void)
this->spectrum_coeff * 3./2.); // incompressibility projection factor
this->velocity->ift();
- // PUT THE INCOMPRESSIBILITY PROJECTION FACTOR sqrt(2/3)
- DEBUG_MSG("L2Norm before: %g\n",
+ // PUT THE INCOMPRESSIBILITY PROJECTION FACTOR sqrt(2/3)
+ DEBUG_MSG("L2Norm before: %g\n",
this->velocity->L2norm(this->kk));
this->kk->template project_divfree<rnumber>(this->velocity->get_cdata());
- DEBUG_MSG("L2Norm after: %g\n",
+ DEBUG_MSG("L2Norm after: %g\n",
this->velocity->L2norm(this->kk));
-
+
this->ps->completeLoop(sqrt(this->dt));
this->iteration++;
return EXIT_SUCCESS;