Optimal prediction for moment models

Optimal prediction for moment models: Crescendo diffusion and reordered equations

Benjamin Seibold Massachusetts Institute of Technology
77 Massachusetts Avenue
Cambridge, MA 02139
http://www-math.mit.edu/~seibold
and  Martin Frank TU Kaiserslautern, Erwin Schrödinger Strasse
D-67663 Kaiserslautern
Germany
http://www.mathematik.uni-kl.de/~frank
Abstract.

A direct numerical solution of the radiative transfer equation or any kinetic equation is typically expensive, since the radiative intensity depends on time, space and direction. An expansion in the direction variables yields an equivalent system of infinitely many moments. A fundamental problem is how to truncate the system. Various closures have been presented in the literature. We want to generally study moment closure within the framework of optimal prediction, a strategy to approximate the mean solution of a large system by a smaller system, for radiation moment systems. We apply this strategy to radiative transfer and show that several closures can be re-derived within this framework, such as , diffusion, and diffusion correction closures. In addition, the formalism gives rise to new parabolic systems, the reordered equations, that are similar to the simplified equations. Furthermore, we propose a modification to existing closures. Although simple and with no extra cost, this newly derived crescendo diffusion yields better approximations in numerical tests.

Key words and phrases:
radiative transfer, method of moments, optimal prediction, diffusion approximation, crescendo diffusion, reordered equations
2000 Mathematics Subject Classification:
85A25, 78M05, 82Cxx
The authors thank Martin Grothaus for helpful suggestions on measures in function spaces. The support by the German Research Foundation and the National Science Foundation is acknowledged. M. Frank was supported by DFG grant KL 1105/14/2. B. Seibold was partially supported by NSF grant DMS–0813648.

1. Introduction

In many fields, macroscopic equations can be derived from mesoscopic kinetic equations. For instance, in the Navier-Stokes and Euler equations, the macroscopic fluid variables, e.g. density and momentum, are moments of the phase space distribution of the Boltzmann equation. Similarly, in the equations of radiative transfer [29], the direction dependent kinetic equations can be transformed into a coupled system for the infinitely many moments. This can be interpreted as considering an (infinite) expansion in the kinetic variable and then taking only finitely many members of this expansion. This usually means that several “coordinates” of the field are neglected.

A common feature of existing closure strategies is that they are based on truncating and approximating the moment equations, and observing to which extent the solution of the approximate system is close to the true solution. The approximations are supported by physical arguments, such as higher moments being small or adjusting instantaneously. In Sect. 2, we derive a moment model for radiative transfer, and in Sect. 3, we outline some classical linear closure strategies.

A more potent strategy would be to have an identity for the evolution of the first moments, and then deriving closures by approximating this identity. Moment models have first been studied from a higher perspective by Levermore [28]. Recent systematic studies of moment systems include also Struchtrup’s order of magnitude method [35] which leads to the R13 equations of gas dynamics [36, 37]. In this paper, we take a similar approach. Extending results of [21], we show that the method of optimal prediction [13, 10, 12, 8, 9] can be applied to the equations of radiative transfer and yields closed systems of finitely many moments. Optimal prediction, outlined in Sect. 4, approximates the mean solution of a large system by a smaller system, by averaging the equations with respect to an underlying probability measure. It can be understood as removing undesired modes, but in an averaged fashion, instead of merely neglecting them.

Optimal prediction has been formulated for Hamiltonian partial differential equations [13, 14], however, without exploiting the full formalism. It has been applied to partial differential equations [15, 2], however, only after reducing them to a system of ordinary differential equations using a Fourier expansion or a semi-discretization step. In addition, most considered examples are Hamiltonian, for which a canonical measure exists. In contrast, here we encounter partial differential equations (in space) after a Fourier expansion (in the angular variable). Hence, the methodology is generalized to semigroups. Furthermore, the radiation system is linear. Using Gaussian measures, the formalism is linear, and it yields an identity in the presence of a memory kernel. Here we restrict ourselves to Gaussian measures with vanishing covariance. We present linear optimal prediction in function spaces, and show various possible approximations of the memory term.

In Sect. 5, we apply linear optimal prediction to the radiation moment system, and derive existing and propose new closure relations. The new formalism allows a better understanding of the errors done due to the truncation of the infinite system. The performance of the newly derived closures is investigated numerically in Sect. 6.

2. Moment Models for Radiative Transfer

The equations of radiative transfer [29] are

 1c∂tI(x,Ω,t)+Ω∇I(x,Ω,t)+(σ(x)+κ(x))I(x,Ω,t)=σ(x)4π∫4πI(x,Ω′,t)dΩ′+q(x,t). (1)

