Skip to content

Add Herschel-Bulkley non-Newtonian viscosity model#1298

Draft
jasontruong2707 wants to merge 11 commits intoMFlowCode:masterfrom
jasontruong2707:feature/non-newtonian-viscosity
Draft

Add Herschel-Bulkley non-Newtonian viscosity model#1298
jasontruong2707 wants to merge 11 commits intoMFlowCode:masterfrom
jasontruong2707:feature/non-newtonian-viscosity

Conversation

@jasontruong2707
Copy link

Summary

  • Implements Herschel-Bulkley (HB) viscosity model with Papanastasiou regularization
  • Supports power-law, Bingham, and HB fluids via per-fluid parameters
  • Replaces static Res_viscous array with dynamic s_compute_re_visc for both Newtonian and non-Newtonian flows
  • Validated against Li et al. (2015) lid-driven cavity benchmarks (Re=500, n=0.5 and n=1.5)

New files

  • src/simulation/m_hb_function.fpp — Core HB viscosity computation
  • src/simulation/m_re_visc.fpp — Dynamic Reynolds number module
  • examples/2D_lid_driven_cavity_nn/ — Validation case
  • examples/2D_poiseuille_nn/ — Validation case

Test plan

  • Builds with gfortran (CPU)
  • 2D lid-driven cavity non-Newtonian cases run without errors
  • Results match Li et al. (2015) benchmarks
  • Existing viscous test cases unaffected (non_newtonian=F by default)

@qodo-code-review
Copy link
Contributor

Review Summary by Qodo

Add Herschel-Bulkley non-Newtonian viscosity model with dynamic Reynolds number computation

✨ Enhancement 🧪 Tests

Grey Divider

Walkthroughs

Description
• 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)
Diagram
flowchart 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
Loading

Grey Divider

File Changes

1. src/simulation/m_time_steppers.fpp ✨ Enhancement +1075/-1061

Integrate non-Newtonian viscosity into time-stepping CFL calculations

• Added use m_re_visc module import for non-Newtonian viscosity computations
• Modified s_compute_dt() subroutine to compute per-phase Reynolds numbers for non-Newtonian
 fluids using s_compute_re_visc() and s_compute_mixture_re()
• Added local variables Re_visc_per_phase to track per-phase viscous Reynolds numbers in the CFL
 computation

src/simulation/m_time_steppers.fpp


2. src/simulation/m_viscous.fpp ✨ Enhancement +1578/-1622

Replace static viscous Reynolds array with dynamic computation

• Removed static Res_viscous array allocation and initialization from
 s_initialize_viscous_module()
• Replaced four instances of inline Reynolds number computation with calls to s_compute_re_visc()
 and s_compute_mixture_re() functions
• Added local variable Re_visc_nn to store per-phase viscous Reynolds numbers
• Added comment explaining that dynamic viscosity computation now handles both Newtonian and
 non-Newtonian cases

src/simulation/m_viscous.fpp


3. src/simulation/m_global_parameters.fpp ✨ Enhancement +20/-2

Add non-Newtonian fluid parameters and detection logic

• Added any_non_newtonian logical flag to track if any fluid is non-Newtonian
• Initialized non-Newtonian fluid parameters (non_newtonian, tau0, K, nn, mu_max,
 mu_min, mu_bulk, hb_m) in fluid_pp structure
• Added detection logic to set any_non_newtonian flag based on fluid properties
• Updated GPU declarations to include any_non_newtonian flag

src/simulation/m_global_parameters.fpp


View more (20)
4. src/simulation/m_data_output.fpp ✨ Enhancement +1984/-1974

Integrate non-Newtonian viscosity into data output stability criteria

• Added use m_re_visc module import for non-Newtonian viscosity computations
• Added Re_visc_per_phase local variable to compute per-phase Reynolds numbers
• Integrated non-Newtonian Reynolds number computation in stability criteria calculation using
 s_compute_re_visc() and s_compute_mixture_re()

src/simulation/m_data_output.fpp


5. src/common/m_derived_types.fpp ✨ Enhancement +8/-0

Extend physical parameters type with non-Newtonian fluid properties

• Extended physical_parameters derived type with eight new fields for non-Newtonian fluid modeling
• Added Herschel-Bulkley model parameters: tau0 (yield stress), K (consistency index), nn
 (flow behavior index)
• Added viscosity bounds: mu_max, mu_min, and mu_bulk for non-Newtonian fluids
• Added hb_m parameter for Papanastasiou regularization

src/common/m_derived_types.fpp


6. .github/workflows/frontier_amd/bench.sh ⚙️ Configuration changes +0/-1

Add Frontier AMD benchmark script symlink

• Created symbolic link pointing to parent directory benchmark script

.github/workflows/frontier_amd/bench.sh


7. src/simulation/m_riemann_solvers.fpp ✨ Enhancement +5272/-5257

Replace static viscosity array with dynamic non-Newtonian computation

• Removed static Res_gs array and replaced with dynamic s_compute_re_visc function calls
• Added per-phase Reynolds number arrays Re_visc_per_phase_L and Re_visc_per_phase_R for
 non-Newtonian support
• Updated GPU parallel loop declarations to include new per-phase Re arrays
• Replaced inline Re computation loops with calls to s_compute_mixture_re for both Newtonian and
 non-Newtonian cases

src/simulation/m_riemann_solvers.fpp


8. src/simulation/m_pressure_relaxation.fpp ✨ Enhancement +279/-308

Remove static viscosity array from pressure relaxation module

• Removed static Res_pr array allocation and initialization
• Simplified s_initialize_pressure_relaxation_module and s_finalize_pressure_relaxation_module
 to no-ops
• Removed viscosity computation from s_correct_internal_energies (now handled dynamically via
 m_re_visc)

src/simulation/m_pressure_relaxation.fpp


9. src/simulation/m_re_visc.fpp ✨ Enhancement +346/-0

New dynamic viscosity module for Newtonian and non-Newtonian fluids

• New module implementing dynamic viscosity computation for both Newtonian and non-Newtonian fluids
• s_compute_velocity_gradients_at_cell computes strain rate tensor components using finite
 differences
• s_compute_re_visc returns per-phase Re_visc = 1/mu using Herschel-Bulkley model for
 non-Newtonian fluids
• s_compute_mixture_re combines per-phase values with volume fractions to compute mixture Reynolds
 number

src/simulation/m_re_visc.fpp


10. src/simulation/m_hb_function.fpp ✨ Enhancement +77/-0

New Herschel-Bulkley non-Newtonian viscosity model implementation

• New module implementing Herschel-Bulkley viscosity model with Papanastasiou regularization
• f_compute_hb_viscosity computes viscosity from yield stress, consistency index, and flow
 behavior index
• f_compute_shear_rate_from_components calculates shear rate magnitude from strain rate tensor
 components

src/simulation/m_hb_function.fpp


11. src/simulation/m_checker.fpp Error handling +27/-0

Add validation for non-Newtonian fluid parameters

• Added s_check_inputs_non_newtonian subroutine to validate non-Newtonian fluid parameters
• Validates that viscosity is enabled, consistency index K > 0, flow behavior index nn > 0, yield
 stress tau0 >= 0
• Validates viscosity bounds and Papanastasiou regularization parameter hb_m > 0

