add back symmetrize and friends

src/fluid_solver.cpp

@@ -173,7 +173,7 @@ fluid_solver<R>::~fluid_solver() \
 template<> \
 void fluid_solver<R>::compute_vorticity() \
 { \
-    /*this->low_pass_Fourier(this->cu, 3, this->kM);*/ \
+    this->low_pass_Fourier(this->cu, 3, this->kM); \
     CLOOP( \
             this->cvorticity[3*cindex+0][0] = -(this->ky[yindex]*this->cu[3*cindex+2][1] - this->kz[zindex]*this->cu[3*cindex+1][1]); \
             this->cvorticity[3*cindex+1][0] = -(this->kz[zindex]*this->cu[3*cindex+0][1] - this->kx[xindex]*this->cu[3*cindex+2][1]); \
@@ -182,20 +182,20 @@ void fluid_solver<R>::compute_vorticity() \
             this->cvorticity[3*cindex+1][1] =  (this->kz[zindex]*this->cu[3*cindex+0][0] - this->kx[xindex]*this->cu[3*cindex+2][0]); \
             this->cvorticity[3*cindex+2][1] =  (this->kx[xindex]*this->cu[3*cindex+1][0] - this->ky[yindex]*this->cu[3*cindex+0][0]); \
             ); \
-    /*this->symmetrize(this->cvorticity, 3);*/ \
+    this->symmetrize(this->cvorticity, 3); \
 } \
  \
 template<> \
 void fluid_solver<R>::compute_velocity(C *vorticity) \
 { \
     double k2; \
-    /*this->low_pass_Fourier(vorticity, 3, this->kM);*/ \
+    this->low_pass_Fourier(vorticity, 3, this->kM); \
     std::fill_n((R*)this->cu, this->cd->local_size*2, 0.0); \
     CLOOP( \
             k2 = (this->kx[xindex]*this->kx[xindex] + \
                   this->ky[yindex]*this->ky[yindex] + \
                   this->kz[zindex]*this->kz[zindex]); \
-            /*if (k2 < this->kM2)*/ \
+            if (k2 < this->kM2) \
             { \
                 this->cu[3*cindex+0][0] = -(this->ky[yindex]*vorticity[3*cindex+2][1] - this->kz[zindex]*vorticity[3*cindex+1][1]) / k2; \
                 this->cu[3*cindex+1][0] = -(this->kz[zindex]*vorticity[3*cindex+0][1] - this->kx[xindex]*vorticity[3*cindex+2][1]) / k2; \
@@ -204,7 +204,7 @@ void fluid_solver<R>::compute_velocity(C *vorticity) \
                 this->cu[3*cindex+1][1] =  (this->kz[zindex]*vorticity[3*cindex+0][0] - this->kx[xindex]*vorticity[3*cindex+2][0]) / k2; \
                 this->cu[3*cindex+2][1] =  (this->kx[xindex]*vorticity[3*cindex+1][0] - this->ky[yindex]*vorticity[3*cindex+0][0]) / k2; \
             } \
-            /*else \
+            else \
             { \
                 this->cu[3*cindex+0][0] = 0.0; \
                 this->cu[3*cindex+1][0] = 0.0; \
@@ -212,13 +212,13 @@ void fluid_solver<R>::compute_velocity(C *vorticity) \
                 this->cu[3*cindex+0][1] = 0.0; \
                 this->cu[3*cindex+1][1] = 0.0; \
                 this->cu[3*cindex+2][1] = 0.0; \
-            }*/ \
+            } \
             ); \
     if (this->cd->myrank == this->cd->rank[0]) \
         std::fill_n((R*)(this->cu), 6, 0.0); \
-    /*this->impose_zero_modes();*/ \
-    /*this->symmetrize(this->cu, 3);*/ \
-    /*this->low_pass_Fourier(this->cu, 3, this->kM);*/ \
+    this->impose_zero_modes(); \
+    this->symmetrize(this->cu, 3); \
+    this->low_pass_Fourier(this->cu, 3, this->kM); \
 } \
  \
 template<> \
@@ -329,7 +329,7 @@ void fluid_solver<R>::omega_nonlin( \
             ); \
     this->symmetrize(this->cu, 3); \
     this->add_forcing(this->cu, 1.0); \
-    /*this->force_divfree(this->cu);*/ \
+    this->force_divfree(this->cu); \
 } \
  \
 template<> \
@@ -409,6 +409,8 @@ int fluid_solver<R>::read(char field, char representation) \
                 FFTW(execute)(*((FFTW(plan)*)this->r2c_vorticity )); \
         } \
         this->low_pass_Fourier(this->cvorticity, 3, this->kM); \
+        this->force_divfree(this->cvorticity); \
+        this->symmetrize(this->cvorticity, 3); \
         return EXIT_SUCCESS; \
     } \
     if ((field == 'u') && (representation == 'c')) \