fixes intialization of runs with different interpolation

cpp/full_code/NSVEparticles.cpp

@@ -98,6 +98,7 @@ int NSVEparticles<rnumber>::initialize(void)
         this->particles_output_writer_mpi->close_file();
     }
 
+    DEBUG_MSG("fs->iteration = %d\n", this->fs->iteration);
     DEBUG_MSG_WAIT(MPI_COMM_WORLD, "about to call particles_system_builder\n");
     this->ps = particles_system_builder(
                 this->fs->cvelocity,              // (field object)

examples/convergence.py

@@ -8,13 +8,16 @@ from TurTLE import DNS, PP
 base_niterations = 256
 test_niterations = 32
 
-nparticles = 1000
+nparticles = 100
 particle_random_seed = 15
-base_particle_dt = 0.6
+base_particle_dt = 0.5
+test_niterations_particles = 128
 
 neighbours_smoothness_list = [
         (1, 1),
-        (1, 0)]
+        (2, 1),
+        (3, 2),
+        (4, 3)]
 
 def main():
     #generate_initial_conditions()
@@ -33,7 +36,7 @@ def generate_initial_conditions():
     # number of MPI processes to use
     nprocesses = 4
     # number of OpenMP threads per MPI process to use
-    nthreads_per_process = 2
+    nthreads_per_process = 1
 
 # 1. Generate quasistationary state to use for initial conditions.
     # create a dns object
@@ -104,7 +107,7 @@ def run_simulations_field():
     # number of MPI processes to use
     nprocesses = 4
     # number of OpenMP threads per MPI process to use
-    nthreads_per_process = 2
+    nthreads_per_process = 1
 
 # 1. Run NSVE for the three resolutions.
     for factor in [1, 2, 4]:
@@ -249,22 +252,24 @@ def run_simulations_particles():
     # number of MPI processes to use
     nprocesses = 4
     # number of OpenMP threads per MPI process to use
-    nthreads_per_process = 2
+    nthreads_per_process = 3
 
 # 1. Run NSVEparticles for a few iterations, to build up rhs values, consistent
 #    with the relevant velocity field.
     # create dns object
-    cc = DNS()
     # launch simulation
