Skip to content
Snippets Groups Projects

Compare revisions

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

Source

Select target project
No results found

Target

Select target project
  • fluiddyn/fluiddyn_papers
1 result
Show changes
Showing
with 2464 additions and 0 deletions
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm, ticker
from util_simuls_regimes import get_sim
from fluidsim import load
from util import (
compute_kf_kb_ko_keta_kd,
customize,
get_path_finer_resol,
save_fig,
)
plt.rcParams["text.usetex"] = True
###
N = 20
Rb = None
proj = False
ratio_one = True
###
path = get_path_finer_resol(N, Rb, proj, ratio_one)
sim = load(path)
fig, axes = plt.subplots(
ncols=2, nrows=2, figsize=(10, 2 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
mean_values = sim.output.get_mean_values(tmin="t_last-2", customize=customize)
Uh2 = mean_values["Uh2"]
epsK = mean_values["epsK"]
Fh = mean_values["Fh"]
proj = sim.params.projection
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
data = sim.output.spectra.load_kzkh_mean(
tmin, key_to_load=["A", "Khd", "Kz", "Khr"]
)
kh = data["kh_spectra"]
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]
KH, KZ = np.meshgrid(kh, kz)
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
EKhr = data["Khr"]
Epolo = EKhd + EKz
Etoro = EKhr
E = Epolo + Etoro + EA
Ee = 2 * np.minimum(Epolo, EA)
Ed = Epolo + EA - Ee
E[E == 0] = 1e-15
levels = np.linspace(0, 1, 51, endpoint=True)
K = np.sqrt(KH**2 + KZ**2)
K[K == 0] = 1e-15
# EA / E (kh, kz)
cs = ax0.contourf(KH, KZ, Ee / E, cmap=cm.binary, levels=levels)
ax0.set_title(r"$E_{equi}/E$", fontsize=16)
# Epolo / E (kh, kz)
cs = ax1.contourf(KH, KZ, Ed / E, cmap=cm.binary, levels=levels)
ax1.set_title(r"$E_{diff}/E$", fontsize=16)
# Etoro / E (kh, kz)
cs = ax2.contourf(KH, KZ, Etoro / E, cmap=cm.binary, levels=levels)
ax2.set_title(r"$E_{toro}/E$", fontsize=16)
data = sim.output.spect_energy_budg.load_mean(tmin=tmin)
kh = data["kh"]
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]
KH, KZ = np.meshgrid(kh, kz)
DA = data["diss_A"]
TA = data["transfer_A"]
TK = data["transfer_Kh"] + data["transfer_Kz"]
K2A = data["conv_K2A"]
DKh = data["diss_Kh"]
DKz = data["diss_Kz"]
DK = DKh + DKz
D = DA + DK
T = TA + TK
levels = np.linspace(-1, 1, 51, endpoint=True)
cs2 = ax3.contourf(
KH,
KZ,
K2A / (D + np.abs(TA) + np.abs(TK) + np.abs(K2A)),
cmap=cm.seismic,
levels=levels,
)
ax3.set_title(r"$\tilde{\mathcal{B}}$", fontsize=16)
# ax0.legend()
th = np.linspace(0, np.pi / 2, 50)
for ax in [ax0, ax1, ax2, ax3]:
ax.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-")
ax.plot(kb * np.sin(th), kb * np.cos(th), color="k", linestyle="dotted")
ax.plot(ko * np.sin(th), ko * np.cos(th), "k--")
a = 3
xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
ax.plot(
xa,
xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dashed",
color="gray",
)
a = 1 / 3
xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
ax.plot(
xa,
xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dotted",
color="gray",
)
ax.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-")
ax.plot(kf * np.sin(th), kf * np.cos(th), linestyle="--", color="orange")
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="--", color="g")
ax.set_xlim([kh[1], 0.8 * max(kh)])
ax.set_ylim([kz[1], 0.8 * max(kh)])
ax.set_xscale("log")
ax.set_yscale("log")
for ax in [ax0, ax2]:
ax.set_ylabel(r"$k_z$", fontsize=16)
ax.set_yticks(
[1e1, 1e2, 1e3, kb, ko],
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_b$", r"$k_O$"],
)
for ax in [ax2, ax3]:
ax.set_xlabel(r"$k_h$", fontsize=16)
ax.set_xticks(
[1e1, 1e2, 1e3, kb, ko],
[r"$10^1$", r"$10^2$", r"$10^3$", r"$k_b$", r"$k_O$"],
)
for ax in [ax0, ax1]:
ax.set_xticks([])
for ax in [ax1, ax3]:
ax.set_yticks([])
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.53, 0.02, 0.35])
cbar = fig.colorbar(cs, cax=cbar_ax)
cbar.set_ticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
# fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.11, 0.02, 0.35])
cbar = fig.colorbar(cs2, cax=cbar_ax)
cbar.set_ticks([-1.0, -0.5, 0.0, 0.5, 1.0])
# fig.tight_layout()
save_fig(
fig,
f"fig_seb_regimes_one_couple_proj{proj}_ratio_one{ratio_one}_N{N}_Rb{Rb}.png",
)
if __name__ == "__main__":
plt.show()
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm, ticker
from util_simuls_regimes import get_sim
import matplotlib.patches as patches
from math import degrees
from fluidsim import load
from util import (
compute_kf_kb_ko_keta_kd,
customize,
get_path_finer_resol,
get_paths,
save_fig,
)
plt.rcParams["text.usetex"] = True
def plot_B(sim, ax):
mean_values = sim.output.get_mean_values(tmin="t_last-2", customize=customize)
R4 = mean_values["R4"]
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 - 3.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
data = sim.output.spect_energy_budg.load_mean(tmin=tmin)
# data = sim.output.spect_energy_budg.load_mean(tmin=t_last-1.2, tmax=t_last-1.0)
kh = data["kh"]
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)
DA = data["diss_A"]
TA = data["transfer_A"]
TK = data["transfer_Kh"] + data["transfer_Kz"]
K2A = data["conv_K2A"]
DKh = data["diss_Kh"]
DKz = data["diss_Kz"]
DK = DKh + DKz
D = DA + DK
T = TA + TK
levels = np.linspace(-1, 1, 51, endpoint=True)
cs = ax.contourf(
KH,
KZ,
K2A / (D + np.abs(TA) + np.abs(TK) + np.abs(K2A)),
cmap=cm.seismic,
levels=levels,
)
th = np.linspace(0, np.pi / 2, 100, endpoint=True)
ax.plot(kb * np.sin(th), kb * np.cos(th), color="k", linestyle="dashed")
ax.plot(
kf_max * np.sin(th),
kf_max * np.cos(th),
color="orange",
linestyle="dashed",
)
a = 1 / 3
xa = np.linspace(delta_kh, a**1.5 * ko, 50, endpoint=True)
ax.plot(
xa,
xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dashed",
color="c",
)
# 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="dashed",
color="magenta",
)
# ax.plot([delta_kh, max(kh)], [delta_kh, max(kh)], "k-")
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="--", color="g")
# 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=2,
)
)
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=2,
)
)
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=2,
)
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=2,
)
ax.set_xlim([delta_kh, 2 * max(kh) / 3])
ax.set_ylim([delta_kh, 2 * max(kh) / 3])
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel(r"$k_h$", fontsize=20)
ax.set_xticks(
[1e1, 1e2],
[r"$10^1$", r"$10^2$"],
fontsize=14,
)
ax.set_ylabel(r"$k_z$", fontsize=20)
ax.set_yticks(
[1e1, 1e2],
[r"$10^1$", r"$10^2$"],
fontsize=14,
)
return cs
Ns = [10, 20, 80]
nbax = 0
css = [None for i in range(6)]
fig, axes = plt.subplots(
ncols=2, nrows=3, figsize=(10, 1.2 * 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]
for N in Ns:
for proj in [False, True]:
# path = get_path_finer_resol(N=N, Rb=None, proj=proj, ratio_one=True)
path = get_paths(
N, Rb=None, nh=640, proj=proj, ratio_one=True, reverse=False
)
sim = load(path)
css[nbax] = plot_B(sim, axs[nbax])
nbax += 1
for ax in [ax0, ax1, ax2, ax3]:
ax.set_xlabel("")
ax.set_xticks([])
for ax in [ax1, ax3, ax5]:
ax.set_ylabel("")
ax.set_yticks([])
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$N=10$", fontsize=20)
ax1.set_title(r"Without vortical modes" + "\n" + r"$N=10$", fontsize=20)
ax2.set_title(r"$N=20$", fontsize=20)
ax3.set_title(r"$N=20$", fontsize=20)
ax4.set_title(r"$N=80$", fontsize=20)
ax5.set_title(r"$N=80$", fontsize=20)
ax0.text(2e0, 1e2, r"$k_{\rm b}$", color="k", fontsize=14)
ax0.text(4e0, 1e2, r"$k_{\eta}$", color="g", fontsize=14)
ax0.text(8e0, 1e2, r"$\chi_{\boldmath{k}} = 1/3$", color="c", fontsize=14)
ax0.text(3.2e1, 1e2, r"$\gamma_{\boldmath{k}} = 1$", color="m", fontsize=14)
fig.tight_layout()
fig.subplots_adjust(bottom=0.15)
# cbar_ax = fig.add_axes([0.88, 0.17, 0.02, 0.65])
cbar_ax = fig.add_axes([0.15, 0.08, 0.8, 0.01])
cbar = fig.colorbar(css[3], cax=cbar_ax, cmap=cm, orientation="horizontal")
cbar.set_label(
r"$\tilde{\mathcal{B}}$",
fontsize=20,
)
cbar.set_ticks([-1.0, -0.5, 0.0, 0.5, 1.0])
cbar.set_ticklabels([r"$-1$", r"$-0.5$", r"$0$", r"$0.5$", r"$1$"], fontsize=14)
fig.subplots_adjust(bottom=0.15, wspace=0.05, hspace=0.16)
save_fig(fig, f"fig_seb_transition_ratio_one.png")
if __name__ == "__main__":
plt.show()
import h5py
import matplotlib.cm
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.collections import LineCollection
from fluidsim import load
from fluidsim.util import load_params_simul, times_start_last_from_path
from util import (
compute_kf_kb_ko_keta_kd,
customize,
formatters,
get_path_finer_resol,
get_paths,
pos_closest_value,
save_fig,
spectral_fit_k_sint,
spectral_fit_kh_kz,
)
cm = matplotlib.cm.get_cmap("inferno", 100)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
def plot_spectra(sim, ax, key="Ee"):
p_oper = sim.params.oper
kmax = p_oper.coef_dealiasing * p_oper.nx * np.pi / p_oper.Lx
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
data = sim.output.spectra.load_kzkh_mean(
tmin, key_to_load=["A", "Khd", "Kz", "Khr"]
)
kh = data["kh_spectra"]
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]
kh, kz = np.meshgrid(kh, kz)
if key == "EA":
spectrum = data["A"]
elif key == "Epolo":
EKhd = data["Khd"]
EKz = data["Kz"]
spectrum = EKhd + EKz
elif key == "Ee":
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
Epolo = EKhd + EKz
spectrum = 2 * np.minimum(EA, Epolo)
elif key == "Ed":
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
Epolo = EKhd + EKz
spectrum = EA + Epolo - 2 * np.minimum(EA, Epolo)
elif key == "Etoro":
spectrum = data["Khr"]
else:
print(f"Don't know key: {key} \n")
exit
k = np.sqrt(kh**2 + kz**2)
k_nozero = k.copy()
k_nozero[k_nozero == 0] = 1e-16
sint = kh / k_nozero
theta = np.arcsin(sint)
"""
kminreg = 1*kf
kmaxreg = 0.5*kb
thetaminreg = np.pi/4
thetamaxreg = np.pi/2
c, ah, az = spectral_fit_kh_kz(spectrum=spectrum, kh=kh, kz=kz, kmin=kminreg, kmax=kmaxreg, thetamin=thetaminreg, thetamax=thetamaxreg)
spectrum_fit = c * (kh**ah) * (kz**az)
"""
"""
c, ak, at = spectral_fit_k_sint(spectrum=spectrum, k=k, sint=sint, kmin=kminreg, kmax=kmaxreg, thetamin=thetaminreg, thetamax=thetamaxreg)
spectrum_fit = c * (k**ak) * (sint**at)
"""
for nz in range(len(kz[:, 0])):
if kz[nz, 0] <= kb and kz[nz, 0] > 0 * kf:
cs = ax.plot(
kh[nz, :],
spectrum[nz, :] * kh[nz, :] ** (5 / 3),
color=cm(kz[nz, 0] / kb),
linestyle="-",
)
"""
spec = spectrum[nz,:]
khs = kh[nz,:]
spec = spec[khs >= kz[nz,0]]
khs = khs[khs >= kz[nz,0]]
cs = ax.plot(khs, spec*khs**2, color=cm(kz[nz,0]/kb), linestyle="-")
"""
"""
k = k.flatten()
kh = kh.flatten()
kz = kz.flatten()
spectrum = spectrum.flatten()
spectrum_fit = spectrum_fit.flatten()
theta = theta.flatten()
for n in range(len(k)):
if k[n] <= kmaxreg and k[n] >= kminreg and theta[n] >= thetaminreg and theta[n] <= thetamaxreg:
cs = ax.scatter(kh[n], spectrum[n]/spectrum_fit[n], color=cm(kz[n]/kb), linestyle="-")
ax.axvline(kminreg, color='m', linestyle='dashed')
ax.axvline(kmaxreg, color='m', linestyle='dotted')
"""
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")
ks = np.array([1.5 * kf, 0.5 * keta])
ax.plot(ks, 10 * ks ** (-2.69 + 5 / 3), "k-")
ax.text(15, 1, r"$\propto k_h^{-2.69}$", fontsize=16)
ax.set_xlim([delta_kh, kmax])
ax.set_ylim([1e-3, 1e1])
ax.set_xscale("log")
ax.set_yscale("log")
ax.grid(True)
return cs
Ns = [10, 20, 80]
nbax = 0
css = [None for i in range(6)]
fig, axes = plt.subplots(
ncols=2, nrows=3, figsize=(10, 3 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ax4 = axes[2, 0]
ax5 = axes[2, 1]
axs = [ax0, ax1, ax2, ax3, ax4, ax5]
for N in Ns:
for proj in [False, True]:
path = get_path_finer_resol(N=N, Rb=None, proj=proj, ratio_one=True)
sim = load(path, hide_stdout=True)
css[nbax] = plot_spectra(sim, axs[nbax], key="Ee")
nbax += 1
for ax in [ax0, ax2, ax4]:
ax.set_ylabel(r"$E_{equi}(k_h, k_z) \times k_h^{5/3}$", fontsize=16)
for ax in [ax4, ax5]:
ax.set_xlabel(r"$k_h$", fontsize=16)
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"$(a)$", fontsize=16)
ax1.set_title(r"Without vortical modes" + "\n" + r"$(b)$", fontsize=16)
ax2.set_title(r"$(c)$", fontsize=16)
ax3.set_title(r"$(d)$", fontsize=16)
ax4.set_title(r"$(e)$", fontsize=16)
ax5.set_title(r"$(f)$", fontsize=16)
norm = matplotlib.colors.Normalize(vmin=0, vmax=1)
fig.tight_layout()
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.28, 0.02, 0.35])
cbar = fig.colorbar(
matplotlib.cm.ScalarMappable(norm=norm, cmap=cm),
cax=cbar_ax,
cmap=cm,
orientation="vertical",
)
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$"], fontsize=14
)
cbar.set_label(r"$k_z/k_b$", fontsize=16)
save_fig(fig, "fig_spectra_kh_kz_ratio_one.png")
if __name__ == "__main__":
plt.show()
import glob
import os
import sys
from pathlib import Path
import h5py
import matplotlib.cm
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.collections import LineCollection
from fluidsim import load
from fluidsim.util import load_params_simul, times_start_last_from_path
from util import (
compute_kf_kb_ko_keta_kd,
customize,
formatters,
get_path_finer_resol,
get_paths,
pos_closest_value,
save_fig,
spectral_fit_k_sint,
spectral_fit_kh_kz,
)
cm = matplotlib.cm.get_cmap("inferno", 100)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
def plot_spectra(sim, ax, key="Ee"):
p_oper = sim.params.oper
kmax = p_oper.coef_dealiasing * p_oper.nx * np.pi / p_oper.Lx
t_start, t_last = sim.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
data = sim.output.spectra.load_kzkh_mean(
tmin, key_to_load=["A", "Khd", "Kz", "Khr"]
)
kh = data["kh_spectra"]
kz = data["kz"]
delta_kh = kh[1]
delta_kz = kz[1]
if key == "EA":
spectrum = data["A"]
elif key == "Epolo":
EKhd = data["Khd"]
EKz = data["Kz"]
spectrum = EKhd + EKz
elif key == "Ee":
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
Epolo = EKhd + EKz
spectrum = 2 * np.minimum(EA, Epolo)
elif key == "Ed":
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
Epolo = EKhd + EKz
spectrum = EA + Epolo - 2 * np.minimum(EA, Epolo)
elif key == "Etoro":
spectrum = data["Khr"]
else:
print(f"Don't know key: {key} \n")
exit
k = np.sqrt(kh**2 + kz**2)
k_nozero = k.copy()
k_nozero[k_nozero == 0] = 1e-16
sint = kh / k_nozero
theta = np.arcsin(sint)
cs = ax.pcolormesh(
kh,
kz,
np.log10(spectrum),
cmap=cm,
vmin=-7.0,
vmax=-1.0,
shading="nearest",
)
th = np.linspace(0, np.pi / 2, 100, endpoint=True)
ax.plot(kb * np.sin(th), kb * np.cos(th), color="k", linestyle="dotted")
ax.plot(ko * np.sin(th), ko * np.cos(th), "k--")
a = 3
xa = np.linspace(delta_kh, a**1.5 * ko, 50, endpoint=True)
ax.plot(
xa,
xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dashed",
color="gray",
)
a = 1 / 3
xa = np.linspace(delta_kh, a**1.5 * ko, 50, endpoint=True)
ax.plot(
xa,
xa * np.sqrt((a**1.5 * ko / xa) ** 0.8 - 1),
linestyle="dotted",
color="gray",
)
ax.plot([delta_kh, max(kh)], [delta_kh, max(kh)], "k-")
ax.plot(kf * np.sin(th), kf * np.cos(th), linestyle="--", color="orange")
ax.plot(keta * np.sin(th), keta * np.cos(th), linestyle="--", color="g")
ax.set_xlim([delta_kh, kmax])
ax.set_ylim([delta_kz, kmax])
ax.set_xscale("log")
ax.set_yscale("log")
ax.grid(True)
return cs
Ns = [10, 20, 80]
nbax = 0
css = [None for i in range(6)]
fig, axes = plt.subplots(
ncols=2, nrows=3, figsize=(10, 3 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ax4 = axes[2, 0]
ax5 = axes[2, 1]
axs = [ax0, ax1, ax2, ax3, ax4, ax5]
for N in Ns:
for proj in [False, True]:
path = get_path_finer_resol(N=N, Rb=None, proj=proj, ratio_one=True)
sim = load(path, hide_stdout=True)
css[nbax] = plot_spectra(sim, axs[nbax], key="Ee")
nbax += 1
for ax in [ax0, ax2, ax4]:
ax.set_ylabel(r"$k_z$", fontsize=16)
for ax in [ax4, ax5]:
ax.set_xlabel(r"$k_h$", fontsize=16)
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"$(a)$", fontsize=16)
ax1.set_title(r"Without vortical modes" + "\n" + r"$(b)$", fontsize=16)
ax2.set_title(r"$(c)$", fontsize=16)
ax3.set_title(r"$(d)$", fontsize=16)
ax4.set_title(r"$(e)$", fontsize=16)
ax5.set_title(r"$(f)$", fontsize=16)
"""
norm = matplotlib.colors.Normalize(vmin=-7, vmax=-3)
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.28, 0.02, 0.35])
cbar = fig.colorbar(matplotlib.cm.ScalarMappable(norm=norm, cmap=cm), cax=cbar_ax, cmap=cm, orientation='vertical')
cbar.set_ticks([-7, -6, -5, -4, -3], position='top')
cbar.set_label(r"$E_{equi}$", fontsize=16)
"""
fig.tight_layout()
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.28, 0.02, 0.35])
cbar = fig.colorbar(css[0], cax=cbar_ax, cmap=cm, orientation="vertical")
cbar.set_label(r"$E_{equi}(k_h, k_z)$", fontsize=16)
save_fig(fig, "fig_spectra_kh_kz_ratio_one_2D.png")
if __name__ == "__main__":
plt.show()
import glob
import os
import sys
from pathlib import Path
import h5py
import matplotlib.cm
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.collections import LineCollection
from util_simuls_regimes import get_sim
from fluidsim import load
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,
formatters,
get_path_finer_resol,
get_paths,
paths_simuls_regimes,
paths_simuls_regimes_proj,
pos_closest_value,
save_fig,
)
cm = matplotlib.cm.get_cmap("inferno", 100)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
print(sys.argv)
letter = sys.argv[-1]
if letter not in "DLOWPU":
letter = "L"
sim = get_sim(letter)
path = paths_simuls_regimes[letter]
sim_proj = get_sim(letter, proj=True)
path_proj = paths_simuls_regimes_proj[letter]
# 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=2, figsize=(10, 2 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ks_mins = [-0.02, 0.18, 0.38, 0.58, 0.78, 0.98]
ks_maxs = [0.02, 0.22, 0.42, 0.62, 0.82, 1.02]
# Standard Navier-Stokes
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)
mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
R2 = mean_values["R2"]
Fh = mean_values["Fh"]
epsK = mean_values["epsK"]
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:
# Get the data
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)
Ed = EA + Epolo - Ee
spectrum = Ee
omega_emp, delta_omega_emp = compute_omega_emp_vs_kzkh(
N, spectrum, kh, kz, omegas
)
KH, KZ = np.meshgrid(kh, kz)
K = (KH**2 + KZ**2) ** 0.5
K_NOZERO = K.copy()
K_NOZERO[K_NOZERO == 0] = 1e-16
omega_disp = N * KH / K_NOZERO
nlb = delta_omega_emp / omega_disp
chi = (K**2 * epsK) ** (1 / 3) / omega_disp
DO_PLOT = [0 for i in range(len(ks_maxs))]
for nh in range(len(kh)):
for nz in range(len(kz)):
k = np.sqrt(kh[nh] ** 2 + kz[nz] ** 2)
sint = kh[nh] / k
omega_waves = N * sint
theta = np.arcsin(sint)
color = cm(k / kb)
for i in range(len(ks_maxs)):
if (
k / kb > ks_mins[i]
and k / kb < ks_maxs[i]
and max(spectrum[nz, nh, :]) >= 1e-16
and DO_PLOT[i] < 3
):
cs = ax0.plot(
(omegas - omega_waves) / N,
spectrum[nz, nh, :] / max(spectrum[nz, nh, :]),
color=color,
)
cs = ax2.plot(
(omegas - omega_waves) / N,
spectrum[nz, nh, :] / E[nz, nh, :],
color=color,
)
DO_PLOT[i] += 1
# Without vortical modes
t_start, t_last = sim_proj.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim_proj, tmin)
mean_values = sim_proj.output.get_mean_values(tmin=tmin, customize=customize)
R2 = mean_values["R2"]
Fh = mean_values["Fh"]
epsK = mean_values["epsK"]
N = sim_proj.params.N
path_spec = sorted(path_proj.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:
# Get the data
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)
Ed = EA + Epolo - Ee
spectrum = Ee
omega_emp, delta_omega_emp = compute_omega_emp_vs_kzkh(
N, spectrum, kh, kz, omegas
)
KH, KZ = np.meshgrid(kh, kz)
K = (KH**2 + KZ**2) ** 0.5
K_NOZERO = K.copy()
K_NOZERO[K_NOZERO == 0] = 1e-16
omega_disp = N * KH / K_NOZERO
nlb = delta_omega_emp / omega_disp
DO_PLOT = [0 for i in range(len(ks_maxs))]
for nh in range(len(kh)):
for nz in range(len(kz)):
k = np.sqrt(kh[nh] ** 2 + kz[nz] ** 2)
sint = kh[nh] / k
omega_waves = N * sint
theta = np.arcsin(sint)
color = cm(k / kb)
for i in range(len(ks_maxs)):
if (
k / kb > ks_mins[i]
and k / kb < ks_maxs[i]
and max(spectrum[nz, nh, :]) >= 1e-16
and DO_PLOT[i] < 3
):
cs = ax1.plot(
(omegas - omega_waves) / N,
spectrum[nz, nh, :] / max(spectrum[nz, nh, :]),
color=color,
linestyle="-",
)
cs = ax3.plot(
(omegas - omega_waves) / N,
spectrum[nz, nh, :] / E[nz, nh, :],
color=color,
linestyle="-",
)
DO_PLOT[i] += 1
# Axis
for ax in [ax0, ax1, ax2, ax3]:
ax.set_xlim([-1, 3])
ax.plot([0, 0], [1e-3, 1e0], "k-")
ax.grid(True)
for ax in [ax0, ax1]:
ax.set_yscale("log")
ax.set_ylim([1e-3, 1e0])
for ax in [ax2, ax3]:
ax.set_ylim([0, 1e0])
for ax in [ax2, ax3]:
ax.set_xlabel(r"$(\omega - \omega_{\bm{k}})/N$", fontsize=16)
for ax in [ax0, ax1]:
ax.set_xticklabels([])
for ax in [ax1, ax3]:
ax.set_yticklabels([])
ax0.set_ylabel(
r"$E_{equi}(k_h, k_z, \omega) / \max\limits_{\omega} ~ E_{equi}(k_h, k_z, \omega)$",
fontsize=16,
)
ax2.set_ylabel(r"$E_{equi}(k_h, k_z, \omega) / E(k_h, k_z, \omega)$", fontsize=16)
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$(a)$", fontsize=16)
ax1.set_title(r"Without vortical modes" + "\n" + r"$(b)$", fontsize=16)
ax2.set_title(r"$(c)$", fontsize=16)
ax3.set_title(r"$(d)$", fontsize=16)
cmap = matplotlib.cm.binary
norm = matplotlib.colors.Normalize(vmin=0, vmax=1)
"""
fig.subplots_adjust(bottom=0.2)
cbar_ax = fig.add_axes([0.2, 0.1, 0.6, 0.02])
cbar = fig.colorbar(matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap), cax=cbar_ax, cmap=cm, orientation='horizontal')
cbar.set_ticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], position='top')
cbar.set_label(r"$k/k_b$", fontsize = 12, rotation=0)
"""
(leg0,) = ax3.plot([-2, -2], [0, 0], color=cm(0.0))
(leg1,) = ax3.plot([-2, -2], [0, 0], color=cm(0.2))
(leg2,) = ax3.plot([-2, -2], [0, 0], color=cm(0.4))
(leg3,) = ax3.plot([-2, -2], [0, 0], color=cm(0.6))
(leg4,) = ax3.plot([-2, -2], [0, 0], color=cm(0.8))
(leg5,) = ax3.plot([-2, -2], [0, 0], color=cm(1.0))
ax3.legend(
[leg1, leg2, leg3, leg4, leg5],
[
r"$k/k_b \simeq 0.2$",
r"$k/k_b \simeq 0.4$",
r"$k/k_b \simeq 0.6$",
r"$k/k_b \simeq 0.8$",
r"$k/k_b \simeq 1.0$",
],
loc="lower right",
fontsize=16,
)
fig.tight_layout()
save_fig(fig, f"fig_spectra_slices_omega_kh_kz_regimes_{letter}.png")
if __name__ == "__main__":
plt.show()
import glob
import os
from pathlib import Path
import h5py
import matplotlib.cm
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.collections import LineCollection
from util_simuls_regimes import get_sim
from fluidsim import load
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,
formatters,
get_path_finer_resol,
get_paths,
paths_simuls_regimes,
paths_simuls_regimes_proj,
pos_closest_value,
save_fig,
)
cm = matplotlib.cm.get_cmap("inferno", 100)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
###
N = 120
Rb = 10
ratio_one = True
###
path = get_path_finer_resol(N, Rb, proj=False, ratio_one=ratio_one)
sim = load(path)
path_proj = get_path_finer_resol(N, Rb, proj=True, ratio_one=ratio_one)
sim_proj = load(path_proj)
# 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=2, figsize=(10, 2 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ks_mins = [-0.02, 0.18, 0.38, 0.58, 0.78, 0.98]
ks_maxs = [0.02, 0.22, 0.42, 0.62, 0.82, 1.02]
# Standard Navier-Stokes
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)
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:
# Get the data
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)
Ed = EA + Epolo - Ee
spectrum = Ee
omega_emp, delta_omega_emp = compute_omega_emp_vs_kzkh(
N, spectrum, kh, kz, omegas
)
KH, KZ = np.meshgrid(kh, kz)
K = (KH**2 + KZ**2) ** 0.5
K_NOZERO = K.copy()
K_NOZERO[K_NOZERO == 0] = 1e-16
omega_disp = N * KH / K_NOZERO
nlb = delta_omega_emp / omega_disp
DO_PLOT = [0 for i in range(len(ks_maxs))]
for nh in range(len(kh)):
for nz in range(len(kz)):
k = np.sqrt(kh[nh] ** 2 + kz[nz] ** 2)
sint = kh[nh] / k
omega_waves = N * sint
theta = np.arcsin(sint)
if nlb[nz, nh] <= 1 / 3 and k <= kb:
color = cm(k / kb)
cs = ax0.scatter(
(omegas - omega_waves) / N,
spectrum[nz, nh, :] / max(spectrum[nz, nh, :]),
color=color,
)
cs = ax2.scatter(
(omegas - omega_waves) / N,
spectrum[nz, nh, :] / E[nz, nh, :],
color=color,
)
# Without vortical modes
t_start, t_last = sim_proj.output.print_stdout.get_times_start_last()
tmin = t_last - 2.0
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim_proj, tmin)
N = sim_proj.params.N
path_spec = sorted(path_proj.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:
# Get the data
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)
Ed = EA + Epolo - Ee
spectrum = Ee
omega_emp, delta_omega_emp = compute_omega_emp_vs_kzkh(
N, spectrum, kh, kz, omegas
)
KH, KZ = np.meshgrid(kh, kz)
K = (KH**2 + KZ**2) ** 0.5
K_NOZERO = K.copy()
K_NOZERO[K_NOZERO == 0] = 1e-16
omega_disp = N * KH / K_NOZERO
nlb = delta_omega_emp / omega_disp
DO_PLOT = [0 for i in range(len(ks_maxs))]
for nh in range(len(kh)):
for nz in range(len(kz)):
k = np.sqrt(kh[nh] ** 2 + kz[nz] ** 2)
sint = kh[nh] / k
omega_waves = N * sint
theta = np.arcsin(sint)
if nlb[nz, nh] <= 1 / 3 and k <= kb:
color = cm(k / kb)
cs = ax1.plot(
(omegas - omega_waves) / N,
spectrum[nz, nh, :] / max(spectrum[nz, nh, :]),
color=color,
linestyle="--",
)
cs = ax3.plot(
(omegas - omega_waves) / N,
spectrum[nz, nh, :] / E[nz, nh, :],
color=color,
linestyle="--",
)
# Axis
for ax in [ax0, ax1, ax2, ax3]:
ax.set_xlim([-1, 3])
ax.plot([0, 0], [1e-3, 1e0], "k-")
ax.grid(True)
for ax in [ax0, ax1]:
ax.set_yscale("log")
ax.set_ylim([1e-3, 1e0])
for ax in [ax2, ax3]:
ax.set_ylim([0, 1e0])
for ax in [ax2, ax3]:
ax.set_xlabel(r"$(\omega - \omega_{\bm{k}})/N$", fontsize=16)
for ax in [ax0, ax1]:
ax.set_xticklabels([])
for ax in [ax1, ax3]:
ax.set_yticklabels([])
ax0.set_ylabel(
r"$E_{equi}(k_h, k_z, \omega) / \max\limits_{\omega} ~ E_{equi}(k_h, k_z, \omega)$",
fontsize=16,
)
ax2.set_ylabel(r"$E_{equi}(k_h, k_z, \omega) / E(k_h, k_z, \omega)$", fontsize=16)
ax0.set_title(r"$(a)$", fontsize=16)
ax1.set_title(r"$(b)$", fontsize=16)
ax2.set_title(r"$(c)$", fontsize=16)
ax3.set_title(r"$(d)$", fontsize=16)
cmap = matplotlib.cm.binary
norm = matplotlib.colors.Normalize(vmin=0, vmax=1)
"""
fig.subplots_adjust(bottom=0.2)
cbar_ax = fig.add_axes([0.2, 0.1, 0.6, 0.02])
cbar = fig.colorbar(matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap), cax=cbar_ax, cmap=cm, orientation='horizontal')
cbar.set_ticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0], position='top')
cbar.set_label(r"$k/k_b$", fontsize = 12, rotation=0)
"""
(leg0,) = ax3.plot([-2, -2], [0, 0], color=cm(0.0))
(leg1,) = ax3.plot([-2, -2], [0, 0], color=cm(0.2))
(leg2,) = ax3.plot([-2, -2], [0, 0], color=cm(0.4))
(leg3,) = ax3.plot([-2, -2], [0, 0], color=cm(0.6))
(leg4,) = ax3.plot([-2, -2], [0, 0], color=cm(0.8))
(leg5,) = ax3.plot([-2, -2], [0, 0], color=cm(1.0))
ax3.legend(
[leg1, leg2, leg3, leg4, leg5],
[
r"$k/k_b \simeq 0.2$",
r"$k/k_b \simeq 0.4$",
r"$k/k_b \simeq 0.6$",
r"$k/k_b \simeq 0.8$",
r"$k/k_b \simeq 1.0$",
],
loc="lower right",
fontsize=16,
)
fig.tight_layout()
save_fig(fig, f"fig_spectra_slices_omega_kh_kz_regimes_ratio_one_N{N}.png")
if __name__ == "__main__":
plt.show()
import sys
import matplotlib.pyplot as plt
import numpy as np
from util_simuls_regimes import get_sim
from fluidsim import load
from util import (
compute_kf_kb_ko_keta_kd,
customize,
get_paths,
pos_closest_value,
save_fig,
)
print(sys.argv)
letter = sys.argv[-1]
if letter not in "DLOWP":
letter = "P"
sim = get_sim(letter)
coef_compensate = 5 / 3
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
if letter == "P":
fig, ax = plt.subplots()
temp = sim.output.spectra.load3d_mean(tmin)
print(temp)
EK = temp["spectra_E"]
EA = temp["spectra_A"]
EKhd = temp["spectra_Khd"]
EKhr = temp["spectra_Khr"]
EKh = EKhr + EKhd
EKz = EK - EKh
Epolo = EKhd + EKz
Etoro = EKhr
k = temp["k"]
ax.plot(
k,
EA * k**coef_compensate,
"b-",
label=(r"$E_{A}$"),
)
ax.plot(
k,
Epolo * k**coef_compensate,
"m-",
label=(r"$E_{polo}$"),
)
ax.plot(
k,
Etoro * k**coef_compensate,
"r-",
label=(r"$E_{toro}$"),
)
ax.legend(fontsize=10, loc="upper left")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel(r"$k$", fontsize=16)
ax.set_ylabel(r"$E(k) ~ k^{5/3}$", fontsize=16)
# ax.text(1.1*kb, 1e-10, r"$k_b$", fontsize=16)
# ax.text(1.1*ko, 1e-10, r"$k_O$", fontsize=16)
ax.set_xlim([k[1], 0.8 * max(k)])
ax.set_ylim(bottom=1e-4, top=1e1)
fig.tight_layout()
save_fig(fig, f"fig_spectra_slices_regime_{letter}_3D.png")
if True:
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_kz = kz[1]
delta_kh = kh[1]
EA = data["A"]
EKhd = data["Khd"]
EKz = data["Kz"]
EKhr = data["Khr"]
Epolo = EKhd + EKz
Etoro = EKhr
EKh = EKhr + EKhd
kstxts = ["kfkb", "kbko", "koketa"]
kss = [0.5 * kb, (kb * ko) ** 0.5, (ko * keta) ** 0.5]
for n in range(3):
ks = kss[n]
kstxt = kstxts[n]
ikz = pos_closest_value(kz, ks)
ikh = pos_closest_value(kh, ks)
print(ks, kz[ikz], ikz)
print(ks, kh[ikh], ikh)
EA_kz = EA[:, ikh]
Epolo_kz = Epolo[:, ikh]
Etoro_kz = Etoro[:, ikh]
EKz_kz = EKz[:, ikh]
EKh_kz = EKh[:, ikh]
EA_kh = EA[ikz, :]
Epolo_kh = Epolo[ikz, :]
Etoro_kh = Etoro[ikz, :]
EKz_kh = EKz[ikz, :]
EKh_kh = EKh[ikz, :]
fig, ax = plt.subplots() # figsize=(10, 3 * 4.5 / 2))
if kstxt == "kfkb":
kstitle = (
r"$k_h=$"
+ f"{kh[ikh]/kb:.1f}"
+ r"$k_b ~ (--)$"
+ " or "
+ r"$k_z=$"
+ f"{kz[ikz]/kb:.1f}"
+ r"$k_b ~ (-)$"
)
elif kstxt == "kbko":
kstitle = (
r"$k_h=$"
+ f"{kh[ikh]/ko:.1f}"
+ r"$k_O ~ (--)$"
+ " or "
+ r"$k_z=$"
+ f"{kz[ikz]/ko:.1f}"
+ r"$k_O ~ (-)$"
)
else:
kstitle = (
r"$k_h=$"
+ f"{kh[ikh]/ko:.1f}"
+ r"$k_O ~ (--)$"
+ " or "
+ r"$k_z=$"
+ f"{kz[ikz]/ko:.1f}"
+ r"$k_O ~ (-)$"
)
ax.plot(
kz,
EA_kz * kz**coef_compensate,
"b--",
label=None,
)
ax.plot(
kz,
Epolo_kz * kz**coef_compensate,
"m--",
label=None,
)
ax.plot(
kz,
Etoro_kz * kz**coef_compensate,
"r--",
label=None,
)
"""
ax.plot(
kz,
EKz_kz * kz ** coef_compensate,
"g--",
)
ax.plot(
kz,
EKh_kz * kz ** coef_compensate,
"y--",
)
"""
ax.plot(
kh,
EA_kh * kh**coef_compensate,
"b-",
label=(r"$E_A$"),
)
ax.plot(
kh,
Epolo_kh * kh**coef_compensate,
"m-",
label=(r"$E_{polo}$"),
)
ax.plot(
kh,
Etoro_kh * kh**coef_compensate,
"r-",
label=(r"$E_{toro}$"),
)
"""
ax.plot(
kh,
EKz_kh * kh ** coef_compensate,
"g-",
label=(
r"$E_{z}$"
),
)
ax.plot(
kh,
EKh_kh * kh ** coef_compensate,
"y-",
label=(
r"$E_{h}$"
),
)
"""
if letter != "P":
if kstxt == "kfkb":
x = [delta_kz, kb]
y = [
5e-3 * x[0] ** (-2 + coef_compensate),
5e-3 * x[1] ** (-2 + coef_compensate),
]
ax.plot(x, y, "k-")
ax.text(
(x[0] * x[1]) ** 0.5,
0.1 * 5e-3 * ((x[0] * x[1]) ** 0.5) ** (-2 + coef_compensate),
r"$k_h^{-2}$",
fontsize=16,
)
ax.legend(fontsize=10, loc="lower left")
if kstxt == "kbko":
x = [kb, ko]
y = [
2e-4 * kb ** (-5 / 3 + coef_compensate),
2e-4 * ko ** (-5 / 3 + coef_compensate),
]
ax.plot(x, y, "k-")
ax.text(
ks,
0.1 * 2e-4 * ks ** (-5 / 3 + coef_compensate),
r"$k_h^{-5/3}$",
fontsize=16,
)
# ax.legend(fontsize=10, loc = "lower left")
if kstxt == "koketa":
x = [delta_kz, kb]
y = [
1e-11 * delta_kz ** (1 + coef_compensate),
1e-11 * kb ** (1 + coef_compensate),
]
ax.plot(x, y, "k-")
ax.text(
(delta_kz * kb) ** 0.5,
0.2
* 1e-11
* ((delta_kz * kb) ** (0.5 * (1 + coef_compensate))),
r"$k_h^1$",
fontsize=16,
)
# ax.legend(fontsize=10, loc = "upper left")
x = [delta_kz, ko]
y = [1e-9 * delta_kz**coef_compensate, 1e-9 * ko**coef_compensate]
ax.plot(x, y, "k--")
ax.text(
(x[0] * x[1]) ** 0.5,
0.07 * 1e-9 * ((kb * ko) ** (0.5 * coef_compensate)),
r"$k_z^0$",
fontsize=16,
)
ax.axvline(kb, color="k", linestyle="dotted")
ax.text(1.1 * kb, 2e-10, r"$k_b$", fontsize=16)
ax.axvline(ko, color="k", linestyle="dashed")
ax.text(1.1 * ko, 2e-10, r"$k_O$", fontsize=16)
# ax.legend(fontsize=10, loc = "lower right")
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel(r"$k_h, k_z$", fontsize=16)
ax.set_ylabel(
r"$E_i ~ k_h^{5/3} ~ (-)$ or $E_i ~ k_z^{5/3} (--)$", fontsize=16
)
# ax.text(1.1*kb, 1e-10, r"$k_b$", fontsize=16)
# ax.text(1.1*ko, 1e-10, r"$k_O$", fontsize=16)
ax.set_xlim([kh[1], 0.8 * max(kh)])
ax.set_ylim(bottom=1e-10, top=1e0)
# plt.legend(loc=2, prop={'size': 10})
ax.set_title(kstitle, fontsize=12)
fig.tight_layout()
save_fig(fig, f"fig_spectra_slices_regime_{letter}_{kstxt}.png")
if __name__ == "__main__":
plt.show()
from util_dataframe import df, df_proj
from util import formatters, tmp_dir
# print(df.columns)
# fmt: off
columns = [
"N", "Rb", #"Re",
"nx", "nz", "k_max*eta",
# "R4",
#"epsK2/epsK", "k_max*lambda",
"Fh", "R2", #"Re_lambda",
#"Uh2", "epsK",
# "Gamma", "I_velocity", "I_dissipation",
]
header = [
"$N$", "$\R_i$", #"$Re_i$",
"$n_h$", "$n_z$",r"$\kmax\eta$",
# "$\R_4$",
#r"$\epsK_2/\epsK$", r"$\kmax\lambda$",
"$F_h$", r"$\R$", #r"$Re_\lambda$"
#"${U_h}^2$", r"$\epsK$",
# "$\Gamma$", "$I_{velo}$", "$I_{diss}$",
]
# fmt: on
column_format = "rrrr|" + 3 * "r"
df.to_latex(
buf=tmp_dir / "table_better_simuls.tex",
columns=columns,
formatters=formatters,
column_format=column_format,
index=False,
header=header,
escape=False,
caption=(
r"List of simulations without poloidal projection (Standard Navier-Stokes)."
),
label="table-better-simuls",
)
df_proj.to_latex(
buf=tmp_dir / "table_better_simuls_proj.tex",
columns=columns,
formatters=formatters,
column_format=column_format,
index=False,
header=header,
escape=False,
caption=(
r"List of simulations with poloidal projection (without vortical modes)."
),
label="table-better-simuls-proj",
)
from util_dataframe import df_proj_ratio_one, df_ratio_one
from util import formatters, tmp_dir
# print(df.columns)
# fmt: off
columns = [
"N",
"nx", "k_max*eta",
#"epsK2/epsK", "k_max*lambda",
"Fh", "R4", #"Re_lambda",
#"Uh2", "epsK",
# "Gamma", "I_velocity", "I_dissipation",
]
header = [
"$N$",
"$n_h=n_z$",r"$\kmax\eta_4$",
#r"$\epsK_2/\epsK$", r"$\kmax\lambda$",
"$F_h$", r"$\R_4$", #r"$Re_\lambda$"
#"${U_h}^2$", r"$\epsK$",
# "$\Gamma$", "$I_{velo}$", "$I_{diss}$",
]
# fmt: on
column_format = "rr|" + 3 * "r"
df_ratio_one.to_latex(
buf=tmp_dir / "table_better_simuls_ratio_one.tex",
columns=columns,
formatters=formatters,
column_format=column_format,
index=False,
header=header,
escape=False,
caption=(
r"List of simulations with aspect ratio one and Standard Navier-Stokes."
),
label="table-better-simuls-ratio-one",
)
df_proj_ratio_one.to_latex(
buf=tmp_dir / "table_better_simuls_proj_ratio_one.tex",
columns=columns,
formatters=formatters,
column_format=column_format,
index=False,
header=header,
escape=False,
caption=(
r"List of simulations with aspect ratio one and without vortical modes."
),
label="table-better-simuls-proj-ratio-one",
)
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_{velo}$', '$I_{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",
)
import os
import sys
from itertools import product
from math import sqrt
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from fluiddyn.util import has_to_be_made as _has_to_be_made
from fluidsim.util import get_dataframe_from_paths, times_start_last_from_path
# import scipy.optimize as optimize
path_base = os.environ["STRAT_TURB_POLO2022"]
path_base_proj = os.environ["STRAT_TURB_POLO_PROJ2022"]
path_base_ratio_one = os.environ["STRAT_WAVES2022"]
paths_all = sorted(Path(path_base).glob("simul_folders/ns3d*"))
paths_all_proj = sorted(Path(path_base_proj).glob("simul_folders/ns3d*"))
paths_all_ratio_one = sorted(
Path(path_base_ratio_one).glob("simul_folders/ns3d*")
)
here = Path(__file__).absolute().parent
tmp_dir = here.parent / "tmp"
tmp_dir.mkdir(exist_ok=True)
def has_to_be_made(name, sources: list):
if not isinstance(name, str):
names = name
return any(has_to_be_made(name, sources) for name in names)
if isinstance(sources, str):
sources = [sources]
sources.append("util.py")
if not any(name.endswith(ext) for ext in (".png", ".tex")):
name += ".png"
return _has_to_be_made(tmp_dir / name, sources, source_dir=here)
height = 3.7
plt.rc("figure", figsize=(1.33 * height, height))
def get_paths(N, Rb, nh, proj=False, ratio_one=False, reverse=False):
str_N = f"_N{N}_"
str_nh = f"_{nh}x{nh}"
if not ratio_one:
str_Rb = f"_Rb{Rb:.3g}_"
str_Rb2 = f"_Rb{Rb}_"
if proj:
paths_tmp = paths_all_proj
else:
paths_tmp = paths_all
paths = [
p
for p in paths_tmp
if str_N in p.name
and (str_Rb in p.name or str_Rb2 in p.name)
and str_nh in p.name
]
else:
str_proj = "_projpoloidal_"
paths_tmp = paths_all_ratio_one
if proj:
paths = [
p
for p in paths_tmp
if str_N in p.name and str_proj in p.name and str_nh in p.name
]
else:
paths = [
p
for p in paths_tmp
if str_N in p.name and str_proj not in p.name and str_nh in p.name
]
paths.sort(key=lambda p: int(p.name.split("x")[1]), reverse=reverse)
return paths[0]
def get_paths_couple(N, Rb, proj=False, ratio_one=False, reverse=False):
str_N = f"_N{N}_"
if not ratio_one:
str_Rb = f"_Rb{Rb:.3g}_"
str_Rb2 = f"_Rb{Rb}_"
if proj:
paths_tmp = paths_all_proj
else:
paths_tmp = paths_all
paths_couple = [
p
for p in paths_tmp
if str_N in p.name and (str_Rb in p.name or str_Rb2 in p.name)
]
else:
str_proj = "_projpoloidal_"
paths_tmp = paths_all_ratio_one
if proj:
paths_couple = [
p for p in paths_tmp if str_N in p.name and str_proj in p.name
]
else:
paths_couple = [
p for p in paths_tmp if str_N in p.name and str_proj not in p.name
]
paths_couple.sort(key=lambda p: int(p.name.split("x")[1]), reverse=reverse)
return paths_couple
def get_path_finer_resol(N, Rb, proj=False, ratio_one=False):
paths_couple = get_paths_couple(N, Rb, proj, ratio_one, reverse=True)
for path in paths_couple:
t_start, t_last = times_start_last_from_path(path)
if t_last > t_start + 1:
return path
def lprod(a, b):
return list(product(a, b))
couples320 = set(
lprod([10, 20, 40], [5, 10, 20, 40, 80, 160])
+ lprod([30], [10, 20, 40])
+ lprod([6.5], [100, 200])
+ lprod([4], [250, 500])
+ lprod([3], [450, 900])
+ lprod([2], [1000, 2000])
+ lprod([0.66], [9000, 18000])
+ [(14.5, 20), (5.2, 150), (2.9, 475), (1.12, 3200), (0.25, 64000)]
)
couples320.add((60, 10))
couples320.add((60, 20))
couples320.add((80, 10))
couples320.add((100, 10))
couples320.add((120, 10))
couples320.remove((40, 160))
# Small Rb
couples320.update(lprod([20], [1, 2]))
couples320.update(lprod([40], [1, 2]))
couples320.update(lprod([80], [0.5, 1]))
N_ratio_one = [10, 20, 50, 80, 120]
has_to_save = "SAVE" in sys.argv
def save_fig(fig, name):
if has_to_save:
print(f"saving file {tmp_dir.name}/{name}")
fig.savefig(tmp_dir / name, dpi=300)
def customize(result, sim):
EKh = result["EKh"]
EKz = result["EKz"]
EK = EKh + EKz
U = sqrt(2 * EK / 3)
nu_2 = sim.params.nu_2
epsK = result["epsK"]
result["name"] = sim.output.name_run
if nu_2 != 0.0:
result["lambda"] = sqrt(U**2 * nu_2 / epsK)
result["Re_lambda"] = U * result["lambda"] / nu_2
result["Rb"] = float(sim.params.short_name_type_run.split("_Rb")[-1])
# else:
# result["lambda"] = INFINITY
# result["Re_lambda"] = INFINITY
# result["Rb"] = INFINITY
result["nx"] = sim.params.oper.nx
result["nz"] = sim.params.oper.nz
result["proj"] = sim.params.projection
def get_customized_dataframe(paths):
df = get_dataframe_from_paths(
paths, tmin="t_last-2", use_cache=1, customize=customize
)
if "nu0.0" in paths[0].name:
columns_old = df.columns.tolist()
# fmt: off
first_columns = [
"N", "nx", "nz", "Fh", "R4", "Gamma",
"lx1", "lx2", "lz1", "lz2", "I_velocity", "I_dissipation"]
# fmt: on
else:
df["Re"] = df.Rb * df.N**2
columns_old = df.columns.tolist()
# fmt: off
first_columns = [
"N", "Rb", "Re", "nx", "nz", "Fh", "R2", "k_max*eta", "epsK2/epsK", "Gamma",
"lx1", "lx2", "lz1", "lz2", "I_velocity", "I_dissipation"]
# fmt: on
columns = first_columns.copy()
for key in columns_old:
if key not in columns:
columns.append(key)
df = df[columns]
return df
def plot(
df,
x,
y,
logx=True,
logy=False,
c=None,
cmap=None,
vmin=None,
vmax=None,
s=None,
ax=None,
):
ax = df.plot.scatter(
x=x,
y=y,
logx=logx,
logy=logy,
c=c,
edgecolors="k",
vmin=vmin,
vmax=vmax,
s=s,
ax=ax,
)
if c is not None:
pc = ax.collections[-1]
if cmap == "seismic":
pc.set_cmap("seismic")
else:
pc.set_cmap("binary")
cbar = plt.colorbar(pc, ax=ax)
if vmin == -1 and vmax == 1:
cbar.set_ticks([-1.0, -0.5, 0.0, 0.5, 1.0])
return ax
N_1couple = 40 # 40
Rb_1couple = 20 # 20
paths_1couple = get_paths_couple(N_1couple, Rb_1couple)
print([p.name for p in paths_1couple])
params_simuls_regimes = {
"D": (40, 2),
"L": (40, 20), # (40, 20),
"O": (10, 80),
"W": (6.5, 200),
"P": (0.66, 18000),
"U": (80, 10),
}
paths_simuls_regimes = {
k: get_path_finer_resol(*params)
for k, params in params_simuls_regimes.items()
}
paths_simuls_regimes = {
k: v for k, v in paths_simuls_regimes.items() if v is not None
}
paths_simuls_regimes_proj = {
k: get_path_finer_resol(*params, proj=True)
for k, params in params_simuls_regimes.items()
}
paths_simuls_regimes_proj = {
k: v for k, v in paths_simuls_regimes_proj.items() if v is not None
}
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.13
R2_limit = 10.0
def compute_kf_kb_ko_keta_kd(sim, tmin):
mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
# data = sim.output.spectra.load_kzkh_mean(tmin)
# kh = data["kh_spectra"]
# delta_kh = kh[1]
# assert delta_kh == 2 * np.pi / sim.params.oper.Lx
delta_kh = 2 * np.pi / sim.params.oper.Lx
if sim.params.oper.nx == sim.params.oper.nz:
kf = 3.5 * delta_kh
else:
kf = 20 * delta_kh
N = sim.params.N
nu = sim.params.nu_2
nu_4 = sim.params.nu_4
Fh = mean_values["Fh"]
R4 = mean_values["R4"]
Uh2 = mean_values["Uh2"]
kb = N / Uh2**0.5
epsK = mean_values["epsK"]
ko = (N**3 / epsK) ** 0.5
if sim.params.nu_2 != 0.0:
R2 = mean_values["R2"]
keta = ko * R2**0.75
kd = (N / nu) ** 0.5
else:
keta = 1e16
kd = 1e16
ketah = ko * (Fh * R4) ** (3 / 10)
keta = min(keta, ketah)
kdh = (N / nu_4) ** 0.25
kd = min(kd, kdh)
return kf, kb, ko, keta, kd
def pos_closest_value(input_list, input_value):
arr = np.asarray(input_list)
i = (np.abs(arr - input_value)).argmin()
return i
def compute_omega_emp_vs_kzkh(
N,
spectrum,
kh_spectra,
kz_spectra,
omegas,
):
KH, KZ = np.meshgrid(kh_spectra, kz_spectra)
K = (KH**2 + KZ**2) ** 0.5
K_NOZERO = K.copy()
K_NOZERO[K_NOZERO == 0] = 1e-16
omega_disp = N * KH / K_NOZERO
omega_emp = np.zeros((len(kz_spectra), len(kh_spectra)))
delta_omega_emp = np.zeros((len(kz_spectra), len(kh_spectra)))
omega_norm = np.zeros((len(kz_spectra), len(kh_spectra)))
# we compute omega_emp first
for io in range(len(omegas)):
omega_emp += omegas[io] * spectrum[:, :, io]
omega_norm += spectrum[:, :, io]
omega_norm[omega_norm == 0] = 1e-16
omega_emp = omega_emp / omega_norm
# then we conpute delta_omega_emp
for io in range(len(omegas)):
delta_omega_emp += ((omegas[io] - omega_disp) ** 2) * spectrum[:, :, io]
delta_omega_emp = (np.divide(delta_omega_emp, omega_norm)) ** 0.5
return omega_emp, delta_omega_emp
def spectral_fit_kh_kz(
spectrum, kh, kz, kmin, kmax, thetamin=0.0, thetamax=np.pi / 2, plot=False
):
F = spectrum.flatten()
KH = kh.flatten()
KZ = kz.flatten()
K = np.sqrt(KH**2 + KZ**2)
SINT = KH / K
KH = KH[F >= 1e-16]
KZ = KZ[F >= 1e-16]
K = K[F >= 1e-16]
SINT = SINT[F >= 1e-16]
F = F[F >= 1e-16]
F = F[K >= kmin]
KH = KH[K >= kmin]
KZ = KZ[K >= kmin]
SINT = SINT[K >= kmin]
K = K[K >= kmin]
F = F[K <= kmax]
KH = KH[K <= kmax]
KZ = KZ[K <= kmax]
SINT = SINT[K <= kmax]
K = K[K <= kmax]
F = F[SINT <= np.sin(thetamax)]
KZ = KZ[SINT <= np.sin(thetamax)]
K = K[SINT <= np.sin(thetamax)]
KH = KH[SINT <= np.sin(thetamax)]
SINT = SINT[SINT <= np.sin(thetamax)]
F = F[SINT >= np.sin(thetamin)]
KH = KH[SINT >= np.sin(thetamin)]
K = K[SINT >= np.sin(thetamin)]
KZ = KZ[SINT >= np.sin(thetamin)]
SINT = SINT[SINT >= np.sin(thetamin)]
def error(params):
c, ah, az = params
res = np.sum(
(np.log(F) - np.log(c) - ah * np.log(KH) - az * np.log(KZ)) ** 2
)
return res
if plot:
plt.figure()
c = 0.01
for ah in np.linspace(-4, 4, 20):
for az in np.linspace(-4, 4, 20):
plt.scatter(
ah, az, c=np.log(error([c, ah, az])), vmin=-3, vmax=10
)
# print(error([c, ah, az]))
plt.show()
initial_guess = [1e-1, -2, 0]
result = optimize.minimize(error, initial_guess, method="Nelder-Mead")
if result.success:
fitted_params = result.x
print("[c, ah, az] = ", fitted_params)
return fitted_params
else:
raise ValueError(result.message)
def spectral_fit_k_sint(
spectrum, k, sint, kmin, kmax, thetamin=0.0, thetamax=np.pi / 2, plot=False
):
F = spectrum.flatten()
K = k.flatten()
SINT = sint.flatten()
K = K[F >= 1e-16]
SINT = SINT[F >= 1e-16]
F = F[F >= 1e-16]
F = F[K >= kmin]
SINT = SINT[K >= kmin]
K = K[K >= kmin]
F = F[K <= kmax]
SINT = SINT[K <= kmax]
K = K[K <= kmax]
F = F[SINT <= np.sin(thetamax)]
K = K[SINT <= np.sin(thetamax)]
SINT = SINT[SINT <= np.sin(thetamax)]
F = F[SINT >= np.sin(thetamin)]
K = K[SINT >= np.sin(thetamin)]
SINT = SINT[SINT >= np.sin(thetamin)]
def error(params):
c, ak, at = params
res = np.sum(
(np.log(F) - np.log(c) - ak * np.log(K) - at * np.log(SINT)) ** 2
)
return res
if plot:
plt.figure()
c = 0.01
for ak in np.linspace(-10, 10, 20):
for at in np.linspace(-10, 10, 20):
plt.scatter(
ak, at, c=np.log10(error([c, ak, at])), vmin=0, vmax=6
)
# print(error([c, ak, at]))
plt.show()
initial_guess = [1e-1, -2, 0]
result = optimize.minimize(error, initial_guess, method="Nelder-Mead")
if result.success:
fitted_params = result.x
print("[c, ah, az] = ", fitted_params)
return fitted_params
else:
raise ValueError(result.message)
import numpy as np
from fluidsim import load
from util import (
N_ratio_one,
compute_kf_kb_ko_keta_kd,
couples320,
get_customized_dataframe,
get_path_finer_resol,
path_base_ratio_one,
)
# Contruct dataframe
def construct_df(proj=False, ratio_one=False):
paths = []
if not ratio_one:
for N, Rb in sorted(couples320):
if N == 2.9:
continue
path = get_path_finer_resol(N, Rb, proj, ratio_one)
if path is not None:
paths.append(path)
else:
for N in N_ratio_one:
path = get_path_finer_resol(N, None, proj, ratio_one)
if path is not None:
paths.append(path)
print(f"Using {len(paths)} simulations")
df = get_customized_dataframe(paths)
if not ratio_one:
df["k_max*lambda"] = df["k_max"] * df["lambda"]
Etoro = []
Epolo = []
E = []
Kmaxeta = []
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)
E.append(e)
Epolo.append(epolo)
Etoro.append(etoro)
if ratio_one:
p_oper = sim.params.oper
kmax = p_oper.coef_dealiasing * p_oper.nx * np.pi / p_oper.Lx
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
kmaxeta = kmax / keta
Kmaxeta.append(kmaxeta)
df["E"] = E
df["Epolo"] = Epolo
df["Etoro"] = Etoro
if ratio_one:
df["k_max*eta"] = Kmaxeta
return df
df = construct_df(proj=False, ratio_one=False)
df_proj = construct_df(proj=True, ratio_one=False)
df_ratio_one = construct_df(proj=False, ratio_one=True)
df_proj_ratio_one = construct_df(proj=True, ratio_one=True)
from util import get_customized_dataframe, paths_simuls_regimes
df = get_customized_dataframe(paths_simuls_regimes.values())
df["letter"] = paths_simuls_regimes.keys()
regimes = [
"Diffusive",
"LAST",
"Optimal",
"Weakly stratified",
"Passive scalar",
]
regimes = {regime[0]: regime for regime in regimes}
print(regimes)
df["regime"] = [regimes[letter] for letter in df["letter"]]
print(df)
from fluidsim import load
from util import (
params_simuls_regimes,
paths_simuls_regimes,
paths_simuls_regimes_proj,
)
_simuls = {}
def get_sim(letter, proj=False):
if proj:
sim = load(paths_simuls_regimes_proj[letter], hide_stdout=True)
else:
sim = load(paths_simuls_regimes[letter], hide_stdout=True)
_simuls[letter] = sim
return sim
[project]
name = "2022strat_polo_waves"
version = "0.1.0"
description = "Article 2022strat_polo_waves"
authors = [
{name = "pierre.augier", email = "pierre.augier@univ-grenoble-alpes.fr"},
]
dependencies = [
"formattex",
"black",
"fluidsim>=0.8.4",
"fluidfft>=0.4.3",
"pyfftw>=0.15.0",
"jinja2>=3.1.4",
]
requires-python = "==3.11.*"
readme = "README.md"
license = {text = "MIT"}
[tool.pdm]
distribution = false