|
MFC
Exascale flow solver
|
|
MFC
Exascale flow solver
|
Reference:
P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Computers & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
References:
P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Computers & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
H. Chen, C. Si, Y. Wu, H. Hu, and Y. Zhu, “Numerical investigation of the effect of equivalence ratio on the propagation characteristics and performance of rotating detonation engine”, Int. J. Hydrogen Energy, Mar. 2023. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.ijhydene.2023.03.190
Reference:
Trojak, W., & Dzanic, T. Positivity-preserving discoutinous spectral element method for compressible multi-species flows. arXiv:2308.02426
A temporally-evolving supersonic reacting shear layer between a hot air stream and an N2-diluted hydrogen stream. The base state comes from a 1-D flamelet solve (Cantera + Pyrometheus + JAX) extruded into 3D by hcid=371. This is the supersonic counterpart to examples/2D_reacting_mixing_layer, which runs the same flamelet machinery at Mc = 0.3.
| Parameter | Value |
|---|---|
| Oxidizer stream | air, X_O2 = 0.21, 500 K |
| Fuel stream | X_H2 = 0.5, balance N2, 300 K |
| Pressure | 101325 Pa |
| Vorticity thickness delta_omega | 1.0e-3 m |
| Convective Mach number Mc | 1.5 |
| Domain | [15, 20, 10] delta_omega in (x, y, z) |
| Grid | 14 pts/delta_omega in x and z, 28 in y, so 210 x 560 x 140 |
| Time step | 1e-9 s, RK3 |
| Numerics | WENO5 (mapped, monotonicity-preserving), HLLC |
| Boundaries | periodic in x and z, ghost-cell extrapolation in y |
| Chemistry | San Diego mechanism, 9 species, unity-Lewis transport |
x is streamwise, y is cross-stream (the flamelet profile axis), z is spanwise. The resolution density follows the temporal mixing-layer DNS of Wang et al. (Combustion and Flame, 2024), with the box reduced to a quarter of theirs in each direction.
The fuel stream is diluted because pure H2's sound speed is over 4x the oxidizer's, so reaching Mc = 1.5 would demand a velocity split of roughly 2650 m/s. X_H2 = 0.5 brings that to about 1394 m/s and moves the stoichiometric mixture fraction off the domain edge.
case.py calls flamelet_ic.py to solve a 1-D flamelet on the cross-stream grid and write it as prim.<n>.00.000000.dat under IC/. hcid=370 would extrude those files uniformly across z and leave the flow z-invariant, so this case uses hcid=371: it scales the file's cross-stream velocity by 1 + 0.5*cos(k_z z) and sets the spanwise component from the result, giving the IC 3D content at step 0. k_z = 2*pi/L_z is taken over the global domain, so the IC does not depend on the MPI decomposition.
The in-plane (x, y) perturbation is baked into the files by perturb_xy, seeded by perturb_seed in case.py. IC/ is gitignored and regenerated on a fresh checkout, so the fixed seed is what keeps the case reproducible. Regeneration is skipped when IC/ already matches the grid and the physical parameters, tracked in .cache_key.json.
The file spacing must match the run grid. A mismatch aborts in pre_process, so delete IC/ and let it regenerate after changing the grid or --scale.
./mfc.sh run examples/3D_reacting_mixing_layer/case.py -n 8
--scale shrinks the grid for cheap runs; --scale 0.05 gives 32^3, which is what the 3D -> Chemistry -> Reacting Mixing Layer regression test uses. --hot runs the full flamelet Newton/BDF solve instead of the default cold mollified profile.
The mechanism ships alongside the case as sandiego.yaml (UC San Diego Combustion Research Group, https://web.eng.ucsd.edu/mae/groups/combustion/mechanism.html).
A 2D flat plate pitching about its leading edge, 0 → 45° on an Eldredge smoothed ramp (a = 21), K = π/8 (case C1), Re_c = 300, Ma 0.2, plate thickness 2.5 % of chord. There is a published measurement to compare against:
Jantzen, Taira, Granlund & Ol, Phys. Fluids 26, 053606 (2014), Fig. 10, 2D panel, curve C1.
C_L = F_y / (½ ρ U² c) = 2 F_y here, with ρ = U = c = 1 and MFC's 2D force being per unit depth.
This case exists to measure the defect in #1849, not to fix it. Whoever does fix it needs a number to move.
NCELL sets how many cells lie across the plate thickness. The physical problem is identical at every level — same chord, same thickness, same domain, same times — so the sweep isolates the resolution requirement from any change of geometry, which a thickness sweep would confound.
SUMMARY=1 python3 case.py prints the grid and step count without running anything.
| cells across the thickness | peak C_L | C_L at t = 4 | rms difference from the reference |
|---|---|---|---|
| 2 | 6.413 | 1.209 | 20.9 % |
| 4 | 4.682 | 1.271 | 33.5 % |
| 8 | 4.458 | 1.111 | 27.9 % |
| 16 | 4.833 | 1.005 | 26.9 % |
| reference | 6.998 | 1.521 | — |
Two separate readings:
The peak converges from four cells up. 4.682, 4.458 and 4.833 across a fourfold refinement — a spread of 8 % with no trend. Two cells is genuinely under-resolved and 35 % out; four is enough. So a thin body does not need ten or more cells before the immersed boundary resolves it, which is worth knowing on its own for cost estimates.
The disagreement with the reference does not shrink. The last column sits at 27–34 % at every resolution and shows no sign of falling as the grid refines. That is the measurement that matters for #1849: it rules out under-resolution as the explanation and points at the force computation itself.
The obvious cheap test — a symmetric body whose true force is exactly zero — does not apply to a pitching plate, whose lift is large and unknown. For that style of check see examples/2D_ibm_force_decomposition (a cylinder, whose lift must be zero) and examples/3D_ibm_neighborhood_radius. This case trades the exact reference for a published one.
The scaling case can exercise both weak- and strong-scaling. It adjusts itself depending on the number of requested ranks.
This directory also contains a collection of scripts used to test strong and weak scaling on OLCF Frontier.
Pass --scaling weak. The --memory option controls (approximately) how much memory each rank should use, in Gigabytes. The number of cells in each dimension is then adjusted according to the number of requested ranks and an approximation for the relation between cell count and memory usage. The problem size increases linearly with the number of ranks.
Pass --scaling strong. The --memory option controls (approximately) how much memory should be used in total during simulation, across all ranks, in Gigabytes. The problem size remains constant as the number of ranks increases.
For example, to run a weak-scaling test that uses ~4GB of GPU memory per rank on 8 2-rank nodes with case optimization, one could:
Reference:
C. W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of Computational Physics 77 (2) (1988) 439–471. doi:10.1016/0021-9991(88)90177-5.
Reference: See Example 4.8.
A.S. Chamarthi, S.H. Frankel, A. Chintagunta, Implicit gradients based novel finite volume scheme for compressible single and multi-component flows, arXiv preprint arXiv:2106.01738 (2021).
Validates the combined non-Newtonian terms of MFC's Herschel-Bulkley viscosity against a closed-form analytic Poiseuille profile: a shear-thinning power law (nn = 0.5 < 1) and a yield stress (tau0 > 0) acting together. The companion examples isolate each effect — 2D_poiseuille_nn / 2D_poiseuille_thickening_nn (power-law only, tau0 = 0) and 2D_bingham_poiseuille_nn (yield only, nn = 1). The signature of a correct general Herschel-Bulkley model is a rigid plug near the centerline (where |tau| < tau0) joined to a shear-thinning sheared profile at the walls.
Single Papanastasiou-regularized Herschel-Bulkley fluid with both a sub-unity flow index and a finite yield stress:
| Parameter | Value | Role |
|---|---|---|
| fluid_pp(1)K | 1.5e-2 | consistency index |
| fluid_pp(1)nn | 0.5 | flow index < 1 -> shear-thinning |
| fluid_pp(1)tau0 | 3.5e-3 | yield stress -> plug half-width y0 = 0.35 H |
| fluid_pp(1)hb_m | 1.0e4 | sharp Papanastasiou yield regularization |
| fluid_pp(1)mu_min/mu_max | 1e-6 / 0.3 | viscosity clamp (rigid plug) |
Driven by g_x = 0.1, rho = 1, pres = 10, giving tau_w = rho*g*H = 1e-2, u_plug ~ 4e-3 and Mach ~1e-3. Channel L_y = 0.2, H = 0.1, no-slip walls, periodic in x. Grid m = 24 (x), n = 63 (y).
The plug viscosity diverges as the shear rate -> 0, so the clamp mu_max sets the plug rigidity. mu_max = 0.3 (~6x the wall effective viscosity) keeps a clear plug while keeping the explicit viscous timestep dt ~ dy^2 rho/mu_max tractable: dt scales as 1/mu_max, so set mu_max just above the physical maximum required.
The shear stress is tau = rho*g*(H - y); the fluid only flows where |tau| > tau0. With tau = tau0 + K*|du/dy|^n and tau_w = rho*g*H > tau0:
plug half-width : y0 = tau0/(rho*g)
sheared region : u(y) = (n/((n+1)*rho*g)) * K^(-1/n) *
[ (tau_w - tau0)^((n+1)/n)
- (rho*g*(H-y) - tau0)^((n+1)/n) ] (|y-H| >= y0)
plug : u_plug = (n/((n+1)*rho*g)) * K^(-1/n) *
(tau_w - tau0)^((n+1)/n) (|y-H| < y0)
(upper half mirrors about y = H). Requires tau_w > tau0 for any flow.
./mfc.sh run examples/2D_herschel_bulkley_poiseuille_nn/case.py -n 2 python examples/2D_herschel_bulkley_poiseuille_nn/compare_analytic.py
The run reaches t_stop = 0.4 in ~5 min on 2 CPU ranks with dt ~ 1e-5.
Relative L2 error vs. the analytic Herschel-Bulkley profile: 3.8% (2-rank run, steady-state drift between the last two saves ~1.1%). The sheared-region momentum balance K*|du/dy|^n + tau0 = rho*g*(H-y) holds to ~1% near the walls. A flat plug forms at the centerline: because the Papanastasiou plug is regularized (not perfectly rigid), the strict >=99% u_max band understates it, but the >=95% u_max plug half-width = 0.0359 = 1.03 y0, matching the analytic y0 = tau0/(rho*g) = 0.035 = 0.35 H.
Reference: See Example 4.13.
A.S. Chamarthi, S.H. Frankel, A. Chintagunta, Implicit gradients based novel finite volume scheme for compressible single and multi-component flows, arXiv preprint arXiv:2106.01738 (2021)., see Example 4.13
Reference:
Chamarthi, A., & Hoffmann, N., & Nishikawa, H., & Frankel S. (2023). Implicit gradients based conservative numerical scheme for compressible flows. arXiv:2110.05461
Reference:
V. A. Titarev, E. F. Toro, Finite-volume WENO schemes for three-dimensional conservation laws, Journal of Computational Physics 201 (1) (2004) 238–260.
Reference:
Woodward, P. (1984). The numerical simulation of two-dimensional fluid flow with strong shocks. Journal of Computational Physics, 54(1), 115–173. https://doi.org/10.1016/0021-9991(84)90140-2
This example case contains an automated convergence test using a 1D, two-component advection case. The case can be run by executing the bash script ./submitJobs.sh in a terminal after enabling execution permissions with chmod +x ./submitJobs.sh and setting the ROOT_DIR and MFC_DIR variables. By default the script runs the case for 6 different grid resolutions with 1st, 3rd, and 5th, order spatial reconstructions. These settings can be modified by editing the variables at the top of the script. You can also run different model equations by setting the ME variable and different Riemann solvers by setting the RS variable.
Once the simulations have been run, you can generate convergence plots with matplotlib by running python3 plot.py in a terminal. This will generate plots of the L1, L2, and Linf error norms and save the results to errors.csv.
Validates the power-law term of MFC's Herschel-Bulkley non-Newtonian viscosity against a closed-form analytic Poiseuille profile. Demonstrates the shear-thinning regime (nn < 1): the velocity profile is blunter (flatter-topped) than a parabola.
Single Papanastasiou-regularized Herschel-Bulkley fluid with no yield stress (tau0 = 0), so the effective viscosity is the pure power law mu = K*gamma_dot^(n-1), clamped to [mu_min, mu_max]:
| Parameter | Value | Role |
|---|---|---|
| fluid_pp(1)K | 2.0e-2 | consistency index |
| fluid_pp(1)nn | 0.7 | flow index < 1 -> shear-thinning |
| fluid_pp(1)tau0 | 0.0 | no yield stress (pure power law) |
| fluid_pp(1)hb_m | 1000.0 | Papanastasiou regularization parameter |
| fluid_pp(1)mu_min/mu_max | 1e-6 / 10.0 | viscosity clamp |
Driven by a constant body acceleration g_x = 8e-2, rho = 1, pres = 10 (sound speed ~3.74), giving u_max ~ 0.011 and Mach ~3e-3 (effectively incompressible). Channel L_y = 0.2, half-height H = L_y/2 = 0.1, no-slip walls at y = 0, L_y, periodic in x.
Fully-developed steady channel flow balances the body force against the shear stress, tau = rho*g*(H - y). With tau = K*|du/dy|^n (power law) the closed form is
u(y) = (n/(n+1)) * (rho*g/K)^(1/n) * ( H^((n+1)/n) - |y - H|^((n+1)/n) )
(n < 1 blunt/flat-topped; n = 1 parabola; n > 1 pointed). For n < 1 the viscosity diverges at the shear-free centerline, so any regularized solver caps it there; the near-wall momentum balance K*|du/dy|^n = rho*g*(H-y) is the cleanest pointwise correctness test and holds regardless of the cap.
./mfc.sh run examples/2D_poiseuille_nn/case.py -n 2 python examples/2D_poiseuille_nn/compare_analytic.py
Relative L2 error vs. the analytic power-law profile: 0.68% (2-rank run, steady state confirmed). The local momentum balance holds to ~0.1% across the channel, and the profile bluntness (mean/peak = 0.706) matches the n = 0.7 theory (n+1)/(2n+1) = 0.708.
Reference:
Bezgin, D. A., & Buhendwa A. B., & Adams N. A. (2022). JAX-FLUIDS: A fully-differentiable high-order computational fluid dynamics solver for compressible two-phase flows. arXiv:2203.13760
Ghia, U., & Ghia, K. N., & Shin, C. T. (1982). High-re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. Journal of Computational Physics, 48, 387-411
Reference:
G. B. Skinner and G. H. Ringrose, “Ignition Delays of a Hydrogen—Oxygen—Argon Mixture at Relatively Low Temperatures”, J. Chem. Phys., vol. 42, no. 6, pp. 2190–2192, Mar. 1965. Accessed: Oct. 13, 2024.
After running the simulation, compare MFC species mass fractions and induction time against a Cantera 0-D ideal-gas reactor reference:
This reads the Silo output, runs an equivalent Cantera reactor, prints the induction times (Skinner et al. / Cantera / (Che)MFC), and saves plots-nD_perfect_reactor-example.png. All dependencies are installed automatically by the MFC toolchain.
Validates the power-law term of MFC's Herschel-Bulkley non-Newtonian viscosity against a closed-form analytic Poiseuille profile, in the shear-**thickening** regime (nn > 1): the velocity profile is more pointed (sharper-topped) than a parabola. Companion to examples/2D_poiseuille_nn (shear-thinning, nn = 0.7).
Single Papanastasiou-regularized Herschel-Bulkley fluid with no yield stress (tau0 = 0), so the effective viscosity is the pure power law mu = K*gamma_dot^(n-1), clamped to [mu_min, mu_max]:
| Parameter | Value | Role |
|---|---|---|
| fluid_pp(1)K | 5.0e-2 | consistency index |
| fluid_pp(1)nn | 1.5 | flow index > 1 -> shear-thickening |
| fluid_pp(1)tau0 | 0.0 | no yield stress (pure power law) |
| fluid_pp(1)hb_m | 1000.0 | Papanastasiou regularization parameter |
| fluid_pp(1)mu_min/mu_max | 1e-6 / 0.035 | viscosity clamp |
Driven by a constant body acceleration g_x = 5e-2, rho = 1, pres = 10 (sound speed ~3.74), giving u_max ~ 0.013 and Mach ~3e-3 (effectively incompressible). Channel L_y = 0.2, half-height H = L_y/2 = 0.1, no-slip walls at y = 0, L_y, periodic in x. Grid m = 24 (x), n = 63 (y).
For n > 1 the maximum physical viscosity is at the wall (highest shear rate), mu_wall = K^(1/n) * (rho*g*H)^((n-1)/n) = 0.0232. mu_max = 0.035 ~ 1.5*mu_wall sits just above that maximum, so the clamp never activates (the analytic profile stays exact everywhere) while keeping the explicit viscous timestep large — the timestep scales as 1/mu_max, so set mu_max just above the physical maximum viscosity.
Fully-developed steady channel flow balances the body force against the shear stress, tau = rho*g*(H - y). With tau = K*|du/dy|^n (power law) the closed form is
u(y) = (n/(n+1)) * (rho*g/K)^(1/n) * ( H^((n+1)/n) - |y - H|^((n+1)/n) )
(n < 1 blunt/flat-topped; n = 1 parabola; n > 1 pointed). For n > 1 the effective viscosity mu = K*gamma_dot^(n-1) -> 0 at the shear-free centerline (rather than diverging as for n < 1), so no regularization cap is needed and the analytic profile is an exact reference everywhere.
./mfc.sh run examples/2D_poiseuille_thickening_nn/case.py -n 2 python examples/2D_poiseuille_thickening_nn/compare_analytic.py
The run reaches t_stop = 0.9 (~2 wall-viscous diffusion times) in ~1 min on 2 CPU ranks with dt = 8.4e-5.
Relative L2 error vs. the analytic power-law profile: 1.46% (2-rank run, steady-state drift between the last two saves 3.9e-3). The local momentum balance K*|du/dy|^n = rho*g*(H-y) holds to ~1.4% across the channel, and the profile bluntness (mean/peak = 0.626) matches the n = 1.5 theory (n+1)/(2n+1) = 0.625, confirming the pointed shear-thickening profile. u_max = 1.288e-2 matches the analytic 1.291e-2, confirming the clamp stays inactive.
Reference:
Coralic, V., & Colonius, T. (2014). Finite-volume Weno scheme for viscous compressible multicomponent flows. Journal of Computational Physics, 274, 95–121. https://doi.org/10.1016/j.jcp.2014.06.003
Validates the yield-stress term of MFC's Herschel-Bulkley non-Newtonian viscosity against a closed-form analytic Poiseuille profile. Demonstrates the Bingham regime (nn = 1, tau0 > 0): a rigid plug of uniform velocity forms near the centerline, where the shear stress falls below the yield stress.
Single Papanastasiou-regularized Herschel-Bulkley fluid with unit flow index, so K = mu is a plain Newtonian consistency and the only non-Newtonian effect is the yield stress:
| Parameter | Value | Role |
|---|---|---|
| fluid_pp(1)K | 5.0e-2 | n = 1 -> plain dynamic viscosity mu |
| fluid_pp(1)nn | 1.0 | flow index = 1 (Bingham, no power-law) |
| fluid_pp(1)tau0 | 4.0e-3 | yield stress -> plug half-width y0 = 0.4 H |
| fluid_pp(1)hb_m | 1.0e4 | sharp Papanastasiou yield regularization |
| fluid_pp(1)mu_min/mu_max | 1e-6 / 1.0 | viscosity clamp (rigid plug) |
Driven by g_x = 0.1, rho = 1, pres = 10, giving tau_w = rho*g*H = 1e-2, u_plug ~ 3.6e-3 and Mach ~1e-3. Channel L_y = 0.2, H = 0.1, no-slip walls, periodic in x.
The shear stress is tau = rho*g*(H - y); the fluid only flows where |tau| > tau0. With n = 1, K = mu, tau_w = rho*g*H > tau0:
plug half-width : y0 = tau0/(rho*g)
sheared region : u(y) = (1/(2*mu*rho*g)) *
[ (tau_w - tau0)^2 - (rho*g*(H-y) - tau0)^2 ] (|y-H| >= y0)
plug : u_plug = (1/(2*mu*rho*g)) * (tau_w - tau0)^2 (|y-H| < y0)
The signature of a correct yield term is the flat plug of uniform velocity within |y - H| < y0. Requires tau_w > tau0 for any flow.
./mfc.sh run examples/2D_bingham_poiseuille_nn/case.py -n 2 python examples/2D_bingham_poiseuille_nn/compare_analytic.py
Relative L2 error vs. the analytic Bingham profile: 2.5% (2-rank run, steady state confirmed). A plug forms at the centerline with half-width matching the analytic y0 = tau0/(rho*g) = 0.4 H.
Validates the immersed boundary (IBM) + non-Newtonian viscosity interaction: the channel's no-slip walls are two rectangular IB slabs instead of domain boundary conditions, so the flow exercises the per-stencil-sample Herschel-Bulkley viscosity mu_eff used by the IBM ghost-point and force machinery (s_compute_viscous_stress_tensor in m_viscous.fpp, consumed by s_compute_ib_forces in m_ibm.fpp). Companion to examples/2D_poiseuille_thickening_nn, which validates the same fluid against the same analytic profile with BC walls.
Domain x in [0, 0.2] (periodic), y in [0, 0.3], grid m = 24, n = 95 (dy = 0.003125). Two rectangle IB slabs (patch_ibgeometry = 3, no-slip) bound the flow gap y in [0.05, 0.25] (half-height H = 0.1, centerline y_c = 0.15, 64 cells across the gap). Each slab extends beyond the domain in x and mostly outside the domain in y, so the only IB surface seen by the flow is its flat gap face; the slab centroids sit just inside the domain (a centroid exactly on the boundary is owned by no rank and its ib_state force record is never written). The domain BCs behind the slabs are no-slip walls (bc_y = -16), which keep the body-forced dead fluid inside the slabs benign. patch_ibmass = 0 so the reported IB force is the pure pressure+viscous volume integration (no bf_x*mass bookkeeping term); ib_state_wrt = T writes it at every save.
Fluid and forcing match the BC-walled template (single Papanastasiou-regularized Herschel-Bulkley fluid, no yield stress):
| Parameter | Value | Role |
|---|---|---|
| fluid_pp(1)K | 5.0e-2 | consistency index |
| fluid_pp(1)nn | 1.5 | flow index > 1 -> shear-thickening |
| fluid_pp(1)tau0 | 0.0 | no yield stress (pure power law) |
| fluid_pp(1)hb_m | 1000.0 | Papanastasiou regularization parameter |
| fluid_pp(1)mu_min/mu_max | 1e-6 / 0.035 | viscosity clamp (~1.5x mu_wall = 0.0232, clamp inactive) |
Driven by g_x = 5e-2, rho = 1, pres = 10 (Mach ~3e-3); cfl_adap_dt with cfl_target = 0.3 to t_stop = 0.9 (~2 wall-viscous diffusion times).
In the gap the steady fully-developed power-law profile is
u(y) = (n/(n+1)) * (rho*g/K)^(1/n) * ( H^((n+1)/n) - |y - y_c|^((n+1)/n) )
and the steady x-force per wall per unit depth is tau_w * L_x with tau_w = rho*g*H = 5e-3.
./mfc.sh run examples/2D_ibm_poiseuille_nn/case.py -n 2 ./build/venv/bin/python3 examples/2D_ibm_poiseuille_nn/compare_analytic.py
(~2.5 min on 2 CPU ranks.) For the n = 1 equivalence check, run the newtonian and nn1 modes into two scratch copies and compare:
for MODE in newtonian nn1; do
mkdir -p build/ibm_nn_equiv/$MODE
cp examples/2D_ibm_poiseuille_nn/case.py build/ibm_nn_equiv/$MODE/
IBM_NN_MODE=$MODE ./mfc.sh run build/ibm_nn_equiv/$MODE/case.py -n 2
done
./build/venv/bin/python3 examples/2D_ibm_poiseuille_nn/check_equivalence.py \
build/ibm_nn_equiv/newtonian build/ibm_nn_equiv/nn1
A — n = 1 Newtonian equivalence (IBM mu_eff degeneracy). A power-law fluid with nn = 1, tau0 = 0, K = 0.02 is analytically the same fluid as a Newtonian one with mu = 0.02. Running both modes with the same fixed dt = 6e-5 to t = 0.3 gives bitwise identical velocity fields (max abs and rel L2 difference 0.0) and bitwise identical IBM-integrated wall forces — the non-Newtonian IBM path reduces exactly to the Newtonian one.
B — analytic power-law profile (n = 1.5). Relative L2 error of the steady x-averaged gap profile vs. the analytic solution: 5.8% with the nominal H = 0.1 (slab faces), 2.6% with H = 0.1016 fitted from the u -> 0 crossings (IBM walls are sharp only to ~half a cell; the fitted walls sit ~dy/2 = 0.0016 outside the faces). Steady-state drift between the last two saves 2.0e-3; profile bluntness (mean/peak) 0.640 vs. the n = 1.5 theory 0.625 (parabola 0.667), confirming the pointed shear-thickening profile. The BC-walled template achieves 1.46% on the same fluid; the extra error is the diffuse-wall representation, not the viscosity model.
C — IBM-integrated wall force. The volume-integrated x-force converges to 8.05e-4 per wall (both walls identical by symmetry; plateaued by t = 0.9) vs. the analytic tau_w*L_x = 1.0e-3 — a ratio of 0.80. The deficit is the known coarseness of the volume-integration force estimator (second-order finite differences of ghost/dead-cell states inside the body), not the viscosity model: in Validation A the same integral is bitwise identical between the Newtonian and non-Newtonian code paths.
Reference:
P. D. Lax, Weak solutions of nonlinear hyperbolic equations and their numerical computation, Communications on pure and applied mathematics 7 (1) (1954) 159–193.
s_apply_boundary_patches (src/pre_process/m_boundary_conditions.fpp) dispatches by dimensionality:
There is no else. A geometry belonging to the other dimensionality falls straight through: the patch is never applied, nothing is printed, and the face silently keeps whatever bc_[xyz] gave it.
That is quiet in the worst way. A nozzle cut into a no-slip wall with a 3D geometry in a 2D case simply stays a solid wall — the run completes, writes output, and the jet has a velocity of exactly zero for all time with no indication anything was ignored.
| before | after | |
|---|---|---|
| GEOMETRY=1 | passes | passes |
| GEOMETRY=3 | passes, then does nothing at run time | patch_bc(1)geometry must be 1 (line segment) in 2D; geometry 3 is never applied |
The 3D direction is symmetric: geometry 1 in a 3D case is equally ignored, and is now equally refused.
The case carries a second initial-condition patch a few cells thick at the inlet, and it is load-bearing: the Dirichlet buffer is filled by pre_process from the initial condition at the boundary face. A domain initialised at rest stores rest in that buffer, and the "inflow" then delivers zero for all time — again with no warning. This is not what the validator change addresses; it is noted here because the two failures look identical from the outside, and knowing that saves working out which one is in play.
Validator only; no source change and no golden files. All 182 example cases still validate.
Qualitative demonstration of MFC's Herschel-Bulkley non-Newtonian viscosity in a recirculating flow. A shear-thinning fluid fills a unit square cavity driven by a moving top lid. Unlike the Poiseuille examples, this case has no closed-form analytic solution — it is a qualitative demonstration of the expected shear-thinning trend (a primary vortex center shifted toward the moving lid, with stronger near-wall velocity gradients relative to the Newtonian case).
Two identical Papanastasiou-regularized Herschel-Bulkley fluids (a two-fluid setup sharing one rheology), pure power law (tau0 = 0):
| Parameter | Value | Role |
|---|---|---|
| fluid_pp(i)K | 1.0e-2 | consistency index; Re_eff = 1/K = 100 at unit shear rate |
| fluid_pp(i)nn | 0.5 | flow index < 1 -> shear-thinning |
| fluid_pp(i)tau0 | 0.0 | no yield stress (pure power law) |
| fluid_pp(i)hb_m | 1000.0 | Papanastasiou regularization parameter |
| fluid_pp(i)mu_min/mu_max | 1e-6 / 1.0 | viscosity clamp |
Unit square [0,1]^2, m = n = 99 (coarse smoke-run grid), all walls no-slip with the top lid moving at bc_yve1 = 0.5. Effective Reynolds number Re_eff = 1/K = 100 at unit shear rate; the conventional lid-based Reynolds number rho*U*L/mu(1) = 0.5/1e-2 = 50 with U = 0.5 and mu(1) = K.
The auto-registered CI test of this example is truncated to 50 time steps and serves as smoke coverage only, not a physics anchor.
Incompressible recirculating cavity flow with the shear-dependent power-law viscosity mu = K*gamma_dot^(n-1). With n < 1 the fluid thins under the strong shear beneath the lid and in the corner boundary layers, while the slow cavity core stays comparatively viscous.
What to look for (qualitative, no analytic match): a primary recirculating vortex whose center, relative to a Newtonian cavity at the same Reynolds number, is shifted toward the moving lid, with stronger near-wall velocity gradients — the expected shear-thinning trend. Do not expect a quantitative error; the committed grid (m = n = 99) is intentionally coarse. For quantitative comparison use m = n = 499 or finer with a longer t_step_stop.
./mfc.sh run examples/2D_lid_driven_cavity_nn/case.py -n 2
Post-process and inspect the velocity / vorticity (omega_wrt(3)) fields; there is no compare_analytic.py for this qualitative case.
Reference: See Example 4.18.
A.S. Chamarthi, S.H. Frankel, A. Chintagunta, Implicit gradients based novel finite volume scheme for compressible single and multi-component flows, arXiv preprint arXiv:2106.01738 (2021).
Reference:
Hillewaert, K. (2013). TestCase C3.5 - DNS of the transition of the Taylor-Green vortex, Re=1600 - Introduction and result summary. 2nd International Workshop on high-order methods for CFD.
This figure shows the isosurface with zero q-criterion.
Reference:
Panchal et. al., A Seven-Equation Diffused Interface Method for Resolved Multiphase Flows, JCP, 475 (2023)
