Regular update 2026-07: remeshing robustness, restart metadata, GPU and build fixes

[chaseshyu]

Summary

Regular maintenance update (2026-07): 21 commits of bug fixes and robustness work accumulated in five groups — remeshing robustness, boundary conditions, restart/output metadata, GPU fixes, and build/tooling.

Remeshing robustness

• Boundary point deletion merged segments incorrectly — delete_points_on_boundary() (renamed delete_points_and_merge_boundary()) already fully removes the points and merges the adjacent boundary segments; the caller ran delete_points() again on the same list, a second deletion pass over an already-compacted array with stale indices. That rewrote the boundary segments of the highest-index nodes, detaching the top surface from the side wall and collapsing the corner into a spurious cliff.
• mesh.max_steiner_factor (new parameter, default 30, 2D only) — caps Triangle's Steiner-point insertion at factor × npoints during remeshing. Bounds the non-terminating encroached-subsegment splitting Triangle enters when two boundary segments meet at a tiny angle (highly distorted mesh at extreme strain), which otherwise spins forever and exhausts memory. Applied only when no area constraint is active, so initial meshing is unaffected; <= 0 disables.
• True BFS in barycentric interpolation — the enclosing-element search after remeshing only expanded one neighbor level and hard-exited at a fixed 256-element cap. Now a level-by-level BFS (capped at 3 layers, unordered_set dedup) with a best-candidate tolerance for nodes numerically on a face, and a nearest-node fallback for boundary nodes that legitimately fall just outside the old mesh. Interior nodes that still fail print a diagnostic instead of dying.
• dt refreshed unconditionally after remeshing — dt is otherwise recomputed only every slow_updates_interval steps, so a remesh that refines the mesh could run up to 10 steps on a stale, too-large dt — enough for an explicit CFL runaway to invert elements. compute_mass consumes var.dt, so the refresh comes first.
• Sediment deposition guarded against degenerate surface facets — a near-zero-area facet (collapsed "cliff" at a convergent trench) made dv_apply/base blow up, placing the marker ~10¹³ m below the domain and aborting the run (exit 168). Deposition on such a facet is skipped for the step (volume accumulator kept, so it deposits after the next remesh regularizes the surface).

Boundary conditions

• Open-sidewall lithostatic traction floored at zero — for wall facets above the datum, ref_pressure options 0/4 extrapolate linearly negative, flipping the support traction into an outward suction that pulls the topographic strip out of the domain. The true exterior there is air/vacuum. Below the datum the clamp is inert, so normal runs are bit-identical.

Restart & output metadata

• Frame files now embed the full .info row (steps, time, dt, walltime, nnode, nelem, nseg) — as FieldData in HDF5 frames and named header scalars in plain binary frames. restart() falls back to the embedded metadata when the .info file is missing, and the new utils/recreate_info.py rebuilds a deleted .info from the frame files alone. Dynearthsol.py gains scan_frames() to read the per-frame metadata directly.
• plastic strain-rate read back at restart — delta_plstrain carries per-element state across steps but was never restored, so restarted runs diverged from continuous ones.
• delta_plstrain reset each step — only the plastic branches write it, so an element that stopped yielding kept its last increment forever; the NN remap then smeared that fossil across the mesh at each remesh into a growing spurious halo. Clearing it at the top of update_stress makes "plastic strain-rate" a true per-step rate. Diagnostic-only: no feedback into the physics.
• Final averaged-fields frame reached past max_time — with is_outputting_averaged_fields the run now continues to the next averaging-aligned step so the last frame is written instead of silently dropped.
• info_display_next_step made consistent at restart, and the non-HDF5 checkpoint scalars converted from one packed 7-double array to named scalars via new BinaryOutput::write_scalar / BinaryInput::read_scaler (matching the existing HDF5 reader API).

GPU (OpenACC)

• Fixed device data management for zmin in apply_vbcs (2D): the reduction result was never made present on the device for the subsequent kernel.
• Fixed 2D regular-mesh generation under GPU builds: create_elem_from_cell / create_rect_node write host-allocated arrays and must not be ACC kernels.

Build & tooling

• macOS/clang: OPENMP_RPATH uses an absolute path instead of @executable_path, so the binary runs from any working directory (e.g. benchmarks-cores subdirectory runs).
• benchmarks-cores/Makefile: new INDIR=1 option runs each case in its own model-<CASE>-<MODELNAME> subdirectory to avoid clobbering output from other runs.
• Unused pseudo-transient (PT) solver remnants (compute_dt_PT, update_velocity_PT, var.dt_PT) commented out.
• CI workflows now clean the configuration they rebuild — clean removes only $(OBJS)/$(EXE), whose names depend on ndims and the OpenACC suffix, so the previous bare make clean only removed default-config (3D, non-GPU) objects. A make ndims=2 step after an earlier 2D build with different flags (e.g. opt=0 openmp=0) silently reused the stale objects, since flags are not encoded in object names. All build steps in basic-build, macos-build, mmg-build, and nvc-build now pass the matching ndims=/openacc=1 to their clean invocation.

Compatibility notes

• Non-HDF5 checkpoint format changed: .chkpt scalars are now individual named entries instead of one packed array, so plain-binary checkpoints written before this branch cannot be restarted by the new binary (HDF5 checkpoints are unaffected). Restart also now expects the plastic strain-rate array in the save frame (present in all existing output).
• New parameter mesh.max_steiner_factor defaults to 30 (documented in examples/defaults.cfg); set <= 0 for the old unlimited behavior.

Testing

• 2D build verified on macOS arm64 (clang++, Homebrew boost/HDF5/libomp).
• The remeshing and restart fixes were each validated on long-running 2D convergent-margin simulations (multi-Myr, thousands of remeshes) where the original failures reproduced; restart determinism was checked bit-exact against continuous runs.

🤖 Generated with Claude Code

[chaseshyu]

Hi @echoi , this PR includes the bugfixes and minor updates splitted from #74. Once this PR has been aggreed and merged, I will reopen #74 of remeshing stress remap feature.