src/simulation/m_checker.fpp


12. src/post_process/m_mpi_proxy.fpp ✨ Enhancement +4/-0

Broadcast non-Newtonian fluid parameters in post-processing MPI proxy

• Added MPI broadcast calls for non-Newtonian fluid properties: non_newtonian, tau0, K, nn,
 mu_max, mu_min, mu_bulk, hb_m

src/post_process/m_mpi_proxy.fpp


13. src/pre_process/m_mpi_proxy.fpp ✨ Enhancement +5/-0

Broadcast non-Newtonian fluid parameters in pre-processing MPI proxy

• Added MPI broadcast calls for non-Newtonian fluid properties: non_newtonian, tau0, K, nn,
 mu_max, mu_min, mu_bulk, hb_m

src/pre_process/m_mpi_proxy.fpp


14. src/simulation/m_mpi_proxy.fpp ✨ Enhancement +4/-0

Broadcast non-Newtonian fluid parameters in simulation MPI proxy

• Added MPI broadcast calls for non-Newtonian fluid properties: non_newtonian, tau0, K, nn,
 mu_max, mu_min, mu_bulk, hb_m

src/simulation/m_mpi_proxy.fpp


15. src/simulation/m_ibm.fpp ✨ Enhancement +9/-1

Support non-Newtonian viscosity in immersed boundary method

• Added conditional logic to use mu_max or K as reference viscosity for non-Newtonian fluids in
 IBM calculations
• Falls back to consistency index K if mu_max is not set

src/simulation/m_ibm.fpp


16. src/post_process/m_global_parameters.fpp ✨ Enhancement +8/-0

Initialize non-Newtonian fluid parameters in post-processing

• Initialize non-Newtonian fluid properties: non_newtonian, tau0, K, nn, mu_max, mu_min,
 mu_bulk, hb_m
• Default values: non_newtonian=.false., nn=1.0, hb_m=1000.0, others zero or default

src/post_process/m_global_parameters.fpp


17. src/pre_process/m_global_parameters.fpp ✨ Enhancement +8/-0

Initialize non-Newtonian fluid parameters in pre-processing

• Initialize non-Newtonian fluid properties: non_newtonian, tau0, K, nn, mu_max, mu_min,
 mu_bulk, hb_m
• Default values: non_newtonian=.false., nn=1.0, hb_m=1000.0, others zero or default

src/pre_process/m_global_parameters.fpp


18. examples/2D_poiseuille_nn/case.py 🧪 Tests +158/-0

Add 2D Poiseuille non-Newtonian validation case

• New validation case for 2D Poiseuille flow with Herschel-Bulkley non-Newtonian fluid
• Pressure-driven channel flow with periodic BCs in streamwise direction and no-slip walls
• Configurable HB parameters: yield stress, consistency index, flow behavior index, regularization
 parameter

examples/2D_poiseuille_nn/case.py


19. examples/2D_lid_driven_cavity_nn/case.py 🧪 Tests +116/-0

Add 2D lid-driven cavity non-Newtonian validation case

• New validation case for 2D lid-driven cavity with Herschel-Bulkley non-Newtonian fluid
• Re=5000, shear-thickening fluid (nn=1.5) with mesh stretching near walls
• Demonstrates non-Newtonian behavior in complex flow with recirculation zones

examples/2D_lid_driven_cavity_nn/case.py


20. toolchain/mfc/params/definitions.py ⚙️ Configuration changes +3/-0

Add non-Newtonian fluid parameter definitions to toolchain

• Added parameter definitions for non-Newtonian fluid properties: non_newtonian (logical), tau0,
 K, nn, mu_max, mu_min, mu_bulk, hb_m (real)
• All parameters registered under viscosity category

toolchain/mfc/params/definitions.py


21. .github/workflows/frontier_amd/submit.sh ⚙️ Configuration changes +0/-1

Add AMD frontier workflow submit script link

• Symbolic link to parent frontier submit script

.github/workflows/frontier_amd/submit.sh


22. .github/workflows/frontier_amd/build.sh ⚙️ Configuration changes +0/-1

Add AMD frontier workflow build script link

• Symbolic link to parent frontier build script

.github/workflows/frontier_amd/build.sh


23. .github/workflows/frontier_amd/test.sh ⚙️ Configuration changes +0/-1

Add AMD frontier workflow test script link

• Symbolic link to parent frontier test script

.github/workflows/frontier_amd/test.sh


Grey Divider

Qodo Logo

@qodo-code-review
Copy link
Contributor

qodo-code-review bot commented Mar 9, 2026

Code Review by Qodo

🐞 Bugs (4) 📘 Rule violations (1) 📎 Requirement gaps (0)

Grey Divider


Action required

1. non_newtonian lacks case validation 📘 Rule violation ⛯ Reliability
Description
New HB parameters (fluid_pp(*)%non_newtonian, tau0, K, nn, mu_*, hb_m) are defined but
not validated in toolchain/mfc/case_validator.py, allowing invalid configurations to pass
toolchain checks. This should be validated in the toolchain since physics constraints are explicitly
enforced in the Fortran input checker.
Code

toolchain/mfc/params/definitions.py[R1061-1063]

+        _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"})
Evidence
PR Compliance ID 6 requires new parameters to be added to toolchain definitions and also to
case_validator.py when constraints apply. The PR adds the non-Newtonian parameters to
definitions.py and enforces constraints in src/simulation/m_checker.fpp, but
case_validator.py's viscosity validation only checks Re(*) and contains no validation for
non_newtonian or HB parameters.

CLAUDE.md
toolchain/mfc/params/definitions.py[1061-1063]
src/simulation/m_checker.fpp[95-116]
toolchain/mfc/case_validator.py[993-1024]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

## Issue description
The PR introduces new non-Newtonian viscosity inputs (`fluid_pp(*)%non_newtonian`, `tau0`, `K`, `nn`, `mu_min`, `mu_max`, `mu_bulk`, `hb_m`) but the toolchain validator does not validate them, so invalid cases can pass toolchain checks and only fail later at runtime.

## Issue Context
Fortran-side validation was added in `src/simulation/m_checker.fpp`, indicating physics constraints apply and should also be enforced in `toolchain/mfc/case_validator.py` per the compliance checklist.

## Fix Focus Areas
- toolchain/mfc/case_validator.py[993-1024]

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


2. IGR uses conservative gradients 🐞 Bug ✓ Correctness
Description
In s_compute_dt, the IGR branch calls s_compute_re_visc with q_cons_ts(1)%vf, but s_compute_re_visc
differentiates q_prim_vf(momxb:), so it treats momenta as velocities and computes an incorrect shear
rate and viscosity. This can yield an incorrect (potentially unstable) time step whenever
any_non_newtonian and igr are enabled.
Code

src/simulation/m_time_steppers.fpp[R748-754]

+                    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)
Evidence
The new dt logic calls s_compute_re_visc with conservative variables when igr is true, but the
viscosity routine explicitly expects primitive variables and computes du/dx etc. from
q_prim_vf(momxb) and q_prim_vf(momxb+1/2), which are velocities in the primitive layout (not
momenta).

