Remeshing stress remap: free-surface pin, surface-relative reference, Deborah blend, untouched-element carry-through

[chaseshyu]

Six commits that change how element stress crosses a remesh. Stress is the one field for which remeshing is not a passive interpolation problem: it carries elastic memory in the cold lid, it must satisfy the free-surface condition exactly, and it sits on a steep lithostatic background that any smoothing operation will corrupt. This PR treats those three facts as design principles instead of accepting the remap error they otherwise produce.

Method: how stress crosses a remesh

DES remaps element stress through a superconvergent patch recovery (SPR) round trip: element → nodal least-squares fit on the old mesh → interpolation to the new nodes → node → element average. To keep the fitted quantity small, the stress is first pressure-centered (the lithostatic reference is added on the old mesh and subtracted at the new element centroids). Every commit in this series is a correction to one link of that chain, each derived from a principle rather than a symptom fix.

flowchart TB
    subgraph OLD["ON THE OLD MESH — spr_elem_to_node"]
        T([remesh triggered]) --> W["P3 : compute blend weight per element<br/>De = (η/G) / Δt_rm , w = smoothstep(log₁₀ De)"]
        W --> CTR["P2 : pressure-center at the ELEMENT CENTROID,<br/>relative to the CURRENT surface:<br/>σ += p_ref(z_eff(centroid)) , z_eff from SurfaceTopo<br/>P4 : record the added p_ref per element"]
        CTR --> SPR1["SPR patch fit: element → nodal σ_n"]
    end

    subgraph MID["MESH REGENERATION + TRANSFER"]
        SPR1 --> NEWM([new mesh built])
        NEWM --> BARY["nodal transfer: σ_n → new nodes<br/>(barycentric interpolation)"]
        NEWM --> NN["NN element transfer: σ, σ_yy, w, p_ref_old<br/>(verbatim copy where the element is unchanged)<br/>P4 : flag unchanged elements (is_changed == 0)"]
    end

    subgraph NEW["ON THE NEW MESH — spr_node_to_elem"]
        BARY --> PIN["P1 : pin the free surface at each TOP NODE:<br/>total σ_zz = 0, zero shear — in the centered variable<br/>σ_n = p_ref(z_eff(node)), and z_eff = 0 at the surface<br/>→ the pin is exact, independent of topography"]
        NN -.NN snapshot.-> UNT
        PIN --> AVG["Step C: nodal → element average<br/>(the SPR-recovered stress, all elements)"]
        AVG --> UNT{"P4 : element<br/>unchanged?"}
        UNT -- "yes (~93%)" --> KEEP["keep pre-remesh stress verbatim:<br/>restore NN snapshot, skip fallback & blend"]
        UNT -- no --> FB{"P1 : surface element left<br/>LESS compressive than<br/>before the remesh?"}
        FB -- yes --> REV["revert to NN stress<br/>(σ and σ_yy atomically)"]
        FB -- no --> BL
        REV --> BL["P3 : Deborah blend<br/>σ = w·σ_NN + (1−w)·σ_SPR"]
        BL --> DEC["P2 : de-center at the NEW ELEMENT CENTROID:<br/>σ −= p_ref(z_eff(new centroid))"]
        KEEP --> DEC2["P4 : de-center with the CARRIED reference:<br/>σ −= recorded p_ref_old (bit-exact round trip)"]
    end

    DEC --> DONE([remapped stress on the new mesh])
    DEC2 --> DONE

Loading

Principles and commits

1. Exact information beats extrapolation — pin the free surface.
The SPR patch fit is one-sided at surface nodes, its least reliable place in the mesh, while the free-surface condition (total σ_zz = 0, zero shear) is known exactly. Pin it before the node→element average, keep an NN-remapped fallback for any surface element the SPR average would leave less compressive than before the remesh (the spurious-tension direction), and stop reallocate_variables from clobbering the remapped stress the fallback needs.

