Graziani F. (editor) Computational Methods in Transport

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

Accurate and Eﬃcient Radiation Transport in Optically Thick Media 265

= −σ



D −

, (23)

giving the formal solution

D(s)=D(0)exp



−





)ds





−





dB(s



)



exp



−









)ds





. (24)

The formal solution shows that the boundary condition, D(0), decays in

a few absorption lengths. Deep in the material σ



varies slowly. In fact both



and dB/ds are constant to ﬁrst order in 

stream

. Equation (24) can then

be integrated. The result is

D(s)=



D(0) +





exp[−σ



s] −



. (25)

It shows that

D(s)=−



(s)

dB(s)

(26)

is the steady-state solution of (24) and that any boundary value of D(0)

decays to it in a few absorption lengths. A result that is correct to second

order in 

stream

was given in [SB03].

Asymptotic Expansion

The relative orders of various terms in the transport equation, (18), in thick

media were estimated in Sect. 2.2. The estimation is valid far from boundary

layers and time transients. To ﬁrst order in 

stream

, there are only two terms

0=−σ



D − Ω·∇B, (27)

giving the solution

D = −



Ω·∇B = −



∂B

∂T

Ω·∇(T

) . (28)

The radiation ﬂux, from (20), is

rad

= −





∞

dν



∂B

∂T



∇(T

)=−

3σ

∇(T

) . (29)

This is the correct diﬀusion limit of the transport equation. We also recovered

the correct deﬁnition of the Rosseland mean opacity, σ

; see [Cas00, MM84].

To ﬁrst order (29) is identical to (25). It conﬁrms the ﬁrst order accuracy of

the diﬀusion ﬂux [Mor00]. Note the utter simplicity of the derivation.

An expansion in higher orders of 

stream

can also be carried out. The

results are similar to those of Morel [Mor00], but they are slightly diﬀerent

and more consistent. A short discussion was given in [SB03].

266 A. Sz˝oke, et al.

3.3 The Diﬀerence Formulation, in LTE, with Scattering

The Scattering Term

The scattering term was displayed in (5). In LTE, the Planck distribution at

the material temperature is stationary and it also satisﬁes detailed balance.

This imposes thermodynamic conditions on the scattering cross sections



(ν



→ ν, Ω·Ω



)B(ν



)



B(ν)

2hν



= σ

(ν → ν



, Ω·Ω



) B(ν)



B(ν



)

2hν

3



, (30)

where the x,t; T dependence of σ

and B has been suppressed. The scat-

tering term expressed in terms of the diﬀerence ﬁeld, D, and the material

temperature, T, is obtained by substituting I = B + D in (5) and then using

(30) to rewrite it in the following form

Q(D)=



∞

dν





4π

dΩ



(ν



→ ν, Ω·Ω



)D(ν



, Ω



)



1 − e

−hν



/kT

1 − e

−hν/kT

D(ν, Ω)

2hν



−



∞

dν





4π

dΩ



(ν → ν



, Ω·Ω



)D(ν, Ω)



1 − e

−hν/kT

1 − e

−hν



/kT

D(ν



, Ω



)

2hν

3



(31)

Regrettably, there was an error in (31) in our original paper, [SB03].

If scattering does not change the radiation energy, e.g. in Thomson scat-

tering, the stimulated emission terms in (5) cancel identically and an isotropic

distribution is stationary under those conditions. In fact, we deﬁne, as usual

J(x,t; ν):=

4π



4π

dΩ I(x,t; ν, Ω) . (32)

Then I = J is stationary and (5) can be written as

mono

(I − J)=



4π

dΩ



(ν, Ω·Ω



)[I(ν, Ω



) − J(ν)]

−



4π

dΩ



(ν, Ω·Ω



)[I(ν, Ω) − J(ν)] . (33)

In the scattering terms, (5), (31), (33), both σ

(ν



→ ν, Ω·Ω



)andI(ν, Ω)

or D(ν, Ω) can be expanded in spherical harmonics [Brn95]. The integrals are