src/simulation/m_time_steppers.fpp[729-755]
src/simulation/m_re_visc.fpp[26-32]
src/simulation/m_re_visc.fpp[64-72]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
When `igr` is enabled, `s_compute_dt` calls `s_compute_re_visc` with `q_cons_ts(1)%vf` (conservative state). `s_compute_re_visc` assumes its first argument is primitive variables and differentiates `q_prim_vf(momxb:...)` as velocities, so the IGR path computes wrong shear rate/viscosity and thus wrong dt.

### Issue Context
`m_re_visc` computes velocity gradients directly from the provided field, without converting momenta to velocities.

### Fix Focus Areas
- src/simulation/m_time_steppers.fpp[748-754]
- src/simulation/m_re_visc.fpp[26-90]

### Suggested fix
- Option A (simplest): In `s_compute_dt`, ensure `q_prim_vf` is valid even when `igr` is true (perform conservative→primitive conversion when `any_non_newtonian` is true, or always for dt computation).
- Option B: Add a conservative-aware path in `s_compute_re_visc` (or a new routine) that computes velocities `u=mom/rho` at neighbor cells before differencing, and use that path from the IGR branch.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


3. Frontier AMD scripts mispath 🐞 Bug ⛯ Reliability
Description
The newly added .github/workflows/frontier_amd/*.sh files contain only a relative path like
"../frontier/build.sh", but the workflows execute them via "bash
.github/workflows/frontier_amd/build.sh" from the checkout directory, so the relative path resolves
outside the repo and fails. This breaks Frontier AMD build/test/bench workflow steps.
Code

.github/workflows/frontier_amd/build.sh[1]

+../frontier/build.sh
Evidence
The workflow invokes the Frontier AMD scripts using a repo-root relative path (or after cd pr),
meaning the scripts’ current working directory is not .github/workflows/frontier_amd/. Because the
scripts use ../frontier/... without anchoring to the script directory, the target path is resolved
relative to the current working directory and points outside the repo.

.github/workflows/frontier_amd/build.sh[1-1]
.github/workflows/test.yml[243-266]
.github/workflows/bench.yml[113-117]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
The new `frontier_amd/*.sh` wrappers are one-line scripts like `../frontier/build.sh`. Workflows run them from the checkout directory (repo root or `cd pr`), so `../frontier/...` resolves outside the repo and the job fails.

### Issue Context
GitHub Actions `run:` commands and `bash &lt;path&gt;` do not change CWD to the script’s directory; relative paths inside the script are resolved from the caller’s CWD.

### Fix Focus Areas
- .github/workflows/frontier_amd/build.sh[1-1]
- .github/workflows/frontier_amd/submit.sh[1-1]
- .github/workflows/frontier_amd/test.sh[1-1]
- .github/workflows/frontier_amd/bench.sh[1-1]

### Suggested fix
Replace each one-liner with a small bash wrapper:
```bash
#!/usr/bin/env bash
set -euo pipefail
DIR=&quot;$(cd -- &quot;$(dirname -- &quot;${BASH_SOURCE[0]}&quot;)&quot; &amp;&amp; pwd)&quot;
exec bash &quot;$DIR/../frontier/build.sh&quot; &quot;$@&quot;
```
(and similarly for submit/test/bench), or adjust paths to `.github/workflows/frontier/...` if you prefer not to compute `DIR`.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools



Remediation recommended

4. IBM non-Newtonian mu wrong 🐞 Bug ✓ Correctness
Description
In m_ibm, for non-Newtonian fluids the “reference viscosity” used in immersed-boundary viscous
forces falls back to fluid_pp%K when mu_max is not explicitly set, but K is not the dynamic
viscosity for HB fluids except for the special case tau0=0 and nn=1. This makes IBM viscous
forces/torques inconsistent with the HB viscosity model for yield-stress and general power-law
cases.
Code

src/simulation/m_ibm.fpp[R1016-1023]

+                if (fluid_pp(fluid_idx)%non_newtonian) then
+                    ! Non-Newtonian: use mu_max as reference viscosity for IBM
+                    if (fluid_pp(fluid_idx)%mu_max < dflt_real .and. &
+                        fluid_pp(fluid_idx)%mu_max > sgm_eps) then
+                        dynamic_viscosities(fluid_idx) = fluid_pp(fluid_idx)%mu_max
+                    else
+                        dynamic_viscosities(fluid_idx) = fluid_pp(fluid_idx)%K
+                    end if
Evidence
IBM computes a scalar dynamic_viscosity by volume-averaging dynamic_viscosities(fluid_idx) and
then uses it in viscous stress computation. For non-Newtonian fluids, dynamic_viscosities is set
to mu_max (if set) else to K, but the HB model defines viscosity as a function of shear rate
involving tau0, K, and nn, so using K alone does not represent viscosity in general.

src/simulation/m_ibm.fpp[1014-1027]
src/simulation/m_ibm.fpp[1067-1072]
src/simulation/m_hb_function.fpp[7-45]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
`m_ibm` needs a scalar “dynamic viscosity” per fluid for immersed boundary viscous forces. For non-Newtonian fluids it uses `mu_max` if set, else falls back to `K`, which is not equal to HB viscosity for general (`tau0`, `nn`) settings.

### Issue Context
HB viscosity is shear-rate dependent: `mu(gdot) = tau0*(1-exp(-m*gdot))/gdot + K*gdot^(nn-1)`.

### Fix Focus Areas
- src/simulation/m_ibm.fpp[1014-1027]
- src/simulation/m_hb_function.fpp[23-49]

### Suggested fix
- Prefer computing a reference viscosity via `f_compute_hb_viscosity(...)` at a documented representative shear rate (e.g., `gdot_ref = 1._wp` or a user-provided parameter), then clamp by `mu_min/mu_max`.
- Alternatively, if IBM must use a bound, require `mu_max` to be explicitly set for non-Newtonian fluids when IBM+viscous is enabled (checker-level validation) and remove the `K` fallback.

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


5. Hardcoded shear-rate floor 🐞 Bug ✓ Correctness
Description
The shear-rate magnitude used by the HB model is hard-clamped to be at least 1e-2, so all true shear
rates below 1e-2 produce identical viscosity inputs. This makes viscosity insensitive in low-shear
regions and can override user intent for hb_m/mu_max tuning below that threshold.
Code

src/simulation/m_hb_function.fpp[R68-74]

+        ! 2*D_ij*D_ij = 2*(D_xx^2+D_yy^2+D_zz^2+2*(D_xy^2+D_xz^2+D_yz^2))
+        shear_rate = sqrt(2._wp*(D_xx*D_xx + D_yy*D_yy + D_zz*D_zz + &
+                                 2._wp*(D_xy*D_xy + D_xz*D_xz + D_yz*D_yz)))
+
+        ! Clamp for numerical safety
+        shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)
+
Evidence
HB viscosity depends directly on shear_rate (including division by shear_rate in the yield term),
but the implementation clamps shear_rate to a minimum of 1e-2 before viscosity is computed,
collapsing all lower-shear behavior to a single effective value.

src/simulation/m_hb_function.fpp[42-45]
src/simulation/m_hb_function.fpp[68-74]

Agent prompt
The issue below was found during a code review. Follow the provided context and guidance below and implement a solution

### Issue description
`f_compute_shear_rate_from_components` clamps `shear_rate &gt;= 1e-2`, which forces a constant effective shear rate in all regions below that threshold. Since HB viscosity is a function of shear rate, this can materially change viscosity in low-shear regions.

### Issue Context
The yield term is well-behaved as `gdot -&gt; 0` in the Papanastasiou model, but computing `(1-exp(-m*gdot))/gdot` can suffer cancellation; that should be handled numerically rather than by a large floor.

### Fix Focus Areas
- src/simulation/m_hb_function.fpp[42-45]
- src/simulation/m_hb_function.fpp[68-74]

### Suggested fix
- Reduce the floor to a true numerical epsilon (e.g., `sgm_eps`-scaled) or make it configurable.
- In `f_compute_hb_viscosity`, compute the yield term with a stable formulation near zero shear, e.g. using `expm1`:
 `yield_term = tau0 * (-expm1(-hb_m_val*shear_rate)) / max(shear_rate, eps)`
 (or use the analytic limit `tau0*hb_m_val` when `shear_rate` is very small).

ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools


Grey Divider

ⓘ The new review experience is currently in Beta. Learn more

Grey Divider

Qodo Logo

Comment on lines +1061 to +1063
_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.

Comment on lines +748 to +754
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.

@coderabbitai
Copy link
Contributor

coderabbitai bot commented Mar 9, 2026

Caution

Review failed

Failed to post review comments

📝 Walkthrough

Walkthrough

Adds 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)
Check name Status Explanation
Title check ✅ Passed The title clearly and specifically describes the main change—adding a Herschel-Bulkley non-Newtonian viscosity model—which is the primary objective of the PR.
Description check ✅ Passed The PR description includes most required sections: Summary explains the implementation and validation approach; New files and Test plan are detailed. However, Type of change and Testing details are missing, and the Checklist items are not explicitly addressed.
Docstring Coverage ✅ Passed Docstring coverage is 100.00% which is sufficient. The required threshold is 80.00%.

✏️ Tip: You can configure your own custom pre-merge checks in the settings.

Tip

Try Coding Plans. Let us write the prompt for your AI agent so you can ship faster (with fewer bugs).
Share your feedback on Discord.


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.

❤️ Share

Comment @coderabbitai help to get the list of available commands and usage tips.

@sbryngelson sbryngelson marked this pull request as draft March 9, 2026 19:38
jasontruong2707 and others added 5 commits March 9, 2026 15:44
Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Format branch before merge to avoid whitespace conflicts in
m_riemann_solvers.fpp. Add new modules to module_categories.json.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
@github-actions
Copy link

Claude Code Review

Head SHA: ab1a4c8
Files changed: 294 (core physics: ~15 files)

Key files: src/simulation/m_hb_function.fpp (new), src/simulation/m_re_visc.fpp (new), src/simulation/m_viscous.fpp, src/simulation/m_riemann_solvers.fpp, src/common/m_derived_types.fpp, src/simulation/m_global_parameters.fpp, toolchain/mfc/params/definitions.py, toolchain/mfc/case_validator.py

Summary

  • Adds Herschel-Bulkley (HB) non-Newtonian viscosity with Papanastasiou regularization via two new modules (m_hb_function, m_re_visc)
  • Replaces the static Res_viscous/Res_gs arrays with a dynamic per-phase s_compute_re_visc called at every cell, every time step
  • All 4 parameter-system locations updated; runtime checks in m_checker.fpp and Python validator
  • Includes a notable IBM torque-calculation correctness fix (moved s_cross_product after viscous force accumulation)
  • Also bundles 2D modal Fourier / 3D spherical harmonic patch geometry additions and various CI infrastructure refactors

Findings

Critical / Correctness

1. Hard-coded shear-rate clamp is unphysical and problem-specific
src/simulation/m_hb_function.fpp, line ~73:

shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)

This clamps γ̇ to [10⁻², 10⁵] unconditionally. For highly viscous or very slow flows the lower bound of 0.01 s⁻¹ is orders of magnitude above realistic values and will artificially cap the yield-stress plateau viscosity. The upper bound likewise truncates high-shear behaviour. These should either be user-configurable (hb_gdot_min / hb_gdot_max) or at minimum documented as known limitations. As-is, a researcher simulating Re ≪ 1 Bingham flow will silently get wrong results.

2. Riemann-solver index mapping for NORM_DIR=2 and NORM_DIR=3 is incorrect
src/simulation/m_re_visc.fpp / m_riemann_solvers.fpp:

The Riemann solver loops over rotated index triplets (j, k, l) where the variable roles are permuted according to norm_dir. For NORM_DIR == 2, the x-normal interface loop uses (k, j, l) as (x, y, z) cell indices passed into s_compute_re_visc, so the gradient calculations in that routine (which index x_cc(j±1), y_cc(k±1), z_cc(l±1)) receive physically wrong coordinates. Specifically for NORM_DIR == 2:

call s_compute_re_visc(q_prim_vf, alpha_L, k, j, l, ...)

but inside s_compute_re_visc the first argument is treated as the x-index (j), the second as y (k), third as z (l). The velocity gradient routine then uses x_cc(k±1) as if it were x spacing, producing wrong shear rates. For Newtonian fluids this code path is not taken (early-exit in the if (.not. any_non_newtonian) branch), so existing tests won't catch this regression.

3. s_compute_re_visc is called with pre-computed gradients only from m_viscous, but NOT from m_riemann_solvers
In m_riemann_solvers.fpp all calls to s_compute_re_visc omit grad_x_vf/grad_y_vf/grad_z_vf. This means the Riemann-solver path always recomputes velocity gradients from scratch via s_compute_velocity_gradients_at_cell, which uses first-order one-sided differences at boundaries. This is inconsistent with the higher-order gradients already available in grad_x_vf at those call sites, and is a performance regression for every non-Newtonian Riemann solve.

Significant

4. Performance regression on Newtonian-only cases
Every call site now calls s_compute_re_visc unconditionally when viscous=T, even when any_non_newtonian=F. The any_non_newtonian flag gates the inner gradient/HB computation, but the per-phase loop and initialization still run. The old static Res_viscous lookup was O(Re_size) per cell; the new code is O(num_fluids) per cell with additional branch overhead. For production Newtonian GPU runs this is a measurable throughput regression. Consider keeping the static array path for the common Newtonian case.

5. fluid_pp struct fields not explicitly GPU-updated after MPI bcast
In s_initialize_global_parameters the new fields (fluid_pp(i)%non_newtonian, tau0, K, etc.) are initialized and MPI-broadcast. The struct as a whole is accessed on GPU via GPU_DECLARE(create='[...]') for scalar flags, but there is no GPU_UPDATE(device='[fluid_pp]') after the broadcast populates the new fields. The existing code updates any_non_newtonian via GPU_UPDATE, but the per-fluid parameters used inside the GPU kernel (fluid_pp(q)%tau0, etc.) may not reach device memory. This could cause silent wrong results or segfaults on GPU builds.

6. IBM non-Newtonian viscosity uses gdot=1 reference point — not local shear rate
src/simulation/m_ibm.fpp, lines ~1015-1022:

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 γ̇. This is acknowledged by the "compute reference viscosity at gdot=1" comment but should at minimum be a documented limitation with a @:PROHIBIT or warning if ib .and. any_non_newtonian.

Minor

7. File header says .f90 not .fpp
src/simulation/m_hb_function.fpp, line 2; src/simulation/m_re_visc.fpp, line 2:

!! @file m_hb_function.f90   ! should be .fpp
!! @file m_re_visc.f90       ! should be .fpp

8. s_compute_velocity_gradients_at_cell — mixed-precision shear components at boundaries
In the boundary fallback for D_xy (and D_xz), the routine only accumulates one of the two symmetric terms when the other direction lacks a neighbor. This results in D_xy ≈ 0.5*(du/dy) instead of the full symmetric average, underestimating shear at grid boundaries.

9. Missing regression test for non-Newtonian viscous flows
The PR validates against examples but no entry in toolchain/mfc/test/cases.py appears to be added for CI testing. Viscous tests exist (confirmed by the viscous_weno5_sgb_acoustic benchmark update), so a 2D non-Newtonian test case should be added to the automated test suite to prevent future regressions.


The IBM torque bug fix, the integral%zmin/zmax copy-paste fix, and the hyper_cleaning wave-speed correction are genuine improvements. The core architecture (per-fluid non_newtonian flag, any_non_newtonian short-circuit, Papanastasiou regularization) is sound. The shear-rate clamp and GPU data coherence issues should be addressed before merge.

Add Herschel-Bulkley parameter descriptions to descriptions.py.
Add fluid materials table entries and detailed section in case.md
with model equation, parameter table, and special cases.
Improve PHYSICS_DOCS entry with fuller explanation.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
@github-actions
Copy link

Claude Code Review

Head SHA: 4d042f5
Files changed: 294 (core Fortran sources: 33 files in `src/`)

Key files:

  • `src/simulation/m_hb_function.fpp` (new) — HB viscosity model
  • `src/simulation/m_re_visc.fpp` (new) — dynamic Re/viscosity computation
  • `src/simulation/m_viscous.fpp` — removes static `Res_viscous`, delegates to `m_re_visc`
  • `src/simulation/m_start_up.fpp` — moves `mytime += dt` after RK
  • `src/simulation/m_acoustic_src.fpp` — removes `t_step` arg, uses `mytime`
  • `src/common/m_derived_types.fpp` — adds HB fields to `physical_parameters`
  • `toolchain/mfc/params/definitions.py` — registers 8 new `fluid_pp` params

Summary:

  • Implements Herschel-Bulkley viscosity with Papanastasiou regularization for non-Newtonian flows
  • Replaces the static `Res_viscous` array with a dynamic per-cell `s_compute_re_visc` routine
  • Fixes a pre-existing `integral%zmin/zmax` init bug and the `mytime` timing relative to the RK stepper
  • No automated regression tests added to `toolchain/mfc/test/cases.py` for the new model

Findings

1. Shear-rate lower clamp is physically wrong — m_hb_function.fpp (line 71)

! Clamp for numerical safety
shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)

A minimum shear-rate of 1e-2 is far too large. In creeping viscous flows (Re ≪ 1), Poiseuille near-wall regions, or any stagnation zone, the physical shear rate is routinely 1e-6 to 1e-10. Clamping to 1e-2 imposes an artificial lower bound on viscosity and introduces large physical errors in exactly the regimes where non-Newtonian models are most important. The Papanastasiou regularization was specifically designed to eliminate the need for this kind of hard clamp. Replace with sgm_eps or at most 1.0e-10_wp.

2. Division by zero in f_compute_hb_viscositym_hb_function.fpp (line 43)

yield_term = tau0*(1._wp - exp_term)/shear_rate

If tau0 > 0 and shear_rate = 0, this divides by zero silently. The function is pure and public, so there's no guard; a caller passing shear_rate = 0 (e.g., a future use path that skips f_compute_shear_rate_from_components) would produce NaN. The correct limit as gdot → 0 via L'Hôpital is tau0 * hb_m. Either add an explicit if (shear_rate < epsilon(1._wp)) then branch returning tau0*hb_m_val + K_val*(epsilon(1._wp)**(nn_val-1._wp)), or add @:ASSERT(shear_rate > 0._wp, "shear_rate must be positive") to fail loudly.

3. No automated tests for non-Newtonian flows

toolchain/mfc/test/cases.py has no additions for the HB model. The PR test plan only covers manual runs of the example cases, which are not wired into ./mfc.sh test. Per the project contract, new physics features should have at least one entry in cases.py with a golden file. Please add a minimal 2D non-Newtonian test (e.g., a single-fluid power-law case nn ≠ 1, tau0 = 0) to provide CI regression coverage.

4. mytime timing change affects ALL simulations — m_start_up.fpp (line ~848)

Moving mytime = mytime + dt from before to after s_tvd_rk is coupled with replacing t_step*dt by mytime in s_acoustic_src_calculations. The acoustic source now evaluates at mytime (end of previous step), whereas before it used t_step*dt (beginning of current step with the old pre-increment). These are the same value with the old code. After this change, at t_step=0, mytime=0.0 (no increment yet), so the first source evaluation is at t=0 rather than t=dt. For cases relying on acoustic pulse timing with delay(ai), this is a one-dt shift. This behavioral change for existing acoustic cases should be called out explicitly in the PR description and validated against a known acoustic result.

5. File header typo — m_hb_function.fpp (line 2)

!! @file m_hb_function.f90

Should be m_hb_function.fpp.


Minor notes

  • integral%zmin/zmax fix (m_global_parameters.fpp ~line 825): the duplicate ymin/ymax initialization corrected to zmin/zmax is a genuine pre-existing bug fix — good catch.
  • s_finalize_viscous_module is now an empty stub with no body and no comment. Consider adding a brief comment explaining the `Res_viscous` removal so future readers understand the intent.
  • D_xy at boundary cells (m_re_visc.fpp ~line 130): at cells where only one of the two cross-derivative terms is available (boundary), only the available half is added. The result is no longer 0.5*(∂u/∂y + ∂v/∂x) but just one of the halves. This asymmetry is minor for interior-dominant physics but could matter in near-wall shear computations.

🤖 Generated with Claude Code

Fix .f90 -> .fpp in Doxygen headers for m_hb_function and m_re_visc.
Add shear-thinning and Bingham regression tests to CI test suite.
Re-merge with upstream/master after ruff PR (MFlowCode#1300) landed.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
@sbryngelson
Copy link
Member

Code Review — Non-Newtonian Viscosity

Thorough investigation of all findings from the automated review. Each item was verified against the source code.


Critical / Correctness

1. Index mapping bug in Riemann solver for NORM_DIR=2,3 ⚠️

src/simulation/m_riemann_solvers.fpp passes rotated loop indices to s_compute_re_visc, but src/simulation/m_re_visc.fpp interprets its j,k,l arguments as unrotated physical (x,y,z) indices and uses x_cc(j), y_cc(k), z_cc(l) for gradient computation.

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 s_compute_re_visc, D_xx is computed using x_cc(j') where j'=k (which spans the y-direction) — wrong grid spacing. Similarly for NORM_DIR=3.

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 m_riemann_solvers.fpp: lines ~464-480 (HLL), ~1268-1284 (HLLC), ~2180-2194 (HLLD), ~2827-2841, ~3266-3280.

Fix options:

  • (A) Pass unrotated physical indices to s_compute_re_visc (simplest)
  • (B) Make s_compute_re_visc accept a norm_dir argument and rotate internally

2. GPU coherence: fluid_pp not updated on device after MPI broadcast ⚠️

src/simulation/m_start_up.fpps_initialize_gpu_vars() (lines ~1231-1280) updates hundreds of variables on the GPU device but does not include fluid_pp. After MPI broadcast populates fluid_pp(i)%tau0, %K, %nn, etc. on the host, these values are never copied to device memory.

GPU kernels in s_compute_re_visc (called from GPU_PARALLEL_LOOP regions in m_riemann_solvers.fpp) read fluid_pp(q)%tau0, fluid_pp(q)%K, etc. from device memory, which contains stale/uninitialized values.

Impact: Silent wrong results on all GPU builds when any_non_newtonian=True.

Fix: Add $:GPU_UPDATE(device='[fluid_pp]') to s_initialize_gpu_vars() in m_start_up.fpp.


Significant

3. Hard-coded shear-rate clamp is not parameterizable

src/simulation/m_hb_function.fpp, line 73:

shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)

The bounds [10⁻², 10⁵] are hard-coded with no corresponding user parameters. For low-Re Bingham flows, the lower bound of 0.01 s⁻¹ is orders of magnitude above realistic values and will artificially cap the yield-stress plateau viscosity. No other MFC code has similar hard-coded clamps.

Recommendation: Either make these user-configurable (fluid_pp(i)%shear_rate_min/max) or document as a known limitation.


4. Missing pre-computed gradients in Riemann solver calls

All calls to s_compute_re_visc from m_riemann_solvers.fpp omit the optional grad_x_vf/grad_y_vf/grad_z_vf arguments. This forces recomputation of velocity gradients from scratch via s_compute_velocity_gradients_at_cell (1st-order one-sided differences at boundaries), even though higher-order gradients are already available at the Riemann solver call sites.

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 gdot=1 reference shear rate

src/simulation/m_ibm.fpp, lines ~1018-1024:

dynamic_viscosities(fluid_idx) = f_compute_hb_viscosity( &
    ..., 1._wp, ...)  ! Fixed shear rate = 1

The 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: s_compute_re_visc is guarded by if (viscous) at all call sites, and the any_non_newtonian flag inside the function gates the expensive gradient/HB computation. The Newtonian-only path is a simple O(Re_size) lookup. No measurable regression.


Minor (fixed in latest push)

7. File headers say .f90 instead of .fpp ✅ Fixed

m_hb_function.fpp line 2 and m_re_visc.fpp line 2 had @file m_*.f90 — corrected to .fpp.

8. Boundary shear rate asymmetry

s_compute_velocity_gradients_at_cell uses one-sided differences at boundaries, so D_xy only includes one of the two symmetric terms when a neighbor is missing. This underestimates shear at grid boundaries. Low impact for interior-dominated flows.

9. Missing regression test ✅ Fixed

Added two non-Newtonian test cases to toolchain/mfc/test/cases.py:

  • Shear-thinning (nn=0.5, tau0=0, K=0.1)
  • Bingham plastic (nn=1, tau0=0.1, K=0.01)

Both run in the single-fluid viscous test block for each dimension.


Summary

# Finding Status Severity
1 Index mapping bug (NORM_DIR=2,3) Needs fix Critical
2 GPU coherence (fluid_pp not on device) Needs fix Critical
3 Hard-coded shear-rate clamp Needs decision Medium
4 Missing pre-computed gradients in Riemann solver Enhancement Low
5 IBM fixed gdot=1 Document limitation Low
6 Newtonian performance No issue
7 File headers .f90→.fpp Fixed Trivial
8 Boundary shear asymmetry Known Low
9 Missing regression test Fixed High

Items 1 and 2 should be addressed before merge — they produce silently wrong results in 2D/3D and GPU builds respectively.

@github-actions
Copy link

Claude Code Review

Head SHA: 6f880a980fa505ae7568e1d3a39a74b1472e908b

Files changed: 28

Key files:

  • src/simulation/m_hb_function.fpp (new)
  • src/simulation/m_re_visc.fpp (new)
  • src/simulation/m_riemann_solvers.fpp
  • src/simulation/m_viscous.fpp
  • src/common/m_derived_types.fpp
  • toolchain/mfc/params/definitions.py
  • toolchain/mfc/case_validator.py
  • examples/2D_lid_driven_cavity_nn/case.py (new)
  • examples/2D_poiseuille_nn/case.py (new)
  • toolchain/mfc/test/cases.py

Summary

  • Adds Herschel-Bulkley (HB) non-Newtonian viscosity with Papanastasiou regularization, supporting power-law, Bingham, and HB fluids via per-fluid parameters
  • Replaces static pre-computed Res_viscous/Res_gs arrays with a dynamic per-cell viscosity routine (s_compute_re_visc)
  • Adds 8 new fluid_pp(i) parameters (non_newtonian, tau0, K, nn, mu_min, mu_max, mu_bulk, hb_m) with Python and Fortran-side validation
  • Includes two validation examples and 2 regression test cases

Findings

🔴 CRITICAL — Missing m_start_up.fpp namelist bindings

The CLAUDE.md parameter checklist requires updating 3 locations in Fortran. This PR updates m_global_parameters.fpp (defaults) and m_mpi_proxy.fpp (MPI broadcast) but does not update m_start_up.fpp in any of the three targets. Without namelist bindings in m_start_up.fpp, the new fluid_pp(i)%non_newtonian, %tau0, %K, etc. parameters are never read from the .inp file. The simulation will always use the Fortran-side defaults regardless of what the user specifies in their case file. This is a silent failure.

Required files (not present in diff):

  • src/simulation/m_start_up.fpp
  • src/pre_process/m_start_up.fpp
  • src/post_process/m_start_up.fpp

🔴 CRITICAL — Hard-coded shear rate clamp overrides user mu_min/mu_max

src/simulation/m_hb_function.fpp, line ~703:

shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)

This unconditionally clamps the shear rate to [0.01, 1e5] before the viscosity is computed, then f_compute_hb_viscosity clamps the resulting mu to [mu_min, mu_max]. The shear rate lower bound of 1e-2 is not physically motivated and will produce incorrect viscosities for flows with gdot < 0.01 (e.g., near-stagnation regions in lid-driven cavity). For a power-law fluid with n < 1, mu → ∞ as gdot → 0, and the user-supplied mu_max is the correct limiter. This clamp should be removed from f_compute_shear_rate_from_components; mu_min/mu_max in f_compute_hb_viscosity are sufficient.


🟠 SIGNIFICANT — Fortran mu_max check is logically dead

src/simulation/m_checker.fpp, line ~509:

@:PROHIBIT(fluid_pp(i)%mu_max < dflt_real .and. &
    fluid_pp(i)%mu_max <= fluid_pp(i)%mu_min, &
    "Non-Newtonian fluid mu_max must be > mu_min when set")

dflt_real is -huge(0._wp) (a large-negative sentinel). fluid_pp(i)%mu_max is initialized to dflt_real, meaning any user-supplied positive value satisfies mu_max > dflt_real, so the first condition mu_max < dflt_real is never true. This means the mu_max > mu_min constraint is never enforced in Fortran. The intended logic appears to be mu_max /= dflt_real .and. mu_max <= mu_min (i.e., mu_max has been set AND is invalid).


🟠 SIGNIFICANT — Performance regression for all-Newtonian cases

The PR removes pre-built static arrays Res_viscous, Res_gs, and Res_pr and replaces them with per-cell calls to s_compute_re_visc. For pure Newtonian cases (any_non_newtonian = .false.), the Newtonian branch of s_compute_re_visc still loops over Re_size(i) to look up fluid_pp(q)%Re(i) values, doing work on every face that was previously pre-computed once at startup. This affects all existing viscous simulations, not just non-Newtonian ones. The simplest fix is a guard: when any_non_newtonian is false, reintroduce the static array path.


🟠 SIGNIFICANT — Double velocity gradient computation in Riemann solver non-Newtonian path

In m_riemann_solvers.fpp, when any_non_newtonian is true and grad_x_vf/grad_y_vf/grad_z_vf are not passed, s_compute_re_visc calls s_compute_velocity_gradients_at_cell (finite differences from q_prim_vf) for every Riemann face. However, m_viscous.fpp already computes these gradients via the provided grad_x_vf etc. arrays. The Riemann solver path does not pass these gradient arrays to s_compute_re_visc, so gradients are computed twice per cell. Passing the pre-computed gradient arrays here would fix both the redundancy and ensure consistency between the viscous flux and the HB viscosity.


🟡 MINOR — L'Hôpital singularity handling not explicit

src/simulation/m_hb_function.fpp, line ~673:

yield_term = tau0*(1._wp - exp_term)/shear_rate

At shear_rate = 0, this is 0/0. The current mitigation relies on the shear rate clamp in f_compute_shear_rate_from_components. If that clamp is removed (as recommended above), a proper if (shear_rate < eps) yield_term = tau0 * hb_m_val branch is needed (from L'Hôpital: limit is tau0 * m).


🟡 MINOR — Re(1) still required alongside non_newtonian = T in examples

Both example cases set both fluid_pp(i)%Re(1) and fluid_pp(i)%non_newtonian = T. The Python validator (check_viscosity) still requires Re(1) to be set when viscous = T, even for non-Newtonian fluids. The non-Newtonian path ignores Re(1) at runtime, so requiring it is misleading. Consider relaxing the Re(1) requirement in check_viscosity when non_newtonian = T for that fluid.


🟡 MINOR — Test coverage is 1D only, examples use num_fluids=2

The regression tests in cases.py (2 new cases) run only in 1D with num_fluids=1 (since if num_fluids == 1). The example cases use num_fluids=2. A 2D test case with num_fluids=2 would catch index mapping bugs in the y-direction Riemann solver path.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
@github-actions
Copy link

Claude Code Review

Head SHA: fcf7646da9bb254f3761fc5d3e99a40e8bf4d601
Files changed: 28 — src/simulation/m_hb_function.fpp, src/simulation/m_re_visc.fpp, src/simulation/m_riemann_solvers.fpp, src/simulation/m_viscous.fpp, src/common/m_derived_types.fpp, src/simulation/m_global_parameters.fpp, src/simulation/m_checker.fpp, src/simulation/m_pressure_relaxation.fpp, toolchain/mfc/case_validator.py, toolchain/mfc/params/definitions.py (+18 more)

Summary

  • Adds Herschel-Bulkley (HB) non-Newtonian viscosity with Papanastasiou regularization via two new modules: m_hb_function.fpp (HB formula) and m_re_visc.fpp (dynamic Re dispatch)
  • Replaces static Res_viscous/Res_gs arrays with per-call s_compute_re_visc, touching all Riemann solver paths and m_viscous.fpp
  • New parameters are properly declared in derived types, default-initialized in m_global_parameters.fpp for all three targets, MPI-broadcast in all three m_mpi_proxy.fpp, and validated in m_checker.fpp + case_validator.py — but see finding Error on READE.md #1
  • Two example cases added; regression tests added in cases.py for shear-thinning and Bingham sub-cases

Findings

🔴 Critical — Missing namelist bindings in m_start_up.fpp (all 3 targets)

The 4-location parameter checklist (per CLAUDE.md and .claude/rules/parameter-system.md) requires:

  1. toolchain/mfc/params/definitions.py
  2. src/*/m_global_parameters.fpp
  3. src/*/m_start_up.fpp ❌ — MISSING
  4. toolchain/mfc/case_validator.py

None of the three m_start_up.fpp files appear in this diff. Without namelist bindings, the Fortran runtime will never read non_newtonian, tau0, K, nn, mu_min, mu_max, mu_bulk, or hb_m from the .inp input file. MPI rank 0 will always hold the defaults (non_newtonian = .false., K = 0, etc.), and MPI_BCAST will then propagate wrong defaults to all ranks. The feature will silently never activate in practice.

Fix required: Add all 8 new fluid_pp(i)% members to the namelist /user_inputs/ declaration in src/simulation/m_start_up.fpp, src/pre_process/m_start_up.fpp, and src/post_process/m_start_up.fpp.


🟠 Correctness — Coordinate-frame mismatch in off-axis Riemann solver calls (m_riemann_solvers.fpp)

s_compute_re_visc(q_prim_vf, alpha, j_arg, k_arg, l_arg, ...) uses j_arg, k_arg, l_arg as physical (x,y,z) indices internally — it accesses x_cc(j_arg), y_cc(k_arg), z_cc(l_arg) and velocity components momxb, momxb+1, momxb+2 accordingly.

In the Riemann solver, the loop runs in a rotated frame where j is the normal direction. For NORM_DIR == 2 the calls are:

call s_compute_re_visc(q_prim_vf, alpha_L, k, j, l, Re_visc_per_phase_L)

Here k (physical y-index) is passed as the x-position argument and j (physical x-index) as the y-position argument. s_compute_velocity_gradients_at_cell will then compute du/dx using cells at x_cc(k ± 1) (wrong coordinate), yielding a physically incorrect shear rate — and therefore incorrect viscosity — on y-normal faces for non-Newtonian fluids. Same issue applies for NORM_DIR == 3.

For the Newtonian code path (any_non_newtonian = .false.), s_compute_re_visc doesn't call s_compute_velocity_gradients_at_cell, so Newtonian cases are unaffected.


🟠 Correctness — s_compute_re_visc called without pre-computed gradients in Riemann solvers

m_viscous.fpp passes pre-computed grad_x_vf/grad_y_vf/grad_z_vf to s_compute_re_visc, which are computed centrally with proper boundary stencils. The Riemann solver calls (in m_riemann_solvers.fpp) omit these optional gradient arguments, so s_compute_velocity_gradients_at_cell re-derives gradients independently using a simpler first/second-order central stencil. This means the viscosity applied in the flux computation uses a different shear rate estimate than the viscous stress term — potential for inconsistency and spurious oscillations.


🟡 Robustness — Hardcoded shear-rate clamp in m_hb_function.fpp line 703

shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)

The floor 1e-2 is a physical parameter, not a numerical epsilon. For nearly-quiescent regions in a Bingham fluid (the classic use case for HB), the true shear rate may be much smaller than 1e-2, and this clamp directly limits the effective yield-stress approximation regardless of hb_m. The upper bound 1e5 could also be restrictive for high-Re problems. Consider exposing these as optional parameters or at minimum documenting them prominently, since they govern the accuracy of the Papanastasiou approximation.


🟡 Robustness — f_compute_hb_viscosity unsafe with shear_rate = 0 (m_hb_function.fpp line 673)

yield_term = tau0*(1._wp - exp_term)/shear_rate   ! division by shear_rate

The function relies on callers guaranteeing shear_rate > 0. The clamp in f_compute_shear_rate_from_components provides this safety for all call sites that go through that helper, but f_compute_hb_viscosity is public and takes an arbitrary shear_rate input. The IBM path (m_ibm.fpp line 732) hard-codes shear_rate = 1._wp as a workaround (which is a physically questionable "reference viscosity at unit shear rate"), confirming the interface is fragile. Consider asserting shear_rate > 0._wp at the top of f_compute_hb_viscosity via @:ASSERT(shear_rate > 0._wp, ...).


🟡 Behavioural — s_finalize_viscous_module GPU memory leak

s_initialize_viscous_module now ends with:

$:GPU_ENTER_DATA(copyin='[is1_viscous,is2_viscous,is3_viscous,iv]')

But s_finalize_viscous_module is now empty (the @:DEALLOCATE(Res_viscous) was removed and nothing replaced it). The matching GPU_EXIT_DATA(delete='[is1_viscous,is2_viscous,is3_viscous,iv]') should be present in finalize. This may be a pre-existing omission, but the PR makes it more visible by emptying the subroutine body.


🔵 Minor — K = 0 default passes Python validator silently

fluid_pp(i)%K defaults to 0._wp in Fortran. When a user sets non_newtonian = T without setting K, the Python validator's check_non_newtonian sees K is None (not yet assigned) and skips the K <= 0 check, so no Python-side error is raised. The Fortran checker in m_checker.fpp will still abort correctly, but the error will be a runtime abort rather than a pre-run validation message. Consider seeding K = 0 in the Python defaults so the validator can catch it.


🔵 Minor — Example case docstring mismatch (2D_lid_driven_cavity_nn/case.py lines 108–110)

The docstring says Re_gen = … = 5000 and lid velocity is described as Re=5000, n=1.5 but nn = 1.5 is shear-thickening. The title at line 106 says Re = 5000 while line 109 uses n = 1.5. This is internally consistent but the docstring description is brief enough to be confusing — the example parameters were validated against Li et al. (2015) for Re=500, but the file is configured for Re=5000. Worth a comment clarifying which benchmark this is.

Generated from 6 test cases (shear-thinning + Bingham, 1D/2D/3D).

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
@github-actions
Copy link

Claude Code Review

Head SHA: 9c15db4
Files changed: 40 (selected paths below)

Key files
  • src/simulation/m_hb_function.fpp (new)
  • src/simulation/m_re_visc.fpp (new)
  • src/simulation/m_riemann_solvers.fpp
  • src/simulation/m_viscous.fpp
  • src/simulation/m_pressure_relaxation.fpp
  • src/simulation/m_ibm.fpp
  • src/common/m_derived_types.fpp
  • toolchain/mfc/params/definitions.py
  • toolchain/mfc/case_validator.py
  • examples/2D_lid_driven_cavity_nn/case.py (new)

Summary

  • Adds Herschel-Bulkley viscosity with Papanastasiou regularization as a per-fluid option, gated behind fluid_pp(i)%non_newtonian
  • Replaces the pre-computed static Res_viscous/Res_gs/Res_pr arrays with a dynamic s_compute_re_visc call at every viscous evaluation site
  • Correctly threads new parameters through MPI broadcast, GPU declare, checker validation, and case validator
  • Two example cases and 6 new regression tests included

Findings

🔴 Critical — Missing namelist bindings in m_start_up.fpp

Per the project's 4-location rule, every new Fortran parameter must be bound in the namelist declaration in src/*/m_start_up.fpp. None of the three target m_start_up.fpp files appear in this diff.

