- Email: [email protected]
Closed linear one-cell functional spatial approximations
Ann. nucl. Energy, Vol. 14, No. 4, pp. 177-191, 1987
0306-4549/87 $3.00+0.00 Copyright © 1987 Pergamon Journals Ltd
Printed in Great Britain. All rights reserved
CLOSED L I N E A R ONE-CELL F U N C T I O N A L SPATIAL APPROXIMATIONS P. NELSON" Department of Applied Mathematics, California Institute of Technology, Pasadena, CA 91125, U.S.A.
(Received 29 July 1986) Abstract--Closed linear one-cell functional (closed LOF) methods are introduced as a class of discrete spatial approximations to the discrete-ordinate equations for anisotropically scattering monoenergetic transport in slab geometry. It is shown that most spatial approximations currently used or suggested for use, for this class of problems, are specific realizations of closed LOF methods. An appropriate concept of consistency is introduced for closed LOF methods. It is shown that a consistent closed LOF method is asymptotically a finite-difference method, in the sense that the approximating equations can be expressed in terms of approximations to the cell-edge fluxes, for sufficiently fine meshes. Further questions related to closed LOF methods are suggested for future study.
l. INTRODUCTION
We write the discrete-ordinate equations for anisotropically scattering monoenergetic transport in slab geometry as / z ~ x ' +a~bi = S~, i = 1. . . . . N,
(la)
where N
Si(x) : =
~ wrkir(X)Or(x)+qi(x).
i'=l
(lb)
Here the p~s [ - 1, 1] ~ {0} are the quadrature points, with associated quadrature weights {w~}, a = a(x) is the scattering kernel and q is the given source function. We shall assume that a and the kit are piecewise constant. Further we shall consider equations (la,b) subject only to boundary conditions representing specified angular fluxes incident upon the slab faces, say qJi(0) = f
for/~i > 0,
(2a)
~bi(a) = f
for #~ < 0,
(2b)
where a > 0. Suppose momentarily that a and the {wi}, {k~r}, {qi} and {f} are all nonnegative, as is the case in most instances of interest. In this situation, if there exists a nonnegative solution of equations (la,b) and (2a,b), as is known (Nelson, 1973) to be the case if the underlying exact transport problem is subcritical and the quadrature rule is sufficiently fine, then that solution must be [see for example, Nelson (1971)] the limit, as n -~ ~ , of the {~/,~n)}, where the latter are the inner iteration approximations defined by 1 ~ i ~ - [ - o - ~ 1 } n+l)
:
S~ n),
i = 1. . . . . N,
(3a)
and N
S} ")(x) = ~ wrk,r (x)~P},")(x) + q, (x).
(3b)
i'=l
t Permanent address: Department of Mathematics, Institute for Numerical Transport Theory, Texas Tech University, P.O. Box 4319, Lubbock, TX 79409, U.S.A. 177
178
P. NELSON h
X I
Xr
Fig. 1. Notation for a typical mesh cell, C.
Here the initial approximations {~b}°)} or {S}°)} are arbitrary, modulo suitable technical conditions, and the approximate angular fluxes {~k}")} also satisfy the boundary conditions (2a,b). The angular coupling inherent in equations (la,b) tends to recede to the background within the inner iteration process, in the sense that at each iterative stage (i.e. each fixed value of n) the approximate angular fluxes are determined by N monodirectional monoenergetic transport equations,
dq,
# dxx + trio = S,
(4)
where S is known. Further, each such ff is specified at only a single point, which is to say that, in mathematical terminology, equation (4) is to be solved subject to initial conditions. The process of solving such an initialvalue problem for equation (4) often is termed ray tracing. Because of the facts described in the preceding paragraph, spatial approximations for equations (1 a,b) and (2a,b) often are described explicitly in their application to the ray-tracing problem only. In Nelson and Zelazny (1986) the class of linear one-cell functional (LOF) methods for ray tracing was introduced, and the consistency, convergence and stability properties of such methods were studied. More specifically, it was shown that, in application to the ray-tracing problem, a consistent LOF method is both convergent and (zero) stable. The LOF methods seem to include very nearly all methods that have been used, or suggested for use, for raytracing calculations as they occur for a discrete spatial aproximation of equations (1) and (2) within the framework of an overall inner iteration scheme. Implementation of an LOF method for equation (4) requires, as detailed in Section 2, the values of certain quantities, the so-called basic linear functionals of the particular method, corresponding to the particular source(s) under consideration. Nelson and Zelazny (1986) assumed that any information regarding S that was required to construct these values was available as needed. In point of fact, within the iterative scheme (3) the relevant source depends, via equation (3b), upon the approximate angular flux generated within the preceding iteration. Therefore, the information required, at each iterative stage, to construct, possibly approximately, the required values of the basic linear functionals must have been generated (and stored) during the preceding iteration. In the present paper we show that this requirement can be met by appending certain closure approximations to an LOF method, and further we make an initial study of the resulting closed LOF methods. Closed LOF methods are defined in Section 2. In Section 3 we show how a variety of specific spatial approximations fall within the class of closed LOF methods. Section 4 is devoted to considerations relating to extension of the concept of consistency, as defined by Nelson and Zelazny (1986) for LOF methods, to encompass the closure approximations, and thus closed LOF methods. The basic result of Section 5 is a theorem to the effect that, under suitable assumptions on the associated basic linear functionals, a consistent closed LOF method asymptotically is a finite-difference method ; that is, we show that, for suitable consistent closed LOF methods and sufficiently fine spatial meshes, the required values of the basic linear functionals can be expressed in terms of the approximate values of the angular flux at the cell edges (mesh points). Some possible topics for further study, relating to closed LOF methods, are briefly outlined in the concluding Section 6.
2. CLOSED LOF METHODS
We begin by defining an LOF method for the ray-tracing problem for equation (4). Consider a typical spatial cell, C = [x~,Xr], as in Fig. 1. For/~ > 0 and given ~k(xt), possibly approximate, an LOF method for this ray-tracing problem is an approximation of the form
I/l(Xr) ~-~ ~(Xl, h ; ~b(xt)) : = z(a(xl+), h, #)~(x,) + ~ bm(a(xt+),h, #)Lm(S [C; xl, h), m=l
(5)
Closed linear one-cell functional methods
179
where dependence upon the data of definition (5) (i.e./~, a and S) has not been notationally included in the argument list for ~, because we usually consider these arbitrary but fixed, but all dependencies have been explicitly indicated on the r.h.s, of expression (5), for clarity. Here z and the bk are ordinary point functions, but for each fixed h > 0 the L,, (" ; xz, h) are linear functionals that are termed the basic linear functionals of the particular method. The notation a(x~+) indicates the right-hand limit of a(x) at the left cell edge. If the mesh is arranged so that material interfaces are cell edges, as is invariably done in practice and we henceforth assume for the remainder of this paper, then this will be the value of tr throughout the cell ; unless a material interface is located at xt, this also will be the same as the left-hand limit, a ( x r ) . The notation 'S I C ' ( = S restricted to C) indicates that the L,(" ; xt, h) are allowed to depend only upon the behavior of S within the cell C. F o r / t < 0 an approximation similar to formula (5) is assumed to hold, with the roles of xt and x~ reversed. By comparing formula (5) with the corresponding exact solution of equation (4). subject to given ~,(xt), it is clear that z(a(x~+), h, /~) should be an approximation to the exponential attenuation factor across C, in direction/~ and that the second term on the r.h.s, of equation (5) should approximate the uncollided flux, emergent from C, in direction #, stemming from sources within C. These considerations will be made precise in our discussion of consistency, in Section 4. We remark that one-step finite-difference methods of the type frequently used for general-purpose computing are instances of LOF methods ; however, in such methods the Lk(SI C; x~, h) invariably are taken as point evaluations of S, whereas other types of functionals often are employed within transport theory, as shown by the examples of the next section. In application to ray tracing for equation (3a), an LOF method has the form
~) + ~" ~0<.+ i.r 1) = z(a(xt+),h, I,,~ i"~,l,(n+ ~ bm(a(xl+),h,#,)L,,,(S}")lC;xt, h). lWi,l
(6a)
m--I
By virtue of equation (3b), and constancy of the kii. within C, the values of the basic linear functionals that appear in equation (6a) can be expressed as N
L,,(S} ") I C;x,,h) = ~ wrk,,,L,,(~b},") I C; x,,h)+ Lm(qt I C;xz, h).
(6b)
i'=1
As the qi are given (usually piecewise constant), we assume that the Lm(qi[ C; x~, h) can be computed satisfactorily, at least approximately. Thus, in order to close the iterative loop, in the sense of permitting the approximate cell-edge fluxes to be computed via equation (6a), it remains only to devise some satisfactory method for approximating the Lm(~l},n) l C ; Xl, h). As we only know that the ~,}")satisfy equation (3a), with n replaced by n - 1, the task just described essentially is that of approximating Lm(¢l C), for any given S, where ff and S are related by equation (4), subject to known ~b(xt) (~9(xr)) for # > 0 (/~ < 0). The form (5) of an LOF approximation suggests the use of closure approximations :
Lm(~l C;xt, h) ~ £,n(SI C;xt, h;~O(xt)) : = "Cm(ff(xt+),h,ll)~(xl)q- ~ bmm,(a(Xz+),h,#)L,n,(S ] C;xt, h), (7) m'--I
where it is to be understood that ~k and S are related as described above. The dependence upon tr and # is not explicitly indicated in the argument list for the £m, as for ff alone, but the dependence upon S is indicated, as consideration of these approximations for varying S is central to the use of closure approximations. A closed LOF method for equations (1) and (2) is simply an LOF method along with closure approximations of the form (7) f o r / t > 0, and corresponding approximations, but with xz replaced by xr and x~+ by xr , for # < 0. In Section 5 we show, in detail, how a closed LOF method provides an effective iteration scheme, analogous to scheme (3), for generating approximations to the solution of the basic underlying problem (1, 2). However, we take as our next task to provide a demonstration that a variety of well-known spatial approximations for problem (1, 2) are specific instances of closed LOF methods.
3. EXAMPLES OF CLOSED LOF METHODS
Throughout this section we denote the approximations to the solution of equation (4) at the left-hand and right-hand cell edge by q~l and ~ , respectively. Thus, for each LOF method considered,
180
P, NELSON
~L := qT(x,,h; ~7,)
(8)
in the notation of definition (5).
Example A ; the diamond-difference method This method traditionally is considered to be defined by the particle-balance equation, h (qTr- @) + ~r~ = S
(9)
q7 = (q7r + @)/2,
(10)
and the supplementary equation
where q7 and ~qare interpreted as approximations to the averages over C of ~Oand S , repectively. (Here, as in the remainder of this section, a denotes the constant value of the cross section within the particular cell under consideration.) These can be solved for qTr,in terms of ~qand @, to yield
1 - ah/2# ~
2h
This clearly has the form of an LOF method, in the notation of formulas (5) and (8), provided we take p = 1 and
L,(S[ C;x,,h) : = h
S ( x ) d x = ~q
(11)
as the sole basic linear functional. If we further take ~ as £ ~, the required approximation to L ~(~ [ C; xt, h), then there results
£ , ( S I C ; x , , h ; ~ , ) = 1 [1 + (2# - ah)/(2# + crh)]@ + 2# +h ah S, _ which clearly has the form (7) of a closure approximation, as required to display the diamond-difference method as a closed L O F method.
Example B; the step-characteristic method (Lathrop, 1969) This method is defined by replacing S in equation (4) by its cell average, formula (11), and then solving the result, subject to presumed known qT(x,) = ~, (for # > 0), for ~ = ~(x), an approximate angular flux defined throughout the cell. This yields 1
~r = exp ( - ah/#)~, + ~ [1 -
exp
( -
ah/#)]~,
which clearly is of the form (5), as required to comprise an L O F method. If, further, the required approximation to L ~(~k I C; xl, h) is taken as L 1(~ [ C; xl, h) -- £ ~, then we find
£,(Sl C ; x , , h ; ~ , ) = ~h which is precisely of the form (7), as required to constitute a closure approximation. Thus the step-characteristic method also can be regarded as a closed L O F method.
Example C ; moment-characteristic methods o f Vaidyanathan (1977a, b, 1979 ; Lee and Vaidyanathan, 1980) These are a generalization of the step-characteristic method. They are defined by again replacing S in equation (4) by an approximation, and solving the result, subject to known qT(x3 = ~t (# > 0), for q7 = qT(x), which is to be interpreted as an approximate angular flux defined throughout the cell. The approximate source used is
Closed linear one-cell functional methods p
181
1
g(x) := y~ Smpm(x), m=0
where the p,.(x) are the cell-based Legendre polynomials of Larsen and Miller (1980) and the p basic linear functionals of the method are defined by
am
=
Lm+ ,(SI C;x~, h) . - 2 mh+ l fc S(x)p,,(x)dx,
m = 0. . . . , p - 1 .
(12)
The resulting value of ~r is ? P if, = ff(xr) = exp ( - ah/#)@ + E S,,_, | exp [-- a(x, -- X')/ItlPm_l (X') dx', m=0
(13)
dC
which clearly has the required form (5) of an LOF method.By using fundamental properties of the Legendre polynomials, a simple recursive formula can be given for the integrals in equation (13), but that is not of real concern here. For these methods the required closure approximations to Lm(~OI C; h) are simply taken as Lm(i~ [ C; h). The latter are perhaps best computed by multiplying equation (4), with S replaced by S and ~ by qT, by Pm 1(x) and integrating the result over C. After some manipulations, there results [(m- 2)/2]
# --(2m--1)[gTr+(--1)m~7,]+2(2m--1) £m := Lm(qTIC ; x , , h ) = ~h
Y~
}
Lm 1 2,(OIC;x,,h)
l=0
+
1
cr
Lm(SlC;x,,h),
m = l . . . . . p,
(14)
where the arguments of £,, are as before and [x] denotes the greatest integer not exceeding x. Upon using equation (13), it is clear that formula (14) can be employed recursively to yield closure approximations of the required form (7). The detailed form of the resulting closure approximations is not essential to our purposes. Clearly the step-characteristic method is the case p = 1 of this class of closed LOF methods.
Example D ; continuous methods (Neta and Victory, 1982, 1983 ; Victory and Ganguly, 1986) In these methods the basic linear functionals again are given by definition (12), and the moment equations (14) again play a major role. The difference from the moment-characteristic methods stems from the fact that now the approximate angular flux within the cell is taken as a polynomial of degree p, p+l
~ ( x ) = ~ Zm(tfflf;x,,h)pm l(X). m-I
This yields p+ I
qT~ : = qT(x~) = ~ Lm(~lC;x,,h),
(15)
m=l
and the result of substituting this into the moment equations (14) can be regarded as a system o f p equations in the p 4-1 unknowns Lm(~[ C; xt, h), with known terms (i.e. 'right-hand side') depending linearly upon ~t and the Lm(S] C; xt, h). This is made into a determined system of equations by requiring that ~(xt) = ~t, p-I-I
(-1)m-~Lm(~[C;x,,h) = ~,, m=l
whence the term 'continuous' in the name of this class of methods. Victory and Ganguly (1986) have shown (Lemma 4 of the cited work) that the resulting system has a unique solution for the Lm(qT] C), and clearly this solution gives expressions of the form (7), with Lm($1 C; x~, h ; @) : = L,, (~ ] C; xt, h), as required for closure approximations. With these values of Lm(qT] C; xt, h), formula (15) then provides the required approximation of the form (5), and thus defines the underlying LOF method. The diamond-difference method is the case p = 1 of the class of continuous methods thus defined.
182
P. NELSON
Example E; the linear discontinuous method (Hill, 1975) The conventional interpretation (Alcouffe et al., 1979) of this method involves a linear approximation to the within-cell angular flux, which is conveniently written in the form ~
X -- X l
qT(x) = 4~(x,+) + ~
~
[q,(x,) - qY(x,+)].
By taking suitably approximating S, and then taking the appropriate moments of equation (4) over C [for details see Alcouffe et al. (1979)], the parameters ~(xl+)and qT(xr), as required to determine ~, are obtained as qT(x,+) -
6 + 4~ qT(x,)+ e 6 + ~.(4+ e) a [ 6 + ~ ( 4 + 8)] [(I +~)S(x,+)-S(xr)]
(16a)
and 6-2~
~(x,) - 6 + ~(4 + ~) qT(x,)+ ~[6 + e(4 + ,1)][3S(x~+) + (3 + e)S(xr)]
(16b)
where : = ~h/,u.
(16c)
Clearly formulas (16b,c) give an approximation to ~O(x,) that has the form (5) of an LOF method. The corresponding linear functionals are the right-hand limit of S at the left cell edge,
L,(SI C;x,,h) := S(x,+), and the value of S at the right cell edge, L2(S I C;xl,
h) :=
S(xr).
The corresponding closure approximations to L~(~,] C;x~,h) and L2(~]C;x~,h) are simply provided by equations (16a) and (16b), respectively. Note that generally qT(xt+) 4: qT(x~), whence cometh the term 'discontinuous' in the name of this method. Hennart (1981) has proposed a modified interpretation of the underlying linear discontinuous LOF scheme, in which the approximate angular flux is a continuous (piecewise) quadratic. It is not clear what the appropriate closure approximations would be, within the framework of this interpretation. A variety of additional specific examples of closed LOF methods could be given. However, those already indicated, plus others that easily can be provided by the interested reader, serve to make the major point of interest here, namely that many spatial approximations to 1-D transport, including nearly all such approximations that are designed for solution by an analog of the inner iteration process, are instances of closed LOF methods. Thus general properties of this class of methods would seem to be worthy of study. 4. CONSISTENCY OF CLOSED LOF METHODS
We begin with the concept of consistency for an LOF method, as introduced by Nelson and Zelazny (1986). For this purpose let, for/z > 0, ~(x ; xl, Y0) denote the angular flux at x~ C, corresponding to incident flux Y0 at xt; i.e. let x -~ ~b(x; xt, yo) be the exact solution of equation (4) that satisfies ~b(xt ; xt, Yo) = Yo- For a given LOF method (5), the corresponding truncation error over C is then defined by 1
~
T(x,, h ; Yo) --- h [~b(x,, h ;Y0) - ~0(x, + h ; x,, Y0)].
(17)
Such a method then is said to be consistent with equation (4) if the limit lira T(xt, h ;Y0) = 0
h~0
holds uniformly for Y0 in an arbitrary bounded set, 0 ~< xt ~< a, and x~ away from material interfaces ; i.e. if, given positive constants e and b, there exists h0 >0, possibly depending on ~ and b (and/~), but not on xt or Y0 (for the given fixed piecewise continuous ~ and S), such that if [Y0] < b,0 < h < h0, 0 ~< x~ ~< a - h , and there is no discontinuity of a or S lying strictly between xl and xt+h, then [T(xt, h;y0)[ < e. The method is
Closed linear one-cell functional methods
183
said to be consistent of order n with equation (4) (or simply order n) if T(xl, h ; yo)/h" is similarly uniformly bounded as h ~ 0. A proof of the following algebraic characterization of consistent LOF methods can be given by using the easily computed analytic form of q/(x ; xt, Y0). Theorem 1 (Nelson and Zelazny, 1986). Suppose o and S are piecewise continuous on [0, a]. Then an L O F method (5) is consistent with equation (4) if, and only if, 1
lim
h~0+ h
P
m-:---~,bm (O(X,+), h,/0Lm (S I C; x,, h) = S(xt+ )/# =
and lim h~ 0+
r(a(x,+), h, ~ ) -
1
h
= - tr(xl+)/It
hold uniformly for 0 ~< xt ~< a, xt away from material interfaces, and # bounded away from zero. In the same m a n n e r we define consistency of a closure approximation (7) as follows: the associated Lm
truncation error over C is 1
Tm(x,,h; yo) = ~[£m(S[ C ;x,,h; y o ) - Lm(~k(" ;x/,yo)l C ;x,,h)], and the closure approximation is consistent with equation (4) if lim T,,(xt, h;yo) = 0,
h~0
again uniformly for bounded Y0, x~ s [0, a], and x~ away from material interfaces. Likewise the closure approximation is consistent of order n with equation (4) (or simply order n) if T(xt, h ;yo)/h" is similarly uniformly bounded as h ~ 0. In the present generality we do not have available an analytic form of L,,(- ; xl, h) to use in establishing an algebraic characterization of consistency for closure approximations analogous to that provided by Theorem 1 for LOF methods. However, all that is really needed to establish such a characterization is a suitable assumption on the behavior of Lm(" ; Xl, h) for small h, as we demonstrate in our next result. For purposes of stating this result, let us introduce the notation BM = { f : f is a real-valued continuous function defined on [0, a] with piecewise continuous first and second derivatives there, such that If<'~(x)l ~< M for i = 0, 1, 2 and x ~ [0, a] such that f ' ( x ) and i f ( x ) exist}, where M is an arbitrary positive number. Theorem 2. Suppose that for any continuous and twice piecewise continuously differentiablef defined on [0, a], there are constants a0, and at such that
Lm ( f l C; x,, h) = a o f (x,+ ) + a, f" (x,)h + o(h), uniformly in f belonging to any BM and xl~ [0, a] lying away from discontinuities o f f ' or f". Then a closure approximation to L,,, of the form (7), is consistent with equation (4), for piecewise continuously differentiable S if, and only if, the limits 1
lim
P
S~ b,,k (a(x~+), h, #)Lk (SI C ; x/, h) = a, S(x~÷)/#
h - 0 + ~ k~,=
(18a)
and lim
h~0+
Zm(a(xt+ ), h, la) - a o h
=
- a oa (xt+)/#
(18b)
184
P. NELSON hold uniformly for x~ in [0, a] and away from material interfaces (i.e. discontinuities of a or S), and [p[ bounded away from zero.
Proof From the assumed form for L,,, we compute
l{
E
Tm(x~,h;yo)=~ [zm(~(xt*),h,l~)-ao]yo+
~ b,,k(a(xr),h,~)Lk(SIC;x~,h)-aiS(xt*)h/#
k=l
] } +o(h) ,
(19)
where for any fixed S the asymptotic estimate is uniform in xz~ [0, a] lying away from material interfaces, [#l bounded away from zero, and Yo lying in any bounded set. (The restriction of y0 to some bounded set, along with piecewise continuous differentiability of S and I~1 being bounded away from zero, insures that the relevant O(';xt, y0) belong to some B~ not depending on Yo or #.) If the subject closure approximation is consistent, then by taking Yo = 0 we see that equation (18a) must hold, in the uniform sense described. Similarly equation (18b) follows, in the appropriate sense, by choosing any other value of Y0. Conversely, if equations (18a,b) hold in the uniform sense described in the statement of the theorem, then equation (19) clearly shows that T,~(xt, h;yo) ~ 0 as h --. 0, in the uniform sense required for consistency. This completes the proof of Theorem 2. A closed LOF method is consistent if the underlying LOF method and all the associated closure approximations are consistent. A consistent closed LOF method has closed order n if (i) each closure approximation (7) has order nm, where (ii) T(xz, h;yo)/h n is bounded as h - ~ 0 , uniformly for Y0 in an arbitrary bounded set and x~e [0, a] lying away from material interfaces, with T(xt, h, yo) defined as was the truncation error T(xt, h, yo) in equation (17), except with each Lm(SI C; xz, h) in the expression (5) for ff(xz, h ;Y0) replaced by Lm(S [C; Xl, h) + O( h"-). Obviously, it is desirable for a closed LOF method to have its closed order no less than the order of the underlying LOF method. Our next result shows that this indeed holds for most of the specific methods introduced in the preceding section. Theorem 3. If ~ is positive everywhere, then the moment characteristic and continuous methods on p basic linear functionals have both their order and closed order equal to 2p.
Proof We refer to Nelson (1987) or to Victory and Ganguly (1986) for proofs of the assertions that the underlying LOF methods have order 2p. Granted this, the argument that they retain these orders as closed LOF methods is as follows. For the moment characteristic methods, the moment equations (14), along with the order of the underlying LOF method, as just described, imply that the closure approximation to Lm is consistent of order 2 p - m . (The hypothesis that tr is positive, along with the earlier assumption that it is piecewise constant, serves to bound tr above zero.) Now note that Taylor's theorem applied to x' ~ exp [-(T(x r--XZ)/lU], along with maxlpm(x)l = O(1)ash ~ 0 , x~C
shows that ceXp [ - O'(Xr- -
Xz)/IA]Pm-
1(x/) dx"
=
O(h") as h ~ 0.
It then follows that if S,,_1 in equation (13) is replaced by S,,-1 =Sm l+O(h 2p "), then the resulting approximation to ~(x,) is also accurate within O(h2p). But this is precisely the statement that the method has closed order 2p. For the continuous methods, the computation of the orders of the closure approximations is exactly the same. In Victory and Ganguly 0986) and Nelson (1987) it is shown that the basic LOF approximation (5) for the continuous methods can be written as
Closed linear one-cell functional methods ~ Pe(e;r) h'-' ~=Ot~+;~__0Sm(-1)"
185
te'm~P'(e;m)]e -(~+')
(2m)!Fp t e . m a - n m! [_ . . . . . . . . .
Qp(e;m)]
'
where e is as in equation (16a), and P,,(e ; m)/Q,,(e; m) is the mth diagonal Pad6 approximant (Baker, 1975) to exp ( - e ) . As such a Pad6 approximant has a Taylor's series, about e = 0, that agrees with that of exp ( - e ) to the first 2m+ 1 terms, and further each Qm(e;m) is bounded as h ~ 0, it follows that the coefficient of S,, in this equality is O(h ~+ ~). Therefore, if S,, is replaced by S,, = Sm+ O(h zp-m- ~), then the resulting approximation to ~r is of order 2p. This shows that the continuous methods have closed order 2p, and completes our proof of Theorem 3. We remark that the calculation in Nelson (1987) of the orders of the moment-characteristic and continuous methods as LOF methods is closely based upon the techniques of Victory and Ganguly (1986) ; however, the latter actually prove a somewhat different (and stronger) result, as will be elaborated in the next section. The following example shows it is possible for a (reasonable) closed LOF method to have its closed order less than the order of the associated LOF method. See also Neta and Victory (1982) for a discussion of this method from a somewhat different point of view.
Example F; the quadratic method of Gopinath et al. (1980) For this method the underlying LOF approximation is the same as that for the quadratic moment-characteristic method, which is to say the case p = 3 of the moment-characterstic methods. Therefore its order as an LOF method is 2p = 6. The closure approximation to L t(~k ] C;x~, h) also is the same as for the quadratic moment-characteristic method, namely that defined by formulas (12)-(14). However, for m = 2, 3, formula (14) is replaced by /22 : = ( ~ r - @ ) / 2
and
//,3
'--
I~r--[~l iLl. 2
Here ~r and £~ are defined by formulas (12)-(14), and the resulting expressions for £2 and £3 are readily seen indeed to have the form (7) of closure approximations. These closure approximations are easily seen to have respective orders 3 and 2. The corresponding values of T, as in part (ii) of the definition of 'closed order' then are both seen, as in proof of Theorem 3, to be O(h 4) as h -~ 0. Thus, this method has closed order of at least 4. Computational experiments by Larsen and Miller (1980) indicate that there are problems for which its (closed) order is exactly 4, and thus less than the order of the underlying LOF method. To this point we have not used the factorization b,,(~, h, p)L,, (S [ C; x~, h) appearing in the approximate source term for an LOF method (5), nor the corresponding factorization in the closure approximations (7). That is, for all considerations up until this point these terms could equally well have been written in unfactored form, say/~,,,(S[ C;x~,h,g, #). We now describe some considerations that make essential use of this factorization. A minor reason for writing these terms in factored form is simply the observation that for all actual methods presently used or proposed, the dependence on ¢rand p splits offas such factors. However, note that dependence on h appears in both factors. The major reason for such factored form is to capture the asymptotic dependence upon h in the {b,,}, and thus isolate it from the functional dependence upon S. Toward that end, let us define a family {L,,(- ; xt, h)} of basic linear functionals to be normalized if, for each m and xte [0, a],
~= o(1),
IIL,.(';xt,h)ll [ ¢ o(1)
as h ~ 0 .
Here and henceforth the norm of such a basic linear functional is that as a linear functional on the space of continuous real-valued functions defined on [0, a]. Note that the families of basic linear functionals for the methods of Section 3 are normalized. This usually is the case, for the natural definitions of the basic linear functionals. Further, it always can be arranged, except possibly in the 'degenerate' case that some Lm(" ; xt, h) has norm zero for all sufficiently small positive h. (This normalization may lead to dependence upon the 'reference points' x~ in the b,,, but this would not affect any of our results.)
186
P. NELSON
*o (=
o) ~'1
x,t (= a) Xll
~./-1
Xi
Xj+I
•17J- 1
ZJ-2
Cj-1
C,
[
Cj
hj-l ~hj---'~
~,-hx "t"
Fig. 2. A mesh of spatial cells, and the associated standard notation. Consider now a consistent LOF method on a normalized family of basic linear functionals. Theorem 1 insures that P ; xt, h) = O(h) as h ~ 0 , uniformly in xte[0, a] lying away from ~=tbm(a(xz+),h,#)Lm(" material interfaces and # bounded away from zero.
(20)
In the ensuing we will wish to use this asymptotic estimate to infer the consequence that
b,,(a(xt+), h,/~) = o(1) as h ~ 0, uniformly in xze [0, a] lying away from material interfaces and p bounded away from zero, for m = 1. . . . ,p. (21) In general this inference is incorrect, as can be seen by easily constructed examples. However, let us define a normalized family {L,,(- ; xt, h)} of basic linear functionals to be asymptotically linearly independent if for any xt~ [0, a] the asymptotic estimate (20), with O(h) replaced by o(1), implies formula (21). If an LOF method has asymptotically linearly independent basic linear functionals, then the desired inference clearly does hold, along with the analogous result for closure approximations (7). We record this as a formal result, for ease of later reference. Lemma 1. Let formula (5) be an LOF method such that the associated family of basic linear functionals is asymptotically linearly independent. If this method is consistent with equations (3), then the coefficients {b,,} satisfy formula (21). If further the basic linear functionals satisfy the hypotheses of Theorem 2, and the associated closure approximations (7) are consistent with equation (4) (so that the closed LOF method is likewise consistent), then the coefficients {bmm,} satisfy bmm,(a(x~+),h, p ) ~ 0 as h ~ O, uniformly in x~E [0, a] lying away from material interfaces and/~ bounded away from zero, for 1 ~< m, m' ~
ARE ASYMPTOTICALLY FINITE-DIFFERENCE METHODS
Closure approximations are defined precisely so that one can carry out an analog of the inner iterations (3) within a discrete spatial approximation. We now wish to consider some of the details of such a discrete inner interation. Thus we must consider a mesh of cells, say ~ = {Ci}, where Cj = [xj, xj+ d, with the remaining notation as shown in Fig. 2. By contrast, the previous discrete considerations have been limited to a single cell. The discrete version of ray tracing consists simply of the explicit (in j ) difference equations P
> O,
(22a)
~O,+l) * Fid+~(n+t)-I-1 ,5 ..F("!. ij = ~ij' ~ ~ ~mq ~mo, ]~i < O,
(22b)
~n,j++]) = .Cij~};+ l, ~_ E b , . o £ ~ ,
Iz,
m=l
m=l
ff}g+o=f
foru,>0,
~}7+')=f
for/ti<0,
(22c)
where z0 : = z(aj, hi, #~),
(23a)
Closed linear one-cell functional methods
187
bmi~ : = b~(aj, hi, #3
(23b)
and aj : = a(x),
xeC ° := (xj,xj+,).
(23c)
E("! Here "~/7~" r U + ° is to be interpreted as an approximation to O("+ °(xj, #3, and --mu as an approximation to L,,(S}")ICj;xj, hj). Equations (22) are the LOF equations (6a), or their counterpart for #i < 0, in application to the particular cell Cj. In consonance with equation (6b) and the interpretation of £~]j, the latter are given in terms of £ ~ j ~ L,.(01") I Cj ; xj, hi) by N
£ ~,~?,j= ~ wr k~,7£~],j + L,,, (q~ ] C ; x j, hi),
(24a)
i'=1
where
kirj : = kii,(x),
x e C °.
(24b)
It is convenient to consider the t~("~ as the beginning of an iteration. If these are known, the E{"!. -,,u can be --m U computed from formulas (24), and then the ",IT(,+ i) from equations (22), subject to the known values of rq ~ ( n + l) i0 = ~O(0,/~i) = f and ws ,FA-+1) = ~(a, #i) = f , for #i > 0 and #t < 0, respectively. The -£(,.+ - m U o, as needed to close the iterative loop, are computed by the closure approximations (7). In application to equation (3a), over the cell Cj, these become p
(,+ 1) _-- ~ ~e -~ij'r,F,(,+ q 1)+ ~ ~hm m ... ij -E~"?.. - m lj,
Pi > O,
(25a)
]~i <
(25b)
m'=l
and p
£(,+ m i j 1) -_- . ~ m i j ' ~,r,(,+ ' ij+ o I al-
E m'=l
h
,.. E("?..
v m m tj - - m U'
O,
where L,,, (S}") I Cj ; xj, hi) has been approximated by _,,,,j,t.(")..as is consistent with the interpretation of the latter. Here we have also introduced an obvious extension of the notation (23). In carrying out this iterative scheme, one presumably would compute the appropriate approximate flux functionals, via the closure approximations (25), more or less simultaneously with the discrete ray-tracing calculations (22). The required values of the approximate source functionals would be computed from formulas (24), as needed. This requires that the p N J values of the approximate flux functionals be stored at the end of each iteration, in order to be available as needed within the next iteration. Thus the number of basic linear functionals, p, is also the amount of storage, per cell and discrete direction, required by a closed LOF method. (For some methods it will be necessary to correct this to account for the fact that a given approximate source functional is used in each of two adjacent cells. In practice the iterative procedure just described is often modified so that the approximate source functionals, {/2~}, are stored at the end of each iteration, for use in the next. The advantage of this is that often these are only weakly dependent on the direction index (t), so that savings in storage might be possible through some type of data compression in this index. For an extreme example, if the scattering and sources are isotropic, then the {/2~]j} are independent of i; more generally, a few ( < N ) discrete Legendre coefficients may suffice. For monoenergetic problems this is not particularly useful as, in any case, the approximate flux functionals, {£~]j} must be temporarily stored during the ray-tracing calculations (22) and (25), in order to permit ultimate computation of the approximate source functionals via formulas (24). However, for multigroup problems this strategy can, particularly for scattering that is only slightly anisotropic, yield a significant reduction in storage requirements. Regardless of such considerations, formulas (22), (24) and (25), along with the boundary conditions inherited from (2), do define an iterative scheme. For certain particular closed LOF methods, including the diamonddifference and step-characteristic methods, conditions are known (Menon and Sahni, 1985 ; Nelson, 1986) that assure convergence of this scheme. If we assume this convergence, then the limiting values (Lmij = niim £~],J ~o etc.) satisfy
188
P. NELSON
~i,j+l = "CiJ~ijq- ~ bmijEmij,
~i • O,
(26a)
rn=l
~i,j = rij~i,j+l J[- ~ b,~u £mij,
]~i < O,
(26b)
m=l N
£,~j = ~ wi,k~rj£,~rj + Lm(qi l C j ; xj, hj),
(26c)
i'=l P
£,,u = z,,~:~u+ ~ bm,-'u/Sm'u, ]~i > 0,
(26d)
rW--I
£mij = Tmij~i,j+ I "JU~ bmm,ijff~m,ij, m'--I
]2 i
(26e)
< O,
ff, o = f ,
/t~ > 0,
(26f1
OiJ=f,
~,<0.
(26g)
and
Regardless of the question of convergence of formulas (22), (24) and (25), one can consider the question of whether approximation (26) has a unique solution and, if so, the manner (if at all) in which q7u approximates 0(xj, #~), £mU approximates Lm(S~ICj;xj, hi) etc. [In fact we could equally well have used the basic LOF and closure approximations to write down equations (26) directly as an approximation to problem (1, 2).] This is the setting in which Victory and Ganguly (1986) carried out their fundamental studies of consistency, along with convergence and stability, for the moment-characteristic and continuous methods. In this sense, the results of the cited work also include that part of Theorem 3 dealing with the closed orders of these methods. The discrete approximation (26), as defined by a closed LOF method, involves three sets of discrete dependent variables, namely the approximate cell-edge fluxes {ffij}, the approximate source functionals {Zmij} a n d the approximate flux functionals {£mU}" In general it does not seem possible to eliminate any of these, in terms of the others, and thereby obtain a system of reduced size. However, for a consistent LOF method, a sufficiently small mesh, and under further reasonable assumptions, it turns out that any two of these sets can be eliminated, thereby leaving a system of equations with the remaining set as the only unknowns. The purpose of the remainder of this section is to provide precise formulations and proofs of these (three) assertions. We begin with : Lemma 2. Suppose that the closed LOF method underlying equations (26) has asymptotically linearly independent basic linear functionals, that the associated closure approximations are consistent with equation (4), and further that the mesh ~ has each material interface (i.e. discontinuity of the piecewise constant functions k or a) coinciding with some cell edge. Then there exists some e such that if II~ II : = l<~j<~J max hi~ < e,
(27)
then the system consisting of equations (26c-e), for m and j ranging between 1 and p or 1 and N, respectively, and fixed j subject to 1 ~< j ~< J, has a unique solution for either £mU, or £,,ij, in terms of the ~ej having #~ > 0 and the ff*.a+l having/z, > 0.
Proof. By substituting equation (26c) into equations (26d,e), we obtain equations of the form p
N
£m,j = Z
E
m'=l
i'=1
amim,,,jLm,,,j+Bm,j,
where the flmU involve the q7u and qTu+ ~, as described in the statement of the theorem, and otherwise known data. By Gerschgorin's theorem, it suffices to show that there exists e such that II # IP < ~ implies
Closed linear one-cell functional methods
189
N
• IA..,m',vl < 1, m'=l
(28)
i'--I
for e a c h j = 1 , . . . , J. But if#i > 0, then the 1.h.s. of inequality (28) is equal to p
N
Z
~ Ibmm,ullwrllk~,jl,
m'=l
i'=l
where bmm,ij = b m m , ( ~ j , h j , ]~i)"
As the w~, and k~rj are uniformly bounded, over all possible combinations of subscripts and meshes, existence of e as desired follows from Lemma 1, along with the requirement that each material interface be a cell edge. For #i < 0 the required estimate (28) is obtained similarly. This completes the proof for the {£mU}' That for the {/Sin0.is similar, which completes our outline of a proof of Lemma 2. If the solution for the {£mU},obtained as described in Lemma 2, is substituted into equations (26a,b), then there results a system of equations in only the approximate cell-edge fluxes, {qTu}. This system is equivalent to system (26), in the sense that the corresponding {/2,,u} and {£mU} can be obtained as described in Lemma 2. This is the essence of the proof of the result encapsulated in the title of this section. A precise detailed statement of this result is : Theorem 4. Suppose the hypotheses of Lemma 2 hold. Then there exists e > 0 such that if expression (27) holds, then there exist (mesh-dependent) coefficients B~evand ~u, 1 ~< i, i' ~< N, 1 ~< j ~< J, such that : (i) if the system (26) has a solution for {~u}, {£,,u} and {£mu}, 1 ~< i ~< N, 1 ~
~,.j+ , = ~,j~,j + E B..j~,j + E B.,j~.j+ , + ~,j, #~ > O, i'
i" + i"
ff~0=f
for#i>0,
i"
tffij=f
l~i<.N,
forp~<0,
l<.j<~J,
(29)
where ~ ( ~ , ) indicates summation over only those i' such that/~r > 0 ( < 0) ; (ii) conversely, if system (29) has a solution for {flu}, for 1 ~< i ~< N and each j, 1 ~< j ~< J, and {/2mu}, {£mU} are calculated from these {flu} as described in Lemma 2, then the resulting values of {qTu}, {£~u} and {£,.u} comprise a solution of the system (26), for 1 ~< i ~< N, 1 ~< j ~< J and 1 ~< m ~
lffq =
k= 1
l
1 "~il
m= 1
bmik£,nik
+
Zi~ f ,
/~, > 0,
(30a)
and lifo =
bmikZmik +
"~il k =j L \ l=j
1
rit f ,
#, < 0.
(30b)
=j
We display this solution explicily, because it has some inherent interest as the discrete counterpart of the explicit solution of equation (4), with kt = #i, S given, and subject to the boundary conditions (2a,b) By means of equations (30), the system for the {£,,u} in terms of the {qTu}, as described in Lemma 2, can be converted into a system involving only the {/2mu}. Likewise, by using equation (26c), one arrives at a system involving only the {/2mU}as unknowns. We formally record the latter result.
190
P. NELSON Theorem 5. Suppose the hypotheses of Lemma 2 hold. Then there exists e > 0 such that if
formula (27) holds, then there exist coefficients C,~,.,ii,jj,and flmis, 1 <~ m, m' <<.p, 1 <~ i, i' <~ N, 1 ~< j, j ' ~< J, depending upon the mesh and the data for the boundary-value problem (1, 2) such that : (i) if the system (26) has a solution for {~o}, {£mij} and {£m~j}, where the indices range as above, then the {£mlj} satisfy p
£mij = E m'=l
N
J
E ECmm'ii'ii'jj'£m'i'J'~-flmij , i'=l j'-I
l<.m<~p,
l<~i<~N,
l<~j<~J;
(31)
(ii) conversely, if equation (31) has a solution for the {£m~S}, and {£mO}, {qT~j}are calculated by equations (26c) and (26a,b) [or equations (30)], respectively, then the resulting values comprise a solution of system (26), for all relevant values of the indices. There is, of course, a similar result for the approximate source functionals, {£,,o}. None of these (asymptotically possible) reduced forms of the system (26) have been exploited in studies of specific closed LOF methods. Rather the standard approach has been to consider the system in the {qTtj} and {£m~j} that results from eliminating the {£mO} by equation (26c). For example, in Nelson (1987) an explicit solution is given for the system in the case of monodirectional (N = 1) continuous methods. This is how the explicit expression for ~r for these methods, as given in the preceding section, was obtained. Victory and Ganguly (1986) even obtain an analytic solution of this system for continuous methods and a general value of N. This solution is a key to their study of convergence, stability and discretization error for the continuous methods.
6. C O N C L U D I N G R E M A R K S
Our primary objective in introducing the class of closed LOF methods is to provide a suitable framework for a general theory of spatial approximations to the monoenergetic discrete-ordinate equations in slab geometry, and possibly even eventually in general multidimensional geometries. An essential step in developing such a theory, and thereby filling out this framework, is to show that a consistent closed LOF method is both stable and convergent to the solution of the underlying continuous problem. This would generalize previous results for a number of specific closed LOF methods. One possible two-step approach to obtaining this desideratum would be as follows. First, extend the basic result of Nelson and Zelazny (1986) to show that a consistent closed LOF method is both stable and convergent for the monoenergetic discrete-ordinate equations subject to initial conditions. Secondly, extend to closed LOF methods the result of Keller and White (1975) to the effect that, for a linear system of differential equations, stability and convergence of a given finitedifference approximation for linear boundary conditions is equivalent to these same properties for initial conditions, provided both the initial and boundary problems are well-posed. A second desirable objective is that of showing convergence of the discrete inner iteration process (22)-(25), for a consistent closed LOF method and sufficiently small mesh. Existing methods (Menon and Sahni, 1985 ; Nelson, 1986) of proof seem essentially limited to methods of relatively low order. The collectively compact operator theory of Anselone (1971) provides a possible point of attack in attempting to obtain such a result. This theory requires an associated 'injection' of the discrete approximate angular fluxes back into the space of continuous angular flux; i.e. each discrete approximation to the angular flux must be associated to a continuous approximation. For the moment-characteristic and continuous methods there is a natural definition of such continuous approximation, although the desired convergence apparently has not even been established for these methods. For general abstract closed LOF methods there is no such 'natural' definition of an associated continuous approximation. Our Theorem 4 suggests the possibility of simply using linear interpolation, although obviously this will not provide an extremely accurate continuous approximation. Finally, it would seem that the clarification of the requirements of a discrete spatial approximation, as provided by the explicit description of closed LOF methods, could provide a stimulus for the development of novel specific methods. Most methods that have been used or proposed are Galerkin methods (i.e. the continuous methods) or 'improved' Galerkin methods (i.e. the moment-characteristic methods), or otherwise involve piecewise polynomial approximations to the source function (S), so that the basic linear functionals are polynomial moments. However, Runge-Kutta methods essentially are closed LOF methods in which the basic linear functionals are point evaluations. These do not seem to have been used in transport theory. Also,
Closed linear one-cell functional methods
191
there would a p p e a r to be some possible advantages (Nelson a n d Zelazny, 1986) to using (closed) linear multicell methods, in analogy to the linear multistep m e t h o d s (e.g. L a m b e r t , 1973) t h a t are widely used for numerical solution of general initial-value probems. Acknowledgemen~This research was partially supported by the U.S. Department of Energy, under Contract No. DEAA03-76-SF00767. REFERENCES
Alcouffe R. E., Larsen E. W., Miller W. F. Jr and Wienke B. R. (1979) Nucl. Sci. Engng 71, 111. Anselone P. M. (1971) Collectively Compact Operator Theory. Prentice-Hall, Englewood Cliffs, N.J. Baker G. A. Jr (1975) Essentials ofPadb Approximants. Academic Press, New York. Gopinath D. V., Natarajan A. and Sundararaman V. (1980) Nucl. Sci. Engng 73, 181. Hennart J. P. (1981) Ann. nucl. Energy 8, 677. Hill T. R. (1975) Report LA-5990-MS. Keller H. B. and White A. B. Jr (1975) S l A M Jl numer. Analysis 12, 791. Lambert J. D. (1973) Computational Methods in Ordinary Differential Equations. Wiley, London. Larsen E. W. and Miller W. F. Jr (1980) Nucl. Sci. Engng 73, 76. Lathrop K. D. (1969) J. comput. Phys. 4, 475. Lee S. M. and Vaidyanathan R. (1980) Nucl. Sci. Engng 76, 1. Menon S. V. G. and Sahni D. C. (1985) Trans. theor, statist. Phys. 14, 353. Nelson P. (1971) J. math. Analysis Applic. 35, 90. Nelson P. (1973) S I A M J1 numer. Analysis" 10, 175. Nelson P. (1986) S l A M Jl numer. Analysis. 23, 1017. Nelson P. (1987) In preparation. Nelson P. and Zelazny R. (1986) Nucl. Sci. Engng 93, 283. Neta B. and Victory H. D. Jr (1982) Numer. Funct. Analysis" Optmizn 5, 85. Neta B. and Victory H. D. Jr (1983) S l A M Jl numer. Analysis 20, 94. Vaidyanathan R. (1977a) Atomkernenergie 29, 30l. Vaidyanathan R. (1977b) Proc. 5th Int. Conf. Reactor Shielding, p. 754. Science Press, Princeton, N.J. Vaidyanathan R. (1979) Nucl. Sci. Engng 71, 46. Victory H. D. Jr and Ganguly R. (1986) S I A M J. numer. Analysis 23, 78.
ANE
14:4-D
0306-4549/87 $3.00+0.00 Copyright © 1987 Pergamon Journals Ltd
Printed in Great Britain. All rights reserved
CLOSED L I N E A R ONE-CELL F U N C T I O N A L SPATIAL APPROXIMATIONS P. NELSON" Department of Applied Mathematics, California Institute of Technology, Pasadena, CA 91125, U.S.A.
(Received 29 July 1986) Abstract--Closed linear one-cell functional (closed LOF) methods are introduced as a class of discrete spatial approximations to the discrete-ordinate equations for anisotropically scattering monoenergetic transport in slab geometry. It is shown that most spatial approximations currently used or suggested for use, for this class of problems, are specific realizations of closed LOF methods. An appropriate concept of consistency is introduced for closed LOF methods. It is shown that a consistent closed LOF method is asymptotically a finite-difference method, in the sense that the approximating equations can be expressed in terms of approximations to the cell-edge fluxes, for sufficiently fine meshes. Further questions related to closed LOF methods are suggested for future study.
l. INTRODUCTION
We write the discrete-ordinate equations for anisotropically scattering monoenergetic transport in slab geometry as / z ~ x ' +a~bi = S~, i = 1. . . . . N,
(la)
where N
Si(x) : =
~ wrkir(X)Or(x)+qi(x).
i'=l
(lb)
Here the p~s [ - 1, 1] ~ {0} are the quadrature points, with associated quadrature weights {w~}, a = a(x) is the scattering kernel and q is the given source function. We shall assume that a and the kit are piecewise constant. Further we shall consider equations (la,b) subject only to boundary conditions representing specified angular fluxes incident upon the slab faces, say qJi(0) = f
for/~i > 0,
(2a)
~bi(a) = f
for #~ < 0,
(2b)
where a > 0. Suppose momentarily that a and the {wi}, {k~r}, {qi} and {f} are all nonnegative, as is the case in most instances of interest. In this situation, if there exists a nonnegative solution of equations (la,b) and (2a,b), as is known (Nelson, 1973) to be the case if the underlying exact transport problem is subcritical and the quadrature rule is sufficiently fine, then that solution must be [see for example, Nelson (1971)] the limit, as n -~ ~ , of the {~/,~n)}, where the latter are the inner iteration approximations defined by 1 ~ i ~ - [ - o - ~ 1 } n+l)
:
S~ n),
i = 1. . . . . N,
(3a)
and N
S} ")(x) = ~ wrk,r (x)~P},")(x) + q, (x).
(3b)
i'=l
t Permanent address: Department of Mathematics, Institute for Numerical Transport Theory, Texas Tech University, P.O. Box 4319, Lubbock, TX 79409, U.S.A. 177
178
P. NELSON h
X I
Xr
Fig. 1. Notation for a typical mesh cell, C.
Here the initial approximations {~b}°)} or {S}°)} are arbitrary, modulo suitable technical conditions, and the approximate angular fluxes {~k}")} also satisfy the boundary conditions (2a,b). The angular coupling inherent in equations (la,b) tends to recede to the background within the inner iteration process, in the sense that at each iterative stage (i.e. each fixed value of n) the approximate angular fluxes are determined by N monodirectional monoenergetic transport equations,
dq,
# dxx + trio = S,
(4)
where S is known. Further, each such ff is specified at only a single point, which is to say that, in mathematical terminology, equation (4) is to be solved subject to initial conditions. The process of solving such an initialvalue problem for equation (4) often is termed ray tracing. Because of the facts described in the preceding paragraph, spatial approximations for equations (1 a,b) and (2a,b) often are described explicitly in their application to the ray-tracing problem only. In Nelson and Zelazny (1986) the class of linear one-cell functional (LOF) methods for ray tracing was introduced, and the consistency, convergence and stability properties of such methods were studied. More specifically, it was shown that, in application to the ray-tracing problem, a consistent LOF method is both convergent and (zero) stable. The LOF methods seem to include very nearly all methods that have been used, or suggested for use, for raytracing calculations as they occur for a discrete spatial aproximation of equations (1) and (2) within the framework of an overall inner iteration scheme. Implementation of an LOF method for equation (4) requires, as detailed in Section 2, the values of certain quantities, the so-called basic linear functionals of the particular method, corresponding to the particular source(s) under consideration. Nelson and Zelazny (1986) assumed that any information regarding S that was required to construct these values was available as needed. In point of fact, within the iterative scheme (3) the relevant source depends, via equation (3b), upon the approximate angular flux generated within the preceding iteration. Therefore, the information required, at each iterative stage, to construct, possibly approximately, the required values of the basic linear functionals must have been generated (and stored) during the preceding iteration. In the present paper we show that this requirement can be met by appending certain closure approximations to an LOF method, and further we make an initial study of the resulting closed LOF methods. Closed LOF methods are defined in Section 2. In Section 3 we show how a variety of specific spatial approximations fall within the class of closed LOF methods. Section 4 is devoted to considerations relating to extension of the concept of consistency, as defined by Nelson and Zelazny (1986) for LOF methods, to encompass the closure approximations, and thus closed LOF methods. The basic result of Section 5 is a theorem to the effect that, under suitable assumptions on the associated basic linear functionals, a consistent closed LOF method asymptotically is a finite-difference method ; that is, we show that, for suitable consistent closed LOF methods and sufficiently fine spatial meshes, the required values of the basic linear functionals can be expressed in terms of the approximate values of the angular flux at the cell edges (mesh points). Some possible topics for further study, relating to closed LOF methods, are briefly outlined in the concluding Section 6.
2. CLOSED LOF METHODS
We begin by defining an LOF method for the ray-tracing problem for equation (4). Consider a typical spatial cell, C = [x~,Xr], as in Fig. 1. For/~ > 0 and given ~k(xt), possibly approximate, an LOF method for this ray-tracing problem is an approximation of the form
I/l(Xr) ~-~ ~(Xl, h ; ~b(xt)) : = z(a(xl+), h, #)~(x,) + ~ bm(a(xt+),h, #)Lm(S [C; xl, h), m=l
(5)
Closed linear one-cell functional methods
179
where dependence upon the data of definition (5) (i.e./~, a and S) has not been notationally included in the argument list for ~, because we usually consider these arbitrary but fixed, but all dependencies have been explicitly indicated on the r.h.s, of expression (5), for clarity. Here z and the bk are ordinary point functions, but for each fixed h > 0 the L,, (" ; xz, h) are linear functionals that are termed the basic linear functionals of the particular method. The notation a(x~+) indicates the right-hand limit of a(x) at the left cell edge. If the mesh is arranged so that material interfaces are cell edges, as is invariably done in practice and we henceforth assume for the remainder of this paper, then this will be the value of tr throughout the cell ; unless a material interface is located at xt, this also will be the same as the left-hand limit, a ( x r ) . The notation 'S I C ' ( = S restricted to C) indicates that the L,(" ; xt, h) are allowed to depend only upon the behavior of S within the cell C. F o r / t < 0 an approximation similar to formula (5) is assumed to hold, with the roles of xt and x~ reversed. By comparing formula (5) with the corresponding exact solution of equation (4). subject to given ~,(xt), it is clear that z(a(x~+), h, /~) should be an approximation to the exponential attenuation factor across C, in direction/~ and that the second term on the r.h.s, of equation (5) should approximate the uncollided flux, emergent from C, in direction #, stemming from sources within C. These considerations will be made precise in our discussion of consistency, in Section 4. We remark that one-step finite-difference methods of the type frequently used for general-purpose computing are instances of LOF methods ; however, in such methods the Lk(SI C; x~, h) invariably are taken as point evaluations of S, whereas other types of functionals often are employed within transport theory, as shown by the examples of the next section. In application to ray tracing for equation (3a), an LOF method has the form
~) + ~" ~0<.+ i.r 1) = z(a(xt+),h, I,,~ i"~,l,(n+ ~ bm(a(xl+),h,#,)L,,,(S}")lC;xt, h). lWi,l
(6a)
m--I
By virtue of equation (3b), and constancy of the kii. within C, the values of the basic linear functionals that appear in equation (6a) can be expressed as N
L,,(S} ") I C;x,,h) = ~ wrk,,,L,,(~b},") I C; x,,h)+ Lm(qt I C;xz, h).
(6b)
i'=1
As the qi are given (usually piecewise constant), we assume that the Lm(qi[ C; x~, h) can be computed satisfactorily, at least approximately. Thus, in order to close the iterative loop, in the sense of permitting the approximate cell-edge fluxes to be computed via equation (6a), it remains only to devise some satisfactory method for approximating the Lm(~l},n) l C ; Xl, h). As we only know that the ~,}")satisfy equation (3a), with n replaced by n - 1, the task just described essentially is that of approximating Lm(¢l C), for any given S, where ff and S are related by equation (4), subject to known ~b(xt) (~9(xr)) for # > 0 (/~ < 0). The form (5) of an LOF approximation suggests the use of closure approximations :
Lm(~l C;xt, h) ~ £,n(SI C;xt, h;~O(xt)) : = "Cm(ff(xt+),h,ll)~(xl)q- ~ bmm,(a(Xz+),h,#)L,n,(S ] C;xt, h), (7) m'--I
where it is to be understood that ~k and S are related as described above. The dependence upon tr and # is not explicitly indicated in the argument list for the £m, as for ff alone, but the dependence upon S is indicated, as consideration of these approximations for varying S is central to the use of closure approximations. A closed LOF method for equations (1) and (2) is simply an LOF method along with closure approximations of the form (7) f o r / t > 0, and corresponding approximations, but with xz replaced by xr and x~+ by xr , for # < 0. In Section 5 we show, in detail, how a closed LOF method provides an effective iteration scheme, analogous to scheme (3), for generating approximations to the solution of the basic underlying problem (1, 2). However, we take as our next task to provide a demonstration that a variety of well-known spatial approximations for problem (1, 2) are specific instances of closed LOF methods.
3. EXAMPLES OF CLOSED LOF METHODS
Throughout this section we denote the approximations to the solution of equation (4) at the left-hand and right-hand cell edge by q~l and ~ , respectively. Thus, for each LOF method considered,
180
P, NELSON
~L := qT(x,,h; ~7,)
(8)
in the notation of definition (5).
Example A ; the diamond-difference method This method traditionally is considered to be defined by the particle-balance equation, h (qTr- @) + ~r~ = S
(9)
q7 = (q7r + @)/2,
(10)
and the supplementary equation
where q7 and ~qare interpreted as approximations to the averages over C of ~Oand S , repectively. (Here, as in the remainder of this section, a denotes the constant value of the cross section within the particular cell under consideration.) These can be solved for qTr,in terms of ~qand @, to yield
1 - ah/2# ~
2h
This clearly has the form of an LOF method, in the notation of formulas (5) and (8), provided we take p = 1 and
L,(S[ C;x,,h) : = h
S ( x ) d x = ~q
(11)
as the sole basic linear functional. If we further take ~ as £ ~, the required approximation to L ~(~ [ C; xt, h), then there results
£ , ( S I C ; x , , h ; ~ , ) = 1 [1 + (2# - ah)/(2# + crh)]@ + 2# +h ah S, _ which clearly has the form (7) of a closure approximation, as required to display the diamond-difference method as a closed L O F method.
Example B; the step-characteristic method (Lathrop, 1969) This method is defined by replacing S in equation (4) by its cell average, formula (11), and then solving the result, subject to presumed known qT(x,) = ~, (for # > 0), for ~ = ~(x), an approximate angular flux defined throughout the cell. This yields 1
~r = exp ( - ah/#)~, + ~ [1 -
exp
( -
ah/#)]~,
which clearly is of the form (5), as required to comprise an L O F method. If, further, the required approximation to L ~(~k I C; xl, h) is taken as L 1(~ [ C; xl, h) -- £ ~, then we find
£,(Sl C ; x , , h ; ~ , ) = ~h which is precisely of the form (7), as required to constitute a closure approximation. Thus the step-characteristic method also can be regarded as a closed L O F method.
Example C ; moment-characteristic methods o f Vaidyanathan (1977a, b, 1979 ; Lee and Vaidyanathan, 1980) These are a generalization of the step-characteristic method. They are defined by again replacing S in equation (4) by an approximation, and solving the result, subject to known qT(x3 = ~t (# > 0), for q7 = qT(x), which is to be interpreted as an approximate angular flux defined throughout the cell. The approximate source used is
Closed linear one-cell functional methods p
181
1
g(x) := y~ Smpm(x), m=0
where the p,.(x) are the cell-based Legendre polynomials of Larsen and Miller (1980) and the p basic linear functionals of the method are defined by
am
=
Lm+ ,(SI C;x~, h) . - 2 mh+ l fc S(x)p,,(x)dx,
m = 0. . . . , p - 1 .
(12)
The resulting value of ~r is ? P if, = ff(xr) = exp ( - ah/#)@ + E S,,_, | exp [-- a(x, -- X')/ItlPm_l (X') dx', m=0
(13)
dC
which clearly has the required form (5) of an LOF method.By using fundamental properties of the Legendre polynomials, a simple recursive formula can be given for the integrals in equation (13), but that is not of real concern here. For these methods the required closure approximations to Lm(~OI C; h) are simply taken as Lm(i~ [ C; h). The latter are perhaps best computed by multiplying equation (4), with S replaced by S and ~ by qT, by Pm 1(x) and integrating the result over C. After some manipulations, there results [(m- 2)/2]
# --(2m--1)[gTr+(--1)m~7,]+2(2m--1) £m := Lm(qTIC ; x , , h ) = ~h
Y~
}
Lm 1 2,(OIC;x,,h)
l=0
+
1
cr
Lm(SlC;x,,h),
m = l . . . . . p,
(14)
where the arguments of £,, are as before and [x] denotes the greatest integer not exceeding x. Upon using equation (13), it is clear that formula (14) can be employed recursively to yield closure approximations of the required form (7). The detailed form of the resulting closure approximations is not essential to our purposes. Clearly the step-characteristic method is the case p = 1 of this class of closed LOF methods.
Example D ; continuous methods (Neta and Victory, 1982, 1983 ; Victory and Ganguly, 1986) In these methods the basic linear functionals again are given by definition (12), and the moment equations (14) again play a major role. The difference from the moment-characteristic methods stems from the fact that now the approximate angular flux within the cell is taken as a polynomial of degree p, p+l
~ ( x ) = ~ Zm(tfflf;x,,h)pm l(X). m-I
This yields p+ I
qT~ : = qT(x~) = ~ Lm(~lC;x,,h),
(15)
m=l
and the result of substituting this into the moment equations (14) can be regarded as a system o f p equations in the p 4-1 unknowns Lm(~[ C; xt, h), with known terms (i.e. 'right-hand side') depending linearly upon ~t and the Lm(S] C; xt, h). This is made into a determined system of equations by requiring that ~(xt) = ~t, p-I-I
(-1)m-~Lm(~[C;x,,h) = ~,, m=l
whence the term 'continuous' in the name of this class of methods. Victory and Ganguly (1986) have shown (Lemma 4 of the cited work) that the resulting system has a unique solution for the Lm(qT] C), and clearly this solution gives expressions of the form (7), with Lm($1 C; x~, h ; @) : = L,, (~ ] C; xt, h), as required for closure approximations. With these values of Lm(qT] C; xt, h), formula (15) then provides the required approximation of the form (5), and thus defines the underlying LOF method. The diamond-difference method is the case p = 1 of the class of continuous methods thus defined.
182
P. NELSON
Example E; the linear discontinuous method (Hill, 1975) The conventional interpretation (Alcouffe et al., 1979) of this method involves a linear approximation to the within-cell angular flux, which is conveniently written in the form ~
X -- X l
qT(x) = 4~(x,+) + ~
~
[q,(x,) - qY(x,+)].
By taking suitably approximating S, and then taking the appropriate moments of equation (4) over C [for details see Alcouffe et al. (1979)], the parameters ~(xl+)and qT(xr), as required to determine ~, are obtained as qT(x,+) -
6 + 4~ qT(x,)+ e 6 + ~.(4+ e) a [ 6 + ~ ( 4 + 8)] [(I +~)S(x,+)-S(xr)]
(16a)
and 6-2~
~(x,) - 6 + ~(4 + ~) qT(x,)+ ~[6 + e(4 + ,1)][3S(x~+) + (3 + e)S(xr)]
(16b)
where : = ~h/,u.
(16c)
Clearly formulas (16b,c) give an approximation to ~O(x,) that has the form (5) of an LOF method. The corresponding linear functionals are the right-hand limit of S at the left cell edge,
L,(SI C;x,,h) := S(x,+), and the value of S at the right cell edge, L2(S I C;xl,
h) :=
S(xr).
The corresponding closure approximations to L~(~,] C;x~,h) and L2(~]C;x~,h) are simply provided by equations (16a) and (16b), respectively. Note that generally qT(xt+) 4: qT(x~), whence cometh the term 'discontinuous' in the name of this method. Hennart (1981) has proposed a modified interpretation of the underlying linear discontinuous LOF scheme, in which the approximate angular flux is a continuous (piecewise) quadratic. It is not clear what the appropriate closure approximations would be, within the framework of this interpretation. A variety of additional specific examples of closed LOF methods could be given. However, those already indicated, plus others that easily can be provided by the interested reader, serve to make the major point of interest here, namely that many spatial approximations to 1-D transport, including nearly all such approximations that are designed for solution by an analog of the inner iteration process, are instances of closed LOF methods. Thus general properties of this class of methods would seem to be worthy of study. 4. CONSISTENCY OF CLOSED LOF METHODS
We begin with the concept of consistency for an LOF method, as introduced by Nelson and Zelazny (1986). For this purpose let, for/z > 0, ~(x ; xl, Y0) denote the angular flux at x~ C, corresponding to incident flux Y0 at xt; i.e. let x -~ ~b(x; xt, yo) be the exact solution of equation (4) that satisfies ~b(xt ; xt, Yo) = Yo- For a given LOF method (5), the corresponding truncation error over C is then defined by 1
~
T(x,, h ; Yo) --- h [~b(x,, h ;Y0) - ~0(x, + h ; x,, Y0)].
(17)
Such a method then is said to be consistent with equation (4) if the limit lira T(xt, h ;Y0) = 0
h~0
holds uniformly for Y0 in an arbitrary bounded set, 0 ~< xt ~< a, and x~ away from material interfaces ; i.e. if, given positive constants e and b, there exists h0 >0, possibly depending on ~ and b (and/~), but not on xt or Y0 (for the given fixed piecewise continuous ~ and S), such that if [Y0] < b,0 < h < h0, 0 ~< x~ ~< a - h , and there is no discontinuity of a or S lying strictly between xl and xt+h, then [T(xt, h;y0)[ < e. The method is
Closed linear one-cell functional methods
183
said to be consistent of order n with equation (4) (or simply order n) if T(xl, h ; yo)/h" is similarly uniformly bounded as h ~ 0. A proof of the following algebraic characterization of consistent LOF methods can be given by using the easily computed analytic form of q/(x ; xt, Y0). Theorem 1 (Nelson and Zelazny, 1986). Suppose o and S are piecewise continuous on [0, a]. Then an L O F method (5) is consistent with equation (4) if, and only if, 1
lim
h~0+ h
P
m-:---~,bm (O(X,+), h,/0Lm (S I C; x,, h) = S(xt+ )/# =
and lim h~ 0+
r(a(x,+), h, ~ ) -
1
h
= - tr(xl+)/It
hold uniformly for 0 ~< xt ~< a, xt away from material interfaces, and # bounded away from zero. In the same m a n n e r we define consistency of a closure approximation (7) as follows: the associated Lm
truncation error over C is 1
Tm(x,,h; yo) = ~[£m(S[ C ;x,,h; y o ) - Lm(~k(" ;x/,yo)l C ;x,,h)], and the closure approximation is consistent with equation (4) if lim T,,(xt, h;yo) = 0,
h~0
again uniformly for bounded Y0, x~ s [0, a], and x~ away from material interfaces. Likewise the closure approximation is consistent of order n with equation (4) (or simply order n) if T(xt, h ;yo)/h" is similarly uniformly bounded as h ~ 0. In the present generality we do not have available an analytic form of L,,(- ; xl, h) to use in establishing an algebraic characterization of consistency for closure approximations analogous to that provided by Theorem 1 for LOF methods. However, all that is really needed to establish such a characterization is a suitable assumption on the behavior of Lm(" ; Xl, h) for small h, as we demonstrate in our next result. For purposes of stating this result, let us introduce the notation BM = { f : f is a real-valued continuous function defined on [0, a] with piecewise continuous first and second derivatives there, such that If<'~(x)l ~< M for i = 0, 1, 2 and x ~ [0, a] such that f ' ( x ) and i f ( x ) exist}, where M is an arbitrary positive number. Theorem 2. Suppose that for any continuous and twice piecewise continuously differentiablef defined on [0, a], there are constants a0, and at such that
Lm ( f l C; x,, h) = a o f (x,+ ) + a, f" (x,)h + o(h), uniformly in f belonging to any BM and xl~ [0, a] lying away from discontinuities o f f ' or f". Then a closure approximation to L,,, of the form (7), is consistent with equation (4), for piecewise continuously differentiable S if, and only if, the limits 1
lim
P
S~ b,,k (a(x~+), h, #)Lk (SI C ; x/, h) = a, S(x~÷)/#
h - 0 + ~ k~,=
(18a)
and lim
h~0+
Zm(a(xt+ ), h, la) - a o h
=
- a oa (xt+)/#
(18b)
184
P. NELSON hold uniformly for x~ in [0, a] and away from material interfaces (i.e. discontinuities of a or S), and [p[ bounded away from zero.
Proof From the assumed form for L,,, we compute
l{
E
Tm(x~,h;yo)=~ [zm(~(xt*),h,l~)-ao]yo+
~ b,,k(a(xr),h,~)Lk(SIC;x~,h)-aiS(xt*)h/#
k=l
] } +o(h) ,
(19)
where for any fixed S the asymptotic estimate is uniform in xz~ [0, a] lying away from material interfaces, [#l bounded away from zero, and Yo lying in any bounded set. (The restriction of y0 to some bounded set, along with piecewise continuous differentiability of S and I~1 being bounded away from zero, insures that the relevant O(';xt, y0) belong to some B~ not depending on Yo or #.) If the subject closure approximation is consistent, then by taking Yo = 0 we see that equation (18a) must hold, in the uniform sense described. Similarly equation (18b) follows, in the appropriate sense, by choosing any other value of Y0. Conversely, if equations (18a,b) hold in the uniform sense described in the statement of the theorem, then equation (19) clearly shows that T,~(xt, h;yo) ~ 0 as h --. 0, in the uniform sense required for consistency. This completes the proof of Theorem 2. A closed LOF method is consistent if the underlying LOF method and all the associated closure approximations are consistent. A consistent closed LOF method has closed order n if (i) each closure approximation (7) has order nm, where (ii) T(xz, h;yo)/h n is bounded as h - ~ 0 , uniformly for Y0 in an arbitrary bounded set and x~e [0, a] lying away from material interfaces, with T(xt, h, yo) defined as was the truncation error T(xt, h, yo) in equation (17), except with each Lm(SI C; xz, h) in the expression (5) for ff(xz, h ;Y0) replaced by Lm(S [C; Xl, h) + O( h"-). Obviously, it is desirable for a closed LOF method to have its closed order no less than the order of the underlying LOF method. Our next result shows that this indeed holds for most of the specific methods introduced in the preceding section. Theorem 3. If ~ is positive everywhere, then the moment characteristic and continuous methods on p basic linear functionals have both their order and closed order equal to 2p.
Proof We refer to Nelson (1987) or to Victory and Ganguly (1986) for proofs of the assertions that the underlying LOF methods have order 2p. Granted this, the argument that they retain these orders as closed LOF methods is as follows. For the moment characteristic methods, the moment equations (14), along with the order of the underlying LOF method, as just described, imply that the closure approximation to Lm is consistent of order 2 p - m . (The hypothesis that tr is positive, along with the earlier assumption that it is piecewise constant, serves to bound tr above zero.) Now note that Taylor's theorem applied to x' ~ exp [-(T(x r--XZ)/lU], along with maxlpm(x)l = O(1)ash ~ 0 , x~C
shows that ceXp [ - O'(Xr- -
Xz)/IA]Pm-
1(x/) dx"
=
O(h") as h ~ 0.
It then follows that if S,,_1 in equation (13) is replaced by S,,-1 =Sm l+O(h 2p "), then the resulting approximation to ~(x,) is also accurate within O(h2p). But this is precisely the statement that the method has closed order 2p. For the continuous methods, the computation of the orders of the closure approximations is exactly the same. In Victory and Ganguly 0986) and Nelson (1987) it is shown that the basic LOF approximation (5) for the continuous methods can be written as
Closed linear one-cell functional methods ~ Pe(e;r) h'-' ~=Ot~+;~__0Sm(-1)"
185
te'm~P'(e;m)]e -(~+')
(2m)!Fp t e . m a - n m! [_ . . . . . . . . .
Qp(e;m)]
'
where e is as in equation (16a), and P,,(e ; m)/Q,,(e; m) is the mth diagonal Pad6 approximant (Baker, 1975) to exp ( - e ) . As such a Pad6 approximant has a Taylor's series, about e = 0, that agrees with that of exp ( - e ) to the first 2m+ 1 terms, and further each Qm(e;m) is bounded as h ~ 0, it follows that the coefficient of S,, in this equality is O(h ~+ ~). Therefore, if S,, is replaced by S,, = Sm+ O(h zp-m- ~), then the resulting approximation to ~r is of order 2p. This shows that the continuous methods have closed order 2p, and completes our proof of Theorem 3. We remark that the calculation in Nelson (1987) of the orders of the moment-characteristic and continuous methods as LOF methods is closely based upon the techniques of Victory and Ganguly (1986) ; however, the latter actually prove a somewhat different (and stronger) result, as will be elaborated in the next section. The following example shows it is possible for a (reasonable) closed LOF method to have its closed order less than the order of the associated LOF method. See also Neta and Victory (1982) for a discussion of this method from a somewhat different point of view.
Example F; the quadratic method of Gopinath et al. (1980) For this method the underlying LOF approximation is the same as that for the quadratic moment-characteristic method, which is to say the case p = 3 of the moment-characterstic methods. Therefore its order as an LOF method is 2p = 6. The closure approximation to L t(~k ] C;x~, h) also is the same as for the quadratic moment-characteristic method, namely that defined by formulas (12)-(14). However, for m = 2, 3, formula (14) is replaced by /22 : = ( ~ r - @ ) / 2
and
//,3
'--
I~r--[~l iLl. 2
Here ~r and £~ are defined by formulas (12)-(14), and the resulting expressions for £2 and £3 are readily seen indeed to have the form (7) of closure approximations. These closure approximations are easily seen to have respective orders 3 and 2. The corresponding values of T, as in part (ii) of the definition of 'closed order' then are both seen, as in proof of Theorem 3, to be O(h 4) as h -~ 0. Thus, this method has closed order of at least 4. Computational experiments by Larsen and Miller (1980) indicate that there are problems for which its (closed) order is exactly 4, and thus less than the order of the underlying LOF method. To this point we have not used the factorization b,,(~, h, p)L,, (S [ C; x~, h) appearing in the approximate source term for an LOF method (5), nor the corresponding factorization in the closure approximations (7). That is, for all considerations up until this point these terms could equally well have been written in unfactored form, say/~,,,(S[ C;x~,h,g, #). We now describe some considerations that make essential use of this factorization. A minor reason for writing these terms in factored form is simply the observation that for all actual methods presently used or proposed, the dependence on ¢rand p splits offas such factors. However, note that dependence on h appears in both factors. The major reason for such factored form is to capture the asymptotic dependence upon h in the {b,,}, and thus isolate it from the functional dependence upon S. Toward that end, let us define a family {L,,(- ; xt, h)} of basic linear functionals to be normalized if, for each m and xte [0, a],
~= o(1),
IIL,.(';xt,h)ll [ ¢ o(1)
as h ~ 0 .
Here and henceforth the norm of such a basic linear functional is that as a linear functional on the space of continuous real-valued functions defined on [0, a]. Note that the families of basic linear functionals for the methods of Section 3 are normalized. This usually is the case, for the natural definitions of the basic linear functionals. Further, it always can be arranged, except possibly in the 'degenerate' case that some Lm(" ; xt, h) has norm zero for all sufficiently small positive h. (This normalization may lead to dependence upon the 'reference points' x~ in the b,,, but this would not affect any of our results.)
186
P. NELSON
*o (=
o) ~'1
x,t (= a) Xll
~./-1
Xi
Xj+I
•17J- 1
ZJ-2
Cj-1
C,
[
Cj
hj-l ~hj---'~
~,-hx "t"
Fig. 2. A mesh of spatial cells, and the associated standard notation. Consider now a consistent LOF method on a normalized family of basic linear functionals. Theorem 1 insures that P ; xt, h) = O(h) as h ~ 0 , uniformly in xte[0, a] lying away from ~=tbm(a(xz+),h,#)Lm(" material interfaces and # bounded away from zero.
(20)
In the ensuing we will wish to use this asymptotic estimate to infer the consequence that
b,,(a(xt+), h,/~) = o(1) as h ~ 0, uniformly in xze [0, a] lying away from material interfaces and p bounded away from zero, for m = 1. . . . ,p. (21) In general this inference is incorrect, as can be seen by easily constructed examples. However, let us define a normalized family {L,,(- ; xt, h)} of basic linear functionals to be asymptotically linearly independent if for any xt~ [0, a] the asymptotic estimate (20), with O(h) replaced by o(1), implies formula (21). If an LOF method has asymptotically linearly independent basic linear functionals, then the desired inference clearly does hold, along with the analogous result for closure approximations (7). We record this as a formal result, for ease of later reference. Lemma 1. Let formula (5) be an LOF method such that the associated family of basic linear functionals is asymptotically linearly independent. If this method is consistent with equations (3), then the coefficients {b,,} satisfy formula (21). If further the basic linear functionals satisfy the hypotheses of Theorem 2, and the associated closure approximations (7) are consistent with equation (4) (so that the closed LOF method is likewise consistent), then the coefficients {bmm,} satisfy bmm,(a(x~+),h, p ) ~ 0 as h ~ O, uniformly in x~E [0, a] lying away from material interfaces and/~ bounded away from zero, for 1 ~< m, m' ~
ARE ASYMPTOTICALLY FINITE-DIFFERENCE METHODS
Closure approximations are defined precisely so that one can carry out an analog of the inner iterations (3) within a discrete spatial approximation. We now wish to consider some of the details of such a discrete inner interation. Thus we must consider a mesh of cells, say ~ = {Ci}, where Cj = [xj, xj+ d, with the remaining notation as shown in Fig. 2. By contrast, the previous discrete considerations have been limited to a single cell. The discrete version of ray tracing consists simply of the explicit (in j ) difference equations P
> O,
(22a)
~O,+l) * Fid+~(n+t)-I-1 ,5 ..F("!. ij = ~ij' ~ ~ ~mq ~mo, ]~i < O,
(22b)
~n,j++]) = .Cij~};+ l, ~_ E b , . o £ ~ ,
Iz,
m=l
m=l
ff}g+o=f
foru,>0,
~}7+')=f
for/ti<0,
(22c)
where z0 : = z(aj, hi, #~),
(23a)
Closed linear one-cell functional methods
187
bmi~ : = b~(aj, hi, #3
(23b)
and aj : = a(x),
xeC ° := (xj,xj+,).
(23c)
E("! Here "~/7~" r U + ° is to be interpreted as an approximation to O("+ °(xj, #3, and --mu as an approximation to L,,(S}")ICj;xj, hj). Equations (22) are the LOF equations (6a), or their counterpart for #i < 0, in application to the particular cell Cj. In consonance with equation (6b) and the interpretation of £~]j, the latter are given in terms of £ ~ j ~ L,.(01") I Cj ; xj, hi) by N
£ ~,~?,j= ~ wr k~,7£~],j + L,,, (q~ ] C ; x j, hi),
(24a)
i'=1
where
kirj : = kii,(x),
x e C °.
(24b)
It is convenient to consider the t~("~ as the beginning of an iteration. If these are known, the E{"!. -,,u can be --m U computed from formulas (24), and then the ",IT(,+ i) from equations (22), subject to the known values of rq ~ ( n + l) i0 = ~O(0,/~i) = f and ws ,FA-+1) = ~(a, #i) = f , for #i > 0 and #t < 0, respectively. The -£(,.+ - m U o, as needed to close the iterative loop, are computed by the closure approximations (7). In application to equation (3a), over the cell Cj, these become p
(,+ 1) _-- ~ ~e -~ij'r,F,(,+ q 1)+ ~ ~hm m ... ij -E~"?.. - m lj,
Pi > O,
(25a)
]~i <
(25b)
m'=l
and p
£(,+ m i j 1) -_- . ~ m i j ' ~,r,(,+ ' ij+ o I al-
E m'=l
h
,.. E("?..
v m m tj - - m U'
O,
where L,,, (S}") I Cj ; xj, hi) has been approximated by _,,,,j,t.(")..as is consistent with the interpretation of the latter. Here we have also introduced an obvious extension of the notation (23). In carrying out this iterative scheme, one presumably would compute the appropriate approximate flux functionals, via the closure approximations (25), more or less simultaneously with the discrete ray-tracing calculations (22). The required values of the approximate source functionals would be computed from formulas (24), as needed. This requires that the p N J values of the approximate flux functionals be stored at the end of each iteration, in order to be available as needed within the next iteration. Thus the number of basic linear functionals, p, is also the amount of storage, per cell and discrete direction, required by a closed LOF method. (For some methods it will be necessary to correct this to account for the fact that a given approximate source functional is used in each of two adjacent cells. In practice the iterative procedure just described is often modified so that the approximate source functionals, {/2~}, are stored at the end of each iteration, for use in the next. The advantage of this is that often these are only weakly dependent on the direction index (t), so that savings in storage might be possible through some type of data compression in this index. For an extreme example, if the scattering and sources are isotropic, then the {/2~]j} are independent of i; more generally, a few ( < N ) discrete Legendre coefficients may suffice. For monoenergetic problems this is not particularly useful as, in any case, the approximate flux functionals, {£~]j} must be temporarily stored during the ray-tracing calculations (22) and (25), in order to permit ultimate computation of the approximate source functionals via formulas (24). However, for multigroup problems this strategy can, particularly for scattering that is only slightly anisotropic, yield a significant reduction in storage requirements. Regardless of such considerations, formulas (22), (24) and (25), along with the boundary conditions inherited from (2), do define an iterative scheme. For certain particular closed LOF methods, including the diamonddifference and step-characteristic methods, conditions are known (Menon and Sahni, 1985 ; Nelson, 1986) that assure convergence of this scheme. If we assume this convergence, then the limiting values (Lmij = niim £~],J ~o etc.) satisfy
188
P. NELSON
~i,j+l = "CiJ~ijq- ~ bmijEmij,
~i • O,
(26a)
rn=l
~i,j = rij~i,j+l J[- ~ b,~u £mij,
]~i < O,
(26b)
m=l N
£,~j = ~ wi,k~rj£,~rj + Lm(qi l C j ; xj, hj),
(26c)
i'=l P
£,,u = z,,~:~u+ ~ bm,-'u/Sm'u, ]~i > 0,
(26d)
rW--I
£mij = Tmij~i,j+ I "JU~ bmm,ijff~m,ij, m'--I
]2 i
(26e)
< O,
ff, o = f ,
/t~ > 0,
(26f1
OiJ=f,
~,<0.
(26g)
and
Regardless of the question of convergence of formulas (22), (24) and (25), one can consider the question of whether approximation (26) has a unique solution and, if so, the manner (if at all) in which q7u approximates 0(xj, #~), £mU approximates Lm(S~ICj;xj, hi) etc. [In fact we could equally well have used the basic LOF and closure approximations to write down equations (26) directly as an approximation to problem (1, 2).] This is the setting in which Victory and Ganguly (1986) carried out their fundamental studies of consistency, along with convergence and stability, for the moment-characteristic and continuous methods. In this sense, the results of the cited work also include that part of Theorem 3 dealing with the closed orders of these methods. The discrete approximation (26), as defined by a closed LOF method, involves three sets of discrete dependent variables, namely the approximate cell-edge fluxes {ffij}, the approximate source functionals {Zmij} a n d the approximate flux functionals {£mU}" In general it does not seem possible to eliminate any of these, in terms of the others, and thereby obtain a system of reduced size. However, for a consistent LOF method, a sufficiently small mesh, and under further reasonable assumptions, it turns out that any two of these sets can be eliminated, thereby leaving a system of equations with the remaining set as the only unknowns. The purpose of the remainder of this section is to provide precise formulations and proofs of these (three) assertions. We begin with : Lemma 2. Suppose that the closed LOF method underlying equations (26) has asymptotically linearly independent basic linear functionals, that the associated closure approximations are consistent with equation (4), and further that the mesh ~ has each material interface (i.e. discontinuity of the piecewise constant functions k or a) coinciding with some cell edge. Then there exists some e such that if II~ II : = l<~j<~J max hi~ < e,
(27)
then the system consisting of equations (26c-e), for m and j ranging between 1 and p or 1 and N, respectively, and fixed j subject to 1 ~< j ~< J, has a unique solution for either £mU, or £,,ij, in terms of the ~ej having #~ > 0 and the ff*.a+l having/z, > 0.
Proof. By substituting equation (26c) into equations (26d,e), we obtain equations of the form p
N
£m,j = Z
E
m'=l
i'=1
amim,,,jLm,,,j+Bm,j,
where the flmU involve the q7u and qTu+ ~, as described in the statement of the theorem, and otherwise known data. By Gerschgorin's theorem, it suffices to show that there exists e such that II # IP < ~ implies
Closed linear one-cell functional methods
189
N
• IA..,m',vl < 1, m'=l
(28)
i'--I
for e a c h j = 1 , . . . , J. But if#i > 0, then the 1.h.s. of inequality (28) is equal to p
N
Z
~ Ibmm,ullwrllk~,jl,
m'=l
i'=l
where bmm,ij = b m m , ( ~ j , h j , ]~i)"
As the w~, and k~rj are uniformly bounded, over all possible combinations of subscripts and meshes, existence of e as desired follows from Lemma 1, along with the requirement that each material interface be a cell edge. For #i < 0 the required estimate (28) is obtained similarly. This completes the proof for the {£mU}' That for the {/Sin0.is similar, which completes our outline of a proof of Lemma 2. If the solution for the {£mU},obtained as described in Lemma 2, is substituted into equations (26a,b), then there results a system of equations in only the approximate cell-edge fluxes, {qTu}. This system is equivalent to system (26), in the sense that the corresponding {/2,,u} and {£mU} can be obtained as described in Lemma 2. This is the essence of the proof of the result encapsulated in the title of this section. A precise detailed statement of this result is : Theorem 4. Suppose the hypotheses of Lemma 2 hold. Then there exists e > 0 such that if expression (27) holds, then there exist (mesh-dependent) coefficients B~evand ~u, 1 ~< i, i' ~< N, 1 ~< j ~< J, such that : (i) if the system (26) has a solution for {~u}, {£,,u} and {£mu}, 1 ~< i ~< N, 1 ~
~,.j+ , = ~,j~,j + E B..j~,j + E B.,j~.j+ , + ~,j, #~ > O, i'
i" + i"
ff~0=f
for#i>0,
i"
tffij=f
l~i<.N,
forp~<0,
l<.j<~J,
(29)
where ~ ( ~ , ) indicates summation over only those i' such that/~r > 0 ( < 0) ; (ii) conversely, if system (29) has a solution for {flu}, for 1 ~< i ~< N and each j, 1 ~< j ~< J, and {/2mu}, {£mU} are calculated from these {flu} as described in Lemma 2, then the resulting values of {qTu}, {£~u} and {£,.u} comprise a solution of the system (26), for 1 ~< i ~< N, 1 ~< j ~< J and 1 ~< m ~
lffq =
k= 1
l
1 "~il
m= 1
bmik£,nik
+
Zi~ f ,
/~, > 0,
(30a)
and lifo =
bmikZmik +
"~il k =j L \ l=j
1
rit f ,
#, < 0.
(30b)
=j
We display this solution explicily, because it has some inherent interest as the discrete counterpart of the explicit solution of equation (4), with kt = #i, S given, and subject to the boundary conditions (2a,b) By means of equations (30), the system for the {£,,u} in terms of the {qTu}, as described in Lemma 2, can be converted into a system involving only the {/2mu}. Likewise, by using equation (26c), one arrives at a system involving only the {/2mU}as unknowns. We formally record the latter result.
190
P. NELSON Theorem 5. Suppose the hypotheses of Lemma 2 hold. Then there exists e > 0 such that if
formula (27) holds, then there exist coefficients C,~,.,ii,jj,and flmis, 1 <~ m, m' <<.p, 1 <~ i, i' <~ N, 1 ~< j, j ' ~< J, depending upon the mesh and the data for the boundary-value problem (1, 2) such that : (i) if the system (26) has a solution for {~o}, {£mij} and {£m~j}, where the indices range as above, then the {£mlj} satisfy p
£mij = E m'=l
N
J
E ECmm'ii'ii'jj'£m'i'J'~-flmij , i'=l j'-I
l<.m<~p,
l<~i<~N,
l<~j<~J;
(31)
(ii) conversely, if equation (31) has a solution for the {£m~S}, and {£mO}, {qT~j}are calculated by equations (26c) and (26a,b) [or equations (30)], respectively, then the resulting values comprise a solution of system (26), for all relevant values of the indices. There is, of course, a similar result for the approximate source functionals, {£,,o}. None of these (asymptotically possible) reduced forms of the system (26) have been exploited in studies of specific closed LOF methods. Rather the standard approach has been to consider the system in the {qTtj} and {£m~j} that results from eliminating the {£mO} by equation (26c). For example, in Nelson (1987) an explicit solution is given for the system in the case of monodirectional (N = 1) continuous methods. This is how the explicit expression for ~r for these methods, as given in the preceding section, was obtained. Victory and Ganguly (1986) even obtain an analytic solution of this system for continuous methods and a general value of N. This solution is a key to their study of convergence, stability and discretization error for the continuous methods.
6. C O N C L U D I N G R E M A R K S
Our primary objective in introducing the class of closed LOF methods is to provide a suitable framework for a general theory of spatial approximations to the monoenergetic discrete-ordinate equations in slab geometry, and possibly even eventually in general multidimensional geometries. An essential step in developing such a theory, and thereby filling out this framework, is to show that a consistent closed LOF method is both stable and convergent to the solution of the underlying continuous problem. This would generalize previous results for a number of specific closed LOF methods. One possible two-step approach to obtaining this desideratum would be as follows. First, extend the basic result of Nelson and Zelazny (1986) to show that a consistent closed LOF method is both stable and convergent for the monoenergetic discrete-ordinate equations subject to initial conditions. Secondly, extend to closed LOF methods the result of Keller and White (1975) to the effect that, for a linear system of differential equations, stability and convergence of a given finitedifference approximation for linear boundary conditions is equivalent to these same properties for initial conditions, provided both the initial and boundary problems are well-posed. A second desirable objective is that of showing convergence of the discrete inner iteration process (22)-(25), for a consistent closed LOF method and sufficiently small mesh. Existing methods (Menon and Sahni, 1985 ; Nelson, 1986) of proof seem essentially limited to methods of relatively low order. The collectively compact operator theory of Anselone (1971) provides a possible point of attack in attempting to obtain such a result. This theory requires an associated 'injection' of the discrete approximate angular fluxes back into the space of continuous angular flux; i.e. each discrete approximation to the angular flux must be associated to a continuous approximation. For the moment-characteristic and continuous methods there is a natural definition of such continuous approximation, although the desired convergence apparently has not even been established for these methods. For general abstract closed LOF methods there is no such 'natural' definition of an associated continuous approximation. Our Theorem 4 suggests the possibility of simply using linear interpolation, although obviously this will not provide an extremely accurate continuous approximation. Finally, it would seem that the clarification of the requirements of a discrete spatial approximation, as provided by the explicit description of closed LOF methods, could provide a stimulus for the development of novel specific methods. Most methods that have been used or proposed are Galerkin methods (i.e. the continuous methods) or 'improved' Galerkin methods (i.e. the moment-characteristic methods), or otherwise involve piecewise polynomial approximations to the source function (S), so that the basic linear functionals are polynomial moments. However, Runge-Kutta methods essentially are closed LOF methods in which the basic linear functionals are point evaluations. These do not seem to have been used in transport theory. Also,
Closed linear one-cell functional methods
191
there would a p p e a r to be some possible advantages (Nelson a n d Zelazny, 1986) to using (closed) linear multicell methods, in analogy to the linear multistep m e t h o d s (e.g. L a m b e r t , 1973) t h a t are widely used for numerical solution of general initial-value probems. Acknowledgemen~This research was partially supported by the U.S. Department of Energy, under Contract No. DEAA03-76-SF00767. REFERENCES
Alcouffe R. E., Larsen E. W., Miller W. F. Jr and Wienke B. R. (1979) Nucl. Sci. Engng 71, 111. Anselone P. M. (1971) Collectively Compact Operator Theory. Prentice-Hall, Englewood Cliffs, N.J. Baker G. A. Jr (1975) Essentials ofPadb Approximants. Academic Press, New York. Gopinath D. V., Natarajan A. and Sundararaman V. (1980) Nucl. Sci. Engng 73, 181. Hennart J. P. (1981) Ann. nucl. Energy 8, 677. Hill T. R. (1975) Report LA-5990-MS. Keller H. B. and White A. B. Jr (1975) S l A M Jl numer. Analysis 12, 791. Lambert J. D. (1973) Computational Methods in Ordinary Differential Equations. Wiley, London. Larsen E. W. and Miller W. F. Jr (1980) Nucl. Sci. Engng 73, 76. Lathrop K. D. (1969) J. comput. Phys. 4, 475. Lee S. M. and Vaidyanathan R. (1980) Nucl. Sci. Engng 76, 1. Menon S. V. G. and Sahni D. C. (1985) Trans. theor, statist. Phys. 14, 353. Nelson P. (1971) J. math. Analysis Applic. 35, 90. Nelson P. (1973) S I A M J1 numer. Analysis" 10, 175. Nelson P. (1986) S l A M Jl numer. Analysis. 23, 1017. Nelson P. (1987) In preparation. Nelson P. and Zelazny R. (1986) Nucl. Sci. Engng 93, 283. Neta B. and Victory H. D. Jr (1982) Numer. Funct. Analysis" Optmizn 5, 85. Neta B. and Victory H. D. Jr (1983) S l A M Jl numer. Analysis 20, 94. Vaidyanathan R. (1977a) Atomkernenergie 29, 30l. Vaidyanathan R. (1977b) Proc. 5th Int. Conf. Reactor Shielding, p. 754. Science Press, Princeton, N.J. Vaidyanathan R. (1979) Nucl. Sci. Engng 71, 46. Victory H. D. Jr and Ganguly R. (1986) S I A M J. numer. Analysis 23, 78.
ANE
14:4-D