Thermosolutal convection of a Navier-Stokes fluid in a rectangle

\[\begin{split} \mathbb{S}_{\mathbf{u},p,c,\theta} \begin{cases} \Omega = [0, \mathcal{A}X] \times [0, X] & \text{aspect ratio } \mathcal{A}=\mathcal{O}(1)\\ c_0(x,y)=y+\mathcal{N}(x,y) & \text{perturbed initial linear concentration} \\ \theta_0(x,y)=1-y+\mathcal{N}(x,y) & \text{perturbed initial linear temperature} \\ \textbf{u}_0=\textbf{0} & \text{static initial velocity} \\ p_0=0 & \text{static initial pressure} \\ c_{\text{D}}(x,y=0)=0 & \text{light lower boundary} \\ c_{\text{D}}(x,y=X)=1 & \text{heavy upper boundary} \\ \theta_{\text{D}}(x,y=0)=1 & \text{hot lower boundary} \\ \theta_{\text{D}}(x,y=X)=0 & \text{cold upper boundary} \\ c_{\text{N}}(x=0,y)=0 & \text{no-solutal-flux on left boundary}\\ c_{\text{N}}(x=\mathcal{A}X,y)=0 & \text{no-solutal-flux on right boundary}\\ \theta_{\text{N}}(x=0,y)=0 & \text{no-thermal-flux on left boundary}\\ \theta_{\text{N}}(x=\mathcal{A}X,y)=0 & \text{no-thermal-flux on right boundary}\\ \textbf{u}_{\text{E}}\vert_{\partial\Omega}=\textbf{0} & \text{no-flow on entire boundary} \\ \rho(c, \theta) = c - \gamma\theta & \text{linear density} \\ \textbf{e}_g=-\textbf{e}_y & \text{vertically downward gravity} \\ \end{cases} \end{split}\]

from lucifex.sim import run
from lucifex.utils.npy_utils import as_index
from lucifex.plt import (
    plot_line, create_animation, plot_colormap,
    save_figure, display_animation, plot_colormap_multifigure,
)
from lucifex.solver import maximum
from py.C31_navier_stokes_thermosolutal import navier_stokes_thermosolutal_rectangle

simulation = navier_stokes_thermosolutal_rectangle(
    aspect=2.0,
    Nx=64,
    Ny=64,
    Ra=5e4,
    gamma=2.0,
    Le=10.0,
    dt_max=0.01, 
)

n_stop = 400
dt_init = 1e-6 
n_init = 10
run(simulation, n_stop=n_stop, dt_init=dt_init, n_init=n_init)

c, theta, rho, u = simulation['c', 'theta', 'rho', 'u']
mesh = c.mesh

title_triplet = lambda time_triplet, tex_triplet: (
    tuple(f'${n}(t={t:.6f})$' for t, n in zip(time_triplet, tex_triplet))
)
series_triplet = (c, theta, rho)
tex_triplet = ('c', '\\theta', '\\rho')

time_slice = slice(0, None, 2)
time_triplet_series = [ti for ti in zip(*[s.time_series[time_slice] for s in series_triplet])]

anim = create_animation(
    plot_colormap_multifigure(n_cols=3, cbars=False),
    cmap=('magma', 'plasma', 'hot'),
)(
    [ui for ui in zip(*[s.series[time_slice] for s in series_triplet])],
    title=[title_triplet(i, tex_triplet) for i in time_triplet_series],
)
anim_path = save_figure('c(x,y,t)_theta(x,y,t)_rho(x,y,t)', return_path=True)(anim)

display_animation(anim_path)

time_index = -1
time_triplet = tuple(u.time_series[time_index] for u in series_triplet)

mfig, *_ = plot_colormap_multifigure(n_cols=3, cbars=True)(
    [u.series[time_index] for u in series_triplet],
    ('hot', 'magma', 'plasma'),
    title=title_triplet(time_triplet, tex_triplet),
)
save_figure(f'c(x,y,t={time_triplet[0]})_theta(x,y,t={time_triplet[1]})_rho(x,y,t={time_triplet[2]})')(fig)

time_indices = as_index(rho.time_series, (0, 0.5, -1), fraction=True)
for i in time_indices:
    fig, ax = plot_colormap(rho.series[i], title=f'$\\rho(t={rho.time_series[i]:.3f})$')
    save_figure(f'rho(x,y,t={rho.time_series[i]:.3f})', thumbnail=(i == -1))(fig)

uMax = [maximum(i) for i in u.series]
fig, ax = plot_line(
    (u.time_series, uMax), 
    x_label='$t$', 
    y_label='$\max_{\\textbf{x}}|\\textbf{u}|$',
)
save_figure('uMax(t)')(fig)