Skip to content
Snippets Groups Projects
save_spectra_1R.py 4.08 KiB
Newer Older
vlabarre's avatar
vlabarre committed
import matplotlib.cm
Pierre Augier's avatar
Pierre Augier committed
import matplotlib.pyplot as plt
vlabarre's avatar
vlabarre committed

plt.rcParams["text.usetex"] = True
Pierre Augier's avatar
Pierre Augier committed
from fluidsim import load

vlabarre's avatar
vlabarre committed
import numpy as np

from util import (
    save_fig,
    couples320,
    get_path_finer_resol,
    customize,
    add_letters,
)
vlabarre's avatar
vlabarre committed
cm = matplotlib.cm.get_cmap("inferno", 100)

Pierre Augier's avatar
Pierre Augier committed
R_i = 20

Ns = sorted(N for (N, R_i_loc) in couples320 if R_i_loc == R_i)

vlabarre's avatar
vlabarre committed
Ns = [N for N in Ns if N < 60]

Pierre Augier's avatar
Pierre Augier committed
paths = [get_path_finer_resol(N, R_i) for N in Ns]
paths = [path for path in paths if path is not None]

vlabarre's avatar
vlabarre committed
fig, axes = plt.subplots(ncols=2, nrows=3, figsize=(10, 3 * 3 * 4.5 / 4))
vlabarre's avatar
vlabarre committed

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]
Pierre Augier's avatar
Pierre Augier committed


for path, N in zip(paths, Ns):
    sim = load(path, hide_stdout=True)
vlabarre's avatar
vlabarre committed
    mean_values = sim.output.get_mean_values(tmin="t_last-2", customize=customize)
Pierre Augier's avatar
Pierre Augier committed
    Uh2 = mean_values["Uh2"]
vlabarre's avatar
vlabarre committed
    kb = N / Uh2**0.5
Pierre Augier's avatar
Pierre Augier committed
    epsK = mean_values["epsK"]
vlabarre's avatar
vlabarre committed
    ko = (N**3 / epsK) ** 0.5
    Fh = mean_values["Fh"]
Pierre Augier's avatar
Pierre Augier committed

    t_start, t_last = sim.output.print_stdout.get_times_start_last()
    tmin = max(10, min(t_start + 2, t_last - 1))
Pierre Augier's avatar
Pierre Augier committed

vlabarre's avatar
vlabarre committed
    data = sim.output.spectra.load_kzkh_mean(
        tmin, key_to_load=["A", "Khd", "Kz", "Khr"]
    )
vlabarre's avatar
vlabarre committed
    kh = data["kh_spectra"]
    kz = data["kz"]
    delta_kh = kh[1]
    delta_kz = kz[1]
Pierre Augier's avatar
Pierre Augier committed
    KH, KZ = np.meshgrid(kh, kz)
Pierre Augier's avatar
Pierre Augier committed
    EA = data["A"]
    EKhd = data["Khd"]
    EKz = data["Kz"]
    EKhr = data["Khr"]

vlabarre's avatar
vlabarre committed
    Epolo = EKhd + EKz
    Etoro = EKhr
Pierre Augier's avatar
Pierre Augier committed
    EKh = EKhr + EKhd

vlabarre's avatar
vlabarre committed
    Epolo_vs_kh = np.sum(Epolo, axis=0) * delta_kz
    Epolo_vs_kz = np.sum(Epolo, axis=1) * delta_kh

    Etoro_vs_kh = np.sum(Etoro, axis=0) * delta_kz
    Etoro_vs_kz = np.sum(Etoro, axis=1) * delta_kh

    EA_vs_kh = np.sum(EA, axis=0) * delta_kz
    EA_vs_kz = np.sum(EA, axis=1) * delta_kh
Pierre Augier's avatar
Pierre Augier committed

vlabarre's avatar
vlabarre committed
    color = cm((np.log10(Fh) - (-2)) / (0 - (-2)))
vlabarre's avatar
vlabarre committed
    # Potential
vlabarre's avatar
vlabarre committed
    (line,) = ax0.plot(
        kh / kb,
        EA_vs_kh[:] * kh[:] ** (5 / 3),
        label=rf"$F_h={Fh:.2f}$",
        color=color,
    )
    ax0.axvline(ko / kb, linestyle="dashed", color=color)
vlabarre's avatar
vlabarre committed
    ax1.plot(kz / kb, EA_vs_kz[:] * kz[:] ** (5 / 3), color=color)
