Short characteristic integration of radiative transfer problems: Formal solution in two-dimensional slabs

.I. Quant. Spcctrosc.

Radiat.

Transfer

Vol. 39, No. I, Pp. 61-79, 1988

Printed in Great Britain

0022-4073/8833.00+ 0.00 Pergamon Journals Ltd

SHORT CHARACTERISTIC INTEGRATION OF RADIATIVE TRANSFER PROBLEMS: FORMAL SOLUTION IN TWO-DIMENSIONAL SLABS PAUL KUNASZ~ and LAWRENCE H. AIER Los Alamos National Laboratory, P.O. Box 1663, Los Alamos, NM 87545, U.S.A.

(Received 12 March 1987)

Abstract-A short characteristic method based on parabolic approximation of the source function is developed and applied to the solution of the two-dimensional radiative transfer problem on Cartesian meshes. The method is significantly faster for the evaluation of multidimensional radiation fields than those currently in use. Convergence as a functional of the grid resolution is discussed and linear and parabolic upwind interpolation are compared.

1. INTRODUCTION

The development of new methods for the solution of multidimensional radiative transfer problems is motivated by the existence of both astrophysical and laboratory applications for which the plane parallel or spherical one-dimensional approximation is inadequate or inappropriate. Unfortunately, the techniques that have been already proposed suffer from computational inefficiency. A review of the methods available 10 yrs ago, and a discussion of the physical effects of multidimensional radiative transfer’” can be found in the references. In this series of papers, we will develop methods which are both accurate and efficient. In the present paper, we will address the solution of the two-dimensional transfer equation when the source function is explicitly known. In subsequent papers, we will extend this work to treat noncoherent scattering and non-LTE two-level atoms, The obvious difficulty that arises as one considers multidimensional geometry is that the number of grid points at which the radiation field is to be found is large. For two-dimensional grids, the number of points scales as N*, where N is a measure of the number of points on a side. To achieve the same physical resolution as for one dimension, therefore, requires a factor of N more computation even if the method scaled only as the number of grid points. In fact, for many computational methods, the scaling of the computational effort is much more severe than just the number of grid points. For example, the obvious method, known as the long characteristic or LC method involves the solution of a separate transfer problem for each grid point. Through each grid point and for each direction needed to resolve the angular variation of the radiation field, one constructs an artificial ray along which the physical variables are estimated by interpolation. The transfer problem is then solved for the ray to determine the radiation field at the grid point being considered. Because the time to solve for the radiation on a single ray scales as N, the time to find the full radiation field scales as N3, which is thus even a factor of N worse than the scaling of the number of grid points. The basic nature of the radiative transfer equation is that photons are propagated downstream along a ray. In order to determine the radiation at a mesh point along a ray in a specified direction, it is necessary to know, first, the intensity of radiation entering the cell at the upstream end of the ray segment, or short characteristic, and, then, adjust this quantity by the amount of radiation absorbed and emitted in the cell. The artificial rays constructed in the long characteristic method automatically provide this information, but at the cost of increased computational effort. jCorrespondence

address: 315 Skylark Way, Boulder, CO 80303, U.S.A. 67

PAUL K~NASZ and LAWRENCEH. AUER

68

Alternatively, for a given direction, if we know the radiation at the upstream cell corners, the value of the upwind intensity may be found by interpolation. This is known as the short characteristic (SC) method, and is the basis of the technique we are proposing for the formal solution of the multidimensional transfer problem. Because in the SC technique, for each point in the grid one only solves the transfer equation across the single neighboring cell, the timing for this method scales only as the number of cells, i.e. as N2. The difficulties inherent in LCs led Mihalas et d4 (hereafter MAM) to apply SC techniques for the solution of radiative transfer in two-dimensional slab geometries. At the time of this work, the most powerful techniques for the solution of one-dimensional transfer problems were those based on difference approximations of the second-order differential equation.s7 The multidimensional work presented in MAM was accordingly based on generalizations of the difference schemes used in the previous work. This, as we shall now see, unfortunately sacrifices some of the timing advantages of the SC method. In the MAM approach, illustrated in Fig. 1, at each point a three-point difference equation is constructed. For the point, 0, this will involve the values at the intersection points, M and P. Because physical values are known only on the grid, the values at M and P are found by implicit interpolation on the nine-point stencil. Unfortunately, the resulting equations of the radiation at a point are intrinsically implicit. That is, even if the source functions were a priori known, in order to determine the value at 0, one has to know the values at the eight neighboring grid points not just the upstream points. More specifically, if there are N, grid points in the x-direction, and N, for the z-direction, when we solve the sparse matrix arising from the transfer equation, in the x-direction first, the timing will scale as NINIL, where L is the number of frequency-angle quadrature points. This is a factor of NfNz poorer than the result predicted for the SC approach. The earlier approach of Cannon8 which applies the Feautrier method in x-y geometry, scales as

