|
PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
|
|
Loading...
Searching...
No Matches
|
PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
|
|
This page describes the pressure-correction solve path used by projection in PICurv.
The correction solve enforces incompressibility through:
\[ \nabla^2 \phi = \frac{1}{\Delta t}\nabla\cdot\mathbf{u}^*. \]
In code terms:
The operator and the projection evaluate the gradient flux through each face with the same stencil, including its one-sided forms beside non-periodic boundaries and solid cells, so the Laplacian the solver inverts is exactly the divergence of the gradient the projection subtracts. Non-periodic faces are homogeneous Neumann. The constant null space of that problem is removed on every level by a callback attached to each level's operator.
The integral of the right-hand side is reported as Poisson Source Imbalance in Continuity_Metrics.log. It is the net volume flux into the domain carried by the momentum step's velocity before the correction, and must stay near zero for the all-Neumann equation to have a solution. It covers every boundary face, unlike that log's Net Flux column, which sums only the inlet and outlet faces.
On its first call for a block, PoissonSolver_Multigrid does the following:
KSP (options prefix ps_) with a multiplicative V-cycle PCMG,and stores the result in the finest level's context. Every later step reuses it: only the right-hand side is formed and the finest-level system solved for Phi. The operator depends only on the grid metrics, the solid field and the boundary types, all fixed for a run. Rebuilding the coarse factorization every step had cost 60-65% of each step on a 0.52M-cell wall-resolved LES duct before the solver was kept. The solver is created at the first step of every process, so --continue and --restart-from with a different rank count build it for the new layout.
After the Poisson solve:
The implementation that rebuilt the solver every step, with its dormant immersed-boundary flux corrections and multi-region null space, is available at commit 53ba654 and earlier.
From solver.yml via picurv_cli/core.py:
poisson_solver.method -> -ps_ksp_type; fgmres (default) or cg. The multigrid preconditioner is not a single fixed linear operator, which only a flexible Krylov method tolerates. Under gmres, lgmres and bcgs - left- or right-preconditioned - the Krylov residual fell to 1e-12 while the true residual stalled near 1e-3, and the projection left a divergence of 2e-4 in a duct; validation refuses all three. cg emits -ps_ksp_norm_type unpreconditioned so it stops on the true residual; fgmres and cg reproduce each other to 1e-14 (poisson-options-2026-09-21). Any other KSP type needs petsc_passthrough_options and is unverified.poisson_solver.absolute_tolerance -> -ps_ksp_atolpoisson_solver.relative_tolerance -> -ps_ksp_rtolpoisson_solver.max_iterations -> -ps_ksp_max_itpoisson_solver.gmres.restart -> -ps_ksp_gmres_restartpoisson_solver.preconditioner.type -> -ps_pc_type; only multigrid is supported todaypoisson_solver.multigrid.levels -> -mg_levelpoisson_solver.multigrid.pre_sweeps -> -mg_pre_itpoisson_solver.multigrid.post_sweeps -> -mg_post_itpoisson_solver.multigrid.semi_coarsening.{i,j,k} -> -mg_i_semi, -mg_j_semi, -mg_k_semipoisson_solver.multigrid.level_solvers.level_N.method -> -ps_mg_levels_N_ksp_typepoisson_solver.multigrid.level_solvers.level_N.preconditioner -> -ps_mg_levels_N_pc_typepoisson_solver.multigrid.level_solvers.level_N.max_it, .rtol, .atol -> -ps_mg_levels_N_ksp_max_it, _ksp_rtol, _ksp_atol (the level_0 forms use the -ps_mg_coarse_ prefix). Before 2026-09-18 these were emitted without the ksp_ prefix, which PETSc left unused, so they had no effect.poisson_solver.multigrid.cycle and .mode are accepted structured keys; current supported values are v and multiplicativepressure_solver remains a legacy alias for poisson_solverpetsc_passthrough_options -> advanced PETSc flags not exposed as structured YAMLFinal option parsing happens in function CreateSimulationContext during context creation.
MG level numbering follows PETSc/PICurv convention: level_0 is the coarsest grid. Larger level numbers are progressively finer. The default MG level preconditioner is block Jacobi (bjacobi) when not specified. Each smoothed level runs pre_sweeps iterations before the coarse-grid correction and post_sweeps after it. When the two counts differ, the post-smoother is a separate PETSc solver: it starts as a copy of the configured pre-smoother (method, preconditioner, tolerances) and reads any further options under -ps_mg_levels_N_up_. A per-level max_it in level_solvers sets the pre-smoother's count only; the post-smoother keeps post_sweeps unless -ps_mg_levels_N_up_ksp_max_it is passed through. Equal counts keep one shared smoother.
Common MG-level preconditioner notes:
jacobi and sor are simple smoother PCs. SOR can use -ps_mg_levels_N_pc_sor_omega <positive-real> through passthrough when needed.ilu and lu are factor PCs. Useful passthrough knobs include -ps_mg_levels_N_pc_factor_levels <nonnegative-integer> and -ps_mg_levels_N_pc_factor_shift_amount <nonnegative-real>.bjacobi owns nested block solves; inspect exact nested PETSc prefixes with -ps_ksp_view before tuning sub-KSP/sub-PC options for a specific PETSc build.Current implementation includes:
monitor.yml -> solver_monitoring.poisson.pic_true_residual.A pressure solve that cannot be completed stops the run. If PETSc reports a non-finite residual (DIVERGED_NANORINF) or a preconditioner it could not build (DIVERGED_PC_FAILED), the solver aborts at that step with the KSP reason, because the projection would otherwise proceed on an unsolved Phi. Stopping at max_iterations is this solve's normal mode and is not reported; any other divergence prints a warning and continues. One confirmed cause of DIVERGED_PC_FAILED is a hierarchy coarsened too far: a duct 9 nodes across with levels: 3 fails at step 1, while the same grid with levels: 2 runs cleanly. The node count at the coarsest level does not predict it on its own - a 5-node grid coarsened to 3 nodes with levels: 2 also runs - so reduce the level count when this reason appears.
If pressure solve quality degrades, check first:
A V-cycle attacks the error at two different scales, with two different tools.
On every level above the coarsest, a smoother runs a fixed number of relaxation sweeps. Relaxation is very good at removing error components whose wavelength is comparable to the local mesh spacing, and almost useless against error that is smooth on that mesh. That is the whole design: each level strips out the high-frequency error it can see, restricts what is left to a coarser mesh where the remaining error looks high-frequency again, and repeats.
At the base of the cycle sits the coarse solve. By construction it receives the error that every smoother above it was blind to: the smooth, long-wavelength, global component. Nothing below it will get another chance at that error, so the coarse solve is expected to remove it essentially exactly, in one shot.
These are categorically different jobs, and PICurv's level naming actively hides the distinction:
| YAML key | PETSc option prefix | Role |
|---|---|---|
level_solvers.level_0 | -ps_mg_coarse_ | Coarse solve at the base of the V-cycle |
level_solvers.level_1 .. level_N | -ps_mg_levels_N_ | Smoothers, coarse to fine |
level_0 looks like just another entry in an evenly spaced list. It is not. It is the one entry that is not a smoother, and configuring it as though it were is the root of the defect described next.
Multigrid here is a preconditioner, not a standalone solver. The outer Krylov method's convergence theory assumes the preconditioner is a fixed linear operator: feeding it the same vector must always return the same vector.
A Krylov method violates that. Its iteration builds a subspace tailored to the vector it was given, and its stopping test fires after however many iterations that particular vector needs. Two different inputs get two different numbers of inner iterations, so the map from input to output is not linear and not even fixed. Put a Krylov method at level_0 and the entire multigrid preconditioner becomes a nonlinear operator.
Pinning ksp_max_it does not fix this. A fixed iteration count still builds a Krylov space out of the input vector, so the operator still depends on its input in a nonlinear way. This is exactly why the classical smoothers - Richardson, Jacobi, Chebyshev, SOR - are the right tools for the smoothing levels: each is a fixed linear operator, applied a fixed number of times.
The outer method PICurv uses, fgmres, is flexible: it tolerates a varying preconditioner well enough to keep constructing a solution. What it cannot do is keep its Arnoldi recurrence consistent with the true residual. FGMRES uses right preconditioning, so under a fixed preconditioner its recurrence-tracked residual equals the true residual b - Ax. Under a varying one, the two quantities separate - and the convergence test reads the tracked one. The solver reports convergence against a number that no longer describes the constraint it exists to enforce.
| Coarse-grid unknowns | Recommended level_0 | Why |
|---|---|---|
| up to ~1e4 | {method: preonly, preconditioner: redundant} | Every rank forms and factors the whole coarse operator with LU. Exact, fixed, linear, no iteration to vary. Cheap because the grid is tiny. |
| larger | {method: preonly, preconditioner: telescope} | Redundant LU on every rank stops being cheap. PCTELESCOPE moves the coarse solve onto a subset of ranks and does a direct solve there, still fixed and linear. |
| last resort | a Krylov method | Only when the coarse grid is genuinely too large to factor. Then set tolerances against the true residual, keep pic_true_residual on, and treat the reported convergence as unverified until you have checked the two norms agree. |
PICurv logs a startup warning when a Krylov ksp_type is configured at level_0. It stays a warning rather than an error because the last-resort case above is legitimate at large scale.
Two constraints bound multigrid.levels from opposite ends.
From below, aim for a coarsest grid of roughly 1e3 to 1e4 unknowns. Smaller wastes a level; larger makes the replicated LU factor expensive on every rank. Because each level removes a factor of eight from a 3-D grid, this scales logarithmically: a grid of ~5M cells wants 5 levels, not 4.
From above, coarsenability. Each level halves an axis as IM -> (IM+1)/2, so IM must stay odd at every level for the coarsening to be exact. IM is the node count, one more than the cell count that case.yml and a grid.gen config state, so an exact ladder needs an even cell count: 128 cells, not 129. Declaring mg_levels in a grid.gen config makes the generator refuse a count that cannot coarsen that far. The chain runs IM_fine = 2 * IM_coarse - 1, giving usable ladders such as
5 -> 9 -> 17 -> 33 -> 65 -> 129 -> 257 4 -> 7 -> 13 -> 25 -> 49 -> 97
An even count still runs, but logs
and proceeds on a slightly misaligned coarse grid. The MPI bound described in the next section applies on top of both of these.
redundant scales indefinitely provided levels are added as the grid is refined. Keeping the level count fixed while refining grows the coarse grid with the fine one, and the replicated LU factor eventually dominates. Adding a level as the grid grows keeps the coarse grid, and therefore the coarse solve, roughly constant.
richardson + bjacobi at level_0 is a fixed linear operator and so is correct, but it is not an exact coarse solve. As the coarse grid grows, Richardson leaves more and more of the smooth error behind and the V-cycle stops being mesh-independent: iteration counts creep up with resolution. It is a reasonable stopgap, not a scalable answer.
DualKSPMonitor (src/logging.c) prints both numbers to <run.runtime_logs>/Poisson_Solver_Convergence_History_Block_0.log when monitor.yml -> solver_monitoring.poisson.pic_true_residual is on:
Unprecond Norm** is PETSc's rnorm, carried by the Krylov recurrence. This is what the convergence test reads.True Norm** is an explicit KSPBuildResidual recomputation of b - Ax.Under a correct configuration the two agree to many digits at every iteration. When they separate, the convergence test is passing on a fiction.
Worked example. A 32x32x192 curved, wall-clustered grid at Re = 40,000 on 8 MPI ranks. Identical case and solver files; only level_solvers.level_0 differs:
level_0 setting | Tracked residual | Recomputed b-Ax | Max divergence | Iterations |
|---|---|---|---|---|
{fgmres, bjacobi} | 4.80e-10 | 1.52e-05 | 1.02e-08 | 14-16 |
{richardson, bjacobi} | 1.67898e-09 | 1.67898e-09 | 5.18e-12 | 11 |
{preonly, redundant} | 8.62828e-10 | 8.62848e-10 | 1.97e-12 | 10-11 |
Read the first row carefully. The tracked residual is the smallest of the three, which is why the defect survived: by the number the solver reports, the Krylov coarse solve looked best. The true residual is five orders of magnitude larger, and the physical consequence follows directly - maximum divergence is 1e-08 rather than 1e-12. The incompressibility constraint the Poisson solve exists to enforce was being violated by six orders of magnitude, silently. The fixed-operator settings are simultaneously more accurate and cheaper, converging in fewer iterations.
The failure was rank-dependent: 6 and 10 ranks were clean, 4 and 8 were broken, with the same configuration. That follows from bjacobi, whose block structure is inherited from the DMDA decomposition, so the coarse preconditioner - and hence how nonlinear the coarse solve behaves - changes with the rank layout. A configuration that looks fine on your development rank count can be wrong on the production one.
Diagnostic sequence:
pic_true_residual and compare the two norms in the Poisson log.level_solvers.level_0 first.<run.runtime_logs>/Continuity_Metrics.log: a decoupled residual shows up as a max divergence orders of magnitude above the solver tolerance you asked for.tests/smoke/run_driven_periodic_regression.sh asserts this invariant at 4 and 10 ranks - two counts that previously disagreed.
multigrid.levels cannot be chosen independently of the rank layout. Every level's DMDA must leave each rank at least stencil_width nodes along every axis, and grid setup requests a stencil width of 3 whenever any axis is periodic (2 otherwise). Coarsening halves the node count per level, so a deep hierarchy spread over many ranks eventually starves the coarsest grid and PETSc aborts during DM creation, before the first timestep:
With levels: L and N cells along an axis, the coarsest grid holds N / 2^(L-1) cells and one more node. Grid setup gives each DMDA one point more than its node count, so the coarsest DMDA holds M = N / 2^(L-1) + 2 points (each level keeps (nodes + 1) / 2 nodes, so this is exact when N divides by 2^(L-1)). PETSc distributes those over P ranks so the smallest rank holds floor(M / P). The constraint is:
floor( (N / 2^(L-1) + 2) / P ) >= stencil_width
Both sides were checked on 2026-10-09 on a 25 x 25 x 97 wall-bounded channel at three levels: the coarsest axis holds 8 points, and 4 ranks on it ran while 5 aborted with Local x-width of domain x 1 is smaller than stencil width s 2.
Worked maxima for a triply periodic box (stencil width 3):
| Cells per axis | Ranks per axis | Total ranks | Max levels |
|---|---|---|---|
| 64 | 2 | 8 | 5 |
| 64 | 4 | 64 | 3 |
| 128 | 4 | 64 | 4 |
| 128 | 8 | 512 | 3 |
| 192 | 6 | 216 | 4 |
Two consequences worth planning around:
grid.da_processors_x/y/z, picurv run checks this before launch, including under --dry-run, for file, grid_gen and programmatic grids, and names the largest rank count each axis can take. When the layout is left to PETSc, the decomposition is not known until runtime, so the first evidence is the aborted job. Set the layout explicitly for any run where it matters.Reducing levels is the usual fix; also delete the now-unused level_solvers.level_N entry for the level that no longer exists. Shallower hierarchies converge more slowly but stay usable. Measure the cost in <run.runtime_logs>/Poisson_Solver_Convergence_History_Block_0.log: if iteration counts grow sharply, prefer lowering the rank count and restoring the deeper hierarchy.
make unit-poisson-rhs covers the module directly:
operator-and-projection-share-one-face-gradient fills every face metric with random, non-orthogonal values, places a solid cell, projects a zero flux with a random Phi, and requires the right-hand side of the result to equal -A Phi on every fluid row to 1e-12. That identity holds for arbitrary metrics only if the operator and the projection use the same face gradient, including its one-sided forms.poisson-solver-multigrid-projects-to-divergence-free perturbs interior face fluxes on a 17-cubed three-level hierarchy, solves and projects, and checks the result is divergence-free.poisson-solver-multigrid-reuses-its-solver checks that a second step reuses the same solver and operator, converges, and adds one header to the convergence log.poisson-solver-multigrid-honours-distinct-sweeps checks the separate pre- and post-smoothing counts.poisson-null-space-removes-the-interior-mean exercises the attached null space.poisson-solver-multigrid-refuses-an-overcoarsened-hierarchy checks that a hierarchy coarsened past what a 9-cubed grid supports stops with the fatal Poisson error rather than projecting on an unsolved Phi.The periodic stencil branches are exercised end to end rather than by a unit test.
Every user-selectable option was run end to end by poisson-options-2026-09-21: one and two levels, full coarsening, one and three sweeps, Chebyshev and Jacobi smoothers, a Krylov coarse solve, a per-level iteration cap, and cg each reproduced the baseline velocity to 7e-15 and gauge-free pressure to 8e-14. The same measurement found gmres, lgmres and bcgs converging their Krylov residual to 1e-12 while the true residual stalled near 1e-3, with left or right preconditioning; those methods are now refused (see 3. YAML Mapping and PETSc Options), and the per-level max_it, rtol and atol keys, which were emitted without PETSc's ksp_ prefix and silently ignored, now reach the level solvers.
The rewritten solver was put through the same exercise on Grace by poisson-option-matrix-2026-10-10: 39 variants, each changing one setting from a baseline. On the laminar square duct (28 variants) these were both methods, gmres.restart 5 and 60, one to five levels, equal and unequal pre/post sweeps including an -ps_mg_levels_N_up_ override, every semi-coarsening axis, Chebyshev, SOR, GMRES, ILU and capped smoothers, Krylov and LU coarse solves, the shipped tolerances, all four monitoring flags, and 1 and 48 ranks. Every variant gave the same pressure gradient to the printed digits and a velocity within 1.4e-9 of the baseline, with divergence at most 5e-11 (8e-8 with the shipped tolerances). Both methods also ran on the curved bent channel, and eight settings on a triply periodic 32^3 decaying-turbulence box, a pure-Neumann problem with a nonzero source every step, agreed with their baselines to 3e-15. With the matrix's 1e-12 tolerances, 17 duct variants hit the 300-iteration cap near steady state without losing agreement, so its iteration counts do not measure cost at production tolerances.
On a turbulent case, hom02-poisson-rewrite-2026-10-10 reran two of the HOM02 decaying isotropic turbulence variants (no model and constant Smagorinsky, 64^3, 452 steps) with the rewritten solver: every decay and spectrum metric against Wray's 512^3 DNS matched the earlier runs to the printed digits, the Poisson solve fell from 1.89 to 0.16 s per step, and each step took 2.4-2.6x less time.
End to end, make smoke-driven-periodic asserts at 4 and 10 ranks that the multigrid coarse solve keeps tracked and true residuals within 1e-4 of each other until both fall below 1e-10 of the step's initial residual (below that the two drift apart in round-off even with an exact coarse solve; a row whose true residual stays high is always compared) and the maximum divergence below 1e-11. Its 32^3 Cartesian fixture does not reproduce the Krylov-level_0 failure above: there the outer FGMRES tolerates even a coarse solve truncated at two iterations, so the check guards the working configuration rather than detecting that defect. The recorded measurement duct-poiseuille-picard-2026-09-18 (see Capability Evidence Matrix) reproduced the analytic axial pressure gradient of laminar square-duct flow at second order on three grids, with divergence at most 5.5e-8.
The rewrite that builds the solver once per run is covered by poisson-persistent-solver-2026-10-10. On unoptimized builds it gave checkpoints and convergence logs bitwise identical to the previous solver, including through --continue at a different rank count and --restart-from. It reproduced the duct measurement above and the Humphrey bend comparison to their printed digits, and it ran the 0.52M-cell wall-resolved LES duct 3.0x faster per step, the Poisson solve 7.0x faster, with the coarse factorization done three times in a 200-step run instead of 600. On the optimized cluster build the two solvers agree to 1e-12 in velocity and a constant pressure offset rather than bitwise.
1.9.8