[chaseshyu]

Why

DES's subduction rules, phase_change_option = 1, hard-code materials 0–7 and serpentinize only
elements with hydrous markers. geoflac's rules (eclogite, mantle
⇄ serpentinite, chlorite-out, the sediment chain) had no DES port. They serpentinize mantle above
subducted wet crust (oceanic crust, sediment or sedimentary rock), which needs a marker-cell index
that follows the slab and that each checkpoint carries. This PR adds the index, VirtualTopology,
and the rule engine on it, phase_change_option = 3, with depth measured below the current surface.

What

31 commits, each message with its why and what; (10) is the 10th, in the order of the Commits tab.
C++: +1886 −101, mostly phasechanges.cxx and mesh.cxx/mesh.hpp. A regression gate for option 3
follows as its own PR.

Part	Commits	Content	Changes runs without option 3
1	1–6	Fixes of code that predates this PR: phase-change input check (1), option 1 without hydration (2), averaged output (3, 4), lumped masses after phase changes (5), compare.py marker z in 3D (6)	3, 4 (output only), 5 (options 1, 101), 6
2	7–9	Two keys renamed (7), depth below the surface without the attenuation table (8), VirtualTopology: a marker-cell index on a background grid, checkpointed (9)	7 (key names)
3	10–11	The rule engine (10): a cell is wet crust above half its markers; examples/subduction.cfg (11)	—
4	12–17	Speed-ups (12–15), no lock (16), marker passes on the GPU (17); bit-identical on CPU	—
5	18–28	Review revisions: the "wet crust" name (18); surface plastic strain decays in every top element (19) and the eruption's leftovers go (20); unset or shared material indices stop at input (21); a restart with another grid cell size rebuilds the grid (22); the grid is rebuilt at each quality check (23), stored in one array (24) and ends below the deepest mantle that can serpentinize (25, 26); phase changes run before update_mesh() (27); a temporary repair of the bottom-snap fold (28)	19 (2D), 21 (input), 27 (options 1, 101), 28 (only a run that would stop)
6	29–30	Second review round: the pressure the rules read is a key (29); the serpentinization depth cap and the wet-crust fraction are keys (30)	—
—	31	The CHANGELOG entry	—

Breaking:

• mat.mattype_crust → mat.mattype_continental_crust, mat.mattype_sed → mat.mattype_sediment;
  an old name stops at input with exit 11 (7).
• mat.mattype_continental_crust and mat.mattype_oceanic_crust default to −1 (unset):
  ic.temperature_option = 1, or an active option-3 rule that reads one, stops at input with
  exit 11 until it is set (21).
• In 2D, surface plastic strain decays in every top element; it skipped oceanic crust, a remnant of
  the mid-ocean-ridge eruption that Enable GPU acceleration and parallelization across Dynearthsol #25 removed. 2D runs with oceanic crust at the surface change (19).
• Averaged output (the default) divides by the steps of its window: about 1 % at the default
  mesh.quality_check_step_interval; with an interval of 1 it is instantaneous (3, 4). Physics
  unchanged.
• Runs with phase_change_option ≠ 0: conversions read the positions at the start of the step,
  and the lumped masses are rebuilt from the converted composition (5, 27). Options 1 and 101
  change from their first conversion; on a static mesh (control.has_moving_mesh = no,
  control.rsf_dtheta_max = 0), from step 10.

For the reviewer

• Feedback wanted: a check of the review revisions (parts 5 and 6), and a physics check of the
  rule table. The three open questions are answered:
  1. Pressure: selectable, mat.phase_change_pressure_option (29). 0 (default) = the element mean
     stress, 1 = ref_pressure() at the depth below the local surface.
  2. Constants: keys, mat.serpentinization_max_depth (150e3; it also ends the grid) and
     mat.wet_crust_fraction_threshold (0.5) (30).
  3. Static-mesh mass rebuild: kept, because density changes with temperature too.
• Start at: phasechanges.cxx, subduction_phase_changes() and build_phase_transition_rules();
  then the VirtualTopology overview in mesh.hpp.
