Governing equations

MILET solves creeping (zero-Reynolds-number) thermo-mechanical flow: a mixed Stokes problem for velocity and pressure with visco-elasto-plastic (VEP) rheology, coupled to an energy equation (Heat transport) and a marker-in-cell material description (Markers). This page covers the Stokes system as implemented in src/physics/Stokes.jl, the buoyancy term from src/physics/Gravity.jl / src/time_stepping/PicardLoop.jl, and the compressible-form variants in src/physics/Compressibility.jl.

Strong form

In the default Boussinesq formulation (formulation = "boussinesq"):

$$-\mathrm{\nabla }\cdot (2{\eta }_{\mathrm{VEP}}\epsilon (\mathbf{u}))+\mathrm{\nabla }p=\rho \mathbf{g}\text{(momentum)}$$$$\mathrm{\nabla }\cdot \mathbf{u}=0\text{(mass)}$$

with $\epsilon (\mathbf{u})=\frac{1}{2}(\mathrm{\nabla }\mathbf{u}+\mathrm{\nabla }{\mathbf{u}}^{\mathsf{T}})$ the symmetric strain-rate tensor and ${\eta }_{\mathrm{VEP}}$ the effective visco-elasto-plastic viscosity (see Rheology). MILET uses a y-positive-down convention, so gravity is gy = +9.81 in the [gravity] TOML section.

Weak form

build_stokes_forms (src/physics/Stokes.jl) assembles the bilinear and linear forms. Multiplying momentum by a test function $\mathbf{v}\in V$, integrating the stress term by parts, and testing continuity with $q\in Q$:

$$a((\mathbf{u},p),(\mathbf{v},q))={\int }_{\mathrm{\Omega }}2{\eta }_{\mathrm{VEP}}\epsilon (\mathbf{u}):\epsilon (\mathbf{v})d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}p(\mathrm{\nabla }\cdot \mathbf{v})d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}q(\mathrm{\nabla }\cdot \mathbf{u})d\mathrm{\Omega }+\sum _{{\mathrm{\Gamma }}_{\mathrm{fs}}}{\int }_{\mathrm{\Gamma }}\gamma (\mathbf{u}\cdot \mathbf{n})(\mathbf{v}\cdot \mathbf{n})d\mathrm{\Gamma }$$$$l((\mathbf{v},q))={\int }_{\mathrm{\Omega }}\rho \mathbf{g}\cdot \mathbf{v}d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}\chi {\mathit{\sigma }}_{\mathrm{old}}:\epsilon (\mathbf{v})d\mathrm{\Omega }+\sum _{{\mathrm{\Gamma }}_t}{\int }_{\mathrm{\Gamma }}\mathbf{t}\cdot \mathbf{v}d\mathrm{\Gamma }$$

The system is discretized with inf-sup stable Q2–Q1 Taylor–Hood elements (fe_order_v = 2, fe_order_p = 1 in the [fe] section; see Discretization). ${\eta }_{\mathrm{VEP}}$ enters as a piecewise-constant CellField projected from the markers (project_markers_to_cells! in src/particles/Projection.jl). assemble_stokes wraps the forms into a Gridap AffineFEOperator.

Free-slip penalty

Walls flagged free_slip (bc_left, bc_right, bc_surface, bc_deep, plus bc_front/bc_back in 3D) contribute the penalty term

$${\int }_{\mathrm{\Gamma }}\gamma (\mathbf{u}\cdot \mathbf{n})(\mathbf{v}\cdot \mathbf{n})d\mathrm{\Gamma },$$

which weakly enforces no-penetration ($\mathbf{u}\cdot \mathbf{n}=0$) while leaving the tangential component free. build_free_slip_bcs sets the penalty scale per wall as

$$\gamma =\mathtt{\text{free\_slip\_penalty}}\cdot \frac{{\eta }_{\mathrm{ref}}}{h_x},$$

where ${\eta }_{\mathrm{ref}}$ is the volume-mean viscosity of the current iteration (floored at eta_min) and $h_x$ the cell width — so the dimensionless config key free_slip_penalty (default 1000) is portable across viscosity regimes and resolutions.

Traction (Neumann) walls

Walls flagged traction impose $\mathit{\sigma }\cdot \mathbf{n}=\mathbf{t}$ via the RHS boundary integral ${\int }_{\mathrm{\Gamma }}\mathbf{t}\cdot \mathbf{v}d\mathrm{\Gamma }$. build_traction_bcs composes the global-frame traction vector from the per-wall config keys t_traction_<wall>_n (normal, positive outward) and t_traction_<wall>_t (tangential), where <wall> is one of left | right | surface | deep. Walls flagged prescribed use strong Dirichlet values instead (src/boundary_conditions/Velocity.jl).

VEP elastic memory source

For Maxwell visco-elasticity the previous-step deviatoric stress enters the RHS as

$${\int }_{\mathrm{\Omega }}\chi {\mathit{\sigma }}_{\mathrm{old}}:\epsilon (\mathbf{v})d\mathrm{\Omega },\chi =1-Z,$$

where $Z=G\mathrm{\Delta }t/({\eta }_{\mathrm{VP}}+G\mathrm{\Delta }t)$ is the visco-elasticity factor and $\chi$ the elastic-memory weight (derivation in Rheology — Maxwell elasticity). ${\mathit{\sigma }}_{\mathrm{old}}$ is a cell-piecewise SymTensorValue field built from the marker-averaged stress components by build_sigma_old_cellfield. picard_stokes! only constructs this term when any cell has $\chi >0$ (i.e. some rock has shear_modulus $G>0$); with $\chi =0$ everywhere the system reduces to pure visco-plastic Stokes.

