diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index 3dbcebd3cf294bf42e8769ae354125d29df117fa..8f19b71862049f2e2db05dc6d8ed5363f413e904 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -62,36 +62,45 @@ class NavierStokes(bfps.code):
         return None
     def write_fluid_stats(self):
         self.fluid_includes += '#include <cmath>\n'
-        self.fluid_variables += ('double stats[3];\n' +
+        self.fluid_includes += '#include "fftw_tools.hpp"\n'
+        self.fluid_variables += ('double stats[4];\n' +
                                  'FILE *stat_file;\n')
         self.fluid_definitions += """
                 //begincpp
                 void do_stats(fluid_solver<float> *fsolver)
                 {
-                    double vel_tmp;
+                    double vel_tmp, val_tmp;
                     fsolver->compute_velocity(fsolver->cvorticity);
                     stats[0] = .5*fsolver->correl_vec(fsolver->cvelocity,  fsolver->cvelocity);
                     stats[1] = .5*fsolver->correl_vec(fsolver->cvorticity, fsolver->cvorticity);
-                    fs->ift_velocity();
-                    stats[2] = sqrt(fs->ru[0]*fs->ru[0] +
-                                    fs->ru[1]*fs->ru[1] +
-                                    fs->ru[2]*fs->ru[2]);
-                    for (ptrdiff_t rindex = 1; rindex < fs->rd->local_size; rindex++)
-                    {
-                        vel_tmp = sqrt(fs->ru[rindex*3+0]*fs->ru[rindex*3+0] +
-                                       fs->ru[rindex*3+1]*fs->ru[rindex*3+1] +
-                                       fs->ru[rindex*3+2]*fs->ru[rindex*3+2]);
+                    fsolver->ift_velocity();
+                    val_tmp = (fsolver->ru[0]*fsolver->ru[0] +
+                               fsolver->ru[1]*fsolver->ru[1] +
+                               fsolver->ru[2]*fsolver->ru[2]);
+                    stats[2] = sqrt(val_tmp);
+                    stats[3] = val_tmp*.5;
+                    //for (ptrdiff_t rindex = 0; rindex < (fsolver->cd->local_size/3)*2; rindex++)
+                    RLOOP_FOR_OBJECT(
+                        fsolver,
+                        val_tmp = (fsolver->ru[rindex*3+0]*fsolver->ru[rindex*3+0] +
+                                   fsolver->ru[rindex*3+1]*fsolver->ru[rindex*3+1] +
+                                   fsolver->ru[rindex*3+2]*fsolver->ru[rindex*3+2]);
+                        stats[3] += val_tmp*.5;
+                        vel_tmp = sqrt(val_tmp);
                         if (vel_tmp > stats[2])
                             stats[2] = vel_tmp;
-                    }
+                        );
+                    stats[3] /= fsolver->normalization_factor;
+                    MPI_Allreduce((void*)(stats + 3), (void*)(&val_tmp), 1, MPI_DOUBLE, MPI_SUM, fsolver->rd->comm);
+                    stats[3] = val_tmp;
                     if (myrank == 0)
                     {
                         fwrite((void*)&fsolver->iteration, sizeof(int), 1, stat_file);
                         fwrite((void*)&t, sizeof(double), 1, stat_file);
-                        fwrite((void*)stats, sizeof(double), 3, stat_file);
+                        fwrite((void*)stats, sizeof(double), 4, stat_file);
                     }
-                    fs->write_spectrum("velocity", fs->cvelocity);
-                    fs->write_spectrum("vorticity", fs->cvorticity);
+                    fsolver->write_spectrum("velocity",  fsolver->cvelocity);
+                    fsolver->write_spectrum("vorticity", fsolver->cvorticity);
                 }
                 //endcpp
                 """
@@ -99,7 +108,8 @@ class NavierStokes(bfps.code):
                                      ('t',         np.float64),
                                      ('energy',    np.float64),
                                      ('enstrophy', np.float64),
-                                     ('vel_max',   np.float64)])
+                                     ('vel_max',   np.float64),
+                                     ('renergy',   np.float64)])
         pickle.dump(
                 self.stats_dtype,
                 open(os.path.join(
@@ -127,12 +137,21 @@ class NavierStokes(bfps.code):
                 strncpy(fs->forcing_type, forcing_type, 128);
                 fs->iteration = iter0;
                 fs->read('v', 'c');
+                fs->ift_vorticity();
+                fs->write('v', 'r');
+                fs->compute_velocity(fs->cvorticity);
+                fs->ift_velocity();
+                fs->write('u', 'r');
                 if (myrank == 0)
                 {
                     sprintf(fname, "%s_stats.bin", simname);
                     stat_file = fopen(fname, "ab");
+                    sprintf(fname, "%s_time_i%.5x", simname, iter0);
+                    FILE *time_file = fopen(fname, "rb");
+                    fread((void*)&t, sizeof(double), 1, time_file);
+                    fclose(time_file);
                 }
-                t = dt*iter0;
+                MPI_Bcast((void*)&t, 1, MPI_DOUBLE, 0, fs->cd->comm);
                 do_stats(fs);
                 //endcpp
                 """
@@ -148,6 +167,10 @@ class NavierStokes(bfps.code):
                 if (myrank == 0)
                 {
                     fclose(stat_file);
+                    sprintf(fname, "%s_time_i%.5x", simname, fs->iteration);
+                    FILE *time_file = fopen(fname, "wb");
+                    fwrite((void*)&t, sizeof(double), 1, time_file);
+                    fclose(time_file);
                 }
                 fs->write('v', 'c');
                 fs->write('u', 'r');
@@ -317,17 +340,26 @@ class NavierStokes(bfps.code):
                 os.path.join(self.work_dir, self.name + '_dtype.pickle'), 'r'))
         return np.fromfile(os.path.join(self.work_dir, simname + '_stats.bin'),
                            dtype = dtype)
+    def generate_initial_condition(self, simname = 'test'):
+        np.array([0.0]).tofile(
+                os.path.join(
+                        self.work_dir, simname + '_time_i00000'))
+        self.generate_vector_field(simname = 'test')
+        for species in range(self.particle_species):
+            self.generate_tracer_state(
+                    simname = simname,
+                    species = species)
+        return None
 
 import subprocess
 import matplotlib.pyplot as plt
 
 def test(opt):
-    if opt.run or opt.clean:
+    if opt.clean:
         subprocess.call(['rm {0}/test_*'.format(opt.work_dir)], shell = True)
         subprocess.call(['rm {0}/*.pickle'.format(opt.work_dir)], shell = True)
         subprocess.call(['rm {0}/*.elf'.format(opt.work_dir)], shell = True)
         subprocess.call(['rm {0}/*version_info.txt'.format(opt.work_dir)], shell = True)
-    if opt.clean:
         subprocess.call(['make', 'clean'])
         return None
     c = NavierStokes(work_dir = opt.work_dir)
@@ -337,7 +369,7 @@ def test(opt):
     c.parameters['nu'] = 5.5*opt.n**(-4./3)
     c.parameters['dt'] = 5e-3
     c.parameters['niter_todo'] = opt.nsteps
-    c.parameters['famplitude'] = 1.
+    c.parameters['famplitude'] = 0.
     c.parameters['nparticles'] = 32
     if opt.particles:
         c.add_particles()
@@ -346,16 +378,22 @@ def test(opt):
     c.write_src()
     c.write_par(simname = 'test')
     if opt.run:
-        c.generate_vector_field(simname = 'test')
-        if opt.particles:
-            c.generate_tracer_state(simname = 'test', species = 0)
-            c.generate_tracer_state(simname = 'test', species = 1)
+        if opt.iteration == 0:
+            c.generate_initial_condition(simname = 'test')
         c.run(ncpu = opt.ncpu,
-              simname = 'test')
+              simname = 'test',
+              iter0 = opt.iteration)
+        subprocess.call([
+            'cp',
+            os.path.join(c.work_dir, 'test_rvelocity_i{0:0>5x}'.format(opt.nsteps)),
+            os.path.join(c.work_dir, 'tmp')])
+        opt.iteration += opt.nsteps
+        c.run(ncpu = opt.ncpu,
+              simname = 'test',
+              iter0 = opt.iteration)
     stats = c.read_stats()
     k, enespec = c.read_spec()
     k, ensspec = c.read_spec(field = 'vorticity')
-    k, k2enespec = c.read_spec(field = 'kvelocity')
 
     # plot energy and enstrophy
     fig = plt.figure(figsize = (12, 6))
@@ -367,6 +405,7 @@ def test(opt):
     a.legend(loc = 'best')
     a = fig.add_subplot(122)
     a.plot(stats['t'], stats['energy'], label = 'energy', color = (0, 0, 1))
+    a.plot(stats['t'], stats['renergy'], label = 'energy', color = (1, 0, 1), dashes = (2, 2))
     a.set_ylabel('energy', color = (0, 0, 1))
     a.set_xlabel('$t$')
     for tt in a.get_yticklabels():
@@ -379,28 +418,18 @@ def test(opt):
     fig.savefig('stats.pdf', format = 'pdf')
 
     # plot spectra
-    spec_skip = 64
-    fig = plt.figure(figsize=(12,4))
-    a = fig.add_subplot(131)
+    spec_skip = 1
+    fig = plt.figure(figsize=(12,6))
+    a = fig.add_subplot(121)
     for i in range(0, enespec.shape[0], spec_skip):
         a.plot(k, enespec[i]['val'], color = (i*1./enespec.shape[0], 0, 1 - i*1./enespec.shape[0]))
     a.set_xscale('log')
     a.set_yscale('log')
     a.set_title('velocity')
-    a = fig.add_subplot(132)
+    a = fig.add_subplot(122)
     for i in range(0, ensspec.shape[0], spec_skip):
         a.plot(k, ensspec[i]['val'],
                color = (i*1./ensspec.shape[0], 0, 1 - i*1./ensspec.shape[0]))
-        a.plot(k, k2enespec[i]['val'],
-               color = (0, 1 - i*1./ensspec.shape[0], i*1./ensspec.shape[0]),
-               dashes = (2, 2))
-    a.set_xscale('log')
-    a.set_yscale('log')
-    a.set_title('vorticity')
-    a = fig.add_subplot(133)
-    for i in range(0, ensspec.shape[0], spec_skip):
-        a.plot(k, k2enespec[i]['val'],
-               color = (i*1./ensspec.shape[0], 0, 1 - i*1./ensspec.shape[0]))
         a.plot(k, k**2*enespec[i]['val'],
                color = (0, 1 - i*1./ensspec.shape[0], i*1./ensspec.shape[0]),
                dashes = (2, 2))
@@ -408,5 +437,19 @@ def test(opt):
     a.set_yscale('log')
     a.set_title('vorticity')
     fig.savefig('spectrum.pdf', format = 'pdf')
+
+    ##check io
+    #fig = plt.figure(figsize = (12, 6))
+    #a = fig.add_subplot(121)
+    #b = np.fromfile('data/test_rvelocity_i00004', dtype = np.float32).reshape(c.parameters['nz'], c.parameters['ny'], c.parameters['nx'], 3)
+    #c = np.fromfile('data/tmp', dtype = np.float32).reshape(c.parameters['nz'], c.parameters['ny'], c.parameters['nx'], 3)
+    #bla = np.where(np.abs(b-c) > 0.01)
+    #print bla[0].shape, 3*32**3
+    #a.imshow((b - c)[0, :, :, 0], interpolation = 'none')
+    #a = fig.add_subplot(222)
+    #a.imshow(b[0, :, :, 0], interpolation = 'none')
+    #a = fig.add_subplot(224)
+    #a.imshow(c[0, :, :, 0], interpolation = 'none')
+    #fig.savefig('io_test.pdf', format = 'pdf')
     return None
 
diff --git a/bfps/code.py b/bfps/code.py
index ddb26aba45c19a791ef441be59a634f56b6d67c0..e18d6611fb1bf02f725d6a0f70a89016f5baf60f 100644
--- a/bfps/code.py
+++ b/bfps/code.py
@@ -62,7 +62,9 @@ class code(base):
                     {
                         strcpy(simname, argv[1]);
                         iter0 = atoi(argv[2]);
-                        std::cerr << "myrank = " << myrank << ", simulation name is " << simname << std::endl;
+                        std::cerr << "myrank = " << myrank <<
+                                     ", simulation name is " << simname <<
+                                     " and iter0 is " << iter0 << std::endl;
                     }
                     read_parameters();
                 //endcpp
diff --git a/makefile b/makefile
index 09f353d41f5413dbcf6503392401bd6dc5c6fc26..b0eada98cf3afa13eec58ae34d7f583f0b23c208 100644
--- a/makefile
+++ b/makefile
@@ -33,6 +33,7 @@ DEFINES = #-DNDEBUG
 CFLAGS  = -Wall \
 		  -O2 \
 		  -g \
+		  -ffast-math \
 		  ${LOCAL_INCLUDE} \
 		  ${FFTW_INCLUDE}
 		  #-pg \
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index 6da9d0471518fb2a69842c9014ddae6ca65d0b6b..9d8618a0c060f1b7a4b43b8016e5a05db6d3b0c2 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -204,44 +204,37 @@ 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 \
-            { \
-                this->cu[3*cindex+0][0] = 0.0; \
-                this->cu[3*cindex+1][0] = 0.0; \
-                this->cu[3*cindex+2][0] = 0.0; \
-                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->symmetrize(this->cu, 3);*/ \
 } \
  \
 template<> \
 void fluid_solver<R>::ift_velocity() \
 { \
+    std::fill_n(this->ru, this->cd->local_size*2, 0.0); \
     FFTW(execute)(*((FFTW(plan)*)this->c2r_velocity )); \
 } \
  \
 template<> \
 void fluid_solver<R>::ift_vorticity() \
 { \
+    std::fill_n(this->rvorticity, this->cd->local_size*2, 0.0); \
     FFTW(execute)(*((FFTW(plan)*)this->c2r_vorticity )); \
 } \
  \
 template<> \
 void fluid_solver<R>::dft_velocity() \
 { \
+    std::fill_n((R*)this->cu, this->cd->local_size*2, 0.0); \
     FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
 } \
  \
 template<> \
 void fluid_solver<R>::dft_vorticity() \
 { \
+    std::fill_n((R*)this->cvorticity, this->cd->local_size*2, 0.0); \
     FFTW(execute)(*((FFTW(plan)*)this->r2c_vorticity )); \
 } \
  \
@@ -293,10 +286,12 @@ void fluid_solver<R>::omega_nonlin( \
     /* compute velocity */ \
     this->compute_velocity(this->cv[src]); \
     /* get fields from Fourier space to real space */ \
-    /*std::fill_n(this->ru, 2*this->cd->local_size, 0);     */ \
-    /*std::fill_n(this->rv[src], 2*this->cd->local_size, 0);*/ \
+    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]));      \
+    this->clean_up_real_space(this->ru, 3); \
+    this->clean_up_real_space(this->rv[src], 3); \
     /* compute cross product $u \times \omega$, and normalize */ \
     R tmpx0, tmpy0, tmpz0; \
     RLOOP ( \
@@ -308,9 +303,11 @@ void fluid_solver<R>::omega_nonlin( \
              this->ru[(3*rindex)+2] = tmpz0 / this->normalization_factor; \
             ); \
     /* go back to Fourier space */ \
+    this->clean_up_real_space(this->ru, 3); \
     FFTW(execute)(*((FFTW(plan)*)this->r2c_velocity )); \
-    this->low_pass_Fourier(this->cu, 3, this->kM); \
     this->symmetrize(this->cu, 3); \
+    this->low_pass_Fourier(this->cu, 3, this->kM); \
+    this->force_divfree(this->cu); \
     /* $\imath k \times Fourier(u \times \omega)$ */ \
     R tmpx1, tmpy1, tmpz1; \
     CLOOP( \
@@ -327,8 +324,8 @@ void fluid_solver<R>::omega_nonlin( \
             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); \
+    this->symmetrize(this->cu, 3); \
     this->force_divfree(this->cu); \
 } \
  \
@@ -336,6 +333,13 @@ template<> \
 void fluid_solver<R>::step(double dt) \
 { \
     double k2, factor0, factor1; \
+    std::fill_n((R*)this->cu, this->cd->local_size*2, 0.0); \
+    std::fill_n((R*)this->ru, this->cd->local_size*2, 0.0); \
+    std::fill_n((R*)this->cv[1], this->cd->local_size*2, 0.0); \
+    std::fill_n((R*)this->cv[2], this->cd->local_size*2, 0.0); \
+    std::fill_n((R*)this->rv[0], this->cd->local_size*2, 0.0); \
+    std::fill_n((R*)this->rv[1], this->cd->local_size*2, 0.0); \
+    std::fill_n((R*)this->rv[2], this->cd->local_size*2, 0.0); \
     this->omega_nonlin(0); \
     CLOOP( \
             k2 = (this->kx[xindex]*this->kx[xindex] + \
@@ -350,9 +354,9 @@ void fluid_solver<R>::step(double dt) \
             this->cv[1][3*cindex+2][1] = (this->cv[0][3*cindex+2][1] + dt*this->cu[3*cindex+2][1])*factor0; \
             ); \
  \
-    this->low_pass_Fourier(this->cv[1], 3, this->kM); \
     this->force_divfree(this->cv[1]); \
     this->omega_nonlin(1); \
+    this->low_pass_Fourier(this->cv[1], 3, this->kM); \
     CLOOP( \
             k2 = (this->kx[xindex]*this->kx[xindex] + \
                   this->ky[yindex]*this->ky[yindex] + \
@@ -367,9 +371,9 @@ void fluid_solver<R>::step(double dt) \
             this->cv[2][3*cindex+2][1] = (3*this->cv[0][3*cindex+2][1]*factor0 + (this->cv[1][3*cindex+2][1] + dt*this->cu[3*cindex+2][1])*factor1)*0.25; \
             ); \
  \
-    this->low_pass_Fourier(this->cv[2], 3, this->kM); \
     this->force_divfree(this->cv[2]); \
     this->omega_nonlin(2); \
+    this->low_pass_Fourier(this->cv[2], 3, this->kM); \
     CLOOP( \
             k2 = (this->kx[xindex]*this->kx[xindex] + \
                   this->ky[yindex]*this->ky[yindex] + \
@@ -382,8 +386,21 @@ void fluid_solver<R>::step(double dt) \
             this->cv[3][3*cindex+1][1] = (this->cv[0][3*cindex+1][1]*factor0 + 2*(this->cv[2][3*cindex+1][1] + dt*this->cu[3*cindex+1][1]))*factor0/3; \
             this->cv[3][3*cindex+2][1] = (this->cv[0][3*cindex+2][1]*factor0 + 2*(this->cv[2][3*cindex+2][1] + dt*this->cu[3*cindex+2][1]))*factor0/3; \
             );  \
-    this->low_pass_Fourier(this->cv[0], 3, this->kM); \
+    this->symmetrize(this->cu, 3); \
     this->force_divfree(this->cv[0]); \
+    this->low_pass_Fourier(this->cv[0], 3, this->kM); \
+    this->write_base("cvorticity", this->cvorticity); \
+    this->read_base("cvorticity", this->cvorticity); \
+        this->low_pass_Fourier(this->cvorticity, 3, this->kM); \
+        this->force_divfree(this->cvorticity); \
+        this->symmetrize(this->cvorticity, 3); \
+    /*this->ift_vorticity(); \
+    this->dft_vorticity(); \
+    CLOOP( \
+            for (int component = 0; component < 3; component++) \
+            for (int imag_part = 0; imag_part < 2; imag_part++) \
+                this->cvorticity[cindex*3+component][imag_part] /= this->normalization_factor; \
+            );*/ \
  \
     this->iteration++; \
 } \
diff --git a/src/fluid_solver_base.cpp b/src/fluid_solver_base.cpp
index 0c0eef566441b6d6daa31f48d799446b350ac1ae..adea136a5e4284fd2c0be25f40b9fdcc08cec9c3 100644
--- a/src/fluid_solver_base.cpp
+++ b/src/fluid_solver_base.cpp
@@ -30,6 +30,13 @@
 #include "fftw_tools.hpp"
 
 
+template <class rnumber>
+void fluid_solver_base<rnumber>::clean_up_real_space(rnumber *a, int howmany)
+{
+//    for (ptrdiff_t rindex = 0; rindex < this->cd->local_size*2; rindex += howmany*(this->rd->subsizes[2]+2))
+//        std::fill_n(a+rindex+this->rd->subsizes[2]*howmany, 2*howmany, 0.0);
+}
+
 template <class rnumber>
 void fluid_solver_base<rnumber>::cospectrum(cnumber *a, cnumber *b, double *spec, const double k2exponent)
 {
@@ -117,6 +124,7 @@ fluid_solver_base<R>::fluid_solver_base( \
     ntmp[0] = ny; \
     ntmp[1] = nz; \
     ntmp[2] = nx/2 + 1; \
+    ntmp[3] = 3; \
     this->cd = new field_descriptor<R>( \
             4, ntmp, MPI_CNUM, this->rd->comm);\
  \
diff --git a/src/fluid_solver_base.hpp b/src/fluid_solver_base.hpp
index bb1ebcad6e594cd79b61420ba150731f89564845..b948ab5d474acf961a471287957c9904bf99821a 100644
--- a/src/fluid_solver_base.hpp
+++ b/src/fluid_solver_base.hpp
@@ -76,6 +76,7 @@ class fluid_solver_base
         void low_pass_Fourier(cnumber *a, int howmany, double kmax);
         void force_divfree(cnumber *a);
         void symmetrize(cnumber *a, int howmany);
+        void clean_up_real_space(rnumber *a, int howmany);
         rnumber correl_vec(cnumber *a, cnumber *b);
         void cospectrum(cnumber *a, cnumber *b, double *spec, const double k2exponent = 0.0);
         void write_spectrum(const char *fname, cnumber *a, const double k2exponent = 0.0);
@@ -124,6 +125,23 @@ class fluid_solver_base
     } \
 }
 
+/* real space loop */
+#define RLOOP_FOR_OBJECT(obj, expression) \
+ \
+{ \
+    ptrdiff_t rindex; \
+    for (int zindex = 0; zindex < obj->rd->subsizes[0]; zindex++) \
+    for (int yindex = 0; yindex < obj->rd->subsizes[1]; yindex++) \
+    { \
+        rindex = (zindex * obj->rd->subsizes[1] + yindex)*(obj->rd->subsizes[2]+2); \
+    for (int xindex = 0; xindex < obj->rd->subsizes[2]; xindex++) \
+        { \
+            expression; \
+            rindex++; \
+        } \
+    } \
+}
+
 /*****************************************************************************/
 
 #endif//FLUID_SOLVER_BASE
diff --git a/test.py b/test.py
index 997fe73886652c270ba2f99e38f1b1cfbb58b254..4192f131d7b223e5c1f781a7c2a9c3577e27db2e 100755
--- a/test.py
+++ b/test.py
@@ -268,6 +268,7 @@ def Kolmogorov_flow_test(opt):
 
 if __name__ == '__main__':
     parser = argparse.ArgumentParser()
+    parser.add_argument('--iteration', type = int, dest = 'iteration', default = 0)
     parser.add_argument('--particles', dest = 'particles', action = 'store_true')
     parser.add_argument('--clean', dest = 'clean', action = 'store_true')
     parser.add_argument('--run', dest = 'run', action = 'store_true')
@@ -276,7 +277,25 @@ if __name__ == '__main__':
     parser.add_argument('-n', type = int, dest = 'n', default = 64)
     parser.add_argument('--wd', type = str, dest = 'work_dir', default = 'data')
     opt = parser.parse_args()
-    test_curl(opt)
-    #NStest(opt)
+    #test_curl(opt)
+    NStest(opt)
     #resize_test(opt)
+    Rdata = np.fromfile(
+            'data/test_rvorticity_i00000',
+            dtype = np.float32).reshape(opt.n,
+                                        opt.n,
+                                        opt.n, 3)
+    tdata = Rdata.transpose(3, 0, 1, 2).copy()
+    tdata.tofile('../vortex/input_split_per_component')
+    stats_vortex = np.loadtxt('../vortex/sim_000000.log')
+    dtype = pickle.load(open('data/NavierStokes_dtype.pickle', 'r'))
+    stats = np.fromfile('data/test_stats.bin', dtype = dtype)
+    fig = plt.figure(figsize = (12, 6))
+    a = fig.add_subplot(121)
+    a.plot(stats['t'], stats['energy'])
+    a.plot(stats_vortex[:, 2], stats_vortex[:, 3], dashes = (2, 4))
+    a = fig.add_subplot(122)
+    a.plot(stats['t'], stats['enstrophy'])
+    a.plot(stats_vortex[:, 2], stats_vortex[:, 9]/2, dashes = (2, 4))
+    fig.savefig('vortex_comparison.pdf', format = 'pdf')