Add Herschel-Bulkley non-Newtonian viscosity model#1298
Add Herschel-Bulkley non-Newtonian viscosity model#1298jasontruong2707 wants to merge 26 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. 4-space indentation in HB modules
|
![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.
📝 WalkthroughWalkthroughThis pull request introduces comprehensive support for Herschel-Bulkley non-Newtonian fluid viscosity modeling throughout the CFD simulation framework. New modules 🚥 Pre-merge checks | ✅ 4 | ❌ 1❌ Failed checks (1 warning)
✅ Passed checks (4 passed)
✏️ Tip: You can configure your own custom pre-merge checks in the settings. 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 ReviewHead SHA: 55ded73 Files changed:
Findings1. Non-Newtonian fluids without
|
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)
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).
|
Persistent review updated to latest commit 3c636c5 |
|
Persistent review updated to latest commit 3c636c5 |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 8
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
src/simulation/m_ibm.fpp (1)
988-1027:⚠️ Potential issue | 🟠 Major | 🏗️ Heavy liftUse stencil-local viscosity for the stress samples.
dynamic_viscosityis computed once at(i,j,k)and then reused fors_compute_viscous_stress_tensor(..., i±1, ...)and(..., j±1, ...). That is no longer correct for Herschel-Bulkley fluids, wheremudepends on the local shear rate. The resulting IB force/torque will be biased anywhere the shear varies across one cell.Compute
museparately at each neighbor location, or move the HB evaluation insides_compute_viscous_stress_tensorso each stencil sample uses its own viscosity.
🧹 Nitpick comments (1)
src/simulation/m_riemann_solvers.fpp (1)
326-336: 🏗️ Heavy liftConsider centralizing rotated→physical index mapping for
s_compute_re_visc.The same mapping logic is duplicated across several solver branches. A shared helper (left/right physical index builder from
norm_dir) would reduce drift and help prevent reintroducing direction-mapping bugs.Also applies to: 1021-1031, 1892-1902, 2514-2524, 2912-2922
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: CHILL
Plan: Pro
Run ID: 5496c2f1-bd8c-429e-bcf7-2c9e4a470ce7
⛔ Files ignored due to path filters (3)
examples/2D_lid_driven_cavity_nn/Lid_driven_cavity_Re_100_n_0.5.pngis excluded by!**/*.pngexamples/2D_lid_driven_cavity_nn/Lid_driven_cavity_Re_100_n_1.5.pngis excluded by!**/*.pngexamples/2D_poiseuille_nn/velocity_profile.pngis excluded by!**/*.png
📒 Files selected for processing (37)
docs/documentation/case.mddocs/module_categories.jsonexamples/2D_lid_driven_cavity_nn/case.pyexamples/2D_poiseuille_nn/case.pysrc/common/m_derived_types.fppsrc/post_process/m_global_parameters.fppsrc/post_process/m_mpi_proxy.fppsrc/pre_process/m_global_parameters.fppsrc/pre_process/m_mpi_proxy.fppsrc/simulation/m_checker.fppsrc/simulation/m_data_output.fppsrc/simulation/m_global_parameters.fppsrc/simulation/m_hb_function.fppsrc/simulation/m_ibm.fppsrc/simulation/m_mpi_proxy.fppsrc/simulation/m_pressure_relaxation.fppsrc/simulation/m_re_visc.fppsrc/simulation/m_riemann_solvers.fppsrc/simulation/m_time_steppers.fppsrc/simulation/m_viscous.fpptests/0CD345C7/golden-metadata.txttests/0CD345C7/golden.txttests/1E64A4D9/golden-metadata.txttests/1E64A4D9/golden.txttests/30713C1C/golden-metadata.txttests/30713C1C/golden.txttests/45DA4729/golden-metadata.txttests/45DA4729/golden.txttests/8B093EC0/golden-metadata.txttests/8B093EC0/golden.txttests/A46BB5B4/golden-metadata.txttests/A46BB5B4/golden.txttoolchain/mfc/case_validator.pytoolchain/mfc/params/definitions.pytoolchain/mfc/params/descriptions.pytoolchain/mfc/params/registry.pytoolchain/mfc/test/cases.py
| Re = 5000, n = 1.5 (shear-thickening) with mesh stretching near walls. | ||
|
|
||
| HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1) | ||
| Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000 | ||
|
|
||
| Mesh stretching: cosh-based clustering near all 4 walls (x_a, x_b, y_a, y_b). |
There was a problem hiding this comment.
Docstring describes mesh stretching, but the case config does not enable it.
The header says wall clustering/stretching is used, but no stretch_* toggles or stretching coefficients are set in the emitted JSON.
📝 Proposed fix
- Re = 5000, n = 1.5 (shear-thickening) with mesh stretching near walls.
+ Re = 5000, n = 1.5 (shear-thickening) on a uniform grid.
@@
- Mesh stretching: cosh-based clustering near all 4 walls (x_a, x_b, y_a, y_b).
+ (No mesh stretching parameters are enabled in this case file.)📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| Re = 5000, n = 1.5 (shear-thickening) with mesh stretching near walls. | |
| HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1) | |
| Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000 | |
| Mesh stretching: cosh-based clustering near all 4 walls (x_a, x_b, y_a, y_b). | |
| Re = 5000, n = 1.5 (shear-thickening) on a uniform grid. | |
| HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1) | |
| Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000 | |
| (No mesh stretching parameters are enabled in this case file.) |
| if (fluid_pp(q)%non_newtonian) then | ||
| ! Non-Newtonian: compute shear mu from HB model | ||
| mu_fluid = f_compute_hb_viscosity(fluid_pp(q)%tau0, fluid_pp(q)%K, fluid_pp(q)%nn, fluid_pp(q)%mu_min, & | ||
| & fluid_pp(q)%mu_max, shear_rate, fluid_pp(q)%hb_m) | ||
| Re_visc_per_phase(q, 1) = 1._wp/mu_fluid | ||
| ! Bulk viscosity | ||
| if (fluid_pp(q)%mu_bulk > 0._wp) then | ||
| Re_visc_per_phase(q, 2) = 1._wp/fluid_pp(q)%mu_bulk | ||
| else | ||
| Re_visc_per_phase(q, 2) = dflt_real | ||
| end if | ||
| else | ||
| ! Newtonian: return Re input values | ||
| $:GPU_LOOP(parallelism='[seq]') | ||
| do i = 1, 2 | ||
| if (Re_size(i) > 0) then | ||
| $:GPU_LOOP(parallelism='[seq]') | ||
| do r = 1, Re_size(i) | ||
| if (Re_idx(i, r) == q) then | ||
| Re_visc_per_phase(q, i) = fluid_pp(q)%Re(i) |
There was a problem hiding this comment.
Don't mix raw Reynolds numbers with 1/mu in Re_visc_per_phase.
s_compute_mixture_re assumes every entry in Re_per_phase is the same physical quantity and harmonically averages them. In this routine the HB branch stores 1/mu, but the Newtonian fallback stores fluid_pp(q)%Re(i) directly. A case with one Newtonian phase and one HB phase will therefore blend incompatible units and produce the wrong mixture viscosity/Re.
Please normalize the Newtonian path to the same basis as the HB path, or block mixed-mode setups during validation.
Also applies to: 298-307
| ! Map rotated (j,k,l) to physical (x,y,z) indices | ||
| #:if NORM_DIR == 1 | ||
| call s_compute_re_visc(q_prim_vf, alpha_L, j, k, l, Re_visc_per_phase_L) | ||
| call s_compute_re_visc(q_prim_vf, alpha_R, j + 1, k, l, Re_visc_per_phase_R) | ||
| #:elif NORM_DIR == 2 | ||
| call s_compute_re_visc(q_prim_vf, alpha_L, k, j, l, Re_visc_per_phase_L) | ||
| call s_compute_re_visc(q_prim_vf, alpha_R, k, j + 1, l, Re_visc_per_phase_R) | ||
| #:else | ||
| call s_compute_re_visc(q_prim_vf, alpha_L, l, k, j, Re_visc_per_phase_L) | ||
| call s_compute_re_visc(q_prim_vf, alpha_R, l, k, j + 1, Re_visc_per_phase_R) | ||
| #:endif | ||
| call s_compute_mixture_re(alpha_L, Re_visc_per_phase_L, Re_L) | ||
| call s_compute_mixture_re(alpha_R, Re_visc_per_phase_R, Re_R) | ||
| end if |
There was a problem hiding this comment.
Initialize viscous Re_L/Re_R for multi-fluid bubble-Euler runs.
Line 2514’s block still computes Re_L/Re_R only under num_fluids == 1. In viscous bubble-Euler cases with num_fluids > 1, those values can remain undefined and later contaminate Re_avg_rs* computations.
| if (.not. igr .or. dummy .or. any_non_newtonian) then | ||
| call s_convert_conservative_to_primitive_variables(q_cons_ts(1)%vf, q_T_sf, q_prim_vf, idwint) | ||
| end if |
There was a problem hiding this comment.
This non-Newtonian CFL path dereferences primitive storage that is never allocated in IGR mode.
q_prim_vf(:)%sf is only allocated in the .not. igr branch of this module, but this condition now converts into and reads from q_prim_vf whenever any_non_newtonian is true. On igr=.true. runs, that makes s_compute_dt write/read unallocated fields before calling s_compute_re_visc.
Either allocate the primitive buffers for the IGR path as well, or reject igr + any_non_newtonian before entering s_compute_dt.
Also applies to: 663-665
| def check_non_newtonian(self): | ||
| """Checks constraints on non-Newtonian fluid parameters""" | ||
| viscous = self.get("viscous", "F") == "T" | ||
| num_fluids = self.get("num_fluids", 1) | ||
|
|
||
| for i in range(1, num_fluids + 1): | ||
| nn_flag = self.get(f"fluid_pp({i})%non_newtonian", "F") == "T" | ||
| if not nn_flag: | ||
| continue | ||
|
|
||
| self.prohibit(not viscous, f"fluid_pp({i})%non_newtonian requires viscous = T") | ||
|
|
||
| K = self.get(f"fluid_pp({i})%K") | ||
| nn = self.get(f"fluid_pp({i})%nn") | ||
| tau0 = self.get(f"fluid_pp({i})%tau0") | ||
| mu_min = self.get(f"fluid_pp({i})%mu_min") | ||
| mu_max = self.get(f"fluid_pp({i})%mu_max") | ||
| hb_m = self.get(f"fluid_pp({i})%hb_m") | ||
|
|
||
| self.prohibit(K is not None and K <= 0, f"fluid_pp({i})%K (consistency index) must be > 0") | ||
| self.prohibit(nn is not None and nn <= 0, f"fluid_pp({i})%nn (flow behavior index) must be > 0") | ||
| self.prohibit(tau0 is not None and tau0 < 0, f"fluid_pp({i})%tau0 (yield stress) must be >= 0") | ||
| self.prohibit(mu_min is not None and mu_min < 0, f"fluid_pp({i})%mu_min must be >= 0") | ||
| self.prohibit(mu_max is not None and mu_min is not None and mu_max <= mu_min, f"fluid_pp({i})%mu_max must be > mu_min") | ||
| self.prohibit(hb_m is not None and hb_m <= 0, f"fluid_pp({i})%hb_m (Papanastasiou parameter) must be > 0") | ||
|
|
There was a problem hiding this comment.
check_non_newtonian() is currently ineffective (not invoked + missing required-field enforcement).
The new check is never called by any validate_* entrypoint, and within the method, core HB fields (K, nn, hb_m) are only range-checked when present, so missing values can pass toolchain validation and fail later.
🔧 Proposed fix
def check_non_newtonian(self):
"""Checks constraints on non-Newtonian fluid parameters"""
viscous = self.get("viscous", "F") == "T"
num_fluids = self.get("num_fluids", 1)
for i in range(1, num_fluids + 1):
nn_flag = self.get(f"fluid_pp({i})%non_newtonian", "F") == "T"
if not nn_flag:
continue
self.prohibit(not viscous, f"fluid_pp({i})%non_newtonian requires viscous = T")
K = self.get(f"fluid_pp({i})%K")
nn = self.get(f"fluid_pp({i})%nn")
tau0 = self.get(f"fluid_pp({i})%tau0")
mu_min = self.get(f"fluid_pp({i})%mu_min")
mu_max = self.get(f"fluid_pp({i})%mu_max")
hb_m = self.get(f"fluid_pp({i})%hb_m")
+ self.prohibit(K is None, f"fluid_pp({i})%K must be set when non_newtonian = T")
+ self.prohibit(nn is None, f"fluid_pp({i})%nn must be set when non_newtonian = T")
+ self.prohibit(hb_m is None, f"fluid_pp({i})%hb_m must be set when non_newtonian = T")
+
self.prohibit(K is not None and K <= 0, f"fluid_pp({i})%K (consistency index) must be > 0")
self.prohibit(nn is not None and nn <= 0, f"fluid_pp({i})%nn (flow behavior index) must be > 0")
self.prohibit(tau0 is not None and tau0 < 0, f"fluid_pp({i})%tau0 (yield stress) must be >= 0")
self.prohibit(mu_min is not None and mu_min < 0, f"fluid_pp({i})%mu_min must be >= 0")
self.prohibit(mu_max is not None and mu_min is not None and mu_max <= mu_min, f"fluid_pp({i})%mu_max must be > mu_min")
self.prohibit(hb_m is not None and hb_m <= 0, f"fluid_pp({i})%hb_m (Papanastasiou parameter) must be > 0") def validate_simulation(self):
"""Validate simulation-specific parameters"""
self.validate_common()
@@
self.check_body_forces()
self.check_viscosity()
+ self.check_non_newtonian()
self.check_mhd_simulation()📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| def check_non_newtonian(self): | |
| """Checks constraints on non-Newtonian fluid parameters""" | |
| viscous = self.get("viscous", "F") == "T" | |
| num_fluids = self.get("num_fluids", 1) | |
| for i in range(1, num_fluids + 1): | |
| nn_flag = self.get(f"fluid_pp({i})%non_newtonian", "F") == "T" | |
| if not nn_flag: | |
| continue | |
| self.prohibit(not viscous, f"fluid_pp({i})%non_newtonian requires viscous = T") | |
| K = self.get(f"fluid_pp({i})%K") | |
| nn = self.get(f"fluid_pp({i})%nn") | |
| tau0 = self.get(f"fluid_pp({i})%tau0") | |
| mu_min = self.get(f"fluid_pp({i})%mu_min") | |
| mu_max = self.get(f"fluid_pp({i})%mu_max") | |
| hb_m = self.get(f"fluid_pp({i})%hb_m") | |
| self.prohibit(K is not None and K <= 0, f"fluid_pp({i})%K (consistency index) must be > 0") | |
| self.prohibit(nn is not None and nn <= 0, f"fluid_pp({i})%nn (flow behavior index) must be > 0") | |
| self.prohibit(tau0 is not None and tau0 < 0, f"fluid_pp({i})%tau0 (yield stress) must be >= 0") | |
| self.prohibit(mu_min is not None and mu_min < 0, f"fluid_pp({i})%mu_min must be >= 0") | |
| self.prohibit(mu_max is not None and mu_min is not None and mu_max <= mu_min, f"fluid_pp({i})%mu_max must be > mu_min") | |
| self.prohibit(hb_m is not None and hb_m <= 0, f"fluid_pp({i})%hb_m (Papanastasiou parameter) must be > 0") | |
| def check_non_newtonian(self): | |
| """Checks constraints on non-Newtonian fluid parameters""" | |
| viscous = self.get("viscous", "F") == "T" | |
| num_fluids = self.get("num_fluids", 1) | |
| for i in range(1, num_fluids + 1): | |
| nn_flag = self.get(f"fluid_pp({i})%non_newtonian", "F") == "T" | |
| if not nn_flag: | |
| continue | |
| self.prohibit(not viscous, f"fluid_pp({i})%non_newtonian requires viscous = T") | |
| K = self.get(f"fluid_pp({i})%K") | |
| nn = self.get(f"fluid_pp({i})%nn") | |
| tau0 = self.get(f"fluid_pp({i})%tau0") | |
| mu_min = self.get(f"fluid_pp({i})%mu_min") | |
| mu_max = self.get(f"fluid_pp({i})%mu_max") | |
| hb_m = self.get(f"fluid_pp({i})%hb_m") | |
| self.prohibit(K is None, f"fluid_pp({i})%K must be set when non_newtonian = T") | |
| self.prohibit(nn is None, f"fluid_pp({i})%nn must be set when non_newtonian = T") | |
| self.prohibit(hb_m is None, f"fluid_pp({i})%hb_m must be set when non_newtonian = T") | |
| self.prohibit(K is not None and K <= 0, f"fluid_pp({i})%K (consistency index) must be > 0") | |
| self.prohibit(nn is not None and nn <= 0, f"fluid_pp({i})%nn (flow behavior index) must be > 0") | |
| self.prohibit(tau0 is not None and tau0 < 0, f"fluid_pp({i})%tau0 (yield stress) must be >= 0") | |
| self.prohibit(mu_min is not None and mu_min < 0, f"fluid_pp({i})%mu_min must be >= 0") | |
| self.prohibit(mu_max is not None and mu_min is not None and mu_max <= mu_min, f"fluid_pp({i})%mu_max must be > mu_min") | |
| self.prohibit(hb_m is not None and hb_m <= 0, f"fluid_pp({i})%hb_m (Papanastasiou parameter) must be > 0") |
When shear_rate == 0, the naive formula produces 0/0 (yield term) and 0^(nn-1) which is Inf for nn<1 (power-law term). Both corrupt viscosity and propagate into Re, fluxes, and dt. Fix: - Yield term: use analytic L'Hopital limit (1-exp(-m*g))/g -> m at g=0, so yield_term = tau0*hb_m when shear_rate <= verysmall - Power-law term: use g_eff = max(shear_rate, verysmall) so nn<1 fluids get a large-but-finite viscosity at rest, bounded by mu_max Addresses reviewer comment from qodo-code-review.
1. check_non_newtonian() was never called from validate_simulation(), so all non-Newtonian parameter constraints were silently skipped. Also add required-field checks for K, nn, hb_m and an IGR prohibition. 2. Re_visc_nn was missing from private(...) in all four GPU_PARALLEL_LOOP calls in m_viscous.fpp, causing thread race conditions on GPU builds. 3. igr + non_newtonian is now prohibited at validation time: q_prim_vf sub-arrays are only allocated in the non-IGR path, so s_compute_dt would dereference unallocated storage with any_non_newtonian=.true.
Line ending with & followed by a continuation line with only & and a comment is non-standard. Cray ftn (ftn-71) rejects it. Move the comment to its own line.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1298 +/- ##
==========================================
- Coverage 64.75% 64.69% -0.07%
==========================================
Files 71 73 +2
Lines 18721 18852 +131
Branches 1551 1559 +8
==========================================
+ Hits 12123 12196 +73
- Misses 5640 5693 +53
- Partials 958 963 +5 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
Automated Code Review — Follow-up Fixes AppliedSummary: 3 issues addressed. Fix branch: sbryngelson/MFC:pr-1298 Fixed
Remaining (not addressed)
|
- Fortran checker: prohibit non-Newtonian + igr combination (IGR does not allocate primitive variable storage needed by HB viscosity) - Fortran checker: require mu_max when nn < 1 (shear-thinning) to bound viscosity at zero shear rate - case_validator.py: matching Python-side shear-thinning mu_max check - cases.py: add weno_Re_flux=T variant to Non-Newtonian test block; generate golden files for 1D/2D/3D shear-thinning+weno_Re_flux
m_time_steppers.fpp: remove invalid trailing & on call to s_compute_moment_of_inertia (Cray ftn rejects blank continuation lines). m_re_visc.fpp: add #:include 'case.fpp'. Without it, MFC_CASE_OPTIMIZATION is undefined (False) during Fypp preprocessing, so AMD case-opt builds take the USING_AMD dimension(3) path while callers that include case.fpp take the dimension(num_fluids) path — causing 'actual argument has fewer elements' errors when num_fluids != 3.
….fpp Cray ftn ftn-7032 rejects derived-type dummy argument arrays (assumed size) inside OpenMP target regions. Pre-extract all non-Newtonian HB parameters from fluid_pp into plain local arrays before the GPU_PARALLEL_LOOP and add them to the copyin clause, eliminating fluid_pp access inside the loop. Also adds case.fpp include to fix MFC_CASE_OPTIMIZATION-gated AMD sizing.
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)