Feature/particle integration test

Adds an in-depth convergence test for particle tracking of a stationary field.

In the below figure "o" stands for "solver order", "n" for interpolation neighbours and "m" for interpolation smoothness.

Notes:

• I also added a classic Runge Kutta solver for particles, to emphasize effects of interpolation. The plot suggests that the implementation is correct (I need to see order 4 convergence to believe the code is correct).
• The number of particles is relatively small now because otherwise the test is too slow in DEBUG mode. There seems to be a separate issue with optimization that I'll need to address.
• The effects of the interpolation are only visible if the particle trajectories are long enough --- i.e. particles need to traverse a few grid cells in order for cell-crossing effects to show up. This is why there are so many iterations...

TurTLE/test/test_particle_integration.py

+import os
+import numpy as np
+import h5py
+import matplotlib.pyplot as plt
+
+from TurTLE import TEST
+
+def main():
+    c = TEST()
+    if not os.path.exists('test.h5'):
+        c.launch([
+            'test_particle_integration',
+            '--np', '4',
+            '--ntpp', '3'])
+
+    data_file = h5py.File(c.simname + '_particles.h5', 'r')
+
+    # to figure out reasonable timestep to satisfy "particle shall not move more than one grid cell"
+    #vv = data_file['particle_o2_n1_m1_f1/velocity/0'][()]
+    #print(np.min(vv), np.max(vv))
+
+
+    f = plt.figure(figsize = (9, 9))
+    a = f.add_subplot(111)
+
+    for ff in [1, 2, 4, 8]:
+        gg = data_file['particle_o2_n2_m1_f{0}/position'.format(ff)]
+        xx = read_trajectory(gg, 5, 7)
+        a.plot(xx[:, 0], xx[:, 2], marker = '.')
+    f.tight_layout()
+    f.savefig('test_particle_integration_trajectory.pdf')
+    plt.close(f)
+
+    f = plt.figure(figsize = (9, 9))
+    a = f.add_subplot(111)
+    dt0 = c.get_data_file()['parameters/dt'][()]
+    dtlist = [dt0, dt0/2, dt0/4]
+    for n, m in [(0, 0),
+                 (1, 1),
+                 (2, 1),
+                 (3, 2)]:
+        for oo in [1, 2, 4]:
+            err_list = []
+            for factor in [1, 2, 4]:
+                xx1 = data_file['particle_o{0}_n{1}_m{2}_f{3}/position/{4}'.format(oo, n, m, factor, factor*8)][()]
+                xx2 = data_file['particle_o{0}_n{1}_m{2}_f{3}/position/{4}'.format(oo, n, m, factor*2, factor*16)][()]
+                diff = np.sqrt(np.sum((xx2 - xx1)**2, axis = -1))
+                err_list.append(np.mean(diff))
+            a.plot(dtlist,
+               err_list,
+               label = 'o{0}_n{1}_m{2}'.format(oo, n, m),
+               dashes = (oo, oo),
+               marker = '.')
+    a.plot(dtlist,
+           np.array(dtlist)**2 * (1e-2),
+           color = (0, 0, 0),
+           dashes = (1, 1),
+           label = '$\\propto \\Delta t^2$')
+    a.plot(dtlist,
+           np.array(dtlist)**3 * (1e-2),
+           color = (0, 0, 0),
+           dashes = (2, 1, 2),
+           label = '$\\propto \\Delta t^3$')
+    a.plot(dtlist,
+           np.array(dtlist)**4 * (1e-2),
+           color = (0, 0, 0),
+           dashes = (3, 5, 3),
+           label = '$\\propto \\Delta t^4$')
+    a.legend(loc = 'upper left')
+    a.set_xlabel('$\\Delta t$')
+    a.set_ylabel('mean trajectory error')
+    a.set_xscale('log')
+    a.set_yscale('log')
+    f.tight_layout()
+    f.savefig('test_particle_integration_evdt.pdf')
+    f.savefig('test_particle_integration_evdt.svg')
+    plt.close(f)
+
+    #f = plt.figure(figsize = (6, 6))
+    #a = f.add_subplot(111)
+    #err_mean_1 = []
+    #err_max_1 = []
+    #n = 1
+    #phi0 = data_file['particle/phi_{0}{1}/0'.format(n  , 1)][()]
+    #phi1 = data_file['particle/phi_{0}{1}/0'.format(n+1, 1)][()]
+    #err = np.abs(phi0 - phi1)
+    #err_mean_1.append(err.mean())
+    #err_max_1.append(err.max())
+    #err_mean_2 = []
+    #err_max_2 = []
+    #for n in [2, 3, 4]:
+    #    phi0 = data_file['particle/phi_{0}{1}/0'.format(n  , 1)][()]
+    #    phi1 = data_file['particle/phi_{0}{1}/0'.format(n+1, 1)][()]
+    #    err = np.abs(phi0 - phi1)
+    #    err_mean_1.append(err.mean())
+    #    err_max_1.append(err.max())
+    #    phi0 = data_file['particle/phi_{0}{1}/0'.format(n  , 2)][()]
+    #    phi1 = data_file['particle/phi_{0}{1}/0'.format(n+1, 2)][()]
+    #    err = np.abs(phi0 - phi1)
+    #    err_mean_2.append(err.mean())
+    #    err_max_2.append(err.max())
+    #neighbour_list = np.array([1, 2, 3, 4]).astype(np.float)
+    #a.plot(neighbour_list, err_mean_1, marker = '.', label = 'm=1, mean')
+    #a.plot(neighbour_list, err_max_1,  marker = '.', label = 'm=1, max')
+    #a.plot(neighbour_list, 1e-4*neighbour_list**(-4),  color = 'black', dashes = (2, 2))
+    #neighbour_list = np.array([2, 3, 4]).astype(np.float)
+    #a.plot(neighbour_list, err_mean_2, marker = '.', label = 'm=2, mean')
+    #a.plot(neighbour_list, err_max_2,  marker = '.', label = 'm=2, max')
+
+    #a.set_xscale('log')
+    #a.set_yscale('log')
+    #a.legend(loc = 'best')
+
+    #f.tight_layout()
+    #f.savefig('test_interpolation_methods_err.pdf')
+    #plt.close(f)
+
+    #data_file.close()
+    return None
+
+def read_trajectory(group, ix, iz):
+    iter_names = group.keys()
+    iterations = np.sort([int(ii) for ii in iter_names])
+
+    xx = []
+    for ii in iterations:
+        xx.append(group['{0}'.format(ii)][iz, ix])
+    return np.array(xx)
+
+if __name__ == '__main__':
+    main()
+