diff --git a/cpp/field.cpp b/cpp/field.cpp
index 7d3f0afee6dcf5b8dc0f7aa8bece40aa88a54444..2f95949f0dd12b3c9a17f6f94fd4c289babd648e 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -2119,10 +2119,13 @@ template <typename rnumber,
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)
+ const double dissipation,
+ const double Lint,
+ const double etaK,
+ const double c_L,
+ const double c_eta,
+ const double coefficient,
+ const int rseed)
{
TIMEZONE("make_gaussian_random_field");
// initialize a separate random number generator for each thread
@@ -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 / (4.*M_PI*k2));
+ output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum / (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 / 3. / (4.*M_PI*k2));
+ output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum / 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 / 9. / (4.*M_PI*k2));
+ output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(spectrum / 9. / (4.*M_PI*k2));
+ }
+ break;
+ }
}
});
output_field->symmetrize();
diff --git a/cpp/field.hpp b/cpp/field.hpp
index f352be510fbbf4f980cef7c5c60446b666834d8d..3c109b29203c0bbd3d6af87313192bc50c3eb406 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -430,10 +430,13 @@ template <typename rnumber,
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.);
+ const double dissipation,
+ const double Lint,
+ const double etaK,
+ const double c_L = 6.78,
+ const double c_eta = 0.40,
+ const double coefficient = 1.,
+ const int rseed = 0);
#endif//FIELD_HPP
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index 3ec76efc03ce3c39ab2151cfaac27346e71ed950..2c9aa904775db6e3685af471e540df894c89e578 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,13 @@ 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->spectrum_dissipation,
+ this->spectrum_Lint,
+ this->spectrum_etaK,
+ this->spectrum_c_L,
+ this->spectrum_c_eta,
+ 3./2., // incompressibility projection factor
+ this->field_random_seed + 10*this->iteration); // 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 +297,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_c_L = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_c_L");
+ this->spectrum_c_eta = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_c_eta");
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");