A Radiation Transfer Solver for Athena using Short Characteristics

0.81 MB·English

by  Shane W. Davis

Draftversion January12,2012
PreprinttypesetusingLATEXstyleemulateapjv.4/12/04
A RADIATION TRANSFER SOLVER FOR ATHENA USING SHORT CHARACTERISTICS
Shane W. Davis1, James M. Stone2, and Yan-Fei Jiang2
Draft version January 12, 2012
ABSTRACT
We describe the implementation of a module for the Athena magnetohydrodynamics (MHD) code
which solves the time-independent, multi-frequency radiative transfer (RT) equation on multidimen-
2 sional Cartesian simulation domains, including scattering and non-LTE effects. The module is based
1 on well-known and well-tested algorithms developed for modeling stellar atmospheres, including the
0 methodofshortcharacteristicstosolvetheRTequation,acceleratedLambdaiterationtohandlescat-
2 teringandnon-LTEeffects, andparallelizationvia domaindecomposition. The module servesseveral

n purposes: it can be used to generate spectra and images, to compute a variable Eddington tensor
a (VET) for full radiation MHD simulations, and to calculate the heating and cooling source terms
J in the MHD equations in flows where radiation pressure is small compared with gas pressure. For

1 the latter case, the module is combined with the standard MHD integrators using operator-splitting:
1 we describe this approach in detail, including a new constraint on the time step for stability due
to radiation diffusion modes. Implementation of the VET method for radiation pressure dominated

] flows is described in a companion paper. We present results from a suite of test problems for both
M
the RT solver itself, and for dynamical problems that include radiative heating and cooling. These
I tests demonstrate that the radiative transfer solution is accurate,and confirmthat the operator split
h. method is stable, convergent,and efficient for problems of interest. We demonstrate there is no need
p to adopt ad-hoc assumptions of questionable accuracy to solve RT problems in concert with MHD:
- the computational cost for our general-purpose module for simple (e.g. LTE grey) problems can be
o
comparable to or less than a single timestep of Athena’s MHD integrators, and only few times more
r
t expensive than that for more general (non-LTE) problems.
s
Subject headings: (magnetohydrodynamics:) MHD methods: numerical radiative transfer
a
[

1 1. INTRODUCTION from optically thick to optically thin regimes, such
v as flux-limited diffusion (e.g. Levermore & Pomraning
Radiation is of fundamental importance for the ther-
2 1981, hereafter FLD). This includes applications such
modynamics of most astrophysical systems. It can be
2 as accretion flows, star formation, neutrino transport in
thedominantsourceofheatingandcoolingofastrophys-
2
supernovae,stellaratmospheresandwinds,cosmological
2 ical plasmas. Even in those systems where it plays a
reionization, and many others. Indeed, there is a large
. minor role in energy transport,it is the dominant mech-
1 and growing list of astrophysical MHD codes that
anism through which we perceive and explore the uni-
0 utilize FLD or a similar prescribed closure relation
verse. Nevertheless, it has often proven difficult to di-
2 (includinge.g.Turner & Stone2001;Bruenn et al.2006;
rectlymodeltheeffectsofradiationaccuratelyinmodern
1 Hayes et al. 2006; Gonz´alez et al. 2007; Krumholz et al.
multidimensional astrophysical(magneto)hydrodynamic
v: (MHD) codes due to both computational expense and 2007a; Gittings et al. 2008; Swesty & Myra 2009;
i conceptual complexity. Commer¸con et al. 2011; van der Holst et al. 2011;
X
Zhang et al. 2011).
Most approaches to adding radiative transfer to dy-
r Numericalmethodsfordirectlysolvingradiativetrans-
namical simulations are based on adopting restrictive
a
fer (RT) have been implemented (e.g. Stone et al.
assumptions or approximations. For example, often
1992; van Noort et al. 2002; Hayes & Norman 2003;
the flow is assumed to be optically thin to radiation
Hubeny & Burrows2007),buttheirapplicationtoastro-
everywhere and for all time, or the radiation field
physicalproblemshasbeensomewhatlimited,especially
is assumed to originate in a small number of point
in full 3D. A notable exception is the progress made in
sources,withthe diffuse emissionfromscatteredorrera-
simulating the atmospheres of the Sun and other cool
diatedphotonsignored(suchasinCosmologicalreioniza-
stars. In the solar physics community, multidimensional
tion problems e.g.Abel & Wandelt 2002; Mellema et al.
MHD simulations of convection with realistic RT have
2006; Rijkhorst et al. 2006; Whalen & Norman 2006;
been performed for decades with increasing sophistica-
Reynolds et al. 2009; Finlator et al. 2009).
tion. (see e.g. Nordlund 1982; Stein & Nordlund 1998;
For problems in which the diffuse emission cannot
Vo¨gler et al. 2005; Heinemann et al. 2007; Hayek et al.
be ignored, the dynamics of the radiation field is often
2010)
treatedbysolvingthe radiationmomentequationsusing
Encouraged by recent work modeling the departure
ad hoc closure prescriptions to handle the transition
of the radiation field from local thermodynamic equi-
1CanadianInstituteforTheoreticalAstrophysics. Toronto,ON librium (LTE) due to the presence of electron scat-
M5S3H4,Canada tering in three-dimensional MHD simulations (see e.g.
2 Department of Astrophysical Sciences, Princeton University, Hayek et al. 2010), we have implemented a general-
Princeton,NJ08544,USA
purpose RT solver in Athena (Stone et al. 2008), based
2
on the methods widely used in the stellar atmospheres In the above, ρ is the gas density, p, v is the fluid
community. Athena is a general purpose astrophysical velocity, and B is the magnetic field. The total stress
MHD code, which is being actively developed and al- tensor T is defined as
ready includes several modules for handling a variety of T=(p+B2/2)I BTB, (5)
physicalprocesses. Effectively,wehavecombinedAthena −
and E is the total (fluid) energy
withamodernstellaratmospherescode. Infact,Athena
already has a RT module that computes the effects of p 1 B2
E = + ρv2+ , (6)
ionization radiation from a single point source on the γ 1 2 2
surrounding gas (Krumholz et al. 2007b). However, this where p is the gas pres−sure and I is identity matrix.
moduleisnotwell-suitedformodelingtheradiationfrom
Thesourcetermonthe righthandsideofequation(4)
diffuse continuum emission.
is the net gain or loss of energy due to radiative heating
The addition of a RT solver to Athena enables three
and cooling and is given (for a static medium) by
goals: (1) it canbe used as a diagnostic toolto compute
self-consistentlyspectraandimagesfromtime-dependent Q =4π ∞χtot(J S )dν. (7)
MHD simulations for direct comparison to astronomical rad ν ν − ν
Z0
observations; (2) it allows us to compute a variable Ed- This is an integral over frequency ν of the difference be-
dington tensor (VET) for the integration of the coupled tween mean intensity J and the total source function
ν
MHD and radiation moment equations (Sekora & Stone S , weighted by the total opacity3 We do not attempt
ν
2010; Jiang et al., submitted to ApJS, hereafter JSD12) to add the corresponding radiation source term to the
forfullradiationMHDsimulationsinregimeswhereboth momentum equation. This limits us to applications in
energy and momentum transport by photons is impor- which radiation pressure is at most a modest fraction of
tant;and(3)itallowsustocomputetheradiationsource gaspressure. AnintegratorforthecoupledMHDandra-
terms in the energy equations and directly couple them diation moment equations based on the one-dimensional
totheMHDintegratortocomputethedynamicsofflows algorithms discussed in Sekora & Stone (2010) has been
where radiation pressure can be ignored. implementedinAthenaandextendedtomultidimensions
This paper focuses on describing our implementation by JSD12. These more advanced techniques are needed
ofmethodstosolvetheRTequation,andthecouplingof to handle the stiff source terms and modified dynamics
the solver with the MHD integrator to compute the ra- in radiation pressure dominated flows.
diationsource term in the energy equation. The compu- In order to compute the energy source term due to
tation of the VET and solution of the radiationmoment radiation,the MHD equationsmustbe supplementedby
equations is described in JSD12. The plan of this work the time-independent equation for RT
is as follows: In Section 2 we summarize the equations
nˆ I =χtot(S I ), (8)
that are solved. In section 3 we describe the detailed ·∇ ν ν ν − ν
implementation of our solver and the iterative methods where Iν is the specific intensity for an angle defined
usedmodeldeviationsfromLTEandhandlecertain(e.g. by the unit vector nˆ. In this work, we consider opaci-
periodic) boundary conditions. In section 4 we describe ties due to scattering χsνc and true absorptionχaνbs, with
how we compute the radiation source terms in the en- χtot = χabs +χsc. It is convenient to define the pho-
ν ν ν
ergy equation and incorporate them into the MHD in- ton destruction probability ǫ = χabs/χtot. The source
ν ν ν
tegration. In Section 5 we present the results of several function is then given by
test problems not only to assess the accuracy of the RT
S =ǫ B +(1 ǫ )J , (9)
ν ν ν ν ν
solver,but alsoto evaluatethe performanceofthe MHD −
whereB isthermalsourcefunction. Themeanintensity
integrator when the energy source terms are included. ν
J is the “zeroth” moment, or average, of I over solid
We summarize our results in Section 6. ν ν
angle
2. MHDEQUATIONSWITHRT 1
J = I (nˆ)dΩ. (10)
ν ν
In this workwe solvethe usualequations of compress- 4π
Z
ibleMHD,including the sourceterminthe energyequa- Whenabsorptiondominatesǫ 1andS B ,but
ν ν ν
tionto accountfor heating and coolingdue to radiation. when scattering dominates ǫ →0 and S →J . Note
ν ν ν
These sourceterms are computed directly froma formal that this expressionassumes th→at scattering→is isotropic.
solutionofthetime-independentRTequation. Thus,the Although this is not strictly true for many scattering
basic equations are continuity processes(e.g. electronscattering),itwill generallybe a
∂ρ good approximationfor problems of interest.
+ (ρv)=0, (1) In addition to J we will also use H and K , the first
∂t ∇· ν ν ν
andsecondmoments,respectively. Theircomponentsare
momentum conservation
given by
∂(ρv)
+ (ρvv+T)=0, (2) 1
∂t ∇· Hνi=4π Iν(nˆ)µidΩ, (11)
the induction equation Z
1
∂B (v B)=0, (3) Kνij=4π Iν(nˆ)µiµjdΩ, (12)
∂t −∇× × Z
and energy conservation 3 Note that χtνot has units of [cm−1]. Throughout this work
we will use χ for quantities with these dimensions and κ = χ/ρ
∂E + (Ev+T v)=Q . (4) for quantities with dimensions of [cm2/g], but will refer to these
∂t ∇· · rad interchangeablyasopacities.
3
where dΩ is the differential of solid angle, and µ nˆ solutionofradiationtransfer. The computation ofRT is
i
≡ ·
xˆ . These moments are related to the radiation energy describedin Section 3 and the interface of the RT solver
i
density E , radiationflux F , and radiation pressure and MHD integrator is described in Section 4. The se-
rad rad
P via the standard definitions quence for a single timestep can summarized as follows:
rad
1) Using the hydrodynamic variables (typically T and
E =4π ∞J dν, (13) ρ)fromtheprevioustimestepasinputs,wecomputeχtot,
rad ν ν
c Z0 ǫν, and Bν, or each frequency in each grid zone.
2)We solveEquation(8) using the methods described
∞
F =4π H dν, (14)
rad ν in Section 3, yielding S and J everywhere in the do-
ν ν
Z0 main.
Prad=4π ∞Kνdν. (15) 3)Using Sν andJν (orHν),we compute the radiation
c Z0 source term Qrad and update equation (4) as described
in Section 4.
Integration of equation (8) over solid angle yields
4) We advance the MHD variables using the standard
Athena integrators.
F =4π ∞χtot(J S )dν. (16)
−∇· rad Z0 ν ν − ν 3. SOLUTIONOFRADIATIONTRANSFER
andprovidesanalternative(differential)formforthera- An extensive literature on the solution of RT for
diationsourceterminequation(4). Thedifferentialform astrophysical problems in multidimensions exists and
tends to perform better in regions where optical depths there are numerous monographs and review articles on
across a gridzone are large, while the integral form is the topic (e.g. Mihalas & Mihalas 1984; Castor 2004;
preferableinregionsoflowopticaldepth. Hence,wewill Carlsson 2008). With this literature to draw from, we
use both expressions, as discussed in section 4. havelargelyadoptedastrategyofimplementingexisting,
We have notbeen forcedto make distinctions between well-developed algorithms. Since there are many differ-
the Eulerian and comoving frame for radiation vari- entmethodswithdifferentstrengthsandweaknesses,the
ables as we have dropped all velocity dependent terms majorchallengeisfinding a methodwhichbestsuits our
in equations (7), (8), and (16). We neglect these terms particular needs. Our most salient constraints include:
because they are negligible for the tests considered in 1) The method needs to be amenable to domain de-
this paper. However, we anticipate solving problems composition since this is the primary algorithm for par-
where the velocity dependent terms may be important allelizing the solution of the MHD equations in Athena.
and can implement terms that are first order in v/c in 2) The method must be able to handle the explicit
our RT solver, where necessary. For consistency with dependence of the source term on J in equation (9) for
ν
the VET solver (JSD12), we will adopt the mix frame problems in which scattering is important (i.e. we must
approach where Iν, its moments, ν, and nˆ are Eu- be able solve non-LTE problems).
lerian frame variables, while opacities and emissivities 3) The method needs to be able to handle (shearing)
are defined in the comoving frame. Derivations of the periodic boundary conditions.
mixed frame equations can be found in Mihalas & Klein 4)Themethodmustberobustandcapableofhandling
(1982), Mihalas & Mihalas (1984), Lowrie et al. (1999), discontinuities in temperature and density which arise
and Hubeny & Burrows (2007). when shocks are present in the flow.
Since we neglect the time derivative of Iν and terms 5) Ideally, the method should be efficient enough that
that are first order in v/c in equation (8), our method for simple problems (e.g. LTE with grey or mean opaci-
is only formally reliable in the static diffusion and free ties), neither the memory constraints nor the total com-
streaming-limits. Specifically,thetimescaleforfluidflow putational time is dominated by the solution of RT.
tf L/v across a characteristic length scale L in the With these considerations in mind, we have
∼
simulationdomainmustbe longerthanthe time ittakes implemented a short-characteristics based solver
for radiation to diffuse tdif L2χtot/c or free-stream (Mihalas et al. 1978; Olson & Kunasz 1987;
∼
tfs L/c acrossthe domain (see e.g. Mihalas & Mihalas Kunasz & Auer 1988). In this method the specific
∼
1984)). Thisissufficientforthetestproblemsconsidered intensityisdiscretizedonasetofraysateachcellcenter
here and should be adequate for many of the problems in the simulation domain. Equation (8) is integrated
ofprimaryinterestto us. Whennecessary,wecanretain along each ray using initial intensities interpolated
terms first order in v/c in equations (7) and (8) and the from neighboring grid zones. Since only neighboring
code will be formally accurate in the dynamic diffusion grid zones are used for this integration, the total
limit (tf .tdif) as well. computational cost (per iteration) scales linearly with
Throughout this work Bν is assumed to correspond the number of gridzones in the domain. This is also
to the Planck function and is a function only of ν and simple to parallelize with domain decomposition as only
gas temperature T. We assume an ideal equation of information from cells on the faces of the neighboring
state with p = ρRT and gas thermal energy density sub-domains need to be passed.
Egas = p/(γ 1). Here R is the gas constant and γ This is in contrast to a long characteristics method
−
is the adiabatic index. The adiabatic sound speed is (e.g.Feautrier1964),whichwouldintegratetheRTequa-
a= γp/ρ. tion along each ray through all gridzones intersected by
The methods for solving the MHD equations without the ray until the edges of the simulation domain are
p
the radiation source term are described in detail in pre- reached. Such a method is generally more computation-
vious publications (Gardiner & Stone 2008; Stone et al. allyexpensivesincecomputationofthe specificintensity
2008;Stone & Gardiner2009)andareunchangedbythe in each gridzone typically requires integrating equation
4
(8)through N1/3gridzones(whereN isthetotalnum- compositionfor parallelization. Hayek et al. (2010) used
∼
ber of gridzone in the domain). It is also more cumber- their code to solve the RT equation including scatter-
some to use with domain decomposition (see, however, ing,inMHDsimulationsofstellaratmospheresonthree-
Heinemann et al. 2006) since it may require the passing dimensional Cartesian grids. Hence, the effectiveness
oflargerblocksofdata,including informationfromnon- of several key aspects of our module have already been
neighboring subdomains. demonstrated in a sophisticated MHD code and applied
Although short characteristic methods are computa- to realistic astrophysicalapplications.
tionally more expedient, they suffer from greater nu-
merical diffusion due to the interpolation that is re- 3.1. Frequency Discretization
quiredtocomputetheintensityinneighboringgridzones
The scheme we have implemented allows for the com-
(Kunasz & Auer 1988). For problems where a few grid-
putation of frequency dependent, grey, or monochro-
zones (or point sources) dominate the total emissivity, a
maticRT.Radiationvariables(moments andspecific in-
short characteristics solver may require very high angu-
tensities) and radiative properties of the fluid such as
larresolutionto accuratelyresolvethe radiationfieldfar
the opacities, thermal source function, and photon de-
from the dominant source. If the angular resolution is
struction parameter are tabulated on a grid of n 1
too low, anomalous structure (e.g. spokes) in the heat- f ≥
discrete frequencies or frequency groups. For flexibility,
ing and cooling rates will emanate from the dominant
the functional form of opacities and emissivities can be
sources (see e.g. Finlator et al. 2009). (In this case the
specifiedvia user-definedfunctions. Ingeneral,the com-
numerical diffusion introduced by interpolation can be
putational cost and memory footprint of problems scale
beneficial.) Instead, the emission from point sources is
linearly with n .
better handled by suitably designed long characteristics f
These frequency bins can simply be discrete frequen-
methods (Abel & Wandelt 2002; Krumholz et al. 2007b,
cieswhenRTisusedtogeneratediagnosticoutputssuch
althoughseealsoRijkhorst et al.2006). Fortheapplica-
as images and spectra. Group mean opacities and emis-
tions of interest in this work (e.g. accretion flows), the
sivities (e.g. Mihalas & Mihalas 1984; Skartlien 2000)
diffuse radiation field dominates. Moreover, even when
and correspondingquadrature weights must be specified
point sources are present, the diffuse radiation field due
when the RT solver is used to compute the radiation
to scattering or re-emission (e.g. HII regions) cannot
sourcetermsorVET.Inthesimplestcase,n =1andan
generally be ignored, and therefore we anticipate such f
appropriate frequency integrated mean opacity is speci-
problems may be accommodated in the future by a hy-
fied.
brid scheme which uses short characteristics for the dif-
Unless otherwise noted, we will drop subscripts de-
fuse emission, and long characteristics for bright point
noting the frequency dependence of radiation variables
sources.
andonlydescribethemonochromaticproblemhereafter.
Non-LTE problems are handled via iteration. For
For the problems under present consideration, there is
each time step the formal solution of the whole do-
no explicit coupling of the specific intensity at different
main is repeated, updating J and S during each
ν ν frequencies so the frequency dependent calculation is a
iteration, until some formal convergence criterion is
trivial generalizationof the monochromatic problem.
met. As discussed below, we implement an accel-
erated (or approximate) lambda iteration (hereafter
3.2. Angular Discretization
ALI) algorithm based on the Gauss-Seidel method
of Trujillo Bueno & Fabiani Bendicho (1995) (hereafter We discretize the specific intensity on both angular
TF95). The TF95 method is efficient for solving non- andspatialgrids. Forone-dimensionalproblems,thedis-
LTE problems because it significantly increases the con- cretization is chosen so that polar angles correspond to
vergence rate without significantly increasing the com- the abscissas for Gaussian quadrature. In multidimen-
putational cost (or memory footprint) per iteration. sions, discretization of the angles proceeds according to
Iteration is also used in LTE problems to handle the algorithm described in Appendix B of Bruls et al.
boundary conditions at the interface of subdomains and (1999),whichisbasedontheprinciplesoftypeAquadra-
for physical periodic boundary conditions at domain ture described in Carlson (1963).
edges. Oneachiterationtheincomingintensityfromthe This method attempts to distribute the rays as evenly
neighboring subdomain is fixed from the previous itera- aspossibleovertheunitsphere,subjecttotheconstraint
tion (or timestep for the first iteration). The outgoing that each octant of the unit sphere is discretized identi-
intensity, which corresponds to the incoming intensity cally. Hence the angle discretization is invariant for 90
◦
in the neighboring subdomain, is then updated and the rotationsaboutthecoordinateaxes. Thisisdesirablebe-
formalsolutionisiteratedtoconvergence. ForLTEprob- cause Athena is designed to be a general purpose code,
lems, this is not the most efficient method for handling and there is often no preferred direction with which to
the subdomain boundaries (see e.g. Heinemann et al. align the angular grid (as in some atmosphere calcula-
2006). For the moment, we are primarily interested in tions). Without this constraint, the result would gener-
non-LTEproblemswhereiterationisrequiredregardless. ally depend on the orientation of boundary and initial
Wegenerallyfindfairlyrapidconvergence(requiringonly conditions relative to the coordinate axes.
afewiterations)formostofourLTEtestproblemswhen The user specifies the number of polar angles n , and
µ
iterations is used, so this is not a significant limitation. the code generates an array of n rays covering the unit
a
In many respects our short-characteristics RT solver sphere. In one dimension, this corresponds to n = n
a µ
is similar to those of van Noort et al. (2002) and rays because of axisymmetry. For multidimensional do-
Hayek et al. (2010) in that both implement ALI to han- mains n =n (n +2) rays. However,in two-dimension
a µ µ
dle deviations from LTE and both utilize domain de- only half of these are unique due to the implied invari-
5
G H I The short characteristic method (Mihalas et al. 1978;
S+ Olson & Kunasz 1987; Kunasz & Auer 1988) has been
discussed previously by several authors. The basic com-
S+ putation step for a single gridzone in both LTE and
non-LTE problems is illustrated in Figure 1 for the two-
dimensional case. Fluid radiative properties and radia-
tion variables (e.g. χtot, B, ǫ, I , J, S) are defined on a
D S0 E F k
radiationgrid. Theverticesofthisgridcorrespondtothe
S0 cell centers of the MHD domain so that fluid radiative
properties are computed directly from the cell centered
I,− S− MHD fluid variables. Generalizing to three dimensional
domains is straight-forward.
At each vertex, the specific intensity I0 at x0 is com-
y k
j A I,− S− B C putedalongeachraynˆk fromx−k tox+k. Forsecond-order
interpolation the intensity is given by
x
Fig. 1.—i Schematic of the RT solution for an individual grid- Ik0=Ik−e(−∆τk−)+Ψ−kSk−+Ψ0kS0+Ψ+kSk+, (20)
zraodnieawtiohnosgercidel.lcIennttehriscocrarseespSo0n,datnodveχr0teaxrEe kinnoawtnwoa-tdiEm,enansidonIa0l where Ψ−k, Ψ0k, and Ψ+k denote interpolation coefficients
vtisieortntoic(be2se0)oc.foStmhinpecugetreiqddu.,atnIht0eiytiisemscuoIsm−tp,buSet±ecdo,maanploduntχged±avdeiaaocinhnotrtearcypoorulraseitsnipognonefqdruotamo- wthheicohptdiecpalenddepotnh itnhteerovpaalsci∆tieτsk−χa−kn,dχ∆0,τak+n.d χ+k through
neighboring grid zones. The red ray intersects rows between ver- The formofthe interpolationcoefficientsΨ−k, Ψ0k, and
ticesAandB(upwind)orHandI(downwind). Hencethevalues Ψ+dependsontheinterpolationmethodused. Thestan-
of S, χ and I at these vertices are used to interpolate S±, χ±, k
andI−. Vertices CandGarealsousedwithquadraticinterpola- dardexpressionsforsecondorderinterpolationarelisted
tion. A similarcase holds forthe blueray, but withinterpolation in equations (7a)-(9c) of Kunasz & Auer (1988). One
performedoncolumnsA-D-GandC-F-I.Theopenandclosedcir- drawback of these expression is that they are subject to
cdlaesshdeedncoutervveeirstiacnesewxtheincshioanreoufstehdefbolruethreayinwtehripcohlainttioenrseocftIs−A.-BT-hCe overshootwhere gradients in Sk± and χ±k are steep. For-
tunately,thesecasescanbehandledwithB´ezier-typein-
row.
terpolationasdescribedinAuer (2003)andHayek et al.
(2010). With B´ezier-type interpolation schemes, one
ance of physical quantities in the third dimension and tcoandeutteilrimzeinaeciofnotvroerlsphooionttsSakcre=prSe0se−nt0i.n5∆tτhk−e(s∂tSan/∂daτr)d0k
na =nµ(nµ+2)/2. second-order expressions. If Skc < min(Sk−,S0) or Skc >
anSaelottgionugsntoµin=vok2ining taheotnwe-od-simtreenasmioanpaplrcoaxlicmulaattiioonn, iins max(Sk−,S0) overshoots are present and alternative ex-
pressions are utilized. Suitable choices follow from set-
wmhaitcehdtbhyetrraadnisafteironalofinegldaosfinegalcehrhaeym. Tishpihsearessiusmapptpioronxiis- ting Skc = Sk− or Skc = S0. Hayek et al. (2010) provide
commonly used to derive analytic solutions, and allows thecorrespondingexpressionsforΨ±k andΨ0kintheirAp-
pendix A. Similar methods are also used to compute in-
theratioofH/J tovarybutkeepstheratioofK/J fixed
at1/3,consistentwiththeEddingtonapproximation. In tervals∆τk± (seee.g. equationA.3ofHayek et al.2010).
two (three) dimensional calculations, choosing nµ = 2 For one-dimensional problems, x−k and x+k correspond
to neighboring grid vertices x and x . Hence,
approximates each quadrant (octant) with a single ray i 1 i+1
−
and also yields Kij = 1/3δijJ. The algorithm is well- Ik− =Ii−1,k,whichwasjustcomputedintheneighboring
defined and unique only for nµ ≤12 (Bruls et al. 1999), zone while Sk±, and χ±k can be computed directly from
corresponding to n =168 (84 in two dimensions). This hydrodynamics variables at x . In multidimensional
a i 1
should not be prohibitive for the problems of interest. problems,x−k andx+k nolonger±correspondto verticesof
For each ray nˆ , we compute a vector of direction
cosines (µ0k,µ1k,µk2k) with µik =nˆk xˆi and quadrature tbheeinrtaedripaotiloantedg.riWd eanimdpvlaermiaebnltesanIdk−,teSstk±b,oatnhdfiχrs±kt-omrduestr
·
weights wk. Then equations (10)-(12) become (linear) and monotonic second order (quadratic) inter-
polationschemes (Auer & Paletou1994). Bothmethods
na 1
− prevent overshoots and enforce positivity of the inter-
J= w I (17)
k k
polants. The choice is particularly relevant for I , as
k
k=0
X second-ordermethodsgenerallyproducemuchlessdiffu-
na−1 sion of the radiation beam. A drawback of second order
H = w I µ (18)
i k k ik interpolation is that it places additional constraints on
Xk=0 theorderinwhichonesweepsthroughgridzonesandthe
na 1 stencil used for the evaluation of I .
− k
Kij= wkIkµikµjk, (19) Consider the two rays depicted in Figure 1. We com-
Xk=0 pute interpolants Sk± and χ±k using known quantities at
vertices of the radiation grid. If row A-B-C or column
where I I(nˆ ).
k ≡ k A-D-Gcorrespondtoghost(boundary)zones,Ik− canbe
computed from the (prescribed) boundary intensities. If
3.3. Implementation of the Short-Characteristics
they are not ghost zones, interpolation can only be per-
Algorithm
formedonzonesinwhichI hasalreadybeencomputed.
k
6
N vertex directly below A (Kunasz & Auer 1988).
These two solutions share common drawbacks. For
parallelization with domain decomposition, only one
ghostzone is needed per gridzone on a subdomain face,
when only vertices A-I are used. Extension of shallow
raysbeyondthisstencilrequiresthepassingofadditional
data and associated bookkeeping. More philosophically,
we feel it is desirable to treat all rays as consistently
as possible. In either of these schemes, RT along some
rayswillbecomputedusingonlyneighboringgridzones,
... while other rays will not. Our preference is to treat all
rays on the same footing.
For this reason, we have decided to switch the order
of the sweep for shallow rays. Athena is implemented
sothat eachsub-gridof the domains has regularspacing
and therefore gridzones with fixed aspect ratio. This
meansthatthedistinctionbetweenraysthatareshallow
and those that are not is equivalent for each grid zone.
y
j However, our definition of a shallow ray depends upon
the direction of the sweep. The blue ray in Figure 1 is
1
1 ... N shallow because we first traverse the grid along rows of
x fixed yj, only moving to yj+1 when intensity has been
i
computed for all gridzones in the row y , as depicted in
Fig. 2.— Progression of the sweep through a two-dimensional j
Figure 2.
grid(domainorsubdomain)whenlinearinterpolationisused. The
forwardsweep(redcurve)firstprogressesparalleltoyˆ,computing Ifwereversethesweepsothatwefirsttraversecolumns
RTonlyforupwardpointingrays(nˆ·yˆ>0). Foreachrow(fixed of fixed x , the blue raywill no longer be shallow,as the
i
yj), one first sweeps parallel to xˆ, computing RT along rays with intensity at G will be computed before it is needed for
anˆn·dxˆc>om0puutnetsilarloeancghrinaygstwheithgrnˆid·xˆed<ge0xuNnt,ilthreeanchrienvgertsheesgdriirdecetdiogne the computation of the intensity at E. In this case the
x0. Thiscontinues untilonereaches gridzone(x0,yN). Theback- redrayisnowashallowrayastheintensityatCwillnot
wardsweep(bluecurves)invertstheforwardsweep,computingRT have been computed before it is needed to compute the
only for downward pointing rays (nˆ·yˆ< 0). In the Gauss-Seidel
intensity at E. Hence, by varying the sweep direction,
method, updated values of Si,j are incorporated into the back-
we can handle all rays and accommodate a quadratic
ward sweep, beginning with SN,N. The three-dimensional case is
astraightforwardgeneralization. interpolation scheme which computes all intensities in a
gridzone (x ,y ) only using intensities from neighboring
i j
gridzones (x ,y ).
i 1 j 1
± ±
If we first sweep along rows of fixed y (as in Figure 2), 3.4. Iterative Methods for Non-LTE Problems
j
I has only been computed at vertices A, B, C, and D.
k Wenowdescribehowwehandlenon-LTEproblemsit-
This means that I is known for all vertices used in the
k eratively. Following common convention we denote the
linearinterpolationofIk−aswellasforquadraticinterpo- angleaveragedformalsolutionoftheRTequation(here-
lation(andanyhigherorderinterpolation)ofrayswhich
after, simply the formal solution) in operator notation
intersect row A-B-C.
as
WerefertoraysthatintersectthecolumnA-D-G,such
J =ΛS. (21)
as the blue one in Figure 1, as “shallow” rays. Shallow
rays are a potential problem for quadratic (and higher Here,Λisalinearoperatorrepresentingthe(discretized)
order) interpolation, since I at G has not been com- formal solution, and J and S are vectors spanning each
k
puted. When quadratic or higher order interpolation is gridzone in the simulation domain. Using equation (9)
desired, such rays can be handled in a number of ways. toeliminate J,oneobtainsanequationfor S intermsof
One possibility is to switch the order of the sweep for B
shallow rays so that it first proceeds in the y direction S =(1 ǫ)Λ[S]+ǫB. (22)
−
along columns of fixed x . In this case I for shallow
i k
Since Λ is a linear operator we can solve for S
rayswillbe knownatverticesA,B,D,andG.Themain
drawback(discussedfurtherinSection3.4below)isthat S =[1 (1 ǫ)Λ] 1[ǫB]. (23)
−
one is unable to implement a Gauss-Seidel iteration for − −
non-LTE problems. If one can invert Λ a formal solution of the non-LTE
One can also construct alternatives by extending the problem follows from solving (23) and obtaining J from
stencilbeyondverticesA-I.Forexample,onecanextend (21). However, for three-dimensional problems Λ is a
shallowrays until they intersectrow A-B-C as shownby very large matrix and not sparsely populated when sys-
thedashedcurveinFigure1. Adrawbackofthissolution tems are far from LTE so its direct inversion is imprac-
isthatitrequiresmodestadditionaleffortforcomputing tical. Therefore, equation (22) is usually solved via iter-
Ψ−k, although this can be alleviated by computing only ation.
on the first iteration and reusing it for subsequent iter- Asimpleiterativeschemeforsolvingequation(22)be-
ations (e.g. Hayek et al. 2010). Alternatively quadratic gins with an initial guess for the source function SN,
interpolation could be preformed using A, D, and the which is then used to compute an improved estimate
7
Sn+1 =(1 ǫ)Λ[Sn]+ǫB. However,this method (often Sn . However, this alone is not sufficient to make
− N,N
referred to as Lambda Iteration) has very poor conver- it a Gauss-Seidel scheme, because the contributions to
gence properties. For practical problems, ALI methods Jn+1 , Jn+1 , and Jn+1 from upward directed
N 1,N N,N 1 N 1,N 1
(Cannon 1973) are commonly used. Rybicki & Hummer ray−s on the for−ward sweep−use−d Sn . These must also
(1991), Hubeny (2003) and TF95 provide useful reviews N,N
be updated using ∆S = Sn+1 Sn and weights
of ALI methods and we refer the reader to these works N,N N,N − N,N
which were saved on the forward sweep. We also up-
for a more in-depth discussion. Here we just summarize
date the outgoing intensities In+1 (since they were also
the basic concepts involved. k
computed using Sn ) as they correspondto the incom-
In ALI methods one solves equation (23) directly, but N,N
using an approximate form Λ which is easier to invert ing intensities in neighboring gridzones. Since the corre-
∗
then the full Λ operator. Since only the approximate spondingweightshavealreadybeencomputedaspartof
Λ is used, iteration is still necessary. Numerous choices the forward sweep, the additional computational cost is
∗
for Λ have been proposed, but it has been argued that very modest.
∗
simply taking the diagonal elements of the full Λ ma- Following the discussion in Section 3.3, we note that
trixrepresentsanear-optimalchoice(Olson et al.1986). feasibility of performing a Gauss-Seidel iteration with
Olson & Kunasz (1987) provide expressions for diago- quadratic interpolation is dependent on the way shal-
nal elements of Λ when short characteristics are used. low rays are handled. Reorienting the sweep for shal-
In each grid zone the change in the source functions low rays so that j is more rapidly varying index, but
∆S =Sn+1 Sn can be written as keepingiasthe rapidlyvaryingindex forremainingrays
i i − i does not allow for an efficient Gauss-Seidel scheme be-
∆S = (1−ǫi)Jin+ǫiBi−Sin, (24) cause some of the necessary Jn+1 (and therefore Sn+1)
i 1 (1 ǫ )Λ are not available when the backward sweep begins.4 In
i ii
− − thelightofthisissue,wehaveimplementedGauss-Seidel
wherethe subscriptienumeratesallgridzonesinthe do-
routines only with linear interpolation. For problems
main.
where quadratic interpolation is preferable, we default
As TF95 discuss, when Jn is exclusively used in equa-
to the Jacobi method (i.e. standard ALI).
tion(24),theALIschemeisequivalenttotheJacobiiter-
We continue the iteration until some convergence cri-
ativemethodforsolvinglinearsystems. TF95showthat
terionis met. Consistentwithpreviouswork,westopit-
one can construct a Gauss-Seidel algorithm by incorpo-
eratingwhenthemaximumrelativechangeinthesource
rating the new values of Jn+1 in equation (24) as these
i′<i function over the whole domain is less than some pre-
become available. Here i′ < i refers to gridzones where scribed threshold δc
J has already been updated. The complexity of devis-
inga Gauss-SeidelalgorithmforRT comesfromthe fact ∆S
i
max | | δ . (25)
tshuabtsetthoefctohmepruatyastinˆon noefesdpetcoifibceincotemnpsiuttyedIi,uksifnogr sSonme (cid:18) Si (cid:19)≤ c
k i′<i
rather than Sn+1 (i.e. old rather than new values of the ForLTEproblemsthatuseiterationtohandleboundary
i′<i
source function). Therefore the contribution from these conditions, S does not change from one iteration to the
particularraystoJn+1mustbecorrectedastheupdated next and we replace Si with Ji in equation (25).
i′<i The choice of δ is clearly an important input to the
values Sn+1 become available. c
i′<i method, but there is no firmly established criterion and
TF95 give a detailed description of how to implement
the optimal choice depends on a number of considera-
suchanalgorithmonaone-dimensionaldomain. Theal-
tionsthatmaybe problemdependent. Since thecompu-
gorithmrequiresstoringamodestamountofdataineach
tationalcostofthe methodgenerallyscaleslinearlywith
gridzone, but very little additional computation. The
thenumberofiterationsperformedandalowerthreshold
convergencerateis improvedbyafactoroftwo,soprob-
leads to more iterations, there is a tradeoff between ac-
lems requiring several iterations gain nearly a factor of
curacy and computational expediency. With the excep-
two decrease in computational effort for only a minor
tion of the uniform temperature non-LTE atmosphere,
increase in code complexity.
the tests presented in section 5 were performed using
When linear interpolation is used, the generaliza- δ =10 5. Increasingδ to 10 4 hada negligibleimpact
tion of their one-dimensional method to two and three- c − c −
on the linear wave tests.
dimensional domains is straight-forward. The two-
Ourexpectationsbasedonthetestswehaveperformed
dimensionalsweepproceedsasdepictedinFigure2. The
sofar arethatfor the problemsofprimaryinterestto us
vertices in the radiation grid correspond to cell centers
(e.g. shearing box simulations of accretion disks) δ
(xi,yj). The sweep generally proceeds with i as the 10 5 10 3 will be sufficient, consistent with studcie∼s
− −
more rapidly varying index. Consider a domain with −
usingsimilarmethods(e.g.Hayek et al.2010). However,
N = N = N for simplicity. In each gridzone (x ,y ),
x y i j weemphasizethattheappropriatechoicewillbeproblem
we first compute the intensity I for all upwarddirected
k dependent andmust be assessedon a case-by-casebasis.
rays (nˆ yˆ > 0) in the forward sweep and then for all
downwakrd· directedrays(nˆ yˆ<0)onthereversesweep. We view the choice of δc in roughly the same terms as
k· we view the choice of grid resolution. One can adopt a
On the reverse sweep, the upper right gridzone
threshold based on previous results and experience, but
(x ,y ) is the firstin which the computationof allnew
N N ultimately one needs to compute the problem using a
intensitiesIn+1 iscompleted. AtthispointJn+1 iscom-
k N,N
pletely specified and we compute Sn+1. From here on, 4Fortwo-dimensionaldomainsonecandeviseanefficientGauss-
N,N Seidelalgorithmthatsweepsdiagonallythroughthegrid,butthis
all subsequent RT computations use Sn+1 rather than implementationdoes notgeneralizetothreedimensions.
N,N
8
range of δ and choose a sufficiently small value such the domain decomposition should not lead to significant
c
that the results are insensitive to the choice. errors.
3.5. Boundary Conditions and Parallelization 4. INTERFACEOFTHERADIATIVETRANSFERSOLVER
TOTHEMHDINTEGRATOR
Boundary conditions and domain decomposition in
There are two regimes in which the effect of radia-
Athena are both implemented for MHD via the use of
tion on the MHD is important. The first is when the
ghost zones, and we implement RT boundary conditions
radiation field is a significant contribution to both the
in an analogous way. The solver computes RT in grid-
energyandmomentumfluxesintheflow. Inthisregime,
zonesonaboundary(domainorsubdomain)inthesame
the radiationsourcetermsinthe MHDequationscanbe
wayasaninteriorgridzone,butusingtheintensitiesand
very stiff, and the equations contain wave modes which
source functions from the ghost zones to compute the
propagate at close to the speed of light. Both of these
relevant integration weights and interpolants. The in-
properties require significant modification to the under-
tensities and source functions in the ghost zones are de-
lying MHD integrators in order to enable accurate and
termined according to prescribed boundary conditions.
stable integration. In JSD12 we describe a method for
In general, boundary conditions for the MHD integra-
this regime based on an extension of the modified Go-
torwillnottranslatedirectly toboundaryconditions for
dunov method of SS10 to multidimensions, with a VET
the RT solver. Different problems with the same MHD
(defined as f = P /E ) computed from a formal so-
boundaryconditionsmayrequiredifferentboundariesfor rad rad
lution of the RT equationusing the module described in
the radiation field. Hence separate boundary conditions
this paper. At each time step, the RT solver computes
must be prescribed when using the RT solver. For the
theradiationfieldasdescribedinSection3,evaluatingK
code test problems presented in section 5, we have im-
and J via equations (17) and (19). We the compute the
plemented two types of boundary conditions specifying
VETusingf =K/J asdescribedinsection3.4ofJSD12.
either fixed incident intensity or periodic intensities on
The second regime is when the radiation pressure can
the boundaries. Otherboundaryconditionscanbespec-
beignored,andtheeffectofradiationisonlythroughthe
ified via user defined functions.
heatingandcoolingsourcetermsintheenergyequations.
Athenarunsonparallelmachinesusingdomaindecom-
In principle, the modified Godunov method adopted in
positionimplementedthroughMPIcalls. TheMHDinte-
JSD12wouldbe anattractive approachfor handling the
grator passes all conserved variables and passive scalars
stiff energy source term that can arise in this regime as
from faces of neighboring subdomains to ghost zones.
well. However, the modified Godunov method requires
The MHD integrator requires either four or five ghost
that one compute the gradientof radiationsource terms
zones for each gridzone on the subdomain face. The
onthe planeofprimitive variables. This inturnrequires
RT solver operates analogously, passing intensities and
analytic expression for the radiation sources in terms of
source functions, but only requires one ghost zone for
the fluid variables. Hence, it is generally not a viable
each gridzone on a subdomain face.
method for problems where the radiation properties are
The main differences between the RT solver and the
complicated functions of frequency and fluid variables,
MHD integrator are the frequency and quantity of data
as may be the case with bound-free and bound-bound
that must be passed. For each frequency bin in every
atomic opacities or Compton scattering.
ghostzoneI mustbepassedforalln rayswithquadratic
a These limitations motivate us to implement an alter-
interpolation or, alternatively, n /2 incoming rays with
a native method to directly compute the radiation source
linear interpolation. For non-LTE problems S and H
term in the fluid energy equation (4) and couple it to
must also be passed. Hence for quadratic interpolation,
the standard MHD integrators. When operating in this
thecodepassesatotalofn (n +1+n )floatingpoint
f a dim mode, we perform the formal solution at the beginning
variables per face gridzone per iteration, where n is
dim of each timestep. We first compute fluid radiation prop-
the number of dimensions in the domain. In contrast,
erties in each gridzone x of the domain. This includes
the MHD integrator typically passes 50 floating point i
variables per face gridzone per times∼tep. For problems the variables χtiot, Bi, and ǫi, which are computed via
user-defined functions of the conserved MHD variables
wheren andn aresmallandfewiterationsarerequired
a f and passive scalars from the previous timesteps. We use
(e.g. an LTE grey problem), the volume of RT data is
these, along with J from the previous time-step, to ini-
therefore comparable to and may even be less than the i
tialize S . Once the formal solution is completed, we
amount of data passed by the MHD integrator. i
accountfor the source function onthe righthandside of
We note that the use of iteration to handle subdo-
equation (4) via an operator split update of E. We first
mainboundaryconditionsmayleadtosomedependence
compute the radiative source function in each zone and
on the number of subdomains that are used. We have
then update the total energy
considered the sensitivity of our results to this issue by
performing most of the tests described in section 5 both ∆Ei =δt(Qrad)i. (26)
with and without domain decomposition. In practice, The standardMHD integrationalgorithmthen proceeds
the converged mean intensities do not differ (relative using this “new” value for E .
i
to the non decomposed domain) by more than ∼ δcJ. WecomputeQradinoneoftwoways,dependingonthe
The sensitivity is highest for problems where the optical characteristic optical depth. We either use the integral
depthacrossanindividualsubdomainisoforderunityor form
smaller, Problems with optically thick subdomains gen- Qint =4πχtot(J S )=4πχabs(J B ), (27)
erally lead to smaller discrepancies. Since we already i i i− i i i− i
or the differential form
choose our convergence criterion to be at a level that
minimizes the impact on our results, this sensitivity to Qdif = 4π H , (28)
i − ∇· i
9
Previouswork(Bruls et al.1999, andreferencestherein)
has demonstrated that the integral form is inaccurate
whentheopticaldepthpergridzoneislarge. Inthiscase
J B B while χabs is large so round-off errors can
i− i ≪ i i
begreatlyamplified. Theintegralform,however,ismore
accurate when χabs∆x .1 (Bruls et al. 1999).
i i
Therefore, we have designed our RT solver to com-
pute eitherformofQ , depending onthe regimeofthe
rad
computation. In most applications of interest, there is a
transition from optically thick to optically thin regions,
so we must specify a criterion for switching between the
differential and integral forms in the same domain. For
thetestproblemsconsideredhere,wefindasimpleswitch
Qint if χtot∆x 1
(Qrad)i = Qidif otheirwisei ≤
(cid:26) i
tobesufficient. Thishastheadvantagethatitisapurely
localcriterion. Usingamethodwhichmoresmoothlyin-
terpolatesbetweenthetworegimes(seee.g.Hayek et al.
2010)didnotimproveperformanceinameasurableway,
but may be preferable for more sophisticated applica-
tions. Fig. 3.— Convergence of the ALI methods on a highly non-
Duetotheexplicitupdate,wemusttakecareinchoos- LTE uniform atmosphere with ǫ = 10−6. The curves show the
maximumrelativedifferencebetweenthenumericallycomputedS
ing a time step. In the absence of radiation the MHD
andtheanalyticsolutionsimpliedbyequation(30)fortheJacobi
integratorchoosesatimestepδtCbasedontheCFLcon- (solid) and Gauss-Seidel (dotted) methods. We compute the nu-
straint. In principle, this time step can be much larger mericalsolutionsusingaone-dimensionaldomainwith9gridzones
perdecadeinopticaldepth.
thantheradiativecoolingtime,whichcouldleadtoobvi-
ouserrors,suchastheenergydensitybecomingnegative.
Asweelaborateuponinsection5.4,onecanderiveagen- tests of the coupled MHD integrator and RT solver in
eralizedCFL conditionforaradiatingfluidbasedonthe fully time dependent calculations. The former are par-
need to resolve the damping time for a non-equilibrium ticularlyusefulforevaluatingtheRTsolversperformance
radiation diffusion mode. This time scale δtrd is gener- onmultidimensionalandnon-LTEproblems. Forthelat-
allymostrestrictivewhentheopticaldepthpergridzone ter, we focus primarily on simpler LTE problems, so we
χtot∆x 1, in which case can compare the simulations result directly to precise
∼
E a analytic or semi-analytic solutions.
gas
δtrd δtC, (29) Further tests of the RT solver as part of the VET
≈ E c
rad method are presented in JSD12.
assuming the adiabatic sound speed a (rather than the
Alfv´en speed) sets the CFL condition. This generalized 5.1. Uniform Temperature Non-LTE Atmosphere
CFL constraint can be quite stringent, requiring short WebeginbysolvingthemonochromaticRTproblemin
time steps and increasing computational costs if either auniformtemperature,one-dimensionalscatteringdom-
a c or Erad & Egas. Hence, many problems will re- inated atmosphere. This test is particularly useful for
≪
quire the use of the VET method described in JSD12, evaluating the RT solver’s performance on a non-LTE
which uses timesteps determined by the standard (non- systemsandevaluatingtheconvergencepropertiesofJa-
radiative) CFL constraint. In practice, we are almost cobi and Gauss-Seidel iterative schemes. We adopt the
always limited to problems with Egas < Erad, so we do two-stream approximation for the RT solver so we can
not attempt to include the radiation momentum source compare directly with analytic solutions based on the
terminequation(2) asitis generallysmallfor problems Eddington approximation. Since we assume a uniform
thatarecomputationallyfeasiblewithoperatorsplitting. opacity κ and temperature T, the analytic solution is
The algorithm described above will, in general, only only a function of optical depth dτ = χdz, the thermal
be first order convergent. Note that we could construct source function B and photon destruction probability ǫ.
a second-order scheme when using Athena’s VL+CT WiththeseassumptionsthemeanintensityJ isgivenby
integrator (Stone & Gardiner 2009), by performing the
operator split update before the corrector step in the e−√3ǫτ
J =B 1 . (30)
predictor-corrector scheme. However, some of the ad- − 1+√ǫ!
vantagesofthesecondorderconvergencewillbelostdue
to the increased diffusivity of the VL+CT relative to Weassumethatχ ρandρincreasesexponentially(but
∝
the CTU+CT scheme (Gardiner & Stone 2008). Hence, keepǫconstant)withdistancefromtheupperboundary,
we have not yet pursued the possibility although it may whichhas no incoming intensity. This providesan expo-
prove to be a useful avenue for future work. nential variation in τ which is well-suited for resolving
the transition from LTE to non-LTE within the atmo-
5. TESTS sphere.
Our test problems fall into roughly two categories: Figure3showstheconvergenceofthetrueerrorofthe
stand-alone tests of the RT solver on fixed domains and numerically derived solutions. This is evaluated as the
10
Fig. 4.— Comparisonofnumerical(solid)andanalytic(dashed)
solutionsofthesourcefunctionmonochromatic,uniformnon-LTE
atmospheresasafunctionofopticaldepth. Eachsetofcurvescor-
respondstoadifferentphotondestructionparameter,runningfrom
ǫ=10−2to10−10fromtoptobottom. Wecomputethenumerical
solutions simulations using cubic three-dimensional domains with
643 gridzones, distributed over 64 MPI subdomains. The optical Fig. 5.— Comparison of Athena (diamonds) and Feautrier
depthvariationisalignedwiththez-axisofthesimulationdomain (solid) emission spectra from the upper boundary of a one di-
andthesolutionisuniforminthehorizontaldirections. mensionalatmosphere. Bothcalculationsassumeisotropicelectron
scatteringandfree-free(Bremsstrahlung)absorptionandemission
for a completely ionized H plasma. The intensities are computed
using the same angular grid corresponding to abscissas of a 16
maximumrelativedifference ∆S /S,with∆S thediffer- point Gauss-Legendre quadrature of the interval (1,1). The plot-
| | ted intensities (from top to bottom) correspond to cosi =0.10,
ence of the numerically derived S from the analytic so-
0.28, 0.46, 0.62, 0.76, and 0.99. The atmospheres have constant
lution. We first consider a one-dimensionaldomain with temperature (106 K) and density which varies exponentially with
ǫ = 10−6, as this gives a highly non-LTE atmosphere distance,risingfrom10−6 g/cm2 attheupper(surface)boundary
and facilitates direct comparisonwith Figure 3 in TF95. to 10−4 g/cm2 and lower boundary. For comparison, we plot the
We initialize the radiationfieldto be inLTEeverywhere correspondingblackbodyat106 Kasadashedcurve.
(J = B). We consider two different iterative schemes:
Jacobi and Gauss-Seidel. As expected, the convergence
rate of the Gauss-Seidel methods is nearly a factor of
two better than Jacobi. We assume nine gridzones per solutiontends to increase as ǫ decreasesand the domain
decade in τ to match TF95 and our convergence rates deviates more strongly from LTE. The accuracy of the
agree reasonably well with those shown their Figure 3. numerical solution improves with increasing resolution,
We have also implemented the successive over- but the number of iterations needed for convergence in-
relaxation (SOR) method of TF95, and find rapid con- creases roughly linearly with resolution. The number
vergence, consistent with that shown in Figure 3 of of iterations required for convergence also increases as ǫ
TF95. WehavetestedSORonbothone-dimensionaland decreases. Hence, greater deviations from LTE require
two-dimensional domains and find that it is an effective a greater number of iterations for convergence, as one
method aslong asallboundary intensities arefixed dur- would expect.
ing iteration. However, if the intensities on one of the Although some RT problems do require explicit fre-
boundaries vary from one iteration step to the next, the quency coupling (e.g. Compton scattering, partial re-
method is generally not stable. For example, instabil- distribution), many problems can be treated in the ap-
ity occurredwhenwe usedperiodic boundary conditions proximation that frequencies are not explicitly coupled.
or when we employed subdomain decomposition. Since Multifrequency problems are then just a series of sin-
most of our primary science goals involve problems that gle frequency calculations and, hence, a straightforward
require the use of periodic boundary conditions or do- generalization of the monochromatic problem. Figure 5
maindecomposition,wedonotconsiderSORagenerally shows the intensity spectrum from a multifrequency cal-
viable method for our work. Nevertheless, it may be an culationdonewithAthenaforauniformtemperatureat-
effective method for a modest sized problemthat can be mosphere. We again assume ρ varies exponentially with
run serially with fixed boundary intensities. distance,rising from10 6 g/cm2 at the upper boundary
−
We next consider the same test problem, but use a to 10 4 g/cm2 and lower boundary. The results plotted
−
cubic three-dimensional domain with n = 2. We align here are for N =256.
µ x
the variation of density with the z axis of the domain Incoming intensity at the upper boundary is assumed
anduse periodic boundaries in the horizontaldirections. to be zero and I = B at the lower boundary. We in-
ν ν
Figure 4 shows a comparison of the numerical and ana- clude isotropic electron scattering opacity and free-free
lytical solutions for various choices of ǫ. The agreement (Bremsstrahlung)emissionandabsorption. Theelectron
betweenthenumericandanalyticsolutionsisquitegood scatteringismodeledasisotropicandthecross-sectionis
overall, but tends to be poorest at low optical depths. the Thomson cross-section. For simplicity free-free pro-
For fixed resolution, the discrepancies with the analytic cesses are computed assuming a Gaunt factor of unity.