Вы находитесь на странице: 1из 148

# 5mm.

## INF5620: Numerical Methods for Partial

Differential Equations
Hans Petter Langtangen

## Simula Research Laboratory, and

Dept. of Informatics, Univ. of Oslo

Nonlinear PDEs

Examples

## −u′′ (x) = f (u), u(0) = uL , u(1) = uR ,

−(α(u)u′ )′ = 0, u(0) = uL , u(1) = uR
∂u
−∇ · [α(u)∇u] = g(x), with u or − α B.C.
∂n

Discretization methods:
standard finite difference methods
standard finite element methods
the group finite element method
We get nonlinear algebraic equations
Solution method: iterate over linear equations

## Nonlinear PDEs – p.3/148

Nonlinear discrete equations; FDM

## Finite differences for −u′′ = f (u):

1
− 2 (ui−1 − 2ui + ui+1 ) = f (ui )
h

## Finite differences for (α(u)u′ )′ = 0:

1
(α(ui+1/2 )(ui+1 − ui ) − α(ui−1/2 )(ui − ui−1 )) = 0
h2

1
2
([α(ui+1 ) + α(ui )](ui+1 − ui ) − [α(ui ) + α(ui−1 )](ui − ui−1 )) = 0
2h
⇒ nonlinear system of algebraic equations

## F (u) = 0 or A(u)u = b Nonlinear PDEs – p.4/148

Nonlinear discrete equations; FEM

## Finite elements for −u′′ = f (u) with u(0) = u(1) = 0

Pn
Galerkin approach: find u = k=1 uk ϕk (x) ∈ V such that
Z 1 Z 1
u′ v ′ dx = f (u)vdx ∀v ∈ V
0 0

## Left-hand side is easy to assemble: v = ϕi

Z 1 X
1
− (ui−1 − 2ui + ui+1 ) = f( uk ϕk (x))ϕi dx
h 0 k

P P
We write u = k uk ϕk instead of u = k ck ϕk since
ck = u(xk ) = uk and uk is used in the finite difference schemes

## Nonlinear PDEs – p.5/148

Nonlinearities in the FEM

Note that X
f( ϕk (x)uk )
k

is a complicated function of u0 , . . . , uN
F.ex.: f (u) = u2
Z !2
1 X
ϕk u k ϕi dx
0 k

## gives rise to a difference representation

h 2 2 2

ui−1 + 2ui (ui−1 + ui+1 ) + 6ui + ui+1
12

## (compare with f (ui ) = u2i in FDM!)

R
Must use numerical integration in general to evaluate f (u)ϕi dx

## Nonlinear PDEs – p.6/148

The group finite element method

## The group finite element method:

X n
X
f (u) = f ( uj ϕj (x)) ≈ f (uj )ϕj
j j=1

R1 R1P
Resulting term: 0
f (u)ϕi dx = 0 j ϕi ϕj f (uj ) gives

X Z 1 
ϕk ϕi dx f (uj )
j 0

## This integral andPformulation also arise from approxmating some

function by u = j uj ϕj
We can write the term as M f (u), where M has rows consisting
of h/6(1, 4, 1), and row i in M f (u) becomes

h
(f (ui−1 ) + 4f (ui ) + f (ui+1 ))
6 Nonlinear PDEs – p.7/148
FEM for a nonlinear coefficient

We now look at

## Using a finite element method (fine exercise!) results in an integral

Z 1 X
α( uk ϕk )ϕ′i ϕ′j dx
0 k

## ⇒ complicated to integrate by hand or symbolically

Linear P1 elements and trapezoidal rule (do it!):

1 1
(α(ui ) + α(ui+1 ))(ui+1 − ui ) − (α(ui−1 ) + α(ui ))(ui − ui−1 ) = 0
2 2

## Nonlinear PDEs – p.8/148

Nonlinear algebraic equations

## FEM/FDM for nonlinear PDEs gives

nonlinear algebraic equations:
(α(u)u′ )′ = 0 ⇒ A(u)u = b
−u′′ = f (u) ⇒ Au = b(u)
In general a nonlinear PDE gives

F (u) = 0

or

F0 (u0 , . . . , uN ) = 0
...
FN (u0 , . . . , uN ) = 0

## Nonlinear PDEs – p.9/148

Solving nonlinear algebraic eqs.

Have
A(u)u − b = 0, Au − b(u) = 0, F (u) = 0
Idea: solve nonlinear problem as a sequence of linear
subproblems
Must perform some kind of linearization
Iterative method: guess u0 , solve linear problems for u1 , u2 , . . .
and hope that
lim uq = u
q→∞

## Nonlinear PDEs – p.10/148

Picard iteration (1)

## Model problem: A(u)u = b

Simple iteration scheme:

A(uq )uq+1 = b, q = 0, 1, . . .

## Must provide (good) guess u0

Termination:
||uq+1 − uq || ≤ ǫu
or using the residual (expensive, req. new A(uq+1 )!)

## ||b − A(uq+1 )uq+1 || ≤ ǫr

Relative criteria:
||uq+1 − uq || ≤ ǫu ||uq ||
or (more expensive)

## ||b − A(uq+1 )uq+1 || ≤ ǫr ||b − A(u0 )u0 || Nonlinear PDEs – p.11/148

Picard iteration (2)

## Model problem: Au = b(u)

Simple iteration scheme:

Auq+1 = b(uq ), q = 0, 1, . . .

Relaxation:

## (may improve convergence, avoids too large steps)

This method is also called Successive substitutions

## Nonlinear PDEs – p.12/148

Newton’s method for scalar equations

## The Newton method for f (x) = 0, x ∈ IR

Given an approximation xq
Approximate f by a linear function at xq :

## q+1 q q+1 q f (xq )

M (x ;x ) = 0 ⇒ x =x − ′ q
f (x )

## Nonlinear PDEs – p.13/148

Newton’s method for systems of equations

## Multi-dimensional Taylor-series expansion:

M (u; uq ) = F (uq ) + J (u − uq ), J ≡ ∇F

∂Fi
Ji,j =
∂uj
Iteration no. q:
solve linear system J (uq )(δu)q+1 = −F (uq )
update: uq+1 = uq + (δu)q+1
Can use relaxation: uq+1 = uq + ω(δu)q+1

## Nonlinear PDEs – p.14/148

The Jacobian matrix; FDM (1)

## Model equation: u′′ = −f (u)

Scheme:
1
Fi ≡ 2
(ui−1 − 2ui + ui+1 ) + f (ui ) = 0
h

∂Fi
Ji,j =
∂uj

Only

## ∂Fi ∂Fi ∂Fi

Ji,i−1 = , Ji,i = , Ji,i+1 = 6= 0
∂ui−1 ∂ui ∂ui+1

⇒ Jacobian is tridiagonal
Nonlinear PDEs – p.15/148
The Jacobian matrix; FDM (2)

1
Fi ≡ 2 (ui−1 − 2ui + ui+1 ) − f (ui ) = 0
h
Derivation:
∂Fi 1
Ji,i−1 = = 2
∂ui−1 h
∂Fi 1
Ji,i+1 = = 2
∂ui+1 h
∂Fi 2
Ji,i = = − 2 + f ′ (ui )
∂ui h
Must form the Jacobian matrix J in each iteration and solve

J δuq+1 = −F (uq )

## and then update

uq+1 = uq + ωδuq+1
Nonlinear PDEs – p.16/148
The Jacobian matrix; FEM

## −u′′ = f (u) on Ω = (0, 1) with u(0) = u(1) = 0 and FEM

Z " #
1 X X
Fi ≡ ϕ′i ϕ′k uk − f( us ϕs )ϕi dx = 0
0 k s