Without namelist bindings, gfortran silently skips unknown namelist variables — meaning the 8 new parameters (non_newtonian, tau0, K, nn, mu_min, mu_max, mu_bulk, hb_m) are never read from the .inp file. Every HB simulation would silently run with defaults (K=0, nn=1, tau0=0), producing Newtonian behavior with no error.

Action required: Add all 8 fields to the &user_inputs namelist in src/simulation/m_start_up.fpp, src/pre_process/m_start_up.fpp, and src/post_process/m_start_up.fpp.


🔴 Critical — definitions.py adds only 3 lines for 8 new parameters

toolchain/mfc/params/definitions.py shows +3 / -0. With 8 new fluid_pp fields, most are missing from the parameter registry. Parameters not registered there won't be written to the .inp namelist file by the Python toolchain, compounding the m_start_up.fpp issue above.


🟡 Hard-coded shear rate clamp (m_hb_function.fpp, line 703)

shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)

The lower bound of 1e-2 is physically significant: at initialization when the flow is at rest, shear_rate = 0 gets clamped to 0.01, producing a non-zero artificial effective viscosity. For nn > 1 (shear-thickening) fluids this is especially large. The bounds are not tied to mu_min/mu_max and are undocumented. Consider deriving the bounds from the viscosity clamps or exposing them as parameters, and document the physical rationale.