vlabarre's avatar
vlabarre committed
    # ax1.plot([1, 10], [3*1, 3*10.**(5 / 3 - 2.3)], "k-")
    ax1.axvline(ko / kb, linestyle="dashed", color=color)
vlabarre's avatar
vlabarre committed

    # Poloidal
vlabarre's avatar
vlabarre committed
    ax2.plot(kh / kb, Epolo_vs_kh[:] * kh[:] ** (5 / 3), color=color)
vlabarre's avatar
vlabarre committed
    ax2.axvline(ko / kb, linestyle="dashed", color=color)
vlabarre's avatar
vlabarre committed
    ax3.plot(kz / kb, Epolo_vs_kz[:] * kz[:] ** (5 / 3), color=color)
vlabarre's avatar
vlabarre committed
    # ax1.plot([1, 10], [3*1, 3*10.**(5 / 3 - 2.3)], "k-")
    ax3.axvline(ko / kb, linestyle="dashed", color=color)
vlabarre's avatar
vlabarre committed
    # Toroidal
vlabarre's avatar
vlabarre committed
    ax4.plot(kh / kb, Etoro_vs_kh[:] * kh[:] ** (5 / 3), color=color)
vlabarre's avatar
vlabarre committed
    ax4.axvline(ko / kb, linestyle="dashed", color=color)
vlabarre's avatar
vlabarre committed
    ax5.plot(kz / kb, Etoro_vs_kz[:] * kz[:] ** (5 / 3), color=color)
vlabarre's avatar
vlabarre committed
    # ax1.plot([1, 10], [3*1, 3*10.**(5 / 3 - 2.3)], "k-")
    ax5.axvline(ko / kb, linestyle="dashed", color=color)
vlabarre's avatar
vlabarre committed

vlabarre's avatar
vlabarre committed
ax0.set_xlabel(r"$k_h / k_{\rm b}$", fontsize=14)
ax1.set_xlabel(r"$k_z / k_{\rm b}$", fontsize=14)
ax2.set_xlabel(r"$k_h / k_{\rm b}$", fontsize=14)
ax3.set_xlabel(r"$k_z / k_{\rm b}$", fontsize=14)
ax4.set_xlabel(r"$k_h / k_{\rm b}$", fontsize=14)
ax5.set_xlabel(r"$k_z / k_{\rm b}$", fontsize=14)
ax0.set_ylabel(r"$E_{\rm A}(k_h) {k_h}^{5/3}$", fontsize=14)
ax1.set_ylabel(r"$E_{\rm A}(k_z) {k_z}^{5/3}$", fontsize=14)
ax2.set_ylabel(r"$E_{\rm polo}(k_h) {k_h}^{5/3}$", fontsize=14)
ax3.set_ylabel(r"$E_{\rm polo}(k_z) {k_z}^{5/3}$", fontsize=14)
ax4.set_ylabel(r"$E_{\rm toro}(k_h) {k_h}^{5/3}$", fontsize=14)
ax5.set_ylabel(r"$E_{\rm toro}(k_z) {k_z}^{5/3}$", fontsize=14)
vlabarre's avatar
vlabarre committed
ax0.legend(fontsize=10)
vlabarre's avatar
vlabarre committed
for ax in [ax0, ax1, ax2, ax3, ax4, ax5]:
Pierre Augier's avatar
Pierre Augier committed
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_ylim(bottom=1e-2, top=4)

vlabarre's avatar
vlabarre committed
for ax in [ax0, ax1, ax2, ax3]:
    ax.set_xlabel(None)
    ax.set_xticks([])

for ax in [ax1, ax3, ax5]:
    kk = np.array([1, 10])
vlabarre's avatar
vlabarre committed
    ax.plot(kk, 0.5 * kk ** (-3 + 5 / 3), "k-")
vlabarre's avatar
vlabarre committed
    ax.text(1e1, 1, r"$\propto k_z^{-2}$", fontsize=12)
vlabarre's avatar
vlabarre committed
    ax.plot(kk, 3 * kk ** (-2 + 5 / 3), "k-")
vlabarre's avatar
vlabarre committed
    ax.text(0.8, 2e-1, r"$\propto k_z^{-3}$", fontsize=12)


Pierre Augier's avatar
Pierre Augier committed
add_letters(fig, "ab")
Pierre Augier's avatar
Pierre Augier committed
fig.tight_layout()
save_fig(fig, "fig_spectra_1R.png")

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