fix nonlinear term computation

I was normalizing many many times because of the previous error, so now
everything seems to make sense, and I get collapse of energy and
enstrophy with those given by vortex

src/fluid_solver.cpp

@@ -271,15 +271,14 @@ void fluid_solver<R>::omega_nonlin( \
     /* compute cross product $u \times \omega$, and normalize */ \
     R tmpx0, tmpy0, tmpz0; \
     RLOOP ( \
-             tmpx0 = (this->ru[(3*rindex)+1]*this->rv[src][(3*rindex)+2] - this->ru[(3*rindex)+2]*this->rv[src][(3*rindex)+1]) / this->normalization_factor; \
-             tmpy0 = (this->ru[(3*rindex)+2]*this->rv[src][(3*rindex)+0] - this->ru[(3*rindex)+0]*this->rv[src][(3*rindex)+2]) / this->normalization_factor; \
-             tmpz0 = (this->ru[(3*rindex)+0]*this->rv[src][(3*rindex)+1] - this->ru[(3*rindex)+1]*this->rv[src][(3*rindex)+0]) / this->normalization_factor; \
+             tmpx0 = (this->ru[(3*rindex)+1]*this->rv[src][(3*rindex)+2] - this->ru[(3*rindex)+2]*this->rv[src][(3*rindex)+1]); \
+             tmpy0 = (this->ru[(3*rindex)+2]*this->rv[src][(3*rindex)+0] - this->ru[(3*rindex)+0]*this->rv[src][(3*rindex)+2]); \
+             tmpz0 = (this->ru[(3*rindex)+0]*this->rv[src][(3*rindex)+1] - this->ru[(3*rindex)+1]*this->rv[src][(3*rindex)+0]); \
              this->ru[(3*rindex)+0] = tmpx0 / this->normalization_factor; \
              this->ru[(3*rindex)+1] = tmpy0 / this->normalization_factor; \
              this->ru[(3*rindex)+2] = tmpz0 / this->normalization_factor; \
             ); \
     /* go back to Fourier space */ \
-    /* TODO: is 0 padding needed here? */ \
     FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
     this->low_pass_Fourier(this->cu, 3, this->kM); \
     this->symmetrize(this->cu, 3); \
@@ -292,12 +291,12 @@ void fluid_solver<R>::omega_nonlin( \
             tmpx1 =  (this->ky[yindex]*this->cu[3*cindex+2][0] - this->kz[zindex]*this->cu[3*cindex+1][0]); \
             tmpy1 =  (this->kz[zindex]*this->cu[3*cindex+0][0] - this->kx[xindex]*this->cu[3*cindex+2][0]); \
             tmpz1 =  (this->kx[xindex]*this->cu[3*cindex+1][0] - this->ky[yindex]*this->cu[3*cindex+0][0]); \
-            this->cu[3*cindex+0][0] = tmpx0 / this->normalization_factor;\
-            this->cu[3*cindex+1][0] = tmpy0 / this->normalization_factor;\
-            this->cu[3*cindex+2][0] = tmpz0 / this->normalization_factor;\
-            this->cu[3*cindex+0][1] = tmpx1 / this->normalization_factor;\
-            this->cu[3*cindex+1][1] = tmpy1 / this->normalization_factor;\
-            this->cu[3*cindex+2][1] = tmpz1 / this->normalization_factor;\
+            this->cu[3*cindex+0][0] = tmpx0; \
+            this->cu[3*cindex+1][0] = tmpy0; \
+            this->cu[3*cindex+2][0] = tmpz0; \
+            this->cu[3*cindex+0][1] = tmpx1; \
+            this->cu[3*cindex+1][1] = tmpy1; \
+            this->cu[3*cindex+2][1] = tmpz1; \
             ); \
     /*this->symmetrize(this->cu, 3);*/ \
     this->add_forcing(this->cu, 1.0); \

test.py

@@ -223,15 +223,17 @@ def convergence_test(opt):
     stats1 = np.fromfile('test1_stats.bin', dtype = dtype)
     stats2 = np.fromfile('test2_stats.bin', dtype = dtype)
     stats_vortex = np.loadtxt('../vortex/sim_000000.log')
+    print stats1
+    print stats2
     fig = plt.figure(figsize = (12,6))
     a = fig.add_subplot(121)
     a.plot(stats1['t'], stats1['energy'])
-    a.plot(stats2['t'], stats2['energy'])
-    a.plot(stats_vortex[:, 2], stats_vortex[:, 3])
+    a.plot(stats2['t'], stats2['energy'], dashes = (3, 3))
+    a.plot(stats_vortex[:, 2], stats_vortex[:, 3], dashes = (2, 4))
     a = fig.add_subplot(122)
     a.plot(stats1['t'], stats1['enstrophy'])
-    a.plot(stats2['t'], stats2['enstrophy'])
-    a.plot(stats_vortex[:, 2], stats_vortex[:, 9]/2)
+    a.plot(stats2['t'], stats2['enstrophy'], dashes = (3, 3))
+    a.plot(stats_vortex[:, 2], stats_vortex[:, 9]/2, dashes = (2, 4))
     fig.savefig('test.pdf', format = 'pdf')
 
     #fig = plt.figure(figsize=(12, 12))