Graziani F. (editor) Computational Methods in Transport

Подождите немного. Документ загружается.

366 D.E. Shumaker and C.S. Woodward

0 0.005 0.01 0.015

Temperature (keV)

time (µs)

−3

−2

−1

Fuel Density (g/cm

)

Fig. 6. Time dependence of radiation and material temperatures and fuel density

at x =0,y =0

however, limit the exchange of energy between radiation and material. The

combination of a large diﬀusion loss rate and a poor coupling of radiation and

material energy results in a lower material temperature. Due to the strong

dependence of the reaction rate on temperature, T

, the lower temperature

outer region can have a signiﬁcantly lower reaction rate. This lower rate

results in incomplete consumption of fuel in this region.

Figure 7 give contours of material temperature, radiation temperature

and fuel density at t =0.0085 µs. At this time, the fuel has been signiﬁcantly

depleted in the central region. The material temperature is high around the

inner edge of the remaining fuel and the central material temperature is held

down by the transfer to the cooler radiation. The radiation is heated by the

cylindrical shell-like region of high material temperature. This heating, along

with diﬀusion, keeps the radiation temperature in the central region ﬂat.

Lastly, the shell-like region of unconsumed fuel is still present at the end of

the simulation.

In order to study the convergence of this problem with respect to reducing

the tolerance, we made several simulations varying RTOL = ATOL from 10

−7

to 10

−10

for maximum orders of integration of 2 and 5. For this series of runs

the numerical statistical counters and parameters deﬁned by,

RTOL = relative tolerance,

MO = maximum order allowed,

NST = time steps,

NNI = nonlinear iterations,

NLI = linear iterations,

RT = run time in seconds,

Implicit Solution of Non-Equilibrium Radiation Diﬀusion 367

0.2 0.4 0.6 0.8

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x (cm)

y (cm)

Radiation Temperature (keV)

0.6 0.8 1 1.2 1.4 1.6 1.8

0.2 0.4 0.6 0.8

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x (cm)

y (cm)

Material Temperature (keV)

1 1.5 2 2.5 3

0.2 0.4 0.6 0.8

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x (cm)

y (cm)

Fuel Density(g/cm

)

0.005 0.01 0.015 0.02

Fig. 7. Radiation temperature, material temperature and fuel density at t =

0.0085 µs

are given in Table 2. Here we see a substantial decrease in the run times due

to using a higher maximum integration order. The time stepping algorithm is

able to use the higher integration order to meet accuracy requirements with

larger time steps resulting in a faster overall runtime for a given accuracy

with the higher maximum integration order. This advantage is seen in savings

for all integration statistics, including the cumulative numbers of linear and

nonlinear iterations required for solution. Lastly, we again see a signiﬁcant

increase in run time to meet accuracy requirements for RTOL values below

−9

Table 2. Statistics for 2D Fusion Fuel Problem

RTOL MO NST NNI NLI RT (sec)

−7

2 17,384 20,807 62,322 1,730

−7

5 14,510 17,479 52,369 1,449

−8

2 45,768 56,333 136,973 4,152

−8

5 35,686 45,797 116,451 3,347

−9

2 113,545 140,788 292,722 9,313

−9

5 81,535 110,791 224,178 7,414

−10

2 266,722 323,349 621,406 20,564

−10

5 172,137 242,944 433,689 15,500

As a measure of accuracy for these runs, we determined the maximum

relative error in material temperatures over the entire 2D grid. This quantity

is shown in Fig. 8. The relative error is computed with respect to the most

highly resolved run, RTOL = 10

−10

, MO = 5. This plot indicates that the

reduction in error is approximately linearly related to RTOL. Note that Fig. 8

also indicates that MO = 5 runs are more accurate than MO = 2 for the

same RTOL. The relative error plots are all similar in that the maximum is

368 D.E. Shumaker and C.S. Woodward

5 6 7 8 9 10

x 10

−3

−8

−7

−6

−5

−4

−3

−2

time (µs)

Max relative error

MO = 5

RTOL = 10

−7

RTOL = 10

−8

RTOL = 10

−9

