|
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 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;geometric, constant_flux, or initial_flux handler.gradient_model.enabled is what would settle it. Wall functions obtain the friction velocity from an iterative root-find, and a matrix-free Jacobian action amplifies that solve's tolerance by 1/h, so a loose inner tolerance shows up as a degraded Jacobian action rather than as an error.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.
The two driven periodic handlers are inside the validated scope. What makes them admissible is that their momentum source is frozen for the whole timestep. That takes two independent gates, both keyed on simCtx->step, because two pieces of state sit between the flux measurement and the applied force:
bulkVelocityCorrection, computed in the handler's PreStep from the field at the start of the step (src/BC_Handlers.c);SimCtx and applied in ComputeDrivenChannelFlowSource() (src/BodyForces.c).Both functions run once per residual evaluation, so gating only the first still leaves the applied force walking toward its target across evaluations. With both gated, every residual evaluation in the Newton solve sees the same body force.
This is worth stating explicitly because it is not automatic. ApplyBoundaryConditions runs the handler PreStep sweep three times per call, and the residual callback calls it on every evaluation - so a handler that recomputed its source in PreStep would let that source drift with the trial vector, silently changing the operator the Newton solve is converging on and invalidating the finite-difference Jacobian action. The per-timestep freeze in PreStep_PeriodicDrivenConstant is what prevents that. The boundary trim (boundaryVelocityCorrection) is deliberately not frozen, but it acts on boundary fluxes rather than on the momentum source; see 5.2 Update cadence, and why it matters.
Because the source is frozen at the start of a step, the flux controller lags the target by one step - first order in dt. Halving dt halves the lag. That is the same lag the Picard solver sees, so results from the two solvers are comparable.
The paired-face checks below still apply: both faces of the driven axis must be PERIODIC and must carry the same handler.
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;PreStep, so a handler carrying per-timestep controller state must guard it against being advanced here; see 1.1 Driven periodic faces;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.
No other finite-difference mode or Jacobian type is implemented; any other value is rejected at validation.
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.
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.Use nonlinear_solver.eisenstat_walker.enabled to select fixed KSP tolerances or PETSc inexact-Newton forcing. The block exposes PETSc versions 1–4 and every EW parameter: initial and maximum relative tolerance, gamma, exponent, safeguard exponent, and safeguard threshold. PETSc may change the effective KSP tolerance at every Newton iteration. See SNESKSPSetUseEW and the SNES manual. Other PETSc SNES/KSP options remain accessible through prefixed petsc_passthrough_options.
Newton–Krylov monitors are enabled under solver_monitoring.momentum (see Configuration Reference: Monitor YAML):
newton_krylov_history -> -mom_nk_pic_monitor: PICurv's own nonlinear and inner-linear iteration histories;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_Linear_History_Block_<b>.log: one row per KSP iteration (step | block | newton | krylov | requested_rtol | reported_residual_norm), including the effective tolerance after any EW update;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 assembles a same-cell 3x3 frozen-coefficient approximation in the current F=-R residual convention. Matrix rows are residual components and columns are the same-cell contravariant velocity components being differentiated; for constraint rows it inserts the exact modern derivative (+1 identity for fixed rows, or +1/-1 for periodic duplicates). Its viscous diagonal carries the same effective viscosity the residual diffuses with, nu + nu_t, using the residual's own face average of the eddy viscosity; omitting the eddy term left the matrix modelling a viscous diagonal smaller than the operator's by the eddy-to-molecular ratio, which on a developed large-eddy simulation is order one or more.
FrozenMomentumJacobian_FaceEddyViscosity() also returns zero on a face lying on a WALL. That branch currently changes no assembled entry, and it does not match the residual. It changes nothing because the eddy viscosity reaches the matrix only through the three diagonal entries, so a row carries only the viscosity of its own axis, and the two coordinates that trigger the branch are exactly the ones ClassifyMomentumRow() classifies as boundary-pinned - those rows take an identity stamp and never assemble a block at all. Both rules restate one boundary fact: the wall-normal flux at a no-slip wall is not an unknown. It does not match the residual because Viscous() substitutes the wall-model eddy viscosity lnu_wall there and falls back to zero only when no wall model is active. Widening the stencil breaks the first of those and exposes the second, because an interior row would then need the eddy viscosity on the wall face to build its off-diagonal. Issue #8 tracks fixing it before any preconditioner with a wider stencil. It still intentionally omits pressure, the eddy-viscosity derivatives with respect to velocity, nonorthogonal viscous cross-couplings, the gradient (Clark) stress, boundary-map derivatives, and body-force derivatives. The eddy viscosity and the Clark stress are omitted for different reasons, and only one of them was a judgement call: the eddy term was restored because it changes the viscous diagonal by the eddy-to-molecular ratio, while the Clark term is higher order and non-diagonal, so representing it would change what kind of matrix this is.
PCPBJACOBI is only the current internal backend mapping; it is not a user-facing numerical model.
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.
| Value | Maps to | Status |
|---|---|---|
frozen_momentum_jacobian | frozen_momentum_jacobian | experimental |
none | none | supported |
Identity. momentum_solver.newton_krylov.preconditioner.model: none (the default) -> -mom_nk_preconditioner_model none.
What it does. Runs the Krylov solve unpreconditioned on the matrix-free Jacobian. No preconditioning matrix is created, assembled, or stored.
When to choose it. The default, and the right starting point: it has no assembly cost, no extra memory, and no approximation of its own to be wrong. Move to frozen_momentum_jacobian only once Krylov iteration counts are demonstrably the bottleneck on your case.
Parameters it owns. None. Setting preconditioner.structure alongside this model is a validation error - none accepts no matrix structure, because there is no matrix.
Interactions. Mutually exclusive with frozen_momentum_jacobian. The Krylov method and its tolerances are configured separately under linear_solver.
Diagnostics. The Krylov iteration count per Newton step is the signal. It is reported by the SNES/KSP monitors described in 7. Monitors and Log Output; a rising count across timesteps is what motivates preconditioning.
Evidence. Integration verified - make unit-newton-krylov exercises this path. Production exercised - turbulent-channel-nk-2026-09-29, retained at 10.1 Retained turbulent-channel campaign (2026-09-29), uses this model on 144 ranks.
Limitations. Iteration counts grow with conditioning, so on stiff or highly stretched grids the unpreconditioned solve can dominate the timestep cost. Supported within the solver's declared scope; one production campaign establishes neither performance scaling nor accuracy for other flows.
Identity. momentum_solver.newton_krylov.preconditioner.model: frozen_momentum_jacobian with preconditioner.structure.type: point_block -> -mom_nk_preconditioner_model frozen_momentum_jacobian and -mom_nk_preconditioner_structure point_block.
What it does. Assembles a separate AIJ matrix holding a same-cell 3x3 frozen-coefficient momentum block in the current residual convention, and uses it as the preconditioning operator for the matrix-free Jacobian. Constraint rows carry the exact current derivative.
When to choose it. When Krylov iteration counts under none are the measured bottleneck. Its benefit depends on the grid, state, timestep, and omitted operator terms, so compare Krylov counts and solve time on the intended case.
Parameters it owns. preconditioner.structure.type, which must be point_block. It is not an independent choice: the model determines it, and any other value is a validation error.
Interactions. Requires structure.type: point_block and is rejected without it. Its viscous diagonal carries the residual's own effective viscosity nu + nu_t, so a turbulence model is represented at the level the diagonal can represent it. A wall model is not: the wall-face eddy viscosity is the one place this matrix does not follow the residual, which is harmless only while the stencil stays at zero (issue #8). The matrix still deliberately omits pressure, the eddy-viscosity derivatives with respect to velocity, nonorthogonal viscous cross-couplings, the gradient (Clark) stress, boundary-map derivatives, and body-force derivatives - so its quality degrades as those terms matter more.
Diagnostics. Krylov iteration counts before and after are the only meaningful diagnostic. PCPBJACOBI appears in PETSc output as the internal backend mapping; it is not a user-facing numerical model and should not be read as one.
Evidence. Integration verified - make unit-newton-krylov covers model/backend/ownership wiring, exact constraint rows, matrix structure and reuse, and serial/MPI application.
Limitations. Experimental, and no performance claim is made. It costs an extra assembled matrix in memory and an assembly per update, and the omitted terms mean it is a same-cell approximation rather than an approximate Jacobian in any global sense.
The Newton–Krylov path has two regression levels and the retained production measurement below:
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.
The repository owner approved promotion of the solver and preconditioner.model: none after this campaign closed the previously recorded production-size gap. tests/tooling/measurement_records.json, record turbulent-channel-nk-2026-09-29, is the durable source for the configuration, measurements, acceptance criterion, provenance hashes, and limitations. Its verdict concerns successful production execution. The DNS comparison below remains exploratory; it does not establish quantitative turbulent-flow accuracy.
The run used commit e37e868821beae3732b97cce8c0c6bc125914b82, PETSc 3.20.3 debug, OpenMPI 4.1.5, and 144 Slurm ranks with a 4 x 4 x 9 decomposition. The grid had 96 x 96 x 256 physical cells in spanwise, wall-normal, and streamwise directions, respectively, spanning 2 pi x 2 x 4 pi, with stretched wall-normal spacing. It used a single block, no-slip walls, periodic spanwise boundaries, constant-flux streamwise forcing, bulk target 1, viscosity 1/2800, and timestep 0.005. Momentum used central convection, matrix-free finite differences, SNES newtonls with backtracking, GMRES, and no preconditioner; pressure used FGMRES and four-level multigrid. LES, wall functions, IBM, and particles were disabled.
The supplied continuation logs cover steps 10,001 through 20,000. Every one of those 10,000 steps reports a converged Newton solve and committed state; maximum logged divergence is 5.26e-11. Logs for startup steps 0 through 10,000 were not available for this audit. The step-20,000 checkpoint contains 1,999 accepted statistics samples with total time weight 49.975, last sample step 19,996, and one restart inside the statistics window. This establishes continuation of the accumulators; it is not a comparison against an uninterrupted run.
The comparison uses the Lee–Moser Re_tau = 182.088 mean profile and velocity fluctuations. Here h = 1; streamwise velocity is the k component and wall-normal velocity is the j component. Plane-averaged checkpoint statistics are compared at equal wall distance y/h. Wall friction is inferred from a no-slip quadratic fit to the first two cell-center means at each wall; the two wall stresses are averaged before taking the square root for u_tau. It is not a recorded solver traction. Skin friction is Cf = 2 u_tau^2 / Ub^2. The reference Cf uses its reported viscosity, friction Reynolds number, and bulk velocity.
| Retained quantity | Channel run | Lee–Moser reference |
|---|---|---|
| Bulk velocity | 0.9999903 | 1 |
| Inferred Re_tau | 202.876 | 182.088 |
| Cf | 0.01049987 | 0.00812326 |
| Peak Reynolds shear, -uv / Ub^2 | 0.00392530 | 0.00295518 |
| Peak Reynolds shear, -uv / u_tau^2 (each flow's own u_tau) | 0.7477 | 0.7276 |
| Centerline velocity / Ub | 1.18153 | 1.16426 |
The run's friction is 29.3% above the reference and its bulk-normalized shear peak is 32.8% above it. Rescaling by each flow's own friction velocity reduces the shear-peak difference to 2.8%; that rescaling does not remove the friction error. The internal total-shear residual has an interior RMS of 1.22% of the inferred wall stress. This balance is a consistency check, not independent DNS validation.
Separate halves reconstructed from the weighted checkpoint accumulators give Re_tau 203.584 then 202.166, volume-averaged temporal TKE 0.00876582 then 0.00899821, and peak -uv/Ub^2 0.00390817 then 0.00394227. Full-window TKE is 0.00911349; variation of the mean between halves adds variance, so it need not lie between the two half-window TKE values. These two blocks do not establish a confidence interval. Spectra were produced from 51 checkpoints, but no DNS spectrum acceptance test was performed.
| Experimental surface exercised | Lifecycle decision from this run |
|---|---|
momentum.newton_krylov, momentum.solver: Newton Krylov, momentum.nk_preconditioner: none | Supported within the declared scope; record the campaign under production evidence |
cluster.scheduling | Remains experimental: Slurm submission and dependent post-processing ran, but live cancellation, graceful cancellation, and accounting reconciliation were not checked |
workspace.asset_lifecycle | Remains experimental: manifest records content-addressed grid/IC materialization through hardlinks; cluster import modes, campaign-scale reuse, reflinks, and remote-backed pruning were not checked |
Statistics, spectra, initial-condition generation, periodic forcing, and restart already have supported records. Asset materialization through hardlinks does not establish that the selectable workspace.input_import_mode: hardlink operation was invoked. The frozen-Jacobian preconditioner, LES, wall functions, storage operations, and parameter sweeps were not exercised by this campaign. No refinement study, solver-strategy comparison, parallel scaling measurement, or fully converged DNS validation is claimed.
| 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