∂Fi
First term of the Jacobian Ji,j = ∂uj :

Z X Z Z 1
∂ 1 1
∂ X ′ ′
ϕ′i ϕ′k uk dx = ϕi ϕk uk dx = ϕ′i ϕ′j dx
∂uj 0 0 ∂uj 0
k k

Second term:
Z 1 X Z 1 X
∂ ′
− f( us ϕs )ϕi dx = − f( us ϕs )ϕj ϕi dx
∂uj 0 s 0 s

P
because when u = s us ϕs ,
∂ ′ ∂u ′ ∂
P ′
∂uj f (u) = f (u) ∂uj = f (u) ∂uj s u s ϕ s = f (u)ϕj Nonlinear PDEs – p.17/148
A 2D/3D transient nonlinear PDE (1)

## PDE for heat conduction in a solid where the conduction depends

on the temperature u:

∂u
̺C = ∇ · [κ(u)∇u]
∂t

## (f.ex. u = g on the boundary and u = I at t = 0)

Stable Backward Euler FDM in time:

un − un−1
= ∇ · [α(un )∇un ]
∆t

with α = κ/(̺C)
n
P n
Next step: Galerkin formulation, where u = j j ϕj is the
u
unknown and un−1 is just a known function

## Nonlinear PDEs – p.18/148

A 2D/3D transient nonlinear PDE (2)

## FEM gives nonlinear algebraic equations:

Fi (un0 , . . . , unN ) = 0, i = 0, . . . , N

where
Z
 n n−1
 n n

Fi ≡ u −u v + ∆tα(u )∇u · ∇v dΩ

## Nonlinear PDEs – p.19/148

A 2D/3D transient nonlinear PDE (3)

Picard iteration:
Use “old” un,q in α(un ) term, solve linear problem for un,q+1 ,
q = 0, 1, . . .
Z
Ai,j = (ϕi ϕj + ∆tα(un,q )∇ϕi · ∇ϕj ) dΩ

Z
bi = un−1 ϕi dΩ

## Newton’s method: need Jacobian,

∂Fi
Ji,j =
∂unj
Z
Ji,j = (ϕi ϕj + ∆t(α′ (un,q )ϕj ∇un,q · ∇ϕi + α(un,q )∇ϕi · ∇ϕj )) dΩ

## Nonlinear PDEs – p.20/148

Iteration methods at the PDE level

## Consider −u′′ = f (u)

Could introduce Picard iteration at the PDE level:

d2 q+1
− 2u = f (uq ), q = 0, 1, . . .
dx

## ⇒ linear problem for uq+1

A PDE-level Newton method can also be formulated
(see the HPL book for details)
We get identical results for our model problem
Time-dependent problems: first use finite differences in time, then
use an iteration method (Picard or Newton) at the time-discrete
PDE level

## Nonlinear PDEs – p.21/148

Continuation methods

## Challenging nonlinear PDE:

∇ · (||∇u||q ∇u) = 0

## For q = 0 this problem is simple

Idea: solve a sequence of problems, starting with q = 0, and
increase q towards a target value
Sequence of PDEs:

## with = 0 < q0 < q1 < q2 < · · · < qm = q

Start guess for ur is ur−1
(the solution of a “simpler” problem)
CFD: The Reynolds number is often the continuation parameter q

Exercise 1

## Derive the nonlinear algebraic equations for the problem

 
d du
α(u) = 0 on (0, 1), u(0) = 0, u(1) = 1,
dx dx

## using a finite difference method and the Galerkin method with P1

finite elements and the Trapezoidal rule for approximating
integrals.

Exercise 2

## For the problem in Exercise 1, use the group finite element

method with P1 elements and the Trapezoidal rule for integrals
and show that the resulting equations coincide with those
obtained in Exercise 1.

Exercise 3

## For the problem in Exercises 1 and 2, identify the F vector in the

system F = 0 of nonlinear equations. Derive the Jacobian J for a
general α(u). Then write the expressions for the Jacobian when
α(u) = α0 + α1 u + α2 u2 .

Exercise 4

## Explain why discretization of nonlinear differential equations by

finite difference and finite element methods normally leads to a
Jacobian with the same sparsity pattern as one would encounter
in an associated linear problem. Hint: Which unknowns will enter
equation number i?

Exercise 5

## Show that if F (u) = 0 is a linear system of equations,

F = Au − b, for a constant matrix A and vector b, then Newton’s
method (with ω = 1) finds the correct solution in the first iteration.

Exercise 6

## The operator ∇ · (α∇u), with α = ||∇u||q , q ∈ IR, and || · || being

the Eucledian norm, appears in several physical problems,
especially flow of non-Newtonian fluids. The quantity ∂α/∂uj is
central when formulating a Newton method, where uj is the
P
coefficient in the finite element approximation u = j uj ϕj . Show
that

||∇u||q = q||∇u||q−2 ∇u · ∇ϕj .
∂uj

Exercise 7

## Consider the PDE

∂u
= ∇ · (α(u)∇u)
∂t
discretized by a Forward Euler difference in time. Explain why this
nonlinear PDE gives rise to a linear problem (and hence no need
for Newton or Picard iteration) at each time level.
Discretize the PDE by a Backward Euler difference in time and
realize that there is a need for solving nonlinear algebraic
equations. Formulate a Picard iteration method for the spatial
PDE to be solved at each time level. Formulate a Galerkin method
for discretizing the spatial problem at each time level. Choose
some appropriate boundary conditions.
Explain how to incorporate an initial condition.

Exercise 8

## Repeat Exercise 7 for the PDE

∂u
̺(u) = ∇ · (α(u)∇u)
∂t

Exercise 9

## For the problem in Exercise 8, assume that a nonlinear Newton

cooling law applies at the whole boundary:

∂u
−α(u) = H(u)(u − uS ),
∂n

## where H(u) is a nonlinear heat transfer coefficient and uS is the

temperature of the surroundings (and u is the temperature). Use a
Backward Euler scheme in time and a Galerkin method in space.
Identify the nonlinear algebraic equations to be solved at each
time level. Derive the corresponding Jacobian.
The PDE problem in this exercise is highly relevant when the
temperature variations are large. Then the density times the heat
capacity (̺), the heat conduction coefficient (α) and the heat
transfer coefficient (H) normally vary with the temperature (u).

Exercise 10

## In Exercise 8, restrict the problem to one space dimension,

choose simple boundary conditions like u = 0, use the group finite
element method for all nonlinear coefficients, apply P1 elements,
use the Trapezoidal rule for all integrals, and derive the system of
nonlinear algebraic equations that must be solved at each time
level.
Set up some finite difference method and compare the form of the
nonlinear algebraic equations.

Exercise 11

## In Exercise 8, use the Picard iteration method with one iteration at

each time level, and introduce this method at the PDE level.
Realize the similarities with the resulting discretization and that of
the corresponding linear diffusion problem.

## Nonlinear PDEs – p.33/148

Shallow water waves

Tsunamis

## Waves in fjords, lakes, or oceans, generated by

slide
earthquake
subsea volcano
asteroid
human activity, like nuclear detonation, or slides generated by oil
drilling, may generate tsunamis
Propagation over large distances
Hardly recognizable in the open ocean, but wave amplitude
increases near shore
Run-up at the coasts may result in severe damage
Giant events: Dec 26 2004 (≈ 300000 killed), 1883 (similar to
2004), 65 My ago (extinction of the dinosaurs)

## Shallow water waves – p.35/148

Norwegian tsunamis

Tromsø

70˚

Bodø

EN
ED
SW
65˚

