- Email: [email protected]
Recent enhancements in MUDPACK , a multigrid software package for elliptic partial differential equations
Recent Enhancements in MUDPACK, a Multigrid Software Package for Elliptic Partial Differential Equations John
C. Adams
Nutional Center for Atmospheric P.O. Box 3000 Boulder, Colorado 80307
Research*
ABSTRACT
Enhancements discussed
1.
and
in a new illustrated
with
version
of the
multigrid
software
package
\IUIII’M:L
are
examples.
INTRODUCTION
The multigrid software package MUDPACK has undergone substantial revision since an earlier version dcscrihed in detail in this journal (see [I]). The changes include the availability of multigrid order solvers, the removal of some constraints
options, the addition of fourthon grid size, and the addition
of hybrid multigrid-direct-method solvers. These constitute a significant enhancement of the original package. The new version replaces and is incompatible with the earlier version. However, most of the details of the discretization, boundary conditions, relaxation methods, etc., contained in [l] are still accurate and are not repeated here. In the next section we briefly review the purpose and contents of MUDPACK. In the last two sections, we discuss each of the major changes and present examples which illustrate the advantages they provide. More atten-
*The National tion.
Center
for Atmospheric
Research
is sponsored
APPLZED MATHEMATZCS AND COMPCJTATZON 43:79-93 0 Elsevier Science Publishing Co., Inc., 1991 655 Avenue of the Americas, New York, NY 10010
by the National
Scienw
Founda-
(1991)
79
0096-3003/91/$03.50
JOHN C. ADAMS
80
tion is paid to supercomputer performance than in the examples given in [l, 21. Applications of the software to “real” problems taken from a variety of disciplines within the atmospheric sciences are given in [5].
2.
PACKAGE MUDPACK
tial equations
DESCRIPTION
software provides efficient solutions to elliptic partial differen(PDEs) using multigrid iteration [7, 8, 10, 121. The package is
user-friendly, is FORTRAN-77-portable, and vectorizes on Cray computers. The codes have been successfully run on a variety of computers and operating systems. Currently there are 27 FORTRAN solvers in the package for generating second- and fourth-order approximations to elliptic PDEs in a variety of forms. These include two- and three-dimensional, real and complex, separable and nonseparable, and self-adjoint equations. The solution regions are rectangular in the sense that the domain of each independent variable must be a bounded interval on the real line. This means that curvilinear coordinate systems such as spherical or cylindrical coordinates are acceptable. The codes are not restricted to Cartesian coordinates. Allowable boundary conditions include any combination of periodic, Dirichlet, or mixed derivative. The continuous problems are automatically discretized using second-order finite-difference formulae on uniform grids superimposed on the solution regions. This is transparent to a user, who only needs to supply the PDE, boundary conditions, and grid size. The discretization includes adjustment of second-order coefficients at coarser grid levels if there are nonzero first-order terms in the PDE which dominate and make it appear hyperbolic. If untreated, these can inhibit convergence. The resulting linear block-tridiagonal systems using multigrid iteration. Vectorized Gauss-Seidel
are approximately solved relaxation schemes based
on alternating points (red-black), lines, or planes (for three-dimensional anisotropic PDEs [ 161) are available [q]. Choice of the appropriate method of relaxation, which can usually be determined by PDE inspection, is crucial in achieving rapid convergence. Block Gaussian elimination can be used to solve coarse-grid equations with two-dimensional nonseparable hybrid solvers. Deferred corrections [13] are used to improve second-order estimates to fourth order. Error control for convergence monitoring is included. Full multigrid cycling (FMG) is available [7] and eliminates the need for initial guesses. Exact minimal work-space requirements are output. Equivalencing can be used to save storage. The package includes extensive documentation and test programs which can be used for testing when installing the codes on new systems.
Recent Developments 3.
SOFTWARE
in
81
MUDPACK
CHANGES
(a) Multigrid Options The new version of MUDPACK includes options for those wishing to implement variants of multigrid iteration. This provides users familiar with multigrid the opportunity to experiment with and fine-tune the underlying algorithms for greater efficiency. Default options are retained for those preferring black-box solvers. This is accomplished with a integer vector, MGOET,
cycles
which sets the kind of cycling (e.g., F [ll], V, W, or more general can be selected), the number of pre- and post-coarse-grid-correction
relaxation sweeps, the method of residual restriction, and the method of prolongation. A listing describing MGOFT, as it appears in all solver documentation files, is contained in the Appendix. The new default multigrid options set W(2,l) cycling, fully weighted residual restriction, and cubic prolongation. These are more robust than V(2,l) cycling with linear prolongation, which were all that was available in the previous version of MUDPACK. The new default options increase computational time by at least 25% pr multigrid cycle in three dimensions and 50% in two dimensions. Nevertheless, they are more reliable and can give better overall efficiency (see Example A). V cycling with linear prolongation can still be selected as a possibly more efficient choice for certain problems. (19) Fourth-Order Fourth-order
Solvers solvers in two and three dimensions
are now available
that
use deferred corrections [13] to improve second-order estimates to fourthorder accuracy. We briefly review this technique in one dimension. The extension to higher dimensions is straightforward. Suppose we wish operator notation):
to solve
the
following
linear
elliptic
equation
(in
lu= f Using standard second-order finite-difference grid (of size h) yields the linear system LU= which we can solve to discretization-level truncation error t=F-Lu
approximations
on a uniform
F, error using multigrid
iteration.
The
82
JOHN C. ADAMS
measures
how closely
tion equations.
the exact continuous
Simple Taylor-series
t =
h’(
cu,,,
solution
arguments +
satisfies
the discretiza-
show that t has the form
dux,,,) + O( h-‘))
where c, d are known coefficients from the elliptic equation. If U satisfies the discretization equations exactly, then it can he used to generate second-order approximations to uXI, and u,,,,. For example, if the uniform grid on [u, b] is a=x(l)
< ...
< ..*
and U(i) is the approximation at x(i), then the following difference formula can be used to approximate the third and fourth partial derivatives of u(x): at x=a:
u,,,
=
u,,,,
=
-5U(
1) + lSU(2)
-24U(3)
+ 14U(4)
-3U(5)
zh”
+ O(h”),
3U(1)-14U(2)+26U(3)-24U(4)+11U(5)-2U(6)
+ O(h”);
h’
at x=a+h:
u,,,
=
-3U(l)
+ lOU(2)
- 12U(3)
+6U(4)
2h3 2U(1)-9U(2)+16U(3)-14U(4)+6U(5)-U(6)
UXX.r, =
h’
at x = x(i),
I.4XXX=
u,,,,
=
where
-U(5)
+ O(h”),
+ 0( h”);
a + h < x < b - h:
-U(i-2)+2u(i-l)-2u(i+l)+u(i+2)
+ 0( h”),
2hR U(i-2)-4U(i-l)+SU(i)-4U(i+l)+U(i+Z)
h”
+ 0( h”).
Recent Developments Similar
difference
in
83
MUDPACK
formulae
are used at x = b - h and x = 17. These
can be
obtained by using the FORTRAN difference-formula generators described in [4]. The necessary difference equations are encoded in the MCDPACK fourthorder solvers. If we denote all of these by the difference operators
S”( II) =
u,,,. +
S“(V) =
u,,,,
O( h’),
+ 0( h’)
and let
T = h"[
cS”( V) + &P(V)] ,
then
T=t+O(h’). The fourth-order truncation-error estimate T is computed and passed down to coarser grids using weighted averaging. Then one full multigrid cycle with the default options is used to solve the correction
equation
LE=-T. Since
E is an O(h”)
approximation
to the exact error
e=u--U=L_1(-t), it follows that V=V+E yields a fourth-order approximation in V. Another related and effective method for generating higher-order approximations with multigrid is r-extrapolation (see [7]). The use of higher-order stencils is investigated in [14]. Examples B, D in the next section illustrate the significant improvements in accuracy that the fourth-order solvers can provide at reasonable computational and storage costs.
84
JOHN C. ADAMS
(c) Grid Size Second- and fourth-order approximations are generated on uniform Iby-m-by-n grids superimposed on boxes in three dimensions or l-by-m grids superimposed on rectangles in two dimensions. The new grid sizes (including boundaries) have the form 1=p2’+1
m=92j+l n=r2k+l
(x-direction
size),
( y-direction
size),
(z-direction
size),
i,j, k > 0 and p, 9, r > 1 are integers. In the earlier version of i = j = k was required. Since p, 9, r should be small for effective error reduction with multigrid iteration, the old constraint was not well suited for asymmetric grid sizes (see Example C). Let G denote the Z-by-m-by-n fine grid. Multigrid iteration is implemented on the ascending chain of grids where
MUDPACK,
G(0) < ...
< -..
where t = max(i, j, k) and each G(s) n(s) grid points given by l(S)
m(s)
=
for s = 0,. .,t
p 2nw(s+i--t,O)
=q2
n(s) =r2
mnx(s
+
1,
+j-t,O) +
1,
max(s+k-t,O)
+
=G,
has Z(s) by m(s)
by
1.
The coarsest grid, G(O), has p + 1 by 9 + 1 by r + 1 points. p, 9, r should all be chosen as small as possible and i, j, k, as large as possible within grid-size constraints for effective error reduction with multigrid iteration. p, 9, and r should each be 2 or a small odd value, since even values greater than 2 can be reduced by increasing i, j, or k. Large values for p (9, r) can reduce the convergence rate even if line relaxation in the x (y, Z) direction is chosen (see Example D). Larger values for p or 9 cause no problem if one of the two-dimensional “hybrid’ solvers (discussed next) is used. The earlier MUDPACK requirement that p (9, r.) must be greater than 2 when line relaxation in the x (y, z) direction is used and the x (y, z) boundary condition is periodic has been removed.
Recent Developments in
85
MUDPACK
(d) Hybrid Soloers The certainty of direct methods is combined with the efficiency of solvers in MUDPACK. multigrid iteration in “hybrid” multigrid-direct-method This has been done for two-dimensional nonseparable elliptic equations. Separable PDEs can be approximated with cyclic reduction [6, 151if a direct method is required. The hybrid solvers use block Gaussian elimination whenever the coarsest grid level is encountered within multigrid cycling. This eliminates the constraint that p, y must be small and provides a natural way to compare solutions from direct and iterative schemes. In the extreme case of I = p + I and rn = y + 1, the hybrid codes become (storage permitting) direct-method solvers. The use of Gaussian elimination requires approximately 2(p + l)(p + l)(y + 1) additional words of storage if periodic boundary conditions are not set in the y direction, and approximately 4(p + l)(p + l)(y + 1) additional words of storage if they are. If there are n = I = m grid points in each direction, then balance between multigrid iteration, which is an O(n”) algorithm, and the direct method, which O(p”) operations for solution on the coarsest grid, is roughly requires achieved when
This holds when
kc-
grid levels
are used
before
log, n
switching
2 to the direct
method.
Choosing
p
approximately equal to n /Zk will achieve rough parity between the direct and iterative parts of the hybrid algorithm. Larger values for p mean the direct method will dominate the computation, while smaller values will only marginally 4.
increase
the cost of multigrid
iteration.
EXAMPLES
Examples A, B, C, and D illustrate the advantages of the revisions discussed in Section 3. All runs were done on a Cray Y-Mpg/S64 computer. The name of the MUDPACK solver used is included in each table. Initial guesses are not provided. The results for the fourth-order solvers assume a second-order solution has been computed. The minimum required storage (in megawords), the exact least-squares error, and the execution time and
JOHN C. ADAMS
86
megaflop rate for one full multigrid cycle are recorded in Examples D. Discretization times and other overhead are not measured.
B, C and
We present full-multigrid (FMG) reEXAXIPLLCIA (Multigrid options). sults when solving a two- and three-dimensional elliptic PDE with different multigrid options. For both problems, we determine the accuracy in the approximation generated after only one full multigrid cycle by computing
where II.11is the Euclidean
norm,
h is the mesh size in each coordinate uh is the exact multigrid cycles),
discrete
solution
direction, (generated
by executing
u is the approximation generated after one full multigrid u0 is the exact PDE solution evaluated on grid h.
numerous
cycle,
Values for E,, less than 2.00 are considered to be acceptable FMG results. We compare results using V(2,l) cycles with linear prolongation and W(2,l) cycles with cubic prolongation. For both problems, point relaxation and fully weighted residual restriction are used. First we consider the two-dimensional cross derivative term on the unit square:
+[(1+x”)(l+yL)11/2-$-
nonseparable
elliptic
xyu(x,y)
Assume u(x, y) is specified at the upper x and lower satisfies the mixed derivative boundary conditions
~+y$-y2u(o,Y)=P(Y)
at
PDE
with a
= r(x,y).
y boundaries
x=0
and
Recent Deoelopments
in
87
MUVPACK
TABLE
1
\11’1>2(:R(.ll'l.O-l~lhll:NSIOYAI.) E/t With V cycles
11 I
iK
1 GiY I I!,?. I %i
With W cycles
141.66
1.69
140.90
0.99
137.16
0.99
130.35
1.00
and
The exact
is used
solution
for testing,
and the results
We next consider
are tabulated
the three-dimensional
in Table
separable
1.
elliptic
PDE
a”u a% a’u ax”+dys+a-“-(“+“+“)U(“,E/,-_)=T(l,y,-) & on the unit cube
is used
grids
with
for testing.
with
Dirichlet
Results
2” + 1 points
boundary
are tabulated
in each
direction
TABLE 41uD3sP
h I
conditions.
in Table
The exact
2 on three-dimensional
for increasing
2
With W cycles
ii7
4.73
z1
7.15
1.01
r64
9.78
1.00
12.45
1.00
14.55
1.00
1 ii% 1 Ez
k.
(THREI~-~IMENSIONAL)
With V cycles
solution
1.04
JOHN C. ADAMS
88 TABLE 3 brun3sr (sEcOND-OKIER) Grid 17X17X17 33X33X33 65X65x65 129X129X129 257X257X257
Storage (Mwords)
Rate (Mflops)
Time (set>
Error
0.013 0.088 0.651 4.995 39.144
26 52 92 141 178
0.01 0.06 0.30 1.63 10.42
0.19E - 2 0.51F.-3 O.l3c-3 0.33E - 4 0.83E-5
For both problems, results are poor using V cycles with linear prolongation and excellent using W cycles with cubic prolongation. See [2] for several examples where FMG results are satisfactory using the less expensive V cycles.
EXAMPLE B (Fourth-order). We approximate the three-dimensional PDE in Example A with second- and fourth-order solvers (see Tables 3 and 4). The default multigrid options and point relaxation are used. Computation and storage of the fourth-order truncation-error estimate account for the reduced megaflop rate and additional expense with the fourth-order solver. If accuracy is measured as a function of computational expense, then the fourth-order scheme is very cost-effective. It achieves more accuracy on the 33 X33 x 33 grid than the second-order solver does on the 257X257 X257 grid. The fourth-order solver requires too much storage (approximately 56 Mwords) for solution on the highest-resolution grid. The error reductions, with each doubling of resolution, by a factor of 4 for second-order and a factor of 16 for fourth-order, verify that discretization-level error has been reached.
TABLE 4 ~cn34s~ (FOURTH-0~13313) Grid 17X17X17 33X33X33 65X65X65 129x129x129
Storage (Mwords)
Rate (Mflops)
Time (set)
0.018 0.124 0.926 7.142
19 35 72 98
0.03 0.16 0.62 3.76
Error 0.26~ 0.13E 0.78~ 0.47E
- 4 -5 - 7 - 8
Recent Developments EXAMPLE
C
(Grid
use more grid points than in the vertical horizontal
on the
than vertical
region
where z() boundary in x and derivative
BY
in MUDPACK
size).
Typically,
three-dimensional
resolution
in solving
solution
= 10 and xc, = y. = 2500, is used to set the right-hand side and conditions and to compute the error. Here we assume u is periodic y, specified at the lower z boundary, and satisfies the mixed condition
The PDE
coefficients
are given
a(x)
default
parameters are shown
models
the equation
0 < x, y < 10,000 km and 0 < z < 10 km. The exact
uz+u=h(x,y)
The
weather
in the horizontal direction, where the scales are larger, direction. To simulate this, we use eight times more
multigrid
options
at
z=lO.
by
=l+sin”(7rx/x,,),
and
are p = y = r = 3 and in Table 5.
point
relaxation
i = j = k +3
TABLE
are
selected.
for increasing
Grid
5
hlUlI3bP storage Grid size 49X49X7 97X97x
13
(Mwords)
Rate (Mflops)
Solution time (set)
size
k. The results
Error
0.057
37
0.058
0.46~ - 4
0.324
76
0.189
0.12E - 4
193X193x25
2.299
128
0.829
0.29~ - 5
395X395x49
17.271
193
4.633
0.77E - 6
JOHN C. ADAMS
90
EXAMPLE D (Hybrid solvers). Here we consider Helmholtz equation in spherical coordinates,
the two-dimensional
on a one-degree grid on the full surface of a sphere of radius 1 (4 and 0 are the longitude and colatitude). ~(6, 0) is specified at the poles 8 = 0,~ and is periodic in 4. For testing purposes, we use the exact solution
and the coefficient
functions o(4,e)
= A(4,e)
=s+sin’ecos”$.
Choosing p = 45, y = 45, i = 4, and j = 3 fits the required 361-by-181 grid exactly. Multigrid is implemented on the sequence of grid sizes 46X 46, 91 X 46, 181 X91, 361 X 181. One FMG cycle with the default multigrid options is executed. Both point relaxation and line relaxation in the 4 direction are used. Multigrid iteration with relaxation and Gaussian elimination on the 46-by-46 coarsest subgrid are compared. The results are shown in Table 6. There is a reduction in the megaflop rate with the hybrid method. Nevertheless, the results indicate that relaxation alone is ineffective in reducing error with the high-resolution coarse grid. In fact, if we use the most robust form of two-dimensional relaxation (lines in both the r$ and the 0 direction) at all grid levels, then 36 multigrid cycles and almost 4 seconds of computer time are required to reach discretization level error. In only one FMG is discretization-level error reached in less than i second, when the
TABLE 6 ~~1~2,~112, ~~~124
Method
Storage (Mwords)
Rate (Mflops)
Solution time (set)
Relax (point) Relax (line 4) Hybrid (point) Hybrid (line 4) Hybrid (fourth)
0.618 1.059 0.815 1.256 1.321
118 111 51 53 58
0.067 0.086 0.229 0.245 0.285
Error 0.73E 0.67~ 0.44E 0.43E 0.47~
-
2 2 5 5 8
Recent Developments in
91
MUDPACK
coarse-grid direct method is used with either point or line relaxation at the higher-resolution grids. For a small additional cost, the hybrid method improves the second-order estimate to a fourth-order approximation.
APPENDIX.
DOCUMENTATION
A listing of the multigrid documentation
FOR
MGOPT
options parameter
file for each second-order
~~:oPI’, as it appears
in the
solver, follows:
MUDPACK
C MGOPT c
AN
C
AMONG
C
A
C
WILL
C
ARE
C
GRID
C
INTEGER
VECTOR
VARIOUS
DEFAULT
SET
BE
OF
OF
5
AND
CORRECTION
ALLOWS
DEFINED
ALGORITHM
THE
USER
IF RGOPTcI)=o
PARAMETERS OTHERWISE
SELECTED.
INTERNALLY
WHICH
OPTIONS.
MULTIGRID
INTERNALLY
SET
LENGTH
MULTIGRID
AS
TO
IS
(CHOSEN FOR PARAMETERS
FOLLOWS:
(SEE
SELECT
INPUT
THEN
ROBUSTNESS) IN MGOPT
BASIC
COARSE
BELOW)
KCYCLE=MGOPTclI
C
=-1
C
=0
IF IF
DEFAULT
CKCYCLE
C
IS TO
F CYCLING
BE
MULTIGRID
WILL
BE
C
=I
IF V
C
=2
IF w CYCLING
CYCLING
OPTIONS
RESET IS
USED
TO
ARE
TO
BE
USED
2)
TO
BE
USED
(THE
LEAST
IS TO
BE
USED
(THE
DEFAULT)
EXPENSIVE
K
CYCLING
LARGER
THAN
PER
CYCLE)
C
C .GE.
3 IF
MORE
GENERAL
(WARNING--VALUES
C
C
THE
C
RESULT
EXECUTION
C
IPRER=MGOPTCZI
TIME
IN THE
C
THE
C
RESIDUAL
NUMBER IS
C
COARSER
GRID
PER
CYCLE
NON-FATAL
OF
IS 1
USED
INCREASE
”
AND
SWEEPS
CYCLING
(DEFAULT
AND
IERROR=-
.‘PRE-RELAXATION LEVEL
BE
2
CONSIDERABLY
ERROR
RESTRICTED
TO
OR
EXECUTED
IS
VALUE
IS
BEFORE
INVOKED 2
AT
WHENEVER
THE
THE NEXT
MGOPTcl)=O)
C C
IPOST=MGOPTc3)
C
THE
NUMBER
C
HAS
BEEN
C
CORRECTION
C
WHENEVER
C C
OF
“POST
INVOKED HAS
AT BEEN
RELAXATION THE
NEXT
”
SWEEPS
COARSER
TRANSFERRED
EXECUTED
GRID
BACK
AFTER
LEVEL
(DEFAULT
AND VALUE
CYCLING
THE IS
RESIDUAL 1
MGOPTcl)=O).
IRESW=MGOPTc4)
=O
IF UNWEIGHTED
C
(WARNING--ORDINARILY
C
WHEN
USED
(IDENTITY)
RESIDUAL THIS
WITHIN
MULTIGRID
OPTION
RESTRICTI?“S GIVES
POOR
ARE
USED
?‘iSULTS
ITERATION)
C C C C
=l
IF DEFAULT
FULLY
WEIGHTED VALUE
WHENEVER
RESIDUAL
RESTRICTIONS
MGOPTcl)=O).
ARE
USED
(THIS
IS
THE
92
JOHN =2
C
IF
HALF
C
(THIS
C
AND,
C
SIMILAR
WEIGHTING
OPTION
IS
USED
REQUIRES
WITH
WITH
LESS
RED/BLACK
POINT
CONVERGENCE
RESIDUAL
TIME
THAN
RESTRICTIONS. FULL
RELAXATION,
RATES.
IT
WEIGHTING
SOMETIMES
SHOULD
C. ADAMS
BE
GIVES
USED
WITH
CAUTION)
C C
INTPOL=MGOPTC5)
C =i
C
MULTILINEAR
IF
C
TRANSFER
C
FROM
PROLONGATION
RESIDUAL COARSE
(INTERPOLATION)
CORRECTIONS
TO
FINE
AND
GRIDS
THE
WITHIN
Is PDE
FULL
USED
TO
APPROXIMATION
MULTIGRID
CYCLING.
C C
3
IF
MULTICUBIC
C
TRANSFER
C
FROM
C
(THIS
PROLONGATION RESIDUAL
COARSE
TO
IS
THE
(INTERPOLATION)
CORRECTIONS FINE
AND
GRIDS
DEFAULT
WITHIN
VALUE
IS
THE
PDE
FULL
WHENEVER
USED
TO
APPROXIMATION
MULTIGRID
CYCLING.
MGOPTCl)=O).
C C
THE
DEFAULT
C
FOR
ROBUSTNESS.
C WILL
VALUES
GIVE
c
ONE
c
C1,2,1,1,1)
IN
GOOD
AN
(2,2,1,1,3) IN
SOME
RESULTS
EARLIER
THE
IN CASES
WITH
VC2,l) LESS
VERSION
OF
VECTOR
MGOPT
CYCLES
WITH
COMPUTATION.
MUDPACK)
WERE
THIS
CAN
BE
CHOSEN
LINEAR
PROLONGATION
CHOICE
SET
WITH
(THE
THE
ONLY
VECTOR
REFERENCES 1
J. Adams, linear
Multigrid
MUDPACK:
elliptic
partial
FORTRAN software
differential
equations,
for the efficient
Appl.
Math.
solution
of
34:113-146
Comput.
(1989). 2
J.
Adams,
FMG
Proceedings
results
of the Fourth
with
the
Copper
multigrid Mountain
software
package
Conference
in
~~UDFAC:K,
SIAM,
on Multigrid,
1989, pp 1-12. 3 4 5
J. Adams,
2.0, NCAR,
J. Adams,
Fortran
subprograms
for finite
difference
26:113-116
(1978).
J. Adams,
R. Garcia, of the
science,
B. Gross,
J, Hack,
multigrid
scientific
software
Partial
J. Comput.
and R. Sweet, differential
of elliptic
V. Pizzo, and E. Ridley,
package
W. Briggs,
partial
Efficient
Fortran
equations,
in
subprograms
Elliptic
D.
A Multigrid
Jespersen,
appear.
Tutorial,
Vectorization
Vieweg,
Multigrid
SIAM,
of Compuer
Soloers
problems,
Math.
and S. Parter,
Philadelphia Programs
1987.
with Applications
to Computa-
1984.
methods
for partial
Anal. 24, Math. Assoc. Amer.,
J. Mandel
for the
Problem
(1977).
tional Fluid Dynamics, Numer.
in atmo-
MUDPACK
(M. Schultz, Ed.), Academic, New York, 1982, pp. 333-390. A. Brandt, Multi-level adaptive solutions to boundary value
W. Gentzsch,
Phys.
in preparation.
solution
8
Differential
Mar. 1990. formula,
D. Haidvogel,
J. Adams, P. Swarztrauber,
9
11
Elliptic
Version
Comp. 31:333-390
10
for Linear
User Document
spheric
7
Software
SCD
Applications 6
Multigrid
~IUDPACK:
Equations,
differential
equations,
in
Stud.
1984.
On the Multigrid
F-cycle,
Appl.
Math.
Comput.,
to
Recent 12
Developments
S. McCormick,
93
in MUDPACK Multigrid
Frontiers
Methods,
Appl. Math. 3, SIAM,
Philadelphia,
1987. 13
V. Pereyra,
Highly
value problems 14
S. Schaffer,
15
P.
16
Higher
Swarztrauber,
(G. Golub,
Ed.),
accurate
numerical
in n dimensions,
order multigrid Fast
Poisson
methods, solvers,
Math. Assoc. Amer.,
C. Thole
and
treatment
of elliptic
U. Trotenberg,
solution
of casilinear
Math. Comp. 24:771-783
Math. Camp. 43:89-115 in
Studies
boundary-
in
Numerical
(1984). Analysis
1985, pp. 319-369.
Basic
3D operators,”
elliptic
(1970).
smoothing
procedures
Appl. Math. Comput.
for the
19:333-345
multigrid (1986).
C. Adams
Nutional Center for Atmospheric P.O. Box 3000 Boulder, Colorado 80307
Research*
ABSTRACT
Enhancements discussed
1.
and
in a new illustrated
with
version
of the
multigrid
software
package
\IUIII’M:L
are
examples.
INTRODUCTION
The multigrid software package MUDPACK has undergone substantial revision since an earlier version dcscrihed in detail in this journal (see [I]). The changes include the availability of multigrid order solvers, the removal of some constraints
options, the addition of fourthon grid size, and the addition
of hybrid multigrid-direct-method solvers. These constitute a significant enhancement of the original package. The new version replaces and is incompatible with the earlier version. However, most of the details of the discretization, boundary conditions, relaxation methods, etc., contained in [l] are still accurate and are not repeated here. In the next section we briefly review the purpose and contents of MUDPACK. In the last two sections, we discuss each of the major changes and present examples which illustrate the advantages they provide. More atten-
*The National tion.
Center
for Atmospheric
Research
is sponsored
APPLZED MATHEMATZCS AND COMPCJTATZON 43:79-93 0 Elsevier Science Publishing Co., Inc., 1991 655 Avenue of the Americas, New York, NY 10010
by the National
Scienw
Founda-
(1991)
79
0096-3003/91/$03.50
JOHN C. ADAMS
80
tion is paid to supercomputer performance than in the examples given in [l, 21. Applications of the software to “real” problems taken from a variety of disciplines within the atmospheric sciences are given in [5].
2.
PACKAGE MUDPACK
tial equations
DESCRIPTION
software provides efficient solutions to elliptic partial differen(PDEs) using multigrid iteration [7, 8, 10, 121. The package is
user-friendly, is FORTRAN-77-portable, and vectorizes on Cray computers. The codes have been successfully run on a variety of computers and operating systems. Currently there are 27 FORTRAN solvers in the package for generating second- and fourth-order approximations to elliptic PDEs in a variety of forms. These include two- and three-dimensional, real and complex, separable and nonseparable, and self-adjoint equations. The solution regions are rectangular in the sense that the domain of each independent variable must be a bounded interval on the real line. This means that curvilinear coordinate systems such as spherical or cylindrical coordinates are acceptable. The codes are not restricted to Cartesian coordinates. Allowable boundary conditions include any combination of periodic, Dirichlet, or mixed derivative. The continuous problems are automatically discretized using second-order finite-difference formulae on uniform grids superimposed on the solution regions. This is transparent to a user, who only needs to supply the PDE, boundary conditions, and grid size. The discretization includes adjustment of second-order coefficients at coarser grid levels if there are nonzero first-order terms in the PDE which dominate and make it appear hyperbolic. If untreated, these can inhibit convergence. The resulting linear block-tridiagonal systems using multigrid iteration. Vectorized Gauss-Seidel
are approximately solved relaxation schemes based
on alternating points (red-black), lines, or planes (for three-dimensional anisotropic PDEs [ 161) are available [q]. Choice of the appropriate method of relaxation, which can usually be determined by PDE inspection, is crucial in achieving rapid convergence. Block Gaussian elimination can be used to solve coarse-grid equations with two-dimensional nonseparable hybrid solvers. Deferred corrections [13] are used to improve second-order estimates to fourth order. Error control for convergence monitoring is included. Full multigrid cycling (FMG) is available [7] and eliminates the need for initial guesses. Exact minimal work-space requirements are output. Equivalencing can be used to save storage. The package includes extensive documentation and test programs which can be used for testing when installing the codes on new systems.
Recent Developments 3.
SOFTWARE
in
81
MUDPACK
CHANGES
(a) Multigrid Options The new version of MUDPACK includes options for those wishing to implement variants of multigrid iteration. This provides users familiar with multigrid the opportunity to experiment with and fine-tune the underlying algorithms for greater efficiency. Default options are retained for those preferring black-box solvers. This is accomplished with a integer vector, MGOET,
cycles
which sets the kind of cycling (e.g., F [ll], V, W, or more general can be selected), the number of pre- and post-coarse-grid-correction
relaxation sweeps, the method of residual restriction, and the method of prolongation. A listing describing MGOFT, as it appears in all solver documentation files, is contained in the Appendix. The new default multigrid options set W(2,l) cycling, fully weighted residual restriction, and cubic prolongation. These are more robust than V(2,l) cycling with linear prolongation, which were all that was available in the previous version of MUDPACK. The new default options increase computational time by at least 25% pr multigrid cycle in three dimensions and 50% in two dimensions. Nevertheless, they are more reliable and can give better overall efficiency (see Example A). V cycling with linear prolongation can still be selected as a possibly more efficient choice for certain problems. (19) Fourth-Order Fourth-order
Solvers solvers in two and three dimensions
are now available
that
use deferred corrections [13] to improve second-order estimates to fourthorder accuracy. We briefly review this technique in one dimension. The extension to higher dimensions is straightforward. Suppose we wish operator notation):
to solve
the
following
linear
elliptic
equation
(in
lu= f Using standard second-order finite-difference grid (of size h) yields the linear system LU= which we can solve to discretization-level truncation error t=F-Lu
approximations
on a uniform
F, error using multigrid
iteration.
The
82
JOHN C. ADAMS
measures
how closely
tion equations.
the exact continuous
Simple Taylor-series
t =
h’(
cu,,,
solution
arguments +
satisfies
the discretiza-
show that t has the form
dux,,,) + O( h-‘))
where c, d are known coefficients from the elliptic equation. If U satisfies the discretization equations exactly, then it can he used to generate second-order approximations to uXI, and u,,,,. For example, if the uniform grid on [u, b] is a=x(l)
< ...
< ..*
and U(i) is the approximation at x(i), then the following difference formula can be used to approximate the third and fourth partial derivatives of u(x): at x=a:
u,,,
=
u,,,,
=
-5U(
1) + lSU(2)
-24U(3)
+ 14U(4)
-3U(5)
zh”
+ O(h”),
3U(1)-14U(2)+26U(3)-24U(4)+11U(5)-2U(6)
+ O(h”);
h’
at x=a+h:
u,,,
=
-3U(l)
+ lOU(2)
- 12U(3)
+6U(4)
2h3 2U(1)-9U(2)+16U(3)-14U(4)+6U(5)-U(6)
UXX.r, =
h’
at x = x(i),
I.4XXX=
u,,,,
=
where
-U(5)
+ O(h”),
+ 0( h”);
a + h < x < b - h:
-U(i-2)+2u(i-l)-2u(i+l)+u(i+2)
+ 0( h”),
2hR U(i-2)-4U(i-l)+SU(i)-4U(i+l)+U(i+Z)
h”
+ 0( h”).
Recent Developments Similar
difference
in
83
MUDPACK
formulae
are used at x = b - h and x = 17. These
can be
obtained by using the FORTRAN difference-formula generators described in [4]. The necessary difference equations are encoded in the MCDPACK fourthorder solvers. If we denote all of these by the difference operators
S”( II) =
u,,,. +
S“(V) =
u,,,,
O( h’),
+ 0( h’)
and let
T = h"[
cS”( V) + &P(V)] ,
then
T=t+O(h’). The fourth-order truncation-error estimate T is computed and passed down to coarser grids using weighted averaging. Then one full multigrid cycle with the default options is used to solve the correction
equation
LE=-T. Since
E is an O(h”)
approximation
to the exact error
e=u--U=L_1(-t), it follows that V=V+E yields a fourth-order approximation in V. Another related and effective method for generating higher-order approximations with multigrid is r-extrapolation (see [7]). The use of higher-order stencils is investigated in [14]. Examples B, D in the next section illustrate the significant improvements in accuracy that the fourth-order solvers can provide at reasonable computational and storage costs.
84
JOHN C. ADAMS
(c) Grid Size Second- and fourth-order approximations are generated on uniform Iby-m-by-n grids superimposed on boxes in three dimensions or l-by-m grids superimposed on rectangles in two dimensions. The new grid sizes (including boundaries) have the form 1=p2’+1
m=92j+l n=r2k+l
(x-direction
size),
( y-direction
size),
(z-direction
size),
i,j, k > 0 and p, 9, r > 1 are integers. In the earlier version of i = j = k was required. Since p, 9, r should be small for effective error reduction with multigrid iteration, the old constraint was not well suited for asymmetric grid sizes (see Example C). Let G denote the Z-by-m-by-n fine grid. Multigrid iteration is implemented on the ascending chain of grids where
MUDPACK,
G(0) < ...
< -..
where t = max(i, j, k) and each G(s) n(s) grid points given by l(S)
m(s)
=
for s = 0,. .,t
p 2nw(s+i--t,O)
=q2
n(s) =r2
mnx(s
+
1,
+j-t,O) +
1,
max(s+k-t,O)
+
=G,
has Z(s) by m(s)
by
1.
The coarsest grid, G(O), has p + 1 by 9 + 1 by r + 1 points. p, 9, r should all be chosen as small as possible and i, j, k, as large as possible within grid-size constraints for effective error reduction with multigrid iteration. p, 9, and r should each be 2 or a small odd value, since even values greater than 2 can be reduced by increasing i, j, or k. Large values for p (9, r) can reduce the convergence rate even if line relaxation in the x (y, Z) direction is chosen (see Example D). Larger values for p or 9 cause no problem if one of the two-dimensional “hybrid’ solvers (discussed next) is used. The earlier MUDPACK requirement that p (9, r.) must be greater than 2 when line relaxation in the x (y, z) direction is used and the x (y, z) boundary condition is periodic has been removed.
Recent Developments in
85
MUDPACK
(d) Hybrid Soloers The certainty of direct methods is combined with the efficiency of solvers in MUDPACK. multigrid iteration in “hybrid” multigrid-direct-method This has been done for two-dimensional nonseparable elliptic equations. Separable PDEs can be approximated with cyclic reduction [6, 151if a direct method is required. The hybrid solvers use block Gaussian elimination whenever the coarsest grid level is encountered within multigrid cycling. This eliminates the constraint that p, y must be small and provides a natural way to compare solutions from direct and iterative schemes. In the extreme case of I = p + I and rn = y + 1, the hybrid codes become (storage permitting) direct-method solvers. The use of Gaussian elimination requires approximately 2(p + l)(p + l)(y + 1) additional words of storage if periodic boundary conditions are not set in the y direction, and approximately 4(p + l)(p + l)(y + 1) additional words of storage if they are. If there are n = I = m grid points in each direction, then balance between multigrid iteration, which is an O(n”) algorithm, and the direct method, which O(p”) operations for solution on the coarsest grid, is roughly requires achieved when
This holds when
kc-
grid levels
are used
before
log, n
switching
2 to the direct
method.
Choosing
p
approximately equal to n /Zk will achieve rough parity between the direct and iterative parts of the hybrid algorithm. Larger values for p mean the direct method will dominate the computation, while smaller values will only marginally 4.
increase
the cost of multigrid
iteration.
EXAMPLES
Examples A, B, C, and D illustrate the advantages of the revisions discussed in Section 3. All runs were done on a Cray Y-Mpg/S64 computer. The name of the MUDPACK solver used is included in each table. Initial guesses are not provided. The results for the fourth-order solvers assume a second-order solution has been computed. The minimum required storage (in megawords), the exact least-squares error, and the execution time and
JOHN C. ADAMS
86
megaflop rate for one full multigrid cycle are recorded in Examples D. Discretization times and other overhead are not measured.
B, C and
We present full-multigrid (FMG) reEXAXIPLLCIA (Multigrid options). sults when solving a two- and three-dimensional elliptic PDE with different multigrid options. For both problems, we determine the accuracy in the approximation generated after only one full multigrid cycle by computing
where II.11is the Euclidean
norm,
h is the mesh size in each coordinate uh is the exact multigrid cycles),
discrete
solution
direction, (generated
by executing
u is the approximation generated after one full multigrid u0 is the exact PDE solution evaluated on grid h.
numerous
cycle,
Values for E,, less than 2.00 are considered to be acceptable FMG results. We compare results using V(2,l) cycles with linear prolongation and W(2,l) cycles with cubic prolongation. For both problems, point relaxation and fully weighted residual restriction are used. First we consider the two-dimensional cross derivative term on the unit square:
+[(1+x”)(l+yL)11/2-$-
nonseparable
elliptic
xyu(x,y)
Assume u(x, y) is specified at the upper x and lower satisfies the mixed derivative boundary conditions
~+y$-y2u(o,Y)=P(Y)
at
PDE
with a
= r(x,y).
y boundaries
x=0
and
Recent Deoelopments
in
87
MUVPACK
TABLE
1
\11’1>2(:R(.ll'l.O-l~lhll:NSIOYAI.) E/t With V cycles
11 I
iK
1 GiY I I!,?. I %i
With W cycles
141.66
1.69
140.90
0.99
137.16
0.99
130.35
1.00
and
The exact
is used
solution
for testing,
and the results
We next consider
are tabulated
the three-dimensional
in Table
separable
1.
elliptic
PDE
a”u a% a’u ax”+dys+a-“-(“+“+“)U(“,E/,-_)=T(l,y,-) & on the unit cube
is used
grids
with
for testing.
with
Dirichlet
Results
2” + 1 points
boundary
are tabulated
in each
direction
TABLE 41uD3sP
h I
conditions.
in Table
The exact
2 on three-dimensional
for increasing
2
With W cycles
ii7
4.73
z1
7.15
1.01
r64
9.78
1.00
12.45
1.00
14.55
1.00
1 ii% 1 Ez
k.
(THREI~-~IMENSIONAL)
With V cycles
solution
1.04
JOHN C. ADAMS
88 TABLE 3 brun3sr (sEcOND-OKIER) Grid 17X17X17 33X33X33 65X65x65 129X129X129 257X257X257
Storage (Mwords)
Rate (Mflops)
Time (set>
Error
0.013 0.088 0.651 4.995 39.144
26 52 92 141 178
0.01 0.06 0.30 1.63 10.42
0.19E - 2 0.51F.-3 O.l3c-3 0.33E - 4 0.83E-5
For both problems, results are poor using V cycles with linear prolongation and excellent using W cycles with cubic prolongation. See [2] for several examples where FMG results are satisfactory using the less expensive V cycles.
EXAMPLE B (Fourth-order). We approximate the three-dimensional PDE in Example A with second- and fourth-order solvers (see Tables 3 and 4). The default multigrid options and point relaxation are used. Computation and storage of the fourth-order truncation-error estimate account for the reduced megaflop rate and additional expense with the fourth-order solver. If accuracy is measured as a function of computational expense, then the fourth-order scheme is very cost-effective. It achieves more accuracy on the 33 X33 x 33 grid than the second-order solver does on the 257X257 X257 grid. The fourth-order solver requires too much storage (approximately 56 Mwords) for solution on the highest-resolution grid. The error reductions, with each doubling of resolution, by a factor of 4 for second-order and a factor of 16 for fourth-order, verify that discretization-level error has been reached.
TABLE 4 ~cn34s~ (FOURTH-0~13313) Grid 17X17X17 33X33X33 65X65X65 129x129x129
Storage (Mwords)
Rate (Mflops)
Time (set)
0.018 0.124 0.926 7.142
19 35 72 98
0.03 0.16 0.62 3.76
Error 0.26~ 0.13E 0.78~ 0.47E
- 4 -5 - 7 - 8
Recent Developments EXAMPLE
C
(Grid
use more grid points than in the vertical horizontal
on the
than vertical
region
where z() boundary in x and derivative
BY
in MUDPACK
size).
Typically,
three-dimensional
resolution
in solving
solution
= 10 and xc, = y. = 2500, is used to set the right-hand side and conditions and to compute the error. Here we assume u is periodic y, specified at the lower z boundary, and satisfies the mixed condition
The PDE
coefficients
are given
a(x)
default
parameters are shown
models
the equation
0 < x, y < 10,000 km and 0 < z < 10 km. The exact
uz+u=h(x,y)
The
weather
in the horizontal direction, where the scales are larger, direction. To simulate this, we use eight times more
multigrid
options
at
z=lO.
by
=l+sin”(7rx/x,,),
and
are p = y = r = 3 and in Table 5.
point
relaxation
i = j = k +3
TABLE
are
selected.
for increasing
Grid
5
hlUlI3bP storage Grid size 49X49X7 97X97x
13
(Mwords)
Rate (Mflops)
Solution time (set)
size
k. The results
Error
0.057
37
0.058
0.46~ - 4
0.324
76
0.189
0.12E - 4
193X193x25
2.299
128
0.829
0.29~ - 5
395X395x49
17.271
193
4.633
0.77E - 6
JOHN C. ADAMS
90
EXAMPLE D (Hybrid solvers). Here we consider Helmholtz equation in spherical coordinates,
the two-dimensional
on a one-degree grid on the full surface of a sphere of radius 1 (4 and 0 are the longitude and colatitude). ~(6, 0) is specified at the poles 8 = 0,~ and is periodic in 4. For testing purposes, we use the exact solution
and the coefficient
functions o(4,e)
= A(4,e)
=s+sin’ecos”$.
Choosing p = 45, y = 45, i = 4, and j = 3 fits the required 361-by-181 grid exactly. Multigrid is implemented on the sequence of grid sizes 46X 46, 91 X 46, 181 X91, 361 X 181. One FMG cycle with the default multigrid options is executed. Both point relaxation and line relaxation in the 4 direction are used. Multigrid iteration with relaxation and Gaussian elimination on the 46-by-46 coarsest subgrid are compared. The results are shown in Table 6. There is a reduction in the megaflop rate with the hybrid method. Nevertheless, the results indicate that relaxation alone is ineffective in reducing error with the high-resolution coarse grid. In fact, if we use the most robust form of two-dimensional relaxation (lines in both the r$ and the 0 direction) at all grid levels, then 36 multigrid cycles and almost 4 seconds of computer time are required to reach discretization level error. In only one FMG is discretization-level error reached in less than i second, when the
TABLE 6 ~~1~2,~112, ~~~124
Method
Storage (Mwords)
Rate (Mflops)
Solution time (set)
Relax (point) Relax (line 4) Hybrid (point) Hybrid (line 4) Hybrid (fourth)
0.618 1.059 0.815 1.256 1.321
118 111 51 53 58
0.067 0.086 0.229 0.245 0.285
Error 0.73E 0.67~ 0.44E 0.43E 0.47~
-
2 2 5 5 8
Recent Developments in
91
MUDPACK
coarse-grid direct method is used with either point or line relaxation at the higher-resolution grids. For a small additional cost, the hybrid method improves the second-order estimate to a fourth-order approximation.
APPENDIX.
DOCUMENTATION
A listing of the multigrid documentation
FOR
MGOPT
options parameter
file for each second-order
~~:oPI’, as it appears
in the
solver, follows:
MUDPACK
C MGOPT c
AN
C
AMONG
C
A
C
WILL
C
ARE
C
GRID
C
INTEGER
VECTOR
VARIOUS
DEFAULT
SET
BE
OF
OF
5
AND
CORRECTION
ALLOWS
DEFINED
ALGORITHM
THE
USER
IF RGOPTcI)=o
PARAMETERS OTHERWISE
SELECTED.
INTERNALLY
WHICH
OPTIONS.
MULTIGRID
INTERNALLY
SET
LENGTH
MULTIGRID
AS
TO
IS
(CHOSEN FOR PARAMETERS
FOLLOWS:
(SEE
SELECT
INPUT
THEN
ROBUSTNESS) IN MGOPT
BASIC
COARSE
BELOW)
KCYCLE=MGOPTclI
C
=-1
C
=0
IF IF
DEFAULT
CKCYCLE
C
IS TO
F CYCLING
BE
MULTIGRID
WILL
BE
C
=I
IF V
C
=2
IF w CYCLING
CYCLING
OPTIONS
RESET IS
USED
TO
ARE
TO
BE
USED
2)
TO
BE
USED
(THE
LEAST
IS TO
BE
USED
(THE
DEFAULT)
EXPENSIVE
K
CYCLING
LARGER
THAN
PER
CYCLE)
C
C .GE.
3 IF
MORE
GENERAL
(WARNING--VALUES
C
C
THE
C
RESULT
EXECUTION
C
IPRER=MGOPTCZI
TIME
IN THE
C
THE
C
RESIDUAL
NUMBER IS
C
COARSER
GRID
PER
CYCLE
NON-FATAL
OF
IS 1
USED
INCREASE
”
AND
SWEEPS
CYCLING
(DEFAULT
AND
IERROR=-
.‘PRE-RELAXATION LEVEL
BE
2
CONSIDERABLY
ERROR
RESTRICTED
TO
OR
EXECUTED
IS
VALUE
IS
BEFORE
INVOKED 2
AT
WHENEVER
THE
THE NEXT
MGOPTcl)=O)
C C
IPOST=MGOPTc3)
C
THE
NUMBER
C
HAS
BEEN
C
CORRECTION
C
WHENEVER
C C
OF
“POST
INVOKED HAS
AT BEEN
RELAXATION THE
NEXT
”
SWEEPS
COARSER
TRANSFERRED
EXECUTED
GRID
BACK
AFTER
LEVEL
(DEFAULT
AND VALUE
CYCLING
THE IS
RESIDUAL 1
MGOPTcl)=O).
IRESW=MGOPTc4)
=O
IF UNWEIGHTED
C
(WARNING--ORDINARILY
C
WHEN
USED
(IDENTITY)
RESIDUAL THIS
WITHIN
MULTIGRID
OPTION
RESTRICTI?“S GIVES
POOR
ARE
USED
?‘iSULTS
ITERATION)
C C C C
=l
IF DEFAULT
FULLY
WEIGHTED VALUE
WHENEVER
RESIDUAL
RESTRICTIONS
MGOPTcl)=O).
ARE
USED
(THIS
IS
THE
92
JOHN =2
C
IF
HALF
C
(THIS
C
AND,
C
SIMILAR
WEIGHTING
OPTION
IS
USED
REQUIRES
WITH
WITH
LESS
RED/BLACK
POINT
CONVERGENCE
RESIDUAL
TIME
THAN
RESTRICTIONS. FULL
RELAXATION,
RATES.
IT
WEIGHTING
SOMETIMES
SHOULD
C. ADAMS
BE
GIVES
USED
WITH
CAUTION)
C C
INTPOL=MGOPTC5)
C =i
C
MULTILINEAR
IF
C
TRANSFER
C
FROM
PROLONGATION
RESIDUAL COARSE
(INTERPOLATION)
CORRECTIONS
TO
FINE
AND
GRIDS
THE
WITHIN
Is PDE
FULL
USED
TO
APPROXIMATION
MULTIGRID
CYCLING.
C C
3
IF
MULTICUBIC
C
TRANSFER
C
FROM
C
(THIS
PROLONGATION RESIDUAL
COARSE
TO
IS
THE
(INTERPOLATION)
CORRECTIONS FINE
AND
GRIDS
DEFAULT
WITHIN
VALUE
IS
THE
PDE
FULL
WHENEVER
USED
TO
APPROXIMATION
MULTIGRID
CYCLING.
MGOPTCl)=O).
C C
THE
DEFAULT
C
FOR
ROBUSTNESS.
C WILL
VALUES
GIVE
c
ONE
c
C1,2,1,1,1)
IN
GOOD
AN
(2,2,1,1,3) IN
SOME
RESULTS
EARLIER
THE
IN CASES
WITH
VC2,l) LESS
VERSION
OF
VECTOR
MGOPT
CYCLES
WITH
COMPUTATION.
MUDPACK)
WERE
THIS
CAN
BE
CHOSEN
LINEAR
PROLONGATION
CHOICE
SET
WITH
(THE
THE
ONLY
VECTOR
REFERENCES 1
J. Adams, linear
Multigrid
MUDPACK:
elliptic
partial
FORTRAN software
differential
equations,
for the efficient
Appl.
Math.
solution
of
34:113-146
Comput.
(1989). 2
J.
Adams,
FMG
Proceedings
results
of the Fourth
with
the
Copper
multigrid Mountain
software
package
Conference
in
~~UDFAC:K,
SIAM,
on Multigrid,
1989, pp 1-12. 3 4 5
J. Adams,
2.0, NCAR,
J. Adams,
Fortran
subprograms
for finite
difference
26:113-116
(1978).
J. Adams,
R. Garcia, of the
science,
B. Gross,
J, Hack,
multigrid
scientific
software
Partial
J. Comput.
and R. Sweet, differential
of elliptic
V. Pizzo, and E. Ridley,
package
W. Briggs,
partial
Efficient
Fortran
equations,
in
subprograms
Elliptic
D.
A Multigrid
Jespersen,
appear.
Tutorial,
Vectorization
Vieweg,
Multigrid
SIAM,
of Compuer
Soloers
problems,
Math.
and S. Parter,
Philadelphia Programs
1987.
with Applications
to Computa-
1984.
methods
for partial
Anal. 24, Math. Assoc. Amer.,
J. Mandel
for the
Problem
(1977).
tional Fluid Dynamics, Numer.
in atmo-
MUDPACK
(M. Schultz, Ed.), Academic, New York, 1982, pp. 333-390. A. Brandt, Multi-level adaptive solutions to boundary value
W. Gentzsch,
Phys.
in preparation.
solution
8
Differential
Mar. 1990. formula,
D. Haidvogel,
J. Adams, P. Swarztrauber,
9
11
Elliptic
Version
Comp. 31:333-390
10
for Linear
User Document
spheric
7
Software
SCD
Applications 6
Multigrid
~IUDPACK:
Equations,
differential
equations,
in
Stud.
1984.
On the Multigrid
F-cycle,
Appl.
Math.
Comput.,
to
Recent 12
Developments
S. McCormick,
93
in MUDPACK Multigrid
Frontiers
Methods,
Appl. Math. 3, SIAM,
Philadelphia,
1987. 13
V. Pereyra,
Highly
value problems 14
S. Schaffer,
15
P.
16
Higher
Swarztrauber,
(G. Golub,
Ed.),
accurate
numerical
in n dimensions,
order multigrid Fast
Poisson
methods, solvers,
Math. Assoc. Amer.,
C. Thole
and
treatment
of elliptic
U. Trotenberg,
solution
of casilinear
Math. Comp. 24:771-783
Math. Camp. 43:89-115 in
Studies
boundary-
in
Numerical
(1984). Analysis
1985, pp. 319-369.
Basic
3D operators,”
elliptic
(1970).
smoothing
procedures
Appl. Math. Comput.
for the
19:333-345
multigrid (1986).