A 2D pot heated from below — Rayleigh–Bénard convection (heat + flow)

Download .py All tutorials

Heat a layer of fluid from below and, above a critical temperature difference, it stops merely conducting: hot fluid becomes buoyant and rises, cold fluid sinks, and the layer organises into rolling convection cells. This is the canonical multiphysics problem — incompressible flow coupled to heat transport — and it falls out of one jno.fem([...]).

The Boussinesq model — three fields, two-way coupling

\[ \begin{aligned} \partial_t u + (u\cdot\nabla)u &= -\nabla p + \mathrm{Pr}\,\nabla^2 u + \mathrm{Pr}\,\mathrm{Ra}\,T\,\hat e_y &&\text{(momentum + buoyancy)}\\ \nabla\cdot u &= 0 &&\text{(incompressible)}\\ \partial_t T + u\cdot\nabla T &= \nabla^2 T &&\text{(heat: advected + diffused)} \end{aligned} \]

Velocity \(u\) (P2), pressure \(p\) (P1) and temperature \(T\) (P1) are three coupled fields. The coupling is genuinely two-way, and that is what makes it multiphysics rather than two solves glued together:

buoyancy \(\mathrm{Pr}\,\mathrm{Ra}\,T\) feeds temperature into the momentum balance — a linear cross term, like the exchange term in the two-temperature example;
advection \(u\cdot\nabla T\) feeds the flow into the heat balance — a product of two different unknowns, so the whole system is nonlinear and routes through the coupled Newton path (the same machinery as the Navier–Stokes convective term \((u\cdot\nabla)u\)).

u, v = d.fem_symbols(value_shape=(2,), order=2)   # velocity (P2)
p, q = d.fem_symbols(order=1)                       # pressure (P1)
T, s = d.fem_symbols(order=1)                        # temperature (P1)

momentum = (u.t·v) + (u·∇)u·v + Pr·∇u:∇v - p·div v - Pr*Ra*T·v_y   # buoyancy: T -> momentum
continuity = q · div u
energy   = T.t·s + (u·∇T)·s + ∇T·∇s                                # advection: u -> heat
fem = jno.fem([momentum, continuity, energy, *walls, *initial_condition])

fem.solve() marches it implicitly (backward-Euler + Newton per step, internally) and returns the forward trajectory; we watch the rolls grow from rest.

A pot: no-slip walls, hot floor, cold lid

no-slip on every wall — u(walls) - 0 (one all-component vector Dirichlet);
a hot floor / cold lid, with the conductive profile held on the walls — T(walls) - (1 - y);
the fluid starts at rest from the conductive state plus a tiny perturbation that seeds the instability — u(initial) - 0 (an all-component vector initial condition) and T(initial) - T0.

