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.cm
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -2,4 +2,3 @@ ...@@ -2,4 +2,3 @@
import matplotlib.cm import matplotlib.cm
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np import numpy as np
...@@ -5,5 +4,4 @@ ...@@ -5,5 +4,4 @@
import numpy as np import numpy as np
from util_simuls_regimes import get_sim
from math import degrees
from util import ( from util import (
...@@ -8,4 +6,3 @@ ...@@ -8,4 +6,3 @@
from util import ( from util import (
compute_kf_kb_ko_keta_kd,
save_fig, save_fig,
...@@ -11,6 +8,4 @@ ...@@ -11,6 +8,4 @@
save_fig, save_fig,
compute_E_waves_vs_kh_kz, customize,
spectra_vs_khomega_slice,
spectra_vs_kzomega_slice,
) )
...@@ -15,7 +10,3 @@ ...@@ -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) cm = matplotlib.cm.get_cmap("inferno", 100)
...@@ -21,8 +12,7 @@ ...@@ -21,8 +12,7 @@
cm = matplotlib.cm.get_cmap("inferno", 100) cm = matplotlib.cm.get_cmap("inferno", 100)
cmbin = matplotlib.cm.get_cmap("binary", 100)
# Latex # Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
...@@ -23,85 +13,10 @@ ...@@ -23,85 +13,10 @@
# Latex # Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
print(sys.argv) def plot_spectra_vs_omegak(sim, ax, key="Ee"):
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):
N = sim.params.N N = sim.params.N
t_start, t_last = sim.output.print_stdout.get_times_start_last() t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0 tmin = t_last - 2.0
...@@ -105,7 +20,6 @@ ...@@ -105,7 +20,6 @@
N = sim.params.N N = sim.params.N
t_start, t_last = sim.output.print_stdout.get_times_start_last() t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0 tmin = t_last - 2.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin) mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
E_waves_omega, E_waves, E_tot, kh, kz, omegas = compute_E_waves_vs_kh_kz(sim) Uh = np.sqrt(mean_values["Uh2"])
ratio_waves = E_waves / E_tot
...@@ -111,44 +25,33 @@ ...@@ -111,44 +25,33 @@
ax.set_xlabel(r"$k_h$", fontsize=20) path = sim.params.path_run
ax.set_xlim([kh[1], max(kh)]) path_spec = sorted(path.glob(f"spatiotemporal/periodogram_[0-9]*.h5"))
ax.set_xticks([25, 50, 75, 100, 125]) assert len(path_spec) == 1, f"Not only 1 periodogram in {path} \n"
ax.set_xticklabels( path_spec = path_spec[0]
[r"$25$", r"$50$", r"$75$", r"$100$", r"$125$"], fontsize=14 with h5py.File(path_spec, "r") as f:
) kh = f["kh_spectra"][:]
ax.set_ylabel(r"$k_z$", fontsize=20) kz = f["kz_spectra"][:]
ax.set_ylim([kz[1], max(kz)]) omegas = f["omegas"][:]
ax.set_yticks([25, 50, 75, 100, 125]) Epot = f["spectrum_A"][:]
ax.set_yticklabels( EKz = f["spectrum_K"][:] - f["spectrum_Khd"][:] - f["spectrum_Khr"][:]
[r"$25$", r"$50$", r"$75$", r"$100$", r"$125$"], fontsize=14 Epolo = f["spectrum_Khd"][:] + EKz
) Etoro = f["spectrum_Khr"][:]
cs = ax.pcolormesh( delta_kh = kh[1]
kh, delta_kz = kz[1]
kz, if key == "Epot":
ratio_waves, spectrum = Epot
cmap=cmbin, elif key == "EK":
vmin=0, spectrum = Epolo + Etoro
vmax=1, elif key == "Epolo":
shading="nearest", spectrum = Epolo
) elif key == "Ee":
f = kh / np.tan(0.3) spectrum = 2 * np.minimum(Epot, Epolo)
ax.plot(kh, f, linestyle="-", linewidth=1, color="orange") elif key == "Ed":
f = kh / np.tan(0.15) spectrum = Epot + Epolo - 2 * np.minimum(Epot, Epolo)
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange") elif key == "Etoro":
f = kh / np.tan(0.075) spectrum = Etoro
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange") elif key == "Etot":
f = kh / np.tan(0.225) spectrum = Epolo + Etoro + Epot
ax.plot(kh, f, linestyle="--", linewidth=1, color="orange") else:
print(f"Don't know key: {key} \n")
# kb exit
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",
)
...@@ -154,20 +57,34 @@ ...@@ -154,20 +57,34 @@
delta_kh = kh[1] KH, KZ = np.meshgrid(kh, kz)
delta_kz = kz[1] K = (KH**2 + KZ**2) ** 0.5
kf_min = sim.params.forcing.nkmin_forcing * delta_kz K_NOZERO = K.copy()
kf_max = sim.params.forcing.nkmax_forcing * delta_kz K_NOZERO[K_NOZERO == 0] = 1e-16
angle = sim.params.forcing.tcrandom_anisotropic.angle omega_disp = N * KH / K_NOZERO
delta_angle = sim.params.forcing.tcrandom_anisotropic.delta_angle
# Forcing omega_disp_grid = np.linspace(0, N, 25, endpoint=True)
ax.add_patch(
patches.Arc( OMEGA_DISP_GRID, OMEGA_GRID = np.meshgrid(omega_disp_grid, omegas)
xy=(0, 0), SPECTRUM = np.zeros(OMEGA_DISP_GRID.shape)
width=2 * kf_max,
height=2 * kf_max, for iod in range(len(omega_disp_grid) - 1):
angle=0, for io in range(len(omegas)):
theta1=90.0 - degrees(angle) - 0.5 * degrees(delta_angle), omega_disp_min = omega_disp_grid[iod]
theta2=90.0 - degrees(angle) + 0.5 * degrees(delta_angle), omega_disp_max = omega_disp_grid[iod + 1]
linestyle="-", IOD = np.where(
color="orange", np.logical_and(
linewidth=1, 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 @@ ...@@ -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 @@ ...@@ -174,44 +97,8 @@
) )
ax.add_patch( ax.set_ylabel(r"$\omega / N$", fontsize=20)
patches.Arc( ax.set_ylim([0, 3])
xy=(0, 0), ax.set_yticks([0, 1, 2, 3])
width=2 * kf_min, ax.set_yticklabels([r"$0$", r"$1$", r"$2$", r"$3$"], fontsize=14)
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,
)
return cs return cs
...@@ -215,4 +102,6 @@ ...@@ -215,4 +102,6 @@
return cs 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( fig, axes = plt.subplots(
...@@ -218,4 +107,4 @@ ...@@ -218,4 +107,4 @@
fig, axes = plt.subplots( 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 @@ ...@@ -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 @@ ...@@ -224,7 +120,14 @@
delta = 0.1 sim = get_sim("S1", proj=False)
sim = get_sim(letter, proj=False) css[0] = plot_spectra_vs_omegak(sim, ax0, key="Etot")
# 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_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 @@ ...@@ -230,5 +133,7 @@
sim = get_sim(letter, proj=True) sim = get_sim("S3", proj=False)
# cs1 = plot_slices_omega_ki_waves(sim, ax1, key="Ee", key_k="kh", ik=3, delta=delta) css[4] = plot_spectra_vs_omegak(sim, ax4, key="Etot")
cs1 = plot_ratio_E_waves_vs_kh_kz(sim, ax1)
sim_proj = get_sim("S3", proj=True)
css[5] = plot_spectra_vs_omegak(sim_proj, ax5, key="Etot")
...@@ -234,6 +139,12 @@ ...@@ -234,6 +139,12 @@
ax1.set_yticklabels([]) for ax in [ax1, ax3, ax5]:
ax1.set_ylabel("") 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) 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) ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
...@@ -237,4 +148,8 @@ ...@@ -237,4 +148,8 @@
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", 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) 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 @@ ...@@ -240,3 +155,5 @@
for ax in [ax0, ax1, ax2, ax3, ax4, ax5]:
ax.plot([0, 1], [0, 1], "w--")
fig.tight_layout() fig.tight_layout()
...@@ -241,5 +158,9 @@ ...@@ -241,5 +158,9 @@
fig.tight_layout() 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 @@ ...@@ -244,11 +165,7 @@
cbar_ax2 = fig.add_axes([0.88, 0.17, 0.02, 0.65]) cbar.ax.set_ylabel(
cbar2 = fig.colorbar(cs1, cax=cbar_ax2, cmap=cmbin, orientation="vertical") r"$\log_{10} E(\sin \theta_{\boldsymbol{k}},\omega)$", fontsize=20
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
) )
...@@ -252,7 +169,9 @@ ...@@ -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__": if __name__ == "__main__":
plt.show() plt.show()
import sys import sys
import matplotlib.cm import matplotlib.cm
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np import numpy as np
...@@ -4,4 +5,4 @@ ...@@ -4,4 +5,4 @@
import numpy as np import numpy as np
from util import save_fig from math import degrees
...@@ -7,4 +8,12 @@ ...@@ -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) cm = matplotlib.cm.get_cmap("inferno", 100)
...@@ -9,8 +18,9 @@ ...@@ -9,8 +18,9 @@
cm = matplotlib.cm.get_cmap("inferno", 100) cm = matplotlib.cm.get_cmap("inferno", 100)
cmbin = matplotlib.cm.get_cmap("binary", 100)
# Latex # Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
...@@ -11,19 +21,23 @@ ...@@ -11,19 +21,23 @@
# Latex # Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
def plot_ratio_waves_vs_Fh(ax, data=df, datab=df_proj): print(sys.argv)
ax.scatter( letter = sys.argv[-1]
data["Fh"],
data["R_waves"], if letter not in "DLOWP":
c=np.log10(data["R2"]), letter = "L"
cmap="inferno",
edgecolors="k",
s=40, normalize = False
marker="o",
label=r"Standard Navier-Stokes", # def plot_spectra(sim, ax, key="Ee"):
vmin=-1,
vmax=4,
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 @@ ...@@ -29,13 +43,52 @@
) )
cs = ax.scatter(
datab["Fh"], if key_k == "kh":
datab["R_waves"], spectra_normalized = spectra_vs_khomega_slice(
c=np.log10(datab["R2"]), E_waves_omega, kh, omegas, ik, normalize=normalize
cmap="inferno", )
edgecolors="k", ax.set_xlabel(r"$k_h$", fontsize=20)
s=40, ax.set_xlim([kh[1], max(kh)])
marker="^", cs = ax.pcolormesh(
label=r"Without vortical modes", kh,
vmin=-1, omegas / N,
vmax=4, 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 @@ ...@@ -41,16 +94,8 @@
) )
ax.set_xlim([1e-3, 20]) ax.set_yticks([0, 1, 2, 3])
ax.set_xticks([1e-3, 1e-2, 1e-1, 1e0, 1e1]) ax.set_yticklabels([r"$0$", r"$1$", r"$2$", r"$3$"], fontsize=14)
ax.set_xticklabels([1e-3, 1e-2, 1e-1, 1e0, 1e1], fontsize=14) cs.cmap.set_bad("k")
ax.set_xscale("log") cs.cmap.set_under("k")
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 return cs
...@@ -54,15 +99,30 @@ ...@@ -54,15 +99,30 @@
return cs return cs
def plot_ratio_waves_vs_R(ax, data=df, datab=df_proj): def plot_ratio_E_waves_vs_kh_kz(sim, ax):
ax.scatter( N = sim.params.N
data["R2"], t_start, t_last = sim.output.print_stdout.get_times_start_last()
data["R_waves"], tmin = t_last - 2.0
c=np.log10(data["Fh"]), kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
cmap="inferno", E_waves_omega, E_waves, E_tot, kh, kz, omegas = compute_E_waves_vs_kh_kz(sim)
edgecolors="k", ratio_waves = E_waves / E_tot
s=40,
marker="o", ax.set_xlabel(r"$k_h$", fontsize=20)
label=r"Standard Navier-Stokes", ax.set_xlim([kh[1], max(kh)])
vmin=-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_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, vmax=1,
...@@ -68,2 +128,24 @@ ...@@ -68,2 +128,24 @@
vmax=1, 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 @@ ...@@ -69,13 +151,22 @@
) )
cs = ax.scatter(
datab["R2"], delta_kh = kh[1]
datab["R_waves"], delta_kz = kz[1]
c=np.log10(datab["Fh"]), kf_min = sim.params.forcing.nkmin_forcing * delta_kz
cmap="inferno", kf_max = sim.params.forcing.nkmax_forcing * delta_kz
edgecolors="k", angle = sim.params.forcing.tcrandom_anisotropic.angle
s=40, delta_angle = sim.params.forcing.tcrandom_anisotropic.delta_angle
marker="^", # Forcing
label=r"Without vortical modes", ax.add_patch(
vmin=-3, patches.Arc(
vmax=1, 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 @@ ...@@ -81,10 +172,14 @@
) )
ax.set_xlim([1e-1, 1e5]) ax.add_patch(
ax.set_xticks([1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5]) patches.Arc(
ax.set_xticklabels([1e-1, 1e0, 1e1, 1e2, 1e3, 1e4, 1e5], fontsize=14) xy=(0, 0),
ax.set_xscale("log") width=2 * kf_min,
ax.set_ylim([0, 1]) height=2 * kf_min,
ax.set_yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0]) angle=0,
ax.set_yticklabels( theta1=90.0 - degrees(angle) - 0.5 * degrees(delta_angle),
[r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1.0$"], fontsize=14 theta2=90.0 - degrees(angle) + 0.5 * degrees(delta_angle),
linestyle="-",
color="orange",
linewidth=1,
)
) )
...@@ -90,8 +185,32 @@ ...@@ -90,8 +185,32 @@
) )
ax.set_xlabel(r"$\mathcal{R}$", fontsize=20) ax.plot(
ax.set_ylabel(r"$\tilde{E}_{\rm wave}$", fontsize=20) [
# ax.grid(True) 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 return cs
fig, axes = plt.subplots( fig, axes = plt.subplots(
...@@ -94,6 +213,6 @@ ...@@ -94,6 +213,6 @@
return cs return cs
fig, axes = plt.subplots( 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 @@ ...@@ -99,7 +218,5 @@
) )
#
ax0 = axes[0] ax0 = axes[0]
ax1 = axes[1] ax1 = axes[1]
...@@ -102,7 +219,10 @@ ...@@ -102,7 +219,10 @@
ax0 = axes[0] ax0 = axes[0]
ax1 = axes[1] 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 @@ ...@@ -108,44 +228,12 @@
# Emphasizing the points (N,Rb) = (40,20) sim = get_sim(letter, proj=True)
df_L = df[df["N"] == 40] # cs1 = plot_slices_omega_ki_waves(sim, ax1, key="Ee", key_k="kh", ik=3, delta=delta)
df_L = df_L[df_L["Rb"] == 20] cs1 = plot_ratio_E_waves_vs_kh_kz(sim, ax1)
df_L_proj = df_proj[df_proj["N"] == 40]
df_L_proj = df_L_proj[df_L_proj["Rb"] == 20] ax1.set_yticklabels([])
ax0.scatter( ax1.set_ylabel("")
df_L["Fh"],
df_L["R_waves"], ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", fontsize=20)
edgecolors="r", ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
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,
)
...@@ -150,8 +238,6 @@ ...@@ -150,8 +238,6 @@
ax0.legend(loc="upper right", fontsize=14) fig.tight_layout()
fig.subplots_adjust(right=0.85)
ax0.set_title(r"$\rm (a)$", fontsize=20)
ax1.set_title(r"$\rm (b)$", fontsize=20)
...@@ -156,12 +242,9 @@ ...@@ -156,12 +242,9 @@
ax1.set_yticklabels("") cbar_ax2 = fig.add_axes([0.88, 0.17, 0.02, 0.65])
ax1.set_ylabel("") 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)
fig.tight_layout() cbar2.set_ticks([0, 0.2, 0.4, 0.6, 0.8, 1])
cbar2.set_ticklabels(
cbar = fig.colorbar(cs0, ax=ax0, orientation="vertical") [r"$0$", r"$0.2$", r"$0.4$", r"$0.6$", r"$0.8$", r"$1$"], fontsize=14
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
) )
...@@ -167,11 +250,4 @@ ...@@ -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") save_fig(fig, "figure15.png")
......
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.cm
...@@ -3,8 +4,5 @@ ...@@ -3,8 +4,5 @@
from matplotlib.path import Path from util import save_fig
from matplotlib.patches import PathPatch
from util import save_fig, Fh_limit, R2_limit
from util_dataframe import df, df_proj from util_dataframe import df, df_proj
...@@ -8,4 +6,7 @@ ...@@ -8,4 +6,7 @@
from util_dataframe import df, df_proj from util_dataframe import df, df_proj
cm = matplotlib.cm.get_cmap("inferno", 100)
# Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
...@@ -11,2 +12,4 @@ ...@@ -11,2 +12,4 @@
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
...@@ -12,6 +15,41 @@ ...@@ -12,6 +15,41 @@
fig, ax = plt.subplots( def plot_ratio_waves_vs_Fh(ax, data=df, datab=df_proj):
ncols=1, nrows=1, figsize=(7.5, 1.5 * 3 * 4.5 / 4), constrained_layout=True 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 @@ ...@@ -16,15 +54,45 @@
cs = ax.scatter( def plot_ratio_waves_vs_R(ax, data=df, datab=df_proj):
df["Fh"], ax.scatter(
df["R2"], data["R2"],
# c=np.log10(df["R_waves"] * (df["Epolo"] + df["EA"]) / (df["Epolo"] + df["EA"] + df["Etoro"])), data["ratio Ewaves"],
# cmap="inferno", c=np.log10(data["Fh"]),
color="k", cmap="inferno",
# edgecolors="k", edgecolors="k",
marker="o", s=40,
# vmin=-2, marker="o",
# vmax=0, label=r"Standard Navier-Stokes",
# s=50 * df_proj["I_dissipation"] * 4 / 3, # 4/3 because of projection vmin=-3,
label="Standard Navier-Stokes", 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 @@ ...@@ -30,2 +98,3 @@
) )
#
...@@ -31,10 +100,51 @@ ...@@ -31,10 +100,51 @@
cs = ax.scatter(
df_proj["Fh"], ax0 = axes[0]
df_proj["R2"], ax1 = axes[1]
color="k",
marker="^", cs0 = plot_ratio_waves_vs_Fh(ax0, data=df, datab=df_proj)
label="Without vortical modes", 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 @@ ...@@ -38,13 +148,5 @@
) )
# Le Reun et al. 2018 ax0.legend(loc="upper right", fontsize=14)
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
...@@ -50,7 +152,5 @@ ...@@ -50,7 +152,5 @@
Fh = EpsK / (Uh**2 * N) ax0.set_title(r"$\rm (a)$", fontsize=20)
R = Re * Fh**2 ax1.set_title(r"$\rm (b)$", fontsize=20)
cs = ax.scatter(Fh, R, color="k", marker="s", label="Le Reun et al. (2018)")
...@@ -55,49 +155,7 @@ ...@@ -55,49 +155,7 @@
# Table 2 of Rodda et al. 2022 ax1.set_yticklabels("")
Lh = 2.0 # 6.0 # m ax1.set_ylabel("")
nu = 1e-6 # m^2.s^-1
F = np.array( fig.tight_layout()
[
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
...@@ -103,20 +161,6 @@ ...@@ -103,20 +161,6 @@
Uxrms = np.array( cbar = fig.colorbar(cs0, ax=ax0, orientation="vertical")
[ cbar.set_ticks([-1, 0, 1, 2, 3, 4])
0.00077, cbar.set_ticklabels(
0.0047, [r"$-1$", r"$0$", r"$1$", r"$2$", r"$3$", r"$4$"], fontsize=14
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,
]
) )
...@@ -122,31 +166,4 @@ ...@@ -122,31 +166,4 @@
) )
Uyrms = np.array( fig.text(0.42, 0.07, r"$\log_{10} \mathcal{R}$", fontsize=20)
[
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)")
...@@ -151,115 +168,9 @@ ...@@ -151,115 +168,9 @@
ax.set_xlim([1e-3, 20]) cbar = fig.colorbar(cs1, ax=ax1, orientation="vertical")
ax.set_xscale("log") cbar.set_ticks([-3, -2, -1, 0, 1])
ax.set_ylim([1e-1, 1e5]) cbar.set_ticklabels([r"$-3$", r"$-2$", r"$-1$", r"$0$", r"$1$"], fontsize=14)
ax.set_yscale("log") fig.text(0.88, 0.07, r"$\log_{10} F_h$", fontsize=20)
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, "figure16.png") 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 @@ ...@@ -34,7 +34,7 @@
marker="o", marker="o",
vmin=0.0, vmin=0.0,
vmax=1.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 @@ ...@@ -80,7 +80,7 @@
ax1.scatter(df_L_proj["Fh"], df_L_proj["R2"], edgecolors="r", facecolors='none', s=60, marker="s", linewidths=2) 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(2e-3, 3e2, r"LAST", fontsize=14, alpha=0.5)
ax0.text(1.5e-1, 1e4, r"Weakly", fontsize=14, alpha=0.5) ax0.text(1.5e-1, 1e4, r"Weakly", fontsize=14, alpha=0.5)
......
...@@ -6,6 +6,10 @@ ...@@ -6,6 +6,10 @@
from util import save_fig 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 plt.rcParams["text.usetex"] = True
fig, axes = plt.subplots( 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 sys
import numpy as np import h5py
import matplotlib.cm
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -4,4 +5,6 @@ ...@@ -4,4 +5,6 @@
import matplotlib.pyplot as plt 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 util_simuls_regimes import get_sim
...@@ -6,8 +9,17 @@ ...@@ -6,8 +9,17 @@
from util_simuls_regimes import get_sim 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 # Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
...@@ -9,8 +21,11 @@ ...@@ -9,8 +21,11 @@
# Latex # Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
cm = matplotlib.cm.get_cmap("inferno", 100)
print(sys.argv) print(sys.argv)
letter = sys.argv[-1] letter = sys.argv[-1]
...@@ -18,24 +33,25 @@ ...@@ -18,24 +33,25 @@
letter = "L" letter = "L"
fig, axes = plt.subplots( def plot_omega_spectra(sim, ax):
ncols=2, nrows=2, figsize=(10, 1.2 * 2 * 3 * 4.5 / 4), constrained_layout=True path = sim.params.path_run
) t_start, t_last = times_start_last_from_path(path)
tmin = t_last - 2
ax0 = axes[0, 0] mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
ax1 = axes[0, 1] params = load_params_simul(path)
ax2 = axes[1, 0] proj = params.projection
ax3 = axes[1, 1] 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"
coef_compensate = 5 / 3 path_spec = path_spec[0]
with h5py.File(path_spec, "r") as f:
# Standard Navier-Stokes kh = f["kh_spectra"][:]
sim = get_sim(letter) kz = f["kz_spectra"][:]
t_start, t_last = sim.output.print_stdout.get_times_start_last() omegas = f["omegas"][:]
tmin = t_last - 2.0 EA = f["spectrum_A"][:]
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin) EKz = f["spectrum_K"][:] - f["spectrum_Khd"][:] - f["spectrum_Khr"][:]
data = sim.output.spectra.load_kzkh_mean(tmin) Epolo = f["spectrum_Khd"][:] + EKz
kh = data["kh_spectra"] Etoro = f["spectrum_Khr"][:]
kz = data["kz"] E = Epolo + Etoro + EA
Ee = 2 * np.minimum(EA, Epolo)
...@@ -41,10 +57,7 @@ ...@@ -41,10 +57,7 @@
data = sim.output.spectra.load_kzkh_mean( EA = np.sum(EA, axis=0)
tmin, key_to_load=["A", "Khd", "Kz", "Khr"] Epolo = np.sum(Epolo, axis=0)
) Etoro = np.sum(Etoro, axis=0)
kh = data["kh_spectra"] Ee = np.sum(Ee, axis=0)
kz = data["kz"] E = np.sum(E, axis=0)
delta_kh = kh[1]
delta_kz = kz[1]
KH, KZ = np.meshgrid(kh, kz)
...@@ -50,19 +63,7 @@ ...@@ -50,19 +63,7 @@
EA = data["A"] EA = np.sum(EA, axis=0)
EKhd = data["Khd"] Epolo = np.sum(Epolo, axis=0)
EKz = data["Kz"] Etoro = np.sum(Etoro, axis=0)
EKhr = data["Khr"] Ee = np.sum(Ee, axis=0)
E = np.sum(E, axis=0)
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
...@@ -68,38 +69,31 @@ ...@@ -68,38 +69,31 @@
ax0.plot( coef_compensate = 0
kh, cs = ax.plot(
Epolo_vs_kh * kh ** (coef_compensate), omegas / N,
"-g", EA * omegas**coef_compensate,
label=r"$E_{\rm polo}(k_h)$", color="b",
) label=r"$E_{\rm pot}(\omega)$",
ax0.plot( )
kh, cs = ax.plot(
Etoro_vs_kh * kh ** (coef_compensate), omegas / N,
"-r", Epolo * omegas**coef_compensate,
label=r"$E_{\rm toro}(k_h)$", color="g",
) label=r"$E_{\rm polo}(\omega)$",
ax0.plot( )
kh, if proj == None:
Epot_vs_kh * kh ** (coef_compensate), cs = ax.plot(
"-b", omegas / N,
label=r"$E_{\rm pot}(k_h)$", Etoro * omegas**coef_compensate,
) color="r",
ax0.plot( label=r"$E_{\rm toro}(\omega)$",
kz, )
Epolo_vs_kz * kz ** (coef_compensate), om = np.array([0.1 * N, N])
"--g", ax.plot(
label=r"$E_{\rm polo}(k_z)$", om / N,
) 1e-4 * (om / N) ** (-2 + coef_compensate),
ax0.plot( "--",
kz, color="gray",
Etoro_vs_kz * kz ** (coef_compensate), label=None,
"--r", )
label=r"$E_{\rm toro}(k_z)$", ax.text(0.5, 1e-3, r"$\omega^{-2}$", color="gray", fontsize=14)
)
ax0.plot(
kz,
Epot_vs_kz * kz ** (coef_compensate),
"--b",
label=r"$E_{\rm pot}(k_z)$",
)
...@@ -105,17 +99,40 @@ ...@@ -105,17 +99,40 @@
sim.output.spect_energy_budg.plot_fluxes( ax.plot(
tmin=tmin, key_k="kh", ax=ax2, plot_conversion=False om / N,
) 1e-5 * (om / N) ** (-3 / 2 + coef_compensate),
for ax in [ax0, ax2]: "-.",
ax.axvline(kb, color="k", linestyle="dotted") color="gray",
ax.axvline(ko, color="k", linestyle="dashed") label=None,
ax.axvline(kf, color="orange", linestyle="dashed") )
ax.set_xlim([kh[1], max(kh)]) ax.text(0.5, 5e-6, r"$\omega^{-3/2}$", color="gray", fontsize=14)
ax.set_xticks([1e1, 1e2, 1e3, kb, ko])
ax.set_xticklabels( # Forcing
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"], angle = sim.params.forcing.tcrandom_anisotropic.angle
fontsize=14, delta_angle = sim.params.forcing.tcrandom_anisotropic.delta_angle
) omega_fmin = N * np.sin(angle - 0.5 * delta_angle)
# sim.output.spect_energy_budg.plot_fluxes(tmin=tmin, key_k="kz", ax=ax4) 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 @@ ...@@ -120,12 +137,5 @@
# Without vortical modes sim = get_sim(letter)
sim = get_sim(letter, proj=True) sim_proj = 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"]
...@@ -131,28 +141,2 @@ ...@@ -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 @@ ...@@ -158,42 +142,7 @@
ax1.plot( # TODO: uncomment this assert
kh, # assert (
Epolo_vs_kh * kh ** (coef_compensate), # sim.params.oper.nx == sim_proj.params.oper.nx
"-g", # ), f"Not the same resolution for simulation Without vortical modes: {sim.params.oper.nx} vs {sim_proj.params.oper.nx}"
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)
...@@ -198,32 +147,6 @@ ...@@ -198,32 +147,6 @@
for ax in [ax0, ax1]: fig, axes = plt.subplots(
ax.set_xscale("log") ncols=2, nrows=1, figsize=(10, 1.2 * 3 * 4.5 / 4), constrained_layout=True
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")
...@@ -229,63 +152,5 @@ ...@@ -229,63 +152,5 @@
""" ax0 = axes[0]
ax.plot( ax1 = axes[1]
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)
...@@ -290,7 +155,5 @@ ...@@ -290,7 +155,5 @@
ax2.set_xlabel(r"$k_h$", fontsize=20) cs0 = plot_omega_spectra(sim, ax0)
ax3.set_xlabel(r"$k_h$", fontsize=20) cs1 = plot_omega_spectra(sim_proj, ax1)
# ax4.set_xlabel(r"$k_z$", fontsize=20)
# ax5.set_xlabel(r"$k_z$", fontsize=20)
...@@ -296,6 +159,7 @@ ...@@ -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) 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) ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
...@@ -299,10 +163,6 @@ ...@@ -299,10 +163,6 @@
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", 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) 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() fig.tight_layout()
...@@ -306,5 +166,7 @@ ...@@ -306,5 +166,7 @@
fig.tight_layout() fig.tight_layout()
save_fig(fig, f"figure5.png") save_fig(fig, f"figure5.png")
...@@ -309,4 +171,5 @@ ...@@ -309,4 +171,5 @@
save_fig(fig, f"figure5.png") save_fig(fig, f"figure5.png")
if __name__ == "__main__": if __name__ == "__main__":
plt.show() plt.show()
import matplotlib.cm import sys
import matplotlib.pyplot as plt
import numpy as np
...@@ -4,8 +2,6 @@ ...@@ -4,8 +2,6 @@
from util import ( import numpy as np
compute_kf_kb_ko_keta_kd, import matplotlib.pyplot as plt
save_fig,
)
from util_simuls_regimes import get_sim from util_simuls_regimes import get_sim
...@@ -9,9 +5,9 @@ ...@@ -9,9 +5,9 @@
from util_simuls_regimes import get_sim 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 # Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
...@@ -13,48 +9,8 @@ ...@@ -13,48 +9,8 @@
# Latex # Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
print(sys.argv)
def plot_spectra(sim, ax, key="Ee", key_k="kh"): letter = sys.argv[-1]
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
...@@ -60,56 +16,5 @@ ...@@ -60,56 +16,5 @@
if key_k == "kh": if letter not in "DLOWPU":
for iz in range(int(np.floor(np.log(len(kz)) / np.log(1.2)))): letter = "L"
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
...@@ -114,6 +19,3 @@ ...@@ -114,6 +19,3 @@
### Figure: Kinetic energy
nbax = 0
css = [None for i in range(4)]
fig, axes = plt.subplots( fig, axes = plt.subplots(
...@@ -119,5 +21,5 @@ ...@@ -119,5 +21,5 @@
fig, axes = plt.subplots( 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] ax0 = axes[0, 0]
...@@ -125,6 +27,95 @@ ...@@ -125,6 +27,95 @@
ax2 = axes[1, 0] ax2 = axes[1, 0]
ax3 = axes[1, 1] 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 @@ ...@@ -129,9 +120,80 @@
for proj in [False, True]: # Without vortical modes
sim = get_sim("L", proj=proj) sim = get_sim(letter, proj=True)
css[nbax] = plot_spectra(sim, axs[nbax], key="EK", key_k="kh") t_start, t_last = sim.output.print_stdout.get_times_start_last()
css[nbax + 2] = plot_spectra(sim, axs[nbax + 2], key="EK", key_k="kz") tmin = t_last - 2.0
nbax += 1 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 @@ ...@@ -136,12 +198,80 @@
for ax in [ax0, ax2]: for ax in [ax0, ax1]:
ax.set_ylabel(r"$E_{\rm kin}(k_h, k_z)$", fontsize=20) ax.set_xscale("log")
ax.set_yticks( ax.set_yscale("log")
[1e-10, 1e-8, 1e-6, 1e-4, 1e-2], ax.set_ylim([1e-3, 1e1])
[r"$10^{-10}$", r"$10^{-8}$", r"$10^{-6}$", r"$10^{-4}$", r"$10^{-2}$"], 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, fontsize=14,
) )
for ax in [ax1, ax3]: for ax in [ax1, ax3]:
ax.set_yticklabels([]) ax.set_yticklabels([])
...@@ -143,10 +273,31 @@ ...@@ -143,10 +273,31 @@
fontsize=14, fontsize=14,
) )
for ax in [ax1, ax3]: for ax in [ax1, ax3]:
ax.set_yticklabels([]) 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) 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) ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
ax2.set_title(r"$\rm (c)$", fontsize=20) ax2.set_title(r"$\rm (c)$", fontsize=20)
ax3.set_title(r"$\rm (d)$", fontsize=20) ax3.set_title(r"$\rm (d)$", fontsize=20)
...@@ -148,7 +299,8 @@ ...@@ -148,7 +299,8 @@
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", 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) ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
ax2.set_title(r"$\rm (c)$", fontsize=20) ax2.set_title(r"$\rm (c)$", fontsize=20)
ax3.set_title(r"$\rm (d)$", 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 @@ ...@@ -154,4 +306,3 @@
norm = matplotlib.colors.Normalize(vmin=-1, vmax=1)
fig.tight_layout() fig.tight_layout()
...@@ -156,32 +307,6 @@ ...@@ -156,32 +307,6 @@
fig.tight_layout() fig.tight_layout()
fig.subplots_adjust(right=0.85) save_fig(fig, f"figure6.png")
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")
if __name__ == "__main__": if __name__ == "__main__":
plt.show() plt.show()
import sys
import matplotlib.cm import matplotlib.cm
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -2,5 +1,4 @@ ...@@ -2,5 +1,4 @@
import matplotlib.cm import matplotlib.cm
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np import numpy as np
...@@ -5,8 +3,9 @@ ...@@ -5,8 +3,9 @@
import numpy as np import numpy as np
from math import degrees from util import (
compute_kf_kb_ko_keta_kd,
from util import compute_kf_kb_ko_keta_kd, save_fig save_fig,
)
from util_simuls_regimes import get_sim from util_simuls_regimes import get_sim
...@@ -17,14 +16,7 @@ ...@@ -17,14 +16,7 @@
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}" plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
print(sys.argv) def plot_spectra(sim, ax, key="Ee", key_k="kh"):
letter = sys.argv[-1]
if letter not in "DLOWP":
letter = "L"
def plot_spectra(sim, ax, key="Ee"):
p_oper = sim.params.oper p_oper = sim.params.oper
kmax = p_oper.coef_dealiasing * p_oper.nx * np.pi / p_oper.Lx 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() t_start, t_last = sim.output.print_stdout.get_times_start_last()
...@@ -37,6 +29,5 @@ ...@@ -37,6 +29,5 @@
kz = data["kz"] kz = data["kz"]
delta_kh = kh[1] delta_kh = kh[1]
delta_kz = kz[1] delta_kz = kz[1]
if key == "EA": if key == "EA":
spectrum = data["A"] spectrum = data["A"]
...@@ -41,5 +32,10 @@ ...@@ -41,5 +32,10 @@
if key == "EA": if key == "EA":
spectrum = data["A"] spectrum = data["A"]
elif key == "EK":
EKhd = data["Khd"]
EKz = data["Kz"]
Etoro = data["Khr"]
spectrum = EKhd + EKz + Etoro
elif key == "Epolo": elif key == "Epolo":
EKhd = data["Khd"] EKhd = data["Khd"]
EKz = data["Kz"] EKz = data["Kz"]
...@@ -62,15 +58,40 @@ ...@@ -62,15 +58,40 @@
print(f"Don't know key: {key} \n") print(f"Don't know key: {key} \n")
exit exit
cs = ax.pcolormesh( if key_k == "kh":
kh, for iz in range(int(np.floor(np.log(len(kz)) / np.log(1.2)))):
kz, nz = int(np.floor(1.2**iz))
np.log10(spectrum), if kz[nz] <= 10 * kb and kz[nz] >= 0.1 * kb:
cmap=cm, cs = ax.plot(
vmin=-9.0, kh,
vmax=-1.0, spectrum[nz, :],
shading="nearest", color=cm(0.5 + 0.5 * np.log10((kz[nz] / kb))),
) linestyle="-",
)
th = np.linspace(0, np.pi / 2, 50) 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 @@ ...@@ -76,93 +97,3 @@
# Forcing ax.set_ylim([1e-10, 1e-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,
)
)
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_xticks( ax.set_xticks(
...@@ -168,5 +99,5 @@ ...@@ -168,5 +99,5 @@
ax.set_xticks( ax.set_xticks(
[1e1, kb, 1e2, ko, 1e3], [1e1, 1e2, 1e3, kb, ko],
[r"$10^1$", r"$k_{\rm b}$", r"$10^2$", r"$k_{\rm O}$", r"$10^3$"], [r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"],
fontsize=14, fontsize=14,
) )
...@@ -171,14 +102,6 @@ ...@@ -171,14 +102,6 @@
fontsize=14, 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_xscale("log")
ax.set_yscale("log") ax.set_yscale("log")
# ax.grid(True) # ax.grid(True)
...@@ -181,7 +104,12 @@ ...@@ -181,7 +104,12 @@
ax.set_xscale("log") ax.set_xscale("log")
ax.set_yscale("log") ax.set_yscale("log")
# ax.grid(True) # 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 return cs
...@@ -185,4 +113,7 @@ ...@@ -185,4 +113,7 @@
return cs return cs
### Figure: Kinetic energy
nbax = 0
css = [None for i in range(4)]
fig, axes = plt.subplots( fig, axes = plt.subplots(
...@@ -188,8 +119,8 @@ ...@@ -188,8 +119,8 @@
fig, axes = plt.subplots( 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] ax0 = axes[0, 0]
ax1 = axes[0, 1] ax1 = axes[0, 1]
ax2 = axes[1, 0] ax2 = axes[1, 0]
ax3 = axes[1, 1] ax3 = axes[1, 1]
...@@ -190,9 +121,7 @@ ...@@ -190,9 +121,7 @@
) )
ax0 = axes[0, 0] ax0 = axes[0, 0]
ax1 = axes[0, 1] ax1 = axes[0, 1]
ax2 = axes[1, 0] ax2 = axes[1, 0]
ax3 = axes[1, 1] ax3 = axes[1, 1]
ax4 = axes[2, 0]
ax5 = axes[2, 1]
...@@ -198,14 +127,4 @@ ...@@ -198,14 +127,4 @@
axs = [ax0, ax1, ax2, ax3, ax4, ax5] axs = [ax0, ax1, ax2, ax3]
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")
...@@ -210,8 +129,9 @@ ...@@ -210,8 +129,9 @@
for ax in [ax0, ax2, ax4]: for proj in [False, True]:
ax.set_ylabel(r"$k_z$", fontsize=20) sim = get_sim("L", proj=proj)
for ax in [ax4, ax5]: css[nbax] = plot_spectra(sim, axs[nbax], key="EK", key_k="kh")
ax.set_xlabel(r"$k_h$", fontsize=20) css[nbax + 2] = plot_spectra(sim, axs[nbax + 2], key="EK", key_k="kz")
nbax += 1
...@@ -216,11 +136,17 @@ ...@@ -216,11 +136,17 @@
for ax in [ax0, ax1, ax2, ax3]: for ax in [ax0, ax2]:
ax.set_xticklabels([]) ax.set_ylabel(r"$E_{\rm kin}(k_h, k_z)$", fontsize=20)
for ax in [ax1, ax3, ax5]: 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([]) ax.set_yticklabels([])
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", 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) ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
ax2.set_title(r"$\rm (c)$", fontsize=20) ax2.set_title(r"$\rm (c)$", fontsize=20)
ax3.set_title(r"$\rm (d)$", fontsize=20) ax3.set_title(r"$\rm (d)$", fontsize=20)
...@@ -221,10 +147,8 @@ ...@@ -221,10 +147,8 @@
ax.set_yticklabels([]) ax.set_yticklabels([])
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", 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) ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
ax2.set_title(r"$\rm (c)$", fontsize=20) ax2.set_title(r"$\rm (c)$", fontsize=20)
ax3.set_title(r"$\rm (d)$", 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 @@ ...@@ -229,19 +153,7 @@
# Add annotations for the lines and forcing region norm = matplotlib.colors.Normalize(vmin=-1, vmax=1)
ax1.text(
4.5,
1.5e1,
r"Forcing",
color="orange",
rotation=0.7 * 180 / np.pi,
fontsize=10,
)
ax1.text(1e1, 1.5e1, r"$\chi_{\mathbf{k}} < 1/3$", color="cyan", fontsize=14)
ax1.text(1e2, 1.5e2, r"$\chi_{\mathbf{k}} < 3$", color="cyan", fontsize=14)
ax1.text(7e1, 5e1, r"$\gamma_{\mathbf{k}} < 1$", color="m", fontsize=14)
fig.tight_layout() fig.tight_layout()
fig.subplots_adjust(right=0.85) fig.subplots_adjust(right=0.85)
...@@ -244,13 +156,15 @@ ...@@ -244,13 +156,15 @@
fig.tight_layout() fig.tight_layout()
fig.subplots_adjust(right=0.85) fig.subplots_adjust(right=0.85)
cbar_ax0 = fig.add_axes([0.88, 0.675, 0.02, 0.25]) cbar_ax = fig.add_axes([0.88, 0.1, 0.02, 0.33])
cbar0 = fig.colorbar(cs0, cax=cbar_ax0, cmap=cm, orientation="vertical") cbar = fig.colorbar(
cbar0.set_label(r"$\log_{10} E_{\rm toro}$", fontsize=20) matplotlib.cm.ScalarMappable(norm=norm, cmap=cm),
cax=cbar_ax,
cmap=cm,
cbar_ax2 = fig.add_axes([0.88, 0.3725, 0.02, 0.25]) orientation="vertical",
cbar2 = fig.colorbar(cs2, cax=cbar_ax2, cmap=cm, orientation="vertical") )
cbar2.set_label(r"$\log_{10} E_{\rm polo}$", fontsize=20) 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 @@ ...@@ -256,14 +170,15 @@
cbar_ax4 = fig.add_axes([0.88, 0.07, 0.02, 0.25]) cbar_ax = fig.add_axes([0.88, 0.565, 0.02, 0.33])
cbar4 = fig.colorbar(cs4, cax=cbar_ax4, cmap=cm, orientation="vertical") cbar = fig.colorbar(
cbar4.set_label(r"$\log_{10} E_{\rm pot}$", fontsize=20) matplotlib.cm.ScalarMappable(norm=norm, cmap=cm),
cax=cbar_ax,
for cbar in [cbar0, cbar2, cbar4]: cmap=cm,
cbar.set_ticks([-9, -7, -5, -3, -1]) orientation="vertical",
cbar.set_ticklabels( )
[r"$-9$", r"$-7$", r"$-5$", r"$-3$", r"$-1$"], fontsize=14 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") save_fig(fig, "figure7.png")
...@@ -266,6 +181,7 @@ ...@@ -266,6 +181,7 @@
save_fig(fig, "figure7.png") save_fig(fig, "figure7.png")
if __name__ == "__main__": if __name__ == "__main__":
plt.show() plt.show()
import sys
import matplotlib.cm
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
...@@ -1,2 +3,3 @@ ...@@ -1,2 +3,3 @@
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np import numpy as np
...@@ -2,7 +5,4 @@ ...@@ -2,7 +5,4 @@
import numpy as np import numpy as np
from matplotlib import cm
import matplotlib.patches as patches
from util_simuls_regimes import get_sim
from math import degrees from math import degrees
...@@ -6,5 +6,7 @@ ...@@ -6,5 +6,7 @@
from math import degrees 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 @@ ...@@ -10,2 +12,5 @@
cm = matplotlib.cm.get_cmap("inferno", 100)
# Latex
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
...@@ -11,3 +16,4 @@ ...@@ -11,3 +16,4 @@
plt.rcParams["text.usetex"] = True plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
...@@ -12,56 +18,5 @@ ...@@ -12,56 +18,5 @@
fig, axes = plt.subplots( print(sys.argv)
ncols=2, nrows=2, figsize=(10, 2 * 3 * 4.5 / 4), constrained_layout=True letter = sys.argv[-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)
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)
...@@ -67,68 +22,5 @@ ...@@ -67,68 +22,5 @@
data = sim.output.spect_energy_budg.load_mean(tmin=tmin) if letter not in "DLOWP":
kh = data["kh"] letter = "L"
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,
)
...@@ -133,44 +25,16 @@ ...@@ -133,44 +25,16 @@
# Without vortical modes def plot_spectra(sim, ax, key="Ee"):
sim = get_sim("L", proj=True) p_oper = sim.params.oper
kmax = p_oper.coef_dealiasing * p_oper.nx * np.pi / p_oper.Lx
mean_values = sim.output.get_mean_values(tmin="t_last-2", customize=customize) t_start, t_last = sim.output.print_stdout.get_times_start_last()
R2 = mean_values["R2"] tmin = t_last - 2.0
Uh2 = mean_values["Uh2"] kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
epsK = mean_values["epsK"] data = sim.output.spectra.load_kzkh_mean(
Fh = mean_values["Fh"] tmin, key_to_load=["A", "Khd", "Kz", "Khr"]
proj = sim.params.projection )
t_start, t_last = sim.output.print_stdout.get_times_start_last() kh = data["kh_spectra"]
tmin = t_last - 2.0 kz = data["kz"]
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin) delta_kh = kh[1]
delta_kz = kz[1]
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
...@@ -176,24 +40,25 @@ ...@@ -176,24 +40,25 @@
# Etoro / E (kh, kz) if key == "EA":
cs = ax1.contourf(KH, KZ, Etoro / E, cmap=cm.binary, levels=levels) spectrum = data["A"]
elif key == "Epolo":
data = sim.output.spect_energy_budg.load_mean(tmin=tmin) EKhd = data["Khd"]
kh = data["kh"] EKz = data["Kz"]
kz = data["kz"] spectrum = EKhd + EKz
delta_kh = kh[1] elif key == "Ee":
delta_kz = kz[1] EA = data["A"]
KH, KZ = np.meshgrid(kh, kz) EKhd = data["Khd"]
DA = data["diss_A"] EKz = data["Kz"]
TA = data["transfer_A"] Epolo = EKhd + EKz
TK = data["transfer_Kh"] + data["transfer_Kz"] spectrum = 2 * np.minimum(EA, Epolo)
K2A = data["conv_K2A"] elif key == "Ed":
DKh = data["diss_Kh"] EA = data["A"]
DKz = data["diss_Kz"] EKhd = data["Khd"]
DK = DKh + DKz EKz = data["Kz"]
D = DA + DK Epolo = EKhd + EKz
T = TA + TK spectrum = EA + Epolo - 2 * np.minimum(EA, Epolo)
levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True) elif key == "Etoro":
spectrum = data["Khr"]
# (Epolo - Epot)/ (Epolo + Epot) (kh, kz) else:
cs2 = ax3.contourf(KH, KZ, Ed, cmap=cm.seismic, levels=levels) print(f"Don't know key: {key} \n")
exit
...@@ -199,25 +64,11 @@ ...@@ -199,25 +64,11 @@
cs = ax.pcolormesh(
th = np.linspace(0, np.pi / 2, 50) kh,
for ax in [ax1, ax3]: kz,
# kb np.log10(spectrum),
ax.plot(kb * np.sin(th), kb * np.cos(th), color="cyan", linestyle="dashed") cmap=cm,
vmin=-9.0,
# Chi_L = 1/3, 3 vmax=-1.0,
a = 1 / 3 shading="nearest",
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",
) )
...@@ -222,11 +73,4 @@ ...@@ -222,11 +73,4 @@
) )
for ax in [ax1, ax3]: th = np.linspace(0, np.pi / 2, 50)
# 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")
...@@ -232,32 +76,9 @@ ...@@ -232,32 +76,9 @@
ax.set_xlabel(r"$k_h$", fontsize=20) # Forcing
ax.set_xticks( kf_min = sim.params.forcing.nkmin_forcing * delta_kz
[1e1, 1e2, 1e3, kb, ko], kf_max = sim.params.forcing.nkmax_forcing * delta_kz
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_{\rm b}$", r"$k_{\rm O}$"], angle = sim.params.forcing.tcrandom_anisotropic.angle
fontsize=14, delta_angle = sim.params.forcing.tcrandom_anisotropic.delta_angle
)
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 # Forcing
ax.add_patch( ax.add_patch(
patches.Arc( patches.Arc(
...@@ -312,6 +133,100 @@ ...@@ -312,6 +133,100 @@
linewidth=1, 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 # Add annotations for the lines and forcing region
ax1.text( ax1.text(
...@@ -324,8 +239,8 @@ ...@@ -324,8 +239,8 @@
) )
ax1.text(1e1, 1.5e1, r"$\chi_{\mathbf{k}} < 1/3$", color="cyan", fontsize=14) 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(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.tight_layout()
fig.subplots_adjust(right=0.85) fig.subplots_adjust(right=0.85)
...@@ -328,13 +243,10 @@ ...@@ -328,13 +243,10 @@
fig.tight_layout() fig.tight_layout()
fig.subplots_adjust(right=0.85) 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_ax0 = fig.add_axes([0.88, 0.675, 0.02, 0.25])
cbar.set_ticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0]) cbar0 = fig.colorbar(cs0, cax=cbar_ax0, cmap=cm, orientation="vertical")
cbar.set_ticklabels( cbar0.set_label(r"$\log_{10} E_{\rm toro}$", fontsize=20)
[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)
...@@ -339,11 +251,10 @@ ...@@ -339,11 +251,10 @@
# fig.subplots_adjust(right=0.85) cbar_ax2 = fig.add_axes([0.88, 0.3725, 0.02, 0.25])
cbar_ax = fig.add_axes([0.88, 0.11, 0.02, 0.36]) cbar2 = fig.colorbar(cs2, cax=cbar_ax2, cmap=cm, orientation="vertical")
cbar = fig.colorbar(cs2, cax=cbar_ax) cbar2.set_label(r"$\log_{10} E_{\rm polo}$", fontsize=20)
cbar.set_ticks([-2 / 3, -1 / 3, 0.0, 1 / 3, 2 / 3])
cbar.set_ticklabels( cbar_ax4 = fig.add_axes([0.88, 0.07, 0.02, 0.25])
[r"$-2/3$", r"$-1/3$", r"$0$", r"$1/3$", r"$2/3$"], fontsize=14 cbar4 = fig.colorbar(cs4, cax=cbar_ax4, cmap=cm, orientation="vertical")
) cbar4.set_label(r"$\log_{10} E_{\rm pot}$", fontsize=20)
cbar.set_label(r"$\tilde{\mathcal{D}}$", fontsize=20)
...@@ -349,5 +260,12 @@ ...@@ -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__": if __name__ == "__main__":
plt.show() plt.show()
...@@ -12,6 +12,6 @@ ...@@ -12,6 +12,6 @@
fig, axes = plt.subplots( 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 @@ ...@@ -16,8 +16,9 @@
) )
ax0 = axes[0] ax0 = axes[0, 0]
ax1 = axes[1] ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
# Standard Navier-Stokes # Standard Navier-Stokes
sim = get_sim("L", proj=False) sim = get_sim("L", proj=False)
...@@ -61,6 +62,8 @@ ...@@ -61,6 +62,8 @@
K = np.sqrt(KH**2 + KZ**2) K = np.sqrt(KH**2 + KZ**2)
K[K == 0] = 1e-15 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) data = sim.output.spect_energy_budg.load_mean(tmin=tmin)
kh = data["kh"] kh = data["kh"]
...@@ -79,12 +82,8 @@ ...@@ -79,12 +82,8 @@
T = TA + TK T = TA + TK
levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True) levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True)
cs0 = ax0.contourf( # (Epolo - Epot)/ (Epolo + Epot) (kh, kz)
KH, cs2 = ax2.contourf(KH, KZ, Ed, cmap=cm.seismic, levels=levels)
KZ,
K2A / (D + np.abs(TA) + np.abs(TK) + np.abs(K2A)),
cmap=cm.seismic,
levels=levels,
)
th = np.linspace(0, np.pi / 2, 50) th = np.linspace(0, np.pi / 2, 50)
...@@ -89,9 +88,11 @@ ...@@ -89,9 +88,11 @@
th = np.linspace(0, np.pi / 2, 50) th = np.linspace(0, np.pi / 2, 50)
for ax in [ax0]: for ax in [ax0, ax2]:
ax.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-") # kb
# Chi_d = 1 ax.plot(kb * np.sin(th), kb * np.cos(th), color="cyan", linestyle="dashed")
a = 1.0
xa = np.linspace(kh[1], kd, 50, endpoint=True) # Chi_L = 1/3, 3
a = 1 / 3
xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
ax.plot( ax.plot(
xa, xa,
...@@ -96,4 +97,4 @@ ...@@ -96,4 +97,4 @@
ax.plot( ax.plot(
xa, xa,
xa * np.sqrt((kd / xa) ** (4 / 3) - 1), xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dotted", linestyle="dotted",
...@@ -99,4 +100,12 @@ ...@@ -99,4 +100,12 @@
linestyle="dotted", 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 @@ ...@@ -101,6 +110,6 @@
) )
for ax in [ax0]: for ax in [ax0, ax2]:
# Chi_nu = 1 (Kolmogorov scale) # Chi_nu = 1 (Kolmogorov scale)
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g") ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g")
...@@ -165,6 +174,8 @@ ...@@ -165,6 +174,8 @@
K = np.sqrt(KH**2 + KZ**2) K = np.sqrt(KH**2 + KZ**2)
K[K == 0] = 1e-15 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) data = sim.output.spect_energy_budg.load_mean(tmin=tmin)
kh = data["kh"] kh = data["kh"]
...@@ -183,12 +194,8 @@ ...@@ -183,12 +194,8 @@
T = TA + TK T = TA + TK
levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True) levels = np.linspace(-2 / 3, 2 / 3, 51, endpoint=True)
cs1 = ax1.contourf( # (Epolo - Epot)/ (Epolo + Epot) (kh, kz)
KH, cs2 = ax3.contourf(KH, KZ, Ed, cmap=cm.seismic, levels=levels)
KZ,
K2A / (D + np.abs(TA) + np.abs(TK) + np.abs(K2A)),
cmap=cm.seismic,
levels=levels,
)
th = np.linspace(0, np.pi / 2, 50) th = np.linspace(0, np.pi / 2, 50)
...@@ -193,9 +200,11 @@ ...@@ -193,9 +200,11 @@
th = np.linspace(0, np.pi / 2, 50) th = np.linspace(0, np.pi / 2, 50)
for ax in [ax1]: for ax in [ax1, ax3]:
ax.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-") # kb
# Chi_d = 1 ax.plot(kb * np.sin(th), kb * np.cos(th), color="cyan", linestyle="dashed")
a = 1.0
xa = np.linspace(kh[1], kd, 50, endpoint=True) # Chi_L = 1/3, 3
a = 1 / 3
xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
ax.plot( ax.plot(
xa, xa,
...@@ -200,4 +209,4 @@ ...@@ -200,4 +209,4 @@
ax.plot( ax.plot(
xa, xa,
xa * np.sqrt((kd / xa) ** (4 / 3) - 1), xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dotted", linestyle="dotted",
...@@ -203,4 +212,12 @@ ...@@ -203,4 +212,12 @@
linestyle="dotted", 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 @@ ...@@ -205,6 +222,6 @@
) )
for ax in [ax1]: for ax in [ax1, ax3]:
# Chi_nu = 1 (Kolmogorov scale) # Chi_nu = 1 (Kolmogorov scale)
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g") ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="dotted", color="g")
...@@ -227,10 +244,13 @@ ...@@ -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_ylabel(None)
ax.set_yticks([]) ax.set_yticks([])
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", 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) ax1.set_title(r"Without vortical modes" + "\n" + r"$\rm (b)$", fontsize=20)
...@@ -231,7 +251,9 @@ ...@@ -231,7 +251,9 @@
ax.set_ylabel(None) ax.set_ylabel(None)
ax.set_yticks([]) ax.set_yticks([])
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$\rm (a)$", 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) 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 @@ ...@@ -237,5 +259,5 @@
for ax in [ax0, ax1]: for ax in [ax0, ax1, ax2, ax3]:
# Forcing # Forcing
ax.add_patch( ax.add_patch(
patches.Arc( patches.Arc(
...@@ -290,6 +312,18 @@ ...@@ -290,6 +312,18 @@
linewidth=1, 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 @@ ...@@ -294,3 +328,4 @@
fig.tight_layout()
fig.subplots_adjust(right=0.85) fig.subplots_adjust(right=0.85)
...@@ -296,7 +331,17 @@ ...@@ -296,7 +331,17 @@
fig.subplots_adjust(right=0.85) fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.17, 0.02, 0.65]) cbar_ax = fig.add_axes([0.88, 0.535, 0.02, 0.36])
cbar = fig.colorbar(cs1, cax=cbar_ax) 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_ticks([-2 / 3, -1 / 3, 0.0, 1 / 3, 2 / 3])
cbar.set_ticklabels( cbar.set_ticklabels(
[r"$-2/3$", r"$-1/3$", r"$0$", r"$1/3$", r"$2/3$"], fontsize=14 [r"$-2/3$", r"$-1/3$", r"$0$", r"$1/3$", r"$2/3$"], fontsize=14
) )
...@@ -299,8 +344,8 @@ ...@@ -299,8 +344,8 @@
cbar.set_ticks([-2 / 3, -1 / 3, 0.0, 1 / 3, 2 / 3]) cbar.set_ticks([-2 / 3, -1 / 3, 0.0, 1 / 3, 2 / 3])
cbar.set_ticklabels( cbar.set_ticklabels(
[r"$-2/3$", r"$-1/3$", r"$0$", r"$1/3$", r"$2/3$"], fontsize=14 [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") 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 @@ ...@@ -16,8 +16,6 @@
from util_simuls_regimes import get_sim 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) cm = matplotlib.cm.get_cmap("inferno", 100)
cmbin = matplotlib.cm.get_cmap("binary", 100) cmbin = matplotlib.cm.get_cmap("binary", 100)
......
from util_dataframe import df, df_proj 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) # print(df.columns)
...@@ -4,12 +35,6 @@ ...@@ -4,12 +35,6 @@
# print(df.columns) # print(df.columns)
""" df["Lh/Lz"] = df["nx"] / df["nz"]
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)
"""
...@@ -14,8 +39,8 @@ ...@@ -14,8 +39,8 @@
columns = ["N", "Rb", "nx", "nz"] columns = ["N", "Rb", "Lh/Lz"]
header = [r"$N$", r"$\R_i= 1/ \nu N^2$", r"$n_h$", r"$n_z$"] header = [r"$N$", r"$\R_i$", r"$L_h/L_z$"]
column_format = "|llll|" column_format = "|lll|"
df.to_latex( df.to_latex(
buf=tmp_dir / "table_params.tex", buf=tmp_dir / "table_params.tex",
...@@ -30,6 +55,4 @@ ...@@ -30,6 +55,4 @@
) )
column_format = "|" + 3 * "l" + "|"
columns = [ columns = [
...@@ -35,4 +58,5 @@ ...@@ -35,4 +58,5 @@
columns = [ columns = [
"nx",
"k_max*eta", "k_max*eta",
"Fh", "Fh",
"R2", "R2",
...@@ -40,4 +64,5 @@ ...@@ -40,4 +64,5 @@
# "Uh2", "epsK", "Gamma", "I_velocity", "I_dissipation", # "Uh2", "epsK", "Gamma", "I_velocity", "I_dissipation",
] ]
column_format = "|" + len(columns) * "l" + "|"
...@@ -43,10 +68,6 @@ ...@@ -43,10 +68,6 @@
header = [ df_tmp = df[columns]
r"$\kmax\eta$",
r"$F_h$", header = ["$n_h$", r"$\kmax\eta$", r"$F_h$", r"$\mathcal{{R}}$"]
r"$\mathcal{{R}}$" df_tmp.to_latex(
# "$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(
buf=tmp_dir / "table_better_simuls.tex", buf=tmp_dir / "table_better_simuls.tex",
...@@ -52,5 +73,5 @@ ...@@ -52,5 +73,5 @@
buf=tmp_dir / "table_better_simuls.tex", buf=tmp_dir / "table_better_simuls.tex",
columns=columns, # columns=columns,
formatters=formatters, formatters=formatters,
column_format=column_format, column_format=column_format,
index=False, index=False,
...@@ -60,12 +81,7 @@ ...@@ -60,12 +81,7 @@
label=None, label=None,
) )
header = [
r"$\kmax\eta$", df_tmp = df_proj[columns]
r"$F_h$", df_tmp.to_latex(
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(
buf=tmp_dir / "table_better_simuls_proj.tex", buf=tmp_dir / "table_better_simuls_proj.tex",
...@@ -71,5 +87,5 @@ ...@@ -71,5 +87,5 @@
buf=tmp_dir / "table_better_simuls_proj.tex", buf=tmp_dir / "table_better_simuls_proj.tex",
columns=columns, # columns=columns,
formatters=formatters, formatters=formatters,
column_format=column_format, column_format=column_format,
index=False, 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 @@ ...@@ -15,9 +15,5 @@
path_base = os.environ["STRAT_TURB_POLO2022"] path_base = os.environ["STRAT_TURB_POLO2022"]
path_base_proj = os.environ["STRAT_TURB_POLO_PROJ2022"] 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 = sorted(Path(path_base).glob("simul_folders/ns3d*"))
paths_all_proj = sorted(Path(path_base_proj).glob("simul_folders/ns3d*")) paths_all_proj = sorted(Path(path_base_proj).glob("simul_folders/ns3d*"))
...@@ -22,8 +18,14 @@ ...@@ -22,8 +18,14 @@
paths_all = sorted(Path(path_base).glob("simul_folders/ns3d*")) paths_all = sorted(Path(path_base).glob("simul_folders/ns3d*"))
paths_all_proj = sorted(Path(path_base_proj).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 here = Path(__file__).absolute().parent
tmp_dir = here.parent / "tmp" tmp_dir = here.parent / "tmp"
...@@ -156,6 +158,6 @@ ...@@ -156,6 +158,6 @@
result["name"] = sim.output.name_run result["name"] = sim.output.name_run
if nu_2 != 0.0: if nu_2 != 0:
result["lambda"] = sqrt(U**2 * nu_2 / epsK) result["lambda"] = sqrt(U**2 * nu_2 / epsK)
result["Re_lambda"] = U * result["lambda"] / nu_2 result["Re_lambda"] = U * result["lambda"] / nu_2
...@@ -160,8 +162,3 @@ ...@@ -160,8 +162,3 @@
result["lambda"] = sqrt(U**2 * nu_2 / epsK) result["lambda"] = sqrt(U**2 * nu_2 / epsK)
result["Re_lambda"] = U * result["lambda"] / nu_2 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 @@ ...@@ -167,3 +164,4 @@
result["Rb"] = float(sim.params.short_name_type_run.split("_Rb")[-1])
result["nx"] = sim.params.oper.nx result["nx"] = sim.params.oper.nx
result["nz"] = sim.params.oper.nz result["nz"] = sim.params.oper.nz
...@@ -168,7 +166,57 @@ ...@@ -168,7 +166,57 @@
result["nx"] = sim.params.oper.nx result["nx"] = sim.params.oper.nx
result["nz"] = sim.params.oper.nz 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): def get_customized_dataframe(paths):
df = get_dataframe_from_paths( df = get_dataframe_from_paths(
...@@ -171,6 +219,6 @@ ...@@ -171,6 +219,6 @@
def get_customized_dataframe(paths): def get_customized_dataframe(paths):
df = get_dataframe_from_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 @@ ...@@ -176,19 +224,13 @@
) )
if "nu0.0" in paths[0].name: df["Re"] = df.Rb * df.N**2
columns_old = df.columns.tolist()
# fmt: off columns_old = df.columns.tolist()
first_columns = [
"N", "nx", "nz", "Fh", "R4", "Gamma", # fmt: off
"lx1", "lx2", "lz1", "lz2", "I_velocity", "I_dissipation"] first_columns = [
# fmt: on "N", "Rb", "Re", "nx", "nz", "Fh", "R2", "k_max*eta", "epsK2/epsK", "Gamma",
else: "lx1", "lx2", "lz1", "lz2", "I_velocity", "I_dissipation"]
df["Re"] = df.Rb * df.N**2 # fmt: on
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() columns = first_columns.copy()
for key in columns_old: for key in columns_old:
...@@ -274,35 +316,6 @@ ...@@ -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 Fh_limit = 0.14
R2_limit = 10.0 R2_limit = 10.0
......
import numpy as np from util import couples320, get_customized_dataframe, get_path_finer_resol
from fluidsim import load
from util import (
couples320,
get_customized_dataframe,
get_path_finer_resol,
compute_E_waves_vs_kh_kz,
)
...@@ -11,6 +2,5 @@ ...@@ -11,6 +2,5 @@
# Contruct dataframe
def construct_df(proj=False, ratio_one=False): def construct_df(proj=False, ratio_one=False):
paths = [] paths = []
for N, Rb in sorted(couples320): for N, Rb in sorted(couples320):
...@@ -32,33 +22,7 @@ ...@@ -32,33 +22,7 @@
df = get_customized_dataframe(paths) df = get_customized_dataframe(paths)
if not ratio_one: df["k_max*lambda"] = df["k_max"] * df["lambda"]
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"]
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)
...@@ -64,36 +28,8 @@ ...@@ -64,36 +28,8 @@
E.append(e) df = df[df["k_max*eta"] > 0.4]
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
return df return df
df = construct_df(proj=False) df = construct_df(proj=False)
df_proj = construct_df(proj=True) df_proj = construct_df(proj=True)
...@@ -94,10 +30,6 @@ ...@@ -94,10 +30,6 @@
return df return df
df = construct_df(proj=False) df = construct_df(proj=False)
df_proj = construct_df(proj=True) 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)