tweaks PP file output

TurTLE/PP.py

@@ -783,9 +783,9 @@ class PP(_code):
         if 'field_dtype' in df.keys():
             # we don't need to do anything, raw binary files are used
             return None
-        last_iteration = df['iteration'][...]
-        cppf = df['parameters/checkpoints_per_file'][...]
-        niter_out = df['parameters/niter_out'][...]
+        last_iteration = df['iteration'][()]
+        cppf = df['parameters/checkpoints_per_file'][()]
+        niter_out = df['parameters/niter_out'][()]
         with h5py.File(os.path.join(self.work_dir, self.simname + '_fields.h5'), 'a') as ff:
             ff.require_group('vorticity')
             ff.require_group('vorticity/complex')
@@ -807,7 +807,13 @@ class PP(_code):
         if self.dns_type == 'resize':
             with h5py.File(os.path.join(self.work_dir, self.pp_parameters['new_simname'] + '.h5'), 'a') as ff:
                 for kk in df['parameters'].keys():
-                    ff['parameters/' + kk] = df['parameters/' + kk][...]
+                    ff['parameters/' + kk] = df['parameters/' + kk][()]
+                ff['parameters/nx'][()] = self.pp_parameters['new_nx']
+                ff['parameters/ny'][()] = self.pp_parameters['new_ny']
+                ff['parameters/nz'][()] = self.pp_parameters['new_nz']
+                if 'iteration' not in ff.keys():
+                    ff['iteration'] = last_iteration
+                    ff['checkpoint'] = df['checkpoint'][()]
         return None
     def launch_jobs(
             self,

documentation/sphinx_static/convergence.py

+import numpy as np
+import h5py
+import matplotlib.pyplot as plt
+
+import TurTLE
+from TurTLE import DNS, PP
+
+base_niterations = 256
+test_niterations = 32
+
+def run_simulations():
+    # change these two values as neded.
+    # number of MPI processes to use
+    nprocesses = 4
+    # number of OpenMP threads per MPI process to use
+    nthreads_per_process = 2
+
+# 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_1x',
+            '--niter_todo', '{0}'.format(base_niterations),
+            '--niter_out', '{0}'.format(base_niterations),
+            '--overwrite-src-parameters',
+            '--kMeta',  '1.0',
+            '--niter_stat', '16',
+            '--checkpoints_per_file', '{0}'.format(16)])
+
+# 2. Generate initial conditions for NSVE runs at 2x and 4x the
+#    base resolution.
+    for factor in [2, 4]:
+        # create postprocess object
+        cc = PP()
+        # create a larger version of the field at iteration 512
+        cc.launch([
+                'resize',
+                '--np', '{0}'.format(nprocesses),
+                '--ntpp', '{0}'.format(nthreads_per_process),
+                '--simname', 'base_1x',
+                '--precision', 'double',
+                '--iter0', '{0}'.format(base_niterations),
+                '--iter1', '{0}'.format(base_niterations),
+                '--new_nx', '{0}'.format(factor*32),
+                '--new_ny', '{0}'.format(factor*32),
+                '--new_nz', '{0}'.format(factor*32),
+                '--new_simname', 'base_{0}x'.format(factor)])
+        # update "checkpoint" file, so we can use "base__x" as "src-simname"
+        # this is needed after "resize"
+        f1 = h5py.File('base_{0}x_checkpoint_0.h5'.format(factor), 'a')
+        f1['vorticity/complex/{0}'.format(base_niterations)] = h5py.ExternalLink(
+            'base_{0}x_fields.h5'.format(factor),
+            'vorticity/complex/{0}'.format(base_niterations))
+        f1.close()
+
+# 3. Generate initial conditions for NSE runs at 1x, 2x and 4x the
+#    base resolution.
+    for factor in [1, 2, 4]:
+        # create postprocess object
+        cc = PP()
+        # compute velocity field at last iteration
+        cc.launch([
+                'get_velocity',
+                '--np', '{0}'.format(nprocesses),
+                '--ntpp', '{0}'.format(nthreads_per_process),
+                '--simname', 'base_{0}x'.format(factor),
+                '--precision', 'double',
+                '--iter0', '{0}'.format(base_niterations),
+                '--iter1', '{0}'.format(base_niterations)])
+        # we don't need to update "checkpoint" file after "get_velocity"
+
+# 4. 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(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', 'nsve_{0}x'.format(factor),
+                '--precision', 'double',
+                '--dtfactor', '0.3',
+                '--kMeta', '{0}'.format(1.0*factor),
+                '--niter_todo', '{0}'.format(test_niterations*factor),
+                '--niter_out', '{0}'.format(test_niterations*factor),
+                '--niter_stat', '{0}'.format(4*factor)])
+
+# 5. 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(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', 'nse_{0}x'.format(factor),
+                '--precision', 'double',
+                '--dtfactor', '0.3',
+                '--kMeta', '{0}'.format(1.0*factor),
+                '--niter_todo', '{0}'.format(test_niterations*factor),
+                '--niter_out', '{0}'.format(test_niterations*factor),
+                '--niter_stat', '{0}'.format(4*factor)])
+    return None
+
+def ift(cfield):
+    ff = np.fft.irfftn(cfield, axes = (0, 1, 2))
+    field = np.transpose(ff, (1, 0, 2, 3)).copy()
+    return field
+
+def compute_curl(
+        cc,
+        cfield,
+        inverse_curl = False):
+    kx = cc.get_data_file()['kspace/kx'][()]
+    ky = cc.get_data_file()['kspace/ky'][()]
+    kz = cc.get_data_file()['kspace/kz'][()]
+    ff = cfield.copy()
+    ff[..., 0] = 1j*(ky[:, None, None]*cfield[..., 2] - kz[None, :, None]*cfield[..., 1])
+    ff[..., 1] = 1j*(kz[None, :, None]*cfield[..., 0] - kx[None, None, :]*cfield[..., 2])
+    ff[..., 2] = 1j*(kx[None, None, :]*cfield[..., 1] - ky[:, None, None]*cfield[..., 0])
+    if inverse_curl:
+        k2 = (kx[None, None,    :]**2 +
+              ky[   :, None, None]**2 +
+              kz[None,    :, None]**2)
+        k2[0, 0, 0] = 1.
+        ff /= k2[:, :, :, None]
+    return ff
+
+def get_velocity(
+        cc,
+        iteration = 0):
+    data_file = h5py.File(
+            cc.get_checkpoint_fname(iteration = iteration),
+            'r')
+    if 'velocity' in data_file.keys():
+        if 'real' in data_file['velocity'].keys():
+            vel = data_file['velocity/real/{0}'.format(iteration)][()]
+        else:
+            cvel = data_file['velocity/complex/{0}'.format(iteration)][()]
+            vel = ift(cvel)
+            #vel *= cc.parameters['nx']*cc.parameters['ny']*cc.parameters['nz']
+    else:
+        assert('complex' in data_file['vorticity'].keys())
+        cvort = data_file['vorticity/complex/{0}'.format(iteration)][()]
+        cvel = compute_curl(cc, cvort, inverse_curl = True)
+        vel = ift(cvel)
+    data_file.close()
+    return vel
+
+def plot_error():
+    c0_list = [DNS(simname = 'nsve_{0}x'.format(factor))
+               for factor in [1, 2, 4]]
+    c1_list = [DNS(simname = 'nse_{0}x'.format(factor))
+               for factor in [1, 2, 4]]
+    for cc in c0_list + c1_list:
+        cc.compute_statistics()
+
+    factor_list = [1, 2, 4]
+
+    error_list = []
+    for ii in range(len(factor_list)):
+        factor = factor_list[ii]
+        c0 = c0_list[ii]
+        c1 = c1_list[ii]
+        df0 = h5py.File(c0.get_checkpoint_fname(iteration = test_niterations*factor), 'r')
+        df1 = h5py.File(c1.get_checkpoint_fname(iteration = test_niterations*factor), 'r')
+        vel0 = get_velocity(c0, iteration = test_niterations*factor)
+        vel1 = get_velocity(c1, iteration = test_niterations*factor)
+        diff = vel1 - vel0
+        vv = np.sqrt(np.sum(vel0**2, axis = 3))
+        dd = np.sqrt(np.sum(diff**2, axis = 3))
+        err = np.max(dd / vv)
+        error_list.append(err)
+        print('maximum error for factor {0} is {1}'.format(
+            factor, err))
+        df0.close()
+        df1.close()
+
+    f = plt.figure(figsize = (4, 4))
+    a = f.add_subplot(111)
+    fl = np.array(factor_list).astype(np.float)
+    a.plot(fl,
+           error_list,
+           marker = '.',
+           label = 'NSVE vs NSE')
+    a.plot(fl, fl**(-3),
+           dashes = (1, 1),
+           color = 'black',
+           label = '$\propto f^{-3}$')
+    a.set_ylabel('relative error')
+    a.set_xlabel('resolution factor $f$')
+    a.set_xscale('log')
+    a.set_yscale('log')
+    a.legend(loc = 'best')
+    f.tight_layout()
+    f.savefig('field_err_vs_dt.pdf')
+    plt.close(f)
+    return None
+
+def main():
+    run_simulations()
+    plot_error()
+    return None
+
+if __name__ == '__main__':
+    main()
+