diff --git a/cpp/full_code/NSVEparticles.cpp b/cpp/full_code/NSVEparticles.cpp
index e9bdfaf760565ba477f8634d34b98b5719143ff2..ba541554ad357d7344b54bc8b57a332bfdf7b2ed 100644
--- a/cpp/full_code/NSVEparticles.cpp
+++ b/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)
diff --git a/examples/convergence.py b/examples/convergence.py
index 9ee1e6c1da18a04d987eabd64d33c8095f305126..9560ab4d7e61f2f0eef7b94635c558f2994cbf9c 100644
--- a/examples/convergence.py
+++ b/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)