Isoviscous convection, Ra = 10⁴ to 10⁶ (2-D)

Basally heated, isoviscous, Boussinesq convection in the spherical annulus, the canonical problem of Hernlund & Tackley (2008). One driver covers the three Rayleigh numbers; the Rayleigh number sets the vigour of the flow, the thickness of the thermal boundary layers and how many plumes the annulus holds, and above about 10⁵ the flow is time dependent.

This page is the "look at the flow" companion of the verification driver benchmarks/hernlund_tackley_annulus.jl. The physics and the set-up are the same; the example writes VTK frames and figures and is shorter by default, the benchmark carries the reference tables and a convergence guard. Quote numbers from the benchmark, as described in Benchmarks.

The model

$$\mathrm{\nabla }\cdot \left(2\eta \dot{\mathit{\epsilon }}\right)-\mathrm{\nabla }p={\rho }_0\alpha (T-T_0)g\hat{\mathbf{r}},\mathrm{\nabla }\cdot \mathbf{v}=0,\frac{\partial T}{\partial t}+\mathbf{v}\cdot \mathrm{\nabla }T=\kappa {\mathrm{\nabla }}^{2}T,$$

with constant $\eta$ chosen from Ra (eta_for_Ra in examples/common_2d.jl), free-slip walls, $T=1300$ K at the CMB and $300$ K at the surface, and the net rotation removed after every solve (cfg.removeNetRotation). The equations and their discretisation are described in Governing equations and Discretization.

parameter	value	where
geometry	spherical annulus, $f=r_i/r_o=0.5495$, $D=1000$ km	ht_shell!
grid	Nr × 8 Nr, uniform in radius (default 64 × 512)	ht_shell!
Ra	${10}^4$, ${10}^5$ or ${10}^6$ (first argument)	eta_for_Ra
$\eta$	$9.7\times {10}^{22}$, $9.7\times {10}^{21}$, $9.7\times {10}^{20}$ Pa s
$\mathrm{\Delta }T$	1000 K
temperature transport	markers, with the grid diffusion increment handed back (Crank-Nicolson, $\theta =0.5$)	Markers
Stokes solver	MUMPS direct	cfg.solverType = "mumps"
initial condition	linear ramp, $m$ warm cosine peaks of 20 K at the CMB, 1 K noise	_initial_temperature_2d in markers/init_markers.jl
dtMax	$5\times {10}^{13}$, ${10}^{13}$, $5\times {10}^{12}$ s	CFL controls once convection starts
run length	0.3 thermal diffusion times, about 9500 Myr	cfg.maxTimeMyr
output	VTK and figures every 2000 steps, diagnostics.csv every 500

Planform

At these Rayleigh numbers the annulus supports several stable planforms, and the seeded wavenumber $m$ (third argument, default 6) decides which one the run settles into. At Ra = 10⁵, Nr = 64, seeding $m=4,6,8$ gives $V_{\mathrm{rms}}$ = 176.2, 158.2, 156.4 and Nu = 9.05, 8.88, 8.40. Quote $m$ with every number.

A short preview run of the driver at Ra = 10⁵, seeded m = 6, 64 × 512 (3000 steps, 0 to 910 Myr of model time, about 0.03 diffusion times; a frame every 60 steps): temperature, cold blue to hot red. The seeded perturbation grows into six hot upwellings within the first 100 Myr; the cold boundary layer then drips between them and the flow settles into a time-dependent pattern of narrow plumes and cold sheets. This is far shorter than the 0.3 diffusion times the driver runs by default, so it shows the onset, not a statistically steady state; take numbers from the full run.

The same run at 0, 56, 188 and 891 Myr.

Running it

Locally, one Rayleigh number at a time:

export GKSwstype=nul                      # headless figure rendering
julia --project=. -t 4 examples/annulus_convection_2d.jl 1e5 64 6

Arguments are Ra Nr m. Shorten a run with MILEVA_MAXTIME_MYR or MILEVA_MAXSTEPS, and redirect it with MILEVA_OUTDIR. Output goes to output/examples/annulus_Ra1e+05_N64_m6/.

On a cluster, hpc/annulus_convection_2d.slurm runs all three Rayleigh numbers as a SLURM array (see Running on HPC):

sbatch hpc/annulus_convection_2d.slurm

At the end the driver prints the mean of ${\mathrm{Nu}}_{\mathrm{top}}$, ${\mathrm{Nu}}_{\mathrm{bot}}$ and $V_{\mathrm{rms}}$ over the last third of the recorded steps, next to the Hernlund & Tackley (2008) Table 1 values for orientation.

What to expect

The numbers below were measured with the benchmark driver.

• Ra = 10⁴: steady convection. With the benchmark driver at Nr = 64, $m=6$, MILEVA measured Nu = 4.14 and $V_{\mathrm{rms}}$ = 35.1 against the published 4.18 and 37.7. The isoviscous Ra = 10⁴ case settles at about 0.3 diffusion times, which is this example's cap; for a converged number use the benchmark driver, which runs to 0.6.

• Ra = 10⁵: time-dependent convection with thinner boundary layers, and more than one stable planform. Seeded with $m=3$, the benchmark driver settles into the three plumes of the published run and measures Nu = 7.70 and $V_{\mathrm{rms}}$ = 161.8 against 7.39 and 160. Seeded with $m=6$ it keeps more, narrower cells and Nu = 8.88, about 20 % higher; see Benchmarks.

• Ra = 10⁶: strongly time dependent flow with many thin plumes; the published annulus values are Nu = 14.4 and $V_{\mathrm{rms}}$ = 640. Use Nr = 128 or more to resolve the boundary layers.

A useful internal check that needs no reference: in a statistically steady state the same power crosses every shell, so the time means of ${\mathrm{Nu}}_{\mathrm{top}}$ and ${\mathrm{Nu}}_{\mathrm{bot}}$ must agree. If they do not, the run has not equilibrated.

Smoke test

A quick check that the driver runs:

MILEVA_MAXSTEPS=6 julia --project=. -t 2 examples/annulus_convection_2d.jl 1e5 16 6

6 steps at 16 × 128 in 41 s of wall time, most of it compilation; temperature and velocity finite throughout. At the default 64 × 512 grid, 40 steps took 0.41 s per step and 2.3 GB.