Skip to content
Snippets Groups Projects
save_nonlinear_ratio_one.py 4.58 KiB
Newer Older
import h5py
vlabarre's avatar
vlabarre committed
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm, ticker
vlabarre's avatar
vlabarre committed
from util_simuls_regimes import get_sim

from fluidsim import load
from fluidsim.util import load_params_simul, times_start_last_from_path

vlabarre's avatar
vlabarre committed
from util import (
vlabarre's avatar
vlabarre committed
    compute_kf_kb_ko_keta_kd,
    compute_omega_emp_vs_kzkh,
vlabarre's avatar
vlabarre committed
    customize,
    get_path_finer_resol,
    save_fig,
vlabarre's avatar
vlabarre committed
)
vlabarre's avatar
vlabarre committed

# Latex
vlabarre's avatar
vlabarre committed
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
vlabarre's avatar
vlabarre committed
cm = cm.get_cmap("inferno", 100)
vlabarre's avatar
vlabarre committed


def plot_nlb(sim, ax):
    t_start, t_last = times_start_last_from_path(path)
    tmin = t_last - 2
vlabarre's avatar
vlabarre committed
    kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
vlabarre's avatar
vlabarre committed
    mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
vlabarre's avatar
vlabarre committed
    Fh = mean_values["Fh"]
    epsK = mean_values["epsK"]
    params = load_params_simul(path)
    nh = nx = params.oper.nx
    proj = params.projection
    N = sim.params.N
    path_spec = sorted(path.glob(f"spatiotemporal/periodogram_[0-9]*.h5"))
    assert len(path_spec) == 1, f"Not only 1 periodogram in {path} \n"
    path_spec = path_spec[0]
    with h5py.File(path_spec, "r") as f:
        kh = f["kh_spectra"][:]
        kz = f["kz_spectra"][:]
        delta_kh = kh[1]
        delta_kz = kz[1]
        omegas = f["omegas"][:]
        EA = f["spectrum_A"][:]
        EKz = f["spectrum_K"][:] - f["spectrum_Khd"][:] - f["spectrum_Khr"][:]
        Epolo = f["spectrum_Khd"][:] + EKz
        Etoro = f["spectrum_Khr"][:]
        E = Epolo + Etoro + EA
vlabarre's avatar
vlabarre committed
        Ee = 2 * np.minimum(EA, Epolo)
vlabarre's avatar
vlabarre committed
        Ed = EA + Epolo - Ee
        spectrum = Epolo + EA
        omega_emp, delta_omega_emp = compute_omega_emp_vs_kzkh(
            N, spectrum, kh, kz, omegas
        )
        KH, KZ = np.meshgrid(kh, kz)
        K = (KH**2 + KZ**2) ** 0.5
        K_NOZERO = K.copy()
        K_NOZERO[K_NOZERO == 0] = 1e-16
        omega_disp = N * KH / K_NOZERO
vlabarre's avatar
vlabarre committed
        chi = (K**2 * epsK) ** (1 / 3) / omega_disp
vlabarre's avatar
vlabarre committed

        cs = ax.pcolormesh(
            kh,
            kz,
vlabarre's avatar
vlabarre committed
            np.log10(delta_omega_emp / omega_disp),
vlabarre's avatar
vlabarre committed
            cmap=cm,
            vmin=-0.5,
            vmax=1.5,
            shading="nearest",
        )

vlabarre's avatar
vlabarre committed
        th = np.linspace(0, np.pi / 2, 100, endpoint=True)
vlabarre's avatar
vlabarre committed
        ax.plot(kb * np.sin(th), kb * np.cos(th), color="k", linestyle="dotted")
        ax.plot(ko * np.sin(th), ko * np.cos(th), "k--")
        a = 3
        xa = np.linspace(delta_kh, a**1.5 * ko, 50, endpoint=True)
        ax.plot(
            xa,
            xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
            linestyle="dashed",
            color="gray",
        )
        a = 1 / 3
        xa = np.linspace(delta_kh, a**1.5 * ko, 50, endpoint=True)
        ax.plot(
            xa,
            xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
            linestyle="dotted",
            color="gray",
        )
vlabarre's avatar
vlabarre committed

        ax.plot([delta_kh, max(kh)], [delta_kh, max(kh)], "k-")
vlabarre's avatar
vlabarre committed
        ax.plot(kf * np.sin(th), kf * np.cos(th), linestyle="--", color="orange")
        ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="--", color="g")
vlabarre's avatar
vlabarre committed
        ax.set_xlim([delta_kh, 2 * max(kh) / 3])
        ax.set_ylim([delta_kh, 2 * max(kh) / 3])
        # ax.set_xscale("lin")
        # ax.set_yscale("lin")
vlabarre's avatar
vlabarre committed
    return cs


Ns = [10, 20, 80]
nbax = 0
css = [None for i in range(6)]

vlabarre's avatar
vlabarre committed
fig, axes = plt.subplots(
    ncols=2, nrows=3, figsize=(10, 3 * 3 * 4.5 / 4), constrained_layout=True
)
vlabarre's avatar
vlabarre committed
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ax4 = axes[2, 0]
ax5 = axes[2, 1]
vlabarre's avatar
vlabarre committed
axs = [ax0, ax1, ax2, ax3, ax4, ax5]
vlabarre's avatar
vlabarre committed

for N in Ns:
    for proj in [False, True]:
        path = get_path_finer_resol(N=N, Rb=None, proj=proj, ratio_one=True)
        sim = load(path)
        css[nbax] = plot_nlb(sim, axs[nbax])
        nbax += 1


for ax in [ax0, ax2, ax4]:
vlabarre's avatar
vlabarre committed
    ax.set_ylabel(r"$k_z$", fontsize=16)
vlabarre's avatar
vlabarre committed
for ax in [ax4, ax5]:
vlabarre's avatar
vlabarre committed
    ax.set_xlabel(r"$k_h$", fontsize=16)
vlabarre's avatar
vlabarre committed
for ax in [ax0, ax1, ax2, ax3]:
vlabarre's avatar
vlabarre committed
    ax.set_xticks([])
for ax in [ax1, ax3, ax5]:
    ax.set_yticks([])

vlabarre's avatar
vlabarre committed
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$(a)$", fontsize=16)
vlabarre's avatar
vlabarre committed
ax1.set_title(r"Without vortical modes" + "\n" + r"$(b)$", fontsize=16)
ax2.set_title(r"$(c)$", fontsize=16)
ax3.set_title(r"$(d)$", fontsize=16)
ax4.set_title(r"$(e)$", fontsize=16)
ax5.set_title(r"$(f)$", fontsize=16)
vlabarre's avatar
vlabarre committed
fig.tight_layout()
vlabarre's avatar
vlabarre committed
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.28, 0.02, 0.35])
cbar = fig.colorbar(css[3], cax=cbar_ax)
cbar.set_ticks([-0.5, 0.0, 0.5, 1.0, 1.5])
vlabarre's avatar
vlabarre committed
cbar.ax.set_ylabel(
    r"$\log_{10}\left(\delta \omega_{\bm{k}}/ \omega_{\bm{k}}\right)$",
    fontsize=16,
)
vlabarre's avatar
vlabarre committed


save_fig(fig, f"fig_nonlinear_ratio_one.png")

if __name__ == "__main__":
    plt.show()