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
Commits on Source (11)
Showing
with 4741 additions and 0 deletions
[settings]
profile = black
known_first_party = fluiddyn,fluidsim
known_local_folder = util
NAME = article
$(NAME).pdf: figures
cd input && latexmk -f -shell-escape -pdf $(NAME).tex && rsync $(NAME).pdf ../$(NAME).pdf
figures:
pdm run python py/make_figures.py SAVE
clean:
rm -f tmp/*.tex
cd input && rm -f *.aux *.fdb_latexmk *.fls *.log *.bak* *.bbl *.blg *.out *Notes.bib
cleantmp:
rm -rf tmp/*
cleanpdf:
rm -f $(NAME).pdf input/$(NAME).pdf
cleanall: clean cleantmp cleanpdf
format:
cd input && pdm run formattex *.tex -i -v
black:
pdm run black -l 82 py
install-deps: sync
sync:
pdm sync --clean
lock:
pdm lock
# Source of a great article
## How to build the document?
1. One needs to download the dataset which is now available in [MyCore
CNRS](https://mycore.core-cloud.net/index.php/s/rpiQSLSoauNCVQE).
1. One needs to define the environment variable `STRAT_TURB_TORO2022` with the
path towards the dataset directory.
Something like this can be added in your `~/.bashrc` file:
```bash
export STRAT_TURB_TORO2022=/data/MyCore/2022strat-turb-toro
```
1. Install a virtual env with `pdm sync --clean` and activate it.
1. The command `make` should produce all Latex tables, figures and finally the
main pdf file `article.pdf`.
Note also than the Latex code can be formatted with the command `make format`.
File added
\documentclass[%
aps,
superscriptaddress,
longbibliography,
12pt,
onecolumn,
a4paper,
% preprint,
% reprint,
linenumbers,
showpacs,
showkeys,
% draft,
amsfonts, amssymb, amsmath]{revtex4-2}
% The result seems nicer with revtex4-2, but
% Debian 9 (Stretch) does not have revtex4-2
\usepackage[utf8]{inputenc}
% \usepackage{epstopdf} % Converts .eps to
% \epstopdfsetup{update}
\usepackage[caption=false]{subfig}
\usepackage{hyperref} % Required for customising links and the PDF*
\linespread{1.05}
\hypersetup{pdfpagemode={UseOutlines},
bookmarksopen=true,
bookmarksopenlevel=0,
hypertexnames=false,
colorlinks=true, % Set to false to disable coloring links
citecolor=blue, % The color of citations
linkcolor=red, % The color of references to document elements (sections, figures, etc)
urlcolor=black, % The color of hyperlinks (URLs)
pdfstartview={FitV},
unicode,
breaklinks=true,
}
\usepackage{graphicx,amssymb,amsmath}
% \usepackage[english]{babel}
\usepackage{grffile}
\usepackage{color}
\usepackage{array}
\usepackage{hhline}
\usepackage[]{algorithm2e}
\usepackage{ulem}
% \usepackage{float}
\usepackage{siunitx}
\sisetup{
inter-unit-product = \ensuremath{{}\!\cdot\!{}},
detect-all,
separate-uncertainty = true,
exponent-product = \times,
space-before-unit = true,
output-decimal-marker = {,},
multi-part-units = brackets,
range-phrase = --,
% allow-number-unit-breaks,
list-final-separator = { et },
list-pair-separator = { et },
abbreviations
}
\setlength{\tabcolsep}{7pt}
\usepackage{booktabs}
\usepackage{float}
\newcommand{\cor}[1]{\textcolor{red}{#1}}
\newcommand{\todo}[1]{\textcolor{red}{TODO: #1}}
\newlength{\figwidth}
\setlength{\figwidth}{120mm}
% \setlength{\figwidth}{0.7\textwidth} % useful in single column
\newcommand{\R}{\mathcal{R}}
\newcommand{\eps}{\varepsilon}
\newcommand{\epsK}{{\varepsilon_{\!\scriptscriptstyle kin}}}
\newcommand{\epsKK}{{\varepsilon_{\!\scriptscriptstyle K 2}}}
\newcommand{\epsKKKK}{{\varepsilon_{\!\scriptscriptstyle K 4}}}
\newcommand{\epsA}{{\varepsilon_{\!\scriptscriptstyle pot}}}
\newcommand{\xx}{\boldsymbol{x}}
\newcommand{\rr}{\boldsymbol{r}}
\newcommand{\kk}{\boldsymbol{k}}
\newcommand{\eek}{\boldsymbol{e}_{\boldsymbol{k}}}
\newcommand{\eeh}{\boldsymbol{e}_h}
\newcommand{\eep}{\boldsymbol{e}_{p\kk}}
\newcommand{\eet}{\boldsymbol{e}_{t\kk}}
\newcommand{\eetheta}{\boldsymbol{e}_\theta}
\newcommand{\eex}{\boldsymbol{e}_x}
\newcommand{\eey}{\boldsymbol{e}_y}
\newcommand{\eez}{\boldsymbol{e}_z}
\newcommand{\cc}{\boldsymbol{c}}
\newcommand{\uu}{\boldsymbol{u}}
\newcommand{\vv}{\boldsymbol{v}}
\newcommand{\vvk}{\hat{\boldsymbol{v}}_{\kk}}
\newcommand{\vvp}{\hat{v}_{p\kk}}
\newcommand{\vvt}{\hat{v}_{t\kk}}
\newcommand{\vvs}{\hat{\boldsymbol{v}}_{s\kk}}
\newcommand{\bbk}{\hat{b}_{\kk}}
\newcommand{\ff}{\boldsymbol{f}}
\newcommand{\ffk}{\boldsymbol{\hat{f}}_{\kk}}
\newcommand{\bomega}{\boldsymbol{\omega}}
\newcommand{\bnabla}{\boldsymbol{\nabla}}
\newcommand{\Dt}{\mbox{D}_t}
\newcommand{\p}{\partial}
\newcommand{\mean}[1]{\langle #1 \rangle}
\newcommand{\epsP}{\varepsilon_{\!\scriptscriptstyle P}}
\newcommand{\epsm}{\varepsilon_{\!\scriptscriptstyle m}}
\newcommand{\CKA}{C_{K\rightarrow A}}
\newcommand{\D}{\mbox{D}}
\newcommand{\diff}{\text{d}}
\newcommand{\bv}{Brunt-V\"ais\"al\"a }
\newcommand{\kmax}{k_{\max}}
\newcommand{\thk}{\theta_{\kk}}
\newcommand{\phk}{\varphi_{\kk}}
\newcommand{\thf}{\theta_f}
\newcommand{\ok}{\omega_{\kk}}
\setlength\parindent{0pt}
\begin{document}
\title{Weak internal gravity waves in stratified flows with and without vortical modes}
\author{Vincent Labarre}
\email[]{vincent.labarre@oca.eu}
\affiliation{Universit\'{e} C\^{o}te d'Azur, Observatoire de la C\^{o}te d'Azur, CNRS,
Laboratoire Lagrange, Nice, France.}
\author{Pierre Augier}
\email[]{pierre.augier@univ-grenoble-alpes.fr}
\affiliation{Laboratoire des Ecoulements G\'eophysiques et Industriels, Universit\'e
Grenoble Alpes, CNRS, Grenoble-INP, F-38000 Grenoble, France}
\author{Giorgio Krstulovic}
\email[]{giorgio.krstulovic@oca.eu}
\affiliation{Universit\'{e} C\^{o}te d'Azur, Observatoire de la C\^{o}te d'Azur, CNRS,
Laboratoire Lagrange, Nice, France.}
\author{Sergey Nazarenko}
\email[]{sergey.nazarenko@unice.fr}
\affiliation{Universit\'{e} C\^{o}te d'Azur, CNRS, Institut de Physique de Nice -
INPHYNI, Nice, France}
\begin{abstract}
In order to test weak wave turbulence, we perform direct numerical simulations of
stratified turbulence without shear modes, and with or without vortical modes. Removing
vortical modes naturally helps to have a better overall balance between poloidal
kinetic energy, involved in internal gravity waves, and potential energy. A
spatio-temporal analysis reveals that removing vortical modes is not enough to
concentrate the energy in temporal scales, meaning that removing vortical modes does
not allow to get that closer to an internal gravity wave turbulence regime. For highly
stratified simulations, waves dominate in a sub-region of the buoyancy range defined by
the non-linearity parameter, as observed in Yokoyama and Takaoka (2019)
\cite{yokoyama_energy-based_2019}. However, conversion between kinetic and potential
energy in that region are not consistent with a wave turbulence regime for all
simulations. We give simple dimensional argument to stress the importance of the aspect
ratio for weak non-linearity to hold. For simulations with aspect ratio one, when the
waves dominated region starts to enclose the forcing region, we observe a transition to
a new regime that shares characteristics with internal gravity wave turbulence.
\todo{Target journal: JFM?}
\end{abstract}
%----------------------------------------------------------------------------------------
% Print the title
\maketitle
%----------------------------------------------------------------------------------------
% ARTICLE CONTENTS
%----------------------------------------------------------------------------------------
\section{Introduction}
\label{sec:introduction}
\input{introduction.tex}
\section{Methods}
\label{sec:methods}
\input{methods.tex}
\section{Results}
\label{sec:results}
\input{results.tex}
\section{Conclusions}
\label{sec:conclusions}
\input{conclusions.tex}
\begin{acknowledgments}
This project was supported by the Simons Foundation through the Simons collaboration on
wave turbulence. Part of the computations have been done on the ``Mesocentre SIGAMM''
machine, hosted by Observatoire de la Cote d'Azur. The authors are grateful to the OPAL
infrastructure from Université Côte d’Azur and the Université Côte d’Azur’s Center for
High-Performance Computing for providing resources and support. This work was granted
access to the HPC/AI resources of IDRIS under the allocation 2022-A0122A13417 made by
GENCI.
\end{acknowledgments}
\appendix
\section{List of simulations}
\input{../tmp/table_better_simuls.tex}
\input{../tmp/table_better_simuls_proj.tex}
\input{../tmp/table_better_simuls_ratio_one.tex}
\input{../tmp/table_better_simuls_proj_ratio_one.tex}
\newpage
%\appendix\section{A great appendix}
%\label{appendix}
\bibliography{biblio}
\end{document}
This diff is collapsed.
In this study, we performed direct numerical simulations of stratified turbulence
without shear modes, and with or without vortical modes. We observed that removing
vortical modes helps to have a better overall balance between poloidal kinetic energy
and potential energy. However, a spatiotemporal analysis showed that removing vortical
mode does not help to concentrate energy in temporal scales. Localisation of energy in
temporal scales is rather controlled by the stratification, i.e. by the Froude number.
Following \cite{yokoyama_energy-based_2019}, we performed a spectral energy budget in
the two-dimensional spectral space. As in \cite{yokoyama_energy-based_2019}, we
observed that the wave dominated region is delimited by the non-linearity parameter
$\chi_{\kk}$, with a threshold $\chi_{\kk} = 1/3$. We also observe a dependence in
$k/k_b$ (where $k_b$ is the buoyancy wave-vector) in our simulations such that the wave
dominated region should lies in $\left\{\kk ~ | ~ \chi_{\kk} <1/3, k < k_b \right\}$.
Consequently, performing numerical simulation in the buoyancy range (with $k_b \eta
\leq 1$) is tempting for observing wave turbulence at a lower numerical cost. However,
the conversion between kinetic energy and potential energy do not necessarily shows
fluctuations around zero in the wave dominated regime, as we expect for a system of
statistically stationary waves.
With a simple dimensional analysis, we showed that the aspect ratio of the numerical
domain (or of the container) is particularly important for the weak non-linearity
assumption to hold. We therefore conducted simulations with aspect ratio one. These
simulations revealed that a transition occurs when the forcing region starts to lie in
the waves dominated region: at this point strong fluctuations of the conversion between
kinetic energy and potential energy appears in the all wave dominated range, which
indicates the presence of statistically stationary waves. This transition has a simple
physical meanings: (i) When the forcing lies in the waves dominated region, injected
energy is efficiently converted into waves. (ii) The forcing does not feed other modes
of motions that could cascade and perturbs waves over the all scales. (iii) We start to
observe a temporal inertial range for waves. It is expected that the same transition
could occurs in others WWT systems, including MHD and rotating flows.
For highly stratified simulations, the spatial spectra of the equipartition energy $E_e
= 2 \min(E_{polo}, E_A)$ deviates from the prediction of waves turbulence theory
\cite{dematteis_downscale_2021, dematteis_origins_2022}, and the Garrett-Munk 1976
model \cite{garrett_internal_1979} where $k_z \leq k_h \leq k_b$. It is interesting to
note that the waves turbulence prediction is derived in the limit $k_z \gg k_h$, which
is the opposite limit that our numerical results, and the
\cite{yokoyama_energy-based_2019}'s ones suggest: the non-linear broadening is shown to
be large for $k_z \gg k_h$, therefore breaking a WWT theory's assumption in this region
of the spectral space. For more, the energy spectra are not separable, contrary to what
it is usually assumed in WWT. This motivate further studies to remove the $k_z \ll k_h$
assumption in the derivation of the kinetic equation of internal gravity waves.
2022strat_polo_waves/input/fig_poloidal-toroidal.png

72.2 KiB

% WWT
As eddies in classical hydrodynamic turbulence, waves in nonlinear systems interact and
transfer conserved quantities along scales in a cascade process. The Weak-Wave
Turbulence (WWT) theory aims to provides a statistical description of the system when
the non-linearity is small \cite{zakharov_kolmogorov_1992,nazarenko_wave_2011,
nazarenko_wave_2015}. The applications of this theory encompass capillary-gravity waves
\cite{falcon_experiments_2022}, Alfvén waves in magnetohydrodynamics (MHD)
\cite{galtier_weak_2000}, inertial and internal waves in rotating stratified fluids
\cite{caillol_kinetic_2000, galtier_weak_2003, medvedev_turbulence_2007}, Kelvin waves
in superfluids \cite{lvov_weak_2010}, elastic plates \cite{during_weak_2006},
gravitational waves \cite{galtier_turbulence_2017}, and density waves in Bose-Einstein
condensates \cite{dyachenko_optical_1992}.
% Weak non-linearity and small non-linear broadening
Two key hypotheses are used in WWT. The first one is weak non-linearity of the
dynamical equations. This condition is often translated in terms of the separation of
spatial scales: the considered scale on one side and the saturation scale on the other
side. The saturation scale is here defined as the scale at which the weak non-linearity
does not holds anymore, inducing wave breaking. This hypothesis is quite similar to the
separation between the integral scale and the dissipative scale in the more famous
three-dimensional hydrodynamic turbulence. The second hypothesis corresponds to a time
scales separation, and is more specific to wave turbulence. It consists to require that
the linear time, given by the wave period $\tau_{L}$, is small compared to the
non-linear time of interactions between waves $\tau_{NL}$ \cite{nazarenko_wave_2011}.
For many physical systems, at least one the two hypotheses is broken in some range of
scales, which reduce or even break the validity of WWT \cite{biven_breakdown_2001}.
Yet, when scales separations are observed in space and time, a weak wave turbulent
range can appear in spatiotemporal spectra of waves' energy. Besides experimental and
numerical difficulties for obtaining scales' separations, testing WWT in isotropic
systems is conceptually simpler. For this reason, many progress have been made in the
experimental and numerical verification of WWT for elastic plates
\cite{miquel_nonstationary_2011, yokoyama_identification_2014}, capillary-gravity waves
\cite{falcon_experiments_2022}, Gross-Pitaevskii equation \cite{zhu_testing_2022}, or
2D acoustic waves \cite{griffin_energy_2022}.
Anisotropic turbulence is, generally speaking, more difficult to investigate due to the
reduce number of symmetries. For example, it has been shown that studying
two-dimensional spatial spectra instead of one-dimensional integrated spectra is
essential to investigate stratified turbulence \cite{yokoyama_energy-based_2019}. Then,
the anisotropy makes the problem multidimensional in Fourier space, which makes the
notion of spectral energy fluxes more difficult than in isotropic turbulence
\cite{yokoyama_energy-flux_2021}. An additional difficultly is that scales separations
required by WWT can also depend on anisotropy. As an example, for linearly stratified
flows, separations of scales require
\begin{equation}
k/k_b \ll 1 ~~~~ \text{and} ~~~~ \tau_{L}/\tau_{NL} = \frac{(\epsK k^2)^{1/3}}{N k_h/k} \ll 1,
\end{equation}
where $k_b = N/U_h$ is the buoyancy wave-vector, $U_h$ is the rms of the horizontal
velocity, $N$ is the \bv frequency, $\epsK$ is the kinetic energy dissipation rate, $k$
is the wave-vector modulus, and $k_h$ is the horizontal wave-vector modulus. Then, due
to the anisotropic dispersion relation of internal gravity waves, time scales
separation will be less valid for small $k_h$, and even impossible for modes with
$k_h=0$. Consequently, separation of times scales can be violated even if a large
separation of spatial scales is observed (i.e. $k/k_b \ll 1$). It is therefore more
difficult to observe signatures of wave turbulence in anisotropic systems like MHD, or
rotating and stratified flows. \\
Stratified flows are not only an interesting conceptual problem, but they are also
important for the understanding of geophysical flows \cite{staquet_internal_2002,
vallis_atmospheric_2017}. In particular, understanding the role of waves in mixing is
required for sub-grid parameterizations in climate models \cite{mackinnon_climate_2017,
gregg_mixing_2018}. It is therefore not surprising that stratified flows received a
particular attention from both the ``strong'' turbulence community
\cite{billant_self-similarity_2001, waite_stratified_2004, waite_stratified_2006,
lindborg_energy_2006, brethouwer_scaling_2007, waite_stratified_2011,
kimura_energy_2012, bartello_sensitivity_2013, brunner-suzuki_upscale_2014,
augier_stratified_2015, maffioli_vertical_2017}, and the WWT community
\cite{caillol_kinetic_2000, lvov_hamiltonian_2001, lvov_weak_2010,
dematteis_downscale_2021, dematteis_origins_2022}. It turns out that internal waves are
effectively important in the dynamics and mixing in stratified flows
\cite{maffioli_signature_2020, lam_partitioning_2020, lam_energy_2021}, and that,
three-wave resonant interactions are responsible for slow, net energy transfers between
different wave-numbers \cite{davis_succession_2020, rodda_experimental_2022}. However,
many questions and issues remain about the applicability of WWT to stratified flows.
A recent study, \cite{yokoyama_energy-based_2019}, showed that a wave dominated region
should lie in the spectral region defined by $\tau_{L} / \tau_{NL} < 1/3$, in agreement
with observations made for MHD \cite{meyrand_direct_2016}. Yet, non-waves structures
like shear modes (purely vertical shear) and vortical modes (vertical vorticity) still
make the observation of a system of weakly interacting internal gravity waves
impossible. The same problem appear in rotating flows, in which the geostrophic modes
plays the role of the non-propagative structure. Experimentalists bypassed this
difficulty to observed weak inertial wave turbulence in a rotating tank, by using
honeycomb grids at top and bottom boundaries to dissipate geostrophic modes efficiently
\cite{brunet_shortcut_2020, monsalve_quantitative_2020}. These works point out that
some wave systems have a tendency to generate non-wave motions that can severely affect
the wave dynamics and should be suppressed in experiments aiming to observe wave
turbulence. In rotating inertial waves these modes are geostrophic flows, while in
internal gravity waves these modes are vertical shear and vertical vorticity
\cite{nazarenko_verifying_2020}. \\
The present study first deals with the existence and properties of a weak internal
gravity waves turbulence regime. To this end, we performed numerical simulations of
stratified turbulence at various stratifications. In the same spirit than experiments
in a rotating tank presented in \cite{brunet_shortcut_2020,
monsalve_quantitative_2020}, we removed shear modes in all of our simulations and the
vortical modes for half of them. The removing of vortical modes is done using the
Craya-Herring spectral decomposition \cite{craya_contribution_1957}. Consequently,
contrary to the removing of geostrophic modes in rotating tank, removing vortical modes
in a stratified flows is hardly doable experimentally. Note that, in order to get rid
of non waves structures, numerical simulations reduced dynamical equations (i.e.
without non-waves structures) of stratified rotating flows in the hydrostatic balance
where already done \cite{lvov_nonlinear_2009}. Despite the simplifications, these
simulations reproduced some key features of oceanic internal-wave spectra: accumulation
of energy at near-inertial waves and realistic frequency and horizontal wave-number
dependencies of spatiotemporal spectra. In the present work, we don't account for
rotation, and don't use the hydrostatic balance approximation. Also, having twins
simulations (i.e. with or without vortical modes) allows to have a better comparison to
understand the role of the vortical modes on the dynamics of stratified flows. In
particular, our study shows that the existence of vortical modes is not the limiting
factor when trying to observe internal gravity wave turbulence. In fact, the
stratification strength and the position of the forcing in spectral space appears to be
more important in that context.
The most commonly used WWT prediction is certainly the scale invariant stationary
solutions to the kinetic equation, which gives the expected spatial waves' energy
spectra in the statistically steady state. The solutions are usually of two types: the
thermodynamic equilibrium solution (Rayleigh-Jeans spectra), and the non equilibrium
solution(s) that are linked to the cascade(s) of the dynamical invariant(s) of the
system among scales (Kolmogorov-Zakharov spectra) \cite{nazarenko_wave_2011}. It is
important to note that these solutions can be considered as valid only if the collision
integral in the waves' kinetic equation converges. The importance of the nonlinear
interactions among internal gravity waves was recognized early, leading to several
derivations of waves' kinetic equations (see \cite{muller_nonlinear_1986,
lvov_resonant_2012}). Derivation of the Kolmogorov-Zakharov in the limit of large
vertical wave-number $k_z \simeq k \gg k_h$ can be found in \cite{caillol_kinetic_2000,
lvov_hamiltonian_2001}. This solution corresponds to an energy cascade, and is given by
\begin{equation}
n(k_h, k_z) \sim k_h^{-7/2} k_z^{-1/2} ~~~~ \text{and} ~~~~ E(k_h, k_z) \sim k_h^{-5/2} k_z^{-3/2}
\end{equation}
respectively for the wave-action spectra and the energy spectra. Yet, it was noted that
this solution comes through a cancellation between oppositely signed divergent
contributions in their respective collision integrals. Therefore, the existence of this
solution is fortuitous. Later, it was shown that power law solutions $n(k_h,k_z) \sim
k_h^{-\alpha_h} ~ k_z^{-\alpha_z}$ have convergent collision integral's contributions
only on the segment $\alpha_h \in ]3, 4[, \alpha_z=0$ \cite{lvov_oceanic_2010,
dematteis_downscale_2021}. The collision integral was then computed numerically for
$\alpha_z=0$, and it was deduced that the only scale invariant stationary solution to
the kinetic equation (zero of the collision integral) was close to
\begin{equation}
n(k_h, k_z) \sim k_h^{-3.69} k_z^{0} ~~~~ \text{and} ~~~~ E(k_h, k_z) \sim k_h^{-2.69} k_z^{-1}.
\end{equation}
It was also shown that, on the segment $\alpha_z =0$, the dominant contributions to the
collision integral corresponds to non-local transfers first identified by McComas
\cite{McComas_resonant_1977}, and by interactions between wave-vectors that are almost
collinear in the horizontal plane \cite{dematteis_downscale_2021}.
The manuscript is organized as follows. In section \ref{sec:methods}, we present our
methodology including a presentation of the code and the simulations. Our results are
presented in section \ref{sec:results}. Subsection \ref{subsec:global} is devoted to
the study of global energy ratio in the control parameters space for simulations at
aspect ratios less than unity. It allows us to observe that WWT is more likely to
occurs at high stratification and without vortical modes. Subsections
\ref{subsec:khkz} and \ref{subsec:khkzomega} deal with the spatiotemporal analysis of a
couple of strongly stratified simulations in order to investigate further the presence
of linear waves in spatial scales. It turns out that aspect small aspect ratio prevents
a WWT regime to occur. For this reason, we performed additional numerical simulations
with aspect ratio one. This new simulations are presented in subsection
\ref{subsec:ratio-one}. We give concluding remarks in section \ref{sec:conclusions}.
% Basics
We start from the three-dimensional Navier-Stokes equations under the Boussinesq
approximation:
\begin{align}
\label{eq:Continuity}
\bnabla \cdot \vv &= 0 \\
\label{eq:Impulsion}
\p_t\vv + (\vv \cdot \bnabla)\vv &= b\boldsymbol{e}_z - \bnabla p +
\nu \nabla^2\vv + \ff, \\
\label{eq:Buoyancy}
\p_t{b} + (\vv \cdot \bnabla)b &= -N^2 v_z + \kappa\nabla^2{b},
\end{align}
where $\rr = (x,y,z)$ represents the three spatial coordinates in the cartesian frame
$(O, \eex, \eey, \eez)$, $\eez$ is the stratification axis, $\vv=(v_x, v_y, v_z)$ is
the velocity, $b$ the buoyancy, $p$ the total kinematic pressure, $N$ the \bv
frequency, $\nu$ is the viscosity, $\kappa$ is the diffusivity, and $\ff$ is the
velocity forcing. We fix the Schmidt number $Sc = \nu/\kappa$ to one. We consider a
periodic domain of horizontal size $L_x = L_y = L_h = 3$. The vertical size of the
domain $L_z$, is varied. We note $(n_x, n_y, n_z)$ the numbers of collocations points
in the three spatial directions, with $n_x = n_x = n_h$.
The Fourier transform of the velocity field $\vvk = (\hat{v}_{x\kk}, \hat{v}_{y\kk},
\hat{v}_{z\kk})$ can be written using the poloidal-toroidal-shear decomposition (see
e.g. \cite{craya_contribution_1957, maffioli_vertical_2017})
\begin{equation}
\vvk = \begin{cases} \vvp \eep + \vvt \eet ~~~~ \text{if} ~ k_h \neq 0, \\
\vvs = \hat{v}_{x \kk} \eex + \hat{v}_{y \kk} \eey ~~~~ \text{if} ~ k_h = 0.
\end{cases}
\end{equation}
where
\begin{equation}
\label{eq:poloidal-toroidal}
\eek = \frac{\kk}{k}, ~~~~
\eep = \frac{\kk \times (\kk \times \eez)}{|\kk \times (\kk \times \eez)|}, ~~~~
\eet = \frac{\kk \times \eez}{|\kk \times \eez|},
\end{equation}
$\vvp$ is the poloidal velocity, $\vvt$ is the toroidal velocity, $\vvs$ is the shear
modes velocity, $\kk=(k_x, k_y, k_z)$ denotes the wave-vector, $k =|\kk| = \sqrt{k_x^2
+ k_y^2 + k_z^2}$ is its norm, and $k_h = \sqrt{k_x^2 + k_y^2}$ is the norm of the
horizontal wave-vector (Figure \ref{fig:poloidal-toroidal}).
\begin{figure}[H]
\centering
\includegraphics[width=0.5\textwidth]{fig_poloidal-toroidal}
\caption{Illustration of the poloidal-toroidal basis $(\eek, \eep, \eet)$ defined by
equations (\ref{eq:poloidal-toroidal}). $\thk$ is the angle between $\eez$ and $\eek$.
$\phk$ is the angle between the horizontal projection of $\kk$ and $\eex$.
\label{fig:poloidal-toroidal}}
\end{figure}
% Description of linear modes and shear modes dynamics
In spectral space, the equations of motion then reads
\begin{align}
\label{eq:StratifiedSpectralPoloidalToroidal}
\begin{cases}
\dot{\hat{v}}_{p\kk} &= - (\widehat{\vv \cdot \bnabla \vv})_{\kk} \cdot \eep - \bbk \sin \thk - \nu k^2 \vvp + \hat{f}_{\kk}, \\
\dot{\hat{v}}_{t\kk} &= - (\widehat{\vv \cdot \bnabla \vv})_{\kk} \cdot \eet - \nu k^2 \vvt, \\
\dot{\hat{b}}_{\kk} &= - (\widehat{\vv \cdot \bnabla b})_{\kk} + N^2 \vvp \sin \thk - \kappa k^2 \bbk
\end{cases}
\end{align}
if $k_h \neq 0$, and
\begin{align}
\label{eq:StratifiedSpectralShearModes}
\begin{cases}
\dot{\hat{\vv}}_{s\kk} &= - (\widehat{\vv \cdot \bnabla \vv_h})_{\kk} - \nu k^2 \vvs, \\
\dot{\hat{b}}_{\kk} &= - (\widehat{\vv \cdot \bnabla b})_{\kk} - \kappa k^2 \bbk
\end{cases}
\end{align}
if $k_h=0$, where $\vv_h = (v_x,v_y,0)$ is the horizontal component of $\vv$. $k_\alpha
\in \delta k_\alpha \mathbb{Z}$ with $\delta k_\alpha = 2\pi /L_\alpha$ for $\alpha
=x,y,z$. We will also note $\delta k_x = \delta k_y = \delta k_h$.
$\left(\widehat{\cdot}\right)_{\kk}$ denotes the Fourier transform. In the
poloidal-toroidal-shear decomposition, it is clear that linear waves live in the
poloidal poloidal velocity and the buoyancy, but not in the toroidal velocity.
The linear, unforced, inviscid dynamics can be written as
\begin{align}
\label{eq:StratifiedSpectralPoloidalToroidalLinear}
&\dot{a}_{\kk}^{(-)} = -i \ok a_{\kk}^{(-)}, ~~~~ \dot{a}_{\kk}^{(0)} = 0, ~~~~ \dot{a}_{\kk}^{(+)} = +i \ok a_{\kk}^{(+)} ~~~~ \text{if} ~~ k_h \neq 0 \\
\label{eq:StratifiedSpectralShearModesLinear}
&\dot{\hat{\vv}}_{s\kk} = \mathbf{0}, ~~~~ \dot{\hat{b}}_{\kk} = 0 ~~~~ \text{if} ~~ k_h = 0.
\end{align}
where $a_{\kk}^{(\mp)} = \left(\vvp \mp i \frac{\bbk}{N} \right)/\sqrt{2\ok}$ are the
waves modes,
\begin{equation}
\ok =N k_h/k = N \sin \thk
\end{equation}
is the pulsation of the waves, $\thk$ is the angle between $\kk$ and the stratification
axis $\eez$, and $a_{\kk}^{(0)} = \vvt$ are the vortical modes (vertical vorticity).
Equations
(\ref{eq:StratifiedSpectralPoloidalToroidalLinear}-\ref{eq:StratifiedSpectralShearModesLinear})
show that both shear modes and vortical modes have zero frequency.
% Forcing, truncations, and hyperviscosity, code
Since we are motivated by forcing internal gravity waves, in which only the poloidal
part is involved and whose dispersion relation is anisotropic, we choose a time
correlated, anisotropic, poloidal velocity forcing $\ffk = \hat{f}_{\kk} \eep$. The
flow is forced at large spatial scales and small angle $\left\{ \kk ~ | ~ 5 \leq
k/\delta k_h \leq 20, ~ |\sin \thk - 0.3| \leq 0.05 \right\}$, meaning that relatively
slow internal waves are forced. The correlation time of the forcing is equal to the
period of the forced waves $T_c = 2\pi / 0.3 N$. It is computed in spectral space such
that the kinetic energy injection rate $P_K$ is constant and equal to unity. The
forcing scheme is described by the Algorithm \ref{algorithm:forcing}. \\
\begin{algorithm}[H]
$t=0$: generate two random numbers $f_{0\kk}, f_{1\kk} \sim \mathcal{P}$ $\forall \kk$ in the forcing region; $t_0 = 0$; \\
\While{$t \leq T$}{
$t \rightarrow t +\delta t$; \\
\If{$t - t_0 \geq T_c$}{
$t_0 \rightarrow t$; \\
$f_{0 \kk} \rightarrow f_{1 \kk}$; \\
genererate $f_{1\kk} \sim \mathcal{P}$
}{}
$\ffk = \left\{ f_{0\kk} - \dfrac{(f_{1\kk} - f_{0\kk})}{2} \left[ \cos\left( \dfrac{\pi(t-t_0)}{T_c} \right) + 1 \right] \right\} ~ \eep$ \\
$\ffk$ are normalized to ensure $P_K = 1$
}
\caption{Forcing scheme. $\mathcal{P}$ is the probability law which has for distribution $p_X(x) = \left\{
\frac{1}{4} ~~ \text{if} ~ Re(x) \in [-1:1] ~ \text{and} ~ Im(x) \in [-1:1], 0 ~~ \text{otherwise} \right\}$. $\delta t$ is the time increment at each time step. \label{algorithm:forcing}}
\end{algorithm}
Some modes are removed in all the present simulations because they cause numerical and
physical problems and/or are not consistent with experiments in which the flow is
bounded with walls. (i) All modes with wavenumber modulus larger than $\kmax = 0.8
(n_h/2) \delta k_h$ are truncated to limit aliasing. (ii) Shear modes are truncated.
(iii) Vertically invariant vertical velocity (internal waves at $\omega = N$) is also
forbidden. We add hyperviscosity and hyperdiffusivity terms of order $4$, $\nu_4 =
\kappa_4$, in many simulations in order to avoid numerical issues. They are fixed
according to the numerical resolution such that the hyperviscous and hyperdifussive
dissipations act at the smallest scales. For most of the simulations, the resolution
increased sufficiently for the hyperviscosity and hypediffusivity to have no
significant impact on measured global quantities, and on spectra in the inertial range.
We measure the turbulent kinetic dissipations $\epsKK$ and $\epsKKKK$ based on both
viscosities, and the total energy dissipation $\epsK = \epsKK + \epsKKKK$. The product
of the maximal wave-vector $\kmax$ with the Kolmogorov scales $\eta \equiv (\nu^3 /
\epsKK)^{1/4}$ and $\eta_4 \equiv (\nu_4^3/\epsK)^{1/10}$ are computed to estimate que
quality of the simulations. Only some simulations at high $N$ remain affected by
hyperdiffusivity. The time advancement is performed using the $4^{th}$ order
Runge-Kutta scheme. The numerical simulations are performed using the pseudo-spectral
solver \texttt{ns3d.strat} from the FluidSim Python package.
\cite{augier_fluiddyn_2019, mohanan_fluidfft_2019, mohanan_fluidsim_2019}.
% Poloidal projection
For this study, we performed two types of simulations. In the first ones, we keep
vortical modes. For each simulation with vortical modes, we made a ``twin'' simulation
with the same values of input parameters but without vortical modes. The goal of this
procedure is to remove two non propagative structures, namely the shear modes and the
vortical modes in the second dataset in order to get closer to a system on internal
gravity waves. The reader is refered to \todo{cite manuscript for poloidal projection}
for more details about the numerical simulations with vortical modes. The simulations
without vortical modes are realised in the very same way.
% Dimensionless numbers and useful quantities
The main physical input parameters are the \bv frequency and the viscosity. Since both
forcing length and energy injection rate are in practice equal or close to unity, we
can take as input parameters the \bv frequency and an input buoyancy Reynolds number
$\R_i = 1/ \nu N^2$. The turbulent non-dimensional numbers are the horizontal turbulent
Froude number and the buoyancy Reynolds numbers that are respectively
\begin{equation}
F_h = \epsK / ({U_h}^2 N), ~~~~ \R = \epsK / (\nu N^2), ~~~~ \R_4 = \frac{\epsK U_h^2}{\nu_4 N^4},
\end{equation}
where $\epsK$ is the average kinetic energy dissipation rate, and $U_h$ is the rms of
the horizontal velocity. We also compute the buoyancy and the Ozmidov wave-vectors $k_b
\equiv N/U_h$ and $k_O \equiv \sqrt{N^3 / \epsK}$. All the quantities presented in this
manuscript are computed from averaging over 2 unit times of the simulations, when the
flow is stationary, at the end of the simulations. A list of the simulations we used in
subsection \ref{subsec:global}, with relevant parameters and physical quantities is
given by tables \ref{table-better-simuls} and \ref{table-better-simuls-proj}. \\
% Simulations with aspect ratio one
For reasons that will be explained in subsection \ref{subsec:ratio-one}, we also
performed additional simulations with aspect ratio one, $\nu=0$, a forcing at large
spatial scales but large angle $\left\{ \kk ~ | ~ 3 \leq k/\delta k_h \leq 3.5, ~ |\thk
- 0.8\pi/2| \leq \pi/40 \right\}$, meaning that relatively fast internal waves are
forced. The correlation time of the forcing is equal to the period of the forced waves
$T_c = 2\pi / \sin(0.8\pi/2) N$. We also changed the de-aliasing coefficient to $2/3$
such that $\kmax = n_h \delta k_h/3$. A list of these simulations is given byt tables
\ref{table-better-simuls-ratio-one} and \ref{table-better-simuls-proj-ratio-one}.
This diff is collapsed.
This diff is collapsed.
import sys
from runpy import run_path
import matplotlib.pyplot as plt
print("import util")
from util import (
has_to_be_made,
has_to_save,
here,
paths_all,
paths_all_proj,
paths_all_ratio_one,
)
def make_fig(name, source=None):
if source is None:
source = f"save_{name}.py"
fig_name = f"fig_{name}.png"
if has_to_be_made(fig_name, source):
print(f"Calling {source} to make {fig_name}")
run_path(str(here / source))
def make_table(name, source=None):
if source is None:
source = f"save_table_{name}.py"
if has_to_be_made(f"table_{name}.tex", source):
print(f"Calling {source} to make {name}")
run_path(str(here / source))
print(paths_all_ratio_one)
make_table("better_simuls")
make_table("better_simuls_ratio_one")
make_fig("ratio_E_vs_Fh")
sys.argv.append("")
sys.argv[-1] = 0
make_fig("seb_regimes_L_projNone", "save_seb_regimes.py")
del sys.argv[-1]
sys.argv.append("")
sys.argv[-1] = 1
make_fig("seb_regimes_L_projpoloidal", "save_seb_regimes.py")
del sys.argv[-1]
make_fig("nonlinear_regimes_L", "save_nonlinear_regimes.py")
make_fig(
"spectra_slices_omega_kh_kz_regimes_L",
"save_spectra_slices_omega_kh_kz_regimes.py",
)
# make_fig("delta_omega_spectra")
make_fig("ratio_E_vs_Fh_ratio_one")
make_fig("seb_transition_ratio_one")
make_fig("nonlinear_ratio_one")
make_fig("spectra_kh_kz_ratio_one")
make_fig("spectra_kh_kz_ratio_one_2D")
if not has_to_save:
plt.show()
import h5py
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 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,
get_path_finer_resol,
save_fig,
)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
cm = cm.get_cmap("inferno", 100)
def plot_nlb(sim, ax):
t_start, t_last = times_start_last_from_path(path)
tmin = t_last - 2
kf, kb, ko, keta, kd = compute_kf_kb_ko_keta_kd(sim, tmin)
mean_values = sim.output.get_mean_values(tmin=tmin, customize=customize)
Fh = mean_values["Fh"]
epsK = mean_values["epsK"]
params = load_params_simul(path)
nh = nx = params.oper.nx
proj = params.projection
N = sim.params.N
path_spec = sorted(path.glob(f"spatiotemporal/periodogram_[0-9]*.h5"))
assert len(path_spec) == 1, f"Not only 1 periodogram in {path} \n"
path_spec = path_spec[0]
with h5py.File(path_spec, "r") as f:
kh = f["kh_spectra"][:]
kz = f["kz_spectra"][:]
delta_kh = kh[1]
delta_kz = kz[1]
omegas = f["omegas"][:]
EA = f["spectrum_A"][:]
EKz = f["spectrum_K"][:] - f["spectrum_Khd"][:] - f["spectrum_Khr"][:]
Epolo = f["spectrum_Khd"][:] + EKz
Etoro = f["spectrum_Khr"][:]
E = Epolo + Etoro + EA
Ee = 2 * np.minimum(EA, Epolo)
Ed = EA + Epolo - Ee
spectrum = Epolo + EA
omega_emp, delta_omega_emp = compute_omega_emp_vs_kzkh(
N, spectrum, kh, kz, omegas
)
KH, KZ = np.meshgrid(kh, kz)
K = (KH**2 + KZ**2) ** 0.5
K_NOZERO = K.copy()
K_NOZERO[K_NOZERO == 0] = 1e-16
omega_disp = N * KH / K_NOZERO
chi = (K**2 * epsK) ** (1 / 3) / omega_disp
cs = ax.pcolormesh(
kh,
kz,
np.log10(delta_omega_emp / omega_disp),
cmap=cm,
vmin=-0.5,
vmax=1.5,
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, 2 * max(kh) / 3])
ax.set_ylim([delta_kh, 2 * max(kh) / 3])
# ax.set_xscale("lin")
# ax.set_yscale("lin")
return cs
Ns = [10, 20, 80]
nbax = 0
css = [None for i in range(6)]
fig, axes = plt.subplots(
ncols=2, nrows=3, figsize=(10, 3 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
ax4 = axes[2, 0]
ax5 = axes[2, 1]
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)
css[nbax] = plot_nlb(sim, axs[nbax])
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_xticks([])
for ax in [ax1, ax3, ax5]:
ax.set_yticks([])
ax0.set_title(r"Standard Navier-Stokes" + "\n" + r"$(a)$", fontsize=16)
ax1.set_title(r"Without vortical modes" + "\n" + r"$(b)$", fontsize=16)
ax2.set_title(r"$(c)$", fontsize=16)
ax3.set_title(r"$(d)$", fontsize=16)
ax4.set_title(r"$(e)$", fontsize=16)
ax5.set_title(r"$(f)$", fontsize=16)
fig.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[3], cax=cbar_ax)
cbar.set_ticks([-0.5, 0.0, 0.5, 1.0, 1.5])
cbar.ax.set_ylabel(
r"$\log_{10}\left(\delta \omega_{\bm{k}}/ \omega_{\bm{k}}\right)$",
fontsize=16,
)
save_fig(fig, f"fig_nonlinear_ratio_one.png")
if __name__ == "__main__":
plt.show()
import glob
import os
import re
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 fluiddyn.util import modification_date
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,
get_path_finer_resol,
get_paths,
paths_simuls_regimes,
paths_simuls_regimes_proj,
pos_closest_value,
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"
sim = get_sim(letter)
path = paths_simuls_regimes[letter]
sim_proj = get_sim(letter, proj=True)
print(sim_proj.params.projection)
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=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]
# Standard Navier-Stokes
t_start, t_last = times_start_last_from_path(path)
tmin = t_last - 2
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"]
params = load_params_simul(path)
nh = nx = params.oper.nx
proj = params.projection
N = sim.params.N
path_spec = sorted(path.glob(f"spatiotemporal/periodogram_[0-9]*.h5"))
assert len(path_spec) == 1, f"Not only 1 periodogram in {path} \n"
path_spec = path_spec[0]
with h5py.File(path_spec, "r") as f:
# List all groups
# print("Keys: %s" % f.keys())
# a_group_key = list(f.keys())[0]
# 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
xb = np.linspace(kh[1], kb, 50, endpoint=True)
for ax in [ax2, ax4]:
ax.plot(xb, np.sqrt(kb**2 - xb**2), color="k", linestyle="dotted")
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",
)
Ee = np.sum(Ee, axis=2)
cs0 = ax0.scatter(
omega_disp.flatten() / N,
omega_emp.flatten() / N,
c=np.log10(Ee.flatten()),
cmap=cm,
vmin=-7.0,
vmax=-3.0,
# shading="nearest",
)
cs2 = ax2.pcolormesh(
kh,
kz,
np.abs((omega_emp - omega_disp)) / N,
cmap=cm,
vmin=0.0,
vmax=1,
shading="nearest",
)
cs4 = ax4.pcolormesh(
kh,
kz,
np.log10(delta_omega_emp / omega_disp),
cmap=cm,
vmin=-0.5,
vmax=1.5,
shading="nearest",
)
# Without vortical modes
t_start, t_last = times_start_last_from_path(path_proj)
tmin = t_last - 2
k_old = kb
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)
print(k_old, " ", kb)
R2 = mean_values["R2"]
Fh = mean_values["Fh"]
params = load_params_simul(path)
nh = nx = params.oper.nx
proj = params.projection
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:
# List all groups
# print("Keys: %s" % f.keys())
# a_group_key = list(f.keys())[0]
# 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
xb = np.linspace(kh[1], kb, 50, endpoint=True)
for ax in [ax1, ax3]:
ax.plot(xb, np.sqrt(kb**2 - xb**2), color="k", linestyle="dotted")
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",
)
Ee = np.sum(Ee, axis=2)
cs1 = ax1.scatter(
omega_disp.flatten() / N,
omega_emp.flatten() / N,
c=np.log10(Ee.flatten()),
cmap=cm,
vmin=-7,
vmax=-3,
# shading="nearest",
)
cs3 = ax3.pcolormesh(
kh,
kz,
np.abs((omega_emp - omega_disp)) / N,
cmap=cm,
vmin=0.0,
vmax=1,
shading="nearest",
)
cs5 = ax5.pcolormesh(
kh,
kz,
np.log10(delta_omega_emp / omega_disp),
cmap=cm,
vmin=-0.5,
vmax=1.5,
shading="nearest",
)
for ax in [ax0, ax1]:
ax.plot([0, 1], [0, 1], "k-")
ax.set_xlim([0, 1])
ax.set_ylim([0, 2])
for ax in [ax2, ax3, ax4, ax5]:
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.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-")
ax.plot(
[kh[1], kf * 0.35],
[kf * (1 - 0.25**2) ** 0.5, kf * (1 - 0.25**2) ** 0.5],
linestyle="--",
color="orange",
)
ax.plot(
[kf * 0.35, kf * 0.35],
[kz[1], kf * (1 - 0.25**2) ** 0.5],
linestyle="--",
color="orange",
)
# ax.plot(xo, kb*np.sqrt(xo/delta_kh), "m--")
for ax in [ax2, ax3]:
ax.set_xticklabels([])
for ax in [ax1, ax3, ax5]:
ax.set_yticklabels([])
for ax in [ax2, ax3, ax4, ax5]:
ax.set_xlabel(r"$k_h$", fontsize=16)
for ax in [ax2, ax4]:
ax.set_ylabel(r"$k_z$", fontsize=16)
for ax in [ax0, ax1]:
ax.set_xlabel(r"$\omega_{\bm{k}}/N$", fontsize=16)
ax0.set_ylabel(r"$\omega_{emp, \bm{k}}/N$", 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)
ax4.set_title(r"$(e)$", fontsize=16)
ax5.set_title(r"$(f)$", fontsize=16)
fig.tight_layout()
# fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.71, 0.02, 0.23])
cbar = fig.colorbar(cs0, cax=cbar_ax)
cbar.set_ticks([-7, -6, -5, -4, -3])
cbar.ax.set_ylabel(r"$\log_{10} E_{equi}$", fontsize=16)
# fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.385, 0.02, 0.23])
cbar = fig.colorbar(cs2, cax=cbar_ax)
cbar.set_ticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
cbar.ax.set_ylabel(r"$|\omega_{emp, \bm{k}}- \omega_{\bm{k}}|/N$", fontsize=16)
# fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.06, 0.02, 0.23])
cbar = fig.colorbar(cs4, cax=cbar_ax)
cbar.set_ticks([-0.5, 0.0, 0.5, 1.0, 1.5])
cbar.ax.set_ylabel(
r"$\log_{10}\left(\delta \omega_{\bm{k}}/ \omega_{\bm{k}}\right)$",
fontsize=16,
)
fig.subplots_adjust(right=0.85, wspace=0.1, hspace=0.4)
save_fig(fig, f"fig_nonlinear_regimes_{letter}.png")
if __name__ == "__main__":
plt.show()
import glob
import os
import re
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 fluiddyn.util import modification_date
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,
get_path_finer_resol,
get_paths,
paths_simuls_regimes,
paths_simuls_regimes_proj,
pos_closest_value,
save_fig,
)
# Latex
plt.rcParams["text.usetex"] = True
plt.rcParams["text.latex.preamble"] = r"\usepackage{bm}"
cm = matplotlib.cm.get_cmap("binary", 100)
###
N = 80
Rb = 20
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)
fig, axes = plt.subplots(
ncols=2, nrows=2, figsize=(10, 2 * 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0, 0]
ax1 = axes[0, 1]
ax2 = axes[1, 0]
ax3 = axes[1, 1]
# Standard Navier-Stokes
t_start, t_last = times_start_last_from_path(path)
tmin = t_last - 2
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"]
params = load_params_simul(path)
nh = nx = params.oper.nx
proj = params.projection
N = sim.params.N
path_spec = sorted(path.glob(f"spatiotemporal/periodogram_[0-9]*.h5"))
assert len(path_spec) == 1, f"Not only 1 periodogram in {path} \n"
path_spec = path_spec[0]
with h5py.File(path_spec, "r") as f:
# List all groups
# print("Keys: %s" % f.keys())
# a_group_key = list(f.keys())[0]
# 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 = Epolo + EA
omega_emp, delta_omega_emp = compute_omega_emp_vs_kzkh(
N, spectrum, kh, kz, omegas
)
KH, KZ = np.meshgrid(kh, kz)
K = (KH**2 + KZ**2) ** 0.5
K_NOZERO = K.copy()
K_NOZERO[K_NOZERO == 0] = 1e-16
omega_disp = N * KH / K_NOZERO
chi = (K**2 * epsK) ** (1 / 3) / omega_disp
xb = np.linspace(kh[1], kb, 50, endpoint=True)
for ax in [ax0, ax2]:
ax.plot(xb, np.sqrt(kb**2 - xb**2), color="k", linestyle="dotted")
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",
)
cs0 = ax0.pcolormesh(
kh,
kz,
np.abs((omega_emp - omega_disp)) / N,
cmap=cm,
vmin=0.0,
vmax=1,
shading="nearest",
)
cs2 = ax2.pcolormesh(
kh,
kz,
np.log10(delta_omega_emp / omega_disp),
cmap=cm,
vmin=-0.5,
vmax=1.5,
shading="nearest",
)
"""
cs2 = ax2.scatter(
chi*(kb/K)**(2/3),
K/kb,
c = np.log10(delta_omega_emp/omega_disp),
cmap=cm,
vmin=-0.5,
vmax=1.5,
)
"""
# Without vortical modes
t_start, t_last = times_start_last_from_path(path_proj)
tmin = t_last - 2
k_old = kb
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)
print(k_old, " ", kb)
# R2 = mean_values["R2"]
Fh = mean_values["Fh"]
epsK = mean_values["epsK"]
params = load_params_simul(path)
nh = nx = params.oper.nx
proj = params.projection
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:
# List all groups
# print("Keys: %s" % f.keys())
# a_group_key = list(f.keys())[0]
# 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
chi = (K**2 * epsK) ** (1 / 3) / omega_disp
xb = np.linspace(kh[1], kb, 50, endpoint=True)
xa = np.linspace(kh[1], a**1.5 * ko, 50, endpoint=True)
for ax in [ax1, ax3]:
ax.plot(xb, np.sqrt(kb**2 - xb**2), color="k", linestyle="dotted")
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",
)
cs1 = ax1.pcolormesh(
kh,
kz,
np.abs((omega_emp - omega_disp)) / N,
cmap=cm,
vmin=0.0,
vmax=1,
shading="nearest",
)
cs3 = ax3.pcolormesh(
kh,
kz,
np.log10(delta_omega_emp / omega_disp),
cmap=cm,
vmin=-0.5,
vmax=1.5,
shading="nearest",
)
"""
cs3 = ax3.scatter(
chi*(kb/K)**(2/3),
K/kb,
c = np.log10(delta_omega_emp/omega_disp),
cmap=cm,
vmin=-0.5,
vmax=1.5,
)
"""
th = np.linspace(0, np.pi / 2, 50)
for ax in [ax0, ax1, ax2, ax3]:
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.plot([kh[1], max(kh)], [kh[1], max(kh)], "k-")
ax.plot(kf * np.sin(th), kf * np.cos(th), linestyle="--", color="orange")
for ax in [ax0, ax1]:
ax.set_xticklabels([])
for ax in [ax1, ax3]:
ax.set_yticklabels([])
for ax in [ax2, ax3]:
ax.set_xlabel(r"$k_h$", fontsize=12)
for ax in [ax0, ax2]:
ax.set_ylabel(r"$k_z$", fontsize=12)
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)
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.53, 0.02, 0.35])
cbar = fig.colorbar(cs0, cax=cbar_ax)
cbar.set_ticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
cbar.ax.set_ylabel(r"$|\omega_{emp, \bm{k}}- \omega_{\bm{k}}|/N$", fontsize=16)
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([-0.5, 0.0, 0.5, 1.0, 1.5])
cbar.ax.set_ylabel(
r"$\log_{10}\left(\delta \omega_{\bm{k}}/ \omega_{\bm{k}}\right)$",
fontsize=16,
)
# fig.tight_layout()
save_fig(fig, f"fig_nonlinear_regimes_ratio_one_N{N}.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 Fh_limit, plot, 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.set_xlim([1e-3, 20])
ax0.set_xscale("log")
ax0.set_ylim([0, 1])
ax0.set_xlabel(r"$F_h$", fontsize=16)
ax0.set_ylabel(r"$E_{toro}/E$", fontsize=16)
ax0.set_title(r"$(a)$", fontsize=16)
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=16)
ax1.set_xlim([1e-3, 20])
ax1.set_xscale("log")
ax1.set_ylim([-1, 1])
ax1.set_yticks([-1, -0.5, 0.0, 0.5, 1.0])
ax1.set_xlabel(r"$F_h$", fontsize=16)
ax1.set_ylabel(r"$\tilde{\mathcal{D}}$", fontsize=16)
ax1.set_title(r"$(b)$", fontsize=16)
ax1.grid(True)
fig.tight_layout()
fig.subplots_adjust(right=0.85, wspace=0.3)
cbar_ax = fig.add_axes([0.88, 0.15, 0.02, 0.7])
cbar = fig.colorbar(cs, cax=cbar_ax, orientation="vertical")
# cbar.set_ticklabels([-1, 0, 1, 2, 3, 4], position="right")
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.88, 0.07, r"$\log_{10}(\mathcal{R})$", fontsize=16)
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_proj_ratio_one, df_ratio_one
from util import Fh_limit, plot, save_fig
df = df_ratio_one
df_proj = df_proj_ratio_one
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]
print(df)
ax0.scatter(
df["Fh"],
df["Etoro"] / df["E"],
c=np.log10(df["R4"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
vmin=-2,
vmax=3,
)
# ax0.set_xlim([1e-3, 20])
ax0.set_xscale("log")
ax0.set_ylim([0, 1])
ax0.set_xlabel(r"$F_h$", fontsize=16)
ax0.set_ylabel(r"$E_{toro}/E$", fontsize=16)
ax0.set_title(r"$(a)$", fontsize=16)
ax0.grid(True)
cs = ax1.scatter(
df["Fh"],
(df["Epolo"] - df["EA"]) / (df["Epolo"] + df["EA"]),
c=np.log10(df["R4"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="o",
vmin=-2,
vmax=3,
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["R4"]),
cmap="inferno",
edgecolors="k",
s=40,
marker="^",
label=r"Without vortical modes",
)
ax1.legend(loc="lower center", fontsize=16)
# ax1.set_xlim([1e-3, 20])
ax1.set_xscale("log")
ax1.set_ylim([-1, 1])
ax1.set_yticks([-1, -0.5, 0.0, 0.5, 1.0])
ax1.set_xlabel(r"$F_h$", fontsize=16)
ax1.set_ylabel(r"$\tilde{\mathcal{D}}$", fontsize=16)
ax1.set_title(r"$(b)$", fontsize=16)
ax1.grid(True)
fig.tight_layout()
fig.subplots_adjust(right=0.85, wspace=0.3)
cbar_ax = fig.add_axes([0.88, 0.15, 0.02, 0.7])
cbar = fig.colorbar(cs, cax=cbar_ax, orientation="vertical")
cbar.set_ticks([-2, -1, 0, 1, 2, 3])
cbar.set_ticklabels(
[r"$-2$", r"$-1$", r"$0$", r"$1$", r"$2$", r"$3$"], fontsize=14
)
# cbar.set_label(r"$\log_{10}(\mathcal{R})$", fontsize = 12, rotation=0)
fig.text(0.88, 0.07, r"$\log_{10}(\mathcal{R}_4)$", fontsize=16)
save_fig(fig, f"fig_ratio_E_vs_Fh_ratio_one.png")
if __name__ == "__main__":
plt.show()
import matplotlib.pyplot as plt
import numpy as np
from util_dataframe import df
from util import (
Fh_limit,
R2_limit,
formatters,
params_simuls_regimes,
plot,
save_fig,
)
key = ["Epolo", "Etoro", "EA"]
EA = df["EA"]
Epolo = df["Epolo"]
Etoro = df["Etoro"]
E = df["E"]
# Plot 1: (EA+Epolo)/E
ax = plot(
df,
"Fh",
"R2",
c=(EA + Epolo) / E,
cmap="binary", # "binary" "seismic"
logy=True,
vmin=0,
vmax=1,
s=35,
)
ax.set_xlim([1e-3, 20])
ax.set_ylim([1e-1, 1e5])
# ax.axvline(Fh_limit, linestyle=":")
# ax.axvline(1.0, linestyle=":")
# Fh_min, Fh_max = ax.get_xlim()
# ax.plot([Fh_min, Fh_limit], [R2_limit, R2_limit], linestyle=":")
ax.set_xlabel(r"$F_h$", fontsize=16)
ax.set_ylabel(r"$\mathcal{R}$", fontsize=16)
ax.set_title(r"$(E_{polo}+E_{A})/E$", fontsize=16)
fig = ax.figure
fig.tight_layout()
save_fig(fig, f"fig_ratio_EAEpolo_E_vs_FhR.png")
fig, axes = plt.subplots(
ncols=2, nrows=1, figsize=(10, 3 * 4.5 / 4), constrained_layout=True
)
ax0 = axes[0]
ax1 = axes[1]
# Plot 2: Etoro/E
ax0 = plot(
df,
"Fh",
"R2",
c=Etoro / E,
cmap="binary", # "binary" "seismic"
logy=True,
vmin=0,
vmax=1,
s=35,
)
ax0.set_xlim([1e-3, 20])
ax0.set_ylim([1e-1, 1e5])
ax0.set_xlabel(r"$F_h$", fontsize=16)
ax0.set_ylabel(r"$\mathcal{R}$", fontsize=16)
ax0.set_title(r"$E_{toro}/E$", fontsize=16)
save_fig(fig, f"fig_ratio_Etoro_E_vs_FhR.png")
# Plot 3: (Epolo-EA)/(Epolo+EA)
ax1 = plot(
df,
"Fh",
"R2",
c=(Epolo - EA) / (Epolo + EA),
cmap="seismic", # "binary" "seismic"
logy=True,
vmin=-1,
vmax=1,
s=35,
)
ax1.set_xlim([1e-3, 20])
ax1.set_ylim([1e-1, 1e5])
ax1.set_xlabel(r"$F_h$", fontsize=16)
ax1.set_ylabel(r"$\mathcal{R}$", fontsize=16)
ax1.set_title(r"$(E_{polo}-E_{A})/(E_{polo}+E_{A})$", fontsize=16)
fig = axes.figure
fig.tight_layout()
save_fig(fig, f"fig_ratio_Epolo-EA_Epolo+EA_vs_FhR.png")
if __name__ == "__main__":
plt.show()
import sys
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm, ticker
from util_simuls_regimes import get_sim
from util import compute_kf_kb_ko_keta_kd, customize, save_fig
plt.rcParams["text.usetex"] = True
print(sys.argv)
proj = bool(sys.argv[-1])
if proj not in [True, False]:
letter = False
sim = get_sim("L", proj)
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)
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]
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(
[kh[1], kf * 0.35],
[kf * (1 - 0.25**2) ** 0.5, kf * (1 - 0.25**2) ** 0.5],
linestyle="--",
color="orange",
)
ax.plot(
[kf * 0.35, kf * 0.35],
[kz[1], kf * (1 - 0.25**2) ** 0.5],
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.tight_layout()
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])
save_fig(fig, f"fig_seb_regimes_L_proj{proj}.png")
if __name__ == "__main__":
plt.show()