MO = 2

RTOL = 10

−7

RTOL = 10

−8

RTOL = 10

−9

Fig. 8. Maximum relative error in material temperature vs time for maximum

order, MO, 2 and 5

during the sharp rise in material temperature. Relative error plots for the

other components of the system, radiation temperature and fuel density, are

similarandarethusnotshown.

5 Conclusions

We have presented fully implicit solutions of a highly nonlinear three equation

system including: (1) radiation energy evolution with diﬀusion and exchange

with material energy; (2) material energy evolution with reaction heating and

exchange with radiation energy; and (3) fuel density evolution with depletion

due to consumption. Our solution method makes use of high order in time

integration methods, inexact Newton-Krylov nonlinear solvers, and multigrid

preconditioning using a Schur complement strategy. We have demonstrated

accurate solution of the fusion heating and fuel evolution by comparison

with an analytic solution for a simpliﬁed problem. Using a 2D problem which

encompasses all the processes included in our model, we have demonstrated

convergence to a highly resolved result. Our test problems have shown that

using a higher order of time integration leads to a more accurate and eﬃcient

method than with lower orders.

Implicit Solution of Non-Equilibrium Radiation Diﬀusion 369

Acknowledgements

This work was performed under the auspices of the U. S. Dept. of Energy

by University of California, Lawrence Livermore National Laboratory under

contract W-7405-ENG-48.

References

1. A. J. Mewadows, Stellar Evolution, Pergamon Press, 1978, p. 97.

2. J. J. Duderstadt, G. A. Moses, Inertial Conﬁnement Fusion, John Wiley and

Sons, 1982, p. 131.

3. D. A. Knoll, W. J. Rider, G. L. Olson, Nonlinear convergence, accuracy, and

time step control in nonequilibrium radiation diﬀusion, J. Quant. Spec. and

Rad. Trans. 70 (1) (2001) 25–36.

4. V. A. Mousseau, D. A. Knoll, W. J. Rider, Physics-based preconditioning and

the Newton–Krylov method for non-equilibrium radiation diﬀusion, J. of Com-

put. Phys. 160 (2000) 743–765.

5. D. J. Mavriplis, Multigrid approaches to non-linear diﬀusion problems on un-

structured meshes, Num. Lin. Alg. with App. 8 (8) (2001) 499–512.

6. L. Stals, Comparison of non-linear solvers for the solution of radiation transport

equations, Elec. Trans. Num. Anal. 15 (2003) 78–93.

7. P. N. Brown, C. S. Woodward, Preconditioning strategies for fully implicit

radiation diﬀusion with material-energy transfer, SIAM J. Sci. Comput. 23 (2)

(2001) 499–516.

8. P. N. Brown, D. E. Shumaker, C. S. Woodward, Fully implicit solution of large-

scale non-equilibrium radiation diﬀusion with high order time integration, J.

Comp. Phys. 204 (2) (2005) 760–783.

9. D. L. Ropp, J. N. Shadid, C. C. Ober, Studies of the accuracy of time integration

methods for reaction-diﬀusion equations, J. Comp. Phys. 194 (2004) 544–574.

10. C. C. Ober, J. N. Shadid, Studies on the accuracy of time integration methods

for the radiation-diﬀusion equations, J. Comp. Phys. 195 (2004) 743–772.

11. G. C. Pomraning, The Equations of Radiation Hydrodynamics, Pergamon, New

York, 1973.

12. R. L. Bowers, J. R. Wilson, Numerical Modeling in Applied Physics and As-

trophysics, Jones and Bartlett, Boston, 1991.

13. E. M. Corey, D. A. Young, A new prototype equation of state data library, Tech.

Rep. UCRL-JC-127698, Lawrence Livermore National Laboratory, Livermore,

CA, submitted to American Physical Society Meeting (1997).

14. J. Wesson, Tokamaks, Clarendon Press, Oxford, 2004, p. 13.

15. D. J. Rose, J. Melville Clark, Plasmas and Controlled Fusion, M.I.T. Press,

1961, p. 322.