then reduced to relaxation equations for the spherical harmonic components

of D,orI − J, respectively. In particular, if the scattering is isotropic, (33)

reduces to

mono

(I − J)=σ

(ν)(I − J) . (34)

Finally, we note that the scattering terms are always proportional to

(ν)=



4π

dΩ



(ν, Ω·Ω



) , (35)

and to the analogous expressions in (5), (31).

Accurate and Eﬃcient Radiation Transport in Optically Thick Media 267

The Full Equations

The full equations in the diﬀerence formulation, in LTE, are obtained by

adding the right-hand side of (31) to the right-hand side of (17). As the

change of the radiation energy caused by scattering comes from the matter,

the integral of (31) over frequency has to be subtracted from (21), in analogy

to (8). They have to be solved together with the conservation equations, (9).

In the more general case they have to be solved together with the full set of

equations of radiation hydrodynamics [Cas00], [MM84, Pom73].

Rather than presenting the formal development of the equations, we will

sketch their simpliﬁed version that gives some insight into the relaxation

behavior of the radiation intensity. We start by rewriting (1) in terms of the

path variable, s, as deﬁned in (22). We also take liberties with the scattering

term; we tacitly assume it to be monochromatic and isotropic, as in (34).

= −σ



(I − B) −σ

(I − J) . (36)

We deﬁne the total extinction coeﬃcient, σ

= σ



+ σ

. Obviously,







= −σ



(I − B) − σ

(I − J) . (37)

In analogy to the diﬀerence formulation, we subtract appropriate terms from

the equation to give



d(I − B)

d(I − J)

= −σ



(I − B) − σ

(I − J) −



−

. (38)

Depending on the relative magnitudes of σ



and σ

, it is easy to see that

asymptotically, in thick media, the equation approaches either the diﬀusion

approximation or the Eddington approximation.

An alternative rearrangement gives

= −σ



(J − B) − (σ



+ σ

)(I − J) . (39)

We now subtract dB/ds from both sides of the equation, to give

d(I − J)

d(J − B)

= −σ



(J − B) − (σ



+ σ

)(I − J) −

. (40)

With the further assumption of dB/ds = 0 and the constancy of σ



and

along the radiation path, it is easy to verify that (40) has the solution

J − B =[J(0) − B(0)]exp[−σ



s] , (41)

268 A. Sz˝oke, et al.

I − J =[I(0) − J(0)]exp[−(σ



+ σ

)s] . (42)

We note at this point that the material relaxation equation, (8), can

always be written as

∂E

mat

∂t

=4π



∞

dνσ



(J − B)+G, (43)

so J − B is a natural variable to consider.

We emphasize again that our derivation of the solution of (40) was highly

simplistic. In particular it did not take into account material relaxation and

the possibility of scattering with change of photon frequency, e.g. Compton

scattering [Cas00]. We hope to report later on further developments along

these lines.

4 Test Problems

A variety of test problems have been investigated in [Bro05] in order to eval-

uate the computational eﬃciency and accuracy of the diﬀerence formulation,

employing the Symbolic Implicit Monte Carlo (SIMC) method. Our goal was

to analyze some simple situations that indicate the potential impact of the

diﬀerence formulation in more complex physics environments. The key issues

are accuracy and eﬃciency for both thick and thin media, and speciﬁcally the

frequent occurrence of media which are thick at one frequency while being

thin at others. In addition to comparing the Monte Carlo solution of the two

formulations of transport, we compared them to analytic, or semi-analytic

solutions where they are available. The main results of that study will be

summarized below.

4.1 Behavior of a Finite Slab Heated from the Outside

A very basic test of a time-dependent thermal transport algorithm is whether

it correctly approaches the steady state solution for a ﬁnite slab immersed

in a heat bath with diﬀerent temperatures on either side. The fourth power

of the temperature should become a linear function of the optical depth in

a gray, thick medium, in steady state. Deviations from such a straight line

are indicative of boundary layers when they occur within a mean free path

or so of a surface, but otherwise indicate serious problems in the numerical