🟡 Performance regression for all Newtonian viscous cases

The PR removes the pre-computed GPU-resident arrays Res_viscous, Res_gs, and Res_pr, replacing them with s_compute_re_visc at every evaluation site. For pure Newtonian flows (any_non_newtonian = .false.), the "all Newtonian" branch in s_compute_re_visc still iterates over fluid_pp (host-side struct) at every cell/face instead of reading from a pre-computed GPU array.

The old path was O(Re_size) array reads from GPU memory; the new path reads from fluid_pp per call. This is a performance regression for all existing viscous Newtonian simulations. Consider keeping the fast static path for the Newtonian case, or confirm via profiling that this overhead is acceptable.


🟡 IBM reference viscosity at shear_rate = 1.0 is undocumented (m_ibm.fpp, line 1728)

dynamic_viscosities(fluid_idx) = f_compute_hb_viscosity( &
    fluid_pp(fluid_idx)%tau0, ..., 1._wp, ...)  ! shear_rate = 1.0

The comment says "compute reference viscosity at gdot = 1" but provides no physical justification. For highly shear-dependent fluids, the viscosity at gdot=1 could be far from the actual wall shear rate, causing incorrect IBM force corrections. This should at least be documented with a physical rationale or an issue filed.


🟢 Positive observations

  • GPU macros ($:GPU_ROUTINE(parallelism='[seq]')) correctly applied to device-called functions in m_hb_function.fpp and m_re_visc.fpp
  • MPI broadcast added consistently across all three targets
  • @:PROHIBIT checks in m_checker.fpp cover all new parameters including the mu_max > mu_min invariant
  • Fypp-based index remapping in m_riemann_solvers.fpp for NORM_DIR == 1/2/3 looks correct
  • Dead viscous Re computation in m_pressure_relaxation.fpp correctly removed
  • any_non_newtonian flag avoids overhead of shear-rate computation in the common Newtonian case

Review generated by Claude Code (claude-sonnet-4-6)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

2 participants