Under pressure!

Fast pressure calculations in the Deft package

Jos van Kan

j.vankan@math.tudelft.nl

Delft Institute of Applied mathematics (DIAM)

Deft package

Delft Flow and Transport

• Navier Stokes equations for incompressible flow on general domains

• Offshoot of ISNaS (Information System Navier Stokes)

Design decisions

• Finite Volume Method

• Rectangular blocks of curvilinear coordinates • Staggered grid

• Time-dependent algorithm

• Pressure correction for the incompressibility condition

Navier Stokes equations

\frac{∂u}{∂t} + ∇ · (uu^{T}) + ∇p = ∇ · T + f
∇ · u = 0

• u velocities

• p pressure

• T stress tensor

• f body forces like gravity Variations in density are not taken into account. p and T are scaled quantities.

Method of lines

\begin{aligned}\frac{du_{h}}{dt} + Gp_{h} & = N(u_{h}) + L(u_{h}) + f_{h} \\ Du_{h} & = 0 \end{aligned}

Time Integration

Implicit time integration, like for instance Crank-Nicolson

u_{h}^{n+1}+ΔtGp_{h}^{n+1/2} =
u_{h}^{n}+\frac{1}{2}Δt (N(u_{h}^{n+1}) + L(u_{h}^{n+1}) + f_{h}^{n+1}) +

\frac{1}{2}Δt (N(u_{h}^{n}) + L(u_{h}^{n}) + f_{h}^{n})

Du_{h}^{n+1} = 0

Matrix Structure

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.7/??

Pressure Correction

u_{h}^{∗} + ΔtGp_{h}^{n−1/2} =
u_{h}^{n}+\frac{1}{2}Δt (N(u_{h}^{∗}) + L(u_{h}^{∗}) + f_{h}^{n+1}) +

\frac{1}{2}Δt (N(u_{h}^{n}) + L(u_{h}^{n}) + f_{h}^{n})

u_{h}^{n+1} − u_{h}^{∗}+ΔtGΔp = 0
−Du_{h}^{∗}+ΔtDGΔp = 0
p^{n+1/2} = p^{n−1/2} + Δp

Pressure Correction

u_{h}^{∗} + ΔtGp_{h}^{n−1/2} =
u_{h}^{n}+\frac{1}{2}Δt (N(u_{h}^{∗}) + L(u_{h}^{∗}) + f_{h}^{n+1}) +

\frac{1}{2}Δt (N(u_{h}^{n}) + L(u_{h}^{n}) + f_{h}^{n})

u_{h}^{n+1} − u_{h}^{∗}+ΔtGΔp = 0
−Du_{h}^{∗}+ΔtDGΔp = 0
p^{n+1/2} = p^{n−1/2} + Δp

p Matrix structure 2D

Example: n × n block, N = n^{2}. • Tridiagonal block matrix of tridiagonal matrices • Bandwidth: O(\sqrt{N})

• Flops LU decomposition: O(N^{2}). (One time only on fixed domains)

• Flops LU backsubstitution: O(N^{3/2})

• Flops matrix vector multiplication: O(N).

p Matrix structure (3D)

Example: n × n × n block, N = n^{3}.

• Tridiagonal block matrix of tridiagonal blockmatrices of tridiagonal matrices

• Bandwidth O(N^{2/3})

• Flops LU decomposition: O(N^{7/3}). (Only once on fixed domains)

• Flops LU backsubstitution O(N^{5/3})

• Flops matrix vector multiplication: O(N)

A Classic: Defect Correction

Solve Ax = b Presets: x^{0} = 0, r^{0} = b − Ax^{0} = b, k = 0 while ∥r^{k}∥_{∞} > ε∥b∥_{∞} do

Solve P c^{k} = r^{k} {P is preconditioner} x^{k+1} = x^{k} + c^{k} r^{k+1} = r^{k} − Ac^{k} k = k + 1

end while

Each iteration requires O(N) flops.

Preconditioners

Classic With A = D − L − U • Jacobi: P = D

• Gauss-Seidel: P = D − L

• Successive overrelaxation: P = (D/ω − L) Modern With A = LU

• Incomplete LU (ILU):

A = L̃Ũ + E, P = L̃Ũ. L̃ and Ũ sparse, usually the same sparsity pattern as A.

• Incomplete Block LU (IBLU).

DC Error Reduction

r^{k+1} = r^{k} − Ac^{k} = (I − AP^{−1})r^{k}
A^{−1}b − x^{k+1} = A^{−1}r^{k+1} = ε^{k+1}
ε^{k+1} = A^{−1}(I − AP^{−1})Aε^{k} = (I − P^{−1}A)ε^{k}

Reduction governed by spectral radius of (I − AP^{−1}). For the Laplacian:

• Jacobi and Gauss-Seidel: 1 − O(h^{2})

• SOR with optimal ω and a whole slew of other conditions: 1 − O(h).

.

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.13/??

Effectiveness of DC

How many iterations to gain a decimal digit?

ε^{n+k} = λ_{1}^{k}ε^{n}
λ_{1}^{k} = 0.1
k log λ_{1} = − log 10
k = −\frac{2.3}{log λ_{1}}
k = O(\frac{1}{1 − λ_{1}})