• fix: preserve surface-element stress through remeshing (SPR pin + NN fallback)
• fix: OpenACC coverage and async ordering for the surface-stress remap —
  device pragmas for the new loops; the NN stress injection is ordered after the async barrier so the old tensor is never freed mid-flight.
• fix: include out-of-plane sigma_yy in the surface-stress remap fallback —
  in plane strain σ_yy is a third of the mean stress; it joins the same compressiveness test and reverts atomically with the in-plane tensor, so an element is never left in a mixed SPR / pre-remesh state.

2. Reference states must follow the current geometry.
Pressure-centering against the fixed datum z = 0 breaks down as soon as topography exists: the "small deviatoric residual" the SPR relies on grows by ρg × relief, and the surface pin lands at the wrong lithostat. SurfaceTopo rebuilds the surface from the current top boundary at every remesh (no persistent state, 2-D and 3-D) and evaluates the reference at the depth below the local surface, with a DCT-I depth-attenuation kernel (e^{−|k|d}, the harmonic-load half-space decay) so short-wavelength relief does not project undamped to depth. On a flat surface the behavior is bit-exact unchanged.

• feat: topography-corrected reference pressure in SPR remeshing (SurfaceTopo)

3. Choose the remap by the physics timescale.
The two transfer operators fail in opposite regimes: SPR's smoothing stabilizes the weak, fast-deforming region (the no-SPR variant cascades and dies) but artificially relaxes elastic stress in the cold lid, whose Maxwell time dwarfs the remesh interval; plain NN preserves memory but also preserves element-scale noise. The Deborah number is exactly the ratio that separates the regimes. Per element $e$, with viscosity $\eta_e$, shear modulus $G_e$ and the time since the previous remesh as the loading timescale,

$$\mathrm{De}_e = \frac{t_{M,e}}{\Delta t_{\mathrm{rm}}},\qquad t_{M,e} = \frac{\eta_e}{G_e},\qquad \Delta t_{\mathrm{rm}} = \max\!\big(t - t_{\mathrm{last\,remesh}},\ \Delta t\big),$$

the blend weight is a smoothstep in $\log_{10}\mathrm{De}$ between the two bounds (remesh_deborah_min = 1, remesh_deborah_max = 100 by default),

$$s_e = \mathrm{clip}\!\left( \frac{\log_{10}\mathrm{De}_e - \log_{10}\mathrm{De}_{\min}} {\log_{10}\mathrm{De}_{\max} - \log_{10}\mathrm{De}_{\min}},\ 0,\ 1\right), \qquad w_e = s_e^{2}\,(3 - 2 s_e),$$

and the remapped stress is the convex combination of the two transfers, applied in the pressure-centered variable (so the blend is exact for the lithostatic background) and after the surface fallback (so a surface element is never made less compressive again):

$$\boldsymbol{\sigma}_e \leftarrow w_e\,\boldsymbol{\sigma}_e^{\mathrm{NN}} + (1 - w_e)\,\boldsymbol{\sigma}_e^{\mathrm{SPR}}, \qquad \sigma_{yy,e} \ \text{likewise in plane strain.}$$

$w_e$ is evaluated on the old mesh — where strain rate, temperature and element markers are still mutually consistent, and before the centering, since the viscosity law reads the stress trace — and rides through the remesh with the other element fields. mesh.remesh_deborah_blend (default on; off = bit-exact pure-SPR path).

• feat: Deborah-number-weighted blend of NN and SPR stress remap at remeshing

4. Identity geometry ⇒ identity remap.
Most elements are not changed by a typical remesh, yet SPR rewrote every one of them at every event. Every other element field already crosses a remesh verbatim when the element is unchanged (inject_field); this commit gives stress the same guarantee. Unchanged elements (strict is_changed == 0 in the NN pass) keep stress and σ_yy bit-exactly: the Step C average is overwritten by the NN snapshot, the element is excluded from the surface fallback and the blend (w·s + (1−w)·s is an identity only in exact arithmetic), and Step C' subtracts the recorded reference added at centering rather than the new-topo value — the attenuation table is global, so the add/subtract pair only cancels bit-exactly with the carried number. A remesh that changes nothing is now a stress no-op.

