Kernel Preserving Multigrid Methods for Convection-Diffusion Equations

We propose a kernel preserving multigrid approach for solving convection-diffusion equations. The multigrid methods use Petrov-Galerkin coarse grid correction and linear interpolation. The restriction operator is constructed by preserving the kernel of the convection-diffusion operator. The construction considers constant and variable coefficient problems as well as cases where the convection term is not known explicitly. For constant convection-diffusion problems, we prove that the resulting Petrov-Galerkin coarse grid correction has small phase errors and the coarse grid matrix is essentially an M-matrix. The algebraic construction of the interpolation and restriction operators is a generalization of the energy minimization approach for elliptic PDEs.

Theoretical results

Denote the iteration matrix of multigrid (2-level) as follows:

where P = prolongation and R = restriction. The Fourier symbol of the coarse grid correction matrix C is:

We are interested in small frequencies |μ|, |ν| ≈ 0. If one ignores the effect of P & R, then

Denote the strength of convection by β. As β → ∞, the coarse and fine grid matrices differ by

Thus, the coarse grid matrix shifts waves of all frequencies by 1 grid point to the right. (Similar, it shifts by 1 grid point to the left if β → -∞. Hence, the coarse grid error is off by 1 grid point, resulting in at best first order accuracy.

If one considers the effect of the kernel preserving restriction, then

Thus the kernel preserving restriction counteracts the shifting caused by the coarse grid matrix.

To summarize, the phase errors of various coarse grid matrices:

	Petrov-Galerkin	Direct discretization	Galerkin
Cross-characteristic components (μ = ν)	O(μπh)^2	O(μπh)	O(μπh)^3
Characteristic components (μ = -ν)	O(μπh)^4	1/2	O(μπh)^2

(Note: Galerkin coarse grid matrix is often unstable for strong convection.)

Stability analysis

Consider the Petrov-Galerkin coarse grid matrix A^H = R A^h P, R = kernel preserving restriction, P = linear interpolation. Assume A^h is an M-matrix (e.g. Scharfetter-Gummel discretization). Then

Thus, for |β| is sufficiently large, then A^H is numerically an M-matrix.

Numerical results

BMG = Blackbox multigrid
EMG(Jac) = our multigrid with Jacobi smoothing
EMG(GS1) = our multigrid with GS smoothing (order first in x (left to right) and then y (bottom to top))
EMG(GS2) = our multigrid with GS smoothing (order first in y (top to bottom) and then x (left to right))
EMG(GS3) = our multigrid with GS smoothing (order first in x (right to left) and then y (top to bottom))
EMG(GS4) = our multigrid with GS smoothing (order first in y (bottom to top) and then x (right to left))

Number of multigrid V-cycles, h = 1/32.

δ	BMG	EMG(Jac)	EMG(GS1)	EMG(GS2)	EMG(GS3)	EMG(GS4)
1	11	11	6	6	8	6
1/2	13	10	5	7	9	7
1/4	20	10	4	6	11	6
1/8	>20	9	3	8	15	8
1/16	>20	9	3	10	18	10

EMG(GS) results are generally worse then those for EMG(Jac), except for EMG(GS1), where the ordering is consistent with the flow direction. The worst results are given by EMG(GS3) whose ordering if opposite to the flow direction. Thus, in contrast with diffusion dominated problems, Gauss-Seidel shows a poorer smoothing effect compared to damped Jacodi, but can still be an effective smoother if the ordering happens to be consistent with the flow direction.

Number of multigrid V(2,2)-cycles, δh = 10^-3.

h	EMG(Jac)	EMG(GS1)	EMG(GS2)	EMG(GS3)	EMG(GS4)
1/16	9	2	6	11	6
1/32	9	2	10	18	10
1/64	9	3	15	28	15
1/128	9	3	12	14	12

Example 2: recirculating flow

a = 4x(x-1)(1-2y), b = -4y(y-1)(1-2x).

Number of multigrid V(2,1)-cycles, ε = 10^-3.

h	EMG(βmidpt)	EMG(βave)	EMG(βweighted)
1/16	6	6	4
1/32	7	7	6
1/64	10	10	8
1/128	12	12	10

EMG(β_midpt), EMG(β_ave), and EMG(β_weighted) are 3 variations of our multigrid methods.

Example 3: Asian option pricing

S = stock price
A = arithmetic average of the stock price over some time interval T
τ = time from the expiration date T
r = interest rate
σ = volatility

Note: (1) There is no diffusion in the A direction. (2) The convection approaches infinity as τ approaches the expiration date T.

Analytical and numerical solutions for an Asian call option with zero strike and the average number of multigrid V(2,2)-cycles over all time steps.

σ	T	Analytic	Numerical	EMG(GS)
0.1	0.25	98.763	98.7604	12
	0.50	97.547	97.5411	9
	1.00	95.175	95.1626	6
0.2	0.25	98.763	98.7574	12
	0.50	97.547	97.5309	9
	1.00	95.175	95.1785	6
0.4	0.25	98.763	98.7673	12
	0.50	97.547	97.5447	9
	1.00	95.175	95.1633	6