N:NzL3. The basic reason for using the difference equation approach in MAM was to achieve second order accuracy. This is required because in the limit of great optical depth and high scattering albedo, very large errors can occur unless the formulae being used are exact for parabolic source functions. The difference equations in MAM are sufficiently accurate. In this paper we will show that parabolic interpolation may also be easily used if the source function is specified explicitly and that, therefore, for the solution of the formal transfer problem in two dimensions one may use direct integration of the transfer equation across a cell. Parabolic interpolation yields an SC method which scales only as the number of grid points, NXNZand yet has sufficient accuracy to maintain the correct asymptotic diffusion limit.

Fig. 1. A ray is shown propagating on a two-dimensional Cartesian grid. MO is a short characteristic for this ray. The points a.. .I and 0 form the 13-point stencil used by the parabolic SC method.

69

Formal solution of 2-D radiative transfer problems 2. THE

DIRECT

INTEGRATION

SHORT

CHARACTERISTIC

METHOD

Transfer equation and formal solution on a short characteristic Following

the notation introduced by MAM, we define a two-dimensional planar region (X 2 0), with angle cosines y and ~1on the x- and z-axes, respectively, the transfer equation may be written as

(Z < 0), on which the transfer equation will be solved. For an intensity stream propagating

Wx)@Zlaz)

+ (Y/x)@Z@x) = S - Z

where Z = Z(x, z, y, p, v). Along a characteristic is

path for a given frequency, the transfer equation

az/as = x(S - Z), where s is the characteristic

coordinate

(1)

(2)

and is measured along the direction of propagation ds = dx/y = dz/p;

also x(x, z, v) is the monochromatic opacity. In Fig. 1, point M is the upwind cell boundary on the short characteristic MO. The intensity at target-point 0 is the sum of its value at M, attenuated by the optical length Ar, on segment MO (called At, on segment OP), and an integral for the cell contribution. This sum is the following formal solution at 0: 44 Zo = ZMexp( - AZ,) +

dtS(t)exp[-(Ar, - t)]. (3) s0 The optical depth step AZ, is the path integral of x(x, z) from M to 0. Since x and S are known everywhere on the spatial mesh, their values at M and P are easily found by interpolation. Then, using values at M, 0 and P, interpolation formulas can be constructed for the evaluation of the integrals for the MO and OP intervals. The advantage of Eq. (3) is that it is analytically correct. Its accuracy is determined by the numerical approximations used in the evaluation. In fact, it is often critical that parabolic (or better) interpolation be used in the computation of the integral. Linear interpolation for S(r) yields a biased value of Z and does not recover the correct asymptotic behavior in optically thick media for nonlinear source functions. This problem is particularly pronounced when the formal solution is being used, as in the next paper of this series, for media in which the source function is dominated by scattering. Fortunately knowledge of the source functions at M, 0, and P permits the construction of a quadratic approximation for s(t).

Discretization

