[ADDITION][FIXES] DYREL VS; Sol Benchmarks and Bug fixes

[aelligp]

What is the purpose of this PR?

• Bug fix
• New feature
• Documentation update
• Solver and benchmark improvements

Summary

This PR adds a 2D variational Stokes implementation of the self-tuned dynamic-relaxation (DYREL) solver, improves the original variational APT solver, adds free-surface stabilization (FSSA), and introduces SolCx/SolKz/SolVi validation benchmarks. It also fixes correctness, conditioning, masking, boundary-condition, non-uniform-grid, and 3D stress-rotation issues uncovered during validation.

Following Larionov, Batty, and Bridson (2017), pressure, stress, and staggered velocity samples are weighted by their rock control-volume fractions, collected in diagonal matrices W_p, W_τ, and W_u. Schematically,

$$R_u = D^T W_\tau \tau - G W_p p - f = 0, \qquad R_p = G^T u = 0.$$

Variables whose constraints reference a zero-weight velocity sample are removed from the active system, eliminating the variational null spaces described in Section 5.3 of the paper.

Variational DYREL solver

• Adds _solve_DYREL! dispatch on ϕ::RockRatio, mirroring the existing variational APT dispatch. DYREL remains 2D because its Gershgorin estimator is currently 2D-only.
• Adds constructors for the bulk response, penalty coefficients, preconditioner, and spectral estimates of the reduced variational space.
• Adds masked divergence, strain-rate, pressure-residual, stress, momentum-residual, and damped velocity-update kernels.
• Adds a ϕ-aware Gershgorin stencil using the same weighted viscosity, penalty, pressure, stress, and FSSA contributions as the matrix-free operator.
• Ensures fully masked non-finite values contribute exact zero instead of Inf*0 = NaN.

Reduced-space consistency

The active pressure rule is

$$p_{ij}\text{ active} \iff W^p_{ij}>0 \quad\text{and every velocity face in }(G^T u)_{ij}\text{ has positive weight}.$$

The same mask is used by the pressure residual, convergence norm, field spans, and active degree-of-freedom count. The empirical lower-eigenvalue estimate used by the DR auto-tuner is likewise restricted to active, boundary-trimmed velocity rows:

$$\widehat\lambda_{\min} =\frac{\left|\Delta u^T(R^{n+1}-R^n)\right|}{\Delta u^T\Delta u}, \qquad \Delta u,R\in V_{\mathrm{active}}.$$

Previously, eliminated cut-cell and boundary entries could contaminate this quotient even though they were absent from the convergence norm.

Powell-Hestenes penalty weighting

For an active pressure sample, W_p G^T u=0 is equivalent to the retained constraint G^T u=0. The numerical Powell-Hestenes coefficient is therefore not itself volume-weighted:

$$\theta_c=\gamma R_p+\Delta P_\psi, \qquad R_u^{(\gamma)}=-G W_p\theta_c.$$

The variational gradient supplies the single physical volume weight. The previous construction weighted γ_eff by ϕ and then applied the weighted gradient, producing approximately ϕ² scaling near the interface:

$$G W_p\Gamma W_pG^T \quad\longrightarrow\quad G W_p\Gamma G^T.$$

The physical bulk viscosity remains η_b = ϕ K_b Δt; only the algorithmic PH penalty is unweighted before the weighted momentum gradient.

Automatic DR tuning and Gershgorin validation

Pseudo-time steps and damping coefficients are selected from

$$\Delta\tau_V=\frac{2\,\mathrm{CFL}_V}{\sqrt{\lambda_{\max}}}, \qquad \beta_V=\frac{2\Delta\tau_V}{2+c_V\Delta\tau_V}, \qquad \alpha_V=\frac{2-c_V\Delta\tau_V}{2+c_V\Delta\tau_V},$$

with λ_max supplied by the D⁻¹A Gershgorin row bound and c_V derived from the reduced-space estimate of λ_min.

A CPU regression assembles the actual matrix-free velocity Jacobian by applying production kernels to unit velocity vectors. For full and cut domains it verifies symmetry to roundoff, full reduced rank, positive definiteness, and that every exact preconditioned absolute row sum is bounded by the stored Gershgorin estimate.

For a controlled horizontal cut, the Jacobi/Gershgorin-preconditioned spectra are:

Grid	Full κ	Cut κ
4²	39.2	44.2
6²	81.7	89.9
8²	137.7	151.5
10²	207.4	229.2

