A 2 h ( u 2 h ) What is in FAS? • As in the linear case, there are two basic possibilities: A 2 h ( u 2 h ) • is determined by discretizing the nonlinear operator A(u) in the same fashion as was employed to A h ( u h ) obtain , except that the coarser mesh spacing is used. A 2 h ( u 2 h ) • is determined from the Galerkin condition 2 h A h ( u h ) I 2 h A 2 h ( u 2 h ) = I h h where the action of the Galerkin product can be captured in an implementable formula. • The first method is usually easier, and more common. 21 of 104
Nonlinear problems: an example • Consider − ∆ u ( x , y ) + γ u ( x , y ) e u ( x , y ) = f ( x , y ) on the unit square [0,1] x [0,1] with homogeneous Dirichlet boundary conditions and a regular Cartesian grid. • Suppose the exact solution is u ( x , y ) = ( x 2 − x 3 ) sin ( 3 π y ) 22 of 104
Discretization of nonlinear example • The operator can be written (sloppily) as 1 − 1 h j u � � h j h j i , u u e f 1 4 1 + γ = − − i i i j , , , � � h 2 1 − � � h h A u ( ) • The relaxation is given by h h A u f ( ( ) ) − i j , i j , h j h j u u ← − i i , , 4 u h j i , j 1 u e + γ ( + ) i , 2 h 23 of 104
FAS and Newton’s method on − ∆ u ( x , y ) + γ u ( x , y ) e u ( x , y ) = f ( x , y ) • FAS γ 1 10 100 1000 convergence factor 0.135 0.124 0.098 0.072 number of FAS cycles 12 11 11 10 • Newton’s Method γ 1 10 100 1000 convergence factor 4.00E-05 7.00E-05 3.00E-04 2.00E-04 number of Newton iterations 3 3 3 4 24 of 104
Newton, Newton-MG, and FAS on − ∆ u ( x , y ) + γ u ( x , y ) e u ( x , y ) = f ( x , y ) • Newton uses exact solve, Newton-MG is inexact Newton with a fixed number of inner V(2,1)-cycles the Jacobian problem, overall stopping criterion 1 0 1 − | r 0 | | | < 2 Outer Inner Method iterations iterations Megaflops Newton 3 1660.6 Newton-MG 3 20 56.4 Newton-MG 4 10 38.5 Newton-MG 5 5 25.1 Newton-MG 10 2 22.3 Newton-MG 19 1 24.6 FAS 11 27.1 25 of 104
Comparing FMG-FAS and FMG-Newton − ∆ u ( x , y ) + γ u ( x , y ) e u ( x , y ) = f ( x , y ) We will do one FMG cycle using a single FAS V(2,1) - cycle as the “solver” at each new level. We then follow that with sufficiently many FAS V(2,1)-cycles as is necessary to obtain ||r|| < 10 -10 . Next, we will do one FMG cycle using a Newton- multigrid step at each new level (with a single linear V(2,1)-cycle as the Jacobian “solver.”) We then follow that with sufficiently many Newton-multigrid steps as is necessary to obtain ||r|| < 10 -10 . 26 of 104
Comparing FMG-FAS and FMG-Newton − ∆ u ( x , y ) + γ u ( x , y ) e u ( x , y ) = f ( x , y ) | r h Cycle | r h | e h Mflops | e h Mflops Cycle | | | | | | | | | | | | FMG-FAS 1.10E-02 2.00E-05 3.1 1.06E-02 2.50E-05 2.4 FMG-Newton FAS V 6.80E-04 2.40E-05 5.4 6.70E-04 2.49E-05 4.1 Newton-MG FAS V 5.00E-05 2.49E-05 7.6 5.10E-05 2.49E-05 5.8 Newton-MG FAS V 3.90E-06 2.49E-05 9.9 6.30E-06 2.49E-05 7.5 Newton-MG FAS V 3.20E-07 2.49E-05 12.2 1.70E-06 2.49E-05 9.2 Newton-MG FAS V 3.00E-08 2.49E-05 14.4 5.30E-07 2.49E-05 10.9 Newton-MG FAS V 2.90E-09 2.49E-05 16.7 1.70E-07 2.49E-05 12.6 Newton-MG FAS V 3.00E-10 2.49E-05 18.9 5.40E-08 2.49E-05 14.3 Newton-MG FAS V 3.20E-11 2.49E-05 21.2 1.70E-08 2.49E-05 16.0 Newton-MG 5.50E-09 2.49E-05 17.7 Newton-MG 1.80E-09 2.49E-05 19.4 Newton-MG 5.60E-10 2.49E-05 21.1 Newton-MG 1.80E-10 2.49E-05 22.8 Newton-MG 5.70E-11 2.49E-05 24.5 Newton-MG 27 of 104
Outline ✔ • Nonlinear Problems • Neumann Boundary Conditions • Anisotropic Problems • Variable Mesh Problems • Variable Coefficient Problems • Algebraic Multigrid 28 of 104
Neumann Boundary Conditions • Consider the (1-d) problem u x f x < x 0 1 − ″ ( ) = ( ) < u 0 u 1 0 ′ ( ) = ′ ( ) = 1 • We discretize this on the interval [0,1] with h = N 1 + grid spacing for j=1,2, …, N+1 . x j = j h • We extend the interval with two ghost points 0 x 1 N+2 -1 0 1 2 3 … j ... N-1 N N+1 29 of 104
We use central differences • We approximate the derivatives with differences, using the ghost points: 0 1 j-1 j j+1 N N+1 N+2 -1 u 2 u u u u − + − − u u − j 1 j j 1 N + 2 N − + 1 1 − u 1 u 0 u x ′ ( ) ≈ ″ ( ) ≈ ′ ( ) ≈ j 2 h 2 h 2 h • Giving the system u 2 u u − + − j 1 j j 1 − + 0 j N 1 ≤ ≤ + f = j 2 h u u u u − − 1 1 N + 2 N − 0 0 = = 2 h 2 h 30 of 104
Eliminating the ghost points • Use the boundary conditions to eliminate , u − 1 u N + 2 u u u u − − N + 2 N u u 1 1 N + 2 = − u − u 1 = 0 N = 0 = 1 2 h 2 h • Eliminating the ghost points in the j=0 and j=N+1 equations gives the (N+2)x(N+2) system of equations: u 2 u u − + − j 1 j j 1 − + 0 j N 1 ≤ ≤ + f = j 2 h 2 u 2 u − + N N 1 2 u 2 u + − 0 1 f = N 1 + f = 2 0 h 2 h 31 of 104
We write the system in matrix form h h h • We can write , where A u f = � 2 − 2 � − 1 2 − 1 � � � � − 1 2 − 1 A h = 1 � � ⋅⋅⋅ ⋅⋅⋅ � � h 2 � � − 1 2 − 1 � � � � − 1 2 � � A h • Note that is (N+2)x(N+2), nonsymmetric, and h the system involves unknowns and at the h u N + 1 u 0 boundaries. 32 of 104
We must consider a compatibility condition • The problem , for and u x f x − ″ ( ) = ( ) 0 < x 1 < with is not well-posed! u 0 u 1 0 ′ ( ) = ′ ( ) = • If u(x) is a solution, so is u(x)+c for any constant c. • We cannot be certain a solution exists. If one does, it must satisfy 1 1 1 u ″( x ) dx = � f ( x ) dx − [ u ′( 1 ) − u ′( 0 ) ] = � f ( x ) dx − � � � � 0 0 0 1 0 = � f ( x ) dx � 0 • This integral compatibility condition is necessary! If f(x) doesn’t satisfy it, there is no solution! 33 of 104
The well-posed system • The compatibility condition is necessary for a solution to exist. In general, it is also sufficient, which can be proven that is a well-behaved 2 ∂ − 2 operator in the space of functions u(x) that have x ∂ zero mean. • Thus we may conclude that if f(x) satisfies the compatibility condition, this problem is well-posed: u x f x < x − ″ ( ) = ( ) 0 1 < u 0 u 1 0 ′ ( ) = ′ ( ) = 1 u x d x 0 ( ) = � � 0 • The last says that of all solutions u(x)+c we choose the one with zero mean. 34 of 104
The discrete problem is not well posed A h • Since all row sums of are zero, 1 h ∈ NS ( A h ) A h • Putting into row-echelon form shows that A h h h hence d i m N S 1 N S ( A s p a n 1 ( ( ) ) = ) = ( ) A h • By the Fundamental Theorem of Linear Algebra, h h T has a solution if and only if f N S A ⊥ ( ( ) ) T T • It is easy to show that h N S A c 1 2 1 1 1 1 2 ( ( ) ) = ( / , , , . . . , , / ) h h • Thus, has a solution if and only if h A u f = T h f c 1 2 1 1 1 1 2 ⊥ ( / , , , . . . , , / ) • That is, N 1 1 h h h f f f 0 + � + = 0 j N 1 + 2 2 j 1 = 35 of 104
We have two issues to consider h h T • Solvability. A solution exists iff f N S A ⊥ ( ( ) ) u h + c 1 h A h u h = f h • Uniqueness. If solves so does u h A h = ( A h ) T NS ( A h ) = NS ( ( A h ) T ) • Note that if then and solvability and uniqueness can be handled together • This is easily done. Multiply 1st and last equations by 1/2, giving � 1 − 1 � − 1 2 − 1 � � � � − 1 2 − 1 h = 1 � � A ⋅⋅⋅ ⋅⋅⋅ � � h 2 � � − 1 2 − 1 � � � � − 1 1 � 36 of 104 �
The new system is symmetric h u h = f h • We have the symmetric system : A f h 0 / 2 � � h u 0 � � � � � 1 − 1 � f h � � � � h u 1 − 1 2 − 1 1 � � � � � � � � � � � � f h − 1 2 − 1 1 h u 2 � � � � � � 2 = ⋅⋅⋅ ⋅⋅⋅ � � � � � � h 2 � � � � � � � � u h − 1 2 − 1 � � � � � � f h N � � N � � � � − 1 1 � u h � � � � � f h N + 1 � � N + 1 / 2 � � � � h • Solvability is guaranteed by ensuring that is f 1 h orthogonal to the constant vector : N 1 + h h h f 1 f 0 , = = j � j 0 = 37 of 104
The well-posed discrete system • The (N+3)x(N+2) system is: u 2 u u − + − j 1 j j 1 − + 0 j N 1 ≤ ≤ + f = j 2 h u 0 − u 1 f 0 = 2 h 2 − u N + u N + 1 f N + 1 = 2 h 2 N 1 + (choose the zero mean solution) h u = 0 i � i 0 = or, more simply h u h = f h A h h u 1 0 , = 38 of 104
Multigrid for the Neumann Problem • We must have the interval endpoints on all grids x h x h h / 2 h x N x N 0 1 2 h x 0 2 h 2 h x x N 2 N 4 / / • Relaxation is performed at all points, including endpoints: h 2 h h v v h f + + + h h j 1 j 1 j − + 2 2 h h h h v v h f v v h f h ← + ← + v ← N 1 N N 1 0 1 0 + + j 2 • We add a global Gram-Schmidt step after relaxation on each level to enforce the zero-mean condition h h v 1 , h h h v v 1 ← − h h 1 1 , 39 of 104
Interpolation must include the endpoints • We use linear interpolation: 1 � � 1 2 1 2 � / / � � � 1 � � 1 2 1 2 � / / � h � � I 2 1 = h � � 1 2 1 2 � / / � � � 1 � � 1 2 1 2 � / / � � � 1 � � 40 of 104
Restriction also treats the endpoints 1 T 2 h h • For restriction, we use , yielding the I = 2 I ( ) h 2 h values 1 1 2 h h h f f f = + 0 0 1 2 4 1 1 1 2 h h h h f f f f = + + j 2 j 1 2 j 2 j 1 − + 4 2 4 1 1 2 h h h f f f = + N 1 N N 1 + + 4 2 41 of 104
The coarse-grid operator • We compute the coarse-grid operator using the Galerkin condition h = 2 h 2 h h A I A I h 2 h j 2h = 0 1 2 h 1 0 ε � 1 − 1 0 � − 1 2 − 1 � � 1 h � � 2 h I 1 0 ε − 1 2 − 1 1 2 h = 2 h 0 � � 2 A ⋅⋅⋅ ⋅⋅⋅ � � 4 h 2 1 1 � � 2 h h 2 h − 1 2 − 1 A I 0 � � ε − h 2 h 0 2 2 � � 2 h 2 h − 1 1 � � 1 1 h 2 h h 2 h I A I ε − h 2 h 0 2 2 4 h 4 h 42 of 104
Coarse-grid solvability h h h • Assuming satisfies , the solvability f f 1 0 , = condition, we can show that theoretically the coarse- 2 h h h 2 h grid problem is also solvable. 2 h h A u I f A v = ( − ) h • To be certain numerical round-off does not perturb solvability, we incorporate a Gram-Schmidt like step 2 h each time a new right-hand side is generated for f the coarse grid: 2 h 2 h f 1 , 2 h 2 h 2 h f f 1 ← − 2 h 2 h 1 1 , 43 of 104
Neumann Problem: an example • Consider the problem , u x 2 x 1 − ″ ( ) = − u 0 u 1 0 < x 0 1 ′ ( ) = ′ ( ) = < 2 3 x x which has as a solution for any c u x c ( ) = − + 2 3 ( c=-1/12 gives the zero mean solution) . grid size average number | e h | r h | | | | | | N conv. factor of cycles 31 6.30E-11 0.079 9.70E-05 9 63 1.90E-11 0.089 2.40E-05 10 127 2.60E-11 0.093 5.90E-06 10 255 3.70E-11 0.096 1.50E-06 10 511 5.70E-11 0.100 3.70E-07 10 1027 8.60E-11 0.104 9.20E-08 10 2047 2.10E-11 0.112 2.30E-08 10 4095 5.20E-11 0.122 5.70E-09 10 44 of 104
Outline ✔ • Nonlinear Problems ✔ • Neumann Boundary Conditions • Anisotropic Problems • Variable Mesh Problems • Variable Coefficient Problems • Algebraic Multigrid 45 of 104
Anisotropic Problems − 1 • All problems considered thus far have had as h 2 the off-diagonal entries. • We consider two situations when the matrix has two different constant on the off-diagonals. These situations arise when – the (2-d) differential equation has constant, but different, coefficients for the derivatives in the coordinate directions – the discretization has constant, but different, mash spacing in the different coordinate directions 46 of 104
We consider two types of anisotropy • Different coefficients on the derivatives u u f − − α = x x y y discretized on a uniform grid with spacing h . • Constant, but different, mesh spacings: 1 h x h = = N h y h x h = y α h x 47 of 104
Both problems lead to the same stencil u 2 u u u 2 u u − + − − + − j 1 k j k j 1 k j k 1 j k j k 1 − , , + , , − , , + + α 2 2 h h − α 1 A h � � = 1 2 2 1 − + α − � � 2 h − α � � u 2 u u u 2 u u − + − − + − j 1 k j k j 1 k j k 1 j k j k 1 − , , + , , − , , + + 2 2 h h � � � � α � � 48 of 104
Why standard multigrid fails − α 1 • Note that has weak connections in the A h � � 1 2 2 1 = − + α − � � 2 h − α y -direction. MG convergence factors degrade as α gets � � small. Poor performance at α = 0.1 . 0 1 • Consider α � 0 . A h � � 1 2 2 1 − + α − � � � 2 h 0 � � • This is a collection of disconnected 1-d problems! • Point relaxation will smooth oscillatory errors in the x - direction (strong connections), but with no connections in the y -direction the errors in that direction will generally be random, and no point relaxation will have the smoothing property in the y-direction. 49 of 104
We analyze weighted Jacobi • The eigenvalues of the weighted Jacobi iteration matrix for this problem are λ i , l = 1 − 2 ω � i π � l π � sin 2 � + α sin 2 � � � � � � � � 1 + α � 2 N � 2 N � � � � 1 1 1 0.9 0.9 0.8 0.8 0.8 0.7 0.7 0.6 0.6 0.6 0.4 0.5 0.5 0.2 0.4 0.4 15 0 0 0.3 0.3 10 0 2 4 6 8 10 12 14 16 0 2 4 6 8 10 12 14 16 i 5 l l 5 10 i 15 0 l k 50 of 104
Two strategies for anisotropy • Semicoarsening Because we expect MG-like convergence for the 1-d problems along lines of constant y , we should coarsen the grid in the x - direction, but not in the y -direction. • Line relaxation Because the the equations are strongly coupled in the x -direction it may be advantageous to solve simultaneously for entire lines of unknowns in the x -direction (along lines of constant y ) 51 of 104
Semicoarsening with point relaxation − α 1 A h � � • Point relaxation on smooths in the x - 1 2 2 1 = − + α − � � 2 h − α � � direction. Coarsen by removing every other y -line. Ω h Ω 2 h • We do not coarsen along the remaining y -lines. • Semicoarsening is not as “fast” as full coarsening. The Ω 2 h number of points on is about half the number of Ω h points on , instead of the usual one-fourth. 52 of 104
Interpolation with semicoarsening • We interpolate in the 1-dimensional way along each line of constant y . • The coarse-grid correction equations are h h 2 h v v v = + 2 j k 2 j k j k , , , 2 h 2 1 h v v + j k j k , + , h h v v = + 2 j 1 k 2 j 1 k + , + , 2 53 of 104
Line relaxation with full coarsening • The other approach to this problem is to do full coarsening, but to relax entire lines (constant y ) of variables simultaneously. A h • Write in block form as D − cI � � − cI D − cI � � A h = � � − cI D − cI � � ⋅⋅⋅ ⋅⋅⋅ − cI � � � � − cI D � � where � 2 + 2 α − 1 α and � c − 1 2 + 2 α − 1 D = 1 = � � 2 ⋅⋅⋅ h � � h 2 � � − 1 2 + 2 α � � 54 of 104
Line relaxation • One sweep of line relaxation consists of solving a tridiagonal system for each line of constant y . h h v h • The kth such system has the form where D v g k = k k v h is the kth subvector of with entries and h h v v ( ) = k j j k , the kth right-hand side subvector is α h h h k h k g f v v ( ) = + ( + ) k j j k j 1 j 1 , , − , + 2 h • Because D is tridiagonal, the kth system can be solved very efficiently. 55 of 104
Why line relaxation works • The eigenvalues of the weighted block Jacobi iteration matrix are 2 i l ω π π � 2 � � 2 � � � 1 s i n s i n λ = − + α i l , � � � � � � 2 N 2 N i π 2 � � � � � � � � 2 s i n + α � � 2 N � � 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.2 0.4 0 0.2 0.2 -0.2 0 15 0 15 10 10 -0.2 5 -0.2 5 0 2 4 6 8 10 12 14 16 l 0 0 i 0 2 4 6 8 10 12 14 16 i l 56 of 104
Semicoarsening with line relaxation • We might not know the direction of weak coupling or it might vary. • Suppose we want a method that can handle either − 1 − α = 1 h = 1 or h � � � � A 2 − α 2 + 2 α A 1 − α − 1 2 + 2 α − 1 � � h 2 � � � h 2 � − 1 � − α � • We could use semicoarsening in the x-direction to h handle and line relaxation in the y-direction to A 1 h take care of . A 2 57 of 104
Semicoarsening with line relaxation • Original grid • The original grid • Coarsening is viewed as a done by deleting stack of every other “pencils.” Line pencil relaxation is used to solve problem along each pencil. 58 of 104
An anisotropic example • Consider with u=0 on the boundaries u u f − − α = x x y y of the unit square, and stencil given by − α 1 A h � � 1 2 2 1 = − + α − � � 2 h − α � � • Suppose so the exact f ( x , y ) = 2 ( y − y 2 ) + 2 α( x − x 2 ) solution is given by u ( x , y ) = ( y − y 2 ) ( x − x 2 ) • Observe that if α is small, the x -direction dominates while if α is large, the y -direction dominates 59 of 104
What is smooth error? • Consider α =0.001 and suppose point Gauss-Seidel is applied to a random initial guess. The error after 50 sweeps appears as: 0 . 1 2 0 . 1 0 . 0 8 0 . 0 6 0.1 0 . 0 4 0 . 0 2 0.08 0 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 Error along line of constant x 0.06 0.04 0.02 0 . 1 2 0 0 . 1 1 0 . 0 8 1 0 . 0 6 0.8 0.5 0.6 0 . 0 4 y 0.4 x 0.2 0 . 0 2 0 0 0 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 0 . 8 0 . 9 1 Error along line of constant y 60 of 104
We experiment with 3 methods • Standard V(2,1)-cycling, with point Gauss-Seidel relaxation, full coarsening, and linear interpolation • Semicoarsening in the x -direction. Coarse and fine grids have the same number of points in the y - direction. 1-d full weighting and linear interpolation are used in the x -direction, there is no y-coupling in the intergrid transfers • Semicoarsening in the x -direction combined with line relaxation in the y -direction. 1-d full weighting and interpolation. 61 of 104
With semicoarsening, the operator must change • To account for unequal mesh spacing, the residual and relaxation operators must use a modified stencil − α � � 2 � � h y � � � � � 2 − 1 2 + 2 α − 1 � � � A = � � 2 � 2 � 2 � � h x h x h y h x � � � � − 1 � � 2 � � h y � � h x • Note that as grids become coarser, grows while remains constant. h y 62 of 104
How do the 3 methods work for various values of α α ? α α • Asymptotic convergence factors: α scheme 1000 100 10 1 0.1 0.01 0.001 1E-04 V(2,1)-cycles 0.95 0.94 0.58 0.13 0.58 0.90 0.95 0.95 0.94 0.99 0.98 0.93 0.71 0.28 0.07 0.07 emicoarsening in x 0.04 0.08 0.08 0.08 0.07 0.07 0.08 0.08 semiC / line relax y -direction strong x -direction strong • Note: semicoarsening in x works well for α < .001 but degrades noticeably even at α = .1 63 of 104
A semicoarsening subtlety • Suppose α is small, so that semicoarsening in x is h x − 2 used. As we progress to coarser grids, gets h y − 2 small but remains constant. h x − 2 • If, on some coarse grid, becomes comparable α h y − 2 to , the problem effectively becomes recoupled in the y -direction. Continued semicoarsening can produce artificial anisotropy, strong in the y -direction. • When this occurs, it is best to stop semicoarsening and continue with full coarsening on any further coarse grids. 64 of 104
Outline ✔ • Nonlinear Problems ✔ • Neumann Boundary Conditions ✔ • Anisotropic Problems • Variable Mesh Problems • Variable Coefficient Problems • Algebraic Multigrid 65 of 104
Variable Mesh Problems • Non-uniform grids are commonly used to accommodate irregularities in problem domains • Consider how we might approach the 1-d problem − u ″( x ) = f ( x ) 0 < x < 1 u ( 0 ) = u ( 1 ) = 0 posed on this grid: x = 1 x = 0 x j + 1 x j − 1 x j x N x 0 66 of 104
We need some notation for the mesh spacing • Let N be a positive integer. We define the spacing interval between and : x j x j + 1 h j + 1 / 2 ≡ x j + 1 − x j j = 0 , 1 , ... , N − 1 h j + 1 / 2 x j + 1 x j − 1 x j x N x 0 67 of 104
We define the discrete differential operator • Using second order finite differences (and plugging through a mess of algebra!) we obtain this discrete representation for the problem: h u j + 1 h ) u j h u j − 1 h h + β j h − β j h h = f j 1 ≤ j ≤ N − 1 − α j + ( α j h = u N h = 0 u 0 • where 2 2 h = h = β j α j h j + 1 / 2 ( h j − 1 / 2 + h j + 1 / 2 ) h j − 1 / 2 ( h j − 1 / 2 + h j + 1 / 2 ) 68 of 104
We modify standard multigrid to accommodate variable spacing • We choose every second fine-grid point as a coarse-grid point h h x 2 j x h x N 0 2 h x 2 h 2 h x j x 0 N 2 / • We use linear interpolation, modified for the h v 2 h spacing. If , then for v h = I 2 h 1 ≤ j ≤ N / 2 − 1 2 h + h 2 j + 1 / 2 v j + 1 2 h h 2 j + 3 / 2 v j h h = v j 2 h v 2 j + 1 v 2 j = h 2 j + 1 / 2 + h 2 j + 3 / 2 69 of 104
We use the variational properties A 2 h to derive restriction and . T 2 h = 1 h ) 2 h A h I 2 h A 2 h = I h h I h 2 ( I 2 h Ω 2 h • This produces a stencil on that is similar, but not identical, to the fine-grid stencil. If the resulting system is scaled by , then ( h j − 1 / 2 + h j + 1 / 2 ) the Galerkin product is the same as the fine-grid stencil. • For 2-d problems this approach can be generalized readily to logically rectangular grids. However, for irregular grids that are not logically rectangular, AMG is a better choice. 70 of 104
Outline ✔ • Nonlinear Problems ✔ • Neumann Boundary Conditions ✔ • Anisotropic Problems ✔ • Variable Mesh Problems • Variable Coefficient Problems • Algebraic Multigrid 71 of 104
Variable coefficient problems • A common difficulty is the variable coefficient problem, given in 1-d by − ( a ( x ) u ′( x ) ) ′ = f ( x ) 0 < x < 1 u ( 0 ) = u ( 1 ) = 0 where a(x) is a positive function on [0,1] • We seek to develop a conservative, or self-adjoint, method for discretizing this problem. • Assume we have available to us the values of a(x) at the midpoints of the grid a j + 1 / 2 ≡ a ( x j + 1 / 2 ) x h x h x h x h x h 0 N j − 1 j 72 of 104 j + 1
We discretize using central differences • We can use second-order differences to approximate the derivatives. To use a grid spacing of h we evaluate a(x)u’(x) at points midway between the gridpoints: ( au ′) | xj + 1 / 2 − ( a u ′) | xj − 1 / 2 + O ( h 2 ) ( a ( x ) u ′( x ) ) ′ ≈ h x j Points used to evaluate (au’ ) ’ at x j x h x h x h x h x h 0 N j − 1 j j + 1 73 of 104
We discretize using central differences • To evaluate we must sample a(x) at the ( au ′) | x j + 1 / 2 point and use second order differences: x j + 1 / 2 u j − u j − 1 u j + 1 − u j ( a u ′ ) | x j − 1 / 2 ≈ a j − 1 / 2 ( a u ′ ) | x j + 1 / 2 ≈ a j + 1 / 2 h h where a j + 1 / 2 ≡ a ( x j + 1 / 2 ) Points used to Points used to evaluate (au’ ) ’ at x j evaluate u’ at x j + 1 / 2 x h x h x h x h x h 0 j − 1 N j j + 1 74 of 104
The basic stencil is given • We combine the differences for u’ and for (au’ ) ’ to obtain the operator u j + 1 − u j − 1 u j − u j − 1 − a j + 1 / 2 + a j − 1 / 2 h h − ( a ( x j ) u ′ ( x j ) ) ′ ( x j ) ≈ h and the problem becomes, for 1 ≤ j ≤ N − 1 1 h 2 ( − a j − 1 / 2 u j − 1 + ( a j − 1 / 2 + a j + 1 / 2 ) u j − a j + 1 / 2 u j + 1 ) = f j u 0 = u N = 0 75 of 104
Coarsening the variable coefficient problem • A reasonable approach is to use a standard multigrid algorithm with linear interpolation, full weighting, and the stencil 1 A 2 h = 2 h 2 h 2 h 2 h ( 2 h ) 2 ( − a j − 1 / 2 a j − 1 / 2 + a j + 1 / 2 − a j + 1 / 2 ) h h a 2 j + 1 / 2 + a 2 j + 3 / 2 where 2 h a j + 1 / 2 = 2 • The same stencil can be obtained via the Galerkin relation 76 of 104
Standard multigrid degrades if a(x) is highly variable • It can be shown that the variable coefficient discretization is equivalent to using standard multigrid with simple averaging on the Poisson problem on a certain variable-mesh spacing. -(a(x)u’)’ = f -u’’(x) = f • But simple averaging won’t accurately represent smooth components if is close to but far from . h h h x 2 j x 2 j + 1 x 2 j + 2 h h h x 2 j x 2 j + 2 x 2 j + 1 2 h 2 h x j x j + 1 77 of 104
One remedy is to apply operator interpolation • Assume that relaxation does not change smooth error, so the residual is approximately zero. Applying at yields h x 2 j + 1 h h h h h h h − a 2 j + 1 / 2 e 2 j + ( a 2 j + 1 / 2 + a 2 j + 3 / 2 ) e 2 j + 1 − a 2 j + 3 / 2 e 2 j + 2 = 0 h 2 • Solving for h e 2 j + 1 2 h + a 2 j + 3 / 2 h h 2 h a 2 j + 1 / 2 e j e j + 1 h e 2 j + 1 = h h a 2 j + 1 / 2 + a 2 j + 3 / 2 78 of 104
Thus, the operator induced interpolation is h = v j 2 h v 2 j 2 h + a 2 j + 3 / 2 h h 2 h a 2 j + 1 / 2 v j v j + 1 h v 2 j + 1 = h h a 2 j + 1 / 2 + a 2 j + 3 / 2 • And, as usual, the restriction and coarse-grid operators are defined by the Galerkin relations 2 h = c ( I 2 h h ) T 2 h A h I 2 h A 2 h = I h h I h 79 of 104
A Variable coefficient example • We use V(2,1) cycle, full weighting, linear interpolation. • We use and a ( x ) = ρ sin ( k π x ) a ( x ) = ρ rand ( k π x ) a ( x ) = ρ rand ( k π x ) a ( x ) = ρ sin ( k π x ) ρ k=3 k=25 k=50 k=100 k=200 k=400 0 0.085 0.085 0.085 0.085 0.085 0.085 0.085 0.25 0.084 0.098 0.098 0.094 0.093 0.083 0.083 0.5 0.093 0.185 0.194 0.196 0.195 0.187 0.173 0.75 0.119 0.374 0.387 0.391 0.39 0.388 0.394 0.85 0.142 0.497 0.511 0.514 0.514 0.526 0.472 0.95 0.191 0.681 0.69 0.694 0.699 0.745 0.672 80 of 104
Outline ✔ • Nonlinear Problems ✔ • Neumann Boundary Conditions ✔ • Anisotropic Problems ✔ • Variable Mesh Problems ✔ • Variable Coefficient Problems • Algebraic Multigrid 81 of 104
Algebraic multigrid: for unstructured-grids ● Automatically defines coarse “grid” ● AMG has two distinct phases: — setup phase: define MG components — solution phase: perform MG cycles ● AMG approach is opposite of geometric MG — fix relaxation (point Gauss-Seidel) — choose coarse “grids” and prolongation, P , so that error not reduced by relaxation is in range(P) — define other MG components so that coarse- grid correction eliminates error in range(P) (i.e., use Galerkin principle) (in contrast, geometric MG fixes coarse grids, then defines suitable operators and smoothers) 82 of 112
AMG has two phases: • Setup Phase Ω m + 1 m – Select Coarse “grids,” 1 2 , = , , . . . I m m 1 2 , = , , . . . – Define interpolation, m + 1 – Define restriction and coarse-grid operators T 1 = m 1 m m m 1 m m + + + I I A I A I = ( ) m m m 1 m 1 + + ● Solve Phase — Standard multigrid operations, e.g., V-cycle, W-cycle, FMG, FAS, etc ● Note: Only the selection of coarse grids does not parallelize well using existing techniques! 83 of 112
AMG fundamental concept: Smooth error = “small” residuals • Consider the iterative method error recurrence 1 − k 1 k + e I Q A e = ( − ) • Error that is slow to converge satisfies Q − 1 − 1 I Q A e e A e 0 ( − ) ≈ ≈ � r 0 ≈ � • More precisely, it can be shown that smooth error satisfies ( 1 ) r D − 1 « e A 84 of 112
AMG uses strong connection to determine MG components • It is easy to show from (1) that smooth error ( 2 ) satisfies A e e « D e e , , • Define i is strongly connected to j by a m a x a 0 1 < − ≥ θ { − } , θ ≤ i j i k k i ≠ • For M-matrices, we have from (2) a e e − − 2 1 � i j i j � � � � « 1 e 2 � 2 a � � � i i i i ≠ j � � � � • implying that smooth error varies slowly in the direction of strong connections 85 of 112
Some useful definitions u i • The set of strong connections of a variable , that is, the variables upon whose values the value of depends, is defined as u i � � S i j a m a x a = : − > θ − i j i j � � j i ≠ � � • The set of points strongly connected to a variable u i S T is denoted: . j j S = { : ∈ } i i • The set of coarse-grid variables is denoted C . • The set of fine-grid variables is denoted F . • The set of coarse-grid variables used to C i interpolate the value of the fine-grid variable u i is denoted . 86 of 112
Choosing the Coarse Grid • Two Criteria j ∈ S i ∈ F – (C1) For each , each point should either be i C in or should be strongly connected to at least one point C i in C – (C2) should be a maximal subset with the property that no two -points are strongly connected to each C other. • Satisfying both (C1) and (C2) is sometimes impossible. We use (C2) as a guide while enforcing (C1). 87 of 112
Selecting the coarse-grid points C-point selected (point with largest “value”) Neighbors of C-point become F- points Next C-point selected (after updating “values”) F-points selected, etc . 88 of 112
Examples: Laplacian Operator 5-pt FD, 9-pt FE (quads), and 9-pt FE (stretched quads) 9-pt FE (stretched quads) 1 1 1 1 − − − − � � � � 1 4 1 1 8 1 − − − − � � � � 1 − 1 1 1 � � − − − � � 1 4 1 − − − 9-pt FE (quads) 5-pt FD � � 2 8 2 � � 1 4 1 − − − � � 89 of 112
Prolongation is based on smooth error, strong connections (from M-matrices) Smooth error is given by: C r a e a e 0 = + ≈ F i i i i i j j � F j N ∈ i Prolongation : i C C e i C , ∈ i � � P e ( ) = e i F F ω , ∈ � i k k i � � k C ∈ i � 90 of 112
Prolongation is based on smooth error, strong connections (from M-matrices) Sets: C i C Strongly connected -pts. D s C F F Strongly connected -pts. i F D w Weakly connected points. i i The definition of smooth error, C C a e a e ≈ − i i i i j j � j i ≠ Gives: F a e a e a e a e ≈ − − − i j j i j j i j j i i i � � � s w j C ∈ j D j D ∈ ∈ i i i Strong C Strong F Weak pts. 91 of 112
Finally, the prolongation weights are defined • In the smooth-error relation, use for weak e e j = i connections. For the strong F -points use : � � � � e a e a = ⁄ j j k k j k � � � � � � k C k C ∈ ∈ � � � � i i yielding the prolongation weights: a a i k k j a + i j � a s k m i � j D ∈ m C ∈ i w = − i j a a + i i i n � w n D ∈ i 92 of 112
AMG setup costs: a bad rap • Many geometric MG methods need to compute prolongation and coarse-grid operators • The only additional expense in the AMG setup phase is the coarse grid selection algorithm • AMG setup phase is only 10-25% more expensive than in geometric MG and may be considerably less than that! 93 of 112
AMG Performance: Sometimes a Success Story • AMG performs extremely well on the model problem (Poisson’s equation, regular grid)- optimal convergence factor (e.g., 0.14) and scalability with increasing problem size. • AMG appears to be both scalable and efficient on diffusion problems on unstructured grids (e.g., 0.1- 0.3). • AMG handles anisotropic diffusion coefficients on irregular grids reasonably well. • AMG handles anisotropic operators on structured and unstructured grids relatively well (e.g., 0.35). 94 of 112
So, what could go wrong? Strong F-F connections: weights are dependent on each other i • For point the value is interpolated from , , e j k 1 k 2 and is needed to make the interpolation weights for e i approximating . j • For point the value is interpolated from , , e i k 2 k 1 and is needed to make the interpolation weights for approximating . e j • It’s an implicit system! k C ∈ 1 i k C ∈ 1 j i j k 4 ∈ C k 3 ∈ C j i k C ∈ 2 i 95 of 112 k C ∈ 2 j
Is there a fix? • A Gauss-Seidel like iterative approach to weight definition is implemented. Usually two passes suffice. But does it work? ● Frequently, it does: Convergence factors for Laplacian, stretched quadrilaterals theta Standard Iterative 0.25 0.47 0.14 x 1 0 y ∆ = ∆ 0.5 0.24 0.14 0.25 0.83 0.82 x 1 0 0 y ∆ = ∆ 0.53 0.5 0.23 96 of 112
AMG for systems • How can we do AMG on systems? A A u f 1 1 1 2 � � � � � � = g v � � � � � � A A 2 1 2 2 � � � � � � • Naïve approach: “Block” approach (block Gauss-Seidel, using scalar AMG to “solve” at each cycle) 1 − u A f A v ← ( ) ( − ) 1 1 1 2 1 − v A g A u ← ( ) ( − ) 2 2 2 1 ● Great Idea! Except that it doesn’t work! (relaxation does not evenly smooth errors in both unknowns) 97 of 112
AMG for systems: a solution • To solve the system problem, allow interaction between the unknowns at all levels: k k k I 0 ( ) A A k 1 + 1 1 1 2 u k � � k � � A I and = = k 1 + � � � � k k k 0 I A A ( ) k 1 + 2 1 2 2 � � � � v • This is called the “unknown” approach. • Results: 2-D elasticity, uniform quadrilateral mesh: m esh spacing 0.125 0.0625 0.03135 0.015625 Convergence factor 0.22 0.35 0.42 0.44 98 of 112
How’s it perform (vol I)? Regular grids, plain, old, vanilla problems • The Laplace Operator: Convergence Time Setup Stencil per cycle Complexity per Cycle Times 5-pt 0.054 2.21 0.29 1.63 5-pt skew 0.067 2.12 0.27 1.52 9-pt (-1,8) 0.078 1.30 0.26 1.83 9-pt (-1,-4,20) 0.109 1.30 0.26 1.83 U U − ε − • Anisotropic Laplacian: x x y y Epsilon 0.001 0.01 0.1 0.5 1 2 10 100 1000 Convergence/cycle 0.084 0.093 0.058 0.069 0.056 0.079 0.087 0.093 0.083 99 of 112
How’s it perform (vol II)? Structured Meshes, Rectangular Domains • 5-point Laplacian on regular rectangular grids Convergence factor (y-axis) plotted against number of nodes (x-axis) 0 .1 6 0 .1 4 6 0 .1 4 0 .1 4 2 0 .1 3 5 0 .1 3 3 0 .1 2 0 .1 0 .0 8 0 .0 8 2 0 .0 6 0 .0 5 5 0 .0 4 0 .0 2 0 0 5 0 0 0 0 0 1 0 0 0 0 0 0 100 of 112
Recommend
More recommend
