LEAST-SQUARES FINITE-ELEMENT SOLUTION OF THE
NEUTRON TRANSPORT EQUATION IN
DIFFUSIVE REGIMES
by
KLAUS JURGEN RESSEL
B. S. (Math), Universitat zu Koln, 1985
B. S. (Physics), Universitat zu Koln, 1986
M. S. (Math), Universitat zu Koln, 1991
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
1994
This thesis for the Doctor of Philosophy
degree by
Klaus Jurgen Ressel
has been approved for the
Graduate School
by
(? //h'fa**-/ /33
Date
Ressel, Klaus Jurgen (Ph.D., Applied Mathematics)
Least-Squares Finite-Element Solution of the Neutron Transport Equation in Diffusive Regimes
Thesis directed by Professor Thomas A. Manteuffel
ABSTRACT
A systematic solution approach for the neutron transport equation is considered that is
based on a Least-Squares variational formulation and includes theory for the existence and
uniqueness of the analytical as well as for the discrete solution, bounds for the discretization
error and guidance for the development of an efficient solver for the resulting discrete system.
In particular, the solution of the transport equation for diffusive regimes is studied.
In these regimes the transport equation is nearly singular and its solution becomes a solution
of a diffusion equation. Therefore, to guarantee an accurate discrete solution, a discretization
of the transport operator is needed that is at the same time a good approximation of a
diffusion operator in diffusive regimes. Only few discretizations are known that have this
property. Also, a Least-Squares discretization with piecewise linear elements in space fails
to be accurate in diffusive regimes, which is shown by means of an asymptotic expansion.
For this reason a scaling transformation is developed that is applied to the trans-
port operator prior to the discretization in order to increase the weight for the important
components of the solution in the Least-Squares functional. Not only for slab geometry but
also for x-y-z geometry it is proven that the resulting Least-Squares bilinear form is contin-
uous and V-elliptic with constants independent of the total cross section and the scattering
cross section. For a variety of discrete spaces this leads to bounds for the discretization
error that stay also valid in diffusive regimes. Thus, the Least-Squares approach in combi-
nation with the scaling transformation represents a general framework for the construction
of discretizations that are accurate in diffusive regimes.
For the discretization with piecewise linear elements in space a multigrid solver in
space was developed that gives V-cycle convergence rates in the order of 0.1 independent of
the size of the total cross section, so that one full multigrid cycle of this algorithm computes
a solution with an error in the order of the discretization error.
This abstract accurately represents the cgnten;
publication.
Signei
Thomas A. Manteuffel
CONTENTS
Acknowledgements ........................................................ v
List of Notation..................................................... vi
CHAPTER
1 INTRODUCTION AND PRELIMINARIES ....................................... 1
1.1 Introduction and Outline......................................... 1
1.1.1 Opening Remarks.............?.......................... 1
1.1.2 Outline ............................................... 1
1.2 Neutron Transport Equation and Diffusion Limit............... 2
1.2.1 Neutron Transport Equation . .............................. 2
1.2.2 Diffusion Limit.......................................... 6
1.3 Previous Work on Numerical Solution ....................... 7
1.4 Least-Squares Approach .......................................... 9
2 SLAB GEOMETRY....................................................... 13
2.1 Problems with Direct Least-Squares Approach.................. . 14
2.2 Scaling Transformation.......................................... 17
2.3 Error Bounds for NondiffUsive Regimes........................ 20
2.3.1 Continuity and V-ellipticity.............................. 20
2.3.2 Error Bounds.............................................. 27
2.4 Continuity and V-ellipticity with respect to a scaled norm... 30
2.5 Error Bounds for Diffusive Regimes........................... 38
. 3 X-Y-Z GEOMETRY .................................................... 46
3.1 Continuity and V-ellipticity ................................... 47
3.2 Spherical Harmonics ........................................... 51
3.3 Error Bounds .................................................. 54
4 MULTIGRID SOLVER AND NUMERICAL RESULTS............................... 60
4.1 Sjv-Flux and Pjv-i-Moment Equations ............................ 60
4.1.1 S^-Flux Equations......................................... 60
4.1.2 Moment Equations.......................................... 62
4.1.3 Least-Squares Discretization of the Flux and Moment Equations 64
4.2 Properties of the Least-Squares Discretization............... 66
4.3 Multigrid Solver ............................................. 75
4.3.1 Sjy Flux Equations........................................ 75
4.3.2 Moment Equations.......................................... 76
5 CONCLUSIONS........................................................ 81
5.1 Summary of Results.............................................. 81
5.2 Recommendations for Future Work................................. 82
BIBLIOGRAPHY............................................................... 83
APPENDIX
A FLUX STENCIL ......................................................... 87
B MOMENT STENCIL ..................................................... 98
ACKNOWLEDGEMENTS
I wish to thank, first and foremost, my advisor, Prof. Tom Manteuffel, for his
academic and financial support, which made this thesis possible. His ready availability to
discuss problems in deep detail, along with his mathematical insight and creativity resulted
always in helpful answers and hints, which formed the foundation of this thesis. Secondly, I
would like to express my appreciations to Prof. Steve McCormick, who has been a consultant
to this work since its beginning and whose expertise in multilevel algorithms has proved
invaluable. I wish also to thank the remainder of my committee, Professors Jan Mandel,
Jim Morel and Tom Russell. Jan Mandel and Tom Russell taught me in their classes the
mathematical theory of Finite-Elements. From the many discussions with Jim Morel during
my stay of 9 months at Los Alamos National Laboratory I gained a lot of insight into
transport problems. I am also grateful to the Center of Nonlinear Studies at Los Alamos
National Laboratory for the financial support and the use of their facilities during this time.
Particular thanks goes to Debbie Wangerin, who guided me through the jungle of rules
imposed by the Graduate School. Moreover, I wish to express my gratitude to Dr. Gerhard
Starke, who proofread most parts of this thesis and made useful comments. Last, but not
least, I would like to thank all members of the Center for Computational Mathematics of
the University of Colorado at Denver and my friends Marian Brezina, Dr. Max Lemke, Dr.
Jim Otto, Radek Tezaur and Dr. Petr Vanek for their friendship and support.
LIST OF NOTATION
For the most part, the following notational conventions are used in this thesis. Par-
ticular usage and exceptions should be clear from the context.
Scalars, Vectors and Sets
x,y,z standard space coordinates.
Zl,Zr left and right boundary of the slab.
9 polar angle with respect to z-axis.
V = cos(0).
Ot total cross section.
0's scattering cross section.
Oa =
This result can be directly extended to the three-dimensional case (Pomraning [48]).
With scaling (1.10), the three-dimensional transport equation becomes
Â£ip(r,Q)
0 V + -(/ P) + eaP
Â£
ip(r, Â£2) = eq(r, Â£2)
and its solution has the diffusion expansion
(1.14)
Vfc il) =
where
equation
- V- JV
offt
1.3 Previous Work on Numerical Solution
The number of actual neutron transport problems that can be solved in closed an-
alytical form is very small. Some examples are presented .in Wing [51, Chapter 6] and in
Duderstadt and Martin [17, Chapter 2]. Therefore, computational methods are required for
the solution of most neutron transport problems. They fall into two broad classes: deter-
ministic and stochastic. Stochastic methods, of which the Monte Carlo method is a chief
example, involves determination of the neutron distribution via random sampling of a large
number of neutrons in the system. As the number of particles in the simulation is increased,
the statistical accuracy of the resulting solution is improved. Consequently, these methods
are often prohibitively expensive, especially when high accuracy is needed. On the con-
trary, deterministic methods involve a discretization that transforms the neutron transport
equation into a finite system of algebraic equations, which can be solved by computers.
Since in this thesis a new deterministic approach is developed, we restrict the fol-
lowing overview, which is far from complete, to previous deterministic methods and begin
with the slab geometry case. For the discretization of. the angle dependence, most frequently
a discrete ordinates (Sn) method (Carlson and Lathrop [13]) is used. This approximation
assumes that the angular dependence of the solution can be expanded in a finite number
of Legendre Polynomials and a set of discrete equations results form a collocation at Gauss
quadrature points. For slab geometry, this discretization is equivalent to a Pjv-i approxi-
mation, which is a spectral Galerkin discretization, using the first N Legendre Polynomials
as basis functions (Lewis and Miller [34, Appendix D]).
7
Because of the property that the analytical solution of the transport equation is
converging in the diffusion limit to a solution of a diffusion equation in the interior of the slab,
the discretization of the spatial dependence is much more difficult. For accuracy reasons, the
discrete solution must have the same property. Therefore, this requires a discretization of the
transport operator that becomes a good approximation of a diffusion operator in diffusive
regimes. By applying the asymptotic expansion technique introduced in Section 1.2, to
the discrete solution, Larsen, Morel, and Miller [32] analyze the behavior of various special
discretizations in the diffusion limit. In the discrete case, the mesh size h has to be considered
as a second parameter besides the parameter e. Therefore, they define in their work the
following two different limits. If, for a fixed mesh size h, the discretization approximates
a diffusion operator in the limit s > 0, then the discretization is said to have the correct
thick diffusion limit. On the other hand, if the mesh size h varies linearly with e in the limit
Â£ 0 and this limit results in a consistent discretization of a diffusion operator, then the
discretization is said to have the correct intermediate diffusion limit.
Since standard finite difference discretizations, such as upwind differencing, fail to
have a correct thick diffusion limit, special discretizations have been developed that behave
correctly in diffusive regimes. Among them are the Diamond difference scheme and the
difference schemes of Lund-Wilson and of Castor, which, according to the analysis of Larsen,
Morel, and Miller [32], give the correct thick diffusion limit for the cell average flux.
Moreover, Finite-Element discretizations have been applied for spatial discretiza-
tion of the neutron transport equation in different ways. The direct Galerkin approach to
the first order integro-diiferential form (1.6) of the transport equation, as considered first
theoretically by Ukai [50] and numerically by Martin [40], does not have the correct behav-
ior in the diffusion limit, except when special discontinuous finite elements are used. The
use of discontinuous finite elements results, for example, in the Linear Discontinuous (LD)
discretization (Alcouffe et al. [2]) and the Modified Linear Discontinuous (MLD) scheme
(Larsen and Morel [33]). MLD has the additional property that with a suitable fine spatial
mesh, it can resolve boundary.layers at exterior boundaries or at interior boundaries between
media with different material cross sections.
Further, a variety of Ritz variational formulations have been proposed (see Kaplan
and Davis [26] for a summary). They have the self-adjoint second-order even-parity4 form
of the transport equation as its Euler equation and lead, independent of the choice of the
discrete Finite-Element space, to correct diffusion limit discretizations. However, the even-
parity form of the transport equation is only valid for nonvacuum regions and becomes very
tedious for anisotropic scattering or anisotropic sources (Lewis and Miller [34, p. 260]).
For the solution of the discrete system, a simple splitting iteration known as source
iteration or transport sweep has been used in the past. Because of the slow convergence of
this iteration for diffusive regimes (convergence factor 1 0(A,-)), the Diffusion Synthetic
Acceleration (DSA) method (Larsen [29]) was developed, which uses a diffusion approxima-
tion to accelerate the source iteration. By spectral analysis, Faber and Manteuffel [18] have
shown why this method is successful for problems with isotropic scattering. However, for
problems with anisotropic scattering, DSA is less effective.
Moreover, multigrid methods have been employed for the solution of discrete neu-
tron transport problems. For the LD scheme, a multigrid algorithm in space was developed
by Barnett, Morel and Harris [4] which proved to be effective even for highly anisotropic
4 It can be shown [35] by a certain transformation that the even-parity form is closely related to the
Least-Squares formulation considered here.
8
problems. For isotropic problems, this algorithm is competitive with DSA, although it uses
an expensive block-smoothing. The multigrid algorithm in space of ManteufFel, McCormick,
Morel, Oliveira and Yang [36] for isotropic problems, discretized in space by the MLD scheme,
employs a special operator induced interpolation and has been ported very efficiently to a
parallel architecture [37]. For anisotropic problems a technique called multigrid-in-angle was
developed by Morel and Manteuffel [44]. This scheme involves a shifted transport sweep to
attenuate the error in the upper half of the moments, so that the remaining error can be
approximated by the solution of a problem discretized in angle based on only the lower half
of the moments. Recursive application of this procedure leads to an isotropic problem on
the coarsest level, which can be solved by a multigrid method in space.
For higher dimensional problems, the discretization of the angle dependence also
becomes a problem. For problems with isolated sources in a strongly absorbing medium,
anomalies in the flux distribution, called ray effects (Lewis and Miller [34, p. 194]), are likely
to arise in combination with a discrete ordinates (Sjv) discretization. The 5jv discretization
causes a loss of rotational invariance, since this discretization transforms the fully rotational
invariant transport equation into a set of coupled equations that are at most invariant under
few discrete rotations. Thus, an azimuthally uniform flux, for example, is approximated by a
set of 6-functions at discrete angles, which can be very poor if the number of discrete angles
is not sufficiently large. One potential remedy is a Pn discretization, which is a spectral
Galerkin method using spherical harmonics as basis functions. This discretization results in
a fully rotational invariant discrete problem. However, the coupling of the discrete equations
is complicated and the treatment of boundary conditions is less straight forward.
As in the one-dimensional case, for higher dimensions the discretization in space
must have the correct behavior in the diffusion limit in order to obtain accurate discrete
solutions for diffusive regimes. The direct extension of the appropriate one-dimensional
discretizations is complicated, however. Bogers, Larsen and Adams [7] have shown that the
linear discontinuous (LD) finite element discretization on rectangles does not yield a correct
diffusion limit discretization, whereas the MLD discretization does. However, the efficient
solution of the discrete system resulting from the MLD discretization is an open problem.
Applying a similar multigrid algorithm, which was developed by Manteuffel et al. [36] for
the one-dimensional case, would require the extension of the operator induced interpolation
to higher dimensions, which is complicated. Morel et al. are in the process of developing a
method for three space dimensions based on the even-parity form of the transport equation
and using a Pjy discretization in angle.
We conclude that an arsenal of highly specialized computational methods exists,
whose design is adapted for particular transport problems. However, there is lack of a general
systematic solution approach that includes existence theory of the analytic and discrete
solution, error bounds for the error of discretization and guidance for the development of an
efficient solver of the resulting discrete problem. Especially for higher dimensional problems,
such an approach seems to be needed.
1.4 Least-Squares Approach
In this section, we introduce a systematic solution approach to the neutron trans-
port equation that relies on a Least-Squares self-adjoint variational formulation of (1.4), and
we summarize the associated standard Finite-Element theory. The Least-Squares approach
can be considered as a systematic solution approach, since it includes theory for the existence
and uniqueness of the analytical and discrete problem, as well as bounds of the discretization
error for a whole class of Finite-Element spaces. Furthermore, this approach will guide the
9
development of a Multilevel Projection Method (McCormick [43]) for the efficient solution
of the resulting discrete system.
A Least-Squares Finite-Element discretization with piecewise linear basis functions
in space directly applied to (1.4) does not have the correct behavior in the diffusion limit
(see Section 2.1). For this reason, the scaling-transformation [P + t(I P)], with parameter
t Â£ IR+ specified later, is applied to the transport operator prior to the discretization:
C := [P + T(I-P)][Q-V + atI-asP]
= P(Q.-Â¥.) + T(I~P){n-V) + T(Tt(I-P) + ((Tt-as)P (1.17)
= P(Q Y) + r(J P)(Q V) + rat(I -P) + aaP,
where in the last equation (at as) was Substituted by the absorption scattering cross section
aa. In this transformed operator, the Least-Squares variational formulation of (1.4) is given
by
minP(V), with F(ip) := f [\Â£ip(r,Q qa(r,Â£i)]2 dQdr, (1.18)
V J J
n
with qs = Sq. The Hilbert space V with underlying norm J| ||v will be specified later.
A necessary condition for ip Â£ V to be a minimizer of the functional F in (1.18)
is that the first variation (Gateaux derivative) of F vanishes at ip for all admissible v Â£ V,
resulting in the problem: find ip Â£ V such that
a(ip,v) := = //
m s1 n s1
qs Cv dFldr Vu Â£ V.
(1.19)
The essential part of the theory is to show that the symmetric bilinear form a(*, *) in (1.19)
is V-elliptic i.e., there exists a constant Ce> 0 such that, for all v Â£ V,
a(v,v)> Ce |M|v, (1-20)
and continuous-, i.e., there exists a constant Cc> 0 such that, for every u, v Â£ V,
|a(u,t>)| < Cc |||k IHk-
(1.21)
The proof of the continuity is straightforward, but the proof of the V-ellipticity is difficult
and tricky.
Denote the standard inner product and associated norm of L2(72. xS1) by
(
-//
u v* dQdr;
n s1
[[(/.[[ := y/(u, u) Vu,vÂ£ L2(TZ x S1),
where v* is the complex conjugate of v. Using (1.20), (1.21) and the assumption that
Qs(l,^Q Â£ L2(1Z x S1), which ensures that the functional
qs Cv dCldr
71 s1
is bounded (|/(-u)| < Cc^2||g3|| ||t>||v ), then the Lcix-Milgram Lemma (Ciarlet and Lions
[16, p. 29]) can be applied. It follows that problem (1.19) is well posed in the sense that its
solution exists, is unique and depends continuously on the data qs. The latter follows from
Ce\\ip\\v < a(ip,ip) = l(ip) < Cc1/2||g5|| ||V-[k,
10
so
pl/2
MW < ~^\ks\\.
For the Least-Squares Finite-Element discretization of (1.19), the Hilbert space V
is replaced by a finite-dimensional subspace Vh C V, and (1.19) becomes: find fa e Vh
such that
a{fa,vh) = l(vh) Vvh E Vh. (1-22)
The existence and uniqueness of a solution fa Â£ Vh of the discrete problem (1.22) follows
again from the Lemma of Lax-Milgram since, as a finite-dimensional space, Vh is a closed
subspace of the Hilbert space V and is, therefore, also a Hilbert space with respect to the
inner product of V restricted to Vh.
By subtracting (1.22) from (1.19), it follows immediately that the error is orthogonal
to Vh with respect to the bilinear form a(-, ):
- fa,Vh) = 0 Vuj, Â£ Vb. (1-23)
The Cauchy-Schwarz inequality and (1.23) lead directly to Ceas Lemma (Brenner
and Scott [8, p.62]):
ai^-fa^-^^^a^-Vh^-Vk) Vvh Â£Vh or
fc~
\W i>h\\v < \ -pr- minJlV-t'hlk,
V vhÂ£Vh
(1.24)
with the use of the V-ellipticity (1.20) and the continuity (1.21). By (1.24), the problem of
finding an estimate of the error is therefore reduced to estimating min ||^> ^>h||y. These
vhÂ£Vh
kinds of estimates are provided by approximation theory for a wide class of spaces Vh. For
example, when we consider for simplicity only a semi-discretization in space where Vh is
formed by piecewise polynomials of degree r, V = Hm(JZ) x L2(Sl), Vh C V, and the exact
solution is in Hr+1(7l) x L2(Sl), it can be proved (Ciarlet and Lions [16, Theorem 16.2])
that
min \\'ip-Vh\\m,o
vkevh
where h is the maximum mesh size of the triangulation of 12, used and
M|m,o
Y, f J \Dav\2dtldr
1/2
\i>\
fc+i,o :=
. M< n s1
1 1/2
Y j j \DaTj>\2dQ,dr
|a|=fc+l si
denote the standard Sobolev norm and semi-norm (Adams [1]), respectively. Here, we use
the standard multi-index notation
D^v :=
d^v
dxPldyfodzP*
\P\ :=EA-
i=l
for j3 := (Pi,P2,p3).
11
For Vh formed by piecewise polynomials of degree r, the combination of (1-24) and
(1.25) results in the overall error bound
U MW
V
The crucial point here is that we have shown V-ellipticity (1.20) and continuity
(1.21) with respect to a weighted norm with constants Ce and Cc independent of crt and
cr0, so that an error bound similar to (1.26) for a discretization in space and in angle holds
independent of the size of crt and cra. Hence, the Least-Squares Finite-Element discretization
of the scaled transport operator with piecewise polynomials of degree r > 1 as basis functions
yields an accurate discrete solution even in the diffusion limit.
12
CHAPTER 2
SLAB GEOMETRY
The Least-Squares approach is applied in this chapter to the one dimensional (slab
geometry) neutron transport equation (1.6). Throughout this chapter, we assume without
loss of generality the following:
1) The total scattering cross section is constant in space, so o-t(z) = crt. This can be
established by the transformation
Z
f(Tt(s)ds
*'=7T---------- (2-1)
f crt(s) ds
z\
The transport equation then becomes
+ rti1 ~ p) +
with
ZT Z j* Zf
*1 = J at(s) ds,
2| Zl Zl
2) The slab has length T, so (zr z\) = 1. If the transformation (2.1) was already
applied, this is directly fulfilled; otherwise, this can be established with the simple
transformation z" = This changes
o" (zr zi)
3) We impose homogeneous (vacuum) boundary, conditions, so g\{n) = 0 and gr{n) = 0
in (1.6). This can be done in the following way. Define
gi(fi) for (j, > 0
gr(fi) for n<0
Then, clearly, if>t(z,fj.) 6 Jf1([z;, zr]) x L2([1,1]), so that Ctpt is well defined and
we can solve the problem Ci})o = q C^b with homogeneous boundary conditions.
The original solution is then given by ip = ipo + Vv
In this chapter, let D := [zi, zr] x [1,1] and let
i
uvd/idz and ||u|| := \/(u, u)
zi -l
denote the standard inner product and associated norm of L2(D).
2.1 Problems with Direct Least-Squares Approach
In the following, we give an explanation as to why a Least-Squares Finite-Element
discretization applied to (1.6) using piecewise linear basis functions in space does not, in
general, yield a correct diffusion limit discretization. We recall that the Least-Squares vari-
ational formulation of (1.6) is given by
Jr 1
min
in F(>!>), with := J J fi) eq(z, fi) dfidx, (2.2)
Z\ -1
and
f d'v
V := < v(z,fi) G L2(D) : fi-^- G L2(D),v(zi,fi) = 0 for ft > 0, v(zr,ft) = 0 for fi < 0
In (2.2) we used the parameterized form (1.11) of the transport equation, since it is better
suited for a diffusion limit analysis.
For the discretization of (2.2), the minimization of the Least-Squares functional is
restricted to a finite-dimensional subspace Vh C V. Without loss of generality, in the follow-
ing analysis for the discretization in angle we use a Pi approximation, which assumes that
the angle dependence of the solution has an expansion into the first two Legendre Polyno-
mials. One reason for this is that a semi discretization only in angle by a Pi approximation
results in a diffusion equation [14, Section 8.3]. Second, the behavior of the discretization in
diffusive regimes, where according to (1.3) the exact solution is nearly independent of angle,
is analyzed here; thus, a Pi approximation allowing a linear dependence in angle is sufficient.
For the discretization in space, we use piecewise polynomials on a partition Th of the slab.
Altogether, this results in the discrete space
Vh := {vh G C(D) : vh(z,fi) =.0(z) + fii(a;), where 0,i G TPr{Th)\
Vh(zi,fi) = 0 for fi > 0, vh(zr,(j,) = 0 for /i < 0 } , ' '
where lPr(Th) denotes the space of piecewise polynomials.of degree < r on the partition Th
of the slab.
By the asymptotic expansion introduced in Section 1.2, the minimizer of the Least-
Squares functional can be characterized as follows.
Theorem 2.1 (Characterization of Least-Squares minimizer)
Let the Least-Squares functional F and the discrete space Vh be given as defined in (2.2)
and (2.3) respectively. Suppose tj>h G Vh minimizes F restricted to Vh. Suppose further
that e < 1 and that iph has. the asymptotic expansion in e given by
i>h{z,y) = eQ{z) +n\{z), with
= #(z) :=
v-0 z/=0
(2.4)
where rjv,5v G lPr(Th) are independent of parameter e for all v. We then have:
(i) 50{z) = 0.
(ii) %(z) = -h{z).
14
(iii) Let
Uh := {Vo Â£ IPr{Th) : Vo(zl) = ^Io(zr) = 0, rjo fulfills (ii) for some i$i Â£ IPr(Th)} .
Then for all tjq Â£Uh:
J ^vWo + Wo dz =
*T
J
qr)o dz.
Z\
Zl
Proof. We first prove (i). Using expansion (2.4) in (1.11) we have
Zij>h = j [Mo] + m'o + ft + Mi] + 0(e),
and, therefore,
with
F(1>k) = Â£ evFv(iph),
v--2
Zr
F-2(i>h) = | J 6o(z) dz,
Zl
Zr
F-ityh) = ^ J%(z)6o(z) + S0(z)61(z)dz,
(2.5)
and Fv(%l>h) independent of e for v > 0. For e < 1, it is possible to bound Jr('0/!) from above
independent of Â£ by
F(i>h)
since tph minimizes F and 0 Â£ Vh. Therefore, we must have F-2(i>h) ~ 0 and F-i(ip) = 0,
since otherwise F(iph) oo in the limit Â£ > 0, which contradicts (2.6). In combination with
(2.5), we conclude that
tfoCO = 0. .
To prove (ii), by virtue of (i) we can restrict the minimization of F to the space
( OO OO
wh := < wh Â£ C(D) : wh(z,n) = ^Â£"^(z) + //^V<5(z);
v=0
V = 1
Wh(zi,fJ.) = 0 for fi > 0, Wh(zr,fj) = 0 for /z < 0
where tjv(z), 8v(z) Â£ IPr(Fh) are independent of Â£ for all v.
A necessary condition for Â£ Wh to minimize F is that the first variation of F
vanishes at rph for all admissible Wh Â£ Wh, that is,
(Zi>h,Zwh) = (eq, Cwhj VwhEWh and V e > 0. (2.7)
15
For wh G Wh we have
Cwh= \p-io++e[iiT)[+n26[->r 1182 +ar]o\
+s2 [ni2 + fi26'2 + /i63 + arn] + 0(e3).
Therefore, (2.7) is equivalent to
ZT 1
J Js \(vWo + Mo) + (?0<5i + dfidz + eli + e212 + 0(e3)
zi -1
2r 1
H2qS[ + aqrjo dfidz + 0(e3),
zx -1
(2.8)
where
zr 1
hJ J H2 (TlWi +
Zl -1
Zr 1
h J J fJ-2 (V0V2 + M2) + (7o53 + M3) + (Mo + V^i) + (^3f?o + ^3<5i)
Zl -1
+rjWi + V1&2 + od>[T)o + M'i + M2 + aMi]
+H46[6[ + a2rjoTfo dfidz.
Since (2.8) holds for all e > 0 and for all Wh G Wh, in particular for Wh = ij>k, it follows that
ZT 1 Zr
0 = J J f!2 (rjon'o + 2r?D6i + Mi) dfidz = | J (rj'0 + 61) dz.
Zl -1 Zl
Thus,
Finally, we prove (iii). Because of (ii), we can restrict minimization of F to the
space Wh := {wh G Wh : rf'0(z) = 61(2)}. The choice w;t 6 Wh in (2.7) will zero out the
0(1) and 0(e) term in (2.8). Comparing the 0(e2) term on the left-hand and right-hand
sides gives
J g ^ivWi + Mi) + (M2 + 8262) + a (Mo + 77o<5i)]
+ r Mi + 2a2rjoTfo dz
0
dz
(2.9)
16
for all qv 6 Wh with v > 0 and for all 6 G Wh with v > 1. From the choice 6[ = 0, ?j0 = 0,
T)[ = rj[ and S2 = S2, we conclude that
ZT
J (jfi+Sz) dz = 0 => rj[ = -S2. (2.10)
ZJ
Substituting (2.10) into (2.9) results in
J (^i^o + rjo^Cj + ^i6i + Z^VoVo dz = J ^qS[ + 2aqrjo dz.
z\ Z\
Choosing <5i = 0, then integration by parts leads to
%r Zr
j a6irj'0 + 2a2T)Q-qQ dz = 2a J qr]o dz,
Zl Z\
which with (ii) and after division by 2a becomes
zr z r
J \vWo + Wo dz = J qr)o dz. (2.11)
Zl Z\
Because of the choice wj, G Wh, equation (2.11) holds for all 770 G Uh. Q
One major implication of Theorem 2.1 is that, when tj(z) and <5(z) are continuous
piecewise linear functions, (ii) can only be fulfilled if rjo is a linear function. Otherwise,
Si = t?0 is a step function, which would not be continuous. Taking the boundary conditions
into account, it follows that rjo = 0. Therefore, Uh = {0}, so that (iii) is a vacant statement
and does not contradict the fact that rjo = 0 is a solution. Consequently, in the diffusion
limit e 0 the discrete minimizer i>h converges to iph = 0, independent of the choice of the
right hand side q. This shows that the Least-Squares Finite-Element discretization of (1.6)
with linear basis elements in space does not give a correct diffusion limit approximation,
except in the case q = 0. For a different way of proving this result, we refer to (Manteuffel
and Ressel [38]).
On the other hand, for piecewise polynomial basis functions of degree r > 1, con-
dition (ii) does not restrict rjo to a linear function. Therefore, Uh contains also nontrivial
functions, so that (iii) implies that rjo is a Galerkin approximation of the diffusion equation
+ = q. Thus, the Least-Squares Finite-Element discretization with piecewise poly-
nomials of degree r > 1 yields a correct discrete diffusion limit solution. However, numerical
results for a discretization in space by piecewise quadratic basis functions show that applying
a scaling transformation (introduced in the next section) prior to the discretization enhances
the accuracy.
2.2 Scaling Transformation
In this section, we introduce a scaling transformation that is applied to the transport
operator prior to the Least-Squares discretization. This scaling transformation plays a key
role in this thesis, since it guarantees the accuracy of the Least-Squares discretization in
17
diffusive regimes even for simple Finite-Element spaces, such as spaces using continuous
piecewise linear elements in space (see Section 2.5).
To motivate the scaling transformation we introduce the moment representation
of the flux. Let Pi(fj) denote the 1-th Legendre polynomial. The normalized Legendre
polynomials pi(fi) := y/2l + 1 Pi(p) form an orthonormal basis of L2{[ 1,1]):
l
\ JPk(p)pi(p)dfi = 6ki, (2.12)
-l
where Ski denotes Kronecker delta, i.e, Ski 1 for k = l and Ski = 0 otherwise. Assuming
that il>{z,n) 6 L2([1,1]) for all z 6 [z;, zr], then ij) has the following expansion (moment
representation) in angle:
CO
HZP) = Mz)pi(p)< (2-13)
1=0
where the Fourier coefficients i{z), which are called moments in neutron transport theory,
are given by
l
j'l>{z,p)pi(p)dP- (2.14)
-l
We see directly that the projection operator P, defined in (1.7), is a projection onto zeor
moments (Pip = o), the operator Pfi is a projection onto first moments (Ppip = and
the operator (IP) is a projection onto moment one and all other higher moments. With the
concept of the moment representation, the diffusion expansion (1.13) leads to the implication
that, in diffusive regimes, only moment zero and one are the important components of the
solution.
Because of Ceas Lemma (1.24), the solution of the Least-Squares discretization
can be viewed as the best approximation to the exact solution in the discrete space Vh with
respect to the norm \/a(-, ) := < Â£,Â£ >. However, the different terms in the operator Â£,
as defined in (1.11), are unbalanced (there are 0( j), 0(1) and 0(e) terms), so that different
components of the approximation error are weighted differently in \/a(-, ). The leading term
of Â£ is ~(I P), which means that the error in the higher moments is weighted in this norm
very strongly in diffusive regimes (very small e), even though this part is not important
according to the diffusion expansion (1.13). On the contrary, the error in moment zero,
which is the important part in diffusive regimes, is hardly measured in the norm y/d(-, ),
since it is weighted by e.
The basic idea is, therefore, to scale equation (1.11), thus changing the weighting in
the norm used in the Least-Squares discretization to determine the best approximation to the
exact solution in the discrete space. Define for r 6 JR+ the following scaling transformation
and its inverse:
S~P + t(I-P), S-1 =P + ^(I-P). (2.15)
After applying the scaling transformation S from the left and dividing by e, equation (1.11)
becomes
1 ~ 1 (9?/} 7*
Cil> := -SCi> = -Sii-Â£ + -^(I-P)4> + aPi> = qs, (2.16)
18
where q, Sq and
l ^ dip \ dip t dip
= -Pn-g- + J P)^.
e dz Â£ dz Â£ dz
Clearly, choosing r = 0(e) will increase the weight for moment zero and reduce the weights
for the higher moments.
Equation (2.16) can be balanced further by a scaling transformation from the right.
Let the domain of operator C in (2.16) be the Hilbert space V. Then we define the space V
by
V := S'-1 V, so that ^ _
v = S_1v for all Â£V and Sv = v for all v E V. \ J
Scaling (2.16) also from the right results in
^ 1 f)tb 7^
CSS~liP = CSiP = -SiiS^f- + -r{I P)i> + aPiP = qs,
Â£ OZ Â£*
(2.18)
where
-SfiS = -(r t2)(Ph + fiP) +
e Â£
For r = 0(e) we have f E 0(1), so that in (2.18) the derivative of moment zero and one and
the moments themselves are weighted equally. Moreover, we point out that the double-scaled
operator CS can be bounded independent of Â£.
In the Least-Squares context, the additional scaling from the right can be avoided,
since
min {CSiP qs, CSip qs) <=> min {Cip qs,Cip qs), (2-19)
which will simplify the boundary conditions and so also the computations. Further, for slab
geometry, because of transformation (2.1) we may assume without loss of generality that
at and parameter Â£ are constant in space. However, for higher dimensional problems, this
cannot be established, so that Â£ = e(r). For inhomogeneous material, e(r) is in general
discontinuous, so that the scaling parameter r, which was chosen to be 0(e), would be
discontinuous. To perform the scaling would then require to prescribe jump conditions in
the scaled solution v across material interfaces.
Therefore, we use the additional scaling from the right only as motivation for the
choice of the scaling transformation and as a tool in the theory in Section 2.4, where we
exploit the nice form of the double scaled operator (2.18). For another way of motivating
the scaling by way of the moment equations, we refer the reader to Manteuffel and Ressel [38].
As outlined in Section 1.4, a necessary condition for ip E V to be a minimizer of
the Least-Squares functional (2.19) is that the first variation vanishes at ip, which results in
the problem: find ip EV such that
a(ip,v) := (Cip,Cv) = (qs,Cv) Vu G V.
(2.20)
For a discretization of problem (2.20), the bilinear form a(-, ) is restricted to a finite di-
mensional subspace Vh C V. In the remaining of this chapter, we analyze the error of this
discretization for various subspaces Vh.
19
2.3 Error Bounds for Nondiffusive Regimes
In this section we establish bounds in an unsealed norm for the discretization error
of the Least-Squares discretization. However, in this norm it is not possible to prove V-
ellipticity and continuity of the bilinear form (2.20) with constants independent of parameters
e and a. In diffusive regimes, where e is very small, these bounds blow up and are therefore
useless. Nevertheless, the bounds for diffusive regimes that are derived in Section 2.5 are
only valid for e < l/\/3, so that the bounds in this section can be used to cover the range
Â£ > 1/v^-
As outlined in Section 1.4, the first step on the way to bounding the error is
to prove V'-ellipticity and continuity of the bilinear form a(-,-) in some norm. From the
view of standard elliptic boundary value problems, the choice V = ff1([z/) zr]) x L2([ 1,1])
(Adams [1]) with the norm
l?,o
//(Â£) +v3ilMlz
seems natural. However, it is easy to see that the bilinear form a(-, ) cannot be bounded
from below in this norm. Let Vk := \/2 sin(&7rz) B(/x) with
B(ji) :=
/~3~ 6+li
Y 26 S
fÂ£tzJL
V 26 6
for (i Â£ [6,0]
for fi Â£ [0, <5]
otherwise
Then, for all k Â£ IN, we have Vk Â£ i?1([z/,zr]) x L2([1,1]) and ||fc||i,o = (kit)2 + 1. Some
simple calculations show that
, . 1 + r2 (kir)262 2r2 + 2s4a2
a[V, V) < ----5----------- + --------7-----
Choosing 6 = then the bilinear form a(-, ) is bounded for all k while ^lim [|vjb||i,o = -
Thus, there is no lower bound for a(-, ) in the norm 11-1]! 0.
The next obvious choice is
Ml|2:=
dv
lYz
+
V := {u Â£ C(D), v(zi,fi) = 0 for (i > 0, v(zr,fi) = 0 for fj, < 0}.
(2.21)
Closure here is with respect to the norm |j|-[[[, so that V is a Hilbert space.
In the following, we bound the Least-Squares discretization error in norm (2.21) for
various Finite-Element spaces.
2.3.1 Continuity and V-ellipticity
Before we establish continuity and V-ellipticity of the bilinear form a(-, ), we sum-
marize some simple properties of our operators.
20
Lemma 2.2 (Properties of P, S and fi-^)
For all u,v G V, we have:
(i) (Pu,v) = {u,Pv)i and {(I P)u, v) = {u, (I P)v);
P2 P\ and (I P)2 = (I P).
Thus P and (J P) are orthogonal projections;
(ii) {Pu,v) = (Pu,Pv)\ and {(I P)u,v) = ((/ P)u,(I P)v}\
(0 IHI <
-Sv
e
for scaling parameter r = e and e < 1.
(iv) ||H|2> - \W -PMl) ;
(v) (y,v)>0; l
Proof.
(i):
Zf 1 1
(Pu, v) = f f Pu v dfidz f Pu f v d\idz
Z\ -1 Z\ -1
Zf zr 1
= f 2Pu Pv dz f Pv f u d\idz
Z l Zl -1
z r 1
= f f u- Pv d^dz = (u, Pv),
Zl -1
and the second identity follows directly from the first. From the definition of P, it
is obvious that P2 = P and, therefore, (I P)2 = (I P).
(ii) : follows immediately from (i).
(iii) : IHI2 = ||Pv||2 + ||(I P)u||2 < fHI^II2 + Wi1 ~ -PHI2 = Il'^SuH2, since e <1.
(iv) :
Zr 1
= l \\Pv\\2 + 2 J J n2Pv (I P)v dfi dz + MI P)v
Zl -1
21
The mixed term can be bounded by the Holder inequality as follows:
zr 1
J J n2Pv (I-P)vdfidz
Zi -1
Zr 1
< J \Pv\ Jn[ji\(I P)v] dfidz
*i -i
Therefore
|fi(I P)V>|2 d/i dz
1/2
IIHf > l \\Pv\\2 ~ ||P*|| ||/i(J P)*|| + MI P)v
IM-IK/-p)||
2
(v): Applying integration by parts with respect to z, we get
Zr 1 1
Zr 1
J J /i^ -v dfidz = J (J, [v2(zr,n) v2(zi,n)] dji-J Jn^-vdfidz-
Z\ 1 1 Z\ 1
Taking into account the boundary conditions for v E V, it follows that
l
fiJzV'} = \ J fI[v2(Zr>fJ')-v2(zhlJ')] dV
-1
= ^|/fiv2(zr,fi) dfi- J^/iv2(zi,ii) dfij >
0.
It now easily follows from the Cauchy-Schwarz inequality that the bilinear form
a(-, ) is continuous in the norm |||'|||, since for any u,v G V
Ku>)l = ||(Â£u,Â£v)| < l|Â£ll ll^ll
<
1 + T
du
dz
, T + Â£ a ii i
+ M
1 + T
dv
+
t + sza
HI
< Ce lllulll IIIHII.
Here Cc '=
step.
(ir)2 _j_ ^zP_gPj I ^ ancj we used the discrete Holder inequality in the last
22
We prove now V-ellipticity of the bilinear form a(-, ) when a ^ 0.
Lemma 2.3 (V-ellipticity for a 0)
Suppose a^O and let r = e^fa. Then there exists Ce > 0 such that, for all v G V,
a(v,v) > Ce ||M||2,
where Ce := min {^, a, a2}.
Proof. We have
a(v, v) = (jCVyjCrV)
1 dv 2 , dv
e2 + a V~PT,
+ -t\\(I-P)vtf + a>\)Pv\\
(2.22)
The second mixed term can be written using (ii) of Lemma 2.2 as
() = ) - p)
According to (v) of Lemma 2.2, the first term here is always positive and the second term
cancels with the first mixed term in (2.22), so that
a(v,v) =
1 dv 2 . dv
e* + a
+ ^\\(I-P)v\\2 + a2\\Pv\\
> Ce IIM
with Ce := min { 77, a, a2}, which proves the lemma.
For the more difficult task of establishing the V-ellipticity of the bilinear form a(-, )
when a = 0, we need the following Poincare-Friedrichs inequality.
Lemma 2.4 (Poincare-Friedrichs Inequality)
For any v G V, we have
\\flV
<
Proof. We have
f dv(s,n) r
torf-
Zl
<
f dv(s, fi)
~Jtai-0
. Z
(2.23)
23
< <
z
1
zi
z r
I
9v(s,n)
ds
dv(s,n)
ds
ds for fj, > 0
ds for fx < 0
Zt / 2>r
< J dz < (zr ~ Zi)1/2 | J |/i
dv
dz
1/2
dz
Taking into account that assumption (zr z/) = 1 implies
IIHI2 <
we obtain the lemma.
We are now in a position to establish
Lemma 2.5 (V'-ellipticity for a = 0)
Suppose that a = 0 and 0 < e < 1, and let
dv
T =
yf'
2 +
52
^ = 7 + V1 ~ P)/*Â£ + 72 (I-P> + Pv-
Then there exists Ce > 0 such that, for any v G V,
a(v>v) > Ce [|MI|2 i
where Ce = ^54.
Proof. Recall that
1 dv
-Pp-z~
Â£ dz Â£
Because of (i) in Lemma 2.2, we have
a(v,v) = {Â£v,Â£ v)
1 / dv dv\ r2 / dv
2r2 / f)v \
<(' P>> V ~ P + IF {(P~ P)^ V ~ P>)
Analyzing the mixed term and using (i) of Lemma 2.2, we see that
[v-^TzV- F}") = ({I~ p)pf?) = (''I' ) (P"S'p
24
The first term is always positive according to (v) of Lemma 2.2. Consider the following
arithmetic-geometric inequality: for any rj E 1R+ and for any a, be 1R, 2ab < -qa2 + We
can thus bound the second term according to
P^pv)<
dv
P^z
ll-Pwll 2
dv
^Tz
Therefore, the bilinear form a(-, ) can be bounded from below by
a(v, v) >
dv 2 T2 , a dv
^z H2 2 U-F>S
+Â£ll (I ~ PM2 ~ ^\\Pv\\2 .
Defining
6 :=
[|Pu||2
l|2
dv I
T;=Jtei
/Sir'
SO
that q
( T)_ Wit
the above bound thus simplifies to
a(v,v) > Ci
dv
%
+ ^2 |M|2 with
Ci = ^2 {j--Jt,7 + t2^1~7^)
ft ?(Â£>-*>-Â£),
(2.24)
By proper choice of 77, we now need only establish that Ci and C2 are positive. Unfortunately,
for large enough 6, C2 will be negative, so we will need in this case to readjust the terms in
(2.24), which we do by way of the Poincare-Friedrichs inequality.
Case 1: S > yÂ§: From the Poincare-Friedrichs inequality (2.23) and (iii) of Lemma 2.2, we
conclude that
dv
Tz
>\\H\2>[^\\Pv\\-MI-P)v\
Since fi E [1,1], then clearly ||/i(J P)v[| < ||(/ P)u||. Therefore,
||P|| ||ai(J PHI > ^ ||P|| ||(I P)VII > 0,
where the last inequality follows from the assumption 8 > y| > | since
(2.25)
i||Pu||2>||(J-PH|2^35>(l-^)
>-i
25
From (2.25), we get
-J= \\Pv\\ MI pm) 2 > ||Pt,|| II(I P)vI
It then follows that
IHJlf
Thus,
We now use (2.26) to rewrite (2.24) as
a(v, v) >[Ci-j
>|C1 + 5?
dv
1 dz
dv
"Tz
+
2e2
dv
Tz
+ c2 ihi
+ l^.+ C2)|M|2
V26e2
Choosing 7] = -j- and using the fact that 6 < 1 results in
T^ 1 /
+Ci =
26s2
and
= ?(?(1-)+,(s-b))^^>
c-Â£*-?(i(1-tJ'(1+S)) + tJ >yi
since t 1/^2 + fÂ§, so r2 (l + <1.
Poo/s O* X / i^ Plinrvoinrr m ---- 0/1
Case 2: <5 < 1|: Choosing 77 = 24s, for Cl and C2 in (2.24) we obtain that
and
' r2 / 25 \ r2 / 12 25\ r2
C2 e4 ^24 j e4 V1 13 24) 26s4 >
Ci = \ (7 (1 24t2) + r2 (1 7)) > ^ > -
since 24r2 = 7^2 < 1 for Â£ < y/ff
(2.26)
26
Thus, altogether we have
a(v, v) > Ct
2
with
which completes the proof.
From continuity and V-ellipticity of the bilinear form a(-, ) it follows directly from
the Lax-Milgram Lemma (see Section 1.4) that problem (2.20) and all of its discretizations
are well posed. The next step is to obtain discretization error bounds for a variety of discrete
subspaces Vh, which is done in the next subsection.
2.3.2 Error Bounds
As outlined in the introduction, continuity and P-ellipticity of the bilinear form
a(-, ) lead directly to Ceas Lemma (1.24). Therefore, bounding the discretization error
||\ij) VViIII is reduced to the problem of bounding min |||?/> /i|||, which is a problem of
vK&Vh
approximation theory and depends on the choice of the finite-dimensional space Vh. Here
we consider two main classes of discrete spaces Vh.
The first consists of spaces with functions that can be expanded into the first N
normalized Legendre polynomials with respect to the direction angle p and are piecewise
polynomials of degree r in z on a partition Th of the slab [z\, zr\. This choice of the finite
dimensional space Vh corresponds to discretization by a spectral method in angle p and
a Finite-Element discretization in space. In transport theory the spectral discretization
in angle with the first N Legendre Polynomials as basis functions is also called a -P/v-i
discretization in angle.
For any f(z) Â£ Hm([zi,zr]), with 1 < m < r+ 1, let HZf{z) denote the interpolant
of f(z) by piecewise polynomials of degree r > 1 on a partition of [z\, zr\. It then can be
shown (Johnson [25, p. 91]) that
\\f(z) nj(z)\\LH[zitZr]) < Chm Wf^\{z)\\L^[zuzA)
[f(z)-nzf(z)]
< Ch!
Â£2([*.,*r])
1 /(m)(z)
(2.27)
where h is the maximum cell-width of the partition and the constant C is independent of h
and /.
Further, for any g(z, fi) Â£ Hm{[zi,zr]) x H2([ 1,1]) we define
JV--1
n;vÂ£f(z,M) := XI with Mz)
r=o
l
\ j 9(z,p)pi(p)dv
-i
(2.28)
as the truncated expansion of g into the first N normalized Legendre Polynomials pi(fi) (see
(2.12) in Section 2.2). Note that the normalized Legendre-polynomials form an orthogonal
basis of L2([ 1,1]) and are the eigenfunctions of the Sturm-Liouville operator (Gottlieb and
Orszag [21, p. 37])
Â£sPi(n) ~
A. \(i dP'(p)
dp _ ' dp
= 1(1+ l)pi(p).
(2.29)
27
Then the error of the truncated expansion can be bounded as follows
Lemma 2.6 (Truncated expansion into Legendre polynomials)
For r > 0, let g(z,[i) 6 Hr([zj,zT]) x H2([ 1,1]) and let 11// be defined as in (2.28). Then
have:
IIIM <
(ii)
<_______
- '('+1)
(iii) For any m
Vro < r and Vz G [zi, zr];
dmg C dmg
dzm 11jv dzm ~ N s dzm
(2.30)
with C independent of g and N.
Proof.
(i): II// is orthogonal projection with respect to the inner product of L2([1,1]), so
||n/n7||i2([_li;l]) < IMIL2([_u]), and hence ||IIjvff|| < ||ff||.
(ii): By. definition (2.29) and integration by parts, we have
jJrHz) =lf
~ 21(1 + 1) / Cs dz P,(M) dp ~ ^21(1
Therefore,
<
/(/+!)
Cs
+ 1)
dmg
Cs
dmg
dzm
Â£2([-Ul)
(2.31)
Sz
(iii): Since the Legendre Polynomials are an orthogonal basis, from (2.31) we obtain
dmg dmg
- II//
dz
dzn
= 2Â£ |^m)WI2<
i2([-l.l]) l=N
Cs
dmg
dzr
Â£2([-Ml) 1^2
For / > 1 we have j < so that the sum can be bounded by
OO
o - | . OO ^ TO ^ y I 4 4
Therefore
with C = yjÂ§,
which proves the lemma.
dmg <9mff ^ C dmg
<9zm N dzm ~ N Ls dzm
28
Theorem 2.7 (Finite-Element in space, spectral discretization in angle)
Suppose that 7ft is a partition of the slab [27, zr\ with maximum mesh size h. Let V be given
as defined in (2.21). Define
{JV1
vh e c(Dy, vh = J2 (*) e JPtVh) for / = o,..., n -1
1=0
vh(z,,n) = 0 for n > 0,
Vh(zr,fj,) = 0 for n < 0
where 2Pr(7/l) denotes the space of piecewise polynomials of degree < r on the partition 7ft.
Suppose 1 < m < r + 1 and let ij> G Vn (Hm([zi, zT]) x H2([-1,1])) be the solution of (2.20)
and iph Â£ Vh be the solution of (2.20) restricted to Vh. Then
VIII \\Csi>\\h0 + C2h
m1
dmip
dz"
Proof. From Ceas Lemma (1.24), we have
WH-MW <
min
VhÂ£Vh
HIV1 vh
<
|||V> Hivn^lll
Now note that |||?;||| < Q. Therefore, by (i) of Lemma2.6, (2.27) and (2.30), we conclude
III$ n^n2v>||| < I]i> nJvV,llij0 + - n^V,)lli,0
which proves the theorem.
The second main class of finite-dimensional spaces considered here are formed by
functions that are piecewise polynomials in space z as well as in angle fi. This choice
corresponds to a Finite-Element discretization in both space and angle with rectangular
elements.
Suppose that 7ft is a partition of the computational domain D = [27,27.] x [1,1]
into rectangles T = [zj,zj+1] x [//, fJ,v+i] of maximum diameter h. To be able to handle the
boundary conditions properly, we assume in addition that (27,0) and (zr, 0) are nodes of the
triangulatiori 7ft. By 7ft we define the discrete space:
dmip
dzm
yh
vheC(D); *ft|T=
0
(2.32)
Vh{zi,n) = 0 for n> 0,
Vh(zr,fi) = 0 for fi < 0
29
For all v Â£ V, let H/,u Â£ Vh denote an interpolant1 of v with respect to the partition 7*. It
can be proved (Ciarlet [16, Theorem 16.2]) that, for v Â£ V n7Tr^"1(Â£)) the following bound
for the interpolation error holds:
11^ ~ < Chr+1~m \v\r+1, (2.33)
where 0 < m < k + 1 and | \h-+'l(d) is the semi-norm of Hk+1(D) (Adams [1]). Combining
Cea's Lemma and (2.33), we get
Theorem 2.8 (Finite-Elements in space and angle)
Let Vh, h be given as defined above. Suppose ip Â£ V D Hk+1(D) is the solution of (2.20)
and let iph Â£ Vh be the solution of (2.20) restricted to Vh defined in (2.32). Then we have:
\U-M\ < \f^chr\ip\Hr+1(D).
Proof. By Ceas Lemma, we need only to bound |||V< II/1i/,|||- Note that, for all v Â£
V H Hr+1(D), IIMII < Hwllj 0 < ||v||ffl(1J). Thus, using (2.33) with m = 1, it follows that
|||-0 Hft'0111 < chr \i>\Hr+^D),
which proves the theorem.
We point out that the error bounds in Theorem 2.7 and Theorem 2.8 depend on
the ratio In the V-ellipticity bounds in Lemma 2.3 and Lemma 2.5, the scaling
parameter r was chosen to be 0(e). Therefore, when e is small, Ce is 0(a) when a: / 0,
while Ce is 0(1) when a = 0. In addition, for r = O(e), the continuity constant Cc is
O(j-), so we have = O(j-), which blows up for diffusive regimes, where e is very small.
However, numerical results show that the Least-Squares discretization of the scaled transport
equation stays accurate in diffusive regimes. Thus, we conclude that the bounds, derived in
this section, are not sharp enough to reflect the accuracy of the Least-Squares discretization
in diffusive regimes. In order to obtain error bounds that do not blow up in diffusive regimes,
it is essential to prove continuity and F-ellipticity of the bilinear form a(-, ) with constants
independent of parameters e and a. This is done in the next section with respect to a scaled
norm.
2.4 Continuity and V-ellipticity with respect to a scaled norm
In this section, which is the central part of this thesis, we prove continuity and
V-ellipticity of the form a(-, ) in (2.20) with constants independent of parameters e and a.
This is the foundation for the bounds in Section 2.5 of the Least-Squares discretization error
that do not blow up for diffusive regimes. Throughout this section, we assume that
t = e and a < 1.
In order to obtain continuity and V-ellipticity with constants independent of e and a, we
use a scaled norm. To motivate its choice, we look at the double-scaled (from left and right)
1For r > 2, there axe many different interpolants, depending on the choice of the support abscissas
and support ordinates on the rectangle, which are not specified here. For an overview of commonly used
interpolants for rectangles, we refer the reader to Ciarlet [16, p. 129].
30
transport operator (2.18). Let V denote the domain of the single-scaled (only from the left)
transport operator (2.16) and V = S 1V the domain of the double-scaled transport operator
(2.18). Defining
Q := jSfiS = (l-e)(Pfi + fiP) + Â£fil, (2-34)
we see that the norm
2^
+ IHI2
(2.35)
for v Â£ V would be a natural choice for bounding the double scaled bilinear form
a(u,v) := (Â£Su,Â£Sv) . (2.36)
However, because of the reasons mentioned in Section 2.2, it is desired to use the single-
scaled transport operator for the computations. Therefore, using the relation v = 5_1 V, we
derive from (2.35) the following norm for v Â£ V:
II
= jSfiSS- -i dv dz
ii 1 dv -Sfi-K- Â£ OZ 2 +
= 1 dv -Sfi-K- Â£ OZ 2 +
+ lls-M
Pv+ -(I-P)v
Â£
(I-P)v
+ \\Pv\\2 =:
\v
We define
(2.37)
V := {v Â£ C(D); v(zi,fj,) = 0 for fi > 0; v(zrifi) = 0 for (i < 0}, (2.38)
where closure is with respect to the norm |] ||y, so that V is a Hilbert space with respect
to the inner product
{u,v)v := /^sfi^,^Sfi^\ + /^(I-P)u,^(I-P)v\ + (Pu,Pv).
At a first view, the norm ll-llv also seems to be useless for error bounds when e > 0.
However, suppose we could show that the error of discretization is bounded by
\W ~$h\\v ,qs,h,N), with lim qst h, N) - 0.
iV oo, h+ 0
Then, in particular, every part of the norm ||^i i>h\\y is bounded by C, so that ||(7 P)(ip iph)\\ <
eC, for example, which means that the error in the higher moments is decreasing with e.
First, we prove continuity of the bilinear from a(-, ) in the norm H'lly. From the
Cauchy-Schwarz inequality, it follows directly that, for any u,v Â£ V,
|a(u,u)| = |(Â£u,Â£u)| < ||Â£u|| H-Ci'H .
31
Employing the discrete Holder inequality and using the assumption a < 1, we obtain
IIjCuII <
Isc-
e *dz
+
-(I-P)u
Â£
+ PM
1 du 2 1,
-Sfi-K- e dz +
+ ll-Pu;
\ 1/2
I2 =V3
< V3
Thus, for all u, v 6 V,
|a(u,v)| < Cc ||||v |M|V
with Cc = 3.
To prove V-ellipticity of the bilinear from a(-, ), we exploit the convenient form
of the double-scaled transport operator and prove first that the double-scaled bilinear from
a(-, ) in (2.36) is V-elliptic. The V-ellipticity of the bilinear from a(-, ) then follows easily
as in Corollary 2.12. In order to prove V-ellipticity of the bilinear from a(-, ), we need the
following lemmas.
Lemma 2.9
For all u, v 6 V, v Â£ V and e. < 1, we have
(i) ^ll^ll e||(J p)v\\ < IMI;
(ii) {QTz'*) -0;
(iii) (Poincare-Friedrichs inequality)
IIHI <
Proof.
dv 1 dv dv
< ^Yz < -SfiS Â£ dz QY
(2.39)
(i): Since \\fi(I P)w|| < ||(I P)v[| and
ZT 1
Zl -1
l
J J /x2 (Pv)2 dfidz = J (Pv)2 J /i2 dfidz = J (Pv)
1 Zl
dz
z r 1
\j J (Pv)2 dtidz=\\\Pv\\2,
Zl -1
then
iihi = iia*^+- p)]v ii > ii^ii c\w p)* i:
>^npii-eii(/-p)%
(ii): From (i) in Lemma 2.2, it follows that (Su,v) = (u,Sv) Vu, v E V. Therefore,
using (2.17) leads to |
(f?s) = (H5!?') =; (^st) = l ('Â£)
where the last inequality follows from (v) of Lemma 2.2.
32
(iii): From the Poincare-Friedrichs inequality (2.23) proved in Lemma 2.4, we have ||^|| <
llpfell- Using (iii) of Lemma 2.2 and the relation (2.17) results in
dv 1 dv 1 dv dv
"T* < -SV7T e dz -SfiS=- e az
The following technical lemma is tedious but simplifies the proof of the major result,
Theorem 2.11.
Lemma 2.10
Suppose 0 < e < Then, for any 6 Â£ [0,1], there exists 6 > 0 such that
H(b, 5) := | y^6 -(1+ Â£)&)' + 46(1-5)+ (l J) b + 6 j < 0.988,
where s := j^1/2 e\/3(l 6)^2j In particular, for 5 < 0.875, we can choose 6 = 0.
Proof. For 6 < 0.875, we choose 6=0 and get
77(0,6) = i |\/46 362 + 6} < 77(0,0.875) < 0.986,
since 77(0,6) is monotone increasing for 6 Â£ [0,1].
For 6 > 0.875, using the assumption e-y/3 < 1, we have
s = [61/2 eV3(l 6)1/2]2 > [61/2 (1 6)1/2]2 = 1 2^6(1 6) =: (3.
Suppose we restrict the choice of 6 to 6 < |6. It then follows that
< (- (i+1) <) < (- (*+1) o s (* C1+f) 0
since s < 1. From (2.40), we conclude that
(s-(1+l)t)"s (i-(i + f)t)2.
Therefore,
(2.40)
H(b, 6) < U J (s (i + Â£) &) + 46(1 6) + (l f'] + > =: 77(6,6).
Simple calculus shows that
36 _ (3-/?) 36(1-6) 3
(3 + /?) (3 + /3)V P ~4
33
minimizes H and that b* > 0 if 6 > 0.875. After tedious but straightforward manipulation,
we have that H(b*, 6) attains its maximum at 5* 0.893 and that H(b*, 6*) < 0.988.
We are now ready to prove the central result of this section.
Theorem 2.11 (F-ellipticity of a(-, ))
Let a(-, ) and || ||~ be given as in (2.36) and (2.35). Suppose that 0 < a < 1, 0 < Â£ < ^=.
Then there exists a constant Ce > 0 such that, for all
a(v, v)
dv
Q-q^ + aPv + (I P)v
>Ce
+ M2)=Ce
(2.41)
where Ce = 0.012, which is independent of e and a.
Proof. We have
d(v, v) =
dv
Q + aPv + (I P)v
oz
+ a2||PS||2 + |](L-i>||2
+2 (Â£)+:2 (^,-p)")
+ a2\\Pvf + \\(I-P)v\\
+2a ( Qyz, v\ + 2(1 a) /q||, (I P)v
(2.42)
For the last term, we may write for any d G [0,1], by using (ii) of Lemma 2.2 and (ii) of
Lemma 2.9:
d (q^, (I P)v} + (1 -d) (Q^z, (I P)v}
= d(Qtd)'d (Qt *>) + (1 d) (Qt V ~ P*)
> -d (pQyz, Pv} + (1 d) ({I P)Q^z, (I P)v} .
Substituting this into (2.42) and bounding the fourth term in (2.42) by (ii) of Lemma 2.9,
34
we get
Z(v,v) > ||g|f|| +a*\\Pv\\* + \\(I-P)v\\*
-2(1 a)d\(PQÂ§,Pv) | 2(1 a)(l i) |((J P)Qf, (I P)*>|
which can be reduced by setting a = 0 to
a(v,v) > ||Qf||2 + ||(7_P)l;||2
-2d | (PQÂ§, Pv) | 2(1 d) |((7 P)QÂ§, (I P)v) |.
Defining
<:=!Â£Â¥, so that (1 ) = IIMC,
(2-43)
12
npqj?ii solhat fl T)_ii(f-p)osir
Htf llelfir '
and using Cauchy-Schwarz inequality, we conclude from (2.43) that
a(v, v) >
Q
dv
dz
+ (1 5)||t)||2 2dV6y/j
-2(i-d)V(T^)V(i^)
tdv
Tz
(2.44)
To maximize the lower bound in (2.44), we divide the region (5,7) E [0,1] x [0,1] into two
triangles and choose d as follows:
:-{i
for 5 + 7 < 1
for 5 + 7 > 1 '
Next, we consider these two cases separately.
Case 1: 5 + 7 < l,d = 1: For any 77 > 0 and any u,v and any norm || ||, the
arithmetic-geometric mean inequality
2IMHMI < 7lMI2 +;HMI2
holds. We thus have
a(v, v) > (1 777)
2
+
1
(2.45)
It remains to choose 77 such that the terms (1 777) and (1 5 ^) in (2.45) can be bounded
by a positive constant from below for all possible 7,5 with 7 + 5 < 1. For 5 < 0.5, we choose
77 so that
(1 77?) = (l 5 0 ,
35
which yields
77 =
6 + y/6*~+467
Applying Lemma 2.10, since 7 + 8 < 1 and therefore, 7 < 1 <5, then we have
71 = \{6 + V^+467} < \ + y/P + 45(1 5)} = H{ 0,8) < 0.988.
Thus, from (2.45), the V-ellipticity of a(-, ) directly follows with Ce > 0.012.
On the other hand, if 8 > 0.5, the second term in (2.45) can become negative. To
keep it positive, we then rewrite (2.45) for any 6 6 [0,1] as follows
a(v,v) > (1 b yr]) Q +b Q + (1 6 -
^dv 2 _dv
Q9! + 6 %
Since 8 > 0.5 and e < l/\/3 we have 0 < ^^=V6 ey/1 6^. We can now use the Poincare-
Friedrichs inequality (2.39) of Lemma 2.9 and inequality (i) of Lemma 2.9 to bound the
second term by
o _ 2
V6-ey/0^S) ||i;|p
dv 2 r
> imi2 >
which results in
a(v, H) > (1 b 777)
+
(-
- b 5
8+-s-----
3 77
where
Again, we choose 77 so that
s := (v^-ev/Ml-^))2-
(1 6 777) = (\ 8 + ^
(2.46)
(2.47)
which yields ___________________
- (6 + |s 8) + J(b+ fs 6)2 + 467
V = 27 '
Next, to attain a positive constant in the lower bound in (2.46), we need to show that for
all possible 8,7 with 8 + 7 < 1, 8 > 0.5, a positive b can be selected so that
1
Gi{b, 8,7) 6 + 777
6 + - 6^ + 467 + fb + 8 - } < 1.
Since Gi(6,8,7) < Gi(6,8,1 6) for 8 + 7 < 1, it then is sufficient to prove
V6e [0.5,1] 36 >0 : Gi(6,8,1 8) < C* < 1.
But this follows immediately from Lemma 2.10 with G* = 0.988 since Gi(6,8,16) = H(6,6).
The V-ellipticity of a(-, ) then follows directly from (2.46) with Ge > 1 0.988 = 0.012.
36
Case 2: 5 + 7 > 1: Setting d = 0 in (2.44) and preceding as in case 1 results in
2
d(v,v)> (1-(1-7)77)
dv
' dz
+1-5-
1-6
For 6 < 0.5, we choose
so that
S+^+4(l-5)(l-7)
V 2(1-7)
Using Lemma 2.10, since 6 + 7 > 1, and therefore, 5 > 1 7, we then have
1
(2.48)
(1 (1 7)77) = ^1 6 -
- 7 > 1, and therefore, <5 > 1 7, w<
(1 -7)V = \ |\A2 + 4(1 5)(1 7) + 5} < | {y<52 + 4(l-5)5 + 5} = H(0,6) < 0.988.
The Vr-ellipticity of a(-, ) then follows with Ce = 1 0.988 = 0.012.
On the other hand, for 6 > 0.5, we introduce, as in case 1, a parameter b [0,1]
and use the Poincare-Friedrichs inequality (2.39) to conclude from (2.48) that
a{v,v) > (1 b (1 7)77)
S+-s-
1-6
f]
(2.49)
with s as defined in (2.47). Again, we first choose 77, so that
(1 6 (1 7)77) = (l 6 + Is -
which yields
O + Is s) + \/(&+!s-<5) +4(1 <5)(1 7)
V~ 2(1-7)
In order to attain a positive constant in the lower bound (2.49), we need to show that, for
all possible 6,7 with 6 + 7 > 1, 6 > 0.5, a positive 6 can be selected so that,
G2(b,6,j) := 6+ (1-7)77
1
2
b+-s-5
2
+ 4(l-)(l-7) +
i + -3*
< C* < 1.
Since G^b, 6,7) < G2(b, 6,15) = H(b, 6) for 5+7 > 1, this follows directly from Lemma 2.10
with C* 0.988. Finally, from (2.49) with Ce := 1 C* = 0.012 the U-ellipticity of a(-, )
follows, which proves the theorem.
From the U-ellipticity of the bilinear form a(-, ) the U-ellipticity of the bilinear
form a(-, ) can be proved as follows.
Corollary 2.12 ( U-ellipticity of a(-, ) )
Let a(-, ) and H-Hy be given as in (2.20) and (2.37). Assume that 0 < a < 1, 0<Â£<
Then there exists a constant Ce > 0 such that, for all dGU,
a(v,v) > Ce \\v\\2v, (2.50)
37
where C& = 0.012, which is independent of a and e.
Proof. By the definition of the norm ||-||K in (2.37) and the relation (2.17) we have =
||i)||p. Therefore, using (2.41) in Theorem 2.11 we obtain for any v 6 V
a(v, v)
Â£v Â£v dfidz
Â£Sv Â£Sv dfidz
= a(v,v) > Ce
KM
2
V >
which proves the corollary
2.5 Error Bounds for Diffusive Regimes
Using continuity and V-ellipticity of the bilinear from a(-, ) in the norm ||-||v with
constants independent of a and e, in the following section we establish discretization error
bounds that do not blow up in the limit e 0. We use the same discrete spaces introduced
in Section 2.3.2.
We first consider discrete spaces with functions that can be expanded into the first
N normalized Legendre polynomials with respect to the direction angle fi (Pn-i discretiza-
tion in angle) and are piecewise polynomials of degree < r in z on a partition 7), of the slab.
To combine Ceas Lemma (1.24) with the interpolation error bounds in (2.27) and (2.30) to
obtain an error bound for this class of discrete spaces, we need the following lemma.
Lemma 2.13 (Bound for commutator [IIjvÂ£ Til at])
Let Â£ be the transport operator as defined in (2.16), Â£s the Sturm-Liouville operator as
defined in (2.29) and II/y the projection operator as defined in (2.28). Suppose IV > 2, m > 0,
and v Â£ V fi Hm+3(D). Then, there exists a constant C > 0 independent of e and a such
that
[BjvT TII/v]
dmv
dzm
c dm+1v
- N2 Lsdzm+1
(2.51)
Proof. Recall that
and that has the moment expansion (see (2.13) in Section 2.2)
dmv
dzm
= Â£ Am\z)pi{p) with = \ ( 9
;=0 Ti
Note that
and
ftin /i , v v
n n p = <^nm p n jf
" dzm 0 dzm
N-1
= n* 4, = <Â£im) = Am)pi = p
1=0
38
Therefore,
[njv^-^Djv] = UN fi-^- fiUN
OZ Oz
Using the relation (see Chapter 4)
P Pi (p) = k-1 p(_ i (fi) + hi pl+1 (fi),
with
h -
l + 1
v/4(/+ 1)2-1
and p~i(p) = 0,
we see that
d dmv
UHIMTzd^ =IIiT
Pi+i
Lf=o
(=0
AT
JV-2
1=0
1=0
On the other hand,
d dmv
dz dzm
=/-x;V+v 1=0
= ^2 <{>\m+1)bi-iPi-i 1=0 N-1 + Â£#! 1=0
(m+l) ,
Thus, in combination with (2.52) we have
dmv
[n^-Â£ Â£Hjv] = (j)^+^ 6jv_i pn-i ~ ^n-i Pn
= bN-i (4>p?+1'*pn-i -
Now we notice that, for any integers k, l > 0,
/ (m+l)
(2.52)
(2.53)
ZT 1 Zr
J (<^;m+1)(z)) j p\{p) dpdz = 2 j (<^m+1)(z)) dz
Zi -1 z I
t)
(m+l)
<
W +1)]2
Cs
3m+1v
3zm+!
where the last inequality follows from (ii) of Lemma 2.6. Therefore, (2.53) can be bounded
39
as follows:
[IIjvÂ£ Â£II^r]
dmv
dz"
N
+
V4AT2 1 \N(N 1) N(N + 1)
Â£s
dm+1v
dzm+1
JV2- 1
Â£s
dm+1v
dzm+1
<
1 8 1
Â£s
dm+1ip
dz*+'
y/3 3 N2
since N^_1 < O x 2, < |, which is valid for N > 2.
We are now ready to prove the following error bound.
Theorem 2.14 (Finite-Element in space, spectral discretization in angle)
Suppose that Th is a partition of the slab [zj, zr\ with maximum mesh size h. Let V be given
as defined in (2.38). Define
Nl
V
ft ___________
Vh
e cm vh=Y, WnOO; (*) e iPr(Th) for / = o,..., n -1;
[=0
Vh(zi,fi) = 0 for // > 0 Vh(zr,fi) = 0 for /i < 0 > ,
where TPr{Th) denotes the space of piecewise polynomials of degree < r on the partition Th-
Suppose 0 < e < 0 < a.< 1 and 1 < m < r + 1. Let ip E VflHm+3(D) be the solution of
(2.20) with right-hand side qs E Hm+2(D). Let iph G Vh be the solution of (2.20) restricted
to Vh. Then
\\lp ihWv
2
+
l(I-P)(lP-lPh)
+\m-Tph)\\2
1/2
<
1
vc,;
Ciu r ,ca
C3hm
dmg,
dzm
hm
N2
Cs
dm+1ip
dzm+l
>
with Ci,C2, C3, C4 independent of a and e. In particular,
(2.54)
Proof.
d(ip iph)
dz
Using Ceas Lemma (1.24)
< eeh, \\(I P)(tp iph)\\ < eeh,
< eft, \\P(ip iph)\\ < eft.
and the U-ellipticity of the bilinear form a(-, ), from
40
Corollary 2.12 we conclude that
HV-V/illv < ~^rVaii> iph) < -7= min \/a{ip vh,ip vh)
VOe vOe vheVh
< ~^r^/a(^ ~ ^z^-Nip, ip TlzllNip) = -J=z \\Â£(ip n^njvVOII (2.55)
< -j= {\\a> crniPW + \\cnNij) rn.n^Vll},
since 11*11^ Â£ Vh.
In order to bound the first term in (2.55), we use (2.51) of Lemma 2.13 and (2.30)
of Lemma 2.6 to get
WW jcnNii>\\ < \\Â£ip n^'vil + \\nNCip Â£nNip\\
< ||g. -nw3s|| +
r dip
(2.56)
jCi r C2
n HÂ£s9sII 7V2
r dip
Cs dz
To bound the second term in (2.55), observe that PHzv = IIzPv for any v Â£ V,
since P operates only on fi. Therefore, denoting T -(I P) + aP, we have
\\Â£U.Nip Â£n.zUNip\\ <
C^SfiJlNip- ^ShUzILniP
+ \\TKNip THzUNip\\
d_
dz
-S(iuNip nz-SfiiLNip
S
+ \\TUniP n.THivV-ll.
Applying (2.27) to bound the last two terms leads to
\\Â£uNip jCnziiNip\\
Qm+l 2
dzm+1 e
S/jJIpjip
+ Chn
dm
dzm
TUniP
< Chm
Iljv
dmip
dzm
CHn
dmip
dzm
where we used the V-ellipticity of the bilinear form a(-, ) in the last step. We proceed
further by using (2.51) of Lemma 2.13 and the fact that Iljv is an orthogonal projection ((i)
of Lemma 2.6) to get
\\Â£uNip Â£uziiniP\\
njvÂ£
dmip
dzm
+
, dmip
- n^Â£]
dmqs
dzm
, C4 hm
+ VCIN*
dm+1ip
Ls dzm+'
(2.57)
41
Finally, substituting (2.56) and (2.57) into (2.55) results in (2.54).
Remark 2.15 (Interpretation of error bound)
In the following, we interpret the error bound (2.54) more closely. For diffusive regimes
(e < 1), the exact solution of the continuous problem has the diffusion expansion2 (see
(1.13) in Section 1.2.2)
=
while the the right-hand side in this case is assumed to have the form
2 = 2o (z) + efiqi(z) + 0(e2).
Therefore, it follows that JZs4> = 0(e) and Â£,sqs = 0(e2), since qs = Sq = Pq + e(I P)q.
Taking this into account, for the error e/, in (2.54) we get
eft
= yB:($0(Â£,)+Â£0W)+A(ai'
dmqs
dzn
+ Caj^O(s)
Thus, the error in the zeroth moment \\P(ip Vft)|| is bounded by 0(hm)+0(e) and the error
in the higher moments ||(I P)(ip ^>ft)|| is bounded by 0(ehm) + 0(e2). In particular, for
diffusive regimes, where e is very small, convergence of the discrete solution is also assured
by the above bound for small N, which is a reasonable choice in this case, since the exact
solution is nearly independent of fx.
Moreover, the bound in Theorem 2.14 directly gives the optimal order of conver-
gence for the spatial discretization without the use of Nitsches trick (Johnson [25, p. 97]).
For example, for piecewise linear elements, which means r = 1, we can choose m = 2 and
get an 0{h2) error bound in the L2-norm under the regularity assumptions on ip required in
the theorem.
On the other hand, if e is close to so that e > hm and s > 1/N, then an error
bound can be obtained more easily, since
Mlv <
i + 4
dv
dz
+
i+4
Ml
2
1,0
Therefore, for any v G VTl (Hm+1([zi, zr] x H2([ 1,1]), Ceals Lemma (1.24) and the bounds
in (2.27) and (2.30) for the interpolation error lead directly to
fC~ f 1\1/2
\H~M\v + IM-n*n^lli.o
i \ 1/2 / /-Y
1 + ?) hr||^||i,o + ^
m1
dmip
(2.58)
dzm
1,0/
However, we point out that this bound will blow up in the diffusion limit e 0 for fixed N
and h.
2The assumption that the exact solution has a diffusion expansion would have simplified the proof of
Theorem 2.14; however, we preferred to present here the more general version.
42
For the case e > ^=, which is not covered by Theorem 2.14, the error bounds in
Section 2.3.2 can be used instead.
The second main class of finite-dimensional spaces considered here are formed by
functions that are piecewise polynomials in space z as well as in angle fi. This choice
corresponds to a Finite-Element discretization in both space and angle. Because Section 2.3.2
contains error bounds for this class of discrete spaces that are valid in non-diffusive regimes,
we concentrate in the following only on error bounds for diffusive regimes. Therefore, we
assume in the following theorem, which combines Ceas Lemma and (2.33), that the exact
solution has a diffusion expansion, which simplifies the proof.
Theorem 2.16 (Finite-Elements in space and angle)
Suppose that 0 < a < 1 and 0 < e < . Let T/,bea partition of the computational domain
D = [zi, zr] x [1,1] into rectangles T= \zj,Zj+1] x [/i,/i+i] of maximum diameter h. To
be able to handle the boundary conditions properly, we assume in addition that (zj, 0) and
(zr, 0) are nodes of the triangulation. Let V be given as in (2.38) and define
W1 := | vh E C(D); vh\T = Â£ Cp,y z? f V T E Th
l 0^/2,7
Vh(zi,ft) = o for /I > 0,
Vh (zr /i) = 0 for \i< 0
Suppose ip E V C\Hr+1(D) is the solution of (2.20) and let iph G Vh be the solution of (2.20)
restricted to Vh. Further, suppose that the diffusion expansion ip(z, p) =
is valid. Then we have
\C<
\\ip- i>h\\v < \r?rChr (|Vlff'-+i(c) + Hr\h-+i(d)) ,
(2.59)
with C independent of e, a, and h.
Proof. From Ceas Lemma (1.24), we have
U-M\v <^feM-nhiP\\v
<
Ipfl^-UhTP)
+
{I P)ti{ip -uhip)
+
_(/ PM Jihip)
1/2
+ m n^voif
(2.60)
<4
-Pn^-{rp Hhip)
e az
+
(/ p)ii-^(ip n h*p)
+
~(i p)(v> n^)
+ ITO-n^)|| .
43
By (2.33) we now bound any of the above four terms separately and use the fact that lift is
as an interpolation operator linear, so that = Ilft0 + ellhR-
Since (z) in the diffusion expansion of if> is independent of angle //, we conclude
that Ilft^z) is also independent of fi. Therefore, Pfi= 0 = so that
ip/^-HftVO
1 d
(2.61)
< ||^ii Ilft^Rd^!^
where the last inequality follows from (2.33).
Using (2.33) for the second term in (2.60) results in
(/-P)/Z^(V<-IIftf)
5: IH HfcVllfri(D) < C hr IVI Hr+1(_D)-
(2.62)
Because and Ilft^ are independent of angle, we have (I P) = 0 = (I P)TLh-
Therefore, the third term in (2.60) can be bounded by
1
(I-P)W- IlftVO
-(i-p)e{R-nhR) <||^-nft^||
< Chr + 1\R\Hr+l(D) < Chr\^>R\Hr+l^D^
(2.63)
since h < (zr z{) = 1.
Similarly, a bound for the last term in (2.60) is given by
\m n^VOII < U n^ll < Chr+1 < Chr WW^y (2.64)
Inserting (2.61), (2.62), (2.63), and (2.64) into (2.60) results in (2.59), which proves the
theorem.
Remark 2.17 (V-ellipticity constant Ce)
Both error bounds (2.54) an
Ce- According to Theorem 2.11 and Corollary 2.12, Ce = 0.012, so 1/Ce = 83.3, which is
fairly large. However, we would like to point out that we simplified the proof of Theorem 2.11
by considering only the worst case a = 0. Without setting a = 0, (2.46) would change to
i(u,u) >
+ (1-6) ||u|| 2cf(l a)VV7
-2(1 d)(l a)y(T35)^(rZ7j
Q
dv
8z
HI.
HI
(2.65)
which clearly shows that the V'-ellipticity constant Ce increases with a. To judge the quanti-
tative behavior, we computed Ce for certain values of a using (2.65). The results are plotted
in Figure 2.1. Already for a = 0.3, 1/Ce drops down to 7.04.
44
Ce vs. alpha
Figure 2.1: V-ellipticity constant Ce and its reciprocal as a function of the absorption
parameter a.
45
CHAPTER 3
X-Y-Z GEOMETRY
In this chapter, we generalize the scaling transformation and the error bounds
for the Least-Squares Finite-Element discretization, from one-dimensional slab geometry to
three dimensions. Since the main focus of this chapter is on diffusive transport problems, in
the following we use the parameterized form (1.14) of the transport operator C. In addition,
we assume that the total cross section Ut is constant in space (
parameter e is constant on the computational domain D := 1Z x S1, where 1Z C IR3 is a
region with sufficiently smooth boundary, for example, of class C1,1 (Grisvard [22, p. 5]),
and S1 denotes the unit sphere. Further, we suppose throughout this chapter that a < 1.
Moreover, in the following we restrict our attention without loss of generality to problems
with vacuum boundary conditions (so g(r,Q) = 0 in (1.4)). Problems with inhomogeneous
boundary conditions can easily be transformed to problems with homogeneous boundary
conditions and different right-hand sides (Oden and Carey [45, p. 27]).
As in the one-dimensional case, a scaling transformation is applied to the transport
operator prior to the Least-Squares discretization to ensure accuracy of the discrete solution
in diffusive regions. In the three-dimensional case, the scaling transformation and its inverse
are given by
S := P + e(I P)] S'1 :=P + j(I-P). (3.1)
They have the same form as in the one-dimensional case with the only difference that r = e
and that the L2-orthogonal projection P onto the space of functions that are independent
of direction vector Q is. now defined by
Pil> := J ip(r,Q) dn. (3.2)
s1
After applying the scaling transformation S from the left to the transport operator C in
(1.14) and dividing by e, the transport equation becomes
W := -SXV = YV> + -(/ P)i> + aPif> = qs,
e e e
(3.3)
with qs := Sq.
Throughout this chapter, we denote the standard inner product and the associated
norm of L2(1Z x S1) by
u v* dQ, dr;
R s1
|]ti|| := y/(u, u) Vu,v G L2(TZ x S1),
where v* is the complex conjugate of v. Further, for u,v E C(D), we define the following
inner product
(u,v)v
II
TZ S1
-50- -Vu -50 Vv +
Â£ Â£
(I-P)u.(I-P)v
+ Pu Pv dQdr,
its associated norm
IMIv := y/(u>u)v =
1 2 1
-SQ Vu Â£ + -(I -P)u
+ IM2
1/2
(3.4)
(3.5)
and the space
V := {v E C(D); u(r, 0) = 0 for r E dlZ, and fi n{r_) < 0}, (3-6)
where the closure is with respect to the norm ||-||y, so that it is a Hilbert space.
As mentioned in the introduction, the Least-Squares variational formulation of
(3.3) is given by (1.19). Our first goal is to show that the bilinear form a(-,-), defined
in (1.19) is continuous and V-elliptic with constants independent of parameter e and a.
From these results will follow not only the well posedness of problem (1.19), as outlined in
the introduction, but also the accuracy of the Least-Squares Finite-Element discretization
applied to it in the diffusion limit.
3.1 Continuity and V-ellipticity
As in Chapter 2, we conclude from the Cauchy-Schwarz inequality that the bilinear
form a(-, ) is continuous, since
|a(u, u)| = |(Cu, Cv)\ < ||Â£u|| ||Â£v||.
We now use the discrete Holder inequality to get
||Â£u|| <
< Vs
SQ Vu
Â£
+
(I-P)u
+ HPul
1 2 l,
Vu Â£ + -Â£(I-P)u
+ Il-Pul
\ 1/2
I2 =vS!
IV I
since a < 1 by assumption. Thus, for any u, v E V,
K,tOI < Ce Hiy |M|V,
(3.7)
with Cc = 3.
For the more difficult part of the proof of the V-ellipticity of a(-, ), we proceed in
the same way as in the one-dimensional case. We first scale (3.3) in addition from the right
by S to get
Â£S$ = vV+ {I P)$ + ocPi> = qs,
47
(3.8)
where -0 := S 1'0. Define the new space V and associated norm by
where
V S-'V, |||$ :=||Q-Yt>f + ||v|
Q := -SQS = (1 e) (PQ + QP) + eQI.
Shortly, we will prove that the bilinear form
a(u, v) := J J CSu LSv dtldr V u, v Â£ V
n s1
(3.9)
is V-elliptic. The K-ellipticity of a(-, ) in (1.19) will then follow in the same way as in
Corollary 2.12.
Before we do this, we first establish the following lemmas, which are generalizations
of Lemma 2.2, Lemma 2.4, Lemma 2.9, and Lemma 2.10.
Lemma 3.1
For all u, v Â£ V, and e < 1, we have
(i) (Pu,v) = (u,Pv); and {(I P)u,v) = {u,(I P)v);
P2 = P; and (7 P)2 = (I P).
Thus P and (I P) are orthogonal projections.
(ii) (Pu, v) = (Pu, Pv); and ((7 P)u, v) = ((7 P)u, (7 P)v);
(iii) ||t)|| <
-5t
e
(iv) ||P||-e||(7-P)w|| < |H|;
(v) (fi Vv, v) > 0;
(vi) (Q Vv, u) > 0.
Proof. The proofs of (i), (ii), (iii), (iv) follow analogously to that in Lemma 2.2 and
Lemma 2.9. To prove (v), we apply the fundamental Greens formula (Ciarlet [16, p.34]) to
get
(Q Vu, v) :
n s1
//S2'
Vv v dQdr
We therefore have
J J v Q Vv dQdr + J J v2 Q n dQds.
n s1 on s1
(Q Vv, v) = i J J v2 Q-n dÂ£lds.
8KS1
48
Splitting the boundary dll x S1 into the parts T+ := {(r, Q) G dll x S1; 0. n(r) >0} and
r~ := {(r, fi) G dll x S1; Q-n(r) < 0}, the boundary integral becomes
//2
97IS1
H n dÂ£lds
vZQ. 21 dfids +
v2Q n dQds
?;2Q n dllds > 0,
since v(r,Q) = 0 for (r, Â£2) G F and v G V. Thus, altogether we obtain
(Â£2 Vu, u) > 0.
To prove (vi), we observe from (i) that S is self-adjoint with respect to ). Con-
sequently, we have
(Q Vv, v)
\ 1 1
SQS Vv,v) = (Q VSv, Sv) = (Â£2 Vv, v) > 0,
where the last inequality follows from (v).
Lemma 3.2 (Poincare-Friedrichs Inequality)
Suppose e < 1 and let 11 C M3 be a bounded domain. Then, for any v G V, we have
||v|| < diam(7Â£) [|Â£2-Vt7|| < diam(7Â£) ||Q-Vu||. (3.10)
Proof. For ri}rk G 11 let
[Li, Lk] := {n : r = 2Ii + (1 s)rk, for 0 < s < 1}
denote the line segment between and rk. Let arbitrary r Ell and Â£2 G S'1 be given. We
define
*1 := min{f G IR [r, r -(- fÂ£2] C 11}
^2 := max {2 G IR : [r, r + tÂ£2] C 11}
n := r + tiQ,] 7*2 := r + t2Q.
Then it is easy to see that Â£2 -2l(2Li) < 0. Taking into account the boundary conditions for
v G V, we therefore have that 9(7^, Â£2) = 0, hence,
r r
v(r, Â£2) = Jd/ds = j 0 Vv ds,
where ds denotes the arc-length differential along the line {r + tU,t G 1R}. Therefore, we
conclude from Holders inequality that
L L2 / U
Ke,2)I < / |Â£2 Vu| ds < J [Â£2 Vu| ds < diam^)1/2 I j |Â£2 Vu|2 ds
ni 211 Vi
\ 1/2
/
49
Applying Fubinis theorem, it follows that
2
J J kfcil)|2 dQdr
Li n s1
< diam(77)2 J J |fi-Vu|2 dtldr.
71 S1
R. S1
From the relation v = Sv and (iii) of Lemma 3.1, we thus have
-SQ.S Vv
e
= diam(77) ||Q Vv
988,
||u[| < diam(77) ||fi Vu|| = diam(77) [|fi VSu|| < diam(77)
which proves the lemma.
Lemma 3.3
Suppose 0 < Â£ < 1. Then, for any 6 Â£ [0,1], there exists 6 > 0 such that
H(b, S) := \ | ^ (5 (1 + |)6) 2 + 45(1 <5) + (l 0 6 + 5 J < 0.
where s := jtf1/2 e(l 6)1(,2J In particular, for 5 < 0.875, we can choose 6 = 0.
Proof. The only difference to the proof of Lemma 2.10 is that now s = j^1/2 e(l <5)1 ^2j
instead of s jV^2 eV3(l 6)1/,2J Therefore, when 6 > 0.875, we use the assump-
tion e < 1 to get s > 1 2yjd(l d) =: (3. Everything else is analogous to the proof of
Lemma 2.10.
We are now in a position to state the central result of this section.
Theorem 3.4 (17-ellipticity of (-, ))
Let a(-, ) and ||-||~ be given as in (3.9) and (3.8). Suppose that 0 < a < 1, 0 < Â£ < 1 and
that the diameter diam(77) of the domain1 77 is 1. Then there exists a constant Ce> 0 such
that, for all v Â£ V,
12
a(v,v)
Vv + aPv + (/ P)u
(3.11)
Ce (||Q Vti||2 + ||n||2) = Ce \\v\\2~,
where Ce = 0.012, which is independent of e and a.
Proof. In the proof of Theorem 2.11, we replace QÂ§% by Q Vv and for P and a(-, ) use
the definitions of this chapter. Then the proof of Theorem 3.4 follows exactly els the proof of
Theorem 2.11, except that the Poincare-Friedrichs inequality (3.10) of Lemma 3.2 and (iv)
of Lemma 3.1 are now used to get
||Q-Vu||2 > [V6 eVI <5]2 ||u||2
1This can be established by a simple transformation of the space coordinates r.
50
Therefore, s in (2.47) is replaced by s =
bound the functions Gi and G2.
[Vtf e-s/1 6
and Lemma 3.3 is applied to
From the V-ellipticity of the bilinear form a(-, ), the V-ellipticity of the bilinear
form a(-, ) follows immediately as in Chapter 2. We summarize this result in the following
corollary.
Corollary 3.5 ( V-ellipticity of a(-, ) )
Let a(-, ) and ||-||F be given as defined in (1.19) and (3.5). Suppose that 0
and that diam(P) = 1. Then there exists a constant Ge > 0 such that, for all v G V,
a(v,v) > Ce ||v||y , (3.12)
where Ce = 0.012, which is independent of a and e.
3.2 Spherical Harmonics
Since a truncated expansion into spherical harmonics (Pjv-approximation) is used
throughout this chapter for the the disretization in angle, we introduce here the spherical
harmonics and summarize important properties that are needed for the error bounds.
Recall that the associated Legendre polynomials are defined for / > 0 and
m = 0,..., l by (Margenau and Murphy [39, p. 106])
jm
(3.13)
where P;(/i) is the (unnormalized) Legendre polynomial of degree l. By the formula of
Rodrigues (Arfken [3, p. .554]) for the Legendre Polynomials given by
this definition becomes
1 j/+m
= (3.14)
Expression (3.14) can be used to extend the definition of P,m(fi) to negative integer values
of m. It follows that P,m(/r) and Pl~Tn(fi) are related by
pfmoo = (3-i5)
The associated Legendre Polynomials satisfy the following recurrence relations (Ar-
fken [3, p. 560]):
pr() = [(Z + mjP^/x) + (Z m + ljPfoOO] , (3.16)
y/T^fiPTO*) = 27^1 [^(^-^i-f1^)]. (3.17)
y/l-fPTijl) = ^-[(Z + mKZ + m-ljP^Oi)
(Z m + 1)(Z m + 2)PI^1(/i)] . (3.18)
51
These recurrence relations, although derived in (Arfken [3]) only for positive integers m,
remain valid for negative values of m. This can be easily checked by substituting in the
relation (3.15) into the left and right parts of the recurrence relations.
Further, the associated Legendre polynomials satisfy the orthogonality relation
-1
(3.19)
Based on the associated Legendre Polynomials the spherical harmonics are de-
fined by (Arfken [3, p. 571])
Y,m(6, p) := (-l)m Ci,m P,m(cos(fl)) eimv,
(3.20)
where
. (2l+l)(l-m)\
\l----n X,----
(/ + m)!
Here, 9 denotes the polar angle with respect to the z-axis, while p denotes the azimuthal
angle about the z-axis is.
The spherical harmonics form an orthonormal basis of L2(51): In particular,
2tr 7r
^: J J Y(9, p) YlYl'*(9, p) sin(0) d9 dp SU' 8mm>,
o o
which, by letting L! := (6, p) and dQ := sm(ff)dedv | can be written as
J-Ylm(a)Y,r'*(tt)dn = 6w6mm.,
(3.21)
where Y'*(Q) denotes the complex conjugate of Y (H) From (3.15) it follows directly
that
Y,-m(Q) = (-l)mY,m*({2). (3.22)
By the definitions of 6 and p, we have
= (f2x, Qy, fiz) = (cos(v?) sin(0), sin(^) sin(0), cos(0))
= (cos(^)a/1 p2, sin(y>)\/l -
(3.23)
with fi := cos(0).
The spherical harmonics satisfy the following recurrence relations.
Lemma 3.6 (Recurrence relations for spherical harmonics)
For all l > 0 and we have,
nxY,m = Pl^Y^t1 l,mC -Pl,-mYÂ£? +
Qy Y,m = i ((-l)A,m^r + + A.-mFTr1 -
QzY = J^mY! + TJ^mY] j,
(3.24)
52
where
&I,m
7l,m
(/ + m + 2)(1 + m + 1)
4(2Z + l)(2/ + 3) Pl,m
(l m)(l m 1)
4(2/ 1)(2Z + 1)
(/ m)(l + m)
(2/-l)(2/ + l)
Vl,m
(l + m + 1)(1 m + 1)
(2/ + 3)(2/+ 1)
Proof. To prove the first recurrence relation, we replace cos(^) by
ei
so that
Using for the first term recurrence relation (3.17) and for the second term relation (3.18),
after simple but tedious calculations we obtain the first recurrence relation. The second
recurrence relation follows in a way similar to the first. For the third relation, recurrence
relation (3.16) is used.
Since the spherical harmonics form an orthonormal basis of i2(S'1), every v Â£
Hr(TV) x i2(51) has an expansion of the form
oo l C
(r,Q)=S E hm(r)Yrm, with ^,m(r)= / w(r, ft)Y,m*01) dSl. (3.25)
/=0 m=l qi
For any v Â£ Hr(TZ) x iT2 (S'1), we define
jv-i ; r
^Nv(r,n):= 53 53 l>m(r)Yr(Q), with ^,m(r) = j v(r, Q)Y;m*(i2) dSl (3.26)
/=0 m= l
S1
as the truncated expansion of v into spherical harmonics. To bound the error of the truncated
expansion, in the following lemma we use the fact that the spherical harmonics are the
eigenfunctions of the Laplacian operator on the unit sphere, so
AnV,m(fi)
sin 9 39 (Sm6dd
__L_____
+ sin2 9 dip2
y,m(fl)
(3.27)
= -l(l + l)Y,m(
for / > 0 and m = / + 1,..., 0,..., l.
Lemma 3.7 (Truncated expansion into spherical harmonics)
Let /3 be any multi-index and recall that D@v := gxfi1dJyfild2?3 Suppose that v(r,Q) Â£
H^(H) X if2(51) and, for JV > 2, let IIjv be defined as in (3.26). Then
\\AnDev\\ for l > 0; -/ < m < l, (3.28)
||v njwll < ^||Anu||, (3.29)
53
with C independent of v and N.
Proof. By the definition of 4>i,m and (3.27), we have
2
1
[/(/ + i)F
(lDfv
1
[/(/+1)]2
2
where we applied integration by parts in the last step and took into account the fact that
the boundary terms vanish. Since ||T(m||L2(5i) = 1) then the Holder inequality implies that
||^m||<^||A.vD'4-
Further, since the spherical harmonics form an orthonormal basis, we have
CO l
CO l
A
||w-hhI2 = E X / \
]=Nm=-li l=Nm=-l
Using (3.28) now results in
lu i -
||v njvv|| <-||An|| X JX^ [/(/+!)] = X[/(/ + l)]2'
CO f
21+1
Since < js and < jf for N > 2, then we can bound the sum in the
following way:
V 2/ + 1 f. 2 f 2
(^-1):
<
N2
so that
||v njvv|| < ||Anv||
3.3 Error Bounds
In this section, we use continuity and V-ellipticity of the bilinear form a(-, ) and
Ceas Lemma to bound the error of a Least-Squares Finite-Element discretization applied to
problem (1.19). The discretizations considered here are based on finite-dimensional spaces
54
with functions that have a truncated expansion into spherical harmonics with respect to
direction angle Q and are piecewise polynomials of degree Iona triangulation of the region
% into tetrahedrons. This choice of finite-dimensional spaces corresponds to a discretization
by a spectral method in direction angle and a Finite-Element discretization in space.
In transport theory, the spectral discretization in angle using spherical harmonics as basis
functions is also called a Pjv-approximation.
Let Th be a triangulation of 1Z into tetrahedrons T of maximum diameter h. For any
u(r, Q) Â£ Hr+1(Tl) X L2{S1), let Il/jt; denote the interpolant of v by piecewise polynomials
of degree r on the triangulation Th- Then, similar to (2.27), it can be shown (Ciarlet [16])
that, for 0 < m < 1,
II nHL,0 < C hr+1~m Mr+1,0, (3.30)
where ||-||m 0 denotes the standard norm of Hm(TZ) x L2 (S1) and | |r+i,o denotes the standard
semi-norm of Hr+1(H) x L2(5'1) (see Section 1.4).
In order to combine Ceas Lemma with (2.29) and (3.30) to obtain a discretization
error bound, we need the following lemmas.
Lemma 3.8 (Bound for commutator [!!//Â£ Â£II/v])
Suppose N >2 and let the operator An be defined as in (3.27). Let v Â£ H1(It) x H2(S1).
Then there exists C > 0 independent of a and e such that
||[IIjvÂ£ Â£IIjv] v|| < |Anv|i,o. (3.31)
Proof. By expansion (3.25), it easily follows that HpfPv = PUpjv and IljvPil Vu =
Pfi VII/vu. Therefore,
[HjvÂ£ Â£11*] = IM- V-fi- VBjv.
(3.32)
Now using the recurrence relation (3.24) in Lemma 3.6, we get
OO /
dv
ox
E E ^(/v^-zv^-r1
1=0 ml
N 1
1=0 m=l
N-2 1
f=0 m=l
Similarly, we have
1=0 m=l
55
so that
m=N
N-l
Ed^N 1 ,m / \rm1 ym+l\
dx (.-N-l.-mljv OCN-l,mYN ) .
(3.33)
m=JV+1
Similarly we get
i(n',0'S_t* Â£}**) =
m= JV+1
3j/
and
N
dv d
nNnz -r qz ^-iijvv = Y
dN,m vm 9N-l,r,
m=N
dz
"Tat,
Vm ST' vm
(m7jy_i / , ^ ^AT-l,mJJV
m=-JV+l
Denoting by |77.| the L2-measure of 1Z it is clear that |]Y;m[| < \J\1Z\. Therefore, by (3.28) in
Lemma 3.7, we can bound (3.33) as follows:
dv 9 a dv
a ox < Asit- ox
2y/W
N(N + 1)
N
Y', @N,m +
m= N
2 VW\
(N-1)N
N-1
E
m=JV+1
The sums can be bounded by
f _ V l(N-m)(N-m-1) A (IV m) lV(2Ar + 1)
2-/ w-m 2^ \/ 4(2A^ 1)(2JV +1) ~ 9f'9Ar-1'\ 9f9/V-1'l -
m=-N
m=N
2(2AT 1) 2(27V 1)
IV
m=N
and
K yk1 l(N + m-l)(N + m) y^1 N + m (2JV 1)JV
4(2JV 1)(2JV H-1) ij+1 2(2* 1, 2(2JV-1)
so that
IljV fij; Clx "T
oa:
In the same way, it follows that
di) d
Djvfiy a--DyIl/yU
dy dy
<
<
*VW\
N
dv
An-r-
ox
N
An
dv
dy
56
and
dv d dv
Hiv Clz- uz Unv dz dz <
VW\
N
N(N +
TT X 7Nm +
' m=-N
We now continue by bounding the sums in the following way:
m*
N M 0 N-1 ____
^7Ar,m 2N 1 + 2N 1 ^ ^N2
-JV m=l
s2j^i+2?n(w-1)w=w
and
so that
V'1 - V'1 l(N + m)(N-m) ^ A .
X, w i,rn X Y (2iV + 1)(2JV 1) ^ TJV>
)VJ-1 m=-JV+l V ^ ^ A > m=N
m=N+1
----IlivU
0Z dz
< 3^
~ AT
An
<9u
3z
which proves the lemma.
Lemma 3.9
Let V and ||-||^ be given as defined in (3.8). Then for all v E V fl (H1( 1Z) x L2(51)):
1% < CIRIi.o.
with C independent of e and a.
Proof. By definition, it follows that
11% < ([(1 <0 {\m V|| + ||n -Pv.o11} + e\u- u||]2 + ||u||2)
1/2
Notice that
lia-vfi|| =
dv dv
dv
^ dv dv dv
< fbr-S- d x + + dz
< v^3 ||u||1)0
since |fij;| <1, |f2y| < 1, and |fi0| < 1. Similarly we have
||Â£2 PVu|| <
nxp
dv
dx
+
n P1
UyPdy
+
n*pl
and
< a/3 ||u||10 ,
||Pfi-Vu||<||Q-Vu||||1)0.
VN-l,m
(3.34)
57
From these bounds, (3.34) follows immediately with C = ^[3%/3]2 + lj = \/28.
Now we are in a position to establish the following error bound.
Theorem 3.10 (Finite Element in space, Pm in angle )
Let 7^ be a triangulation of 1Z into tetrahedrons of maximal diameter h. Suppose 0 < a < 1,
0 < e < ^=, and diam(77) < 1. Let V be given as defined in (3.6) and let Vh be defined by
Vh:=LhV: ft(r,Q) = f IPr(Th)\ ,
l 1=0 m=-l J
where fPr(7ft) denotes the space of piecewise polynomials of degree < r on the triangulation
7ft. Let il> eVn (Hr+1(Tl) x H2(S1)) be the solution of (1.19) with qs e L2{11) x H^S1)
and let tph 6 Vh be the solution of (1.19) restricted to Vh. Further assume that ip has the
diffusion expansion ip(r,Q) = Q). Then
U Mv (llA^s|| + |Antf|i,o)
fc,
+\ 7T C2 hr (|<Â£|r+i,o + |^ir|r+i,o),
V '-'e
with Ci and C2 independent of a and e.
Proof. By Ceas Lemma, we have
(3.35)
H-M\v
<\PÂ£ (\w n^v-llv + ||n^ nftnjvVll)
V Og
(3.36)
The first term can be bounded by using the V-ellipticity of a(-, ), (3.29) of Lemma 3.7 and
(3.31) of Lemma 3.8 in the following way:
W-nNnv < -4= ll^(^-nJvV,)ll
<
Vcl
1
-4= (\\C1> - + IIPjvjC Â£11*] tf[|) (3.37)
V^e
<
c
c,
We U l|An9i|1 + NMh
To bound the second term in (3.36), we use that, by definition, ||i||v = \\t>\\y and
that S~1EhUNip = TlhS-'ENip = Ilftl!*^ Therefore,
\\ENip -EhEN^\\v = E-M-ip- lift II* V
< C
n*v> nftiLvV
1,0
58
where the last inequality follows from (3.34) in Lemma 3.9. We now use (3.30) and the
diffusion expansion of ip to get
nNip nhn^
< chr n^ = or
V 1 Â£ n /
1,0
r+1,0/
since Plljv =: IIjvP.
Remark 3.11 (Nondiffusive regimes)
In Theorem 3.10, it is assumed that the analytical solution has a diffusion expansion in
order to get an error bound in (3.38) with a constant that is independent of parameter e.
For regimes where the diffusion expansion is not valid, j is of moderate size, so that there
is no need for an error bound that is independent of e. Therefore, in this case, (3.38) can
simply be bounded by
so that the overall bound becomes
lr+1,0 '
59
CHAPTER 4
MULTIGRID SOLVER AND NUMERICAL RESULTS
According to the theory derived in our earlier chapters, the Least-Squares approach
yields accurate discrete solutions, even for diffusive regimes. In this chapter we confirm this
very efficiently by a full multigrid solver.
The following tests are restricted to the one-dimensional transport problem (1.6).
For the discretization in angle, a Pjv-approximation is used, which is a spectral method
using the first N Legendre polynomials as basis functions. For the discretization in space,
slab (((j>f(z) G JPr(Th) with r = 1 or r = 2). Defining m := dim(IPr(7)l)) 1 and letting
{rjk(z), k = 0,1,..., m} be a Finite-Element basis of IPr(T/t) (in the case r = 1 for example,
r]k(z) are the usual hat functions), then the discrete space is given by
for all / = 0,..., N 1 and all k = 0,..., m.
For the computations in this chapter, we simplify the discrete problem (4.2) further
by using a Gauss quadrature formula to approximate all integrations over angle fi, resulting
in a Least-Squares discretization of either the 5V-flux equations or the moment equations.
These two semi-discrete forms of the transport problem are introduced in the next section.
result by numerical tests and demonstrate that the resulting discrete system can be solved
we employ a Finite-Element discretization with linear or quadratic basis functions. To be
more precise, we recall that the analytical solution has the moment expansion
For the discretization, we truncate the sum
N-1
h(z,fj.) = Â£>?(*) Mm)
(4.2)
4.1 Sjy-Flux and P/v-i-Moment Equations
4.1.1 S'jv-Flux Equations
Let
Recall that, for a function v E VN we can use a Gauss quadrature formula with N support
abscissas {fi\,..., and weights {wi,..., wjy} to write
(4.3)
since this quadrature with N support abscissas is exact for polynomials of degree < 2N 1
(Stoer and Bulirsch [49, p. 153]).
In the 5jv-discretization, the flux is only computed for the discrete set of angles
{fix,..., Hn}, so that the unknowns are given by the vector
fp :=
/ VfoMi) \
\ ^(z,hn)
By collocating at the Gauss points and approximating the operator P by the sum in (4.3),
the following SV-flux equations for ip can be derived from the transport equation (1.6):
Hip
M + CTt{lN ~ R) + VaR
where
q :=
( 1 ^
, I := , :=
\ q(z,m) / l 1 )
=q,
( u l
\ &N
(4.4)
M :=diag(/ii,...,//jv), i?:=WT.
Further, for v E VN we note that the scaling transformation S = P + t(I P) in this
notation becomes
Sv = Snv := [R + t(Ijv R)] u,
with In denoting the N x N identity matrix. Therefore, the scaled Sjv-flux equations are
given by
ILip := StfILip =
3
SnM + T(Tt(lN ~ R) + VaR
02
i> = q.>
(4.5)
with q := Sjv
In., 0
V>(z() =
9i{v-i) \
9i(v%) /
and
Q,Ii
ip{zr) -
/ griUK+x) \
V 9r{m) J
(4.6)
respectively, where In denotes the ^ X y identity matrix.
61
4.1.2 Moment Equations
In order to derive the moment equations from the flux equations, we note that, for
V> G VN,
N-l
tP=^2Mz)pi(v)- (4-7)
1=0
The moments are given by
1 } N
Mz) = o / $(z>P)Pi{p)dP = VHz>N)Pi(Pi) (4-8)
-'i
where in the last step we again used the Gauss quadrature formula. Defining
<Â£(z) := ((^0(z),..., ^jv_.!(2:))t and the matrices T and II as
PDij := Pi-i(Pj), == diag(wi,...,uN),
then the respective relationships (4.7) and (4.8) between the flux ij> and the moments can
be written in matrix vector notation as
l = TQ, and V>=TtÂ£. (4.9)
In the following lemma, we summarize simple properties of the matrices R and T that we
need for the derivation of the moment equations.
Lemma 4.1 (Properties of R and T)
We have:
(i) wTl 1, so R2 = R and R(In R) = 0;
(ii) Rt fl = QR;
(iii) TQTJ = IN and TJTi2 = IN]
(iv) Tfil = Tw = (1,0,..., 0)T and utTt = (1,0,..., 0);
(v) Leting :*-* 3 the n
r o bo 0
&0 0 h 0
TÂ£2MTt = 0 h 0 h 0 B
NxN
(vi)
r i 0 .. 0 1
0 0 .. . 0
TQRTt
_ 0 0 .. 0 NxN
62
Proof.
N . 1
(i) : wTl = E = 2 / 1 dp = 1. Therefore, R2 = luTluT = l(wTl)wT = luT = R.
i=i -i
(ii) : RtQ ~ wl_TÂ£l = u>uiT = filwT = QR.
N 1
(iii) : (TOTt) = Â£ Pi-i(Pk)ukPj-i{Pk) = \ f pi-i(p) pj-i(p) dp, = 6i:j.
,J *=i -l
Therefore, TQTt = I, so Tr nonsingular => 3C such that TtC = I
=> TQTJC = TQ => C = TO => TtTO = I.
N N 1
(iv) : (TO),. = J2 Pi-i(Vj)uj = E pi-i(Vj)ujPo(pj) = \ f pi-i(p)p0(p) dp = 6a.
j=i i=i -i
(v) : The unnormalized Legendre polynomials Pi(p) satisfy the recursion (Arfken [3, p.
540])
ppM = + ^^fl+iO*)- (4-10)
Since the normalized Legendre polynomials are given by pi(p) := y/21 + 1 Pi(p), from
(4.10) we have
PPi(p) = bi-ipi^i(p) + bipl+i(p,) (4.11)
with bi-i := \y-J2_1 Using (4.11), we then have
(MTT)i,j = ViPj-ifai) = h3-iPi-z(Pi) + bj-iPj(Pi)-
Therefore,
N .
[TQMTr]. = Y^Pi-i(Pk)uk (bj-2Pj-2(Pk) + bj-ipj(pk))
k=\
1 1
= bj-Jpi-i{p)pj-2{p) dp + bj-v^ Jpi-i(p)pj(p) dp
-1 -1
bj2^i j'1 T bj 15* j y -(_ 1.
(vi): mRTr = (TO 1) (wtTt) = Tu utTt = (1,0,..., 0)T(1,0,..., 0), where we used
(iv) in the last step.
Multiplying the unsealed flux operator IL by TT from the right and by TO from
63
the left and using Lemma 4.1 gives
TQILTt = TQMTt ^ +
( ' 1 0 . . 0 ' \ ' 1 0 . . 0 '
Bk+' In - 0 0 . . 0 + CTfl 0 0 . . 0
0 o , 0. ) _ 0 0 . 0.
(4.12)
~Btz +
=: JM,
where 1M is the unsealed moment operator. Therefore, multiplying (4.5) by TO from the
left and using (4.9) results in the following scaled moment equations:
TOZLV> = TnsNTTmmTr
= (TÂ£lRTr + tTOTt tTQRT7)
/ 1 T \
0 + tIn - 0
V 0 0 )
JMj>_
(4.13)
IM$h =: TMj>_ = TO<^
Again using relation (4.9), it follows from (4.6) that the boundary conditions for the moment
equations are given by
la, 0
2
TT{Zl) =
( 9l(j*l) \
\ gi(fJ-K) /
0 ,Ia\TTl(zr) =
/ gr(x+1)
V 9i(vn)
(4.14)
We conclude that the SV-flux equations and the P/v-moment equations are equiv-
alent semi-discrete forms of the transport equation. The difference between these two sets
of equations is that the non-derivative part in the flux equations is fully coupled, while the
derivative part is decoupled. For the moment equations, the reverse is true.
4.1.3 Least-Squares Discretization of the Flux and Moment Equa-
tions
After deriving the Â£V-flux and moment equations we return to the discrete P/v-
problem (4.2). Using a Gauss quadrature formula with weights {wi,...,wjv} and points
64
{fii,..., fix} to approximate the integration over angle n results in
Zv N
/ pi)] (fij) [Â£bki,(z, fi)] (Hj)dz
1 i=1
ZI N
I i=i
(4.15)
We note that on the left-hand side, the approximation of the integration over angle jx by the
Gauss quadrature formula is exact as long as / < N 1: since then Cbk,i is a polynomial
in fj, of degree l + 1, while Lijih is a polynomial in fi of degree N, then the product is
a polynomial in fi of degree < 2N 1, for which the Gauss quadrature formula with N
support abscissas is exact. Therefore, only in the equations for i, k = 0,..., m} must
we introduce an error on the left-hand side by approximating the integration over angle. On
the other hand, the same argument shows that the right-hand side is represented exactly
by the Gauss quadrature formula as long as qs(z,fi) has an expansion into the first N 2
Legendre polynomials.
With the notation introduced in Section 4.1.1 we have
( [Liph{z,ii)\ (m)
V [Abh(z,n)]{nN)
= and
[Cbkti(ztfi)](fxi) >
[Â£bki,(z,n)}(/iN)
JLr}k(z)tJ+1,
where tj+l denotes the (/+l)-nth column of the matrixTT defined in Section 4.1.2. Denoting
by (,')jftN the standard Euclidean inner product of 1RN, (4.15) then becomes
Z r
j (MLh,ILT)k(z)ti+1)mN dz
Z\
ZT
= J dz
Zl
(4.16)
for all k G {0,1,..., m} and / G {0,1,..., N 1}. Since the columns of TJ span 1RN, then
we can substitute {t^,.. by the canonical basis {e^ ..., ew} of IRN and we recognize
that (4.16) is a Least-Squares discretization of the SW-flux equations using the discrete space
' / Vh{z,fl i) > < N m
jl Hft = >v-h = J2Y,vik,ik(z)z-j
\ vh(z,fxN) y O II il II
(4.17)
This is the space of iV-vector functions whose components are piecewise linear (for r = 1)
or piecewise quadratic (for r = 2) polynomials on the partition Tk of the slab.
65
Using (4.9) and (iii) of Lemma 4.1, we can rewrite (4.16) as
J (pÂ£TTlh,TTTnÂ£TTTntf+1r)k(z)}iRN dz
*1
Zr
= J(^,TTmmTTmtJ+1r,k(z))MN dz
Zi
for all k G {0,..., m} and / G {0,..TV 1}, which by (4.13) and (iii) of Lemma 4.1 is
equivalent to
2r
j (Mlh,Mem{z))^ dz
(4.18)
Zr
= J (TQ.qs,Meir)k{z))^ dz.
Hi
This is a Least-Squares discretization of the moment equations using the discrete space
' f &%(z) > m
jl o 1* d) i-h . for l = 0,..., N 1
^ N-l(z) J k=0 J
All computations in the following sections are based either on discrete problem
(4.16) or (4.18).
4.2 Properties of the Least-Squares Discretization
In this section, we use the results of numerical experiments to observe properties
of the Least-Squares discretization.
The results plotted in Figure 4.1 and Figure 4.2 demonstrate the accuracy of the
Least-Squares discretization in combination with the scaling transformation for diffusive
regimes. The test problem we chose here is the same one used by Larsen et al. in [32].
The exact solution of the corresponding diffusion equation is (z) = 3/2z2 + 15z, which is
plotted in solid in Figure 4.1 and Figure 4.2. The scalar flux 0 '= Piph of the solution iph
of the Least-Squares discretization of the scaled transport equation using piecewise linear
elements in space is shown by the crosses. For the problem in Figure 4.1, where the absorption
cross section is zero, we used r = 1/of = e2 as the scaling parameter, which gives a higher
accuracy than the scaling with r = e. An explanation of this result is given in the analysis
presented in (Manteuffel and Ressel [38]). For the test problem in Figure 4.2, where cra ^ 0,
the scaling parameter was chosen to be r = \/(Ta/crt = e^/a. Taking into account the fact
that the mesh size of 1.25 is order of magnitudes larger than l/
both cases the results are very accurate.
Further, we mention that the Least-Squares discretization without the scaling trans-
formation results in the zero solution for both cases as indicated by the asterisks in Figure 4.1
and Figure 4.2. This outcome confirms the asymptotic analysis in Theorem 2.1, according to
which the scalar flux of the Least-Squares discretization of the unsealed transport equation
with piecewise linear basis elements in space is a straight line connecting the values at the
boundary in diffusive regimes.
Moreover, the asymptotic analysis in Theorem 2.1 asserts that the Least-Squares
discretization of the unsealed transport equation using piecewise polynomials of degree > 2 in
space has the correct diffusion limit. This too is supported by the observed maximum errors
for a Least-Squares discretization of the unsealed transport equation with piecewise quadratic
elements in space, which we list in Table 4.1. However, using the scaling transformation in
combination with the Least-Squares discretization with quadratic elements in space achieves
dramatically better accuracy in the discrete solution.
For piecewise linear elements in space, the error bound in Theorem 2.14 indicates
an 0(h2) behavior of the Least-Squares discretization error for a sufficient smooth solutions.
To analyze the order of the Least-Squares discretization numerically, we used a problem with
smooth exact solution sin(7rz). We then computed the discrete L2-error of the Least-Squares
discretization with linear elements in space for a sequence of grids that were created from the
coarsest grid by halving the mesh size from one to another grid. Table 4.2 depicts the ratio
of these errors for each two consecutive grids. The value of approximately 4 of this quotient
confirms numerically an 0(h2) behavior of the discretization error for linear elements.
The solution of the transport equation is physically a density distribution and
should therefore always be positive. The Least-Squares discretization has the drawback
that it does not in general guarantee a positive solution. This is shown by the example in
Figure 4.3, where the exact solution of the corresponding diffusion equation is again plotted as
a solid line and the discrete Least-Squares solution is depicted by the crosses. Of course, this
boundary layer can be resolved by refinement of the mesh, as shown in Figure 4.4. However,
in the region [2,10], the solution is nearly constant, so that a refinement makes sense only
in the region around the boundary layer. Therefore, the aim is to use adaptive refinement,
which can be combined very naturally with a full multigrid solver (McCormick [42]). One
easy criterion for determining the area of further refinement would then be to check where
the solution is negative. Of course, this has to be combined with more sophisticated criteria
that compare the solution of consecutive grids, for example.
Besides having the correct diffusion limit, a discretization for transport problems
must satisfy the extra condition to resolve, with a suitable fine spatial mesh, interior bound-
ary layers between media with different material cross sections. To test numerically if the
Least-Squares discretization meets these extra conditions, we used the test problem from
(Larsen and Morel [33]), which is given in Figure 4.5. The solid solution plotted in Fig-
ure 4.5 is computed by a Least-Squares discretization using 50 cells in both [0,1] and [1,11].
This solution approximates the exact solution plotted in (Larsen and Morel [33]) fairly well.
We see further that the boundary layer is not resolved fully when the mesh spacing for the
Least-Squares discretization is too coarse (crosses in Figure 4.5). In addition, the Least-
Squares solution itself indicates an error by becoming negative. Again adaptive refinement
would be an appropriate remedy.
67
Scalar Flux
dib
/i^- + 100(J-P)V> = 0.01
oz
-0(0, //) = 0 for fi > 0
-0(10, //) = 0 for fj, < 0
(e = 0.01, a = 0.0)
Solution of corresponding diffusion equation: cf>(x) = 3/2 x2 + 15s.
Figure 4.1: Scalar Flux of exact (solid) and Least-Squares solution with scaling transfor-
mation (crosses) and without (asteriks).
68
Scalar Flux
(e = 0.01, a = 1.0)
Solution of corresponding diffusion equation: = 3/2 x'2 + 15*.
Figure 4.2: Scalar Flux of exact (solid) and Least-Squares solution with scaling transfor-
mation (crosses) and without (asteriks).
69
Table 4.1: Comparison of maximum error for scaled and unsealed Least-Squares discretiza-
tion with piecewise quadratic elements in space.
OL 0.0 OL 1.0
l/h scaled unsealed scaled unsealed
4 1.8-10-4 5.2-10-2 1.2-10-4 3.2-102
8 1.2-10-5 1.2-10-2 7.7-10-6 7.7-10-3
16 ^3 1^ O 1 3 1 O i-H CO 4.8-10-7 1.9-103
32 4.8-10-8 1 o T( 00 3.0-10-8 4.6-10-4
64 3.0-10-9 1.8-10-4 1.8-109 1.1-10-4
128 1.8-1010 3.8-105 1.1-1010 2.2-10-5
256 1.3-10-11 6.3-10-6 7.9-10-12 3.8-106
Test problem:
[ftjl + vtil-P) + (z,n) =q for (z, y) 6 [0,1] x [1,1]
< -0(0, /z) = 0 for fi > 0 > ,
V^IjZ-O = 0 for fj, < 0
where = 1000.0, (ra = q := (in cos(irz) +
Exact Solution: = sin(Trz), Number of Moments: N = 2.
70
Table 4.2: Order of Least-Squares discretization for linear elements in space.
Â£ = 1.0
Â£ = 0.001
a = 0.0
a = 1.0
a = 0.0
||e2fc||3
ll^lh
o= 1.0
IlCah [|2
II gfc. II a
IMI2
2h
kid
limits
IMh
C2h
l5ll
IKIh
4
8
16
32
64
128
256
512
1024
2048
4096
8192
8.5-10-2
2.1 -102
5:4-10-3
1.3- 10-3
3.4- 10-4
8.5- 10-5
2.1-10-5
5.3- 10"6
1.3- 10-6
3.3- 10-7
8.2-10-8
2,0-10-
4.0
3.9
4.1
3.8
4.0
4.0
3.9
4.0
3.9
4.0
4.1
2.610
6.7- 10'
1.7- 10'
4.2-10'
1.0-10
2.6-10'
6.6-10
1.6-10
4.1-10'
1.0-10
2.5- 10
6.5- 10
3.9
3.9
4.0
4.2
3.8
3.9
4.1
3.9
4.1
4.0
3.8
1.5-10-2
3.8-10-3
9.7- 10-4
2.4- 10-4
6.1-10-5
1.5- 105
3.8- 10-6
9.5- 10-7
2.4- 10-7
6;o-io-8
1.4- 10-
3.5- 10-9
3.9
3.9
4.0
3.9
4.0
3.9
4.0
3.9
4.0
4.2
4.0
1.2-10-2
2.9- 10-3
7.5- 10-4
1.8-10-4
4.7- 10-5
1.2-10-5
2.9- 10-6
7.3-10-7
1.8- 10-7
4.6- 10-8
l.l-lO-8
2.8-10-9
4.1
3.9
4.1
3.8
3.9
4.1
3.9
4.0
3.9
4.2
3.9
Test problem:
[a*If +
V>(0,//)
q for (z,n) G [0,1] x [-1,1]
0 for fi > 0
0 for n < 0
where at = 1, cra = q := fiir cos(wz) +
Exact Solution: ij>(z,n) = sin(7r;z), Number of Moments: N = 4.
71
Scalar Flux
(e = 0.01,. a = 10.0)
Solution of corresponding diffusion equation:
= 1^e20v^0.+ (X + e2ojo_i) 6~VZ-
Mesh size: h ^, Moments: N = 4.
Figure 4.3: Example of violation of the positivity property by the Least-Squares discretiza-
tion
72
Scalar Flux
(e.= 0.01, a = 10.0)
Solution of corresponding diffusion equation:
Mesh size: h = , Moments: N = 4.
Figure 4.4: Refinement resolves boundary layer.
73
Scalar Flux
( dib
+ a-tip Pcrsip = 0.0
= 1 for fi > 0
. ip{ll,/j,) = 0 for // < 0 ,
2 0 < z < 1 , r\_/ 0 0 < z < 1
100 1 < z < 11 an
Number of Moments: N = 16. The solid-line solution is computed by the Least-Squares
discretization using 50 cells in both [0,1] and in [1,11].
Figure 4.5: Behavior of the Least-Squares discretization in the case of interior boundary
layers.
74
4.3 Multigrid Solver
In this section we describe the multigrid solvers, that were developed for solving
the problems resulting from a Least-Squares discretization of Sjy-flux (4.16) and moment
(4.18) equations with piecewise linear elements in space. We refer the reader who is not
familiar with multigrid methods to (Briggs [6]) for an introduction and to (Hackbusch [24])
and (McCormick [41]) for more advanced topics.
Essential for the efficiency of a multigrid solver is the proper choice of its compo-
nents, mainly the intergrid transfer operators, coarse grid problems, and relaxation schemes.
The choice of the first two components is naturally given by the Least-Squares variational
formulation: the sequence of discrete spaces V\ C V2 C C Vj = Vh determines the coarse
grid problems since they are just the restriction of the variational problem to these discrete
subspaces; the prolongation operator, which is a mapping from a coarse grid to the next
finer grid in the grid sequence, is formed directly by composing the isomorphisms between
the discrete spaces and their corresponding coordinate spaces with the injection mapping
between 14_i and 14 (Bramble [5]), (McCormick [43]); and the restriction operators, which
are mappings from a finer grid to the next coarser grid, are just the adjoints of the prolon-
gation operators. Therefore, the only multigrid components that need to be chosen here are
the sequence of discrete spaces and the relaxation.
No relaxation scheme is currently in use for transport problems that smooths the
error in angle and in space simultaneously. Thus, instaed of devising a multigrid scheme that
coarsens simultaneously in space and in angle, we consider first applying the multilevel-in-
angle technique of (Morel and Manteuffel [44]), which is based on a shifted source relaxation
scheme. After reducing the degrees of freedom in angle, a multigrid method in space is used
to solve the remaining discrete problem. Thus, here we consider only the development of a
multigrid solver in space.
For the discrete subspaces, we Use the Finite-Element spaces with linear basis ele-
ments on increasingly finer partitions (halving the cells) of the slab.
4.3.1 Sn Flux Equations
The stencil that results from a Least-Squares discretization of the Sn flux equations
(4.16) with these Finite-Element spaces is given in Appendix A and shows full coupling in
angle. This suggests the use of a line relaxation in angle, which updates all angles for a
given spatial point simultaneously. The matrix that must be inverted for each spatial point
for this scheme is of the form (see Appendix A)
:= (aifi + a2QM + az&M2) + (ci MwuTM + c2ljujJ) .
The first part is diagonal, and the second has the rank 2 factorization
(ci MujjjT M + c2uxjJ) =
Thus, A{ can be cheaply inverted by the Sherman-Morrison formula (Golub and Van Loan [20,
p. 51]). Our computational tests showed essentially no differences in the error reduction and
smoothing properties of this line relaxation scheme for various different orderings of the spa-
tial points. To save computational, we thus use this line relaxation scheme in a red-black
fashion, since then the residual after one relaxation sweep is zero at the black points and
75
not need not to be computed for the restriction to the next coarser grid. This scheme is also
more amenable to advanced computer aechitecture efficiency.
The convergence factors for this multigrid algorithm 1 listed in Table 4.3, are com-
puted in the following way. A problem with zero source term and and whose exact solution
equal is zero is used in combination with a randomly generated initial iterate. Then 30 multi-
grid cycles are performed and the convergence rate is computed from the geometric average
of the per-cycle reduction factors of the last 20 cycles. We thus reduce the influence of the
initial iterate on convergence and observe what tends to be the worst-case factors. Here we
study the (1, l)-V-cycle, which uses one relaxation before and one after coarse grid correc-
tion. Observed factors for
coefficient a. Factors for (2,l)-V-cycles are also included. Such factors are sufficient to get
a solution with an error on the order of the discretization error by one full multigrid cycle,
as demonstrated by the results in Table 4.4. The additional V cycle on the finest level
10, performed subsequent to the full multigrid cycle, is reducing the error only by a small
amount. Thus, we can conclude, that the error after the full multigrid cycle is completed is
already on the order of the discretization error.
4.3.2 Moment Equations
The Stencil for the Least-Squares discretization of the moment equation (4.18) is
given in Appendix B. In the interior of the computational domain, it is a 15-point stencil that
connects the neighboring spatial points and the two higher and two lower moments. At the
spatial boundary, however, the stencil couples all moments. Therefore, we use a line moment
relaxation, that updates all moments simultaneously for a given spatial point all moments
simultaneously. Since the efficiency of the smoothing again is observed to be independent
from the relaxation ordering, as in the Sn flux case, we use a red-black ordering of the lines.
The convergence factors for this multigrid algorithm, are listed in Table 4.5. For very large
values of at, this multigrid solver is more stable with regard to roundoff errors than the
multigrid solver for the Sn flux equations. Even for values of 106, we get (1, l)-V-cycle
convergence factors of order 0.1. Again, these convergence factors are sufficient to get a
solution with an error on the order of the discretization error by one single full multigrid
cycle, as demonstrated in Table 4.6.
1This algorithm was implemented in C++ and a special array class was designed for this purpose.
76
Table 4.3: Multigrid convergence factors for solving the flux equations.
(l,l)-V-cycle
O'* a = 1.0 a = 0.5 a = 0.25 a = 0.1 a 0.0
Id)0- 0.088 0.085 0.087 0.118 0.169
101 0.082 0.083 0.083 0.110 0.136
102 0.052 0.052 0.053 0.106 0.130
103 0.088 0.091 0.088 0.105 0.130
104 0.091 0.091 0.091 0.105 0.130
105 0.092 0.092 0.092 0.105 0.130
106 0.090 0.092 0.092 0.102 0.133
(2,l)-V-cycle
10 0.053 0.050 0.053 0.105 0.155
101 0.047 0.047 0.047 0.082 0.104
102 0.019 0.024 0.024 0.077 0.097
103 0.020 0.021 0.021 0.076 0.096
104 0.020 0.022 0.022 0.076 0.096
105 0.020 0.011 0.023 0.076 0.096
10s 0.019 0.023 0.018 0.077 0.099
Test problem:
\p-jh +
V>(0,/z)
= 0 for (z,fi) 6 [0,1] x [1,1]
= 0 for n > 0
= 0 for fi < 0
where aa =
Exact Solution: iJ>(z,(j,) = 0.
Initial Iterate: randomly generated grid function.
Mesh size: h =
Number of Moments: N = 8.
77
Table 4.4: Full Multigrid (l,l)-V-Cycle convergence factors for solving the SW-flux equa-
tions.
Test problem:
+
ip(0,n)
= q for Â£ [0,1] x [1,1]
= 0 for fj, > 0
= 0 for (i < 0
where cra = q := /iircos(Trz) + aa sin(7rz), Exact Solution: ip(z, /j.) = sin(7rz), Number of
Moments: N = 4.
78
Table 4.5: Multigrid convergence factors for solving the moment equations.
(1,1 -V-cycle
a 1.0 a = 0.5 a = 0.25 a = 0.1 o o II S
To0- 0.052 0.086 0.083 0.118 0.169
101 0.091 0.092 0.091 0.117 0.136
102 0.056 0.056 0.071 0.106 0.131
103 0.092 0.093 0.092 0.105 0.127
104 0.095 0.094 0.094 0.106 0.129
105 0.095 0.094 0.093 0.107 0.130
106 0.095 0.092 0.092 0.107 0.130
107 0.095 0.092 0.092 0.107 0.130
10 0.095 0.092 0.092 0.107 0.130
109 0.095 0.094 0.092 0.107 0.130
1010 0.095 0.094 0.092 0.106 0.130
(2,1' -V-cycle
0-( a= 1.0 a = 0.5 or = 0.25 a = 0.1 a 0.0
0.074 0.051 0.054 0.105 0.155
101 0.055 0.055 0.055 0.082 0.104
102 0.025 0.025 0.039 0.077 0.097
10 0.023 0.026 0.042 0.076 0.096
104 0.023 0.023 0.042 0.076 0.096
105 0.023 0.023 0.042 0.076 0.096
106 0.023 0.023 0.042 0.076 0.095
107 0.023 0.023 0.042 0.076 0.095
10 0.023 0.023 0.042 0.076 0.095
109 0.023 0.023 0.042 0.076 0.095
1010 0.023 0.023 0.042 0.076 0.095
Test problem:
[l*-jfc +
V>(0,/r)
= 0 for (z, fi) 6 [0,1] x [-1,1]
= 0 for fi > 0
= 0 for n < 0
where cra = ^.
Exact Solution: i>(z, fi) = 0.
Initial Iterate: randomly generated grid function.
Mesh size: h =
Number of Moments: N = 8.
79
Table 4.6: Full Multigrid (l,l)-V-Cycle convergence factors for solving the moment equa-
tions.
Test problem:
[Vlh + (z, fl)
V(O^)
0(1. A*)
= q for (z, fi) G [0,1] x [1,1]
= 0 for fi > 0
= 0 for n < 0
where cra = q := fnrcos(irz) + cra sin(7r.z), Exact Solution: ip(z, fi) = sin(7rz), Number of
Moments: N = 4.
80
CHAPTER 5
CONCLUSIONS
5.1 Summary of Results
In this thesis, we have studied a systematic Least-Squares approach to the neutron
transport equatio. The Least-Squares formulation converts the first-order transport problem
into a self-adjoint variational form, which makes it accessible to the standard Finite-Element
theory. Essential for this theory is the V-ellipticity and the continuity of the variational form,
which leads directly to the existence and uniqueness of the analytic and discrete solutions
and to bounds for the discretization error for a variety of different discrete spaces. Moreover,
the variational formulation guides in a natural way the development of a multigrid solver for
the resulting discrete problem. However, due to special properties of the transport equation,
the Least-Squares approach is less straightforward than it first appears.
In this thesis, we focused on neutron transport problems in diffusive regimes. In
these regimes, the transport equation is singularly perturbed and its solution tends to a
solution of a diffusion equation. Therefore, to guarantee an accurate discrete solution, a
discretization for the transport operator is needed, that becomes a good approximation of
the diffusion operator in diffusive regimes. Only a few conventional discretization schemes
are known to have this property.
By an asymptotic expansion, we show in Theorem 2.1 for slab geometry that a
Least-Squares discretization with piecewise linear elements in space fails to be accurate in
diffusive regimes. The choice of linear elements in space will for any right-hand side always
result in a straight line connecting the prescribed values at the boundary for the principal
part of the solution, which is independent of direction angle fi. Numerical tests confirm this
behavior. On the other hand, we prove in Theorem 2.1 that, if piecewise polynomials of
degree > 2 are used, then the principal part of the discrete Least-Squares solution becomes
a Galerkin approximation to the correct diffusion equation in diffusive regimes. This means
that the Least-Squares discretization will be accurate in this case. Numerical tests with
piecewise quadratic elements again confirm this result.
Because of Ceas Lemma, the Least-Squares discretization can be viewed as the best
approximation to the exact solution in the discrete space with respect to the Least-Squares
norm ||Â£-||, where C denotes the transport operator. In diffusive regimes, the different terms
in the transport operator become totally unbalanced, which means that different parts of
the solution are weighted much differently by the Least-Squares norm. With P denoting the
L2-orthogonal projection onto the space of functions that are independent of direction angle,
it is clear that the Least-Squares norm in diffusive regimes hardly measures the components
Pip of the solution ip, although this is the main component for these regimes. The idea
is therefore to scale the transport operator prior to the Least-Squares discretization, with
the effect of changing the weighting in the Least-Squares norm. Clearly, the scaling from
the left by S = P + r(I P) with r = 0(1/
solution component Pip. Numerical tests show that a Least-Squares discretization of the
scaled transport equation, even for piecewise linear elements in space, yields an accurate
solution in diffusive regimes. Moreover, they show for piecewise quadratic elements that the
scaling transformation dramatically increases accuracy.
The major part of this thesis is devoted to proving that the Least-Squares dis-
cretization in combination with the scaling transformation S gives for a variety of simple
Finite-Element spaces always accurate discrete solutions, even in diffusive regimes. As men-
tioned above, essential for bounding the error is the V-ellipticity and the continuity of the
Least-Squares form with respect to some norm. It is easy to show that the scaled Least-
Squares form cannot be bounded from below by a standard Sobolev norm. Therefore, the
first obvious choice in the one-dimensional case is the norm (ll^fe-ll2 + II II2)1/2- With re-
spect to this norm, we prove V-ellipticity and continuity of the scaled Least-Squares bilinear
form and derive error bounds for discrete spaces that use piecewise polynomials in space
and piecewise polynomials or Legendre polynomials in angle as basis functions. However,
since the V-ellipticity and the continuity constants for this norm depend on at and aa, these
bounds blow up in diffusive regimes. To prove the V-ellipticity and continuity with constants
independent of cr* and aa, we use a scaled norm. Based on the V-ellipticity and continuity
with constants independent of at and aa, we obtain discretization error bounds for the same
discrete spaces mentioned above, with constants independent of at and aa. Thus, these
bounds stay valid also in diffusive regimes. This result is generalized to three-dimensional
x-y-z geometry for discrete space that use piecewise polynomials as basis functions in space
and spherical harmonics as basis functions in angle.
We conclude that the Least-Squares approach in combination with the scaling trans-
formation represents a general framework for finding discretizations for the transport equa-
tion that are accurate in diffusive regimes. Further, it naturally guides naturally the develop-
ment of an efficient multigrid solver for the resulting discrete system. This is demonstrated
in this thesis for slab geometry and piecewise linear elements. The developed multigrid solver
for this discrete problem has convergence factors on the order of 0.1, so that one full multigrid
cycle of this algorithm computes a solution with an error on the order of the discretization
error.
5.2 Recommendations for Future Work
Our numerical results show that, when simple discrete spaces in space are used,
refinement is needed in order to resolve boundary layers. Therefore, the aim for the future
would be to combine the full multigrid solver with adaptive refinement. On the other hand,
with the V-ellipticity and the continuity given, it seems fairly straightforward to establish
error bounds for more complicated discrete spaces that can better resolve boundary layers,
including those of exponential or hierarchical type.
Furthermore, generalization of the scaling technique to anisotropic transport prob-
lems suggests itself.
82
BIBLIOGRAPHY
[1] R.A. Adams, Sobolev Spaces, Academic Press, 1975.
[2] R.E. Alcouffe, E.W. Larsen, W.F. Miller and B.R. Wienke, Computational
Efficiency of Numerical Methods for the Multigroup, Discrete Ordinates Neutron Trans-
port Equations: The Slab Geometry Case, Nuclear Science and Engineering 71, pp.
111-127, 1979.
[3] G.B. Arfken, Mathematical Methods for Physicists, second edition, Academy Press,
New York, 1971.
[4] A. Barnett, J.E. Morel and D.R. Harris, A Multigrid Acceleration Method for
the One-Dimensional Sn Equations with Anisotropic Scattering, Nuclear Science and
Engineering 102, pp. 1-21, 1989.
[5] J.H. Bramble, Multigrid Methods, Pitman Research Notes in Mathematics Series 294,
Longman Scientific and Technical, Essex, 1993.
[6] W.L. Briggs, A Multigrid Tutorial, SIAM, Philadelphia, 1987.
[7] C. Borgers, E.W. Larsen and M. L. Adams, The Asymptotic Diffusion Limit of a
Linear Discontinuous Discretization of a Two-Dimensional Linear Transport Equation,
Journal of Computational Physics 98, pp. 285-300, 1992.
[8] S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element Methods,
Texts in applied mathematics, Springer Verlag Inc., New York, 1994.
[9] E. Broda, Ludwig Boltzmann. Mensch. Physiker. Philosoph., Franz Deuticke Verlags-
gesellschaft m.b.H., Wien, 1986.
[10] Z. Cai, R. Lazarov, T.A. Manteuffel and S.F. McCormick, First-Order System
Least Squares for Partial Differential Equations: Part I, SIAM J. Numer. Anal., Vol.
31, 1994.
[11] Z. Cai, T.A. Manteuffel and S.F. McCormick, First-Order System Least Squares
for Partial Differential Equations: Part II, submitted to SIAM J. Numer. Anal., March
1994.
[12] Z. Cai, T.A. Manteuffel and S.F. McCormick, First-Order System Least-Squares
for the Stokes Equation, submitted to SIAM J. Numer. Anal., June 1994.
[13] B.G. Carlson and K.D. Lathrop, Transport Theory The Method of Discrete Or-
' dinates, in Computing Methods in Reactor Physics, (H. Greenspan, C.N. Kelber, and
D. Okrent, eds.), Gordon and Breach, New York, p. 166, 1968.
[14] K.M. Case and P.F. Zweiffel, Linear Transport Theory, Addison-Wesley Publishing
Company, Reading, Massachusetts, 1967.
[15] C. CERCIGNANI, The Boltzmann Equation and Its Applications, Applied Mathematical
Sciences, Vol. 67, Springer-Verlag, New York, 1988.
[16] P.G. ClARLET AND J.L. Lions, Handbook of Numerical Analysis, v. II, Finite Element
Methods, Elsevier Science Publishers B. V. North-Holland, Amsterdam, 1991.
[17] J.J. Duderstadt and W.R Martin, Transport Theory, John Wiley & Sons, New
York, 1978.
[18] V. Faber AND T.A. Manteuffel, Neutron Transport from the Viewpoint of Linear
Algebra, Transport Theory, Invariant Imbedding and Integral Equations, (Nelson, Faber,
Manteuffel, Seth, and White, eds.), Lecture Notes in Pure and Applied Mathematics,
115, pp. 37-61, Marcel-Decker, April 1989.
[19] K.O. Friedrichs, Asymptotic Phenomena in Mathematical Physics, Bull. Am. Math.
Soc., 61, pp. 485-504, 1955.
[20] G.H. Golub and C.F. Van Loan, Matrix Computations, second edition, The Johns
Hopkins University Press, Baltimore, 1989.
[21] D. Gottlieb and S.A. Orszag, Numerical Analysis of Spectral Methods: Theory and
Applications, Regional Conference Series in Applied Mathematics, SIAM, Philadelphia,
1977.
[22] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman Advanced Publishing
Program, Boston, 1985.
[23] G.J. Habetler and B.J. Matkowsky, Uniform Asymptotic Expansion in Transport
Theory with Small Free Paths, and the Diffusion Approximation, Journal of Mathemat-
ical Physics 16, No. 4, pp. 846-854, April 1975.
[24] W. Hackbusch, Multi-Grid Methods and Applications, Springer, Berlin, 1985.
[25] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element
Method, Cambridge University Press, Cambridge, 1990.
[26] S. Kaplan and J.A. Davis, Canonical and Involutory Transformations of the Varia-
tional Problems of Transport Theory, Nucl. Sci. Eng., 28, pp. 166-176, 1967.
[27] J.R. Lamarsh, Introduction to Nuclear Reactor Theory, Addison-Wesley Publishing
Company, Inc., Reading, Massachusetts, 1965.
84
[28] E.W. Larsen, Diffusion Theory as an Asymptotic Limit of Transport Theory for Nearly
Critical Systems with Small Mean Free Path, Annals of Nuclear Energy, Vol. 7, pp. 249-
255.
[29] E.W. Larsen, Diffusion-Synthetic Acceleration Method for Discrete Ordinates Prob-
lems, Transport Theory and Statistical Physics, 13, pp. 107-126, 1984.
[30] E.W. Larsen, The Asymptotic Diffusion Limit of Discretized Transport Problems,
Nuclear Science and Engineering 112, pp. 336-346, 1992.
[31] E.W. Larsen and J.B. Keller, Asymptotic Solution of Neutron Transport Problems
for Small Mean Free Paths, J. Math. Phys., Vol. 15, No. 1, pp. 75-81, January 1974.
[32] E.W. Larsen, J.E. Morel, and W.F. Miller, Asymptotic Solutions of Numerical
Transport Problems in Optically Thick, Diffusive Regimes, J. Comp. Phys., 69, pp.
283-324, 1987.
[33] E.W. Larsen and J.E. Morel, Asymptotic Solutions of Numerical Transport Prob-
lems in Optically Thick Diffusive Regimes II, J. Comp. Phys. 83, (1989), p. 212.
[34] E.E. Lewis and W.F. Miller, Computational Methods of Neutron Transport, John
Wiley & Sons, New York, 1984.
[35] T.A. Manteuffel, unpublished personal notes on even-parity.
[36] T.A. Manteuffel, S.F. McCormick, J.E. Morel, S. Oliveira and G. Yang,
A Fast Multigrid Solver for Isotropic Transport Problems, submitted to SIAM J. Sci.
Comp., to appear.
[37] T.A Manteuffel, S.F. McCormick, J.E. Morel, S. Oliveira and G. Yang, A
parallel Version of a Multigrid Algorithm for Isotropic Transport Equations, submitted
to SIAM J. Sci. and Stat. Comp. 15, No 2, pp. 474-493, March 1994.
[38] T.A. Manteuffel and K.J. Ressel, Multilevel Methods for Transport Equations
in Diffusive Regimes, Proceedings of the Copper Mountain Conference on Multigrid
Methods, April 5-9, 1993.
[39] H. Margenau AND G.M. Murphy, The Mathematics of Physics and Chemistry, sec-
ond edition, D. Van Nostrand Company, Inc., Princeton, 1968.
[40] W.R. Martin, The Application of the Finite Element Method to the Neutron Transport
Equation, Ph.D. Thesis, Nuclear Engineering Department, The University of Michigan,
Ann Arbor, Michigan, 1976.
[41] S.F. McCormick, Multigrid Methods, Frontiers in Applied Mathematics 3, SIAM,
Philadelphia, 1987.
85
[42] S.F. McCormick, Multilevel Adaptive Methods for Partial Differential Equations,
Frontiers in Applied Mathematics, SIAM, Philadelphia, 1989.
[43] S.F. McCormick, Multilevel Projection Methods for Partial Differential Equations,
SIAM, Philadelphia, 1992.
[44] J.E. Morel and T.A. Manteuffel, An Angular Multigrid Acceleration Technique
for the Sn Equations with Highly Forward-Peaked Scattering, Nuclear Science and En-
gineering, 107, pp. 330-342, 1991.
[45] J.T. Oden and G.F. Carey, Finite Elements, Mathematical Aspects, Volume IV,
Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1983.
[46] R.K. Osborn, S. Yip, The Foundations of Neutron Transport Theory, Gordon and
Breach, Science Publishers, Inc., New York, 1966.
[47] A.I. Pehlivanov, G.F. Carey and R.D. Lazarov, Least-Squares Mixed Finite
Elements for Second-order Elliptic Problems, SIAM J. Numer. Anal., Vol. 31, No. 5, pp
1368-1377, October 1994.
[48] G.C. Pomraning, Diffusive Limits for Linear Transport Equations, Nuclear Science
and Engineering 112, pp. 239-255, 1992.
[49] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, second edition, Texts
in applied mathematics, Springer Verlag, New York, 1993.
[50] S. Ukai, Solution of Multi-dimensional Neutron Transport Equations by Finite Element
Methods, Journal of Nuclear Science and Technology, 9(6), pp. 366-373, 1972.
[51] G.M.Wing, An Introduction to Transport Theory, John Wiley and Sons, Inc., New
York, 1962.
86
APPENDIX A
FLUX STENCIL
In this section, we derive the stencil for the Least-Squares discretization of the Sn
flux equations (4.16) for piecewise linear elements. We assume that the slab is partitioned
into zi zq < zi << zm = zr and denote by h* := z*, Zk-i the cell width of cell k.
We are then looking for a discrete solution in the form
m N m
h(*)=EE^ = Y.Mz)>
fc=0 j = l
Jfe=0
where := ..., ^fc,jv)T and
with
Vizi+1) 0 elsewhere
Zh-z hk ~j
2r,*jW II l
= (0, ,o)T.
Plugging (A.l) into (4.16) gives then
(A.1)
Â£
fc-1 Zk-1 k-1Zk-l
for all
(i,j) e {(i, i) : 1 < Z < rn, 1 < i < ivj U j(0, j) : 1 < j < j} U |(m, j) : y < j < ivj .
Since j has support only in cell i and cell i + 1, then (A.2) is equivalent to
zi *i+1
J (aiL>1Lnr,ij)]Rrd* + J {QIL>IL!h,i+ij)Kd*
(*i)
(*2)
Z% zi+l
f (qq ,ILt] \ dz + [ (dq ,TLr\1 \ dz
J \ ** -r,hJ/lRN J \ is A*+WjRW
(A.3)
Zi-1
(*3)
(*4)
In the following, we consider the terms (*1) to (*4) separately.
To (*1):
Applying the substitution z = z z,-_i, we have
Zi
j (^,%r..)^dz =
Zi-1
dz
where h hi and
i :=
-r,J
Z
:= h**
with ipl := V'fc_1 and i/,r := }P_k Then it follows that
11% j = \ T ~ T)] Afe.j + ^ [r<7t(J lwT) + cralwT] ze;-
1L = ^ [lwT + r(J lwT)] M ^ ^ [r7t(J lwT) + o-alwT] - z) + rz
88
so that
J(niLt,ILrLri.)]RNdz
o
h [lwT + t(I lwT)] M (ifrr - , i [luj + t(I lwT)] Me^J
+y ^ jr [r
+Y [i^T + r(J ^-T)] M (^r ^j) X ['r<7<(/ i^T) + ^lw7] e_j ^
+y |j<7t(J lwT) + <7alwT] ipj, ^ [rat(I lwT) +
+y [Tr, ^ [r
= l ( lM^T + ^{1 -hLX)]M (r - , e;)^
+\ (Q KMIwt + rVtM(J lwT)] (3^ + V-J ^ ^
(fi [^aIwT + rVt(I lwT)] M (jfr. 3^) ,Â£,
+/, ^SJ KliST + rV?(7 laT)] (=t + ==) .S,)^
Consider all possible j and recall that h = h{ and that T are the corresponding values
in cell i which we denote by cra i, crt i, respectively r,-, and that 3^ = ^Pi_l,'Pr = 3^. We then
89
get the following contribution from (*1):
1
\(t? 1) Mum7M t?Q.M2]
fli
[W 0-o.i) (MuuJ wwtM)]
+ J - Ti*i,d ^T] j 4-1
+ [(i T?) MuujJM + rfQ,M2]
[(Ti(TtiQM]
+T t7"2*?,.'0 + (ah Tiali) f i-
To (*2);
Applying the substitution z = z Zi, we have
2i+l
-2*
fc+l
J (QIL&>1L?hj)KNdz>
where now (/i := hj+i) and
v>
h z
hi := *
Therefore,
ILrij . = [lwT + r(7 lwT)] Mej + ^ [r<7t(7 lu/r) + cralwT] (h z) e_j
K, = ^ [lwT + r(7 lwT)] M ij^j + ^ [r
90
so that
J(tlJL^lL^^dz
o
h [lwT + t(I lwT)] M ^ ^ ^ [lwT + t(I lwT)] Me^
+y [r
+y [JmT + T'U lwT)] M (r V>() i [r
+y -i^-T) + o'alw1"] t_v ^ [r<7t(J- lwT) + o-alwT]
+y [r
= x (Q [m^T + t2m(/-^ M (t- -&) .s#)^
[
+ ^ (fi KiwT + T2r - ,5#')^
+4 (si KlMT + *VU luT)] (t + &) ,
Consider all possible j and recall that h = h8-+1 and that (cra,at,T) are the corresponding
values in cell i + 1, denoted by (ca,i+i,
the following contribution from (*2):
91
(tt- K1 T?+i) MuultM + r?+1QM2]
+ ^ [('r/+l<7*,i+l O'a.i+i) (MywT + uwJM) 27f+1o-M+1ftM]
+ ^Yi W+lO'M+l^ + K,i+1 T?+1Or2<+1) WWT] | Â£.
+ (T- [ft+i -1)M<^TM 7i2+1fiM2]
[(r2+1trtii+i o-a.i+i) (Mwut wwtM)]
+^r1 N+i^+i0 + Ki+1 ^+1^+1) ^T] } Â£+!
To (*3): We have
where
^(z) = SjvÂ£(z) = IwTÂ£(z) + r(7 lwT)g(z),
/ \
(*) := :
\ q(z,Vn) )
Therefore,
Zi
/ (nt.-Oa.u)*, *
Zt-l
*i
= j (filwTÂ£, [hjT + r(I lwT)] Me^)^ dz
Zi-1
+ j (filwTg, [Tat(I lwT) +
Zi- 1
+ J (nT(I-lwT)q,[lur +T(I-lwT)]MejT]rti)JRNdz
Zi
+ J 1wt)Â£, [r
92