Skip to content
Snippets Groups Projects
save_figure9.py 9.22 KiB
Newer Older
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
import matplotlib.patches as patches
from util_simuls_regimes import get_sim

from math import degrees

from util import compute_kf_kb_ko_keta_kd, customize, save_fig

plt.rcParams["text.usetex"] = True


fig, axes = plt.subplots(
    ncols=2, nrows=2, figsize=(10, 2 * 3 * 4.5 / 4), constrained_layout=True
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]

# Standard Navier-Stokes
sim = get_sim("L", proj=False)

mean_values = sim.output.get_mean_values(tmin="t_last-2", customize=customize)
R2 = mean_values["R2"]
Uh2 = mean_values["Uh2"]
epsK = mean_values["epsK"]
Fh = mean_values["Fh"]
proj = sim.params.projection
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)

data = sim.output.spectra.load_kzkh_mean(
    tmin, key_to_load=["A", "Khd", "Kz", "Khr"]
)
kh = data["kh_spectra"]
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]

# Forcing
kf_min = sim.params.forcing.nkmin_forcing * delta_kz
kf_max = sim.params.forcing.nkmax_forcing * delta_kz
angle = sim.params.forcing.tcrandom_anisotropic.angle
delta_angle = sim.params.forcing.tcrandom_anisotropic.delta_angle

KH, KZ = np.meshgrid(kh, kz)
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
EKhr = data["Khr"]
Epolo = EKhd + EKz
Etoro = EKhr
E = Epolo + Etoro + EA
Ee = 2 * np.minimum(Epolo, EA)
Ed = (Epolo - EA) / (Epolo + EA)
E[E == 0] = 1e-15
levels = np.linspace(0, 1, 51, endpoint=True)
K = np.sqrt(KH**2 + KZ**2)
K[K == 0] = 1e-15

# Etoro / E (kh, kz)
cs = ax0.contourf(KH, KZ, Etoro / E, cmap=cm.binary, levels=levels)

data = sim.output.spect_energy_budg.load_mean(tmin=tmin)
kh = data["kh"]
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]
KH, KZ = np.meshgrid(kh, kz)
DA = data["diss_A"]
TA = data["transfer_A"]
TK = data["transfer_Kh"] + data["transfer_Kz"]
K2A = data["conv_K2A"]
DKh = data["diss_Kh"]
DKz = data["diss_Kz"]
DK = DKh + DKz
D = DA + DK
T = TA + TK
levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True)

# (Epolo - Epot)/ (Epolo + Epot) (kh, kz)
cs2 = ax2.contourf(KH, KZ, Ed, cmap=cm.seismic, levels=levels)


th = np.linspace(0, np.pi / 2, 50)
for ax in [ax0, ax2]:
    # kb
    ax.plot(kb * np.sin(th), kb * np.cos(th), color="cyan", linestyle="dashed")

    # Chi_L = 1/3, 3
    a = 1 / 3
    xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
    ax.plot(
        xa,
        xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
        linestyle="dotted",
        color="cyan",
    )
    a = 3
    xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
    ax.plot(
        xa,
        xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
        linestyle="dotted",
        color="cyan",
for ax in [ax0, ax2]:
    # Chi_nu = 1 (Kolmogorov scale)
    ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g")

    ax.set_xlim([kh[1], 0.8 * max(kh)])
    ax.set_ylim([kz[1], 0.8 * max(kh)])
    ax.set_xscale("log")
    ax.set_yscale("log")

    ax.set_xlabel(r"$k_h$", fontsize=20)
    ax.set_xticks(
        [1e1, 1e2, 1e3, kb, ko],
        [r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
        fontsize=14,
    )
    ax.set_ylabel(r"$k_z$", fontsize=20)
    ax.set_yticks(
        [1e1, 1e2, 1e3, kb, ko],
        [r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
        fontsize=14,
    )


# Without vortical modes
sim = get_sim("L", proj=True)

mean_values = sim.output.get_mean_values(tmin="t_last-2", customize=customize)
R2 = mean_values["R2"]
Uh2 = mean_values["Uh2"]
epsK = mean_values["epsK"]
Fh = mean_values["Fh"]
proj = sim.params.projection
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)

data = sim.output.spectra.load_kzkh_mean(
    tmin, key_to_load=["A", "Khd", "Kz", "Khr"]
)
kh = data["kh_spectra"]
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]

# Forcing
kf_min = sim.params.forcing.nkmin_forcing * delta_kz
kf_max = sim.params.forcing.nkmax_forcing * delta_kz
angle = sim.params.forcing.tcrandom_anisotropic.angle
delta_angle = sim.params.forcing.tcrandom_anisotropic.delta_angle

KH, KZ = np.meshgrid(kh, kz)
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
EKhr = data["Khr"]
Epolo = EKhd + EKz
Etoro = EKhr
E = Epolo + Etoro + EA
Ee = 2 * np.minimum(Epolo, EA)
Ed = (Epolo - EA) / (Epolo + EA)
E[E == 0] = 1e-15
levels = np.linspace(0, 1, 51, endpoint=True)
K = np.sqrt(KH**2 + KZ**2)
K[K == 0] = 1e-15

# Etoro / E (kh, kz)
cs = ax1.contourf(KH, KZ, Etoro / E, cmap=cm.binary, levels=levels)

data = sim.output.spect_energy_budg.load_mean(tmin=tmin)
kh = data["kh"]
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]
KH, KZ = np.meshgrid(kh, kz)
DA = data["diss_A"]
TA = data["transfer_A"]
TK = data["transfer_Kh"] + data["transfer_Kz"]
K2A = data["conv_K2A"]
DKh = data["diss_Kh"]
DKz = data["diss_Kz"]
DK = DKh + DKz
D = DA + DK
T = TA + TK
levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True)