Trondheim

AY
RW Stockholm
Oslo
NO
Bergen

60˚

10˚

15˚

## Circules: Major incidents, > 10 killed; Triangles: Selected smaller

incidents; Square: Storegga (5000 B.C.)

## Shallow water waves – p.36/148

Tsunamis in the Pacific

## Scenario: earthquake outside Chile, generates tsunami, propagating

at 800 km/h accross the Pacific, run-up on densly populated coasts in
Japan; Shallow water waves – p.37/148
Selected events; slides

Loen 1905 40m 61
Tafjord 1934 62m 41
Loen 1936 74m 73
Storegga 5000 B.C. 10m(?) ??
Vaiont, Italy 1963 270m 2600
Litua Bay, Alaska 1958 520m 2
Shimabara, Japan 1792 10m(?) 15000

## Shallow water waves – p.38/148

Selected events; earthquakes etc.

## location year strength run-up dead

Thera 1640 B.C. volcano ? ?
Thera 1650 volcano ? ?
Lisboa 1755 M=9 ? 15(?)m ?000
Portugal 1969 M=7.9 1m
Amorgos 1956 M=7.4 5(?)m 1
Krakatao 1883 volcano 40 m 36 000
Flores 1992 M=7.5 25 m 1 000
Nicaragua 1992 M=7.2 10 m 168
Sumatra 2004 M=9 50 m 300 000

## The selection is biased wrt. European events; 150 catastrophic

tsunami events have been recorded along along the Japanese coast in
modern times.

## Tsunamis: no. 5 killer among natural hazards

Shallow water waves – p.39/148
Why simulation?

## Increase the understanding of tsunamis

Assist warning systems
Assist building of harbor protection (break waters)
Recognize critical coastal areas (e.g. move population)
Hindcast historical tsunamis (assist geologists/biologists)

Problem sketch

η(x,y,t)
z
y
x

H(x,y,t)

## Assume wavelength ≫ depth (long waves)

Assume small amplitudes relative to depth
Appropriate approx. for many ocean wave phenomena

## Shallow water waves – p.41/148

Mathematical model

PDEs:
 
∂η ∂ ∂ ∂H
= − (uH) − (vH) −
∂t ∂x ∂y ∂t
∂u ∂η
= − , x ∈ Ω, t > 0
∂t ∂x
∂v ∂η
= − , x ∈ Ω, t > 0
∂t ∂y

## η(x, y, t) : surface elevation

u(x, y, t) and v(x, y, t) : horizontal (depth averaged) velocities
H(x, y) : stillwater depth (given)
Boundary conditions: either η, u or v given at each point
Initial conditions: all of η, u and v given

Primary unknowns

## Discretization: finite differences

Staggered mesh in time and space
⇒ η, u, and v unknown at different points:

ℓ ℓ+ 21 ℓ+ 12
ηi+ 1 1, ui,j+ 1 , vi+ 1 ,j+1
2 ,j+ 2 2 2

ℓ+ 21
vi+ 1 ,j+1
2

ℓ+ 12 ℓ+ 21
ui,j+ 1 ℓ
r ui+1,j+ 1
2 ηi+ 1
,j+ 1 2
2 2

ℓ+ 21
vi+ 1 ,j
2
Shallow water waves – p.43/148
A global staggered mesh

## Widely used mesh in computational fluid dynamics (CFD)

Important for Navier-Stokes solvers
Basic idea:
centered differences in time and space

## Shallow water waves – p.44/148

Discrete equations; η

∂η ∂ ∂
= − (uH) − (vH)
∂t ∂x ∂y
1 1 1
at (i + , j + , ℓ − )
2 2 2
ℓ− 12
[Dt η = −Dx (uH) − Dy (vH)]i+ 1 ,j+ 1
2 2

1 h ℓ ℓ−1
i 1 h ℓ− 21 ℓ− 21
i
ηi+ 1 ,j+ 1 − ηi+ 1 1 = − (Hu)i+1,j+ 1 − (Hu)i,j+ 1
∆t 2 2 2 ,j+ 2 ∆x 2 2

1 h ℓ− 12 ℓ− 21
i
− (Hv)i+ 1 ,j+1 − (Hv)i+ 1 ,j
∆y 2 2

## Shallow water waves – p.45/148

Discrete equations; u

∂u ∂η 1
= − at (i, j + , ℓ)
∂t ∂x 2

[Dt u = −Dx η]i,j+ 1
2

1 h ℓ+ 12 ℓ− 12
i 1 h ℓ ℓ
i
u 1 − ui,j+ 1 = − ηi+ 1 ,j+ 1 − ηi− 1 1
∆t i,j+ 2 2 ∆x 2 2 2 ,j+ 2

## Shallow water waves – p.46/148

Discrete equations; v

∂v ∂η 1
= − at (i + , j, ℓ)
∂t ∂y 2

[Dt v = −Dy η]i+ 1
2,j

1 h ℓ+ 12 ℓ− 1
i 1 h ℓ ℓ
i
vi+ 1 ,j − vi+ 12,j = η 1 1 − ηi+ 1 ,j− 1
∆t 2 2 ∆y i+ 2 ,j+ 2 2 2

## Shallow water waves – p.47/148

Complicated costline boundary

## Saw-tooth approximation to real boundary

Successful method, widely used
Warning: can lead to nonphysical waves

## Shallow water waves – p.48/148

Relation to the wave equation

## Eliminate u and v (easy in the PDEs) and get

∂2η
= ∇ · [H(x, y)∇η]
∂t2

## Eliminate discrete u and v

⇒ Standard 5-point explicit finite difference scheme for discrete η
(quite some algebra needed, try 1D first)

## Shallow water waves – p.49/148

Stability and accuracy

## Centered differences in time and space

Truncation error, dispersion analysis: O(∆x2 , ∆y 2 , ∆t2 )
Stability as for the std. wave equation in 2D:
s
− 12 1
∆t ≤ H 1 1
∆x2 + ∆y 2

(CFL condition)
If H const, exact numerical solution is possible for
one-dimensional wave propagation

## Shallow water waves – p.50/148

Verification of an implementation

## How can we verify that the program works?

Compare with an analytical solution
(if possible)
Check that basic physical mechanisms are reproduced in a
qualitatively correct way by the program

## Shallow water waves – p.51/148

Tsunami due to a slide

## Surface elevation ahead of the slide, dump behind

Initially, negative dump propagates backwards
The surface waves propagate faster than the slide moves
Shallow water waves – p.52/148
Tsunami due to faulting

## The sea surface deformation reflect the bottom deformation

Velocity of surface waves (H ∼ 5 km): 790 km/h
Velocity of seismic waves in the bottom: 6000–25000 km/h
Shallow water waves – p.53/148
Tsunami approaching the shore

p
The velocity of a tsunami is gH(x, y, t).

The back part of the wave moves at higher speed ⇒ the wave
becomes more peak-formed
Deep water (H ∼ 3 km): wave length 40 km, height 1 m
Shallow water waves – p.54/148
Tsunamis experienced from shore

## As a fast tide, with strong currents in fjords

A wall of water approaching the beach

Wave breaking: the top has larger effective depth and moves faster
than the front part (requires a nonlinear PDE)

## Shallow water waves – p.55/148

Convection-dominated flow

## Convection-dominated flow – p.56/148

Typical transport PDE

## Transport of a scalar u (heat, pollution, ...) in fluid flow v:

∂u
+ v · ∇u = α∇2 u + f
∂t

## Convection (change of u due to the flow): v · ∇u

Diffusion (change of u due to molecular collisions): α∇2 u
Common case: convection ≫ diffusion → numerical difficulties
Important dimensionless number: Peclet number Pe

