added vorticity field to ou process

cpp/full_code/ornstein_uhlenbeck_process.cpp

@@ -28,6 +28,12 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
         nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
     *this->ou_field = 0.0;
     this->ou_field->dft(); //have it complex all zero TODO:Ask Chichi
+
+    this->ou_field_vort = new field<rnumber,be,THREE>(
+        nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
+    *this->ou_field_vort = 0.0;
+    this->ou_field_vort->dft();
+
     this->B = new field<rnumber,be,THREExTHREE>(
         nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
 
@@ -48,6 +54,7 @@ ornstein_uhlenbeck_process<rnumber,be>::~ornstein_uhlenbeck_process()
     TIMEZONE("ornstein_uhlenbeck_process::~ornstein_uhlenbeck_process");
     delete this->kk;
     delete this->ou_field;
+    delete this->ou_field_vort;
     delete this->B;
 }
 
@@ -81,6 +88,9 @@ void ornstein_uhlenbeck_process<rnumber,be>::step(
 
         }
     });
+
+    this->calc_ou_vorticity();
+
     // this->ou_field->symmetrize(); ???
 }
 
@@ -244,6 +254,30 @@ void ornstein_uhlenbeck_process<rnumber,be>::add_to_field_gaussian(
   );
 }
 
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::calc_ou_vorticity(void)
+{
+  this->kk->CLOOP_K2(
+              [&](ptrdiff_t cindex,
+                  ptrdiff_t xindex,
+                  ptrdiff_t yindex,
+                  ptrdiff_t zindex,
+                  double k2){
+      if (k2 <= this->kk->kM2)
+      {
+          this->ou_field_vort->cval(cindex,0,0) = -(this->kk->ky[yindex]*this->ou_field->cval(cindex,2,1) - this->kk->kz[zindex]*this->ou_field->cval(cindex,1,1));
+          this->ou_field_vort->cval(cindex,0,1) =  (this->kk->ky[yindex]*this->ou_field->cval(cindex,2,0) - this->kk->kz[zindex]*this->ou_field->cval(cindex,1,0));
+          this->ou_field_vort->cval(cindex,1,0) = -(this->kk->kz[zindex]*this->ou_field->cval(cindex,0,1) - this->kk->kx[xindex]*this->ou_field->cval(cindex,2,1));
+          this->ou_field_vort->cval(cindex,1,1) =  (this->kk->kz[zindex]*this->ou_field->cval(cindex,0,0) - this->kk->kx[xindex]*this->ou_field->cval(cindex,2,0));
+          this->ou_field_vort->cval(cindex,2,0) = -(this->kk->kx[xindex]*this->ou_field->cval(cindex,1,1) - this->kk->ky[yindex]*this->ou_field->cval(cindex,0,1));
+          this->ou_field_vort->cval(cindex,2,1) =  (this->kk->kx[xindex]*this->ou_field->cval(cindex,1,0) - this->kk->ky[yindex]*this->ou_field->cval(cindex,0,0));
+      }
+      else
+          std::fill_n((rnumber*)(this->ou_field_vort->get_cdata()+3*cindex), 6, 0.0);
+  }
+  );
+  this->ou_field_vort->symmetrize();
+}
 
 template class ornstein_uhlenbeck_process<float,FFTW>;
 template class ornstein_uhlenbeck_process<double,FFTW>;

cpp/full_code/ornstein_uhlenbeck_process.hpp

@@ -22,6 +22,7 @@ class ornstein_uhlenbeck_process{
         double ou_energy_amplitude;
 
         field<rnumber,be,THREE> *ou_field;
+        field<rnumber,be,THREE> *ou_field_vort;
         field<rnumber,be,THREExTHREE> *B;
         // field<rnumber,be,ONE> *number_of_modes;
         kspace<be,SMOOTH> *kk;
@@ -61,6 +62,8 @@ class ornstein_uhlenbeck_process{
         void add_to_field_replace(
           field<rnumber,be,THREE> *src);
 
+        void calc_ou_vorticity(void);
+
 
 };