Markers

MILEVA carries composition and material history on Lagrangian markers (tracers) that move through the Eulerian grid. In the 2-D default configuration markers also carry temperature; see Temperature transport for the alternative of advecting temperature on the grid. The code lives in markers/: init_markers.jl (storage and initial placement), advect_markers.jl (Runge-Kutta advection), markers_to_grid.jl (projection), interp_grid_to_markers.jl (hand-back), subgrid_diffusion.jl and reseed_markers.jl.

Storage

Markers (markers/init_markers.jl) is a struct of plain vectors, one per field, with the position in polar or spherical coordinates:

Field	Meaning
r, th, ph	radius (m), angle or colatitude, longitude (3-D only)
T	temperature (K)
d	grain size (m)
C, pile	continent and thermochemical-pile fractions
basalt, depletion	basalt fraction and melt depletion
eps_acc	accumulated strain
erupt_time	time basalt was last deposited at the surface
hpe	heat-producing elements relative to primitive mantle
volatile	passive volatile concentration (ppm), no feedback
w	persistent weight: the area (2-D) or volume (3-D) the marker stands for
cstage, cmox, cmred, root0	per-parcel state for driver hooks, inert unless a hook writes them

MARKER_VECTOR_FIELDS lists every per-marker vector. Any code that culls, subsets or appends markers must touch all of them together, because the hot loops run with @inbounds and a field left behind desynchronises the set.

Left: the 17 per-marker vectors of MARKER_VECTOR_FIELDS, each a Float64 vector, so a 2-D marker costs 136 bytes (3-D adds the longitude). Right: memory of the marker set at the default density, 9 markers per cell in 2-D ($N_t=4N_r$) and 8 per cell on both Yin-Yang patches in 3-D ($N_{\theta }=N_r$). A 2-D run at $N_r=64$ carries $1.5\times {10}^5$ markers in 20 MB; the 3-D set passes 16 GB near $N_r=128$. Script: docs/scripts/fig_numerics_markers_sections.jl.

Initial placement and weights

init_markers places cfg.markersPerCell markers per cell, rounded to a square $n^2$ in 2-D (default 9 gives 9) and a cube $n^3$ in 3-D (default 9 gives 8). Each cell is split into an $n\times n$ (or $n^3$) sub-grid and one marker is placed in each sub-cell at a random offset within its central 80 %, so coverage is uniform with no lattice artefacts. Placement uses the stretched radial faces (grid.rF, grid.drVec).

Each marker starts with weight $w=V_{\text{cell}}/n_{\text{pc}}$, where the cell measure is exact for the geometry (cell_area_2d in markers/init_markers.jl, which calls cell_measure in utils/geometry.jl; $\int r^{m}dr\mathrm{\Delta }\theta$ in 2-D, the spherical cell volume in 3-D). Every composition inventory in MILEVA is $\sum _m{w}_m{f}_m$, so this weight is the measure the whole marker scheme conserves against.

Markers from init_markers in an 11-cell (31°) sector of a $16\times 128$ annulus with the default 9 markers per cell, coloured by their initial temperature. Each cell holds one marker per sub-cell at a random offset, so coverage is even without a visible lattice.

Advection

advect_markers! integrates $d\mathbf{x}/dt=\mathbf{V}(\mathbf{x})$ with classical fourth-order Runge-Kutta (cfg.advectionScheme = "RK4", default) or the midpoint rule ("RK2"). In 2-D, with $\mathbf{V}=(V_r,V_{\theta })$,

$${\mathbf{k}}_s=(V_r,\frac{V_{\theta }}{r}){|}_{{\mathbf{x}}_s},{\mathbf{x}}_{n+1}={\mathbf{x}}_n+\frac{\mathrm{\Delta }t}{6}({\mathbf{k}}_1+2{\mathbf{k}}_2+2{\mathbf{k}}_3+{\mathbf{k}}_4),$$

where the stages are evaluated at ${\mathbf{x}}_n$, ${\mathbf{x}}_n+\frac{\mathrm{\Delta }t}{2}{\mathbf{k}}_1$, ${\mathbf{x}}_n+\frac{\mathrm{\Delta }t}{2}{\mathbf{k}}_2$ and ${\mathbf{x}}_n+\mathrm{\Delta }t{\mathbf{k}}_3$. The angle advances with the angular velocity $V_{\theta }/r$ and is wrapped into $[0,2\pi )$; in 3-D $\dot{\phi }=V_{\phi }/(r\sin \theta )$. Velocities are interpolated bilinearly (trilinearly in 3-D) from their own staggered positions (interp_velocity, interp_velocity_3d in utils/interp_velocity.jl), stretched-grid aware through searchsortedfirst on grid.rF and grid.r. Final radii are clamped 10 m inside the walls.

In 3-D, markers that leave the Yin patch's parametric rectangle after the step are converted into the other patch's coordinates with the Kageyama-Sato transform (_yinyang_exchange_markers!); the very rare marker that lands outside both rectangles at a corner is clamped back inside and counted.

