Rayleigh-Bénard convection of a Darcy fluid in a porous rectangle

\[\begin{split} \mathbb{S}_{\psi,c} \begin{cases} \Omega = [0, \mathcal{A}X] \times [0, X] & \text{aspect ratio } \mathcal{A}=\mathcal{O}(1)\\ c_0(x,y)=1-y+\mathcal{N}(x,y) & \text{perturbed initial linear concencentration} \\ c_{\text{D}}(x,y=0)=1 & \text{light lower boundary} \\ c_{\text{D}}(x,y=X)=0 & \text{heavy upper boundary} \\ c_{\text{N}}(x=0,y)=0 & \text{no-flux on left boundary}\\ c_{\text{N}}(x=\mathcal{A}X,y)=0 & \text{no-flux on right boundary}\\ \psi_{\text{D}}\vert_{\partial\Omega}=0 & \text{no-penetration on entire boundary} \\ \phi = 1 & \text{constant porosity} \\ \mathsf{D} = \mathsf{I} & \text{constant isotropic dispersion}\\ \mathsf{K} = \mathsf{I} & \text{constant isotropic permeability}\\ \mu = 1 & \text{constant viscosity} \\ \rho(c) = -c & \text{linear density} \\ \textbf{e}_g=-\textbf{e}_y & \text{vertically downward gravity} \\ \end{cases} \end{split}\]

import numpy as np
from lucifex.fdm import AB2, CN
from lucifex.sim import run
from lucifex.mesh import as_grid_mesh
from lucifex.fem import cross_section_grid
from lucifex.utils.npy_utils import as_index
from lucifex.plt import plot_colormap, plot_line, create_animation, save_figure, display_animation
from py.C11_darcy_rayleigh_benard import darcy_rayleigh_benard_rectangle

simulation = darcy_rayleigh_benard_rectangle(
    aspect=2.0,  
    Nx=64,
    Ny=64,
    cell='quadrilateral', 
    scaling='advective',
    Ra=500.0, 
    c_ampl=1e-3, 
    c_freq=(12, 8), 
    c_seed=(123, 456), 
    D_adv=AB2,
    D_diff=CN,
    diagnostic=True,
)

n_stop = 200
dt_init = 1e-6
n_init = 5
run(simulation, n_stop=n_stop, dt_init=dt_init, n_init=n_init)

c = simulation['c']

INFO:root:running build_ext

INFO:root:building 'libffcx_forms_7d56d224dfc767a25f4745d7b4a5a89b2c2c261d' extension

INFO:root:clang -Wno-unused-result -Wsign-compare -Wunreachable-code -DNDEBUG -fwrapv -O2 -Wall -fPIC -O2 -isystem /Users/George/miniconda3/envs/lucifex/include -fPIC -O2 -isystem /Users/George/miniconda3/envs/lucifex/include -I/Users/George/miniconda3/envs/lucifex/lib/python3.10/site-packages/ffcx/codegeneration -I/Users/George/miniconda3/envs/lucifex/include/python3.10 -c libffcx_forms_7d56d224dfc767a25f4745d7b4a5a89b2c2c261d.c -o ./libffcx_forms_7d56d224dfc767a25f4745d7b4a5a89b2c2c261d.o -O2 -g0

INFO:root:clang -bundle -undefined dynamic_lookup -Wl,-rpath,/Users/George/miniconda3/envs/lucifex/lib -L/Users/George/miniconda3/envs/lucifex/lib -Wl,-rpath,/Users/George/miniconda3/envs/lucifex/lib -L/Users/George/miniconda3/envs/lucifex/lib ./libffcx_forms_7d56d224dfc767a25f4745d7b4a5a89b2c2c261d.o -o ./libffcx_forms_7d56d224dfc767a25f4745d7b4a5a89b2c2c261d.cpython-310-darwin.so

INFO:root:running build_ext

INFO:root:building 'libffcx_forms_d81b7186c0ce67076dc0c39477a813ca213f4d1f' extension

INFO:root:clang -Wno-unused-result -Wsign-compare -Wunreachable-code -DNDEBUG -fwrapv -O2 -Wall -fPIC -O2 -isystem /Users/George/miniconda3/envs/lucifex/include -fPIC -O2 -isystem /Users/George/miniconda3/envs/lucifex/include -I/Users/George/miniconda3/envs/lucifex/lib/python3.10/site-packages/ffcx/codegeneration -I/Users/George/miniconda3/envs/lucifex/include/python3.10 -c libffcx_forms_d81b7186c0ce67076dc0c39477a813ca213f4d1f.c -o ./libffcx_forms_d81b7186c0ce67076dc0c39477a813ca213f4d1f.o -O2 -g0

INFO:root:clang -bundle -undefined dynamic_lookup -Wl,-rpath,/Users/George/miniconda3/envs/lucifex/lib -L/Users/George/miniconda3/envs/lucifex/lib -Wl,-rpath,/Users/George/miniconda3/envs/lucifex/lib -L/Users/George/miniconda3/envs/lucifex/lib ./libffcx_forms_d81b7186c0ce67076dc0c39477a813ca213f4d1f.o -o ./libffcx_forms_d81b7186c0ce67076dc0c39477a813ca213f4d1f.cpython-310-darwin.so

time_slice = slice(0, None, 2)
titles = [f'$c(t={t:.3f})$' for t in c.time_series[time_slice]]

anim = create_animation(
    plot_colormap,
    colorbar=False,
)(c.series[time_slice], title=titles)
anim_path = save_figure(f'{c.name}(x,y,t)', return_path=True)(anim)

display_animation(anim_path)

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

grid_mesh = as_grid_mesh(c.mesh)
x, y = grid_mesh.axes

y_index = 8
y_value = y[y_index]
cx_series, _ = cross_section_grid(c.series, 'y', y_value)
cxt = np.array([cx.value for cx in cx_series]).T

title = f'$c(y={y_value:.3f})$'
fig, ax = plot_colormap(
    (x, c.time_series, cxt), 
    aspect='auto',
    x_label='$x$',
    y_label='$t$',
    title=title,
    colorbar=(0, 1),
)
save_figure(title.strip('$'))(fig)

uRMS, uMinMax = simulation['uRMS', 'uMinMax']
uMax = uMinMax.sub(1)

fig, ax = plot_line(
    (uRMS.time_series, uRMS.value_series),
    x_label='$t$',
    y_label='$\mathrm{rms}(\\textbf{u})$',
)
save_figure('uRMS(t)')(fig)

fig, ax = plot_line(
    (uMax.time_series, uMax.value_series),
    x_label='$t$',
    y_label='$\max_{\\textbf{x}}|\\textbf{u}|$',
)
save_figure('uMax(t)')(fig)

dt, dtCFL = simulation['dt', 'dtCFL']

fig, ax = plot_line(
    [(dt.time_series, dt.value_series), (dtCFL.time_series, dtCFL.value_series)],
    x_label='$t$',
    legend_labels=['$\Delta t$', '$\Delta t_{\mathbf{u}}$'],
)
ax.set_yscale('log')
save_figure('dt(t)')(fig)

cMinMax = simulation['cMinMax']
cMin, cMax = cMinMax.split()

fig, ax = plot_line(
    (cMax.time_series, cMax.value_series),
    x_label='$t$',
    y_label='$\max_{\\textbf{x}}(c)$',
)
save_figure('cMax(t)')(fig)

fig, ax = plot_line(
    (cMin.time_series, cMin.value_series),
    x_label='$t$',
    y_label='$\min_{\\textbf{x}}(c)$',
)
save_figure('cMin(t)')(fig)