Direct Numerical Simulations of Multiphase Flows

1
Numerical Methods for Elliptic Equations
Instructor Hong G. Im University of Michigan
Fall 2005
2
Outline
Solution Methods for Elliptic Equations

• Direct Methods
• Cramers rule
• Gaussian elimination
• Tridiagonal matrix algorithm (Thomas algorithm)
• Block tridiagonal matrix algorithm
• Iterative Methods
• Jacobi, Gauss-Seidel iteration
• Successive overrelaxation (SOR)
• Multigrid method

3
A model equation for elliptic problems
Poisson equation (Laplace equation if S 0).
On the boundaries (BC)
Dirichlet
Neumann
4
Examples of Elliptic Equations
Steady conduction equation
2-D stream function equation
Projection method (Step 2)
5
Solving the Poisson equation
j1 j j?1
i?1 i i1
Applying central differencing
for which the modified equation becomes
6
If
Sparse matrix only 5 non-zero entries in each row
7
Ultimately, the difference form of the Poisson
equation boils down to solving for
Hence,
Direct method - Solving inverse matrix
directly (Cramers rule) - Inverting matrix
more cleverly (Gaussian elimination) - Other
(L-U decomposition, Thomas algorithm)
8
Direct Method - 1
Cramers Rule

• Check elementary linear algebra book
• Extremely time consuming for large matrices
• Operation count
• For example, 11?11 matrix
• Using 1 Gigaflops (109 flops) computer

9
Direct Method - 2
Gaussian Elimination

• Pivoting rearranging equations to put the
  largest
• coefficient on the main diagonal.
• Eliminate the column below main diagonal.
• Repeat until the last equation is reached.
• Back-substitution
• Special case tri-diagonal matrix - Thomas
  algorithm

?
10
Iterative Method
Discretized Poisson Equation
Rearranging for
Jacobi
Gauss-Seidel
SOR
11
Convergence 1-D Problem (1)
A One Dimensional Example
An equation in the form
can be solved by iterative procedure
for which convergence is achieved if
or
12
Convergence 1-D Problem (2)
The above figure suggests that
for convergence
13
Convergence Multi Dimensional (1)
For multi-dimensional problems
It can be shown that the space of can be
spanned by eigenvectors
such that
Hence,
14
Convergence Multi Dimensional (2)
Therefore, the problem
is equivalent to
or
yielding a convergence condition
15
Convergence Multi Dimensional (3)
General iterative procedure
Let
An iterative scheme is constructed as
For example,
Jacobi
Gauss-Seidel
16
Convergence Multi Dimensional (4)
Requirements

• should be invertible.
• Iteration should converge, i.e.

Define error at n-th iteration
17
Convergence Multi Dimensional (5)
Condition for convergence
which requires
This is achieved if the modulus of the
eigenvalues are less than unity. Therefore, the
convergence condition becomes
Spectral radius of convergence
Eigenvalues of matrix
18
Convergence Jacobi Method (1)
Jacobi Iteration Method for Poisson
Equation (2nd order central difference with
uniform mesh)
and using a discrete analog of separation of
variables, it can be shown that the eigenvalues
of are
Therefore, and the Jacobi method
converges.
19
Convergence Jacobi Method (2)
For a large matrix
largest eigenvalues for
Thus, for a large matrix, is only
slightly less than unity. ? Very slow convergence
20
Convergence Gauss-Seidel (1)
Gauss-Seidel Iteration Method for Poisson
Equation (2nd order central difference with
uniform mesh)
A1
A2
21
Convergence Gauss-Seidel (2)
It can be shown that the eigenvalues of
matrix are simply square of the eigenvalues of
Jacobi method
Thus, Gauss-Seidel method is twice as fast as the
Jacobi method.
22
Successive Overrelaxation - 1
Successive Overrelaxation
Consider the Gauss-Seidel method
If Gauss-Seidel is an attempt to change the
solution as
Can we accelerate the change by a parameter?
23
Successive Overrelaxation - 2
Hence, SOR first uses Gauss-Seidel to compute
intermediate solution,
or
Then accelerate the next iteration solution
Combining the two steps
Note that
24
Successive Overrelaxation - 3
Example Poisson Equation (2nd order CS, uniform
mesh)
G-S
Combining
25
Successive Overrelaxation - 4
Convergence of SOR
From
Eliminating and solving for
Convergence depends on the eigenvalues of
26
Successive Overrelaxation - 5
For the discretized Poisson operator, it can be
shown that
where is an eigenvalue of the Jacobi
matrix
Note that if
(Gauss-Seidel)
Minimum occurs at
Typically
27
Successive Overrelaxation - 6
For problems with irregular geometry and
non-uniform mesh, must be found by trial
and error.
Typical Comparison Chart
Ferziger, J. H., Numerical Method for Engineering
Application (1981)
28
Coloring Scheme (Red Black)
In large computer application (vector or parallel
platform),
SOR faces difficulties in using constantly
updated values.
Remedy Two separate grid system (red black)
do i 1, nx, 2
enddo
do i2, nx, 2
enddo
29
Successive Line Overrelaxation (SLOR) - 1
Line Relaxation Method (Line Gauss-Seidel Method)
Old
New
New
Old
Adding one more coupling (E)
New
Old
New
New
New
New
? Thomas algorithm
30
Successive Line Overrelaxation (SLOR) - 2
SLOR Line Relaxation Overrelaxation
Apply line relaxation for intermediate solution
and then overrelax
which is no more complicated than line
relaxation.
31
Successive Line Overrelaxation (SLOR) - 3
Notes on SLOR