|v · ∇u| V U/L VL
Pe = 2
∼ 2
=
|α∇ u| αU/L α

## V : characteristic velocity v, L: characteristic length scale, α:

diffusion constant, U : characteristic size of u

## Convection-dominated flow – p.57/148

The transport PDE for fluid flow

## The fluid flow itself is governed by Navier-Stokes equations:

∂v 1
+ v · v = − ∇p + ν∇2 v + f
∂t ̺

∇·v =0
Important dimensionless number: Reynolds number Re

convection |v · ∇v| V 2 /L VL
Re = = ∼ =
diffusion |ν∇2 v| αV /L2 ν

## Convection-dominated flow – p.58/148

A 1D stationary transport problem

## Assumption: no time, 1D, no source term

2 ′ ′′ ′ ′′ α
v · ∇u = α∇ u → vu = αu → u = ǫu , ǫ=
v

## ǫ small: boundary layer at x = 1

Standard numerics (i.e. centered differences) will fail!
Cure: upwind differences

## Convection-dominated flow – p.59/148

Notation for difference equations (1)

Define
uni+ 1 ,j,k − uni− 1 ,j,k
[Dx u]ni,j,k ≡ 2 2

h
with similar definitions of Dy , Dz , and Dt
Another difference:

uni+1,j,k − uni−1,j,k
[D2x u]ni,j,k ≡
2h

Compound difference:

1 
[Dx Dx u]ni n n n
= 2 ui−1 − 2ui + ui+1
h

## Convection-dominated flow – p.60/148

Notation for difference equations (1)

## One-sided forward difference:

n n
u − u
[Dx+ u]ni ≡ i+1 i
h

uni − uni−1
[Dx− u]ni ≡
h

[Dx Dx u = −f ]i

## Convection-dominated flow – p.61/148

Centered differences

## ui+1 − ui−1 ui−1 − 2ui + ui+1

=ǫ , i = 2, . . . , n − 1
2h h2
u1 = 0, un = 1
or
[D2x u = ǫDx Dx u]i
Analytical solution:
1 − ex/ǫ
u(x) =
1 − e1/ǫ

## Convection-dominated flow – p.62/148

Numerical experiments (1)

n=20, epsilon=0.1
1.2
1 centered
exact
0.8
0.6
0.4
u(x)

0.2
0
-0.2
-0.4
-0.6
-0.8
0 0.2 0.4 0.6 0.8 1
x

## Convection-dominated flow – p.63/148

Numerical experiments (2)

n=20, epsilon=0.01
1.2
1 centered
exact
0.8
0.6
0.4
u(x)

0.2
0
-0.2
-0.4
-0.6
-0.8
0 0.2 0.4 0.6 0.8 1
x

## Convection-dominated flow – p.64/148

Numerical experiments (3)

n=80, epsilon=0.01
1.2
1 centered
exact
0.8
0.6
0.4
u(x)

0.2
0
-0.2
-0.4
-0.6
-0.8
0 0.2 0.4 0.6 0.8 1
x

## Convection-dominated flow – p.65/148

Numerical experiments (4)

n=20, epsilon=0.001
1.2
1 centered
exact
0.8
0.6
0.4
u(x)

0.2
0
-0.2
-0.4
-0.6
-0.8
0 0.2 0.4 0.6 0.8 1
x

## Convection-dominated flow – p.66/148

Numerical experiments; summary

## The solution is not monotone if h > 2ǫ

The convergence rate is h2
(expected since both differences are of 2nd order)
provided h ≤ 2ǫ
Completely wrong qualitative behavior for h ≫ 2ǫ

Analysis

## Can find an analytical solution of the discrete problem (!)

Method: insert ui ∼ β i and solve for β

1 + h/(2ǫ)
β1 = 1, β2 =
1 − h/(2ǫ)

## cf. HPL app. A.4.4

Complete solution:
ui = C1 β1i + C2 β2i
Determine C1 and C2 from boundary conditions

β2i − β2
ui = n
β2 − β2

Important result

1 + h/(2ǫ)
⇒ <0 ⇒ h > 2ǫ
1 − h/(2ǫ)

## Must require h ≤ 2ǫ for ui to have the same qualitative property as

u(x)
This explains why we observed oscillations in the numerical
solution

## Convection-dominated flow – p.69/148

Upwind differences

Problem:

## Use a backward difference, called upwind difference, for the u′

term:
ui − ui−1 ui−1 − 2ui + ui+1
=ǫ 2
, i = 2, . . . , n − 1
h h

u1 = 0, un = 1
The scheme can be written as

## Convection-dominated flow – p.70/148

Numerical experiments (1)

n=20, epsilon=0.1
1.2
1 upwind
exact
0.8
0.6
0.4
u(x)

0.2
0
-0.2
-0.4
-0.6
-0.8
0 0.2 0.4 0.6 0.8 1
x

## Convection-dominated flow – p.71/148

Numerical experiments (2)

n=20, epsilon=0.01
1.2
1 upwind
exact
0.8
0.6
0.4
u(x)

0.2
0
-0.2
-0.4
-0.6
-0.8
0 0.2 0.4 0.6 0.8 1
x

## Convection-dominated flow – p.72/148

Numerical experiments; summary

## The solution is always monotone, i.e., always qualitatively correct

The boundary layer is too thick
The convergence rate is h
(in agreement with truncation error analysis)

Analysis

## Analytical solution of the discrete equations:

ui = β i ⇒ β1 = 1, β2 = 1 + h/ǫ

ui = C1 + C2 β2i
Using boundary conditions:

β2i − β2
ui = n
β2 − β2

## Convection-dominated flow – p.74/148

Centered vs. upwind scheme

Truncation error:
centered is more accurate than upwind
Exact analysis:
centered is more accurate than upwind when centered is stable
(i.e. monotone ui ), but otherwise useless
ǫ = 10−6 ⇒ 500 000 grid points to make h ≤ 2ǫ
Upwind gives best reliability, at a cost of a too thick boundary layer

## Convection-dominated flow – p.75/148

An interpretation of the upwind scheme

## ui − ui−1 ui−1 − 2ui + ui+1

h h2
or
[Dx− u = ǫDx Dx u]i
can be rewritten as

## ui+1 − ui−1 h ui−1 − 2ui + ui+1

= (ǫ + )
2h 2 h2
or
h
[D2x u = (ǫ + )Dx Dx u]i
2
Upwind = centered + artificial diffusion (h/2)

## Convection-dominated flow – p.76/148

Finite elements for the model problem

Galerkin formulation of

## ui+1 − ui−1 ui−1 − 2ui + ui+1

=ǫ 2
, i = 2, . . . , n − 1
2h h

u1 = 0, un = 1
or
[D2x u = ǫDx Dx u]i
Stability problems when h > 2ǫ

## Convection-dominated flow – p.77/148

Finite elements and upwind differences

## How to construct upwind differences in a finite element context?

One possibility: add artificial diffusion (h/2)

′ h ′′
u (x) = (ǫ + )u (x), x ∈ (0, 1), u(0) = 0, u(1) = 1
2

## Can be solved by a Galerkin method

Another, equivalent strategy: use of perturbed weighting functions

## Convection-dominated flow – p.78/148

Perturbed weighting functions in 1D

Take
wi (x) = ϕi (x) + τ ϕ′i (x)
or alternatively written

## where v is the standard test function in a Galerkin method

Use this wi or w as test function for the convective term u′ :
Z 1 Z 1 Z 1
u′ wdx = u′ vdx + τ u′ v ′ dx
0 0 0

