Add Herschel-Bulkley non-Newtonian viscosity model#1298
Add Herschel-Bulkley non-Newtonian viscosity model#1298jasontruong2707 wants to merge 20 commits intoMFlowCode:masterfrom
Conversation
Review Summary by QodoAdd Herschel-Bulkley non-Newtonian viscosity model with dynamic Reynolds number computation
WalkthroughsDescription• Implements Herschel-Bulkley non-Newtonian viscosity model with Papanastasiou regularization for power-law, Bingham, and HB fluids • Replaces static viscous Reynolds arrays (Res_viscous, Res_gs, Res_pr) with dynamic computation via new m_re_visc module • Adds per-fluid non-Newtonian parameters: yield stress (tau0), consistency index (K), flow behavior index (nn), viscosity bounds (mu_max, mu_min, mu_bulk), and regularization parameter (hb_m) • Integrates non-Newtonian viscosity computation into time-stepping CFL calculations, Riemann solvers, and data output stability criteria • Adds validation for non-Newtonian fluid parameters in input checker • Extends MPI communication across pre-processing, simulation, and post-processing modules to broadcast non-Newtonian properties • Supports non-Newtonian viscosity in immersed boundary method calculations • Includes two validation cases: 2D Poiseuille flow and 2D lid-driven cavity with non-Newtonian fluids • Adds non-Newtonian parameter definitions to toolchain configuration • Maintains backward compatibility with Newtonian flows (non-Newtonian disabled by default) Diagramflowchart LR
A["Input Parameters<br/>non_newtonian, tau0, K, nn"] --> B["m_hb_function<br/>Herschel-Bulkley<br/>Viscosity Model"]
B --> C["m_re_visc<br/>Dynamic Re Computation"]
C --> D["Time Stepping<br/>CFL Calculation"]
C --> E["Riemann Solvers<br/>Viscous Fluxes"]
C --> F["Data Output<br/>Stability Criteria"]
G["Velocity Gradients"] --> C
H["Strain Rate Tensor"] --> B
File Changes1. src/simulation/m_time_steppers.fpp
|
Code Review by Qodo
1. non_newtonian lacks case validation
|
![qodo-code-review[bot]](https://avatars.githubusercontent.com/in/484649?s=60&v=4){kind=small}
| _r(f"{px}non_newtonian", LOG, {"viscosity"}) | ||
| for a in ["tau0", "K", "nn", "mu_max", "mu_min", "mu_bulk", "hb_m"]: | ||
| _r(f"{px}{a}", REAL, {"viscosity"}) |
This comment was marked as outdated.
This comment was marked as outdated.
Sorry, something went wrong.
| if (any_non_newtonian) then | ||
| if (igr) then | ||
| call s_compute_re_visc(q_cons_ts(1)%vf, alpha, j, k, l, Re_visc_per_phase) | ||
| else | ||
| call s_compute_re_visc(q_prim_vf, alpha, j, k, l, Re_visc_per_phase) | ||
| end if | ||
| call s_compute_mixture_re(alpha, Re_visc_per_phase, Re) |
This comment was marked as outdated.
This comment was marked as outdated.
Sorry, something went wrong.
|
Important Review skippedDraft detected. Please check the settings in the CodeRabbit UI or the ⚙️ Run configurationConfiguration used: defaults Review profile: CHILL Plan: Pro Run ID: You can disable this status message by setting the Use the checkbox below for a quick retry:
📝 WalkthroughWalkthroughAdds Herschel–Bulkley non‑Newtonian support: eight new fields to the derived type physical_parameters (non_newtonian, tau0, K, nn, mu_min, mu_max, mu_bulk, hb_m); new modules m_hb_function and m_re_visc for HB viscosity and viscous/Re computations; input validation via s_check_inputs_non_newtonian; MPI broadcasts and GPU updates extended to the new fields; time‑stepper dt computation updated to use per‑phase viscous terms; pressure‑relaxation internals simplified; two example case scripts added; and CI workflow scripts reference shared frontier scripts. 🚥 Pre-merge checks | ✅ 3✅ Passed checks (3 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. ✨ Finishing Touches🧪 Generate unit tests (beta)
Tip 💬 Introducing Slack Agent: The best way for teams to turn conversations into code.Slack Agent is built on CodeRabbit's deep understanding of your code, so your team can collaborate across the entire SDLC without losing context.
Built for teams:
One agent for your entire SDLC. Right inside Slack. Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
Claude Code ReviewHead SHA: ab1a4c8 Key files: Summary
FindingsCritical / Correctness1. Hard-coded shear-rate clamp is unphysical and problem-specific shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)This clamps 2. Riemann-solver index mapping for NORM_DIR=2 and NORM_DIR=3 is incorrect The Riemann solver loops over rotated index triplets call s_compute_re_visc(q_prim_vf, alpha_L, k, j, l, ...)but inside 3. Significant4. Performance regression on Newtonian-only cases 5. 6. IBM non-Newtonian viscosity uses dynamic_viscosities(fluid_idx) = f_compute_hb_viscosity(..., 1._wp, ...)The IBM force/torque integration evaluates viscosity at a fixed shear rate of 1 s⁻¹ rather than the actual local shear rate at each IB ghost-cell. This will give incorrect viscous forces on immersed boundaries in non-Newtonian flows with spatially varying Minor7. File header says !! @file m_hb_function.f90 ! should be .fpp
!! @file m_re_visc.f90 ! should be .fpp8. 9. Missing regression test for non-Newtonian viscous flows The IBM torque bug fix, the |
Claude Code ReviewHead SHA: 4d042f5 Key files:
Summary:
Findings1. Shear-rate lower clamp is physically wrong —
|
Code Review — Non-Newtonian ViscosityThorough investigation of all findings from the automated review. Each item was verified against the source code. Critical / Correctness1. Index mapping bug in Riemann solver for NORM_DIR=2,3
For NORM_DIR=2 (y-normal faces), the Riemann solver calls: call s_compute_re_visc(q_prim_vf, alpha_L, k, j, l, Re_visc_per_phase_L)Inside Impact: Wrong shear rate → wrong viscosity at y-normal and z-normal Riemann interfaces in 2D/3D non-Newtonian simulations. Most severe on anisotropic grids. Newtonian fluids are unaffected (early exit before gradient computation). Affected locations in Fix options:
2. GPU coherence:
GPU kernels in Impact: Silent wrong results on all GPU builds when Fix: Add Significant3. Hard-coded shear-rate clamp is not parameterizable
shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)The bounds Recommendation: Either make these user-configurable ( 4. Missing pre-computed gradients in Riemann solver calls All calls to Impact: Performance regression (redundant gradient computation) and slight inconsistency with the viscous flux gradients. Not a correctness bug. 5. IBM non-Newtonian viscosity uses fixed
dynamic_viscosities(fluid_idx) = f_compute_hb_viscosity( &
..., 1._wp, ...) ! Fixed shear rate = 1The IBM force/torque calculation evaluates viscosity at a constant shear rate of 1 s⁻¹ rather than the local shear rate at each IB ghost cell. The comment "compute reference viscosity at gdot=1" acknowledges this but doesn't explain the physical reasoning. Impact: Incorrect viscous forces on immersed boundaries in flows with spatially varying shear rate. Low severity if IBM + non-Newtonian is not a primary use case. 6. Newtonian performance regression ✅ No issue Verified: Minor (fixed in latest push)7. File headers say
8. Boundary shear rate asymmetry
9. Missing regression test ✅ Fixed Added two non-Newtonian test cases to
Both run in the single-fluid viscous test block for each dimension. Summary
Items 1 and 2 should be addressed before merge — they produce silently wrong results in 2D/3D and GPU builds respectively. |
Claude Code ReviewHead SHA: Files changed: 28 Key files:
Summary
Findings🔴 CRITICAL — Missing
|
Claude Code ReviewHead SHA: Summary
Findings🔴 Critical — Missing namelist bindings in
|
Claude Code ReviewHead SHA: 9c15db4 Key files
Summary
Findings🔴 Critical — Missing namelist bindings in
|
Claude Code ReviewHead SHA: 4635109 Summary:
Finding 1 — CRITICAL: Missing
|
sbryngelson
force-pushed
the
feature/non-newtonian-viscosity
branch
from
5647466 to
d2fc161
Compare
Claude Code ReviewIncremental review from: `1ab3e1b86357970ac502583f37157b7a177ab44b` New findings since last Claude review: Finding — HIGH: `weno_Re_flux requires viscous` validation silently droppedFile: `toolchain/mfc/case_validator.py` While fixing the previous review's Finding 8 (duplicate `weno_Re_flux` check), the commit that added the standalone `check_non_newtonian` function also removed the only surviving copy of the constraint from `check_viscosity`: ```python These 3 lines were removed from check_viscosity:weno_Re_flux = self.get("weno_Re_flux", "F") == "T" There is no Fortran runtime check for this constraint. A case with `weno_Re_flux=T` and `viscous=F` will now pass `./mfc.sh validate` without error, but will silently produce incorrect results at runtime (or crash). Fix: restore those 3 lines to `check_viscosity`, after the per-fluid loop. |
sbryngelson
force-pushed
the
feature/non-newtonian-viscosity
branch
from
d2fc161 to
b8d0fe0
Compare
|
<!-- claude-review: thread=primary; reviewed_sha=1d178c437ed1cb04046e1197fd1b81be70a365d7; mode=incremental --> |
|
Claude Code Review Incremental review from: cc6373c New findings since last Claude review:
|
sbryngelson
force-pushed
the
feature/non-newtonian-viscosity
branch
6 times, most recently
from
e01e168 to
b7a3b4d
Compare
sbryngelson
force-pushed
the
feature/non-newtonian-viscosity
branch
from
b7a3b4d to
1be6ef3
Compare
… to c_fast from old merge)
… merge artifact from PR)
Resolve conflicts: - m_data_output.fpp: keep master's ib_state_unit and ccfl removal - m_ibm.fpp: keep PR's non-Newtonian branch with master's Re(1) > 0 fix
- module_categories.json: keep m_hb_function, m_re_visc and add m_collisions - m_global_parameters.fpp: use eqn_idx%adv/E naming (master) - m_time_steppers.fpp: keep IB comment from master - cases.py: keep all 4 example skip entries (both branches) - m_ibm.fpp: remove old E_idx/momxb copyin; use eqn_idx%E; drop duplicate s_cross_product (now in m_helper.fpp) - m_data_output.fpp: take master s_write_ib_state_file subroutines; shorten docstrings; add q_T_sf param - m_riemann_solvers.fpp: keep NORM_DIR-aware s_compute_re_visc calls (non-Newtonian fix)
Claude Code ReviewPR #1298 — Add Herschel-Bulkley non-Newtonian viscosity model Critical Issues1. Hard-coded shear rate clamp overrides user viscosity bounds shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)The lower bound of 2. IBM force computation uses a fixed non-physical shear rate dynamic_viscosities(fluid_idx) = f_compute_hb_viscosity(..., 1._wp, ...)Viscosity at immersed boundary cells is evaluated at Important Issues3. @:PROHIBIT(fluid_pp(i)%mu_max < dflt_real .and. fluid_pp(i)%mu_max <= fluid_pp(i)%mu_min, ...
4. Mixed Newtonian/non-Newtonian fluids produce dimensionally inconsistent Re When
5.
Suggestions
Strengths
|
…D_xz/D_yz boundary fallback, weno_Re_flux validation
…ectives GPU_LOOP macros expand to OpenMP directives, which are not permitted inside pure subroutines (nvfortran S-0155). Remove pure from s_compute_re_visc and s_compute_mixture_re, which each contain multiple GPU_LOOP macros. The pure qualifier was semantically correct but syntactically invalid for older nvfortran (24.3, 24.5, 24.11) in OpenMP target offload builds.
The shear rate clamp removal (fix MFlowCode#2) slightly changes viscosity in low-shear-rate boundary cells, producing different output that the old golden did not capture. This regenerates the golden using the corrected code with --no-gpu on Phoenix (matching the gfortran CPU builds used by GitHub CI).
2 participants
Summary
Implements Herschel-Bulkley (HB) non-Newtonian viscosity with Papanastasiou regularization, covering power-law, Bingham plastic, and general HB fluids. The effective viscosity at each cell is computed dynamically from the local strain-rate tensor rather than a pre-cached Reynolds number, so viscosity varies spatially in every time step.
The implementation is backward-compatible: all new parameters default to Newtonian behavior (
non_newtonian = F), and the performance-criticalany_non_newtonianflag gates the shear-rate computation so purely Newtonian runs are unaffected.Physics model
The Herschel-Bulkley model with Papanastasiou regularization:
where
γ̇ = √(2 Dᵢⱼ Dᵢⱼ)is the second invariant of the strain-rate tensor. The Papanastasiou term avoids the singularity atγ̇ → 0analytically, removing the need for any hard clamp on the shear rate.Special cases:
τ₀ = 0, n = 1: Newtonian withμ = Kτ₀ = 0, n ≠ 1: Power-law fluidτ₀ > 0, n = 1: Bingham plasticτ₀ > 0, n ≠ 1: Herschel-BulkleyNondimensionalization: All HB parameters (
K,τ₀,mu_min,mu_max,mu_bulk) must be provided in MFC's nondimensional units, i.e., scaled byρ_ref · U_ref · L_ref. At any shear rate,1/μ_effin these units equals the local effective Reynolds number, consistent with howfluid_pp(i)%Re(1)is used for Newtonian fluids. The mixture formulaRe_mix = 1 / Σᵢ(αᵢ / Re_eff_i)is the correct harmonic-mean viscosity mixing rule and applies uniformly to Newtonian, non-Newtonian, and mixed fluid systems.New parameters (
fluid_pp(i)%...)non_newtonianFKdflt_realnndflt_realtau0dflt_realhb_mdflt_realmu_mindflt_realmu_maxdflt_realmu_bulkdflt_realAll parameters are registered in all four required locations:
params/definitions.py,m_global_parameters.fpp(all three targets),m_start_up.fpp(all three targets), andcase_validator.py.New files
src/simulation/m_hb_function.fppTwo pure GPU_ROUTINE functions:
f_compute_hb_viscosity(tau0, K, nn, mu_min, mu_max, shear_rate, hb_m)— evaluates the Papanastasiou-regularized HB formula, clamps to[mu_min, mu_max]when those are explicitly set (checked againstdflt_realsentinel, not> 0)f_compute_shear_rate_from_components(D_xx, D_yy, D_zz, D_xy, D_xz, D_yz)— computesγ̇ = √(2 Dᵢⱼ Dᵢⱼ)from strain-rate tensor components; no artificial clamp (Papanastasiou handles near-zero analytically)src/simulation/m_re_visc.fppThree pure GPU_ROUTINE subroutines:
s_compute_velocity_gradients_at_cell(q_prim_vf, j, k, l, D_xx, ...)— finite-difference strain-rate tensor at a single cell. Uses 2nd-order central differences in interior; in boundary cells, each of the six tensor components (D_xx,D_yy,D_zz,D_xy,D_xz,D_yz) independently falls back to whatever terms have valid neighbors, accumulating partial contributions rather than zeroing the whole component. This matches the D_xy boundary treatment throughout.s_compute_re_visc(q_prim_vf, alpha, j, k, l, Re_visc_per_phase, [grad_x_vf, grad_y_vf, grad_z_vf])— fills a(num_fluids, 2)array with per-phase1/μ_eff(shear, bulk) for the local cell. Optional pre-computed gradient arrays are used when available (viscous flux path); otherwise falls back tos_compute_velocity_gradients_at_cell. The Newtonian path (whenany_non_newtonian = F) early-exits and readsfluid_pp(q)%Re(i)directly, preserving pre-PR performance.s_compute_mixture_re(alpha, Re_visc_per_phase, Re_mix)— harmonic-mean mixture:Re_mix(i) = 1 / Σ_q(α_q / Re_visc_per_phase(q,i)), skipping entries atdflt_real(inactive phases)examples/2D_lid_driven_cavity_nn/Shear-thinning (
n = 0.5) and shear-thickening (n = 1.5) lid-driven cavity cases validated against Li et al. (2015), Re = 500.examples/2D_poiseuille_nn/Power-law Poiseuille flow with analytical velocity profile comparison.
Changes to existing modules
src/common/m_derived_types.fppAdds eight new fields to
physical_parameters:non_newtonian(logical),tau0,K,nn,mu_min,mu_max,mu_bulk,hb_m(allreal(wp)).src/simulation/m_global_parameters.fppany_non_newtonianlogical flag (set at startup, GPU_DECLARE'd)fluid_pparray so HB parameters are resident on devicefluid_ppafter MPI broadcast so device copy is coherentsrc/simulation/m_viscous.fppRemoves the pre-allocated
Res_viscous(2, Re_size_max)array and its GPU data region. Each GPU parallel loop now callss_compute_re_visc+s_compute_mixture_reper cell, using pre-computedgrad_{x,y,z}_vfto avoid redundant finite differences. Theany_non_newtonianflag gates the shear-rate computation so Newtonian runs are unaffected.src/simulation/m_riemann_solvers.fppRemoves
Res_gscached array. All three Riemann solvers (HLL, HLLC, HLLD) now calls_compute_re_visc+s_compute_mixture_reper face to compute the local viscous Re. The optional gradient arguments are not passed here (gradients are recomputed from cell-centered primitives via finite differences), which is consistent with the existing first-order viscous flux reconstruction in the Riemann solver path.src/simulation/m_time_steppers.fpps_compute_dtconverts to primitive variables whenany_non_newtonianis true before the CFL loop, then callss_compute_re_viscto get the local effective viscosity for the viscous CFL constraint.src/simulation/m_pressure_relaxation.fppRemoves
Res_prcached array; init/finalize are now empty stubs.src/simulation/m_ibm.fppPreviously approximated IBM viscous forces using
μ_eff(γ̇ = 1)as a spatially uniform reference viscosity — physically wrong for non-Newtonian fluids where viscosity varies with local shear rate. Now:1/Re(1)for Newtonian fluids (correct, since their viscosity is constant)any_non_newtonian, callss_compute_velocity_gradients_at_cellat the current cell(i, j, k)to get the actual local strain-rate tensor, thenf_compute_shear_rate_from_components, thenf_compute_hb_viscosityper non-Newtonian fluiddynamic_viscosityis accumulated with the correct per-cell, per-fluid viscosity before computing the viscous stress divergencesrc/simulation/m_checker.fppAdds
s_check_inputs_non_newtonianwith@:PROHIBITguards:K > 0,nn > 0,tau0 ≥ 0,mu_min ≥ 0,hb_m > 0. Themu_maxconstraint usesf_is_default(fluid_pp(i)%mu_max)(fromm_helper_basic) to correctly detect whether the user has set the parameter, rather than the incorrectmu_max < dflt_realsentinel comparison (which is never true sincedflt_real = -1e6).src/simulation/m_data_output.fppAdds output of per-fluid HB viscosity fields when
any_non_newtonianis true.toolchain/mfc/params/definitions.pyRegisters all eight new
fluid_ppparameters with types, defaults, tags, and documentation strings.toolchain/mfc/case_validator.pycheck_non_newtonian(): Python-side validation of HB parameter ranges (K > 0, nn > 0, tau0 ≥ 0, mu_min ≥ 0, mu_max > mu_min when set, hb_m > 0)weno_Re_flux requires viscousconstraint tocheck_viscosity()(had been inadvertently removed during refactoring)src/{pre_process,post_process}/m_global_parameters.fppandm_start_up.fppNew parameters declared and bound in namelists for all three executables.
Regression tests
Four new test UUIDs in
toolchain/mfc/test/cases.py, each with golden files:0CD345C7non_newtonian=T, K=0.1, nn=0.5, tau0=0, hb_m=1001E64A4D9non_newtonian=T, K=0.1, nn=1.5, tau0=0, hb_m=10030713C1Cnon_newtonian=T, K=0.01, nn=1, tau0=0.1, hb_m=10045DA4729non_newtonian=T, K=0.1, nn=0.7, tau0=0.05, mu_min=1e-4, mu_max=10Test plan
non_newtonian=Fby default)