-    factor = 2
-    cc.launch([
+    factor = 1
+    for neighbours, smoothness in neighbours_smoothness_list:
+        interp_name = 'I{0}O{1}'.format(2*neighbours+2, 2*smoothness+1)
+        cc = DNS()
+        cc.launch([
             'NSVEparticles',
             '-n', '{0}'.format(factor*32),
             '--np', '{0}'.format(nprocesses),
             '--ntpp', '{0}'.format(nthreads_per_process),
             '--src-simname', 'base_{0}x'.format(factor),
             '--src-iteration', '{0}'.format(base_niterations),
-            '--simname', 'nsvep_base',
+            '--simname', 'nsvep_base_' + interp_name,
             '--precision', 'double',
             '--dtfactor', '{0}'.format(base_particle_dt/4),
             '--kMeta', '{0}'.format(1.0*factor),
@@ -273,6 +278,8 @@ def run_simulations_particles():
             '--niter_stat', '{0}'.format(1),
             '--nparticles', '{0}'.format(nparticles),
             '--niter_part', '{0}'.format(1),
+            '--tracers0_neighbours', '{0}'.format(neighbours),
+            '--tracers0_smoothness', '{0}'.format(smoothness),
             '--cpp_random_particles', '{0}'.format(particle_random_seed)])
 
 # 2. Prepare initial conditions
@@ -282,15 +289,15 @@ def run_simulations_particles():
             df = h5py.File('nsvep_{0}x_{1}_checkpoint_0.h5'.format(factor, interp_name), 'w')
             # field
             df['vorticity/complex/{0}'.format(8*factor)] = h5py.ExternalLink(
-                'nsvep_base_checkpoint_0.h5',
+                'nsvep_base_{0}_checkpoint_0.h5'.format(interp_name),
                 'vorticity/complex/16')
             # particles
             df['tracers0/state/{0}'.format(8*factor)] = h5py.ExternalLink(
-                    'nsvep_base_checkpoint_0.h5',
+                    'nsvep_base_{0}_checkpoint_0.h5'.format(interp_name),
                     'tracers0/state/16')
             # copy rhs
             source_file = h5py.File(
-                    'nsvep_base_checkpoint_0.h5', 'r')
+                    'nsvep_base_{0}_checkpoint_0.h5'.format(interp_name), 'r')
             rhs = source_file['tracers0/rhs/16'][()]
             if factor == 4:
                 rhs[0] = source_file['tracers0/rhs/16'][0]
@@ -324,7 +331,7 @@ def run_simulations_particles():
             # launch simulation
             cc.launch([
                     'NSVEparticles',
-                    '-n', '{0}'.format(64),
+                    '-n', '{0}'.format(32),
                     '--np', '{0}'.format(nprocesses),
                     '--ntpp', '{0}'.format(nthreads_per_process),
                     '--simname', 'nsvep_{0}x_I{1}O{2}'.format(
@@ -333,9 +340,9 @@ def run_simulations_particles():
                         2*smoothness+1),
                     '--precision', 'double',
                     '--dtfactor', '{0}'.format(base_particle_dt/factor),
-                    '--kMeta', '{0}'.format(2.0),
-                    '--niter_todo', '{0}'.format(test_niterations*factor+8*factor),
-                    '--niter_out', '{0}'.format(test_niterations*factor+8*factor),
+                    '--kMeta', '{0}'.format(1.0),
+                    '--niter_todo', '{0}'.format(test_niterations_particles*factor+8*factor),
+                    '--niter_out', '{0}'.format(test_niterations_particles*factor+8*factor),
                     '--niter_stat', '{0}'.format(4*factor),
                     '--nparticles', '{0}'.format(nparticles),
                     '--niter_part', '{0}'.format(2*factor),
@@ -366,9 +373,6 @@ def plot_error_particles():
             final_iter_x2 = c2.get_data_file()['iteration'][()]
             f1 = h5py.File(c1.simname + '_checkpoint_0.h5', 'r')
             f2 = h5py.File(c2.simname + '_checkpoint_0.h5', 'r')
-            print('in error plot')
-            print(final_iter_x1 - 8*factor)
-            print(final_iter_x2 - 8*factor*2)
             x1 = f1['tracers0/state/{0}'.format(final_iter_x1)][()]
             x2 = f2['tracers0/state/{0}'.format(final_iter_x2)][()]
             diff = np.sqrt(np.sum((x2-x1)**2, axis = -1))
@@ -410,18 +414,20 @@ def plot_traj_particles():
         for factor in [1, 2, 4]:
             cc = DNS(simname = 'nsvep_{0}x_'.format(factor) + interp_name)
             cc.compute_statistics(iter0 = 8*factor)
-            tt = 14
+            tt = 3
             xx = read_trajectory(cc, tt, iter0 = 8*factor)
-            #a.plot(xx[:, 0], label = '$f = {0}$'.format(factor), marker = 'x')
-            a.plot(xx[:, 0] - xx[:, 0].mean(),
-                   xx[:, 1] - xx[:, 1].mean(),
-                   label = '$f = {0}$'.format(factor), marker = 'x')
+            a.plot(xx[:, 0],
+                   xx[:, 1],
+                   label = interp_name + '$f = {0}$'.format(factor),
+                   marker = 'x',
+                   dashes = (4/factor, 3, 4/factor),
+                   alpha = 0.2)
+            #a.plot(xx[:, 0], xx[:, 1], dashes = (4/factor, 3, 4/factor), alpha = 0.2)
     a.legend(loc = 'best')
     a.set_xlabel('$x$ [code units]')
     a.set_ylabel('$y$ [code units]')
     f.tight_layout()
     f.savefig('trajectories.pdf')
-    f.savefig('trajectories.png')
     f.savefig('trajectories.svg')
     plt.close(f)
     return None
@@ -430,7 +436,7 @@ def read_trajectory(cc, traj_index, iter0):
     cc.parameters['niter_part'] = cc.get_data_file()['parameters/niter_part'][()]
     df = cc.get_particle_file()
     xx = []
-    for ii in range(iter0, cc.get_data_file()['iteration'][()], cc.parameters['niter_part']):
+    for ii in range(iter0, cc.get_data_file()['iteration'][()]+1, cc.parameters['niter_part']):
         xx.append(df['tracers0/position/{0}'.format(ii)][traj_index])
     df.close()
     return np.array(xx)