16. T. J. Dolan, Fusion Research Vol. 1, Principles, Pergamon Press, 1980, p. 29.

17. G. D. Byrne, A. C. Hindmarsh, PVODE, an ODE solver for parallel computers,

Int. J. High Perf. Comput. Appl. 13 (1999) 354–365.

18. A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E.

Shumaker, C. S. Woodward, SUNDIALS: Suite of nonlinear and diﬀeren-

tial/algebraic equation solvers, ACM Trans. Math. Softw. 31 (3) (2005) 363–

396.

370 D.E. Shumaker and C.S. Woodward

19. P. N. Brown, G. D. Byrne, A. C. Hindmarsh, VODE: A Variable-Coeﬃcient

ODE Solver, SIAM J. Sci. Stat. Comput. 10 (5) (1989) 1038–1051.

20. K. R. Jackson, R. Sacks-Davis, An alternative implementation of variable step-

size multistep formulas for stiﬀ ODEs, ACM Trans. Math. Software 6 (1980)

295–318.

21. Y. Saad, M. H. Schultz, GMRES: A generalized minimal residual algorithm for

solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (3) (1986)

856–869.

22. S. F. Ashby, R. D. Falgout, A parallel multigrid preconditioned conjugate gra-

dient algorithm for groundwater ﬂow simulations, Nuclear Science and Engi-

neering 124(1) (1996) 145–159.

23. S. Schaﬀer, A semi-coarsening multigrid method for elliptic partial diﬀerential

equations with highly discontinuous and anisotropic coeﬃcients, SIAM J. Sci.

Comp. 20 (1) (1998) 228–242.

24. P. N. Brown, R. D. Falgout, J. E. Jones, Semicoarsening multigrid on distrib-

uted memory machines, SIAM J. Sci. Stat. Comput. 21 (5) (2000) 1823–1834.

25. J. E. Jones, C. S. Woodward, Newton-Krylov-multigrid solvers for large-scale,

highly heterogeneous, variably saturated ﬂow problems, Advances in Water

Resources 24 (2001) 763–774.

26. W. F. Briggs, V. E. Henson, S. F. McCormick, A Multigrid Tutorial, 2nd.

Edition, SIAM, Philadelphia, PA, 2000.

Part IV

Mathematics and Computer Science

Transport Approximations

in Partially Diﬀusive Media

Guillaume Bal

Department of Applied Physics and Applied Mathematics, Columbia University,

New York NY, 10027, USA

gb2030@columbia.edu

Abstract. This paper concerns the analysis of approximations of transport equa-

tions in diﬀusive media. Firstly, we consider a variational formulation for the ﬁrst-

order transport equation that has the correct diﬀusive behavior in the limit of small

mean free paths. The associated bilinear form is shown to be coercive on a classical

Hilbert space in transport theory with a constant of coercivity independent of the

mean free path. This allows us to obtain the diﬀusion approximation of transport

as an orthogonal projection onto a subspace of functions that are independent of

the angular variable. Similarly, projections onto functions that are independent of

the angular variable only in subsets of the full domain can be interpreted as a

transport-diﬀusion coupling method. Convergence results based on averaging lem-

mas and error estimates are presented. Secondly, we address the problem of ex-

tended non-scattering layers or ﬁlaments surrounded by highly scattering media

and derive generalized diﬀusion equations to model transport in such geometries.

1 Introduction

The solution of linear transport equations of Boltzmann type describing the

phase-space density of particles, is of interest in many applications ranging

from neutron transport in nuclear reactor physics [14] and photon trans-

port in human tissues to wave propagation in highly heterogeneous me-

dia [10,18,31]. Among the many deterministic methods available [23], several

numerical schemes have been obtained recently by using a variational frame-

work. Variational methods have a long history in transport theory, both in the

analytical and numerical approaches of source term and eigenvalue transport

problems. A very short and incomplete list of bibliographical references in-

cludes [8,11,14,30,32]. Although variational methods that involve symmetric

positive deﬁnite forms have been known for a long time in the even-parity for-

mulation of transport (see e.g. [1,23,27]), their derivations for the ﬁrst-order