• Check closely:
  1. antigorite_stable() bounds both serpentinite rules; its low-pressure point (650 °C, 0.2 GPa)
     replaces geoflac's (670 °C, 0.6 GPa) in the mantle rule. With geoflac's two fields, mantle at
     642–730 °C between the low-pressure lines serpentinizes and dehydrates back at every call; one
     field leaves no pressure and temperature where both rules apply.
  2. The grid is rebuilt at each quality check, right after update_mesh(): markers are added or
     removed only there and in a remesh, which rebuilds it too, so the grid indexes every marker.
     An assert stops the run (exit 60) if a marker count or id changes between rebuilds. A remesh
     step builds the grid twice.
  3. (28) is temporary: the full fix, bc.has_winkler_bottom_restore, comes in a later release.
     It acts only when MMG returns a mesh the existing tangle check would stop on.
  4. The OpenACC passes of (17) have no data clauses: they rely on mem:managed, which every
     openacc=1 build sets.

Try it

• Default: phase_change_option = 0. With 3, a rule runs when its mat.has_* switch is on; its
  source and product must then name distinct, set materials, or the run stops at input (exit 11).
  Serpentinization (on) needs mat.mattype_serpentinite and mat.mattype_sediment;
  eclogitization (on) needs mat.mattype_oceanic_crust and mat.mattype_eclogite.
• Keys that change it:
  • mat.serpentinization_thickness (10 km): raise it for a thicker serpentinite band.
  • mesh.vt_size (0.5 × mesh.resolution): lower it for a sharper slab outline, while cells still
    hold markers.
  • mat.serpentinization_wet_crust_margin (−1, off): set it ≥ 0 to drop false positives under a
    slab bend.
  • mat.phase_change_pressure_option (0, the element mean stress): set it to 1 for the
    reference lithostat that the phase diagrams use.
  • mat.serpentinization_max_depth (150 km) and mat.wet_crust_fraction_threshold (0.5):
    untuned.
• Start from: examples/subduction.cfg (2D, usemmg=1; 5 Myr in about 22 min on 8 threads of a
  TACC Vista node); it lists every option-3 key, most with a comment.

Numerical effect

• Effect: intended. Option 3 is new: (12)–(17) keep its results (CPU; on a GPU, (17) matches
  (16)); part 5 changes them as the review asked; part 6 keeps them bit-exact at the defaults of
  its new keys. Runs without option 3 change only as Breaking lists; option 0 physics is bit-exact but for the 2D surface plastic strain of (19).
• Reproduce: make ndims=2 usemmg=1, then in examples/ OMP_NUM_THREADS=8 ../dynearthsol2d subduction.cfg: frame 100 (5 Myr) holds 3292 serpentinite (material 4) and 392 eclogite
  (material 5) markers (g++ 14.2.0 on Vista). Another compiler gives other counts.
• Cost: phase changes take 2.79 % of a 21 min 34 s run of examples/subduction.cfg on 8 threads
  of a TACC Vista node; (27) cut the time per step by 3 %. On a busy Apple M3 Max with 6 threads,
  the same run takes 28 min 56 s, and phase changes take 4.43 %.

Cost per part, and the host

From the end-of-run time summary (also in the .manifest): 5.78 M steps, 220 remeshes, 1294 s in
all. Option 3 runs every 10 steps.

Part	Time	Share of the run	Per call
Phase changes, both marker passes	36.1 s	2.79 %	0.062 ms, 578 000 calls
Remesh with MMG, including its grid rebuild	7.2 s	0.55 %	33 ms, 220 remeshes
Output, including the grid in checkpoints	6.1 s	0.47 %	60 ms, 101 frames
Grid rebuild at each quality check	no timer of its own: inside Compute (96.19 %)	—	5780 calls

Host, from the run's .manifest: a TACC Vista CPU node (NVIDIA Grace, Arm Neoverse-V2; Rocky Linux
9.7; g++ 14.2.0), 8 OpenMP threads pinned to one 8-core block (OMP_PLACES=cores,
OMP_PROC_BIND=close); other runs used other cores of the node. (27): 214 → 207 µs per step, two
repeats each, on the same node.

