Hyperelastic reference-map implementation is currently non-operational and internally inconsistent

[sbryngelson]

@mrodrig6 @ChrisZYJ — this audit found several independent problems in the current hyperelastic path. The proposed workflow is to remove the exposed, non-operational capability in one PR and add it back through a stacked restoration PR after the implementation and validation below are complete.

Summary

At audited commit 8821173c20b645f1869b0914f4a64673d6021455, the hyperelastic implementation should not be treated as a working user-facing capability:

1. hyper_model is initialized to dflt_int (-100) but is absent from the parameter registry and generated namelist. Therefore neither constitutive routine is callable from a normal MFC case. The routine nevertheless consumes the derived elastic-energy field immediately afterward.
2. m_hyperelastic.fpp is structurally 3D-only but validation permits 2D. The kernel unconditionally uses all nine gradient components and the y/z finite-difference arrays; in 2D, tensor_size == 5 while entries through 9 are accessed.
3. The tensor contraction computes the wrong Cauchy–Green tensor. With the current storage order, the inverted tensor is F^T; the existing row-dot-row contraction produces F^T F (right Cauchy–Green), while the spatial Cauchy stress requires b = F F^T (left Cauchy–Green).
4. Stress and the derived elastic-energy coefficient are included in sys_size and therefore read by the generic RHS flux divergence, but HLLC does not assign fluxes for those derived entries. HLL does not assign the reference-map flux either.
5. s_Mooney_Rivlin_cauchy_solver is algebraically identical to the neo-Hookean routine and has no second-invariant contribution or independent material coefficient.
6. The HLLC elastic-energy flux uses vel_[LR](dir_idx(i)) in one part of the expression but vel_[LR](i) in its star-state correction, so transverse-direction fluxes use inconsistent velocity components.
7. The HLLC energy initialization includes both left and right hyperelastic energies only when both G_L and G_R are nonzero. Each side must be gated independently.
8. If det(grad xi) <= verysmall, the constitutive projection silently retains stale stress/energy instead of reporting loss of map admissibility.
9. There are no hyperelastic regression cases or golden-test metadata in the repository.

Intended provenance and constitutive decision

The reference-map approach itself is appropriate. The project-specific formulation is documented most directly by Carcana Barbosa, Rodriguez Jr., and Yang (2024):

• F = (grad xi)^(-1)
• b = F F^T
• sigma_dev = (G/J) dev(b)
• conservative reference-map transport d_t(rho xi) + div(rho u tensor xi) = 0

The same closure is reported in the MFC 5.0 paper. Therefore the immediate tensor fix should implement b = F F^T, not change the overall modeling family.

There is a separate constitutive question that should be decided explicitly with the model authors. Kamrin, Rycroft, and Nave's compressible neo-Hookean example uses the isochoric invariant J^(-2/3) tr(b) and obtains a deviatoric factor proportional to G J^(-5/3) dev(b), whereas the MFC/Barbosa formulation uses G J^(-1) dev(b). The restoration PR should not silently change this exponent; it should preserve the published MFC closure unless the authors confirm a different intended strain-energy split.

Concrete correction for the tensor contraction

The current derivative layout is

! [xi_x,x xi_y,x xi_z,x;
!  xi_x,y xi_y,y xi_z,y;
!  xi_x,z xi_y,z xi_z,z] = (grad xi)^T

After inversion, tensora stores F^T. Thus b = F F^T = transpose(tensora)*tensora:

tensorb(1) = tensora(1)**2 + tensora(4)**2 + tensora(7)**2
tensorb(5) = tensora(2)**2 + tensora(5)**2 + tensora(8)**2
tensorb(9) = tensora(3)**2 + tensora(6)**2 + tensora(9)**2
tensorb(2) = tensora(1)*tensora(2) + tensora(4)*tensora(5) + tensora(7)*tensora(8)
tensorb(3) = tensora(1)*tensora(3) + tensora(4)*tensora(6) + tensora(7)*tensora(9)
tensorb(6) = tensora(2)*tensora(3) + tensora(5)*tensora(6) + tensora(8)*tensora(9)

A rotated uniaxial stretch is the essential regression test: it distinguishes F F^T from F^T F, unlike diagonal stretches and simple identity checks.

Recommended restoration scope

1. Register hyper_model; default to neo-Hookean and reject unsupported values. Do not expose Mooney–Rivlin until it has a real constitutive law and parameters.
2. Initially constrain the restored implementation to 3D Cartesian HLLC, matching the actual kernel. Add 2D/axisymmetric support separately.
3. Correct b = F F^T.
4. Keep only rho xi as evolved hyperelastic state. Store derived stress and elastic energy as auxiliary/projected fields rather than conservation equations. If a minimal transitional implementation retains those slots, initialize every flux and re-project after every Runge–Kutta stage.
5. Compute the HLLC reference-map flux as the upwind reference map multiplied by the actual mixture mass flux; add an explicit HLL passive-scalar flux only if HLL is supported.
6. Fix the rotated-direction velocity indices and gate left/right energies independently.
7. Fail clearly or apply a documented admissibility policy when det(grad xi) <= 0.
8. Add CPU bounds-checking tests and analytic tests for identity, rigid rotation, simple shear, rotated uniaxial stretch, uniform reference-map advection, and determinant failure. Add GPU parity before declaring the restored capability ready.

References

• Kamrin, Rycroft, and Nave, Reference map technique for finite-strain elasticity and fluid–solid interaction, JMPS 60 (2012), DOI: https://doi.org/10.1016/j.jmps.2012.06.003, open PDF: https://www.math.mcgill.ca/jcnave/papers/elasticity.pdf
• Carcana Barbosa, Rodriguez Jr., and Yang, Numerical simulations of inertial microbubble collapse at a gel–water interface with finite elasticity and phase change, CTR Summer Program 2024: https://web.stanford.edu/group/ctr/ctrsp24/iv07_CARCANA.pdf
• Wilfong et al., MFC 5.0: An exascale many-physics flow solver, CPC 322 (2026), DOI: https://doi.org/10.1016/j.cpc.2026.110055, open PDF: https://benwilfong.com/papers/wilfong-mfc-26.pdf

AI-assistance disclosure

This issue was prepared with AI assistance. The assistant inspected the cited source files and generated namelist at the audited commit, compared the tensor layout and constitutive formulas against the references above, and ran repository searches for parameter registration and test coverage. The conclusions and proposed fixes should receive maintainer/model-author review before merge.

[sbryngelson]

Implementation workflow is now linked:

• Remove non-operational hyperelasticity capability #1675 is the draft removal PR against master. It deletes the exposed hyperelastic path completely while preserving hypoelasticity.
• Restore constrained neo-Hookean reference-map hyperelasticity #1677 is the draft stacked add-back PR against agent/remove-hyperelasticity. It restores a constrained neo-Hookean, 3D Cartesian, HLLC-only path and applies the audited fixes.

The add-back PR preserves the published MFC/Barbosa/Rodriguez closure tau = (G/J) dev(b), corrects b = F F^T, registers and constrains hyper_model, removes the pseudo-Mooney-Rivlin and bespoke pre-stress paths, defines the missing HLLC fluxes, and adds a deforming 3D shear golden canary (C06DAA3C) with nonzero off-diagonal stress. The shear state distinguishes F F^T from F^T F just as a rotated stretch would.

CPU debug precheck/build and both the new hyperelastic canary and an existing hypoelastic continuum-damage regression pass. GPU parity remains unchecked and is called out explicitly in #1677.

After #1675 merges, #1677 should be retargeted to master.