Anchored Dual-pass HLLD for Hypoelasticity (+ HLLC and interface-consistent HLL)

docs/documentation/case.md

@@ -515,7 +515,8 @@ See @ref equations "Equations" for the mathematical models these parameters cont
 | `flux_lim`                 | Integer | Flux limiter for post-process: [1] minmod; [2] MUSCL; [3] OSPRE; [4] SUPERBEE |
 | `ic_eps`                   | Real    | Interface compression threshold (default: 1e-4) |
 | `ic_beta`                  | Real    | Interface compression sharpness parameter (default: 1.6) |
-| `riemann_solver`           | Integer | Riemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact*; [4] HLLD	(only for MHD) |
+| `riemann_solver`           | Integer | Riemann solver algorithm: [1] HLL*; [2] HLLC; [3] Exact*; [4] HLLD (MHD or hypoelasticity) |
+| `hll_u_interface`          | Logical | HLL Method 2 (u-interface) for volume fraction advection (default F) |
 | `low_Mach`                 | Integer | Low Mach number correction for HLLC Riemann solver: [0] None; [1] Pressure (\cite Chen22); [2] Velocity (\cite Thornber08)	 |
 | `avg_state`	               | Integer | Averaged state evaluation method: [1] Roe average*; [2] Arithmetic mean  |
 | `wave_speeds`              | Integer | Wave-speed estimation: [1] Direct (\cite Batten97); [2] Pressure-velocity* (\cite Toro09)	 |
@@ -536,6 +537,9 @@ See @ref equations "Equations" for the mathematical models these parameters cont
 | `surface_tension`          | Logical | Activate surface tension |
 | `viscous`                  | Logical | Activate viscosity |
 | `hypoelasticity`           | Logical | Activate hypoelasticity* |
+| `riemann_hypo_ADC`         | Logical | Enable hypo anti-diffusion correction for HLLC/HLLD (default F) |
+| `ADC_kappa`                | Real    | ADC sensor scaling parameter (default 1.0) |
+| `hypo_hll_interface_rhs`   | Logical | HLL uses interface-consistent hypo RHS (default F) |
 | `igr`                      | Logical | Enable solution via information geometric regularization (IGR) \cite Cao24 |
 | `igr_order`                | Integer | Order of reconstruction for IGR [3,5] |
 | `alf_factor`               | Real    | Alpha factor for IGR entropic pressure (default 10) |