The same run on an Apple M3 Max (12 performance + 4 efficiency cores; macOS 26.6; Apple clang 21),
6 OpenMP threads, not pinned: 28 min 56 s, 5.89 M steps, 279 µs per step. Phase changes took 4.43 %
(0.131 ms per call), remeshing 0.47 % (200 remeshes) and output 0.38 %. Other jobs shared the Mac
during the run: the 1-minute load average was 14–20 on its 16 cores (mean 16), so an idle M3 Max
runs faster. The step and remesh counts differ from Vista's because another compiler gives other
round-off.

Checks of the second review round (29)–(30)

Check	Result
(29) alone and the tip, ndims 2/3, usemmg=1 hdf5=1, Apple clang; the tip also with usemmg=0	builds
Tip on the follow-up PR's gate (2D), against pins of (28)'s binary	pass; all compared fields bit-exact. Two new items check pressure option 1: the hot case at the edge of the antigorite field, with the pressure recomputed as the lithostat (15 band markers, 82 serpentinite, as option 0), and gravity off, where the lithostat is 0 (no serpentinite; option 0 makes 133)
The two new items against broken builds	a build that negates the depth leaves 133 wedge mantle markers; a build whose option 1 still reads the mean stress makes 133 serpentinite; each fails its item
A pressure option other than 0 or 1, a depth cap ≤ 0, a fraction outside [0, 1)	each stops at input, exit 11
examples/subduction.cfg to 1 Myr, pressure option 0 against 1 (Apple M3 Max, 6 threads)	serpentinite 704 and 741, eclogite 45 and 23 markers; one pair of runs, which diverge

Not run for (29)–(30): g++, nvc++ openacc=1, and the 3D gate.

Checks of the review revisions (18)–(28)

Check	Result
Each commit, ndims 2/3, usemmg=1 hdf5=1, Apple clang ((25) in 3D together with (26); (28) in 3D with g++ 14.2.0 only)	builds
Each commit with its own CPU binaries, on the follow-up PR's gate (2D), now also with restarts at another grid cell size, input checks and the grid's depth limit	pass; option 0 cases (test-tiny, test-topo, test-core, test-rect-tiny) bit-exact but at (19), which changes 2D surface plastic strain as intended (test-rect-tiny: two top elements; test-core); option-3 results re-pinned only where the commit intends a change
Tip, g++ 14.2.0 on Vista: both gates (2D and 3D)	pass
(25), (27), (28): examples/subduction.cfg to 5 Myr on Vista, against the commit before	(25): serpentinite −27 % at 250 kyr and −7 % at 5 Myr, as the grid rows now sit at fixed depths below the surface; (27): within a few % up to 3.5 Myr; (28): the run without it stops at 3.7 Myr (exit 42), with it reaches 5 Myr, and a restart from 3.5 Myr matches the fresh run exactly

Checks of the first round (1)–(17)

Check	Result
Every commit, ndims 2/3, usemmg=1 hdf5=1, Apple clang and gcc 8.5	builds
(10)–(17), nvc++ openacc=1, ndims 2/3	builds
(1)–(9): test-rect-tiny	bit-exact against master
(10)–(17), each with its own CPU binaries, on the follow-up PR's gate. 2D: test-rect-tiny without option 3; a small subduction model with its grid checked against the markers in every checkpoint, restarts after a remesh, with the crust margin on and from a conversion step, rule checks, a hot case at the edge of the antigorite field, and gravity off to check the pressure source. 3D: grid checks and a restart across a remesh	pass; all compared fields bit-exact against one reference set
(17) against (16) on a GPU (RTX A5000): examples/subduction.cfg to 100 and 500 kyr	bit-exact frames

Demo: examples/subduction.cfg to 5 Myr (g++ 14.2.0 on Vista, 8 threads, 22 min)