After each Stokes solve, update_marker_stresses! applies the marker-level stress update (Gerya 2010, eq. 12.13)

$${\mathit{\sigma }}_{\mathrm{new}}=2{\eta }_{\mathrm{VEP}}\dot{\epsilon }(\mathbf{u})+\chi {\mathit{\sigma }}_{\mathrm{old}}$$

followed by an explicit Jaumann co-rotation increment $\mathrm{\Delta }t(\mathbf{W}\mathit{\sigma }-\mathit{\sigma }\mathbf{W})$ with spin tensor $\mathbf{W}=\frac{1}{2}(\mathrm{\nabla }\mathbf{u}-\mathrm{\nabla }{\mathbf{u}}^{\mathsf{T}})$; the 3D variant _update_marker_stresses_3d! carries all six deviatoric components and the full three-component spin. recover_invariants! computes cell-mean ${\dot{\epsilon }}_{II}=\sqrt{\frac{1}{2}\epsilon :\epsilon }$ and ${\sigma }_{II}=2\eta {\dot{\epsilon }}_{II}$ for diagnostics and for the next nonlinear iteration.

Buoyancy: Boussinesq perturbation form

buoyancy_cellfield (src/physics/Gravity.jl) builds the body-force CellField used in $l$. By default the perturbation form is used:

$${\mathbf{f}}_{\mathrm{buoy}}=(\rho -\overline{\rho })\mathbf{g},$$

where $\overline{\rho }$ is the volume-weighted mean cell density of the current Picard iteration (a plain arithmetic mean would bias toward refined regions on AMR meshes and produce spurious buoyancy). The hydrostatic component $\overline{\rho }\mathbf{g}$ is absorbed into a reference pressure that is dropped from the system entirely — this prevents direct-solver pivoting on the singular saddle point from leaking the hydrostatic gradient into the velocity field (observed as ${10}^4\times$ overestimated plume velocities before the fix; see the commentary in src/time_stepping/PicardLoop.jl).

Options:

Config key	Default	Effect
body_force_full	false	true → use full $\rho \mathbf{g}$ (I2ELVIS/FD convention), no subtraction
buoyancy_from_T	false	true → build $\rho (T_h)={\rho }_0(1-\alpha (T_h-T_{\mathrm{ref}}))$ from the continuous FE temperature (resolves thermal boundary layers that the cell-constant marker projection under-drives at high Ra); $T_{\mathrm{ref}}$ = T_ref_density
gravity.model	"vertical_const"	"vertical_const" (constant $(g_x,g_y[,g_z])$ from gx,gy,gz) or "radial" ($-g_{\mathrm{mag}}\hat{\mathbf{r}}$ toward the origin, magnitude g_mag) — gravity_fn in src/physics/Gravity.jl

ALA / TALA continuity

Setting formulation = "ala" or "tala" (in the [rheology] section) replaces incompressibility with the anelastic mass balance

$$\mathrm{\nabla }\cdot ({\rho }_{\mathrm{adi}}\mathbf{u})=0,$$

which build_stokes_forms expands by the product rule into the test equation

$$-q({\rho }_{\mathrm{adi}}\mathrm{\nabla }\cdot \mathbf{u}+\mathbf{u}\cdot \mathrm{\nabla }{\rho }_{\mathrm{adi}}).$$

picard_stokes! constructs the 1-D adiabatic reference profile (adiabatic_density_profile in src/physics/Compressibility.jl)

$${\rho }_{\mathrm{adi}}(y)={\rho }_{\mathrm{surface}}{e}^{y/H}$$

(exponential with depth, $y$ down) from the config keys rho_ref_ala (default 3300) and the scale height H_scale (default 2.9e6 m), and passes it as the ρ_adi_cell kwarg to assemble_stokes. TALA (truncated ALA) uses the same continuity but drops the dynamic-pressure contribution to density. Valid values of formulation are checked by ala_continuity_modification: boussinesq | ala | tala | isothermal_compression. A linearized adiabatic temperature profile helper (adiabatic_temperature_profile, $T(y)=T_s(1+\alpha gy/c_p)$) is also provided.

Nonlinear (Picard) iteration

The viscosity depends on strain rate, $\eta =\eta ({\dot{\epsilon }}_{II})$, so the Stokes operator is nonlinear. picard_stokes! (src/time_stepping/PicardLoop.jl) iterates:

1. project markers → cell fields ($\eta$, $\rho$, $\chi$, ${\mathit{\sigma }}_{\mathrm{old}}$),

2. assemble and solve the linear Stokes system with frozen ${\eta }^{(k)}$,

3. recover ${\dot{\epsilon }}_{II}^{(k+1)}$,

4. test $\Vert {\mathbf{u}}^{(k+1)}-{\mathbf{u}}^{(k)}\Vert /\Vert {\mathbf{u}}^{(k+1)}\Vert <\mathtt{\text{picard\_tol}}$.

Controls in [fe]: n_picard (default 3), picard_tol (default 1e-3), nonlinear_method = "picard" | "newton" and newton_damping ($\alpha \in (0,1]$ under-relaxation; $\alpha =1$ recovers Picard). Sparse patterns, factorizations, and multigrid hierarchies are cached across iterations — see Linear solvers and Geometric multigrid.

Related pages

• Rheology — how ${\eta }_{\mathrm{VEP}}$, $Z$, $\chi$ are computed

• Heat transport — the coupled energy equation

• Two-phase flow — the three-field $(u_s,p_f,p_c)$ extension

• Free surface — sticky-air and ALE treatments

• Parameter file — full key reference