Add Herschel-Bulkley non-Newtonian viscosity model

[sbryngelson]

Summary

Adds per-fluid Herschel–Bulkley (HB) non-Newtonian viscosity with Papanastasiou regularization to the simulation target. The effective viscosity is computed from the local strain-rate tensor every step, so it varies in space and time. Fully backward-compatible: all new parameters default to Newtonian behavior, and when no fluid is non-Newtonian every code path is bitwise-identical to today.

This is a clean reimplementation on current master, replacing the stale #1298 (which could no longer be merged after master refactored the Riemann reconstruction buffers and IB types). It reuses the validated HB formula from that work and routes the shear-dependent viscosity through MFC's existing viscous machinery rather than recomputing gradients — which structurally avoids the rotated→physical index-mapping bug that affected #1298.

What it can model

Effective viscosity:

μ(γ̇) = (τ₀/γ̇)(1 − e^(−m·γ̇)) + K·γ̇^(n−1),   γ̇ = √(2 Dᵢⱼ Dᵢⱼ)

Regime	Parameters
Newtonian	tau0=0, nn=1 (μ = K)
Power-law (shear-thinning / -thickening)	tau0=0, nn≠1
Bingham plastic	nn=1, tau0>0
General Herschel–Bulkley	tau0>0, nn≠1

Parameters (per fluid)

fluid_pp(i)%: non_newtonian (enable), K (consistency), nn (flow index), tau0 (yield stress), hb_m (Papanastasiou regularization), mu_min/mu_max (viscosity clamps; mu_max required), mu_bulk (optional bulk viscosity). All HB parameters are nondimensional (scaled by ρ_ref U_ref L_ref), so 1/μ equals the local effective Reynolds number — consistent with how Re(1) is used for Newtonian fluids. See the new docs/documentation/case.md section.

Scope / integration

Shear-dependent viscosity is injected into MFC's existing viscous paths, all gated behind a global any_non_newtonian flag:

• Cartesian & cylindrical viscous flux (HLL riemann_solver=1, HLLC =2) — via the shared source-flux routines, reusing the interface strain tensor and its existing physical indexing.
• Immersed boundaries — per-stencil-sample viscosity inside the stress-tensor routine (each sample uses its own local shear rate).
• Viscous CFL/dt — bounded by mu_max for a guaranteed-stable timestep.
• Validator guards (fail fast at ./mfc.sh validate): requires viscous, K, nn, mu_max; tau0>0 requires hb_m; incompatible with igr; model_eqns ∈ {2,3}; riemann_solver ∈ {1,2}.

Validation

Check	Result
Power-law Poiseuille vs analytic profile	0.68% L2, correct shear-thinning bluntness, momentum balance to 0.1%
Bingham/HB Poiseuille vs analytic plug-flow	2.5% L2, rigid plug at y₀ = τ₀/(ρg)
GPU build (nvfortran / OpenACC, V100)	compiles
GPU runtime (NN cases on V100 vs CPU goldens)	4/4 pass
Viscous (63) + IBM (48) + NN (4) regression	pass, zero changes to existing goldens
Newtonian invariance	bitwise-identical when any_non_newtonian = .false.
Case-optimized build	compiles

Examples added: 2D_poiseuille_nn (power-law + analytic check), 2D_bingham_poiseuille_nn (yield-stress + plug-flow check), 2D_lid_driven_cavity_nn (Li et al. 2015).

Known limitations (honest)

• Non-Newtonian supported only with HLL/HLLC (riemann_solver 1, 2); LF/HLLD are guarded off.
• Cylindrical effective shear rate uses the grid-direction strain components (curvature corrections to γ̇ not yet included).
• Cylindrical, IBM, multi-fluid-NN-mixture, and 3D paths are covered by shared-code reuse + Newtonian invariance + regression goldens, not by independent analytic benchmarks. Shear-thickening (n>1) has a regression golden but no analytic gate.

Opening as draft for review.

[github-actions]

Claude Code Review

Head SHA: f1a27cf

Files changed:

• 23
• src/simulation/m_hb_function.fpp
• src/simulation/m_global_parameters.fpp
• src/simulation/m_riemann_solvers.fpp
• src/simulation/m_viscous.fpp
• src/simulation/m_time_steppers.fpp
• src/simulation/m_data_output.fpp
• src/common/m_derived_types.fpp
• toolchain/mfc/case_validator.py
• toolchain/mfc/params/definitions.py
• toolchain/mfc/test/cases.py

Findings

NaN in f_compute_hb_viscosity when tau0 = 0 and hb_m is not set

File: src/simulation/m_hb_function.fpp, lines 32–36

if (shear_rate <= verysmall) then
    yield_term = tau0*hb_m_val
else
    yield_term = tau0*(1._wp - exp(-hb_m_val*shear_rate))/shear_rate
end if

When tau0 = 0 (the default, set in all three m_global_parameters.fpp init blocks) and hb_m is not provided by the user (hb_m_val = dflt_real, a large-negative sentinel), the else branch evaluates:

• -hb_m_val * shear_rate = +|dflt_real| * shear_rate → large positive → exp(+∞) = +Inf
• 1._wp - Inf = -Inf
• tau0 * (-Inf) = 0.0 * (-Inf) → NaN (IEEE 754)

The shear_rate <= verysmall branch is safe (0 * dflt_real = 0, finite × finite), but any cell with nonzero shear rate silently produces NaN for all downstream viscosity and stress computations.

The validator (toolchain/mfc/case_validator.py, new check_non_newtonian method) only requires hb_m when tau0 > 0, leaving the tau0 = 0 (pure power-law) path unprotected. All provided examples happen to set hb_m even with tau0 = 0, but a user following the validator's rules literally can omit it and trigger this path.

Fix option A — guard the yield term in f_compute_hb_viscosity:

if (tau0 == 0._wp) then
    yield_term = 0._wp
else if (shear_rate <= verysmall) then
    yield_term = tau0*hb_m_val
else
    yield_term = tau0*(1._wp - exp(-hb_m_val*shear_rate))/shear_rate
end if

Fix option B — require hb_m unconditionally in check_non_newtonian:

self.prohibit(self.get(f"fluid_pp({i})%hb_m") is None,
              f"fluid_pp({i})%non_newtonian requires fluid_pp({i})%hb_m to be set")

Option A is safer because it defends the function regardless of validator coverage.

[sbryngelson]

Thanks for the review — responses to both findings:

Finding 1 (mu_bulk silently non-functional) — correct, good catch. hb_mu_bulk was indeed unreachable for non-Newtonian fluids: they're absent from Re_idx(2,…) (they set Re(1)=1/K but not Re(2)), so the bulk branch in s_compute_mixture_inv_re never processed them, and the CFL block only bounded Re(1). The documented knob was a silent no-op with zero test coverage. Rather than wire up an untested bulk-viscosity path, I've forbidden it in the case validator — setting fluid_pp(i)%mu_bulk on a non-Newtonian fluid is now rejected ("not yet supported"), and the parameter is marked reserved in case.md. Bulk non-Newtonian viscosity can be a focused follow-up.

Finding 2 (riemann_solver guard short-circuits when the key is absent) — not a live issue. It's already covered upstream: check_riemann_solver does prohibit(riemann_solver is None, "riemann_solver must be specified…"), so any case omitting riemann_solver fails validation regardless of the non-Newtonian check. The is not None short-circuit therefore can't be reached with riemann_solver actually absent — there's no path where a non-Newtonian case runs with an unspecified solver. Left as-is.

[Copilot]

Pull request overview

Adds a per-fluid Herschel–Bulkley (HB) non-Newtonian viscosity option (with Papanastasiou regularization) to the simulation target, integrating shear-dependent viscosity into existing viscous flux / stress machinery, plus toolchain validation and example/benchmark cases.

Changes:

• Introduces HB utilities (m_hb_function) for viscosity evaluation, shear-rate computation, and mixture effective Reynolds-number computation.
• Threads HB-aware viscosity into key simulation paths (viscous stresses/fluxes, dt/CFL stability reporting, and MPI parameter broadcast), gated by any_non_newtonian.
• Extends toolchain parameter inventory/validation and adds/updates non-Newtonian example cases + test goldens.