In these equations, the radiative intensity can be viewed as the the number of particles at time , position , traveling into direction . Equation (1) is a mesoscopic phase space equation, modeling absorption and emission (-term), scattering (-term), and containing a source term . Due to the large number of unknowns, a direct numerical simulation of (1) is very costly. Often times only the lowest moments of the intensity with respect to the direction are of interest. Moment models attempt to approximate (1) by a coupled system of moments.

For the sake of notational simplicity, we consider a slab geometry. However, all methods presented here can be easily generalized. Consider a plate that is finite along the -axis and infinite in the and directions. The system is assumed to be invariant under translations in and and in rotations around the -axis. In this case the radiative intensity can only depend on the scalar variable and on the azimuthal angle between the -axis and the direction of motion. Furthermore, we select units such that . The system becomes

 ∂tI(x,μ,t)+μ∂xI(x,μ,t)+(σ(x)+κ(x))I(x,μ,t)=σ(x)2∫1−1I(x,μ′,t)dμ′+q(x,t) (2)

with , , and . The system is supplied with boundary conditions that either prescribe ingoing characteristics or are periodic boundary conditions

 {I(a,μ,t)|μ>0=0I(b,μ,t)|μ<0=0or{I(a,μ,t)=I(b,μ,t)Ix(a,μ,t)=Ix(b,μ,t)

and initial conditions

 I(x,μ,0)=˚I(x,μ).

Under very general assumptions, this problem admits a unique solution , i.e. at each time is square-integrable in both spatial and angular variables [16].

Macroscopic approximations can be derived using angular basis functions. Commonly used are Legendre polynomials [17, 7], since they form an orthogonal basis on . Define the moments

 Il(x,t)=∫1−1I(x,μ,t)Pl(μ)dμ.

Multiplying (2) with and integrating over gives

 ∂tIk(x,t)+∂x∫1−1μPk(μ)I(x,μ,t)dμ+(σ(x)+κ(x))Ik(x,t)=(σ(x)I0(x,t)+2q(x,t))δk,0.

Using the recursion relation for the Legendre polynomials leads to

 ∂tIk+∂x∞∑l=0bklIl={κ(2q−I0),k=0−(κ+σ)Ik,k>0

where . This is an infinite tridiagonal system for

 ∂tIk+bk,k−1∂xIk−1+bk,k+1∂xIk+1=−ckIk+qk, (3)

of first-order partial differential equations. Using the (infinite) matrix notation

 B=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝011302325035370⋱⋱⋱⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,C=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝κκ+σκ+σ⋱⎞⎟ ⎟ ⎟ ⎟ ⎟⎠andq=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝2κq0⋮⋮⎞⎟ ⎟ ⎟ ⎟ ⎟⎠ (4)

we can write (3) as

 ∂tI+B⋅∂xI=−C⋅I+q. (5)

The infinite moment system (3) is equivalent to the transfer equation (2), since we represented an function in terms of its Fourier components. In order to admit a numerical computation, (3) has to be approximated by a system of finitely many moments , i.e. all modes with are not considered. Mathematically, it is evident that such a truncation approximates the full system, since the decay faster than , due to

 ∥I(x,⋅,t)∥2L2(−1,1)=∞∑l=02l+12Il(x,t)2<∞.

3. Moment Closure

In order to obtain a closed system, in the equation for , the dependence on has to be eliminated. A question of fundamental interest is how to close the moment system, i.e. by what to replace the dependence on . In the following, we name three types of linear closure approaches: the closure, higher-order diffusion (correction) approximations, and the simplified closure.

3.1. Pn closure

The simplest closure, the so-called closure [6] is to truncate the sequence , i.e.  for . The physical argument is that if the system is close to equilibrium, then the underlying particle distribution is uniquely determined by the lowest-order moments. This can be justified rigorously by an asymptotic analysis of Boltzmann’s equation [5].

3.2. Diffusion correction closures

The classical diffusion closure is defined for . We assume to be quasi-stationary and neglect for , thus the equations read

 ∂tI0+∂xI1 =−κI0+q0 13∂xI0 =−(κ+σ)I1.

Solving the second equation for and inserting it into the first equation yields the diffusion approximation

 ∂tI0−∂x13(σ+κ)∂xI0=−κI0+q0. (6)

A new hierarchy of approximations, denoted diffusion correction or modified diffusion closure, has recently been proposed by Levermore [28]. In slab geometry, it can be derived in the following way: We assume that for . Contrary to the closure, the -st moment is assumed to be quasi-stationary. Setting yields the algebraic relation

 IN+1=−1κ+σN+12N+3∂xIN,

which, substituted into the equation for , yields an additional diffusion term for the last moment:

 ∂tIN+bN,N−1∂xIN−1−1κ+σϑN∂xxIN={κ(2q−I0),N=0−(κ+σ)IN,N>0 (7)

where . For this closure becomes the classical diffusion closure (6).

3.3. Simplified Pn closure

Other higher order diffusion approximations exist, such as the so-called simplified () equations. For the steady case, these have been derived in an ad hoc fashion [22, 23, 24] and have subsequently been substantiated via asymptotic analysis [26] and via a variational approach [38, 4]. In the time-dependent case, the simplified equations have been derived by asymptotic analysis [20]. They are a system of parabolic partial differential equations. The equations are identical to the diffusion approximation. The equations are

 ∂tϕ =13σt∇2x[ϕ+2ϕ2−ζ]−σaϕ+q ∂tϕ2 =13σt∇2x[215αϕ+1121αϕ2]−13ασtϕ2 ∂tζ =13σt∇2x[ϕ+2ϕ2+(125(1−α)−1)ζ]−σaϕ+q−σtζ.

In these equations, is a free parameter. To obtain a well-posed parabolic system, it should be set [20] to . There are two particularly sensible choices: For , the equations reduce to the steady-state equations, whereas for , they are asymptotically close to the time-dependent equations. In the following, we set to .

These equations can be approximated by a quasi-steadiness assumption. This gives the simplified simplified () approximation

 ∂tϕ =13σt∇2x[ϕ+2ϕ2]−σaϕ+q ∂tϕ2 =13σt∇2x[215αϕ+1121αϕ2]−13ασtϕ2.

3.4. Other types of closures

Various nonlinear approximations exist in the literature, most prominently flux-limited diffusion [27] and minimum entropy methods [31, 1, 18, 39, 19].

4. Optimal Prediction

Optimal prediction, as introduced by Chorin, Hald, Kast, Kupferman et al. [13, 10, 12, 8, 9] seeks the mean solution of a time-dependent system, when only part of the initial data is known. It is assumed that a measure on the phase space is available. Fundamental interest lies in nonlinear systems, for which the mean solution decays to a thermodynamical equilibrium. Applications include molecular dynamics [34] and finance [32]. The formalism has been developed in detail [10, 11, 9] for dynamical systems

 ˙x(t)=R(x(t)) ,x(0)=˚x. (8)

Let the vector of unknowns be split into a part of interest , for which the initial conditions are known, and a part to be “averaged out” , for which the initial conditions are not known or not of relevance. In addition, let a probability measure be given (for Hamiltonian systems this could be the grand canonical distribution ). Given the known initial conditions , the measure induces a conditioned measure for the remaining unknowns. An average of a function with respect to is the conditional expectation

 Pu=E[u|^x]=∫u(^x,~x)f(^x,~x)d~x∫f(^x,~x)d~x. (9)

It is an orthogonal projection with respect to the inner product , which is defined by the expectation . Let denote the solution of (8), for the initial conditions . Then optimal prediction seeks for the mean solution

 (10)

A possible, but computationally expensive approach to approximate (10) is by Monte-Carlo sampling, as presented in [12]. Optimal prediction formulates a smaller system for . First order optimal prediction [15] constructs this system by applying the conditional expectation (9) to the original equation’s (8) right hand side. For Hamiltonian systems, the arising system is again Hamiltonian [9].

An approximate formula for the mean solution can be derived by applying the Mori-Zwanzig formalism [30, 40] in a version for conditional expectations [10] to the Liouville equation for (8). It yields an integro-differential equation, that involves the first order right hand side, plus a memory kernel. First order optimal prediction can be interpreted as the crude approximation of dropping the memory term. Various better approximations have been presented in [11, 12, 8]. For nonlinear systems, optimal prediction remains an approximation, even if the memory term is computed exactly.

4.1. Linear Optimal Prediction

Consider a linear system of evolution equations (such as (5))

 ∂tu=Ru ,u(0)=˚u, (11)

where is a differential operator. Let the unknowns and the operator be split

 u=[^u~u]andR=[^^R^~R~^R~~R].

Lebesgue measures cannot be generalized to infinite dimensions. Instead, as in [25, 3], we define measures on function spaces as follows. Let be a space of test functions (e.g. the Schwartz space) and let be the corresponding space of distributions, such that is a Gelfand triple [3]. Measures on vector valued function spaces are defined via a vector valued Gelfand triple .

Let be a vector-valued test function, a vector-valued distribution, and a matrix-valued distribution. Then the pair-dot-product

 ⟨f,u⟩SS′=∑j⟨fj,uj⟩SS′∈R

is a number, and the pair-matrix-vector-product

 ⟨f,Σ⟩SS′=⎛⎜ ⎜ ⎜⎝∑j⟨fj,Σj1⟩SS′⋮∑j⟨fj,Σjn⟩SS′⎞⎟ ⎟ ⎟⎠∈Rn

is a vector.

We define a Gaussian measure by its Fourier transform

 ∫exp(i⟨f,u⟩SS′)dλ(u)=exp(−∥⟨f,Σ⟩SS′∥222+i⟨f,m⟩SS′), (12)

where the variance is a matrix-valued function and the expectation value is a vector-valued distribution.

Here we assume a measure with vanishing covariance between the different moments (i.e. diagonal ). As derived in [21], this leads to the simple linear projections

 Px=Exand(I−P)x=Fx

with

Considering the solution operators and , which we assume to be well posed, the Mori-Zwanzig formalism [30, 40] yields the identity

 ∂t(etR)=RetR+etRFRF+∫t0K(t−s)esRds, (13)

where

 R=RE (14)

is the projected differential operator, and

 K(t)=etRFRFRE (15)

is a memory kernel for the dynamics. The projected operator represents the mean solution in the context of optimal prediction. Its evolution is given by

 ∂t(etRE)=RetRE+∫t0K(t−s)esREds. (16)

Compared to (13), the middle term has canceled, since . The projected differential operators, the projected right hand side (14) and the memory kernel (15) have the block form

Consequently, in (16) the evolution of the resolved moments is decoupled from the unresolved moments. However, the resolved part of the memory kernel involves the orthogonal dynamics , whose solution is as complex as the full system (11). In the following, we derive new closures, by approximating the memory term in two steps. First, the integral in (16) is approximated using a quadrature rule. This is a general step that is independent of the specific form of the differential operator . We present various approximations in Sect. 4.2. Second, the exact orthogonal dynamics are approximated by a simpler evolution operator. Good approximations employ the specific structure of . We present such an approximation in Sect. 5.

4.2. Approximations of the Memory Integral

The simplest approximation of the memory integral is to drop it. This yields the first order optimal prediction system

 ∂tu=Ru, (17)

which in some applications yields satisfying results. However, in most cases better results can be expected when approximating the memory term more accurately. This approach is generally denoted second order optimal prediction. Assume for the moment that the memory kernel be exactly accessible. The memory term in (16) uses the solution at all previous times. In principle, one could store the solution at all previous times, and thus approximate the integral very accurately (as presented in [33]). However, such approaches are typically very costly in computation time and memory. In addition, for partial differential equations, the stability of such integro-differential equations (or delay equations, after approximating the integral), is a very delicate aspect. Therefore, here we focus on approximation that only use the current state of the system. We present three approximations.

• Piecewise constant quadrature: Assume the memory kernel decays with a characteristic time scale . Then a piecewise constant quadrature rule for yields

 ∫t0K(s)u(t−s)ds≈∫τ0K(0)u(t)ds=τK(0)u(t). (18)
• Truncated piecewise constant quadrature: For , the memory integral in (18) cannot reach over the full length . Hence, it is reasonable to approximate

 ∫t0K(s)u(t−s)ds≈min{τ,t}K(0)u(t). (19)
• Modified piecewise constant quadrature: Using a piecewise constant quadrature rule for only, yields

 ∫t0K(s)u(t−s)ds≈∫t0K(s)u(t)ds=∫t0K(s)dsu(t). (20)

5. Application to the Radiation Moment System

We now turn our attention to the infinite moment system (3). For consistency with the notation used in Sect. 4.1, we denote the (infinite) vector of moments by . In addition, we neglect the forcing term, since it is unaffected by any truncation. The radiation system (3) can be written as

 ∂tu=Ru, (21)

where the differential operator involves the (infinite) matrices (4). We consider a splitting of the moments , where contains the moments resolved. The radiative intensity is always the starting moment. The other moments, however, could be reordered. The system splits into blocks

 u=[^u~u],B=[^^B^~B~^B~~B],andC=[^^C001τI],

where , which we also assume to be the characteristic time scale of the system ( is a time scale since we have set ). The optimal prediction system (16) involves the block operators

 RE=−[^^B∂x+^^C0~^B∂x0],RF=−[0^~B∂x0~~B∂x+1τI],andRFRE=⎡⎣^~B~^B∂xx0~~B~^B∂xx+1τ~^B∂x0⎤⎦.

Here the upper left block of is the projected right hand side, which yields the first order optimal prediction system. The operator describes the orthogonal dynamics system

 {∂t^u=−^~B∂x~u∂t~u=−~~B∂x~u−1τ~u. (22)

Its solution operator is contained in the memory kernel (15). The equation for in (22) is similar to the original full system (21), however, with one fundamental simplification. The decay matrix is a multiple of the identity, a fact that may be employed for fast approximations of the orthogonal dynamics. We shall pursue this approach in future work. For now, we construct a simple approximate system by neglecting the advection in the unresolved variables, i.e. . The approximate orthogonal dynamics operator

 G=−[0^~B∂x01τI]≈RF

has a simple solution

 etG[^u~u]=[^u−τ(1−e−t/τ)^~B∂x~ue−t/τ~u].

Hence, the resolved component of the approximate memory kernel is given by

 ˆˆetGRFRE =^~B~^B∂xx−τ(1−e−t/τ)^~B∂x(~~B~^B∂xx+1τ~^B∂x) =e−t/τ^~B~^B∂xx−τ(1−e−t/τ)^~B~~B~^B∂xxx.

In the following, we work out the specific expressions for two different types of ordering of the moments.

5.1. Standard ordering

We consider and . With respect to this ordering, the projected right hand side is a simple truncated version of the full right hand side. The triple matrix product vanishes . And the double matrix product equals

 ^~B~^B=⎛⎜ ⎜⎝0…0⋮⋱⋮0…ϑN⎞⎟ ⎟⎠,

where . Hence, with the above simplified orthogonal dynamics, the optimal prediction integro-differential equation (16) equals

 (23)

Due to the special structure of , the memory term only modifies the -th equation, thus it is in fact a closure. The various approximations to the memory integral, presented in Sect. 4.2, yield the following modifications to the equation of the -th moment:

• First order optimal prediction: When dropping the integral term, the resulting system (23) becomes exactly the closure. Note that by considering measures with non-vanishing covariances between the moments, all other linear closures can be derived [21].

• Piecewise constant quadrature: The integral term in (23) is approximated by

 τϑN∂xxuN(t). (24)

A comparison with (7) reveals that we have re-derived Levermore’s diffusion correction closure [28]. In the case , it is equivalent to the classical diffusion approximation.

• Truncated/modified piecewise constant quadrature: The approximations (19), respectively (20) lead to the integral approximation

 f(t)ϑN∂xxuN(t). (25)

Compared to the diffusion correction closure (24), the constant is replaced by a time dependent function , that increases gradually towards its maximum value . The diffusion is turned on gradually over time. Hence, we denote these types of approximations crescendo diffusion (), respectively crescendo diffusion correction () closure. More specifically, we call the approximation truncated crescendo diffusion (closure) when it is based on (19). In this case . Similarly, the approximation based on (20), that leads to , is called modified crescendo diffusion (closure).

The new crescendo closures (25) introduce an explicit time dependence in the diffusion term. The physical rationale is that at , the state of the system (at least the unknowns of interest) is known exactly. Information is lost as time evolves, due to the approximation. This implies that crescendo diffusion closures shall only be used if the initial conditions are known with higher accuracy than the approximation error due to truncation. It is one way of circumventing an initial layer in time. In Sect. 6.1, we investigate the quality of these new approximations numerically.

Although here we have assumed spatially homogeneous coefficients, we expect that the equations can be adapted to the space-dependent case in analogy to diffusion theory. Specifically, if and are space dependent, we define , and replace by . The validity of this approximation will be addressed in future research.

5.2. Reordered equations

In the following, we consider a reordering of the unknowns, given by the reordered vector , where is a permutation matrix that leaves as the first component. This reordering is inspired by Gelbard’s original derivation [22, 23, 24] of the steady simplified equations, in which the unknowns are the even-order moments and the odd-order moments are algebraically eliminated. The reordered system is

 ∂tv=−B∂xv−Cv,

where we have the permuted advection matrix . The decay term remains unchanged, since . Of particular interest is the following even-odd reordering: For a given number , consider the following permutation of the unknowns

 v=Pu=(u0,u2,…,u2N,u1,u3,…,u2N+1,u2N+2,…)T,

and consider the lowest even moments as of interest. The resulting reordered advection matrix has the following structure (here for ):

 [^^B^~B~^B~~B]=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣12/53/54/95/91/32/33/74/75/116/116/137/137/15⋱⋱⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦.

The solid lines indicate separations between and . For the even-odd reordering, we have . As for standard ordering, the triple matrix product vanishes . The double matrix product becomes more interesting than for standard ordering. For any choice of , it is observed to be positive definite. In addition, lower order matrices are contained as submatrices in higher-order ones. With the simplified orthogonal dynamics, the optimal prediction integro-differential equation (16) equals

 ∂t^v(t)=−^^C^v(t)+∫t0e−s/τD∂xx^v(t−s)ds. (26)

The various approximations to the memory integral, presented in Sect. 4.2, now become:

• First order optimal prediction: Dropping the integral term in (26) yields a mere exponential decay for each moment

 ∂t^v=−^^C^v,

which is obviously a horrible approximation. For even-odd reordering, a non-trivial approximation to the integral term is strictly required.

• Piecewise constant quadratures: System (26) is approximated by

 ∂t^v=−^^C^v+f(t)D∂xx^v, (27)

where the coefficient function may be constant or time-dependent, depending on the integral approximation used. When based on (18), we have . In this case, we denote system (27) the reordered equations (). When the approximation (19) is used, we have , and speak of the truncated crescendo equations. Similarly, when (20), we have , and speak of the modified crescendo equations. In the case , these systems become the corresponding diffusion closures (6).

The system (27) has strong similarity with the time dependent equations [20]. These are derived via asymptotic analysis. Although the unknowns are not directly comparable, a comparison of the respective diffusion matrices is compelling. For the choice of the diffusion matrix is

 D=⎛⎜⎝0.33330.6667−0.33330.13330.523800.33330.66670.2000⎞⎟⎠.

In comparison, the matrix is

 D=⎛⎜⎝0.33330.666700.13330.52380.342900.19050.5065⎞⎟⎠.

Observe that the respective upper submatrices are identical. The other entries differ, which can be expected, since the third unknown in the equations is not the fourth moment, as it is for the equations. Assuming the third variable in be quasi-steady leads to the system [20], that turns out to be identical to the system. Hence, similarly as both classical and new diffusion approximations can be derived from optimal prediction with standard ordering, even-odd ordering allows the derivation of existing and new parabolic systems that approximate the radiative transfer equations. In Sect. 6.2, we investigate the quality of these new systems numerically.

6. Numerical Results

As a test case, we consider the 1D slab geometry system (2) on (periodic boundary conditions) with , and no source . The initial conditions are and . The numerical solution is obtained using a second order staggered grid finite difference method with resolution and . Any diffusive terms are implemented by two Crank-Nicolson half steps.

6.1. Diffusion Correction Approximations

The time evolution of the radiative intensity for various system sizes is shown in Fig. 1 for , Fig. 4 for , and Fig. 5 for . In each case, the thick solid curve is the solution to the full radiative transfer equation, obtained using a closure, which has converged in (). The dotted function is the classical closure (Sect. 3.1). The dashed curve is the diffusion correction closure (7), which for equals the diffusion closure (6). The dash-dotted curve and the thin solid curve are the two versions of the crescendo diffusion correction (25), as derived in Sect. 5.1. The dash-dotted curve denotes the sharp time dependence, while the solid curve denotes the exponential time dependence.

For , the closure gives a sole exponential decay of the initial conditions, which is a very bad approximation. The diffusion closure is significantly better than , but too diffusive. Both versions of crescendo diffusion remedy this problem fairly well, and yield good approximations, especially for short times. Note that here the truncated crescendo diffusion reaches its full diffusivity at .

For , the solution splits into peaks which move sideways with different velocities. This is due to the fact that the model is a linear system of hyperbolic equations with fixed wave speeds. With increasing, the height of the peaks goes to zero, and the solutions approach the true solution. The diffusive closures smear out the highest considered moment , thus making it more uniform and reducing its influence on the second-highest moment . Hence, qualitatively an -th order diffusion correction solution lies “between” the and the solution. The results shown in Fig. 4 and Fig. 5 visualize this effect. In particular, one can observe that the classical diffusion correction is too close to , i.e. the diffusion in the -th moment is too large. In contrast, the crescendo diffusion correction approximations ameliorate this effect. The peak height is reduced, and the approximate solutions are closer to the truth.

The quality of the various approximations can be quantified by considering the time evolution of the -error of the approximate radiative intensity with respect to the true solution

 e(t)=∫10uapprox0(x,t)−utrue0(x,t)dx.

The error of the various approximations is shown in Fig. 4 for , in Fig. 4 for , and in Fig. 8 for . One can observe that the errors become smaller as increases. Within each figure, the closure yields the largest error at the final time . In contrast, the diffusion (correction) closures show smaller errors at . However, for short times, the error is increased compared to the closure. This initial layer in time is particularly evident for , shown in Fig. 4.

Both versions of crescendo diffusion (correction) remove the initial layer of the diffusion closure, and the following increase in error is significantly smaller than with the closure. In terms of accuracy, no significant difference can be observed between the two versions.

6.2. Parabolic Approximations

Fig. 8 shows a similar comparison for parabolic type approximations, as derived in Sect. 5.2. The thick solid curve is the true solution. The dotted function is the classical diffusion approximation, which equals the system. The dashed curve shows the system (27), and the dash-dotted curve shows the system. For increasing order, the reordered equations converge to some limiting solution that is not equal to the true solution. While this is easily explained by the fact that the odd moments are not resolved, it implies that these systems are mainly useful for small values of . For comparison, the thin solid curve shows the system (with . One can observe that both the diffusion approximations and the equations are too diffusive. The equations are far from perfect, but they do yield better approximations.

The time evolution of the errors is shown in Fig. 8. For the considered examples, the errors obtained by the equations are slightly smaller than the errors obtained by the diffusion closure or by the .

7. Conclusions and Outlook

We have applied the method of optimal prediction to the infinite moment systems that model radiative transfer. An integro-differential identity for the evolution of a finite number of moments is obtained. The memory kernel involves an orthogonal dynamics evolution operator. Simple approximations of this operator and of the integral have been presented, and various approximate systems been derived. While traditionally closures had been derived using physical arguments or by asymptotic analysis, the optimal prediction formalism provides a very different strategy in which closures are derived from a mathematical identity.

Using this methodology, existing closures can be re-derived, such as classical closures, diffusion and diffusion correction closures. In addition, new closures have been obtained. Fundamentally new aspects are crescendo diffusion correction closures, that modify classical diffusion correction closures by turning on the diffusion gradually with time. In the context of optimal prediction, this explicit time dependence has a natural interpretation as loss of information. While crescendo diffusion comes at no additional cost, numerical tests indicate that the quality of the results is generally improved.

In addition, the application of the formalism to a reordered version of the moment system yields approximations of parabolic nature, which we denote reordered equations. Some versions of the existing simplified equations can be re-derived, and other systems can be obtained. A pleasing feature of of these new equations is that the diffusion matrices are always positive definite.

Future research directions will focus on better approximations of the orthogonal dynamics operator, as well as of the memory integral. Approximations of the orthogonal dynamics can be based on the fact that higher order moments possess a uniform decay rate. The memory integral can be approximated using the solution at former times, thus leading to delay equations. While such systems pose numerical challenges, it may be very rewarding to significantly improve approximations for very few considered moments. Another limitation of the presented analysis lies in its linearity. Nonlinear approximations, such as flux-limited diffusion [27] and minimum entropy closures [31, 1, 18] would require nonlinear projections. While the optimal prediction formalism is not restricted to linear projections, a precise formulation in the case of partial differential equations is yet to be presented.

References

• [1] A. M. Anile, S. Pennisi, and M. Sammartino, A thermodynamical approach to Eddington factors, J. Math. Phys. 32 (1991), 544–550.
• [2] J. Bell, A. J. Chorin, and W. Crutchfield, Stochastic optimal prediction with application to averaged Euler equations, Proc. 7th Nat. Conf. CFD, 2000, pp. 1–13.
• [3] Y. M. Berezansky and Y. G. Kondratiev, Spectral methods in infinite-dimensional analysis, Kluwer Academic Publishers, Dordrecht, 1995.
• [4] P. S. Brantley and E. W. Larsen, The simplified approximation, Nucl. Sci. Eng. 134 (2000), 1.
• [5] R. E. Caflisch, The fluid dynamic limit of the nonlinear Boltzmann equation, Commun. Pure Appl. Math. 33 (1980), 651–666.
• [6] S. Chandrasekhar, On the radiative equilibrium of a stellar atmosphere, Astrophys. J. 99 (1944), 180.
• [7] by same author, Radiative transfer, Dover, 1960.
• [8] A. J. Chorin, Conditional expectations and renormalization, Multiscale Model. Simul. 1 (2003), 105–118.
• [9] A. J. Chorin and O. H. Hald, Stochastic tools in mathematics and science, Springer, 2006.
• [10] A. J. Chorin, O. H. Hald, and R. Kupferman, Optimal prediction and the Mori-Zwanzig representation of irreversible processes, Proc. Natl. Acad. Sci. USA 97 (2000), 2968–2973.
• [11] by same author, Non-Markovian optimal prediction, Monte Carlo Meth. Appl. 7 (2001), 99–109.
• [12] A. J. Chorin, O. H. Hald, and Kupferman R., Optimal prediction with memory, Physica D 166 (2002), 239–257.
• [13] A. J. Chorin, A. P. Kast, and Kupferman R., Optimal prediction of underresolved dynamics, Proc. Natl. Acad. Sci. USA 95 (1998), 4094–4098.
• [14] by same author, Unresolved computation and optimal predictions, Comm. Pure Appl. Math. 52 (1998), 1231–1254.
• [15] by same author, On the prediction of large-scale dynamics using unresolved computations, Contemp. Math. AMS 238 (1999), 53–75.
• [16] R. Dautray and J. L. Lions, Mathematical analysis and numerical methods for science and technology, sixth ed., Springer, Paris, 1993.
• [17] B. Davison, Neutron transport theory, Clarendon Press, Oxford, 1958.
• [18] B. Dubroca and J. L. Feugeas, Entropic moment closure hierarchy for the radiative transfer equation, C. R. Acad. Sci. Paris Ser. I 329 (1999), 915–920.
• [19] M. Frank, B. Dubroca, and A. Klar, Partial moment entropy approximation to radiative heat transfer, J. Comput. Phys. 218 (2006), 1–18.
• [20] M. Frank, A. Klar, E. W. Larsen, and S. Yasuda, Time-dependent simplified pn approximation to the equations of radiative transfer, J. Comput. Phys. 226 (2007), 2289–2305.
• [21] M. Frank and B. Seibold, Optimal prediction for radiative transfer: A new perspective on moment closure, submitted (2008).
• [22] E. M. Gelbard, Applications of spherical Harmonics method to reactor problems, Tech. Report WAPD-BT-20, Bettis Atomic Power Laboratory, 1960.
• [23] by same author, Simplified spherical harmonics equations and their use in shielding problems, Tech. Report WAPD-T-1182, Bettis Atomic Power Laboratory, 1961.
• [24] by same author, Applications of the simplified spherical Harmonics equations in spherical geometry, Tech. Report WAPD-TM-294, Bettis Atomic Power Laboratory, 1962.
• [25] T. Hida, H.-H. Kuo, J. Potthoff, and L. Streit, White noise. An infinite dimensional calculus, Kluwer Academic Publisher, Dordrecht, Boston, London, 1993.
• [26] E. W. Larsen, J. E. Morel, and J. M. McGhee, Asymptotic derivation of the multigroup and simplified equations with anisotropic scattering, Nucl. Sci. Eng. 123 (1996), 328.
• [27] C. D. Levermore, Relating Eddington factors to flux limiters, J. Quant. Spectrosc. Radiat. Transfer 31 (1984), 149–160.
• [28] by same author, Transition regime models for radiative transport, Presentation at IPAM: Grand challenge problems in computational astrophysics workshop on transfer phenomena, 2005.
• [29] M. F. Modest, Radiative heat transfer, second ed., Academic Press, 1993.
• [30] H. Mori, Transport, collective motion and Brownian motion, Prog. Theor. Phys. 33 (1965), 423–455.
• [31] I. Müller and T. Ruggeri, Rational extended thermodynamics, second ed., Springer, New York, 1993.
• [32] P. Okunev, A fast algorithm for computing expected loan portfolio tranche loss in the Gaussian factor model, Tech. report, Lawrence Berkeley National Laboratory, 2005.
• [33] B. Seibold, Non-Markovian optimal prediction with integro-differential equations, Tech. report, Lawrence Berkeley National Laboratory, 2001.
• [34] by same author, Optimal prediction in molecular dynamics, Monte Carlo Methods Appl. 10 (2004), no. 1, 25–50.
• [35] H. Struchtrup, Stable transport equations for rarefied gases at high orders in the Knudsen number, Phys. Fluids 16 (2004), 3921–3934.
• [36] H. Struchtrup and M. Torrilhon, Regularization of Grad’s 13 moment equations: Derivation and linear analysis, Phys. Fluids 15 (2003), 2668–2680.
• [37] by same author, Higher-order effects in rarefied channel flows, Phys. Rev. E 78 (2008), 046301.
• [38] D. I. Tomasevic and E. W. Larsen, The simplified approximation, Nucl. Sci. Eng. 122 (1996), 309–325.
• [39] R. Turpault, M. Frank, B. Dubroca, and A. Klar, Multigroup half space moment appproximations to the radiative heat transfer equations, J. Comput. Phys. 198 (2004), 363–371.
• [40] R. Zwanzig, Nonlinear generalized langevin equations, J. Stat. Phys. 9 (1973), no. 3, 215–220.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters