Newer
Older
import matplotlib.pyplot as plt
import numpy as np
from fluidsim import load
from fluidsim.util import load_params_simul, times_start_last_from_path
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
def plot_nlb(sim, ax):
t_start, t_last = times_start_last_from_path(path)
tmin = t_last - 2
mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
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
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
cmap=cm,
vmin=-0.5,
vmax=1.5,
shading="nearest",
)
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",
)
ax.plot([delta_kh, max(kh)], [delta_kh, 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([delta_kh, 2 * max(kh) / 3])
ax.set_ylim([delta_kh, 2 * max(kh) / 3])
# ax.set_xscale("lin")
# ax.set_yscale("lin")
return cs
Ns = [10, 20, 80]
nbax = 0
css = [None for i in range(6)]
fig, axes = plt.subplots(
ncols=2, nrows=3, figsize=(10, 3 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ax4 = axes[2, 0]
ax5 = axes[2, 1]
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]:
ax.set_xticks([])
for ax in [ax1, ax3, ax5]:
ax.set_yticks([])
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$(a)$", fontsize=16)
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)
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])
cbar.ax.set_ylabel(
r"$\log_{10}\left(\delta \omega_{\bm{k}}/ \omega_{\bm{k}}\right)$",
fontsize=16,
)