Add Herschel-Bulkley non-Newtonian viscosity model#1298
Add Herschel-Bulkley non-Newtonian viscosity model#1298jasontruong2707 wants to merge 23 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)
… merge artifact from PR)
Resolve conflicts: - m_data_output.fpp: keep master's ib_state_unit and ccfl removal - m_ibm.fpp: keep PR's non-Newtonian branch with master's Re(1) > 0 fix
- module_categories.json: keep m_hb_function, m_re_visc and add m_collisions - m_global_parameters.fpp: use eqn_idx%adv/E naming (master) - m_time_steppers.fpp: keep IB comment from master - cases.py: keep all 4 example skip entries (both branches) - m_ibm.fpp: remove old E_idx/momxb copyin; use eqn_idx%E; drop duplicate s_cross_product (now in m_helper.fpp) - m_data_output.fpp: take master s_write_ib_state_file subroutines; shorten docstrings; add q_T_sf param - m_riemann_solvers.fpp: keep NORM_DIR-aware s_compute_re_visc calls (non-Newtonian fix)
Claude Code ReviewPR #1298 — Add Herschel-Bulkley non-Newtonian viscosity model Critical Issues1. Hard-coded shear rate clamp overrides user viscosity bounds shear_rate = min(max(shear_rate, 1.0e-2_wp), 1.0e5_wp)The lower bound of 2. IBM force computation uses a fixed non-physical shear rate dynamic_viscosities(fluid_idx) = f_compute_hb_viscosity(..., 1._wp, ...)Viscosity at immersed boundary cells is evaluated at Important Issues3. @:PROHIBIT(fluid_pp(i)%mu_max < dflt_real .and. fluid_pp(i)%mu_max <= fluid_pp(i)%mu_min, ...
4. Mixed Newtonian/non-Newtonian fluids produce dimensionally inconsistent Re When
5.
Suggestions
Strengths
|
…D_xz/D_yz boundary fallback, weno_Re_flux validation
…ectives GPU_LOOP macros expand to OpenMP directives, which are not permitted inside pure subroutines (nvfortran S-0155). Remove pure from s_compute_re_visc and s_compute_mixture_re, which each contain multiple GPU_LOOP macros. The pure qualifier was semantically correct but syntactically invalid for older nvfortran (24.3, 24.5, 24.11) in OpenMP target offload builds.
The shear rate clamp removal (fix MFlowCode#2) slightly changes viscosity in low-shear-rate boundary cells, producing different output that the old golden did not capture. This regenerates the golden using the corrected code with --no-gpu on Phoenix (matching the gfortran CPU builds used by GitHub CI).
|
Persistent review updated to latest commit 3c636c5 |
| exp_term = exp(-hb_m_val*shear_rate) | ||
| yield_term = tau0*(1._wp - exp_term)/shear_rate | ||
| power_law_term = K_val*(shear_rate**(nn_val - 1._wp)) | ||
|
|
There was a problem hiding this comment.
1. Hb zero-shear nans 🐞 Bug ≡ Correctness
f_compute_hb_viscosity divides by shear_rate and evaluates shear_rate**(nn-1) without guarding shear_rate==0, which can produce NaN/Inf viscosities and propagate into Re/fluxes/dt. shear_rate can be exactly 0 for uniform velocity gradients or when gradient computation returns zeros at bounds, making this failure likely in real cases (e.g., initial rest).
Agent Prompt
### Issue description
`f_compute_hb_viscosity` can compute NaN/Inf when `shear_rate` is 0 (division by zero and `0**(nn-1)`), which will then corrupt `Re_visc_per_phase` and mixture viscosity.
### Issue Context
`shear_rate` is computed from velocity gradients and can legitimately be exactly zero in quiescent regions, or be forced to zero when gradient computation returns zeros at bounds.
### Fix Focus Areas
- src/simulation/m_hb_function.fpp[29-46]
### Implementation notes
- Add a small-shear branch, e.g. if `shear_rate <= eps`:
- use the analytic limit for the Papanastasiou yield term: `(1-exp(-m*g))/g -> m` so `yield_term = tau0*hb_m_val`
- avoid `0**(nn-1)` by using `g_eff = max(shear_rate, eps)` for the power-law term, or compute a bounded value consistent with `mu_max` if set.
- Consider using a numerically stable form for `1-exp(-x)` at small `x` (e.g. `-expm1(-x)` if available) to reduce cancellation.
- Ensure the final `mu` is finite and positive before returning (and before computing `1/mu`).
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
|
Persistent review updated to latest commit 3c636c5 |
| use m_derived_types !< Definitions of the derived types | ||
| use m_global_parameters !< Definitions of the global parameters | ||
|
|
||
| implicit none | ||
|
|
||
| private; public :: f_compute_hb_viscosity, f_compute_shear_rate_from_components |
There was a problem hiding this comment.
1. 4-space indentation in hb modules 📘 Rule violation ⚙ Maintainability
New Fortran modules use 4-space indentation instead of the required 2-space indentation, reducing consistency and readability across the codebase. This violates the project formatting convention for new/modified code.
Agent Prompt
## Issue description
New Fortran code added in this PR uses 4-space indentation, but the project requires 2-space indentation for new/modified Fortran sources.
## Issue Context
This affects the newly-added non-Newtonian viscosity implementation modules and should be fixed to maintain consistent formatting across the codebase.
## Fix Focus Areas
- src/simulation/m_hb_function.fpp[1-70]
- src/simulation/m_re_visc.fpp[1-312]
ⓘ Copy this prompt and use it to remediate the issue with your preferred AI generation tools
![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.) |
| exp_term = exp(-hb_m_val*shear_rate) | ||
| yield_term = tau0*(1._wp - exp_term)/shear_rate | ||
| power_law_term = K_val*(shear_rate**(nn_val - 1._wp)) | ||
|
|
||
| mu = yield_term + power_law_term | ||
| mu = max(mu, mu_min_val) | ||
| if (mu_max_val > dflt_real) mu = min(mu, mu_max_val) |
There was a problem hiding this comment.
Guard the zero-shear limit before dividing by shear_rate.
This kernel will hit undefined arithmetic when shear_rate is zero: line 40 becomes 0/0, and line 41 can overflow for shear-thinning inputs (nn < 1). That path is reachable immediately in the new quiescent Poiseuille case, so this can seed NaNs/Infs on the first step.
Suggested fix
- exp_term = exp(-hb_m_val*shear_rate)
- yield_term = tau0*(1._wp - exp_term)/shear_rate
- power_law_term = K_val*(shear_rate**(nn_val - 1._wp))
+ if (shear_rate <= tiny(1._wp)) then
+ yield_term = tau0*hb_m_val
+
+ if (nn_val < 1._wp) then
+ power_law_term = huge(1._wp)
+ else if (nn_val > 1._wp) then
+ power_law_term = 0._wp
+ else
+ power_law_term = K_val
+ end if
+ else
+ exp_term = exp(-hb_m_val*shear_rate)
+ yield_term = tau0*(1._wp - exp_term)/shear_rate
+ power_law_term = K_val*(shear_rate**(nn_val - 1._wp))
+ end if📝 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.
| exp_term = exp(-hb_m_val*shear_rate) | |
| yield_term = tau0*(1._wp - exp_term)/shear_rate | |
| power_law_term = K_val*(shear_rate**(nn_val - 1._wp)) | |
| mu = yield_term + power_law_term | |
| mu = max(mu, mu_min_val) | |
| if (mu_max_val > dflt_real) mu = min(mu, mu_max_val) | |
| if (shear_rate <= tiny(1._wp)) then | |
| yield_term = tau0*hb_m_val | |
| if (nn_val < 1._wp) then | |
| power_law_term = huge(1._wp) | |
| else if (nn_val > 1._wp) then | |
| power_law_term = 0._wp | |
| else | |
| power_law_term = K_val | |
| end if | |
| else | |
| exp_term = exp(-hb_m_val*shear_rate) | |
| yield_term = tau0*(1._wp - exp_term)/shear_rate | |
| power_law_term = K_val*(shear_rate**(nn_val - 1._wp)) | |
| end if | |
| mu = yield_term + power_law_term | |
| mu = max(mu, mu_min_val) | |
| if (mu_max_val > dflt_real) mu = min(mu, mu_max_val) |
| 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
| $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, rho_visc, gamma_visc, pi_inf_visc, alpha_visc_sum , & | ||
| & alpha_visc, alpha_rho_visc, Re_visc, tau_Re]') |
There was a problem hiding this comment.
Privatize the new viscosity temporaries in these GPU kernels.
These loops now write subroutine-scope temporaries like Re_visc_nn (and in the 3D blocks also rho_visc, gamma_visc, pi_inf_visc, and alpha_visc_sum) without listing them in private(...). On GPU builds that lets threads clobber each other and makes the cylindrical stress terms non-deterministic.
Suggested fix
- $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, rho_visc, gamma_visc, pi_inf_visc, alpha_visc_sum , &
- & alpha_visc, alpha_rho_visc, Re_visc, tau_Re]')
+ $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, rho_visc, gamma_visc, pi_inf_visc, alpha_visc_sum, &
+ & alpha_visc, alpha_rho_visc, Re_visc, Re_visc_nn, tau_Re]')
...
- $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, rho_visc, gamma_visc, pi_inf_visc, alpha_visc_sum , &
- & alpha_visc, alpha_rho_visc, Re_visc, tau_Re]')
+ $:GPU_PARALLEL_LOOP(collapse=3, private='[i, j, k, l, rho_visc, gamma_visc, pi_inf_visc, alpha_visc_sum, &
+ & alpha_visc, alpha_rho_visc, Re_visc, Re_visc_nn, tau_Re]')
...
- $:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_visc, alpha_rho_visc, Re_visc, tau_Re]')
+ $:GPU_PARALLEL_LOOP(collapse=3, private='[rho_visc, gamma_visc, pi_inf_visc, alpha_visc_sum, &
+ & alpha_visc, alpha_rho_visc, Re_visc, Re_visc_nn, tau_Re]')
...
- $:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_visc, alpha_rho_visc, Re_visc, tau_Re]')
+ $:GPU_PARALLEL_LOOP(collapse=3, private='[rho_visc, gamma_visc, pi_inf_visc, alpha_visc_sum, &
+ & alpha_visc, alpha_rho_visc, Re_visc, Re_visc_nn, tau_Re]')As per coding guidelines, "Declare private(...) on loop-local variables in GPU kernels to avoid unintended data sharing between threads."
Also applies to: 170-171, 258-258, 348-348
| CMake v3.26.5 on atl1-1-02-002-23-1.pace.gatech.edu | ||
|
|
||
| C : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc) | ||
| Fortran : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran) | ||
|
|
||
| PRE_PROCESS : OFF | ||
| SIMULATION : OFF | ||
| POST_PROCESS : OFF | ||
| SYSCHECK : OFF | ||
| DOCUMENTATION : ON | ||
| ALL : OFF | ||
|
|
||
| MPI : ON | ||
| OpenACC : OFF | ||
| OpenMP : OFF | ||
|
|
||
| Fypp : /storage/home/hcoda1/6/sbryngelson3/scratch/MFC-ruff/build/venv/bin/fypp | ||
| Doxygen : /storage/home/hcoda1/6/sbryngelson3/.local/bin/doxygen | ||
|
|
||
| Build Type : Release | ||
|
|
||
| Configuration Environment: | ||
|
|
||
| CC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc | ||
| CXX : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/g++ | ||
| FC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran | ||
| OMPI_CC : | ||
| OMPI_CXX : | ||
| OMPI_FC : | ||
|
|
||
| simulation: | ||
|
|
||
| CMake Configuration: | ||
|
|
||
| CMake v3.26.5 on atl1-1-02-002-28-2.pace.gatech.edu | ||
|
|
||
| C : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc) | ||
| Fortran : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran) | ||
|
|
||
| PRE_PROCESS : OFF | ||
| SIMULATION : ON | ||
| POST_PROCESS : OFF | ||
| SYSCHECK : OFF | ||
| DOCUMENTATION : OFF | ||
| ALL : OFF | ||
|
|
||
| MPI : ON | ||
| OpenACC : OFF | ||
| OpenMP : OFF | ||
|
|
||
| Fypp : /storage/scratch1/6/sbryngelson3/MFC-ruff/build/venv/bin/fypp | ||
| Doxygen : | ||
|
|
||
| Build Type : Release | ||
|
|
||
| Configuration Environment: | ||
|
|
||
| CC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc | ||
| CXX : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/g++ | ||
| FC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran | ||
| OMPI_CC : | ||
| OMPI_CXX : | ||
| OMPI_FC : | ||
|
|
||
| syscheck: | ||
|
|
||
| CMake Configuration: | ||
|
|
||
| CMake v3.26.5 on atl1-1-02-002-23-1.pace.gatech.edu | ||
|
|
||
| C : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc) | ||
| Fortran : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran) | ||
|
|
||
| PRE_PROCESS : OFF | ||
| SIMULATION : OFF | ||
| POST_PROCESS : OFF | ||
| SYSCHECK : ON | ||
| DOCUMENTATION : OFF | ||
| ALL : OFF | ||
|
|
||
| MPI : ON | ||
| OpenACC : OFF | ||
| OpenMP : OFF | ||
|
|
||
| Fypp : /storage/home/hcoda1/6/sbryngelson3/scratch/MFC-ruff/build/venv/bin/fypp | ||
| Doxygen : | ||
|
|
||
| Build Type : Release | ||
|
|
||
| Configuration Environment: | ||
|
|
||
| CC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc | ||
| CXX : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/g++ | ||
| FC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran | ||
| OMPI_CC : | ||
| OMPI_CXX : | ||
| OMPI_FC : | ||
|
|
||
| pre_process: | ||
|
|
||
| CMake Configuration: | ||
|
|
||
| CMake v3.26.5 on atl1-1-02-002-28-2.pace.gatech.edu | ||
|
|
||
| C : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc) | ||
| Fortran : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran) | ||
|
|
||
| PRE_PROCESS : ON | ||
| SIMULATION : OFF | ||
| POST_PROCESS : OFF | ||
| SYSCHECK : OFF | ||
| DOCUMENTATION : OFF | ||
| ALL : OFF | ||
|
|
||
| MPI : ON | ||
| OpenACC : OFF | ||
| OpenMP : OFF | ||
|
|
||
| Fypp : /storage/scratch1/6/sbryngelson3/MFC-ruff/build/venv/bin/fypp | ||
| Doxygen : | ||
|
|
||
| Build Type : Release | ||
|
|
||
| Configuration Environment: | ||
|
|
||
| CC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc | ||
| CXX : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/g++ | ||
| FC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran | ||
| OMPI_CC : | ||
| OMPI_CXX : | ||
| OMPI_FC : | ||
|
|
||
| post_process: | ||
|
|
||
| CMake Configuration: | ||
|
|
||
| CMake v3.26.5 on atl1-1-02-002-28-2.pace.gatech.edu | ||
|
|
||
| C : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc) | ||
| Fortran : GNU v12.3.0 (/usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran) | ||
|
|
||
| PRE_PROCESS : OFF | ||
| SIMULATION : OFF | ||
| POST_PROCESS : ON | ||
| SYSCHECK : OFF | ||
| DOCUMENTATION : OFF | ||
| ALL : OFF | ||
|
|
||
| MPI : ON | ||
| OpenACC : OFF | ||
| OpenMP : OFF | ||
|
|
||
| Fypp : /storage/scratch1/6/sbryngelson3/MFC-ruff/build/venv/bin/fypp | ||
| Doxygen : | ||
|
|
||
| Build Type : Release | ||
|
|
||
| Configuration Environment: | ||
|
|
||
| CC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gcc | ||
| CXX : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/g++ | ||
| FC : /usr/local/pace-apps/spack/packages/linux-rhel9-x86_64_v3/gcc-11.3.1/gcc-12.3.0-ukkkutsxfl5kpnnaxflpkq2jtliwthfz/bin/gfortran | ||
| OMPI_CC : | ||
| OMPI_CXX : | ||
| OMPI_FC : |
There was a problem hiding this comment.
Redact host/user/path details from committed golden metadata.
This snapshot commits internal hostnames and user-identifying absolute paths (for example at Line 13, Line 29, Line 36, Line 47, Line 63, Line 70). Please sanitize these fields in generated golden metadata (e.g., <host>, <path>, <user>) before committing to avoid infrastructure/PII leakage while keeping reproducibility signals.
| 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.68% -0.08%
==========================================
Files 71 73 +2
Lines 18721 18842 +121
Branches 1551 1559 +8
==========================================
+ Hits 12123 12188 +65
- Misses 5640 5691 +51
- Partials 958 963 +5 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
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)