The all-component vector u=0 Dirichlet and initial condition both exercise multifield BC/IC paths that this work completed (a vector field's BC/IC no longer has to be spelled out per component).

The result

A rectangular pot, temperature shown hot (red) at the floor to cold (blue) at the lid, with the
computed velocity arrows forming counter-rotating convection rolls — hot plumes rising, cold plumes
sinking.

Temperature is the colour (hot floor → cold lid); the arrows are the computed velocity. From rest, the conductive layer is unstable: small perturbations grow into rising hot plumes and sinking cold ones, settling into steady counter-rotating rolls. Both fields are the actual finite-element solution — nothing painted in.

What to notice

One coupled solve. Velocity, pressure and temperature are assembled and Newton-solved together; the buoyancy and advection cross-terms are just ordinary terms in the weak form.
Convection genuinely onsets. The script starts from rest (max|u| = 0) and asserts the flow grows and carries hot fluid upward — the instability is computed, not imposed.
Generality by reuse. This is the nonlinear coupled-multifield path (Navier–Stokes + two temperatures, combined): vector + scalar fields, mixed P2/P1 orders, transient, on a real domain.

Full script

"""12 - A 2D "pot heated from below": Rayleigh-Benard convection (heat + fluid flow, fully coupled).

Heat a layer of fluid from below and it does not just conduct -- above a critical temperature
difference the hot fluid becomes buoyant, rises, cold fluid sinks, and the layer breaks into rolling
convection cells. This is the **Boussinesq** model: incompressible Navier-Stokes with a buoyancy body
force proportional to temperature, two-way coupled to an advection-diffusion equation for the heat.

    du/dt + (u.grad)u = -grad p + Pr lap u + Pr*Ra * T e_y      (momentum + buoyancy)
    div u = 0                                                    (incompressible)
    dT/dt + u.grad T  = lap T                                    (heat: advected + diffused)

Three coupled fields -- velocity ``u`` (P2), pressure ``p`` (P1), temperature ``T`` (P1) -- written as
a single ``jno.fem([...])``. The coupling is genuinely two-way: **buoyancy** ``Pr*Ra*T`` feeds heat
into the momentum balance (a linear cross term), and **advection** ``u.grad T`` feeds the flow into the
heat balance (a product of two *different* unknowns -> nonlinear). The whole system routes through the
coupled nonlinear Newton path; ``fem.solve()`` marches it in time internally -- backward-Euler + Newton
at the domain's time step -- and we watch the rolls grow from rest.

Boundary/initial conditions for the pot: **no-slip** walls (``u=0``), a **hot floor / cold lid** with
the conductive profile held on the walls (``T = 1 - y``), and the fluid starting **at rest** from the
conductive state plus a tiny perturbation that seeds the instability. The animation is the *computed*
temperature with the *computed* velocity arrows -- nothing painted in.
"""

import os

os.environ["MPLBACKEND"] = "Agg"
os.environ.setdefault("XLA_PYTHON_CLIENT_MEM_FRACTION", "0.5")

import jax

jax.config.update("jax_enable_x64", True)  # the assembler builds in float64

import matplotlib.tri as mtri  # noqa: E402  (used by the convective-flux verification below)
import numpy as np  # noqa: E402

import jno  # noqa: E402

Pr, Ra = 1.0, 1.0e4  # Prandtl, Rayleigh (Ra >> Ra_c ~ 1708 -> vigorous convection)
Lx, Ly = 2.0, 1.0  # a wide-ish pot -> a pair of counter-rotating rolls
dt, nsteps, nframes = 0.009, 26, 13  # integration step / count -> fem.solve() reads dt from the domain
#                                      time grid below (stop ~when the rolls establish, no static tail)

d = jno.Shape.rect(0, 0, Lx, Ly, size=0.11).domain(time=(0.0, nsteps * dt, nsteps + 1))
u, v = d.fem_symbols(value_shape=(2,), names=("u", "v"), order=2)  # P2 velocity
p, q = d.fem_symbols(names=("p", "q"), order=1)  # P1 pressure
T, sT = d.fem_symbols(names=("T", "sT"), order=1)  # P1 temperature
xi, yi, ti = d.variable("interior", split=True)
xb, yb, _ = d.variable("boundary", split=True)
ci = d.variable("initial", split=True)
ub, vb = u.bind(x=xi, y=yi, t=ti), v.bind(x=xi, y=yi, t=ti)
pb, qb = p.bind(x=xi, y=yi, t=ti), q.bind(x=xi, y=yi, t=ti)
Tb, sb = T.bind(x=xi, y=yi, t=ti), sT.bind(x=xi, y=yi, t=ti)

ux, uy, vx, vy = ub[0], ub[1], vb[0], vb[1]
uxx, uxy, uyx, uyy = ub.x[0], ub.y[0], ub.x[1], ub.y[1]  # grad-then-index: d u_i / d x_j
vxx, vxy, vyx, vyy = vb.x[0], vb.y[0], vb.x[1], vb.y[1]
momentum = (
    (ub.t[0] * vx + ub.t[1] * vy)  # du/dt
    + ((ux * uxx + uy * uxy) * vx + (ux * uyx + uy * uyy) * vy)  # (u.grad)u  -- nonlinear
    + Pr * (uxx * vxx + uxy * vxy + uyx * vyx + uyy * vyy)  # Pr grad u : grad v
    - pb * (vxx + vyy)  # -p div v
    - Pr * Ra * Tb * vy  # buoyancy: temperature -> momentum
)
continuity = qb * (uxx + uyy)  # div u = 0
energy = Tb.t * sb + (ux * Tb.x + uy * Tb.y) * sb + (Tb.x * sb.x + Tb.y * sb.y)  # dT/dt + u.grad T - lap T

Tcond = 1.0 - ci[1] / Ly  # conductive profile: 1 (hot) at the floor, 0 (cold) at the lid
T0 = Tcond + 0.05 * jno.np.sin(2 * np.pi * ci[0] / Lx) * jno.np.sin(np.pi * ci[1] / Ly)  # seed the rolls
fem = jno.fem(
    [
        momentum,
        continuity,
        energy,
        u(xb, yb) - 0.0,  # no-slip walls (all-component vector Dirichlet)
        T(xb, yb) - (1.0 - yb / Ly),  # hot floor / cold lid, conductive profile on the walls
        p.pin(),  # gauge-fix: remove the pressure null space
        u(*ci) - 0.0,  # start at rest (all-component vector initial condition)
        T(*ci) - T0,
    ]
)
assert fem.is_transient and not fem.is_linear, "Boussinesq convection is transient + nonlinear"
off = fem.offsets
print(f"\n2D Rayleigh-Benard pot (Ra={Ra:g}, Pr={Pr:g}): dofs={fem.dofs}, steps={nsteps}")

# fem.solve() marches the coupled nonlinear system itself -- backward Euler + Newton per step over the
# domain's time grid -- and returns a differentiable trace node; evaluate the forward trajectory (one
# full DOF vector per step, row 0 = the rest initial state) through a minimal crux.
sol = fem.solve()
traj = np.asarray(jno.core([sol.mse]).eval([sol]))  # (nsteps + 1, dofs)
save_every = max(1, nsteps // nframes)
frames = traj[::save_every]  # subsample rows for the animation (row 0 = rest, then every save_every)

pts_v = np.asarray(fem.field_points[0])  # P2 velocity nodes
pts_T = np.asarray(fem.field_points[2])  # P1 temperature nodes
vel = frames[:, off[0] : off[1]].reshape(frames.shape[0], -1, 2)
Tf = frames[:, off[2] :]  # temperature frames
tris = np.asarray(d.built_mesh.cells_dict["triangle"])
triT = mtri.Triangulation(pts_T[:, 0], pts_T[:, 1], tris)
umax0, umaxF = float(np.abs(vel[0]).max()), float(np.abs(vel[-1]).max())
# convective vertical heat flux <u_y * T> (T interpolated to the velocity nodes): positive for
# Rayleigh-Benard convection -- hot fluid rises (u_y>0, T high), cold sinks (u_y<0, T low).
Tvel = np.asarray(mtri.LinearTriInterpolator(triT, Tf[-1])(pts_v[:, 0], pts_v[:, 1]))
conv_flux = float(np.nanmean(vel[-1, :, 1] * Tvel))
print(f"  convection onset: max|u| {umax0:.3f} (rest) -> {umaxF:.2f}  |  convective heat flux <u_y T> = {conv_flux:+.3f}")
assert umaxF > 1.0 and conv_flux > 0.0, "expected convection to develop with upward heat transport"