Restore constrained neo-Hookean reference-map hyperelasticity

[sbryngelson]

Description

This is the add-back half of the stacked workflow proposed in #1674. It is intentionally based on the removal branch in #1675, so this PR shows the hyperelastic capability being added back from a clean baseline. After #1675 merges, this PR should be retargeted to master.

The restored implementation is deliberately narrower than the prior interface: neo-Hookean only, 3D Cartesian only, finite-difference gradients enabled, and HLLC only. It preserves the project-specific Barbosa/Rodriguez-family closure already used by MFC, tau = (G/J) dev(b), rather than silently replacing it with the different J^(-5/3) closure in Kamrin et al.

The repair:

• registers hyper_model in the generated parameter/namelist system and accepts only hyper_model = 1;
• rejects unsupported dimensionality, coordinates, gradient configuration, model equations, and Riemann solvers during case validation;
• removes the undocumented spherical pre_stress initializer;
• corrects the tensor contraction to compute b = F F^T from the stored F^T layout;
• removes the duplicate pseudo-Mooney-Rivlin routine;
• initializes and consistently converts the derived stress and elastic-energy slots;
• gives every hyperelastic state slot a defined HLLC flux and uses the HLLC mass flux for the conservative reference-map transport;
• fixes direction indexing in the elastic energy flux and treats the left/right elastic energies independently;
• adds an always-run, deforming 3D shear canary with nonzero off-diagonal elastic stresses.

The audit, concrete equations, provenance discussion, references, and remaining architectural caveats are in #1674.

Closes #1674.

Type of change (delete unused ones)

• Bug fix
• Refactor
• Documentation

Testing

• ./mfc.sh precheck
• ./mfc.sh build -t pre_process simulation post_process --no-mpi --debug -j 4
• ./mfc.sh test --from C06DAA3C --to C06DAA3C --no-mpi --debug -j 1
• ./mfc.sh test --from 32A3A936 --to 32A3A936 --no-mpi --debug -j 1

The new C06DAA3C golden case is a 3D shear problem. It evolves the reference map away from the identity and produces nonzero off-diagonal elastic stresses, so it exercises the corrected b = F F^T path and the HLLC elastic/reference-map fluxes.

Checklist

• I added or updated tests for new behavior
• I updated documentation if user-facing behavior changed

See the developer guide for full coding standards.

GPU changes (expand if you modified src/simulation/)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

GPU validation remains to be done.

AI assistance disclosure

This audit and patch were AI-assisted. The technical claims were checked against the repository source, CPU debug builds/runs, the new regression output, and the references collected in #1674. Human review is especially requested for the intended constitutive closure and the transitional treatment of derived stress/energy slots.

Review requested from @mrodrig6 and @ChrisZYJ.

[sbryngelson]

Follow-up correctness audit: this should remain a draft

This is an AI-assisted follow-up audit of the add-back at 4aa4488. The constitutive core is substantially improved, but I do not think the current PR is yet a correct end-to-end restoration.

Parts that appear correct

• With the storage layout used in m_hyperelastic.fpp, the corrected contraction computes b = F F^T, not F^T F.
• Retaining tau = (G/J) dev(b) is consistent with the Barbosa/Rodriguez-family formulation already used by MFC; it should not be silently replaced by the different isochoric closure in Kamrin et al.
• Conservative reference-map transport as rho xi, using the mass flux, is conceptually sound.
• Removing the unsupported pseudo-Mooney-Rivlin path and registering/constraining hyper_model are good changes.

1. Reference-map boundary conditions are currently incorrect

The new C06DAA3C test uses BC_GHOST_EXTRAP = -3 on all six faces. The generic routine copies the nearest interior value into every ghost cell. For the undeformed map xi_x = x, second-order centered differentiation at the first cell therefore gives

(d xi_x/dx)_0 = (xi_1 - xi_-1)/(2 dx) = 1/2

rather than one. On an x face this creates

F = diag(2,1,1),  J = 2

and a nonzero elastic stress in an initially undeformed body. Thus the golden test starts with numerical pre-strain near its boundaries.

This should not be addressed by rejecting most MFC boundary conditions. The completed capability needs a dedicated reference-map boundary layer covering every applicable MFC BC:

• Periodic (-1): the reference map is affine-periodic. Copy the tangential components normally, but offset the normal component by the reference-domain length. For example, at an x wrap,

  xi_x(-j,k,l)  = xi_x(m-(j-1),k,l) - Lx_ref
  xi_x(m+j,k,l) = xi_x(j-1,k,l)     + Lx_ref

  The same offset must be applied at global periodic MPI wraps, not at ordinary internal rank boundaries.

• Reflective/symmetry (-2) and characteristic/slip walls (-5, -15): reflect the normal reference coordinate about the reference-coordinate wall location,

  xi_n(ghost) = 2*Xi_wall - xi_n(mirror)

  and mirror the tangential reference-map components evenly. The current unconditional -xi_n(mirror) is only valid when Xi_wall = 0; it is wrong at an upper wall in a general domain.

• Ghost/Riemann extrapolation and outflow families (-3, -4, -8, -9, -10, -12): use affine-preserving extrapolation for every xi component. For second order,

  xi_i(-j)  = (j+1)*xi_i(0) - j*xi_i(1)
  xi_i(m+j) = (j+1)*xi_i(m) - j*xi_i(m-1)

  This preserves the identity map and any other affine deformation exactly. If fourth-order reference-map gradients are supported, the ghost construction must be raised consistently or one-sided fourth-order derivative stencils must be used.