The slab at 0-4 Myr: the mesh with each element filled by its main material, isotherms, the 3 km
VirtualTopology grid, its wet-crust cells (orange: more than half oceanic crust, sediment or
sedimentary rock) and the cells the 5 km search reaches. As the slab sinks to about 110 km, mantle
serpentinizes just above its wet-crust cells, along the slab top and in the wedge corner, and
eclogite forms at its tip.

At 5 Myr. Top left: the same map. Right: a zoom on the slab top with the markers, and its plastic
strain with the same isotherms. Bottom, through the run: the element count with each remesh's
change (the jump at 3.7 Myr is the remesh (28) repairs), the 220 remeshes by trigger (bottom drift
early, then low quality and tiny elements), and the serpentinite and eclogite markers.

Demo before the review revisions: examples/subduction.cfg to 5 Myr (Apple clang, 6 threads, 20 min)

The slab at 0-4 Myr: the mesh with each element filled by its main material, isotherms, the 3 km
VirtualTopology grid, its crust cells (orange: more than half oceanic crust or sediment) and the
cells the 5 km search reaches. As the slab sinks to about 110 km, mantle serpentinizes just above
its crust cells, along the slab top and in the wedge corner, and eclogite forms at its tip.

At 5 Myr. Top left: the same map. Right: a zoom on the slab top with the markers, and its plastic
strain with the same isotherms. Bottom, through the run: the element count with each remesh's
change, the 206 remeshes by trigger (bottom drift early, then low quality and tiny elements), and
the serpentinite and eclogite markers.

New or retired names

• Added: option 3 and its 16 keys (list below); a grid record in option-3 checkpoints (a corrupt one
  stops a restart with exit 22; a checkpoint without one still restarts, and a restart with another
  grid cell size warns and rebuilds the grid); a PhaseChg: line in the time summary and
  phase_change_* fields in the manifest (option 3). No new exit codes.
• Renamed or retired: the two keys under Breaking. --help drops option 2, the mid-ocean-ridge
  rule that Enable GPU acceleration and parallelization across Dynearthsol #25 removed, and marks its four keys unused; they are still accepted (1).

The 16 keys of option 3 with defaults, and the four unused keys

• Materials, −1 = unset (a rule on that reads an unset one stops at input):
  mat.mattype_serpentinite, mat.mattype_eclogite, mat.mattype_hydrous_mantle,
  mat.mattype_sediment_rock, mat.mattype_metasediment, mat.mattype_schist.
• Switches: mat.has_serpentinization (yes), mat.has_eclogitization (yes),
  mat.has_hydrous_mantle_dehydration (no), mat.has_sediment_metamorphism (no).
• Pressure: mat.phase_change_pressure_option (0, the element mean stress; 1, ref_pressure()).
• Grid: mat.serpentinization_thickness (10e3 m), mat.serpentinization_wet_crust_margin (−1),
  mat.serpentinization_max_depth (150e3 m), mat.wet_crust_fraction_threshold (0.5),
  mesh.vt_size (0.5). It counts mat.mattype_oceanic_crust, mat.mattype_sediment and
  mat.mattype_sediment_rock (when set) as wet crust.
• Unused since Enable GPU acceleration and parallelization across Dynearthsol #25, still accepted: mat.mattype_depleted_mantle, mat.mattype_partial_melting_mantle,
  mat.mattype_mor_extrusion, mat.convert_rate_oceanic_crust.

Checklist

• Follows CONTRIBUTING.md; builds for ndims=2 and ndims=3:
  yes, every commit.
• Entry under ## [Unreleased] in CHANGELOG.md: yes (31).
• New parameters registered in input.cxx and documented in examples/defaults.cfg: yes.
• Test under tests/functional/ or benchmarks-cores/: in the follow-up PR.

AI disclosure

Claude Opus 5.5 co-authored every commit, ran the checks and subagent reviews, restructured the
history and drafted this text.

🤖 Generated with Claude Code

[chaseshyu]

Hi @tan2,

Thank you for the careful review.
It made the feature more correct and efficient in some way.

Here is what changed:

Run time. The subduction.cfg header now names the machine. The 5 Myr run takes 21 min 34 s on TACC Vista (NVIDIA Grace, g++ 14.2.0, 8 threads), and 28 min 56 s on an Apple M3 Max with 6 threads. The Mac was busy with other jobs during the run, so an idle one is faster.

Following are the revision regarding each comment:

1. Naming. "Crust" is now "wet crust" everywhere, and the parameter is now serpentinization_wet_crust_margin.

2. Stale grid. The grid is now rebuilt at every quality check. Markers change only at those steps, so the grid is never out of date for long. This let me delete the removed-id set and the unlabelled-marker fallback. Restarts are still bit-exact.

3. Restart. If vt_size or resolution changes, the restart now warns and rebuilds the grid.

4. Memory and speed.

• The grid now uses one flat array instead of one vector per cell.
• The grid now stops below the deepest mantle that can serpentinize.
• phase_changes() now runs before update_mesh(). This removes the extra mass rebuild and makes each step about 3 % faster. I did not add the #pragma acc wait, because the two functions before it already end with one.

5. Input checks. mattype_continental_crust and mattype_oceanic_crust now default to −1. An unset index, or sediment that uses the oceanic-crust index, now stops the run at input. mattype_mantle is pushed only once. I also removed an old oceanic-crust exception in the surface plastic-strain decay, which was left over from the removed mid-ocean-ridge code.

Not in your review. The 5 Myr runs hit a known crash, where the bottom of the mesh folds. I added a temporary repair. It acts only when the mesh folds, so models that run today give the same results. The full fix, bc.has_winkler_bottom_restore, comes in a later release.

Not changed: the SurfaceTopo cache and the parallel build_index(). I measured both on Vista with examples/subduction.cfg, 1 Myr:

Part	Share of run time
SurfaceTopo builds	0.08 %
Grid rebuilds	0.17 %
Whole phase-change rule	1.77 %

In a small 3D test model (about 0.39 M markers), both are about 0.1 %. The grid is rebuilt only once per quality-check interval. Its build time and the step time both grow with the number of markers, so the share should stay this small for large models. Each change would save less than 0.2 % of the run, and a cache could also change results. If a large 3D model shows otherwise, I will speed up the slowest part of the build.

The PR description has the full checks.

AI disclosure: I wrote this reply and the revision with help from Claude Code (Anthropic's Claude Opus 5.5). I reviewed and tested all the changes before I submitted them.

[tan2]

Answers to your open questions

1. Pressure. Please make it selectable. For example, add mat.phase_change_pressure_option, with 0 = element mean stress (the current default) and 1 = ref_pressure(param, -depth), using the surface-relative depth get_ZPT already computes. The cost is the same either way: one multiply or a PREM table lookup per marker. Options 1 and 101 already
   use ref_pressure (the older get_ZPT overload), and the phase-diagram boundaries are in lithostatic terms, so option 1 is the closer match to geoflac and the papers. Mean
   stress lets tectonic over- and underpressure shift the boundaries at the slab top.

2. The constants. Yes, please make them input parameters: mat.serpentinization_max_depth (default 150e3, which also ends the grid) and mat.wet_crust_fraction_threshold
   (default 0.5). Both are labelled "untuned", and comment 1 shows how sensitive the result is to the grid.

3. The static-mesh mass rebuild. No, keep it as it is. The masses depend on more than the composition. The density
   in MatProps::rho() changes with temperature (thermal expansion), so a rebuild every 10 steps isn't a no-op even when nothing converts. Skipping it would leave stale masses in the run.

[chaseshyu]

Hi @tan2,

Thank you. I made both changes and added those three avariables:

mat.phase_change_pressure_option: 0 (default) = the element mean stress; 1 = ref_pressure() at the depth below the local surface.
mat.serpentinization_max_depth: (default 150e3; it also ends the grid)
mat.wet_crust_fraction_threshold: (default 0.5).

The static-mesh mass rebuild stays as it is. With the defaults, results are bit-identical. The PR description has the checks.

[chaseshyu]

Hi @tan2,

Because of your previous requested changes, this PR might need your approval to be merged.

This PR is ready to merge