The cut interface adds about 10% in this case; the dominant O(N²) trend is the expected limitation of diagonal preconditioning.

Free-surface stabilization

The vertical residual contains the implicit Kaus et al. (2010) correction

$$R_y^{\mathrm{FSSA}}=V_y\,\partial_y(\rho g\phi)\,\theta\,\Delta t, \qquad \theta=1.$$

Accordingly:

• the DYREL Gershgorin estimator includes its magnitude in the V_y diagonal and row bound;
• ordinary and variational DYREL rebuild D_y, λ_max,y, Δτ_y, α_y, and β_y before the first DR window and during periodic tuning;
• the variational estimator uses the same ϕ-weighted buoyancy derivative as the residual;
• γ_eff receives an FSSA-dependent floor so the penalty diagonal cannot become negligible relative to FSSA.

Original variational APT solver

The original non-DR solver now follows the same reduced-space convergence rules. Its residuals are active-space RMS values,

$$E_{u_d}=\frac{\|M_{u_d}R_{u_d}\|_2}{\sqrt{N_{u_d}}}, \qquad E_p=\frac{\|M_pR_p\|_2}{\sqrt{N_p}},$$

where M_p is the isvalid_c(ϕ) mask, M_{u_d} selects retained staggered velocity rows, and each N is the corresponding global active degree-of-freedom count. This replaces full-grid normalization, which made a fixed tolerance depend on rock fraction and resolution. It also fixes the 3D continuity metric, which divided the L2 norm by N_p instead of sqrt(N_p).

For free_surface=true, the APT vertical update now includes the FSSA diagonal in its local pseudo-time denominator:

$$V_y^{k+1}=V_y^k+\frac{\eta\,d\tau}{\eta_{\tau,y}}R_y^k, \qquad \eta_{\tau,y}=\eta_\tau+\eta\,d\tau\left|\Delta t\,\partial_y(\rho g\phi)\right|.$$

Because η dτ has dimensions of length squared, the added term has viscosity units. It is identically zero when FSSA is disabled and bounds the pseudo-time amplification caused by a sharp free-surface density jump.

Matched sticky-air comparison

The deterministic plume harness can run ordinary sticky air or reduced variational Stokes with identical geometry, timestep, penalty, tolerances, initial fields, and residual normalization. Without FSSA:

Grid	Sticky-air iterations	Variational iterations
20²	2,600	1,900
40²	3,300	2,350
80²	4,750	3,950

Both formulations converge at every resolution. A permanent 40² regression requires both solves to converge and prevents the reduced solve from becoming more than 25% slower than sticky air.

Other changes

• Adds SolCx, SolKz, and SolVi DYREL miniapps and plain-versus-variational (ϕ ≡ 1) equivalence tests.
• Uses local-viscosity penalty scaling so γ_eff/η = O(γfact) across heterogeneous materials.
• Aligns residual normalization with the APT value-span convention and returns err, iter, and converged diagnostics.
• Prevents the inner DYREL pressure residual from stalling the velocity solve target.
• Warns when inner or outer iteration limits are exhausted without convergence.
• Adds allocation-free MPI reductions for weighted and masked residual norms.
• Fixes a no-slip bottom ghost update, all-air free-surface phase-ratio handling, phase-ratio clamping, degenerate Gershgorin diagonals/moduli, non-uniform-grid topography lookup, and 3D stress-rotation interpolation/vorticity handling.

Validation

• test_dyrel_kernels: 69/69
• core DYREL tests: 67/67
• RockRatio tests: 78/78
• variational free-surface plume, including matched sticky air: 12/12
• Sol benchmark plain/ϕ≡1 equivalence: 18/18
• full non-MPI CPU suite: 1420/1422; the two errors were test-environment failures because ExactFieldSolutions was unavailable to SolCx/SolKz visualization tests
• all MPI tests passed, including 2D and 3D shear-band cases

The latest focused APT rerun is currently blocked before test execution by an external GeoParams version mismatch: the branch references ThreePhase_Density, which is absent from the locally resolved GeoParams. Static parsing and git diff --check pass. CUDA/AMDGPU execution was not available on the development machine, so GPU validation remains for CI or a GPU-enabled runner.

Checklist

• New behavior has focused regression coverage
• Affected miniapps and benchmarks are updated
• Public API remains backward compatible
• Reduced operator, damping estimates, and convergence masks are mutually consistent
• CPU and MPI validation completed before the external dependency mismatch
• CUDA/AMDGPU validation completed on a GPU runner