Reviewed changes

Copilot reviewed 46 out of 48 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
toolchain/mfc/test/cases.py	Adds NN smoke/regression cases and skips slow NN examples in CI example sweep.
toolchain/mfc/params/descriptions.py	Adds human-readable descriptions for new fluid_pp(i)%* NN parameters.
toolchain/mfc/params/definitions.py	Registers new NN parameters in the toolchain parameter registry.
toolchain/mfc/params_tests/negative_tests.py	Adds negative validator tests for NN constraints.
toolchain/mfc/case_validator.py	Adds check_non_newtonian constraints and hooks it into simulation validation.
tests/EE8D38CF/golden-metadata.txt	Adds golden metadata for NN test generation.
tests/E18624FD/golden-metadata.txt	Adds golden metadata for NN test generation.
tests/A582B270/golden-metadata.txt	Adds golden metadata for NN-related test.
tests/78EB6879/golden-metadata.txt	Adds golden metadata for NN test generation.
tests/5304E59F/golden-metadata.txt	Adds golden metadata for NN-related test.
tests/122B0BF3/golden-metadata.txt	Adds golden metadata for NN test generation.
tests/00C20EE8/golden-metadata.txt	Adds golden metadata for NN-related test.
src/simulation/m_viscous.fpp	Uses HB utilities to compute axis-cell effective viscosity/Re when NN is enabled; updates IBM stress tensor to use local NN viscosity.
src/simulation/m_time_steppers.fpp	Adjusts dt/CFL computation to use NN viscosity bounds when enabled.
src/simulation/m_riemann_solvers.fpp	Extends viscous source-flux routines to compute local NN effective Re at interfaces.
src/simulation/m_mpi_proxy.fpp	Broadcasts new NN fluid parameters and flags across ranks.
src/simulation/m_hb_function.fpp	New module providing HB viscosity formula, shear-rate, and mixture effective-Re routines (GPU-routine compatible).
src/simulation/m_global_parameters.fpp	Adds any_non_newtonian and per-fluid HB parameter arrays; populates and GPU-updates them at startup.
src/simulation/m_data_output.fpp	Updates stability-criteria output to account for NN viscosity bounds when enabled.
src/pre_process/m_mpi_proxy.fpp	Broadcasts NN fluid parameters/flag in pre_process.
src/pre_process/m_global_parameters.fpp	Initializes new NN fields in fluid_pp defaults.
src/post_process/m_mpi_proxy.fpp	Broadcasts NN fluid parameters/flag in post_process.
src/post_process/m_global_parameters.fpp	Initializes new NN fields in fluid_pp defaults.
src/common/m_derived_types.fpp	Extends physical_parameters with NN fields (flag + HB/clamp parameters).
examples/2D_poiseuille_thickening_nn/README.md	Documents shear-thickening power-law Poiseuille validation case.
examples/2D_poiseuille_thickening_nn/compare_analytic.py	Adds analytic comparison script for shear-thickening Poiseuille.
examples/2D_poiseuille_thickening_nn/case.py	Adds shear-thickening power-law Poiseuille input case.
examples/2D_poiseuille_nn/README.md	Documents shear-thinning power-law Poiseuille validation case.
examples/2D_poiseuille_nn/compare_analytic.py	Adds analytic comparison script for shear-thinning Poiseuille.
examples/2D_poiseuille_nn/case.py	Adds shear-thinning power-law Poiseuille input case.
examples/2D_lid_driven_cavity_nn/README.md	Documents qualitative lid-driven cavity NN example.
examples/2D_lid_driven_cavity_nn/case.py	Adds lid-driven cavity NN case.
examples/2D_herschel_bulkley_poiseuille_nn/README.md	Documents general HB (yield + power-law) Poiseuille validation case.
examples/2D_herschel_bulkley_poiseuille_nn/compare_analytic.py	Adds analytic comparison script for general HB Poiseuille.
examples/2D_herschel_bulkley_poiseuille_nn/case.py	Adds general HB Poiseuille input case.
examples/2D_bingham_poiseuille_nn/README.md	Documents Bingham (yield-stress) Poiseuille validation case.
examples/2D_bingham_poiseuille_nn/compare_analytic.py	Adds analytic comparison script for Bingham Poiseuille.
examples/2D_bingham_poiseuille_nn/case.py	Adds Bingham Poiseuille input case.
docs/references.bib	Adds the Papanastasiou (1987) reference entry.
docs/module_categories.json	Registers m_hb_function under Physics Models for documentation tooling.
docs/documentation/case.md	Documents the new NN viscosity inputs, formula, and usage notes.

