Radiation Transport for Explosive Outflows: Opacity Regrouping
DD RAFT VERSION O CTOBER
16, 2014
Preprint typeset using L A TEX style emulateapj v. 5/2/11
RADIATION TRANSPORT FOR EXPLOSIVE OUTFLOWS: OPACITY REGROUPING R YAN
T. W
OLLAEGER , AND D ANIEL R. VAN R OSSUM Department of Nuclear Engineering & Engineering Physics, University of Wisconsin, Madison 1500 Engineering Drive, 410 ERB, Madison, WI, 53706;[email protected] and Flash Center for Computational Science, Department of Astronomy & Astrophysics, University of Chicago, Chicago, IL, 60637; daan@flash.uchicago.edu
Draft version October 16, 2014
ABSTRACTImplicit Monte Carlo (IMC) and Discrete Diffusion Monte Carlo (DDMC) are methods used to stochasticallysolve the radiative transport and diffusion equations, respectively. These methods combine into a hybridtransport-diffusion method we refer to as IMC-DDMC. We explore a multigroup IMC-DDMC scheme that,in DDMC, combines frequency groups with sufficient optical thickness. We term this procedure “opacityregrouping”. Opacity regrouping has previously been applied to IMC-DDMC calculations for problems in whichthe dependence of the opacity on frequency is monotonic. We generalize opacity regrouping to non-contiguousgroups and implement this in
SuperNu , a code designed to do radiation transport in high-velocity outflowswith non-monotonic opacities. We find that regrouping of non-contiguous opacity groups generally improvesthe speed of IMC-DDMC radiation transport. We present an asymptotic analysis that informs the nature ofthe Doppler shift in DDMC groups and summarize the derivation of the Gentile-Fleck factor for modifiedIMC-DDMC. We test
SuperNu using numerical experiments including a quasi-manufactured analytic solution,a simple ten-group problem, and the W7 problem for Type Ia supernovae. We find that the opacity regrouping isnecessary to make our IMC-DDMC implementation feasible for the W7 problem and possibly Type Ia supernovasimulations in general. We compare the bolometric light curves and spectra produced by the
SuperNu and
PHOENIX radiation transport codes for the W7 problem. The overall shape of the bolometric light curves are ingood agreement, as are the spectra and their evolution with time. However, for the numerical specifications weconsidered, we find that the peak luminosity of the light curve calculated using
SuperNu is ∼
10% less thanthat calculated using
PHOENIX . Subject headings: methods: numerical radiative transfer stars: evolution supernovae: general INTRODUCTIONType Ia supernovae (SNe Ia) are the explosions of Carbon-Oxygen (C-O) white dwarf stars. In the most widely studiedmodel of SNe Ia, a C-O white dwarf approaching the Chan-drasekhar mass releases energy from nuclear fusion that ex-ceeds gravitational binding energy of the star, causing thestar to explode (Branch & Khokhlov 1995). The resultinghigh-velocity outflow becomes ballistic in a matter of minutes,and thereafter expands homologously. During this expansion,gamma rays from the radioactive decay of Ni heat the out-flow, causing it to radiate, with a peak luminosity that canexceed the host galaxy of the supernova.The majority of observed SNe Ia have similar peak lumi-nosities and spectra (Hillebrandt & Niemeyer 2000). The lightcurves of most SNe Ia obey a peak luminosity-width relation-ship (Phillips 1993). As a result, the light curve data for SNeIa may be fit to a template, enabling its peak luminosity, andtherefore its relative distance, to be determined. Consequently,SNe Ia are important “standard candles” for measuring cosmicdistances and the expansion rate of the universe, and their usefor these purposes led to the discovery of dark energy [see,e.g., Riess et al. (1998); Perlmutter et al. (1999)].Given the significance of SNe Ia in galaxy formation andevolution (Scannapieco et al. 2008) and in nucleosynthesis, aswell as in cosmology, much research has been done to under-stand how model parameters affect the observable propertiesof these events; for example, the connection between explosionasymmetry and anomalies in luminosity (Calder et al. 2004;Kromer & Sim 2009). Other research efforts have focused ongenerating methods, algorithms, and codes that can adequatelytreat the physics of SNe Ia, along with other hydrodynamic and radiative events in astrophysics. Mihalas & Mihalas (1984,pp. 128,144,160) derived the equations of relativistic fluidflow. Castor (2004, pp. 41,49) describes standard Lagrangianand Eulerian methods to solving hydrodynamic problems. The
FLASH code (Fryxell et al. 2000; Calder et al. 2002) provides ameans of solving the Euler equations for compressive, reactivehydrodynamics with nuclear reactions.Radiation transport in Type Ia SNe is a complex problemboth theoretically and practically. From the theoretical per-spective, photons may interact with millions of spectral linesin a heterogeneous material that has multiple ionization states[see, e.g. van Rossum (2012)]. A photon may see an opticallythin environment in one location of the outflow and subse-quently redshift into resonance with a line opacity elsewhere.Such situations provide a challenge to Local ThermodynamicEquilibrium (LTE) calculations, and especially, Nonlocal Ther-modynamic Equilibrium (NLTE) calculations. There is alsothe question of the leading-order behavior of the radiation atdifferent time scales in the presence of material fluid. Lowrieet al. (2001) make the distinction between the radiation timescale and the fluid time scale as a means of preserving correctrelativistic principles in first-order comoving transport.From the practical perspective, high-fidelity Type Ia SNesimulations are generally seen to be demanding in memoryand algorithm efficiency (Baron & Hauschildt 2007). Foran end-to-end simulation, one needs to couple a progenitorexplosion-phase hydrodynamic simulation to the beginning ofthe homologous-expansion phase, and then appropriately treatradiation transport in the latter (Seitenzahl et al. 2013; Longet al. 2013). Numerical simulations of the full evolution ofthe supernova, regardless of the particular explosion model, a r X i v : . [ a s t r o - ph . H E ] O c t involve a large range of densities, temperatures, length scales,time scales, and physical phenomena.Codes can apply transport theory to the homologous-expansion phase of Type Ia supernovae to synthesize lightcurves and spectra. Broadly speaking, transport calcula-tions may be performed deterministically with some subsetof matrix-solution techniques or stochastically with random-number sampling. The stochastic approach gives terms inthe transport equation a probabilistic interpretation; this givesrise to “particles” with sampled properties that can be manip-ulated and tallied to solve the transport equation. Commonmethods of computational transport described by Lewis &Miller (1993) include: discrete ordinates (Lewis & Miller1993, pp. 116,156), integral transport (Lewis & Miller 1993,p. 208), multigroup (Lewis & Miller 1993, p. 61), and finite ele-ments (Adams 2001). The listed methods may be implementedin or in conjunction with Monte Carlo (MC) or deterministicschemes (Urbatsch et al. 1999); the resulting scheme might bedeemed a composite method.Several radiation transport codes have been developed andapplied to the W7 model of Nomoto et al. (1984) and to SNIa models generally. Deterministic codes include PHOENIX ,a code based on the iterative, short characteristic method(Hauschildt & Baron 1999; Baron & Hauschildt 2007; Ol-son & Kunasz 1987). Recently, van Rossum (2012) extended
PHOENIX to be able to calculate self-consistently the tem-poral evolution of the SN Ia outflow. Hauschildt & Wehrse(1991) investigate a discrete ordinates method that incorporatesrelativistic effects to be able to treat explosive outflow. TheMC codes
SEDONA of Kasen et al. (2006), the code of Lucy(2005), and the
ARTIS code of Kromer & Sim (2009) solvemulti-dimensional, time-dependent radiation transport in ho-mologous outflow. Kasen et al. (2006) and Kromer & Sim(2009) solve multifrequency transport by applying the Solobevapproximation (Castor 2004, p. 122) to line transport.Monte Carlo in the context of a velocity field has the fa-vorable property that particles (which are also referred to aspackets) may be tracked in one inertial (lab) frame and interactwith the fluid in the comoving frame. A particle may have itsproperties converted to the comoving frame, updated accord-ing to the interaction, and converted back to the lab frame ifthe particle history is not discontinued. Kasen et al. (2006)applies MC iteratively within a time step to obtain convergedelectron temperatures while Kromer & Sim (2009) find thecontribution of MC iteration to be insignificant if small timesteps are chosen.Instead of treating the temperature structure iteratively orexplicitly, there exist transport methods that are made fullyimplicit (N’Kaoua 1991; Brooks 1989) or semi-implicit (Fleck& Cummings 1971; Carter & Forest 1973) through time dis-cretization of the material equation(s) and adjustment of MonteCarlo interpretations (Densmore & Larsen 2004). To ourknowledge, these methods have not been extensively examinedfor application in the SN Ia problem.Implicit Monte Carlo (IMC) is a stochastic method that maybe applied to solve the time-dependent, nonlinear radiationtransport equations (Fleck & Cummings 1971; Fleck & Can-field 1984). Of the implicit methods referenced towards theend of the preceding paragraph, IMC is quite possibly the sim-plest to implement. The IMC method is made semi-implicitthrough a non-dimensional quantity, referred to as the Fleckfactor, that converts a portion of absorption and reemission to instantaneous “effective scattering.” By introducing effectivescattering, the Fleck factor stabilizes large-time-step radiationtransport calculations that might otherwise suffer significantnon-physical temperature fluctuations (Fleck & Cummings1971). However Larsen & Mercier (1987) demonstrate thatIMC may still be prone to spurious temperature fluctuationfor large time steps and derive a sufficient but not necessaryconstraint on time step size to prevent non-physical behav-ior, which they call the “Maximum Principle” (MP). Recentextensions have been made to IMC that mitigate the patholo-gies associated with the MP [see, e.g., McClarren & Urbatsch(2009); Gentile (2011); McClarren & Urbatsch (2012)].IMC may suffer in performance when effective scatteringdominates over other particle processes. Performance may beimproved for calculations having significant physical or effec-tive scattering by combining IMC with either a deterministicor stochastic diffusion method. Stochastic methods includeRandom Walk (RW) (Fleck & Canfield 1984), Implicit MonteCarlo Diffusion (IMD) (Gentile 2001; Cleveland et al. 2010),and Discrete Diffusion Monte Carlo (DDMC) [see, e.g., Dens-more et al. (2007, 2008, 2012)]. The methods listed have beenhybridized with IMC and applied to both grey and multifre-quency or multigroup problems. Additionally, each methodmay benefit IMC by replacing small-mean-free-path particleprocesses with large diffusion processes. The larger diffusionsteps of the RW method developed by Fleck & Canfield (1984)place a diffusive particle isotropically on the surface of a sphereof several mean free paths in radius centered at the particle’sinitial position. This sphere must be bounded by the spatialgrid that stores the material properties (Fleck & Canfield 1984).Hence, histories in diffusive domains near cell boundaries willnot have sufficiently large displacement spheres; this is foundto limit the increase in IMC efficiency (Densmore et al. 2008).DDMC and IMD differ from RW by discretizing the diffu-sion equation in space; after some algebra, the resulting termsare given a Monte Carlo interpretation (Densmore et al. 2007;Gentile 2001). The discretization implies that a DDMC parti-cle position within a spatial cell is ambiguous (Wollaeger et al.2013). IMD discretizes the diffusion equation in time whileDDMC keeps particle time continuous. Continuous particletime precludes causal ambiguity for each particle (Densmoreet al. 2007).The hybridization of IMC and DDMC, referred to as IMC-DDMC, has been investigated in multigroup problems (Dens-more et al. 2012; Abdikamalov et al. 2012; Wollaeger et al.2013). In each of the IMC-DDMC implementations, thereis a mean-free-path threshold that dictates whether or not acell and group of the spatial and wavelength grids is amenableto diffusion theory. Densmore et al. (2012) investigate a hy-brid for monotonic opacity dependence on frequency that ap-plies grey DDMC in a “large” lower group below a frequencythreshold and multifrequency or multigroup IMC above thefrequency threshold. Abdikamalov et al. (2012) describe ageneral multigroup IMC-DDMC scheme for application toneutrino transport in the presence of a fluid; this makes themethod velocity dependent. Wollaeger et al. (2013) delin-eate a velocity-dependent method for photons that reconciles Note the Fleck factor is not a directly tunable parameter but followsnaturally from linearizing the thermal transport equations within each timestep. Roughly speaking, time steps that result in the deposition of a radiationenergy density that is greater than or of order the material energy density maycause an IMC simulation to become unstable; hence the pathology dependson the evolution of the radiation field (Gentile 2011).
IMC-DDMC to high-velocity, homologous Lagrangian grids.Here, we present some extensions to the particular IMC-DDMC method described by Wollaeger et al. (2013). Theextensions are opacity regrouping (Densmore et al. 2012)and the Gentile-Fleck factor (Gentile 2011). We implementthese features in the IMC-DDMC radiation transport code,
SuperNu (Wollaeger et al. 2013). We first briefly discussthe thermal radiation transport equations. Then we applyan asymptotic analysis to the continuous, comoving trans-port equation on an interior of a frequency domain and in aboundary layer of a frequency domain; this clarifies wherethe DDMC redshift scheme is generally applicable. We sum-marize standard IMC, the Gentile-Fleck factor modified IMCscheme (Gentile 2011), and the hybrid IMC-DDMC equations.Next, we discuss IMC-DDMC processes and a scheme forcombining groups that have DDMC into larger groups to in-crease computational efficiency. The groups belong to thesame spatial cell and must all have opacities that make thecell sufficiently optically thick; this is an optimization sinceeffective scattering for particles in either of the original groupsis reduced (Densmore et al. 2012). We term this optimization“opacity regrouping.” Opacity regrouping was first impliedby Densmore et al. (2012) with a low-frequency DDMC groupadaptively adding or subtracting adjacent IMC groups basedon the mean free path threshold. Moreover, the extensionof the optimization to strongly non-monotonic opacity wasanticipated by Densmore et al. (2012). Recently, an opacityregrouping procedure for non-contiguous groups was imple-mented by Cleveland & Gentile (2014) for Hybrid ImplicitMonte Carlo Diffusion (HIMCD); in addition to improvingcode performance, their approach addresses the effects of tele-portation error (Fleck & Canfield 1984) with new methodcoupling criteria. In addition to the IMC-DDMC mean freepath threshold, τ D , we introduce an additional mean free paththreshold, τ L , that determines regroupable DDMC groups.We investigate the effect of changing regrouping parameterson a simple ten-group problem and the one-dimensional W7problem presented by Nomoto et al. (1984). Additionally,we explore the effect of a modified Fleck factor, presentedby Gentile (2011), on mitigating erroneous fluctuations in thetemperature profile in the W7 test problem.This article is organized as follows. In Section 2, we discussthe approximations to the radiation transport and fluid equa-tions assumed in our code. In Section 3, we perform an asymp-totic analysis which indicates a potential source of discrepancybetween full multigroup IMC with a discretized Doppler shiftcorrection and continuous-frequency IMC in a multigroup ma-terial setting. In Section 4, we describe the Gentile-Fleck factorused in some numerical results and we summarize the IMC-DDMC equations. Additionally, we write the equations foropacity regrouping. In Section 5, we write the formulae usedto regroup subsets of groups. In Section 6, we describe IMC-DDMC particle processes including the opacity regroupingand DDMC redshift schemes. In Section 7, we present somecalculations that highlight the advantages of the Gentile-Fleckfactor and opacity regrouping and demonstrate the applicationof SuperNu to SNe Ia. In Section 7.1, combining the tech-niques of Oberkampf & Roy (2010) and Gentile (2011), weuse a simple quasi-manufactured transport solution for high-velocity outflow to verify the Gentile-Fleck factor’s ability tomitigate spurious overheating. In Section 7.2, we demonstratethe improved performance that using DDMC opacity regroup-ing produces for the multigroup outflow problems presentedby Wollaeger et al. (2013). Finally, in Section 7.3, we explore the application of IMC-DDMC with opacity regrouping andthe Gentile-Fleck factor to the W7 problem. We also inves-tigate the effects of group opacities that are a composite ofRosseland-like and Planck-like opacities. RADIATION AND FLUID EQUATIONSWe review the underlying theory of the IMC-DDMC schemetested. Following Pomraning (1973) and Castor (2004), termsin the comoving fluid frame are subscripted with 0. The ther-mal equation of radiation transport in the lab frame is (Sz˝oke& Brooks 2005; Abdikamalov et al. 2012) c ∂I ν ∂t + ˆΩ · ∇ I ν + σ ν,a I ν = σ ν,a B ν − σ ν,s I ν + (cid:90) π (cid:90) ∞ νν (cid:48) σ s ( (cid:126)r, ν (cid:48) → ν, ˆΩ (cid:48) → ˆΩ) I ν (cid:48) ( (cid:126)r, ˆΩ (cid:48) , t ) dν (cid:48) d Ω (cid:48) (1)where c is the speed of light, t is time, (cid:126)r is the spatial coor-dinate, ˆΩ is unit direction, ν is frequency, σ a,ν is absorptionopacity, σ s,ν is scattering opacity, σ s ( (cid:126)r, ν (cid:48) → ν, ˆΩ (cid:48) → ˆΩ) isdifferential scattering opacity, I ν is the radiation intensity, and B ν is the thermal emission source. The first order comovingform of Eq. (1) is (Castor 2004, p. 111) (cid:32) · (cid:126)Uc (cid:33) c DI ,ν Dt + ˆΩ ·∇ I ,ν − ν c ˆΩ ·∇ (cid:126)U ·∇ ν ˆΩ I ,ν + 3 c ˆΩ · ∇ (cid:126)U · ˆΩ I ,ν = σ ,ν ,a ( B ,ν − I ,ν ) − σ ,ν ,s I ,ν + (cid:90) π (cid:90) ∞ ν ν (cid:48) σ ,s ( (cid:126)r, ν (cid:48) → ν , ˆΩ (cid:48) · ˆΩ ) I ,ν (cid:48) dν (cid:48) d Ω (cid:48) , (2)where (cid:126)r is an Eulerian spatial coordinate, (cid:126)U is the velocityfield, and we have used Castor’s notation to denote the photoncomoving momentum derivative with ∇ ν ˆΩ . The homologousflow equation is (Kasen et al. 2006) (cid:126)r = (cid:126)U t , (3)Equation (3) allows for some simplification to material- radia-tion coupling. The Lagrangian momentum and energy equa-tions, respectively, are ρ D (cid:126)UDt + ∇ P = − (cid:126)g , (4)and C v DTDt + P ∇ · (cid:126)U = − g (0) , (5)where ρ is density, P is fluid pressure, T is fluid temperature, C v is heat capacity per unit volume, and ( g (0) , (cid:126)g ) is a radiationenergy-momentum coupling 4-vector. Following the justifica-tion provided by Kasen et al. (2006) and van Rossum (2012),we neglect P . For the time scales and physical specificationsof interest, much more energy is in the radiation field than thematerial. Incorporating Eq. (3) and P = 0 into Eqs. (4) and (5)yields C v DTDt = (cid:90) π (cid:90) ∞ σ ,ν ,a ( I ,ν − B ,ν ) dν d Ω + (cid:90) π (cid:90) ∞ σ ,ν ,s I ,ν dν d Ω − (cid:90) π (cid:90) ∞ (cid:90) π (cid:90) ∞ ν ν (cid:48) σ ,s ( (cid:126)r, ν (cid:48) → ν , ˆΩ (cid:48) · ˆΩ ) I ,ν (cid:48) dν (cid:48) d Ω (cid:48) dν d Ω = − g (0)0 ,a − g (0)0 ,s , (6)where g (0)0 ,a and g (0)0 ,s are absorption and scattering contribu-tions to the comoving radiation-material coupling, respectively.Equation (6) is similar in form to the material equation pre-sented by Sz˝oke & Brooks (2005) but with a Lagrangian tem-poral derivative. DOPPLER SHIFT GROUP EDGE ANALYSISMonte Carlo particles may be tracked by either discretegroups or continuous values in frequency space. In the contextof relativistic velocity, Doppler shift has an important effecton the radiation intensity’s interaction with a group structure.When considering how to track particles through phase space,it is informative to consider approaches to sustaining consis-tency between multigroup transport and multigroup diffusion.Specifically IMC may have particle frequency tracked andupdated continuously in a multigroup setting through explicitchanges in reference frame. In contrast, a DDMC particlewavelength is essentially unknown within a group since aDDMC particle step in theory replaces multiple correspondingIMC collision steps. Hence, each time a continuous frequencyvalue is needed from a DDMC particle, it must be sampledfrom a subgroup distribution (Densmore et al. 2012). DDMCparticles may be tracked with continuous frequencies or wave-lengths but the values then merely serve as a label for the sur-rounding group. Consequently, multigroup IMC may simulatethe frequency derivative in Eq. (2) exactly while the DDMCscheme described by Wollaeger et al. (2013) can not exactlysimulate the frequency derivative. We perform an asymptoticanalysis for frequency-dependent, semi-relativistic, comovingtransport with the simplification of homologous outflow be-fore considering a group grid that is constant in the comovingframe along with the upwind redshift approximation (Mihalas& Mihalas 1984, p. 475). A group edge of an optically thickregion of frequency is treated in a manner analogous to spatialboundary layers (Habetler & Matkowsky 1975; Malvagi &Pomraning 1991). Incorporating Eq. (3) in Eq. (2), c ∂I ,ν ∂t + ˆΩ · ∇ I ,ν + σ ,ν I ,ν − ν ct ∂I ,ν ∂ν + (cid:126)rct · ∇ I ,ν + 3 ct I ,ν = j ,ν , (7)where σ ,ν = σ ,ν ,a + σ ,ν ,s is isotropic, j ,ν is the totalsource due to scattering and external sources, and the ˆΩ · (cid:126)U /c term multiplying the Lagrangian derivative has been neglected.Following prior authors (Habetler & Matkowsky 1975; Mal-vagi & Pomraning 1991), we introduce a parameter, ε (cid:28) ,and make the following scalings: c → c/ε , σ ,ν → σ ,ν /ε , σ ,ν ,a → εσ ,ν ,a , ω = ( ν − ν b ) /ε m , q → εq , where ν b is a frequency at boundary b in frequency space and q is theexternal or thermal source in j ,ν . The value m is a numberintroduced to control the amount of variation in intensity with respect to frequency. If ∂I ,ν /∂ω is O(1), then ∂I ,ν /∂ν isO( /ε m ). Incorporating the scalings into Eq. (7), ε c ∂I ,ν ∂t + ε ˆΩ · ∇ I ,ν + σ ,ν I ,ν − ε − m ct ν ∂I ,ν ∂ω + ε ct (cid:126)r · ∇ I ,ν + 3 ε ct I ,ν = εj ,ν , (8)and assuming isotropic elastic scattering, εj ,ν = ε q π + (cid:0) σ ,ν − ε σ ,ν ,a (cid:1) π (cid:90) π I ,ν d Ω (cid:48) . (9)For our purposes, we need only consider m ∈ { , } for aninterior group solution ( m = 0 ) and a frequency boundarylayer solution ( m = 1 ). The intensity may then be decomposedas I ,ν = I i + I b (Malvagi & Pomraning 1991) where I i isthe interior frequency solution and I b is the boundary layerfrequency solution. Moreover, all solutions may be expandedas a power series in ε , I ( i,b ) = (cid:80) ∞ k =0 I ( k )( i,b ) ε k . Additionally,we constrain lim ω →∞ I b = 0 ; this constraint is analogous tothe spatial boundary layer constraint of Malvagi & Pomraning(1991) where the value ω would instead correspond to distanceaway from a surface along a normal vector.To ensure validity of the stated scalings, we demonstrate theresulting interior solution is the diffusion approximation to thesemi-relativistic moment equations presented by Castor (2004,p. 113). The interior intensity is subsequently used alongwith the boundary layer to obtain the desired result. Setting m = 0 and incorporating the power series in ε , Eq. (8) may beseparated into O( ε ), O( ε ), and O( ε ) equations: I (0) i = φ (0) i π (10)for O( ε ), I (1) i = φ (1) i π − π ˆΩ σ ,ν · ∇ φ (0) i (11)for O( ε ), and I (2) i = 14 π (cid:20) qσ ,ν + φ (2) i − σ ,ν ,a σ ,ν φ (0) i + ν ctσ ,ν ∂φ (0) i ∂ν − (cid:126)rctσ ,ν · ∇ φ (0) i − ctσ ,ν φ (0) i − cσ ,ν ∂φ (0) i ∂t − ˆΩ σ ,ν · ∇ (cid:32) φ (1) i − ˆΩ σ ,ν · ∇ φ (0) i (cid:33)(cid:35) (12)for O( ε ), where Eq. (10) has been used in Eq. (11) andEqs. (10) and (11) have been used in Eq. (12). The values φ ( k ) i = (cid:82) π I ( m ) i d Ω are the ε power series coefficients forscalar intensity. Integrating Eq. (12) over comoving solid an-gle, c ∂φ (0) i ∂t − ∇ · (cid:18) σ ,ν ∇ φ (0) i (cid:19) + σ ,ν ,a φ (0) i − ν ct ∂φ (0) i ∂ν + (cid:126)rct · ∇ φ (0) i + 3 ct φ (0) i = q . (13)With some manipulation (by reverting (cid:126)r/t to (cid:126)U and /t to ∇ · (cid:126)U ), Eq. (13) can be seen to be the diffusion approximationto the zeroth-moment, frequency-dependent transport equationpresented by Castor (2004, p. 113) under the assumptions ofisotropic, elastic scattering in the comoving frame and homol-ogous flow.Next we set m = 1 and asymptotically analyze the frequencyboundary. In the domain examined, the optically thick regionwill be at higher frequency, or ω > . Applying the ε powerseries again, the O( ε ), O( ε ), and O( ε ) equations for I b are I (0) b = φ (0) b π (14)for O( ε ), I (1) b = 14 π (cid:32) φ (1) b − ˆΩ σ ,ν · ∇ φ (0) b + ν b ctσ ,ν ∂φ (0) b ∂ω (cid:33) (15)for O( ε ), and I (2) b = 14 π (cid:20) φ (2) b − σ ,ν ,a σ ,ν φ (0) b + ν b ctσ ,ν ∂φ (1) b ∂ω − (cid:126)rctσ ,ν · ∇ φ (0) b − ctσ ,ν φ (0) b − cσ ,ν ∂φ (0) b ∂t − ˆΩ σ ,ν · ∇ (cid:32) φ (1) b − ˆΩ σ ,ν · ∇ φ (0) b (cid:33)(cid:35) (16)for O( ε ), where φ ( k ) b = (cid:82) π I ( m ) b d Ω . The term ∂φ (0) b /∂ω =0 from integration of Eq. (16); this is an important result forthe remainder of the derivation and has been used in Eqs (15)and (16). If Eq. (16) is integrated, closure for φ (0) b is notobtained. In particular, ∂φ (1) b /∂ω persists. The O( ε ) solutionin terms of I (1 , , b is c ∂I (1) b ∂t + ˆΩ · ∇ I (2) b + σ ,ν I (3) b − ν b ct ∂I (2) b ∂ω + (cid:126)rct · ∇ I (1) b + 3 ct I (1) b = σ ,ν π φ (3) b − σ ,ν ,a π φ (1) b . (17)To obtain an equation for φ (1) b , Eq. (16) may be incorporatedinto the second and fourth terms on the left hand side ofEq. (17) and the overall result may be integrated in Ω . Uponintegration of ˆΩ · ∇ I (2) b , values in Eq. (16) that are even in ˆΩ vanish. Upon integration of ∂I (2) b /∂ω , values in Eq. (16) thatare odd in ˆΩ vanish. Fortunately, any terms with ∂φ (0) b /∂ω vanish as well. The result is c ∂φ (1) b ∂t − ∇ · (cid:18) σ ,ν ∇ φ (1) b (cid:19) + σ ,ν ,a φ (1) b − ν b ct ∂∂ω (cid:32) ν b ctσ ,ν ∂φ (1) b ∂ω (cid:33) + 3 ct φ (1) b − ν b ct ∂φ (2) b ∂ω + (cid:126)rct · ∇ φ (1) b = 0 . (18) The first and fourth terms in Eq. (18) together resemble adiffusion equation in frequency space. The system of equationsis still not closed, but Eq. (17) along with Eq. (18) imply ∂∂ω (cid:32) ν b ctσ ,ν ∂φ (1) b ∂ω (cid:33) = 0 . (19)Taking σ ,ν = σ ,ν b , Eq. (19) solves to φ (1) b = ctσ ,ν ν b A ω + A , (20)where A and A are constant in ω . But lim ω →∞ φ (1) b = 0 ,so φ (1) b = A = A = 0 . With ∂φ (1) b /∂ω = 0 , integration ofEq. (16) yields c ∂φ (0) b ∂t − ∇ · (cid:18) σ ,ν ∇ φ (0) b (cid:19) + σ ,ν ,a φ (0) b + (cid:126)rct · ∇ φ (0) b + 3 ct φ (0) b = 0 . (21)Equation (21) indicates the leading-order boundary layer so-lution has no Doppler correction term when ∂I b /∂ν variesstrongly (or m = 1 ). Summing Eqs. (13) and (21), c ∂φ (0)0 ,ν ∂t − ∇ · (cid:18) σ ,ν ∇ φ (0)0 ,ν (cid:19) + σ ,ν ,a φ (0)0 ,ν − ν ct ∂φ (0) i ∂ν + (cid:126)rct · ∇ φ (0)0 ,ν + 3 ct φ (0)0 ,ν = q . (22)where φ ,ν = φ (0) i + φ (0) b is the uniformly valid leading-ordersolution. If the interior solution of the upper frequency rangeis constant in frequency, then c ∂φ (0)0 ,ν ∂t − ∇ · (cid:18) σ ,ν ∇ φ (0)0 ,ν (cid:19) + σ ,ν ,a φ (0)0 ,ν + (cid:126)rct · ∇ φ (0)0 ,ν + 3 ct φ (0)0 ,ν = q . (23)The Doppler correction is removed from the leading-orderscalar intensity equation in the range of frequencies ν > ν b when the leading-order interior solution is constant in fre-quency. In a piecewise-constant multigroup setting with high-contrast opacities, the intensity can vary significantly betweengroups and might be treated as constant within groups. In-tegration of Eq. (23) over a group interval does not producecoupling between groups.We now extend the analysis to problems with an inelasticscattering component. The extension is a model that servesto provide theoretical evidence that group discretization mayhave a nontrivial effect on problems with real or effective in-elastic scattering (such as those solved with IMC). Densmore(2011) asymptotically analyzes the effect of treating some ab-sorption and re-emission as instantaneous effective scatteringwhile treating the remainder explicitly with a linear spatialsampling distribution. We draw an analogy here between elas-tic scattering, which preserves ν , and IMC effective scattering,which preserves (cid:126)r . To complete the analogy, inelastic scatter-ing redistributes ν while IMC effective absorption/emissionredistributes (cid:126)r . We now generalize Eq. (7) to include a ped-agogical model of inelastic scattering in the diffusive upperfrequency range. This inelastic scattering component is meantto emulate effective scattering in IMC within one group. Werewrite Eq. (7) as c ∂I ,ν ∂t + ˆΩ · ∇ I ,ν + σ ,ν I ,ν − ν ct ∂I ,ν ∂ν + (cid:126)rct · ∇ I ,ν + 3 ct I ,ν = q π +14 π (1 − χ ) σ s φ ,ν + 14 π χσ s p s ( ν ) φ ,g (24)where χ ∈ [0 , is a elastic/inelastic splitting parame-ter, p s ( ν ) is a probability density function, σ s is a fre-quency independent scattering opacity coefficient, and φ ,g = (cid:82) ν t ν b φ ,ν dν . The value ν t is the upper bound of the diffu-sive region. Constraining (cid:82) ν t ν b p s ( ν ) dν = 1 , the integralof the total scattering source term over frequency is σ s φ ,g .Considering Eq. (24) implies (cid:90) π (cid:90) ν t ν b ν ν (cid:48) σ ,s ( (cid:126)r, ν (cid:48) → ν , ˆΩ (cid:48) · ˆΩ ) I ,ν (cid:48) dν (cid:48) d Ω (cid:48) =14 π (1 − χ ) σ s φ ,ν + 14 π χσ s p s ( ν ) φ ,g , (25)a consistent differential scattering opacity is σ ,s ( (cid:126)r, ν (cid:48) → ν , ˆΩ (cid:48) · ˆΩ ) = σ s π (cid:20) (1 − χ ) δ ( ν − ν (cid:48) ) + χ ν (cid:48) ν p s ( ν ) (cid:21) , (26)where δ ( ν − ν (cid:48) ) is the Dirac distribution. Thus the totalscattering opacity is σ ,ν ,s = σ s (cid:20) (1 − χ ) + χν (cid:90) ν t ν b p s ( ν (cid:48) ) ν (cid:48) dν (cid:48) (cid:21) . (27)Furthermore, we define a secondary distribution, ˜ p s ( ν ) = (cid:18)(cid:90) ν t ν b p s ( ν (cid:48) ) ν (cid:48) dν (cid:48) (cid:19) − p s ( ν ) ν , (28)which is shown below to be the O( ε ) and O( ε ) frequencydependence of scalar intensity. We define φ i,g and φ b,g asthe interior and boundary scalar intensity group integratedcontributions to to the diffusive range. Applying the scalingswith m = 0 , considering the interior solution, and setting φ ,g = (cid:80) ∞ k =0 φ ( k )0 ,g ε k = (cid:80) ∞ k =0 ( φ ( k ) i,g + φ ( k ) b,g ) ε k , the O( ε ),O( ε ), and O( ε ) equations for intensity are I (0) i = 14 π ˜ p s ( ν ) φ (0) i,g , (29) I (1) i = 14 π ˜ p s ( ν ) (cid:32) φ (1) i,g − ˆΩ σ ,ν ,s · ∇ φ (0) i,g (cid:33) , (30) and I (2) i = 14 π (cid:20) qσ ,ν ,s + σ s σ ,ν ,s [(1 − χ ) φ (2) i + χp s ( ν ) φ (2) i,g ] − σ ,ν ,a σ ,ν ,s ˜ p s ( ν ) φ (0) i,g + ν φ (0) i,g ctσ ,ν ,s ∂ ˜ p s ∂ν − ˜ p s ( ν ) ctσ ,ν ,s (cid:126)r · ∇ φ (0) i,g − p s ( ν ) ctσ ,ν ,s φ (0) i,g − cσ ,ν ,s ∂ (˜ p s ( ν ) φ (0) i,g ) ∂t − ˜ p s ( ν ) σ ,ν ,s ˆΩ · ∇ (cid:32) φ (1) i,g − ˆΩ σ ,ν ,s · ∇ φ (0) i,g (cid:33)(cid:35) (31)respectively. Integration of Eq. (31) gives a correct form of thecomoving diffusion equation. Additionally, Eq. (31) indicatesthe Doppler coupling in the diffusion region is dependent onthe inelastic scattering profile. The scattering profile deter-mines the leading interior solution. For m = 1 , the O( ε ) andO( ε ) equations are I (0) b = 14 π ˜ p s ( ν b ) φ (0) b,g , (32) I (1) b = 14 π (cid:18) σ s σ ,ν ,s [(1 − χ ) φ (1) b + χp s ( ν b ) φ (1) b,g ] − ˆΩ σ ,ν ,s · ∇ φ (0) b + ν b ctσ ,ν ,s ∂φ (0) b ∂ω (cid:33) , (33)respectively, where it is assumed the inelastic probability den-sity does not vary strongly in the boundary layer. This as-sumption may be more clearly expressed as an Taylor ex-pansion of p s ( ν ) around ν b at a point in the boundary layer: p s ( ν ) = p s ( ν b ) + εω∂p s ( ν b ) /∂ν . Equation (32) is frequencyindependent; so ∂φ (0) b /∂ω = 0 . Integration of Eq. (33) oversolid angle yields φ (1) b = ˜ p s ( ν b ) φ (1) b,g . (34)Equation (34) implies ∂φ (1) b /∂ω = 0 . Invocation of ∂φ (2) b /∂ω equation was not needed to obtain Eq. (34). The O( ε ) bound-ary layer equation is I (2) b = 14 π (cid:20) σ s σ ,ν ,s ((1 − χ ) φ (2) b + χp s ( ν b ) φ (2) b,g ) − σ ,ν ,a σ ,ν ,s φ (0) b − (cid:126)rctσ ,ν ,s · ∇ φ (0) b − ctσ ,ν ,s φ (0) b − cσ ,ν ,s ∂φ (0) b ∂t − ˆΩ σ ,ν ,s · ∇ (cid:32) φ (1) b − ˆΩ σ ,ν ,s · ∇ φ (0) b (cid:33)(cid:35) . (35)Equation (35) gives a diffusion equation, c ∂φ (0) b ∂t − ∇ · (cid:18) σ ,ν ,s ∇ φ (0) b (cid:19) + σ ,ν ,a φ (0) b + (cid:126)rct ·∇ φ (0) b + 3 ct φ (0) b = χσ s (cid:18) p s ( ν b ) φ (2) b,g − p s ( ν b )˜ p s ( ν b ) φ (2) b (cid:19) , (36)which has an inelastic scattering source from the O( ε ) scalarflux. Finally, integrating Eq. (35) over Ω , differentiating theresult with respect to ω , and using ∂φ (0) b /∂ω = ∂φ (1) b /∂ω = 0 yields ∂φ (2) b ∂ω = (1 − χ )(1 − χ ) + χν b (cid:82) ν t ν b p s ( ν (cid:48) ) ν (cid:48) dν (cid:48) ∂φ (2) b ∂ω . (37)If χ = 0 , scattering is entirely elastic and Eq. (37) is self-consistent. Otherwise, Eq. (37) is solved with ∂φ (2) b /∂ω = 0 (this may be seen from differentiation of Eq. (36) with respectto ω as well). The uniformly valid diffusion equation is c ∂φ (0)0 ,ν ∂t − ∇ · (cid:18) σ ,ν ∇ φ (0)0 ,ν (cid:19) + σ ,ν ,a φ (0)0 ,ν − ν φ (0) i,g ct ∂ ˜ p s ∂ν + (cid:126)rct · ∇ φ (0)0 ,ν + 3 ct φ (0)0 ,ν = q + χσ s (cid:18) p s ( ν ) φ (2)0 ,g − p s ( ν )˜ p s ( ν ) φ (2)0 ,ν (cid:19) . (38)where we have made use of σ ,ν = σ ,ν ,s + O( ε ). Pho-ton number density is proportional to φ ,ν /ν . Setting ˜ φ = φ (0)0 ,ν /ν gives an equation for number density in thecomoving frame: c ∂ ˜ φ∂t − ∇ · (cid:18) σ ,ν ∇ ˜ φ (cid:19) + σ ,ν ,a ˜ φ − φ (0) i,g ct ∂ ˜ p s ∂ν + (cid:126)rct · ∇ ˜ φ + 3 ct ˜ φ = qν + χσ s (cid:18) p s ( ν ) ν φ (2)0 ,g − (cid:90) ν t ν b p s ( ν (cid:48) ) ν (cid:48) dν (cid:48) φ (2)0 ,ν (cid:19) . (39)Integration of Eq. (39) causes the inelastic scattering termon the right hand side to vanish. Consequently, the Dopplercorrection is again dependent on the interior solution but nowalso on the scattering distribution, ˜ p s . If ˜ p s = 1 / ( ν t − ν b ) ,then the comoving photon number density diffusion equationhas no Doppler correction term.The boundary layer solutions do not provide Doppler cor-rections in the sense described by Castor (2004, p. 112). Wethus focus on the Doppler correction that the interior solutionprovides at the group boundary. Additionally, sufficient inelas-ticity in collisions, or χ ∼ O(1) in Eq. (24), makes the Dopplercorrection dependent on the redistribution profile.To obtain the upwind approximation for Doppler shift in allgroups, the transport equation may first be group integrated.We define a frequency grid in the comoving frame with G groups: ν G +1 / < . . . < ν / . Integrating Eq. (7) over acomoving group, g , yields c ∂I ,g ∂t + ˆΩ · ∇ I ,g + σ ,g I ,g + 4 ct I ,g − ct ( ν g − / I ,ν g − / − ν g +1 / I ,ν g +1 / )+ (cid:126)rct ·∇ I ,g = j ,g , (40)where I ,g = (cid:82) ν g − / ν g +1 / I ,ν dν , σ ,g = (cid:82) ν g − / ν g +1 / σ ,ν I ,ν dν / (cid:82) ν g − / ν g +1 / I ,ν dν , and j ,g = (cid:82) ν g − / ν g +1 / j ,ν dν . In practice, σ ,g , might be computedwith an approximation since the exact value is dependenton the solution. Alternatively, one could define the opacityas piecewise constant in frequency. Applying the upwindapproximation to the edge frequency-dependent intensityterms yields (Mihalas & Mihalas 1984, p. 475) c ∂I ,g ∂t + ˆΩ · ∇ I ,g + σ ,g I ,g + 4 ct I ,g + ν g +1 / ct ∆ ν g I ,g + (cid:126)rct · ∇ I ,g = j ,g + ν g − / ct I ,g − ∆ ν g − , (41)where ∆ ν g = ν g − / − ν g +1 / . The upwind approximationmay be extended trivially to find the multigroup form of Eq. (2).The fifth term on the left hand side and the second term onthe right hand side of Eq. (41) are responsible for couplinggroups through Doppler shifting. If the group coupling terms inEq. (41) are removed, then the result describes grey multigrouptransport in the context of homologous outflow. If Eq. (41) issolved with a grey MC transport scheme that includes expan-sion effects (through frame transformations and spatial gridexpansion), then a stochastic interpretation must be given tothe Doppler shift group coupling terms. The diffusion equationcorresponding to Eq. (41) may be found by integrating Eq. (41)over comoving angle and applying Fick’s Law, c ∂φ ,g ∂t − ∇ · (cid:18) σ ,g ∇ φ ,g (cid:19) + σ ,g φ ,g + 4 ct φ ,g + ν g +1 / ct ∆ ν g φ ,g + (cid:126)rct · ∇ φ ,g = 4 πj ,g + ν g − / ct ∆ ν g − φ ,g − , (42)where opacities have been assumed piecewise constant in fre-quency. The Doppler correction terms in Eqs. (41) and (42)can be interpreted as “Doppler shift opacities”, where samplingthe value ν g +1 / /ct ∆ ν g would induce a particle to transitionfrom group g to group g + 1 . If an IMC particle samples aDoppler shift event, the particle’s frequency will be updated toan adjacent group.Instead of assuming a fully grouped approach, we imple-ment a Doppler shift scheme in IMC-DDMC that more closelyemulates continuous frequency transport in the presence ofpiecewise constant opacities. We make the constraint in ourcode that inelastic redistribution at the subgroup level is uni-form, or p s ( ν ) = 1∆ ν g . (43)Considering Eqs. (28), and (29): ˜ p s ∼ /ν ∼ φ (0) i , and theDoppler correction in Eq. (38) and (39) satisfies − ν φ (0) i,g ct ∂ ˜ p s ∂ν = 1 ct φ (0) i . (44)Since the equations for scalar flux in the frequency boundarylayer have no Doppler correction, we assume I b = 0 ; theinterior radiation field thus account for all radiation in the dif-fusive frequency region. Then the entire radiation field has theDoppler correction. Consequently, incorporating Eq. (44) intoEq. (38), neglecting higher order scattering terms, assumingpiecewise constant opacities and integrating over the grouprange yields c ∂φ (0)0 ,g ∂t − ∇ · (cid:18) σ ,g ∇ φ (0)0 ,g (cid:19) + σ a,g φ (0)0 ,g + (cid:126)rct · ∇ φ (0)0 ,g + 4 ct φ (0)0 ,g = q g . (45)Equation (45) is Eq. (42) without upwind Doppler shift terms.We infer that the degree of elasticity (in our model χ ) is impor-tant to how DDMC groups redshift to other groups, particu-larly when DDMC emulates continuous frequency transport.In order to have Eq. (45) represent grey diffusion for the caseof one group, we limit Doppler shift of particles to adjacentgroups for problems with inelastic-dominant collisions, or χ ∼ O(1). Such a constraint should emulate IMC for problems withinelastic-dominant collisions. Assuming a non-zero veloc-ity field exists and inelastic opacity is large with respect to ν g +1 / /ct ∆ ν g , IMC particles would have their frequenciesredistributed many times before streaming to the edge of agroup; this may greatly reduce the occurrence of Doppler shiftbetween groups in IMC. In Section 6, we describe a DDMCDoppler shift scheme that takes into account the degree ofinelasticity in collisions. MULTIGROUP IMC-DDMC EQUATIONSEquation (6) is amenable to the semi-implicit time differencedescribed by Fleck & Cummings (1971). Moreover, the semi-implicit discretization procedure may be applied on Eqs. (2)and (6) to obtain IMC equations for the comoving frame. Themultigroup form of Eq. (6) is C v DTDt = G (cid:88) g =1 (cid:90) π σ a,g I ,g d Ω − cσ P aT − g (0)0 ,s (46)where σ a,g is comoving grouped absorption opacity, σ P iscomoving Planck opacity, and we have compressed the nota-tion of the inelastic scattering contribution since it is a mate-rial source with a treatment described by Fleck & Cummings(1971). Introducing a parameter β = 4 aT /C v and integratingEq. (46) over a time step gives ( aT ) n +1 − ( aT ) n = (cid:90) t n +1 t n β (cid:32) G (cid:88) g =1 (cid:90) π σ a,g I ,g d Ω − cσ P aT − g (0)0 ,s (cid:33) dt , (47)where a value subscripted with n implies evaluation at thebeginning of a time step indexed by n . IMC is made semi-implicit and linear within a time step by setting β = β n , σ a,g = σ a,g,n , and σ P = σ P,n (Fleck & Cummings 1971;Fleck & Canfield 1984). Additionally, setting ∆ t n ¯ I ,g = (cid:82) t n +1 t n I ,g dt , ∆ t n [ αT n +1 + (1 − α ) T n ] = (cid:82) t n +1 t n T , and ∆ t n ¯ g (0)0 ,s = (cid:82) t n +1 t n g (0)0 ,s dt gives aT n +1 − aT n = β n ∆ t n G (cid:88) g =1 (cid:90) π σ a,g,n ¯ I ,g d Ω − c ∆ t n β n σ P,n [ αaT n +1 + (1 − α ) aT n ] − ∆ t n β n ¯ g (0)0 ,s (48) where ∆ t n = t n +1 − t n and α ∈ [0 , is the standard IMCtime centering parameter. With Eq. (48), an expression maybe found for αaT n +1 + (1 − α ) aT n that excludes T n +1 . In-troducing the Fleck factor, f n = 11 + αβ n c ∆ t n σ P,n , (49)the time centered aT is (Abdikamalov et al. 2012) αaT n +1 + (1 − α ) aT n =1 cσ P,n (1 − f n ) G (cid:88) g =1 (cid:90) π σ a,g,n ¯ I ,g d Ω + f n aT n − cσ P,n (1 − f n )¯ g (0)0 ,s . (50)By replacing ¯ I ,g with I ,g , the thermal emission source termfor a group g in the comoving transport equation may be ap-proximated as σ a,g,n B ,g = 14 π caT σ a,g,n b g,n =(1 − f n ) σ a,g,n b g πσ P G (cid:88) g (cid:48) =1 (cid:90) π σ a,g (cid:48) ,n ¯ I ,g (cid:48) d Ω + σ a,g,n b g,n πσ P,n f n acT n − (1 − f n ) σ a,g,n b g,n πσ P,n ¯ g (0)0 ,s . (51)Equations (47)-(51) are not the only way to semi-implicitlydiscretize the temperature equation in time. Moreover, in cer-tain circumstances it may be appropriate to apply differentapproximations in order to avoid problematic IMC errors. Inparticular, Larsen & Mercier (1987) derive a “Maximum Prin-ciple” for IMC that supplies a sufficient but not necessaryupper bound on time step sizes. It follows from their analysisthat IMC is not guaranteed to give a physical result for anypossible numerical setup. If IMC numerical parameters are ill-conditioned, spurious temperature oscillations and overheatingmay occur (McClarren & Urbatsch 2012). Gentile (2011) per-forms a similar discretization but linearly expands opacity and aT from their values at n to values at n + 1 . Despite severeapproximations (Gentile 2011), the result is a modified Fleckfactor that adapts to the state of the radiation field. Instead ofexpanding material quantities in T , an alternative approachto obtaining the result of Gentile (2011) is to make a changeof variables in the time derivative similar to that of Fleck &Cummings (1971). Defining E ∗ = 1 c ∆ t n ¯ σ P (cid:90) t n +1 t n G (cid:88) g =1 (cid:90) π σ a,g I ,g d Ω dt , (52)where ¯ σ P is time centered, Equation (46) may be stated as σ P ˜ β DDt [ σ P ( aT − E ∗ )]= G (cid:88) g =1 (cid:90) π σ a,g I ,g d Ω − cσ P aT − g (0)0 ,s , (53)where ˜ β = 1 C v (cid:20) aT + ( aT − E ∗ ) 1 σ P ∂σ P ∂T (cid:21) . (54)Evaluating σ P ˜ β on the left hand side of Eq. (53) at the be-ginning of a time step, integrating Eq. (53) with respect totime, setting (cid:82) t n +1 t n σ P aT = ∆ t n [ ασ P,n +1 aT n +1 + (1 − α ) σ P,n aT n ] , setting ¯ σ P = ασ P,n +1 + (1 − α ) σ P,n , and set-ting Λ a,n = σ P,n ( aT n − E ∗ ) give Λ a,n +1 − Λ a,n = c ∆ t n σ P,n ˜ β n (cid:16) − α Λ a,n +1 − (1 − α )Λ a,n − ¯ g (0)0 ,s (cid:17) . (55)Defining the Gentile-Fleck factor as ˜ f n = 11 + α ˜ β n c ∆ t σ P,n , (56)The time centered emission term is found to be ασ P,n +1 aT n +1 + (1 − α ) σ P,n aT n =˜ f n σ P,n aT n − (1 − ˜ f n )¯ g (0)0 ,s + ¯ σ P (cid:18) − σ P,n ¯ σ P ˜ f n (cid:19) E ∗ . (57)The next simplification is σ P,n / ¯ σ P in the last term on theright hand side of Eq. (57). By incorporating Eq. (52) for E ∗ , Eq. (52) may be a substitute for the emission term in thecomoving thermal transport equation. The value ˜ f n may beinterpreted in the same manner as f n to control the amount ofeffective scattering and absorption in IMC. Unfortunately, theform of ˜ β n allows ˜ f n to be negative. Gentile (2011) constrains ˜ f n ∈ [0 , by setting ˜ β n =1 C v (cid:34) aT n + max (cid:32) ( aT n − E ∗ ) 1 σ P,n ∂σ P ∂T (cid:12)(cid:12)(cid:12)(cid:12) T n , (cid:33)(cid:35) (58)Additionally, E ∗ is estimated with the tallied radiation en-ergy density from time step n − . Equations (56) and (58)are the exact same equations for the modified Fleck factorderived by Gentile (2011). If the Planck opacity decreaseswith temperature and the radiation temperature is higher thanthe material temperature, then ˜ β n > β n and ˜ f n < f n . FromEq. (58), it is evident that ˜ f n ≤ f n and the Gentile-Fleck factoralways increases effective scattering over the standard Fleckfactor (Gentile 2011). Unfortunately, the cost of more stabilityin IMC temperature update is a decrease in IMC efficiency.However, hybridizing IMC with a diffusion scheme mitigatesthe added cost (Gentile 2011).It remains to assess whether or not such a modification toIMC is needed for problems like the W7 SN Ia describedby Nomoto et al. (1984). The grey form of the MaximumPrinciple of Larsen & Mercier (1987) is ∆ t n (cid:20) ac sup T L
Quasi-Manufactured Verification
Our first problem is a test of the Gentile-Fleck factor us-ing a quasi-manufactured solution (Oberkampf & Roy 2010)for grey transport in a high-velocity outflow. Here, a quasi-manufactured radiation transport solution has an assumed, ormanufactured, radiation energy density profile and, in contrast,a material temperature that is solved for using the manufac-tured radiation energy density and the material equation. Themanufactured source term is incorporated in the radiation trans-port equation to counter redshift and preserve the constancy ofthe manufactured radiation energy density. For the numericalregime considered, we obtain a positive definite source that issimple to implement. The quasi-manufactured solution pro-vides a benchmark demonstrating that the Gentile-Fleck factor(or modified Fleck factor) provides better accuracy relativeto the standard Fleck factor. Specifically, the Gentile-Fleckfactor decreases effective absorption, which mitigates potentialviolations of the IMC Maximum Principle (Larsen & Mercier1987).Equation (57) is implemented approximately (Gentile2011) in an optimized form since computing the deriva-tive of opacity with respect to temperature may be com-putationally expensive. We use ∂σ P,j,n /∂T ≈ ( σ P,j,n − ( ρ j,n /ρ j,n − ) σ P,j,n − ) / ( T j,n − T j,n − ) for n ≥ , and ∂σ P,j, /∂T ≈ ( σ P,j ((1 + ε ) T j, ) − σ P,j, ) / ( εT j, ) where ε is a user defined parameter. The source term from the manufac-turing is positive-definite and yields a solution with non-trivialtime dependence. Gentile (2001) provides an analytic solutionto a spatially independent problem that is used as a benchmarkfor modified IMC in static material. The opacity is propor-tional to T − , implying that increasing temperature reducesemission. The manufacturing and outflow are an extension ofthe solution, but we find our analytic result somewhat simplerin form. Assuming pure absorption, integrating the comovingtransport equation (Eq. (2)) over frequency, and assuming nospatial dependence yields ∂E∂t + 4 t E = cσ ( T )( aT − E ) + S m , (77)and C v ∂T∂t = cσ ( T )( E − aT ) , (78)where E is radiation energy density and S m is the manufac-tured source. The heat capacity C v = ρc v and the opacityis σ ( T ) = κρT , (79)where c v and κ are constants. We manufacture the radiationfield as constant and solve Eq. (78) to obtain a transcendentalexpression for temperature and time. The manufactured sourcemay then be found from S m = 4 t E + C v ∂T∂t (80)by adding Eqs. (77) and (78). It is clear from Eq. (80) that amonotonically increasing temperature over all time ensuresa positive definite source. This should be the case when T ( E/a ) / . Fortunately a low initialtemperature and high initial radiation field is the setup thatinduces the overheating pathology in standard IMC. Followingthe approach of Gentile (2011), Eq. (78) may be re-expressedas (cid:18) ( E/a ) TE/a − T − T (cid:19) ∂T∂t = acκc v (81)where conveniently, ρ cancels through division of σ ( T ) by C v .Equation (81) yields (cid:114) Ea ln (cid:32) [ (cid:112) E/a + T ][ (cid:112) E/a − T ][ (cid:112) E/a − T ][ (cid:112) E/a + T ] (cid:33) −
12 ( T − T ) = acκc v ( t − t ) , (82)where t and T are the initial time and material temperature,respectively. For material and radiation properties of interest,Eq. (82), indicates long equilibration time between the fields.Specifically, for an initial radiation temperature of 1.70 × K,an initial material temperature of 1.16 × K, a specific heatcapacity of 9.3 × erg/K/g, and κ = 1 . × cm K /g,the characteristic equilibrium time is on the order of 10 seconds. These numbers are borrowed or adapted from Gentile(2011). If the scope of simulation time is much smaller, itmay safely be assumed that T , T (cid:28) ( E/a ) / = T r forthe numbers given. When the material temperature and initialtemperature are much smaller than the radiation temperature,Eq. (82) may be approximated by T ( t ) = T r (cid:34) acκc v T r ( t − t ) + (cid:18) T T r (cid:19) (cid:35) / . (83)From Eq. (80), the time integrated manufactured radiationsource is approximately ∆ t n S m,n = 4 t n E ∆ t n + C v,n ( T n +1 − T n ) , (84)for small time steps. Equation (84) is positive definite whenEq. (83) is used ( T ( t n ) = T n ).We construct a problem that induces a “temperature flip”pathology in standard IMC or DDMC. In the first time step,standard IMC-DDMC causes over deposition; this results inthe radiation energy density and material temperature respec-tively dropping and increasing abruptly despite the more grad-ual nature of the actual solution. Given the strong inversedependence of opacity on temperature, emission abruptly be-comes low, causing the material temperature to remain too highfor time spans of interest. Gentile (2011) demonstrates thisIMC pathology in the context of static material. Our problemconsists of a homologous outflow over 10 spatial cells with amaximum speed of 10 cm/s. The material temperature is uni-formly initialized to 1.16 × K and the radiation temperatureis initialized to the manufactured value of 1.70 × K. Start-ing from an expansion time of 2 days, we compute the MC re-sults over a 10th of a millisecond, or t ∈ [2 , . × − ] days. We test both 100 and 1000 time steps in the time spangiven. The source, Eq (84), is applied uniformly across the 10spatial cells. The density is uniform over the spatial domainwith a total constant mass of M = 1 × g. Additionally, κ = 1 . × cm K /g, c v = 9 . × erg/K/g.Similar to findings of Gentile (2011), for this test problemit is found that modified pure IMC is very inefficient; the Gentile-Fleck factor increases effective scattering in IMC toa large extent relative to the standard Fleck factor in IMC.Since grey DDMC does not model effective scattering explic-itly , we test the Gentile-Fleck factor in DDMC; this approachis similar to the use of RW by Gentile (2011) to acceleratea test calculation. In Figure 1, analytic material temperatureis calculated with Eq. (83). The MC temperatures are ob-tained by implementing the manufactured source, Eq. (84),with Eq. (83) used to evaluate T n and T n +1 . For the MC re-sults, the average of the temperature profiles are taken overthe 10 spatial cells (temperature change from cell to cell isinsignificant, however). Figure 1a has material and radiationtemperature results of IMC and DDMC with the standard Fleckfactor, and the quasi-manufactured solution versus time. InFig. 1a, both the IMC and DDMC solutions suffer the “tem-perature flip” error, in which material temperature becomesnon-physically higher than radiation temperature in the firsttime step. Figure 1b has material and radiation temperatureresults for DDMC with the modified Fleck factor using 100 (de-noted “Large ∆ t ”) and 1000 time steps. Results demonstratethe “temperature flip” error is avoided for DDMC modifiedwith the Gentile-Fleck factor. Increasing the number of timesteps from 100 to 1000 further improves agreement towards thequasi-manufactured solution. We conclude that the overheat-ing pathology in IMC and DDMC can occur in high-velocityflows and that the Gentile-Fleck factor mitigates the overheat-ing error in high-velocity outflow. However, the ability of theGentile-Fleck factor to correct the error is apparently limited,since in the early time steps the material temperature becomestoo high while the radiation temperature drops too low relativeto the analytic solutions.7.2. Ten Group Outflow Test
With 10 group, spherical Heaviside source, outflow prob-lems described by Wollaeger et al. (2013), we test the effectof opacity regrouping in IMC-DDMC for simple yet highlystructured opacities. Specifically, we demonstrate the utilityof regrouping non-contiguous groups for radiation transportin a high-velocity fluid with astrophysical properties. The ap-proach is described in Section 5 for LTE transport. The formof the opacities is meant to only allow for significant codespeed-up when opacities for non-adjacent frequency intervalscan be regrouped. With opacity regrouping allowed for non-contiguous group intervals, a DDMC particle has a probabilityof being in any group that satisfies the regrouping criteria; thisgeneralization improves speed without significant detrimentto accuracy relative to the non-opacity-regrouped (non-OR)results for the numerical specifications considered.The problems consist of a homologous outflow with a maxi-mum outer speed of U max = 10 cm/s. The time domain ofthe problem is t ∈ [2 , days. The temperature of the domainis uniformly initialized to 1.16 × K. There is a uniformradiation source density of × /t n erg/cm /s for | (cid:126)U | ∈ [0 , . U max ] . The source is uniform in frequency as well. Thetotal mass is set to × g equally divided amongst spatialcells. The heat capacity is C v = 2 × ρ erg/cm /K. Thegroups are spaced logarithmically from . × − cm to . × − cm in wavelength with g = 1 being the lowestwavelength group. The opacity in cm − (with ρ in g/cm ) is σ a,g = (cid:26) . ρ , g = 2 k − . × − m ρ , g = 2 k , (85)where k ∈ { . . . } and m is set to 4 or 7 (Wollaeger et al.4 −4 Time (s) T e m pe r a t u r e ( K ) Analytic TAnalytic T r Standard IMC TStandard IMC T r Standard DDMC TStandard DDMC T r | (a) −4 Time (s) T e m pe r a t u r e ( K ) Analytic TAnalytic T r DDMC T (Large ∆ t)DDMC T r (Large ∆ t) |DDMC TDDMC T r (b) Figure 1.
For the quasi-manufactured problem described in Section 7.1, wecompare standard IMC, standard DDMC, and DDMC with the modified Fleckfactor against analytic solutions. In Fig. 1a: analytic (solid), standard IMC(dashed), and standard DDMC (dash-dotted) material (T) and radiation (T r )temperatures for the 1000 time step case. The IMC and DDMC results agreevery closely but are both wrong. The IMC (dashed light blue) and DDMC(dash-dotted yellow) radiation temperatures are closer to the analytic materialtemperature (solid blue) than the analytic radiation temperature (solid green).Inversely, the IMC (dashed red) and DDMC (dash-dotted purple) materialtemperatures are closer to the analytic radiation temperature than the analyticmaterial temperature. This pathology indicates the standard Fleck factor isinsufficient for the time step sizes used. In Fig. 1b: material (T) and radiation(T r ) temperatures from the analytic solution (solid), DDMC with the modifiedFleck factor and 100 time steps (dashed, “Large ∆ t”), and DDMC with themodified Fleck factor and 1000 time steps (dash-dotted). The modified Fleckfactor prevents the radiation and material temperatures from “flipping” (seeFig. 1a). Moreover, a decrease in time step size causes further correction ofthe MC solutions towards the analytic solution. m , we use 50 uniform spatial cells,128 uniform time steps, 0 initial particles, and 100,000 sourceparticles per time step. For all the IMC-DDMC calculationspresented, τ L = τ D = 3 mean free paths.Considering the m = 4 disparity, Fig. 2 has radiation en-ergy densities and material temperatures for IMC, non-ORIMC-DDMC, and opacity-regrouped IMC-DDMC; opacity regrouping is not apparently a significant detriment to thesesolutions. In Fig. 3, the L error for the spectra (in erg/s) ofnon-OR and opacity-regrouped IMC-DDMC relative to IMCincrease while DDMC is dominant and subsequently decreaseas outer cells transition to IMC. The DDMC approximationfor the lab frame spectral tally becomes steadily less accuraterelative to the IMC tally as the cells become optically thin. Theinfluence of opacity regrouping in the m = 7 case is similar tothat of the m = 4 case. In other words the conclusions fromFigs. 2 and 3 hold for the m = 7 case.We also incorporate a regrouping cutoff index, g c , as anexperimental parameter. For a group g that meets the regroup-ing criteria, only groups in the neighborhood g ± g c with anumber of mean free paths for inelastic collisions greater than τ L may have their properties used to accelerate the diffusionof particles in g . For t ∈ [2 , . , we test solution speed versusthe cutoff group displacement for different regrouping cutoffs, g c ∈ { . . . } . Table 1 has times of IMC-DDMC for each g c value along with the time for IMC. All times presented are forsimulation on one core. Table 1
Run Times for First 64 Time Steps of Heaviside Problem with g c ∈ { . . . } with 1 Core (minutes) Method g c m = 4 m = 7 IMC - 202.23 505.71HMC 0 23.11 45.09HMC 1 19.60 37.62HMC 2 5.80 6.74HMC 3 5.74 6.71HMC 4 5.72 6.31HMC 5 5.79 6.44HMC 6 5.80 6.41HMC 7 5.79 6.47HMC 8 5.83 6.42HMC 9 5.76 6.54HMC 10 5.81 6.64
From Table 1, it is evident that regrouping only adjacentgroups provides no significant speed up in computation due tothe highly non-monotonic structuring of opacity versus group.However, when the regrouping cutoff parameter, g c , is set to 2,there is a significant reduction in computational cost.For the problems considered in this section, opacity regroup-ing in IMC-DDMC is seen to be a large computational ad-vantage without large cost of accuracy to important quantities(spectra and temperatures). For different problems, the controlparameters for opacity regrouping may need to be adjustedto maintain good agreement with IMC. To balance efficiencywith solution accuracy, adaptive regrouping parameters mightbe considered. However, for the calculations in the followingsection, opacity regrouping is constrained to τ D = τ L = 3 with g c set to the number of groups.7.3. W7 Tests
We now turn to the W7 problem described by Nomoto et al.(1984) and solved by several authors [see, e.g. Kasen et al.(2006); Kromer & Sim (2009); van Rossum (2012)]. TheW7 problem consists of simulating radiative transfer in a onedimensional model of Type Ia supernovae. The W7 specifica-tions include density and mass fractions for elements up to Ni5 Velocity (cm/s) R ad i a t i on E ne r g y D en s i t y ( e r g / c m ) IMC (3.5 days)Non−OR HMC (3.5 days) |HMC (3.5 days)IMC (5 days)Non−OR HMC (5 days)HMC (5 days) (a) Velocity (cm/s) M a t e r i a l T e m pe r a t u r e ( K ) IMC (3.5 days)Non−OR HMC (3.5 days) |HMC (3.5 days)IMC (5 days)Non−OR HMC (5 days)HMC (5 days) (b) −10 −8 −6 −4 −2 Wavelength (cm) F l u x ( e r g / s ) IMC (3.5 days)Non−OR HMC (3.5 days) |HMC (3.5 days)IMC (5 days)Non−OR HMC (5 days)HMC (5 days) (c)
Figure 2.
Radiation energy density, material temperature, and grouped spec-tra of IMC (solid), non-opacity-regrouped (non-OR) IMC-DDMC (dashed, g c = 0 ), and opacity-regrouped IMC-DDMC (dot-solid, g c = 10 ) at 3.5 and5 days for the 10 group, outflow problem with a spherical Heaviside sourcedescribed in Section 7.2. Radiation energy density and material temperatureare plotted versus fluid velocity and spectra are plotted at group wavelengthcenters. The opacity is described by Eq. (85) with m = 4 . In Figs. 2a, 2b,and 2c, radiation energy density, material temperature, group spectra, respec-tively. For this problem, the IMC-DDMC results with opacity regroupingshow good agreement with the non-OR results. L S pe c t r u m E rr o r Non−OR HMC |HMC
Figure 3.
Non-opacity-regrouped (non-OR) IMC-DDMC (solid), and opacityregrouped IMC-DDMC (dashed) L error versus time step of group spectrarelative to pure IMC for radiation escaping the outermost cell of the 10 group,Heaviside source problem described in Section 7.2. The opacity is describedby Eq. (85) with m = 4 . The spectral error for both non-OR and opacity-regrouped IMC-DDMC progressively increases relative to IMC until pureIMC is applied in the outer cells. Error from opacity regrouping appearsinsignificant relative to error from DDMC. on a velocity grid. The radial outflow speed at the outer bound-ary is ∼
7% of the speed of light. In the free-expansion phaseof the supernova radioactive decay of Ni heats the fluid andcauses it to radiate in the UV, visible and infrared ranges ofthe spectrum. For this problem, we apply the modified Fleckfactor, tested in Section 7.1, and opacity regrouping, tested inSection 7.2. Additionally, we test different calculations of thegrouped opacity by introducing uniform subgroups for eachgroup. Despite the physical and algorithmic complexities ofthe opacity, IMC-DDMC yields light curves and spectra thatare in good agreement with those of
PHOENIX for the numeri-cal specifications considered. Moreover, the total computationtimes are on the order of hours (see Table 2).For IMC-DDMC, a method that in our formulation requiresa group structure, the W7 problem has the difficulty of requir-ing many groups for accurate spectra. Specifically, we findthat the number of groups required to achieve a resolved lightcurve is on the order of thousands. While IMC-DDMC iseasily extensible to 2 and 3 spatial dimensions in theory, stor-ing ∼ SuperNu is run on 192 coreson the Cray XE6 supercomputer Beagle at the ComputationInstitute of the University of Chicago.6In each time step, the opacity per group is computed usinga subgroup structure to allow for non-trivial opacity profileweighting. Opacity contributions to each group include bound-bound (bb), bound-free (bf), and free-free (ff) transitions. Un-less otherwise specified, groups are spaced logarithmicallywhile subgroups are treated uniformly. Additionally, there is agrey scattering opacity that is isotropic in the comoving framecalculated as (Castor 2004, p. 161) σ s = 8 π n e − (cid:18) e − m e − c (cid:19) , (86)where e − is electron charge, n e − is electron number density,and m e − is electron mass in cgs units. With mass fractionsknown a priori and given the assumption of LTE, the Saha-Boltzmann equations are used to obtain the excitation densitiesfor each atom in the W7 model (Mihalas & Mihalas 1984,p. 49). To calculate opacity, we introduce a subgrid for eachgroup g with index g g ∈ { . . . G g } . Values for bb opacitiesare calculated from oscillator strength data for each atomicspecies (Kurucz 1994). Furthermore, it is assumed that a lineis entirely included in the subgroup its line center is located.So (Mihalas & Mihalas 1984, pp. 329-332), σ a,g g ,bb =1∆ λ g g (cid:88) s (cid:88) i (cid:88) i (cid:48) >i (cid:18) π ( e − ) m e − c (cid:19) f i,i (cid:48) ,s λ i,i (cid:48) ,s c × [Θ( λ i,i (cid:48) ,s − λ g g − / ) − Θ( λ i,i (cid:48) ,s − λ g g +1 / )] × n i,s (1 − e hckTλi,i (cid:48) ,s ) , (87)where σ a,g g ,bb is the bb contribution to subgroup g g , f i,i (cid:48) ,s is the non-dimensional oscillator strength from state i to i (cid:48) of species s , λ i,i (cid:48) ,s is the wavelength center of the line cor-responding to the i → i (cid:48) transition, n i,s is the total densityof species s occupying state i , and the Θ are Heaviside stepfunctions constraining the sum to opacity profiles centered inthe subgroup. The bound-free opacities are tabulated accord-ing to the analytic fit prescription of Verner et al. (1996). Weapproximate the bf opacity, σ a,g g ,bf , of the subgroup as thevalue of the fit at the center wavelength in the subgroup. Theff opacities, σ a,g g ,ff , are computed with tabulated Gaunt fac-tors based on the work of Sutherland (1998) and are similarlyevaluated in the subgroup. The total absorption opacity forsubgroup g g is σ a,g g = σ a,g g ,bb + σ a,g g ,bf + σ a,g g ,ff (Mihalas& Mihalas 1984, p. 332). The total group opacity may thenbe averaged in some manner over the sub group contributions.We introduce an opacity mixing control parameter α σ ∈ [0 , to linearly combine reciprocal (“Rosseland type”) and directaverages of opacity. Averages of reciprocal opacity may prefer-entially weight lower opacity. For instance, Rosseland opacityis lower than Planck opacity. For some weight function, w ( λ ) ,the group absorption opacity is calculated as σ a,g = (1 − α σ ) 1 w g G g (cid:88) g g σ a,g g (cid:90) λ g +1 / λ g − / w ( λ ) dλ + α σ w g (cid:80) G g g g σ − a,g g (cid:82) λ g +1 / λ g − / w ( λ ) dλ , (88)where w g = (cid:82) λ g +1 / λ g − / w ( λ ) dλ . For a uniform weight function, Eq. (88) simplifies to σ a,g = (1 − α σ ) 1 G g G g (cid:88) g g σ a,g g + α σ G g (cid:80) G g g g /σ a,g g . (89)If LTE is considered, the weight function might be set to thenormalized Planck function; in this case Eq. (88) is a mix ofgrouped Planck and Rosseland opacities.For the W7 tests discussed, gamma ray energy depositionprofiles and the initial material and radiation temperaturesare borrowed from the PHOENIX code (Hauschildt 1992;Hauschildt & Baron 1999; Hauschildt & Baron 2004; vanRossum 2012). We estimate and apply a nominal value ofheat capacity of C v = 2 . × ρ erg/K/cm from Pinto &Eastman (2000) to compute the Fleck factor and update thematerial temperature. It has been found that changing C v bya factor of 3 does not change temperatures and spectra; theinsignificance of C v is attributable to the disparity of energystorage between the radiation and material fields. In the W7problem, the Fleck factor is found to be very small in IMC andIMC-DDMC. Consequently, even a modest group resolutionin IMC causes effective scattering to dominate particle pro-cesses. For the W7 tests attempted, it is apparently unfeasibleto use pure IMC, non-OR IMC-DDMC, or even IMC-DDMCwhere opacity regrouping is limited to adjacent groups. For a100 group W7 simulation with groups logarithmically spacedfrom × − cm to . × − cm, 64 velocity cells spaceduniformly from 0 cm/s to . × cm/s, a time domainof t ∈ [40 , days post explosion with 0.25 day time steps,250,000 initial particles, and 250,000 source particles per timestep, neither IMC nor non-OR IMC-DDMC completed thesimulation with 192 cores and a wall time of 40 hours each.In contrast, fully opacity-regrouped IMC-DDMC ( g c = 100 )completed the same problem with 24 cores in 1018.9 seconds.For the scope of this paper, we focus our attention to opacity-regrouped IMC-DDMC simulations.Our first W7 test problems explore the effect of differentgroup opacity averaging and group resolution. Specifically,Eq. (89) is implemented. The problems considered have 225,400, 625, and 1024 groups, 20 subgroups per group, and anopacity mixing parameter α σ ∈ { . , . , . , . , . } . Eachcalculation has 64 velocity cells uniformly spaced from 0cm/s to . × cm/s, 248 uniform time steps for t ∈ [2 , days, 250,000 initial radiation particles, 250,000 sourceparticles generated per time step, τ D = τ L = 3 mean freepaths, and the opacity-regrouped neighborhoods span the entireset of groups ( g c = G ). Absolute bolometric magnitudes arecalculated with M bol = 4 . − . (cid:18) L . × (cid:19) , (90)where L is luminosity in erg/s. The luminosities are computedby tallying lab frame particle energies escaping the domainand dividing by time step size. Figures 4a, 4b, 4c, and 4d havelight curves calculated with Eq. (90) for G = 225 , G = 400 , G = 625 , and G = 1024 , respectively, and a fixed number ofsubgroups, G g = 20 . Similarly Figs. 5a, 5b, 5c, and 5d havespectra at 20 days post explosion calculated with Eq. (90) for G = 225 , G = 400 , G = 625 , and G = 1024 , respectively,and G g = 20 . For the group resolutions presented, the α σ =1 . case does not appear to converge at the same rate as theother results. In other words, the α σ = 1 . case for Eq. (89)produces more sensitivity in brightness and spectrum versus7course group resolutions. As the mixing parameter is increasedtowards 1, the opacity calculation applies more reciprocalaveraging. Since reciprocal averaging favors smaller subgroupopacity values, it is expected that larger α σ yield earlier andbrighter light curves. Despite producing unrealistic light curvesfor α σ ≈ , α σ may be calibrated between 0 and ∼ . tomake simulations with low or modest group numbers emulatehigh-resolution simulations.Table 2 has computation times for each curve. Timing resultsfor the problem described are for 24 cores. With source particlenumbers kept constant, simulation time scales sub-linearlywith increasing group number. Table 2
Total Run Times for Opacity-Regrouped HMC W7 with 24 Cores (hours) G \ α σ We now examine the effect of the Gentile-Fleck factor, orEqs. (54) and (58) along with the optimization described inthe last paragraph of Section 7.1, on W7 temperatures. Fig-ure 6 has spectra and material temperature profiles shown atday 3 and 32 post-explosion for the W7 problem describedwith α σ = 0 . and G = 225 . At early times ( t (cid:46) days),both IMC-DDMC and modified IMC-DDMC yield outer-celltemperature fluctuations for the numerical specifications con-sidered. The fluctuations are different between the standardand modified methods. Consequently, the application of theGentile-Fleck factor in IMC-DDMC uncovers some uncer-tainty in early spectra. At later times ( t (cid:38) days), theGentile-Fleck factor yields consistently smoother material tem-perature profiles than the standard Fleck factor. However, thespectra at later times are not significantly affected by the fluc-tuations in the outer-cell temperatures because that region isoptically thin at that point.Finally, we compare the results of SuperNu and
PHOENIX for the W7 problem in LTE. We find that the light curve gen-erated by
SuperNu is systematically ∼
10% dimmer at peakthan the light curve generated by
PHOENIX for various timestep and group resolutions. For controlled testing, groupedopacities have been introduced into
PHOENIX . The multigroupcomputations have no opacity mixing, or α σ = 0 . Figure 7 has500 group light curve results from PHOENIX and
SuperNu along with a standard, high-resolution (30,000 wavelengthpoints)
PHOENIX light curve. From inspection of Fig. 7b, itis worth noting that the luminosities of multigroup
PHOENIX and
SuperNu have similar early rising light curves. Thismeans that the different diffusion treatments in the two codesare in good agreement. The standard
PHOENIX light curverises earlier than the multigroup
PHOENIX light curve, as ex-pected. This effect can be emulated in low group resolutionsimulations using the opacity mixing parameter (see Figure4). Increasing α σ from 0 to ∼ . has a similar effect on thelight curve shape as increasing the resolution to convergence.Figure 8 has spectra at 10, 20, and 40 days post-explosion forthe 500 group SuperNu and high-resolution
PHOENIX simu-lations. Despite differences in magnitudes, the time evolutionof the light curves and the shapes of the spectra are in goodagreement. The codes use the same atomic data but the EOS and opacity routines are different; these factors may accountfor some differences in the luminosities and spectra.Resolving the sources of the 10-15% discrepancy will re-quire more in-depth code-to-code comparisons which is workin progress but beyond the scope of this paper. Having per-formed time step and group resolution tests, we also planto perform resolution tests on the spatial grid. It is possi-ble the codes have different convergence properties with gridresolution. In particular, the standard leakage opacity at IMC-DDMC spatial method interfaces may underpredict particletransmission across cell surfaces when DDMC interface cellsare optically thick (Densmore et al. 2007). Densmore et al.(2006) performs an emissivity based derivation to generalizethe standard IMC-DDMC boundary condition and improvethe emission from DDMC to IMC at spatial interfaces. If in-creased grid resolution in
SuperNu increases the luminosity,then the alternate boundary condition presented by Densmoreet al. (2006) may increase the absolute bolometric magnitudeof the light curve at the current 64 cell resolution. We haveperformed preliminary tests with an emissivity based boundarycondition and find a ∼
2% increase in the absolute bolometricmagnitude at peak; despite this modest change, exploring theeffects of increasing the spatial resolution may be revealing.Apart from grid resolution, EOS, opacities, and transport meth-ods, there may be other important reasons for the observeddifferences. CONCLUSIONS AND FUTURE WORKWe have incorporated techniques to mitigate overheatingerrors and combine DDMC groups with high opacity in theIMC-DDMC code,
SuperNu . In Section 6, we describedan approach to Doppler shift DDMC particles. The Dopplershift scheme accounts for the effect of inelastic collisionswith uniform subgroup redistribution. Following Abdikamalovet al. (2012), the Doppler shift scheme is operator split fromthe diffusion scheme; it does not conflict with the opacityregrouping process.We found that opacity regrouping is needed in IMC-DDMCto make the W7 problem feasible; the optimization mitigatescomputational cost in performing the multidimensional calcu-lation. Additionally, we have described and tested an approachto treating the opacity that involves refining the wavelengthgrid to subgroups.In Section 7.1 we used the Gentile-Fleck factor to mitigatean overheating pathology in the presence of strong outflow.The MC results are benchmarked against a quasi-manufacturedsolution. In Section 7.2, we treated structured multigroup prob-lems with IMC-DDMC to test the effect of non-contiguousopacity regrouping. For the problem presented, opacity re-grouping significantly improves efficiency without a significantcost of accuracy in the temperatures and spectra. In Section 7.3,we tested IMC-DDMC with opacity regrouping and subgroup-ing on the W7 problem. We also compared light curves andspectra for the W7 test problem calculated using
SuperNu and
PHOENIX for a similar set-up. We modified
PHOENIX tobe able to use multigroup opacities, which enabled us to domore controlled code-to-code comparisons. The light-curverise times given by multigroup
PHOENIX and
SuperNu arein good agreement for the same group resolution. We findsatisfactory agreement in the shape of the spectra. However,there exists a ∼ SuperNu and
PHOENIX in the luminosity of the light curve around and afterpeak that is currently not fully understood. Time step resolu-tion tests indicate the light curves compared between codes8 A b s o l u t e B o l o m e t r i c M agn i t ude HMC, G=225, α σ =0.0HMC, G=225, α σ =0.3HMC, G=225, α σ =0.5HMC, G=225, α σ =0.8HMC, G=225, α σ =1.0 (a) A b s o l u t e B o l o m e t r i c M agn i t ude HMC, G=400, α σ =0.0HMC, G=400, α σ =0.3HMC, G=400, α σ =0.5HMC, G=400, α σ =0.8HMC, G=400, α σ =1.0 (b) A b s o l u t e B o l o m e t r i c M agn i t ude HMC, G=625, α σ =0.0HMC, G=625, α σ =0.3HMC, G=625, α σ =0.5HMC, G=625, α σ =0.8HMC, G=625, α σ =1.0 (c) A b s o l u t e B o l o m e t r i c M agn i t ude HMC, G=1024, α σ =0.0HMC, G=1024, α σ =0.3HMC, G=1024, α σ =0.5HMC, G=1024, α σ =0.8HMC, G=1024, α σ =1.0 (d) Figure 4.
Opacity-regrouped IMC-DDMC W7 bolometric light curves for opacity mixing α σ ∈ { . , . , . } (solid) and α σ ∈ { . , . } (dashed; so solidand dashed curves alternate versus α σ ) and a fixed number of subgroups, G g = 20 . Equation (89) has been applied for opacity mixing. Light curves are calculatedby tallying particles that have escaped the spatial (velocity) domain per time step and applying Eq. (90). In Figs. 4a, 4b, 4c, and 4d group resolutions are G = 225 , G = 400 , G = 625 , and G = 1024 , respectively. As expected, peak luminosity is earlier and brighter for opacity mixing that favors reciprocal averaging sincesmaller subgroup opacity values are favored. Values of α σ close to 1 are not realistic as opacities of strong absorption lines are more and more neglected. Theopacity mixing parameter can be used to calibrate simulations with modest group resolution to emulate the diffusion characteristics of equivalent high-resolutionsimulations. are converged in time. For certain spatial grid resolutions,DDMC may underpredict spatial leakage of diffusion particlesto IMC (Densmore et al. 2006, 2007). Consequently, spatialgrid resolution tests of SuperNu may be informative.We plan to extend our code to multiple dimensions. TheIMC-DDMC method is simple to extend to two and threedimensions for simple grid geometries. The challenges in per-forming multidimensional simulations of SN Ia light curvesand spectra with IMC-DDMC lies in optimization and mem-ory requirements. In addition to spatial geometry, we planto investigate methods and algorithms that further mitigatespurious temperature spikes due to the Maximum Principle orMC noise. ACKNOWLEDGEMENTSWe would like to thank Donald Lamb, Gregory Moses, andCarlo Graziani for supporting and guiding this work. Wewould like to thank Donald Lamb for the constructive recom- mendations and suggestions. We especially thank our referee,Ernazar Abdikamalov, for the valuable recommendations thatimproved this paper. This research was supported in partby the NSF under grant AST-0909132, and by NIH throughresources provided by the Computation Institute and the Bi-ological Sciences Division of the University of Chicago andArgonne National Laboratory, under grant S10 RR029030-01.This work is supported in part at the University of Chicago bythe National Science Foundation under grant PHY-0822648for the Physics Frontier Center ”Joint Institute for NuclearAstrophysics” (JINA).
REFERENCESAbdikamalov, E., Burrows, A., Ott, C. D., Loffler, F., O’Connor, E., Dolence,J. C., & Schnetter, E. 2012, ApJ, 755, 111Adams, M. L. 2001, Nucl. Sci. Eng., 137Atzeni, S., & ter Vehn, J. M. 2004, The Physics of Inertial Fusion (OxfordUniversity Press) Wavelength [ang] F l u x [ e r g / s / c m ] HMC, G=225, α σ =0.0 (20.25 days)HMC, G=225, α σ =0.3 (20.25 days)HMC, G=225, α σ =0.5 (20.25 days) (a) Wavelength [ang] F l u x [ e r g / s / c m ] HMC, G=400, α σ =0.0 (20.25 days)HMC, G=400, α σ =0.3 (20.25 days)HMC, G=400, α σ =0.5 (20.25 days) (b) Wavelength [ang] F l u x [ e r g / s / c m ] HMC, G=625, α σ =0.0 (20.25 days)HMC, G=625, α σ =0.3 (20.25 days)HMC, G=625, α σ =0.5 (20.25 days) (c) Wavelength [ang] F l u x [ e r g / s / c m ] HMC, G=1024, α σ =0.0 (20.25 days)HMC, G=1024, α σ =0.3 (20.25 days)HMC, G=1024, α σ =0.5 (20.25 days) (d) Figure 5.
Opacity-regrouped IMC-DDMC W7 spectra for opacity mixing α σ = 0 . , . , . (dotted, dashed, and solid, respectively) and a fixed number ofsubgroups, G g = 20 . Equation (89) has been applied for opacity mixing. Spectra are calculated by tallying escaping particles energies per group per time anddividing by group wavelength range. Data are plotted at group centers. In Figs. 5a, 5b, 5c, and 5d group resolutions are G = 225 , G = 400 , G = 625 , and G = 1024 , respectively. Locations of peaks and troughs amongst the different opacity mixings presented appear consistent. For λ ∈ [2000 , , radiationtransmission is larger for larger values of α σ .Baron, E., & Hauschildt, P. H. 2007, A&A, 468, 255Branch, D., & Khokhlov, A. 1995, Physics Reports, 256, 53Brooks, E. D. 1989, J. Comput. Phys., 83Buchler, J. R. 1983, JQSRT, 30, 395Calder, A. C., Plewa, T., Vladimirova, N., Lamb, D. Q., & Truran, J. W. 2004,Astrophysical Journal, LettersCalder, A. C., et al. 2002, Astrophysical Journal, Supplement, 143, 201Carter, L. L., & Forest, C. A. 1973, lA-5038, Los Alamos National LaboratoryCastor, J. I. 2004, Radiation Hydrodynamics (Cambridge University Press)Cleveland, M. A., & Gentile, N. 2014, Transport Theory and StatisticalPhysics, 1Cleveland, M. A., Gentile, N. A., & Palmer, T. S. 2010, J. Comput. Phys., 229,5707Densmore, J. D. 2011, J. Comput. Phys., 230, 1116Densmore, J. D., Davidson, G., & Carrington, D. B. 2006, Ann. Nucl. Energy,33, 583Densmore, J. D., Evans, T. M., & Buksas, M. W. 2008, Nucl. Sci. Eng., 159, 1Densmore, J. D., & Larsen, E. W. 2004, J. Comput. Phys., 199, 175Densmore, J. D., Thompson, K. G., & Urbatsch, T. J. 2012, J. Comput. Phys.,231, 6925Densmore, J. D., Urbatsch, T. J., Evans, T. M., & Buksas, M. W. 2007, J.Comput. Phys., 222, 485Fleck, Jr., J. A., & Canfield, E. H. 1984, J. Comput. Phys., 54, 508Fleck, Jr., J. A., & Cummings, J. D. 1971, J. Comput. Phys., 8, 313Fryxell, B., et al. 2000, ApJS, 131, 273 Gamezo, V. N., Khokhlov, A. M., Oran, E. S., Chtchelkanova, A. Y., &Rosenberg, R. O. 2003, Science, 299, 77Gentile, N. A. 2001, J. Comput. Phys., 172, 543—. 2011, J. Comput. Phys., 230Habetler, G. J., & Matkowsky, B. J. 1975, J. Math. Phys., 16, 846Hauschildt, P. H. 1992, JQSRT, 47, 433Hauschildt, P. H., & Baron, E. 1999, Journal of Computational and AppliedMathematics, 109Hauschildt, P. H., & Baron, E. 2004, A&A, 417, 317Hauschildt, P. H., & Wehrse, R. 1991, JQSRT, 46Hillebrandt, W., & Niemeyer, J. 2000, ARA&A, 38, 191Kasen, D., Thomas, R. C., & Nugent, P. 2006, ApJ, 651, 366Kromer, M., & Sim, S. A. 2009, Mon. Not. R. Astron. Soc., 398Kurucz, R. L. 1994Larsen, E. W., & Mercier, B. 1987, J. Comput. Phys., 71Lewis, E. E., & Miller, Jr., W. F. 1993, Computational Methods of NeutronTransport (American Nuclear Society)Long, M., et al. 2013, ArXiv e-printsLowrie, R. B., Mihalas, D., & Morel, J. E. 2001, JQRST, 69, 291Lucy, L. B. 2005, A&A, 429, 19Malvagi, F., & Pomraning, G. C. 1991, J. Math. Phys., 32, 805McClarren, R. G., Holloway, J. P., & Brunner, T. A. 2008a, JQSRT, 109, 389McClarren, R. G., Lowrie, R. B., Prinja, A. K., & Morel, J. E. 2008b, JQSRT,109, 2590McClarren, R. G., & Urbatsch, T. J. 2009, J. Comput. Phys., 228, 5669 Velocity [cm/s] M a t e r i a l T e m pe r a t u r e [ K ] HMC (6 days)Mod HMC (6 days) (a) Wavelength [ang] F l u x [ e r g / s / c m ] HMC (6 days)Mod HMC (6 days) (b) M a t e r i a l T e m pe r a t u r e [ K ] HMC (32 days)Mod HMC (32 days) (c) Wavelength [ang] F l u x [ e r g / s / c m ] HMC (32 days)Mod HMC (32 days) (d)
Figure 6.
IMC-DDMC (solid) and Gentile-Fleck factor modified IMC-DDMC (dashed) material temperatures (left) and spectra (right) at days 6 and 32post-explosion for the W7 problem described in Section 7.3 with G = 225 and α σ = 0 . . At early time ( t (cid:46) days), the Gentile-Fleck factor slightly modifiesthe spectrum which is sensitive to its effect on the fluctuations in the outer regions of the ejecta. At late times in the W7 expansion ( t (cid:38) days), the Gentile-Fleckfactor consistently mitigates temperature fluctuations in the outer cells. Despite the continued temperature fluctuations in the outer cells for standard IMC-DDMCat later times, the difference in spectra at late times is no longer significant.McClarren, R. G., & Urbatsch, T. J. 2012, in Transactions of the AmericanNuclear SocietyMcKinley, M. S., Brooks, E. D., & Sz˝oke, A. 2003, J. Comput. Phys., 189,330Mihalas, D., & Mihalas, B. W. 1984, Foundations of RadiationHydrodynamics (Oxford University Press)N’Kaoua, T. 1991, SIAM J. Stat. Comput., 12, 505Nomoto, K., Thielemann, F., & Yokoi, K. 1984, ApJ, 286, 644Oberkampf, W. L., & Roy, C. J. 2010, Verification and Validation in ScientificComputing (Cambridge University Press)Olson, G. L., & Kunasz, P. B. 1987, JQSRT, 38Perlmutter, S. 2003, Physics Today, 53Perlmutter, S., et al. 1999, ApJ, 517, 565Petschek, A. 1990, Supernovae (Springer-Verlag)Phillips, M. M. 1993, ApJ, 413, L105Pinto, P. A., & Eastman, R. G. 2000, ApJ, 530, 744Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics(Pergamon Press) Riess, A. G., et al. 1998, AJ, 116, 1009Scannapieco, C., Tissera, P. B., White, S. D. M., & Springel, V. 2008,MNRAS, 389Seitenzahl, I. R., et al. 2013, MNRAS, 429, 1156Su, B., & Olson, G. L. 1999, JQSRT, 62, 279Sutherland, R. S. 1998, MNRAS, 300Sz˝oke, A., & Brooks, E. D. 2005, JQSRT, 91, 95Urbatsch, T. J., Morel, J. E., & Gulick, J. C. 1999, in Proc. Int. Conf.Mathematics and Computation, Reactor Physics, and EnvironmentAnalysis in Nuclear Applicationsvan Rossum, D. R. 2012, ApJ, 756, 31Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465Warsa, J. S., & Densmore, J. D. 2010, Nucl. Sci. Eng., 166, 36Wollaeger, R. T., van Rossum, D. R., Graziani, C., Couch, S. M., Jordan,G. C., Lamb, D. Q., & Moses, G. A. 2013, ApJS, 209 A b s o l u t e B o l o m e t r i c M a g n i t u d e SuperNu G=500Phoenix G=500Phoenix high-res (a) A b s o l u t e B o l o m e t r i c M a g n i t u d e SuperNu G=500Phoenix G=500Phoenix high-res (b)
Figure 7.
SuperNu (blue), with multigroup
PHOENIX (green), and standard
PHOENIX (red) light curves.
SuperNu and multigroup
PHOENIX apply 500groups and directly averaged group opacity, or an α σ = 0 mix. PHOENIX is run in LTE for consistency with
SuperNu . There exists a systematic difference of ∼ PHOENIX light curve rises earlier than multigroup
PHOENIX . This is due tothe high resolution that enables windows of lower opacity through which diffusion is enhanced. Diffusion at low group resolutions can be simulated and calibratedusing the opacity mixing parameter α σ (see Figure 4). F l u x [ e r g / s / c m ] (a) F l u x [ e r g / s / c m ] (b) F l u x [ e r g / s / c m ] (c) F l u x [ e r g / s / c m ] (d) Figure 8.
SuperNu (blue) with 500 groups and standard
PHOENIX (green) spectra for the W7 problem at 10, 20, 30, and 40 days post-explosion. In Fig. 8a,the difference in flux is partly attributable to the earlier rise of the
PHOENIX high-resolution luminosity (see Fig. 7). In Fig. 8b, the W7 supernova is near peakluminosity; resolving the discrepancy in flux requires further code-to-code comparison. In Fig. 8d, the flux of
PHOENIX is not systematically larger than
SuperNu .Around day 40, the high-resolution
PHOENIX light curve is at a lower luminosity than the 500 group