[wilfonba]

This pull request adds an all-Mach projection solver path to MFC. The method supports the 5-equation model with viscosity, surface tension, body forces, and heat conduction, modeling fluids with ideal or stiffened gas equations of state. Figures and comments are contained in "details" blocks to keep the comment from getting too long.

Method

The numerical method can be summarized as

1. Compute $\mathbf{q}^*$ from $\mathbf{q}^n$ by advancing the advective portion of the five-equation system via explicit transport and find the advected pressure

$$\begin{equation*} p^a = p^n + \Delta t \mathbf{u}^n \cdot \nabla p^n \end{equation*}$$

2. Solve for the pressure $p^{n+1}$ at $t^{n+1}$ via

$$\begin{equation*} p^{n+1} + \rho^n \left(c^2\right)^n\Delta t^2 \nabla \cdot \left( \frac{\nabla p^{n+1}}{\nabla \rho^{n+1}}\right) = p^a - \rho^n \left(c^2)\right)^n\nabla \cdot \mathbf{u}^* \end{equation*}$$

using multigrid-preconditioned flexible conjugate gradients.
3. Correct the momentum and energy equations with

$$\begin{align*} \left(\rho\mathbf{u}\right)^{n+1} &= (\rho \mathbf{u})^* - \Delta t \nabla p^{n+1} \\\\\ \rho E^{n+1} &= \rho E^* - \Delta t \nabla \cdot \left(p^{n+1} \mathbf{u}^{n+1}\right) \end{align*}$$

The viscous, capillary, and thermal stresses are treated explicitly. Conduction enters the pressure equation as $\partial_t p += -\nabla \cdot \mathbf{q}/\Gamma$. Future PRs will incorporate split time-stepping of the diffusion terms to take advantage of the advective time step allowed by this method. I am currently testing immersed boundary and chemistry features in a development branch and will include them in future pull requests.

Validation

(2D) Couette flow at Re = 100, Mach = 0.01

Simulation ran on 1 A100 GPU using the `examples/2D_couette_projection/case.py` example. The error $\max|u - u_\text{exact}|$ falls from 8.0e-4 at t = 1 to 2.3e-5 at t = 30.

(2D) Poiseuille flow at Re = 100, Mach = 0.01

Simulation ran on 1 A100 GPU using the `examples/2D_poiseuille_projection/case.py` example. The error $\max|u - u_\text{exact}|$ is ~2.4e-4 at all times.

(2D) Lid driven cavity at Re = 5000

Simulation ran on 1 A100 GPU using the `examples/2D_lid_driven_cavity_projection/case.py` example. The centerline velocities at t = 200 match those of Ghia et al. closely.

(2D) Oscillating ellipsoidal bubble

Simulation ran on 1 A100 GPU using the `examples/2D_oscillating_bubble_projection/case.py` example. The surface-tension-induced oscillation frequency matches linear theory.

(2D) Laplace pressure jump

Simulation ran on 1 A100 GPU using a copy of the `examples/2D_laplace_pressure_jump/case.py` example with the projection method enabled. The pressure jump error converges with grid resolution, and spurious currents remain small.

(2D) Rayleigh-Benard convection

Simulation ran on 1 A100 GPU using the `examples/2D_rayleigh_benard_projection/case.py` example. The steady-state Nusselt number at the hot and cold plates at Ra = 4e4 and Pr = 7 matches that of Waleffe et al. (Phys. Fluids 27, 051702, 2015, Fig. 1: ~3.7 at L/H = 2).

(2D) Rising bubble

Simulations ran on 1 A100 GPU using the `examples/2D_rising_bubble/case.py` example. The simulation uses body forces, viscosity, and surface tension with the projection solver, and comparisons with established data from other solvers are used. Case 1 is a density ratio of 10 and Case 2 is a density ratio of 1000.

(3D) Taylor-Green Vortex

Simulations performed on 1-4 A100 GPUs using the `examples/3D_TaylorGreenVortex_projection/case.py` example. The kinetic energy, kinetic energy dissipation rates, and entropy converge to the spectral DNS results of Dairay, Lamballais, Laizet & Vassilicos, J. Comput. Phys. 337, 2017.

Convergence

(1D) 5-eqn

Convergence is measured using a modified copy of the `examples/1D_euler_convergence/case.py` example. The acoustic wave reduces spatial and temporal accuracy to second and first order, respectively.

(1D) 5-eqn + Viscosity

Convergence is measured using a modified copy of the `examples/1D_euler_convergence/case.py` example. The acoustic wave reduces spatial and temporal accuracy to second and first order, respectively. Viscosity damps the acoustic waves, so the temporal convergence is flatter than in the previous case.