## The new term τ u′ v ′ is the weak formulation of an artificial diffusion

term τ u′′ v
With τ = h/2 we then get the upwind scheme

## Convection-dominated flow – p.79/148

Optimal artificial diffusion

## [θDx− u + (1 − θ)D2x u = ǫDx Dx u]i

Is there an optimal θ?
Yes, for
h 2ǫ
θ(h/ǫ) = coth −
2ǫ h
we get exact ui (i.e. u exact at nodal points)
Equivalent artificial diffusion τo = 0.5hθ(h/ǫ)
Exact finite element method: w(x) = v(x) + τo v ′ (x) for the
convective term u′

## Convection-dominated flow – p.80/148

Multi-dimensional problems

Model problem:
v · ∇u = α∇2 u
or written out:
∂u ∂u
vx + vy = α∇2 u
∂x ∂y
Non-physical oscillations occur with centered differences or
Galerkin methods when the left-hand side terms are large
Remedy: upwind differences
Downside: too much diffusion
Important result: extra stabilizing diffusion is needed only in the
streamline direction, i.e., in the direction of v = (vx , vy )

## Convection-dominated flow – p.81/148

Streamline diffusion

## Idea: add diffusion in the streamline direction

Isotropic physical diffusion, expressed through a diffusion tensor :

d X
X d
∂2u
αδij = α∇2 u
i=1 j=1
∂xi ∂xj

## αδij is the diffusion tensor (here: same in all directions)

Streamline diffusion makes use of an anisotropic diffusion tensor
αij :
Xd X d  
∂ ∂u vi vj
αij , αij = τ 2
i=1 j=1
∂x i ∂xj ||v||

function

## Convection-dominated flow – p.82/148

Perturbed weighting functions (1)

## Consider the weighting function

w = v + τ ∗ v · ∇v
R
for the convective (left-hand side) term: w v · ∇u dΩ
This expands to
Z Z
vv · ∇udΩ + τ ∗ v · ∇u v · ∇vdΩ

The latterP
term can be viewed as the Galerkin formulation of (write
v · ∇u = i ∂u/∂xi etc.)

Xd X d  
∂ ∗ ∂u
τ vi vj
i=1 j=1
∂x i ∂xj

## Convection-dominated flow – p.83/148

Perturbed weighting functions (2)

## ⇒ Streamline diffusion can be obtained by perturbing the weighting

function
Common name: SUPG
(streamline-upwind/Petrov-Galerkin)

Consistent SUPG

## Why not just add artificial diffusion?

Why bother with perturbed weighting functions?
In standard FEM (method of weighted residuals),
Z
L(u)wΩ = 0

## the exact solution is a solution of the FEM equations (it fulfills

L(u))
This no longer holds if we
add an artificial diffusion term (∼ h/2)
use different weighting functions on different terms
Idea: use consistent SUPG
no artificial diffusion term
same (perturbed) weighting function
applies to all terms
Convection-dominated flow – p.85/148
A step back to 1D

## Let us try to use

w(x) = v(x) + τ v ′ (x)
on both terms in u′ = ǫu′′ :
Z 1 Z 1
′ ′ ′
(u v + (ǫ + τ )u v )dx + τ v ′′ u′ dx = 0
0 0

## Problem: last term contains v ′′

Remedy: drop it (!)
Justification: v ′′ = 0 on each linear (P1) element
Drop 2nd-order derivatives of v in 2D/3D too
Consistent SUPG is not so consistent...

Choosing τ ∗

## Choosing τ ∗ is a research topic

Many suggestions
Two classes:
τ∗ ∼ h
τ ∗ ∼ ∆t (time-dep. problems)
Little theory

## Convection-dominated flow – p.87/148

A test problem (1)

y
du/dn=0 or u=0
1
y = x tan θ + 0.25

v
u=1 du/dn=0
θ

0.25
u=0

u=0 1 x

## Convection-dominated flow – p.88/148

A test problem (2)

Methods:
1. Classical SUPG:
Brooks and Hughes: "A streamline upwind/Petrov-Galerkin finite
element formulation for advection domainated flows with particular
emphasis on the incompressible Navier-Stokes equations", Comp.
Methods Appl. Mech. Engrg., 199-259, 1982.

∗ v · ∇u
w = v + τ v · ∇v + τ̂ ∇u
||∇u||2

was proposed in
Hughes, Mallet and Mizukami: "A new finite element formulation
for computational fluid dynamics: II. Beyond SUPG", Comp.
Methods Appl. Mech. Engrg., 341-355, 1986.

## Convection-dominated flow – p.89/148

Galerkin’s method

−0.65
0

1 1
Z 0

X
Convection-dominated flow – p.90/148
SUPG

−0.65
0

1 1
Z 0

X
Convection-dominated flow – p.91/148
Time-dependent problems

Model problem:
∂u
+ v · ∇u = ǫ∇2 u
∂t
Can add artificial streamline diffusion term
Can use perturbed weighting function

w = v + τ ∗ v · ∇v

on all terms
How to choose τ ∗ ?

## Convection-dominated flow – p.92/148

Taylor-Galerkin methods (1)

## Idea: Lax-Wendroff + Galerkin

Model equation:
∂u ∂u
+U =0
∂t ∂x
Lax-Wendroff: 2nd-order Taylor series in time,
 n  2
n
∂u 1 ∂ u
un+1 = un + ∆t + ∆t2
∂t 2 ∂t2

## Replace temporal by spatial derivatives,

∂ ∂
= −U
∂t ∂x

Result:
 n  n
∂u 1 2 2 ∂2u
un+1 n
= u − U ∆t + U ∆t
∂x 2 ∂x2 Convection-dominated flow – p.93/148
Taylor-Galerkin methods (2)

## We can write the scheme on the form

 2
n
+ ∂u 1 2 ∂ u
Dt u + U = U ∆t 2
∂x 2 ∂x

## ⇒ Forward scheme with artificial diffusion

Lax-Wendroff: centered spatial differences,

1 2
[δt+ u + U D2x u = U ∆tDx Dx u]ni
2

## Alternative: Galerkin’s method in space,

1 2
[δt+ u + U D2x u = U ∆tDx Dx u]ni
2

## provided that we lump the mass matrix

This is Taylor-Galerkin’s method
Convection-dominated flow – p.94/148
Taylor-Galerkin methods (3)

In multi-dimensional problems,

∂u
+ v · ∇u = 0
∂t

we have

= −v · ∇
∂t
and (∇ · v = 0)

2 Xd X d  
∂ ∂ ∂
= ∇ · (vv · ∇) = v v
r s
∂t2 r=1 s=1
∂x r ∂xs

## This is streamline diffusion with τ ∗ = ∆t/2:

1
[Dt+ u + v · ∇u = ∆t∇ · (vv · ∇u)]n
2
Convection-dominated flow – p.95/148
Taylor-Galerkin methods (4)

## Can use the Galerkin method in space

(gives centered differences)
The result is close to that of SUPG, but τ ∗ is diferent
⇒ The Taylor-Galerkin method points to τ ∗ = ∆t/2 for SUPG in
time-dependent problems

## Convection-dominated flow – p.96/148

Solving linear systems

## Solving linear systems – p.97/148

The importance of linear system solvers

## PDE problems often (usually) result in linear systems of algebraic

equations
Ax = b
Special methods utilizing that A is sparse is much faster than
Gaussian elimination!
Most of the CPU time in a PDE solver is often spent on solving
Ax = b
⇒ Important to use fast methods

## Solving linear systems – p.98/148

Example: Poisson eq. on the unit cube (1)