• feat: keep pre-remesh stress verbatim in elements the remesh did not change

Benchmark: EVP rifting, five remap variants, 800 kyr

Same protocol as before (derived from examples/rifting-2d.cfg, natural remeshing, has_output_during_remeshing = yes so every event writes a before/after frame pair at the same physical time. Variants:

	stress remap across a remesh
A	nearest-neighbor only, no SPR (the two spr_* calls disabled)
B	SPR, fixed-datum reference (force_topo off at both SPR sites)
C	SPR, topography-corrected reference (this branch, blend off, untouched off)
D	+ Deborah blend (blend on, untouched off)
E	+ untouched-element carry-through (branch defaults)

What the method decides, where

(figure slot: figE_method_anatomy.png)

The implementation made visible, at one representative remesh event (Δt_rm = 19 kyr): (a) the reconstructed blend-weight field splits the domain exactly along the regime boundary — area-weighted mean w = 0.96 in the near-surface band (cold lid keeps NN memory), 0.000 in the low-viscosity band (pure SPR keeps the stabilizing smoothing), 0.30 over the whole domain. (b) the changed/unchanged mask of the same event: 93 % of elements are left geometrically unchanged — the untouched carry-through keeps all of them bit-exact — and the changed ones are the rift interior plus the bottom-boundary strip, which the floor re-discretization rebuilds at every event. (c) that 91–97 % unchanged fraction holds at every event of every variant: a typical remesh touches less than a tenth of the mesh, which is why principle 4 removes most of the remap footprint.

Setup and dynamical stability

(figure slot: figE_setup_stability.png)

The benchmark develops rift topography and a deep low-viscosity band under natural remeshing. A no longer crashes — the branch's remeshing-robustness fixes turn the old 25 kyr triangulation failed death into a graceful permanent cascade (381 remesh events in 87 kyr, low-viscosity band 280 → 3300+ elements, stopped deliberately) — but the physics verdict against NN-only is unchanged. All four SPR variants ran quietly to the 800 kyr time limit with no terminal cascade.

Per-event remap error

RMS over the common window t ≤ 779 kyr (MPa):

variant	ΔP domain	Δτ_II domain	ΔP surf (<15 km)	Δτ_II surf	ΔP low-η band	Δτ_II low-η
A (died 87 kyr)	5.95	0.93	~0	~0	7.47	1.16
B	11.39	11.02	7.38	5.00	1.94	0.11
C	9.94	9.73	5.15	4.24	1.70	0.11
D	9.16	9.10	1.74	1.62	1.87	0.10
E	1.73	1.13	1.21	0.85	1.97	0.42

(figure slot: figE_error_series.png)

Each principle removes the part of the error it targets: topography correction (B→C) −30 % near-surface ΔP; Deborah blend (C→D) −66 % near-surface (the cold lid keeps NN); untouched carry-through (D→E) cuts the whole-domain error ~5× by eliminating the deep error band entirely. In the low-viscosity band — genuinely remeshed, De ≪ 1 — all SPR variants including E are statistically identical in ΔP, so the stabilizing smoothing is preserved exactly where it is needed; E's weak-band Δτ_II is elevated (0.42 vs 0.11 MPa) but bounded — it oscillates and decays, band size stays normal, no A-style runaway.

Where the error lives: one event, and all events

(figure slot: figE_remap_maps.png — element-resolved ΔP at the ~500 kyr event of each variant, after − before at the same physical time)

(figure slot: figE_remap_avg_maps.png — RMS ΔP over ALL remesh events on a common grid)

B and C pay a saturated two-band error at the lithosphere base across the full domain width at every event, plus near-surface fabric; D cleans the lid but keeps the deep bands; E is near-zero everywhere except the actively-remeshed rift interior. A's error is a ±40 MPa element-scale checkerboard confined to the deep weak band — while its surface is usually re-triangulated identically, which is why a surface metric alone underestimates the NN-only failure.

Bias, not noise

(figure slot: figE_remap_mean_maps.png — SIGNED mean ΔP over all events)

Keeping the sign separates the two failure modes. A's signed mean is ≈ 0 despite its saturated RMS: pure zero-mean noise, amplified by the EVP feedback. B/C/D's band dipole survives averaging at nearly full amplitude — every remesh pumps pressure the same direction across the brittle–ductile transition (SPR smoothing the same curvature every time). A same-direction bias accumulates where noise would cancel; this is the mechanism behind SPR's cold-lid stress relaxation, and E removes it.

Where the peak velocity lives

(figure slot: figE_vmax_locations.png — peak-|v| node of every frame over the RMS error background)

E's cost accounting: ~50 % more remesh events (33 vs 22) and recurring post-remesh velocity transients (~15× v_bc, decaying, never cascading, quietening after 600 kyr). Locating the per-frame peak velocity shows these transients live at the bottom wall/floor corner — not at the rift or distributed through the interior; baseline frames peak in the rift upwelling column exactly like B/C/D. The anatomy figure explains why: the floor strip is re-discretized at every remesh (panel b), so its elements always take remapped values against untouched neighbours above, under the Winkler support. The transient is a localized boundary artifact with an evident refinement path (exclude boundary-adjacent elements from the carry-through), not an inherent instability of the method.

Rift architecture

(figure slot: figE_topography_localization.png — topography evolution + fault patterns)

Same graben, conjugate fault sets and ~25 km localization depth in B/C/D/E; shoulder heights at 800 kyr agree to ±2 % with no fidelity ordering (single realization per variant, so shoulder height is realization scatter, not a remap signal).

(figure slot: figE_stressII_comparison.png — deviatoric stress II in the same rift window at 800 kyr)

The stress state each remap leaves behind, as a field: B and C carry the smoothest τ_II — the lid's stress a soft blur, the fault stress-shadows smeared away — the cumulative fingerprint of SPR smoothing at every event. D is crisper; E retains the most structure: distinct low-stress fault shadows in the lid on both flanks, element-scale texture, and the sharpest high-stress core and axial stress shadow, with the same large-scale architecture as the others — fidelity, not divergence.

Testing

• 2-D and 3-D build cleanly; tests/functional/2d-ep-irregular.cfg and 3d-evp-regular.cfg pass with remeshing exercised.
• Both new flags are bit-exact off, and with the flags on the output frames before the first remesh are field-identical to the flags-off run, with only post-remesh frames differing — each option engages exactly at remeshing and nowhere else.
• Restart is bit-deterministic across post-restart remeshes with the blend on (exercising the checkpointed last_remesh_time) and off.
• All four surviving benchmark variants ran to the 800 kyr time limit with no terminal cascade — the deaths reported in the earlier revision of this benchmark (541/758 kyr) were removed by the branch's robustness fixes, independent of remap choice.

Earlier four-variant benchmark (superseded by the five-variant rerun above; kept for the record)

Variant labels in these figures follow the old scheme (A fixed-datum,
B topo-corrected, C NN-only, D blend), which maps to the new B, C, A, D.

🤖 Generated with Claude Code

[Copilot]

Pull request overview

This PR consolidates a set of general-purpose fixes aimed at improving remeshing stability and accuracy (notably stress remap behavior near the free surface and under evolving topography), strengthening interpolation robustness, and polishing build/benchmark tooling.

Changes:

• Refresh dt unconditionally after remeshing to prevent CFL instability on refined meshes.
• Improve remeshing stress handling via NN stress carry-through + SPR surface pinning, and introduce a topography-aware reference pressure (SurfaceTopo) for SPR centering/restoration.
• Replace the interpolation fallback with a true (layer-capped) BFS and add tooling improvements (macOS OpenMP rpath, benchmark isolation via INDIR).

Reviewed changes

Copilot reviewed 12 out of 12 changed files in this pull request and generated 6 comments.

Show a summary per file

File	Description
remeshing.cxx	Recomputes dt unconditionally after remesh; integrates remap flow relying on preserved stress data.
parameters.hpp	Comments out unused pseudo-transient dt_PT state.
nn-interpolation.cxx	NN-remaps var.stress so SPR can fall back for surface elements.
Makefile	Adjusts macOS OpenMP rpath handling when OPENMP_ROOT_DIR is relative.
geometry.hpp	Declares SurfaceTopo and removes/limits PT compute_dt_PT declarations.
geometry.cxx	Implements SurfaceTopo and uses it in SPR centering/restoration plus free-surface pin + surface fallback.
fields.hpp	Comments out unused update_velocity_PT declaration.
fields.cxx	Prevents var.stress from being reallocated (to preserve NN-remapped stress for SPR fallback); comments out PT velocity update.
dynearthsol.cxx	Comments out unused PT dt_PT initialization.
brc-interpolation.cxx	Replaces one-level neighbor expansion with a capped BFS for enclosing-element search.
benchmarks-cores/Makefile	Adds INDIR option to run benchmarks in per-case subdirectories.
bc.cxx	Clamps open-sidewall lithostatic traction at zero to avoid unphysical suction above datum.

[chatgpt-codex-connector]

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: b8b810918c

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

• Open a pull request for review
• Mark a draft as ready
• Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

[chatgpt-codex-connector]

💡 Codex Review

Here are some automated review suggestions for this pull request.

Reviewed commit: 934863d6d7

ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

• Open a pull request for review
• Mark a draft as ready
• Comment "@codex review".

If Codex has suggestions, it will comment; otherwise it will react with 👍.

Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".

[echoi]

This PR is a mixture of one major improvement (in my opinion) in remeshing and other bugfixes. I really appreciate the detailed benchmark for four stress remap variants, which demonstrates the clear merit of the Deborah number-based approach. It could suppress explosive perturbations right after remeshing, a long-time pathology the FLAC family has suffered, by weighting the NN and the SPR remapping: large maxwell time -> NN, small maxwell time -> SPR.

Before approving this PR, I just want to understand some details better.

mesh.remesh_deborah_blend, now on by default (= no restores the bit-identical pure-SPR path); per element σ = w·σ_NN + (1−w)·σ_SPR with w a smoothstep in the local Deborah number De = (η/G)/Δt_remesh, so the cold lid keeps NN (elastic stress memory) and the weak zone keeps SPR (smoothing). Both inputs pressure-centered, so the blend is exact for the lithostat.

1. what does "pressure-centered" mean here?
2. How was this expression for weighting factor derived?

            const double t_maxwell = var.mat->visc(e) / var.mat->shearm(e);
            double t = (std::log10(t_maxwell / dt_remesh) - lde0) / (lde1 - lde0);
            t = std::min(std::max(t, 0.0), 1.0);
            (*var.spr_blend_weight)[e] = t * t * (3.0 - 2.0 * t);  // smoothstep

3. I can see the De's min/max can be set by a user through remesh_deborah_min and remesh_deborah_max. I don't have a very good sense of when to touch them and what values to use.

[chaseshyu]

1. what does "pressure-centered" mean here?

It means stress tenser is subtracted by reference pressure. Magnitude is similar to deviatoric stress but it's not.

2. How was this expression for weighting factor derived?

            const double t_maxwell = var.mat->visc(e) / var.mat->shearm(e);
            double t = (std::log10(t_maxwell / dt_remesh) - lde0) / (lde1 - lde0);
            t = std::min(std::max(t, 0.0), 1.0);
            (*var.spr_blend_weight)[e] = t * t * (3.0 - 2.0 * t);  // smoothstep

The equation here is for gentle transition of weighting.

3. I can see the De's min/max can be set by a user through remesh_deborah_min and remesh_deborah_max. I don't have a very good sense of when to touch them and what values to use.

De < remesh_deborah_min -> pure SPR
De > remesh_deborah_max -> pure NN
range between min and max -> the sharpness of transition.

[echoi]

@chaseshyu I'm sorry I still don't understand this statement: "Both inputs pressure-centered, so the blend is exact for the lithostat." My question is more like, why do pressure-centered (i.e., removed) quantities need to be considered in the first place, and what does "exact for the lithostat" means?

The curve looks nice but the functional form x^2(3-2x) seems arbitrary without further information. That's what I was asking about.

De < remesh_deborah_min -> pure SPR
De > remesh_deborah_max -> pure NN
range between min and max -> the sharpness of transition.

I got it. Still, when should a user consider a different range instead of the default, 2? What about the same range for min max, e.g., 10^3 and 10^5? Would it be physically meaningful? If not, why is it even allowed to happen?

[chaseshyu]

@chaseshyu I'm sorry I still don't understand this statement: "Both inputs pressure-centered, so the blend is exact for the lithostat." My question is more like, why do pressure-centered (i.e., removed) quantities need to be considered in the first place, and what does "exact for the lithostat" means?

Because the background pressure is huge, we use pressure-centered stress to calculate NN and SPR to reduce precision and interpolation errors. In the code context, the analytical reference pressure will be added back after blending the pressure-centered stress. "Exact for the lithostat" just means this process preserves the background pressure field perfectly, without introducing numerical artifacts into the following main loop kernel.

The curve looks nice but the functional form x^2(3-2x) seems arbitrary without further information. That's what I was asking about.

Yes, it is an arbitrary function to switch between NN and SPR.

De < remesh_deborah_min -> pure SPR
De > remesh_deborah_max -> pure NN
range between min and max -> the sharpness of transition.

I got it. Still, when should a user consider a different range instead of the default, 2? What about the same range for min max, e.g., 10^3 and 10^5? Would it be physically meaningful? If not, why is it even allowed to happen?

In Deborah number De = (η/G)/Δt_remesh, η/G is Maxwell relaxation time, which is the timescale over which viscous flow relaxes elastic stress. Δt_remesh is the time elapsed since the last remesh. if De >> 1, we should keep elastic stress. if De << 1, we should relax it. To me, it acts more like a switch: if the velocity of low-viscosity materials becomes chaotic after remeshing, I will increase the min bound. If the near-surface stress diffuses away after remeshing, I will decrease the max bound. I don't know if this really answers your question... XD

[echoi]

we use pressure-centered stress to calculate NN and SPR

Sorry, this is news to me... :(

In Deborah number De = (η/G)/Δt_remesh, η/G is Maxwell relaxation time, which is the timescale over which viscous flow relaxes elastic stress. Δt_remesh is the time elapsed since the last remesh. if De >> 1, we should keep elastic stress. if De << 1, we should relax it. To me, it acts more like a switch: if the velocity of low-viscosity materials becomes chaotic after remeshing, I will increase the min bound. If the near-surface stress diffuses away after remeshing, I will decrease the max bound. I don't know if this really answers your question... XD

Sounds like De =1 is special because it's used as a reference in your reasoning. It seems to me to mean that the minimum need not be a parameter a user can change. I think I need more time to fully understand the whole picture, but your approach makes sense. I just hope your "work plan" is documented somehow.

[chaseshyu]

Sounds like De =1 is special because it's used as a reference in your reasoning. It seems to me to mean that the minimum need not be a parameter a user can change. I think I need more time to fully understand the whole picture, but your approach makes sense. I just hope your "work plan" is documented somehow.

Sorry for the messy commits. I think it’s a great idea to document the work plan. Does it sound like a concept/background of PRs and commits? Maybe I can document it as an issue before the PR or commit it in the doc folder in the PR. That way, reviewers can get a better idea of the whole picture. What do you think? Any suggestions?

[chatgpt-codex-connector]

You have reached your Codex usage limits for code reviews. You can see your limits in the Codex usage dashboard.

[chaseshyu]

Hi @echoi, this PR#74 has been updated to focus on stress remeshing and is ready for review.