diff --git a/cpp/full_code/NSVE_field_stats.cpp b/cpp/full_code/NSVE_field_stats.cpp
index 0969175cc75530e2dad2c3c5dd9e6a0449416ed0..7a4f1312ff15e291321058dbcc0907b1d681f374 100644
--- a/cpp/full_code/NSVE_field_stats.cpp
+++ b/cpp/full_code/NSVE_field_stats.cpp
@@ -101,6 +101,33 @@ int NSVE_field_stats<rnumber>::read_current_cvorticity(void)
}
return EXIT_SUCCESS;
}
+template <typename rnumber>
+int NSVE_field_stats<rnumber>::read_arbitrary_cvorticity(int arbitrary_iteration)
+{
+ TIMEZONE("NSVE_field_stats::read_arbitrary_cvorticity");
+ this->vorticity->real_space_representation = false;
+ if (this->bin_IO != NULL)
+ {
+ char itername[16];
+ sprintf(itername, "i%.5x", arbitrary_iteration);
+ std::string native_binary_fname = (
+ this->simname +
+ std::string("_cvorticity_") +
+ std::string(itername));
+ this->bin_IO->read(
+ native_binary_fname,
+ this->vorticity->get_cdata());
+ }
+ else
+ {
+ this->vorticity->io(
+ this->simname + std::string("_fields.h5"),
+ "vorticity",
+ arbitrary_iteration,
+ true);
+ }
+ return EXIT_SUCCESS;
+}
template <typename rnumber>
int NSVE_field_stats<rnumber>::finalize(void)
diff --git a/cpp/full_code/NSVE_field_stats.hpp b/cpp/full_code/NSVE_field_stats.hpp
index 28a2376f17ac2ac837cbacac828cd91572bb3a17..0939ba03b60e03f66073d1713d5c3905d9110dd6 100644
--- a/cpp/full_code/NSVE_field_stats.hpp
+++ b/cpp/full_code/NSVE_field_stats.hpp
@@ -59,6 +59,7 @@ class NSVE_field_stats: public postprocess
virtual int finalize(void);
int read_current_cvorticity(void);
+ int read_arbitrary_cvorticity(int iter);
};
#endif//NSVE_FIELD_STATS_HPP
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
index 2c680914c3b13042ec50f2f33e0743ba677cced6..653c2540f01ee18661346b04945c6cade09450d2 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -4,7 +4,7 @@
#include <cassert>
#include "scope_timer.hpp"
#include <algorithm>
-
+#include <chrono>
template <class rnumber,field_backend be>
@@ -16,6 +16,7 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
double ou_kmin,
double ou_kmax,
double ou_energy_amplitude,
+ double ou_gamma_factor,
double DKX,
double DKY,
double DKZ,
@@ -28,7 +29,7 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
this->ou_field = new field<rnumber,be,THREE>(
nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
*this->ou_field = 0.0;
- this->ou_field->dft();
+ this->ou_field->dft();
this->ou_field_vort = new field<rnumber,be,THREE>(
nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
@@ -44,21 +45,29 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
this->ou_kmin_squ = pow(ou_kmin,2);
this->ou_kmax_squ = pow(ou_kmax,2);
this->ou_energy_amplitude = ou_energy_amplitude;
+ this->ou_gamma_factor = ou_gamma_factor;
this->epsilon = pow((this->ou_energy_amplitude/this->kolmogorov_constant), 3./2.);
- assert(this->kk->kM2 >= this->ou_kmax_squ);
+ assert(this->kk->kM2 >= this->ou_kmax_squ);
gen.resize(omp_get_max_threads());
+ long now;
+
+ if (this->ou_field->clayout->myrank == 0){
+ now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
+ }
+ MPI_Bcast(&now,1,MPI_LONG,0,this->ou_field->comm);
+
for(int thread=0;thread<omp_get_max_threads();thread++)
{
- int current_seed =
+ long current_seed =
this->ou_field->clayout->myrank*omp_get_max_threads() +
- thread;
+ thread+now;
gen[thread].seed(current_seed);
}
-
+
this->initialize_B();
}
@@ -179,10 +188,13 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void)
{
- //add some logic here TODO
- for(int i=0; i<1000; i++)
+ double ou_kmin = sqrt(this->ou_kmin_squ);
+ double tau = 1.0/this->gamma(ou_kmin);
+ double dt = tau/1000.;
+
+ for(int i=0; i<2000; i++)
{
- this->step(2e-3);
+ this->step(dt);
}
}
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp
index a1aeaa98dd948ca39b2c96177341890b345c9953..5a1288c929afac4cad4a3a6f36ff5b2a0436b1ca 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.hpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp
@@ -19,6 +19,7 @@ class ornstein_uhlenbeck_process{
double ou_kmin_squ;
double ou_kmax_squ;
double ou_energy_amplitude;
+ double ou_gamma_factor;
double kolmogorov_constant = 2;
double epsilon;
@@ -38,6 +39,7 @@ class ornstein_uhlenbeck_process{
double ou_kmin,
double ou_kmax,
double ou_energy_amplitude,
+ double ou_gamma_factor,
double DKX = 1.0,
double DKY = 1.0,
double DKZ = 1.0,
@@ -52,8 +54,7 @@ class ornstein_uhlenbeck_process{
inline double gamma(double kabs)
{
- return this->kolmogorov_constant*sqrt(this->kolmogorov_constant/
- this->ou_energy_amplitude);
+ return this->ou_gamma_factor*pow(kabs,2./3.)*sqrt(this->ou_energy_amplitude);
}
diff --git a/cpp/full_code/ou_vorticity_equation.hpp b/cpp/full_code/ou_vorticity_equation.hpp
index e36292c2271616794815d43d084810b9daeb3fc9..8f8e110c64ca29bd1a5039a96f7cc8e632b5fc7b 100644
--- a/cpp/full_code/ou_vorticity_equation.hpp
+++ b/cpp/full_code/ou_vorticity_equation.hpp
@@ -24,6 +24,7 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
double ou_kmin,
double ou_kmax,
double ou_energy_amplitude,
+ double ou_gamma_factor,
double DKX = 1.0,
double DKY = 1.0,
double DKZ = 1.0,
@@ -34,7 +35,7 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
this->ou = new ornstein_uhlenbeck_process<rnumber,be>(
NAME,
nx, ny, nz,
- ou_kmin, ou_kmax, ou_energy_amplitude,
+ ou_kmin, ou_kmax, ou_energy_amplitude,ou_gamma_factor,
DKX, DKY, DKZ, FFTW_PLAN_RIGOR);
this->ou_forcing_type = "replace";
}