Pressure solvers

The discrete pressure Poisson equation

$$Lp=WMF(u)$$

enforces divergence freeness at every Runge-Kutta stage (see Spatial and temporal discretization). Since this is the only globally coupled (and thus most expensive) part of a time step, choosing an appropriate solver matters:

• psolver_spectral: FFT-based solver. The fastest option, but requires a uniform grid with periodic boundary conditions in all directions.

• psolver_transform: FFT/DCT-based solver for uniform grids with periodic and Dirichlet boundary conditions.

• psolver_tridiagonal: FFT/tri-diagonal solver for channel-like setups (one wall-bounded direction, periodic otherwise). The periodic directions must be uniform, but the wall-normal direction may be stretched. Works on the GPU.

• psolver_direct: sparse direct solver (factorize once, solve every step). Works for all grids and boundary conditions, but the factorization only supports Float64 on the CPU, and memory usage can be prohibitive for large 3D problems. On CUDA arrays, loading CUDSS.jl makes this solver use a GPU factorization instead.

• psolver_cg: matrix-free conjugate gradient solver. Works on all backends and precisions; the go-to fallback when the direct solver does not apply.

• psolver_cg_AMGX: conjugate gradient solver with algebraic multigrid preconditioning from NVIDIA AMGX (requires CUDA and AMGX.jl).

The default (default_psolver) selects the spectral solver for uniform periodic setups, the tri-diagonal solver for channel-like setups, and the direct solver otherwise.

Preconditioning the matrix-free CG solver

psolver_cg is matrix-free: it only applies the Laplacian stencil (laplacian!), so it runs on any backend. It currently uses a diagonal (Jacobi) preconditioner, which parallelizes trivially but does not improve the mesh-dependent conditioning. GPU-friendly upgrades, roughly in order of implementation effort:

• Polynomial preconditioning (Chebyshev or truncated Neumann series): applies the same Laplacian kernel a few times per iteration; no setup phase, no extra storage, fully matrix-free.

• Fast-solver preconditioning: use an FFT-based solver (psolver_spectral, psolver_transform, or psolver_tridiagonal) for a nearby constant-coefficient problem as the preconditioner. One FFT round-trip per iteration, and typically a mesh-independent iteration count when the grid is a smooth deformation of a uniform one.

• Geometric multigrid: matrix-free smoothers (damped Jacobi/Chebyshev) plus coarsening kernels; mesh-independent convergence for all supported grids and boundary conditions.

• Algebraic multigrid: needs the assembled matrix; on CUDA this already exists as psolver_cg_AMGX.

Incomplete factorizations (IC(0)/ILU) are not a good fit: the triangular solves are inherently sequential and perform poorly on GPUs.

API

IncompressibleNavierStokes.amgx_setup Function

Initializes AMGX, all needed objects are returned in a named tuple. Needs to be followed by amgx_close after use.

Becomes available using AMGX.

source

IncompressibleNavierStokes.close_amgx Function

Close all objects created by amgx_setup.

Becomes available using AMGX.

source

IncompressibleNavierStokes.default_psolver Method

default_psolver(
    setup
) -> Union{IncompressibleNavierStokes.var"#psolve!#123", IncompressibleNavierStokes.var"#psolve!#162", IncompressibleNavierStokes.var"#psolve!#98"{Bool}}

Get default Poisson solver from setup.

source

IncompressibleNavierStokes.poisson! Method

poisson!(psolver, f) -> Any

Solve the Poisson equation for the pressure (in-place version).

source

IncompressibleNavierStokes.poisson Method

poisson(psolver, f) -> Any

Solve the Poisson equation for the pressure with right hand side f. For periodic and no-slip BC, the sum of f should be zero.

Differentiable version.

source

IncompressibleNavierStokes.project! Method

project!(u, setup; psolver, p)

Project velocity field onto divergence-free space (in-place version).

source

IncompressibleNavierStokes.project Method

project(u, setup; psolver)

Project velocity field onto divergence-free space (differentiable version).

source

IncompressibleNavierStokes.psolver_cg Method

psolver_cg(
    setup;
    abstol,
    reltol,
    maxiter,
    preconditioner
) -> IncompressibleNavierStokes.var"#psolve!#105"{_A, _B, _C, IncompressibleNavierStokes.var"#laplace_diag#103"{ndrange, workgroupsize}} where {_A, _B, _C, ndrange, workgroupsize}

Conjugate gradients iterative Poisson solver.

source

IncompressibleNavierStokes.psolver_cg_AMGX Function

Poisson solver using conjugate gradient method from AMGX.

Becomes available using AMGX.

source

IncompressibleNavierStokes.psolver_direct Method

psolver_direct(
    setup
) -> IncompressibleNavierStokes.var"#psolve!#98"{Bool}

Create direct Poisson solver using an appropriate matrix decomposition.

source

IncompressibleNavierStokes.psolver_spectral Method

psolver_spectral(
    setup
) -> IncompressibleNavierStokes.var"#psolve!#123"

Create spectral Poisson solver from setup.

source

IncompressibleNavierStokes.psolver_transform Method

psolver_transform(
    setup
) -> IncompressibleNavierStokes.var"#psolve!#118"{_A, NamedTuple{(:uhat, :w, :winv, :perm, :perminv), var"#s182"}} where {_A, var"#s182"<:NTuple{5, Any}}

Create FFT/DCT Poisson solver from setup. This solver does FFT in periodic directions, and DCT in Dirichlet directions. Only works on uniform grids, with periodic/dirichlet BC. If there is only Periodic BC, then psolver_spectral is faster, and uses half as much memory.

Warning: This does transform one dimension at a time. With Float32 precision, this can lead to large errors.

source

IncompressibleNavierStokes.psolver_tridiagonal Method

psolver_tridiagonal(
    setup;
    dir
) -> IncompressibleNavierStokes.var"#psolve!#162"

FFT/tri-diagonal Poisson solver for channel-like setups: one wall-bounded direction (DirichletBC on both sides) and periodic boundary conditions in the other direction(s). FFTs in the periodic directions decouple the Poisson equation into an independent tri-diagonal system along the wall-normal direction for each Fourier mode. The systems are solved in a batched Thomas-algorithm kernel (one mode per thread), so the solver works on the GPU as well.

The periodic directions require uniform grid spacing, but the wall-normal direction may be arbitrarily stretched. This makes this solver a fast direct method for channel flows on stretched wall-normal grids, where psolver_transform (fully uniform grids) does not apply.

The wall-normal direction is inferred from the boundary conditions. It can also be chosen explicitly with dir (it must still be wall-bounded).

source