Boltzmann equation are more recent [7,24–26].

Among the variational formulations involving a symmetric positive def-

inite form over a Hilbert space deﬁned in [7, 25], some enjoy the property

that the constant of coercivity is independent of the transport mean free

path, which measures (locally) the main distance between successive inter-

actions of the particles with the underlying media. This property ensures

that discretizations of solutions of the variational problem are accurate in

374 G. Bal

the transport regime, characterized by mean free paths comparable to the

typical length of variation of the geometrical components of the transport

equation, as well as in the diﬀusive regime, characterized by much smaller

mean free paths. For references on the approximation of transport by diﬀu-

sion, see e.g. [12, 19].

We consider in this paper the SAAF (Self-adjoint angular ﬂux) method

presented in [26] and characterized explicitly in [7] as a variational method

equivalent to the solution of the ﬁrst-order transport equation. We show that

the symmetric bilinear form associated to the variational problem has indeed

a constant of coercivity independent of the mean free path over a Hilbert

space that is also deﬁned independently of the mean free path.

One of the main advantages of symmetric variational formulations is that

approximations of the solutions can be obtained by orthogonal projection

(with respect to the bilinear form and onto a subspace of the considered

Hilbert space). This is analyzed in detail in [24, 25] in the framework of ﬁrst-

order system least squares (FOSLS). In this paper we show that the diﬀusion

approximation of transport may be obtained as such an orthogonal projection

(this already appears in [26]) and show that the error between the transport

and the diﬀusion equations are linear in the mean free path (and quadratic

in certain circumstances) as expected. The use of the variational formulation

in conjunction with averaging lemmas allows us to obtain that the above

error converges to zero in the limit of small mean free paths even when the

limiting solution is not very regular. The advantage of the proposed method

is that there is no need to construct any scaling transformation as in [24] or

adding any terms accounting for boundary conditions as in [25] as both can

be deduced directly from the variational formulation.

In many applications, the mean free path will be small (compared to vari-

ations of the geometrical environment) in certain areas but not necessarily

in the whole domain. A possible numerical approach consists then of solving

the diﬀusion approximation where it is valid and the transport equation else-

where. The main diﬃculty is then to couple both models at their common

interface. Such a coupling was considered in the even-parity formulation of

transport in [5]. See [16, 22, 33, 34] for additional references on the coupling

problem. We consider here the transport-diﬀusion coupling as the orthogonal

projection of the exact transport solution onto the subspace of the underly-

ing Hilbert space of functions that depend only on the spatial variable in the

diﬀusive domain but depend on both the spatial and angular variables in the

rest of the domain. We also show the convergence of the coupled transport-

diﬀusion solution (with appropriate error estimates) to the exact transport

solution when diﬀusion is valid on the domain treated by the diﬀusion equa-

tion. It is interesting to observe that the ﬁrst-order transport equation is not

exactly satisﬁed on the transport area (as opposed to the case where the full

domain is modeled by transport as mentioned earlier; see [7]). Rather the

coupling involves a second-order transport equation in the transport area.

Transport Approximations in Partially Diﬀusive Media 375

This exempliﬁes the inﬂuence of the boundary and interface conditions in

the variational formulation of ﬁrst-order transport.

Many more details on the full discretization of the transport equation

using variational approaches can be found in [24–26]. We do not consider the

discretization problem here except for a brief remark in Sect. 2.7.

For general geometries of embedded domains where the diﬀusion approx-

imation does not hold, the coupling mentioned above may be a useful alter-

native to solving the full transport equations. In certain situations however,

a more macroscopic model can be derived. Initiated by the analysis of clear

layers in optical tomography [13], generalized diﬀusion models have been

derived in [4, 6] to account for the propagation of photons along straight

lines in non-scattering clear layers. In this paper we generalize the analysis

to non-scattering ﬁlaments in three dimensional geometries and obtain new

generalized diﬀusion models. Possible applications in which similar models

may be useful include γ ray propagation in astrophysics and radiation in

atmospheric clouds.

The rest of the paper is structured as follows. Section 2 presents the vari-