−∇2 u = f on an n = q × q × q grid
FDM/FEM result in Ax = b system
FDM: 7 entries pr. row in A are nonzero
FEM: 7 (tetrahedras), 27 (trilinear elms.), or 125 (triquadratic
elms.) entries pr. row in A are nonzero
A is sparse (mostly zeroes)
Fraction of nonzeroes: Rq −3
(R is nonzero entries pr. row)
Important to work with nonzeroes only!

## Solving linear systems – p.99/148

Example: Poisson eq. on the unit cube (2)

## Compare Banded Gaussian elimination (BGE) versus Conjugate

Work in BGE: O(q 7 ) = O(n2.33 )
Work in CG: O(q 3 ) = O(n) (multigrid; optimal),
for the numbers below we use incomplete factorization
preconditioning: O(n1.17 )
n = 27000:
CG 72 times faster than BGE
BGE needs 20 times more memory than CG
n = 8 million:
CG 107 times faster than BGE
BGE needs 4871 times more memory than CG

## Solving linear systems – p.100/148

Classical iterative methods

Ax = b, A ∈ IRn,n , x, b ∈ IRn .

Split A: A = M − N
Write Ax = b as
M x = N x + b,
and introduce an iteration

M xk = N xk−1 + b, k = 1, 2, . . .

## Systems M y = z should be easy/cheap to solve

Different choices of M correspond to different classical iteration
methods:
Jacobi iteration
Gauss-Seidel iteration
Successive Over Relaxation (SOR)
Symmetric Successive Over Relaxation (SSOR)
Solving linear systems – p.101/148
Convergence

M xk = N xk−1 + b, k = 1, 2, . . .
The iteration converges if G = M −1 N has its largest eigenvalue,
̺(G), less than 1
Rate of convergence: R∞ (G) = − ln ̺(G)
To reduce the initial error by a factor ǫ,

||x − xk || ≤ ǫ||x − x0 ||

one needs
− ln ǫ/R∞ (G)
iterations

## Solving linear systems – p.102/148

Some classical iterative methods

Split: A = L + D + U
L and U are lower and upper triangular parts, D is A’s diagonal
Jacobi iteration: M = D (N = −L − U )
Gauss-Seidel iteration: M = L + D (N = −U )
SOR iteration: Gauss-Seidel + relaxation
SSOR: two (forward and backward) SOR steps
Rate of convergence R∞ (G) for −∇2 u = f in 2D with u = 0 as
BC:
Jacobi: πh2 /2
Gauss-Seidel: πh2
SOR: π2h
SSOR: > πh
SOR/SSOR is superior (h vs. h2 , h → 0 is small)

## Solving linear systems – p.103/148

Jacobi iteration

M =D
Put everything, except the diagonal, on the rhs
2D Poisson equation −∇2 u = f :

## ui,j−1 + ui−1,j + ui+1,j + ui,j+1 − 4ui,j = −h2 fi,j

Solve for diagonal element and use old values on the rhs:

## 1 k−1 k−1 k−1 k−1


uki,j = 2
ui,j−1 + ui−1,j + ui+1,j + ui,j+1 + h fi,j
4

for k = 1, 2, . . .

## Solving linear systems – p.104/148

Relaxed Jacobi iteration

## Dx∗ = (−L − U )xk−1 + b

Set
xk = ωx∗ + (1 − ω)xk−1
weighted mean of xk−1 and xk if ω ∈ (0, 1)

## Solving linear systems – p.105/148

Relation to explicit time stepping

∂u
α = ∇2 u + f
∂t

## by an explicit forward scheme until ∂u/∂t ≈ 0, provided

ω = 4∆t/(αh)2
Stability for forward scheme implies ω ≤ 1
In this example: ω = 1 best (⇔ largest ∆t)
Forward scheme for t → ∞ is a slow scheme, hence Jacobi
iteration is slow

## Solving linear systems – p.106/148

Gauss-Seidel/SOR iteration

M =L+D
For our 2D Poisson eq. scheme:

1 k k−1 k−1

uki,j = u k 2
+ ui−1,j + ui+1,j + ui,j+1 + h fi,j
4 i,j−1

i.e. solve for diagonal term and use the most recently computed
values on the right-hand side
SOR is relaxed Gauss-Seidel iteration:
compute x∗ from Gauss-Seidel it.
set xk = ωx∗ + (1 − ω)xk−1
ω ∈ (0, 2), with ω = 2 − O(h2 ) as optimal choice
Very easy to implement!

## Solving linear systems – p.107/148

Symmetric/double SOR: SSOR

## SSOR = Symmetric SOR

One (forward) SOR sweep for unknowns 1, 2, 3, . . . , n
One (backward) SOR sweep for unknowns n, n − 1, n − 2, . . . , 1
M can be shown to be
  −1  
1 1 1 1
M= D+L D D+U
2−ω ω ω ω

## Notice that each factor in M is diagonal or lower/upper triangular

(⇒ very easy to solve systems M y = z)

## Solving linear systems – p.108/148

Status: classical iterative methods

## Jacobi, Gauss-Seidel/SOR, SSOR are too slow for paractical PDE

computations
The simplest possible solution method for −∇2 u = f and other
stationary PDEs in 2D/3D is to use SOR
Classical iterative methods converge quickly in the beginning but
slow down after a few iterations
Classical iterative methods are important ingredients in multigrid
methods

## Solving linear systems – p.109/148

Ax = b, A ∈ IRn,n , x, b ∈ IRn .

## Use a Galerkin or least-squares method to solve a linear system

(!)
Idea: write
k
X
xk = xk−1 + αj q j
j=1

## αj : unknown coefficients, q j : known vectors

Compute the residual:

k
X
r k = b − Axk = rk−1 − αj Aq j
j=1

## Solving linear systems – p.110/148

Galerkin

Residual:
k
X
r k = b − Axk = r k−1 − αj Aq j
j=1

(r k , q i ) = 0

Galerkin’s method (r ∼ R, q j ∼ Nj , αj ∼ uj ):

(rk , q i ) = 0, i = 1, . . . , k

## (·, ·): Eucledian inner product

Result: linear system for αj ,

k
X
(Aq i , q j )αj = (rk−1 , q i ), i = 1, . . . , k
j=1

## Solving linear systems – p.111/148

Least squares

Residual:
k
X
r k = b − Axk = r k−1 − αj Aq j
j=1

(rk , rk ) = 0
∂αi

## Least squares: minimize (r k , rk )

Result: linear system for αj :

k
X
(Aq i , Aq j )αj = (rk−1 , Aq i ), i = 1, . . . , k
j=1

## Solving linear systems – p.112/148

The nature of the methods

In iteration k: seek xk in a k-dimensional vector space Vk
Basis for the space: q 1 , . . . , q k
Use Galerkin or least squares to compute the (optimal)
approximation xk in Vk
Extend the basis from Vk to Vk+1 (i.e. find q k+1 )

## Solving linear systems – p.113/148

Extending the basis

## Vk is normally selected as a so-called Krylov subspace:

Vk = span{r 0 , Ar 0 , . . . , Ak−1 r 0 }

## Alternatives for computing q k+1 ∈ Vk+1 :

k
X
q k+1 = rk + βj q j
j=1
k
X
q k+1 = Aqk + βj q j
j=1

## The first dominates in frequently used algorithms – only that

choice is used hereafter
How to choose βj ?

## Solving linear systems – p.114/148

Orthogonality properties

## Bad news: must solve a k × k linear system for αj in each

iteration (as k → n the work in each iteration approach the work of
solving Ax = b!)
The coefficient matrix in the αj system:

(Aq i , q j ), (Aq i , Aq j )

## Idea: make the coefficient matrices diagonal

