wip move to turtle addons

cpp/full_code/NSVE_ou_forcing.cpp

+#include <string>
+#include "NSVE_ou_forcing.hpp"
+#include "fftw_tools.hpp"
+
+
+
+template <typename rnumber>
+int NSVE_ou_forcing<rnumber>::read_parameters(void)
+{
+  this->NSVE<rnumber>::read_parameters();
+  return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_ou_forcing<rnumber>::initialize(void)
+{
+  this->NSVE<rnumber>::initialize();
+  this->fs = new ou_vorticity_equation<rnumber,FFTW>(
+                this->simname.c_str(),
+                this->nx, this->ny, this->nz,
+                this->dkx, this->dky, this->dkz,
+                this->ou_kmin, this->ou_kmax, this->ou_energy_amplitude,
+                fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+
+  return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_ou_forcing<rnumber>::step(void){
+  this->NSVE<rnumber>::step();
+  //TODO STEP OF OU????
+  return EXIT_SUCCESS;
+}
+
+template class NSVE_ou_forcing<float>;
+template class NSVE_ou_forcing<double>;

cpp/full_code/NSVE_ou_forcing.hpp

+#ifndef NSVE_OU_FORCING_HPP
+#define NSVE_OU_FORCING_HPP
+
+#include <cstdlib>
+#include "base.hpp"
+#include "ou_vorticity_equation.hpp"
+#include "NSVE.hpp"
+
+template <typename rnumber>
+class NSVE_ou_forcing: public NSVE<rnumber>
+{
+    public:
+
+      double ou_kmax;
+      double ou_kmin;
+      double ou_energy_amplitude;
+
+      ou_vorticity_equation<rnumber,FFTW> *fs;
+
+
+      NSVE_ou_forcing(
+              const MPI_Comm COMMUNICATOR,
+              const std::string &simulation_name):
+          NSVE<rnumber>(
+                  COMMUNICATOR,
+                  simulation_name),
+          ou_kmax(5.0), ou_kmin(2.0), ou_energy_amplitude(1.0){}
+
+      ~NSVE_ou_forcing();
+
+      int read_parameters(void);
+      int step(void);
+      int initialize(void);
+
+};
+
+#endif

cpp/full_code/ou_vorticity_equation.cpp

@@ -11,3 +11,68 @@ ou_vorticity_equation<rnumber,be>::~ou_vorticity_equation()
 {
   delete this->ou;
 }
+
+template <class rnumber, field_backend be>
+void ou_vorticity_equation<rnumber, be>::omega_nonlin(int src)
+{
+  DEBUG_MSG("ou_vorticity_equation::omega_nonlin(%d)\n", src);
+  assert(src >= 0 && src < 3);
+  this->compute_velocity(this->v[src]);
+
+  this->add_ou_forcing();
+  // TIME STEP?!?!?!!?! TODO
+
+  /* get fields from Fourier space to real space */
+  this->u->ift();
+  this->rvorticity->real_space_representation = false;
+  *this->rvorticity = this->v[src]->get_cdata();
+  this->rvorticity->ift();
+  /* compute cross product $u \times \omega$, and normalize */
+  this->u->RLOOP(
+              [&](ptrdiff_t rindex,
+                  ptrdiff_t xindex,
+                  ptrdiff_t yindex,
+                  ptrdiff_t zindex){
+      rnumber tmp[3];
+      for (int cc=0; cc<3; cc++)
+          tmp[cc] = (this->u->rval(rindex,(cc+1)%3)*this->rvorticity->rval(rindex,(cc+2)%3) -
+                     this->u->rval(rindex,(cc+2)%3)*this->rvorticity->rval(rindex,(cc+1)%3));
+      for (int cc=0; cc<3; cc++)
+          this->u->rval(rindex,cc) = tmp[cc] / this->u->npoints;
+  }
+  );
+  /* go back to Fourier space */
+  //this->clean_up_real_space(this->ru, 3);
+  this->u->dft();
+  this->kk->template dealias<rnumber, THREE>(this->u->get_cdata());
+  /* $\imath k \times Fourier(u \times \omega)$ */
+  this->kk->CLOOP(
+              [&](ptrdiff_t cindex,
+                  ptrdiff_t xindex,
+                  ptrdiff_t yindex,
+                  ptrdiff_t zindex){
+      rnumber tmp[3][2];
+      {
+          tmp[0][0] = -(this->kk->ky[yindex]*this->u->cval(cindex,2,1) - this->kk->kz[zindex]*this->u->cval(cindex,1,1));
+          tmp[1][0] = -(this->kk->kz[zindex]*this->u->cval(cindex,0,1) - this->kk->kx[xindex]*this->u->cval(cindex,2,1));
+          tmp[2][0] = -(this->kk->kx[xindex]*this->u->cval(cindex,1,1) - this->kk->ky[yindex]*this->u->cval(cindex,0,1));
+          tmp[0][1] =  (this->kk->ky[yindex]*this->u->cval(cindex,2,0) - this->kk->kz[zindex]*this->u->cval(cindex,1,0));
+          tmp[1][1] =  (this->kk->kz[zindex]*this->u->cval(cindex,0,0) - this->kk->kx[xindex]*this->u->cval(cindex,2,0));
+          tmp[2][1] =  (this->kk->kx[xindex]*this->u->cval(cindex,1,0) - this->kk->ky[yindex]*this->u->cval(cindex,0,0));
+      }
+      for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+          this->u->cval(cindex, cc, i) = tmp[cc][i];
+  }
+  );
+  //this->add_forcing(this->u, this->v[src]); TODO not necessary for ou?
+  this->kk->template force_divfree<rnumber>(this->u->get_cdata());
+  this->u->symmetrize();
+}
+
+template <class rnumber, field_backend be>
+void ou_vorticity_equation<rnumber, be>::add_ou_forcing()
+{
+  if (this->ou_forcing_type == "replace"){
+    this->ou->add_to_field_replace(this->u);
+  }
+}

cpp/full_code/ou_vorticity_equation.hpp

@@ -13,6 +13,7 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
 {
     public:
 
+      std::string ou_forcing_type;
 
       ornstein_uhlenbeck_process<rnumber,be> *ou;
       ou_vorticity_equation(
@@ -35,10 +36,15 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
                 nx, ny, nz,
                 ou_kmin, ou_kmax, ou_energy_amplitude,
                 DKX, DKY, DKZ, FFTW_PLAN_RIGOR);
+              this->ou_forcing_type = "replace";
             }
 
         ~ou_vorticity_equation();
 
+        void omega_nonlin(int src);
+
+        void add_ou_forcing(void);
+
 };