Left: 34 markers carried by advect_markers! (RK4) for 600 steps at the default Courant number through one fixed two-cell Stokes flow on a $64\times 256$ annulus; orange dots are the starting points. Right: rms position error of 400 markers after 1.2 Myr against an RK4 reference with 4096 steps. The midpoint rule converges at second order; RK4 is about 20 times more accurate at the default cflSafety, and its order is reduced below four by the piecewise-bilinear velocity interpolation. Script: docs/scripts/fig_numerics_markers_sections.jl.

Projection to the grid

markers_to_grid scatters marker values to cell centres with bilinear (2-D) or trilinear (3-D) weights, scaled by each marker's weight:

$$f_c=\frac{\sum _m{\hat{w}}_m{\varphi }_c({\mathbf{x}}_m)f_m}{\sum _m{\hat{w}}_m{\varphi }_c({\mathbf{x}}_m)},{\hat{w}}_m=w_m/\bar{w}.$$

The $w_m$ factor makes the projection a measure average rather than a marker-number average, which matters because reseeding leaves weights that can span a factor of 60 or more. Dividing by the mean weight keeps the accumulated sums of order one, so the empty-cell threshold keeps its meaning. A node that receives no marker weight falls back to a supplied field (the conductive or adiabatic profile for temperature). The accumulated weight per node, un-normalised, is kept in MarkerScatterWS.W for the hand-back below. Threads accumulate into private buffers that are summed afterwards.

A smooth field $f=1500+300\cos 4\theta \sin (\pi z/D)$ K set on the markers of init_markers and projected with markers_to_grid on a $64\times 256$ annulus. Left: the rms difference from $f$ at the cell centres, 0.34 K at the default 9 markers per cell, falling roughly as $n^{-1/2}$ from 0.68 K at 4 to 0.19 K at 36. Right: the error at 9 markers per cell, mostly cell-scale sampling noise with a weak imprint of the field's curvature. Script: docs/scripts/fig_numerics_markers_sections.jl.

Handing the temperature increment back

When temperature lives on markers, each step projects marker temperature to the grid, solves advection-free diffusion and heating there, and returns only the increment $\mathrm{\Delta }{T}_c=T_c^{\text{new}}-T_c^{\text{proj}}$ to the markers. Doing that conservatively takes two pieces.

Measure-scaled delivery. Plain bilinear interpolation of $\mathrm{\Delta }T$ is not the adjoint of the weighted projection: $\sum _m{w}_m\mathrm{\Delta }T({\mathbf{x}}_m)\ne \sum _c\mathrm{\Delta }{T}_c{V}_c$. delivery_scale(grid, cfg, ws) returns $V_c/W_c$ per node, with $W_c$ the accumulated marker measure from the last projection, and main_SphCon.jl delivers $\mathrm{\Delta }{T}_c{V}_c/W_c$. That makes

$$\sum _m{w}_m\mathrm{\Delta }{T}_m=\sum _c\mathrm{\Delta }{T}_c{V}_c$$

an identity: the markers gain exactly the energy the grid solver computed. The factor is clamped to $[1/s_{\text{max}},s_{\text{max}}]$ with $s_{\text{max}}$ = cfg.deliveryScaleMax (default 2). The cap keeps a sparsely populated wall cell (for example $W_c/V_c=0.185$) from handing its few markers several times the cell's increment. Energy withheld by the cap shows up in the budget residual (see Energy budget), not in the velocity. Marker temperatures are then clamped to the same bounds the grid solver uses.

Sub-grid diffusion. A grid increment moves a cell's mean temperature but leaves the spread of marker temperatures inside the cell alone, so that spread has no dissipation channel and can grow into cell-scale buoyancy noise. subgrid_temperature_correction! (Gerya & Yuen 2003) relaxes each marker toward the grid temperature at its position on the local diffusion time,

$$\delta T_m=(T_{\text{grid}}({\mathbf{x}}_m)-T_m)(1-e^{-d\mathrm{\Delta }t/{\tau }_m}),{\tau }_m=\frac{1}{\kappa (2/\mathrm{\Delta }{r}^2+2/(r\mathrm{\Delta }\theta {)}^2)},$$

with $d$ = cfg.subgridDiffusivity (default 1; 0 disables). The heat this moves, $\sum w_m\delta T_m$ over the markers of each cell divided by the cell measure, is subtracted from the grid increment before delivery, so the step still conserves energy exactly. This correction is 2-D only.

tests/test_marker_energy.jl asserts the closure: one marker temperature step conserves energy to $3\times {10}^{-15}$ in 2-D (both geometries, sub-grid diffusion on) and $5\times {10}^{-11}$ in 3-D. tests/test_quiescence.jl records the price: exact conservation turns a uniform heat source into about 4 % lateral structure in the warming (converging with $N_r$), where the old unscaled delivery showed 0.2 % and lost energy; the sub-grid relaxation is the intended damping of that noise.

