still broken.

still debugging. at the moment, fields are set to zero a bunch of times,
but it doesn't seem to help

src/fluid_solver.cpp

@@ -262,11 +262,15 @@ void fluid_solver<R>::omega_nonlin( \
 { \
     assert(src >= 0 && src < 3); \
     /* compute velocity */ \
-    /*std::fill_n((R*)this->cu, 2*this->cd->local_size, 0);*/ \
+    std::fill_n((R*)this->cu, 2*this->cd->local_size, 0); \
+    this->low_pass_Fourier(this->cv[src], 3, this->kM); \
     this->compute_velocity(this->cv[src]); \
+    this->low_pass_Fourier(this->cu, 3, this->kM); \
     /* get fields from Fourier space to real space */ \
-    FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \
-    FFTW(execute)(*((FFTW(plan)*)this->vc2r[src]));     \
+    /*std::fill_n(this->ru, 2*this->cd->local_size, 0);     */ \
+    /*std::fill_n(this->rv[src], 2*this->cd->local_size, 0);*/ \
+    FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity ));  \
+    FFTW(execute)(*((FFTW(plan)*)this->vc2r[src]));      \
     /* compute cross product $u \times \omega$, and normalize */ \
     R tmpx0, tmpy0, tmpz0; \
     RLOOP ( \
@@ -279,9 +283,10 @@ void fluid_solver<R>::omega_nonlin( \
             ); \
     /* go back to Fourier space */ \
     /* TODO: is 0 padding needed here? */ \
+    std::fill_n((R*)this->cu, 2*this->cd->local_size, 0); \
     FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
     this->low_pass_Fourier(this->cu, 3, this->kM); \
-    this->symmetrize(this->cu, 3); \
+    /*this->symmetrize(this->cu, 3);*/ \
     /* $\imath k \times Fourier(u \times \omega)$ */ \
     R tmpx1, tmpy1, tmpz1; \
     CLOOP( \
@@ -298,7 +303,7 @@ void fluid_solver<R>::omega_nonlin( \
             this->cu[3*cindex+1][1] = tmpy1 / this->normalization_factor;\
             this->cu[3*cindex+2][1] = tmpz1 / this->normalization_factor;\
             ); \
-    this->symmetrize(this->cu, 3); \
+    /*this->symmetrize(this->cu, 3);*/ \
     this->add_forcing(this->cu, 1.0); \
 } \
  \

src/fluid_solver_base.cpp

@@ -245,6 +245,7 @@ void fluid_solver_base<R>::symmetrize(C *data, const int howmany) \
     C *buffer; \
     buffer = FFTW(alloc_complex)(howmany*this->cd->sizes[1]); \
     ptrdiff_t yy; \
+    ptrdiff_t tindex; \
     int ranksrc, rankdst; \
     for (yy = 1; yy < this->cd->sizes[0]/2; yy++) { \
         ranksrc = this->cd->rank[yy]; \
@@ -271,13 +272,15 @@ void fluid_solver_base<R>::symmetrize(C *data, const int howmany) \
         if (this->cd->myrank == rankdst) \
         { \
             for (ii = 1; ii < this->cd->sizes[1]; ii++) \
-                for (cc = 0; cc < howmany; cc++) { \
+                for (cc = 0; cc < howmany; cc++) \
+                { \
                     (*((data + howmany*((this->cd->sizes[0] - yy - this->cd->starts[0])*this->cd->sizes[1] + ii)*this->cd->sizes[2]) + cc))[0] = \
                         (*(buffer + howmany*(this->cd->sizes[1]-ii)+cc))[0]; \
                     (*((data + howmany*((this->cd->sizes[0] - yy - this->cd->starts[0])*this->cd->sizes[1] + ii)*this->cd->sizes[2]) + cc))[1] = \
                        -(*(buffer + howmany*(this->cd->sizes[1]-ii)+cc))[1]; \
                 } \
-            for (cc = 0; cc < howmany; cc++) { \
+            for (cc = 0; cc < howmany; cc++) \
+            { \
                 (*((data + cc + howmany*(this->cd->sizes[0] - yy - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2])))[0] =  (*(buffer + cc))[0]; \
                 (*((data + cc + howmany*(this->cd->sizes[0] - yy - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2])))[1] = -(*(buffer + cc))[1]; \
             } \
@@ -288,12 +291,15 @@ void fluid_solver_base<R>::symmetrize(C *data, const int howmany) \
     /* put asymmetric data to 0 */\
     /*if (this->cd->myrank == this->cd->rank[this->cd->sizes[0]/2]) \
     { \
-        for (cc = 0; cc < howmany; cc++) \
+        tindex = howmany*(this->cd->sizes[0]/2 - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2]; \
+        for (ii = 0; ii < this->cd->sizes[1]; ii++) \
         { \
-            data[cc + howmany*(this->cd->sizes[0]/2 - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2]][0] = 0.0; \
-            data[cc + howmany*(this->cd->sizes[0]/2 - this->cd->starts[0])*this->cd->sizes[1]*this->cd->sizes[2]][1] = 0.0; \
+            std::fill_n((R*)(data + tindex), howmany*2*this->cd->sizes[2], 0.0); \
+            tindex += howmany*this->cd->sizes[2]; \
         } \
-    }*/ \
+    } \
+    tindex = howmany*(); \
+    std::fill_n((R*)(data + tindex), howmany*2, 0.0);*/ \
 } \
  \
 template<> \

test.py

@@ -301,5 +301,6 @@ if __name__ == '__main__':
     a.plot(ycoord, amp*np.sin(ycoord), dashes = (2, 2))
     a.plot(ycoord, vfin[0, :, 0, 2])
     a.plot(ycoord, -np.cos(ycoord), dashes = (2, 2))
+    a.grid()
     fig.savefig('ux_vs_y.pdf', format = 'pdf')