Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • fluiddyn/fluiddyn_papers
1 result
Show changes
Showing
with 2464 additions and 0 deletions
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm, ticker
from util_simuls_regimes import get_sim
from fluidsim import load
from util import (
compute_kf_kb_ko_keta_kd,
customize,
get_path_finer_resol,
save_fig,
)
plt.rcParams["text.usetex"] = True
###
N = 20
Rb = None
proj = False
ratio_one = True
###
path = get_path_finer_resol(N, Rb, proj, ratio_one)
sim = load(path)
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]
mean_values = sim.output.get_mean_values(tmin="t_last-2", customize=customize)
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]
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 - Ee
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
# EA / E (kh, kz)
cs = ax0.contourf(KH, KZ, Ee / E, cmap=cm.binary, levels=levels)
ax0.set_title(r"$E_{equi}/E$", fontsize=16)
# Epolo / E (kh, kz)
cs = ax1.contourf(KH, KZ, Ed / E, cmap=cm.binary, levels=levels)
ax1.set_title(r"$E_{diff}/E$", fontsize=16)
# Etoro / E (kh, kz)
cs = ax2.contourf(KH, KZ, Etoro / E, cmap=cm.binary, levels=levels)
ax2.set_title(r"$E_{toro}/E$", fontsize=16)
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(-1, 1, 51, endpoint=True)
cs2 = ax3.contourf(
KH,
KZ,
K2A / (D + np.abs(TA) + np.abs(TK) + np.abs(K2A)),
cmap=cm.seismic,
levels=levels,
)
ax3.set_title(r"$\tilde{\mathcal{B}}$", fontsize=16)
# ax0.legend()
th = np.linspace(0, np.pi / 2, 50)
for ax in [ax0, ax1, ax2, ax3]:
ax.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-")
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(kh[1], 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(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="gray",
)
ax.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-")
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")
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")
for ax in [ax0, ax2]:
ax.set_ylabel(r"$k_z$", fontsize=16)
ax.set_yticks(
[1e1, 1e2, 1e3, kb, ko],
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_b$", r"$k_O$"],
)
for ax in [ax2, ax3]:
ax.set_xlabel(r"$k_h$", fontsize=16)
ax.set_xticks(
[1e1, 1e2, 1e3, kb, ko],
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_b$", r"$k_O$"],
)
for ax in [ax0, ax1]:
ax.set_xticks([])
for ax in [ax1, ax3]:
ax.set_yticks([])
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.53, 0.02, 0.35])
cbar = fig.colorbar(cs, cax=cbar_ax)
cbar.set_ticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
# fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.11, 0.02, 0.35])
cbar = fig.colorbar(cs2, cax=cbar_ax)
cbar.set_ticks([-1.0, -0.5, 0.0, 0.5, 1.0])
# fig.tight_layout()
save_fig(
fig,
f"fig_seb_regimes_one_couple_proj{proj}_ratio_one{ratio_one}_N{N}_Rb{Rb}.png",
)
if __name__ == "__main__":
plt.show()
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.