That is,
Galerkin: (Aq i , q j ) = 0 for i 6= j
Least squares: (Aq i , Aq j ) = 0 for i 6= j
Use βj to enforce orthogonality of q i

## Solving linear systems – p.115/148

Formula for updating the basis vectors

Define
hu, vi ≡ (Au, v) = uT Av
and
[u, v] ≡ (Au, Av) = uT AT Av
Galerkin: require A-orthogonal qj vectors, which then results in

hrk , q i i
βi = −
hq i , q i i

results in
[rk , q i ]
βi = −
[q i , q i ]

Simplifications

(rk−1 , q k )
αk =
hq k , q k i

## and αi = 0 for i < k (!):

xk = xk−1 + αk q k

## That is, hand-derived formulas for αj

Least squares:
(rk−1 , Aq k )
αk =
[q k , q k ]
and αi = 0 for i < k

Symmetric A

## If A is symmetric (AT = A) and positive definite (positive eigenvalues

⇔ y T Ay > 0 for any y 6= 0), also βi = 0 for i < k
⇒ need to store q k only
(q 1 , . . . , q k−1 are not used in iteration k)

## Solving linear systems – p.118/148

Summary: least squares algorithm

## given a start vector x0 ,

compute r 0 = b − Ax0 and set q 1 = r 0 .
for k = 1, 2, . . . until termination criteria are fulfilled:
αk = (r k−1 , Aq k )/[q k , q k ]
xk = xk−1 + αk q k
r k = rk−1 − αk Aq k
if A is symmetric then
βk = [rk , q k ]/[q k , q k ]
q k+1 = r k − βk q k
else
βj = [r k , q j ]/[q j , q j ], j = 1, . . . , k
k
Pk
q k+1 = r − j=1 βj q j
The Galerkin-version requires A to be symmetric and positive
definite and results in the famous Conjugate Gradient method

## Solving linear systems – p.119/148

Truncation and restart

## Problem: need to store q 1 , . . . , q k

Much storage and computations when k becomes large
Truncation: work with a truncated sum for xk ,

k
X
xk = xk−1 + αj q j
j=k−K+1

## where a possible choice is K = 5

Small K might give convergence problems
Restart: restart the algorithm after K iterations
(alternative to truncation)

## Solving linear systems – p.120/148

Family of methods

## Generalized Conjugate Residual method

= least squares + restart
Orthomin method
= least squares + truncation
= Galerkin + symmetric and positive definite A
Conjugate Residuals method
= Least squares + symmetric and positive definite A
Many other related methods: BiCGStab, Conjugate Gradients
Squared (CGS), Generalized Minimum Residuals (GMRES),
Minimum Residuals (MinRes), SYMMLQ
All of these are easily called in Diffpack

Convergence

## Conjugate Gradient-like methods converge slowly (but usually

faster than SOR/SSOR)
To reduce the initial error by a factor ǫ,

1 2√
ln κ
2 ǫ

## iterations are needed, where κ is the condition number:

largest eigenvalue of A
κ=
smalles eigenvalue of A

## κ = O(h−2 ) when solving 2nd-order PDEs

(incl. elasticity and Poisson eq.)

Preconditioning

M −1 Ax = M −1 b

## solve it with a Conjugate Gradient-like method and construct M

such that
1. κ = O(1) ⇒ M ≈ A (i.e. fast convergence)
2. M is cheap to compute
3. M is sparse (little storage)
4. systems M y = z (occuring in the algorithm due to
M −1 Av-like products) are efficiently solved (O(n) op.)
The preconditioning business: find a good balance between 1-4

## Solving linear systems – p.123/148

Classical methods as preconditioners

## Idea: “solve” M y = z by one iteration with a classical iterative

method (Jacobi, SOR, SSOR)
Jacobi preconditioning: M = D (diagonal of A)
No extra storage as M is stored in A
No extra computations as M is a part of A
Efficient solution of M y = z
But: M is probably not a good approx to A
⇒ poor quality of this type of preconditioners?
Conjugate Gradient method + SSOR preconditioner is widely
used

## Solving linear systems – p.124/148

M as a factorization of A

M = LU

## where L and U are lower and upper triangular matrices resp.

Implications:
1. M = A (κ = 1): very efficient preconditioner!
2. M is not cheap to compute
(requires Gaussian elim. on A!)
3. M is not sparse (L and U are dense!)
4. systems M y = z are not efficiently solved
(O(n2 ) process when L and U are dense)

## Solving linear systems – p.125/148

M as an incomplete factorization of A

## New idea: compute sparse L and U

How? compute only with nonzeroes in A
bU
⇒ Incomplete factorization, M = L b 6= LU
M is not a perfect approx to A
M is cheap to compute and store
(O(n) complexity)
M y = z is efficiently solved (O(n) complexity)
This method works well - much better than SOR/SSOR
preconditioning

How to compute M

## Run through a standard Gaussian elimination, which factors A as

A = LU
Normally, L and U have nonzeroes where A has zeroes
Idea: let L and U be as sparse as A
Compute only with the nonzeroes of A
Such a preconditioner is called Incomplete LU Factorization, ILU
Option: add contributions outside A’s sparsity pattern to the
diagonal, multiplied by ω
Relaxed Incomplete Factorization (RILU): ω > 1
Modified Incomplete Factorization (MILU): ω = 1
See algorithm C.3 in the book

## Solving linear systems – p.127/148

Numerical experiments

## Two test cases:

−∇2 u = f on the unit cube and FDM
−∇2 u = f on the unit cube and FEM
Diffpack makes it easy to run through a series of numerical
experiments, using multiple loops, e.g.,
sub LinEqSolver_prm
set basic method = { ConjGrad & MinRes }
ok
sub Precond_prm
set preconditioning type = PrecRILU
set RILU relaxation parameter = { 0.0 & 0.4 & 0.7 & 1.0 }
ok

## Solving linear systems – p.128/148

Test case 1: 3D FDM Poisson eq.

Equation: −∇2 u = 1
Boundary condition: u = 0
7-pt star standard finite difference scheme
Grid size: 20 × 20 × 20 = 8000 points and 20 × 30 × 30 = 27000
points
Source code:
\$NOR/doc/Book/src/linalg/LinSys4/
All details in HPL Appendix D
Input files:
\$NOR/doc/Book/src/linalg/LinSys4/experiments
Solver’s CPU time written to standard output

## Solving linear systems – p.129/148

Jacobi vs. SOR vs. SSOR

## n = 203 = 8000 and n = 303 = 27000

Jacobi: not converged in 1000 iterations
SOR(ω = 1.8): 2.0s and 9.2s
SSOR(ω = 1.8): 1.8s and 9.8s
Gauss-Seidel: 13.2s and 97s
SOR’s sensitivity to relax. parameter ω:
1.0: 96s, 1.6: 23s, 1.7: 16s, 1.8: 9s, 1.9: 11s
SSOR’s sensitivity to relax. parameter ω:
1.0: 66s, 1.6: 17s, 1.7: 13s, 1.8: 9s, 1.9: 11s
⇒ relaxation is important,
great sensitivity to ω

## Compare Conjugate Residuals with Conjugate Gradients

Or: least squares vs. Galerkin
MinRes: not converged in 1000 iterations
⇒ ConjGrad is clearly faster than the best SOR/SSOR
MinRes: 0.7s and 4s
The importance of preconditioning grows as n grows

## Solving linear systems – p.131/148

Different preconditioners

## ILU, Jacobi, SSOR preconditioners (ω = 1.2)

