# HG changeset patch
# User calpe  <miguel.calpe.linares@gmail.com>
# Date 1538639034 -7200
#      Thu Oct 04 09:43:54 2018 +0200
# Node ID 14d6d55c07bd845dad26f17fd9de3763e9b00f7d
# Parent  e6653ef320349320e4422161eb970d036b8faec5
Scripts to compute and plot buoyancy Reynolds and horizontal Froude.

diff --git a/scripts/ns2d.strat/compute_anisotropy.py b/scripts/ns2d.strat/compute_anisotropy.py
new file mode 100644
--- /dev/null
+++ b/scripts/ns2d.strat/compute_anisotropy.py
@@ -0,0 +1,82 @@
+"""
+compute_anisotropy.py
+=====================
+1/10/2018
+
+Function to compute the anisotropy of a simulation.
+
+"""
+from glob import glob
+
+import h5py
+import numpy as np
+
+def _compute_array_times_from_path(path_simulation):
+    """
+    Compute array with times from path simulation.
+
+    Parameters
+    ----------
+    path_simulation : str
+      Path of the simulation.
+    """
+    times_phys_files = []
+
+    paths_phys_files = glob(path_simulation + "/state_phys_t*")
+    for path in paths_phys_files:
+        if not "=" in path.split("state_phys_t")[1].split(".nc")[0]:
+            times_phys_files.append(
+                float(path.split("state_phys_t")[1].split(".nc")[0]))
+        else:
+            continue
+    return np.asarray(times_phys_files)
+
+def compute_anisotropy(path_simulation, tmin=None):
+    """
+    It computes the anisotropy of a simulation.
+
+    Parameters
+    ----------
+    path_simulation : str
+      Path of the simulation.
+
+    tmin : float
+      Lower limit to compute the time average.
+      By default it takes the last 10 files.
+    """
+
+    # Print out
+    res_out = float(path_simulation.split("NS2D.strat_")[1].split("x")[0])
+    gamma_str = path_simulation.split("_gamma")[1].split("_")[0]
+    if gamma_str.startswith("0"):
+        gamma_out = float(gamma_str[0] + "." + gamma_str[1])
+    else:
+        gamma_out = float(gamma_str)
+
+    print("Compute anisotropy nx = {} and gamma {}..".format(res_out, gamma_out))
+
+    # Compute index start average time imin.
+    times = _compute_array_times_from_path(path_simulation)
+    dt_state_phys = np.median(np.diff(times))
+
+    if not tmin:
+        nb_files = 10
+        tmin = np.max(times) - (nb_files * dt_state_phys)
+    imin = np.argmin(abs(times - tmin))
+
+    # Compute anisotropy
+    anisotropies = []
+    for path in glob(path_simulation + "/state_phys_t*")[imin:]:
+        with h5py.File(path, "r") as f:
+            ux = f["state_phys"]["ux"].value
+            uz = f["state_phys"]["uy"].value
+        anisotropies.append(np.mean(ux**2 / (ux**2 + uz**2)))
+
+    return np.mean(anisotropies)
+
+if __name__ == "__main__":
+   path_simulation = ("/fsnet/project/meige/2015/15DELDUCA/DataSim/" +
+                      "sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma1_2018-08-14_10-01-22")
+
+   anisotropy = compute_anisotropy(path_simulation)
+   print("anisotropy = {}".format(anisotropy))
diff --git a/scripts/ns2d.strat/compute_flow_features.py b/scripts/ns2d.strat/compute_flow_features.py
new file mode 100644
--- /dev/null
+++ b/scripts/ns2d.strat/compute_flow_features.py
@@ -0,0 +1,92 @@
+"""
+compute_flow_features.py
+========================
+28/09/2018
+
+"""
+import os
+import numpy as np
+import h5py
+
+from glob import glob
+import matplotlib.pyplot as plt
+
+from fluidsim import load_params_simul
+
+# Argparse arguments
+nx = 3840
+n_files_average = 50 # Number of files to perform the time average
+
+# Create paths
+path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
+
+if nx == 1920:
+    directory = "sim{}_no_shear_modes".format(nx)
+elif nx == 3840 or nx == 7680:
+    directory = "sim{}_modif_res_no_shear_modes".format(nx)
+else:
+    raise ValueError(".")
+
+path_simulations = sorted(glob(os.path.join(path_root, directory, "NS2D*")))
+
+gammas = []
+anisotropies_gammas = []
+ratio_dissipations = []
+for ipath, path in enumerate(path_simulations):
+    params = params = load_params_simul(path)
+
+    # Add gamma to list gammas
+    gamma_str = path.split("_gamma")[1].split("_")[0]
+    if gamma_str.startswith("0"):
+        gammas.append(float(gamma_str[0] + "." + gamma_str[1]))
+    else:
+        gammas.append(float(gamma_str))
+
+    # Compute time average ratio ux**2 / uy**2 (anisotropy)
+    print("Computing anisotropy for gamma {}...".format(gammas[ipath]))
+    path_phys_files = glob(path + "/state_phys_t*")
+    anisotropies = []
+    for path_file in path_phys_files[-n_files_average:]:
+        with h5py.File(path_file, "r") as f:
+            ux = f["state_phys"]["ux"].value
+            uz = f["state_phys"]["uy"].value
+            anisotropies.append(np.mean(ux**2) / np.mean(uz**2))
+    anisotropies_gammas.append(np.mean(anisotropies))
+
+    # Compute ratio D(k_x)/epsilon
+    print("Computing ratio dissipation for gamma {}...".format(gammas[ipath]))
+    with h5py.File(path + "/spect_energy_budg.h5", "r") as f:
+        kx = f['kxE'].value
+        kz = f['kyE'].value
+        dset_dissEKu_kx = f['dissEKu_kx']
+        dset_dissEKv_kx = f['dissEKv_kx']
+        dset_dissEA_kx = f['dissEA_kx']
+
+        delta_kx = kx[1] - kx[0]
+        delta_kz = kz[1] - kz[0]
+
+        dissEK_kx = (dset_dissEKu_kx[-n_files_average:] + \
+                     dset_dissEKv_kx[-n_files_average:])
+        dissEA_kx = dset_dissEA_kx[-n_files_average:]
+        dissE_kx = (dissEK_kx + dissEA_kx).mean(0)
+        D_kx = dissE_kx.cumsum() * delta_kx
+
+        # Compute k_fx
+        k_fx = (np.sin(params.forcing.tcrandom_anisotropic.angle) *
+                params.forcing.nkmax_forcing * max(delta_kx, delta_kz))
+        ik_fx = np.argmin(abs(kx - k_fx ))
+
+        # Compute ratio
+        ratio_dissipations.append(D_kx[ik_fx] / D_kx[-1])
+
+fig1, ax1 = plt.subplots()
+ax1.set_xlabel(r"$\gamma$")
+ax1.set_ylabel(r"$U_x^2/U_z^2$")
+ax1.plot(gammas, anisotropies_gammas, 'ro')
+
+fig2, ax2 = plt.subplots()
+ax2.set_xlabel(r"$\gamma$")
+ax2.set_ylabel(r"$D(k_{fx})/D(k_{x, max})$")
+ax2.plot(gammas, ratio_dissipations, 'bo')
+
+plt.show()
diff --git a/scripts/ns2d.strat/compute_ratio_dissipation.py b/scripts/ns2d.strat/compute_ratio_dissipation.py
new file mode 100644
--- /dev/null
+++ b/scripts/ns2d.strat/compute_ratio_dissipation.py
@@ -0,0 +1,57 @@
+"""
+compute_ratio_dissipation.py
+============================
+02/10/2018
+"""
+
+import h5py
+import numpy as np
+
+from fluidsim import load_params_simul
+
+def compute_ratio_dissipation(path_simulation, tmin=None):
+    """
+    Compute ratio dissipation from path simulation.
+    """
+    # Print out
+    res_out = float(path_simulation.split("NS2D.strat_")[1].split("x")[0])
+    gamma_str = path_simulation.split("_gamma")[1].split("_")[0]
+    if gamma_str.startswith("0"):
+        gamma_out = float(gamma_str[0] + "." + gamma_str[1])
+    else:
+        gamma_out = float(gamma_str)
+
+    print("Compute dissipation nx = {} and gamma {}..".format(res_out, gamma_out))
+
+    # Load object parameters
+    params = load_params_simul(path_simulation)
+
+    with h5py.File(path_simulation + "/spect_energy_budg.h5", "r") as f:
+        times = f["times"].value
+        kx = f["kxE"].value
+        kz = f['kyE'].value
+        dset_dissEKu_kx = f['dissEKu_kx']
+        dset_dissEKv_kx = f['dissEKv_kx']
+        dset_dissEA_kx = f["dissEA_kx"]
+
+        # Compute itmin time average
+        if not tmin:
+            nb_files = 10
+            dt = np.median(np.diff(times))
+            tmin = np.max(times) - (nb_files * dt)
+        itmin = np.argmin(abs(times - tmin))
+
+        # Compute dissipation curve
+        delta_kx = np.median(np.diff(kx))
+        delta_kz = np.median(np.diff(abs(kz)))
+        dissEK_kx = dset_dissEKu_kx[-itmin:] + dset_dissEKv_kx[-itmin:]
+        dissEA_kx = dset_dissEA_kx[-itmin:]
+        dissE_kx = (dissEK_kx + dissEA_kx).mean(0)
+        D_kx = dissE_kx.cumsum() * delta_kx
+
+    # Compute k_fx
+    k_fx = (np.sin(params.forcing.tcrandom_anisotropic.angle) *
+            params.forcing.nkmax_forcing * max(delta_kx, delta_kz))
+    ik_fx = np.argmin(abs(kx - k_fx ))
+
+    return D_kx[ik_fx] / D_kx[-1]
diff --git a/scripts/ns2d.strat/compute_reynolds_froude.py b/scripts/ns2d.strat/compute_reynolds_froude.py
new file mode 100644
--- /dev/null
+++ b/scripts/ns2d.strat/compute_reynolds_froude.py
@@ -0,0 +1,107 @@
+"""
+compute_reynolds_froude.py
+==========================
+1/10/2018
+
+"""
+
+import h5py
+import numpy as np
+
+from fluidsim import load_params_simul
+
+def _compute_epsilon_from_path(path_simulation, tmin=None):
+    """
+    Computes the mean dissipation from tmin
+    """
+    # Load data dissipation
+    with open(path_simulation + "/spatial_means.txt", "r") as f:
+        lines = f.readlines()
+
+    lines_t = []
+    lines_epsK = []
+
+    for il, line in enumerate(lines):
+        if line.startswith("time ="):
+            lines_t.append(line)
+        if line.startswith("epsK ="):
+            lines_epsK.append(line)
+
+    nt = len(lines_t)
+    t = np.empty(nt)
+    epsK = np.empty(nt)
+
+    for il in range(nt):
+            line = lines_t[il]
+            words = line.split()
+            t[il] = float(words[2])
+
+            line = lines_epsK[il]
+            words = line.split()
+            epsK[il] = float(words[2])
+
+    # Compute start time average itmin
+    dt_spatial = np.median(np.diff(t))
+
+    if not tmin:
+        nb_files = 100
+        tmin = np.max(t) - (dt_spatial * nb_files)
+
+    itmin = np.argmin((abs(t - tmin)))
+
+    return np.mean(epsK[itmin:], axis=0)
+
+def _compute_lx_from_path(path_simulation, tmin=None):
+    """
+    Compute horizontal length from path using appendix B. Brethouwer 2007
+    """
+
+    # Load parameters from simulation
+    params = load_params_simul(path_simulation)
+
+    with h5py.File(path_simulation + "/spectra1D.h5", "r") as f:
+        times_spectra = f["times"].value
+        kx = f["kyE"].value
+        spectrum1Dkx_EK_ux = f["spectrum1Dkx_EK_ux"].value
+
+    # Compute time average spectra
+    dt = np.median(np.diff(times_spectra))
+
+    if not tmin:
+        nb_files = 100
+        tmin = np.max(times_spectra) - (dt * nb_files)
+    itmin = np.argmin(abs(times_spectra - tmin))
+
+    spectrum1Dkx_EK_ux = np.mean(spectrum1Dkx_EK_ux[-itmin:, :], axis=0)
+
+    # Remove modes dealiased
+    ikxmax = np.argmin(abs(kx - (np.max(kx) * params.oper.coef_dealiasing)))
+    kx = kx[:ikxmax]
+    spectrum1Dkx_EK_ux = spectrum1Dkx_EK_ux[:ikxmax]
+    delta_kx = np.median(np.diff(kx))
+
+    # Compute horizontal length scale Brethouwer 2007
+    return (np.sum(spectrum1Dkx_EK_ux * delta_kx) /
+            np.sum(kx * spectrum1Dkx_EK_ux * delta_kx))
+
+def compute_buoyancy_reynolds(path_simulation, tmin=None):
+    """
+    Compute the buoyancy Reynolds number.
+    """
+    params = load_params_simul(path_simulation)
+    epsK = _compute_epsilon_from_path(path_simulation, tmin=tmin)
+    lx = _compute_lx_from_path(path_simulation, tmin=tmin)
+
+    F_h = ((epsK / lx**2)**(1/3)) * (1 / params.N)
+
+    eta_8 = (params.nu_8 ** 3 / epsK)**(1/22)
+    Re_8 = (lx/eta_8)**(22/3)
+
+    R_b = Re_8 * F_h**8
+    return F_h, Re_8, R_b
+
+# path_simulation = "/fsnet/project/meige/2015/15DELDUCA/DataSim/sim1920_no_shear_modes/NS2D.strat_1920x480_S2pix1.571_F07_gamma02_2018-08-14_09-59-55"
+# F_h, Re_8, R_b = compute_buoyancy_reynolds(path_simulation)
+# print("F_h", F_h)
+# print("Re_8", Re_8)
+# print("R_b", R_b)
diff --git a/scripts/ns2d.strat/make_table_parameters.py b/scripts/ns2d.strat/make_table_parameters.py
new file mode 100644
--- /dev/null
+++ b/scripts/ns2d.strat/make_table_parameters.py
@@ -0,0 +1,101 @@
+"""
+make_table_parameters.py
+=========================
+28/09/2018
+
+"""
+import os
+import h5py
+import numpy as np
+
+from glob import glob
+
+from fluidsim import load_sim_for_plot
+
+# Argparse arguments
+nx = 1920
+MAKE_TABLE = False
+
+# Parameters script
+n_files_tmean = 100
+
+# Create path
+path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
+
+if nx == 1920:
+    directory = "sim1920_no_shear_modes"
+
+path_simulations = sorted(glob(os.path.join(path_root, directory, "NS2D*")))
+
+if MAKE_TABLE:
+    path_table = ("/home/users/calpelin7m/" +
+                  "Phd/docs/Manuscript/buoyancy_reynolds_table_n{}.tex".format(nx))
+
+    to_print = ("\\begin{table}[h]\n"
+                "\\centering \n"
+                "\\begin{tabular}{cccc} \n"
+                "\\toprule[1.5pt] \n" + \
+                "\\bm{$\gamma$} & \\bm{$F_h$} & \\bm{$Re_8$} & \\bm{$\mathcal{R}$} \\\\ \n"
+                "\\midrule\ \n")
+
+for path in path_simulations:
+
+    # Load object simulations
+    sim = load_sim_for_plot(path)
+
+    # Compute gamma from path
+    gamma_str = path.split("_gamma")[1].split("_")[0]
+    if gamma_str.startswith("0"):
+        gamma_table = float(gamma_str[0] + "." + gamma_str[1])
+    else:
+        gamma_table = float(gamma_str)
+
+    # Compute mean kinetic dissipation
+    dict_spatial = sim.output.spatial_means.load()
+
+    times = dict_spatial["t"]
+    epsK_tot = dict_spatial["epsK_tot"]
+    epsK_tmean = np.mean(epsK_tot[-n_files_tmean:], axis=0)
+    print("epsilon", epsK_tmean)
+    # Compute horizontal scale as Appendix B. Brethouwer (2007)
+    path_spectra = path + "/spectra1D.h5"
+
+    with h5py.File(path_spectra, "r") as f:
+        times_spectra = f["times"].value
+        kx = f["kyE"].value
+        spectrum1Dkx_EK_ux = f["spectrum1Dkx_EK_ux"].value
+
+    spectrum1Dkx_EK_ux = np.mean(spectrum1Dkx_EK_ux[-100:, :], axis=0)
+    ## Remove modes with dealiasing
+    ikxmax = np.argmin(abs(kx - sim.oper.kxmax_dealiasing))
+    kx = kx[:ikxmax]
+    spectrum1Dkx_EK_ux = spectrum1Dkx_EK_ux[:ikxmax]
+    delta_kx = sim.oper.deltakx
+    lx = (np.sum(spectrum1Dkx_EK_ux * delta_kx) /
+          np.sum(kx * spectrum1Dkx_EK_ux * delta_kx))
+    print("lx", lx)
+    # Compute eta_8
+    eta_8 = (sim.params.nu_8 ** 3 / epsK_tmean)**(1/22)
+    print("eta_8", eta_8)
+
+    # Compute Re_8
+    Re_8 = (lx/eta_8)**(22/3)
+    print("Re_8", Re_8)
+
+    # Compute horizontal Froude
+    F_h = ((epsK_tmean / lx**2)**(1/3)) * (1/sim.params.N)
+    print("F_h", F_h)
+
+    # Reynolds buoyancy 8
+    Rb8 = Re_8 * F_h**8
+    print("Rb8", Rb8)
+
+    if MAKE_TABLE:
+        to_print += ("{} & {:.4f} & {:.4e} & {:.4e} \\\\ \n".format(
+            gamma_table, F_h, Re_8, Rb8))
+
+if MAKE_TABLE:
+    with open(path_table, "w") as f:
+        to_print += ("\\end{tabular} \n"
+                    "\\end{table}")
+        f.write(to_print)
diff --git a/scripts/ns2d.strat/plot_reynolds_froude.py b/scripts/ns2d.strat/plot_reynolds_froude.py
new file mode 100644
--- /dev/null
+++ b/scripts/ns2d.strat/plot_reynolds_froude.py
@@ -0,0 +1,61 @@
+"""
+plot_reynolds_froude.py
+========================
+1/10/2018
+
+Makes plot buoyancy reynolds Vs Froude.
+"""
+
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+from glob import glob
+
+from compute_anisotropy import compute_anisotropy
+from compute_ratio_dissipation import compute_ratio_dissipation
+from compute_reynolds_froude import compute_buoyancy_reynolds
+from fluiddyn.output.rcparams import set_rcparams
+
+
+# Create path simulations
+path_root = "/fsnet/project/meige/2015/15DELDUCA/DataSim"
+directories = ["sim1920_no_shear_modes", "sim1920_modif_res_no_shear_modes"]
+paths_simulations = []
+for directory in directories:
+    paths_simulations += sorted(glob(os.path.join(path_root, directory, "NS2D*")))
+
+froudes = []
+reynoldsb = []
+anisotropies = []
+dissipations = []
+
+set_rcparams(fontsize=14, for_article=True)
+
+fig, ax = plt.subplots()
+ax.set_xlabel(r"$F_h$")
+ax.set_ylabel(r"$\mathcal{R}$")
+ax.set_xscale("log")
+ax.set_yscale("log")
+fig.text(0.8, 4e-7, r"$\frac{D(k_{fx})}{D(k_x)}$", fontsize=16)
+
+for path in paths_simulations:
+    F_h, Re_8, R_b = compute_buoyancy_reynolds(path)
+    anisotropy = compute_anisotropy(path)
+    dissipation = compute_ratio_dissipation(path)
+
+    froudes.append(F_h)
+    reynoldsb.append(R_b)
+    anisotropies.append(anisotropy)
+    dissipations.append(dissipation)
+
+    print("F_h", F_h)
+    print("Re_8", Re_8)
+    print("R_b", R_b)
+
+areas = 500 * np.asarray(anisotropies)**2
+scatter = ax.scatter(froudes, reynoldsb, s=areas, c=dissipations, alpha=0.7)
+ax.scatter(0.7, 1e-4, s=500 * np.asarray(0.5)**2, c="red")
+ax.text(0.64, 1e-5, "isotropy", fontsize=14, color="r")
+fig.colorbar(scatter)
+plt.show()