replace vorticity modes as well

cpp/full_code/ou_vorticity_equation.cpp

@@ -12,6 +12,85 @@ ou_vorticity_equation<rnumber,be>::~ou_vorticity_equation()
   delete this->ou;
 }
 
+template <class rnumber,
+          field_backend be>
+void ou_vorticity_equation<rnumber, be>::step(double dt)
+{
+    DEBUG_MSG("vorticity_equation::step\n");
+    TIMEZONE("vorticity_equation::step");
+    this->ou->add_to_field_replace(this->cvorticity, "vorticity");
+    *this->v[1] = 0.0;
+    this->omega_nonlin(0);
+    this->kk->CLOOP_K2(
+                [&](ptrdiff_t cindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex,
+                    double k2){
+        if (k2 <= this->kk->kM2)
+        {
+            double factor0;
+            factor0 = exp(-this->nu * k2 * dt);
+            for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+                this->v[1]->cval(cindex,cc,i) = (
+                        this->v[0]->cval(cindex,cc,i) +
+                        dt*this->u->cval(cindex,cc,i))*factor0;
+        }
+    }
+    );
+    this->impose_forcing(this->v[1], this->v[0]);
+
+    this->omega_nonlin(1);
+    this->kk->CLOOP_K2(
+                [&](ptrdiff_t cindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex,
+                    double k2){
+        if (k2 <= this->kk->kM2)
+        {
+            double factor0, factor1;
+            factor0 = exp(-this->nu * k2 * dt/2);
+            factor1 = exp( this->nu * k2 * dt/2);
+            for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+                this->v[2]->cval(cindex, cc, i) = (
+                        3*this->v[0]->cval(cindex,cc,i)*factor0 +
+                        ( this->v[1]->cval(cindex,cc,i) +
+                         dt*this->u->cval(cindex,cc,i))*factor1)*0.25;
+        }
+    }
+    );
+    this->impose_forcing(this->v[2], this->v[0]);
+
+    this->omega_nonlin(2);
+    // store old vorticity
+    *this->v[1] = *this->v[0];
+    this->kk->CLOOP_K2(
+                [&](ptrdiff_t cindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex,
+                    double k2){
+        if (k2 <= this->kk->kM2)
+        {
+            double factor0;
+            factor0 = exp(-this->nu * k2 * dt * 0.5);
+            for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+                this->v[3]->cval(cindex,cc,i) = (
+                        this->v[0]->cval(cindex,cc,i)*factor0 +
+                        2*(this->v[2]->cval(cindex,cc,i) +
+                           dt*this->u->cval(cindex,cc,i)))*factor0/3;
+        }
+    }
+    );
+    this->impose_forcing(this->v[0], this->v[1]);
+
+    this->kk->template force_divfree<rnumber>(this->cvorticity->get_cdata());
+    this->cvorticity->symmetrize();
+    this->iteration++;
+}
+
+
 template <class rnumber, field_backend be>
 void ou_vorticity_equation<rnumber, be>::omega_nonlin(int src)
 {

cpp/full_code/ou_vorticity_equation.hpp

@@ -42,6 +42,7 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
         ~ou_vorticity_equation();
 
         void omega_nonlin(int src);
+        void step(double dt);
 
         void add_ou_forcing_velocity(void);
         void add_ou_forcing_vorticity(void);

cpp/vorticity_equation.hpp

@@ -88,7 +88,7 @@ class vorticity_equation
 
         /* solver essential methods */
         virtual void omega_nonlin(int src);
-        void step(double dt);
+        virtual void step(double dt);
         void impose_zero_modes(void);
 
         /** \brief Method that computes force and adds it to the right hand side of the NS equations.