MinRes:
Jacobi: not conv., SSOR: 11.4s, ILU: 4s
Jacobi: 4.8s, SSOR: 2.8s, ILU: 2.7s
Sensitivity to relax. parameter in SSOR, with ConjGrad as solver:
1.0: 3.3s, 1.6: 2.1s, 1.8: 2.1s, 1.9: 2.6s
Sensitivity to relax. parameter in RILU, with ConjGrad as solver:
0.0: 2.7s, 0.6: 2.4s, 0.8: 2.2s, 0.9: 1.9s,
0.95: 1.9s, 1.0: 2.7s
⇒ ω slightly less than 1 is optimal,
RILU and SSOR are equally fast (here)

## Solving linear systems – p.132/148

Test case 2: 3D FEM Poisson eq.

## Equation: −∇2 u = A1 π 2 sin πx + 4A2 π 2 sin 2πy + 9A3 π 2 sin 3πz

Boundary condition: u known
ElmB8n3D and ElmB27n3D elements
Grid size: 21 × 21 × 21 = 9261 nodes and 31 × 31 × 31 = 29791
nodes
Source code:
\$NOR/doc/Book/src/fem/Poisson2
All details in HPL Chapter 3.2 and 3.5
Input files:
\$NOR/doc/Book/src/fem/Poisson2/linsol-experiments
Solver’s CPU time available in casename-summary.txt

## Solving linear systems – p.133/148

Jacobi vs. SOR vs. SSOR

## n = 9261 and n = 303 = 29791, trilinear and triquadratic elms.

Jacobi: not converged in 1000 iterations
SOR(ω = 1.8): 9.1s and 81s, 42s and 338s
SSOR(ω = 1.8): 47s and 248s, 138s and 755s
Gauss-Seidel: not converged in 1000 iterations
SOR’s sensitivity to relax. parameter ω:
1.0: not conv., 1.6: 200s, 1.8: 83s, 1.9: 57s
(n = 29791 and trilinear elements)
SSOR’s sensitivity to relax. parameter ω:
1.0: not conv., 1.6: 212s, 1.7: 207s, 1.8: 245s, 1.9: 435s
(n = 29791 and trilinear elements)
⇒ relaxation is important,
great sensitivity to ω

## Compare Conjugate Residuals with Conjugate Gradients

Or: least squares vs. Galerkin
MinRes: not converged in 1000 iterations
9261 vs 29791 unknowns, trilinear elements
⇒ ConjGrad is clearly faster than the best SOR/SSOR!
MinRes: 5s and 28s
ILU prec. has a greater impact when using triquadratic elements
(and when n grows)

## Solving linear systems – p.135/148

Different preconditioners

## ILU, Jacobi, SSOR preconditioners (ω = 1.2)

MinRes:
Jacobi: 68s., SSOR: 57s, ILU: 28s
Jacobi: 19s, SSOR: 14s, ILU: 16s
Sensitivity to relax. parameter in SSOR, with ConjGrad as solver:
1.0: 17s, 1.6: 12s, 1.8: 13s, 1.9: 18s
Sensitivity to relax. parameter in RILU, with ConjGrad as solver:
0.0: 16s, 0.6: 15s, 0.8: 13s, 0.9: 12s, 0.95: 11s, 1.0: 16s
⇒ ω slightly less than 1 is optimal,
RILU and SSOR are equally fast (here)

## Solving linear systems – p.136/148

More experiments

Convection-diffusion equations:
\$NOR/doc/Book/src/app/Cd/Verify
Files: linsol_a.i etc as for LinSys4 and Poisson2
Elasticity equations:
\$NOR/doc/Book/src/app/Elasticity1/Verify
Files: linsol_a.i etc as for the others
Run experiments and learn!

## Solving linear systems – p.137/148

Multigrid methods

Multigrid methods are the most efficient methods for solving linear
systems
Multigrid methods have optimal complexity O(n)
Multigrid can be used as stand-alone solver or preconditioner
Multigrid applies a hierarchy of grids
Multigrid is not as robust as Conjugate Gradient-like methods and
incomplete factorization as preconditioner, but faster when it
works
Multigrid is complicated to implement
Diffpack has a multigrid toolbox that simplifies the use of multigrid
dramatically

## Solving linear systems – p.138/148

The rough ideas of multigrid

## Observation: e.g. Gauss-Seidel methods are very efficient during

the first iterations
High-frequency errors are efficiently damped by Gauss-Seidel
Low-frequence errors are slowly reduced by Gauss-Seidel
higher frequency
Repeat the procedure
On the coarsest grid: solve the system exactly
Transfer the solution to the finest grid
Iterate over this procedure

## Solving linear systems – p.139/148

Damping in Gauss-Seidel’s method (1)

j+1 + h fj

i :

j+1

is a pseudo time

## Solving linear systems – p.140/148

Damping in Gauss-Seidel’s method (2)

## Can find eℓj with techniques from Appendix A.4:

X
eℓj = Ak exp (i(kjh − ω̃ℓ∆t))
k

## or (easier to work with here):

X
eℓj = Ak ξ ℓ exp (ikjh), ξ = exp (−iω̃∆t)
k

## Inserting a wave component in the scheme:

exp (ikh) 1
ξ = exp (−iω̃∆t) = , |ξ| = √
2 − exp (−ikh) 5 − 4 cos kh

## Solving linear systems – p.141/148

Gauss-Seidel’s damping factor

1
|ξ| = √ , p = kh ∈ [0, π]
5 − 4 cos p

0.9

0.8

0.7

0.6

0.5

0.4

p

## Small p = kh ∼ h/λ: low frequency (relative to the grid) and small

damping
Large (→ π) p = kh ∼ h/λ: high frequency (relative to the grid)
and efficient damping
Solving linear systems – p.142/148
More than one grid

## From the previous analysis: error components with high frequency

are quickly damped
p is increased by a factor of 2, i.e., not so high-frequency waves
on the h grid is efficiently damped by Gauss-Seidel on the h′ grid
Repeat the procedure
On the coarsest grid: solve by Gaussian elimination
Interpolate solution to a finer grid, perform Gauss-Seidel
iterations, and repeat until the finest grid is reached

## Solving linear systems – p.143/148

Transferring the solution between grids

## From fine to coarser: restriction

simple restriction

weighted restriction

## fine grid function

1 2 3 4 5 q-1
1 2 3 4 5 6 7 8 9 q

## coarse grid function

1 2 3 4 5 q−1
1 2 3 4 5 6 7 8 9 q

Smoothers

## The Gauss-Seidel method is called a smoother when used to

damp high-frequency error components in multigrid
Other smoothers: Jacobi, SOR, SSOR, incomplete factorization
No of iterations is called no of smoothing sweeps
Common choice: one sweep

## Solving linear systems – p.145/148

A multigrid algorithm

Perform smoothing (pre-smoothing)
Restrict to coarser grid
Repeat the procedure (recursive algorithm!)
On the coarsest grid: solve accurately
Prolongate to finer grid
Perform smoothing (post-smoothing)
One cycle is finished when reaching the finest grid again
Can repeat the cycle
Multigrid solves the system in O(n) operations

V- and W-cycles

q γq =1 γq =2
4

smoothing

## Solving linear systems – p.147/148

Multigrid requires flexible software

## Many ingredients in multigrid:

pre- and post-smoother
no of smoothing sweeps
solver on the coarsest level
cycle strategy
restriction and prolongation methods
how to construct the various grids?
There are also other variants of multigrid (e.g. for nonlinear
problems)
The optimal combination of ingredients is only known for simple
model problems (e.g. the Poisson eq.)
In general: numerical experimentation is required!
(Diffpack has a special multigrid toolbox for this)
{"found_math": false, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0}