- Email: [email protected]
Numerical solution of continuity equations
Computer Physics Communications 12 (1976) 67—79 © North-Holland Publishing Company
NUMERICAL SOLUTION OF CONTINUITY EQUATIONS J.P. BORIS Plasma Physics Division, Naval Research Laboratory, Washington, D.C. 20375, USA
This paper compares and contrasts the six general methods available for solving the time dependent continuity equation numericaily. Since finite-difference techniques naturally arise as the most generally applicable solution methods, consideration is given to solving the major remaining difficulties encountered in the general application of lagrangian and eulerian finite differences. In the former, triangular mesh ceils with free reconnections are suggested as a complicated but serviceable solution to the lagrangian grid distortion problem. In the latter, nonlinear monotonic difference schemes are introduced as the natural way to avoid the eulerian positivity problem. Flux-corrected transport techniques, the most developed and efficient of the positive (monotone) finite-difference approaches, are described with examples.
1. The continuity equation
~fpd3r_fpv. This paper reviews the different numerical techniques which have been developed to solve the continuity equation ap/at =
—
V
(conservation form).
pV
(la)
In section 2 the six different available methods are compared and categorized. Within a framework provided by six requirements which ought to be satisfied by an ideal continuity equation algorithm, the strengths and weaknesses of the different approaches are analyzed. Finite-difference techniques have particularly attractive properties in this framework, thus section 3 is devoted to solving, at least partially, the remaining outstanding difficulties which finite-difference techniques suffer. Section 4 is devoted to flux-corrected transport (FCT) solutions of the nonlinear continuity equation which preserve the positivity property of the real equations and yet have relatively low numerical diffusion. In addition to eq. (1a), there are two other ways to write the continuity equation which get reflected in some of the numerical solution techniques: ap
dp V Vp ~j~convection
and
=
—
pV V,
compression (convection form)
(lb)
region
dA.
(lc)
region boundary (integral form)
The convective term V Vp displayed explicity in eq. .
.
.
(1b) gives eqs. (Ia)—(l c) their intrinsically hyperbolic form and causes really severe problems numerically. The compression term —pV V is sometimes absent, as in the Liouvifie equation, in which case it is often called the convection or the advection equation. The continuity equation appears in almost all descriptions of dynamic physical systems simply because it expresses two of the general principles in physics, conservation and causality. The integral form given by eq.(lc) states that the amount of material (f pd3r)in a given region (conservation) can only change due to flows of material (p V dA) through the boundary. In other words, any material currently at a given point could only have gotten there by leaving someplace else (causality). Continuity equations underlay compressible and incompressible fluid dynamics, hydrodynamics. plasma physics (Vlasov equation and MHD moment equations) and evenequation quantumdisplays mechanics. The continuity the positivity property that a quantity being transported will never .
turn negative anywhere in a reasonable flow field if it was everywhere positive to start with. This positivity
68
J.P. Boris/Numerical solution of continuity equations
property expresses in a continuum way the intrinsic corpuscular nature of matter. Thus matter cannot be removed from a region which is devoid of matter to begin with. This duality, in which matter obeys both particle and fluid-like equations on microscopic and macroscopic scales, respectively, has its ramifications for numerical solution techniques as well. The first three of the techniques described below take advantage of the underlying discrete basis of the continuity equat~ionwhile the last three aim more directly at solving the partial differential equation itself, The continuity equation is hyperbolic because it has first derivatives in both space and time and yet really is quite different from wave equations which are also hyperbolic. Therefore the solution techniques can also differ quite significantly. In the continuity equation the quantity being propagated appears in both the time and space derivative terms and could be described as “self”-hyperbolic. Generally the characteristics (flow) are in one direction only. Wave equations usually occur in pairs with two oppositely directed characteristics and the two quantities being propagated (say E and B in electromagnetics) appear in the time derivative of one equation and in the space derivative of the other. Thus wave equations might best be classed’ as “mutually” hyperbolic. There is no positivity to speak of and the most useful conservation property is quadratic rather than linear. Because of these distinctions, much of what is said here about continuity equations will not apply simply to wave equations.
2. Numerical solution techniques This survey considers only conservative solution techniques although nonconservative techniques are adequate for some problems. In most situations the conservation properties of the physics have to be mirrored in the numerical method or else nonphysical instabilities or unacceptable secular errors will result. Although positivity is probably as important as conservation, not all of the generally applied methods have recognized this fact to date and thus the survey cannot be limited to positive methods. It ought to be, however, since the most common and annoying problems which arise in the solution of coupled systems of continuity equations arise because the quantity convected becomes falsely negative somewhere or becaLse
numerical oscillations occur in the vicinity of steep gradients and shocks. The analysis is also limited to discussing the basic methods even though there are a host of hybrid techniques combining two or more of the basic approaches. The number of such combinations is staggering and even the definitions of the basic methods are in some dispute. Therefore no large useful purpose is served by trying to classify and analyze all of the possible combinations. Rather the general properties of each are described and the pros and cons considered. It is left intentionally to the reader to select the best method or hybrid of methods for his own particular application. As a framework for evaluating the algorithms, I present six general requirements which an ideal numerical algorithm for solving the continuity equation should satisfy to be generally useful. These requirements are given in what I estimate to be their order of importance. An ideal algorithm should: (1) be linearly stable for all cases of interest; (2) mirror conservation properties of the physics; (3) ensure the positivity property when appropriate; (4) be reasonably accurate; (5) be computationally efficient; and (6) be independent of specific properties of one particular application. This last requirement is a restatement of the facts that we are seeking generally useful methods and that all hybrid schemes cannot be discussed. As with the methods themselves, some of the requirements are related. Stability, conservation, and positivity generally relate to accuracy, for example, and yet no one of these can guarantee any of the others. On the other hand the requirements can also be partially contradictory. Accuracy and efficiency often tug in opposite directions, for example. Table I show the six methods I wish to discuss, the order being chosen to range from most discrete (or particulate) in nature to most continuous (or global). In their simplest form quasi-particle methods attempt to solve eq. (1) by simulating the microscopic physics of the particles which make up the fluid. If a large number N of quasi-particles are initialized in the calculation at positions {X 1(O)} at time t = 0, their later positions can be found at time t by integrating the simple orbit equation forward dX.(t)/dt = v.(t) (2)
J.P. Boris/Numerical solution of continuity equations Table 1 Possible solution techniques
______________
• Quasi-particle methods (1) Collisionless particles, stars and plasmas (2) Collisional particles for fluids • Characteristic methods • Lagrangian finite difference mcthods • Eulerian finite difference methods (1) Explicit versus implicit (2) Order versus accuracy •The finite element method • The spectral method
assuming that the velocity of each particle is known as a function of time. Eq. (2) is an ordinary differential equation and is usually solved easily with a high degree of accuracy. When each quasi-particle has a velocity which is determined solely as a function of the particle’s position at any time, i.e. ~) V1r~ V1X~‘ ~ ‘ —
‘
the quasi-particles cannot pass each other, except as a result of finite time-step errors [or perhaps singularities in V(X, t)]. The model is then appropriate for collisional fluids. When several quasi-particles can occupy essentially the same location in space and yet be moving with significantly different velocities, the model is more appropriately used to describe “collisionless” systems such as plasmas on the microscopic scale or stars in self-gravitating systems on a somewhat grander scale. Such “collisionless” systems of course undergo “collisions” when the impact parameters of two real particles (stars) becomes sufficiently small but such occurrences are rare. When the force law for the interacting “collisionless” particles is inverse square, it can be shown that the major component of the force is coherent and long range. Thus, most solutions of this type are seeking statistical and average properties of the flow and the individual particle orbits are relatively unimportant. The collisional quasi-particle methods which have been tried to date include PlC [1], GAP [2], and VORTEX [3] and more recent attempts of the same genre by workers at Los Alamos [4]. These methods introduce what amounts to a finite-difference grid for some of the derived physical quantities such as pressure while retaining the particle’s mass, momentum, and energy as particle variables. The collisionless quasi-
69
particle methods are reviewed by Birdsall and Langdon [5] and in the Proceedings of the Fourth Conference on the Numerical Simulation of Plasma [61.The .
original applications seem to stem from work by Hockney and Buneman [7] The term quasi-particle has been used because it is almost never possible to carry in the computer memory as many particles as are involved in the actual physical problem of interest. It is not unusual to simulate systems of lO~particles (galaxies) to 1019 particles (I cc of dense plasma) using 103_.106 quasi-particles. In fact, when the number of quasi-particles exceeds ~ or so, it is not usually possible to treat the forces by calculating a direct interaction law. A potential equation is interjected between the source terms (charge or mass) and the corresponding long range fields (electromagnetic or gravitational) and solved on a finite-difference grid in either configuration or transform space. The quasi-particle methods are generally stable, and conserveconservation though where they of should both(almost energy and by definition) momentumalis difficult in the collisionless models. These methods can also guarantee positivity since it is just as hard to remove a quasi-particle from an already empty region of space as it is a real particle. The methods can even be made reasonable general and flexible. Where these quasi-particle methods fall down is in efficiency and accuracy. Because the statistics and hence the smoothness of the computed solutions depend on the number of quasi-particles which the user can afford, many applications of interest are beyond the particle approach. Long-time, turbulent, collisional flows are unbearably expensive and collisionless “plasma” simulations usually break down for calculations which are run longer than a collision time for the quasi-particles. Some workers like to think of plasma simulation as an unavoidably collisional approximation to the collisionless Vlasov equation. The Vlasov equation,in reality, is a collisionless approximation to a real collisional plasma. Characteristic methods, which take advantage of the physical content of the specific system of eqautions being solved, are not entirely counter to our requirement (6) because different characteristic methods can often be formulated for different sets of governing equations. The equations of ideal MHD [8] and ideal compressible flow [9] for example, can be formulated and solved by characteristics and these solutions are often the best that can be obtained because the discon,
70
J.P. Boris/Numerical solution of continuity equations
tinuities which arise naturally in these physical systems are followed individually like particles in the previous methods and the continuous fluid behavior between the characteristics of the discontinuities is easy to represent. Nevertheless, characteristic methods are not generally applicable even though conservation and positivity are relatively easy to ensure. The presence of even one diffusion term in the governing equations usually abrogates the use of characteristic methods because the characteristics of the model usually cease to exist, A further drawback to characteristic methods is their complexity and relative inefficiency for multidimensional flows. The characteristic interfaces are complicated lines in two dimensions and surfaces in three dimensions and these can become progressively more and more difficult to follow as the flows proceed in time. This relative inefficiency, which is not troublesome in simple flows where a few characteristics may be sufficient to describe a whole flow, becomes far more serious when the new parallel and pipeline computers are contemplated. Since the characteristics, like the quasi-particles, must each be considered individually and since aspects of random access to computer memory as well as complicated logic are involved, the future of these methods on the newer computers is correspondingly dim. We come now to a discussion of finite differences. Unlike the quasi-particle and characteristic methods, finite-difference solutions of the continuity equation are based on an approximation to eq. (1) in which the continuous flow parameters p(x, t) and V(x, t) are discretized to a finite set of representative values ~p1},
classic text is Richtmyer and Morton [10] but the literature is virtually limitless. Although a detailed classification of possible finite difference schemes is also probably purposeless (there are two many hybrids) and may be impossible, a few of the possible choices should be discussed. The first of these is the distinction between lagrangian and eulerian methods. A lagrangian finite-difference scheme [11—13] is one in which the location of the grid points is allowed to move along with the fluid. This approach is useful because the convection term is transformed away in the moving system and hence the positivity property is easily ensured as well as conservation. The trouble with this approach is that complicated flows in two or three dimensions will rather quickly distort the moving mesh. Much of the accuracy is lost, often before useful information can be extracted from the calculation. Eulerian methods, on the other hand, keep the grid fixed in position and hence distortion of the mesh can not occur. Since fluid is now flowing across a stationary and usually rather coarse mesh, truncation errors occur which are secular in time and which can play havoc ~yithaccuracy and positivity in a long run [l0,14—16]. Nevertheless, these problems with both the eulerian and lagrangian approaches can be overcome as described in the next section. Here I wish to point out that the lagrangian finite-difference approaches, with the mesh distortion problem cleared up, are not completely general because aspects of eulerian flow appear in realistic coupled fluid systems which are being treated on a lagrangian grid.
defined at N distinct points in space called grid (or mesh) points. The desired continuous functions are assumed to be known only at the mesh points and a whole culture of techniques have been developed to predict future values of fp~}at some time t given ~p1}at t = 0. The term finite differences refers to the fact that the spatial derivatives (and usually the time derivative well)e.g. are approximated using only the grid point as values ~o• —i—-—, (4)
a /at =
—
ax X
__!__
1
X1+i
~i i
The attraction of these finite-difference methods is their richness and simplicity; no one need use the same finite-difference scheme twice in a row. The
Consider the ideal compressible fluid equations. V V (Sa) p apV/at = V- pVV— VP (Sb) —
—
ac/at =
—
V eV— V PV,
(Sc)
where 2 (Sd) e = [P/(7 1)] + ~-pV Ifwe define V~= [(c +P)/e] Vas a new velocity which varies from Vto about ‘yVand which is always parallel .
—
to
v, the energy continuity equation can be written as
ac/at =
—
and is in
V CV*
(5d’)
exactly the same form as eq. (5a) for the mass
J.P. Boris/Numerical solution of continuity equations
density. Since V V’~generally, a mesh moving at V is lagrangian for ~nassdensity, but not for energy density while a mesh moving at V’ is lagrangian for the energy density but not for the mass density. This means that full advantage can never be taken of a lagrangian grid; some aspects of a realistic calculation are always eulerian. Thus the eulerian grid problems with finite differences have to be controlled or solved for lagrangian, as well as eulerian, approaches to be fully viable, The second choice which I wish to consider is the time integration technique, in particular the neverending controversy of choosing an explicit versus an implicit scheme. Implicit schemes for the convective derivatives involve, in the worst cases, gigantic matrix equation solutions and in the best cases the inversion of a tridiagonal matrix. Thus, their operation count and hence speed is appreciably worse than obtained using corresponding explicit schemes. The advantage of implicit schemes professed to outweight the dual disadvantages of program complexity and operation count per timestep is the ability to take much longer finite-difference timesteps without exciting numerical modes of instability. The answer lies in the physics (mathemetics) of the specific problem being solved, When the physical phenomenon of interest varies appreciably on a timescale r, no calculation with & r can reasonably be claimed to reproduce the phenomenon accurately. In such cases computationally efficient explicit methods are called for with L~t
71
forward differenced beyond second-order centered differencing tend to damp the undesired high frequency sound waves strongly. Thus, sonic energy from the whistle gets deposited as heat near the whistle by a conservative, forward-differences implicit scheme. A centered scheme, or quite near centered, would allow most of the sonic energy to propagate away as sound without heating the local air. The strongly dispersed short wavelength sonic components would generally stay in the vicinity of the whistle far too long, however, and might eventually deposit their energy near the whistle anyway. The best solution for such a quandary is usually either to rescale the problem so the two timescales are not so disparate or to reformulate the problem entirely so the two timescales are separated entirely. In the whistle problem this would involve a multiple timescale analysis in which the high frequency oscillations are averaged out. Such a reformulation is not always available in which case neither the implicit nor the explicit methods can be judged fully satisfactory. Both finite-element and spectral methods utilize continuum concepts to the limit in solving eq. (1). The functions of interest, p(X, t) and V(X, t), are expanded in a set of useful basis functions whose expansion coefficients are now determined from a complicated set of coupled, usually nonlinear, but ordinary differential equations. In the case of finite elements [17] the basis functions are localized to a mesh cell or two and the techniques of matching expansion function boundary conditions at cell boundaries reminds one very much of using splines. The spectral methods generally employ a series of global basis functions [18] but otherwise the general characteristic is a projection to recover ordinary differential equations which are implicit in one way or another. The finite element methods acquire their implicitness by coupling of spatially adjacent basis functions via projection when evaluating the left-hand side a/at of eq. (1). The spectral methods become effectively implicit when the huge convolutions are needed to evaluate, for example, V(~V/OX) where the velocity field is represented as a Fourier expansion. The disadvantages of these two approaches are clear: (1) they are slow per timestep compared with corncomparable sized finite-difference methods for general problems; (2) they clearly are ill-suited (spectral methods more so than finite element) to problems displaying
72
J.P. Boris/Numerical solution of continuity equations
highly localized effects such as shocks; and (3) they are generally complicated and expensive to program, particularly so for nonlinear terms, The advantages are greater accuracy by far then simple difference or particle methods for smooth well-behaved flows and oscillations, The spectral methods have been described as “infinite-order” accurate and have been characterized as “optimal” for certain finite-difference evaluations [16]. Whether or not the potential increased accuracy outweight the drawbacks for real problems using either finite elements or spectral methods depends very strongly on the specific problem being solved. While the degree to which these two classes of techniques satisfy or do not satisfy our requirements (5) and (6) may be subject to legitimate debate, there is no doubt that the positivity property of the continuity equation and accuracy cannot be guaranteed when expanding in a time-independent set of basis functions, The Gibbs phenomenon [16] represents a minimal, or lower-bound error, which is only reduced by increased resolution. In a very real philosophical sense an uncertainty principle is at work. The numerical error can only be reduced by increasing the number of discrete degrees of freedom in the representation. Before giving my own evaluations, I would like to quote three conclusions, in a similar context, from another source [19] (1) The spectral method has no stability problems but is much more complicated and slower than generalized difference methods. (2) It is doubtful whether the finite-element method, based on piece-wise polynomials, can compete with the above methods. (3) If difference methods are used, they should be at least fourth-order accurate. While (I agree basically with these remarks, they do not really encompass the quasi-particle schemes nor do they adequately reflect a very important piece of personal experience. Whenever a theory is to be tested, or an algorithm, or a new mathematical technique, attention rather quickly turns to the crucial yet simple conservation equations of ideal compressible flow, eqs. (5a)— (Sd). Now I know that finite-difference methods have solved the transient Rankine—Hugoniot shock problem adequately, and I even know that the various ‘f~vors of quasi-particle methods can do the same vital problem
creditably if enough particles are used. I do not know of any calculation, even in one spatial dimension, using either a finite element or a spectral method which has correctly solved for an ideal gas compressible shock. Until such calculations become common place and computationally attractive, finite-difference methods would seem to have inside track, Therefore, the next two sections take a closer look at the remaining deficiencies in lagrangian and eulerian finite-difference techniques.
3. Improving finite-difference techniques There is undoubtedly some merit in trying to boost the performance of quasi-particle methods on the one hand and the basis-function expansion methods on the other toward the level obtained by finite differences. But it is a low risk—high return investment to patch up the obvious failings of the front runners finite differences. In the case of lagrangian methods, the major outstanding problems arise from secularly unattractive distortions of the grid which wreck calculations of interesting flows quite quickly. In the case of —
. . Fig. 1. A Schematic illustration o! the use of a two-dimensional grid of triangles to represent a iarge amplitude, separating, free-surface gravity wave.
J.P. Boris/Numerical solution of continuity equations
eulerian methods, the major outstanding weakness in a huge class of problems of real interest is the need for a large artificial damping (numerical diffusion) to fill in what would otherwise be pits of “negative density” in the calculated profiles. Since the “eulerian” positivity problem is encountered even in Lagrangian calculations for many situations, it demands the greater share of attention. The distortion problem for lagrangian grids can be handled adequately by using a finite-difference mesh of triangles rather than rectangles (tetrahedrons hi three dimensions). Fig. 1 shows the representational advantage of using triangles for a breaking gravity wave. This advantage, cribbed from engineers and original finite element buffs who have been using meshes of triangles for years, was adopted to finite differences ofincompressible flow [12,20]. Thenoticeableimprovements in performance arise for two reasons: first, the general connectivity premitted by triangles allows smooth representation of much more complicated shapes than can be treated smoothly by an equal number of rectangular cells, and secondly since the number of lines meeting at a vertex is not fixed, it grid distortions. can be varied during a calculation to prevent severe While this approach may not be the only method of dealing with severe grid distortions (see, for example, ref. [13]), it is the most physical and the least dissipative technique employed so far. Fritts [21] has recently employed lagrangian triangular grid techniques in studying nonlinear aspects of free-surface waves including strongly sheared flows. Ten to fifteen full cycles of the waves can be integrated reversibly without significant deterioration of the solution and the reconnection procedure in the presence of strong shear is reversible (i.e. nondissipative). Fig. 2 shows a plot of the wave natural period versus grid size for four different resolution calculations of the same physical problem. The second-order accuracy over long times is demonstrated by the parabolic form of the curve and extrapolation to zero stepsize gives improved accuracy as expected. An added advantage accrues from using triangular cells. Since it is difficult, if not impossible, to label the vertices unambiguously as either even or odd, the nonlinear mesh-separation instabilities, which plague low-dissipation rectangular cell techniques, seem to be absent or damped with triangular cells.
~
73
Bt ~.OO4
,
~
,
.
8t~.OO4~
,
/ ,
/
-
8m~.OO2
,
.
.,
~
_____
I
~i*
~‘~‘~‘;;~‘ ~.i
‘I
ii
ft .
tt~Y*ltttt~i ~
~l.i,% ~
— ~—.- ‘
tti
i
~__L__.~
~ ~—~---.~
-C
-
Hg, 2a. 138
-
136
r
~
-412732 + 00020
134
132
-
130
-
128 126
-
I
I
.1
2
.3
Fig. 2b. Fig. 2. Wave period versus normalized numerical grid size for four different resolution calculations of a finite depth freesurface wave.
The remaining drawbacks to the triangular grid approach are the complexity of the programming required and the fact that the operations which are rerequired seem to hinge strongly on linked lists, random access, and sequential processing. With several pipeline and vector processors working and several more soon to be in operation, efforts toward the “vectorization” for parallel processing of these triangular-grid techni-
74
J.P. Boris/Numerical solution of continuity equations
ques would seem to be well worthwhile. The major remaining problem with eulerian techniques for solving continuity equations is the competition between accuracy and positivity [19,22] wherein fully second order techniques cannot maintain positivity Consider the rather general three-point approximation to eq. (la)
SIDED
~
WE +Vj+i/2(p~l—PiQ)—i)fl/2(p~—~1),
where ei+i/2
~+i/2/~~+i/2’
1~11~>~/+1/2.
2
cycle
00
I
=
60
-
I
~—
—
X (cell no.)
~—~
—
X (cell no.)
---~
21
•~--
L x Ice); ~)
O~cell
~)
r;
2 -~
X (cell no)
00
—-~
•-~
i 0
._
____________
2
FROG
-L
H
I
--
~
-___________________
LEAP-
L_
— ~
I....:
-
P -
s—
2t
cycle 800
2
-
~
(6)
~.
100 0
----~
-i
X (cell no.)
100 -
(7a)
The upper limit arises from the explicit diffusion timestep condition [10] while the lower limit is the Lax— Wendroff damping [23]. Unfortunately positivity is only ensured linearly when 1 61+i/2 (7b) >i the first-order upstream-centered scheme result, We appear to be caught between a rock and a hard place here but the escape route is signaled in the preceding sentence by the word “linearly”. By relaxing the linearity implied by eq. (6) and letting the diffusion coefficients be nonlinear functionals of the flow velocities {c/+i/2}, we can hope to reduce the integrated dissipation below the rather ghastly limit (7b) and yet retain sufficient dissipation near steep gradients to ensure positivity. A literature is beginning to form about these “monotonic difference schemes” [22,24—26] since the dilemma of accuracy versus positivity in eulerian difference schemes can be resolved in no other way. Several different approaches are possible and the area is still largely unexplored. Boris and Book introduced this basic nonlinear approach with the techniques of flux-corrected transport [16, 23—25]. The rewards were quickly seen. Fig. 3 shows a comparison of four common difference schemes solving the standard square wave problem [14,24]. The effects of excess numerical damping in the donor-cell treatment (upstream centered first order), and of excess dispersion in the leapfrog ,
ONE-
Eq. (6) is in finite-
difference conservative form with whole indices re presenting cell centers and half indices indicating cell interfaces. The additional numerical diffusion terms with diffusion coefficients {~/+i/2} have to be added to ensure positivity The stability of eq (6) is ensured at least roughly, when ~ >v
f r 20
SHASTA 2~ ~
2~
(FCT)
.1 -
X (cell no.) —~
E 100 0
x Ice II
no.)
.
..
100
Fig. 3. Comparison of four difference schemes solving the square wave problem. The merits of flux-corrected transport relative to the other methods are clear.
and Lax—Wendroff treatments are clearly visible. Dispersion manifests itself as a trail or projection of oscillations in the computed solution near discontinuities and sharp gradients of the “correct” solution. VanLeer [27] through geometric arguments, and Harten [26] by a more analytic approach, have developed other monotonic finite-difference schemes for avoiding the accuracy— positivity dilema, Clearly the study of monotonic schemes will be fruitful because of the central role that eulerian finite-difference techniques naturally play. ,
4. Flux-corrected transport algorithms The first, and so far the most developed and used, of the monotonic schemes is flux-corrected transport (FCT). The calculation in fig. 3 was performed by the first FCT Algorithm SHASTA [24] and had an error about four or five times smaller than the simple linear methods also shown. The damping was second order as
J.P. Boris/Numerical solution of continuity equations
were the relative phase errors [161. This algorithm was used widely in the two-dimensional ideal compressible fluid dynamics code SHAS2D [28]. The basic techniques was quickly generalized to cylindrical and spherical systems, to lagrangian as well as fixed eulerian grids, and was applied to a number of one-, two-, and threedimensional problems [29,30]. More recent work has been devoted toward extending the basic nonlinear flux-correction techniques to convection algorithms other than SHASTA and toward discovering an “optimum” FCT algorithm, This latter effort [16] came to several conclusions: (1) There is an error inherent in finite-difference techniques which arises solely from the finite resolution of the grid and which does not vanish even when both numerical dispersion and numerical dissipation vanish. This error is two to three times smaller than achieved with the original FCT algorithms. (2) In local algorithms where dispersion and diffusion errors are present also, the residual dispersion errors seem to be more severe. Therefore, the optimal FCT
non for choosing the reduction factors of the antidiffusive fluxes is that the antidiffused solution will exhibit no new maxima or minima where the diffused solution had none. Although the limit (7b) represents the minimum amount to diffusion needed for stability, FCT algorithms generally use a larger zero-order diffusion [16] because it has been found that the correct choice of the ~v
1+112} within the monotonic and stable range will reduce convective phase errors from second to fourth order as suggested by Kreiss. Since the antidiffusion
___________________________________________ =
6
1 ~
~
CYLINDRICAL DIAPHRAGM PROBLEM [pI 4 I
SEC
-
HOCK
S 4
cell)
-
CONTACT DISCONTINUITY
)~3S cells)
-
algorithm was found to be one which has the diffusion coefficients = + and hence which drove
dispersion errors from second to fourth order [19]. (3) The best of the generally useful FCT variants was almost twice as accurate as the original FCT algorithm and within 50% of the nonzero minimal error, This is roughly six to eight times more accurate than the old methods. (4) Implicit convection algorithms, in addition to being somewhat slower than the local explicit methods to compute, also have bad nonlinear properties. Information transfer occurs faster than fluid characteristics hence shocks are relatively poorly treated. Since the latest FCT algorithms have eliminated
75
2
-
I
-
—200 CELLS )ir .-.
01)
100 CELLS 50 CELLS
C0
05
10 r—
15
20
-
Fig. 4. A cylindrical diaphragm problem using FASTiD. Three different resolution are shown to demonstrate convergence of the solution. FCT generally requires shock widths of about i.5 cells.
________________________________________ CYLINDRICAL DIAPHRAGM PROBLEM [~l = 41
I ‘0.8 SEC 4
,Ta-,,..3~
roughly 95% of the removable error and the removable error that remains is barely half of the irreducible error, it is natural to haye turned next toward optimization in speed, flexibility, and generality [31]. In FCT algorithms, the basic convective transport algorithm is augmented with a strong enough linear diffusion to ensure positivity at the expense of excess smoothing. Since the amount of diffusion which has been added is known, FCT then performs a conservative antidiffusion step to remove the diffusion in excess of the stability limit. However, the antidiffusive fluxes are effectively multiplied by a coefficient which ranges from zero to unity to preserve monotonicity. The crite-
SHOCK-/.
L
2
-
I
-
—200 CELLS (Sr ‘.01 n,CC 100 CELLS 000
C0.3 Fig. 5.
o~,~/\.~OTAC
-
DISCONTINUITY
-
I 06
0.7
St ‘.002)
50 CELLS I 04
A blow-up of the
I 0.5
fig. 4 calculation at 0.8 sec just before
the rebounding shock reaches the contact discontinuity. This figure demonstrates that the calculations differ from resolution but not due to accumulated secular truncation errors because the shock and the contact discontinuity are the same location in all three calculations.
76
J.P. Boris/Numerical solution of continuity equations
FPSST2D LPtSER SHELL
~
CC
I
I I
I -
I
I
I
/
I
// CC
C
______________________________________________
—
-
-
-
I
-
I
C
I
CII
—
CC
—
CI
Ii
I
~
-
I
1
-
~C’~
C
__*
.1
/ ppD
X
NRL. Fig. 6a.
can also chosen correspondingly larger, no real price is exacted for this improvement in phase properties. The most recent efforts [31] have taken advantage of this fact to generate minimum-operation-count FCT algorithms for cartesian, cylindrical, and spherical
in FORTRAN on the Texas Instruments Advanced Scientific Computer at the Naval Research Laboratory. The execution time per continuity equation per grid point is roughly 1.3 psec. A complete 2D calculation on a system of 200 grid pts X 200 grid pts requires
coordinate systems with stationary (eulerian) and
roughly 2—2.5 sec per timestep depending on extra
moveable (lagrangian) grid systems. Because these algorithms, and in particular the nonlinear flux-correction formula, were very carefully designed, they are fully “vectorizeable” for pipeline and parallel processing and have been implement in all generality
physics and boundary conditions incorporated in the problem. Figs. 4 and 5 show two aspects of a 1 D calculation performed by the code FAST1D on the ASC. The problem chosen is the “Lapidus” problem [32] in
J.P. Boris/Numerical solution of continuity equations
77
FPtST2D LPSSER SHELL
CC
C
I
C
-~
2.O8*1O~
L
/ /
-
J~•~l.i4.J~’ C’~C
II) C
CC
p I
IC
C
C~\C~1.ILwIiWW1IJl
C
C
-
C~
~
/
y -,
-
PPD
x—
_____________
Fig. 6b.
which eqs. (Sa)—(5d) are solved in cylindrical coordinates with ‘y = 1.4. The diaphragm is originally at r = 1.0 in fig. 4 and bursts at t = 0.0 sec. The density solution is shown at 0.6 sec in fig. 4 just after the shock has reached the axis and rebounded, Three different resolution calculations are overlapped to show the convergence. Even the calculation with 50 cells is at least as accurate as the original Payne and Lapidus solutions with 200 cells. The width of the contact discontinuity is about 3.5 cells while the shock is smeared over only 1.5 cells without noticeable
over- or undershooting. The calculations are performed without any added artificial viscosity. The monotonicity control provided by FCT on each of the continuity equations separately is adequate to ensure stability and accuracy as shown. Figs. 6(a)—(c) show the evolution of Rayleigh— Taylor instability in the implosion of a laser pellet shell 30 pm thick. The calculation was performed using the FAST2D code on the ASC and the thermal conductivity was set to zero to show the full nonlinear evolution of the instability. An Eulerian “sliding-rezone” grid was
78
J.P. Boris/Numerical solution of continuity equations
LASER SHELL
FAST2D
-
3,BS*10’
I
_____ -- CC~-)(
~.f)llI 1,lIj1i~I
—.
-
CC
C
//
I
CCCCCC
•-
C__CC I~C
C
)C
C
C
IC
CC
C
______________________
I~
I
CCC
C C
-
-
111111)
•I
I
I
—
I
I
III
C
CC
C
CC
C
—
I
I
)~I~
._,~.1
I
I
IC,..
CnIC
I
—
C
I
C
:
C
ii
I
)
/
till)
i
PPD
X—
NRL Fig. 6c.
Fig. 6. A series of three plots from a FAST2D calcuiation of a laser pellet shell being imploded (toward the left) by a strong constant 3. (a) pressure (high temperature— on the right. The thermal conductivity is set to zero to show the nonlinear evolution of the Rayleigh— After nsec theThe beginnings of theinRayleigh—Taylor canshell be seen thethick linearwith phase the back sideofof3 the Taylor1.05 instability. time is given the upper right ininstability seconds. The is 30inj.sm an on initial density g/cmshell and a sound/shock wave can be seen traveling inward along the “top” of the shell, (b) After 2.08 nsec the instability is beginning to enter the nonlinear phase and the sheil has moved to the left less than a single shell thickness. Lagrangian calculations would still be acceptable here. (c) After 3.85 nsec the shell has moved roughly two shell thickness and the Rayleigh—Taylor instability is fully developed. Usual lagrangian treatments would be breaking down here because of severe grid distortion.
used with 200 X 200 grid points. Strong deterioration of the shell is apparent by 3.85 nsec but breakthrough has not yet occurred. The jetting of material off the
backside is severe enough by this time that grid distortion would obviate the usual lagrangian solution techniques. Fig. 6(a) shows the linear phase of the instability and fig. 6(b) shows the onset of the nonlinear regime.
J.P. Boris/Numerical solution of continuity equations
79
Acknowledgement
[13] R. Chan, J. Comput. Phys. i7(3) (i974) 3ii—331. See also Hirt, Amsden and Cook, J. Comput. Phys. 14 (1974)
This work was sponsored in part by the Energy Research and Development Administration, Contract AEC-AT(49-20) and by the Office of Naval Research, Contract RRO14-03-02.
1141
References [1] F.H. Harlow, in: Methods in Computational Physics, eds. Alder, Fernbach and Rotenberg, Vol. 3 (Academic Press, New York, 1964) p. 319. [21 B. Marder, “GAP-A PlC Type Fluid Code”, Univ. of California, LASL Scientific Report. [3] J.P. Christiansen, J. Comput. Phys. 13(3) (1973) 363—379. [41 K.A. Taggart, R.L. Morse, R.L. McCrory and R.N. Remund, Bull. Amer. PhYSC Soc. 20(10) (1975) 1378. [5] C.K. Birdsail and D. Fuss, J. Chem. Phys. 3 (i969) 494. See also A.B. Langdon and B.F. Lasinski, in: Methods of Computational Physics, Vol. 16 (Academic Press, New York, 1976). [6] J.P. Boris and R.A. Shanny (Eds.), Proc. Fourth Conf. on Numerical Simulation of Plasmas, 2—3 Nov. 1970, U.S. GPO #0851—00059. See also R.L. Morse, in: Methods in Computational Physics, eds. Alder, Fernbach, and Rotenberg, Vol. 9 (Academic Press, New York, 1970). [7] R.W. Hockney, J. Assoc. Comput. Mach i2 (1965) 95. See also R.W. Hockney, Phys. Fluids 9 (1966) 1826 and 0. Buneman, J. Comput. Phys. 1(1967) 517. [8] A.G. Kuiikovskiy and GA. Lyubimov, Magnetohydrodynamics (Addison-Wesley, Palo Alto, 1965). [9] G. Moretti and M. Abbett, AIAA J. 4(i2) (i966). [10] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Interscience, J. Wiley and Sons, New York, 1967). [ii] C. Brennan and A.K. Whitney, Unsteady free surface flows; solutions employing the iagrangian description of the motion, 8th Symposium on Naval Hydrodynamics, Pasadena, Aug. 1970. [12] W.P. Crowley, Flag, a free-lagrange method for numerically simulating hydrodynamics flows in two dimensions, Proc. Second International Conference on Numerical Methods in Fluid Dynamics (Springer-Verlag, 1971).
MG. Wurtele Tellus XIII (3) (1961) 379—391. [15] P.J. Roache, Computational Fluid Dynamics (Hermosa Publishers, Albuquerque, 1972). [16] J.P. Boris and D.L. Book, J. Comput. Phys. 20 (1976) 000—000 [FCT IIII. [171 G. Strang and G. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, New Jersey, i973). [18] S.A. Orszag and M. Israel, Numerical flow simulation by spectral methods, in Proc. of the Symp. Numerical Models of Ocean Circulation, National Academy of Sciences (USGPO, Washington, D.C., 1972). [i9] H.O. Kreiss, A. comparison of numerical methods used in atmospheric and oceanographic applications, in Proc. of the Symp. Numerical Models of Ocean Circulation, National Academy of Sciences (USGPO, Washington, D.C., 1972). See also in the same volume B. Wendroff, Problems of accuracy with conventional finite-difference methods, and G. Fix, A survey of numerical methods for selected problems in continuum mechanics. [20] J.P. Boris, M.J. Fritts and Klaus Ham, Free surface hydrodynamics using a lagrangian triangular mesh, Proc. First Int. COflf. Numerical Ship Hydrodynamics (NBS, Gaitherburg, Md. Oct., i975). [21] M.J. Fritts, SAl Report SAI-76-528-WA, March (i976). [221 B. VanLeer, J. Comput. Phys. 14 (1974) 36 1—370. [23] J.P. Boris and D.L. Book, in: Methods in Computational Physics, Vol. i6 (Academic Press, New York, 1976) ch. 11. [24] J.P. Boris and DL. Book, J. Comput. Phys. 11(1973) 38ff [FCT I]. [25] DL. Book, J.P. Boris and K.H. Ham, J. Comput. Phys. 18(3) (1975) 248—283 [FCT II]. [26] A. Harten, CIMS Report C00-3077-50, June (1974). [27] B. VanLeer, MUSCL, a new approach to numerical gas dynamics, Proc. of the Second European Conf. Computational Physics, Garching, 27—30 Apr. (i976). [28] J.P. Boris, NRL Memorandum Report 2542 (1972). [29] J.P. Boris, J,H, Gardner and S. Zalesak, NRL Memorandum Report 308i (1975). [30] R.A. Chevalier and J.H. Gardner, Astrophys. J. 192 (1974) 427. [31] J.P. Boris, NRL Memorandum Report 3237 (1976) to be published. [32] A. Lapidus, J. Comput. Phys. 8 (1971) i06—i18.
NUMERICAL SOLUTION OF CONTINUITY EQUATIONS J.P. BORIS Plasma Physics Division, Naval Research Laboratory, Washington, D.C. 20375, USA
This paper compares and contrasts the six general methods available for solving the time dependent continuity equation numericaily. Since finite-difference techniques naturally arise as the most generally applicable solution methods, consideration is given to solving the major remaining difficulties encountered in the general application of lagrangian and eulerian finite differences. In the former, triangular mesh ceils with free reconnections are suggested as a complicated but serviceable solution to the lagrangian grid distortion problem. In the latter, nonlinear monotonic difference schemes are introduced as the natural way to avoid the eulerian positivity problem. Flux-corrected transport techniques, the most developed and efficient of the positive (monotone) finite-difference approaches, are described with examples.
1. The continuity equation
~fpd3r_fpv. This paper reviews the different numerical techniques which have been developed to solve the continuity equation ap/at =
—
V
(conservation form).
pV
(la)
In section 2 the six different available methods are compared and categorized. Within a framework provided by six requirements which ought to be satisfied by an ideal continuity equation algorithm, the strengths and weaknesses of the different approaches are analyzed. Finite-difference techniques have particularly attractive properties in this framework, thus section 3 is devoted to solving, at least partially, the remaining outstanding difficulties which finite-difference techniques suffer. Section 4 is devoted to flux-corrected transport (FCT) solutions of the nonlinear continuity equation which preserve the positivity property of the real equations and yet have relatively low numerical diffusion. In addition to eq. (1a), there are two other ways to write the continuity equation which get reflected in some of the numerical solution techniques: ap
dp V Vp ~j~convection
and
=
—
pV V,
compression (convection form)
(lb)
region
dA.
(lc)
region boundary (integral form)
The convective term V Vp displayed explicity in eq. .
.
.
(1b) gives eqs. (Ia)—(l c) their intrinsically hyperbolic form and causes really severe problems numerically. The compression term —pV V is sometimes absent, as in the Liouvifie equation, in which case it is often called the convection or the advection equation. The continuity equation appears in almost all descriptions of dynamic physical systems simply because it expresses two of the general principles in physics, conservation and causality. The integral form given by eq.(lc) states that the amount of material (f pd3r)in a given region (conservation) can only change due to flows of material (p V dA) through the boundary. In other words, any material currently at a given point could only have gotten there by leaving someplace else (causality). Continuity equations underlay compressible and incompressible fluid dynamics, hydrodynamics. plasma physics (Vlasov equation and MHD moment equations) and evenequation quantumdisplays mechanics. The continuity the positivity property that a quantity being transported will never .
turn negative anywhere in a reasonable flow field if it was everywhere positive to start with. This positivity
68
J.P. Boris/Numerical solution of continuity equations
property expresses in a continuum way the intrinsic corpuscular nature of matter. Thus matter cannot be removed from a region which is devoid of matter to begin with. This duality, in which matter obeys both particle and fluid-like equations on microscopic and macroscopic scales, respectively, has its ramifications for numerical solution techniques as well. The first three of the techniques described below take advantage of the underlying discrete basis of the continuity equat~ionwhile the last three aim more directly at solving the partial differential equation itself, The continuity equation is hyperbolic because it has first derivatives in both space and time and yet really is quite different from wave equations which are also hyperbolic. Therefore the solution techniques can also differ quite significantly. In the continuity equation the quantity being propagated appears in both the time and space derivative terms and could be described as “self”-hyperbolic. Generally the characteristics (flow) are in one direction only. Wave equations usually occur in pairs with two oppositely directed characteristics and the two quantities being propagated (say E and B in electromagnetics) appear in the time derivative of one equation and in the space derivative of the other. Thus wave equations might best be classed’ as “mutually” hyperbolic. There is no positivity to speak of and the most useful conservation property is quadratic rather than linear. Because of these distinctions, much of what is said here about continuity equations will not apply simply to wave equations.
2. Numerical solution techniques This survey considers only conservative solution techniques although nonconservative techniques are adequate for some problems. In most situations the conservation properties of the physics have to be mirrored in the numerical method or else nonphysical instabilities or unacceptable secular errors will result. Although positivity is probably as important as conservation, not all of the generally applied methods have recognized this fact to date and thus the survey cannot be limited to positive methods. It ought to be, however, since the most common and annoying problems which arise in the solution of coupled systems of continuity equations arise because the quantity convected becomes falsely negative somewhere or becaLse
numerical oscillations occur in the vicinity of steep gradients and shocks. The analysis is also limited to discussing the basic methods even though there are a host of hybrid techniques combining two or more of the basic approaches. The number of such combinations is staggering and even the definitions of the basic methods are in some dispute. Therefore no large useful purpose is served by trying to classify and analyze all of the possible combinations. Rather the general properties of each are described and the pros and cons considered. It is left intentionally to the reader to select the best method or hybrid of methods for his own particular application. As a framework for evaluating the algorithms, I present six general requirements which an ideal numerical algorithm for solving the continuity equation should satisfy to be generally useful. These requirements are given in what I estimate to be their order of importance. An ideal algorithm should: (1) be linearly stable for all cases of interest; (2) mirror conservation properties of the physics; (3) ensure the positivity property when appropriate; (4) be reasonably accurate; (5) be computationally efficient; and (6) be independent of specific properties of one particular application. This last requirement is a restatement of the facts that we are seeking generally useful methods and that all hybrid schemes cannot be discussed. As with the methods themselves, some of the requirements are related. Stability, conservation, and positivity generally relate to accuracy, for example, and yet no one of these can guarantee any of the others. On the other hand the requirements can also be partially contradictory. Accuracy and efficiency often tug in opposite directions, for example. Table I show the six methods I wish to discuss, the order being chosen to range from most discrete (or particulate) in nature to most continuous (or global). In their simplest form quasi-particle methods attempt to solve eq. (1) by simulating the microscopic physics of the particles which make up the fluid. If a large number N of quasi-particles are initialized in the calculation at positions {X 1(O)} at time t = 0, their later positions can be found at time t by integrating the simple orbit equation forward dX.(t)/dt = v.(t) (2)
J.P. Boris/Numerical solution of continuity equations Table 1 Possible solution techniques
______________
• Quasi-particle methods (1) Collisionless particles, stars and plasmas (2) Collisional particles for fluids • Characteristic methods • Lagrangian finite difference mcthods • Eulerian finite difference methods (1) Explicit versus implicit (2) Order versus accuracy •The finite element method • The spectral method
assuming that the velocity of each particle is known as a function of time. Eq. (2) is an ordinary differential equation and is usually solved easily with a high degree of accuracy. When each quasi-particle has a velocity which is determined solely as a function of the particle’s position at any time, i.e. ~) V1r~ V1X~‘ ~ ‘ —
‘
the quasi-particles cannot pass each other, except as a result of finite time-step errors [or perhaps singularities in V(X, t)]. The model is then appropriate for collisional fluids. When several quasi-particles can occupy essentially the same location in space and yet be moving with significantly different velocities, the model is more appropriately used to describe “collisionless” systems such as plasmas on the microscopic scale or stars in self-gravitating systems on a somewhat grander scale. Such “collisionless” systems of course undergo “collisions” when the impact parameters of two real particles (stars) becomes sufficiently small but such occurrences are rare. When the force law for the interacting “collisionless” particles is inverse square, it can be shown that the major component of the force is coherent and long range. Thus, most solutions of this type are seeking statistical and average properties of the flow and the individual particle orbits are relatively unimportant. The collisional quasi-particle methods which have been tried to date include PlC [1], GAP [2], and VORTEX [3] and more recent attempts of the same genre by workers at Los Alamos [4]. These methods introduce what amounts to a finite-difference grid for some of the derived physical quantities such as pressure while retaining the particle’s mass, momentum, and energy as particle variables. The collisionless quasi-
69
particle methods are reviewed by Birdsall and Langdon [5] and in the Proceedings of the Fourth Conference on the Numerical Simulation of Plasma [61.The .
original applications seem to stem from work by Hockney and Buneman [7] The term quasi-particle has been used because it is almost never possible to carry in the computer memory as many particles as are involved in the actual physical problem of interest. It is not unusual to simulate systems of lO~particles (galaxies) to 1019 particles (I cc of dense plasma) using 103_.106 quasi-particles. In fact, when the number of quasi-particles exceeds ~ or so, it is not usually possible to treat the forces by calculating a direct interaction law. A potential equation is interjected between the source terms (charge or mass) and the corresponding long range fields (electromagnetic or gravitational) and solved on a finite-difference grid in either configuration or transform space. The quasi-particle methods are generally stable, and conserveconservation though where they of should both(almost energy and by definition) momentumalis difficult in the collisionless models. These methods can also guarantee positivity since it is just as hard to remove a quasi-particle from an already empty region of space as it is a real particle. The methods can even be made reasonable general and flexible. Where these quasi-particle methods fall down is in efficiency and accuracy. Because the statistics and hence the smoothness of the computed solutions depend on the number of quasi-particles which the user can afford, many applications of interest are beyond the particle approach. Long-time, turbulent, collisional flows are unbearably expensive and collisionless “plasma” simulations usually break down for calculations which are run longer than a collision time for the quasi-particles. Some workers like to think of plasma simulation as an unavoidably collisional approximation to the collisionless Vlasov equation. The Vlasov equation,in reality, is a collisionless approximation to a real collisional plasma. Characteristic methods, which take advantage of the physical content of the specific system of eqautions being solved, are not entirely counter to our requirement (6) because different characteristic methods can often be formulated for different sets of governing equations. The equations of ideal MHD [8] and ideal compressible flow [9] for example, can be formulated and solved by characteristics and these solutions are often the best that can be obtained because the discon,
70
J.P. Boris/Numerical solution of continuity equations
tinuities which arise naturally in these physical systems are followed individually like particles in the previous methods and the continuous fluid behavior between the characteristics of the discontinuities is easy to represent. Nevertheless, characteristic methods are not generally applicable even though conservation and positivity are relatively easy to ensure. The presence of even one diffusion term in the governing equations usually abrogates the use of characteristic methods because the characteristics of the model usually cease to exist, A further drawback to characteristic methods is their complexity and relative inefficiency for multidimensional flows. The characteristic interfaces are complicated lines in two dimensions and surfaces in three dimensions and these can become progressively more and more difficult to follow as the flows proceed in time. This relative inefficiency, which is not troublesome in simple flows where a few characteristics may be sufficient to describe a whole flow, becomes far more serious when the new parallel and pipeline computers are contemplated. Since the characteristics, like the quasi-particles, must each be considered individually and since aspects of random access to computer memory as well as complicated logic are involved, the future of these methods on the newer computers is correspondingly dim. We come now to a discussion of finite differences. Unlike the quasi-particle and characteristic methods, finite-difference solutions of the continuity equation are based on an approximation to eq. (1) in which the continuous flow parameters p(x, t) and V(x, t) are discretized to a finite set of representative values ~p1},
classic text is Richtmyer and Morton [10] but the literature is virtually limitless. Although a detailed classification of possible finite difference schemes is also probably purposeless (there are two many hybrids) and may be impossible, a few of the possible choices should be discussed. The first of these is the distinction between lagrangian and eulerian methods. A lagrangian finite-difference scheme [11—13] is one in which the location of the grid points is allowed to move along with the fluid. This approach is useful because the convection term is transformed away in the moving system and hence the positivity property is easily ensured as well as conservation. The trouble with this approach is that complicated flows in two or three dimensions will rather quickly distort the moving mesh. Much of the accuracy is lost, often before useful information can be extracted from the calculation. Eulerian methods, on the other hand, keep the grid fixed in position and hence distortion of the mesh can not occur. Since fluid is now flowing across a stationary and usually rather coarse mesh, truncation errors occur which are secular in time and which can play havoc ~yithaccuracy and positivity in a long run [l0,14—16]. Nevertheless, these problems with both the eulerian and lagrangian approaches can be overcome as described in the next section. Here I wish to point out that the lagrangian finite-difference approaches, with the mesh distortion problem cleared up, are not completely general because aspects of eulerian flow appear in realistic coupled fluid systems which are being treated on a lagrangian grid.
defined at N distinct points in space called grid (or mesh) points. The desired continuous functions are assumed to be known only at the mesh points and a whole culture of techniques have been developed to predict future values of fp~}at some time t given ~p1}at t = 0. The term finite differences refers to the fact that the spatial derivatives (and usually the time derivative well)e.g. are approximated using only the grid point as values ~o• —i—-—, (4)
a /at =
—
ax X
__!__
1
X1+i
~i i
The attraction of these finite-difference methods is their richness and simplicity; no one need use the same finite-difference scheme twice in a row. The
Consider the ideal compressible fluid equations. V V (Sa) p apV/at = V- pVV— VP (Sb) —
—
ac/at =
—
V eV— V PV,
(Sc)
where 2 (Sd) e = [P/(7 1)] + ~-pV Ifwe define V~= [(c +P)/e] Vas a new velocity which varies from Vto about ‘yVand which is always parallel .
—
to
v, the energy continuity equation can be written as
ac/at =
—
and is in
V CV*
(5d’)
exactly the same form as eq. (5a) for the mass
J.P. Boris/Numerical solution of continuity equations
density. Since V V’~generally, a mesh moving at V is lagrangian for ~nassdensity, but not for energy density while a mesh moving at V’ is lagrangian for the energy density but not for the mass density. This means that full advantage can never be taken of a lagrangian grid; some aspects of a realistic calculation are always eulerian. Thus the eulerian grid problems with finite differences have to be controlled or solved for lagrangian, as well as eulerian, approaches to be fully viable, The second choice which I wish to consider is the time integration technique, in particular the neverending controversy of choosing an explicit versus an implicit scheme. Implicit schemes for the convective derivatives involve, in the worst cases, gigantic matrix equation solutions and in the best cases the inversion of a tridiagonal matrix. Thus, their operation count and hence speed is appreciably worse than obtained using corresponding explicit schemes. The advantage of implicit schemes professed to outweight the dual disadvantages of program complexity and operation count per timestep is the ability to take much longer finite-difference timesteps without exciting numerical modes of instability. The answer lies in the physics (mathemetics) of the specific problem being solved, When the physical phenomenon of interest varies appreciably on a timescale r, no calculation with & r can reasonably be claimed to reproduce the phenomenon accurately. In such cases computationally efficient explicit methods are called for with L~t
71
forward differenced beyond second-order centered differencing tend to damp the undesired high frequency sound waves strongly. Thus, sonic energy from the whistle gets deposited as heat near the whistle by a conservative, forward-differences implicit scheme. A centered scheme, or quite near centered, would allow most of the sonic energy to propagate away as sound without heating the local air. The strongly dispersed short wavelength sonic components would generally stay in the vicinity of the whistle far too long, however, and might eventually deposit their energy near the whistle anyway. The best solution for such a quandary is usually either to rescale the problem so the two timescales are not so disparate or to reformulate the problem entirely so the two timescales are separated entirely. In the whistle problem this would involve a multiple timescale analysis in which the high frequency oscillations are averaged out. Such a reformulation is not always available in which case neither the implicit nor the explicit methods can be judged fully satisfactory. Both finite-element and spectral methods utilize continuum concepts to the limit in solving eq. (1). The functions of interest, p(X, t) and V(X, t), are expanded in a set of useful basis functions whose expansion coefficients are now determined from a complicated set of coupled, usually nonlinear, but ordinary differential equations. In the case of finite elements [17] the basis functions are localized to a mesh cell or two and the techniques of matching expansion function boundary conditions at cell boundaries reminds one very much of using splines. The spectral methods generally employ a series of global basis functions [18] but otherwise the general characteristic is a projection to recover ordinary differential equations which are implicit in one way or another. The finite element methods acquire their implicitness by coupling of spatially adjacent basis functions via projection when evaluating the left-hand side a/at of eq. (1). The spectral methods become effectively implicit when the huge convolutions are needed to evaluate, for example, V(~V/OX) where the velocity field is represented as a Fourier expansion. The disadvantages of these two approaches are clear: (1) they are slow per timestep compared with corncomparable sized finite-difference methods for general problems; (2) they clearly are ill-suited (spectral methods more so than finite element) to problems displaying
72
J.P. Boris/Numerical solution of continuity equations
highly localized effects such as shocks; and (3) they are generally complicated and expensive to program, particularly so for nonlinear terms, The advantages are greater accuracy by far then simple difference or particle methods for smooth well-behaved flows and oscillations, The spectral methods have been described as “infinite-order” accurate and have been characterized as “optimal” for certain finite-difference evaluations [16]. Whether or not the potential increased accuracy outweight the drawbacks for real problems using either finite elements or spectral methods depends very strongly on the specific problem being solved. While the degree to which these two classes of techniques satisfy or do not satisfy our requirements (5) and (6) may be subject to legitimate debate, there is no doubt that the positivity property of the continuity equation and accuracy cannot be guaranteed when expanding in a time-independent set of basis functions, The Gibbs phenomenon [16] represents a minimal, or lower-bound error, which is only reduced by increased resolution. In a very real philosophical sense an uncertainty principle is at work. The numerical error can only be reduced by increasing the number of discrete degrees of freedom in the representation. Before giving my own evaluations, I would like to quote three conclusions, in a similar context, from another source [19] (1) The spectral method has no stability problems but is much more complicated and slower than generalized difference methods. (2) It is doubtful whether the finite-element method, based on piece-wise polynomials, can compete with the above methods. (3) If difference methods are used, they should be at least fourth-order accurate. While (I agree basically with these remarks, they do not really encompass the quasi-particle schemes nor do they adequately reflect a very important piece of personal experience. Whenever a theory is to be tested, or an algorithm, or a new mathematical technique, attention rather quickly turns to the crucial yet simple conservation equations of ideal compressible flow, eqs. (5a)— (Sd). Now I know that finite-difference methods have solved the transient Rankine—Hugoniot shock problem adequately, and I even know that the various ‘f~vors of quasi-particle methods can do the same vital problem
creditably if enough particles are used. I do not know of any calculation, even in one spatial dimension, using either a finite element or a spectral method which has correctly solved for an ideal gas compressible shock. Until such calculations become common place and computationally attractive, finite-difference methods would seem to have inside track, Therefore, the next two sections take a closer look at the remaining deficiencies in lagrangian and eulerian finite-difference techniques.
3. Improving finite-difference techniques There is undoubtedly some merit in trying to boost the performance of quasi-particle methods on the one hand and the basis-function expansion methods on the other toward the level obtained by finite differences. But it is a low risk—high return investment to patch up the obvious failings of the front runners finite differences. In the case of lagrangian methods, the major outstanding problems arise from secularly unattractive distortions of the grid which wreck calculations of interesting flows quite quickly. In the case of —
. . Fig. 1. A Schematic illustration o! the use of a two-dimensional grid of triangles to represent a iarge amplitude, separating, free-surface gravity wave.
J.P. Boris/Numerical solution of continuity equations
eulerian methods, the major outstanding weakness in a huge class of problems of real interest is the need for a large artificial damping (numerical diffusion) to fill in what would otherwise be pits of “negative density” in the calculated profiles. Since the “eulerian” positivity problem is encountered even in Lagrangian calculations for many situations, it demands the greater share of attention. The distortion problem for lagrangian grids can be handled adequately by using a finite-difference mesh of triangles rather than rectangles (tetrahedrons hi three dimensions). Fig. 1 shows the representational advantage of using triangles for a breaking gravity wave. This advantage, cribbed from engineers and original finite element buffs who have been using meshes of triangles for years, was adopted to finite differences ofincompressible flow [12,20]. Thenoticeableimprovements in performance arise for two reasons: first, the general connectivity premitted by triangles allows smooth representation of much more complicated shapes than can be treated smoothly by an equal number of rectangular cells, and secondly since the number of lines meeting at a vertex is not fixed, it grid distortions. can be varied during a calculation to prevent severe While this approach may not be the only method of dealing with severe grid distortions (see, for example, ref. [13]), it is the most physical and the least dissipative technique employed so far. Fritts [21] has recently employed lagrangian triangular grid techniques in studying nonlinear aspects of free-surface waves including strongly sheared flows. Ten to fifteen full cycles of the waves can be integrated reversibly without significant deterioration of the solution and the reconnection procedure in the presence of strong shear is reversible (i.e. nondissipative). Fig. 2 shows a plot of the wave natural period versus grid size for four different resolution calculations of the same physical problem. The second-order accuracy over long times is demonstrated by the parabolic form of the curve and extrapolation to zero stepsize gives improved accuracy as expected. An added advantage accrues from using triangular cells. Since it is difficult, if not impossible, to label the vertices unambiguously as either even or odd, the nonlinear mesh-separation instabilities, which plague low-dissipation rectangular cell techniques, seem to be absent or damped with triangular cells.
~
73
Bt ~.OO4
,
~
,
.
8t~.OO4~
,
/ ,
/
-
8m~.OO2
,
.
.,
~
_____
I
~i*
~‘~‘~‘;;~‘ ~.i
‘I
ii
ft .
tt~Y*ltttt~i ~
~l.i,% ~
— ~—.- ‘
tti
i
~__L__.~
~ ~—~---.~
-C
-
Hg, 2a. 138
-
136
r
~
-412732 + 00020
134
132
-
130
-
128 126
-
I
I
.1
2
.3
Fig. 2b. Fig. 2. Wave period versus normalized numerical grid size for four different resolution calculations of a finite depth freesurface wave.
The remaining drawbacks to the triangular grid approach are the complexity of the programming required and the fact that the operations which are rerequired seem to hinge strongly on linked lists, random access, and sequential processing. With several pipeline and vector processors working and several more soon to be in operation, efforts toward the “vectorization” for parallel processing of these triangular-grid techni-
74
J.P. Boris/Numerical solution of continuity equations
ques would seem to be well worthwhile. The major remaining problem with eulerian techniques for solving continuity equations is the competition between accuracy and positivity [19,22] wherein fully second order techniques cannot maintain positivity Consider the rather general three-point approximation to eq. (la)
SIDED
~
WE +Vj+i/2(p~l—PiQ)—i)fl/2(p~—~1),
where ei+i/2
~+i/2/~~+i/2’
1~11~>~/+1/2.
2
cycle
00
I
=
60
-
I
~—
—
X (cell no.)
~—~
—
X (cell no.)
---~
21
•~--
L x Ice); ~)
O~cell
~)
r;
2 -~
X (cell no)
00
—-~
•-~
i 0
._
____________
2
FROG
-L
H
I
--
~
-___________________
LEAP-
L_
— ~
I....:
-
P -
s—
2t
cycle 800
2
-
~
(6)
~.
100 0
----~
-i
X (cell no.)
100 -
(7a)
The upper limit arises from the explicit diffusion timestep condition [10] while the lower limit is the Lax— Wendroff damping [23]. Unfortunately positivity is only ensured linearly when 1 61+i/2 (7b) >i the first-order upstream-centered scheme result, We appear to be caught between a rock and a hard place here but the escape route is signaled in the preceding sentence by the word “linearly”. By relaxing the linearity implied by eq. (6) and letting the diffusion coefficients be nonlinear functionals of the flow velocities {c/+i/2}, we can hope to reduce the integrated dissipation below the rather ghastly limit (7b) and yet retain sufficient dissipation near steep gradients to ensure positivity. A literature is beginning to form about these “monotonic difference schemes” [22,24—26] since the dilemma of accuracy versus positivity in eulerian difference schemes can be resolved in no other way. Several different approaches are possible and the area is still largely unexplored. Boris and Book introduced this basic nonlinear approach with the techniques of flux-corrected transport [16, 23—25]. The rewards were quickly seen. Fig. 3 shows a comparison of four common difference schemes solving the standard square wave problem [14,24]. The effects of excess numerical damping in the donor-cell treatment (upstream centered first order), and of excess dispersion in the leapfrog ,
ONE-
Eq. (6) is in finite-
difference conservative form with whole indices re presenting cell centers and half indices indicating cell interfaces. The additional numerical diffusion terms with diffusion coefficients {~/+i/2} have to be added to ensure positivity The stability of eq (6) is ensured at least roughly, when ~ >v
f r 20
SHASTA 2~ ~
2~
(FCT)
.1 -
X (cell no.) —~
E 100 0
x Ice II
no.)
.
..
100
Fig. 3. Comparison of four difference schemes solving the square wave problem. The merits of flux-corrected transport relative to the other methods are clear.
and Lax—Wendroff treatments are clearly visible. Dispersion manifests itself as a trail or projection of oscillations in the computed solution near discontinuities and sharp gradients of the “correct” solution. VanLeer [27] through geometric arguments, and Harten [26] by a more analytic approach, have developed other monotonic finite-difference schemes for avoiding the accuracy— positivity dilema, Clearly the study of monotonic schemes will be fruitful because of the central role that eulerian finite-difference techniques naturally play. ,
4. Flux-corrected transport algorithms The first, and so far the most developed and used, of the monotonic schemes is flux-corrected transport (FCT). The calculation in fig. 3 was performed by the first FCT Algorithm SHASTA [24] and had an error about four or five times smaller than the simple linear methods also shown. The damping was second order as
J.P. Boris/Numerical solution of continuity equations
were the relative phase errors [161. This algorithm was used widely in the two-dimensional ideal compressible fluid dynamics code SHAS2D [28]. The basic techniques was quickly generalized to cylindrical and spherical systems, to lagrangian as well as fixed eulerian grids, and was applied to a number of one-, two-, and threedimensional problems [29,30]. More recent work has been devoted toward extending the basic nonlinear flux-correction techniques to convection algorithms other than SHASTA and toward discovering an “optimum” FCT algorithm, This latter effort [16] came to several conclusions: (1) There is an error inherent in finite-difference techniques which arises solely from the finite resolution of the grid and which does not vanish even when both numerical dispersion and numerical dissipation vanish. This error is two to three times smaller than achieved with the original FCT algorithms. (2) In local algorithms where dispersion and diffusion errors are present also, the residual dispersion errors seem to be more severe. Therefore, the optimal FCT
non for choosing the reduction factors of the antidiffusive fluxes is that the antidiffused solution will exhibit no new maxima or minima where the diffused solution had none. Although the limit (7b) represents the minimum amount to diffusion needed for stability, FCT algorithms generally use a larger zero-order diffusion [16] because it has been found that the correct choice of the ~v
1+112} within the monotonic and stable range will reduce convective phase errors from second to fourth order as suggested by Kreiss. Since the antidiffusion
___________________________________________ =
6
1 ~
~
CYLINDRICAL DIAPHRAGM PROBLEM [pI 4 I
SEC
-
HOCK
S 4
cell)
-
CONTACT DISCONTINUITY
)~3S cells)
-
algorithm was found to be one which has the diffusion coefficients = + and hence which drove
dispersion errors from second to fourth order [19]. (3) The best of the generally useful FCT variants was almost twice as accurate as the original FCT algorithm and within 50% of the nonzero minimal error, This is roughly six to eight times more accurate than the old methods. (4) Implicit convection algorithms, in addition to being somewhat slower than the local explicit methods to compute, also have bad nonlinear properties. Information transfer occurs faster than fluid characteristics hence shocks are relatively poorly treated. Since the latest FCT algorithms have eliminated
75
2
-
I
-
—200 CELLS )ir .-.
01)
100 CELLS 50 CELLS
C0
05
10 r—
15
20
-
Fig. 4. A cylindrical diaphragm problem using FASTiD. Three different resolution are shown to demonstrate convergence of the solution. FCT generally requires shock widths of about i.5 cells.
________________________________________ CYLINDRICAL DIAPHRAGM PROBLEM [~l = 41
I ‘0.8 SEC 4
,Ta-,,..3~
roughly 95% of the removable error and the removable error that remains is barely half of the irreducible error, it is natural to haye turned next toward optimization in speed, flexibility, and generality [31]. In FCT algorithms, the basic convective transport algorithm is augmented with a strong enough linear diffusion to ensure positivity at the expense of excess smoothing. Since the amount of diffusion which has been added is known, FCT then performs a conservative antidiffusion step to remove the diffusion in excess of the stability limit. However, the antidiffusive fluxes are effectively multiplied by a coefficient which ranges from zero to unity to preserve monotonicity. The crite-
SHOCK-/.
L
2
-
I
-
—200 CELLS (Sr ‘.01 n,CC 100 CELLS 000
C0.3 Fig. 5.
o~,~/\.~OTAC
-
DISCONTINUITY
-
I 06
0.7
St ‘.002)
50 CELLS I 04
A blow-up of the
I 0.5
fig. 4 calculation at 0.8 sec just before
the rebounding shock reaches the contact discontinuity. This figure demonstrates that the calculations differ from resolution but not due to accumulated secular truncation errors because the shock and the contact discontinuity are the same location in all three calculations.
76
J.P. Boris/Numerical solution of continuity equations
FPSST2D LPtSER SHELL
~
CC
I
I I
I -
I
I
I
/
I
// CC
C
______________________________________________
—
-
-
-
I
-
I
C
I
CII
—
CC
—
CI
Ii
I
~
-
I
1
-
~C’~
C
__*
.1
/ ppD
X
NRL. Fig. 6a.
can also chosen correspondingly larger, no real price is exacted for this improvement in phase properties. The most recent efforts [31] have taken advantage of this fact to generate minimum-operation-count FCT algorithms for cartesian, cylindrical, and spherical
in FORTRAN on the Texas Instruments Advanced Scientific Computer at the Naval Research Laboratory. The execution time per continuity equation per grid point is roughly 1.3 psec. A complete 2D calculation on a system of 200 grid pts X 200 grid pts requires
coordinate systems with stationary (eulerian) and
roughly 2—2.5 sec per timestep depending on extra
moveable (lagrangian) grid systems. Because these algorithms, and in particular the nonlinear flux-correction formula, were very carefully designed, they are fully “vectorizeable” for pipeline and parallel processing and have been implement in all generality
physics and boundary conditions incorporated in the problem. Figs. 4 and 5 show two aspects of a 1 D calculation performed by the code FAST1D on the ASC. The problem chosen is the “Lapidus” problem [32] in
J.P. Boris/Numerical solution of continuity equations
77
FPtST2D LPSSER SHELL
CC
C
I
C
-~
2.O8*1O~
L
/ /
-
J~•~l.i4.J~’ C’~C
II) C
CC
p I
IC
C
C~\C~1.ILwIiWW1IJl
C
C
-
C~
~
/
y -,
-
PPD
x—
_____________
Fig. 6b.
which eqs. (Sa)—(5d) are solved in cylindrical coordinates with ‘y = 1.4. The diaphragm is originally at r = 1.0 in fig. 4 and bursts at t = 0.0 sec. The density solution is shown at 0.6 sec in fig. 4 just after the shock has reached the axis and rebounded, Three different resolution calculations are overlapped to show the convergence. Even the calculation with 50 cells is at least as accurate as the original Payne and Lapidus solutions with 200 cells. The width of the contact discontinuity is about 3.5 cells while the shock is smeared over only 1.5 cells without noticeable
over- or undershooting. The calculations are performed without any added artificial viscosity. The monotonicity control provided by FCT on each of the continuity equations separately is adequate to ensure stability and accuracy as shown. Figs. 6(a)—(c) show the evolution of Rayleigh— Taylor instability in the implosion of a laser pellet shell 30 pm thick. The calculation was performed using the FAST2D code on the ASC and the thermal conductivity was set to zero to show the full nonlinear evolution of the instability. An Eulerian “sliding-rezone” grid was
78
J.P. Boris/Numerical solution of continuity equations
LASER SHELL
FAST2D
-
3,BS*10’
I
_____ -- CC~-)(
~.f)llI 1,lIj1i~I
—.
-
CC
C
//
I
CCCCCC
•-
C__CC I~C
C
)C
C
C
IC
CC
C
______________________
I~
I
CCC
C C
-
-
111111)
•I
I
I
—
I
I
III
C
CC
C
CC
C
—
I
I
)~I~
._,~.1
I
I
IC,..
CnIC
I
—
C
I
C
:
C
ii
I
)
/
till)
i
PPD
X—
NRL Fig. 6c.
Fig. 6. A series of three plots from a FAST2D calcuiation of a laser pellet shell being imploded (toward the left) by a strong constant 3. (a) pressure (high temperature— on the right. The thermal conductivity is set to zero to show the nonlinear evolution of the Rayleigh— After nsec theThe beginnings of theinRayleigh—Taylor canshell be seen thethick linearwith phase the back sideofof3 the Taylor1.05 instability. time is given the upper right ininstability seconds. The is 30inj.sm an on initial density g/cmshell and a sound/shock wave can be seen traveling inward along the “top” of the shell, (b) After 2.08 nsec the instability is beginning to enter the nonlinear phase and the sheil has moved to the left less than a single shell thickness. Lagrangian calculations would still be acceptable here. (c) After 3.85 nsec the shell has moved roughly two shell thickness and the Rayleigh—Taylor instability is fully developed. Usual lagrangian treatments would be breaking down here because of severe grid distortion.
used with 200 X 200 grid points. Strong deterioration of the shell is apparent by 3.85 nsec but breakthrough has not yet occurred. The jetting of material off the
backside is severe enough by this time that grid distortion would obviate the usual lagrangian solution techniques. Fig. 6(a) shows the linear phase of the instability and fig. 6(b) shows the onset of the nonlinear regime.
J.P. Boris/Numerical solution of continuity equations
79
Acknowledgement
[13] R. Chan, J. Comput. Phys. i7(3) (i974) 3ii—331. See also Hirt, Amsden and Cook, J. Comput. Phys. 14 (1974)
This work was sponsored in part by the Energy Research and Development Administration, Contract AEC-AT(49-20) and by the Office of Naval Research, Contract RRO14-03-02.
1141
References [1] F.H. Harlow, in: Methods in Computational Physics, eds. Alder, Fernbach and Rotenberg, Vol. 3 (Academic Press, New York, 1964) p. 319. [21 B. Marder, “GAP-A PlC Type Fluid Code”, Univ. of California, LASL Scientific Report. [3] J.P. Christiansen, J. Comput. Phys. 13(3) (1973) 363—379. [41 K.A. Taggart, R.L. Morse, R.L. McCrory and R.N. Remund, Bull. Amer. PhYSC Soc. 20(10) (1975) 1378. [5] C.K. Birdsail and D. Fuss, J. Chem. Phys. 3 (i969) 494. See also A.B. Langdon and B.F. Lasinski, in: Methods of Computational Physics, Vol. 16 (Academic Press, New York, 1976). [6] J.P. Boris and R.A. Shanny (Eds.), Proc. Fourth Conf. on Numerical Simulation of Plasmas, 2—3 Nov. 1970, U.S. GPO #0851—00059. See also R.L. Morse, in: Methods in Computational Physics, eds. Alder, Fernbach, and Rotenberg, Vol. 9 (Academic Press, New York, 1970). [7] R.W. Hockney, J. Assoc. Comput. Mach i2 (1965) 95. See also R.W. Hockney, Phys. Fluids 9 (1966) 1826 and 0. Buneman, J. Comput. Phys. 1(1967) 517. [8] A.G. Kuiikovskiy and GA. Lyubimov, Magnetohydrodynamics (Addison-Wesley, Palo Alto, 1965). [9] G. Moretti and M. Abbett, AIAA J. 4(i2) (i966). [10] R.D. Richtmyer and K.W. Morton, Difference Methods for Initial Value Problems (Interscience, J. Wiley and Sons, New York, 1967). [ii] C. Brennan and A.K. Whitney, Unsteady free surface flows; solutions employing the iagrangian description of the motion, 8th Symposium on Naval Hydrodynamics, Pasadena, Aug. 1970. [12] W.P. Crowley, Flag, a free-lagrange method for numerically simulating hydrodynamics flows in two dimensions, Proc. Second International Conference on Numerical Methods in Fluid Dynamics (Springer-Verlag, 1971).
MG. Wurtele Tellus XIII (3) (1961) 379—391. [15] P.J. Roache, Computational Fluid Dynamics (Hermosa Publishers, Albuquerque, 1972). [16] J.P. Boris and D.L. Book, J. Comput. Phys. 20 (1976) 000—000 [FCT IIII. [171 G. Strang and G. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Englewood Cliffs, New Jersey, i973). [18] S.A. Orszag and M. Israel, Numerical flow simulation by spectral methods, in Proc. of the Symp. Numerical Models of Ocean Circulation, National Academy of Sciences (USGPO, Washington, D.C., 1972). [i9] H.O. Kreiss, A. comparison of numerical methods used in atmospheric and oceanographic applications, in Proc. of the Symp. Numerical Models of Ocean Circulation, National Academy of Sciences (USGPO, Washington, D.C., 1972). See also in the same volume B. Wendroff, Problems of accuracy with conventional finite-difference methods, and G. Fix, A survey of numerical methods for selected problems in continuum mechanics. [20] J.P. Boris, M.J. Fritts and Klaus Ham, Free surface hydrodynamics using a lagrangian triangular mesh, Proc. First Int. COflf. Numerical Ship Hydrodynamics (NBS, Gaitherburg, Md. Oct., i975). [21] M.J. Fritts, SAl Report SAI-76-528-WA, March (i976). [221 B. VanLeer, J. Comput. Phys. 14 (1974) 36 1—370. [23] J.P. Boris and D.L. Book, in: Methods in Computational Physics, Vol. i6 (Academic Press, New York, 1976) ch. 11. [24] J.P. Boris and DL. Book, J. Comput. Phys. 11(1973) 38ff [FCT I]. [25] DL. Book, J.P. Boris and K.H. Ham, J. Comput. Phys. 18(3) (1975) 248—283 [FCT II]. [26] A. Harten, CIMS Report C00-3077-50, June (1974). [27] B. VanLeer, MUSCL, a new approach to numerical gas dynamics, Proc. of the Second European Conf. Computational Physics, Garching, 27—30 Apr. (i976). [28] J.P. Boris, NRL Memorandum Report 2542 (1972). [29] J.P. Boris, J,H, Gardner and S. Zalesak, NRL Memorandum Report 308i (1975). [30] R.A. Chevalier and J.H. Gardner, Astrophys. J. 192 (1974) 427. [31] J.P. Boris, NRL Memorandum Report 3237 (1976) to be published. [32] A. Lapidus, J. Comput. Phys. 8 (1971) i06—i18.