Merge branch 'feature/correlation-time' into develop

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)

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

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);
   }
 }
 

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);
         }
         
 

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";
             }