- Email: [email protected]
The stability and accuracy of EPIC algorithms
Computer Physics Communications 44 (1987) 73—82 North-Holland, Amsterdam
73
THE STABILITY AND ACCURACY OF EPIC ALGORITHMS James W. EASTWOOD UKAEA, Cu/ham Laboratory, Abingdon, Oxon 0X14 3DB, UK (Euralom /UKAEA Fusion Association) Received 1 September 1986; in final form 2 February 1987
In the absence of particle discretisation effects (quadrature errors) EPIC algorithms give accurate and unconditionally stable schemes for integrating kinematic fluid and MHD equations. If particle discretisation is introduced as a consistently applied trapezium rule quadrature, then unconditional stability is retained. However, it is shown that not all quadratures lead to stable schemes. General results on linear stability are given. An analysis of dispersive and dissipative effects of finite particle number is presented, and the efficacy of the algorithms is illustrated by kinematics and dynamic test problems. The kinematic test with discrete particles showed unexpected stability and accuracy for large Courant numbers.
1. Infroduction EPIC algorithms [1—3]arise from optimal particle-mesh methods originally developed for collisionless phase fluids. In high dimensional phase space and shocked fluid flows, the Lagrangian sample points (particles) are the natural principal information carriers. However, in low (~ 3) dimensional smooth flows, the mesh can become a more cost effective alternative. Here, we concentrate on the latter case, where particles are treated either as conceptual devices or as quadrature points in evaluating nasty integrals. Which alternative is appropriate depends on the choice of integration method. In this paper we show that Galerkin EPIC schemes remain stable if uniform quadrature is applied consistently to the projection equation. If the quadrature points are shifted with respect to the mesh, then unstable schemes can be found. Unstable schemes may be interpreted as applying too much ‘antidiffusion’ through the mass matrix. The discussions presented are given for the one dimensional case. Generalisation to higher dimensions is straightforward.
2. Zero quadrature errors Let
at
(1)
ax
The Galerkin EPIC scheme advances the discrete approximation f W,,(x)f~,at each timestep, ~t, according to the prescription [f~(x)~(x)
=
dx]f7~t
[f~(x)~(x~)
dx}f~
(2)
where we treat x as a Lagrangian coordinate; x=x(x0, t)x0+~. We restrict our attention to uniform unit node spacing and uniform advection. The basis functions W then have the displacement invariance property =
OO1O-4655/87/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
W0(x —p)
W(x —p)
(3)
74
J. W. Eastwood / Stability and accuracy of EPIC algorithms
and eq. (2) reduces to ~ ~~4) ~,
I(q—p, 0)f~~~~t=I(q_p, flf where the sum over p is implied and J(q-p, fi =JW(x)W(x_
(p-q+fl)
dx (5)
setting f~o ~ and Fourier transforming eq. (4) gives the dispersion relation —
(6)
Jf(k,0) (k, fi e j~o.~t~_
ployed. A displacement ~ = (~, ~) may be represented by a sequence (~,0), (0, i.,) or to second order by a similar time symmetrised sequence of x, y, y, x displacements, etc. A general element connectivity is best tackled using quadrature. 3. Uniform quadrature The simplest example of quadrature, which is readily generalised to higher dimensions, is obtained when integrals in eq. (2) are replaced by sums (trapezoidal rule): ~ W(xi)J~(xi)~jjt±At
and the linear stability condition T(k,E)
(7)
=
[~:
xi)~(xoi)~If.
(12)
I(k, 0) If we let (8)
H(X)=JW(x)W(x_X)dx then its Fourier transform is
(9) ñ(k)= J’~’(k)I2 The set of nodal values {I(p, ~)} can be represented by the generalised function I~(x)=ffl(x—flH(x)DI(k, = iri(-~—)~
fi
~
(10)
*
~i
(11)
i(k,o) thus demonstrating unconditional stability for all choices of W. k~ is the aliased wavenumber, k~= k 2lTn. The integrals appearing in eq. (2) are only easy to evaluate in one dimension. However, if a urnform mesh is employed in two or more dimensions, then fractional timestepping can be em—
S(q—p, O)f,,’~’=S(q,p, E)f,
(13)
yielding the stability criterion
~(k, fi
(14)
~(k,O)
where the transform, ifi and convolution (*) definitions follow those used in ref. [4]. Evaluating the convolution on the r.h.s. of eq. (6) yields i(k, ~)
There are M = 1/~ uniformly spaced points in each element (“particles per cell”) where positions x, = (i + )6 and offset ~ [0, 1]. Following the procedure used for the integral case, we may write eq. (12) as
where SDS. Theorem The projection scheme, eq. (12) is stable for piecewise linear functions W~(x)if = 0 is used in evaluating the mass matrix appearing on the l.h.s. Proof Let H(X) =J~ffl(Mx_)W(x)W(x_X) DH(k)=W(k)
~ m=
dx
~7(k_2.rmM)e_i2~Tm~ -
oo
(15) where ifi and transform definition follow the no-
75
J. W. Eastwood / Stability and accuracy of EPIC algorithms
tation of ref. [4]. Then we may represent the set { S( p, E)} by the generalised function S~(x) = ffl(x
—
flH(x) 00
§(k, ~) =
n
~ =
—00
ñ(k
—
2’rrn)
21T~~.
ei
(16)
For piecewise linear functions, W(x) (= if I x I then 1 Ix I else 0) ~ W(k) = sin2(k/2~r),so H ~ 0 for all k, and therefore —
§(k, E) = ~ 0I~(k— 2Trn) e ~0H(k—2irn)
S(k,0)
~1
q.e.d.
= 0 in eq. (12) is good choice: It gives a stable scheme whose numerical dissipation tends to zero as ~t tends to zero. In practical terms, this implies that one uniformly distributes the ephemeral particles at the new timelevel, and computes where they came from.
4. Linear dispersion The linear dispersion relation for the EPIC scheme is given by the Fourier—Laplace transform of eq. (13) with ~ = constant: e”~’= ~(k, fl/~(k, 0),
Corollary If the mass matrix on the l.h.s. of (12) is lumped, then scheme remains stable.
(17)
where ~(k,
(18)
~)=>S(~±p)e”.
p Proof Lumping the l.h.s. of eq. (12) reduces it to and the stability condition to S(k, ~) I ~ 1. But -
S(k,0)=
2+82
~
1_82
+
~
t~X~k~1
q.e.d
If lumping is used, or if * 0 on the r.h.s. of eq. (12), then we have schemes which will diffuse!, even in the limit ~t 0 (cf. the Lax finite difference scheme [5]). For example, setting ~ = 0, = 0 on the l.h.s. and = 1/2 on the r.h.s. of eq. (12) gives -~
+
(1
~62
)(~+~ 2~ ±i;~)
=~ + (2±82)
-
-
2f~~fp1)
showing quadrature errors introduce an 2/i~t.Whilst thiseffective may be
diffusion coefficient cc 6 it is undesirable in low of value in shocked flows, dissipation smooth flows. The symmetry of this example shows that if ~ 0 on the l.h.s. of eq. (12), then unstable schemes can result. In particular, = 1/2 on the l.h.s. and = 0 on the r.h.s. gives a negative diffusion coefficient cc 8 2/st as ~t-’0. We conclude from the results of this section that if quadrature is used, then uniformly setting
We restrict our attention here to those cases where the same set of quadrature points are used to evaluate both sides of eq. (13). 4.1. One point per element If we take one quadrature point per element, located at the nodes (8 = 1, = 0) then S(x)
(1—Ixi
~
0
IxI~1
otherwise
and for 0 ~ ~ 1
~(k, fl=1—~(1 —cos k)+i~ sink.
(19)
It follows from eq. (19) that this crude quadrature is equivalent to lumping the l.h.s. of eq. (2). If this quadrature were only applied to the l.h.s. of eq. (2), we find the dispersion relation e~”~’= ~
~) = 1—~(1—cosk)±i~sink I( k, 0) 1 — ~sin2k/2
which gives Ie~~~t 2 ~ 1 for all ~, implying unconditional instability. Placing the one quadrature point per element at the centre of the element (6 = 1, = 1/2) gives ~,
.S~(k,~)
.s~’(k,o)
1+ i2~sink 1±cosk
76
J. W. Eastwood
/ Stability and
which again is unstable since Ie~~t I ~ 1. As noted in the previous section, lumping or using quadrature points with offset = 0 gives stable schemes: For the case, with lumping we find
Ie”~’I2=I5~(k,fl12 =
~[(i
±cos
k)2 ± 4~2
sin2kl ~ 1.
4.2. Two points per element Taking quadrature points, 6 = 1/2, = 1/2 in each element gives S(k, ~) piecewise linear in ~, with minima at ~ = 0 and 1/2 and maxima at = 1/4 and 3/4. It therefore follows from eq. (17) that this scheme is unconditionally unstable. Lumping the mass matrix (giving dispersion relation e”~’= S(k, fi) stabiises, since max(~,k) I S(k, ~) I = 1. Similarly, taking quadrature points at 0, 1/2 to evaluate the mass matrix yields scheme stable for all Both of these stabilised schemes are strongly dissipative: In the limit ~t 0, we find that for fixed element sizes, ca, = f(k )/z~t. Taking quadrature points at 0 and 1/2 (6 = 1/2, =0) gives lower dissipation, unconditionally stable scheme. For this case, we find
~.
—*
accuracy of EPIC algorithms
a = 0 and 1 recovers respectively the unstable and stable one point schemes, and a = 1/2 gives the two point trapezium rule discussed above. Theorem For the set of schemes defined by eq. (21) i) a ~ 1/2 gives unconditional stability; ii) the least dissipative unconditionally stable scheme has a = 1/2 (i.e. the trapezoidal rule); lii) the least dissipative scheme stable for Courant number ~ ~ ~maX s~1/2 has a = ~max Proof By symmetry and periodicity, we need only consider the case 0 ~ ~ ~ 1/2. In this case, we find S(1—fl= [1±2C+a(2~—1)]/4
S(s)
r[1_a(2~_1)]/2
S(i-i-fl= [1—2~+a(2~—1)}/4
S(p±fl=0, p is integer *0, ±1. Hence ei~~t= ~(k, ~ ~(k, 0)
—‘+E[
].
2(k/2) —2a+2icot(k/2) a + cot (22)
S(k, fl=~[3—2~±(1+2flcosk±i4~sink] 1/2
Rearranging the stability condition
and hence
IAI2=~(1’~ flu
2
~(k,o)
1 —2 + i4 cot(k/2) e”~’=1±~I 1±2cot2(k/2)
l_
1 ]
implying that Ie”~’I ~ 1, and in the limit zXt 2 ±2cot2(k/2)’
gives
I
cot2(k/2)(E —~
0
(20)
—*f(k), since ~cc z~t. This example is an instance of the trapezium rule applied with two points per element. A more general 8 = 1/2, = 0 two point quadrature is given by the linear combination of the node centred (S 0) and element centred (S1) one point schemes i.e. ~,
S
aS
(21) 0 + (1
—
—
a)
~
a2(1
—
fi.
This inequality is satisfied for all k E [—ii, IT] and ~ E [0, 1/2) provided a ~ 1/2, thus establishing (i). If a < 1/2, then the scheme is stable for Courant numbers ~ ~ Emax = a. In the absence of numerical dissipation I A I 2 = 1, and I A 2 becomes smaller as dissipation increases. Damping increases with wavenumber; ~IA I 2/ak ~ 0 if a ~ 1. The displacement (or Courant number) Eex at which dissipation is maximum is = ex
a(a±cot2O) 2(a2±cot2O)
J. W. Easiwood
/ Stability and
and at that ~ =
4.3. Trapezium rule schemes
AI2= 1 + a21 tan2O Thus, the least dissipative scheme has the smallest a compatible with stability: For unconditional stability, this is a = 1/2, and for stability when a=~m~.
~Emax,1t~tS
Setting a = ~ in eq. (22) shows that the conditionally stable scheme recovers the property of the zero quadrature error case in that w 1 0 as i~t 0. —~
0. 5
0. 0 0.5
Results of evaluating eq. (17) numerically for different quadrature point spacings 8 = 1/M and offsets c are shown in figs. 1—5. Plotted are the real and imaginary parts of the frequency, Wr and ~ versus wavenumber k, for wavenumbers in the prirncpal zone [4]. The differential problem (eq. (1)) has a dispersion relation (0r~t/1TE
=
k/IT,
co~L~t/iT~ =
0.
—~
LU
0.0
77
accuracy of EPIC algorithms
k
1. 0
The deviation of the curves shown in the figures 0~r and ca, give measures from these lines respectively of theforphase errors and numerical damping. Good numerical schemes generally have small damping for wavenumbers with small phase errors, and large damping where phase errors are large [2]. Economy of storage in multidimensional calculations favour schemes with small phase errors up to large wavenumbers, say k = IT/2 (which corresponds to a wavelength of four node spacings). Fig. 1 shows the dependence of the dispersion curves on the number, M, of quadrature points per element. Fig. la shows that once M exceeds 1, phase errors closely follow those of the integral evaluation limit (M = no). Fig. lb shows dissipalion decreasing as M increases.
¶7
1.0
1.0
0. 5 0. 5
0. 0 0.5
-~
¶7
1.0
0. 0 1.0
Fig. 1. Dispersion curves, (a) real part and (b) imaginary part of frequency versus wavenumber for the linear basis function EPIC scheme using trapezium rule integration with =0. Numbers labelling curves are values of M. Wavenumber k ~TT corresponds to awavelength of two elements. The curve labelled cc is the integral case.
Fig. 2. Numerical dissipation for M = 2, = 1/2 evaluation of assignment (r.h.s. of eq. (12)) and evaluation of mass matrix (l.h.s. of eq. (12)) 1) with M 2, = 1/2, 2) with M = 2, 0 and 3) lumped. Curve 4 is the integral case.
J. W. Eastwood
78
/ Stability and
accuracy of EPIC algorithms
1.0
1.0 CA)~
~f
11~
0. 5
0. 5
0. 0
0. 0 0.5
1.0
1.0
0.5
IT
11
Fig. 3. Scaling of numerical dissipation with ~ t for linear EPIC using exact integration (eq. (2)), showing w, —. 0 as —~ 0. The number labelling curves are Courant number, ~.
Fig. 2 compares dissipation for different M = 2 cases. Curve 4 the curve is the integral limit (labelled no in fig. 1). Curve 1 evaluates both sides of eq. (12) with e = 1/2, M = 2, confirming that this choice of quadrature is numerically unstable (w
<0). This may be interpreted as applying too much ‘antidiffusion’ by inverting the mass matrix. Curve 2 shows that using = 0, M = 2 on the 1.h.s. and M= 2, = 1/2 on the r.h.s. of eq. (12) gives stable antidiffusion, but with slightly larger dis1
1.0
Fig. 5. Scaling of dissipation with ~t, showing ~ increasing as decreases if different quadrature is used on each side of eq. (12)[M = 4, c = 0 on l.h.s, M= 4, =1/2 on r.h.s.J
~ t
sipation than obtains when M = 2, = 0 is used consistently (cf. fig. 1). Curve 3 is for no antidiffusion. The mass matrix on the l.h.s. of eq. (12) is lumped to give a unit matrix, and the result, as expected, is strong numerical dissipation. Figs. 3—5 illustrate the change in dissipation with timestep. Fig. 3 shows ~ 0 as zit 0 in the absence of quadrature errors. The numbers labelling the curves are Courant numbers, ~. Figs. 4 and 5 show the same cases using M = 4 and = 0 for mass matrix evaluation and respectively M = 4, = 0 and M = 4, = 1/2 for evaluation of the r.h.s. Fig. 4 shows w~ constant as z~t 0 and fig. 5 shows w~ no as zXt 0. In all three cases, the scaling of ca~is for fixed node spacing; If i~t 0 and z~ 0, then the differential limit can be recovered. The real parts of the frequency corresponding to curves in figs. 2—5 differ little from those given in fig. I a. —b
—*
—~
—~
0. 5
,0
—~
—~
-~
—‘
5. Kinematic test problem
0. 0 0. 5
k 11
I 0
Fig. 4. Scaling of numerical dissipation for linear EPIC using trapezium rule integration with = 0, M= 4, showing w, -. constant as i~t-~ 0.
used Wetoextend compare here the the integral kinematic form of testEPIC problem (eq. (2)) with conventional finite difference schemes [1—3].A Gaussian hump is transported through a unit length periodic box by a velocity field u(x) that provides both compression and expansion.
J. W. Eastwood
/ Stability and accuracy of EPIC algorithms
79
We take u =2 sin 2i~x,N = 100 nodes in a spatial period and a Gaussian density hump placed initially at x = —0.2 near the velocity maximum (p(x, t = 0) = exp[—[x + 1/5]2/02], a = 0.1). The evolutionary equation being modelled is the continuity equation —
ap
a
(23)
We consider three approximations to eq. (23). The integral EPIC scheme using linear trial and basis functions
f
W1(x)p(x, t+~t)dx =
f
Wq(x)p(xo,
t)
the
dx0,
~24)
where x = x0 + The r.h.s. of eq. (24) is a piecewise cubic, and the l.h.s. becomes a tridiagonal matrix of stencil, [1,4,1]/6 multiplying nodal values of p. Given ~, each timestep comprises i) evaluating the r.h.s. of eq. (24), and then inverting the mass matrix to obtain new nodal values
Fig. 6. Initial and final profiles of a Gaussian convected three times round a unit period by a sinusoidal velocity profile using integral linear basis function EPIC scheme, eq. (2). A Courant number of 0.25 was used.
~.
(pt+at; p(l, N)). The second case assumes uniform quadrature at the new timelevel. Setting x(x) = ~(x0), we approximate eq. (24) by ~W~(x1)p(x~,
t+i~t)6
=~Wq(xj)p(xi—~i, t)~8,
ax
(25)
evaluation. In the absence of numerical errors, the two curves would be identical. After the 2000 steps in this test, the errors are smaller than are achieved by state-of-the-art finite difference schemes [2,3]. Fig. 7 illustrates the effect of quadrature errors. As expected, dispersion remains small but dissipation increases as the number of quadrature points per cell is reduced. Figs. 6 and 7 are for a Courant number ~ = 0.25 at the position of maximum velocity. This value was chosen to afford direct comparison with re-
where Wq(Xj) are ‘triangle’ functions [4], p and X are piecewise linear, ax0/ax 1 a~Jax,x, = i8, i integer and 6 = 1/M. The third alternative is to use the trapezium rule on the r.h.s. of eq. (24) at the old timelevel, and treat the Lh.s. as in eq. (25), giving —
~Wq(xj)p(xj, =
t+I~t)8
~ W.~(x~ + E(x01))p(x01,
r)6,
(26)
where ~ is assume linear in x0 over each element. Figs. 6—8 give results of the test problem. Fig. 6 shows the remarkably low dispersion and numerical dissipation achieved by EPIC with integral
Fig. 7. As fig. 6, but using uniform quadrature at the new timelevel, with = 0 and M as labelling the curves.
80
J. W. Eastwood
/ Stability and
accuracy of EPICalgorithms
where v=v(x, t±L~t), v 0=v(x,0, t), x=x0+E and ~ = v01~st.The EPIC projection equation for eq. (28) is
{J(w1~w1+R’~t~-~-~ ~ 8x
=
ax / dx}vi
fWk(x)vodx,
(29)
where v(x) = ~,W1(x)v,, etc. The algorithm for advancing v by one timestep using eq. (29) is Fig. 8. The same problem as fig. 7, but with = 0, M = 4 and Courant numbers as labelling the curves.
a) predict x = x0 + ~ b) find Mk = f Wk(x)vo dx; c) solve Av = M.
sults for other schemes presented in ref. [2]. However, EPIC schemes with = 0 are unconditionally stable, so we are free to choose larger timesteps. Increasing L~t reduces the number of steps to complete the three periods, so reduces computational costs and numerical dissipation. Fig. 8 illustrates this for M = 4, = 0 and three different Courant numbers. Note that these results are for non-uniform velocity u(x). For uniform advection and integer Courant number, numerical dissipation would be zero.
If required, a correction for ~ could be applied, although we find the simple prediction ~ = v0~tis sufficient. Elements of matrix A are given by evaluating the integral on the l.h.s. of eq. (29). A is tridiagonal with stencil ~[1—a,4-i-2a,1—a];
a—6R’L~t.
The conservative EPIC 2/2) scheme obtained the by andis following writing v av/ax as a/ax(v above discretisation steps. In this instance
and 6. Dynamical test problem We have applied three variations of the linear basis function EPIC scheme to Burgers’ equation 8v
at
—
+
av 8x v—
2 1 a 0
(27)
—
a2
—
ax
=
Jw~(x)vodxo.
(30)
R
=
on the unit period interval with initial conditions O = sin 2 ITX. The three cases are i) nonconservative, ii) conservative and iii) nonconservative with quadrature in space (ie with discrete particle effects). In each case the diffusion term on the r.h.s. of eq. (27) is treated fully implicitly. Integrating the nonconservative form, eq. (27) over a small time interval ~t gives O
{f(wkw±R-’z~t.~±f~ aw 1\
R’L~t—~’— = v
2
ax
0,
(28)
Case iii) differs from case i) in that integrals in eq. (29) are replaced by sums with = 0 (cf. eq. (12)). The discretisation may be interpreted as using particles which are uniformly distributed at the new timelevel. Case i) and ii) give almost identical results. For large Reynolds number, R, the sinusoid rapidly becomes a “sawtooth’, then decays on the resistive timescale. Fig. 9 shows results for case i). At R = 100, v are accurately resolved and smooth. At higher R wiggles develop, but unlike many finite difference schemes, which are nonlinearly unstable
J. W. Eastwood
R
R
100
=
/ Stability and
1 000
accuracy of EPIC algorithms
81
at such R due to aliasing errors [6—8], EPIC remains stable. Finite time singularities are not possible in the limit iXt —‘0 for the EPIC schemes discussed here, since in that limit it follows from the projection equation that d/dtfv2dx < 0. The appearance of wiggles is the Gibbs phenomenon, caused by insufficient resolution. Appropriate sub-grid-scale physics models can be used with EPIC (as with other schemes) to remove these wiggles.
R
=
10 000
R
= 100 000
The extra dissipation (cf. fig. 1) in case iii) suppress some of the wiggles. More interestingly, the robustness of the particle advection allows much larger timesteps to be used. Fig. 10 shows the solution for the test problem remains accurate as well as stable for Courant numbers up to four!
7. Final remarks Fig. 9. Sequences showing reduction in amplitude as time increases of solutions to Burgers’ equation. Courant number = 0.5, N = 100 nodes and output every 100~t were used. Reynolds numbers, R, are as indicated.
C
=
0. S
C=2
C=1
C=4
Fig. 10. Sequence showing evolution of solutions to Burgers’ equation using case iii) of section 6, with M =4 particle per element, predictor x = v L~t, R = ion, N = 100 and Courant numbers as shown. Curves are drawn at intervals of 40/C i~t.
The Ephemeral Particle-in-Cell(EPIC) method gives stable, low dissipation and dispersion schemes for integrating the hyperbolic and nearly hyperbolic equations arising, for example, in fluid and particle simulations. This paper has addressed the question of stability and accuracy of EPIC algorithms. The discussion has been presented in one spatial dimension but the results apply equally well in two and three dimensions. In section 2, the stability of EPIC algorithms for any choice of basis/trial function W was established in the integral limit. Quadrature (which may be interpreted as discrete particle effects) can in certain circumstances destroy the stability. However, in section 3 it was shown that unconditionally stable schemes can be constructed for any W and number of particles (quadrature points) per element. Specific cases considered in section 4 identify the unstable cases as those with too much antidiffusion. The least diffusive schemes with uniformly placed particles are those where at least one particle per element coincided with the element nodes. The kinematic and dynamic test problems confirm this, and show that choosing particle positions uniform at the new timelevel is the most effective option. Unconditional stability allows much larger
82
J. W Ea.stwood
/ Stability and
accuracy of EPIC algorithms
timesteps for kinematic equations than is possible with conventional finite differencing in time. One surprising feature of the dynamical test problem was that it also also remained stable. Whether this is generally true or is a special property of the test case has not yet been established.
[4] R.W. Hackney and J.W. Eastwood, Computer Simulation using Particles (McGraw-Hill, New York, 1981). [5]P.J. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, 1976) p. 242. [6] J.W. Eastwood and W. Arter, Phys. Rev. Lett. 57 (1986)
References
[7] W. Arter and J.W. Eastwood, The Effect of Aliasing Error upon Numerical Solutions of the Hydrodynamic and Magnetohydrodynamic Equations, to appear Transp. Tb. Stat. Phys. (1987). [8] J.W. Eastwood and W. Arter, Cuiham Preprint CLM-P-782 (1986) to appear IMA J. Numer. Anal.
[1] J.W. Eastwood, Particle Methods, in: Astrophysical Radiation Hydrodynamics, eds. K.-H. Winkler and M.L. Norman (Reidel, Dordrecht, 1986). [2] J.W. Eastwood and W. Arter, in: Numerical Methods for Fluid Dynamics II, eds. K.W. Morton and M.J. Baines (Oxford Univ. Press, Oxford, 1986) p. 581.
[3]J.W. Eastwood, Comput. Phys. Commun. 43 (1986) 89.
2528.
73
THE STABILITY AND ACCURACY OF EPIC ALGORITHMS James W. EASTWOOD UKAEA, Cu/ham Laboratory, Abingdon, Oxon 0X14 3DB, UK (Euralom /UKAEA Fusion Association) Received 1 September 1986; in final form 2 February 1987
In the absence of particle discretisation effects (quadrature errors) EPIC algorithms give accurate and unconditionally stable schemes for integrating kinematic fluid and MHD equations. If particle discretisation is introduced as a consistently applied trapezium rule quadrature, then unconditional stability is retained. However, it is shown that not all quadratures lead to stable schemes. General results on linear stability are given. An analysis of dispersive and dissipative effects of finite particle number is presented, and the efficacy of the algorithms is illustrated by kinematics and dynamic test problems. The kinematic test with discrete particles showed unexpected stability and accuracy for large Courant numbers.
1. Infroduction EPIC algorithms [1—3]arise from optimal particle-mesh methods originally developed for collisionless phase fluids. In high dimensional phase space and shocked fluid flows, the Lagrangian sample points (particles) are the natural principal information carriers. However, in low (~ 3) dimensional smooth flows, the mesh can become a more cost effective alternative. Here, we concentrate on the latter case, where particles are treated either as conceptual devices or as quadrature points in evaluating nasty integrals. Which alternative is appropriate depends on the choice of integration method. In this paper we show that Galerkin EPIC schemes remain stable if uniform quadrature is applied consistently to the projection equation. If the quadrature points are shifted with respect to the mesh, then unstable schemes can be found. Unstable schemes may be interpreted as applying too much ‘antidiffusion’ through the mass matrix. The discussions presented are given for the one dimensional case. Generalisation to higher dimensions is straightforward.
2. Zero quadrature errors Let
at
(1)
ax
The Galerkin EPIC scheme advances the discrete approximation f W,,(x)f~,at each timestep, ~t, according to the prescription [f~(x)~(x)
=
dx]f7~t
[f~(x)~(x~)
dx}f~
(2)
where we treat x as a Lagrangian coordinate; x=x(x0, t)x0+~. We restrict our attention to uniform unit node spacing and uniform advection. The basis functions W then have the displacement invariance property =
OO1O-4655/87/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
W0(x —p)
W(x —p)
(3)
74
J. W. Eastwood / Stability and accuracy of EPIC algorithms
and eq. (2) reduces to ~ ~~4) ~,
I(q—p, 0)f~~~~t=I(q_p, flf where the sum over p is implied and J(q-p, fi =JW(x)W(x_
(p-q+fl)
dx (5)
setting f~o ~ and Fourier transforming eq. (4) gives the dispersion relation —
(6)
Jf(k,0) (k, fi e j~o.~t~_
ployed. A displacement ~ = (~, ~) may be represented by a sequence (~,0), (0, i.,) or to second order by a similar time symmetrised sequence of x, y, y, x displacements, etc. A general element connectivity is best tackled using quadrature. 3. Uniform quadrature The simplest example of quadrature, which is readily generalised to higher dimensions, is obtained when integrals in eq. (2) are replaced by sums (trapezoidal rule): ~ W(xi)J~(xi)~jjt±At
and the linear stability condition T(k,E)
(7)
=
[~:
xi)~(xoi)~If.
(12)
I(k, 0) If we let (8)
H(X)=JW(x)W(x_X)dx then its Fourier transform is
(9) ñ(k)= J’~’(k)I2 The set of nodal values {I(p, ~)} can be represented by the generalised function I~(x)=ffl(x—flH(x)DI(k, = iri(-~—)~
fi
~
(10)
*
~i
(11)
i(k,o) thus demonstrating unconditional stability for all choices of W. k~ is the aliased wavenumber, k~= k 2lTn. The integrals appearing in eq. (2) are only easy to evaluate in one dimension. However, if a urnform mesh is employed in two or more dimensions, then fractional timestepping can be em—
S(q—p, O)f,,’~’=S(q,p, E)f,
(13)
yielding the stability criterion
~(k, fi
(14)
~(k,O)
where the transform, ifi and convolution (*) definitions follow those used in ref. [4]. Evaluating the convolution on the r.h.s. of eq. (6) yields i(k, ~)
There are M = 1/~ uniformly spaced points in each element (“particles per cell”) where positions x, = (i + )6 and offset ~ [0, 1]. Following the procedure used for the integral case, we may write eq. (12) as
where SDS. Theorem The projection scheme, eq. (12) is stable for piecewise linear functions W~(x)if = 0 is used in evaluating the mass matrix appearing on the l.h.s. Proof Let H(X) =J~ffl(Mx_)W(x)W(x_X) DH(k)=W(k)
~ m=
dx
~7(k_2.rmM)e_i2~Tm~ -
oo
(15) where ifi and transform definition follow the no-
75
J. W. Eastwood / Stability and accuracy of EPIC algorithms
tation of ref. [4]. Then we may represent the set { S( p, E)} by the generalised function S~(x) = ffl(x
—
flH(x) 00
§(k, ~) =
n
~ =
—00
ñ(k
—
2’rrn)
21T~~.
ei
(16)
For piecewise linear functions, W(x) (= if I x I then 1 Ix I else 0) ~ W(k) = sin2(k/2~r),so H ~ 0 for all k, and therefore —
§(k, E) = ~ 0I~(k— 2Trn) e ~0H(k—2irn)
S(k,0)
~1
q.e.d.
= 0 in eq. (12) is good choice: It gives a stable scheme whose numerical dissipation tends to zero as ~t tends to zero. In practical terms, this implies that one uniformly distributes the ephemeral particles at the new timelevel, and computes where they came from.
4. Linear dispersion The linear dispersion relation for the EPIC scheme is given by the Fourier—Laplace transform of eq. (13) with ~ = constant: e”~’= ~(k, fl/~(k, 0),
Corollary If the mass matrix on the l.h.s. of (12) is lumped, then scheme remains stable.
(17)
where ~(k,
(18)
~)=>S(~±p)e”.
p Proof Lumping the l.h.s. of eq. (12) reduces it to and the stability condition to S(k, ~) I ~ 1. But -
S(k,0)=
2+82
~
1_82
+
~
t~X~k~1
q.e.d
If lumping is used, or if * 0 on the r.h.s. of eq. (12), then we have schemes which will diffuse!, even in the limit ~t 0 (cf. the Lax finite difference scheme [5]). For example, setting ~ = 0, = 0 on the l.h.s. and = 1/2 on the r.h.s. of eq. (12) gives -~
+
(1
~62
)(~+~ 2~ ±i;~)
=~ + (2±82)
-
-
2f~~fp1)
showing quadrature errors introduce an 2/i~t.Whilst thiseffective may be
diffusion coefficient cc 6 it is undesirable in low of value in shocked flows, dissipation smooth flows. The symmetry of this example shows that if ~ 0 on the l.h.s. of eq. (12), then unstable schemes can result. In particular, = 1/2 on the l.h.s. and = 0 on the r.h.s. gives a negative diffusion coefficient cc 8 2/st as ~t-’0. We conclude from the results of this section that if quadrature is used, then uniformly setting
We restrict our attention here to those cases where the same set of quadrature points are used to evaluate both sides of eq. (13). 4.1. One point per element If we take one quadrature point per element, located at the nodes (8 = 1, = 0) then S(x)
(1—Ixi
~
0
IxI~1
otherwise
and for 0 ~ ~ 1
~(k, fl=1—~(1 —cos k)+i~ sink.
(19)
It follows from eq. (19) that this crude quadrature is equivalent to lumping the l.h.s. of eq. (2). If this quadrature were only applied to the l.h.s. of eq. (2), we find the dispersion relation e~”~’= ~
~) = 1—~(1—cosk)±i~sink I( k, 0) 1 — ~sin2k/2
which gives Ie~~~t 2 ~ 1 for all ~, implying unconditional instability. Placing the one quadrature point per element at the centre of the element (6 = 1, = 1/2) gives ~,
.S~(k,~)
.s~’(k,o)
1+ i2~sink 1±cosk
76
J. W. Eastwood
/ Stability and
which again is unstable since Ie~~t I ~ 1. As noted in the previous section, lumping or using quadrature points with offset = 0 gives stable schemes: For the case, with lumping we find
Ie”~’I2=I5~(k,fl12 =
~[(i
±cos
k)2 ± 4~2
sin2kl ~ 1.
4.2. Two points per element Taking quadrature points, 6 = 1/2, = 1/2 in each element gives S(k, ~) piecewise linear in ~, with minima at ~ = 0 and 1/2 and maxima at = 1/4 and 3/4. It therefore follows from eq. (17) that this scheme is unconditionally unstable. Lumping the mass matrix (giving dispersion relation e”~’= S(k, fi) stabiises, since max(~,k) I S(k, ~) I = 1. Similarly, taking quadrature points at 0, 1/2 to evaluate the mass matrix yields scheme stable for all Both of these stabilised schemes are strongly dissipative: In the limit ~t 0, we find that for fixed element sizes, ca, = f(k )/z~t. Taking quadrature points at 0 and 1/2 (6 = 1/2, =0) gives lower dissipation, unconditionally stable scheme. For this case, we find
~.
—*
accuracy of EPIC algorithms
a = 0 and 1 recovers respectively the unstable and stable one point schemes, and a = 1/2 gives the two point trapezium rule discussed above. Theorem For the set of schemes defined by eq. (21) i) a ~ 1/2 gives unconditional stability; ii) the least dissipative unconditionally stable scheme has a = 1/2 (i.e. the trapezoidal rule); lii) the least dissipative scheme stable for Courant number ~ ~ ~maX s~1/2 has a = ~max Proof By symmetry and periodicity, we need only consider the case 0 ~ ~ ~ 1/2. In this case, we find S(1—fl= [1±2C+a(2~—1)]/4
S(s)
r[1_a(2~_1)]/2
S(i-i-fl= [1—2~+a(2~—1)}/4
S(p±fl=0, p is integer *0, ±1. Hence ei~~t= ~(k, ~ ~(k, 0)
—‘+E[
].
2(k/2) —2a+2icot(k/2) a + cot (22)
S(k, fl=~[3—2~±(1+2flcosk±i4~sink] 1/2
Rearranging the stability condition
and hence
IAI2=~(1’~ flu
2
~(k,o)
1 —2 + i4 cot(k/2) e”~’=1±~I 1±2cot2(k/2)
l_
1 ]
implying that Ie”~’I ~ 1, and in the limit zXt 2 ±2cot2(k/2)’
gives
I
cot2(k/2)(E —~
0
(20)
—*f(k), since ~cc z~t. This example is an instance of the trapezium rule applied with two points per element. A more general 8 = 1/2, = 0 two point quadrature is given by the linear combination of the node centred (S 0) and element centred (S1) one point schemes i.e. ~,
S
aS
(21) 0 + (1
—
—
a)
~
a2(1
—
fi.
This inequality is satisfied for all k E [—ii, IT] and ~ E [0, 1/2) provided a ~ 1/2, thus establishing (i). If a < 1/2, then the scheme is stable for Courant numbers ~ ~ Emax = a. In the absence of numerical dissipation I A I 2 = 1, and I A 2 becomes smaller as dissipation increases. Damping increases with wavenumber; ~IA I 2/ak ~ 0 if a ~ 1. The displacement (or Courant number) Eex at which dissipation is maximum is = ex
a(a±cot2O) 2(a2±cot2O)
J. W. Easiwood
/ Stability and
and at that ~ =
4.3. Trapezium rule schemes
AI2= 1 + a21 tan2O Thus, the least dissipative scheme has the smallest a compatible with stability: For unconditional stability, this is a = 1/2, and for stability when a=~m~.
~Emax,1t~tS
Setting a = ~ in eq. (22) shows that the conditionally stable scheme recovers the property of the zero quadrature error case in that w 1 0 as i~t 0. —~
0. 5
0. 0 0.5
Results of evaluating eq. (17) numerically for different quadrature point spacings 8 = 1/M and offsets c are shown in figs. 1—5. Plotted are the real and imaginary parts of the frequency, Wr and ~ versus wavenumber k, for wavenumbers in the prirncpal zone [4]. The differential problem (eq. (1)) has a dispersion relation (0r~t/1TE
=
k/IT,
co~L~t/iT~ =
0.
—~
LU
0.0
77
accuracy of EPIC algorithms
k
1. 0
The deviation of the curves shown in the figures 0~r and ca, give measures from these lines respectively of theforphase errors and numerical damping. Good numerical schemes generally have small damping for wavenumbers with small phase errors, and large damping where phase errors are large [2]. Economy of storage in multidimensional calculations favour schemes with small phase errors up to large wavenumbers, say k = IT/2 (which corresponds to a wavelength of four node spacings). Fig. 1 shows the dependence of the dispersion curves on the number, M, of quadrature points per element. Fig. la shows that once M exceeds 1, phase errors closely follow those of the integral evaluation limit (M = no). Fig. lb shows dissipalion decreasing as M increases.
¶7
1.0
1.0
0. 5 0. 5
0. 0 0.5
-~
¶7
1.0
0. 0 1.0
Fig. 1. Dispersion curves, (a) real part and (b) imaginary part of frequency versus wavenumber for the linear basis function EPIC scheme using trapezium rule integration with =0. Numbers labelling curves are values of M. Wavenumber k ~TT corresponds to awavelength of two elements. The curve labelled cc is the integral case.
Fig. 2. Numerical dissipation for M = 2, = 1/2 evaluation of assignment (r.h.s. of eq. (12)) and evaluation of mass matrix (l.h.s. of eq. (12)) 1) with M 2, = 1/2, 2) with M = 2, 0 and 3) lumped. Curve 4 is the integral case.
J. W. Eastwood
78
/ Stability and
accuracy of EPIC algorithms
1.0
1.0 CA)~
~f
11~
0. 5
0. 5
0. 0
0. 0 0.5
1.0
1.0
0.5
IT
11
Fig. 3. Scaling of numerical dissipation with ~ t for linear EPIC using exact integration (eq. (2)), showing w, —. 0 as —~ 0. The number labelling curves are Courant number, ~.
Fig. 2 compares dissipation for different M = 2 cases. Curve 4 the curve is the integral limit (labelled no in fig. 1). Curve 1 evaluates both sides of eq. (12) with e = 1/2, M = 2, confirming that this choice of quadrature is numerically unstable (w
<0). This may be interpreted as applying too much ‘antidiffusion’ by inverting the mass matrix. Curve 2 shows that using = 0, M = 2 on the 1.h.s. and M= 2, = 1/2 on the r.h.s. of eq. (12) gives stable antidiffusion, but with slightly larger dis1
1.0
Fig. 5. Scaling of dissipation with ~t, showing ~ increasing as decreases if different quadrature is used on each side of eq. (12)[M = 4, c = 0 on l.h.s, M= 4, =1/2 on r.h.s.J
~ t
sipation than obtains when M = 2, = 0 is used consistently (cf. fig. 1). Curve 3 is for no antidiffusion. The mass matrix on the l.h.s. of eq. (12) is lumped to give a unit matrix, and the result, as expected, is strong numerical dissipation. Figs. 3—5 illustrate the change in dissipation with timestep. Fig. 3 shows ~ 0 as zit 0 in the absence of quadrature errors. The numbers labelling the curves are Courant numbers, ~. Figs. 4 and 5 show the same cases using M = 4 and = 0 for mass matrix evaluation and respectively M = 4, = 0 and M = 4, = 1/2 for evaluation of the r.h.s. Fig. 4 shows w~ constant as z~t 0 and fig. 5 shows w~ no as zXt 0. In all three cases, the scaling of ca~is for fixed node spacing; If i~t 0 and z~ 0, then the differential limit can be recovered. The real parts of the frequency corresponding to curves in figs. 2—5 differ little from those given in fig. I a. —b
—*
—~
—~
0. 5
,0
—~
—~
-~
—‘
5. Kinematic test problem
0. 0 0. 5
k 11
I 0
Fig. 4. Scaling of numerical dissipation for linear EPIC using trapezium rule integration with = 0, M= 4, showing w, -. constant as i~t-~ 0.
used Wetoextend compare here the the integral kinematic form of testEPIC problem (eq. (2)) with conventional finite difference schemes [1—3].A Gaussian hump is transported through a unit length periodic box by a velocity field u(x) that provides both compression and expansion.
J. W. Eastwood
/ Stability and accuracy of EPIC algorithms
79
We take u =2 sin 2i~x,N = 100 nodes in a spatial period and a Gaussian density hump placed initially at x = —0.2 near the velocity maximum (p(x, t = 0) = exp[—[x + 1/5]2/02], a = 0.1). The evolutionary equation being modelled is the continuity equation —
ap
a
(23)
We consider three approximations to eq. (23). The integral EPIC scheme using linear trial and basis functions
f
W1(x)p(x, t+~t)dx =
f
Wq(x)p(xo,
t)
the
dx0,
~24)
where x = x0 + The r.h.s. of eq. (24) is a piecewise cubic, and the l.h.s. becomes a tridiagonal matrix of stencil, [1,4,1]/6 multiplying nodal values of p. Given ~, each timestep comprises i) evaluating the r.h.s. of eq. (24), and then inverting the mass matrix to obtain new nodal values
Fig. 6. Initial and final profiles of a Gaussian convected three times round a unit period by a sinusoidal velocity profile using integral linear basis function EPIC scheme, eq. (2). A Courant number of 0.25 was used.
~.
(pt+at; p(l, N)). The second case assumes uniform quadrature at the new timelevel. Setting x(x) = ~(x0), we approximate eq. (24) by ~W~(x1)p(x~,
t+i~t)6
=~Wq(xj)p(xi—~i, t)~8,
ax
(25)
evaluation. In the absence of numerical errors, the two curves would be identical. After the 2000 steps in this test, the errors are smaller than are achieved by state-of-the-art finite difference schemes [2,3]. Fig. 7 illustrates the effect of quadrature errors. As expected, dispersion remains small but dissipation increases as the number of quadrature points per cell is reduced. Figs. 6 and 7 are for a Courant number ~ = 0.25 at the position of maximum velocity. This value was chosen to afford direct comparison with re-
where Wq(Xj) are ‘triangle’ functions [4], p and X are piecewise linear, ax0/ax 1 a~Jax,x, = i8, i integer and 6 = 1/M. The third alternative is to use the trapezium rule on the r.h.s. of eq. (24) at the old timelevel, and treat the Lh.s. as in eq. (25), giving —
~Wq(xj)p(xj, =
t+I~t)8
~ W.~(x~ + E(x01))p(x01,
r)6,
(26)
where ~ is assume linear in x0 over each element. Figs. 6—8 give results of the test problem. Fig. 6 shows the remarkably low dispersion and numerical dissipation achieved by EPIC with integral
Fig. 7. As fig. 6, but using uniform quadrature at the new timelevel, with = 0 and M as labelling the curves.
80
J. W. Eastwood
/ Stability and
accuracy of EPICalgorithms
where v=v(x, t±L~t), v 0=v(x,0, t), x=x0+E and ~ = v01~st.The EPIC projection equation for eq. (28) is
{J(w1~w1+R’~t~-~-~ ~ 8x
=
ax / dx}vi
fWk(x)vodx,
(29)
where v(x) = ~,W1(x)v,, etc. The algorithm for advancing v by one timestep using eq. (29) is Fig. 8. The same problem as fig. 7, but with = 0, M = 4 and Courant numbers as labelling the curves.
a) predict x = x0 + ~ b) find Mk = f Wk(x)vo dx; c) solve Av = M.
sults for other schemes presented in ref. [2]. However, EPIC schemes with = 0 are unconditionally stable, so we are free to choose larger timesteps. Increasing L~t reduces the number of steps to complete the three periods, so reduces computational costs and numerical dissipation. Fig. 8 illustrates this for M = 4, = 0 and three different Courant numbers. Note that these results are for non-uniform velocity u(x). For uniform advection and integer Courant number, numerical dissipation would be zero.
If required, a correction for ~ could be applied, although we find the simple prediction ~ = v0~tis sufficient. Elements of matrix A are given by evaluating the integral on the l.h.s. of eq. (29). A is tridiagonal with stencil ~[1—a,4-i-2a,1—a];
a—6R’L~t.
The conservative EPIC 2/2) scheme obtained the by andis following writing v av/ax as a/ax(v above discretisation steps. In this instance
and 6. Dynamical test problem We have applied three variations of the linear basis function EPIC scheme to Burgers’ equation 8v
at
—
+
av 8x v—
2 1 a 0
(27)
—
a2
—
ax
=
Jw~(x)vodxo.
(30)
R
=
on the unit period interval with initial conditions O = sin 2 ITX. The three cases are i) nonconservative, ii) conservative and iii) nonconservative with quadrature in space (ie with discrete particle effects). In each case the diffusion term on the r.h.s. of eq. (27) is treated fully implicitly. Integrating the nonconservative form, eq. (27) over a small time interval ~t gives O
{f(wkw±R-’z~t.~±f~ aw 1\
R’L~t—~’— = v
2
ax
0,
(28)
Case iii) differs from case i) in that integrals in eq. (29) are replaced by sums with = 0 (cf. eq. (12)). The discretisation may be interpreted as using particles which are uniformly distributed at the new timelevel. Case i) and ii) give almost identical results. For large Reynolds number, R, the sinusoid rapidly becomes a “sawtooth’, then decays on the resistive timescale. Fig. 9 shows results for case i). At R = 100, v are accurately resolved and smooth. At higher R wiggles develop, but unlike many finite difference schemes, which are nonlinearly unstable
J. W. Eastwood
R
R
100
=
/ Stability and
1 000
accuracy of EPIC algorithms
81
at such R due to aliasing errors [6—8], EPIC remains stable. Finite time singularities are not possible in the limit iXt —‘0 for the EPIC schemes discussed here, since in that limit it follows from the projection equation that d/dtfv2dx < 0. The appearance of wiggles is the Gibbs phenomenon, caused by insufficient resolution. Appropriate sub-grid-scale physics models can be used with EPIC (as with other schemes) to remove these wiggles.
R
=
10 000
R
= 100 000
The extra dissipation (cf. fig. 1) in case iii) suppress some of the wiggles. More interestingly, the robustness of the particle advection allows much larger timesteps to be used. Fig. 10 shows the solution for the test problem remains accurate as well as stable for Courant numbers up to four!
7. Final remarks Fig. 9. Sequences showing reduction in amplitude as time increases of solutions to Burgers’ equation. Courant number = 0.5, N = 100 nodes and output every 100~t were used. Reynolds numbers, R, are as indicated.
C
=
0. S
C=2
C=1
C=4
Fig. 10. Sequence showing evolution of solutions to Burgers’ equation using case iii) of section 6, with M =4 particle per element, predictor x = v L~t, R = ion, N = 100 and Courant numbers as shown. Curves are drawn at intervals of 40/C i~t.
The Ephemeral Particle-in-Cell(EPIC) method gives stable, low dissipation and dispersion schemes for integrating the hyperbolic and nearly hyperbolic equations arising, for example, in fluid and particle simulations. This paper has addressed the question of stability and accuracy of EPIC algorithms. The discussion has been presented in one spatial dimension but the results apply equally well in two and three dimensions. In section 2, the stability of EPIC algorithms for any choice of basis/trial function W was established in the integral limit. Quadrature (which may be interpreted as discrete particle effects) can in certain circumstances destroy the stability. However, in section 3 it was shown that unconditionally stable schemes can be constructed for any W and number of particles (quadrature points) per element. Specific cases considered in section 4 identify the unstable cases as those with too much antidiffusion. The least diffusive schemes with uniformly placed particles are those where at least one particle per element coincided with the element nodes. The kinematic and dynamic test problems confirm this, and show that choosing particle positions uniform at the new timelevel is the most effective option. Unconditional stability allows much larger
82
J. W Ea.stwood
/ Stability and
accuracy of EPIC algorithms
timesteps for kinematic equations than is possible with conventional finite differencing in time. One surprising feature of the dynamical test problem was that it also also remained stable. Whether this is generally true or is a special property of the test case has not yet been established.
[4] R.W. Hackney and J.W. Eastwood, Computer Simulation using Particles (McGraw-Hill, New York, 1981). [5]P.J. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, 1976) p. 242. [6] J.W. Eastwood and W. Arter, Phys. Rev. Lett. 57 (1986)
References
[7] W. Arter and J.W. Eastwood, The Effect of Aliasing Error upon Numerical Solutions of the Hydrodynamic and Magnetohydrodynamic Equations, to appear Transp. Tb. Stat. Phys. (1987). [8] J.W. Eastwood and W. Arter, Cuiham Preprint CLM-P-782 (1986) to appear IMA J. Numer. Anal.
[1] J.W. Eastwood, Particle Methods, in: Astrophysical Radiation Hydrodynamics, eds. K.-H. Winkler and M.L. Norman (Reidel, Dordrecht, 1986). [2] J.W. Eastwood and W. Arter, in: Numerical Methods for Fluid Dynamics II, eds. K.W. Morton and M.J. Baines (Oxford Univ. Press, Oxford, 1986) p. 581.
[3]J.W. Eastwood, Comput. Phys. Commun. 43 (1986) 89.
2528.