[sbryngelson]

Comprehensive Multi-Agent Review (round 2, HEAD 466c8917)

Four specialist review agents (general code review, test coverage, silent-failure audit, comment/doc accuracy) ran against the full feature diff. All prior-cycle fixes were re-verified as holding. Every finding below was independently audited against the code before posting.

Verdict

Code correctness: clean. The general code review found no Critical or Important defects: strain-rate assembly is correct at all three sites (Cartesian, cylindrical, IBM), the Newtonian path is bitwise-unchanged, precision/memory/MPI/GPU usage and the case-optimization Fypp guards all check out. The silent-failure audit confirmed no path exists by which a validated non-Newtonian case silently receives Newtonian or zero viscosity, and the mu_max CFL bound is airtight against the clamp logic.

The remaining findings are at the validator boundary and in docs — all cheap to fix, none in the numerical plumbing.

Important (should fix before merge)

1. Leftover fabricated citation — examples/2D_lid_driven_cavity_nn/case.py:9-14. The bogus "Li et al. (2015), Phys. Rev. E 85, 016710" reference was removed from the README in 5790b6c7 but survives in the case.py docstring (PRE 85, 016710 is a 2012 thermal-LBM paper by different authors). Delete and match the README's qualitative framing.
2. Self-contradictory cavity README — examples/2D_lid_driven_cavity_nn/README.md:7 says shear-thinning gives "a blunter primary vortex and reduced near-wall velocity gradients"; line 35 says "vortex center shifted toward the moving lid, with stronger near-wall velocity gradients." The line-35 version is the physically correct one (shear-thinning thins the lid boundary layer); the intro (introduced by the round-1 fix) should be rewritten to match.
3. check_non_newtonian silently skips all validation when num_fluids is unset — toolchain/mfc/case_validator.py:946-948 early-returns on num_fluids is None, while the sibling check_viscosity defaults it to 1. A model_eqns = 1 case (the one model where num_fluids may be omitted) with non_newtonian = 'T' then validates cleanly and runs fully inviscid — including skipping the very guard ("requires model_eqns 2 or 3") that would have rejected it. Default to 1 like check_viscosity.
4. Stray HB parameters are silently ignored when non_newtonian is not set — the per-fluid loop continues past fluids without the flag, so a user who sets K/nn/tau0/mu_max but drops the non_newtonian = 'T' line gets a clean Newtonian run with no warning. The validator's own convention elsewhere is to reject this pattern ("sigma is set but surface_tension is not enabled", num_ibs/Bx0 likewise). Add the matching prohibit.

Suggestions