solution. Teleportation errors result in a wrong slope and cause curvature

in the interior solution [MBS03]; for this reason zone thickness were limited

to one mean free path in our piecewise constant treatment of the material

temperature.

In addition to the issues noted above, the time dependent approach to

steady state oﬀers the opportunity to check the correctness of the implemen-

tation of interior source terms as well as initial and boundary conditions. That

Accurate and Eﬃcient Radiation Transport in Optically Thick Media 269

the steady state temperature is independent of time step, reﬂecting implicit

behavior of the time integration, can also be checked. When the spectrum

of the radiation ﬁeld is examined, even for a grey opacity, the correctness of

the frequency sampling algorithm can be checked. Agreement between the

standard and diﬀerence formulation is non-trivial due to the diﬀerent nature

of the source terms.

In the presentation of our computational results below, we provide a clear

demonstration of rigorous agreement between the two formulations for trans-

port, in terms of their approach to steady state, along with a measure of

increase in computational eﬃciency for the diﬀerence formulation. The mag-

nitude of this increase in computational eﬃciency, in the form of greatly re-

duced Monte Carlo noise as the optical thickness of the problem is increased,

is somewhat surprising even to the authors who were prospecting for it.

In the ﬁrst set of simulations we considered a ﬁnite slab heated from the

left side, with an open boundary on the right, allowing the radiation ﬂowing

through the slab to enter free space. The slab was composed of a uniform,

static material having a frequency-independent (gray) opacity. We calculated

the time dependence of the temperature and of the radiation ﬁeld after a 1

keV black-body source on the left side of the slab is turned on at time t =0.

(1keV≈ 1.2 10

7 ◦

K.)

During the time dependent execution of the problem, a thermal wave,

also known as Marshak wave, sweeps the problem domain and the solution

then approaches steady state. We compared the solutions provided by the two

formulations and their relative noise, for identical problem run times, in order

to obtain a measure of the accuracy and relative computational eﬃciency of

the two formulations for transport, under conditions that the Monte Carlo

portion of the code dominates execution time.

Four instances of this problem are presented below. The slab is composed

of a uniform material having a frequency independent (gray) opacity of 0.1, 1,

10, and 100 mean free paths per cm respectively. The slab is 10 cm thick so the

four opacities correspond to total optical depths of 1, 10, 100 and 1000 mean