(1D) 5-eqn + Conduction

Convergence is measured using a modified copy of the `examples/1D_conduction_convergence/case.py` example.

(2D) 5-eqn + Viscosity

Convergence is measured using a modified copy of the `examples/2D_TaylorGreenVortex/case.py` example.

(2D) 5-eqn + Conduction

Convergence is measured using a 2D extension of the `examples/1D_conduction_convergence/case.py` example.

(1D) Control - Projection solver with no acoustics

When a case contains pure advection, div u = 0, and the projection step doesn't have to correct anything. In this case, we recover 5th-order spatial and 3rd-order temporal convergence.

Conservation

(2D, inviscid) Mass, momentum, and energy conservation

The following figure shows mass, momentum, and energy conservation for a 2D inviscid TGV simulation run with the explicit (dashed) and projection (solid) solvers. The projection solver conserves mass slightly better than the explicit solver despite the larger time-step. Momentum conservation is nearly equivalent for both. The projection solver loses more energy than the explicit solver, but the loss converges to zero at 3rd order with resolution.

(3D, viscous) Mass, momentum, and energy conservation

Mass, momentum, and energy conservation are evaluated using the 3D Taylor green vortex simulation. Momentum conservation is floating-point exact. Mass conservation is close to floating-point exact, and converges with resolution. Energy is not conserved in this test because viscous heating is not turned into thermal energy. Instead, it matches the falling kinetic energy expected in a TGV.

Scaling and performance

Weak Scaling (single fluid ideal gas)

Weak scaling is performed using a tiled 3D Taylor-Green Vortex at Mach 0.01 with 300^3 grid cells per GCD on OLCF Frontier. The code was compiled with AMDFlang AFAR drop 23.2.1, and GPU-direct MPI was used where possible. The figure shows, from top to bottom, the wall time per step, the average number of pressure solve iterations, and the speedup over an explicit solve. Different lines and colors correspond to different iterative solver parameters; the "Proj. Trunc. V-Cycle" line (purple + triangles) uses the defaults and performs best. The pressure solve iteration count is constant across problem sizes for this approach, and speedup exceeds 40x through 9400 nodes.

Strong Scaling (single fluid ideal gas) (coming soon)

Strong scaling was performed on OLCF Frontier with AMD Flang AFAR 24.3.0 using a 1260^3. Strong scaling isn't as important with this solver because you can run a larger problem size per device with far fewer time steps and still get the answer faster than the explicit solve.

Projection vs. Explicit Cost (single fluid ideal gas)

The following table compares the per-step and wall-clock time costs of a Mach 0.01 TGV for the explicit and projection solvers on a single A100 GPU at increasing problem sizes. The projection solver has a higher per-step cost in all cases, but it takes far fewer time steps to reach the same simulation time. Speedups reach up to 51x for large problem sizes.

N	solver	dt	steps to t = 20 tC	s / step	run time	grind (ns/gp/eq/rhs)
100³	projection	0.000834	637	0.082	0.9 min	4.55
100³	explicit	8.26e-06	64,299	0.0274	29.4 min	1.52
100³	explicit / projection			0.33×	34×	0.33×
200³	projection	0.000417	1,274	0.362	7.7 min	2.52
200³	explicit	4.13e-06	128,598	0.177	6.3 h	1.23
200³	explicit / projection			0.49×	49×	0.49×
300³	projection	0.000278	1,910	1.12	35.8 min	2.31
300³	explicit	2.75e-06	192,896	0.581	31.1 h	1.20
300³	explicit / projection			0.52×	52×	0.52×
350³	projection	0.000238	2,229	1.84	1.1 h	2.39
350³	explicit	2.36e-06	225,046	0.923	2.4 days	1.20
350³	explicit / projection			0.50×	51×	0.50×

Projection vs. explicit cost (two fluids ideal+stiffened gas) (coming soon)

The following table compares the per-step and wall-clock time costs of a contrived example in which a water droplet moves diagonally in a tri-periodic box and rotates about its center axis. This provides a near-constant Mach number and enough dynamics for the pressure projection solve to affect runtime. The projection solver takes time steps ~140x as large as the explicit solver, but is limited by the capillary cfl and not the advective CFL.

N	solver	dt	steps per crossing	s / step	run time	grind (ns/gp/eq/rhs)
100³	projection	3.81e-06	210	0.0848	0.3 min	3.14
100³	explicit	2.74e-08	29,192	0.04	19.5 min	1.48
100³	explicit / projection			0.47×	66×	0.47×
200³	projection	1.9e-06	420	0.421	2.9 min	1.95
200³	explicit	1.37e-08	58,383	0.251	4.1 h	1.16
200³	explicit / projection			0.60×	83×	0.60×
300³	projection	1.27e-06	630	1.37	14.4 min	1.88
300³	explicit	9.14e-09	87,574	0.819	19.9 h	1.12
300³	explicit / projection			0.60×	83×	0.60×

