diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index fb6369b2e4ce676936a5aeae4b10b5e0a151ea2a..bc7f9b84c9b960b5c7bd01cf2817c83d61911ec7 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -742,7 +742,10 @@ class DNS(_code):
opt.nz > opt.n):
opt.n = min(opt.nx, opt.ny, opt.nz)
print("Warning: '-n' parameter changed to minimum of nx, ny, nz. This affects the computation of nu.")
- self.parameters['dt'] = (opt.dtfactor / opt.n)
+ if self.dns_type in ['kraichnan_field']:
+ self.parameters['dt'] = opt.dtfactor * 0.5 / opt.n**2
+ else:
+ self.parameters['dt'] = (opt.dtfactor / opt.n)
self.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
# check value of kMax
kM = opt.n * 0.5
diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
index 3337ff956c6cc00fea6d81730197861da1a95cb1..82f5751ebb168d5160bc54c7857fddfe24a9fc7c 100644
--- a/TurTLE/test/test_Gaussian_field.py
+++ b/TurTLE/test/test_Gaussian_field.py
@@ -3,6 +3,7 @@
import numpy as np
from scipy import trapz
from scipy.stats import norm
+from scipy.integrate import quad
import h5py
import sys
import time
@@ -18,11 +19,13 @@ except:
def main():
c = TEST()
# size of grid
- n = 128
+ n = 256
slope = -5./3.
- k_cutoff = 16.
- coeff = 1.
- bin_no = 129
+ k_cutoff = 30.
+ 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
+ bin_no = 100
rseed = int(time.time())
c.launch(
@@ -35,25 +38,29 @@ def main():
'--ntpp', '1',
'--wd', './',
'--histogram_bins', str(bin_no),
- '--max_velocity_estimate', '2.',
+ '--max_velocity_estimate', '8.',
'--spectrum_slope', str(slope),
'--spectrum_k_cutoff', str(k_cutoff),
'--spectrum_coefficient', str(coeff),
'--field_random_seed', str(rseed)] +
sys.argv[1:])
+ plot_stuff(c.simname)
- df = h5py.File(c.simname + '.h5', 'r')
+def plot_stuff(simname):
+ df = h5py.File(simname + '.h5', 'r')
for kk in ['spectrum_slope',
'spectrum_k_cutoff',
'spectrum_coefficient',
'field_random_seed',
'histogram_bins']:
print(kk, df['parameters/' + kk][...])
+ coeff = df['parameters/spectrum_coefficient'][()]
+ k_cutoff = df['parameters/spectrum_k_cutoff'][()]
+ slope = df['parameters/spectrum_slope'][()]
+ bin_no = df['parameters/histogram_bins'][()]
+
f = plt.figure()
-
- # EVERYTHING SHOULD BE READ FROM FILE BECAUSE WE CANT BE SURE THE CODE RAN AGAIN
-
# test spectrum
a = f.add_subplot(121)
kk = df['kspace/kshell'][...]
@@ -74,8 +81,8 @@ def main():
hist_vel = df['statistics/histograms/velocity'][0, :, :3]
f_vel = hist_vel / np.sum(hist_vel, axis=0, keepdims=True).astype(float) / velbinsize
- print('Energy integral: {}'.format(trapz(energy[1:-2], kk[1:-2])))
- print('Energy sum: {}'.format(np.sum(energy[1:-2])))
+ print('Energy analytically: {}'.format(total_energy))
+ print('Energy sum: {}'.format(np.sum(energy*df['kspace/dk'][()])))
print('Moment sum: {}'.format(df['statistics/moments/velocity'][0,2,3]/2))
print('Velocity variances: {}'.format(trapz(vel[:,None]**2*f_vel, vel[:,None], axis=0)))
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
index 2b6926da80e9d5036208bc0bc8b6d7438338103e..7c2a2f746685be066df8c4e79a3a8c3bd9f5b696 100644
--- a/cpp/full_code/Gauss_field_test.cpp
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -73,27 +73,26 @@ int make_gaussian_random_field(
{
if (k2 > 0)
switch(fc)
- {
+ {
case ONE:
{
- output_field->cval(cindex,0) = coefficient * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
- output_field->cval(cindex,1) = coefficient * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
+ 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++)
- { // In the slope, a factor of 1/2 compensates the k^2, a factor of 1/2 goes from the energy to the Fourier coefficient
- // and the 2*pi*k^2 divides out the surface of the shell
- output_field->cval(cindex,cc,0) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, (slope - 2)/4.) * exp(- sqrt(k2)/k_cutoff/2.);
- output_field->cval(cindex,cc,1) = coefficient / 2 / M_PI * rdist(rgen[omp_get_thread_num()]) * pow(k2, (slope - 2)/4.) * exp(- sqrt(k2)/k_cutoff/2.);
- } // TODO: ADJUST OTHER FIELD TYPES
+ {
+ 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) = coefficient * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
- output_field->cval(cindex,cc,ccc,1) = coefficient * rdist(rgen[omp_get_thread_num()]) * pow(k2, slope/4.) * exp(- sqrt(k2)/k_cutoff/2.);
+ 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;
}
@@ -174,7 +173,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(), true);
+ this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata());
/// initialize statistics file.
hid_t stat_group, stat_file;
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index 61192d7ed35b7f6ece154b019cd3814178e9374e..fb65329b776d117a036cd71fa133a2d2e3463cfa 100644
--- a/cpp/full_code/kraichnan_field.cpp
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -102,20 +102,14 @@ int kraichnan_field<rnumber>::initialize(void)
"position/0");
this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
- this->slope = -5./3.;
- this->k_cutoff = 20.; //k_max is approximately N/2 in a N^3 simulation
- const double energy = 1.;
+ this->spectrum_slope = -5./3.;
+ this->spectrum_k_cutoff = 20.; //k_max is approximately N/2 in a N^3 simulation
- // Here we rescale the field to match the desired energy. The divisor
- // is the integral of the power law spectrum with exponential cutoff.
- // The prefactor sqrt(3/2) compensates the projection to incompressibity.
- this->field_coefficient = sqrt(3./2.* energy /
- (1.10369));
- // 0.93058 for k_cutoff = 10
- // 1.10369 for k_cutoff = 20
+ // Projection factor
+ this->spectrum_coeff = 3./2.;
DEBUG_MSG("Coefficient: %g\n",
- this->field_coefficient);
+ this->spectrum_coeff);
return EXIT_SUCCESS;
}
@@ -129,32 +123,19 @@ int kraichnan_field<rnumber>::step(void)
this->kk,
this->velocity,
this->iteration,
- this->slope,
- this->k_cutoff,
- this->field_coefficient);
+ this->spectrum_slope,
+ this->spectrum_k_cutoff,
+ this->spectrum_coeff * 3./2.); // incompressibility projection factor
this->velocity->ift();
- // PUT THE INCOMPRESSIBILITY PROJECTION FACTOR sqrt(2/3)
+ // 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: %g\n",
+ DEBUG_MSG("L2Norm after: %g\n",
this->velocity->L2norm(this->kk));
- // What does the ->template do?
- // ADJUST TOTAL ENERGY, SLOPE AND CUTOFF
- // we actually want to propagate with sqrt(dt) instead of dt
-
- DEBUG_MSG("Iteration: %d\n",
- this->iteration);
- DEBUG_MSG("Field entries: %g\n",
- this->velocity->rval(0,0));
- DEBUG_MSG("Field entries: %g\n",
- this->velocity->rval(0,1));
- DEBUG_MSG("Field entries: %g\n",
- this->velocity->rval(0,2));
- DEBUG_MSG("Field entries: %g\n",
- this->velocity->rval(15,2));
- // In principle this should be sqrt(dt), but to match the NS energy,
- // we absorb this into the velocity field.
- this->ps->completeLoop(this->dt);
+
+ this->ps->completeLoop(sqrt(this->dt));
this->iteration++;
return EXIT_SUCCESS;
}
diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp
index a6be715f6ba6bdcec9e178207478db516eb6c149..e9db021bf72da8b6ade9cfeef266717c23b07f20 100644
--- a/cpp/full_code/kraichnan_field.hpp
+++ b/cpp/full_code/kraichnan_field.hpp
@@ -59,9 +59,9 @@ class kraichnan_field: public direct_numerical_simulation
/* other stuff */
kspace<FFTW, SMOOTH> *kk;
field<rnumber, FFTW, THREE> *velocity;
- double slope;
- double k_cutoff;
- double field_coefficient;
+ double spectrum_slope;
+ double spectrum_k_cutoff;
+ double spectrum_coeff;
/* other stuff */
std::unique_ptr<abstract_particles_system<long long int, double>> ps;