adds alternative computations of R_LL for comparison

examples/correlation_length/test.py

@@ -3,6 +3,7 @@ from TurTLE import DNS, PP
 import os
 import h5py
 import matplotlib.pyplot as plt
+import numpy as np
 
 def read_correlation(
         group,
@@ -15,6 +16,111 @@ def read_correlation(
     Rij[1:] /= 2
     return Rij
 
+def compute_stat_possibilities(
+        cc,
+        iteration):
+    post_file = h5py.File(cc.simname + '_post.h5', 'r')
+    group = post_file['get_3D_correlations/correlations']
+    Rijx = read_correlation(group, components = '(0, 0, x)')
+    Rijy = read_correlation(group, components = '(0, y, 0)')
+    Rijz = read_correlation(group, components = '(z, 0, 0)')
+    post_file.close()
+    data_file = cc.get_data_file()
+    stat_index = iteration//cc.parameters['niter_stat']
+    # read latest slice data
+    cc.statistics['correlations/energy'] = data_file['statistics/moments/velocity'][
+              stat_index, 2, 3] / 2
+    cc.statistics['correlations/Uint'] = (2*(cc.statistics['correlations/energy']/3))**0.5
+    cc.statistics['correlations/R_LL'] = (Rijx[:, 0, 0] + Rijy[:, 1, 1] + Rijz[:, 2, 2])
+    cc.statistics['correlations/rho_LL'] = cc.statistics['correlations/R_LL'] / (2*cc.statistics['correlations/energy'])
+    cc.statistics['correlations/x'] = cc.get_coord('x')[:cc.parameters['nx']//2]
+    dx = cc.statistics['correlations/x'][1] - cc.statistics['correlations/x'][0]
+    cc.statistics['correlations/Lint'] = np.trapz(
+            cc.statistics['correlations/rho_LL'], dx = dx)
+    E_k = 0.5*(data_file['statistics/spectra/velocity_velocity'][stat_index, :, 0, 0]
+           + data_file['statistics/spectra/velocity_velocity'][stat_index, :, 1, 1]
+           + data_file['statistics/spectra/velocity_velocity'][stat_index, :, 2, 2])
+    # normalization factor is (4 pi * shell_width * kshell^2) / (nmodes in shell * dkx*dky*dkz)
+    norm_factor = (4*np.pi*cc.statistics['dk']*cc.statistics['kshell']**2) / (cc.parameters['dkx']*cc.parameters['dky']*cc.parameters['dkz']*cc.statistics['nshell'])
+    E_k_normalized = E_k * norm_factor
+    kval = np.round(cc.statistics['kshell'])
+    kval_normalized = cc.statistics['kshell']
+    # See Pope equations 6.211 and 6.216
+    for key in ['E11_k1_6.211_sum',
+                'E11_k1_6.211_trapz',
+                'E11_k1_6.211_normalized_sum',
+                'E11_k1_6.211_normalized_trapz',
+                'E11_k1_6.216_sum',
+                'E11_k1_6.216_trapz',
+                'E11_k1_6.216_normalized_sum',
+                'E11_k1_6.216_normalized_trapz']:
+        cc.statistics['correlations/' + key] = np.zeros_like(E_k)
+    for kk in range(kval.shape[0]-2):
+        qq = kval[kk] * cc.statistics['correlations/x']
+        cc.statistics['correlations/E11_k1_6.211_sum'][kk] = (
+                (2.0 / np.pi)
+              * np.sum(cc.statistics['correlations/R_LL'] * np.cos(qq))*dx)
+        cc.statistics['correlations/E11_k1_6.211_trapz'][kk] = (
+                (2.0 / np.pi)
+              * np.trapz(cc.statistics['correlations/R_LL'] * np.cos(qq), dx = dx))
+        qq = kval_normalized[kk] * cc.statistics['correlations/x']
+        cc.statistics['correlations/E11_k1_6.211_normalized_sum'][kk] = (
+                (2.0 / np.pi)
+              * np.sum(cc.statistics['correlations/R_LL'] * np.cos(qq))*dx)
+        cc.statistics['correlations/E11_k1_6.211_normalized_trapz'][kk] = (
+                (2.0 / np.pi)
+              * np.trapz(cc.statistics['correlations/R_LL'] * np.cos(qq), dx = dx))
+        prefactor = 3 # because R_LL = 3 * R_11
+        cc.statistics['correlations/E11_k1_6.216_trapz'][kk] = prefactor*np.trapz(
+                (E_k[kk+1:-2] / kval[kk+1:-2]) * (1 - kval[kk]**2 / kval[kk+1:-2]**2),
+                dx = cc.statistics['dk'])
+        cc.statistics['correlations/E11_k1_6.216_sum'][kk] = prefactor*np.sum(
+                (E_k[kk+1:-2] / kval[kk+1:-2]) * (1 - kval[kk]**2 / kval[kk+1:-2]**2))
+        cc.statistics['correlations/E11_k1_6.216_normalized_trapz'][kk] = prefactor*np.trapz(
+                (E_k_normalized[kk+1:-2] / kval_normalized[kk+1:-2])
+                * (1 - kval_normalized[kk]**2 / kval_normalized[kk+1:-2]**2),
+                kval_normalized[kk+1:-2])
+        cc.statistics['correlations/E11_k1_6.216_normalized_sum'][kk] = prefactor*np.sum(
+                (E_k_normalized[kk+1:-2] / kval_normalized[kk+1:-2])
+                * (1 - kval_normalized[kk]**2 / kval_normalized[kk+1:-2]**2))
+    for key in ['R_LL_from_spectrum_sum',
+                'R_LL_from_spectrum_trapz',
+                'R_LL_from_normalized_spectrum_sum',
+                'R_LL_from_normalized_spectrum_trapz']:
+        cc.statistics['correlations/' + key] = np.zeros_like(
+            cc.statistics['correlations/R_LL'])
+        cc.statistics['correlations/' + key][0] = np.sum(E_k)*2
+    for ii in range(1, cc.statistics['correlations/R_LL'].shape[0]):
+        rr = cc.statistics['correlations/x'][ii]
+        qq = kval[1:-2] * rr
+        cc.statistics['correlations/R_LL_from_spectrum_sum'][ii] = 2 * 3 * np.sum(
+                E_k[1:-2]
+                * (np.sin(qq) - qq * np.cos(qq))
+                / qq**3)
+        cc.statistics['correlations/R_LL_from_spectrum_trapz'][ii] = 2 * 3 * np.trapz(
+                E_k[1:-2]
+                * (np.sin(qq) - qq * np.cos(qq))
+                / qq**3,
+                dx = cc.statistics['dk'])
+        qq = kval_normalized[1:-2] * rr
+        cc.statistics['correlations/R_LL_from_normalized_spectrum_sum'][ii] = 2 * 3 * np.sum(
+                E_k_normalized[1:-2]
+                * (np.sin(qq) - qq * np.cos(qq))
+                / qq**3)
+        cc.statistics['correlations/R_LL_from_normalized_spectrum_trapz'][ii] = 2 * 3 * np.trapz(
+                E_k_normalized[1:-2]
+                * (np.sin(qq) - qq * np.cos(qq))
+                / qq**3,
+                kval_normalized[1:-2])
+    kshell = cc.statistics['kshell']
+    cc.statistics['correlations/Lint_from_spectrum'] = (
+            (np.pi /
+              (2*cc.statistics['correlations/Uint']**2)) *
+             np.trapz(E_k_normalized[:-2] /
+                      np.maximum(1, kshell[:-2]),
+                      kshell[:-2]))
+    return None
+
 def main():
     simname = 'fbM_N0128_kMeta1.5_tf2'
     if not os.path.exists(simname + '.h5'):
@@ -28,6 +134,38 @@ def main():
             '--iter0', '8192',
             '--iter1', '8192'])
 
+    c = DNS(work_dir = './',
+            simname = simname)
+    c.compute_statistics()
+    compute_stat_possibilities(c, 8192)
+
+    f = plt.figure()
+    a = f.add_subplot(111)
+    a.plot(c.statistics['correlations/x'],
+           c.statistics['correlations/R_LL'])
+    for kk in ['R_LL_from_spectrum_sum',
+               'R_LL_from_spectrum_trapz',
+               'R_LL_from_normalized_spectrum_sum',
+               'R_LL_from_normalized_spectrum_trapz']:
+        a.plot(c.statistics['correlations/x'],
+               c.statistics['correlations/' + kk])
+    f.savefig('test.pdf')
+    plt.close(f)
+    return None
+
+def main0():
+    simname = 'fbM_N0128_kMeta1.5_tf2'
+    if not os.path.exists(simname + '.h5'):
+        print('please download ' + simname + ' data from shared folder')
+    if not os.path.exists(simname + '_post.h5'):
+        c = PP()
+        c.launch([
+            'get_3D_correlations',
+            '--simname', simname,
+            '--full_snapshot_output', '1',
+            '--iter0', '8192',
+            '--iter1', '8192'])
+
     data_file = h5py.File(simname + '_post.h5', 'r')
     group = data_file['get_3D_correlations/correlations']
     f = plt.figure()