Jacobi and Gauss-Seidel O(h^{−2}) iterations. SOR O(h^{−1}) iterations. In 2D:

• Jacobi and GS O(N^{2}) flops, worse than LU • SOR (O(N^{3/2}) flops, order equal to LU

In 3D:

• Jacobi and GS O(N^{5/3}) flops, order equal to LU • SOR O(N^{4/3}) flops, better than LU

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.14/??

Effectiveness of DC

Jacobi and Gauss-Seidel O(h^{−2}) iterations. SOR O(h^{−1}) iterations. In 2D:

• Jacobi and GS O(N^{2}) flops, worse than LU • SOR (O(N^{3/2}) flops, order equal to LU

In 3D:

• Jacobi and GS O(N^{5/3}) flops, order equal to LU • SOR O(N^{4/3}) flops, better than LU

Convergence properties

Convergence properties

Damped Jacobi, 10 iterations Chebyshev10, λ_{0} = 0.1 Chebyshev10, λ_{0} = 0.03 Chebyshev10, λ_{0} = 0.01

Gradient Methods

• "Best" polynomial on the fine structure of the spectrum

• Gradient Methods are always better than Defect Correction

• Irregular convergence behaviour

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.16/??

Multigrid

Defect Correction: very effective on a part of the spectrum.

The eigenspace of the spectral interval (0.2, 1) is virtually reduced to 0 in a few iterations.

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.17/??

1D example

Consider −\frac{d^{2}u}{dx^{2}} = f, u(0) = u(1) = 0 Discretize into N intervals): Au = f Eigenvalues of 1 − P^{−1}A are

λ_{k} = 1 − sin^{2} \frac{kπ}{2N} , k = 1, …, N − 1.

Corresponding eigenvectors

v_{kj} = sin \frac{kjπ}{N}, k, j = 1, …, N − 1.

Eigenvalues close to 1 correspond to smooth eigenvectors, also for the Laplacian in 2 and 3D.

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.18/??

Rough and smooth spectrum

• Rough part of the spectrum: defect correction, smoother in MG speak.

• Smooth part of the spectrum: solve problem on coarser grid and interpolate. Coarse grid correction in MG speak.

Restriction and prolongation

Fine grid correction: A_{h}c_{h} = r_{h} Coarse grid correction A_{H}c_{H} = r_{H} Transfer operators:

• P_{hH}: prolongation from coarse to fine grid. Interpolation usually.

• R_{Hh}: restriction from fine grid to coarse grid.

R = P^{T} in symmetric problems. The coarse grid correction becomes: R_{Hh}A_{h}P_{hH}c_{H} = R_{Hh}r_{h}

Two Grid Algorithm

Presets: u_{h}^{0}, r_{h}^{0} = f_{h} − Au_{h}^{0} u_{h}^{prs} = S(u_{h}^{0}, b, A, n_{0}){Presmoothing} r_{H} = R_{Hh}r_{h} Solve A_{H}c_{H} = r_{H} u_{h}^{cgc} = u_{h}^{prs} + P_{hH}c_{H} {Coarse Grid Correction} u_{h}^{pos} = S(u_{h}^{cgc}, b, A, n_{1}){Postsmoothing}

Two Grid Algorithm

Presets: u_{h}^{0}, r_{h}^{0} = f_{h} − Au_{h}^{0} u_{h}^{prs} = S(u_{h}^{0}, b, A, n_{0}){Presmoothing} r_{H} = R_{Hh}r_{h} Solve A_{H}c_{H} = r_{H} u_{h}^{cgc} = u_{h}^{prs} + P_{hH}c_{H} {Coarse Grid Correction} u_{h}^{pos} = S(u_{h}^{cgc}, b, A, n_{1}){Postsmoothing}

Two Grid Algorithm

Multi Grid Algorithm

Presets: u_{h}^{0}, r_{h}^{0} = f_{h} − Au_{h}^{0} u_{h}^{prs} = S(u_{h}^{0}, b, A, n_{0}){Presmoothing} r_{H} = R_{Hh}r_{h} Solve A_{H}c_{H} = r_{H} u_{h}^{cgc} = u_{h}^{prs} + P_{hH}c_{H} {Coarse Grid Correction} u_{h}^{pos} = S(u_{h}^{cgc}, b, A, n_{1}){Postsmoothing} Require: A_{ℓ+1} = R_{ℓ+1,ℓ}A_{ℓ}P_{ℓ,ℓ+1} have been

calculated on all levels MGRecursive (A_{ℓ}, r_{ℓ}, c_{ℓ}, ℓ) if ℓ < p then

c_{ℓ} = S(0, r_{ℓ}, A_{ℓ}, n_{0}){Presmoothing} r_{ℓ+1} = R_{ℓ+1,ℓ}(r_{ℓ} − A_{ℓ}c_{ℓ}) {Calculate coarse grid residual} call MGRecursive (A_{ℓ+1}, r_{ℓ+1}, c_{ℓ+1}, ℓ + 1) c_{ℓ} = c_{ℓ} + P_{ℓ,ℓ+1}c_{ℓ+1} {Coarse grid correction} c_{ℓ} = S(c_{ℓ}, r_{ℓ}, A_{ℓ}, n_{1}) {Postsmoothing}

else Solve A_{p}c_{p} = r_{p}{Direct solution on coarsest level}

end if

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.22/??

Require: A_{ℓ+1} = R_{ℓ+1,ℓ}A_{ℓ}P_{ℓ,ℓ+1} have been calculated on all levels MGRecursive (A_{ℓ}, r_{ℓ}, c_{ℓ}, ℓ) if ℓ < p then

c_{ℓ} = S(0, r_{ℓ}, A_{ℓ}, n_{0}){Presmoothing} r_{ℓ+1} = R_{ℓ+1,ℓ}(r_{ℓ} − A_{ℓ}c_{ℓ}) {Calculate coarse grid residual} call MGRecursive (A_{ℓ+1}, r_{ℓ+1}, c_{ℓ+1}, ℓ + 1) c_{ℓ} = c_{ℓ} + P_{ℓ,ℓ+1}c_{ℓ+1} {Coarse grid correction} c_{ℓ} = S(c_{ℓ}, r_{ℓ}, A_{ℓ}, n_{1}) {Postsmoothing}

else Solve A_{p}c_{p} = r_{p}{Direct solution on coarsest level}

end if

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.22/??

Require: A_{ℓ+1} = R_{ℓ+1,ℓ}A_{ℓ}P_{ℓ,ℓ+1} have been calculated on all levels MGRecursive (A_{ℓ}, r_{ℓ}, c_{ℓ}, ℓ) if ℓ < p then

c_{ℓ} = S(0, r_{ℓ}, A_{ℓ}, n_{0}){Presmoothing} r_{ℓ+1} = R_{ℓ+1,ℓ}(r_{ℓ} − A_{ℓ}c_{ℓ}) {Calculate coarse grid residual} call MGRecursive (A_{ℓ+1}, r_{ℓ+1}, c_{ℓ+1}, ℓ + 1) c_{ℓ} = c_{ℓ} + P_{ℓ,ℓ+1}c_{ℓ+1} {Coarse grid correction} c_{ℓ} = S(c_{ℓ}, r_{ℓ}, A_{ℓ}, n_{1}) {Postsmoothing}

else Solve A_{p}c_{p} = r_{p}{Direct solution on coarsest level}

end if

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.22/??

Require: A_{ℓ+1} = R_{ℓ+1,ℓ}A_{ℓ}P_{ℓ,ℓ+1} have been calculated on all levels MGRecursive (A_{ℓ}, r_{ℓ}, c_{ℓ}, ℓ) if ℓ < p then

c_{ℓ} = S(0, r_{ℓ}, A_{ℓ}, n_{0}){Presmoothing} r_{ℓ+1} = R_{ℓ+1,ℓ}(r_{ℓ} − A_{ℓ}c_{ℓ}) {Calculate coarse grid residual} call MGRecursive (A_{ℓ+1}, r_{ℓ+1}, c_{ℓ+1}, ℓ + 1) c_{ℓ} = c_{ℓ} + P_{ℓ,ℓ+1}c_{ℓ+1} {Coarse grid correction} c_{ℓ} = S(c_{ℓ}, r_{ℓ}, A_{ℓ}, n_{1}) {Postsmoothing}

else Solve A_{p}c_{p} = r_{p}{Direct solution on coarsest level}

end if

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.22/??

The MG miracle

• Spectrum of I − P^{−1}A bounded away from 1 uniformly in h.

• Number of iterations does not depend on h.

• The workload is theoretically O(N) flops.

But how big is the multiplicative constant going to be?

Robust Blackbox

Wishlist:

• Good smoother under various circumstances (anisotropy, stretched and skew cells)

• Arbitrary number of points in either direction

Smoothers tested

• (Alternating) damped line Jacobi, 1 or 2 postsmoothing steps, no presmoothing

Smoothers tested

• (Alternating) damped line Jacobi, 1 or 2 postsmoothing steps, no presmoothing

• Incomplete Block LU decomposition, 1 postsmoothing step, no presmoothing

Comparison

Jos van Kan (DIAM), CASA talk February 15, 2006 – p.26/??

Recursion

• Line Jacobi is recursive per line and can be massively parallelized, especially in 3D.

• What about IBLU? Classic IBLU is fully recursive.

Divide and Conquer

The inversion of an n × n tridiagonal matrix can be executed in^{2} log n non recursive steps.

(I + LD^{−1} + UD^{−1})(D − L − U) = D − LD^{−1}U − UD^{−1}L − LD^{−1}L − UD^{−1}U

= D_{1} − L_{1} − U_{1}

Bandwith is doubled in this operation.

Incomplete Block Div and Conq

• Use the same formula, interpreted als blocks • Use incomplete versions of LD^{−1}U etc.

• D, L and U are (block)diagonals, consisting of tridiagonal blocks.

• We need 7 diagonals of D^{−1}

• Calculate from the productform of the inverse, keeping only 7 diagonals in^{2} log n steps

Recursion uses O(^{2} log n) steps as claimed.