[drcole17]

Adds a per-fluid temperature-dependent shear viscosity, the Arrhenius (Andrade) liquid law

$$\mu_i(T) = \exp\left(C_i + \frac{D_i}{T}\right),$$

selected with fluid_pp(i)%visc_model = 1 and coefficients fluid_pp(i)%visc_c, fluid_pp(i)%visc_d. The default visc_model = 0 leaves the constant-Re(1) path untouched; no existing golden changes.

Approach

• Temperature. μ is evaluated at the thermal-equilibrium mixture temperature from f_mixture_temperature (m_eos), the same closure heat conduction (Fourier heat conduction in the energy equation #1906) uses. No new temperature model; like k_therm it requires cv > 0 on every fluid.
• Riemann states. In HLL and HLLC (5- and 6-equation paths), right after s_compute_interface_reynolds fills Re_L/Re_R, Re_L(1)/Re_R(1) are replaced by 1/Σ α_i μ_i(T_L/R). Fluids that are not Arrhenius keep 1/Re_i(1) in the sum. Everything downstream (the Re_avg_rs*_vf face average, the viscous source flux) is unchanged, so the face viscosity follows the reconstructed temperature.
• Time step. The same override is applied to the cell Re(1) before s_compute_dt_from_cfl (m_time_steppers) and in both stability-criterion paths in m_data_output, so the viscous CFL limit sees the local μ(T) rather than the unused placeholder Re(1).
• Exponent cap. C + D/T is capped at 50 (visc_exp_arg_cap in m_constants) and T is floored at sgm_eps, so μ stays finite in a transient or under-resolved cold state.

Why this layout (Herschel-Bulkley, #1545)

The parameter plumbing follows HB: per-fluid fluid_pp fields, gathered at init into device arrays (visc_models, visc_cs, visc_ds) and an any_arrhenius_visc flag that are GPU_DECLARE/GPU_UPDATEd next to HB's, with a matching single @:DEALLOCATE. The μ evaluation is a GPU_ROUTINE(seq) function in its own module, m_arrhenius_viscosity, like m_hb_function. visc_model is an enum parameter with CONSTRAINTS names, generating visc_model_constant/visc_model_arrhenius.

The hook point differs from HB on purpose. HB's μ depends on the strain rate, which only exists in m_viscous, so HB overrides there. μ(T) is a pointwise function of each Riemann state, so it overrides where Re_L/Re_R are formed. That reuses the existing face averaging and viscous flux with no new gradient code. For the time step, HB substitutes its bounded mu_max; here the cell value of μ(T) is exact, so that is what's used.

Verification

Analytic: examples/2D_Couette_Variable_Viscosity

Plane Couette flow between a fixed cold wall (300 K) and a sliding hot wall (400 K), with μ(T₀)/μ(T₁) ≈ 3.5, Re ≈ 20, Ma ≈ 0.1, Pr = 1, Br ≈ 0.014. The reference is the exact steady solution of the coupled momentum/energy BVP (scipy.solve_bvp, ~1e-10). MFC runs with model_eqns = 3, HLLC, upstream k_therm conduction and isothermal no-slip walls (bc_y = -16). It starts from rest, so the curved profile is not baked into the initial condition. Single rank, gfortran 15.2.

Ny	L2(u)/U	Linf(u)/U	L2(T)/ΔT	unsteadiness	steps
33	1.3268e-4	1.8302e-4	1.582e-6	1.4e-9	26752
65	3.4067e-5	4.8511e-5	3.927e-7	5.4e-10	56186
97	1.5310e-5	2.1997e-5	1.752e-7	4.6e-11	125125

The observed L2 order in velocity is 2.003. Mid-gap velocity is u(H/2) = 0.3519 U, versus 0.5 U for constant μ. That ~30% deficit is the μ(T) signal.

This feature was first developed on a fork with its own conduction implementation. Against that version's results:

• Velocity errors agree to 2e-5 relative and step counts are identical.
• L2(T) differs by −0.49% at Ny = 33, −2e-5 at Ny = 65 and −6e-5 at Ny = 97.

The temperature difference is an absolute 8e-9 of ΔT at the coarsest grid, and it vanishes with refinement. That is consistent with the switch to #1906's conduction discretization. Velocity, the quantity μ(T) acts on, is unaffected.

python3 validate.py writes the profile, convergence and pointwise-error figures and summary.json. The example is in casesToSkip: its analytic density initial condition compiles into pre_process, so it needs a case-specific build. Test-suite coverage comes from the golden below.

Golden test

New 6D29F444 ("Viscous -> visc_model=1"): 1D, model_eqns = 2, HLLC. The patches start at T = 2.5, 2.5, 2.0, so μ spans e⁹–e¹⁰ around the base 1/Re(1) = 1e4. The golden was generated with gfortran 15.2; an ifx build of the same case agrees to a maximum relative difference of 2e-14.

Suite and toolchain

• ./mfc.sh test --select-enforce (CPU, gfortran 15.2, local mpirun): 723 passed, 0 failed, 39 skipped (762 total). Coverage selection picked the full suite because macro and codegen inputs changed. All 71 Viscous and conduction tests pass, and no existing golden was modified.
• 9 new TestArrheniusViscosity validator tests pass. test_fluid_pp_parameter_set_unchanged is updated for the three new fluid_pp fields.
• ./mfc.sh precheck: 7/7 pass.

Notes for reviewers

• Validator (check_visc_model). It requires all of:

  • viscous = T
  • model_eqns = 2 or 3
  • riemann_solver = 1 or 2
  • cv > 0 on every fluid
  • stiffened-gas or ideal-gas EOS
  • visc_c set, and visc_d set and nonzero

  It rejects visc_c/visc_d without visc_model = 1, and rejects the model combined with non_newtonian, cyl_coord, bubbles_euler, ib, igr or chemistry.

• Re(1) must still be set for an Arrhenius fluid, because viscous requires it, as it does for non_newtonian. Its value is unused for that fluid's shear viscosity. Re(2) (bulk) stays constant.

• weno_Re_flux = F in the example. It is explained in its README. For a single fluid with smooth μ, weno_Re_flux leaves the 2Δy mode of the wall-normal velocity gradient undamped, and with no advection nothing else damps it. The compact central gradient does.

• run_suite.sh is force-added past the examples/**/*.sh ignore rule, as examples/1D_convergence/submitJobs.sh is.

Known limitations

• Isothermal walls. The wall-face viscosity comes from the reflected ghost state, which is the adjacent interior temperature, not Twall. This is documented in case.md; the Couette case still converges at second order.
• Prohibited combinations.
  • non_newtonian: both models would own Re(1).
  • cyl_coord: the cylindrical viscous source in m_riemann_state forms its own Re_s, and it is not overridden or validated here.
  • Riemann solvers other than HLL/HLLC: the override is only in those two.
• Mixture temperature. It is the single equilibrium temperature from the partial-density-weighted heat capacity Σ α_iρ_i c_v,i n_i, so a fluid without cv would drop out of the heat capacity and bias T high. The validator requires cv on every fluid instead.
• Only the liquid law. Sutherland's law (gases), and the differentially heated cavity benchmark that uses it, are a natural follow-up.
• Not verified locally. GPU (OpenACC/OpenMP offload) builds and multi-rank runs of the new paths are exercised only by CI. The new code uses only GPU_* macros and follows the HB device-array pattern.

Size

About 200 lines of source, toolchain and docs. The physics is the 52-line m_arrhenius_viscosity module plus five 2–4 line call sites. The rest is plumbing the toolchain requires:

• the three new fluid_pp fields in three targets' defaults and the parameter registry
• the device-array gather
• the validator rule and its physics-docs entry (68 lines)
• the case.md rows and paragraph

Tests (validator unit tests, golden) and the example are additional.

---

Contribution Policy

Written with Claude Code. The physics is verified against a closed-form solution, as above. It is a port of a feature developed and validated on a fork, rebuilt onto upstream conduction (#1906).

Acknowledgement

• I confirm this PR meets the above expectations and reflects my own understanding and real-world context.

Left unticked deliberately: that attestation is the author's to make, not the tool's.

🤖 Generated with Claude Code

https://claude.ai/code/session_01VtsGVvCfsDJz5dmtSk9dKt

[github-actions]

Claude Code Review

Head SHA: 0191d97

Files changed:

• 25
• docs/documentation/case.md
• docs/module_categories.json
• src/common/m_constants.fpp
• src/common/m_derived_types.fpp
• src/common/m_global_parameters_common.fpp
• src/post_process/m_global_parameters.fpp
• src/pre_process/m_global_parameters.fpp
• src/simulation/m_arrhenius_viscosity.fpp
• src/simulation/m_global_parameters.fpp
• src/simulation/m_riemann_solver_hllc.fpp

Findings:

• high src/common/m_constants.fpp: visc_model_constant and visc_model_arrhenius are used but never defined. m_global_parameters_common.fpp adds them to its use m_constants, only: list. The pre_process, post_process and simulation m_global_parameters.fpp files assign fluid_pp(i)%visc_model = visc_model_constant. m_arrhenius_viscosity.fpp and the simulation m_global_parameters.fpp compare against visc_model_arrhenius. The only m_constants.fpp change in the diff is visc_exp_arg_cap. A grep of src/ finds no definition of either name, including in the post-change context file. The m_derived_types.fpp comment ("visc_model_* in m_constants") also points to constants that do not exist. This will not compile in any of the three targets. Suggested fix: add integer, parameter :: visc_model_constant = 0 and integer, parameter :: visc_model_arrhenius = 1, matching the CONSTRAINTS["names"] in definitions.py. High confidence. The changes may have been left out of the commit.
• high src/pre_process/m_global_parameters.fpp, src/post_process/m_global_parameters.fpp: neither file imports visc_model_constant. The post_process file has use m_constants, only: format_silo, precision_single. Even with the constants defined in m_constants, these two files need them added to their only: lists (or the lists dropped), or they will not compile. High confidence for post_process. For pre_process the use m_constants line was not located in the context file, so check it there.
• medium src/simulation/m_riemann_solver_hll.fpp (the HLL and HLLC hunks): the Arrhenius override sits inside if (viscous) after s_compute_interface_reynolds. It is applied only when any_arrhenius_visc is set. It relies on Re_L(1) and Re_R(1) being the shear Reynolds number. With weno_Re_flux = F that is the case. Whether the override runs for every path that fills Re_L and Re_R (for example the HLLC branches at the two hunks in m_riemann_solver_hllc.fpp) depends on the unchanged surrounding code. The validator does not reject Arrhenius with weno_Re_flux = T or other viscous-flux paths. Medium-low confidence; worth confirming against the viscous flux code.
• medium src/simulation/m_arrhenius_viscosity.fpp: f_compute_arrhenius_shear_re takes alpha and alpha_rho as dimension(num_fluids), or dimension(3) under AMD without case optimization. Callers pass alpha_L and alpha_rho_L, which are declared elsewhere. The m_data_output.fpp and m_time_steppers.fpp callers pass alpha and alpha_rho from their local declarations. If any of those is sized differently, such as num_fluids_max, the AMD path breaks or gets an explicit-shape mismatch. This hunk is at least consistent with the neighbouring s_compute_stability code, but it was not verified. Medium-low confidence.
• medium src/simulation/m_global_parameters.fpp: visc_models, visc_cs and visc_ds are allocated with @:ALLOCATE and deallocated. The scalar any_arrhenius_visc is declared with GPU_DECLARE(create=...) and updated with GPU_UPDATE(device=...), which is consistent. visc_cs and visc_ds hold dflt_real for non-Arrhenius fluids, but they are only read when visc_models(fl) == visc_model_arrhenius, so that is fine. No problem found. Listed for completeness only.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_arrhenius_viscosity.fpp	32	+32
src/simulation/m_global_parameters.fpp	812	+17
src/simulation/m_riemann_solver_hllc.fpp	1321	+13
src/simulation/m_riemann_solver_hll.fpp	614	+7
src/simulation/m_data_output.fpp	1584	+6
src/simulation/m_time_steppers.fpp	914	+4
src/common/m_derived_types.fpp	486	+3
src/post_process/m_global_parameters.fpp	417	+3
src/pre_process/m_global_parameters.fpp	490	+3
src/common/m_constants.fpp	90	+1
src/common/m_global_parameters_common.fpp	251	+1

Directory	Lines	Diff
common	10531	+5
pre_process	5037	+3
simulation	28610	+79
post_process	3502	+3
total	47680	+90

[codecov]

Codecov Report

❌ Patch coverage is 59.45946% with 15 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.88%. Comparing base (3ef5e1e) to head (0e1527e).

Files with missing lines	Patch %	Lines
src/simulation/m_riemann_solver_hllc.fpp	0.00%	5 Missing ⚠️
src/simulation/m_data_output.fpp	25.00%	2 Missing and 1 partial ⚠️
src/simulation/m_riemann_solver_hll.fpp	0.00%	3 Missing ⚠️
src/simulation/m_arrhenius_viscosity.fpp	75.00%	1 Missing and 1 partial ⚠️
src/simulation/m_time_steppers.fpp	0.00%	1 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2101      +/-   ##
==========================================
- Coverage   61.89%   61.88%   -0.01%     
==========================================
  Files          86       87       +1     
  Lines       22888    22927      +39     
  Branches     3372     3378       +6     
==========================================
+ Hits        14166    14188      +22     
- Misses       6220     6234      +14     
- Partials     2502     2505       +3     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.