• Exact eigenvalues are unknown.
• To ensure convergence,
• Converges approximately twice as fast as
  Gauss-Seidel.
• May be faster than pointwise SOR, but each
  iteration
• takes longer with Thomas algorithm.
• Improved convergence is due to the direct effect
  of the
• boundary condition in each row.

32
Alternating-Direction Implicit - 1
ADI for elliptic equation is analogous to ADI in
parabolic equation
In discrete form
and take it to the limit to obtain the steady
solution.
33
Alternating-Direction Implicit - 2
ADI for
is written as
or
34
Alternating-Direction Implicit - 3
Convergence of ADI

• Iteration parameter usually
  varies with iteration
• For example, Wachspress

lower bound eigenvalue
upper bound eigenvalue

• Comparison with SOR is difficult
• ADI can be efficient if appropriate parameters
  are found.

35
Boundary Conditions for Iterative Method
Dirichlet conditions are easily implemented. For
Neumann condition, the simplest approach is
(1st order)
Update interior points
and then set
? This generally does not converge.
Instead, incorporate BC directly into the
equations
36
Multigrid Method - 1
Basic Idea
Suppose we want to solve a 1-D elliptic equation
An iterative process is analogous to solving an
unsteady problem (recall ADI)
And take it to the limit
37
Multigrid Method - 2
Analytic Solution Behavior
The equation can be solved by Fourier series
and assume that can be expanded as
Substituting into the equation,
38
Multigrid Method - 3
Solving for each
Therefore,
Rate of convergence
High wave number modes damps out faster.
39
Multigrid Method - 4
For iterative method to solve elliptic
equations, there are high wave number and low
wave number components of errors that need to be
damped out.
Consider a domain discretized by
points
kL2?
kH2?N
40
Multigrid Method - 5
For explicit time step, stability condition
yields
If the error at various wave number decays at
where
?
Short-wave errors decay much faster
41
Multigrid Method - 6
t 0
t 5?t
t 500?t
42
Multigrid Method - 7
Multigrid Method the Idea
- A low wave number component on a fine grid
becomes a high wave number component on a
coarse grid - Use a coarse grid system to
converge low wave number component of the
solution rapidly - Map it onto the fine grid
system to converge high wave number component.
43
Multigrid Method - 8
Suppose we are solving a Laplace equation
and iterate until residual becomes smaller than
tolerance
Define
Converged solution Solution after n-th
iteration Correction
44
Multigrid Method - 9
Since
and
we have
Coarse grid correction
Steps Fine grid solution ? Interpolation
onto coarse grid (restriction) ? Coarse grid
correction ? Interpolated onto fine grid
(prolongation) ? Converge on fine grid
45
Multigrid Method - 10
Two-Level Multigrid Method - (nx1, ny1)
(nx/21, ny/21)

• Do n iterations on the fine grid (n 34) using
  G-S
• 2. If then stop.
• 3. Interpolate the residual onto the
  coarse grid
• (restriction (injection if nx/2))
• 4. Iterations on the coarse grid
• 5. Interpolate the correction onto
  the fine grid
• (prolongation)
• 6. Go to 1.

46
Multigrid Method - 11
Details on Prolongation
Coarse Grid Fine grid
47
Multigrid Method - 12
The process can be extended to multi-level
grids (nx1,ny1), (nx/21,ny/21),
(nx/41,ny/41), , (3,3)
Various Strategies
Fine
Coarse
exact solution
V-Cycle
W-Cycle
48
Multigrid Method - 13
Example 4-Level V-Cycle
Fine ? Coarse

• Iterate on Grid 1 for n times
• 2. Restriction by injection to Grid 2
  (Save )
• Iterate on Grid 2
• 
  (Save )
• Restriction by injection to Grid 3
• Iterate on Grid 3
• (Save )
• 6. Restriction by injection to Grid 4
• 7. Iterate on Grid 4

49
Multigrid Method - 14
Coarse ? Fine

• Prolongate from Grid 4 to Grid 3
• 2. Iterate on Grid 3
  (saved)
• Prolongate from Grid 3 to Grid 2
• Iterate on Grid 2
  (saved)
• Prolongate from Grid 2 to Grid 1
• Iterate on Grid 1
• If then stop. Else, repeat the
  entire cycle.

50
Multigrid Method - 15

• Test Problem (Tannehill, p. 170)
• Laplace equation in a square domain
• Dirichlet conditions on four boundaries
• 5 Levels of resolution
• 99, 1717, 3333, 6565, 129129
• 4 Methods
• 1. Conventional Gauss-Seidel (GS)
• 2. Gauss-Seidel-SOR with optimal ? (GSopt )
• 3. Multigrid with 2-level grids (MG2)
• 4. Multigrid with maximum levels of grid (MGMAX)

51
Multigrid Method - 16
52
Multigrid Method - 17
The multigrid method is the fastest
technique currently available for general
elliptic equations, and are said to compare
favorably with the fastest direct solvers for
special problems. MUDPACK http//www.scd.ucar.edu
/css/software/mudpack/ FISHPACK http//www.scd.uca
r.edu/css/software/fishpack/