ational formulation of the ﬁrst-order transport equation. General symmetric

scattering operators are considered in the two-dimensional model only where

the decomposition over spherical harmonics is particularly simple as shown

in Sect. 2.1. The variational formulation in the diﬀusive regime and its main

properties are presented in Sect. 2.2. The diﬀusion approximation is derived

in Sect. 2.3 and analyzed in the spherical harmonics expansion in Sect. 2.4.

Convergence results and error estimates are given in Sects. 2.5 and 2.6. The

transport-diﬀusion coupling is addressed in Sect. 3. Seen as an orthogonal

projection in Sect. 3.1 the resulting local equations are show in Sect. 3.2.

Finally error estimates are given in Sect. 3.3. The derivation of generalized

diﬀusion models to account for thin non-scattering inclusions is addressed in

Sect. 4.

2 Variational Formulation for Transport

We start with the steady-state transport equation written in the form

ω ·∇u + Gu = q, X = Ω × V

u = g, Γ

−

= {(x, ω) ∈ ∂Ω × V, ω · ν(x) < 0} .

(1)

Here u(x, ω) is the particle density in phase space, q(x, ω) is a volume source

term, g(x, ω) is a boundary source term, Ω ⊂ R

for n =2orn = 3 is the

spatial domain, ν(x) is its outward unit normal at x ∈ ∂Ω,andV = S

n−1

the unit sphere in the monogroup-approximation of transport. As usual Γ

−

denotes the set of incoming conditions. The operator G is deﬁned as

Gu(x, ω)=σ(x)u(x, ω) −



n−1

k(x, ω



− ω)u(x, ω



)dµ(ω



) . (2)

376 G. Bal

The integration measure is the normalized Lebesgue measure so that

µ(S

n−1

) = 1. The total absorption σ(x) and the scattering coeﬃcient

k(x, ω



− ω) are positive bounded functions. The scattering coeﬃcient sat-

isﬁes k(x, v)=k(x, −v) and is suﬃciently small so that G is a symmetric

positive deﬁnite operator on L

n−1

) with bounded inverse G

−1

We now consider the variational formulation presented in [7] and men-

tioned in [26]. We ﬁrst recast the equation as

−1

(ω ·∇u)+u = G

−1

(q),X

u = g, Γ

−

(3)

We now multiply the ﬁrst equation above by ω·∇v for a smooth test function

v(x, ω) and integrate over X to get:



−1

(ω ·∇u)ω ·∇vdp−



ω ·∇uv dp +



∂Ω×V

uvω · ν dq



−1

(q)ω ·∇vdp.

Here, dp = dxdµ(ω)anddq = dσ(x)dµ(ω), where dσ(x)isthesurfacemea-

sure on ∂Ω. Upon using the transport equation (1) we recast the above

equality as ﬁnding u ∈ W such that

a(u, v)=L(v), ∀v ∈ W, (4)

where

a(u, v)=





−1

(ω ·∇u)ω ·∇v + Guv



dp +



uvω · ν dq ,

L(v)=





−1

(q)ω ·∇v + qv



dp +



−

gv|ω · ν| dq

W =



u ∈ L

(X), ω ·∇u ∈ L

(X),u

∈ L

(Γ

; |ω · ν|dq)



(5)

The set Γ

= {(x, ω) ∈ ∂Ω × V, ω · ν(x) > 0}. It is easy to verify that

with our assumptions on G, a(u, v) is a continuous, coercive, and symmetric

bilinear form on the Hilbert space W (equipped with its natural norm) [12]

and that L is continuous in the same sense. The above equation (4) thus

admits a unique solution by the Lax-Milgram theory.

2.1 Harmonic Decomposition

It is interesting to analyze the above variational formulation using the spher-

ical harmonics decomposition. We concentrate on the two dimensional case

n = 2 so that the velocity space V = S

is the unit circle. Generalizations to

n = 3 using “classical” spherical harmonics are straightforward though more

tedious computationally [25,26]. In this setting, directions are parameterized

by ω =(cosθ, sin θ)for0≤ θ<2π. We identify u(ω) ≡ u(θ).