add impose_zero_modes. fixes NANs

src/fluid_solver.cpp

@@ -25,6 +25,25 @@
 #include "fftw_tools.hpp"
 
 
+template <class rnumber>
+void fluid_solver<rnumber>::impose_zero_modes()
+{
+    this->cu[0][0] = 0.0;
+    this->cu[1][0] = 0.0;
+    this->cu[2][0] = 0.0;
+
+    this->cv[0][0][1] = 0.0;
+    this->cv[0][1][1] = 0.0;
+    this->cv[0][2][1] = 0.0;
+
+    this->cv[1][0][1] = 0.0;
+    this->cv[1][1][1] = 0.0;
+    this->cv[1][2][1] = 0.0;
+
+    this->cv[2][0][1] = 0.0;
+    this->cv[2][1][1] = 0.0;
+    this->cv[2][2][1] = 0.0;
+}
 /*****************************************************************************/
 /* macro for specializations to numeric types compatible with FFTW           */
 
@@ -169,6 +188,7 @@ void fluid_solver<R>::omega_nonlin( \
             this->cu[cindex*3+1][1] = (this->kz[zindex]*this->cv[src][cindex*3+0][0] - this->kx[xindex]*this->cv[src][cindex*3+2][0]) / k2; \
             this->cu[cindex*3+2][1] = (this->kx[xindex]*this->cv[src][cindex*3+1][0] - this->ky[yindex]*this->cv[src][cindex*3+0][0]) / k2; \
             ); \
+    this->impose_zero_modes(); \
     /* get fields from Fourier space to real space */ \
     FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \
     FFTW(execute)(*((FFTW(plan)*)this->vc2r[src])); \
@@ -187,12 +207,12 @@ void fluid_solver<R>::omega_nonlin( \
     /* $\imath k \times DFT(u \times \omega)$ */ \
     R tmpx1, tmpy1, tmpz1; \
     CLOOP( \
-            tmpx0 = (this->ky[yindex]*this->cu[cindex*3+2][1] - this->kz[zindex]*this->cu[cindex*3+1][1]); \
-            tmpy0 = (this->kz[zindex]*this->cu[cindex*3+0][1] - this->kx[xindex]*this->cu[cindex*3+2][1]); \
-            tmpz0 = (this->kx[xindex]*this->cu[cindex*3+1][1] - this->ky[yindex]*this->cu[cindex*3+0][1]); \
-            tmpx1 = (this->ky[yindex]*this->cu[cindex*3+2][0] - this->kz[zindex]*this->cu[cindex*3+1][0]); \
-            tmpy1 = (this->kz[zindex]*this->cu[cindex*3+0][0] - this->kx[xindex]*this->cu[cindex*3+2][0]); \
-            tmpz1 = (this->kx[xindex]*this->cu[cindex*3+1][0] - this->ky[yindex]*this->cu[cindex*3+0][0]); \
+            tmpx0 = -(this->ky[yindex]*this->cu[cindex*3+2][1] - this->kz[zindex]*this->cu[cindex*3+1][1]); \
+            tmpy0 = -(this->kz[zindex]*this->cu[cindex*3+0][1] - this->kx[xindex]*this->cu[cindex*3+2][1]); \
+            tmpz0 = -(this->kx[xindex]*this->cu[cindex*3+1][1] - this->ky[yindex]*this->cu[cindex*3+0][1]); \
+            tmpx1 =  (this->ky[yindex]*this->cu[cindex*3+2][0] - this->kz[zindex]*this->cu[cindex*3+1][0]); \
+            tmpy1 =  (this->kz[zindex]*this->cu[cindex*3+0][0] - this->kx[xindex]*this->cu[cindex*3+2][0]); \
+            tmpz1 =  (this->kx[xindex]*this->cu[cindex*3+1][0] - this->ky[yindex]*this->cu[cindex*3+0][0]); \
             this->cu[cindex*3+0][0] = tmpx0;\
             this->cu[cindex*3+1][0] = tmpy0;\
             this->cu[cindex*3+2][0] = tmpz0;\

src/fluid_solver.hpp

@@ -74,6 +74,7 @@ class fluid_solver:public fluid_solver_base<rnumber>
 
         void omega_nonlin(int src);
         void step(double dt);
+        void impose_zero_modes();
 };
 
 #endif//FLUID_SOLVER

src/fluid_solver_base.cpp

@@ -106,6 +106,13 @@ fluid_solver_base<R>::fluid_solver_base( \
         else \
             this->knullz[i] = true; \
     } \
+    this->kM = this->kMx; \
+    if (this->kM < this->kMy) this->kM = this->kMy; \
+    if (this->kM < this->kMz) this->kM = this->kMz; \
+    this->kM2 = this->kM * this->kM; \
+    this->dk = this->dkx; \
+    if (this->dk > this->dky) this->dk = this->dky; \
+    if (this->dk > this->dkz) this->dk = this->dkz; \
 } \
  \
 template<> \
@@ -133,14 +140,18 @@ R fluid_solver_base<R>::correl_vec(C *a, C *b) \
             k2 = (this->kx[xindex]*this->kx[xindex] + \
                   this->ky[yindex]*this->ky[yindex] + \
                   this->kz[zindex]*this->kz[zindex]); \
-            factor = (xindex == 0) ? 1 : 2; \
-            val_process += factor * ((*(a + cindex))[0] * (*(b + cindex))[0] + \
-                                     (*(a + cindex))[1] * (*(b + cindex))[1]); \
+            if (k2 < this->kM2) \
+            { \
+                factor = (xindex == 0) ? 1 : 2; \
+                val_process += factor * ((*(a + cindex))[0] * (*(b + cindex))[0] + \
+                                         (*(a + cindex))[1] * (*(b + cindex))[1]); \
+            } \
             );\
     MPI_Allreduce( \
             (void*)(&val_process), \
             (void*)(&val), \
             1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD); \
+    /*return R(val / (this->rd->full_size / 3)) / (this->rd->full_size / 3);*/ \
     return R(val); \
 }