Left: the energy a grid increment hands to the markers, relative to the energy the grid computed, for five increments on a $64\times 256$ annulus. On the markers of init_markers (circles) measure-scaled delivery is exact to ${10}^{-13}$ while plain bilinear delivery misses by up to ${10}^{-3}$. On a marker set stirred for 600 steps and reseeded every 5 (triangles), 23 of 16384 nodes have a factor beyond deliveryScaleMax = 2, mostly in the wall rows; the cap then withholds up to $4\times {10}^{-3}$ of a boundary-layer increment, which is what budget_resid records, while plain delivery is wrong node by node and its total errs in both directions. Right: subgrid_temperature_correction! applied to markers carrying 50 K of random scatter about a smooth field; the scatter decays as $e^{-dt/\tau }$, with $\tau$ = 28 Myr for the mid-depth cells. Script: docs/scripts/fig_numerics_markers_sections.jl.

Reseeding

Divergent flow empties cells and convergent flow crowds them. reseed_markers! runs every cfg.remeshEvery steps (default 5, when cfg.useAdaptiveRemesh is true, the default) after advection, and is exactly conservative for every carried field:

• Death. A cell with more than cfg.reseedNmax markers (default 36) keeps the first reseedNmax; the removed weight is shared equally among the survivors, and the survivors' values are shifted so that $\sum wf$ over the cell is unchanged. Bounded fields stay inside their bounds by redistributing any overshoot among the other survivors.

• Birth. A cell with fewer than cfg.reseedNmin markers (default 4) gets new markers at jittered positions. Weight is redistributed, not created: all markers of the cell share the old total equally. New markers start at the weighted cell mean and the cell is then shifted as above.

• Empty cells. Only a cell with no marker at all creates weight, with values taken from the grid fields passed in. This is the one non-conservative operation left in the marker scheme.

Births inherit the cell-averaged basalt and depletion passed in by the caller, so reseeding does not dilute composition. In 2-D, cell membership is decided by one shared helper, marker_cell_index_2d, from the cell faces, so reseeding and the composition binning in the physics routines agree.

The initial markers of a $64\times 256$ annulus carried for 1500 steps through one fixed two-cell flow. Left: without reseeding 1860 cells fall below reseedNmin and 112 empty completely; with reseed_markers! every 5 steps no cell ends below 4. Of the 300 reseeding calls, 11 met an empty cell; in the others the total weight, $\sum wT$ and $\sum w{f}_{\text{basalt}}$ changed by at most $3\times {10}^{-16}$ relative. Right: reseeding redistributes weight, so the weights spread over a factor of about 300, which is why the projection weights each marker by $w_m$. Script: docs/scripts/fig_numerics_markers_sections.jl.

Markers in 3-D

In 3-D, _init_markers_3d creates a Yin and a Yang copy of the marker set at the same parametric coordinates. Advection (_advect_3d!) and projection (markers_to_grid for MutableGrid3D) use the Yin patch fields, and the patch exchange converts coordinates when a marker leaves the parametric rectangle. Yang's temperature is advected on its own grid (see Temperature transport).

3-D markers

The 3-D marker path is less exercised than the 2-D one. The sub-grid diffusion step is 2-D only, and the delivery cap withholds about 7 % of the increment at Yin's edge nodes, which have half-cell support. The 3-D example advects temperature on the grid instead.

Left: the bottom layer of the 3-D marker set from init_markers on a $4\times 12\times 36$ patch, 13824 markers per patch at 8 per cell, every second one drawn: the Yin copy (black) and the Yang copy at the same parametric positions, mapped into the global frame with yy_yang_to_yin (blue). Right: after shifting the Yin markers 30° in longitude, the 1536 that left the patch rectangle (orange, Yin coordinates) are converted by _yinyang_exchange_markers! into Yang coordinates (blue); mapping them back reproduces their positions to $9\times {10}^{-16}$ rad. Script: docs/scripts/fig_numerics_markers_sections.jl.

Where in the code

Topic	File	Function
Storage, initial placement	markers/init_markers.jl	Markers, init_markers, MARKER_VECTOR_FIELDS
Advection	markers/advect_markers.jl	advect_markers!, _advect_2d!, _advect_3d!, _yinyang_exchange_markers!
Velocity interpolation	utils/interp_velocity.jl	interp_velocity, interp_velocity_3d
Projection	markers/markers_to_grid.jl	markers_to_grid, MarkerScatterWS
Hand-back	markers/interp_grid_to_markers.jl	interp_grid_to_markers!, delivery_scale
Sub-grid diffusion	markers/subgrid_diffusion.jl	subgrid_temperature_correction!
Reseeding	markers/reseed_markers.jl	reseed_markers!, marker_cell_index_2d
Cell measures	utils/geometry.jl, markers/init_markers.jl	cell_measure, cell_area_2d, cell_volume_3d
Tests	tests/	test_markers.jl, test_marker_energy.jl, test_projection_noise.jl, test_quiescence.jl, test_depletion_projection.jl