diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 1c83e8a4bd6bb84302892fd36092cb6b60b313c7..3aea17babb57b56e02f9bdb5a3635caff925cbf5 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -725,9 +725,11 @@ class DNS(_code):
if dns_type == 'kraichnan_field':
pars['output_velocity'] = int(1)
pars['field_random_seed'] = int(1)
- pars['spectrum_slope'] = float(-5./3)
- pars['spectrum_k_cutoff'] = float(16)
- pars['spectrum_coefficient'] = float(0.1)
+ pars['spectrum_dissipation'] = float(0.027)
+ pars['spectrum_Lint'] = float(1.3)
+ pars['spectrum_etaK'] = float(0.3)
+ pars['spectrum_large_scale_const'] = float(6.78)
+ pars['spectrum_small_scale_const'] = float(0.40)
if dns_type == 'NSVE_Stokes_particles':
pars['initial_field_amplitude'] = float(0.0)
pars['initial_particle_vel'] = float(0.05)
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 110e51e90cbadde1bd6ac7eca0a4564b95c8a310..2c2124defa5e26f46e2df2afd4dc126fb0381add 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -133,9 +133,11 @@ class TEST(_code):
if dns_type == 'Gauss_field_test':
pars['histogram_bins'] = int(129)
pars['max_velocity_estimate'] = float(8)
- pars['spectrum_slope'] = float(-5./3)
- pars['spectrum_k_cutoff'] = float(16)
- pars['spectrum_coefficient'] = float(1)
+ pars['spectrum_dissipation'] = float(0.027)
+ pars['spectrum_Lint'] = float(1.3)
+ pars['spectrum_etaK'] = float(0.3)
+ pars['spectrum_large_scale_const'] = float(6.78)
+ pars['spectrum_small_scale_const'] = float(0.40)
return pars
def get_kspace(self):
kspace = {}
@@ -189,6 +191,7 @@ class TEST(_code):
if self.dns_type in ['Gauss_field_test']:
vec_spectra_stats.append('velocity')
vec_spectra_stats.append('k*velocity')
+ vec_spectra_stats.append('k2*velocity')
vec4_rspace_stats.append('velocity')
tens_rspace_stats.append('velocity_gradient')
scal_rspace_stats.append('velocity_divergence')
diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
index 00b1692e084cbb346a3a19d3f046cdba5cb1d01b..314168049fc5afad29df57e13e8c0b14c480ffd1 100644
--- a/TurTLE/test/test_Gaussian_field.py
+++ b/TurTLE/test/test_Gaussian_field.py
@@ -4,6 +4,7 @@ import numpy as np
from scipy import trapz
from scipy.stats import norm
from scipy.integrate import quad
+from scipy.optimize import root
import h5py
import sys, os
import time
@@ -18,12 +19,14 @@ except:
def main():
# size of grid
- n = 1024
- slope = -5./3.
- k_cutoff = 64.
- func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
- total_energy = quad(func, 0.6, k_cutoff*8)[0]
- coeff = 1./total_energy
+ n = 256
+ # from a 1024^3 simulation
+ dissipation = 0.37127850066981344
+ Lint = 0.9029294339535114 # probably the wrong (inconsistent with Pope) definition of Lint
+ etaK = 0.003730966576882254
+
+ c_L, c_eta = calculate_constants(epsilon=dissipation, Lint=Lint, etaK=etaK)
+
bin_no = 100
rseed = int(time.time())
simname = 'Gaussianity_test'
@@ -43,25 +46,202 @@ def main():
'--wd', './',
'--histogram_bins', str(bin_no),
'--max_velocity_estimate', '8.',
- '--spectrum_slope', str(slope),
- '--spectrum_k_cutoff', str(k_cutoff),
- '--spectrum_coefficient', str(coeff),
+ '--spectrum_dissipation', str(dissipation),
+ '--spectrum_Lint', str(Lint),
+ '--spectrum_etaK', str(etaK),
+ '--spectrum_large_scale_const', str(c_L),
+ '--spectrum_small_scale_const', str(c_eta),
'--field_random_seed', str(rseed)] +
sys.argv[1:])
- plot_stuff(simname, total_energy = total_energy)
+ plot_stuff(simname)
return None
-def plot_stuff(simname, total_energy = 1.):
+def main_premultiplied_spectra_test():
+ n_arr = [64, 128, 256, 512]
+ #n_arr = [64]
+ run_no = 3
+ dk = 1.
+ diss = 0.007036634455033392
+ energy = 0.0417604575121591
+ etaK = 0.2157937040813811
+ Lint = energy**(3./2.)/diss
+
+ c_L, c_eta = calculate_constants(epsilon=diss, Lint=Lint, etaK=etaK)
+
+ bin_no = 100
+
+ i = 1
+ if not os.path.exists('conv_test_n64_run0.h5'):
+ for n in n_arr:
+ for run in range(run_no):
+ simname = f'conv_test_n{n}_run{run}'
+
+ rseed = i
+ i += 1
+ c = TEST()
+ opt = c.launch(
+ ['Gauss_field_test',
+ '--nx', str(n),
+ '--ny', str(n),
+ '--nz', str(n),
+ '--dkx', str(dk),
+ '--dky', str(dk),
+ '--dkz', str(dk),
+ '--simname', simname,
+ '--np', '4',
+ '--environment', 'short',
+ '--minutes', '60',
+ '--ntpp', '1',
+ '--wd', './',
+ '--histogram_bins', str(bin_no),
+ '--max_velocity_estimate', '8.',
+ '--spectrum_dissipation', str(diss),
+ '--spectrum_Lint', str(Lint),
+ '--spectrum_etaK', str(etaK),
+ '--spectrum_large_scale_const', str(c_L),
+ '--spectrum_small_scale_const', str(c_eta),
+ '--field_random_seed', str(rseed)] +
+ sys.argv[1:])
+ else:
+ print('Energy = {}'.format(energy))
+
+ Q = quad(lambda k : k**2*E(k, epsilon=diss, etaK=etaK, Lint=Lint, c_L=c_L, c_eta=c_eta),
+ 0., np.inf)[0]/15.
+ print('Q = {}'.format(Q))
+ P = quad(lambda k : k**4*E(k, epsilon=diss, etaK=etaK, Lint=Lint, c_L=c_L, c_eta=c_eta),
+ 0., np.inf)[0]/105.
+ print('P = {}'.format(P))
+
+ for n in n_arr:
+ for run in range(run_no):
+ simname = f'conv_test_n{n}_run{run}'
+ df = h5py.File(simname + '.h5', 'r')
+
+ kk = df['kspace/kshell'][...]
+
+ phi_ij = df['statistics/spectra/velocity_velocity'][0]
+ energy_spec = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
+ k2spec_trace = (
+ df['statistics/spectra/k*velocity_k*velocity'][..., 0, 0]
+ + df['statistics/spectra/k*velocity_k*velocity'][..., 1, 1]
+ + df['statistics/spectra/k*velocity_k*velocity'][..., 2, 2])/2
+ k4spec_trace = (
+ df['statistics/spectra/k2*velocity_k2*velocity'][..., 0, 0]
+ + df['statistics/spectra/k2*velocity_k2*velocity'][..., 1, 1]
+ + df['statistics/spectra/k2*velocity_k2*velocity'][..., 2, 2])/2
+ print('----------------------')
+ print('n = {}, run = {}'.format(n, run))
+ print('Energy = {}'.format(np.sum(energy_spec)))
+ print('Q = {}'.format(np.sum(k2spec_trace)/15.))
+ print('P = {}'.format(np.sum(k4spec_trace)/105.))
+
+
+
+def main_deltak_test():
+ n_arr = [64, 128, 256, 512]
+ dk_arr = [1., 0.5, 0.25, 0.125]
+ diss = 0.007036634455033392
+ energy = 0.0417604575121591
+ etaK = 0.2157937040813811
+ Lint = energy**(3./2.)/diss
+
+ c_L, c_eta = calculate_constants(epsilon=diss, Lint=Lint, etaK=etaK)
+
+ bin_no = 100
+
+ if not os.path.exists('conv_test_n64_run0.h5'):
+ for n, dk in zip(n_arr, dk_arr):
+ for run in range(3):
+ simname = f'conv_test_n{n}_run{run}'
+
+ rseed = int(time.time())
+ c = TEST()
+ opt = c.launch(
+ ['Gauss_field_test',
+ '--nx', str(n),
+ '--ny', str(n),
+ '--nz', str(n),
+ '--dkx', str(dk),
+ '--dky', str(dk),
+ '--dkz', str(dk),
+ '--simname', simname,
+ '--np', '4',
+ '--environment', 'short',
+ '--minutes', '60',
+ '--ntpp', '1',
+ '--wd', './',
+ '--histogram_bins', str(bin_no),
+ '--max_velocity_estimate', '8.',
+ '--spectrum_dissipation', str(diss),
+ '--spectrum_Lint', str(Lint),
+ '--spectrum_etaK', str(etaK),
+ '--spectrum_large_scale_const', str(c_L),
+ '--spectrum_small_scale_const', str(c_eta),
+ '--field_random_seed', str(rseed)] +
+ sys.argv[1:])
+ else:
+ plt.figure()
+ k = np.geomspace(1e-1, 4e1, 1000)
+ plt.plot(k, E(k, epsilon=diss, Lint=Lint, etaK=etaK, c_L=c_L, c_eta=c_eta))
+ for n in n_arr:
+ run = 0
+ simname = f'conv_test_n{n}_run{run}'
+ df = h5py.File(simname + '.h5', 'r')
+
+ kk = df['kspace/kshell'][...]
+ phi_ij = df['statistics/spectra/velocity_velocity'][0] / df['kspace/dk'][()]
+ energy_spec = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
+ plt.scatter(kk[1:-2], energy_spec[1:-2], label=f'n={n}, r={run}')
+ plt.xscale('log')
+ plt.yscale('log')
+ plt.legend()
+ plt.savefig('spectra_convergence.pdf')
+
+
+def calculate_constants(*, epsilon, Lint, etaK):
+ sol = root(optimize_func,
+ x0 = (6.78, 0.40),
+ args = (epsilon, Lint, etaK),
+ method='hybr',
+ options={'eps':0.001, 'factor':0.1})
+ assert sol.success
+ return sol.x
+
+def optimize_func(constants, diss, Lint, etaK):
+ energy = (Lint*diss)**(2./3.)
+ energy_model = quad(
+ lambda k : E(k, epsilon=diss, Lint=Lint, etaK=etaK, c_L=constants[0], c_eta=constants[1]),
+ 0, np.inf)[0]
+ nu = (etaK**4*diss)**(1./3.)
+ diss_model = quad(
+ lambda k : 2*nu*k**2*E(k, epsilon=diss, Lint=Lint, etaK=etaK, c_L=constants[0], c_eta=constants[1]),
+ 0, np.inf)[0]
+ return (energy_model - energy, diss_model - diss)
+
+def E(k, *, epsilon, C=1.5, Lint, etaK, c_L=6.78, c_eta=0.4, beta=5.2):
+ return C*epsilon**(2./3.)*k**(-5./3.)*f_L(k*Lint, c_L=c_L)*f_eta(k*etaK, beta=beta, c_eta=c_eta)
+
+def f_L(x, *, c_L):
+ return (x/(x**2 + c_L)**(1./2.))**(11./3.)
+
+def f_eta(x, *, beta, c_eta):
+ return np.exp(-beta*((x**4 + c_eta**4)**(1./4.) - c_eta))
+
+def plot_stuff(simname):
df = h5py.File(simname + '.h5', 'r')
- for kk in ['spectrum_slope',
- 'spectrum_k_cutoff',
- 'spectrum_coefficient',
+ for kk in ['spectrum_dissipation',
+ 'spectrum_Lint',
+ 'spectrum_etaK',
+ 'spectrum_large_scale_const',
+ 'spectrum_small_scale_const',
'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'][()]
+ dissipation = df['parameters/spectrum_dissipation'][()]
+ Lint = df['parameters/spectrum_Lint'][()]
+ etaK = df['parameters/spectrum_etaK'][()]
+ c_L = df['parameters/spectrum_large_scale_const'][()]
+ c_eta = df['parameters/spectrum_small_scale_const'][()]
bin_no = df['parameters/histogram_bins'][()]
f = plt.figure()
@@ -73,7 +253,10 @@ def plot_stuff(simname, total_energy = 1.):
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.plot(kk[1:-2],
+ E(kk[1:-2],
+ epsilon=dissipation, Lint=Lint, etaK=etaK, c_eta=c_eta, c_L=c_L),
+ ls='--', c='C0')
a.set_xscale('log')
a.set_yscale('log')
# test isotropy
@@ -85,7 +268,6 @@ def plot_stuff(simname, total_energy = 1.):
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 analytically: {}'.format(total_energy))
print(np.sum(energy*np.arange(len(energy))**2))
print('Energy sum: {}'.format(np.sum(energy*df['kspace/dk'][()])))
print('Moment sum: {}'.format(df['statistics/moments/velocity'][0,2,3]/2))
@@ -108,18 +290,18 @@ def plot_stuff(simname, total_energy = 1.):
df['statistics/moments/velocity_gradient'][0, 2].mean()))
print('----------- k2-premultiplied spectrum -----------')
- k2func = lambda k, k_c=k_cutoff, s=slope : k**(2+s)*np.exp(-k/k_c)
- k2sum_analytic = quad(k2func, 0, k_cutoff*20)[0]
- print('Analytically: {}'.format(k2sum_analytic))
+ #k2func = lambda k, k_c=k_cutoff, s=slope : k**(2+s)*np.exp(-k/k_c)
+ #k2sum_analytic = quad(k2func, 0, k_cutoff*20)[0]
+ #print('Analytically: {}'.format(k2sum_analytic))
k2spec_trace = (
df['statistics/spectra/k*velocity_k*velocity'][..., 0, 0]
+ df['statistics/spectra/k*velocity_k*velocity'][..., 1, 1]
+ df['statistics/spectra/k*velocity_k*velocity'][..., 2, 2])
- print('Energy sum: {}'.format(np.sum(k2spec_trace*df['kspace/dk'][()])/2./coeff))
+ #print('Energy sum: {}'.format(np.sum(k2spec_trace*df['kspace/dk'][()])/2./coeff))
df.close()
return None
if __name__ == '__main__':
- main()
+ main_premultiplied_spectra_test()
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 7d3f0afee6dcf5b8dc0f7aa8bece40aa88a54444..2fd79ce4fe7d25f2f7f8bdae0a500ee40e0a8330 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -2120,8 +2120,11 @@ 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 dissipation,
+ const double Lint,
+ const double etaK,
+ const double c_L,
+ const double c_eta,
const double coefficient)
{
TIMEZONE("make_gaussian_random_field");
@@ -2152,30 +2155,43 @@ int make_gaussian_random_field(
const ptrdiff_t zindex,
const 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;
+ if (k2 > 0){
+ // Simple spectrum, diverges at zero
+ // double spectrum = coefficient * pow(k2, slope/2.) * exp(-sqrt(k2)/k_cutoff);
+
+ // Pope spectrum, Pope (2000), ''Turbulent flows'', p. 232
+ double C = 1.5;
+ double beta = 5.2;
+ double f_L = pow(sqrt(k2)*Lint/sqrt(k2*pow(Lint, 2) + c_L), 11./3.);
+ double f_eta = exp(-beta*(pow(pow(k2, 2)*pow(etaK, 4) + pow(c_eta, 4), 1./4.) - c_eta));
+ double spectrum = coefficient*C*pow(dissipation, 2./3.)*pow(k2, -5./6.)*f_L*f_eta;
+ // TODO: What about Delta k?
+ switch(fc)
+ {
+ case ONE:
+ {
+ output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / (4.*M_PI*k2));
+ output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / (4.*M_PI*k2));
+ break;
+ }
+ case THREE:
+ for (int cc = 0; cc<3; cc++)
+ {
+ // factor 3 compensates the trace between the spectral tensor and the energy spectrum
+ output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 3. / (4.*M_PI*k2));
+ output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 3. / (4.*M_PI*k2));
+ }
+ break;
+ case THREExTHREE:
+ for (int cc = 0; cc<3; cc++)
+ for (int ccc = 0; ccc<3; ccc++)
+ {
+ // factor 9 compensates the trace between the spectral tensor and the energy spectrum
+ output_field->cval(cindex,cc,ccc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 9. / (4.*M_PI*k2));
+ output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum * kk->dkx*kk->dky*kk->dkz / 9. / (4.*M_PI*k2));
+ }
+ break;
+ }
}
});
output_field->symmetrize();
@@ -2358,8 +2374,11 @@ template int make_gaussian_random_field<
kspace<FFTW, SMOOTH> *kk,
field<float, FFTW, ONE> *output_field,
const int rseed,
- const double slope,
- const double k_cutoff,
+ const double dissipation,
+ const double Lint,
+ const double etaK,
+ const double c_L,
+ const double c_eta,
const double coefficient);
template int make_gaussian_random_field<
double,
@@ -2369,8 +2388,11 @@ template int make_gaussian_random_field<
kspace<FFTW, SMOOTH> *kk,
field<double, FFTW, ONE> *output_field,
const int rseed,
- const double slope,
- const double k_cutoff,
+ const double dissipation,
+ const double Lint,
+ const double etaK,
+ const double c_L,
+ const double c_eta,
const double coefficient);
template int make_gaussian_random_field<
float,
@@ -2380,8 +2402,11 @@ template int make_gaussian_random_field<
kspace<FFTW, SMOOTH> *kk,
field<float, FFTW, THREE> *output_field,
const int rseed,
- const double slope,
- const double k_cutoff,
+ const double dissipation,
+ const double Lint,
+ const double etaK,
+ const double c_L,
+ const double c_eta,
const double coefficient);
template int make_gaussian_random_field<
double,
@@ -2391,7 +2416,10 @@ template int make_gaussian_random_field<
kspace<FFTW, SMOOTH> *kk,
field<double, FFTW, THREE> *output_field,
const int rseed,
- const double slope,
- const double k_cutoff,
+ const double dissipation,
+ const double Lint,
+ const double etaK,
+ const double c_L,
+ const double c_eta,
const double coefficient);
diff --git a/cpp/field.hpp b/cpp/field.hpp
index f352be510fbbf4f980cef7c5c60446b666834d8d..feb19c804766470197e37f66f8405d8b3bfb9c9f 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -431,8 +431,11 @@ 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 dissipation = 0.3,
+ const double Lint = 1.,
+ const double etaK = 0.01,
+ const double c_L = 6.78,
+ const double c_eta = 0.40,
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 17f2c2023dd858cc657328cdc56f515b06452a71..2a045dfc038a5710c9165fe0e2f892781bfa1a39 100644
--- a/cpp/full_code/Gauss_field_test.cpp
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -82,9 +82,11 @@ int Gauss_field_test<rnumber>::read_parameters()
H5F_ACC_RDONLY,
H5P_DEFAULT);
this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "/parameters/max_velocity_estimate");
- this->spectrum_slope = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_slope");
- this->spectrum_k_cutoff = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_k_cutoff");
- this->spectrum_coefficient = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_coefficient");
+ this->spectrum_dissipation = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_dissipation");
+ this->spectrum_Lint = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_Lint");
+ this->spectrum_etaK = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_etaK");
+ this->spectrum_large_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_large_scale_const");
+ this->spectrum_small_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_small_scale_const");
this->random_seed = hdf5_tools::read_value<int>(parameter_file, "/parameters/field_random_seed");
H5Fclose(parameter_file);
return EXIT_SUCCESS;
@@ -97,15 +99,19 @@ int Gauss_field_test<rnumber>::do_work(void)
/// initialize vector Gaussian random field
make_gaussian_random_field(
- this->kk,
- this->vector_field,
- this->random_seed,
- this->spectrum_slope,
- this->spectrum_k_cutoff,
- this->spectrum_coefficient);
+ this->kk,
+ this->vector_field,
+ this->random_seed,
+ this->spectrum_dissipation,
+ this->spectrum_Lint,
+ this->spectrum_etaK,
+ this->spectrum_large_scale_const,
+ this->spectrum_small_scale_const,
+ 3./2. //incompressibility projection factor
+ );
/// 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;
@@ -153,6 +159,12 @@ int Gauss_field_test<rnumber>::do_work(void)
"k*velocity_k*velocity",
0,
2.);
+ this->kk->template cospectrum<rnumber, THREE>(
+ this->vector_field->get_cdata(),
+ stat_group,
+ "k2*velocity_k2*velocity",
+ 0,
+ 4.);
/// compute basic statistics of Gaussian field
this->vector_field->compute_stats(
diff --git a/cpp/full_code/Gauss_field_test.hpp b/cpp/full_code/Gauss_field_test.hpp
index 855b19279c27da7a058b5c010ef2197308d317ef..b21df6699def682fb1ae3f62c6b14a37442cf3cb 100644
--- a/cpp/full_code/Gauss_field_test.hpp
+++ b/cpp/full_code/Gauss_field_test.hpp
@@ -46,9 +46,11 @@ class Gauss_field_test: public test
/* parameters that are read in read_parameters */
double max_velocity_estimate;
- double spectrum_slope;
- double spectrum_k_cutoff;
- double spectrum_coefficient;
+ double spectrum_dissipation;
+ double spectrum_Lint;
+ double spectrum_etaK;
+ double spectrum_large_scale_const;
+ double spectrum_small_scale_const;
int random_seed;
/* other stuff */
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index 3ec76efc03ce3c39ab2151cfaac27346e71ed950..b38956cf51e438e7a9d83fc89323cea1ee746b43 100644
--- a/cpp/full_code/kraichnan_field.cpp
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -101,9 +101,6 @@ int kraichnan_field<rnumber>::initialize(void)
"position/0");
this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
- DEBUG_MSG("Coefficient: %g\n",
- this->spectrum_coefficient);
-
this->generate_random_velocity();
// is this the first iteration?
if ((this->iteration == 0) and (this->output_velocity == 1))
@@ -127,10 +124,14 @@ int kraichnan_field<rnumber>::generate_random_velocity(void)
make_gaussian_random_field(
this->kk,
this->velocity,
- this->iteration, // not an ideal choice because resulting field sequence depends on MPI/OpenMP configuration. See note below
- this->spectrum_slope,
- this->spectrum_k_cutoff,
- this->spectrum_coefficient * 3./2.); // incompressibility projection factor
+ this->field_random_seed + 10*this->iteration,
+ this->spectrum_dissipation,
+ this->spectrum_Lint,
+ this->spectrum_etaK,
+ this->spectrum_large_scale_const,
+ this->spectrum_small_scale_const,
+ 3./2.); // incompressibility projection factor
+ // not an ideal choice because resulting field sequence depends on MPI/OpenMP configuration. See note below
// this->velocity is now in Fourier space
// project divfree, requires field in Fourier space
// Note on the choice of random seed:
@@ -297,9 +298,11 @@ int kraichnan_field<rnumber>::read_parameters(void)
this->tracers0_integration_steps = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_integration_steps");
this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
- this->spectrum_slope = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_slope");
- this->spectrum_k_cutoff = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_k_cutoff");
- this->spectrum_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_coefficient");
+ this->spectrum_dissipation = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_dissipation");
+ this->spectrum_Lint = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_Lint");
+ this->spectrum_etaK = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_etaK");
+ this->spectrum_large_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_large_scale_const");
+ this->spectrum_small_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_small_scale_const");
this->field_random_seed = hdf5_tools::read_value<int>(parameter_file, "parameters/field_random_seed");
this->output_velocity = hdf5_tools::read_value<int>(parameter_file, "parameters/output_velocity");
this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "parameters/max_velocity_estimate");
diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp
index d7fb215efcdfa57a6ac22c03717526839e1ab016..dc434c2d1f9d819e958960e22bb989b8de6063c1 100644
--- a/cpp/full_code/kraichnan_field.hpp
+++ b/cpp/full_code/kraichnan_field.hpp
@@ -60,9 +60,11 @@ class kraichnan_field: public direct_numerical_simulation
/* other stuff */
kspace<FFTW, SMOOTH> *kk;
field<rnumber, FFTW, THREE> *velocity;
- double spectrum_slope;
- double spectrum_k_cutoff;
- double spectrum_coefficient;
+ double spectrum_dissipation;
+ double spectrum_Lint;
+ double spectrum_etaK;
+ double spectrum_large_scale_const;
+ double spectrum_small_scale_const;
double max_velocity_estimate;
int field_random_seed;