• Delete the unreachable hb_mu_bulk branch in s_compute_mixture_inv_re (m_hb_function.fpp:82-85). Three agents converged on this independently: an NN fluid can never appear in Re_idx(2,:) (the validator forbids both mu_bulk and Re(2) for NN fluids), so the branch is dead by construction, reads as if bulk works, and would still be a silent no-op if the prohibit were ever relaxed.
• Two cheap high-value regression cases for the fast suite: (a) a mixed two-fluid case (one NN + one Newtonian) — the alpha-weighted mixed accumulation in s_compute_mixture_inv_re is the only core constitutive plumbing with zero CI execution, and it sits on the Re_idx indirection most exposed to future index refactors; (b) a tau0 > 0 fast case — the yield/Papanastasiou branch is currently CI-pinned only at the example tolerance (1e-3 ≈ 25% of that case's peak momentum), so a modest mis-scaling of the regularization would pass.
• Don't count the cavity golden as a physics anchor. Its 50-step truncation leaves every NN-sensitive field value below the 1e-3 example tolerance — it is crash/smoke coverage only. (The two Poiseuille example goldens are genuine anchors: they run full t_stop toward a contracting steady state, and the fast-suite nn=0.5 vs nn=1.5 goldens differ at ~4500x the 1e-12 tolerance floor — a swapped exponent or dropped term fails loudly.)
• Document the new validator constraints in case.md: the Re(2)-forbidden rule, the positivity requirements, model_eqns ∈ {2,3}, and the igr incompatibility are enforced but undocumented.
• NN face alphas bypass mpp_lim treatment (m_riemann_solvers.fpp:~4186/~4398): the NN path averages raw cell-centered volume fractions where the Newtonian path uses the reconstructed/limited face alphas; with mpp_lim = T and a sharp interface, alpha undershoot can locally collapse the face viscosity to ~0 via the sgm_eps guard. Interface-localized and bounded — worth a clamp or a documented-limitation note.
• Minor docs: carry the "at unit shear rate" qualifier for the cavity's Re_eff = 1/K = 100 into the README (conventional lid-based Re ≈ 50); drop the in-PR development-history comments ("the previous mu_max = 10 ...") from the thickening README and HB case.py; the negative-test framework asserts "any error" rather than the targeted message, so the riemann_solver=5 NN test may pass via an unrelated constraint.

Strengths

• Newtonian invariance is structurally enforced (single any_non_newtonian gate) and was re-verified empirically — zero changes to the 63 pre-existing viscous goldens.
• Two-tier test design judged adequate: tight-tolerance (1e-12) synthetic anchors for the constitutive formula + loose-tolerance full-physics example anchors, covering both weno_Re_flux viscous-flux code paths, the adaptive-dt bound, and the yield branch.
• Since all 7 NN tests live in the standard suite, the GPU CI lanes now exercise the NN device routines automatically — the manual V100 verification is no longer the only line of defense.
• Validation provenance is strong: four closed-form analytic gates (power-law thinning 0.68%, thickening 1.46%, Bingham plug flow 2.5%, general HB 3.8% relative L2, with plug positions matching τ₀/(ρg)).

Review generated by parallel specialist agents; each Important finding above was manually audited against the source before posting.

[sbryngelson]

All findings from the round-2 review are now addressed, in four commits:

• 83b10784 (validator): check_non_newtonian defaults num_fluids to 1 instead of silently skipping (Important Develop #3); HB parameters set on a fluid without non_newtonian = 'T' are now rejected, matching the validator's sigma/num_ibs convention (Important fix N-D problem with newer compilers #4). The riemann_solver=5 negative test now isolates the NN guard, and two new negative tests cover the added checks (30/30 pass).
• 0f86fe58 (Fortran cleanup): the unreachable hb_mu_bulk branch and its device-array plumbing are deleted (the mu_bulk parameter remains reserved and validator-rejected); NN interface volume fractions are clamped to [0,1] at both flux sites, closing the mpp_lim under/overshoot scenario. Verified golden-neutral: Viscous 63/63 and all NN cases pass with zero changes to existing goldens.
• 747a0ab2 (tests): the two suggested fast-suite cases added with goldens — a tau0 > 0 case pinning the yield/Papanastasiou branch at the 1e-12 tolerance (1D+2D), and a mixed two-fluid NN+Newtonian case exercising the alpha-weighted mixture loop that previously had zero CI execution. NN regression block is now 7 cases.
• b4a3acdb (docs): the leftover fabricated citation in the cavity case.py docstring is removed (Important Error on READE.md #1); the README contradiction is resolved in favor of the physically correct statement — stronger near-wall gradients (Important Update m_global_parameters.f90 #2); Re_eff = 100 is qualified as the unit-shear-rate value (lid-based Re ≈ 50) and the truncated CI run is labeled smoke coverage; case.md now documents all enforced validator constraints (Re(2) rejection, positivity, model_eqns ∈ {2,3}, igr, flag-required-for-HB-params); in-PR development-history notes dropped.

Deliberately not addressed (disclosed scope limits in the PR description): cylindrical/IBM/3D NN paths have no independent benchmark beyond Newtonian invariance + shared-code reuse, and the compare_analytic.py scripts remain non-asserting documentation of the one-time analytic validation.