Tutorial: convergence test¶
An essential property of numerical methods is their accuracy. With standard approaches to differential equations, one may use generic convergence tests to characterize algorithms and/or implementations.
As an example usage of TurTLE, we show here how to characterize convergence for the field solver, or the AdamsBashforth particle solver. The purpose of this tutorial is to introduce TurTLE users to nontrivial workflows that involve keeping track of a number of simulations and prescribing their initial conditions.
The field solver test is a fairly simple pursuit, showing basic operations with the turtle launcher.
The particle solver test is also minimal within the methodological constraints imposed by the use of the Adams Bashforth integration method, leading to a relatively increased complexity. The example shows how to account for a number of technical details: full control over multi-step method initialization and full control of fluid solver forcing parameters.
Navier-Stokes solution¶
TurTLE implements two distinct solvers of the Navier-Stokes equations: one for their vorticity formulation, and one for the more common velocity formulation. While there are inescapable differences in the two implementations, they both use the same time-stepping algorithm (a 3rd order RungeKutta method).
We test convergence properties for a “statistically stationary initial condition”, since the smoothness properties of the field will influence convergence properties in general. As an example of the variety of questions that we can address with TurTLE, we also perform a direct comparison of the NSVE (vorticity solver) and NSE (velocity solver) results.
A basic Python script that performs this convergence test is provided under examples/convergence/fields_temporal.py, and it is discussed below. In brief, simulations are performed using the same initial conditions, but different time steps. Error estimates are computed as differences between the different solutions. The figure shows the L2 norm of these differences, normalized by the root-mean-squared velocity magnitude. Note that the actual value of this “relative error” is irrelevant (it depends on total integration time, spatial resolution, etc), but the plot shows that (1) solutions converge with the 3rd order of the timestep and (2) the systematic difference between NSE and NSVE behaves quite differently from the two individual error estimates.
Particle tracking errors¶
There is a more complex mechanism that generates numerical errors for particle tracking in this setting. In particular, the accuracy of the field itself is relevant, as well as the accuracy of the interpolation method coupled to the ODE solver. For instance the default particle tracking solver is a 4th order Adams Bashforth method in TurTLE, which seems at first unnecessary since the fluid only provides a 3rd order solution. The default interpolation method is that of cubic splines, which should in principle limit ODE error convergence to 2nd order (because convergence relies on existence of relevant Taylor expansion terms, and interpolated fields only have a finite number of continuous derivatives). This discussion is briefly continued below, but we will not do it justice in this tutorial.
We present below a Python script that performs a particle tracking convergence test, provided as examples/convergence/particles_temporal.py. The figure shows trajectory integration errors, in a fairly simple setting, for two interpolation methods: cubic splines (I4O3) and 5th order splines on a kernel of 8 points (I8O5). Furthermore, we show both the maximum and the mean errors for the same set of trajectories. While both mean errors seem to converge with the 3rd order of the timestep (consistent with the error of the velocity field), the maximum error obtained for cubic splines converges with the second order — which is the expected behavior for particle tracking in a field that only has one continuous derivative.
The results would merit a longer and more careful analysis outside the scope of this particular tutorial, so we limit ourselves to noting that (1) different ODEs will have different numerical properties (i.e. rotating particles, or inertial particles) and (2) interpolation errors may in general grow with the Reynolds number, since interpolation errors depend on values of gradients (however the distribution of extreme gradient values also changes with Reynolds number, so the overall behavior is complex).
Python code for fields¶
We present here a Python script that performs simple convergence tests with TurTLE. The script itself is provided as the file examples/convergence/fields_temporal.py.
The preamble imports standard packages, TurTLE itself, and a few simple tools from examples/convergence/fields_base.py (not discussed here):
import numpy as np
import h5py
import matplotlib.pyplot as plt
import TurTLE
from TurTLE import DNS, PP
from fields_base import *
As explained above, we should test convergence for a statistically stationary regime. The base simulation will use a real-space grid of \(128 \times 128 \times 128\) points, and experience shows that reasonable quasistationarity is reached for this problem-size in a few hundred iterations (here base_niterations). For the convergence test itself, we use a relatively small number of iterations (i.e. test_niterations), otherwise deterministic chaos will dominate the difference between different solutions (rather than numerical integration errors).
base_niterations = 512
test_niterations = 32
grid_size = 128
We distinguish between three steps of this script. First, we generate the statistically quasistationary regime. Second, we run the fluid solver. Finally, we “analyze” the results (i.e. we generate simple error plots).
def main():
generate_initial_conditions()
run_simulations_field(
err_vs_t = False)
plot_error_field(
err_vs_t = False)
return None
The err_vs_t parameter may be used if one is interested in looking at the difference between the NSE and NSVE solutions as a function of time. It turns on a fine-grained output of the fields in run_simulations_field and it turns on the corresponding plot in plot_error_field.
generate_initial_conditions performs a few calls to C++ executable codes generated with TurTLE. In brief, it
generates a vorticity field in a quasistationary Navier-Stokes regime (using the NSVE solver).
generates the corresponding velocity field needed for the NSE solver
We use the launch method of the DNS and PP Python classes — the wrappers around the C++ direct_numerical_simulation and postprocess functionality. This means that the current Python code may be translated to a shell script directly, only by writing out the arguments (rather than placing them in a list of Python strings). The full function is given below:
def generate_initial_conditions():
# change these two values as needed.
# number of MPI processes to use
nprocesses = 8
# number of OpenMP threads per MPI process to use
nthreads_per_process = 1
# 1. Generate quasistationary state to use for initial conditions.
# create a dns object
c0 = DNS()
# launch the simulation
c0.launch([
'NSVE',
'-n', '{0}'.format(grid_size),
'--np', '{0}'.format(nprocesses),
'--ntpp', '{0}'.format(nthreads_per_process),
'--precision', 'double',
'--src-simname', 'B32p1e4',
'--src-wd', TurTLE.data_dir,
'--src-iteration', '0',
'--simname', 'base',
'--niter_todo', '{0}'.format(base_niterations),
'--niter_out', '{0}'.format(base_niterations),
'--overwrite-src-parameters',
'--kMeta', '1.5',
'--fk0', '2',
'--fk1', '3',
'--niter_stat', '16',
'--checkpoints_per_file', '{0}'.format(64)])
# 3. Generate initial conditions for NSE runs.
# create postprocess object
cc = PP()
# launches code to compute velocity field at last iteration
cc.launch([
'get_velocity',
'--np', '{0}'.format(nprocesses),
'--ntpp', '{0}'.format(nthreads_per_process),
'--simname', 'base',
'--precision', 'double',
'--iter0', '{0}'.format(base_niterations),
'--iter1', '{0}'.format(base_niterations)])
return None
run_simulations_field runs six DNS (3 resolutions for each of NSE and NSVE solvers):
def run_simulations_field(err_vs_t = False):
# change these two values as needed.
# number of MPI processes to use
nprocesses = 8
# number of OpenMP threads per MPI process to use
nthreads_per_process = 1
if err_vs_t:
niter_out_factor = 1
else:
niter_out_factor = test_niterations
# 1. Run NSVE for the three resolutions.
for factor in [1, 2, 4]:
# create dns object
cc = DNS()
# launch simulation
cc.launch([
'NSVE',
'-n', '{0}'.format(grid_size),
'--np', '{0}'.format(nprocesses),
'--ntpp', '{0}'.format(nthreads_per_process),
'--src-simname', 'base',
'--src-iteration', '{0}'.format(base_niterations),
'--simname', 'nsve_{0}x'.format(factor),
'--precision', 'double',
'--dtfactor', '{0}'.format(0.5 / factor),
'--niter_todo', '{0}'.format(test_niterations*factor),
'--niter_out', '{0}'.format(niter_out_factor*factor),
'--niter_stat', '{0}'.format(factor)])
# 2. Run NSE for the three resolutions.
for factor in [1, 2, 4]:
# create dns object
cc = DNS()
# launch simulation
cc.launch([
'NSE',
'-n', '{0}'.format(grid_size),
'--np', '{0}'.format(nprocesses),
'--ntpp', '{0}'.format(nthreads_per_process),
'--src-simname', 'base',
'--src-iteration', '{0}'.format(base_niterations),
'--simname', 'nse_{0}x'.format(factor),
'--precision', 'double',
'--dtfactor', '{0}'.format(0.5 / factor),
'--niter_todo', '{0}'.format(test_niterations*factor),
'--niter_out', '{0}'.format(niter_out_factor*factor),
'--niter_out', '{0}'.format(test_niterations*factor),
'--niter_stat', '{0}'.format(factor)])
return None
Once the runs are finished successfully, the code calls plot_error_field. Here we first read the output of the six different runs, and we compute their statistics, since the call to compute_statistics will also perform some sanity checks on the output. We then compute the error analysis for each of the two solvers, and also compare the two solutions. The get_velocity and compute_vector_field_distance methods are defined in examples/convergence/fields_base.py.
def plot_error_field(err_vs_t = False):
factor_list = [1, 2, 4]
c0_list = [DNS(simname = 'nsve_{0}x'.format(factor))
for factor in factor_list]
c1_list = [DNS(simname = 'nse_{0}x'.format(factor))
for factor in factor_list]
cl = [c0_list, c1_list]
# sanity checks
for cc in c0_list + c1_list:
cc.compute_statistics()
# errors for individual solvers
error_list = [[], [], []]
for ii in range(len(factor_list)-1):
factor = factor_list[ii]
for jj in [0, 1]:
vel1 = get_velocity(cl[jj][ii ], iteration = test_niterations*factor)
vel2 = get_velocity(cl[jj][ii+1], iteration = 2*test_niterations*factor)
dd = compute_vector_field_distance(vel1, vel2, figname = cl[jj][ii].simname + '_vs_' + cl[jj][ii+1].simname)
error_list[jj].append(dd['L2_rel'])
# comparisons of two solutions
for ii in range(len(factor_list)):
factor = factor_list[ii]
vel1 = get_velocity(cl[0][ii], iteration = test_niterations*factor)
vel2 = get_velocity(cl[1][ii], iteration = test_niterations*factor)
dd = compute_vector_field_distance(vel1, vel2)
error_list[2].append(dd['L2_rel'])
f = plt.figure(figsize = (4, 4))
a = f.add_subplot(111)
a.plot(factor_list[:len(error_list[0])],
error_list[0],
marker = '.',
dashes = (2, 3),
label = 'NSVE')
a.plot(factor_list[:len(error_list[1])],
error_list[1],
marker = '.',
dashes = (3, 5),
label = 'NSE')
a.plot(factor_list,
error_list[2],
marker = '.',
label = 'NSVE vs NSE')
fl = np.array(factor_list).astype(np.float)
for ee in [2, 3]:
a.plot(fl[:2], error_list[0][0] * fl[:2]**(-ee),
dashes = (ee, ee),
color = 'black',
label = '$\propto f^{{-{0}}}$'.format(ee),
zorder = -ee)
a.set_ylabel('relative error')
a.set_xlabel('resolution factor $f$')
a.set_xscale('log')
a.set_yscale('log')
a.legend(loc = 'best', fontsize = 6)
f.tight_layout()
f.savefig('err_vs_dt_field.pdf')
f.savefig('err_vs_dt_field.svg')
plt.close(f)
if err_vs_t:
t = range(test_niterations+1)
err_vs_t = [[], [], []]
# comparisons of two solvers
for ii in range(len(factor_list)):
factor = factor_list[ii]
for tt in t:
vel1 = get_velocity(cl[0][ii], iteration = tt*factor)
vel2 = get_velocity(cl[1][ii], iteration = tt*factor)
dd = compute_vector_field_distance(vel1, vel2)
err_vs_t[ii].append(dd['L2_rel'])
# distance between NSVE and NSE as a function of time
f = plt.figure(figsize = (4, 4))
a = f.add_subplot(111)
a.plot(t, err_vs_t[0], label = 'f = 1')
a.plot(t, err_vs_t[1], label = 'f = 2')
a.plot(t, err_vs_t[2], label = 'f = 4')
a.set_yscale('log')
a.legend(loc = 'best', fontsize = 6)
f.tight_layout()
f.savefig('err_vs_t_field.pdf')
f.savefig('err_vs_t_field.svg')
plt.close(f)
return None
Finally, the script ends with a standard call to main:
if __name__ == '__main__':
main()
Python code for particles¶
We present here a Python script that performs a particle tracking convergence test with TurTLE, in the current default configuration (a particles_system object coupled with an NSVE fluid solver).
The preamble imports standard packages, and TurTLE itself:
import numpy as np
import h5py
import matplotlib.pyplot as plt
import TurTLE
from TurTLE import DNS
We then set a few basic parameters, including a list of interpolation parameters. In brief, the spline interpolations used in TurTLE are parametrized by the size \(I\) of the interpolation kernel and the order \(O\) of the polynomial. In this example code we emphasize two equivalences:
\(I = 2 n + 2\) where \(n\) is the number of neighbouring points (in one orientation) required to interpolate the field within one cell.
\(O = 2 m + 1\) where \(m\) is the order of highest continuous derivative of the interpolating polynomial.
More details on the interpolation method can be found in [lalescu2010jcp].
base_niterations = 256
nparticles = 1000
particle_random_seed = 15
base_particle_dt = 0.5
test_niterations_particles = 128
neighbours_smoothness_list = [
(1, 1),
(3, 2)]
We distinguish between three main parts of this script: generation of initial conditions with statistically stationary regime, computation of numerical solutions, analysis of results (plot commands):
def main():
generate_initial_conditions()
run_simulations_particles()
plot_error_particles()
plot_traj_particles()
return None
generate_initial_conditions consists of a single launch of an NSVE run:
def generate_initial_conditions():
# change these two values as needed.
# number of MPI processes to use
nprocesses = 8
# number of OpenMP threads per MPI process to use
nthreads_per_process = 1
# 1. Generate quasistationary state to use for initial conditions.
# create a dns object
c0 = DNS()
# launch the simulation
c0.launch([
'NSVE',
'-n', '32',
'--np', '{0}'.format(nprocesses),
'--ntpp', '{0}'.format(nthreads_per_process),
'--precision', 'double',
'--src-simname', 'B32p1e4',
'--src-wd', TurTLE.data_dir,
'--src-iteration', '0',
'--simname', 'base',
'--niter_todo', '{0}'.format(base_niterations),
'--niter_out', '{0}'.format(base_niterations),
'--overwrite-src-parameters',
'--kMeta', '1.',
'--niter_stat', '16',
'--checkpoints_per_file', '{0}'.format(16)])
return None
We then define the function that runs the particle simulations, run_simulations_particles. Because we use an Adams-Bashforth method, we need to know the values of the particle velocity (i.e. the right-hand-side of the ODE) at previous steps. By default TurTLE accomplishes this by using the Euler integration scheme for “time step 0”, then using a 2nd order Adams-Bashforth, etc, until the desired number of right-hand-side values is stored. Thus the code first integrates trajectories for 8 iterations, and the output will contain the additional information. This happens because such information is crucial for checkpointing — i.e. stopping/restarting the simulation at arbitrary iterations without affecting the numerical solution.
To put it simply, for the simulation “nsvep_base” the dataset “tracers0/rhs/0” contains only zeros, and TurTLE ignores them. However, the “tracers0/rhs/8” dataset for simulation “nsvep_1x” contains meaningful values (that are used because we request that iter0 is different from 0). In order to generate consistent values of the “rhs” datasets, we first run a simulation at the highest resolution (smallest time-step), and then we extract the appropriate snapshots — see also the second part of the trajectory plot function.
Finally, the launch method is called in an unusual way for TurTLE. Firstly, we demand a DNS that starts at an iteration different from 0. Secondly, we manually copy the forcing parameters from the source simulation; since the source cannot be specified explicitly, the DNS class cannot search for the source parameters itself. Thirdly, I will note that one may in fact break this call by modifying the iter0 and checkpoints_per_file values, since we have explicitly placed iter0 in the checkpoint 0 file. While a full discussion of the iter0 and checkpoint subtleties is beyond the scope of this tutorial, I will mention that it was a conscious choice to use a simple iteration to checkpoint correspondence, since we don’t actually need anything smarter outside of extreme scenarios that need to be addressed individually anyway due to other complexities.
def run_simulations_particles():
# change these two values as needed.
# number of MPI processes to use
nprocesses = 8
# number of OpenMP threads per MPI process to use
nthreads_per_process = 1
# 1. Run NSVEparticles for a few iterations, to build up rhs values, consistent
# with the relevant velocity field.
factor = 1
for neighbours, smoothness in neighbours_smoothness_list:
interp_name = 'I{0}O{1}'.format(2*neighbours+2, 2*smoothness+1)
# create dns object
cc = DNS()
# launch simulation
cc.launch([
'NSVEparticles',
'-n', '{0}'.format(factor*32),
'--np', '{0}'.format(nprocesses),
'--ntpp', '{0}'.format(nthreads_per_process),
'--src-simname', 'base',
'--src-iteration', '{0}'.format(base_niterations),
'--simname', 'nsvep_base_' + interp_name,
'--precision', 'double',
'--dtfactor', '{0}'.format(base_particle_dt/4),
'--kMeta', '{0}'.format(factor*1.0),
'--niter_todo', '{0}'.format(20),
'--niter_out', '{0}'.format(1),
'--niter_stat', '{0}'.format(1),
'--nparticles', '{0}'.format(nparticles),
'--niter_part', '{0}'.format(1),
'--tracers0_neighbours', '{0}'.format(neighbours),
'--tracers0_smoothness', '{0}'.format(smoothness),
'--tracers0_integration_steps', '4',
'--cpp_random_particles', '{0}'.format(particle_random_seed),
'--checkpoints_per_file', '{0}'.format(16)])
# 2. Prepare initial conditions
for neighbours, smoothness in neighbours_smoothness_list:
interp_name = 'I{0}O{1}'.format(2*neighbours+2, 2*smoothness+1)
for factor in [1, 2, 4]:
df = h5py.File('nsvep_{0}x_{1}_checkpoint_0.h5'.format(factor, interp_name), 'w')
# field
field_iteration = 16
df['vorticity/complex/{0}'.format(8*factor)] = h5py.ExternalLink(
'nsvep_base_{0}_checkpoint_0.h5'.format(interp_name),
'vorticity/complex/{0}'.format(field_iteration))
# particles
df['tracers0/state/{0}'.format(8*factor)] = h5py.ExternalLink(
'nsvep_base_{0}_checkpoint_0.h5'.format(interp_name),
'tracers0/state/16')
# copy rhs
source_file = h5py.File(
'nsvep_base_{0}_checkpoint_0.h5'.format(interp_name), 'r')
rhs = source_file['tracers0/rhs/16'][()]
for tt in range(1, rhs.shape[0]):
ii = 16 - tt*4//factor + 1
#print(factor, tt, ii)
rhs[tt] = source_file['tracers0/rhs/{0}'.format(16 - tt*4//factor + 1)][1]
df['tracers0/rhs/{0}'.format(8*factor)] = rhs
df.close()
# 3. Run NSVEparticles
for factor in [1, 2, 4]:
# Test different interpolations.
# The loop iterates over number of neighbours and smoothness.
# The simulation names are defined based on "kernel size" and "order
# of polynomial".
# We do this to emphasize the one-to-one correspondence between
# neighbours and kernel size, and smoothness and order of polynomial.
for neighbours, smoothness in neighbours_smoothness_list:
interp_name = 'I{0}O{1}'.format(2*neighbours+2, 2*smoothness+1)
source_file = h5py.File('nsvep_base_{0}.h5'.format(interp_name), 'r')
# create dns object
cc = DNS()
# launch simulation
cc.launch([
'NSVEparticles',
'-n', '{0}'.format(source_file['parameters/nx'][()]),
'--np', '{0}'.format(nprocesses),
'--ntpp', '{0}'.format(nthreads_per_process),
'--simname', 'nsvep_{0}x_{1}'.format(factor, interp_name),
'--precision', 'double',
'--dtfactor', '{0}'.format(base_particle_dt/factor),
'--niter_todo', '{0}'.format(test_niterations_particles*factor+8*factor),
'--niter_out', '{0}'.format(test_niterations_particles*factor+8*factor),
'--niter_stat', '{0}'.format(1),
## ensure fluid physics is the same
'--dealias_type', '{0}'.format(source_file['parameters/dealias_type'][()]),
'--energy', '{0}'.format(source_file['parameters/energy'][()]),
'--famplitude', '{0}'.format(source_file['parameters/famplitude'][()]),
'--fk0', '{0}'.format(source_file['parameters/fk0'][()]),
'--fk1', '{0}'.format(source_file['parameters/fk1'][()]),
'--fmode', '{0}'.format(source_file['parameters/fmode'][()]),
'--friction_coefficient', '{0}'.format(source_file['parameters/friction_coefficient'][()]),
'--injection_rate', '{0}'.format(source_file['parameters/injection_rate'][()]),
'--nu', '{0}'.format(source_file['parameters/nu'][()]),
'--forcing_type', '{0}'.format(str(source_file['parameters/forcing_type'][()], 'ascii')),
## particle params
'--nparticles', '{0}'.format(nparticles),
'--niter_part', '{0}'.format(1),
'--tracers0_neighbours', '{0}'.format(neighbours),
'--tracers0_smoothness', '{0}'.format(smoothness),
'--tracers0_integration_steps', '4',
'--cpp_random_particles', '0', # turn off C++ particle initialization
'--checkpoints_per_file', '{0}'.format(16)],
iter0 = 8*factor)
source_file.close()
return None
To compute the trajectory integration error, we proceed as explained above, and compare each DNS with its corresponding double-the-resolution DNS:
def plot_error_particles():
f = plt.figure()
a = f.add_subplot(111)
factor_list = [1, 2]
factor_list = np.array(factor_list).astype(np.float)
for neighbours, smoothness in neighbours_smoothness_list:
interp_name = 'I{0}O{1}'.format(2*neighbours+2, 2*smoothness+1)
err_max_list = []
err_mean_list = []
for factor in factor_list:
factor = int(factor)
c1 = DNS(simname = 'nsvep_{0}x_'.format(int(factor)) + interp_name)
c2 = DNS(simname = 'nsvep_{0}x_'.format(int(factor)*2) + interp_name)
for cc in [c1, c2]:
cc.parameters['niter_part'] = cc.get_data_file()['parameters/niter_part'][()]
c1.compute_statistics(iter0 = 8*factor)
c2.compute_statistics(iter0 = 8*factor*2)
final_iter_x1 = c1.get_data_file()['iteration'][()]
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')
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))
err_max_list.append(np.max(diff))
err_mean_list.append(np.mean(diff))
f1.close()
f2.close()
err_max_list = np.array(err_max_list)
err_mean_list = np.array(err_mean_list)
a.plot(factor_list, err_max_list, marker = '.', label = 'max error ' + interp_name)
a.plot(factor_list, err_mean_list, marker = '.', label = 'mean error ' + interp_name)
for ee in [2, 3, 4]:
a.plot(factor_list,
(1e-4)*factor_list**(-ee),
color = 'black',
dashes = (ee, ee),
label = '$\propto f^{{-{0}}}$'.format(ee))
a.set_xscale('log')
a.set_yscale('log')
a.set_xlabel('resolution factor f')
a.set_ylabel('absolute trajectory error [DNS units]')
a.legend(loc = 'best', fontsize = 7)
f.tight_layout()
f.savefig('err_vs_dt_particle.pdf')
f.savefig('err_vs_dt_particle.svg')
plt.close(f)
return None
For the beginner’s convenience, we also provide a method to plot individual particle trajectories. Note that the second part of this function is in fact a sanity check for the initial conditions: the plot confirms that the data required for the multi-step method (explicitly placed in the checkpoint file) is consistent with the simulation parameters.
def plot_traj_particles():
###########################################################################
f = plt.figure()
a = f.add_subplot(111)
for neighbours, smoothness in neighbours_smoothness_list:
interp_name = 'I{0}O{1}'.format(2*neighbours+2, 2*smoothness+1)
for factor in [1, 2, 4]:
cc = DNS(simname = 'nsvep_{0}x_'.format(factor) + interp_name)
cc.compute_statistics(iter0 = 8*factor)
tt = 3
xx = read_trajectory(cc, tt, iter0 = 8*factor)
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', fontsize = 7)
a.set_xlabel('$x$ [code units]')
a.set_ylabel('$y$ [code units]')
f.tight_layout()
f.savefig('trajectories.pdf')
f.savefig('trajectories.svg')
plt.close(f)
###########################################################################
###########################################################################
traj_index = 3
for neighbours, smoothness in neighbours_smoothness_list:
interp_name = 'I{0}O{1}'.format(2*neighbours+2, 2*smoothness+1)
cc = DNS(simname = 'nsvep_base_' + interp_name)
cc.compute_statistics()
xx_base = read_trajectory(cc, traj_index, iter0 = 0, dset_name = 'velocity')
iter0 = 16
tt_base = np.array(range(iter0, iter0 + xx_base.shape[0]))*cc.parameters['dt']
f = plt.figure()
ax_counter = 1
for factor in [1, 2, 4]:
a = f.add_subplot(310 + ax_counter)
ax_counter += 1
a.plot(tt_base,
xx_base[:, 2],
label = 'base_' + interp_name,
marker = '+',
linewidth = 0.5,
alpha = 0.2)
cc = DNS(simname = 'nsvep_{0}x_'.format(factor) + interp_name)
cc.compute_statistics(iter0 = 8*factor)
xx = read_trajectory(cc, traj_index, iter0 = 8*factor, dset_name = 'velocity')
tt = np.array(range(8*factor, xx.shape[0]+8*factor))*cc.parameters['dt']
a.plot(tt,
xx[:, 2],
label = interp_name + '$f = {0}$'.format(factor),
marker = 'x',
dashes = (4/factor, 3, 4/factor),
alpha = 0.2)
df = h5py.File(cc.simname + '_checkpoint_0.h5', 'r')
xx = df['tracers0/rhs/{0}'.format(8*factor)][()]
for rhs_index in [1, 2, 3]:
a.scatter((8*factor - rhs_index)*cc.parameters['dt'],
xx[rhs_index, traj_index, 2],
marker = '*',
zorder = -10,
alpha = 0.5,
color = 'black')
df.close()
a.legend(loc = 'best', fontsize = 7)
a.set_xlabel('$t$ [code units]')
a.set_ylabel('$v_y$ [code units]')
f.tight_layout()
f.savefig('velocities_{0}.pdf'.format(interp_name))
plt.close(f)
###########################################################################
return None
def read_trajectory(
cc,
traj_index,
iter0,
dset_name = 'position'):
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'][()]+1, cc.parameters['niter_part']):
xx.append(df['tracers0/' + dset_name + '/{0}'.format(ii)][traj_index])
df.close()
return np.array(xx)
Finally, the script ends with a standard call to main:
if __name__ == '__main__':
main()