|
MFC
Exascale flow solver
|
|
MFC
Exascale flow solver
|
Verifies the Fourier heat conduction term div(k grad T) on the energy equation against the exact Laplacian of a sinusoidal temperature field.
Ideal gas, so T = p / ((Gamma - 1) * rho * cv). Setting
rho(x) = RHO0 / (1 + A sin(2 pi x / L))
at uniform pressure gives exactly
T(x) = T0 (1 + A sin(2 pi x / L)), T0 = P0 / ((Gamma - 1) RHO0 cv).
At t = 0 the velocity is zero and the pressure is uniform, so every Euler flux vanishes and the entire right-hand side is conduction:
d(rho E)/dt = k d2T/dx2 = -k T0 A (2 pi / L)^2 sin(2 pi x / L).
One time step therefore measures the conduction term in isolation, up to O(dt) time error and O(dx^2) space error. compare_analytic.py forms (rho E(dt) - rho E(0)) / dt from restart_data/lustre_{0,1}.dat and compares it to that expression.
for nx in 50 100 200 400; do
NX=$nx ./mfc.sh run examples/1D_conduction_convergence/case.py -n 1
NX=$nx ./build/venv/bin/python3 examples/1D_conduction_convergence/compare_analytic.py
done
| NX | L2 error | relative | order |
|---|---|---|---|
| 50 | 9.174206e-06 | 1.314570e-03 | – |
| 100 | 2.293085e-06 | 3.285757e-04 | 2.000 |
| 200 | 5.775688e-07 | 8.275971e-05 | 1.989 |
| 400 | 1.442740e-07 | 2.067299e-05 | 2.001 |
| 800 | 4.396086e-08 | 6.299142e-06 | 1.715 |
Second order through NX=400. The drop at NX=800 is the O(dt) time-splitting floor, which by then is comparable to the spatial error; refining dt restores second order.
The result is unchanged on 2 ranks (-n 2), which exercises the MPI temperature halo exchange.
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:
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.
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.
Axisymmetric twin of examples/1D_conduction_convergence. It verifies the radial part of div(k grad T) on a cylindrical grid, and in particular the axis cell k = 0, which has no face pair and is therefore the one cell fed by s_compute_conduction_axis_source.
Ideal gas, so T = p / ((Gamma - 1) * rho * cv). Setting
rho(r) = RHO0 / (1 + A (1 - r^2/R^2) + B (1 - r^4/R^4))
at uniform pressure gives exactly
T(r) = T0 (1 + A (1 - r^2/R^2) + B (1 - r^4/R^4)), T0 = P0 / ((Gamma - 1) RHO0 cv),
whose cylindrical Laplacian is finite on the axis:
(1/r) d/dr (r dT/dr) = -4 A T0 / R^2 - 16 B T0 r^2 / R^4.
At t = 0 the velocity is zero and the pressure is uniform, so every Euler flux and every cylindrical geometric source vanishes and the entire right-hand side is conduction:
d(rho E)/dt = k (1/r) d/dr (r dT/dr).
compare_analytic.py forms (rho E(dt) - rho E(0)) / dt from restart_data/lustre_{0,1}.dat and compares it to that expression.
B = 0 (the default) makes the exact answer a constant: every radial cell, axis included, must return the same number. That is the point of the profile — a sign error or a stray factor in the axis cell is unmissable against a flat field, where a sinusoid would hide it. B > 0 makes the answer vary with r and turns the same case into an ordinary convergence test.
src/pre_process/m_grid.f90 gives the axisymmetric grid a half-width cell at r = 0: cell 0 spans [0, h/2] with center h/4, and every cell above it has width h. Two cells therefore behave differently from the bulk:
# constant exact answer: the axis-cell check
for nr in 25 50 100 200 400; do
NR=$nr ./mfc.sh run examples/2D_axisym_conduction_convergence/case.py -n 1
NR=$nr ./build/venv/bin/python3 examples/2D_axisym_conduction_convergence/compare_analytic.py
done
# r-dependent exact answer: the second-order convergence table
for nr in 25 50 100 200; do
NR=$nr BAMP=0.4 ./mfc.sh run examples/2D_axisym_conduction_convergence/case.py -n 1
NR=$nr BAMP=0.4 ./build/venv/bin/python3 examples/2D_axisym_conduction_convergence/compare_analytic.py
done
Constant exact answer (B = 0, exact = -4.000000e-02 in every cell):
| NR | axis cell | axis rel err | cell 1 rel err | max bulk rel err |
|---|---|---|---|---|
| 25 | -4.00000033e-02 | -8.274e-08 | +3.125e-02 | 8.274e-08 |
| 50 | -4.00000033e-02 | -8.274e-08 | +3.125e-02 | 1.193e-06 |
| 100 | -4.00000033e-02 | -8.274e-08 | +3.125e-02 | 3.413e-06 |
| 200 | -4.00000033e-02 | -8.274e-08 | +3.125e-02 | 7.854e-06 |
| 400 | -3.99999589e-02 | +1.027e-06 | +3.125e-02 | 1.007e-05 |
The axis cell matches the interior to the O(dt) time-splitting floor at every resolution: for a quadratic profile the discretization is exact in space, so nothing else is left. Before the distance weighting was added the axis cell read -5.400e-02, 35% high and independent of NR, which is what this table is built to expose.
Second-order convergence (B = 0.4, r-dependent exact answer):
| NR | axis rel err | order | bulk L2 rel err | order |
|---|---|---|---|---|
| 25 | 1.171e-03 | – | 9.623e-04 | – |
| 50 | 2.856e-04 | 2.04 | 2.348e-04 | 2.03 |
| 100 | 7.062e-05 | 2.02 | 5.799e-05 | 2.02 |
| 200 | 1.757e-05 | 2.01 | 1.448e-05 | 2.00 |
The axis cell converges at the same second order as the bulk. Cell 1 sits at 1.19e-02 at every resolution, as expected from the shared face-gradient form described above.
Reference:
Panchal et. al., A Seven-Equation Diffused Interface Method for Resolved Multiphase Flows, JCP, 475 (2023)
Reference:
Chamarthi, A., & Hoffmann, N., & Nishikawa, H., & Frankel S. (2023). Implicit gradients based conservative numerical scheme for compressible flows. arXiv:2110.05461
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.
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
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.
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
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.
A temporally-evolving supersonic reacting shear layer between a hot air stream and an N2-diluted hydrogen stream. The base state comes from Cantera-only stream mixing (or a counterflow flame with --hot) 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. --hot runs the full Cantera counterflow flame solve instead of the default cold mollified profile. flame_strain_rate in case.py sets the nominal inlet strain rate (default 100/s). The hot profile is mapped by mixture fraction onto the prescribed shear layer; it replaces the former scalar-dissipation-matched flamelet initialization. See Thermochemistry implementation for the model and validation details.
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).
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
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
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.
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.
s_open_probe_files appends whenever D/probe*_prim.dat already exists:
That is correct when a run is being continued. It is wrong when one is being started over, which is what happens every time a case is re-run in place after a parameter change. The second run's rows land on top of the first's, nothing in the file marks the join, and the time column simply resets partway down. A reader sees one monotonic series and is silently wrong.
There is no warning, no header and no separator. The two runs need not even share a grid — a case whose resolution changed between runs produces a file whose first half was recorded at different probe locations.
| rows in D/probe1_prim.dat | |
|---|---|
| before | 40 — two runs of 20, spliced |
| after | 20 |
And the time column, read straight through, before the fix:
Append only when continuing: t_step_start > 0, or n_start > 0 under cfl_dt. A fresh start replaces the file, which is what every other output MFC writes already does. s_open_com_files had the same pattern and gets the same treatment.
This produced three separate wrong numbers in one project before it was noticed. The worst was a jet whose probe files held a t = 15 run of 4,962 rows followed by a t = 40 run of 13,233 — on different grids. Read together they manufactured a velocity drop of 0.99 U_j in a single sample, which was time running backwards at the seam and was diagnosed as a physical instability first.
A related symptom is louder and easier to spot: if the probe output format changes between runs, the column count changes partway down the file and numpy.loadtxt refuses it outright ("the number of columns changed from 11 to 18"). That one at least announces itself. The time reset does not.
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:
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
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:
V. A. Titarev, E. F. Toro, Finite-volume WENO schemes for three-dimensional conservation laws, Journal of Computational Physics 201 (1) (2004) 238–260.
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.
Third member of the conduction verification set, after examples/1D_conduction_convergence (Cartesian) and examples/2D_axisym_conduction_convergence (radial, plus the axis cell). Both of those run at p = 0 and therefore have no azimuthal direction at all. This case covers it.
In 3D cylindrical mode (grid_geometry == 3) the third coordinate is the azimuth and dz is in radians, so the azimuthal term carries two metric factors, (1/r^2) d2T/dtheta2: one for the gradient and one for the divergence. m_rhs.fpp divides the conduction flux difference by dz(l) alone, so both factors have to come out of grid_spacing in m_conduction.fpp, which is why that direction uses y_cc(y)**2 * (z_cc(z+1) - z_cc(z)).
Ideal gas, so T = p / ((Gamma - 1) * rho * cv). Setting
rho(r, theta) = RHO0 / (1 + A (r/R)^2 cos(theta))
at uniform pressure gives exactly
T(r, theta) = T0 (1 + A (r/R)^2 cos(theta)), T0 = P0 / ((Gamma - 1) RHO0 cv),
which is single-valued and smooth on the axis. Its cylindrical Laplacian splits into
(1/r) d/dr (r dT/dr) = +4 A T0 cos(theta) / R^2 (radial half) (1/r^2) d2T/dtheta2 = -1 A T0 cos(theta) / R^2 (azimuthal half)
and the sum, 3 A T0 cos(theta) / R^2, is independent of r. At t = 0 the velocity is zero and the pressure uniform, so every Euler flux and every cylindrical geometric source vanishes and the whole right-hand side is conduction:
(rho E(dt) - rho E(0)) / dt = k grad^2 T + O(dt).
The radial half is reproduced exactly by the discretization: the 3D cylindrical r-grid is uniform (m_grid.f90 only half-cells the axis for grid_geometry == 2), T is quadratic in r, and both the two-point face gradient and the two-point average of F are exact for that. Subtracting it therefore isolates the azimuthal half, and compare_analytic.py reports it ring by ring. Dropping either metric factor multiplies the azimuthal half by r^2, which no refinement removes.
Boundaries: periodic in x and in theta, so the azimuthal direction has no boundary error at all. Excluded from the bulk error are cell k = 0, which reads the ghost across the axis (the documented non-converging axis error), and the two outermost cells: k = NR-1 reads the mirrored wall ghost where dT/dr is not actually zero, and k = NR-2 picks up an O(dt) share of that through the later Runge-Kutta stages.
for np_ in 32 64 128 256; do
NP=$np_ ./mfc.sh run examples/3D_cyl_azimuthal_conduction_convergence/case.py -n 1
NP=$np_ ./build/venv/bin/python3 examples/3D_cyl_azimuthal_conduction_convergence/compare_analytic.py
done
NR = 32, NX = 32, exact = 7.5e-01 cos(theta). The azimuthal truncation error of a centered second difference is -dtheta^2/12 of the azimuthal half, i.e. dtheta^2/36 of the total, and nothing else contributes:
| NP | bulk L2 rel err | order | predicted dtheta^2/36 |
|---|---|---|---|
| 32 | 1.070e-03 | – | 1.071e-03 |
| 64 | 2.676e-04 | 2.00 | 2.677e-04 |
| 128 | 6.689e-05 | 2.00 | 6.693e-05 |
| 256 | 1.670e-05 | 2.00 | 1.673e-05 |
Azimuthal half recovered ring by ring, as a fraction of its exact value (NP = 64, so the expected value is 1 - dtheta^2/12 = 0.999197 at every radius):
| r | 0.094 | 0.531 | 1.031 | 1.844 |
|---|---|---|---|---|
| got/exact | 0.999199 | 0.999197 | 0.999197 | 0.999197 |
Flat across a 20x span in r, i.e. a 380x span in r^2. Refining NR at fixed NP leaves the bulk error unchanged (2.676e-04 at NR = 32, 2.656e-04 at NR = 64), confirming that the radial half contributes no error and that the table above is a pure azimuthal measurement.
Without the r^2 in grid_spacing the same NP = 64 run gives a bulk L2 of 3.542e-01 and a ring-by-ring azimuthal ratio of 0.0088, 0.282, 1.063, 3.397 at those four radii – exactly r^2, and independent of resolution.
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. 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: 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).
ib_neighborhood_radius is a count of rank hops. A rank keeps an immersed-boundary patch only while the patch's centroid lies inside its own subdomain grown outward by that many hops (s_get_neighbor_bounds, f_neighborhood_ranks_own_location), and drops it otherwise. So the radius has to be large enough that the body is reachable within that many hops in every direction — and the hop that needs the most is the one stepping across the thinnest rank.
When the radius is not set in the case file, MFC chooses it from the body's half-extent divided by a rank width. The width it used was assembled the wrong way round:
Each rank reports its widest extent, and the minimum is taken over ranks. On a decomposition where every rank is long in one direction and thin in another, the reported width is the long one, the radius comes out too small, and ranks that should have kept the patch drop it.
A thin plate in a long, narrow channel: 20 chords by 8 by 6, with 400 x 50 x 50 cells chosen so MFC's topology search settles on 16 x 2 x 2 at 64 ranks. That is an ordinary shape for a wake, a jet or a channel — the flow direction resolved far more finely than the cross-stream ones — and it makes the ranks strongly anisotropic:
| direction | ranks | extent per rank |
|---|---|---|
| x | 16 | 1.250 |
| y | 2 | 4.000 |
| z | 2 | 3.000 |
The plate is 1.0 x 2.4 x 0.1, so s_get_ib_bound (geometry 9, the cuboid's half-diagonal) returns 1.3010. Crossing that at 1.250 per hop needs ceil(1.1 * 1.3010 / 1.250) = 2 hops. The old width of 4.000 gives ceil(1.1 * 1.3010 / 4.000) = 1.
and read the line MFC prints at start-up:
| printed radius | |
|---|---|
| before | 1 |
| after | 2 |
Both runs complete; the case is 1 M cells and takes a few minutes on two CPU nodes. SUMMARY=1 python3 case.py prints the half-extent, the rank extents and the arithmetic above without running anything.
The same two production grids that motivated the fix, measured from their own lustre_*_cb.dat:
| case | topology | old width | old radius | new width | new radius |
|---|---|---|---|---|---|
| gust encounter, 128 ranks | 16 x 2 x 4 | 2.051 | 1 | 1.052 | 2 |
| flapping wing, 128 ranks | 8 x 4 x 4 | 1.745 | 1 | 1.027 | 2 |
Both pick 1 where 2 is required.
The new width is never larger than the old one, so the chosen radius never decreases: the change can only make the neighbourhood more conservative, at the cost of more hops in the force reduction. Cases that set ib_neighborhood_radius explicitly are untouched, and so is any decomposition whose ranks are close to cubic, where the widest and narrowest extents coincide — which is why a uniform grid with a balanced topology shows no difference.
Reference:
Trojak, W., & Dzanic, T. Positivity-preserving discoutinous spectral element method for compressible multi-species flows. arXiv:2308.02426
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.
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
