diff --git a/bfps/NavierStokes.py b/bfps/NavierStokes.py
index a075074fb4233e842a9055f74f074dec532dc694..6c0c8b3a5b7c37a21482755344f8548323d7a442 100644
--- a/bfps/NavierStokes.py
+++ b/bfps/NavierStokes.py
@@ -51,6 +51,19 @@ class NavierStokes(_fluid_particle_base):
                 simname = simname,
                 dtype = fluid_precision,
                 use_fftw_wisdom = use_fftw_wisdom)
+        self.parameters['nu'] = 0.1
+        self.parameters['fmode'] = 1
+        self.parameters['famplitude'] = 0.5
+        self.parameters['fk0'] = 1.5
+        self.parameters['fk1'] = 3.0
+        self.parameters['forcing_type'] = 'linear'
+        self.parameters['histogram_bins'] = 256
+        self.parameters['max_velocity_estimate'] = 1.0
+        self.parameters['max_vorticity_estimate'] = 1.0
+        self.parameters['QR2D_histogram_bins'] = 64
+        self.parameters['max_trS2_estimate'] = 1.0
+        self.parameters['max_Q_estimate'] = 1.0
+        self.parameters['max_R_estimate'] = 1.0
         self.file_datasets_grow = """
                 //begincpp
                 std::string temp_string;
@@ -718,6 +731,83 @@ class NavierStokes(_fluid_particle_base):
         with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r+') as ofile:
             kspace = self.get_kspace()
             nshells = kspace['nshell'].shape[0]
+            for k in ['velocity', 'vorticity']:
+                time_chunk = 2**20//(8*3*self.parameters['nx'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/xlines/' + k,
+                                     (1, self.parameters['nx'], 3),
+                                     chunks = (time_chunk, self.parameters['nx'], 3),
+                                     maxshape = (None, self.parameters['nx'], 3),
+                                     dtype = self.dtype,
+                                     compression = 'gzip')
+                time_chunk = 2**20//(8*3*3*nshells)
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/spectra/' + k + '_' + k,
+                                     (1, nshells, 3, 3),
+                                     chunks = (time_chunk, nshells, 3, 3),
+                                     maxshape = (None, nshells, 3, 3),
+                                     dtype = np.float64,
+                                     compression = 'gzip')
+                time_chunk = 2**20//(8*4*10)
+                time_chunk = max(time_chunk, 1)
+                a = ofile.create_dataset('statistics/moments/' + k,
+                                     (1, 10, 4),
+                                     chunks = (time_chunk, 10, 4),
+                                     maxshape = (None, 10, 4),
+                                     dtype = np.float64,
+                                     compression = 'gzip')
+                time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('statistics/histograms/' + k,
+                                     (1,
+                                      self.parameters['histogram_bins'],
+                                      4),
+                                     chunks = (time_chunk,
+                                               self.parameters['histogram_bins'],
+                                               4),
+                                     maxshape = (None,
+                                                 self.parameters['histogram_bins'],
+                                                 4),
+                                     dtype = np.int64,
+                                     compression = 'gzip')
+            for s in range(self.particle_species):
+                time_chunk = 2**20 // (8*3*
+                                       self.parameters['nparticles']*
+                                       self.parameters['tracers{0}_integration_steps'.format(s)])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
+                                     (1,
+                                      self.parameters['tracers{0}_integration_steps'.format(s)],
+                                      self.parameters['nparticles'],
+                                      3),
+                                     maxshape = (None,
+                                                 self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                 self.parameters['nparticles'],
+                                                 3),
+                                     chunks =  (time_chunk,
+                                                self.parameters['tracers{0}_integration_steps'.format(s)],
+                                                self.parameters['nparticles'],
+                                                3),
+                                     dtype = np.float64)
+                time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
+                time_chunk = max(time_chunk, 1)
+                ofile.create_dataset(
+                    '/particles/tracers{0}/velocity'.format(s),
+                    (1,
+                     self.parameters['nparticles'],
+                     3),
+                    chunks = (time_chunk, self.parameters['nparticles'], 3),
+                    maxshape = (None, self.parameters['nparticles'], 3),
+                    dtype = np.float64)
+                if self.parameters['tracers{0}_acc_on'.format(s)]:
+                    ofile.create_dataset(
+                        '/particles/tracers{0}/acceleration'.format(s),
+                        (1,
+                         self.parameters['nparticles'],
+                         3),
+                        chunks = (time_chunk, self.parameters['nparticles'], 3),
+                        maxshape = (None, self.parameters['nparticles'], 3),
+                        dtype = np.float64)
             if self.QR_stats_on:
                 time_chunk = 2**20//(8*3*self.parameters['histogram_bins'])
                 time_chunk = max(time_chunk, 1)
diff --git a/bfps/_fluid_base.py b/bfps/_fluid_base.py
index 31bd4d318ed239191bc63d1af9c524c36a954aa7..4e1a6c6c7824e44c656124fda1e7374666bd21ad 100644
--- a/bfps/_fluid_base.py
+++ b/bfps/_fluid_base.py
@@ -70,19 +70,6 @@ class _fluid_particle_base(_code):
         self.parameters['niter_out'] = 1024
         self.parameters['nparticles'] = 0
         self.parameters['dt'] = 0.01
-        self.parameters['nu'] = 0.1
-        self.parameters['fmode'] = 1
-        self.parameters['famplitude'] = 0.5
-        self.parameters['fk0'] = 1.5
-        self.parameters['fk1'] = 3.0
-        self.parameters['forcing_type'] = 'linear'
-        self.parameters['histogram_bins'] = 256
-        self.parameters['max_velocity_estimate'] = 1.0
-        self.parameters['max_vorticity_estimate'] = 1.0
-        self.parameters['QR2D_histogram_bins'] = 64
-        self.parameters['max_trS2_estimate'] = 1.0
-        self.parameters['max_Q_estimate'] = 1.0
-        self.parameters['max_R_estimate'] = 1.0
         self.fluid_includes = '#include "fluid_solver.hpp"\n'
         self.fluid_variables = ''
         self.fluid_definitions = ''
@@ -403,83 +390,6 @@ class _fluid_particle_base(_code):
             for k in kspace.keys():
                 ofile['kspace/' + k] = kspace[k]
             nshells = kspace['nshell'].shape[0]
-            for k in ['velocity', 'vorticity']:
-                time_chunk = 2**20//(8*3*self.parameters['nx'])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('statistics/xlines/' + k,
-                                     (1, self.parameters['nx'], 3),
-                                     chunks = (time_chunk, self.parameters['nx'], 3),
-                                     maxshape = (None, self.parameters['nx'], 3),
-                                     dtype = self.dtype,
-                                     compression = 'gzip')
-                time_chunk = 2**20//(8*3*3*nshells)
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('statistics/spectra/' + k + '_' + k,
-                                     (1, nshells, 3, 3),
-                                     chunks = (time_chunk, nshells, 3, 3),
-                                     maxshape = (None, nshells, 3, 3),
-                                     dtype = np.float64,
-                                     compression = 'gzip')
-                time_chunk = 2**20//(8*4*10)
-                time_chunk = max(time_chunk, 1)
-                a = ofile.create_dataset('statistics/moments/' + k,
-                                     (1, 10, 4),
-                                     chunks = (time_chunk, 10, 4),
-                                     maxshape = (None, 10, 4),
-                                     dtype = np.float64,
-                                     compression = 'gzip')
-                time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('statistics/histograms/' + k,
-                                     (1,
-                                      self.parameters['histogram_bins'],
-                                      4),
-                                     chunks = (time_chunk,
-                                               self.parameters['histogram_bins'],
-                                               4),
-                                     maxshape = (None,
-                                                 self.parameters['histogram_bins'],
-                                                 4),
-                                     dtype = np.int64,
-                                     compression = 'gzip')
-            for s in range(self.particle_species):
-                time_chunk = 2**20 // (8*3*
-                                       self.parameters['nparticles']*
-                                       self.parameters['tracers{0}_integration_steps'.format(s)])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset('particles/tracers{0}/rhs'.format(s),
-                                     (1,
-                                      self.parameters['tracers{0}_integration_steps'.format(s)],
-                                      self.parameters['nparticles'],
-                                      3),
-                                     maxshape = (None,
-                                                 self.parameters['tracers{0}_integration_steps'.format(s)],
-                                                 self.parameters['nparticles'],
-                                                 3),
-                                     chunks =  (time_chunk,
-                                                self.parameters['tracers{0}_integration_steps'.format(s)],
-                                                self.parameters['nparticles'],
-                                                3),
-                                     dtype = np.float64)
-                time_chunk = 2**20 // (8*3*self.parameters['nparticles'])
-                time_chunk = max(time_chunk, 1)
-                ofile.create_dataset(
-                    '/particles/tracers{0}/velocity'.format(s),
-                    (1,
-                     self.parameters['nparticles'],
-                     3),
-                    chunks = (time_chunk, self.parameters['nparticles'], 3),
-                    maxshape = (None, self.parameters['nparticles'], 3),
-                    dtype = np.float64)
-                if self.parameters['tracers{0}_acc_on'.format(s)]:
-                    ofile.create_dataset(
-                        '/particles/tracers{0}/acceleration'.format(s),
-                        (1,
-                         self.parameters['nparticles'],
-                         3),
-                        chunks = (time_chunk, self.parameters['nparticles'], 3),
-                        maxshape = (None, self.parameters['nparticles'], 3),
-                        dtype = np.float64)
             ofile.close()
         return None