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 1990 additions and 2015 deletions
import sys
import h5py
import matplotlib.cm
import matplotlib.pyplot as plt
......@@ -2,4 +2,3 @@
import matplotlib.cm
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
......@@ -5,5 +4,4 @@
import numpy as np
from math import degrees
from util_simuls_regimes import get_sim
from util import (
......@@ -8,4 +6,3 @@
from util import (
compute_kf_kb_ko_keta_kd,
save_fig,
......@@ -11,6 +8,4 @@
save_fig,
compute_E_waves_vs_kh_kz,
spectra_vs_khomega_slice,
spectra_vs_kzomega_slice,
customize,
)
......@@ -15,7 +10,3 @@
)
from util_simuls_regimes import get_sim
from util_dataframe import df, df_proj, df_ratio_one, df_proj_ratio_one
cm = matplotlib.cm.get_cmap("inferno", 100)
......@@ -21,8 +12,7 @@
cm = matplotlib.cm.get_cmap("inferno", 100)
cmbin = matplotlib.cm.get_cmap("binary", 100)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
......@@ -23,85 +13,10 @@
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
print(sys.argv)
letter = sys.argv[-1]
if letter not in "DLOWP":
letter = "L"
normalize = False
# def plot_spectra(sim, ax, key="Ee"):
def plot_slices_omega_ki_waves(sim, ax, key="Ee", key_k="kh", ik=10, delta=1e-1):
N = sim.params.N
E_waves_omega, E_waves, E_tot, kh, kz, omegas = compute_E_waves_vs_kh_kz(
sim, delta=delta
)
if key_k == "kh":
spectra_normalized = spectra_vs_khomega_slice(
E_waves_omega, kh, omegas, ik, normalize=normalize
)
ax.set_xlabel(r"$k_h$", fontsize=20)
ax.set_xlim([kh[1], max(kh)])
cs = ax.pcolormesh(
kh,
omegas / N,
np.log10(spectra_normalized.transpose()),
cmap=cm,
vmin=-9,
vmax=-4,
shading="nearest",
)
k = (kh**2 + kz[ik] ** 2) ** 0.5
omega_disp = kh / k
ax.plot(kh, omega_disp, "k-")
ax.plot(kh, omega_disp + 0.5 * delta * omega_disp, "b-")
ax.plot(kh, omega_disp - 0.5 * delta * omega_disp, "b-")
ax.text(10, 2.5, rf"$k_z={kz[ik]:.1f}$", color="w", fontsize=14)
elif key_k == "kz":
spectra_normalized = spectra_vs_kzomega_slice(
E_waves_omega, kz, omegas, ik, normalize=normalize
)
ax.set_xlabel(r"$k_z$", fontsize=20)
ax.set_xlim([kz[1], max(kz)])
cs = ax.pcolormesh(
kz,
omegas / N,
np.log10(spectra_normalized.transpose()),
cmap=cm,
vmin=-9,
vmax=-4,
shading="nearest",
)
k = (kh[ik] ** 2 + kz**2) ** 0.5
omega_disp = kh[ik] / k
ax.plot(kz, omega_disp, "k-")
ax.plot(kz, omega_disp + 0.5 * delta * omega_disp, "b-")
ax.plot(kz, omega_disp - 0.5 * delta * omega_disp, "b-")
ax.text(15, 2.5, rf"$k_h={kh[ik]:.1f}$", color="w", fontsize=14)
ax.set_ylabel(r"$\omega / N$", fontsize=20)
ax.set_ylim([0, 3])
ax.set_xticks([25, 50, 75, 100, 125])
ax.set_xticklabels(
[r"$25$", r"$50$", r"$75$", r"$100$", r"$125$"], fontsize=14
)
ax.set_yticks([0, 1, 2, 3])
ax.set_yticklabels([r"$0$", r"$1$", r"$2$", r"$3$"], fontsize=14)
cs.cmap.set_bad("k")
cs.cmap.set_under("k")
return cs
def plot_ratio_E_waves_vs_kh_kz(sim, ax):
def plot_spectra_vs_omegak(sim, ax, key="Ee"):
N = sim.params.N
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
......@@ -105,7 +20,6 @@
N = sim.params.N
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)
E_waves_omega, E_waves, E_tot, kh, kz, omegas = compute_E_waves_vs_kh_kz(sim)
ratio_waves = E_waves / E_tot
mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
Uh = np.sqrt(mean_values["Uh2"])
......@@ -111,44 +25,33 @@
ax.set_xlabel(r"$k_h$", fontsize=20)
ax.set_xlim([kh[1], max(kh)])
ax.set_xticks([25, 50, 75, 100, 125])
ax.set_xticklabels(
[r"$25$", r"$50$", r"$75$", r"$100$", r"$125$"], fontsize=14
)
ax.set_ylabel(r"$k_z$", fontsize=20)
ax.set_ylim([kz[1], max(kz)])
ax.set_yticks([25, 50, 75, 100, 125])
ax.set_yticklabels(
[r"$25$", r"$50$", r"$75$", r"$100$", r"$125$"], fontsize=14
)
cs = ax.pcolormesh(
kh,
kz,
ratio_waves,
cmap=cmbin,
vmin=0,
vmax=1,
shading="nearest",
)
f = kh / np.tan(0.3)
ax.plot(kh, f, linestyle="-", linewidth=1, color="orange")
f = kh / np.tan(0.15)
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange")
f = kh / np.tan(0.075)
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange")
f = kh / np.tan(0.225)
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange")
# kb
th = np.linspace(0, np.pi / 2, 50)
ax.plot(kb * np.sin(th), kb * np.cos(th), linestyle="dashed", color="c")
# 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="c",
)
path = sim.params.path_run
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"][:]
omegas = f["omegas"][:]
Epot = f["spectrum_A"][:]
EKz = f["spectrum_K"][:] - f["spectrum_Khd"][:] - f["spectrum_Khr"][:]
Epolo = f["spectrum_Khd"][:] + EKz
Etoro = f["spectrum_Khr"][:]
delta_kh = kh[1]
delta_kz = kz[1]
if key == "Epot":
spectrum = Epot
elif key == "EK":
spectrum = Epolo + Etoro
elif key == "Epolo":
spectrum = Epolo
elif key == "Ee":
spectrum = 2 * np.minimum(Epot, Epolo)
elif key == "Ed":
spectrum = Epot + Epolo - 2 * np.minimum(Epot, Epolo)
elif key == "Etoro":
spectrum = Etoro
elif key == "Etot":
spectrum = Epolo + Etoro + Epot
else:
print(f"Don't know key: {key} \n")
exit
......@@ -154,20 +57,34 @@
delta_kh = kh[1]
delta_kz = kz[1]
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
# 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,
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
omega_disp_grid = np.linspace(0, N, 25, endpoint=True)
OMEGA_DISP_GRID, OMEGA_GRID = np.meshgrid(omega_disp_grid, omegas)
SPECTRUM = np.zeros(OMEGA_DISP_GRID.shape)
for iod in range(len(omega_disp_grid) - 1):
for io in range(len(omegas)):
omega_disp_min = omega_disp_grid[iod]
omega_disp_max = omega_disp_grid[iod + 1]
IOD = np.where(
np.logical_and(
omega_disp >= omega_disp_min, omega_disp < omega_disp_max
)
)
for i in range(len(IOD[0])):
SPECTRUM[io, iod] += spectrum[IOD[0][i], IOD[1][i], io]
cs = ax.pcolormesh(
OMEGA_DISP_GRID / N,
OMEGA_GRID / N,
np.log10(SPECTRUM),
cmap=cm,
vmin=-6,
vmax=-3,
shading="nearest",
)
......@@ -173,2 +90,8 @@
)
ax.set_xlabel(r"$\sin \theta_{\boldsymbol{k}}$", fontsize=20)
ax.set_xlim([0, 1])
ax.set_xticks([0, 0.25, 0.5, 0.75, 1])
ax.set_xticklabels(
[r"$0$", r"$0.25$", r"$0.5$", r"$0.75$", r"$1$"], fontsize=14
)
......@@ -174,44 +97,8 @@
)
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,
)
ax.set_ylabel(r"$\omega / N$", fontsize=20)
ax.set_ylim([0, 3])
ax.set_yticks([0, 1, 2, 3])
ax.set_yticklabels([r"$0$", r"$1$", r"$2$", r"$3$"], fontsize=14)
return cs
......@@ -215,4 +102,6 @@
return cs
# assert sim.params.oper.nx == sim_proj.params.oper.nx, f"Not the same resolution for simulation Without vortical modes: {sim.params.oper.nx} vs {sim_proj.params.oper.nx}"
fig, axes = plt.subplots(
......@@ -218,4 +107,4 @@
fig, axes = plt.subplots(
ncols=2, nrows=1, figsize=(10, 1.2 * 3 * 4.5 / 4), constrained_layout=True
ncols=2, nrows=3, figsize=(10, 3 * 3 * 4.5 / 4), constrained_layout=True
)
......@@ -220,6 +109,13 @@
)
ax0 = axes[0]
ax1 = axes[1]
css = [None for i in range(6)]
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ax4 = axes[2, 0]
ax5 = axes[2, 1]
......@@ -224,7 +120,14 @@
delta = 0.1
sim = get_sim(letter, proj=False)
# cs0 = plot_slices_omega_ki_waves(sim, ax0, key="Ee", key_k="kh", ik=3, delta=delta)
cs0 = plot_ratio_E_waves_vs_kh_kz(sim, ax0)
sim = get_sim("S1", proj=False)
css[0] = plot_spectra_vs_omegak(sim, ax0, key="Etot")
sim_proj = get_sim("S1", proj=True)
css[1] = plot_spectra_vs_omegak(sim_proj, ax1, key="Etot")
sim = get_sim("S2", proj=False)
css[2] = plot_spectra_vs_omegak(sim, ax2, key="Etot")
sim_proj = get_sim("S2", proj=True)
css[3] = plot_spectra_vs_omegak(sim_proj, ax3, key="Etot")
......@@ -230,5 +133,7 @@
sim = get_sim(letter, proj=True)
# cs1 = plot_slices_omega_ki_waves(sim, ax1, key="Ee", key_k="kh", ik=3, delta=delta)
cs1 = plot_ratio_E_waves_vs_kh_kz(sim, ax1)
sim = get_sim("S3", proj=False)
css[4] = plot_spectra_vs_omegak(sim, ax4, key="Etot")
sim_proj = get_sim("S3", proj=True)
css[5] = plot_spectra_vs_omegak(sim_proj, ax5, key="Etot")
......@@ -234,6 +139,12 @@
ax1.set_yticklabels([])
ax1.set_ylabel("")
for ax in [ax1, ax3, ax5]:
ax.set_ylabel("")
ax.set_yticklabels([])
for ax in [ax0, ax1, ax2, ax3]:
ax.set_xlabel("")
ax.set_xticklabels([])
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)
......@@ -237,4 +148,8 @@
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 (d)$", fontsize=20)
ax4.set_title(r"$\rm (e)$", fontsize=20)
ax5.set_title(r"$\rm (f)$", fontsize=20)
......@@ -240,3 +155,5 @@
for ax in [ax0, ax1, ax2, ax3, ax4, ax5]:
ax.plot([0, 1], [0, 1], "w--")
fig.tight_layout()
......@@ -241,5 +158,9 @@
fig.tight_layout()
fig.subplots_adjust(right=0.85)
# fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.19, 0.02, 0.63])
cbar = fig.colorbar(css[0], cax=cbar_ax)
cbar.set_ticks([-3, -4, -5, -6])
cbar.set_ticklabels([r"$-3$", r"$-4$", r"$-5$", r"$-6$"], fontsize=16)
......@@ -244,11 +165,7 @@
cbar_ax2 = fig.add_axes([0.88, 0.17, 0.02, 0.65])
cbar2 = fig.colorbar(cs1, cax=cbar_ax2, cmap=cmbin, orientation="vertical")
cbar2.set_label(r"$\tilde{E}_{\rm wave}(k_h,k_z)$", fontsize=20)
cbar2.set_ticks([0, 0.2, 0.4, 0.6, 0.8, 1])
cbar2.set_ticklabels(
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1$"], fontsize=14
cbar.ax.set_ylabel(
r"$\log_{10} E(\sin \theta_{\boldsymbol{k}},\omega)$", fontsize=20
)
......@@ -252,7 +169,9 @@
)
save_fig(fig, "figure14.png")
fig.subplots_adjust(right=0.85, wspace=0.05, hspace=0.15)
save_fig(fig, f"figure14.png")
if __name__ == "__main__":
plt.show()
import sys
import matplotlib.cm
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
......@@ -4,4 +5,4 @@
import numpy as np
from util import save_fig
from math import degrees
......@@ -7,4 +8,12 @@
from util_dataframe import df, df_proj
from util import (
compute_kf_kb_ko_keta_kd,
save_fig,
compute_E_waves_vs_kh_kz,
spectra_vs_khomega_slice,
spectra_vs_kzomega_slice,
)
from util_simuls_regimes import get_sim
cm = matplotlib.cm.get_cmap("inferno", 100)
......@@ -9,8 +18,9 @@
cm = matplotlib.cm.get_cmap("inferno", 100)
cmbin = matplotlib.cm.get_cmap("binary", 100)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
......@@ -11,19 +21,23 @@
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
def plot_ratio_waves_vs_Fh(ax, data=df, datab=df_proj):
ax.scatter(
data["Fh"],
data["R_waves"],
c=np.log10(data["R2"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
label=r"Standard Navier-Stokes",
vmin=-1,
vmax=4,
print(sys.argv)
letter = sys.argv[-1]
if letter not in "DLOWP":
letter = "L"
normalize = False
# def plot_spectra(sim, ax, key="Ee"):
def plot_slices_omega_ki_waves(sim, ax, key="Ee", key_k="kh", ik=10, delta=1e-1):
N = sim.params.N
E_waves_omega, E_waves, E_tot, kh, kz, omegas = compute_E_waves_vs_kh_kz(
sim, delta=delta
)
......@@ -29,13 +43,52 @@
)
cs = ax.scatter(
datab["Fh"],
datab["R_waves"],
c=np.log10(datab["R2"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
label=r"Without vortical modes",
vmin=-1,
vmax=4,
if key_k == "kh":
spectra_normalized = spectra_vs_khomega_slice(
E_waves_omega, kh, omegas, ik, normalize=normalize
)
ax.set_xlabel(r"$k_h$", fontsize=20)
ax.set_xlim([kh[1], max(kh)])
cs = ax.pcolormesh(
kh,
omegas / N,
np.log10(spectra_normalized.transpose()),
cmap=cm,
vmin=-9,
vmax=-4,
shading="nearest",
)
k = (kh**2 + kz[ik] ** 2) ** 0.5
omega_disp = kh / k
ax.plot(kh, omega_disp, "k-")
ax.plot(kh, omega_disp + 0.5 * delta * omega_disp, "b-")
ax.plot(kh, omega_disp - 0.5 * delta * omega_disp, "b-")
ax.text(10, 2.5, rf"$k_z={kz[ik]:.1f}$", color="w", fontsize=14)
elif key_k == "kz":
spectra_normalized = spectra_vs_kzomega_slice(
E_waves_omega, kz, omegas, ik, normalize=normalize
)
ax.set_xlabel(r"$k_z$", fontsize=20)
ax.set_xlim([kz[1], max(kz)])
cs = ax.pcolormesh(
kz,
omegas / N,
np.log10(spectra_normalized.transpose()),
cmap=cm,
vmin=-9,
vmax=-4,
shading="nearest",
)
k = (kh[ik] ** 2 + kz**2) ** 0.5
omega_disp = kh[ik] / k
ax.plot(kz, omega_disp, "k-")
ax.plot(kz, omega_disp + 0.5 * delta * omega_disp, "b-")
ax.plot(kz, omega_disp - 0.5 * delta * omega_disp, "b-")
ax.text(15, 2.5, rf"$k_h={kh[ik]:.1f}$", color="w", fontsize=14)
ax.set_ylabel(r"$\omega / N$", fontsize=20)
ax.set_ylim([0, 3])
ax.set_xticks([25, 50, 75, 100, 125])
ax.set_xticklabels(
[r"$25$", r"$50$", r"$75$", r"$100$", r"$125$"], fontsize=14
)
......@@ -41,16 +94,8 @@
)
ax.set_xlim([1e-3, 20])
ax.set_xticks([1e-3, 1e-2, 1e-1, 1e0, 1e1])
ax.set_xticklabels([1e-3, 1e-2, 1e-1, 1e0, 1e1], fontsize=14)
ax.set_xscale("log")
ax.set_ylim([0, 1])
ax.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax.set_yticklabels(
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1.0$"], fontsize=14
)
ax.set_xlabel(r"$F_h$", fontsize=20)
ax.set_ylabel(r"$\tilde{E}_{\rm wave}$", fontsize=20)
# ax.grid(True)
ax.set_yticks([0, 1, 2, 3])
ax.set_yticklabels([r"$0$", r"$1$", r"$2$", r"$3$"], fontsize=14)
cs.cmap.set_bad("k")
cs.cmap.set_under("k")
return cs
......@@ -54,15 +99,30 @@
return cs
def plot_ratio_waves_vs_R(ax, data=df, datab=df_proj):
ax.scatter(
data["R2"],
data["R_waves"],
c=np.log10(data["Fh"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
label=r"Standard Navier-Stokes",
vmin=-3,
def plot_ratio_E_waves_vs_kh_kz(sim, ax):
N = sim.params.N
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)
E_waves_omega, E_waves, E_tot, kh, kz, omegas = compute_E_waves_vs_kh_kz(sim)
ratio_waves = E_waves / E_tot
ax.set_xlabel(r"$k_h$", fontsize=20)
ax.set_xlim([kh[1], max(kh)])
ax.set_xticks([25, 50, 75, 100, 125])
ax.set_xticklabels(
[r"$25$", r"$50$", r"$75$", r"$100$", r"$125$"], fontsize=14
)
ax.set_ylabel(r"$k_z$", fontsize=20)
ax.set_ylim([kz[1], max(kz)])
ax.set_yticks([25, 50, 75, 100, 125])
ax.set_yticklabels(
[r"$25$", r"$50$", r"$75$", r"$100$", r"$125$"], fontsize=14
)
cs = ax.pcolormesh(
kh,
kz,
ratio_waves,
cmap=cmbin,
vmin=0,
vmax=1,
......@@ -68,2 +128,24 @@
vmax=1,
shading="nearest",
)
f = kh / np.tan(0.3)
ax.plot(kh, f, linestyle="-", linewidth=1, color="orange")
f = kh / np.tan(0.15)
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange")
f = kh / np.tan(0.075)
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange")
f = kh / np.tan(0.225)
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange")
# kb
th = np.linspace(0, np.pi / 2, 50)
ax.plot(kb * np.sin(th), kb * np.cos(th), linestyle="dashed", color="c")
# 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="c",
)
......@@ -69,13 +151,22 @@
)
cs = ax.scatter(
datab["R2"],
datab["R_waves"],
c=np.log10(datab["Fh"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
label=r"Without vortical modes",
vmin=-3,
vmax=1,
delta_kh = kh[1]
delta_kz = kz[1]
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
# 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,
)
)
......@@ -81,10 +172,14 @@
)
ax.set_xlim([1e-1, 1e5])
ax.set_xticks([1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5])
ax.set_xticklabels([1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5], fontsize=14)
ax.set_xscale("log")
ax.set_ylim([0, 1])
ax.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax.set_yticklabels(
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1.0$"], fontsize=14
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,
)
)
......@@ -90,8 +185,32 @@
)
ax.set_xlabel(r"$\mathcal{R}$", fontsize=20)
ax.set_ylabel(r"$\tilde{E}_{\rm wave}$", fontsize=20)
# ax.grid(True)
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,
)
return cs
fig, axes = plt.subplots(
......@@ -94,6 +213,6 @@
return cs
fig, axes = plt.subplots(
nrows=1, ncols=2, figsize=(10, 1.2 * 3 * 4.5 / 4), constrained_layout=True
ncols=2, nrows=1, figsize=(10, 1.2 * 3 * 4.5 / 4), constrained_layout=True
)
......@@ -99,7 +218,5 @@
)
#
ax0 = axes[0]
ax1 = axes[1]
......@@ -102,7 +219,10 @@
ax0 = axes[0]
ax1 = axes[1]
cs0 = plot_ratio_waves_vs_Fh(ax0, data=df, datab=df_proj)
cs1 = plot_ratio_waves_vs_R(ax1, data=df, datab=df_proj)
delta = 0.1
sim = get_sim(letter, proj=False)
# cs0 = plot_slices_omega_ki_waves(sim, ax0, key="Ee", key_k="kh", ik=3, delta=delta)
cs0 = plot_ratio_E_waves_vs_kh_kz(sim, ax0)
......@@ -108,44 +228,12 @@
# Emphasizing the points (N,Rb) = (40,20)
df_L = df[df["N"] == 40]
df_L = df_L[df_L["Rb"] == 20]
df_L_proj = df_proj[df_proj["N"] == 40]
df_L_proj = df_L_proj[df_L_proj["Rb"] == 20]
ax0.scatter(
df_L["Fh"],
df_L["R_waves"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax0.scatter(
df_L_proj["Fh"],
df_L_proj["R_waves"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax1.scatter(
df_L["R2"],
df_L["R_waves"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax1.scatter(
df_L_proj["R2"],
df_L_proj["R_waves"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
sim = get_sim(letter, proj=True)
# cs1 = plot_slices_omega_ki_waves(sim, ax1, key="Ee", key_k="kh", ik=3, delta=delta)
cs1 = plot_ratio_E_waves_vs_kh_kz(sim, ax1)
ax1.set_yticklabels([])
ax1.set_ylabel("")
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)
......@@ -150,8 +238,6 @@
ax0.legend(loc="upper right", fontsize=14)
ax0.set_title(r"$\rm (a)$", fontsize=20)
ax1.set_title(r"$\rm (b)$", fontsize=20)
fig.tight_layout()
fig.subplots_adjust(right=0.85)
......@@ -156,12 +242,9 @@
ax1.set_yticklabels("")
ax1.set_ylabel("")
fig.tight_layout()
cbar = fig.colorbar(cs0, ax=ax0, orientation="vertical")
cbar.set_ticks([-1, 0, 1, 2, 3, 4])
cbar.set_ticklabels(
[r"$-1$", r"$0$", r"$1$", r"$2$", r"$3$", r"$4$"], fontsize=14
cbar_ax2 = fig.add_axes([0.88, 0.17, 0.02, 0.65])
cbar2 = fig.colorbar(cs1, cax=cbar_ax2, cmap=cmbin, orientation="vertical")
cbar2.set_label(r"$\tilde{E}_{\rm wave}(k_h,k_z)$", fontsize=20)
cbar2.set_ticks([0, 0.2, 0.4, 0.6, 0.8, 1])
cbar2.set_ticklabels(
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1$"], fontsize=14
)
......@@ -167,11 +250,4 @@
)
fig.text(0.42, 0.07, r"$\log_{10} \mathcal{R}$", fontsize=20)
cbar = fig.colorbar(cs1, ax=ax1, orientation="vertical")
cbar.set_ticks([-3, -2, -1, 0, 1])
cbar.set_ticklabels([r"$-3$", r"$-2$", r"$-1$", r"$0$", r"$1$"], fontsize=14)
fig.text(0.88, 0.07, r"$\log_{10} F_h$", fontsize=20)
save_fig(fig, "figure15.png")
......
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm
......@@ -3,8 +4,5 @@
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from util import save_fig, Fh_limit, R2_limit
from util import save_fig
from util_dataframe import df, df_proj
......@@ -8,4 +6,7 @@
from util_dataframe import df, df_proj
cm = matplotlib.cm.get_cmap("inferno", 100)
# Latex
plt.rcParams["text.usetex"] = True
......@@ -11,2 +12,4 @@
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
......@@ -12,6 +15,41 @@
fig, ax = plt.subplots(
ncols=1, nrows=1, figsize=(7.5, 1.5 * 3 * 4.5 / 4), constrained_layout=True
)
def plot_ratio_waves_vs_Fh(ax, data=df, datab=df_proj):
ax.scatter(
data["Fh"],
data["ratio Ewaves"],
c=np.log10(data["R2"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
label=r"Standard Navier-Stokes",
vmin=-1,
vmax=4,
)
cs = ax.scatter(
datab["Fh"],
datab["ratio Ewaves"],
c=np.log10(datab["R2"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
label=r"Without vortical modes",
vmin=-1,
vmax=4,
)
ax.set_xlim([1e-3, 20])
ax.set_xticks([1e-3, 1e-2, 1e-1, 1e0, 1e1])
ax.set_xticklabels([1e-3, 1e-2, 1e-1, 1e0, 1e1], fontsize=14)
ax.set_xscale("log")
ax.set_ylim([0, 1])
ax.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax.set_yticklabels(
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1.0$"], fontsize=14
)
ax.set_xlabel(r"$F_h$", fontsize=20)
ax.set_ylabel(r"$\tilde{E}_{\rm wave}$", fontsize=20)
# ax.grid(True)
return cs
......@@ -16,15 +54,45 @@
cs = ax.scatter(
df["Fh"],
df["R2"],
# c=np.log10(df["R_waves"] * (df["Epolo"] + df["EA"]) / (df["Epolo"] + df["EA"] + df["Etoro"])),
# cmap="inferno",
color="k",
# edgecolors="k",
marker="o",
# vmin=-2,
# vmax=0,
# s=50 * df_proj["I_dissipation"] * 4 / 3, # 4/3 because of projection
label="Standard Navier-Stokes",
def plot_ratio_waves_vs_R(ax, data=df, datab=df_proj):
ax.scatter(
data["R2"],
data["ratio Ewaves"],
c=np.log10(data["Fh"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
label=r"Standard Navier-Stokes",
vmin=-3,
vmax=1,
)
cs = ax.scatter(
datab["R2"],
datab["ratio Ewaves"],
c=np.log10(datab["Fh"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
label=r"Without vortical modes",
vmin=-3,
vmax=1,
)
ax.set_xlim([1e-1, 1e5])
ax.set_xticks([1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5])
ax.set_xticklabels([1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5], fontsize=14)
ax.set_xscale("log")
ax.set_ylim([0, 1])
ax.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax.set_yticklabels(
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1.0$"], fontsize=14
)
ax.set_xlabel(r"$\mathcal{R}$", fontsize=20)
ax.set_ylabel(r"$\tilde{E}_{\rm wave}$", fontsize=20)
# ax.grid(True)
return cs
fig, axes = plt.subplots(
nrows=1, ncols=2, figsize=(10, 1.2 * 3 * 4.5 / 4), constrained_layout=True
)
......@@ -30,2 +98,3 @@
)
#
......@@ -31,10 +100,51 @@
cs = ax.scatter(
df_proj["Fh"],
df_proj["R2"],
color="k",
marker="^",
label="Without vortical modes",
ax0 = axes[0]
ax1 = axes[1]
cs0 = plot_ratio_waves_vs_Fh(ax0, data=df, datab=df_proj)
cs1 = plot_ratio_waves_vs_R(ax1, data=df, datab=df_proj)
# Emphasizing the points (N,Rb) = (40,20)
df_L = df[df["N"] == 40]
df_L = df_L[df_L["Rb"] == 20]
df_L_proj = df_proj[df_proj["N"] == 40]
df_L_proj = df_L_proj[df_L_proj["Rb"] == 20]
ax0.scatter(
df_L["Fh"],
df_L["ratio Ewaves"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax0.scatter(
df_L_proj["Fh"],
df_L_proj["ratio Ewaves"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax1.scatter(
df_L["R2"],
df_L["ratio Ewaves"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax1.scatter(
df_L_proj["R2"],
df_L_proj["ratio Ewaves"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
......@@ -38,13 +148,5 @@
)
# Le Reun et al. 2018
Re = np.array(
[332.0, 1279.0, 2049.0, 3673.0, 980.0, 1550.0, 2950.0, 500.0, 720.0]
)
N = np.array([1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 2.0, 4.0, 4.0])
Uh = np.array([4.0, 4.8, 4.4, 4.4, 3.7, 5.9, 5.4, 2.3, 3.3])
Uh *= 1e-3
EpsK = np.array([5.0, 3.0, 2.1, 1.8, 1.3, 4.0, 3.0, 1.0, 1.5])
EpsK *= 1e-8
ax0.legend(loc="upper right", fontsize=14)
......@@ -50,7 +152,5 @@
Fh = EpsK / (Uh**2 * N)
R = Re * Fh**2
cs = ax.scatter(Fh, R, color="k", marker="s", label="Le Reun et al. (2018)")
ax0.set_title(r"$\rm (a)$", fontsize=20)
ax1.set_title(r"$\rm (b)$", fontsize=20)
......@@ -55,49 +155,7 @@
# Table 2 of Rodda et al. 2022
Lh = 2.0 # 6.0 # m
nu = 1e-6 # m^2.s^-1
F = np.array(
[
0.14,
0.41,
0.68,
0.29,
0.54,
0.68,
0.36,
0.39,
0.41,
0.41,
0.68,
0.67,
0.68,
0.67,
0.68,
]
)
A = np.array(
[9.0, 4.0, 2.0, 9.0, 4.0, 3.0, 9.0, 9.0, 7.5, 9.0, 4.0, 5.0, 5.0, 9.0, 7.5]
)
A *= 1e-2 # m
N = np.array(
[
0.44,
0.45,
0.45,
0.36,
0.45,
0.45,
0.36,
0.26,
0.45,
0.36,
0.45,
0.26,
0.45,
0.26,
0.45,
]
) # rad.s^-1
Uf = F * A * N
ax1.set_yticklabels("")
ax1.set_ylabel("")
fig.tight_layout()
......@@ -103,20 +161,6 @@
Uxrms = np.array(
[
0.00077,
0.0047,
0.0058,
0.0062,
0.0097,
0.0127,
0.0045,
0.008,
0.011,
0.0054,
0.0095,
0.0056,
0.014,
0.0125,
0.0133,
]
cbar = fig.colorbar(cs0, ax=ax0, orientation="vertical")
cbar.set_ticks([-1, 0, 1, 2, 3, 4])
cbar.set_ticklabels(
[r"$-1$", r"$0$", r"$1$", r"$2$", r"$3$", r"$4$"], fontsize=14
)
......@@ -122,31 +166,4 @@
)
Uyrms = np.array(
[
0.00064,
0.0032,
0.0028,
0.0046,
0.0053,
0.0093,
0.0044,
0.0077,
0.011,
0.0053,
0.0086,
0.0064,
0.015,
0.0084,
0.0141,
]
)
Uh = Uf
EpsK = Uh**3 / Lh
Fh = EpsK / (Uh**2 * N)
R = EpsK / (nu * N**2)
cs = ax.scatter(Fh, R, color="k", marker="*", label="Rodda et al. (2022)")
fig.text(0.42, 0.07, r"$\log_{10} \mathcal{R}$", fontsize=20)
......@@ -151,115 +168,9 @@
ax.set_xlim([1e-3, 20])
ax.set_xscale("log")
ax.set_ylim([1e-1, 1e5])
ax.set_yscale("log")
ax.set_xlabel(r"$F_h$", fontsize=20)
ax.set_ylabel(r"$\mathcal{R}$", fontsize=20)
ax.set_xticks([1e-6, 1e-4, 1e-2, 1e0, 1e2])
ax.set_xticklabels(
[r"$10^{-6}$", r"$10^{-4}$", r"$10^{-2}$", r"$10^{0}$", r"$10^{2}$"],
fontsize=14,
)
ax.set_yticks([1e-5, 1e-3, 1e-1, 1e1, 1e3, 1e5])
ax.set_yticklabels(
[
r"$10^{-5}$",
r"$10^{-3}$",
r"$10^{-1}$",
r"$10^{1}$",
r"$10^{3}$",
r"$10^{5}$",
],
fontsize=14,
)
Fh_min, Fh_max = ax.get_xlim()
ax.axvline(1.0, linestyle="-", color="k")
# ax.axhline(1e0, linestyle="-", color="k")
Fh = np.array([Fh_min, 1])
ax.plot(Fh, [1, 1], linestyle="-", color="k")
ax.plot(
Fh,
Fh ** (6 / 5),
linestyle="-",
color="g",
label=r"$\mathcal{R} = F_h^{6/5}$",
)
ax.plot(Fh, Fh, linestyle="-", color="r", label=r"$\mathcal{R} = F_h$")
ax.plot(
Fh,
Fh ** (2 / 3),
linestyle="-",
color="m",
label=r"$\mathcal{R} = F_h^{2/3}$",
)
Fh = np.array([Fh_min, Fh_max])
ax.plot(
Fh,
Fh**2,
linestyle="-",
color="b",
label=r"$Re = \mathcal{R} F_h^{-2} = 1$",
)
ax.plot(
Fh,
500 * Fh**2,
linestyle="--",
color="b",
label=r"$Re = \mathcal{R} F_h^{-2} = 500$",
)
Fh_0 = 7e-5
Fh_1 = 3e-3
R_min = 1e-5
path = Path([[1, 1], [Fh_0, R_min], [Fh_1, R_min]])
patch = PathPatch(path, facecolor="none")
ax.add_patch(patch)
Fh_on_grid, R_on_grid = np.meshgrid(
np.logspace(-5, 0, 400), np.logspace(-5, 0, 400)
)
def func_color(F, R):
return np.log10(R / F**2)
im = plt.pcolormesh(
Fh_on_grid,
R_on_grid,
func_color(Fh_on_grid, R_on_grid),
cmap=plt.cm.Greys,
vmin=0.1 * func_color(Fh_1, R_min),
vmax=1.2 * func_color(Fh_0, R_min),
clip_path=patch,
clip_on=True,
zorder=0,
)
ax.text(3e0, 1e3, r"$Re > 1$", fontsize=14, color="b")
ax.text(1e1, 2e1, r"$Re <1$", fontsize=14, color="b")
"""
ax.text(1e-5, 10 ** (-1.5), r"$k_{\rm d} > k_{\rm b}$", fontsize=14, color="r")
ax.text(7e-2, 10 ** (-3.5), r"$k_{\rm d} < k_{\rm b}$", fontsize=14, color="r")
ax.text(1e-5, 1e-1, r"$k_{\eta} > k_{\rm b}$", fontsize=14, color="m")
ax.text(7e-2, 1e-3, r"$k_{\eta}< k_{\rm b}$", fontsize=14, color="m")
ax.text(1e-5, 1e-2, r"$\chi_{\rm max} > 1$", fontsize=14, color="g")
ax.text(7e-2, 1e-4, r"$\chi_{\rm max} < 1$", fontsize=14, color="g")
"""
ax.legend(loc="upper left", fontsize=10)
fig.tight_layout()
# fig.subplots_adjust(right=0.85, wspace=0.1)
cbar = fig.colorbar(cs1, ax=ax1, orientation="vertical")
cbar.set_ticks([-3, -2, -1, 0, 1])
cbar.set_ticklabels([r"$-3$", r"$-2$", r"$-1$", r"$0$", r"$1$"], fontsize=14)
fig.text(0.88, 0.07, r"$\log_{10} F_h$", fontsize=20)
save_fig(fig, "figure16.png")
......
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from util import save_fig, Fh_limit, R2_limit
from util_dataframe import df, df_proj
plt.rcParams["text.usetex"] = True
fig, ax = plt.subplots(
ncols=1, nrows=1, figsize=(7.5, 1.5 * 3 * 4.5 / 4), constrained_layout=True
)
cs = ax.scatter(
df["Fh"],
df["R2"],
# c=np.log10(df["R_waves"] * (df["Epolo"] + df["EA"]) / (df["Epolo"] + df["EA"] + df["Etoro"])),
# cmap="inferno",
color="k",
# edgecolors="k",
marker="o",
# vmin=-2,
# vmax=0,
# s=50 * df_proj["I_dissipation"] * 4 / 3, # 4/3 because of projection
label="Standard Navier-Stokes",
)
cs = ax.scatter(
df_proj["Fh"],
df_proj["R2"],
color="k",
marker="^",
label="Without vortical modes",
)
# Le Reun et al. 2018
Re = np.array(
[332.0, 1279.0, 2049.0, 3673.0, 980.0, 1550.0, 2950.0, 500.0, 720.0]
)
N = np.array([1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 2.0, 4.0, 4.0])
Uh = np.array([4.0, 4.8, 4.4, 4.4, 3.7, 5.9, 5.4, 2.3, 3.3])
Uh *= 1e-3
EpsK = np.array([5.0, 3.0, 2.1, 1.8, 1.3, 4.0, 3.0, 1.0, 1.5])
EpsK *= 1e-8
Fh = EpsK / (Uh**2 * N)
R = Re * Fh**2
cs = ax.scatter(Fh, R, color="k", marker="s", label="Le Reun et al. (2018)")
# Table 2 of Rodda et al. 2022
Lh = 2.0 # 6.0 # m
nu = 1e-6 # m^2.s^-1
F = np.array(
[
0.14,
0.41,
0.68,
0.29,
0.54,
0.68,
0.36,
0.39,
0.41,
0.41,
0.68,
0.67,
0.68,
0.67,
0.68,
]
)
A = np.array(
[9.0, 4.0, 2.0, 9.0, 4.0, 3.0, 9.0, 9.0, 7.5, 9.0, 4.0, 5.0, 5.0, 9.0, 7.5]
)
A *= 1e-2 # m
N = np.array(
[
0.44,
0.45,
0.45,
0.36,
0.45,
0.45,
0.36,
0.26,
0.45,
0.36,
0.45,
0.26,
0.45,
0.26,
0.45,
]
) # rad.s^-1
Uf = F * A * N
Uxrms = np.array(
[
0.00077,
0.0047,
0.0058,
0.0062,
0.0097,
0.0127,
0.0045,
0.008,
0.011,
0.0054,
0.0095,
0.0056,
0.014,
0.0125,
0.0133,
]
)
Uyrms = np.array(
[
0.00064,
0.0032,
0.0028,
0.0046,
0.0053,
0.0093,
0.0044,
0.0077,
0.011,
0.0053,
0.0086,
0.0064,
0.015,
0.0084,
0.0141,
]
)
Uh = Uf
EpsK = Uh**3 / Lh
Fh = EpsK / (Uh**2 * N)
R = EpsK / (nu * N**2)
cs = ax.scatter(Fh, R, color="k", marker="*", label="Rodda et al. (2022)")
# Brethouwer et al. 2007 Scaling analysis and simulation of strongly stratified turbulent flows
Fh = 1e-2 * np.array(
[
1.5,
1.6,
0.23,
0.53,
1.2,
1.5,
2.3,
4.2,
0.23,
0.45,
0.81,
1.5,
2.7,
2.6,
0.16,
0.39,
0.75,
0.93,
1.6,
]
)
R = np.array(
[
0.29,
0.69,
0.058,
0.21,
0.75,
1.75,
2.84,
9.3,
0.11,
0.40,
1.09,
2.97,
9.3,
15.6,
0.13,
0.47,
1.57,
4.2,
9.6,
]
)
cs = ax.scatter(Fh, R, color="k", marker="x", label="Brethouwer et al. (2007)")
# Waite and Bartello 2004 Stratified turbulence dominated by vortical motion
# Force vortical modes
# Waite and Bartello 2006 Stratified turbulence generated by internal gravity waves
# I don't see how to extract R from data
# Lindborg 2006 The energy cascade in a strongly stratified fluid
# Used a different viscosities on horizontal and vertical...
# Waite 2011 Stratified turbulence at the buoyancy scale
# Used a different viscosities on horizontal and vertical...
# Waite 2013 Potential enstrophy in stratified turbulence
# Excite only vortical modes
Fh = 1e-2 * np.array([2.2, 1.0, 0.43, 0.046, 2.1, 1.0, 0.45, 0.044, 2.1])
R = np.array([1.8, 0.48, 0.13, 0.002, 3.6, 0.94, 0.26, 0.0041, 1.8])
cs = ax.scatter(Fh, R, color="k", marker="4", label="Waite (2013)")
# Lam et al. 2021 Energy balance and mixing between waves and eddies in stably stratified turbulence
Fh = 1e-2 * np.array([3.5, 2.3, 1.3, 0.7, 0.4, 0.14, 0.045, 1.4, 0.95, 0.22])
R = np.array([11, 5, 1.8, 0.9, 0.5, 0.1, 0.01, 0.7, 0.35, 0.04])
cs = ax.scatter(Fh, R, color="k", marker="+", label="Lam et al. (2021)")
ax.set_xlim([1e-3, 20])
ax.set_xscale("log")
ax.set_ylim([1e-1, 1e5])
ax.set_yscale("log")
ax.set_xlabel(r"$F_h$", fontsize=20)
ax.set_ylabel(r"$\mathcal{R}$", fontsize=20)
ax.set_xticks([1e-6, 1e-4, 1e-2, 1e0, 1e2])
ax.set_xticklabels(
[r"$10^{-6}$", r"$10^{-4}$", r"$10^{-2}$", r"$10^{0}$", r"$10^{2}$"],
fontsize=14,
)
ax.set_yticks([1e-5, 1e-3, 1e-1, 1e1, 1e3, 1e5])
ax.set_yticklabels(
[
r"$10^{-5}$",
r"$10^{-3}$",
r"$10^{-1}$",
r"$10^{1}$",
r"$10^{3}$",
r"$10^{5}$",
],
fontsize=14,
)
Fh_min, Fh_max = ax.get_xlim()
ax.axvline(1.0, linestyle="-", color="k")
# ax.axhline(1e0, linestyle="-", color="k")
Fh = np.array([Fh_min, 1])
ax.plot(Fh, [1, 1], linestyle="-", color="k")
ax.plot(
Fh,
Fh ** (6 / 5),
linestyle="-",
color="g",
label=r"$\mathcal{R} = F_h^{6/5}$",
)
ax.plot(Fh, Fh, linestyle="-", color="r", label=r"$\mathcal{R} = F_h$")
ax.plot(
Fh,
Fh ** (2 / 3),
linestyle="-",
color="m",
label=r"$\mathcal{R} = F_h^{2/3}$",
)
Fh = np.array([Fh_min, Fh_max])
ax.plot(
Fh,
Fh**2,
linestyle="-",
color="b",
label=r"$Re = \mathcal{R} F_h^{-2} = 1$",
)
ax.plot(
Fh,
500 * Fh**2,
linestyle="--",
color="b",
label=r"$Re = \mathcal{R} F_h^{-2} = 500$",
)
Fh_0 = 7e-5
Fh_1 = 3e-3
R_min = 1e-5
path = Path([[1, 1], [Fh_0, R_min], [Fh_1, R_min]])
patch = PathPatch(path, facecolor="none")
ax.add_patch(patch)
Fh_on_grid, R_on_grid = np.meshgrid(
np.logspace(-5, 0, 400), np.logspace(-5, 0, 400)
)
def func_color(F, R):
return np.log10(R / F**2)
im = plt.pcolormesh(
Fh_on_grid,
R_on_grid,
func_color(Fh_on_grid, R_on_grid),
cmap=plt.cm.Greys,
vmin=0.1 * func_color(Fh_1, R_min),
vmax=1.2 * func_color(Fh_0, R_min),
clip_path=patch,
clip_on=True,
zorder=0,
)
ax.text(3e0, 1e3, r"$Re > 1$", fontsize=14, color="b")
ax.text(1e1, 2e1, r"$Re <1$", fontsize=14, color="b")
"""
ax.text(1e-5, 10 ** (-1.5), r"$k_{\rm d} > k_{\rm b}$", fontsize=14, color="r")
ax.text(7e-2, 10 ** (-3.5), r"$k_{\rm d} < k_{\rm b}$", fontsize=14, color="r")
ax.text(1e-5, 1e-1, r"$k_{\eta} > k_{\rm b}$", fontsize=14, color="m")
ax.text(7e-2, 1e-3, r"$k_{\eta}< k_{\rm b}$", fontsize=14, color="m")
ax.text(1e-5, 1e-2, r"$\chi_{\rm max} > 1$", fontsize=14, color="g")
ax.text(7e-2, 1e-4, r"$\chi_{\rm max} < 1$", fontsize=14, color="g")
"""
ax.legend(loc="upper left", fontsize=10)
fig.tight_layout()
# fig.subplots_adjust(right=0.85, wspace=0.1)
save_fig(fig, "figure17.png")
if __name__ == "__main__":
plt.show()
......@@ -34,7 +34,7 @@
marker="o",
vmin=0.0,
vmax=1.0,
s=50 * df_proj["I_dissipation"] * 4 / 3, # 4/3 because of projection
s=50 * df_proj["I_dissipation"],
)
......@@ -80,7 +80,7 @@
ax1.scatter(df_L_proj["Fh"], df_L_proj["R2"], edgecolors="r", facecolors='none', s=60, marker="s", linewidths=2)
"""
#ax0.text(2e-3, 1e3, r"Strongly", fontsize=14, alpha=0.5)
# ax0.text(2e-3, 1e3, r"Strongly", fontsize=14, alpha=0.5)
ax0.text(2e-3, 3e2, r"LAST", fontsize=14, alpha=0.5)
ax0.text(1.5e-1, 1e4, r"Weakly", fontsize=14, alpha=0.5)
......
......@@ -6,6 +6,10 @@
from util import save_fig
# The two following lines are used to do a figure for reviews # TODO: erase it when the manuscript is accepted for publication
# df = df[df["k_max*eta"] > 1.0].copy()
# df_proj = df_proj[df_proj["k_max*eta"] > 1.0].copy()
plt.rcParams["text.usetex"] = True
fig, axes = plt.subplots(
......
import sys
import h5py
import matplotlib.cm
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
from math import degrees
from util_simuls_regimes import get_sim
from fluidsim.util import load_params_simul, times_start_last_from_path
from util import (
compute_kf_kb_ko_keta_kd,
compute_omega_emp_vs_kzkh,
customize,
paths_simuls_regimes,
paths_simuls_regimes_proj,
save_fig,
)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
cm = matplotlib.cm.get_cmap("inferno", 100)
print(sys.argv)
letter = sys.argv[-1]
if letter not in "DLOWPU":
letter = "L"
def plot_omega_spectra(sim, ax):
path = sim.params.path_run
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)
params = load_params_simul(path)
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"][:]
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
Ee = 2 * np.minimum(EA, Epolo)
EA = np.sum(EA, axis=0)
Epolo = np.sum(Epolo, axis=0)
Etoro = np.sum(Etoro, axis=0)
Ee = np.sum(Ee, axis=0)
E = np.sum(E, axis=0)
EA = np.sum(EA, axis=0)
Epolo = np.sum(Epolo, axis=0)
Etoro = np.sum(Etoro, axis=0)
Ee = np.sum(Ee, axis=0)
E = np.sum(E, axis=0)
coef_compensate = 0
cs = ax.plot(
omegas / N,
EA * omegas**coef_compensate,
color="b",
label=r"$E_{\rm pot}(\omega)$",
)
cs = ax.plot(
omegas / N,
Epolo * omegas**coef_compensate,
color="g",
label=r"$E_{\rm polo}(\omega)$",
)
if proj == None:
cs = ax.plot(
omegas / N,
Etoro * omegas**coef_compensate,
color="r",
label=r"$E_{\rm toro}(\omega)$",
)
om = np.array([0.1 * N, N])
ax.plot(
om / N,
1e-4 * (om / N) ** (-2 + coef_compensate),
"--",
color="gray",
label=None,
)
ax.text(0.5, 1e-3, r"$\omega^{-2}$", color="gray", fontsize=14)
ax.plot(
om / N,
1e-5 * (om / N) ** (-3 / 2 + coef_compensate),
"-.",
color="gray",
label=None,
)
ax.text(0.5, 5e-6, r"$\omega^{-3/2}$", color="gray", fontsize=14)
# Forcing
angle = sim.params.forcing.tcrandom_anisotropic.angle
delta_angle = sim.params.forcing.tcrandom_anisotropic.delta_angle
omega_fmin = N * np.sin(angle - 0.5 * delta_angle)
omega_fmax = N * np.sin(angle + 0.5 * delta_angle)
ax.axvline(omega_fmin / N, color="orange", linestyle="dashed")
ax.axvline(omega_fmax / N, color="orange", linestyle="dashed")
ax.set_xlabel(r"$\omega/N$", fontsize=20)
ax.set_xscale("log")
ax.set_xticks([1e-1, 1e0])
ax.set_xticklabels([r"$10^{-1}$", r"$10^{0}$"], fontsize=14)
ax.set_xlim([min(omegas) / N, max(omegas) / N])
ax.set_ylabel(r"$E(\omega)$", fontsize=20)
ax.set_yscale("log")
ax.set_yticks([1e-6, 1e-5, 1e-4, 1e-3, 1e-2])
ax.set_yticklabels(
[
r"$10^{-6}$",
r"$10^{-5}$",
r"$10^{-4}$",
r"$10^{-3}$",
r"$10^{-2}$",
],
fontsize=14,
)
ax.legend(loc="upper right", fontsize=14)
return cs
sim = get_sim(letter)
sim_proj = get_sim(letter, proj=True)
# TODO: uncomment this assert
# assert (
# sim.params.oper.nx == sim_proj.params.oper.nx
# ), f"Not the same resolution for simulation Without vortical modes: {sim.params.oper.nx} vs {sim_proj.params.oper.nx}"
fig, axes = plt.subplots(
ncols=2, nrows=1, figsize=(10, 1.2 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0]
ax1 = axes[1]
cs0 = plot_omega_spectra(sim, ax0)
cs1 = plot_omega_spectra(sim_proj, ax1)
ax1.set_yticklabels([])
ax1.set_ylabel("")
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)
fig.tight_layout()
save_fig(fig, f"figure4.png")
if __name__ == "__main__":
plt.show()
import sys
import numpy as np
import h5py
import matplotlib.cm
import matplotlib.pyplot as plt
......@@ -4,4 +5,6 @@
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
from math import degrees
from util_simuls_regimes import get_sim
......@@ -6,8 +9,17 @@
from util_simuls_regimes import get_sim
from util import save_fig, compute_kf_kb_ko_keta_kd, customize
from fluidsim.util import load_params_simul, times_start_last_from_path
from util import (
compute_kf_kb_ko_keta_kd,
compute_omega_emp_vs_kzkh,
customize,
paths_simuls_regimes,
paths_simuls_regimes_proj,
save_fig,
)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
......@@ -9,8 +21,11 @@
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
cm = matplotlib.cm.get_cmap("inferno", 100)
print(sys.argv)
letter = sys.argv[-1]
......@@ -18,24 +33,25 @@
letter = "L"
fig, axes = plt.subplots(
ncols=2, nrows=2, figsize=(10, 1.2 * 2 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
coef_compensate = 5 / 3
# Standard Navier-Stokes
sim = get_sim(letter)
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)
kh = data["kh_spectra"]
kz = data["kz"]
def plot_omega_spectra(sim, ax):
path = sim.params.path_run
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)
params = load_params_simul(path)
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"][:]
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
Ee = 2 * np.minimum(EA, Epolo)
......@@ -41,10 +57,7 @@
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 = np.sum(EA, axis=0)
Epolo = np.sum(Epolo, axis=0)
Etoro = np.sum(Etoro, axis=0)
Ee = np.sum(Ee, axis=0)
E = np.sum(E, axis=0)
......@@ -50,19 +63,7 @@
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
EKhr = data["Khr"]
Epolo = EKhd + EKz
Etoro = EKhr
EKh = EKhr + EKhd
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
Epot_vs_kh = np.sum(EA, axis=0) * delta_kz
Epot_vs_kz = np.sum(EA, axis=1) * delta_kh
EA = np.sum(EA, axis=0)
Epolo = np.sum(Epolo, axis=0)
Etoro = np.sum(Etoro, axis=0)
Ee = np.sum(Ee, axis=0)
E = np.sum(E, axis=0)
......@@ -68,38 +69,31 @@
ax0.plot(
kh,
Epolo_vs_kh * kh ** (coef_compensate),
"-g",
label=r"$E_{\rm polo}(k_h)$",
)
ax0.plot(
kh,
Etoro_vs_kh * kh ** (coef_compensate),
"-r",
label=r"$E_{\rm toro}(k_h)$",
)
ax0.plot(
kh,
Epot_vs_kh * kh ** (coef_compensate),
"-b",
label=r"$E_{\rm pot}(k_h)$",
)
ax0.plot(
kz,
Epolo_vs_kz * kz ** (coef_compensate),
"--g",
label=r"$E_{\rm polo}(k_z)$",
)
ax0.plot(
kz,
Etoro_vs_kz * kz ** (coef_compensate),
"--r",
label=r"$E_{\rm toro}(k_z)$",
)
ax0.plot(
kz,
Epot_vs_kz * kz ** (coef_compensate),
"--b",
label=r"$E_{\rm pot}(k_z)$",
)
coef_compensate = 0
cs = ax.plot(
omegas / N,
EA * omegas**coef_compensate,
color="b",
label=r"$E_{\rm pot}(\omega)$",
)
cs = ax.plot(
omegas / N,
Epolo * omegas**coef_compensate,
color="g",
label=r"$E_{\rm polo}(\omega)$",
)
if proj == None:
cs = ax.plot(
omegas / N,
Etoro * omegas**coef_compensate,
color="r",
label=r"$E_{\rm toro}(\omega)$",
)
om = np.array([0.1 * N, N])
ax.plot(
om / N,
1e-4 * (om / N) ** (-2 + coef_compensate),
"--",
color="gray",
label=None,
)
ax.text(0.5, 1e-3, r"$\omega^{-2}$", color="gray", fontsize=14)
......@@ -105,17 +99,40 @@
sim.output.spect_energy_budg.plot_fluxes(
tmin=tmin, key_k="kh", ax=ax2, plot_conversion=False
)
for ax in [ax0, ax2]:
ax.axvline(kb, color="k", linestyle="dotted")
ax.axvline(ko, color="k", linestyle="dashed")
ax.axvline(kf, color="orange", linestyle="dashed")
ax.set_xlim([kh[1], max(kh)])
ax.set_xticks([1e1, 1e2, 1e3, kb, ko])
ax.set_xticklabels(
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
fontsize=14,
)
# sim.output.spect_energy_budg.plot_fluxes(tmin=tmin, key_k="kz", ax=ax4)
ax.plot(
om / N,
1e-5 * (om / N) ** (-3 / 2 + coef_compensate),
"-.",
color="gray",
label=None,
)
ax.text(0.5, 5e-6, r"$\omega^{-3/2}$", color="gray", fontsize=14)
# Forcing
angle = sim.params.forcing.tcrandom_anisotropic.angle
delta_angle = sim.params.forcing.tcrandom_anisotropic.delta_angle
omega_fmin = N * np.sin(angle - 0.5 * delta_angle)
omega_fmax = N * np.sin(angle + 0.5 * delta_angle)
ax.axvline(omega_fmin / N, color="orange", linestyle="dashed")
ax.axvline(omega_fmax / N, color="orange", linestyle="dashed")
ax.set_xlabel(r"$\omega/N$", fontsize=20)
ax.set_xscale("log")
ax.set_xticks([1e-1, 1e0])
ax.set_xticklabels([r"$10^{-1}$", r"$10^{0}$"], fontsize=14)
ax.set_xlim([min(omegas) / N, max(omegas) / N])
ax.set_ylabel(r"$E(\omega)$", fontsize=20)
ax.set_yscale("log")
ax.set_yticks([1e-6, 1e-5, 1e-4, 1e-3, 1e-2])
ax.set_yticklabels(
[
r"$10^{-6}$",
r"$10^{-5}$",
r"$10^{-4}$",
r"$10^{-3}$",
r"$10^{-2}$",
],
fontsize=14,
)
ax.legend(loc="upper right", fontsize=14)
return cs
......@@ -120,12 +137,5 @@
# Without vortical modes
sim = get_sim(letter, proj=True)
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
data = sim.output.spectra.load1d_mean(tmin)
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
data = sim.output.spectra.load_kzkh_mean(tmin)
kh = data["kh_spectra"]
kz = data["kz"]
sim = get_sim(letter)
sim_proj = get_sim(letter, proj=True)
......@@ -131,28 +141,2 @@
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
EKh = EKhr + EKhd
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
Epot_vs_kh = np.sum(EA, axis=0) * delta_kz
Epot_vs_kz = np.sum(EA, axis=1) * delta_kh
......@@ -158,42 +142,7 @@
ax1.plot(
kh,
Epolo_vs_kh * kh ** (coef_compensate),
"-g",
label=r"$E_{\rm polo}(k_h)$",
)
ax1.plot(
kh,
Epot_vs_kh * kh ** (coef_compensate),
"-b",
label=r"$E_{\rm pot}(k_h)$",
)
ax1.plot(
kz,
Epolo_vs_kz * kz ** (coef_compensate),
"--g",
label=r"$E_{\rm polo}(k_z)$",
)
ax1.plot(
kz,
Epot_vs_kz * kz ** (coef_compensate),
"--b",
label=r"$E_{\rm pot}(k_z)$",
)
tmp = sim.output.spect_energy_budg.plot_fluxes(
tmin=tmin, key_k="kh", ax=ax3, plot_conversion=False
)
for ax in [ax1, ax3]:
ax.axvline(kb, color="k", linestyle="dotted")
ax.axvline(ko, color="k", linestyle="dashed")
ax.axvline(kf, color="orange", linestyle="dashed")
ax.set_xlim([kh[1], max(kh)])
ax.set_xticks([1e1, 1e2, 1e3, kb, ko])
ax.set_xticklabels(
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
fontsize=14,
)
# tmp = sim.output.spect_energy_budg.plot_fluxes(tmin=tmin, key_k="kz", ax=ax5)
# TODO: uncomment this assert
# assert (
# sim.params.oper.nx == sim_proj.params.oper.nx
# ), f"Not the same resolution for simulation Without vortical modes: {sim.params.oper.nx} vs {sim_proj.params.oper.nx}"
......@@ -198,32 +147,6 @@
for ax in [ax0, ax1]:
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_ylim([1e-3, 1e1])
ax.set_xlabel(r"$k_h, k_z$", fontsize=20)
N = sim.params.N
mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
Uh2 = mean_values["Uh2"]
Uh = np.sqrt(Uh2)
epsK = mean_values["epsK"]
# Plot k^-2 and k^-3
k = np.array([2e1, 8e2])
ax.plot(
k,
4.0 * (k / k[0]) ** (-2.0 + coef_compensate),
"--",
color="gray",
label=None,
)
ax.text(5e2, 3.0, r"$\propto k_i^{-2}$", fontsize=14, color="gray")
ax.plot(
k,
4.0 * (k / k[0]) ** (-3.0 + coef_compensate),
"-.",
color="gray",
label=None,
)
ax.text(5e2, 1e-2, r"$\propto k_i^{-3}$", fontsize=14, color="gray")
fig, axes = plt.subplots(
ncols=2, nrows=1, figsize=(10, 1.2 * 3 * 4.5 / 4), constrained_layout=True
)
......@@ -229,63 +152,5 @@
"""
ax.plot(
k,
0.5 * (epsK**(2/3)) * k ** (-5/3 + coef_compensate),
"-.",
color="gray",
label=None,
)
ax.plot(
k,
0.02 * N * Uh * k ** (-5/3 + coef_compensate),
"-.",
color="gray",
label=None,
)
"""
for ax in [ax2, ax3]:
ax.set_ylim([-0.2, 1.2])
lines = [line for line in ax.get_lines()]
for l in range(3):
lines[l].remove()
labels = [
r"$\Pi_{\rm kin}(k_h)/\varepsilon$",
r"$\Pi_{\rm pot}(k_h)/\varepsilon$",
r"$\varepsilon_{\rm kin}(k_h)/\varepsilon$",
r"$\varepsilon_{\rm pot}(k_h)/\varepsilon$",
]
ax.legend(loc="upper left", fontsize=10, labels=labels)
for ax in [ax2]:
ax.set_yticks([-0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2])
ax.set_yticklabels(
[
r"$-0.2$",
r"$0$",
r"$0.2$",
r"$0.4$",
r"$0.6$",
r"$0.8$",
r"$1.0$",
r"$1.2$",
],
fontsize=14,
)
for ax in [ax1, ax3]:
ax.set_yticklabels([])
ax.set_ylabel("")
ax0.set_yticks([1e-3, 1e-2, 1e-1, 1e0, 1e1])
ax0.set_yticklabels(
[r"$10^{-3}$", r"$10^{-2}$", r"$10^{-1}$", r"$10^{0}$", r"$10^{1}$"],
fontsize=14,
)
ax0.set_ylabel(r"$E_{\rm 1D} \times k_i^{5/3}$", fontsize=20)
for ax in [ax0, ax1]:
ax.legend(loc="lower center", fontsize=10)
# ax.grid(True)
ax0 = axes[0]
ax1 = axes[1]
......@@ -290,7 +155,5 @@
ax2.set_xlabel(r"$k_h$", fontsize=20)
ax3.set_xlabel(r"$k_h$", fontsize=20)
# ax4.set_xlabel(r"$k_z$", fontsize=20)
# ax5.set_xlabel(r"$k_z$", fontsize=20)
cs0 = plot_omega_spectra(sim, ax0)
cs1 = plot_omega_spectra(sim_proj, ax1)
......@@ -296,6 +159,7 @@
ax2.set_ylabel(r"$\Pi(k_h)/ \varepsilon$", fontsize=20)
# ax4.set_ylabel(r"$\Pi(k_z)/ \varepsilon$", fontsize=20)
ax1.set_yticklabels([])
ax1.set_ylabel("")
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)
......@@ -299,10 +163,6 @@
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 (d)$", fontsize=20)
# ax4.set_title(r"$\rm (e)$", fontsize=20)
# ax5.set_title(r"$\rm (f)$", fontsize=20)
fig.tight_layout()
......@@ -306,5 +166,7 @@
fig.tight_layout()
save_fig(fig, f"figure5.png")
......@@ -309,4 +171,5 @@
save_fig(fig, f"figure5.png")
if __name__ == "__main__":
plt.show()
import matplotlib.cm
import matplotlib.pyplot as plt
import numpy as np
import sys
......@@ -4,8 +2,6 @@
from util import (
compute_kf_kb_ko_keta_kd,
save_fig,
)
import numpy as np
import matplotlib.pyplot as plt
from util_simuls_regimes import get_sim
......@@ -9,9 +5,9 @@
from util_simuls_regimes import get_sim
cm = matplotlib.cm.get_cmap("inferno", 100)
from util import save_fig, compute_kf_kb_ko_keta_kd, customize
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
......@@ -13,48 +9,8 @@
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
def plot_spectra(sim, ax, key="Ee", key_k="kh"):
p_oper = sim.params.oper
kmax = p_oper.coef_dealiasing * p_oper.nx * np.pi / p_oper.Lx
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]
if key == "EA":
spectrum = data["A"]
elif key == "EK":
EKhd = data["Khd"]
EKz = data["Kz"]
Etoro = data["Khr"]
spectrum = EKhd + EKz + Etoro
elif key == "Epolo":
EKhd = data["Khd"]
EKz = data["Kz"]
spectrum = EKhd + EKz
elif key == "Ee":
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
Epolo = EKhd + EKz
spectrum = 2 * np.minimum(EA, Epolo)
elif key == "Ed":
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
Epolo = EKhd + EKz
spectrum = EA + Epolo - 2 * np.minimum(EA, Epolo)
elif key == "Etoro":
spectrum = data["Khr"]
else:
print(f"Don't know key: {key} \n")
exit
print(sys.argv)
letter = sys.argv[-1]
......@@ -60,56 +16,5 @@
if key_k == "kh":
for iz in range(int(np.floor(np.log(len(kz)) / np.log(1.2)))):
nz = int(np.floor(1.2**iz))
if kz[nz] <= 10 * kb and kz[nz] >= 0.1 * kb:
cs = ax.plot(
kh,
spectrum[nz, :],
color=cm(0.5 + 0.5 * np.log10((kz[nz] / kb))),
linestyle="-",
)
ks = np.array([0.1 * kf, kb])
ax.plot(ks, 2e-1 * ks ** (-2), "k-")
ax.text(np.sqrt(0.1 * kf * kb), 2e-3, r"$\propto k_h^{-2}$", fontsize=14)
ks = np.array([0.1 * kf, kb])
ax.plot(ks, 1e-10 * ks ** (1), "k-")
ax.text(np.sqrt(0.1 * kf * kb), 3e-10, r"$\propto k_h^{1}$", fontsize=14)
ax.set_xlim([delta_kh, kmax])
ax.set_xlabel(r"$k_h$", fontsize=20)
elif key_k == "kz":
for ih in range(int(np.floor(np.log(len(kh)) / np.log(1.2)))):
nh = int(np.floor(1.2**ih))
if kh[nh] <= 10 * kb and kh[nh] >= 0.1 * kb:
cs = ax.plot(
kz,
spectrum[:, nh],
color=cm(0.5 + 0.5 * np.log10((kh[nh] / kb))),
linestyle="-",
)
ks = np.array([0.3 * kf, ko])
ax.plot(ks, 5e-9 * ks ** (0), "k-")
ax.text(np.sqrt(0.3 * kf * ko), 7e-10, r"$\propto k_z^{0}$", fontsize=14)
ks = np.array([kb, ko])
ax.plot(ks, 1e4 * ks ** (-4), "k-")
ax.text(np.sqrt(kb * ko), 1e-4, r"$\propto k_z^{-4}$", fontsize=14)
ax.set_xlim([delta_kz, kmax])
ax.set_xlabel(r"$k_z$", fontsize=20)
ax.set_ylim([1e-10, 1e-1])
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_xscale("log")
ax.set_yscale("log")
# ax.grid(True)
ax.axvline(kb, color="k", linestyle="dotted")
ax.axvline(ko, color="k", linestyle="dashed")
ax.axvline(kf, color="orange", linestyle="dashed")
ax.axvline(keta, color="g", linestyle="dashed")
return cs
if letter not in "DLOWPU":
letter = "L"
......@@ -114,6 +19,3 @@
### Figure: Kinetic energy
nbax = 0
css = [None for i in range(4)]
fig, axes = plt.subplots(
......@@ -119,5 +21,5 @@
fig, axes = plt.subplots(
ncols=2, nrows=2, figsize=(10, 2 * 3 * 4.5 / 4), constrained_layout=True
ncols=2, nrows=2, figsize=(10, 1.2 * 2 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
......@@ -125,6 +27,95 @@
ax2 = axes[1, 0]
ax3 = axes[1, 1]
axs = [ax0, ax1, ax2, ax3]
coef_compensate = 5 / 3
# Standard Navier-Stokes
sim = get_sim(letter)
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)
kh = data["kh_spectra"]
kz = data["kz"]
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
EKh = EKhr + EKhd
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
Epot_vs_kh = np.sum(EA, axis=0) * delta_kz
Epot_vs_kz = np.sum(EA, axis=1) * delta_kh
ax0.plot(
kh,
Epolo_vs_kh * kh ** (coef_compensate),
"-g",
label=r"$E_{\rm polo}(k_h)$",
)
ax0.plot(
kh,
Etoro_vs_kh * kh ** (coef_compensate),
"-r",
label=r"$E_{\rm toro}(k_h)$",
)
ax0.plot(
kh,
Epot_vs_kh * kh ** (coef_compensate),
"-b",
label=r"$E_{\rm pot}(k_h)$",
)
ax0.plot(
kz,
Epolo_vs_kz * kz ** (coef_compensate),
"--g",
label=r"$E_{\rm polo}(k_z)$",
)
ax0.plot(
kz,
Etoro_vs_kz * kz ** (coef_compensate),
"--r",
label=r"$E_{\rm toro}(k_z)$",
)
ax0.plot(
kz,
Epot_vs_kz * kz ** (coef_compensate),
"--b",
label=r"$E_{\rm pot}(k_z)$",
)
sim.output.spect_energy_budg.plot_fluxes(
tmin=tmin, key_k="kh", ax=ax2, plot_conversion=False
)
for ax in [ax0, ax2]:
ax.axvline(kb, color="k", linestyle="dotted")
ax.axvline(ko, color="k", linestyle="dashed")
ax.axvline(kf, color="orange", linestyle="dashed")
ax.set_xlim([kh[1], max(kh)])
ax.set_xticks([1e1, 1e2, 1e3, kb, ko])
ax.set_xticklabels(
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
fontsize=14,
)
# sim.output.spect_energy_budg.plot_fluxes(tmin=tmin, key_k="kz", ax=ax4)
......@@ -129,9 +120,80 @@
for proj in [False, True]:
sim = get_sim("L", proj=proj)
css[nbax] = plot_spectra(sim, axs[nbax], key="EK", key_k="kh")
css[nbax + 2] = plot_spectra(sim, axs[nbax + 2], key="EK", key_k="kz")
nbax += 1
# Without vortical modes
sim = get_sim(letter, proj=True)
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
data = sim.output.spectra.load1d_mean(tmin)
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
data = sim.output.spectra.load_kzkh_mean(tmin)
kh = data["kh_spectra"]
kz = data["kz"]
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
EKh = EKhr + EKhd
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
Epot_vs_kh = np.sum(EA, axis=0) * delta_kz
Epot_vs_kz = np.sum(EA, axis=1) * delta_kh
ax1.plot(
kh,
Epolo_vs_kh * kh ** (coef_compensate),
"-g",
label=r"$E_{\rm polo}(k_h)$",
)
ax1.plot(
kh,
Epot_vs_kh * kh ** (coef_compensate),
"-b",
label=r"$E_{\rm pot}(k_h)$",
)
ax1.plot(
kz,
Epolo_vs_kz * kz ** (coef_compensate),
"--g",
label=r"$E_{\rm polo}(k_z)$",
)
ax1.plot(
kz,
Epot_vs_kz * kz ** (coef_compensate),
"--b",
label=r"$E_{\rm pot}(k_z)$",
)
tmp = sim.output.spect_energy_budg.plot_fluxes(
tmin=tmin, key_k="kh", ax=ax3, plot_conversion=False
)
for ax in [ax1, ax3]:
ax.axvline(kb, color="k", linestyle="dotted")
ax.axvline(ko, color="k", linestyle="dashed")
ax.axvline(kf, color="orange", linestyle="dashed")
ax.set_xlim([kh[1], max(kh)])
ax.set_xticks([1e1, 1e2, 1e3, kb, ko])
ax.set_xticklabels(
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
fontsize=14,
)
# tmp = sim.output.spect_energy_budg.plot_fluxes(tmin=tmin, key_k="kz", ax=ax5)
......@@ -136,12 +198,80 @@
for ax in [ax0, ax2]:
ax.set_ylabel(r"$E_{\rm kin}(k_h, k_z)$", fontsize=20)
ax.set_yticks(
[1e-10, 1e-8, 1e-6, 1e-4, 1e-2],
[r"$10^{-10}$", r"$10^{-8}$", r"$10^{-6}$", r"$10^{-4}$", r"$10^{-2}$"],
for ax in [ax0, ax1]:
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_ylim([1e-3, 1e1])
ax.set_xlabel(r"$k_h, k_z$", fontsize=20)
N = sim.params.N
mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
Uh2 = mean_values["Uh2"]
Uh = np.sqrt(Uh2)
epsK = mean_values["epsK"]
# Plot k^-2 and k^-3
k = np.array([2e1, 8e2])
ax.plot(
k,
4.0 * (k / k[0]) ** (-2.0 + coef_compensate),
"--",
color="gray",
label=None,
)
ax.text(5e2, 3.0, r"$\propto k_i^{-2}$", fontsize=14, color="gray")
ax.plot(
k,
4.0 * (k / k[0]) ** (-3.0 + coef_compensate),
"-.",
color="gray",
label=None,
)
ax.text(5e2, 1e-2, r"$\propto k_i^{-3}$", fontsize=14, color="gray")
"""
ax.plot(
k,
0.5 * (epsK**(2/3)) * k ** (-5/3 + coef_compensate),
"-.",
color="gray",
label=None,
)
ax.plot(
k,
0.02 * N * Uh * k ** (-5/3 + coef_compensate),
"-.",
color="gray",
label=None,
)
"""
for ax in [ax2, ax3]:
ax.set_ylim([-0.2, 1.2])
lines = [line for line in ax.get_lines()]
for l in range(3):
lines[l].remove()
labels = [
r"$\Pi_{\rm kin}(k_h)/\varepsilon$",
r"$\Pi_{\rm pot}(k_h)/\varepsilon$",
r"$\varepsilon_{\rm kin}(k_h)/\varepsilon$",
r"$\varepsilon_{\rm pot}(k_h)/\varepsilon$",
]
ax.legend(loc="upper left", fontsize=10, labels=labels)
for ax in [ax2]:
ax.set_yticks([-0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2])
ax.set_yticklabels(
[
r"$-0.2$",
r"$0$",
r"$0.2$",
r"$0.4$",
r"$0.6$",
r"$0.8$",
r"$1.0$",
r"$1.2$",
],
fontsize=14,
)
for ax in [ax1, ax3]:
ax.set_yticklabels([])
......@@ -143,10 +273,31 @@
fontsize=14,
)
for ax in [ax1, ax3]:
ax.set_yticklabels([])
ax.set_ylabel("")
ax0.set_yticks([1e-3, 1e-2, 1e-1, 1e0, 1e1])
ax0.set_yticklabels(
[r"$10^{-3}$", r"$10^{-2}$", r"$10^{-1}$", r"$10^{0}$", r"$10^{1}$"],
fontsize=14,
)
ax0.set_ylabel(r"$E_{\rm 1D} \times k_i^{5/3}$", fontsize=20)
for ax in [ax0, ax1]:
ax.legend(loc="lower center", fontsize=10)
# ax.grid(True)
ax2.set_xlabel(r"$k_h$", fontsize=20)
ax3.set_xlabel(r"$k_h$", fontsize=20)
# ax4.set_xlabel(r"$k_z$", fontsize=20)
# ax5.set_xlabel(r"$k_z$", fontsize=20)
ax2.set_ylabel(r"$\Pi(k_h)/ \varepsilon$", fontsize=20)
# ax4.set_ylabel(r"$\Pi(k_z)/ \varepsilon$", fontsize=20)
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 (d)$", fontsize=20)
......@@ -148,7 +299,8 @@
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 (d)$", fontsize=20)
# ax4.set_title(r"$\rm (e)$", fontsize=20)
# ax5.set_title(r"$\rm (f)$", fontsize=20)
......@@ -154,4 +306,3 @@
norm = matplotlib.colors.Normalize(vmin=-1, vmax=1)
fig.tight_layout()
......@@ -156,32 +307,6 @@
fig.tight_layout()
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.1, 0.02, 0.33])
cbar = fig.colorbar(
matplotlib.cm.ScalarMappable(norm=norm, cmap=cm),
cax=cbar_ax,
cmap=cm,
orientation="vertical",
)
cbar.set_ticks([-1, -0.5, 0.0, 0.5, 1])
cbar.set_ticklabels([r"$-1$", r"$-0.5$", r"$0.0$", r"$0.5$", r"$1$"], fontsize=14)
cbar.set_label(r"$\log_{10} (k_h/k_{\rm b})$", fontsize=20)
cbar_ax = fig.add_axes([0.88, 0.565, 0.02, 0.33])
cbar = fig.colorbar(
matplotlib.cm.ScalarMappable(norm=norm, cmap=cm),
cax=cbar_ax,
cmap=cm,
orientation="vertical",
)
cbar.set_ticks([-1, -0.5, 0.0, 0.5, 1])
cbar.set_ticklabels([r"$-1$", r"$-0.5$", r"$0.0$", r"$0.5$", r"$1$"], fontsize=14)
cbar.set_label(r"$\log_{10} (k_z/k_{\rm b})$", fontsize=20)
save_fig(fig, "figure6.png")
save_fig(fig, f"figure6.png")
if __name__ == "__main__":
plt.show()
import sys
import matplotlib.cm
import matplotlib.pyplot as plt
......@@ -2,5 +1,4 @@
import matplotlib.cm
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
......@@ -5,8 +3,9 @@
import numpy as np
from math import degrees
from util import compute_kf_kb_ko_keta_kd, save_fig
from util import (
compute_kf_kb_ko_keta_kd,
save_fig,
)
from util_simuls_regimes import get_sim
......@@ -17,14 +16,7 @@
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
print(sys.argv)
letter = sys.argv[-1]
if letter not in "DLOWP":
letter = "L"
def plot_spectra(sim, ax, key="Ee"):
def plot_spectra(sim, ax, key="Ee", key_k="kh"):
p_oper = sim.params.oper
kmax = p_oper.coef_dealiasing * p_oper.nx * np.pi / p_oper.Lx
t_start, t_last = sim.output.print_stdout.get_times_start_last()
......@@ -37,6 +29,5 @@
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]
if key == "EA":
spectrum = data["A"]
......@@ -41,5 +32,10 @@
if key == "EA":
spectrum = data["A"]
elif key == "EK":
EKhd = data["Khd"]
EKz = data["Kz"]
Etoro = data["Khr"]
spectrum = EKhd + EKz + Etoro
elif key == "Epolo":
EKhd = data["Khd"]
EKz = data["Kz"]
......@@ -62,15 +58,40 @@
print(f"Don't know key: {key} \n")
exit
cs = ax.pcolormesh(
kh,
kz,
np.log10(spectrum),
cmap=cm,
vmin=-9.0,
vmax=-1.0,
shading="nearest",
)
th = np.linspace(0, np.pi / 2, 50)
if key_k == "kh":
for iz in range(int(np.floor(np.log(len(kz)) / np.log(1.2)))):
nz = int(np.floor(1.2**iz))
if kz[nz] <= 10 * kb and kz[nz] >= 0.1 * kb:
cs = ax.plot(
kh,
spectrum[nz, :],
color=cm(0.5 + 0.5 * np.log10((kz[nz] / kb))),
linestyle="-",
)
ks = np.array([0.1 * kf, kb])
ax.plot(ks, 2e-1 * ks ** (-2), "k-")
ax.text(np.sqrt(0.1 * kf * kb), 2e-3, r"$\propto k_h^{-2}$", fontsize=14)
ks = np.array([0.1 * kf, kb])
ax.plot(ks, 1e-10 * ks ** (1), "k-")
ax.text(np.sqrt(0.1 * kf * kb), 3e-10, r"$\propto k_h^{1}$", fontsize=14)
ax.set_xlim([delta_kh, kmax])
ax.set_xlabel(r"$k_h$", fontsize=20)
elif key_k == "kz":
for ih in range(int(np.floor(np.log(len(kh)) / np.log(1.2)))):
nh = int(np.floor(1.2**ih))
if kh[nh] <= 10 * kb and kh[nh] >= 0.1 * kb:
cs = ax.plot(
kz,
spectrum[:, nh],
color=cm(0.5 + 0.5 * np.log10((kh[nh] / kb))),
linestyle="-",
)
ks = np.array([0.3 * kf, ko])
ax.plot(ks, 5e-9 * ks ** (0), "k-")
ax.text(np.sqrt(0.3 * kf * ko), 7e-10, r"$\propto k_z^{0}$", fontsize=14)
ks = np.array([kb, ko])
ax.plot(ks, 1e4 * ks ** (-4), "k-")
ax.text(np.sqrt(kb * ko), 1e-4, r"$\propto k_z^{-4}$", fontsize=14)
ax.set_xlim([delta_kz, kmax])
ax.set_xlabel(r"$k_z$", fontsize=20)
......@@ -76,93 +97,3 @@
# 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
# 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,
)
# 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",
)
# Chi_d = 1
a = 1.0
xa = np.linspace(kh[1], kd, 50, endpoint=True)
ax.plot(
xa,
xa * np.sqrt((kd / xa) ** (4 / 3) - 1),
linestyle="dotted",
color="magenta",
)
# Chi_nu = 1 (Kolmogorov scale)
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g")
ax.set_ylim([1e-10, 1e-1])
ax.set_xticks(
......@@ -168,5 +99,5 @@
ax.set_xticks(
[1e1, kb, 1e2, ko, 1e3],
[r"$10^1$", r"$k_{\rm b}$", r"$10^2$", r"$k_{\rm O}$", r"$10^3$"],
[1e1, 1e2, 1e3, kb, ko],
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
fontsize=14,
)
......@@ -171,14 +102,6 @@
fontsize=14,
)
ax.set_yticks(
[1e1, kb, 1e2, ko, 1e3],
[r"$10^1$", r"$k_{\rm b}$", r"$10^2$", r"$k_{\rm O}$", r"$10^3$"],
fontsize=14,
)
ax.set_xlim([delta_kh, kmax])
ax.set_ylim([delta_kz, kmax])
ax.set_xscale("log")
ax.set_yscale("log")
# ax.grid(True)
......@@ -181,7 +104,12 @@
ax.set_xscale("log")
ax.set_yscale("log")
# ax.grid(True)
ax.axvline(kb, color="k", linestyle="dotted")
ax.axvline(ko, color="k", linestyle="dashed")
ax.axvline(kf, color="orange", linestyle="dashed")
ax.axvline(keta, color="g", linestyle="dashed")
return cs
......@@ -185,4 +113,7 @@
return cs
### Figure: Kinetic energy
nbax = 0
css = [None for i in range(4)]
fig, axes = plt.subplots(
......@@ -188,8 +119,8 @@
fig, axes = plt.subplots(
ncols=2, nrows=3, figsize=(10, 3 * 3 * 4.5 / 4), constrained_layout=True
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]
......@@ -190,9 +121,7 @@
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ax4 = axes[2, 0]
ax5 = axes[2, 1]
......@@ -198,14 +127,4 @@
axs = [ax0, ax1, ax2, ax3, ax4, ax5]
sim = get_sim(letter)
cs0 = plot_spectra(sim, ax0, key="Etoro")
cs2 = plot_spectra(sim, ax2, key="Epolo")
cs4 = plot_spectra(sim, ax4, key="EA")
sim = get_sim(letter, proj=True)
cs1 = plot_spectra(sim, ax1, key="Etoro")
cs3 = plot_spectra(sim, ax3, key="Epolo")
cs5 = plot_spectra(sim, ax5, key="EA")
axs = [ax0, ax1, ax2, ax3]
......@@ -210,8 +129,9 @@
for ax in [ax0, ax2, ax4]:
ax.set_ylabel(r"$k_z$", fontsize=20)
for ax in [ax4, ax5]:
ax.set_xlabel(r"$k_h$", fontsize=20)
for proj in [False, True]:
sim = get_sim("L", proj=proj)
css[nbax] = plot_spectra(sim, axs[nbax], key="EK", key_k="kh")
css[nbax + 2] = plot_spectra(sim, axs[nbax + 2], key="EK", key_k="kz")
nbax += 1
......@@ -216,11 +136,17 @@
for ax in [ax0, ax1, ax2, ax3]:
ax.set_xticklabels([])
for ax in [ax1, ax3, ax5]:
for ax in [ax0, ax2]:
ax.set_ylabel(r"$E_{\rm kin}(k_h, k_z)$", fontsize=20)
ax.set_yticks(
[1e-10, 1e-8, 1e-6, 1e-4, 1e-2],
[r"$10^{-10}$", r"$10^{-8}$", r"$10^{-6}$", r"$10^{-4}$", r"$10^{-2}$"],
fontsize=14,
)
for ax in [ax1, ax3]:
ax.set_yticklabels([])
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 (d)$", fontsize=20)
......@@ -221,10 +147,8 @@
ax.set_yticklabels([])
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 (d)$", fontsize=20)
ax4.set_title(r"$\rm (e)$", fontsize=20)
ax5.set_title(r"$\rm (f)$", fontsize=20)
......@@ -229,19 +153,7 @@
# 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)
norm = matplotlib.colors.Normalize(vmin=-1, vmax=1)
fig.tight_layout()
fig.subplots_adjust(right=0.85)
......@@ -244,13 +156,15 @@
fig.tight_layout()
fig.subplots_adjust(right=0.85)
cbar_ax0 = fig.add_axes([0.88, 0.675, 0.02, 0.25])
cbar0 = fig.colorbar(cs0, cax=cbar_ax0, cmap=cm, orientation="vertical")
cbar0.set_label(r"$\log_{10} E_{\rm toro}$", fontsize=20)
cbar_ax2 = fig.add_axes([0.88, 0.3725, 0.02, 0.25])
cbar2 = fig.colorbar(cs2, cax=cbar_ax2, cmap=cm, orientation="vertical")
cbar2.set_label(r"$\log_{10} E_{\rm polo}$", fontsize=20)
cbar_ax = fig.add_axes([0.88, 0.1, 0.02, 0.33])
cbar = fig.colorbar(
matplotlib.cm.ScalarMappable(norm=norm, cmap=cm),
cax=cbar_ax,
cmap=cm,
orientation="vertical",
)
cbar.set_ticks([-1, -0.5, 0.0, 0.5, 1])
cbar.set_ticklabels([r"$-1$", r"$-0.5$", r"$0.0$", r"$0.5$", r"$1$"], fontsize=14)
cbar.set_label(r"$\log_{10} (k_h/k_{\rm b})$", fontsize=20)
......@@ -256,14 +170,15 @@
cbar_ax4 = fig.add_axes([0.88, 0.07, 0.02, 0.25])
cbar4 = fig.colorbar(cs4, cax=cbar_ax4, cmap=cm, orientation="vertical")
cbar4.set_label(r"$\log_{10} E_{\rm pot}$", fontsize=20)
for cbar in [cbar0, cbar2, cbar4]:
cbar.set_ticks([-9, -7, -5, -3, -1])
cbar.set_ticklabels(
[r"$-9$", r"$-7$", r"$-5$", r"$-3$", r"$-1$"], fontsize=14
)
cbar_ax = fig.add_axes([0.88, 0.565, 0.02, 0.33])
cbar = fig.colorbar(
matplotlib.cm.ScalarMappable(norm=norm, cmap=cm),
cax=cbar_ax,
cmap=cm,
orientation="vertical",
)
cbar.set_ticks([-1, -0.5, 0.0, 0.5, 1])
cbar.set_ticklabels([r"$-1$", r"$-0.5$", r"$0.0$", r"$0.5$", r"$1$"], fontsize=14)
cbar.set_label(r"$\log_{10} (k_z/k_{\rm b})$", fontsize=20)
save_fig(fig, "figure7.png")
......@@ -266,6 +181,7 @@
save_fig(fig, "figure7.png")
if __name__ == "__main__":
plt.show()
import sys
import matplotlib.cm
import matplotlib.pyplot as plt
......@@ -1,2 +3,3 @@
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
......@@ -2,7 +5,4 @@
import numpy as np
from matplotlib import cm
import matplotlib.patches as patches
from util_simuls_regimes import get_sim
from math import degrees
......@@ -6,5 +6,7 @@
from math import degrees
from util import compute_kf_kb_ko_keta_kd, customize, save_fig
from util import compute_kf_kb_ko_keta_kd, save_fig
from util_simuls_regimes import get_sim
......@@ -10,2 +12,5 @@
cm = matplotlib.cm.get_cmap("inferno", 100)
# Latex
plt.rcParams["text.usetex"] = True
......@@ -11,3 +16,4 @@
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
......@@ -12,56 +18,5 @@
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)
print(sys.argv)
letter = sys.argv[-1]
......@@ -67,68 +22,5 @@
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,
)
if letter not in "DLOWP":
letter = "L"
......@@ -133,44 +25,16 @@
# 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
def plot_spectra(sim, ax, key="Ee"):
p_oper = sim.params.oper
kmax = p_oper.coef_dealiasing * p_oper.nx * np.pi / p_oper.Lx
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]
......@@ -176,24 +40,25 @@
# 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)
if key == "EA":
spectrum = data["A"]
elif key == "Epolo":
EKhd = data["Khd"]
EKz = data["Kz"]
spectrum = EKhd + EKz
elif key == "Ee":
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
Epolo = EKhd + EKz
spectrum = 2 * np.minimum(EA, Epolo)
elif key == "Ed":
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
Epolo = EKhd + EKz
spectrum = EA + Epolo - 2 * np.minimum(EA, Epolo)
elif key == "Etoro":
spectrum = data["Khr"]
else:
print(f"Don't know key: {key} \n")
exit
......@@ -199,25 +64,11 @@
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",
cs = ax.pcolormesh(
kh,
kz,
np.log10(spectrum),
cmap=cm,
vmin=-9.0,
vmax=-1.0,
shading="nearest",
)
......@@ -222,11 +73,4 @@
)
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")
th = np.linspace(0, np.pi / 2, 50)
......@@ -232,32 +76,9 @@
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
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
# Forcing
ax.add_patch(
patches.Arc(
......@@ -312,6 +133,100 @@
linewidth=1,
)
# 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",
)
# Chi_d = 1
a = 1.0
xa = np.linspace(kh[1], kd, 50, endpoint=True)
ax.plot(
xa,
xa * np.sqrt((kd / xa) ** (4 / 3) - 1),
linestyle="dotted",
color="magenta",
)
# Chi_nu = 1 (Kolmogorov scale)
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g")
ax.set_xticks(
[1e1, kb, 1e2, ko, 1e3],
[r"$10^1$", r"$k_{\rm b}$", r"$10^2$", r"$k_{\rm O}$", r"$10^3$"],
fontsize=14,
)
ax.set_yticks(
[1e1, kb, 1e2, ko, 1e3],
[r"$10^1$", r"$k_{\rm b}$", r"$10^2$", r"$k_{\rm O}$", r"$10^3$"],
fontsize=14,
)
ax.set_xlim([delta_kh, kmax])
ax.set_ylim([delta_kz, kmax])
ax.set_xscale("log")
ax.set_yscale("log")
# ax.grid(True)
return cs
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]
axs = [ax0, ax1, ax2, ax3, ax4, ax5]
sim = get_sim(letter)
cs0 = plot_spectra(sim, ax0, key="Etoro")
cs2 = plot_spectra(sim, ax2, key="Epolo")
cs4 = plot_spectra(sim, ax4, key="EA")
sim = get_sim(letter, proj=True)
cs1 = plot_spectra(sim, ax1, key="Etoro")
cs3 = plot_spectra(sim, ax3, key="Epolo")
cs5 = plot_spectra(sim, ax5, key="EA")
for ax in [ax0, ax2, ax4]:
ax.set_ylabel(r"$k_z$", fontsize=20)
for ax in [ax4, ax5]:
ax.set_xlabel(r"$k_h$", fontsize=20)
for ax in [ax0, ax1, ax2, ax3]:
ax.set_xticklabels([])
for ax in [ax1, ax3, ax5]:
ax.set_yticklabels([])
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 (d)$", fontsize=20)
ax4.set_title(r"$\rm (e)$", fontsize=20)
ax5.set_title(r"$\rm (f)$", fontsize=20)
# Add annotations for the lines and forcing region
ax1.text(
......@@ -324,8 +239,8 @@
)
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)
ax1.text(7e1, 5e1, r"$\gamma_{\mathbf{k}} < 1$", color="m", fontsize=14)
fig.tight_layout()
fig.subplots_adjust(right=0.85)
......@@ -328,13 +243,10 @@
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)
cbar_ax0 = fig.add_axes([0.88, 0.675, 0.02, 0.25])
cbar0 = fig.colorbar(cs0, cax=cbar_ax0, cmap=cm, orientation="vertical")
cbar0.set_label(r"$\log_{10} E_{\rm toro}$", fontsize=20)
......@@ -339,11 +251,10 @@
# 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)
cbar_ax2 = fig.add_axes([0.88, 0.3725, 0.02, 0.25])
cbar2 = fig.colorbar(cs2, cax=cbar_ax2, cmap=cm, orientation="vertical")
cbar2.set_label(r"$\log_{10} E_{\rm polo}$", fontsize=20)
cbar_ax4 = fig.add_axes([0.88, 0.07, 0.02, 0.25])
cbar4 = fig.colorbar(cs4, cax=cbar_ax4, cmap=cm, orientation="vertical")
cbar4.set_label(r"$\log_{10} E_{\rm pot}$", fontsize=20)
......@@ -349,5 +260,12 @@
save_fig(fig, f"figure8.png")
for cbar in [cbar0, cbar2, cbar4]:
cbar.set_ticks([-9, -7, -5, -3, -1])
cbar.set_ticklabels(
[r"$-9$", r"$-7$", r"$-5$", r"$-3$", r"$-1$"], fontsize=14
)
save_fig(fig, "figure8.png")
if __name__ == "__main__":
plt.show()
......@@ -12,6 +12,6 @@
fig, axes = plt.subplots(
ncols=2, nrows=1, figsize=(10, 1.2 * 3 * 4.5 / 4), constrained_layout=True
ncols=2, nrows=2, figsize=(10, 2 * 3 * 4.5 / 4), constrained_layout=True
)
......@@ -16,8 +16,9 @@
)
ax0 = axes[0]
ax1 = axes[1]
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)
......@@ -61,6 +62,8 @@
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"]
......@@ -79,12 +82,8 @@
T = TA + TK
levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True)
cs0 = ax0.contourf(
KH,
KZ,
K2A / (D + np.abs(TA) + np.abs(TK) + np.abs(K2A)),
cmap=cm.seismic,
levels=levels,
)
# (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)
......@@ -89,9 +88,11 @@
th = np.linspace(0, np.pi / 2, 50)
for ax in [ax0]:
ax.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-")
# Chi_d = 1
a = 1.0
xa = np.linspace(kh[1], kd, 50, endpoint=True)
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,
......@@ -96,4 +97,4 @@
ax.plot(
xa,
xa * np.sqrt((kd / xa) ** (4 / 3) - 1),
xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dotted",
......@@ -99,4 +100,12 @@
linestyle="dotted",
color="magenta",
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",
)
......@@ -101,6 +110,6 @@
)
for ax in [ax0]:
for ax in [ax0, ax2]:
# Chi_nu = 1 (Kolmogorov scale)
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g")
......@@ -165,6 +174,8 @@
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"]
......@@ -183,12 +194,8 @@
T = TA + TK
levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True)
cs1 = ax1.contourf(
KH,
KZ,
K2A / (D + np.abs(TA) + np.abs(TK) + np.abs(K2A)),
cmap=cm.seismic,
levels=levels,
)
# (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)
......@@ -193,9 +200,11 @@
th = np.linspace(0, np.pi / 2, 50)
for ax in [ax1]:
ax.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-")
# Chi_d = 1
a = 1.0
xa = np.linspace(kh[1], kd, 50, endpoint=True)
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,
......@@ -200,4 +209,4 @@
ax.plot(
xa,
xa * np.sqrt((kd / xa) ** (4 / 3) - 1),
xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dotted",
......@@ -203,4 +212,12 @@
linestyle="dotted",
color="magenta",
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",
)
......@@ -205,6 +222,6 @@
)
for ax in [ax1]:
for ax in [ax1, ax3]:
# Chi_nu = 1 (Kolmogorov scale)
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g")
......@@ -227,10 +244,13 @@
)
for ax in [ax1]:
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)
......@@ -231,7 +251,9 @@
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)
......@@ -237,5 +259,5 @@
for ax in [ax0, ax1]:
for ax in [ax0, ax1, ax2, ax3]:
# Forcing
ax.add_patch(
patches.Arc(
......@@ -290,6 +312,18 @@
linewidth=1,
)
fig.tight_layout()
# 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)
......@@ -294,3 +328,4 @@
fig.tight_layout()
fig.subplots_adjust(right=0.85)
......@@ -296,7 +331,17 @@
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.17, 0.02, 0.65])
cbar = fig.colorbar(cs1, cax=cbar_ax)
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
)
......@@ -299,8 +344,8 @@
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{B}}$", fontsize=20)
cbar.set_label(r"$\tilde{\mathcal{D}}$", fontsize=20)
save_fig(fig, f"figure9.png")
......
import numpy as np
import matplotlib.pyplot as plt
import scipy
from fluiddyn.io.redirect_stdout import stdout_redirected
from fluidsim.solvers.ns3d.strat.solver import Simul
from util import save_fig
params = Simul.create_default_params()
params.output.HAS_TO_SAVE = False
params.output.sub_directory = "tmp4fig"
params.forcing.tcrandom.time_correlation = 1.0
params.oper.nx = nx = 12
params.oper.ny = params.oper.nz = nx
params.init_fields.type = "noise"
params.init_fields.noise.length = params.oper.Lz / 2
params.init_fields.noise.velo_max = 1.0
params.forcing.nkmin_forcing = 2
params.forcing.nkmax_forcing = nx // 4
params.forcing.enable = True
params.forcing.type = "tcrandom"
params.forcing.normalized.type = "2nd_degree_eq"
# this actually produces a cleaner forcing...
# params.forcing.normalized.which_root = "positive"
params.forcing.key_forced = "vx_fft"
with stdout_redirected():
sim = Simul(params)
dt = 0.01
t_end = 50
freq_sample = 1 / dt
nt = int(t_end / dt)
fx_vs_time = np.empty([nx**3, nt])
times = []
# effect of the forcing on the state
epsilon = 2e-3
dissipation_coef = 1 - 1e-3
for idx in range(nt):
time = idx * dt
times.append(time)
sim.time_stepping.t = time
sim.forcing.compute()
forcing_spect = sim.forcing.get_forcing()
sim.state.state_spect = (
dissipation_coef * sim.state.state_spect + epsilon * forcing_spect
)
fx_fft = forcing_spect.get_var("vx_fft")
fx = sim.oper.ifft(fx_fft).flatten()
fx_vs_time[:, idx] = fx
times = np.array(times)
freqs, periodograms = scipy.signal.periodogram(
fx_vs_time, freq_sample, scaling="density"
)
periodogram = periodograms.mean(0)
periodogram[0] = np.nan
fig, (ax0, ax1) = plt.subplots(2)
indices = list(range(0, 80, 10))
tmin_plot = 20
tmax_plot = tmin_plot + 10
# tmax_plot = max(times)
# tmin_plot = tmax_plot - 10
index_tmax = abs(times - tmax_plot).argmin()
index_tmin = abs(times - tmin_plot).argmin()
times = times[index_tmin : index_tmax + 1]
fx_vs_time = fx_vs_time[indices, index_tmin : index_tmax + 1]
for ix in range(fx_vs_time.shape[0]):
ax0.plot(times, fx_vs_time[ix])
ax0.set_xlabel("time / correlation time")
ax0.set_title("forcing for few spatial points")
ax0.set_xlim((tmin_plot, tmax_plot))
low, high = ax0.get_ylim()
bound = max(abs(low), abs(high))
ax0.set_ylim(-bound, bound)
ax1.loglog(freqs, periodogram)
ax1.set_xlabel("frequency / correlation frequency")
ax1.set_ylabel("PDS forcing")
ax1.set_xlim(right=10)
ax1.set_ylim(bottom=1e-4)
fig.tight_layout()
save_fig(fig, "forcing_vs_time.png")
from curses import keyname
import matplotlib.pyplot as plt
import numpy as np
from util_dataframe import df, df_proj
from util import save_fig
plt.rcParams["text.usetex"] = True
fig, axes = plt.subplots(
ncols=2, nrows=1, figsize=(10, 4.5), constrained_layout=True
)
ax0 = axes[0]
ax1 = axes[1]
ax0.scatter(
df["Fh"],
df["Etoro"] / df["E"],
c=np.log10(df["R2"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
vmin=-1,
vmax=4,
)
ax0.scatter(
df_proj["Fh"],
df_proj["Etoro"] / df_proj["E"],
c=np.log10(df_proj["R2"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
vmin=-1,
vmax=4,
)
ax0.set_xlim([1e-3, 20])
ax0.set_xticks([1e-3, 1e-2, 1e-1, 1e0, 1e1])
ax0.set_xticklabels([1e-3, 1e-2, 1e-1, 1e0, 1e1], fontsize=14)
ax0.set_xscale("log")
ax0.set_ylim([0, 1])
ax0.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax0.set_yticklabels(
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1.0$"], fontsize=14
)
ax0.set_xlabel(r"$F_h$", fontsize=20)
ax0.set_ylabel(r"$E_{\rm toro}/E$", fontsize=20)
ax0.set_title(r"$(a)$", fontsize=20)
ax0.grid(True)
cs = ax1.scatter(
df["Fh"],
(df["Epolo"] - df["EA"]) / (df["Epolo"] + df["EA"]),
c=np.log10(df["R2"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
vmin=-1,
vmax=4,
label=r"Standard Navier-Stokes",
)
ax1.scatter(
df_proj["Fh"],
(df_proj["Epolo"] - df_proj["EA"]) / (df_proj["Epolo"] + df_proj["EA"]),
c=np.log10(df_proj["R2"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
label=r"Without vortical modes",
)
ax1.legend(loc="lower center", fontsize=14)
ax1.set_xlim([1e-3, 20])
ax1.set_xticks([1e-3, 1e-2, 1e-1, 1e0, 1e1])
ax1.set_xticklabels([1e-3, 1e-2, 1e-1, 1e0, 1e1], fontsize=14)
ax1.set_xscale("log")
ax1.set_ylim([-1, 1])
ax1.set_yticks([-1, -0.5, 0.0, 0.5, 1.0])
ax1.set_yticklabels(
[r"$-1$", r"$-0.5$", r"$0.0$", r"$0.5$", r"$1.0$"], fontsize=14
)
ax1.set_xlabel(r"$F_h$", fontsize=20)
ax1.set_ylabel(r"$\tilde{\mathcal{D}}$", fontsize=20)
ax1.set_title(r"$(b)$", fontsize=20)
ax1.grid(True)
# Emphasizing the points (N,Rb) = (40,20)
df_L = df[df["N"] == 40]
df_L = df_L[df_L["Rb"] == 20]
df_L_proj = df_proj[df_proj["N"] == 40]
df_L_proj = df_L_proj[df_L_proj["Rb"] == 20]
ax0.scatter(
df_L["Fh"],
df_L["Etoro"] / df_L["E"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax0.scatter(
df_L_proj["Fh"],
df_L_proj["Etoro"] / df_L_proj["E"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax1.scatter(
df_L["Fh"],
(df_L["Epolo"] - df_L["EA"]) / (df_L["Epolo"] + df_L["EA"]),
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax1.scatter(
df_L_proj["Fh"],
(df_L_proj["Epolo"] - df_L_proj["EA"])
/ (df_L_proj["Epolo"] + df_L_proj["EA"]),
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
fig.tight_layout()
fig.subplots_adjust(right=0.85, wspace=0.3)
cbar_ax = fig.add_axes([0.88, 0.16, 0.02, 0.75])
cbar = fig.colorbar(cs, cax=cbar_ax, orientation="vertical")
cbar.set_ticks([-1, 0, 1, 2, 3, 4])
cbar.set_ticklabels(
[r"$-1$", r"$0$", r"$1$", r"$2$", r"$3$", r"$4$"], fontsize=14
)
# cbar.set_label(r"$\log_{10}(\mathcal{R})$", fontsize = 12, rotation=0)
fig.text(0.85, 0.07, r"$\log_{10} \mathcal{R}$", fontsize=20)
save_fig(fig, f"fig_ratio_E_vs_Fh.png")
if __name__ == "__main__":
plt.show()
from curses import keyname
import matplotlib.pyplot as plt
import numpy as np
from util_dataframe import df, df_proj
from util import save_fig
plt.rcParams["text.usetex"] = True
fig, axes = plt.subplots(
ncols=2, nrows=1, figsize=(10, 4.5), constrained_layout=True
)
df = df[df["N"] == 20]
df_proj = df_proj[df_proj["N"] == 20]
ax0 = axes[0]
ax1 = axes[1]
ax0.scatter(
df["R2"],
df["Etoro"] / df["E"],
c=np.log10(df["Fh"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
vmin=-3,
vmax=0.5,
)
ax0.scatter(
df_proj["R2"],
df_proj["Etoro"] / df_proj["E"],
c=np.log10(df_proj["Fh"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
vmin=-3,
vmax=0.5,
)
ax0.set_xlim([1e-1, 1e4])
ax0.set_xticks([1e-1, 1e-0, 1e1, 1e3, 1e3, 1e4])
ax0.set_xticklabels([1e-1, 1e-0, 1e1, 1e3, 1e3, 1e4], fontsize=14)
ax0.set_xscale("log")
ax0.set_ylim([0, 1])
ax0.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax0.set_yticklabels(
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1.0$"], fontsize=14
)
ax0.set_xlabel(r"$\mathcal{R}$", fontsize=20)
ax0.set_ylabel(r"$E_{\rm toro}/E$", fontsize=20)
ax0.set_title(r"$(a)$", fontsize=20)
ax0.grid(True)
cs = ax1.scatter(
df["R2"],
(df["Epolo"] - df["EA"]) / (df["Epolo"] + df["EA"]),
c=np.log10(df["Fh"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
vmin=-3,
vmax=0.5,
label=r"Standard Navier-Stokes",
)
ax1.scatter(
df_proj["R2"],
(df_proj["Epolo"] - df_proj["EA"]) / (df_proj["Epolo"] + df_proj["EA"]),
c=np.log10(df_proj["Fh"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
label=r"Without vortical modes",
)
ax1.legend(loc="lower center", fontsize=14)
ax1.set_xlim([1e-1, 1e4])
ax1.set_xticks([1e-1, 1e-0, 1e1, 1e3, 1e3, 1e4])
ax1.set_xticklabels([1e-1, 1e-0, 1e1, 1e3, 1e3, 1e4], fontsize=14)
ax1.set_xscale("log")
ax1.set_ylim([-1, 1])
ax1.set_yticks([-1, -0.5, 0.0, 0.5, 1.0])
ax1.set_yticklabels(
[r"$-1$", r"$-0.5$", r"$0.0$", r"$0.5$", r"$1.0$"], fontsize=14
)
ax1.set_xlabel(r"$\mathcal{R}$", fontsize=20)
ax1.set_ylabel(r"$\tilde{\mathcal{D}}$", fontsize=20)
ax1.set_title(r"$(b)$", fontsize=20)
ax1.grid(True)
# Emphasizing the points (N,Rb) = (40,20)
df_L = df[df["N"] == 40]
df_L = df_L[df_L["Rb"] == 20]
df_L_proj = df_proj[df_proj["N"] == 40]
df_L_proj = df_L_proj[df_L_proj["Rb"] == 20]
ax0.scatter(
df_L["R2"],
df_L["Etoro"] / df_L["E"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax0.scatter(
df_L_proj["R2"],
df_L_proj["Etoro"] / df_L_proj["E"],
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax1.scatter(
df_L["R2"],
(df_L["Epolo"] - df_L["EA"]) / (df_L["Epolo"] + df_L["EA"]),
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
ax1.scatter(
df_L_proj["R2"],
(df_L_proj["Epolo"] - df_L_proj["EA"])
/ (df_L_proj["Epolo"] + df_L_proj["EA"]),
edgecolors="r",
facecolors="none",
s=60,
marker="s",
linewidths=2,
)
fig.tight_layout()
fig.subplots_adjust(right=0.85, wspace=0.3)
cbar_ax = fig.add_axes([0.88, 0.16, 0.02, 0.75])
cbar = fig.colorbar(cs, cax=cbar_ax, orientation="vertical")
cbar.set_ticks([-3, -2, -1, 0])
cbar.set_ticklabels([r"$-3$", r"$-2$", r"$-1$", r"$0$"], fontsize=14)
# cbar.set_label(r"$\log_{10}(\mathcal{R})$", fontsize = 12, rotation=0)
fig.text(0.85, 0.07, r"$\log_{10} F_h$", fontsize=20)
save_fig(fig, f"fig_ratio_E_vs_R.png")
if __name__ == "__main__":
plt.show()
......@@ -16,8 +16,6 @@
from util_simuls_regimes import get_sim
from util_dataframe import df, df_proj, df_ratio_one, df_proj_ratio_one
cm = matplotlib.cm.get_cmap("inferno", 100)
cmbin = matplotlib.cm.get_cmap("binary", 100)
......
from util_dataframe import df, df_proj
from util import formatters, tmp_dir
from util import tmp_dir
def formatter_R(v):
if v % 1 == 0 or v >= 100:
return f"{v:.0f}"
else:
return f"{v:.1f}"
def formatter_N(v):
if v % 1 == 0:
return f"{v:.0f}"
elif v < 10:
return f"{v:.2f}"
else:
return f"{v:.1f}"
formatters = {
"N": formatter_N,
"Rb": formatter_R,
"k_max*eta": lambda v: f"{v:.2f}",
"k_max*lambda": lambda v: f"{v:.2f}",
"epsK2/epsK": lambda v: f"{v:.2f}",
"Fh": lambda v: f"{v:.2e}",
"R2": formatter_R,
"R4": lambda v: f"{v:.2e}",
"Re_lambda": formatter_R,
"Re": formatter_R,
"Lh/Lz": lambda v: f"{v:.0f}",
}
# print(df.columns)
......@@ -4,12 +35,6 @@
# print(df.columns)
"""
df_all = pd.concat([df, df_proj], axis=0)
df_all = df_all.fillna(value="no")
print(df_all)
grouped = df_all.groupby("proj")
print(grouped)
"""
df["Lh/Lz"] = df["nx"] / df["nz"]
......@@ -14,8 +39,8 @@
columns = ["N", "Rb", "nx", "nz"]
header = [r"$N$", r"$\R_i= 1/ \nu N^2$", r"$n_h$", r"$n_z$"]
column_format = "|llll|"
columns = ["N", "Rb", "Lh/Lz"]
header = [r"$N$", r"$\R_i$", r"$L_h/L_z$"]
column_format = "|lll|"
df.to_latex(
buf=tmp_dir / "table_params.tex",
......@@ -30,6 +55,4 @@
)
column_format = "|" + 3 * "l" + "|"
columns = [
......@@ -35,4 +58,5 @@
columns = [
"nx",
"k_max*eta",
"Fh",
"R2",
......@@ -40,4 +64,5 @@
# "Uh2", "epsK", "Gamma", "I_velocity", "I_dissipation",
]
column_format = "|" + len(columns) * "l" + "|"
......@@ -43,10 +68,6 @@
header = [
r"$\kmax\eta$",
r"$F_h$",
r"$\mathcal{{R}}$"
# "$Re_i$", "$\R_4$", r"$\epsK_2/\epsK$", r"$\kmax\lambda$",
# r"$Re_\lambda$", "${U_h}^2$", r"$\epsK$", "$\Gamma$", "$I_{\rm kin}$", "$I_{\rm diss}$",
]
df.to_latex(
df_tmp = df[columns]
header = ["$n_h$", r"$\kmax\eta$", r"$F_h$", r"$\mathcal{{R}}$"]
df_tmp.to_latex(
buf=tmp_dir / "table_better_simuls.tex",
......@@ -52,5 +73,5 @@
buf=tmp_dir / "table_better_simuls.tex",
columns=columns,
# columns=columns,
formatters=formatters,
column_format=column_format,
index=False,
......@@ -60,12 +81,7 @@
label=None,
)
header = [
r"$\kmax\eta$",
r"$F_h$",
r"$\mathcal{{R}}$"
# "$Re_i$", "$\R_4$", r"$\epsK_2/\epsK$", r"$\kmax\lambda$",
# r"$Re_\lambda$", "${U_h}^2$", r"$\epsK$", "$\Gamma$", "$I_{\rm kin}$", "$I_{\rm diss}$",
]
df_proj.to_latex(
df_tmp = df_proj[columns]
df_tmp.to_latex(
buf=tmp_dir / "table_better_simuls_proj.tex",
......@@ -71,5 +87,5 @@
buf=tmp_dir / "table_better_simuls_proj.tex",
columns=columns,
# columns=columns,
formatters=formatters,
column_format=column_format,
index=False,
......
from util_dataframe_simuls_regimes import df
from util import formatters, tmp_dir
# fmt: off
columns = [
"letter", "regime",
"N", "Rb",
'nx', 'nz',
'k_max*eta', 'epsK2/epsK',
'Fh', 'R2',
# 'R4',
# 'Uh2', 'epsK',
# 'Gamma', 'I_velocity', 'I_dissipation',
]
header = [
"", "regime",
"$N$", "$\R_i$",
'$n_x$', '$n_z$',
r'$\kmax\eta$', r'$\epsK_2/\epsK$',
'$F_h$', r'$\R_2$',
# r'$\R_4$',
# '${U_h}^2$', r'$\epsK$',
# '$\Gamma$', '$I_{\rm kin}$', '$I_{\rm diss}$',
]
# fmt: on
df.to_latex(
buf=tmp_dir / "table_simuls_regimes.tex",
columns=columns,
formatters=formatters,
index=False,
header=header,
escape=False,
caption=(r"5 simulations representative of different regimes."),
label="table-simuls-regimes",
)
......@@ -15,9 +15,5 @@
path_base = os.environ["STRAT_TURB_POLO2022"]
path_base_proj = os.environ["STRAT_TURB_POLO_PROJ2022"]
path_base_ratio_one = os.environ["STRAT_WAVES2022"]
paths_all = sorted(Path(path_base).glob("simul_folders/ns3d*"))
paths_all_proj = sorted(Path(path_base_proj).glob("simul_folders/ns3d*"))
......@@ -22,8 +18,14 @@
paths_all = sorted(Path(path_base).glob("simul_folders/ns3d*"))
paths_all_proj = sorted(Path(path_base_proj).glob("simul_folders/ns3d*"))
paths_all_ratio_one = sorted(
Path(path_base_ratio_one).glob("simul_folders/ns3d*")
)
try:
path_base_ratio_one = os.environ["STRAT_WAVES2022"]
except KeyError:
print("warning: STRAT_WAVES2022 environment variable not set")
else:
paths_all_ratio_one = sorted(
Path(path_base_ratio_one).glob("simul_folders/ns3d*")
)
here = Path(__file__).absolute().parent
tmp_dir = here.parent / "tmp"
......@@ -156,6 +158,6 @@
result["name"] = sim.output.name_run
if nu_2 != 0.0:
if nu_2 != 0:
result["lambda"] = sqrt(U**2 * nu_2 / epsK)
result["Re_lambda"] = U * result["lambda"] / nu_2
......@@ -160,8 +162,3 @@
result["lambda"] = sqrt(U**2 * nu_2 / epsK)
result["Re_lambda"] = U * result["lambda"] / nu_2
result["Rb"] = float(sim.params.short_name_type_run.split("_Rb")[-1])
# else:
# result["lambda"] = INFINITY
# result["Re_lambda"] = INFINITY
# result["Rb"] = INFINITY
......@@ -167,3 +164,4 @@
result["Rb"] = float(sim.params.short_name_type_run.split("_Rb")[-1])
result["nx"] = sim.params.oper.nx
result["nz"] = sim.params.oper.nz
......@@ -168,7 +166,57 @@
result["nx"] = sim.params.oper.nx
result["nz"] = sim.params.oper.nz
result["proj"] = sim.params.projection
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_start + 2
data = sim.output.spectra.load1d_mean(tmin, verbose=False)
kz = data["kz"]
delta_kz = kz[1]
EKx = data["spectra_vx_kz"].sum() * delta_kz
EKy = data["spectra_vy_kz"].sum() * delta_kz
EKz = data["spectra_vz_kz"].sum() * delta_kz
result["EK"] = EKx + EKy + EKz
EKhr = data["spectra_Khr_kz"].sum() * delta_kz
EKhd = data["spectra_Khd_kz"].sum() * delta_kz
EKz = data["spectra_vz_kz"].sum() * delta_kz
result["Epolo"] = EKhd + EKz
result["Etoro"] = EKhr
# Get spatiotemporal spectra
path_run = Path(sim.output.path_run)
paths_spec = sorted(path_run.glob("spatiotemporal/periodogram_[0-9]*.h5"))
if not paths_spec:
return
path_spec = paths_spec[-1]
with h5py.File(path_spec, "r") as f:
kh = f["kh_spectra"][:]
kz = f["kz_spectra"][:]
omegas = f["omegas"][:]
EA = f["spectrum_A"][:]
Epolo = f["spectrum_K"][:] - f["spectrum_Khr"][:]
Eequi = 2 * np.minimum(EA, Epolo)
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 = sim.params.N * KH / K_NOZERO
delta = 0.1
E_waves = np.zeros((len(kz), len(kh)))
for index_omega, omega in enumerate(omegas):
weight = gaussian_weight(omega, omega_disp, delta * omega_disp)
E_waves += weight * Eequi[:, :, index_omega]
dk2_dom = kh[1] * kz[1] * omegas[1]
result["E_waves"] = np.sum(np.nan_to_num(E_waves)) * dk2_dom
result["E_waves_norm"] = np.sum(EA + Epolo) * dk2_dom
def get_customized_dataframe(paths):
df = get_dataframe_from_paths(
......@@ -171,6 +219,6 @@
def get_customized_dataframe(paths):
df = get_dataframe_from_paths(
paths, tmin="t_last-2", use_cache=1, customize=customize
paths, tmin="t_start+2", use_cache=1, customize=customize
)
......@@ -176,19 +224,13 @@
)
if "nu0.0" in paths[0].name:
columns_old = df.columns.tolist()
# fmt: off
first_columns = [
"N", "nx", "nz", "Fh", "R4", "Gamma",
"lx1", "lx2", "lz1", "lz2", "I_velocity", "I_dissipation"]
# fmt: on
else:
df["Re"] = df.Rb * df.N**2
columns_old = df.columns.tolist()
# fmt: off
first_columns = [
"N", "Rb", "Re", "nx", "nz", "Fh", "R2", "k_max*eta", "epsK2/epsK", "Gamma",
"lx1", "lx2", "lz1", "lz2", "I_velocity", "I_dissipation"]
# fmt: on
df["Re"] = df.Rb * df.N**2
columns_old = df.columns.tolist()
# fmt: off
first_columns = [
"N", "Rb", "Re", "nx", "nz", "Fh", "R2", "k_max*eta", "epsK2/epsK", "Gamma",
"lx1", "lx2", "lz1", "lz2", "I_velocity", "I_dissipation"]
# fmt: on
columns = first_columns.copy()
for key in columns_old:
......@@ -274,35 +316,6 @@
}
def formatter_R(v):
if v % 1 == 0 or v >= 100:
return f"{v:.0f}"
else:
return f"{v:.1f}"
def formatter_N(v):
if v % 1 == 0:
return f"{v:.0f}"
elif v < 10:
return f"{v:.2f}"
else:
return f"{v:.1f}"
formatters = {
"N": formatter_N,
"Rb": formatter_R,
"k_max*eta": lambda v: f"{v:.2f}",
"k_max*lambda": lambda v: f"{v:.2f}",
"epsK2/epsK": lambda v: f"{v:.2f}",
"Fh": lambda v: f"{v:.2e}",
"R2": formatter_R,
"R4": lambda v: f"{v:.2e}",
"Re_lambda": formatter_R,
"Re": formatter_R,
}
Fh_limit = 0.14
R2_limit = 10.0
......
import numpy as np
from fluidsim import load
from util import (
couples320,
get_customized_dataframe,
get_path_finer_resol,
compute_E_waves_vs_kh_kz,
)
from util import couples320, get_customized_dataframe, get_path_finer_resol
......@@ -11,6 +2,5 @@
# Contruct dataframe
def construct_df(proj=False, ratio_one=False):
paths = []
for N, Rb in sorted(couples320):
......@@ -32,33 +22,7 @@
df = get_customized_dataframe(paths)
if not ratio_one:
df["k_max*lambda"] = df["k_max"] * df["lambda"]
Etoro = []
Epolo = []
E = []
Ratio_E_waves = []
for path in paths:
sim = load(path, hide_stdout=True)
t_start, t_last = sim.output.print_stdout.get_times_start_last()
# nh = sim.params.oper.nx
tmin = t_last - 2.0
data = sim.output.spectra.load1d_mean(tmin)
# data = sim.output.spectra.loadkzkh_mean(tmin, key_to_load = "Khr")
# kx = data["kx"]
kz = data["kz"]
delta_kz = kz[1]
EKx_kz = data["spectra_vx_kz"] * delta_kz
EKy_kz = data["spectra_vy_kz"] * delta_kz
EKz_kz = data["spectra_vz_kz"] * delta_kz
EKhd_kz = data["spectra_Khd_kz"] * delta_kz
EKhr_kz = data["spectra_Khr_kz"] * delta_kz
EA_kz = data["spectra_A_kz"] * delta_kz
eK = np.sum(EKx_kz + EKy_kz + EKz_kz)
eA = np.sum(EA_kz)
e = eK + eA
epolo = np.sum(EKz_kz + EKhd_kz)
etoro = np.sum(EKhr_kz)
df["k_max*lambda"] = df["k_max"] * df["lambda"]
df["E"] = df["EA"] + df["EK"]
df["ratio Ewaves"] = df["E_waves"] / df["E_waves_norm"]
......@@ -64,36 +28,8 @@
E.append(e)
Epolo.append(epolo)
Etoro.append(etoro)
nh = sim.params.oper.nx
nz = sim.params.oper.nz
Gamma = nz / nh
path_spec = sorted(path.glob(f"spatiotemporal/periodogram_[0-9]*.h5"))
if len(path_spec) == 1:
(
E_waves_omega,
E_waves,
E_tot,
kh,
kz,
omegas,
) = compute_E_waves_vs_kh_kz(sim)
E_waves = np.nan_to_num(E_waves)
E_waves_spatiotemporal = np.sum(E_waves)
E_tot_spatiotemporal = np.sum(E_tot)
ratio_E_waves = E_waves_spatiotemporal / E_tot_spatiotemporal
else:
ratio_E_waves = np.nan
Ratio_E_waves.append(ratio_E_waves)
df["E"] = E
df["Epolo"] = Epolo
df["Etoro"] = Etoro
df["R_waves"] = Ratio_E_waves
df = df[df["k_max*eta"] > 0.4]
return df
df = construct_df(proj=False)
df_proj = construct_df(proj=True)
......@@ -94,10 +30,6 @@
return df
df = construct_df(proj=False)
df_proj = construct_df(proj=True)
df_ratio_one = construct_df(proj=False, ratio_one=True)
df_proj_ratio_one = construct_df(proj=True, ratio_one=True)