Skip to content
Snippets Groups Projects
save_figure17.py 6.96 KiB
Newer Older
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
vlabarre's avatar
vlabarre committed
# 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()