fix RLOOP bug

this seems to fix the laminar solution to Kolmogorov flow problem.

src/fluid_solver.cpp

@@ -286,7 +286,7 @@ void fluid_solver<R>::omega_nonlin( \
     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( \
@@ -303,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.hpp

@@ -113,7 +113,7 @@ class fluid_solver_base
     for (int zindex = 0; zindex < this->rd->subsizes[0]; zindex++) \
     for (int yindex = 0; yindex < this->rd->subsizes[1]; yindex++) \
     { \
-        rindex = (zindex * this->rd->subsizes[1] + yindex)*this->rd->subsizes[2]; \
+        rindex = (zindex * this->rd->subsizes[1] + yindex)*(this->rd->subsizes[2]+2); \
     for (int xindex = 0; xindex < this->rd->subsizes[2]; xindex++) \
         { \
             expression; \

test.py

@@ -138,6 +138,7 @@ class stat_test(code):
                 fs->write('v', 'c');
                 fs->write('v', 'r');
                 fs->write('u', 'r');
+                fs->write('u', 'c');
                 delete fs;
                 //endcpp
                 """
@@ -301,6 +302,33 @@ 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.set_xticks(np.linspace(.0, 2*np.pi, 9))
+    a.set_xlim(.0, 2*np.pi)
     a.grid()
     fig.savefig('ux_vs_y.pdf', format = 'pdf')
 
+    fig = plt.figure(figsize=(12, 12))
+    tmp = np.fromfile('test_cvelocity_i{0:0>5x}'.format(stats.shape[0]-1),
+                      dtype = np.complex64).reshape(opt.n, opt.n, opt.n/2+1, 3)
+    a = fig.add_subplot(321)
+    a.plot(np.sum(np.abs(tmp), axis = (1, 2)))
+    a.set_yscale('log')
+    a = fig.add_subplot(323)
+    a.plot(np.sum(np.abs(tmp), axis = (0, 2)))
+    a.set_yscale('log')
+    a = fig.add_subplot(325)
+    a.plot(np.sum(np.abs(tmp), axis = (0, 1)))
+    a.set_yscale('log')
+    a = fig.add_subplot(322)
+    tmp = np.fromfile('test_cvorticity_i{0:0>5x}'.format(stats.shape[0]-1),
+                      dtype = np.complex64).reshape(opt.n, opt.n, opt.n/2+1, 3)
+    a.plot(np.sum(np.abs(tmp), axis = (1, 2)))
+    a.set_yscale('log')
+    a = fig.add_subplot(324)
+    a.plot(np.sum(np.abs(tmp), axis = (0, 2)))
+    a.set_yscale('log')
+    a = fig.add_subplot(326)
+    a.plot(np.sum(np.abs(tmp), axis = (0, 1)))
+    a.set_yscale('log')
+    fig.savefig('cvort.pdf', format = 'pdf')
+