# (Epolo - Epot)/ (Epolo + Epot) (kh, kz)
cs2 = ax3.contourf(KH, KZ, Ed, cmap=cm.seismic, levels=levels)


th = np.linspace(0, np.pi / 2, 50)
for ax in [ax1, ax3]:
    # kb
    ax.plot(kb * np.sin(th), kb * np.cos(th), color="cyan", linestyle="dashed")

    # Chi_L = 1/3, 3
    a = 1 / 3
    xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
    ax.plot(
        xa,
        xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
        linestyle="dotted",
        color="cyan",
    )
    a = 3
    xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
    ax.plot(
        xa,
        xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
        linestyle="dotted",
        color="cyan",
for ax in [ax1, ax3]:
    # Chi_nu = 1 (Kolmogorov scale)
    ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g")

    ax.set_xlim([kh[1], 0.8 * max(kh)])
    ax.set_ylim([kz[1], 0.8 * max(kh)])
    ax.set_xscale("log")
    ax.set_yscale("log")

    ax.set_xlabel(r"$k_h$", fontsize=20)
    ax.set_xticks(
        [1e1, 1e2, 1e3, kb, ko],
        [r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
        fontsize=14,
    )
    ax.set_ylabel(r"$k_z$", fontsize=20)
    ax.set_yticks(
        [1e1, 1e2, 1e3, kb, ko],
        [r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
        fontsize=14,
    )


for ax in [ax0, ax1]:
    ax.set_xlabel(None)
    ax.set_xticks([])
for ax in [ax1, ax3]:
    ax.set_ylabel(None)
    ax.set_yticks([])


ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", fontsize=20)
ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
ax2.set_title(r"$\rm (c)$", fontsize=20)
ax3.set_title(r"$\rm (c)$", fontsize=20)
for ax in [ax0, ax1, ax2, ax3]:
    # Forcing
    ax.add_patch(
        patches.Arc(
            xy=(0, 0),
            width=2 * kf_max,
            height=2 * kf_max,
            angle=0,
            theta1=90.0 - degrees(angle) - 0.5 * degrees(delta_angle),
            theta2=90.0 - degrees(angle) + 0.5 * degrees(delta_angle),
            linestyle="-",
            color="orange",
            linewidth=1,
        )
    )
    ax.add_patch(
        patches.Arc(
            xy=(0, 0),
            width=2 * kf_min,
            height=2 * kf_min,
            angle=0,
            theta1=90.0 - degrees(angle) - 0.5 * degrees(delta_angle),
            theta2=90.0 - degrees(angle) + 0.5 * degrees(delta_angle),
            linestyle="-",
            color="orange",
            linewidth=1,
        )
    )
    ax.plot(
        [
            kf_min * np.sin(angle - 0.5 * delta_angle),
            kf_max * np.sin(angle - 0.5 * delta_angle),
        ],
        [
            kf_min * np.cos(angle - 0.5 * delta_angle),
            kf_max * np.cos(angle - 0.5 * delta_angle),
        ],
        linestyle="-",
        color="orange",
        linewidth=1,
    )
    ax.plot(
        [
            kf_min * np.sin(angle + 0.5 * delta_angle),
            kf_max * np.sin(angle + 0.5 * delta_angle),
        ],
        [
            kf_min * np.cos(angle + 0.5 * delta_angle),
            kf_max * np.cos(angle + 0.5 * delta_angle),
        ],
        linestyle="-",
        color="orange",
        linewidth=1,
    )


# Add annotations for the lines and forcing region
ax1.text(
    4.5,
    1.5e1,
    r"Forcing",
    color="orange",
    rotation=0.7 * 180 / np.pi,
    fontsize=10,
)
ax1.text(1e1, 1.5e1, r"$\chi_{\mathbf{k}} < 1/3$", color="cyan", fontsize=14)
ax1.text(1e2, 1.5e2, r"$\chi_{\mathbf{k}} < 3$", color="cyan", fontsize=14)
# ax1.text(7e1, 5e1, r"$\gamma_{\mathbf{k}} < 1$", color="m", fontsize=14)
fig.tight_layout()
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.535, 0.02, 0.36])
cbar = fig.colorbar(cs, cax=cbar_ax)
cbar.set_ticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
cbar.set_ticklabels(
    [r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1.0$"], fontsize=14
)
cbar.set_label(r"$E_{\rm toro}/E$", fontsize=20)


# fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.11, 0.02, 0.36])
cbar = fig.colorbar(cs2, cax=cbar_ax)
cbar.set_ticks([-2 / 3, -1 / 3, 0.0, 1 / 3, 2 / 3])
cbar.set_ticklabels(
    [r"$-2/3$", r"$-1/3$", r"$0$", r"$1/3$", r"$2/3$"], fontsize=14
)
cbar.set_label(r"$\tilde{\mathcal{D}}$", fontsize=20)

save_fig(fig, f"figure9.png")

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