3-D thermal convection in a box

Bottom-heated, isoviscous convection in a unit cube at Ra = 10⁵ with free-slip walls: cold sheets and hot plumes grow from a perturbed conductive state into time-dependent flow.

Run it with:

julia -t auto --project=. main.jl cases/convection3d.toml

Full case file — cases/convection3d.toml:

# ══════════════════════════════════════════════════════════════════════════════
# 3D thermal convection in a unit cube, Ra = 10⁵, isoviscous, bottom heated,
# free-slip walls.
# Starts from a linear conductive profile with a multi-mode perturbation; cold
# sheets and hot plumes grow and develop into time-dependent flow.
# Box 1×1×1, dimensionless (η = ρ0 = g = κ = ΔT = 1, α = Ra).
# ══════════════════════════════════════════════════════════════════════════════

[model]
case = "convection3d"

dim   = 3
xsize = 1.0
ysize = 1.0
zsize = 1.0
nx    = 20
ny    = 20
nz    = 20
mnx   = 2
mny   = 2
mnz   = 2

gx = 0.0
gy = 1.0          # +y is down; Ra enters through α (rock alpha = Ra)
gz = 0.0
pinit = 0.0

eta_min = 1.0
eta_max = 1.0

maxxystep = 0.4
maxtkstep = 0.1
maxtmstep = 5.0e-4
ntot      = 500

do_stokes        = true
do_heat          = true
do_advect        = true
do_adiabatic     = false
do_shear_heating = false
do_friction      = false
do_phase_change  = false
do_melting       = false
do_free_surface  = false
do_radiogenic    = false

fe_order_v = 2
fe_order_p = 1
fe_order_T = 1
n_picard   = 1
picard_tol = 1e-3

solver_type  = "mumps"
krylov_tol   = 1e-6
krylov_maxit = 200
gmg_levels   = 2

output_dir   = "output/convection3d"
output_every = 10
str_min      = 1e-18

T_top         = 0.0
T_bot         = 1.0
T_ref_density = 0.0
rho_ref_heat  = 1.0
body_force_full = false
marker_T_mode = "pic_increment"
buoyancy_from_T = true

bc_left    = "free_slip"
bc_right   = "free_slip"
bc_surface = "free_slip"
bc_deep    = "free_slip"
bc_front   = "free_slip"
bc_back    = "free_slip"

trench_x = 0.0
plate_thickness = 0.0
plate_age = 1.0
continent_thickness = 0.0
air_thick = 0.0
x_plate_left = 0.0
x_plate_right = 0.0
vx_plate = 0.0
vx_plate_right = 0.0
air_rock_id = 1
surface_smooth = 0.0

[[rock]]
nu0      = 1.0
dh       = 0.0
dv       = 0.0
ss       = 0.0
mm       = 1.0
ro0      = 1.0
alpha    = 1.0e5     # Ra (ρ0 = g = κ = η = ΔT = H = 1)
beta     = 0.0
cp       = 1.0
kt       = 1.0
ht       = 0.0
a0       = 0.0
b0       = 0.0
a1       = 0.0
b1       = 0.0
e0       = 0.0
e1       = 1.0
nu_min   = 1.0
nu_max   = 1.0
immobile = false
G        = 0.0