free paths. The speciﬁc heat of the material is a constant 0.1 jerk/(keV cm

where jerk is an energy unit (1 jerk = 10

ergs), and temperature is measured

in energy units of kT = 1 keV. The slab is initially at a temperature of

0.01 keV, with the radiation ﬁeld being a Planckian in equilibrium with this

temperature. All four problems used a time step of 0.2 sh, where 1 sh =

−8

sec. The problem 1 mean free path thick was run to 20 sh in order to

get close to steady state, as was the problem 10 mean free paths thick. The

problem 100 mean free paths thick was run to 40 sh in order to approach

steady state. The problem 1000 mean free paths heated up very slowly due

to the diﬀusive nature of the solution, requiring 320 sh in order to suitably

approach steady state. The problems 1 and 10 mean free paths thick employed

20 zones, while the thicker problems employed zones one mean free path thick

in order to prevent teleportation error from inﬂuencing the results, and to

270 A. Sz˝oke, et al.

prevent anomalous performance results for the diﬀerence formulation. Equal

thickness zones were used everywhere. Geometric zoning has a role only if

piecewise linear treatment of the material temperature is available to remove

teleportation error.

The speed of light is 300 cm/sh. This provides 6 traversals of the slab for a

time step of 0.2 sh, reﬂecting on the stability of the method. Relaxing the need

for implicit treatment of the source terms was one of the hopes of the authors,

given that implicit treatment requires the solution of a non-linear system of

equations for each time step. Our experiences in this regard, documented

in [Daf05], were made even more diﬃcult by the T

term for thermal emission.

Explicit treatment of the ∂B/∂x source terms was abandoned as a result, but

may be worth revisiting in mixed physics applications where the time step

size is limited for reasons outside of transport physics.

We would like to note that the Monte Carlo solution for the two transport

formulations have entirely diﬀerent requirements for spatial importance sam-

pling if uniform statistical noise, as a function of position, is to be obtained.

In the standard formulation most of the computational eﬀort is spent com-

puting the balance between emission and absorption that produces the local

equilibrium black body ﬁeld. As a result, a scheme that samples particles

with a uniform density in space produces a relatively ﬂat statistical noise

across the slab. The number of source particles born in each zone, during

each time step, is proportional to the thickness of the zone. The statistical

properties of the Monte Carlo solution for the diﬀerence formulation are in

sharp contrast to this. For thicker problems, the statistical noise for the dif-

ference formulation is miniscule near the hot side and increases very rapidly

towards the cold boundary. Suitable importance sampling could ﬂatten out

this growth in noise. An exposition of this topic is beyond the scope of this

paper.

In the results shown below, the unit of source sampling was one ∂B/∂t

particle, and one ∂B/∂x particle pair, for each zone, or zone interface, respec-

tively, for the diﬀerence formulation. We have not attempted to tune the ratio

between ∂B/∂t and ∂B/∂x source particles, nor have we attempted to tune

the relative importance of source sampling across the volume of the problem.

In the standard formulation, the unit of source sampling is one thermally

emitted particle per zone.

In Fig. 1 (a) we show the steady state solution for the material temper-

ature, as a function of position in the slab, for the problem instance that is

one mean free path thick. The average of 100 randomly seeded runs using

the standard formulation, and the average of the same number of instances

using the diﬀerence formulation, is shown by the boxes and the × symbols,

respectively. The two formulations are in perfect agreement. The standard

deviation of the results for the standard formulation, as well as the same

statistic for the diﬀerence formulation, both multiplied by 400, are plotted

using diamonds and triangles, respectively. The Monte Carlo source parti-

Accurate and Eﬃcient Radiation Transport in Optically Thick Media 271

Fig. 1. Temperature distribution in a slab in steady state. The slab is heated from

the left by a 1 keV black body, and it radiates freely on the right. The total optical

depth (OD) of the slab is 1 in (a), 10 in (b), 100 in (c) and 1000 in (d). The

standard deviation of the standard formulation is denoted by diamonds and that of

the diﬀerence formulation by triangles. Note the change in their relative scale with

optical depth. The noise in the standard formulation increases dramatically with

optical depth, while in the diﬀerence formulation it does not

cle counts were selected for equal run times for the two formulations, and

were high enough that the cost of Monte Carlo transport dominated execu-

tion time. The relative performance of the methods, then, is the ratio of the

squares of the standard deviations. There is a small advantage in favor of the

diﬀerence formulation, except at the very right hand side of the slab. Note

that the optical thickness of each zone is only 0.05 optical depths.

In Fig. 1 (b) we show the results for a slab that is 10 mean free paths

thick. Again, we obtain perfect agreement for the two formulations when the

average of the equilibrium material temperature is examined. In order to see

the standard deviation for the two formulations on the same plot, we now

272 A. Sz˝oke, et al.

have to apply a scale factor of 100 for the standard formulation and 1000 for

the diﬀerence formulation. The trend of growth for the standard deviation

for the diﬀerence formulation, as the temperature gradient increases in the

slab from left to right, is becoming apparent.

In Fig. 1 (c) we show the results for a slab that is 100 mean free paths

thick. There are 100 zones in this problem in order to avoid teleportation

error, with every ﬁfth point plotted using a symbol. The other points are

included in the solid lines drawn. The growth in noise for the standard for-

mulation is now becoming visible in the plot. We apply a scale factor of only

10 to see the standard deviation for the standard formulation. The growing

trend of the standard deviation for the diﬀerence formulation, as one tra-

verses the slab from left to right is now quite clear, but the computational

advantage for the diﬀerence formulation is large.

In Fig. 1 (d) we show the results for a slab that is 1000 mean free paths

thick. There are 1000 zones in this problem, a requirement to avoid signiﬁ-

cant teleportation error, with every 50th zone plotted using a symbol. The

disagreement between the two formulations appearing towards the left side

of the slab are statistical ﬂuctuations. The computational advantage of the

diﬀerence formulation is extreme.

At this point, the reader may note that the results for the diﬀerence

formulation have the appearance of smooth curves, regardless of the optical

thickness of the problem. The 1/

√

N noise behavior typical of a Monte Carlo

solution is still present, it is just that the amplitude of the noise is small and

remains small as the optical thickness of the problem increases.

Evidence of the scaling behavior of the computational advantage is shown

in the next ﬁgure. As noted above, the number of Monte Carlo particles for

each optical thickness is set so that the diﬀerence formulation and standard

formulation exhibit the same run time, and that the computer time is dom-

inated by the Monte Carlo, so that the standard deviation is scaling like

one over the square root of the computational work. Under these conditions,

the relative computational eﬃciency of the two methods is the ratio of the

squares of their standard deviations. In Fig. 2 we show the computational

advantage exhibited in the standard deviation of the material temperature,

as a function of position, for the four optical thicknesses. The computational

advantage scales, roughly, as the square of the optical depth of the problem.

Finally, in Fig. 3 we show the penetration of the thermal wave into the

material, 1000 optical depths thick, at an intermediate time. As in Fig. 1 (d),

the material temperature is plotted vs. distance. Reﬁning the time step and

zone size demonstrates that the solution shown is fully converged. This ﬁgure

illustrates the gains that the diﬀerence formulation provides, clearly visible

to the reader.

Accurate and Eﬃcient Radiation Transport in Optically Thick Media 273

Fig. 2. The relative computational advantage of the diﬀerence formulation com-

pared to that of the standard formulation, plotted as a function of position for

various optical depths of the slab. In the 1 OD case, each zone is only 1/20 OD

thick, nevertheless, the diﬀerence formulation is better than the standard one ex-

cept near the surface of the slab on the right hand side. There is a sharp decrease

in computational advantage where the temperature gradient is large. This could

easily have been remedied by spatial importance sampling of the source particles

4.2 Comparison to Analytic Diﬀusion Solution

for a Linearized Problem

In order to check the correctness of our numerical implementation, we have

compared results with the diﬀusion solution appearing in [SO96]. This solu-

tion requires that the material energy take the form

mat

αT

(44)

or, equivalently, a speciﬁc heat of the form

∂E

mat

∂T

= αT

. (45)

The purpose of this form for the material energy is to remove the non-

linearity, T

, that otherwise prevents an analytic solution. The resulting ana-

lytic solution can then be used to check for correct convergent behavior in the

274 A. Sz˝oke, et al.

Fig. 3. Thermal wave (Marshak wave) penetrating a uniform, gray slab of 1000 OD

at an early time, 40 sh. The standard formulation gives a noisy temperature proﬁle

whereas that of the diﬀerence formulation is many orders of magnitude smoother.

The slight diﬀerence in the position of the leading edge is a statistical ﬂuctuation

for the standard formulation

numerical simulation. The behavior of the speciﬁc heat at T = 0, however,

makes things quite fragile unless the numerical code itself is transformed to

handle things in linearized form. To some extent, this defeats the purpose of

checking the original code, much of which had to be modiﬁed to resolve the

issue.

For this problem, α and the cross section were chosen so that the values

in the tables of analytical results appearing in [SO96] could be used directly.

In Fig. 4 we plot the material temperature produced by our Monte Carlo

solution of the standard and diﬀerence formulation, along with data from

Su and Olson’s analytical solution corresponding to a late time, τ = 10. We

obtain good agreement where analytical data are available. This is expected

because for τ = 10 the diﬀusion approximation assumed by Su and Olson is

valid. In Fig. 5 we plot the material temperature for the same problem at a

much earlier time, along with data from Su and Olson’s analytical solution

for a τ = .01. In this case, not surprisingly, there is sharp disagreement