• No-slip walls (-16): distinguish a fixed/clamped solid wall from a slip/symmetry wall. The normal reference component reflects about Xi_wall. Tangential components require the prescribed wall reference map for a clamped or moving wall; simply copying them imposes a zero-normal-gradient condition and is not a general no-slip solid boundary.

• Subsonic/supersonic inflow and Dirichlet (-6, -7, -11, -17): incoming material needs prescribed reference-map data. Guessing xi=x at every time is not sufficient for uniform advection (xi=x-u t) and cannot represent a prescribed incoming deformation. Extend the BC input/buffer interface with reference-map data, preferably a boundary reference origin plus F_in (or equivalently grad_xi_in) and its time dependence. Dirichlet bc_buffers must populate these slots explicitly rather than relying on whatever values happen to occupy the expanded sys_size.

• Null (-13): retain its existing no-op meaning, but do not advertise it as a usable physical hyperelastic boundary if the finite-difference stencil would read unfilled ghosts.

• Axis (-14): this is not applicable to the PR's stated 3D Cartesian-only scope. Literal support for the axis BC requires cylindrical/axisymmetric reference-map gradients, tensor components, geometric sources, and constitutive terms; it is a separate model extension, not merely another ghost-fill case. It should remain explicitly tracked rather than silently accepted.

A useful implementation pattern is to run the existing fluid BC routine and then call one centralized s_apply_reference_map_bc(...) that overwrites only eqn_idx%xi according to the BC code. This avoids scattering slightly different xi rules through every generic boundary primitive.

2. The transverse HLLC elastic-energy correction is structurally wrong

The current correction contains, for transverse components,

(s_S - u_t) * tau_nt / (s_K - u_t)

Here s_S and s_K are wave-normal speeds while u_t is a transverse velocity. An otherwise admissible state can have u_t = s_K, producing an artificial singularity. Also, s_S is not the transverse star velocity. Replacing only the denominator with s_K-u_n removes one singularity but does not provide the missing transverse star state.

The safest first restoration is to use the two-wave HLL flux for the complete conservative hyperelastic state:

F_HLL = (s_R F_L - s_L F_R + s_L s_R (U_R-U_L))/(s_R-s_L)

with physical fluxes

F_mom = rho u_n u + p e_n - tau e_n
F_E   = (E+p)u_n - u . (tau e_n)
F_rhoxi = rho u_n xi

MFC's HLL implementation already has partial hyperelastic wave-speed logic, but its stress-aware momentum and energy branches and all rho xi state fluxes need to be completed. This is more diffusive than HLLC but avoids inventing shear-star relations.

If HLLC is retained, the elastic branch should instead be rewritten from the Rankine-Hugoniot relations. Let P_K = p_K - tau_nn,K. A consistent reduced three-wave assumption is u_t* = u_t and tau_nt* = tau_nt, giving

rho_K* = rho_K (S_K-u_n,K)/(S_K-S*)
P_K*   = P_K + rho_K (S_K-u_n,K)(S*-u_n,K)
E_K*   = ((S_K-u_n,K)E_K - P_K u_n,K + P_K* S*)/(S_K-S*)
F_E,K* = (E_K*+P_K*)S* - sum_t u_t,K tau_nt,K

This has no S_K-u_t denominator and satisfies the assumed three-wave jump relations. It still does not resolve the physical shear waves; resolving those requires a larger elastic HLLC/HLLD-type wave model.

Until that work is complete, validation should at least require direct wave speeds (wave_speeds = 1). The pressure-based branch currently omits the elastic changes used by the direct branch.

3. Invalid maps must not silently lose elasticity

m_hyperelastic.fpp zeros stress and elastic energy before checking det(grad xi), and leaves them zero if the determinant is invalid. That silently turns a failed solid cell into a locally fluid-like cell.

Use a device-safe invalid-map counter/flag, record the minimum determinant, reduce it across MPI ranks, and terminate with a clear diagnostic after the kernel. If regularization is desired instead, it needs to be an explicit documented model with diagnostics, not implicit zeroing.

4. The current golden test is an execution canary, not a correctness test

Because its ghost extrapolation already strains the identity map, C06DAA3C cannot validate the mechanics. Add analytic assertions for:

• identity, including every supported boundary type: F=I, J=1, tau=0, elastic energy zero;
• rigid rotation: b=I, zero stress;
• simple shear with every component checked analytically;
• rotated stretch, which distinguishes F F^T from F^T F;
• affine-periodic maps across serial and MPI periodic wraps;
• uniform inflow/advection with the time-dependent reference map;
• x/y/z rotations of the same Riemann state;
• the formerly singular case u_t = S_L;
• invalid determinant producing the intended failure;
• model-equation and CPU/GPU parity for every configuration claimed as supported.

Until these points are addressed, this PR should remain draft and should not close #1674.

References

• Carcana Barbosa, Rodriguez, Yang et al., CTR Summer Program 2024: https://web.stanford.edu/group/ctr/ctrsp24/iv07_CARCANA.pdf
• MFC equation documentation: https://mflowcode.github.io/documentation/equations.html
• Rodriguez et al., JCP 2019: https://doi.org/10.1016/j.jcp.2018.10.035

This comment and audit were AI-assisted; the findings above are based on direct inspection of the current PR source and the cited formulations.