The spatial mesh consists of an x-grid {xi, i = 1, ZV,},and a z-grid (zj, j = 1, N,}; and we adopt the notation X = x(N,), and Z = z(N,). The distribution of space points may be uniform, or may have geometrically increasing intervals, with or without mesh symmetry about the midplanes. To specify the mean intensity, we must, in general, consider the full angular range - 1 < y 6 1, and - 1 < 1( G 1. If the slab properties and incident intensity are symmetric about the midplane, we have Z(x,z,y,cL,v)=Z(X--x,z,

-y,/J,v),

and we may limit our range to 0 < y < 1, and - 1 G ZJ< 1, while continuing to treat the full space mesh. For all models treated in Sec. 3, we have used discrete angle pairs (y, p) from set B of Carlson. For the examples presented in this paper, the line profile function 4, is a constant-width Doppler profile, and the quadrature sum for J is done by trapezoidal integration on a set of frequency points with 0.5 Doppler-units intervals. If we let k be the index for summation over all combinations of angle and frequency points and let W, be the combined quadrature weight, then the scattering integral

J(x,z)=&

s 4%

do

co dv4 (x, z, 0, v )Z(x, z, a, v) s -cc

PAUL KUNASZ and LAWRENCEH. AVER

70

may be approximated

by (4) k=l

In these equations, we permit both frequency and angle dependence in the profile function, which is required for treating media with different mass motion. Evaluation of the formal solution The intensity value Ihl is found by interpolation on the surrounding upwind intensities. The solution is then advanced by using Eq. (3) to compute the contribution as a ray crosses the cell to the downward vertices. In this way the mesh points are successively filled in with intensity values generated from previously computed upwind intensity values. For rectangular orthogonal grids, the upwind values will always be known as long as one evaluates the nodes in an order which works away from the boundaries, where the incident radiation field has been specified. Whether one computes the values along the x or z axis is immaterial as long as the computational order along the axis moves away from the boundary condition. As we will now show parabolic interpolation for the upwind intensity is easily done on orthogonal grids even if the spacing irregular; however, if the grid is distorted, only linear interpolation is convenient. Although, as we will demonstrate in the next section, parabolic interpolation is usually superior on rectangular grids, we will also present results for linear upwind interpolation in order to justify its use for distorted multidimensional grids and in cases where the source function gradients are so large that parabolic interpolation leads to nonphysical negative values. Parabolic interpolation for ZMcan be usually be done with the three points which are members of the rectangular nine-point stencil centered at target point 0. However, parabolic interpolation may require the use of two points on the rectangular stencil, and one external point, since the only usable members of the nine-point stencil are those on which the intensity has been found, i.e. those on the current and upwind levels. As long as we follow the prescription of always evaluating the intensities in the order which goes away from the boundary conditions, three upwind points will always be available, except in a row or column immediately adjacent to a boundary where incident radiation is being specified. When all possibilities and ray directions are considered, the interpolation stencil for explicit calculation with parabolic fits is seen to be the 13-point stencil shown in Fig. 1. By contrast, in the implicit procedure of MAM, all interpolations were done on the nine-point stencil because all nine points were always available. Which meshlines are cut by ray MPO, the position of points M and P, and the pathlengths As, and Asp, are determined by the signs and relative sizes of y and p, and the relative lengths Ax and AZ, along MO and OP. To systematize the notation we observe the following conventions: we shall use the upper (lower) signs on subscripts involving the index i if y > 0 (y < 0), and the upper (lower) sign on subscripts involving the index j if p > 0 (p < 0). We set aM =

IY IAZ~*I,~/ICL I&I/Z

and bM = l/aM where Azj, 1/2= Izj+I - zj I and A+ i/z = IXiki- XiIf Then, if aM -C 1, point M lies on a horizontal meshline and

AsM=AZ~+I~/IPI; whereas, if bH < 1, M is on a vertical meshline and

Formal solution

For the physical properties formally as

of 2-D radiative

of the slab, the parabolic i+l V,

=

transfer problems

interpolation

71

formula can be expressed

j+l

C

i’zi__l

mi.,j.Vi.,j.

C

(5)

j’cj-1

where V denotes any physical variable, and the indices (i,j) refer to point 0. At most, only three of the m-values will be non-zero, and mu is always zero. If 1, is to be obtained by linear interpolation, we have, for uM < 1, Iivi= (l -

dzi-

I,jlrI +

I

%l;:,j+

@-4

and, for &,,,< 1, 1, = (1 - b,)Zi,

I,j_

1 +

b,Zi;

I,j.

WI

For parabolic interpolation, the standard expressions are used. For the cell distribution, we must evaluate the integral 4 C, = exp( - AT,) S(t)exp(t) dt. s0 In general,

If Ci,j is evaluated by linear interpolation

of S(t) on the segment MO, then

ye = Uo(Ar,) - U,(Ar,)IAr,,

(74

‘J’o= U,(k,)lAr~,

(7b)

Yp=O.

(7c)

and

If parabolic interpolation

is used on the segment MOP, then

‘y, = Uo(Ar,) + [U&W

‘yo= Kb, + bWd&J

- (Ar, +

~A~,)U,(AZ~)II[AZ~(AZ~ + Az,)],

- W&,VI)1/&&,

(84 VW

and

where U,(x) = 1 - exp(-x)

Pa)

U,(x) = x - 1 + exp( -x)

= x - Uo(x),

W-4

and U*(x) = x2 - 2x + 2 - 2 exp( -x)

= x2 - 2U,(x).

PC)

3. RESULTS

The upwind interpolation intrinsic to the SC technique introduces angular dispersion into the solution for the propagation of radiation; thus, the method can never be exact, no matter how Eq. (3) is integrated. In this section, we will address the global accuracy of the SC method. The two questions of particular importance are (1) how does the size of the dispersion vary as a function of the grid resolution and (2) what are the relative merits of parabolic and linear interpolation of the upwind intensities. To these ends, we concentrate on the magnitude and distribution of error, rather than on the distribution of the radiation field itself. The reader interested in the phenomenology of line transfer in a 2-D slab is referred to MAM.

PAUL KUNASZ and LAWRENCEH. ALER

72

There are several free parameters which control the spatial meshes and the angle- and frequency-quadratures. We define the meshes by their starting values (x = z = 0), their final values (X = X > 0, z = 2 < O), and the sizes of the first steps off the boundaries. From these values and the constraint of symmetry about the midplane, a constant step-size ratio, r, or r,, from boundary to midplane, is found. We also consider uniform meshes, with constant intervals (rX= r, = 1). For these, the asymptotic error term for each interpolation process is better by one order. When N, = N,, as in many models treated here, the notation n is used for the common value. The profile function is a constant-width Gaussian. Frequency integration is trapezoidal with intervals of 0.5 Doppler widths. The angle quadrature is from set B of Carlson,g with 1 or 3 angles per octant, a number which we will denote by N,. Error is evaluated on the mesh by comparison with an exact analytic solution, or a highly accurate numerical solution found by integration on the long rays constructed by extending to the slab boundaries the rays emanating from each mesh point in the direction of the angle quadrature. If Q represents Z or J the error E is given by E = 200 /, where Q, is the exact or accurate value. E is percent relative error. All models have left-right symmetry and no background opacity or emissivity. Unless otherwise specified, X and 2 take the values 1, and -2, respectively. Other aspect ratios were treated, with results very similar to those below. Cell integration is done by parabolic interpolation of the source function in f-space, while optical depth t is obtained from a parabolic fit of the opacity x(x, z). By contrast, in MAM, in which solar and stellar applications are of greater interest, linear interpolation was done on log[x(x, z)]. Our models fall into two main classes. In case (I), a unidirectional narrow beam, is incident on the slab and propagates in vacuum across the space mesh. This is the well known searchlight beam problem, and it isolates the effect of 2-D upwind interpolation on streaming radiation. The interpolation error is manifested as diffusion of intensity out of the beam and into the surrounding space. We study the convergence of the computed radiation field, i.e. the reduction of the numerical error with mesh refinement, and compare the errors for linear and parabolic upwind interpolation. In case (2), a diffuse radiatic,’ field is generated by isotropic emission from smoothly distributed sources in the interior of the slab. We divide the second class into models that are optically thin and thick, and, again, concentrate on convergence of the radiation field, and comparison of the two interpolation schemes. Searchlight

beam

Convergence. The upwind interpolation involved in the SC method causes the intensity of an upwind point to be distributed over all its downwind points. The result, as noted, is the dispersion of any unresolved feature. It has been argued that because increasing the grid resolution corresponds to an increased number of upwind interpolations. As we will now show, this is fortunately not the case. Dispersion of a feature is most clearly demonstrated by considering the propagation of a beam across a vacuum. The ray has sharp edges and the solution is analytic: a beam with no attenuation across the grid. Further, because the integral in Eq. (3) is analytically zero in each cell, any error in the solution is attributable to the effect of upwind interpolation. In Fig. 2, we show the results of the SC method for the propagation of such a beam across a grid which is 1 x 1.5 in size using linear interpolation. Figure 2(a) shows a half-tone plot for the case where N, = N, = 16 with the ray moving at 45” with respect to the vertical. It should be noted that this means Ax # AZ, if the steps were equal, for 45” there would be no dispersion in this particular case. The beam is not sharp, partly because on any finite grid a sharp edge reduces to a step function which coupled with the interpolation used for the plotting leads to a fuzzy edge, and partly because the effect of the dispersion which makes the beam wider as it moves across the grid. Both these effects are dramatically reduced in Fig. 2(b), where we have used the SC method but with N, = N, = 64. We will discuss the rate of this convergence numerically in the next section. Comparison of linear and parabolic upwind interpolation. We compare the effects of linear and parabolic interpolation for a hard-edged beam in Fig. 3 and Table 1. The models assumed here

Formal solution of 2-D radiative transfer problems

Fig. 2. Gray-scale plots for a searchlight beam through a vacuum computed by using 16 x 16 and 64 x 64 mesh points.

have dimensions X = 1 and Z = - 1.5, with the searchlight beam entering through the bottom boundary from x = 0 to 0.2, at an angle of 33” to the z-axis, i.e. such that the upper edge of the beam traverses the slab from lower left to upper right. The mesh cells are constant and square, which implies the relationship (NZ - l)/(NX - 1) = IZ/Xl. NX and NZ are defined for three meshes as: 11 x 15, 32 x 49, and 53 x 79. This implies the linear cell dimensions: 0.1, 0.05, 0.031, and 0.019, respectively. In Fig. 3, we compare three solutions: linear, parabolic, and exact, for each of the four meshes. A trough of negative intensity values is clearly seen at the edge of the beam in the parabolic case, and general reproduction of the exact solution seems little better than in the linear case. As mesh resolution increases, the trough narrows on the scale of the model, but, as expected for the case of a hard-edged beam, becomes no less deep. These visual impressions are confirmed in Table 1, in which maximum absolute, average, and root-mean-square errors on the mesh are tabulated in percent for linear and parabolic interpolation. For the parabolic case, two

74

PAUL KIJNASZ and LAWRENCEH. AIJER

Fig. 3. Elevation plots of a hard-edged searchlight beam, computed with linear (top) and parabolic (middle) upwind interpolation; the exact comparison beam (bottom) is also shown, for uniform meshes of increasing resolution. The hard edges of the beam cannot be resolved by refinement of the mesh, and they induce negative overshoot for the parabolic case.

Table 1. The maximum, average, and r.m.s. percent errors in the intensity are shown for the searchlight beam in an empty slab with X= 1 and Z= -1.5. CELL SIZE

!

0.1

1

0.05

/

0.03

1

0.01

Em

67

67

51

50

Ea Er

292

1;

166

1:

Em EF

74 2:

74 1:

65 144

64 1:

Linear

1 Parabolic

The mesh cells are squares, and their sizes are given in the top row. For parabolic interpolation, the largest negative intensity Ineg and the percentage of intensities are also tabulated. All results except Znegare rounded to the nearest integer value.

15

Formal solution of 2-D radiative transfer problems Table 2. The same model is used as in Table 1 but with a logarithmic mesh and boundary steps set equal to 0.01. The step-size ratios are as given. MESH SIZE

12x15

21x31

33x49

53x79

1.14 1.09

1.05 1.03

2.32

1.33

I.7

1 .a0

1.20

Em Ea Er

94

73

64

55

b02

283

169

155

89

55

67

A?

179

70 5 16 -0.14 15

-0.13 19

IX

1

Linear

beg % neg

-0.16 9

-0.08 11

Parabolic 132

additional quantities are tabulated: most negative intensity, and the fraction, in percent, of mesh points on which the intensity is negative. The numbers in Table 1 speak well for the use of linear upwind interpolation for radiative transfer in media with radiative discontinuities. Similar conclusions can probably be drawn for material discontinuities treated by linear interpolation in the cell integration. Though parabolic interpolation is almost as economical as linear, it is not tractable on 2-D and 3-D Lagrangian meshes, thus the importance of the comparative error which we find here for orthogonal meshes. In Table 2, results are presented for the same model, with steps increasing geometrically, from Ar x = At, = 0.01. The column for cell size, now variable across the mesh, is replaced by a column for mesh size. We see that the logarithmic mesh yields larger errors than the uniform mesh, for coarse meshes, but similar errors for the fine ones. Parabolic interpolation improves somewhat over linear in the transition from uniform to logarithmic meshes. Much of the error presented in Tables 1 and 2 and in Figs. 2 and 3 is due to the hard edges of the beam, i.e. its discontinuous spatial distribution. In order to study resolvable streaming radiation, we assume a beam with a Gaussian spatial distribution. The exact intensity is given by Z(x, z) = exp[ - (D/W)2], where D is the perpendicular distance from the beam (even if the perpendicular line segment lies partly outside the slab), and W is 0.20. The center of the Gaussian beam enters the slab at the midpoint of the bottom boundary, and is directed toward the upper right corner. For the model assumed in Table 3, X = 1,Z = - 1S, and all meshes are square except 21 x 36, which is slightly rectangular. One should note that the errors shown in Table 3 decrease as would be expected analytically as a function of the grid resolution. For linear interpolation the decrease is linearly proportional to the number of grid points (i.e. proportional to ~/AZ). For parabolic interpolation the convergence is quadratic. This leads to the significant advantage of parabolic interpolation on fine meshes. Table 3. Percent errors as in Table 1 but the source of radiation is an incident Gaussian beam with half-width = 0.2 MESH SIZE

11 x16

21x36

31 x46

51x76

Em :;

37.1 10.6 4.6

25.6 6.8 2.5

18.4 4.7 1.9

12.3 3.0 1.2

Em :;

23.1 3.0 6.8

7.7 0.7 2.0

3.4 0.4 0.9

1.2 0.1 0.3

‘min % Neg

-0.08 18.0

-0.01 17.0

-0.00 13.0

-0.00 12.0

Linear

Parabolic

PAUL KUNASZand LAWRENCEH. AUER

76

Table 4. The same model is used as in Table 3, with logarithmic meshes having AX = AZ = 0.01 and step-size ratios rx and r, as given. MESH

11x16

SIZE

21x36

2.32 1.80 7.3 20.1 48.0 6.6 16.7 -0.04 8.0

1.05 1.04

1.16

1.33 1.16

51 .o

lmin % Neg

31 x46

1.10

35.4 3.9 12.7

23.2 2.9 7.9

27.6 2.0 6.8 -0.103 13.0

11.1 0.90 2.7

13.7 1.59 4.24

Linear I

2.5 0.23 0.64

-0.03 13.0

Parabolic

-0.0008 12.0

For the Gaussian beams the results of parabolic interpolation clearly fit the exact solution better than those of linear interpolation, except for the negativity near high curvature, and the ripple of positive values next to the negative trough. These effects are significant only on the coarsest mesh. Despite the ripples due to parabolic overshoot, the quantitative comparison in Table 3 favors the parabolic scheme by a factor of 10 in all three measures of error, for the finest mesh. In Table 4, results are given for the model of Table 3, for a logarithmic mesh, with Ax = AZ = 0.01. The parabolic upwind scheme, though still more accurate, does not have quite as much advantage as in the case of the uniform mesh. In both cases, the advantage increases sharply with mesh refinement. D@ise

sources

Optically-thin slabs. In Table 5, the error in 1 is presented for a homogeneous slab with x = S = 1. The meshes are uniform, and N, = 1. Because of the fact that J is a composite of the monochromatic frequency components J,, and because the various frequency components resolve the grid differently, we present results for the far wing of the line, where the opacity is nearly zero in Table 6 as well. The far wing contributes almost nothing to 7, but illustrates the optically thin

Table 5. Maximum, average, and r.m.s. errors in J as functions of mesh size are shown for a homogeneous slab with X = 1.0, and 2 = -2.0, x(x, z) = S(x, z) = I, and N, = 1. The meshes are uniform. N

Em Ea Er

( 11

(

9.65 2.31 3.16

19

7.01 1.39 1.92

(

29

( 39

5.51 0.93 1.33

1

4.68 0.69 1.02

59

1

79

3.19 0.35 0.56

::CZ 0.71

( 151 / 201

2.26 0.18 0.32

1.95 0.14 0.26

Linear

Table 6. Maximum, average, and r.m.s. errors in monochromatic mean intensity are shown in the far wing of the line, as functions of mesh size, for the slab and meshes of Table 5. N

Ea Er Em

1

11

1

19

1

29

1

39

1

59

1

79

/

151

2.44 3.72 13.0

1.42 2.26 9.66

0.93 1.56 7.70

0.69 1.20 6.59

0.46 0.85 5.31

0.34 0.66 4.57

&Xl: Linear 3.27 1

7.55 1.41 2.06

4.83 0.72 1.08

3.49 0.43 0.68

2.80 0.30 0.48

X 0.31

1.68 0.12 0.23

1.17 0.11 0.14

Parabolic

Formal solution of 2-D radiative transfer problems Table 7. Maximum errors are shown in f for the model of Table 5, for a logarithmic space mesh, with As, = Ar, = 0.03; r, and rr are as given. N

7

11

19

29

39

rx rz

3.49 5.21

1.63 2.05

1.14 1.31

1.03 1.12

1.00 1.06

Em Em

8.01 7.65

5.55 3.50

3.56 2.54

3.82 2.83

3.68 2.77

Linear Parabolic

limiting case. Several conclusions can be drawn from the data. (i) Convergence is present over a broad range of mesh sizes; (ii) the error for parabolic fits is roughly a factor of two smaller than that for linear fit; (iii) rates of convergence for line center and line wing are roughly the same, though the degree of error in the wing is greater, and (iv) unlike the searchlight, the asymptotic rates of convergence for J are Er z O(dxO.‘*) and

EP w O(dx0.66),

where dx is the step-size of the mesh. The improvement of EP over E, suggests that, for some problems, higher order interpolation might produce better results than obtainable with parabolic interpolation. Whether the increased complexity of such fits, including the logic of reaching farther off stencil for target points, would be worth the savings in numbers of mesh points needed to achieve a desirable accuracy, has yet to be determined. For the present we will rely on the linear scaling of the method to economically reduce error. In Table 7, maximum errors are presented for the homogeneous model of Tables 5 and 6 with logarithmic mesh spacings, with AZ = 0.03. N increases from 7 to a value large enough that step-size ratios approach unity. Parabolic interpolation is more accurate but not by a large factor. Optically-thick slabs. The models of Table 8 have x = lo* and S = 1, with N, = 1 and Ar x = Ar, = 0.03. For the constant source function, errors in J in the interior are typically small due to the effect of saturation. That is, for optically thick mesh cells, the upwind intensity value, and whatever interpolation errors it carries, are unimportant because the upwind intensity is highly attenuated at the cell center (point 0 in Fig. 1). The frequency components of the radiation field in the wing of the thick line behave exactly as the thin monochromatic components examined above, because they are governed by the same equation. However, they are heavily deweighted in the scattering integral by convolution with the profile function. All measures of error decrease monotonically with mesh refinement. The two interpolation schemes yield roughly the same value for N = 21, with the parabolic scheme providing smaller average and r.m.s. errors for finer meshes. Except for the coarsest mesh, maximum errors for the two schemes are very close. All models up to this point have had zero error in cell integration due to the functional simplicity of x and S; the entire error has been due to upwind interpolation along mesh lines. In Table 9, this constraint is relaxed. Results are given for a model with x(x, z) = lo* and Table 8. Maximum, average, and r.m.s. errors in J as functions of mesh size are shown for a homoeeneous slab with Y = lo*. S = 1. N. = 1. and AT = Ar, = 0.03. The step-size ratios,‘;, and’ r,, are as&given. N

79

rx b

1.15 1.18

Em Ea Er

0.53 0.07 Linear 0.13 1

E; Er

1.29 0.17 0.30

I

0.96 0.09 0.16

I

0.77 0.05 0.10

78

PAUL KLJNASZand LAWRENCE H. AUER

Table 9. Same model as for Table 8, with S(x, z) = sin(nx/X)sin(nz/Z). 15

21

29

13.0 6.43

4.09 6.29

3.14 5.56 3.65

3.16 1.66 1.98

N

9.50

4.46

39

49

3.45 2.21

2.51 1.45

1.11 1.96

1.25

0.90 1.73 1.03

0.89 1.95 1.12

0.60 1.48 0.79

0.49 1.32 0.66

0.43 1.27 0.61

2.42

1.60

59

1 Linear

Parabolic

S(x, z) = sin(xx/X)sin(xz//Z). Other parameters are as in Table 8. The fact that S(x, z) is a transcendental function implies that cell integration with linear or parabolic interpolation of S(x, z) is not exact. Comparing to the results for S = 1, we find larger errors in the interior of the slab, due to the fact that the solution do not saturate to a constant value. The parabolic scheme for upwind interpolation remains superior, especially for coarse meshes. As with S = 1, the method converges monotonically. In Table 10, results are presented for x(x, z) = exp(5z/Z) and S(x, z) = sin(nx/X)sin(xz/Z), other parameters being the same as for Tables 8 and 9. Larger errors are found than for constant opacity. The large maximum error for linear upwind interpolation is greatly reduced by parabolic interpolation as mesh resolution increases. At this point it becomes interesting to vary the interpolation scheme for x(x, z), and for S(x, z) in the cell-contribution integral, Eq. (3). In Table 10, both fits are parabolic. In Table 11, both are linear for the first four meshes, and only the interpolation of x(x, z) is linear for the last three. Surprisingly, the errors are not sensitive to the interpolation scheme for x(x, z). Even more surprising is the fact that when both x and S are interpolated linearly, the parabolic upwind scheme yields the best results of all. This must be due to lack of “proper” mesh resolution for the sinusoidal source function in the central region of the slab, where mesh steps are large. We define “proper” mesh resolution here, as that which allows the higher-order scheme to prevail, knowing that for a sufficiently tight mesh it will prevail. Timing For a spatial mesh of 41 x 41 points, the processing time required on a Cray la computer for

Table 10. The same model is used as for Table 8, with S(x, z) = sin(rrx/X)sin @z/Z) and x(x, z) = exp(5z/Z). 15

21

29

54.6 a.87 13.9

31.6 4.70 6.84

F% 4.37

42.2 13.3 17.1

18.5 5.57 6.92

9.10 2.85 3.52

39

4.12 1.55 1.84

49 13.5 1.56 2.42

Linear

2.79 1.03 1.23

Parabolic

i Table 11. The same model is used as for Table 10. N

15

15.0 Em 36.9 18.6

EF Em Ea Er

20.4 4.84 7.34

=I= 21

29

39

15

21

29

21.7 6.85 9.09

19.5 3.07 5.26

15.4 2.55 3.46

57.2 0.87 14.3

32.6 4.71 6.95

2.70 Linear 4.41 1 23.6

11.3

7.53 1.11 2.01

4.81 0.64 1.07

47.3 12.7 16.8

19.9 5.57 7.06

9.53 2.89 3.62

%

Parabolic

[For the first four meshes, linear interpolation is used for both x(x, z) and S(x, z). For the last three, a linear fit is used for x(x, z) and a parabolic fit for S(x, z).

Formal solution of 2-D radiative transfer problems

19

each frequency-angle point is 0.01 s. Since the scaling is direct with the number of angle-frequency points and with the number of spatial mesh points, we may use the timing formula: T (s) = 0.01N3$&N,/41*. The reader may at first suspect that the evaluation of the exponential in Eq. (9a) will dominate the processing. For a vector computer this is far from true, since one can fill in a string of exponential values as a vectorized process. For the Cray la, the dominant step in the formal solution is the sequential (non-vector) evaluation of the specific intensity over the mesh points of the slab. The process can be vectorized over the frequency-angle set, at the cost of greatly increased storage, but this move is advisable only if the size of the frequency-angle set is sufficiently large relative to the crossover point for vector strings. The crossover point is the point at which saving due to vectorization outweighs the setup time for the loop. For a Cray la, this size is about 10 elements. 4. CONCLUSIONS

The fundamental advantage of the short characteristic method over the long characteristic method is speed. Its only known disadvantage is the inability to track collimated radiation fields without producing lateral dispersion within the same direction of flow. The strength of this source of numerical error is examined by computation of the searchlight beam. We have found that the short characteristic method converges the problem of the hard-edged searchlight beam propagating through a vacuum, with the caveat that if parabolic upwind interpolation is used, negative intensities will occur on the mesh points closest to the beam. Maximum, average, and r.m.s. errors over the space mesh depend only weakly on the choice of upwind interpolation scheme, or whether the mesh has constant or geometrically increasing intervals. For searchlight beams with resolvable edges, errors averaged over the mesh decrease as the square of the step-size for parabolic interpolation, and as a linear function, for linear interpolation, for uniform meshes. For logarithmic meshes, the same relationship holds more loosely. For both types of meshes, the upwind parabolic scheme prevails as the mesh is refined, but more so for the uniform mesh. For optically thin, homogeneous slabs with internally distributed sources, convergence for uniform meshes is somewhat slower than first order for linear interpolation, and faster for parabolic. If the physical gradients are not well sampled on the mesh, parabolic interpolation can give rise to negative intensities, and even negative values of J. However, in a much wider selection of models than those presented here, we have observed negative J-values only in connection with discontinuous distributions of incident intensity. If it is important to avoid negative overshoot in the vicinity of a discontinuity, either the mesh must be arranged to handle the problem, the discontinuous source component must be treated separately, or linear upwind interpolation (or, in the case of material discontinuities, linear fits for cell integration) must be used. These three approaches are all superior to the oft-confessed practice of fixing up the negative values by zeroing them after they are found. The price paid for the positivity of linear upwind interpolation is, of course, an asymptotic truncation error one order lower than that of parabolic interpolation. REFERENCES 1. H. P. Jones, Astrophys. J. 185, 183 (1973). 2. H. P. Jones and A. Skumanich, Astrophys. J. 185, 183 (1973). 3. A. L. Crosbie and T. L. Linsenbardt, JQSRT 19, 257 (1978). 4. D. Mihalas, L. H. Auer, and B. W. Mihalas, Asrrophys. J. 220, 1001 (1978). 5. P. Feautrier, C. r. Acad. Sci. Paris 258, 3184 (1962). 6. L. H. Auer, JQSRT 16, 931 (1976). 7. G. B. Rybicki, JQSRT 11, 589 (1971). 8. C. J. Cannon, Asrrophys. J. 161, 255 (1970). 9. B. G. Carlson, in Methods of Computational Physics (Edited by B. Alder, S. Fembach, and M. Rottenberg). Academic Press, New York (1963).