[github-actions]

Claude Code Review

Head SHA: b03afff

Files changed:

• 29
• .github/CODEOWNERS
• cmake/MFCTargets.cmake
• src/common/include/macros.fpp
• src/common/m_constants.fpp
• src/common/m_derived_types.fpp
• src/common/m_mpi_common.fpp
• src/simulation/m_data_output.fpp
• src/simulation/m_multigrid.fpp
• src/simulation/m_multigrid_exchange.fpp
• src/simulation/m_projection.fpp

Findings:

• medium toolchain/mfc/params/definitions.py: The check_projection_simulation hunk in case_validator.py drops proj_tol <= 0 ("proj_tol must be positive"). The replacement is "proj_tol": {"min": 0}. schema.py:103 rejects only value < min, so proj_tol = 0 now passes validation, where it was rejected before. A zero relative tolerance cannot be reached by the PCG solve, and the solve then runs to proj_max_iters every stage. Fix: keep the positivity check in case_validator.py. Alternatively, if the constraint schema has an exclusive bound, use it. The proj_max_acfl move to {"min": 0} is equivalent. Confidence: fairly high.
• medium toolchain/mfc/test/cases.py: The hunk removes the if len(dimInfo[0]) > 1: guard above cases.append(define_case_d(stack, "2 MPI Ranks", ...)) in alter_projection, and the append stays at the 12-space indent of the capillary block. The "2 MPI Ranks" case now runs only if the enclosing if of the capillary cases is true. If that condition is not len(dimInfo[0]) > 1, the 2-rank projection test silently disappears for some dimensions. It could also run for 1D. CLAUDE.md says the ppn=2 test is what catches missing-broadcast and halo bugs in this PR. Fix: confirm the enclosing condition, or restore the explicit guard. Confidence: low to medium; I could only see the hunk context.
• medium src/simulation/m_time_steppers.fpp: s_propagate_immersed_boundaries is rewritten, not just refactored. The velocity blend flips from c1*step_vel + c2*vel to c1*vel + c2*step_vel. The position and angle updates now use the pre-blend stage velocity. Previously they used the already-updated velocity. The two-way branch changes from moving_ibm == 2 to >= 2. The new blend ordering matches the q_cons blend (c1*current + c2*stored), so it is probably the intended RK-consistency fix. It still changes results for every moving-IB case, and the PR has nothing on it: no mention in the diff beyond the code and no test or golden update tied to it. Fix: split it into its own PR, or add a note and verification. Also confirm that nothing uses moving_ibm > 2. Confidence: medium that this is a deliberate change; the risk is the missing verification.
• medium src/simulation/m_muscl.fpp: The hunk adds if (slopeL < 0._wp) slope = -slope to SUPERBEE. The recent master commit d9fa909 ("Restore the slope sign in the SUPERBEE limiter") appears to do the same thing. If this branch is rebased on top of that commit, the sign is applied twice, which cancels the fix. Fix: rebase and confirm the line appears only once. Confidence: low; this depends on the merge base.

Checked in step 4: the multigrid level/exchange hunks (buffer layout, send/recv offsets, req(12) sizing, self-seam source layer, bottom-level cutoff arithmetic), the ACC_SETUP_mg_levels macro, and the s_projection_apply signature change against its caller. I found no defects there.

[github-actions]

Lines of Code

File	Lines	Diff
src/simulation/m_projection.fpp	860	+860
src/simulation/m_multigrid.fpp	654	+654
src/simulation/m_multigrid_bottom.fpp	275	+275
src/simulation/m_multigrid_exchange.fpp	194	+194
src/simulation/m_multigrid_levels.fpp	81	+81
src/simulation/m_projection_helpers.fpp	52	+52
src/simulation/m_rhs.fpp	2010	+42
src/simulation/m_data_output.fpp	1605	+27
src/simulation/m_time_steppers.fpp	934	+24
src/common/m_derived_types.fpp	505	+22
src/common/m_constants.fpp	105	+16
src/common/include/macros.fpp	143	+15
src/simulation/m_global_parameters.fpp	809	+14
src/common/m_mpi_common.fpp	1504	+10
src/simulation/m_sim_helpers.fpp	265	+10
src/simulation/m_start_up.fpp	1274	+3

Directory	Lines	Diff
common	10589	+63
simulation	30767	+2236
total	49889	+2299