|
PICurv 0.1.0
A Parallel Particle-In-Cell Solver for Curvilinear LES
|
This page documents PICurv's matrix-free Newton–Krylov momentum solver: what it solves, how it is configured, how to read its convergence output, and the residual-purity invariant that makes it work. It is one of two momentum-solution approaches in PICurv; see Momentum Solver Implementations for how it compares to the dual-time Picard–Jameson solver and how to select between them.
The Newton–Krylov solver advances the implicit momentum update by solving the nonlinear momentum residual directly with PETSc's SNES, using a matrix-free (Jacobian-free) Krylov linearization. It is implemented in MomentumSolver_NewtonKrylov (src/momentum_newton_krylov.c) and is selected with strategy.momentum_solver: "Newton Krylov".
Version one is deliberately narrow and validates its inputs up front (MomentumSolver_NewtonKrylov rejects anything outside this set):
--continue); the first solved restart step uses BDF1 because only the checkpoint state is available;Nvert must be fluid everywhere);Within that scope it is a drop-in alternative to the dual-time Picard–Jameson solver and shares the same fractional-step projection, BDF time discretization, boundary system, and pressure solve.
Each physical timestep advances the contravariant velocity Ucont by solving the discrete momentum residual to zero:
\[ F(\mathbf{U}) \;=\; -\,\mathrm{RHS}_{\text{spatial}}(\mathbf{U}) \;+\; \frac{a_0}{\Delta t}\,\mathbf{U} \;-\; (\text{BDF history terms}) \;=\; 0, \]
where RHS_spatial is the convective + viscous + source assembly (ComputeRHS) and the time term is the BDF discretization added by ComputeTotalResidual —
Order selection is centralized in MomentumUsesBDF2()/MomentumBDFCoefficient(), shared with the momentum stability estimate. The Newton solver does not change the residual arithmetic, BDF coefficients, conservation-outlet formulas, or the number of boundary passes; it only changes how the resulting nonlinear system is solved.
The four layers have distinct roles. \(F(U)=0\) is the unchanged nonlinear equation. The Jacobian operator approximates \(dF/dU\) and is the operator GMRES uses for the Newton correction. The preconditioning matrix is a cheaper mathematical approximation to that operator. Finally, the PETSc PC is an internal algorithm that approximately applies the inverse of the supplied preconditioning matrix.
The solver builds a per-step SNES (MomentumSolver_NewtonKrylov):
SNESNEWTONLS** with a backtracking (bt) line search;MatCreateSNESMF, whose action is the finite-difference directional derivative \(J\mathbf{v} \approx [F(\mathbf{X}+h\mathbf{v}) - F(\mathbf{X})]/h\) (MatMFFDComputeJacobian); this is always authoritative;KSPGMRES** uses PCNONE;PCPBJACOBI.Because the Jacobian action is a finite difference of the residual, the residual must be a deterministic function of the trial vector X (see 5. Deterministic Cartesian Seeding (Why It Is Required)). GMRES uses PETSc's default classical Gram–Schmidt orthogonalization; modified Gram–Schmidt is not required and is not enabled.
The Newton loop is: evaluate F(X) → form Krylov solve of J dX = -F → line search along dX → repeat until an SNES convergence test fires. On convergence the solution is committed into Ucont; on failure the entry state is restored (rollback) and the physical step is reported as not converged (simCtx->mom_last_converged).
One residual evaluation (MomentumNewtonKrylov_FormResidual) performs, in order:
X into the global Ucont;Ucont and refresh local lUcont ghosts;Ucat from the current lUcont, finalize periodic Ucat, and refresh lUcat ghosts;F = -Rhs, then a constrained-row pass that replaces every non-independent row (fixed boundary-normal, homogeneous dummy/tangential, and periodic-duplicate rows) with an explicit algebraic equation so the matrix-free operator has no zero Jacobian rows.The three internal boundary passes are unchanged and remain necessary: each pass refreshes the Cartesian state after a boundary correction so the next pass sees a consistent field.
This is the invariant that makes the matrix-free solve correct, and it is easy to break by "simplifying" the residual, so it is documented explicitly.
Every evaluation of F(X) must start from velocity fields derived from that same X. The conservation-outlet handler reads the cell-centered Cartesian velocity lUcat during its first boundary sweep (it measures the uncorrected outlet flux and builds the outlet profile from it). If lUcat were left over from a previous residual or matrix-free evaluation, F(X) would depend on that hidden state — two evaluations at the same X could differ, and the finite-difference Jacobian action would be inconsistent.
The velocity-state relationship is:
So the residual seeds them in exactly this dependency order before the first boundary pass:
Important subtleties, all captured in the source comment above the seed:
Contra2Cart() alone is not sufficient: it rebuilds the interior of the global Ucat but does not refresh lUcat (nor lUcont), and the outlet reads the local ghosted lUcat.ApplyBoundaryConditions() runs after each handler sweep, so it prepares passes two and three — it cannot prepare the very first outlet read of pass one.SynchronizePeriodicCellFields("Ucat") must run before the ghost scatter so periodic duplicate planes are finalized consistently (it is a no-op when no direction is periodic).Removing or shortening this sequence reintroduces a history-dependent residual and invalidates the Newton directions. A permanent regression guards it (Section 10).
Select the solver and (optionally) tune its PETSc controls. Omitted fields keep the defaults established in src/momentum_newton_krylov.c and PETSc.
The point-block alternative replaces the preconditioner block with:
jacobian.type: finite_difference means that the complete deterministic nonlinear residual \(F(U)\) is differentiated numerically. jacobian.finite_difference.mode: matrix_free means PETSc evaluates directional products on demand and does not assemble the Jacobian. Finite difference is the construction type; matrix free is one mode within that type.
The planned second finite-difference mode is colored_sparse, which would assemble a sparse numerical Jacobian using coloring. It is future syntax and is rejected today:
Planned frozen-momentum approximations are a different Jacobian type, not storage modes of finite difference. Their future forms are:
or:
Those forms are also rejected today.
The Jacobian fields map to -mom_nk_jacobian_type finite_difference and -mom_nk_jacobian_fd_mode matrix_free. The preconditioner fields map to application-owned -mom_nk_preconditioner_* selectors. They do not emit a user-selected PETSc PC type. The released linear_solver.preconditioner.type: none spelling is accepted as a deprecated compatibility alias; it conflicts with a non-none new model. Field-by-field mappings and validation rules (nonnegative tolerances and positive iteration/restart counts) are the authoritative configuration reference in Solver Reference, section 4. The complete annotated template is examples/master_template/master_solver.yml.
The planned operator-intent correspondence with legacy selectors is:
| Modern configuration | Legacy operator coverage |
|---|---|
finite_difference / matrix_free + none | top-level -imp 4 |
finite_difference / matrix_free + frozen point block | -imp 5 -imp_type 2 |
finite_difference / colored_sparse | future -imp 5 -imp_type 1 |
frozen momentum / diagonal | future -imp 5 -imp_type 3 |
frozen momentum / full_sparse | future -imp 5 -imp_type 4 |
This table maps numerical operator intent only. It does not reproduce legacy defects, lifecycle, mutable-residual behavior, or historically unknown PETSc runtime options.
Three configuration layers interact, in increasing precedence:
momentum_solver.newton_krylov.*) — the supported surface.petsc_passthrough_options** — raw PETSc options applied last. A raw PC type must match the backend derived by the preconditioning engine or setup fails with an explicit incompatibility error.The tolerances above are a reasonable starting point. Interpretation:
nonlinear_solver.absolute_tolerance stops Newton when the nonlinear residual norm falls below it — the primary physical convergence gate.nonlinear_solver.relative_tolerance stops Newton relative to the initial residual norm.linear_solver.relative_tolerance controls how tightly each inner GMRES solve is converged; a loose 1e-6 inexact-Newton setting is typical and cheap.Newton–Krylov monitors are enabled under solver_monitoring.momentum (see Configuration Reference: Monitor YAML):
newton_krylov_history -> -mom_nk_pic_monitor: PICurv's own per-iteration nonlinear-norm history;snes_monitor -> -mom_nk_snes_monitor, snes_converged_reason -> -mom_nk_snes_converged_reason;ksp_monitor -> -mom_nk_ksp_monitor, ksp_converged_reason -> -mom_nk_ksp_converged_reason.Independently of PETSc monitors, the solver writes structured rank-zero logs into log_dir:
Momentum_Solver_Newton_Krylov_History_Block_<b>.log: one row per Newton iteration (step | block | newton | nonlinear_norm);Momentum_Solver_Newton_Krylov_Summary_Block_<b>.log: one row per physical step (the mathematical Jacobian and preconditioner selections, SNES reason, Newton iterations, residual evaluations, Krylov iterations, initial/final norm, and whether the result was committed or rolled back).A healthy solve on the validated duct case shows the nonlinear norm dropping by several orders of magnitude in about two Newton iterations with accepted line search lambda = 1.
Do not treat all non-convergence the same — the SNES/KSP reason identifies the failure class:
CONVERGED_FNORM_ABS / CONVERGED_FNORM_RELATIVE**: success (absolute or relative nonlinear tolerance met).DIVERGED_MAX_IT**: hit snes_max_it without meeting a tolerance — usually under-resolved inner solves or too tight a nonlinear tolerance for the timestep; loosen nonlinear_solver.relative_tolerance or reduce dt.DIVERGED_LINEAR_SOLVE**: an inner GMRES solve failed to converge — inspect ksp_converged_reason; raise linear_solver.max_iterations or gmres.restart, or loosen linear_solver.relative_tolerance.DIVERGED_LINE_SEARCH**: the backtracking line search could not find a sufficient decrease — typically a poor Newton direction. In this solver that most often means the residual was not deterministic (a broken Cartesian seed, Section 5); it should not occur with the shipped residual.dt.Troubleshooting workflow: enable snes_monitor + snes_converged_reason + ksp_converged_reason, reproduce on a short run, and classify by reason before changing tolerances. If you observe DIVERGED_LINE_SEARCH or non-repeatable nonlinear norms, suspect residual determinism (Section 5) rather than the Krylov settings.
The Jacobian interface owns creation, registration, update, naming, and cleanup of the finite-difference/matrix-free operator. The preconditioner model interface only describes a matrix structure and inserts interior physical coefficients. The common engine owns matrix creation/preallocation, repeated zeroing and assembly, constraint and periodic rows, PETSc backend selection, alias/ownership tracking, and cleanup.
The optional frozen-momentum/point-block matrix is a separate AIJ matrix. For physical rows it reproduces the audited same-cell 3x3 block from the reachable legacy mode-2 approximation, with its sign reversed to match the modern residual convention; for constraint rows it inserts the exact modern derivative (+1 identity for fixed rows, or +1/-1 for periodic duplicates). It intentionally omits pressure, LES/RANS viscosity values and derivatives, nonorthogonal viscous cross-couplings, boundary-map and IBM derivatives, and body-force derivatives.
The legacy face inverse Jacobians are arithmetic averages of neighboring cell inverse Jacobians, and its transverse metric terms average four squared metric-vector norms. The modern directional fields IAj/JAj/KAj directly invert separately constructed face determinants, while ICsi/IEta/IZet, JCsi/JEta/JZet, and KCsi/KEta/KZet come from directional face-center constructions. They are not algebraically equivalent to the legacy expressions on a general curvilinear grid: in particular, an average of squared norms is not the squared norm of an average. The explicit legacy formulas are therefore retained here. Directional fields remain potential inputs to a future, separately specified modern preconditioner model.
The point-block coefficients have an independent test-only legacy transcription covering nonuniform metrics and velocities, every block entry, both BDF coefficients, constraint rows, and one- and multi-rank assembly. No performance claim is made. PCPBJACOBI is only the current internal backend mapping; it is not a user-facing numerical model and is not a historically proven legacy setting.
Future additions are localized as follows: add a Jacobian type/mode beside MomentumNewtonJacobian_Create/Update; add a coefficient provider through MomentumPreconditionerModelOps; add matrix metadata through MomentumPreconditionerDescription; and add a validated structure-to-PETSc backend mapping in MomentumPreconditionerEngine_Create. No placeholder modes are exposed before their implementations exist.
The Newton–Krylov path is covered at two levels:
unit-newton-krylov, part of make check): constraint-row Jacobian structure, matrix-free vs direct differencing, preconditioning-engine model/backend/ownership wiring, small solve/rollback, and residual repeatability. The conservation-outlet conditioned-row derivative test doubles as a seed-removal detector: removing the deterministic Cartesian seed makes that row's self-derivative revert to the decoupled artifact and the test fails.make unit-momentum-newton-boundary-fixedpoint, one and four ranks): on the production-sized straight duct it advances a real physical step 1, then verifies (i) residual purity at the step-2 state (immediate and after real MFFD products), (ii) a complete step-2 solve with the default classical Gram–Schmidt, (iii) that the converged three-pass solution also zeros the 24-pass outlet residual, and (iv) clean pressure projection.Validated behavior on that case: convergence in about two Newton iterations from the true projected step-1 state, identical results with classical and modified Gram–Schmidt, and divergence-free projection, on both one and four ranks.
| Aspect | Newton–Krylov | Dual-Time Picard–Jameson |
|---|---|---|
| Linearization | true Newton (matrix-free Jv) | Picard fixed-point / pseudo-time smoothing |
| Inner solve | PETSc SNES + GMRES | staged Jameson RK pseudo-time |
| Main controls | SNES/KSP tolerances, GMRES restart | pseudo-CFL, pseudo-iterations |
| Maturity | newer, narrow validated scope (Section 1) | established, broadly exercised |
| Failure surface | SNES/KSP convergence reasons | pseudo-CFL rollback / rejection |
See Dual-Time Picard Jameson RK Momentum Solver for the Picard–Jameson solver and Momentum Solver Implementations for selection guidance.
1.9.8