ISSN 09655425, Computational Mathematics and Mathematical Physics, 2015, Vol. 55, No. 8, pp. 1276–1289. © Pleiades Publishing, Ltd., 2015. Original Russian Text © V.T. Zhukov, N.D. Novikova, O.B. Feodoritova, 2015, published in Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki, 2015, Vol. 55, No. 8, pp. 1305–1319.
In blessed memory of Professor A.P. Favorskii
On the Solution of Evolution Equations Based on Multigrid and Explicit Iterative Methods V. T. Zhukov, N. D. Novikova, and O. B. Feodoritova Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, Miusskaya pl. 4, Moscow, 125047 Russia email:
[email protected],
[email protected],
[email protected] Received February 26, 2015
Abstract—Two schemes for solving initial–boundary value problems for threedimensional parabolic equations are studied. One is implicit and is solved using the multigrid method, while the other is explicit iterative and is based on optimal properties of the Chebyshev polynomials. In the explicit iter ative scheme, the number of iteration steps and the iteration parameters are chosen as based on the approximation and stability conditions, rather than on the optimization of iteration convergence to the solution of the implicit scheme. The features of the multigrid scheme include the implementation of the intergrid transfer operators for the case of discontinuous coefficients in the equation and the adaptation of the smoothing procedure to the spectrum of the difference operators. The results pro duced by these schemes as applied to model problems with anisotropic discontinuous coefficients are compared. DOI: 10.1134/S0965542515080151 Keywords: threedimensional parabolic equations, anisotropic discontinuous coefficients, multigrid method, explicit iterative scheme with Chebyshev parameters.
1. INTRODUCTION Below, we analyze the efficiency of two schemes for solving initial–boundary value problems for three dimensional parabolic equations. The need for efficient solutions of such problems is caused by their abundance in mathematical models. Increased interest in this problem is motivated by the development of parallel computations: codes for supercomputers with a large number of processors are required. Such codes must ensure scalable simulation on grids with billion (and more) nodes. By applying a spatial discretization method, the original nonstationary problem can be reduced to a system of ordinary differential equations. We consider two schemes for integrating such a system. One is implicit twolevel and is solved by applying the multigrid algorithm. Hereafter, it is referred to as the MM scheme (i.e., a scheme based on the multigrid method). A detailed presentation of the multigrid algo rithm, which is a version of Fedorenko’s method [1–3], can be found in [4–8]. The other scheme is the explicit iterative LIM one [9], which is designed on the basis of the scheme from [10] for solving parabolic equations and relies on Chebyshev polynomials of special structure. Below, we analyze the performance of the MM and LIM schemes as applied to evolution problems with discontinuous coefficients. 2. FORMULATION OF THE PROBLEM Consider an initial–boundary value problem for the parabolic equation
∂u + Lu = f , ∂t
r ∈ G,
(1)
where r = (t, x, y, z) ∈ G = [t 0;T ] × Ω, Ω ⊂ R 3 , and [t 0;T ] is a given time interval. A boundary condition (Dirichlet, Neumann, or Robin; see [8] for details) is set on the lateral boundary of G, and an initial con dition is specified at t = t 0 . In (1) L is a linear elliptic selfadjoint positive semidefinite operator, f (r ) is a given function, and u(r ) is the sought function. The input data are assumed to ensure the existence and uniqueness of a solution of the required smoothness. 1276
ON THE SOLUTION OF EVOLUTION EQUATIONS
1277
Various generalizations are possible, but the main results will be described for the case
Lu = −div( κ grad u) + a0u
(2)
with piecewise smooth nonnegative functions κ(r ) and a0(r ). For simplicity, let Ω be a threedimensional right parallelepiped. In [t 0,T ] × Ω we introduce a grid Ω h,τ = Ωh × Ωτ, where Ω τ = {t j , 0 ≤ j ≤ J , t J = T } is a grid in time with a variable step τ = t j +1 − t j > 0 and Ω h = {x n ∈ Ω, 0 ≤ n ≤ N } is a Cartesian grid that is nonuniform in each coordinate direction in Ω and depends on a parameter h (mesh size) characterizing the mean cell size. The space U h of functions defined on Ωh and the inner product and the norm in U h are defined in the standard manner (see [8]). The value of the grid function u at the j th time level t j is denoted by u j . On the subspace U h , we define a difference operator Lh approximating the operator L with the second order of accuracy for smooth functions in view of the boundary condition. To be more specific, it can be assumed that the grid functions are given at the nodes of the grid Ωh and Lh is constructed with the help of sevenpoint finitevolume discretization. The operator Lh is selfadjoint, and its eigenvalues λ are nonnegative and belong to the real interval [λ min; λ max ]. The basic information concerning the spectra of simplest elliptic difference operators can be found in [11]. Assume that the estimates for the boundaries λ min ≥ 0 and λ mах of the spectrum of Lh are known; in some cases, they can be calculated (see [6]). 3. FEATURES OF SOLUTIONS TO EVOLUTION PROBLEMS The semidiscrete approximation of problem (1) is written in operator form (taking into account the boundary conditions in the definition of the operator Lh ):
∂u + L u = f . h ∂t The transition from the level t j to the upper level t j +1 = t j + τ can be implemented in different manners. Consider the explicit scheme
u j +1 − u j + Lhu j = f j , τ which imposes a severe restriction on the time step, and the fully implicit scheme
(3)
u j +1 − u j + Lhu j +1 = f j , τ which can be written as a system of linear equations
(4)
( I + τLh ) u j +1 = u j + τf j , where I is the difference identity operator. This system can be written as
Ahuh = g h,
(5)
where Ah = I + τLh is an N × N matrix, uh is the desired grid function, and g h is a given grid function. In contrast to stationary problems, a new aspect in evolution problems is the presence of moving fronts in the solution, which are generated, for example, in heat transfer in a domain consisting initially of dif ferent parts: a cold subdomain with low thermal conductivity and a hot subdomain with high thermal con ductivity. In stationary computations, there are no thermal fronts, since the temperature in the domains becomes steadystate over an infinite time. Another aspect in the solution of evolution problems concerns the accuracy of time integration. There are no many integration schemes to choose from. Splitting schemes over spatial variables are worse paral lelized and can be unsuitable in some cases. Explicit schemes are universal, but they are used extremely rarely because of the wellknown time step constraint. The application of an implicit scheme does not guarantee that the time evolution of the solution is described to high accuracy and, additionally, they lead to systems of linear algebraic equations of form (5), which are difficult to solve for multidimensional prob lems. These systems are usually solved by an iterative method, but it should be kept in mind that the num ber of iterations must satisfy a stability condition. A necessary constraint on the number of iterations is given by the Gelfand–Lokutsievskii theorem [12, p. 275–282] (see also [13, p. 793]). The class of schemes COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
1278
ZHUKOV et al.
falling under this constraint includes iterative schemes equivalent to some twolevel explicit scheme with an extended stencil of influence (see Section 4). For this reason, the LIM scheme appears as a competitor of the MM scheme. The Gelfand–Lokutsievskii theorem holds for LIM and, additionally, it takes into account a wellknown characteristic property of difference approximations to elliptic operators stated by Lusternik in [14]: The first smooth eigenelements approximate the eigenfunctions of a differential opera tor, while the last eigenelements are spurious and depend basically on the properties of the grid and the grid approximation. This implies the wellknown principle that the time evolution the first, lowfrequency components of the discrete solution has to be computed correctly, while, for highfrequency components, it is sufficient to ensure their boundedness for stability. The LIM scheme satisfies this principle. According to the LIM scheme, the computation of the solution at the upper time level consists of executing a partic ular number of explicit iterations with parameters that are the roots of Chebyshev polynomials of the first kind. Chebyshev parameters began to be used in solving evolution equations in the 1980, but some of the basic ideas were formed in the 1950. Once again, we emphasize that the principal idea of the LIM scheme differs from the usual Chebyshev acceleration of iterations: the choice of the number of iterations and iter ation parameters is determined by the approximation and stability conditions rather than by optimizing the convergence of iterations to the solution of an implicit scheme. For parabolic equations, the multigrid method develops in the following direction. After writing implicit scheme (4), system (5) is solved for the unknown function at the upper level by applying the mul tigrid algorithm. The multigrid algorithm uses standard Chebyshev iterations [11] to solve the equations on the coarsest grid and at smoothing stages. As a smoother, we can use a modification of the LIM scheme that ensures the smoothing of the residual in the highfrequency domain of the spectrum of Lh . The fact is that the LIM scheme not only restricts the growth of the harmonic amplitudes to ensure stability, but also guarantees a substantial decrease in the highfrequency components, which agrees with the needs of the multigrid method. In addition to powerful smoothers, algorithmic implementations of the intergrid transfer operators have been constructed in problemdependent form [8], and we analyze their perfor mance in evolution problems. Finally, the basic algorithmic elements of the multigrid method ensure its effective parallel implementation (see [6–8]). For the MM scheme, a new resource of improving its effi ciency is studied, namely, in the case of smoother adaptation, the lower boundary of the highfrequency spectrum obtained at the first time step can be used as a constant at all the subsequent time steps. The complexity of solving parabolic equations is characterized by the “parabolic” Courant number cou = O(τ h − 2 ). For threedimensional problems, the computational cost of the MM scheme at a single −3 time step is O(h ) and, for a fixed grid, is independent of the Courant number cou. However, the constant −3 in O(h ) depends on the resolution of the implicit scheme. In the threedimensional case, the amount of
computations in the LIM scheme is O problem.
(
)
τh −2h −3 and the constant in this estimate depends only on the
Advantages of the MM scheme are that it performs in any range of Courant numbers and is asymptot ically optimal (the number of multigrid iterations required for achieving the prescribed accuracy is inde pendent of h). These advantages are especially pronounced at high Courant numbers. For relatively small cou ≈ 10 2 −10 4 , the LIM scheme becomes competitive (in terms of amount of computations and the accu racy, see Section 5). A distinctive feature of LIM that makes it more competitive is the absence of tuning parameters. The LIM scheme is described in detail in [9, 15]. A brief description is given in the following section. 4. THE LIM LOCAL ITERATION SCHEME In this scheme, the transition operator from the lower level t j to the upper level t j +1 = t j + τ is deter mined by a Chebyshev operator polynomial F p (Lh ) of special structure [10] of algebraic degree p = ⎡π τλ max + 1⎤ , ⎢4 ⎥
(6)
where ⎢⎡ x ⎥⎤ denotes the least integer greater than or equal to x. The polynomial F p is a firstkind Chebyshev polynomial [16] that deviates least from zero on the interval [λ 0; λ max ] with the normalization condition COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
ON THE SOLUTION OF EVOLUTION EQUATIONS
1279
F p (− 1 τ) = 1. Here, λ max is the maximum eigenvalue of the operator Lh and the left end point λ 0 of the interval is determined by the parameters p and λ max using the formulas (see [9, 10])
λ 0 = λ max
z1 − 1 ∈ [−1/ τ, 0] , z1 + 1
z1 = cos π . 2p
The polynomial F p can easily be expressed in terms of the standard Chebyshev polynomial: m= p
F p(λ) = H p(λ) H p(− 1 τ),
H p(λ) =
∏(a
m
− λ) ≡ T p (z1 − (z1 + 1)λ λ max ) ,
(7)
m=1
where T p(z) = cos( p arccos z) is a Chebyshev polynomial of degree p on the interval [–1; 1], to which the interval [λ 0; λ max ] is mapped under the substitution z = z1 − (z1 + 1)λ/ λ max . To construct F p , it is sufficient to set the parameters
am =
λ max (z1 − β m), 1 + z1
m = 1, ..., p ,
(8)
which are expressed in terms of the zeros β m of the optimal Chebyshev polynomial T p , which are ordered for stability [11, 17], but with the preservation of β1 = z1 = cos(0.5π/ p):
⎧ ⎫ β m ∈ K p = ⎨cos 2i − 1 π, i = 1, 2, ..., p⎬ . 2p ⎩ ⎭ The transition from the grid function u j to u j +1 is executed in ν = 2 p − 1 explicit steps
ym =
1 [u + τb y m−1 − τ(L y m−1 − f )], j m h h 1 + τbm
(9)
m = 1, …, ν ,
(10)
where y 0 = u j is the initial approximation and the result is the grid function at the upper level: u j +1 = y ν ; the set of iteration parameters bm is determined as {b1, b2, …, b2 p−1} ≡ {a1, a2, …, a p , a2, a3, …, a p}. For the error δ j = u j − u(t j , x, y, z), the transition from t j to the upper level t j +1 = t j + τ can be written as δ j +1 = W δ j , where the leveltolevel transition operator is a rational function: W = (I − F p2 )(I + τLh )−1 . It is convenient to explain the LIM scheme in terms of the “decrement” function
(11)
1 − F p2(λ) , (12) 1 + τλ which describes the damping of the harmonic amplitudes over a single time step; here, λ ∈ [0; λ max ] is a real interval containing the spectrum of Lh . On the initial segment of the spectrum, the function ρLIМ(λ) approximates the decrement function of the exact leveltolevel transition operator exp(−τ Lh ). Note that the transition operator of implicit scheme (4) has the form (I + τLh )−1, and which of the two approximations of the exact operator exp(−τ Lh ) is better remains an open question and depends on the time step and a particular problem. The LIM scheme requires an estimate for the maximum eigenvalue λ max of Lh , which is obtained by applying the Gershgorin circle theorem (see [18]). The important inequality 0 ≤ ρLIM(λ) ≤ (1 + τλ)−1 holds. In the computations, the parameters bm in (10) can be used in reverse order. Then the last iteration cor responds to an explicit scheme, which agrees with the predictor–corrector structure of difference schemes based on conservation laws. Note also that the lower boundary of the operator spectrum does not partic ipate in the definition of the algorithm. Note that the LIM scheme satisfies the Gelfand–Lokutsievskii theorem mentioned in Section 2, 2 according to which, if the scheme converges as h → 0 (assuming that τ/ h = const ), the necessary num 0.5 ber of iterations is estimated from below by popt = ( 0.5τλ max ) . The LIM scheme executes ν = 2 p − 1
ρLIM(λ) =
COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
1280
ZHUKOV et al. 1.0
LIM(p) Implicit scheme LIM(2p) Exact values Explicit scheme
0.5
0 0
0.05
0.10
0.15
0.20
0.25
Fig. 1. Decrement functions of the schemes on the initial spectral segment.
explicit steps, which is twice as large as the necessary number. In the base scheme LI [10], the number of explicit steps is close to the optimal one according to formula (6), but the solution produced by this scheme can exhibit noticeable spurious oscillations. The scheme with an “optimal” number of iterations popt determines the leveltolevel transition operator as a Chebyshev polynomial of form (7) with param eters am = (1 − β m ) λ max 2, m = 1, …, p . Such a scheme degenerates into a trivial explicit one (see [15]) with the spatial step increased by popt times. In other words, as τ grows, the stability of the “optimal” scheme is preserved, but its accuracy is worsened, since zero coefficients correspond to interior nodes of the differ ence influence stencil in such a scheme. The LI and LIM schemes are free of this shortcoming and, according to the results of [9] concerning the study of the decrement and source functions of these schemes, LI and LIM (especially, the latter) are efficient as applied to evolution problems. The MM scheme has optimality typical of elliptic equations: for fixed τ and prescribed accuracy of the solution to system (5), as h → 0 , the number of multigrid iterations in implicit scheme (4) is independent of h. This advantage, which is substantial in time marching to a steady state, is less important in evolution problems, where the growth of the ratio τ h 2 as h → 0 has to be bounded. The LIM scheme was designed for evolution problems and is nearly optimal in its class of methods in the sense of the Gelfand– Lokutsievskii theorem. Additionally, it has an important quality: for τ ≤ τ exp satisfying the stability condi tion τλ max ≤ 2 for the explicit scheme LIM automatically passes into the explicit scheme ensuring mini mum computational costs at a single time step. For illustrative purposes, Fig. 1 shows the decrement functions of the schemes on the initial spectral segment corresponding to smooth modes. The horizontal axis represents the initial portion [0; 0.25] of the normalized spectral segment [0; λ max h 2] = [0;12]. The decrements are shown for the exact scheme −1 exp(−τλ) (solid line), the implicit scheme (1 + τλ) (dashed line), LIM with the “standard” value p = 10 (dashdotted) and doubled p (circles), and for the explicit scheme ρ = 1 − τλ (dashed line with squares). The graph of the explicit scheme is the tangent to the graph of exp(−τλ) at the origin. This graph goes beyond the figure, and its ordinate at the endpoint of the spectral segment is equal to –150. It can be seen that the decrement functions of the implicit scheme and LIM for doubled p nearly coincide. For the stan dard value of p, the polynomial F p reaches extremal values equal to ± 1. At these points, the function −1 ρLIM(λ) vanishes. This value can be closer to the exact value exp(−τλ) than the value (1 + τλ) of the decrement function for the implicit scheme. For doubled p, the polynomial F p reaches extremal val −1 ues equal to ±π 24. Therefore, the maximum difference between ρLIM(λ) and (1 + τλ) is less than (π/24) 2 ≈ 1.7 × 10 − 2 ; hence, their graphs are indistinguishable. 5. BASIC ELEMENTS OF THE MULTIGRID METHOD Let us briefly describe the basic elements of the multigrid method (MM). To solve implicit scheme (4) represented as a system of grid equations (5), it is used in a form close to Fedorenko’s original method [1–3]. The number of grid levels in the MM algorithm is usually rather small (five or six). Every multigrid itera COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
ON THE SOLUTION OF EVOLUTION EQUATIONS
1281
tion step consists in the transition from a fine grid to one of the next level up to the coarsest grid and back (the socalled Vcycle). It is convenient to describe the basic elements of the method in a doublegrid rep resentation, when there are only two grids: fine h and coarse H. Then the iterating operator of the multi grid method has the form [19]:
Q = S p(I − PAH−1RAh )S p .
(13)
Here, Ah and AH are the operators on the fine and coarse grids, respectively; P and R are the intergrid transfer operators from the coarse to fine grid (interpolation P) and from the fine to coarse grid (projection R); and S p is a smoothing operator with p pre and postsmoothing steps. In the algorithm under consider ation, the operator AH is constructed via rediscretization, i.e., an approximation of the original equation with homogeneous boundary conditions is written on H. The righthand side of the system of equations on H is specified as the projection of the residual rh = g h − Ahuh onto H with the help of the operator R . The efficiency of the algorithm depends on the multigrid triad: the intergrid transfer operators, the algo rithm for solving the coarsest grid equations, and the smoothing operator S p . The algorithmic structures of the intergrid transfer operators P and R = P* are given in [8], where, along with a trilinear interpolation operator, an interpolation operator P based on an approximate solu tion of local difference boundary value problems, is presented, which is critically important for the case of discontinuous coefficients. Such operators P and R = P* are called problemdependent (see [19]). On the coarsest grid, equations of the form AH y = g H with a residual g H on the righthand side are solved using the standard Chebyshev iterative method [11] with iteration parameters ω j : y j = y j −1 − ω j ( AH y j −1 − g H ),
j = 1, …, p,
(14)
where y is the initial approximation, usually, y ≡ 0 . The number p of iterations is determined by the condition of achieving the prescribed accuracy according the criterion || rp|| < ε|| r0|| , where r0 and r p = g H − AH y p are the initial and terminal residuals. An estimate for p has the form (see [11]) 0
0
p = p(ε, η) ≈
(
),
ln ε −1 + ε −2 − 1 ln ρ
ρ=
1+ η , 1− η
η=
λ min , λ max
(15)
where λ min and λ max are estimated minimum and maximum eigenvalues of the difference operator AH . The smoothing operator S p specifies the transition from an iterative approximation v to a “smoother” approximation v new . Two smoothers S p = S p( Ah ) were constructed and examined in [4–8]. The first is an operator Chebyshev polynomial of degree p, while the second smoother is a rational function of form (11). Each smoother is a selfadjoint operator, and its eigenvalues are the values of the function S p(λ) on the spectrum of Lh . The smoothers serve to suppress the residual components corresponding to the highfre quency spectral range [λ*min; λ max ], where λ*min is a formal boundary separating the low and highfrequency spectral ranges. Note that the parameters λ*min and λ max are different for each grid level. In some cases, the interface between the spectral ranges can be estimated. For example, in the isotropic threedimensional case, λ*min = λ max 6. The interface can be corrected in the course of multigrid iterations. Such an adapta tion algorithm is under experimental verification (see [7, 8]) and is described below in this section. Recall the basic information about smoothers (see [8]). When a smoother is applied, the number of steps is chosen so as to suppress the highfrequency residual components the prescribed number of times ε −1 (usually, ε = 0.5 ). The computations are continued until the norm of the entire residual of the iterative approximation is decreased properly. Chebyshev polynomial smoother. Consider explicit iterations of form (14) with a Chebyshev polynomial F p that deviates least from zero on the highfrequency spectral segment [λ*min; λ max ]. Let p = p(ε, η) be defined by formulas (15) with a number p and prescribed accuracy ε = ε smooth . Then the error components on the −1 interval [λ*min; λ max ] are decreased by ε > 1 times uniformly over this spectral range and the smoother’s decrement function ρCheb(λ) = F p(λ) takes the values ±ε at the extremal points of the polynomial on [λ*min; λ max ]. COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
1282
ZHUKOV et al.
LIM smoother. Define p = p(ε, η) by the formula −1 p = ⎡π (ε − 1) η + 1⎤ ⎢4 ⎥
(16)
with η = λ*min λ max and prescribed accuracy ε = ε smooth . Then, to reduce the error by ε interval [λ*min; λ max ], we use the smoothing operator
−1
> 1 times on the
S ν = (I − G p2 )(I + τAh )−1,
(17)
where ν = 2 p − 1 and G p(λ) is a firstkind Chebyshev polynomial that deviates least from zero on the inter val [λ 0; λ max ] with the normalization G p (− 1 τ) = 1. Set z1 = cos (0.5π p). Then the left endpoint of the interval is determined by the parameters p and λ max using the formulas
(
λ 0 = λ max (z1 − 1) (z1 + 1) ∈ [− 1 τ ; 0],
)
τ = 16 p π − 1 λ max . 2
2
(18)
On the interval [0; λ max ], we have | G p(λ)| ≤ 1 and | Sν(λ)| ≤ 1. The polynomial G p is expressed in terms of the optimal polynomial T p by analogy with formulas (7), but, according to (18), the parameter τ is regarded as a function of p. The spectrum of smoothing operator (17) is ρLIM(λ) = (1 − G p2(λ))(1 + τλ)−1.
(19)
It holds that 0 ≤ ρLIM(λ) ≤ (1 + τλ)−1. On the right end of the spectrum, the highfrequency components are damped rapidly: ρLIM(λ max ) ≤ π 2 (16 p 2 ). This property of LIM ensures noticeable advantages in problems with anisotropic discontinuous coefficients, so the LIM smoother is basically used in the exam ples considered below. For the LIM smoother, the transition from the nonsmooth approximation v h on the current grid to the smoother approximation v new is written by analogy with (10). Adaptation of smoothers to the spectrum of difference operators. Adaptation is performed by improving the lower boundary λ*min of the highfrequency spectral segment [λ*min, λ max ] for each grid level. Assume that, in the course of multigrid iterations, after p smoothing steps with some spectral boundaries λ*min and λ max corresponding to the current grid level, the relation δ = || rp ||/|| r0|| is obtained for the initial and termi nal residuals. If δ < 1, for the Chebyshev smoother, inverting formulas of form (15), we find the new interface 2
⎛ ρ − 1⎞ λ*min = ⎜ 1 ⎟ λ max , ⎝ ρ1 + 1⎠
(
)
ρ1 = δ −1 + δ −2 − 1
1/ p
.
(20)
For the LIM smoother, using (18) and (19), the value of λ*min is updated using the formula
λ*min = π 2 (δ −1 − 1)λ max . 16p 2
(21)
If δ ≥ 1 , it is recommended using λ*min = 0.1λ max . The boundaries of the highfrequency spectrum are thus updated in the case of deviations from the prescribed smoothing quality ε smooth and the multigrid iteration rate ρ = ε 2smooth expected on the basis of the analysis of doublegrid representation (13). Adaptation can increase and reduce the number of smoothing steps at every grid level. The boundary λ min of the operator AH on the coarsest grid can be improved in a similar manner. In the computations, the adaptation of the smoothers improves the quality of the MM scheme, but the adaptation algorithm has many tuning parameters and further experimental verification is required, which goes beyond the scope of this paper. 6. NUMERICAL EXAMPLES Consider certain problems. Recall that the Courant number is defined as cou = 0.5τλ max . For diffusion equation (1) with operator (2) in a unit cube with k ≡ 1 and a0 = 0 , in the case of a uniform grid with N 3 nodes, this Courant number is equal to 6 τ N 2 . COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
ON THE SOLUTION OF EVOLUTION EQUATIONS
1283
1.0 1 2 3
0.5
0 0.45
0.50
Fig. 2. Profiles of the approximate solutions produced by (1) LIM with τ = 0.125 (solid), (2) MM with τ = 0.25 (dash dotted), and (3) MM with τ = 0.0625 (dashed).
The tuning parameters of the MM scheme include the maximum number l of grid levels; the relative accuracies ε MG and ε C of multigrid iterations and the solution on the coarsest grid, respectively, which characterize the reduction in the residual norm; the initial interface in the spectrum; the adaptation regime; the smoothing factor ε smooth characterizing the reduction in the residual norm at the smoothing stage; and the type of the interpolation operator (trilinear or problemdependent interpolation). Unless otherwise stated, we mean the following standard options: l = 5 if the number of steps in each direction is N d < 1024, and l = 6 otherwise, ε MG = 10 −5 , ε C = 10 −5, ε smooth = 0.5 , the LIM smoother with adaptation at the initial spectral interface λ*min = λ max 6 for all grid levels, and problemdependent intergrid transfer operators. In the case of adaptation, we analyze the efficiency of using the improved spectral interface at all subsequent time steps. Problem 1. Consider problem (1) with operator (2) in a unit cube with Dirichlet boundary conditions in the case of a0 ≡ 0 and f ≡ 0 . The diffusion coefficient and the initial function are defined as ⎧k , κ(t, x, y, z) = ⎨ 1 ⎩k2,
x < 0.5, x ≥ 0.5,
⎧0, u(0, x, y, z) = ⎨ ⎩1,
x < 0.5, x ≥ 0.5.
The initial function is used to determine boundary values, which are extended to the entire interval 0 = t 0 ≤ t ≤ T = 1. Using this problem as an example, we discuss the features of the schemes. For this pur pose, we set a small jump in the diffusion coefficient, setting k1 = 1 and k2 = 10 −3. Consider the problem 3 on a grid with 128 nodes at τ = 0.0625, 0.125 , 0.25 . Even the least value τ = 0.0625 cannot be treated as small, since it is associated with the large Courant number cou = 0. 5τλ max ~ 6 000 . Figure 2 shows the profiles of the approximate solutions at the time t = T along the Ox axis at y = 0.75. For the LIM scheme, the plots are indistinguishable, so the plot for τ = 0.125 is presented. The LIM scheme with τ = 0.25 is executed 45 times faster than the explicit scheme. For the MM scheme, the plots are shown for τ = 0.0625, 0.25. As the accuracy of the multigrid method is increased, the plots remain unchanged, which means that the solution of the twolevel implicit scheme has been reached. Figure 2 shows that the accuracy of the LIM scheme is higher than that of the implicit scheme. The cause is that the decrement function of the LIM scheme on the initial spectral segment nearly coincides with the exact function and, on the rest of the spectrum, is close to the decrement function of the implicit scheme. A single time step of the multigrid method with accuracy ε MG = 10 −5 requires 6–11 multigrid iter ations, and the total number of steps (each being approximately equivalent to a single step of the explicit scheme) is 12976 and 2720 for τ = 0.0625 and 0.25, respectively. In this computation, we used the LIM smoother and adaptation; moreover, the spectral interface λ*min found at the first time step according COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
1284
ZHUKOV et al. (a)
1
5E05
1.0
(b)
2 0
0.5 −5E05 0 0.4
0.5
0.001
0.42
0.44
0.46
Fig. 3. Profiles of the approximate solutions produced by (а) LIM and MM with operator P and (b) MM with operator P (dashed) and P (solid).
to (21) was used at all the subsequent time steps. As a result, the number of multigrid iterations is reduced (from 11 at the first step to 5–6 at the subsequent time steps). The deterioration of the accuracy of the multigrid method to ε MG = 10 −3 at τ = 0.0625 significantly (by 25 times) accelerates the computations, but the accuracy becomes unacceptable. At τ = 0.25, the 3 LIM scheme requires half as much time as the MM scheme. However, on a grid with 256 nodes, this advantage disappears. Note that, for the MM scheme, the number of multigrid iterations does not grow with mesh refinement (since the method is optimal). Accordingly, the CPU time in the transition to the finer grid increases by 8 times. For the LIM scheme, according to Section 3, the number of iterations grow as τλ max , i.e., in this case, linear growth is observed with respect to the number N of grid nodes in a single direction (since λ max = O (N 2 )). Therefore, the MM scheme is preferable in the case of time marching to a steady state. Problem 2. Consider problem (1) in the right parallelepiped Ω = ((x, y, z): 0 ≤ x, y ≤ 1, 0 ≤ z ≤ 10) with Dirichlet boundary conditions (determined by the initial function) specified on the left and forward faces of Ω and with homogeneous Neumann conditions set on the other faces. Suppose that a0 = 0 is a constant, while the diffusion coefficient, the righthand side, and the initial data are discontinuous:
x < 0.5 or y < 0.5, x ≥ 0.5 and y ≥ 0.5, ⎧0, x < 0.5 or y < 0.5, u(0, x, y, z) ≡ f (t, x, y, z) = ⎨ ⎩1, x ≥ 0.5 and y ≥ 0.5. ⎧k , κ(t, x, y, z) = ⎨ 1 ⎩k2,
(22)
In (22), k1 = 10 −6 and k2 = 10 3. For a preliminary analysis, we use a small grid: 64 × 64 × 16. The time step is τ = 0.05. The computation conditions are the same as in Problem 1. The accuracy of the multigrid method is ε MG = 10 −4 . If the faces x = 1 and y = 1 of the parallelepiped are interpreted as planes of symmetry in view of the boundary condition set on them, then it can be assumed that a heated medium with high thermal conduc tivity is enclosed by a nonheatconducting cold medium (insulator). The insulator gets warmer, but slowly. This problem is difficult for the MM scheme. First, there is noticeable anisotropy: the estima x y tion of the maximum eigenvalues along the Ox, Oy , and Oz axes yields λ max = λ max ≈ 4 × 106 and z 3 λ max = 2.5 × 10 . In the case of simple smoothers of the Gauss–Seidel type(point or linear), such anisot ropy leads to stagnation of multigrid iterations. Second, the diffusion coefficient has a large jump discon tinuity at the interface of two media. The discrete problem is characterized by a large Courant number ( ≈10 6 ). However, the computation of the evolution of the solution with such a step is meaningful: in the zone of sharp variations in the solution, the local Courant number is close to 1. Accordingly, this problem is not very difficult for the LIM scheme. The profiles of the approximate solution near the material interface after two time steps are shown in Fig. 3. For the LIM and MM schemes (with problemdependent intergrid transfer operators), they coin COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
ON THE SOLUTION OF EVOLUTION EQUATIONS
1285
Table 1. Implicit scheme based on the multigrid method First time step Type of intergrid transfer operators P ˜ P
Second time step
number of MM iterations
number of smoothing steps
number of MM iterations
number of smoothing steps
11
18294
4
13288
11
12318
2
6660
cide (see Fig. 3a). In this problem, we see that the solution involves fronts. As was noted above, in problems with discontinuous coefficients, intergrid transfer operators based on trilinear interpolation frequently lead to negative phenomena (divergence, oscillations). In the given problem with such operators, there is no divergence due to the good adaptation of the smoother, but the solution exhibits spurious oscillations, as can be seen in Fig. 3b, which shows the profiles of the approximate solutions near the discontinuity for linear interpolation operators (denoted by P, dashed line) and for problemdependent operators (denoted by P, solid line). It can be seen that the problemdependent operators ensure good accuracy. Table 1 presents numerical results produced by the MM scheme with adaptation at the first two time steps. At the first step in the course of the first 5–6 iterations, the spectral interface is determined accord ing to (21). In the course of further iterations, the convergence rate reaches an asymptotic regime with ρ ≈ 0.25 in accordance with the given smoothing factor ε smooth = 0.5 . The found spectral interface λ*min is used at the next time step, so the number of MM iterations is reduced. Table 1 shows, for example, that a single time step in the MM scheme (starting from the second) takes 6660 smoothing steps. In this case, the explicit scheme would require almost 2 million steps. The LIM scheme takes 2843 steps on the given small grid (4021 steps on a twice finer grid) and is comparable in effi ciency with the MM scheme. Problem 3. In the conditions of Problem 1, consider discontinuous diffusion coefficient and initial data: ⎧k , κ =⎨ 1 ⎩k2,
x < 0.5 x ≥ 0.5
or y < 0.5, and y ≥ 0.5,
⎧0, u(0, x, y, z) = ⎨ ⎩1,
x < 0.5 or y < 0.5, x ≥ 0.5 and y ≥ 0.5,
where k1 = 1 and k2 = 100. To analyze the problem, we use a grid with 643 nodes and the integration inter val 0 = t 0 ≤ t ≤ T = 0.001. The computations are performed with τ = T/8 and τ = T. The value τ = T = 0.001 cannot be treated as small, since, on the given grid, this step is associated with a sufficiently high Courant number (~1600 ). The computation conditions are the same as in Problem 1. The accuracy of the multigrid method is ε MG = 10 −5 . Figure 4 shows the profiles of the approximate solutions. It can be seen that the accuracy depends on τ. For large τ = T , the LIM scheme ensures better accuracy than MM, which is caused by the better de scription of the evolution of the basic smooth modes. For this reason, the computations with varying h and δ = k2 / k1 ( k2 is increased at k1 = 1) are performed at τ = T/8 = 0.000125. A single time step is executed with ε MG = 10 −7 . Below are the characteristics of computational costs for δ = 100 depending on the num ber N d of nodes in a single direction on a grid with N d3 nodes. Grid, Nd
Number of MM iterations
Number of smoothing steps
64 128 256 512
6 10 10 11
444 777 880 966
Note that, in this problem, the spectrum is anisotropic: under double mesh refinement, the relative spectral interface decreases by 4 times. In the course of the first few multigrid iterations, this interface is determined and the convergence rate reaches its asymptotic value ρ = 0.25. Due to the dependence of the COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
1286
ZHUKOV et al. 1 LIМ τ = T/8 2 LIМ τ = T 3 MM τ = T/8 4 MM τ = T
2.0
1.5
0
0.25
0.50
0.75
1.00
Fig. 4. Profiles of the approximate solutions at t = T produced by (1), (2) LIM and (3), (4) MM with τ = T/8 and τ = T.
spectral interface on h, the computational costs also depend (weakly) on h. In the absence of jump discon tinuities, there is no such dependence. Below are the characteristics of computational costs on a grid with 1283 nodes as functions of the jump δ in the diffusion coefficient. δ = k2/k1
Number of MM iterations
Number of smoothing steps
1 100 103 104 106
11 10 9 8 5
90 777 2150 6132 38530
As usual, the spectral interface is determined at the first few multigrid iterations; then the convergence rate reaches its asymptotic value ρ = 0.25. The adaptation of the boundary of the highfrequency domain to the spectral anisotropy of the differ ence operators requires certain computational costs. Eventually, the amount of computations grows, but is no faster than the predicted growth rate (in isotropic problems without adaptation, the growth of com putational costs is directly proportional to δ ). Problem 4. Consider problem (1) with operator (2) in the cube Ω = [− 0.5;0.5]3 with Neumann bound ary conditions on the time interval 0 ≤ t ≤ 0.01. Let κ = 1 and the initial function be u(0, x, y, z ) = 1. The absorption coefficient a0(t, x, y, z) is specified as a discontinuous function of the cylindrical radius r =
x 2 + y 2 as follows: a0 = 0 for r > 0.25, a0 = 10 6 for 0.125 ≤ r ≤ 0.25 , and a0 = 1 for r < 0.125. This
function can be interpreted as zones with different absorbing properties. 3
This problem is solved on a grid with 64 nodes. Note that, in general, the grid is not matched with the discontinuity surface of a0 . Then the domain 0.125 ≤ r ≤ 0.25 is covered by only eight grid cells along the radius, so we consider two grid levels. The maximum eigenvalue of the difference operator D = −Δ h, which is spectrally equivalent to the implicit scheme operator Ah is estimated as λ 0max = 5 × 10 4 , while the pres ence of strong absorption determines the upper boundary of Ah as λ max = 1.05 × 10 6 . Using the indicated information, the initial relative spectral interface η can be specified by the “spectral tuning” formula η = λ 0max ( λ max 6) . Then λ*min = ηλ max . The value η = λ max 6 is referred to as the standard tuning, which is optimal for the strictly isotropic case. COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
ON THE SOLUTION OF EVOLUTION EQUATIONS
1287
0.2 1 LIМ(p) 2 LIМ(2p) 3 MM 4 Explicit
0.1
0 −0.2
0
0.2
Fig. 5. Profiles of the approximate solutions produced by LIМ(p), LIМ(2p), MM, and the explicit scheme.
As a “reference,” we use the approximate solution produced by the explicit scheme with step −7 τ = 2 × 10 (50 thousand steps, 7.5 hours of CPU time). The explicit scheme with maximum step −7 τ = 5 × 10 takes 20 thousand steps over 3 hours. The computation based on the LIM scheme with −6 τ = 2 × 10 yields the polynomial degree p = 2. The total CPU time is 1 hour. Pictorially, the results of these three computations coincide. Let us present the numerical results obtained with τ = 2 × 10 − 3 , which corresponds to cou ≈ 10 . Figure 5 shows the profiles of the reference and approximate solutions at the final time produced by the LIM scheme with standard p specified by (6) (in this case, p = 36), the LIM scheme with doubled p, and the MM scheme. These profiles are denoted as LIМ(p), LIМ(2p), MM, and Explicit and are depicted by solid, dashed, short dashed, and dashdotted lines, respectively. The profiles are shown only for the central domain. Near the boundary of the domain, the differences in the plots are smaller. The MM scheme is applied with LIM smoother without adaptation but with spectral tuning, which ensures a good approximation to the exact spectral interface. The other parameters are two grid levels, the smoothing factor ε smooth = 0.5 , the relative accuracy ε MG = 10 −5, of multigrid iterations, and the relative accuracy ε C = 10 −5 of the solution on the coarsest grid. The time costs of LIМ(p), LIМ(2p), and MM are 18, 23, and 50 sec, respectively. Figure 5 shows that the graph of the LIM scheme with a standard param eter p is the closest to the reference solution. On the entire computation interval, the total costs of the LIM scheme are equivalent to 355 steps of the explicit scheme. The computation based on the MM scheme requires six multigrid iterations; moreover, smoothing on the basic finest grid takes 900 steps overall. Recall that every such step is equivalent in terms of costs to a single step of the explicit scheme. When the accuracy of the MM scheme decreases, for exam ple, to ε MG = 10 −3 , the smoothing procedure takes 480 steps and the graphs of the approximate solutions nearly coincide. For ε MG = 10 −1, a single multigrid iteration is executed, smoothing on the basic grid requires 150 steps, and the numerical accuracy is roughly 1% worse than in the preceding two cases. This computation takes 25% less time than that with the LIM scheme. Let us discuss the efficiency of smoothers in the MM scheme. The LIM smoother performs well in all regimes. Smoothing with adaptation leads to CPU time savings of 5 and 10% in the case of standard and spectral tuning, respectively, in specifying the initial spectral interface. The Chebyshev smoother with given parameters performs only in the regime with adaptation, but the computation time becomes twice as much as in the case of the LIM smoother, since the former leads to growing numbers of multigrid iter ations and smoothing steps at every multigrid iteration. The causes of this were discussed in [6–8]. Note that the solution of the given problem with the MM scheme involved intergrid transfer operators based on linear interpolation. 3
COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
1288
ZHUKOV et al.
Table 2. Comparison of the MM and LIM schemes on a 5123 CPU time, s Number of processors 1 8 64 512
MM
LIM
273 59.3 9.3 1.3
812 221 27.1 3.5
The numerical results show that the tuning parameters for the MM scheme should be chosen experi mentally. For the LIM scheme, there is nothing to choose, since there are no tuning parameters in this scheme. A common task for the two schemes is the estimation of the accuracy of time integration. It may happen that the step τ is determined by a complete problem that involves the parabolic equation as part. Then we have to perform additional computations with step τ 2 to estimate the accuracy and, if the required accuracy is not reached, τ continues to be reduced. As τ decreases, the competitiveness of the LIM scheme rises. The scalability (parallel efficiency) of the LIM and MM schemes has been repeatedly tested in com putations. For illustrative purposes, we give an example. Problem 5. Consider Problem 1 with diffusion coefficient k ≡ 1 and discontinuous initial function u(0, x, y, z) = sign x . The boundary data are determined by the exact solution u(t, x, y, z) = erf (0.5x t ). The computation conditions are the same as in Problem 1. A single time step is taken with τ = 0.05. The com 3 putations were performed on the K100 supercomputer [20]. The scalability was tested on a grid with 512 nodes on 1, 8, 64, and 512 processors of K100. The computation time decreased proportionally to the number of processors (see Table 2). For the given τ, the MM scheme was found to have a better absolute computation time. 7. CONCLUSIONS The multigrid and explicit iterative methods for solving threedimensional evolution problems were studied and further developed on the basis of the optimal properties of Chebyshev polynomials. An implicit scheme for a parabolic equation was studied. The scheme is solved by applying the multigrid method with intergrid transfer operators in problemdependent form and with the adaptation of the lower boundary of the highfrequency domain to the spectrum of the difference operators. The adaptation is performed in the course of multigrid iterations and improves the efficiency of the method, which is espe cially pronounced in the case of evolution problems with anisotropic discontinuous coefficients. The explicit iterative scheme with Chebyshev parameters was examined as a competitor of the multi grid scheme. An analysis of the schemes suggests that the MM scheme is preferable at large values of the parabolic Courant number, while the LIM scheme becomes competitive at relatively low Courant num bers. For evolution problems, the latter scheme can outperform MM even at large Courant numbers if the local Courant number in highgradient solution regions is close to 1. The computations show that the MM and LIM schemes ensure high performance on evolution prob lems with anisotropic and/or discontinuous coefficients. The experience of using these schemes suggests that they are well scalable and make it possible to overcome the difficulties associated with the develop ment of ultra parallel computations and the achievement of exaFLOP speeds in the future. ACKNOWLEDGMENTS This work was supported by the Russian Foundation for Basic Research, project no. 14210025. REFERENCES 1. R. P. Fedorenko, “A relaxation method for solving elliptic difference equations,” Comput. Math. Math. Phys. 1 (4), 1092–1096 (1962). 2. R. P. Fedorenko, “Iterative methods for elliptic difference equations,” Russ. Math. Surv. 28 (2), 129–195 (1973). COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015
ON THE SOLUTION OF EVOLUTION EQUATIONS
1289
3. R. P. Fedorenko, Introduction to Computational Physics (Mosk. Fiz.Tekh. Inst., Moscow, 1994) [in Russian]. 4. V. T. Zhukov, N. D. Novikova, and O. B. Feodoritova, Preprint No. 28, IPM RAN (Keldysh Inst. of Applied Mathematics, Russian Academy of Sciences, Moscow, 2014); http://library.keldysh.ru/preprint.asp?id=201428. 5. V. T. Zhukov, N. D. Novikova, and O. B. Feodoritova, Preprint No. 85, IPM RAN (Keldysh Inst. of Applied Mathematics, Russian Academy of Sciences, Moscow, 2014); http://library.keldysh.ru/preprint.asp?id=201485. 6. V. T. Zhukov, N. D. Novikova, and O. B. Feodoritova, “Parallel multigrid method for solving elliptic equations,” Mat. Model. 26 (1), 55–68 (2014). 7. V. T. Zhukov, N. D. Novikova, and O. B. Feodoritova, “Multigrid method for anisotropic diffusion problems based on adaptive Chebyshev smoothers,” Mat. Model. 26 (9), 126–140 (2014). 8. V. T. Zhukov, N. D. Novikova, and O. B. Feodoritova, “Multigrid method for elliptic equations with anisotropic discontinuous coefficients,” Comput. Math. Math. Phys. 55 (7), (2015). 9. V. T. Zhukov, “On explicit methods for the time integration of parabolic equations,” Math. Model. Comput. Simul. 3 (3), 311–332 (2010). 10. V. O. Lokutsievskii and O. V. Lokutsievskii, “On the numerical solution of boundary value problems for para bolic equations,” Dokl. Akad. Nauk SSSR 291 (3), 540–544 (1986). 11. A. A. Samarskii and E. S. Nikolaev, Solution Methods for Grid Equations (Nauka, Moscow, 1978). 12. I. M. Gel’fand and O. V. Lokutsievskii, “On difference schemes for solving the heat equation,” in Introduction to the Theory of Difference Schemes, by S. K. Godunov and V. S. Ryaben’kii (Fizmatgiz, Moscow, 1962) [in Rus sian]. 13. K. I. Babenko, Fundamentals of Numerical Analysis (NITs RKhD, Moscow, 2002) [in Russian]. 14. L. A. Lyusternik, “On difference approximations of the Laplace operator,” Usp. Mat. Nauk 9 (2), 3–66 (1954). 15. A. S. Shvedov and V. T. Zhukov, “Explicit iterative difference schemes for parabolic equations,” Russ. J. Numer. Anal. Math. Model. 13 (2), 133–148 (1998). 16. P. L. Chebyshev, Issues of least quantities related to approximate representations of functions (St. Petersburg, 1899), Vol. 1 [in Russian]. 17. V. I. Lebedev and S. A. Finogenov, “The order of choice of iteration parameters in the cyclic Chebyshev iteration method,” USSR Comput. Math. Math. Phys. 11 (2) 155–170 (1971). 18. F. R. Gantmacher, The Theory of Matrices (Chelsea, New York, 1959; Nauka, Moscow, 1966). 19. U. Trottenberg, C. W. Oosterlee, and A. Schuller, Multigrid (Academic, New York, 2001). 20. http://www.kiam.ru/MVS/resourses/k100.html.
Translated by I. Ruzanova
COMPUTATIONAL MATHEMATICS AND MATHEMATICAL PHYSICS
Vol. 55
No. 8
2015