@@ -569,7 +573,7 @@ Tangential velocities require viscosity, `weno_avg = T`, and `bc_[x,y,z]%%end =
 The difference of the two models is assessed by (\cite Schmidmayer20).
 Note that some code parameters are only compatible with 5-equation model.
 
-- `alt_soundspeed` activates the source term in the advection equations for the volume fractions, \f$K\nabla\cdot \underline{u}\f$, that regularizes the speed of sound in the mixture region when the 5-equation model is used.
+- `alt_soundspeed` activates the source term in the advection equations for the volume fractions, \f$K\nabla\cdot \underline{u}\f$, that regularizes the speed of sound in the mixture region when the 5-equation model is used. Requires exactly two fluid components (\f$K\f$ is a two-fluid closure).
 The effect and use of the source term are assessed by \cite Schmidmayer20.
 
 - `adv_n` activates the direct computation of number density by the Riemann solver instead of computing number density from the void fraction in the method of classes.
@@ -615,7 +619,11 @@ Setting `muscl_eps = 0` gives textbook limiter behavior where limiters activate
 
 - `riemann_solver` specifies the choice of the Riemann solver that is used in simulation by an integer from 1 through 4.
 `riemann_solver = 1`, `2`, and `3` correspond to HLL, HLLC, and Exact Riemann solver, respectively (\cite Toro09).
-`riemann_solver = 4` is only for MHD simulations. It resolves 5 of the full seven-wave structure of the MHD equations (\cite Miyoshi05).
+`riemann_solver = 4` is the HLLD solver for MHD or hypoelasticity simulations. For MHD it resolves 5 of the full seven-wave structure of the MHD equations (\cite Miyoshi05).
+
+- `hll_u_interface`: Selects between two HLL discretizations of volume fraction advection (`riemann_solver = 1`):
+  - **Default** (``'F'``): \f$\partial_t \alpha_k + u\,\partial_x \alpha_k = 0\f$
+  - **u-interface** (``'T'``, consistent with HLLC): \f$\partial_t \alpha_k + \partial_x(\alpha_k\, u) = \alpha_k\,\partial_x u\f$
 
 - `low_Mach` specifies the choice of the low Mach number correction scheme for the HLLC Riemann solver. `low_Mach = 0` is default value and does not apply any correction scheme. `low_Mach = 1` and `2` apply the anti-dissipation pressure correction method (\cite Chen22) and the improved velocity reconstruction method (\cite Thornber08). This feature requires `model_eqns = 2` or `3`. `low_Mach = 1` works for `riemann_solver = 1` and `2`, but `low_Mach = 2` only works for `riemann_solver = 2`.
 
@@ -635,7 +643,13 @@ This option requires `weno_Re_flux` to be true because cell boundary values are
 
 - `viscous` activates viscosity when set to ``'T'``. Requires `Re(1)` and `Re(2)` to be set.
 
-- `hypoelasticity` activates elastic stress calculations for fluid-solid interactions. Requires `G` to be set in the fluid material's parameters.
+- `hypoelasticity` activates elastic stress calculations for fluid-solid interactions. Requires `G` to be set in `fluid_pp`, and `fd_order` to be set (1, 2, or 4). Compatible with HLL (`riemann_solver = 1`), HLLC (`riemann_solver = 2`), and HLLD (`riemann_solver = 4`). The Riemann solver choice determines how the elastic stress source term \f$\mathbf{S}^e\f$ is discretized:
+  - **HLL**: finite-difference velocity gradient (default), or interface-consistent velocity gradient when ``hypo_hll_interface_rhs = 'T'`` (matches HLLC).
+  - **HLLC**: interface-consistent velocity gradient from the Riemann solution.
+  - **HLLD**: dual-pass approach resolving the elastic wave structure. Requires 2D+ and exactly 2 fluid components.
+  - With hypoelastic HLLD, use characteristic boundary conditions (`bc_{x,y,z}%%beg/end` between -5 and -12) only where the domain boundary remains in a fluid region throughout the simulation. Keep material interfaces and solids out of the nearby cells used by the boundary calculation. A material present only at the numerical volume-fraction floor is allowed, even if it has nonzero `G`. These boundaries treat acoustic waves only. `alt_soundspeed` is supported.
+
+- `riemann_hypo_ADC`: Enables anti-diffusion correction (ADC) for hypoelastic HLLC or HLLD. Blends the HLLC/HLLD flux locally toward the more diffusive HLL flux where a jump sensor (total normal stress, tangential velocity, and tangential stress) detects a strong contact or material interface, improving robustness and reducing interfacial overshoots. This trades contact and shear resolution for robustness where the sensor activates. `ADC_kappa` (default 1.0) scales the reference jump the sensor is normalized by; smaller values blend more toward HLL (more diffusive and robust). Off by default and recommended off: in MFC's tested hypoelastic cases, including strong shock-interface interactions, HLLD is stable and accurate without ADC. The option is a robustness fallback intended for regimes with intense shocks and intense shear at material interfaces.
 
 #### Boundary Condition Patches {#boundary-condition-patches}
 

docs/documentation/equations.md

@@ -874,16 +874,33 @@ Four-state solver resolving the contact discontinuity. Star-state satisfies:
 
 Iterative exact Riemann solver.
 
-#### HLLD (`riemann_solver = 4`, MHD only)
+#### HLLD (`riemann_solver = 4`, MHD or hypoelasticity)
 
-Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contact waves (\cite Miyoshi05). The Riemann fan is divided by outer wave speeds \f$S_L\f$, \f$S_R\f$, Alfven speeds \f$S_L^*\f$, \f$S_R^*\f$, and a middle contact \f$S_M\f$:
+**MHD HLLD.** Seven-state solver for ideal MHD resolving fast magnetosonic, Alfven, and contact waves (\cite Miyoshi05). The Riemann fan is divided by outer wave speeds \f$S_L\f$, \f$S_R\f$, Alfven speeds \f$S_L^*\f$, \f$S_R^*\f$, and a middle contact \f$S_M\f$:
 
 \f[S_M = \frac{(S_R - u_R)\rho_R u_R - (S_L - u_L)\rho_L u_L - p_{T,R} + p_{T,L}}{(S_R - u_R)\rho_R - (S_L - u_L)\rho_L}\f]
 
 \f[S_L^* = S_M - \frac{|B_x|}{\sqrt{\rho_L^*}}, \qquad S_R^* = S_M + \frac{|B_x|}{\sqrt{\rho_R^*}}\f]
 
 where \f$p_T = p + |\mathbf{B}|^2/2\f$ is the total (thermal + magnetic) pressure. Continuity of normal velocity and total pressure is enforced across the Riemann fan.
 
+**Hypoelastic HLLD.** Shares the five-wave Riemann fan structure with MHD HLLD, but differs substantially in formulation. Elastic shear waves play the role of Alfven waves, with total pressure \f$p_T = p - \tau_{nn}\f$ and shear wave impedance \f$C = \hat{\rho}(\hat{G} + \hat{\tau}_{nn})\f$. However, the non-conservative nature of the multi-component elastic stress equations — particularly at material interfaces where \f$G\f$ is discontinuous — requires a dedicated treatment distinct from the MHD solver.
+
+The solver uses an **anchored dual-pass** formulation. Each Riemann solve at interface \f$j{+}1/2\f$ requires cell-centered quantities for the non-conservative products; these are taken as cell-centered values from the cell that the flux will update, rather than from the reconstructed interface states:
+
+\f[\mathbf{F}_{j+1/2}^{\text{left}} = \text{HLLD}\!\left(\mathbf{q}_{j+1/2}^L,\;\mathbf{q}_{j+1/2}^R;\;\mathbf{q}_{\text{cell},\,j}\right), \qquad \mathbf{F}_{j-1/2}^{\text{right}} = \text{HLLD}\!\left(\mathbf{q}_{j-1/2}^L,\;\mathbf{q}_{j-1/2}^R;\;\mathbf{q}_{\text{cell},\,j}\right)\f]
+
+\f[\frac{d\mathbf{U}_j}{dt} = \frac{1}{\Delta x}\!\left(\mathbf{F}_{j-1/2}^{\text{right}} - \mathbf{F}_{j+1/2}^{\text{left}}\right)\f]
+
+This formulation enables HLLD to be used with the non-conservative terms that affect the eigenstructure of the quasi-linear Jacobian. Volume fraction advection is built into the HLLD flux rather than treated as a separate non-conservative step.
+
+#### Developer Notes on Non-Conservative Term Implementation
+
+For details on how the Riemann solvers discretize non-conservative terms (volume fraction advection, Kapila \f$K\,\nabla\!\cdot\!\mathbf{u}\f$, and hypoelastic velocity gradients) and hand off interface quantities to the RHS, see the notes in `misc/dev_notes/`:
+
+- `HLL_HLLC_non_conservative_terms_derivations.md` — Mathematical derivation of HLL Method 1 (alpha-interface), Method 2 (u-interface), and HLLC transport traces, including the \f$S_M\zeta_K\f$ star-branch construction and ADC blending.
+- `Riemann_and_RHS_source_terms_explanations.md` — Code dataflow: which arrays carry what between m_riemann_solvers and m_rhs, overloading of flux_src, the three NC advection modes (adv_src_mode_alpha_iface, adv_src_mode_vel_iface, adv_src_mode_none), and the nc_iface_vel second export channel.
+
 ### 15.3 Time Integration
 
 **Source:** `src/simulation/m_time_steppers.fpp`

docs/module_categories.json

@@ -11,6 +11,7 @@
             "m_riemann_solver_hll",
             "m_riemann_solver_lf",
             "m_riemann_solver_hllc",
+            "m_riemann_solver_hypo_hlld",
             "m_muscl",
             "m_variables_conversion",
             "m_thinc"

examples/1D_sodHypo/case.py

@@ -48,6 +48,7 @@
             "bc_x%end": -3,
             # Turning on Hypoelasticity
             "hypoelasticity": "T",
+            "fd_order": 4,
             # Formatted Database Files Structure Parameters
             "format": "silo",
             "precision": "double",

examples/2D_axisym_hypo_hlld/case.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+import json
+
+# 2D axisymmetric hypoelastic case with HLLD solver.
+# Epoxy rectangular inclusion in water hit by focused acoustic pulse.
+# Water: gamma=4.4, pi_inf=6.0e8, G=0
+# Epoxy: gamma=4.4, pi_inf=2.4e9, rho=1180, G=1.5e9
+
+gamma_w, pi_inf_w, rho_w = 4.4, 6.0e8, 1000.0
+gamma_e, pi_inf_e, rho_e = 4.4, 2.4e9, 1180.0
+
+config = {
+    "run_time_info": "T",
+    # Computational Domain
+    "x_domain%beg": 0.0,
+    "x_domain%end": 1.0,
+    "y_domain%beg": 0.0,
+    "y_domain%end": 1.0,
+    "cyl_coord": "T",
+    "m": 99,
+    "n": 99,
+    "p": 0,
+    "dt": 1.5e-6,
+    "t_step_start": 0,
+    "t_step_stop": 400,
+    "t_step_save": 100,
+    # Simulation Algorithm
+    "num_patches": 2,
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "num_fluids": 2,
+    "mpp_lim": "F",
+    "mixture_err": "F",
+    "time_stepper": 3,
+    "weno_order": 5,
+    "weno_eps": 1.0e-16,
+    "mapped_weno": "T",
+    "null_weights": "F",
+    "mp_weno": "F",
+    "riemann_solver": 4,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "bc_x%beg": -6,
+    "bc_x%end": -6,
+    "bc_y%beg": -2,
+    "bc_y%end": -6,
+    # Output
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "rho_wrt": "T",
+    "parallel_io": "T",
+    # Hypoelasticity
+    "hypoelasticity": "T",
+    "fd_order": 4,
+    # Patch 1: Water background
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": 0.5,
+    "patch_icpp(1)%y_centroid": 0.5,
+    "patch_icpp(1)%length_x": 1.0,
+    "patch_icpp(1)%length_y": 1.0,
+    "patch_icpp(1)%vel(1)": 0.0,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%pres": 1e05,
+    "patch_icpp(1)%tau_e(1)": 0.0,
+    "patch_icpp(1)%alpha_rho(1)": rho_w * (1.0 - 1e-8),
+    "patch_icpp(1)%alpha(1)": 1.0 - 1e-8,
+    "patch_icpp(1)%alpha_rho(2)": rho_e * 1e-8,
+    "patch_icpp(1)%alpha(2)": 1e-8,
+    # Patch 2: Epoxy rectangular inclusion
+    "patch_icpp(2)%alter_patch(1)": "T",
+    "patch_icpp(2)%geometry": 3,
+    "patch_icpp(2)%x_centroid": 0.6,
+    "patch_icpp(2)%y_centroid": 0.2,
+    "patch_icpp(2)%length_x": 0.2,
+    "patch_icpp(2)%length_y": 0.2,
+    "patch_icpp(2)%smoothen": "T",
+    "patch_icpp(2)%smooth_patch_id": 1,
+    "patch_icpp(2)%smooth_coeff": 2.0,
+    "patch_icpp(2)%vel(1)": 0.0,
+    "patch_icpp(2)%vel(2)": 0.0,
+    "patch_icpp(2)%pres": 1e05,
+    "patch_icpp(2)%tau_e(1)": 0.0,
+    "patch_icpp(2)%alpha_rho(1)": rho_w * 1e-8,
+    "patch_icpp(2)%alpha(1)": 1e-8,
+    "patch_icpp(2)%alpha_rho(2)": rho_e * (1.0 - 1e-8),
+    "patch_icpp(2)%alpha(2)": 1.0 - 1e-8,
+    # Acoustic source (axisymmetric focused)
+    "acoustic_source": "T",
+    "num_source": 1,
+    "acoustic(1)%support": 6,
+    "acoustic(1)%loc(1)": 0.1,
+    "acoustic(1)%loc(2)": 0.0,
+    "acoustic(1)%pulse": 2,
+    "acoustic(1)%npulse": 1,
+    "acoustic(1)%mag": 1.0,
+    "acoustic(1)%foc_length": 0.8,
+    "acoustic(1)%aperture": 0.8,
+    "acoustic(1)%gauss_sigma_time": 4e-5,
+    "acoustic(1)%delay": 2e-4,
+    # Fluids Physical Parameters
+    "fluid_pp(1)%gamma": 1.0 / (gamma_w - 1.0),
+    "fluid_pp(1)%pi_inf": gamma_w * pi_inf_w / (gamma_w - 1.0),
+    "fluid_pp(1)%G": 0.0,
+    "fluid_pp(2)%gamma": 1.0 / (gamma_e - 1.0),
+    "fluid_pp(2)%pi_inf": gamma_e * pi_inf_e / (gamma_e - 1.0),
+    "fluid_pp(2)%G": 1.5e9,
+}
+
+print(json.dumps(config, indent=4))

examples/2D_hypo_hlld/case.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+import json
+
+# 2D hypoelastic case with HLLD solver.
+# Epoxy square in water hit by focused acoustic pulse.
+# Water: gamma=4.4, pi_inf=6.0e8, G=0
+# Epoxy: gamma=4.4, pi_inf=2.4e9, rho=1180, G=1.5e9
+
+gamma_w, pi_inf_w, rho_w = 4.4, 6.0e8, 1000.0
+gamma_e, pi_inf_e, rho_e = 4.4, 2.4e9, 1180.0
+
+config = {
+    "run_time_info": "T",
+    # Computational Domain
+    "x_domain%beg": 0.0,
+    "x_domain%end": 1.0,
+    "y_domain%beg": 0.0,
+    "y_domain%end": 1.0,
+    "m": 99,
+    "n": 99,
+    "p": 0,
+    "dt": 1.5e-6,
+    "t_step_start": 0,
+    "t_step_stop": 400,
+    "t_step_save": 100,
+    # Simulation Algorithm
+    "num_patches": 2,
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "num_fluids": 2,
+    "mpp_lim": "F",
+    "mixture_err": "F",
+    "time_stepper": 3,
+    "weno_order": 5,
+    "weno_eps": 1.0e-16,
+    "mapped_weno": "T",
+    "null_weights": "F",
+    "mp_weno": "F",
+    "riemann_solver": 4,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "bc_x%beg": -6,
+    "bc_x%end": -6,
+    "bc_y%beg": -6,
+    "bc_y%end": -6,
+    # Output
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "rho_wrt": "T",
+    "parallel_io": "T",
+    # Hypoelasticity
+    "hypoelasticity": "T",
+    "fd_order": 4,
+    # Patch 1: Water background
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": 0.5,
+    "patch_icpp(1)%y_centroid": 0.5,
+    "patch_icpp(1)%length_x": 1.0,
+    "patch_icpp(1)%length_y": 1.0,
+    "patch_icpp(1)%vel(1)": 0.0,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%pres": 1e05,
+    "patch_icpp(1)%tau_e(1)": 0.0,
+    "patch_icpp(1)%alpha_rho(1)": rho_w * (1.0 - 1e-8),
+    "patch_icpp(1)%alpha(1)": 1.0 - 1e-8,
+    "patch_icpp(1)%alpha_rho(2)": rho_e * 1e-8,
+    "patch_icpp(1)%alpha(2)": 1e-8,
+    # Patch 2: Epoxy square
+    "patch_icpp(2)%alter_patch(1)": "T",
+    "patch_icpp(2)%geometry": 3,
+    "patch_icpp(2)%x_centroid": 0.6,
+    "patch_icpp(2)%y_centroid": 0.5,
+    "patch_icpp(2)%length_x": 0.2,
+    "patch_icpp(2)%length_y": 0.2,
+    "patch_icpp(2)%smoothen": "T",
+    "patch_icpp(2)%smooth_patch_id": 1,
+    "patch_icpp(2)%smooth_coeff": 2.0,
+    "patch_icpp(2)%vel(1)": 0.0,
+    "patch_icpp(2)%vel(2)": 0.0,
+    "patch_icpp(2)%pres": 1e05,
+    "patch_icpp(2)%tau_e(1)": 0.0,
+    "patch_icpp(2)%alpha_rho(1)": rho_w * 1e-8,
+    "patch_icpp(2)%alpha(1)": 1e-8,
+    "patch_icpp(2)%alpha_rho(2)": rho_e * (1.0 - 1e-8),
+    "patch_icpp(2)%alpha(2)": 1.0 - 1e-8,
+    # Acoustic source (2D focused arc)
+    "acoustic_source": "T",
+    "num_source": 1,
+    "acoustic(1)%support": 5,
+    "acoustic(1)%loc(1)": 0.1,
+    "acoustic(1)%loc(2)": 0.5,
+    "acoustic(1)%pulse": 2,
+    "acoustic(1)%npulse": 1,
+    "acoustic(1)%dir": 1.0,
+    "acoustic(1)%mag": 1.0,
+    "acoustic(1)%foc_length": 0.4,
+    "acoustic(1)%aperture": 0.75,
+    "acoustic(1)%gauss_sigma_time": 4e-5,
+    "acoustic(1)%delay": 2e-4,
+    # Fluids Physical Parameters
+    "fluid_pp(1)%gamma": 1.0 / (gamma_w - 1.0),
+    "fluid_pp(1)%pi_inf": gamma_w * pi_inf_w / (gamma_w - 1.0),
+    "fluid_pp(1)%G": 0.0,
+    "fluid_pp(2)%gamma": 1.0 / (gamma_e - 1.0),
+    "fluid_pp(2)%pi_inf": gamma_e * pi_inf_e / (gamma_e - 1.0),
+    "fluid_pp(2)%G": 1.5e9,
+}
+
+print(json.dumps(config, indent=4))