added feature to load field into ou process

cpp/full_code/ornstein_uhlenbeck_process.cpp

@@ -186,18 +186,54 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
 
 
 template <class rnumber, field_backend be>
-void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void)
+void ornstein_uhlenbeck_process<rnumber, be>::setup_field(int iteration,field<rnumber, be, THREE> *src)
 {
-  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(dt);
+  if (iteration == 0){
+    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(dt);
+    }
+          
+  }
+  else {
+    this->set_from_field(src);
   }
 }
 
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber, be>::set_from_field(
+  field<rnumber, be, THREE> *src)
+{
+  assert(src->real_space_representation==false);
+
+  this->kk->CLOOP_K2(
+          [&](ptrdiff_t cindex,
+            ptrdiff_t xindex,
+            ptrdiff_t yindex,
+            ptrdiff_t zindex,
+            double k2){
+
+              if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ){
+
+                for(int cc=0; cc < 3; cc++){
+                  for(int imag=0; imag < 2; imag++){
+                    this->ou_field_vort->cval(cindex,cc,imag) =  src->cval(cindex,cc,imag);
+                  }
+                }
+              }
+
+      }
+
+  );
+
+  this->calc_ou_velocity();
+
+}
+
 template <class rnumber, field_backend be>
 void ornstein_uhlenbeck_process<rnumber, be>::strip_from_field(
   field<rnumber, be, THREE> *src)
@@ -288,5 +324,41 @@ void ornstein_uhlenbeck_process<rnumber,be>::calc_ou_vorticity(void)
   this->ou_field_vort->symmetrize();
 }
 
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::calc_ou_velocity(void)
+{
+  this->kk->CLOOP_K2(
+              [&](ptrdiff_t cindex,
+                  ptrdiff_t xindex,
+                  ptrdiff_t yindex,
+                  ptrdiff_t zindex,
+                  double k2){
+    if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ)
+      {
+            this->ou_field->cval(cindex,0,0) =
+                 -(this->kk->ky[yindex]*this->ou_field_vort->cval(cindex,2,1) - this->kk->kz[zindex]*this->ou_field_vort->cval(cindex,1,1)) / k2;
+
+            this->ou_field->cval(cindex,0,1) = 
+                 (this->kk->ky[yindex]*this->ou_field_vort->cval(cindex,2,0) - this->kk->kz[zindex]*this->ou_field_vort->cval(cindex,1,0)) / k2;
+
+            this->ou_field->cval(cindex,1,0) =
+                 -(this->kk->kz[zindex]*this->ou_field_vort->cval(cindex,0,1) - this->kk->kx[xindex]*this->ou_field_vort->cval(cindex,2,1)) / k2;
+
+            this->ou_field->cval(cindex,1,1) = 
+                 (this->kk->kz[zindex]*this->ou_field_vort->cval(cindex,0,0) - this->kk->kx[xindex]*this->ou_field_vort->cval(cindex,2,0)) / k2;
+
+            this->ou_field->cval(cindex,2,0) =
+                 -(this->kk->kx[xindex]*this->ou_field_vort->cval(cindex,1,1) - this->kk->ky[yindex]*this->ou_field_vort->cval(cindex,0,1)) / k2;
+
+            this->ou_field->cval(cindex,2,1) = 
+                 (this->kk->kx[xindex]*this->ou_field_vort->cval(cindex,1,0) - this->kk->ky[yindex]*this->ou_field_vort->cval(cindex,0,0)) / k2;
+      }
+      else
+          std::fill_n((rnumber*)(this->ou_field_vort->get_cdata()+3*cindex), 6, 0.0);
+    }
+  );
+  this->ou_field->symmetrize();
+}
+
 template class ornstein_uhlenbeck_process<float,FFTW>;
 template class ornstein_uhlenbeck_process<double,FFTW>;

cpp/full_code/ornstein_uhlenbeck_process.hpp

@@ -60,7 +60,7 @@ class ornstein_uhlenbeck_process{
 
         void initialize_B(void);
 
-        void let_converge(void);
+        void setup_field(int iteration,field<rnumber, be, THREE> *src);
 
         void add_to_field_replace(
           field<rnumber,be,THREE> *src, std::string uv);
@@ -68,7 +68,11 @@ class ornstein_uhlenbeck_process{
         void strip_from_field(
           field<rnumber, be, THREE> *src);
 
+        void set_from_field(
+          field<rnumber, be, THREE> *src);
+
         void calc_ou_vorticity(void);
+        void calc_ou_velocity(void);
 
 
 };