A DG-IMEX method for two-moment neutrino transport: Nonlinear solvers for neutrino-matter coupling
M. Paul Laiu, Eirik Endeve, Ran Chu, J. Austin Harris, O. E. Bronson Messer
DD RAFT VERSION F EBRUARY
4, 2021Typeset using L A TEX default style in AASTeX63
A DG-IMEX method for two-moment neutrino transport: Nonlinear solvers for neutrino-matter coupling ∗ M. P
AUL L AIU , E IRIK E NDEVE ,
1, 2 R AN C HU , J. A
USTIN H ARRIS , AND
O. E. B
RONSON M ESSER
3, 4, 21
Multiscale Methods Group, Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA Department of Physics and Astronomy, University of Tennessee Knoxville, Knoxville, TN 37996, USA National Center for Computational Sciences, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA Physics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA (Received —; Revised —; Accepted —)
Submitted to ApJSABSTRACTNeutrino-matter interactions play an important role in core-collapse supernova (CCSN) explosions as theycontribute to both lepton number and/or four-momentum exchange between neutrinos and matter, and thusact as the agent for neutrino-driven explosions. Due to the multiscale nature of neutrino transport in CCSNsimulations, an implicit treatment of neutrino-matter interactions is desired, which requires solutions of couplednonlinear systems in each step of the time integration scheme. In this paper we design and compare nonlineariterative solvers for implicit systems with energy coupling neutrino-matter interactions commonly used in CCSNsimulations. Specifically, we consider electron neutrinos and antineutrinos, which interact with static matterconfigurations through the Bruenn 85 opacity set. The implicit systems arise from the discretization of a non-relativistic two-moment model for neutrino transport, which employs the discontinuous Galerkin (DG) methodfor phase-space discretization and an implicit-explicit (IMEX) time integration scheme. In the context of thisDG-IMEX scheme, we propose two approaches to formulate the nonlinear systems — a coupled approachand a nested approach. For each approach, the resulting systems are solved with Anderson-accelerated fixed-point iteration and Newton’s method. The performance of these four iterative solvers has been compared onrelaxation problems with various degree of collisionality, as well as proto-neutron star deleptonization problemswith several matter profiles adopted from spherically symmetric CCSN simulations. Numerical results suggestthat the nested Anderson-accelerated fixed-point solver is more efficient than other tested solvers for solvingimplicit nonlinear systems with energy coupling neutrino-matter interactions.
Keywords:
Computational methods (1965), Core-collapse supernovae (304), Radiative transfer simulations(1967), Supernova neutrinos (1666) INTRODUCTIONCore-collapse supernovae (CCSNe), the explosive deaths of massive stars, are to a large extent neutrino driven. About 99%of the gravitational potential energy released in a core-collapse event ( ∼ erg) is radiated away by neutrinos, which alsoact as a driver for the expulsion of matter. Near the end of a massive star’s life (i.e. a star with mass exceeding about tensolar masses) its iron core, which does not produce energy by nuclear burning, is held up against gravity by the pressure fromdegenerate electrons. However, the mass of the iron core continues to increase due to silicon burning on its surface. Oncethe mass of the iron core reaches about . solar masses (the Chandrasekhar limit), the electron degeneracy pressure becomesinsufficient in balancing gravity, and the core collapses in on itself. The collapse proceeds until the central rest mass densityexceeds nuclear matter densities ( > g cm − ), when the matter equation of state (EoS) stiffens to halt the collapse and Corresponding author: M. Paul [email protected] ∗ This manuscript has been authored, in part, by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE).The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclu-sive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US gov-ernment purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan( http://energy.gov/downloads/doe-public-access-plan ). a r X i v : . [ a s t r o - ph . H E ] F e b L AIU ET AL .a shock wave is launched into the collapsing outer core. The outward-propagating shock wave loses energy by dissociatingiron nuclei into free nucleons. Furthermore, electron capture on nucleons in the hot matter behind the shock produces copiousamounts of neutrinos and antineutrinos, which can escape the system once the shock reaches sufficiently low densities (about g cm − ). The combination of neutrino emission and iron dissociation weakens the shock, which eventually stalls at a radiusof about − km from the center of the star. At this point, the region below the shock can be divided into the cooling layer and the gain layer , separated by the gain surface, where neutrino heating and cooling balance. In the cooling layer, extendingfrom the surface of the proto-neutron star to the gain surface, there is net energy loss by neutrino emission. At lower densities, inthe heating layer extending from the gain surface to the shock surface, there is net heating by neutrinos emanating from below;see diagram in Figure 1. In the neutrino reheating explosion mechanism (Bethe & Wilson 1985), energy deposition by neutrinosin the heating layer revives the supernova shock wave to disrupt the massive star in a CCSN explosion. This basic descriptionis supported by recent numerical simulations (e.g., Lentz et al. (2015); Melson et al. (2015); Burrows et al. (2020)), but thedetails are more complicated: The CCSN explosion emerges from a complex interplay between between neutrino transport andhydrodynamic (or magnetohydrodynamic) processes, playing out within a curved spacetime (see, e.g., Janka (2012); Burrows(2013); Hix et al. (2014); M¨uller (2016) for reviews). r&ee < Trot = 4 X lO' 51 M R
10 km 2 ms - 2 ergs (1) 1.4M0 [22]. Thus rather small initial rotation periods, < 2 ms, for newly formed neutron stars are required if enough magnetic energy is to potentially arise to power up the typical supernova. These small rotation periods are at variance with the calculated rotational periods of the magnetized cores of supernova progenitors [23], and the extrapolated periods of newly formed neutron stars (> 10 ms). Both of these constraints are "soft" (stellar evolutionary calculations with rotation and magnetic fields are not ab initio, and we have not yet observed a newly formed neutron star), but if they hold then the MHD mechanism will only be relevant to a subset of core collapse supemovae. However, the observations of magnetars, long-duration gamma-ray bursts, and hints of highly coUimated material in some supernova remnants suggests a subclass of events that are magnetically driven.
Neutrino Heating Mechanism
SASIi
FIGURE 1.
Entropy and velocity configuration snapshots of the Model llM_Sym_32R.
The neutrino heating mechanism has a long pedigree extending back to the seminal paper of Colgate and White [24] and its more modem incamation [25]. Following the collapse and bounce of the inner core of massive star at the endpoint of its normal thermonuclear evolution, the shock launched at core bounce stalls in the outer core, losing energy (and therefore post-shock pressure) to nuclear dissociation and electron neutrino losses. Within a short time ('-^
50 ms) a thermodynamic profile is established as illustrated in Figure 1 in which infalling matter encountering the outward flow of 595
Downloaded 08 Feb 2010 to 128.219.49.9. Redistribution subject to AIP license or copyright; see http://proceedings.aip.org/proceedings/cpcr.jsp
Figure 1.
Schematic diagram of the neutrino reheating phase during a CCSN explosion. (This figure is originally from Bruenn et al. (2009),and is reproduced here with permission from AIP Publishing.) Neutrinos emanate from the neutrinosphere — a proxy for the surface of theproto-neutron star (gray shaded region) — and the cooling layer above it. Neutrino heating mediated by neutrino-matter interactions in the gainlayer between the gain radius and the shock provides energy to reinvigorate the stalled supernova shock. During shock revival, the flow belowthe shock may be modulated by neutrino-driven convection and the standing accretion shock instability (SASI).
Multiple interaction processes contribute to lepton number and/or four-momentum exchange between neutrinos and matterin the CCSN explosion. Neutrino emission caused by electron capture on nucleons and the inverse process of neutrino capture(absorption) dominate and contribute to both lepton number and four-momentum exchange. During stellar core collapse, electroncapture on nuclei is critical to the collapse dynamics; e.g., Hix et al. (2003). Neutrino scattering on electrons and nucleons (andnuclei during collapse) are also major contributors to the neutrino opacity. In the seminal work of Bruenn (1985), neutrino-nucleon scattering is treated as isoenergetic (the neutrino changes direction but not energy in its encounter with the nucleon),while full four-momentum exchange is considered for neutrino-electron scattering (NES). In later works (e.g., Thompson et al.(2000); Lentz et al. (2012a); M¨uller et al. (2012); Burrows et al. (2018)), it has been demonstrated that neutrino-nucleon scatteringcontributes to neutrino-matter thermalization in a nontrivial way. (Without modern nuclear electron capture rates, NES also helpsshape conditions in the core prior to the formation of the supernova shock (Mezzacappa & Bruenn 1993; Lentz et al. 2012a).)Neutrino pair creation and annihilation, via electron-positron pairs (Bruenn 1985) or nucleon-nucleon bremsstrahlung (Hannestad
OLVERS FOR NEUTRINO - MATTER COUPLING
3& Raffelt 1998), is yet another significant process in CCSN neutrino transport. In particular, nucleon-nucleon bremsstrahlung isa dominant source of µ and τ neutrinos. Several studies (e.g., Lentz et al. (2012a); Just et al. (2018); Burrows et al. (2018)) haveinvestigated the impact of various neutrino opacities on the CCSN explosion mechanism. Although it is difficult to pin down anexact set of necessary interactions, there is a consensus view that processes that couple globally in neutrino momentum-spaceand/or across neutrino species (e.g., NES and pair processes) must be included for a realistic description.It is computationally expensive to include neutrino transport with satisfactory realism in simulations of CCSNe. First ofall, neutrinos interact relatively weakly with matter in the gain region, which demands a kinetic description. In addition, thehydrodynamics must be modeled without imposed spatial symmetries and with adequate resolution to capture processes thatshape the explosion (see e.g., M¨uller (2020)). The additional requirement of including momentum-space coupling neutrino-matter interactions makes realistic, large scale simulations a computational challenge. For example, if momentum space isdiscretized with N p points, a simple evaluation of the collision operator at one spatial location requires O ( N p ) operations, asopposed to O ( N p ) operations for the emission, absorption and isoenergetic scattering operators. In current three-dimensionalsupernova models, the dimensionality of the neutrino transport problem is reduced by adopting one- or two-moment approachesthat retain the neutrino energy dimension of momentum-space (Just et al. 2015; Skinner et al. 2019; Bruenn et al. 2020), but thisdoes not completely alleviate the computational challenge of including the critical energy coupling interactions.Because the neutrino mean free path can be much shorter than the spatial resolution afforded in parts of the computationaldomain, an implicit treatment of neutrino-matter interactions is desired. In spherically symmetric models (Rampp & Janka 2002;Liebend¨orfer et al. 2004) and multi-dimensional models employing the so-called ray-by-ray approximation (Bruenn et al. 2020),the full set of transport equations is commonly solved with implicit time integration. So far, fully implicit time integration has notbeen the method of choice for truly multi-dimensional neutrino transport (but see Sumiyoshi & Yamada (2012) for an exception).Instead, implicit-explicit (IMEX) methods (Ascher et al. 1997; Pareschi & Russo 2005) have received more attention (see, e.g.,O’Connor (2015); Just et al. (2015); Kuroda et al. (2016); Skinner et al. (2019); Chu et al. (2019)). In the IMEX approach,collisions are treated with implicit methods, while phase-space advection is treated with explicit methods (implicit integrationfor the momentum-space advection terms is also used (Kuroda et al. 2016)). Explicit integration for phase-space advection isadvantageous because it avoids solving a distributed, sparse system of nonlinear equations, and, since neutrino-matter interactionsare completely local in space, the implicit part is embarrassingly parallel. Moreover, the characteristic wave speeds associatedwith the transport and hydrodynamics equations are not too dissimilar in relativistic systems, and it is not clear that the expectedadditional cost of solving the full transport equation with implicit time integration and a larger time step will pay off.Energy coupling neutrino-matter interactions still dominate the computational cost of IMEX-based neutrino transport schemesemploying a spectral two-moment approach, where the neutrino energy domain is discretized with N ε points. The cost per spatialpoint is expected to scale roughly as O ( N pε ) , where the power p is between two and three. The lower bound ( p = 2 ) is motivatedby the expected cost of simply evaluating the energy coupling operators, while the p = 3 scaling is due to the cost of invertingthe Jacobian matrix resulting from an implicit solution algorithm based on Newton’s method. In addition, the neutrino-matterinteraction rates depend nonlinearly on the local thermodynamic properties of the matter, so that in a fully implicit approach(Kuroda et al. 2016), the local computational cost scales linearly with the number of iterations needed to reach convergence.Several studies have investigated approximations that are motivated by potential gains in computational expediency. To alle-viate the cost of including energy coupling neutrino-matter interactions, Just et al. (2015) investigated the effect of evaluatingthe collision terms for energy coupling interactions with matter conditions taken from the known state at high densities (above g cm − ), and using explicit integration at lower densities, and found results that were practically identical to a run witha more implicit treatment. Another approach to circumvent stiffness induced by neutrino-matter interactions is to artificiallyreduce the rates at high densities and use explicit time integration throughout the computational domain (Thompson et al. 2003;O’Connor 2015; Burrows et al. 2018; Just et al. 2018). We also note a simplifying approach to pair processes involving approxi-mating pair annihilation partners by local equilibrium distributions, which essentially renders this interaction local in momentum-space (O’Connor 2015). While some approximations have shown to work well in a limited set of comparisons, they do introduceuncertainties and limit the applicability of the algorithm. Partially for these reasons, we do not follow the approach of alteringthe kernel or collision operator based on physical insight into the problem. Instead we seek to develop fully implicit solvers forenergy coupling neutrino-matter interactions for CCSN simulations and related applications.This paper details the base neutrino transport algorithms implemented in the toolkit for high-order neutrino-radiation hydro-dynamics ( thornado ), which is being developed for simulations of CCSNe and related problems using the discontinuousGalerkin (DG) method for phase-discretization. Specifically, the phase-space discretization in thornado is based on the nodal https://github.com/endeve/thornado, http://dx.doi.org/10.13139/OLCF/1735948 L AIU ET AL .DG method (see, e.g., Hesthaven & Warburton (2008)). The original DG method was developed by Reed & Hill (1973) forsolving neutron transport problems. Since then, it has been extended to the Runge-Kutta discontinuous Galerkin framework forsolving more general hyperbolic partial differential equations (see, e.g., Cockburn & Shu (1998); Cockburn et al. (1990); Cock-burn et al. (1989); Cockburn & Shu (1989, 1991); Cockburn & Shu (2001) for early developments). For more recent developmentsof DG methods, see Shu (2016) and references therein. DG methods are particularly attractive for transport problems since theyrecover the correct asymptotic behavior in the so-called diffusion limit (Larsen & Morel 1989; Adams 2001). They can alsobe easily applied to problems with curvilinear coordinates — necessary when solving general relativistic problems (Teukolsky2016). However, the DG method has so far not been applied to neutrino transport in CCSN models (but see Radice et al. (2013);Endeve et al. (2015); Chu et al. (2019) for applications in simplified settings). The DG method approximates solutions withpiecewise local polynomials, and tracks the evolution of coefficients associated with the polynomial expansion, thus, they areoften referred to as modal DG methods. The nodal DG method uses a particular interpolating polynomial to construct the approx-imation, which allows it to track the evolution of nodal values at the interpolation points. This special polynomial approximationresults in a simple projection operator from the target function to the polynomial space, which enables straightforward parallelimplementation of the nodal DG method (Kl¨ockner et al. 2009). The nodal DG method has been used in various applications,including solving kinetic equations (e.g., Xiong et al. (2015); Juno et al. (2018)).In this paper we design and evaluate nonlinear solvers for neutrino-matter interactions in a two-moment model for neutrinotransport within the IMEX framework. This DG-IMEX method is essentially the same as that described by Chu et al. (2019),but is extended to include curvilinear spatial coordinates, multiple neutrino species, more realistic interactions, and coupling to amaterial background governed by a nuclear EoS. Specifically, we consider electron neutrinos and antineutrinos, and develop non-linear solvers for the opacity set of Bruenn (1985), which includes emission and absorption due to electron and neutrino captureon nucleons and nuclei, isoenergetic neutrino scattering off nucleons and nuclei, inelastic NES, and neutrino-antineutrino pairproduction/annihilation from electron-positron pairs. We use the SFHo EoS (Steiner et al. 2013) in the numerical experiments.The microphysics (neutrino opacities and EoS) has been tabulated by the W
EAK L IB library , which also provides routines foraccess and manipulation (e.g., interpolation and differentiation) of tabulated microphysics data.Several nonlinear solver strategies are considered for the neutrino-matter coupling problem. As a baseline for comparison weconsider Newton’s method, which can offer rapid convergence to the solution if the initial guess is sufficiently close and theobjective function sufficiently regular. However, the necessity of using approximate derivatives due to tabulated opacity kernelsto form the Jacobian matrix can hamper the convergence speed. In addition, the construction of the Jacobian from tabulated dataand the solution of a dense linear system for each iteration is computationally expensive. We consider fixed-point iteration asan alternative to Newton’s method. With fixed-point iteration, Jacobian matrix constructions and dense linear system solutionsare not necessary, but the rate of convergence can be slow. We employ Anderson acceleration to improve the convergence rateof the fixed-point method. This acceleration technique was first proposed by Anderson (1965) for solving integral equations,and has been used to accelerate fixed-point solutions in several applications, including solving radiation-diffusion equations (Anet al. 2017), flow problems (Lott et al. 2012), nuclear reactor simulations (Hamilton et al. 2016), as well as a variety of nonlinearproblems (Walker & Ni 2011). Anderson acceleration speeds up the convergence of standard fixed-point iterations by taking anextrapolation step based on the recent iterates that aims to minimize the residual of the new iterate. The convergence propertiesof Anderson acceleration were analyzed in Toth & Kelley (2015); Kelley (2018); Evans et al. (2020), which include (i) globalconvergence on linear problems under the standard contraction assumption, (ii) local convergence on nonlinear problems underassumptions similar to the standard ones for local convergence of Newton’s method, and (iii) improved local convergence ratewhen applied on linearly converging fixed-point iterations.Although the Anderson-accelerated fixed-point algorithm is generally faster than Newton’s method in the cases investigatedin this paper, we find that the computational cost associated with reevaluating the neutrino opacities in each iteration remainsrelatively high. This observation motivates a nested approach, where the matter quantities are updated in an outer iteration loop,outside an inner iteration loop where the radiation field is iterated to convergence while the matter state is fixed. We considertwo nested iteration schemes, both based on Anderson-accelerated fixed-point iteration in the outer loop. In the first case theinner loop solve is based on Newton’s method. In this nested approach, the Jacobian associated with the inner solve can becomputed analytically since the nonlinear functional is considered independent of the matter state. However, the dense linearsystem to be solved is of similar size as for the fully coupled Newton method. In the second case, the inner loop solve is basedon Anderson-accelerated fixed-point iteration. We find that the nested iteration schemes require fewer opacity evaluations andare more efficient than the fully coupled schemes. We also find that the nested scheme with inner Newton iterations requires https://github.com/starkiller-astro/weaklib, http://dx.doi.org/10.13139/OLCF/1735948 OLVERS FOR NEUTRINO - MATTER COUPLING thornado matures, we intend to account for this physics, and document on the performance in future publications. To the bestof our knowledge, the present paper documents the most advanced application of the DG method to neutrino transport.This paper is organized as follows: in Section 2 we present the mathematical model we adopt for neutrino transport; inSection 3 we introduce the DG-IMEX scheme implemented in thornado ; the nonlinear solvers are detailed in Section 4; andin Section 5 we present numerical results, where we compare the performance of the solvers on (1) relaxation to equilibrium and(2) proto-neutron star deleptonization. We summarize our findings and draw conclusions in Section 6. NEUTRINO/ANTINEUTRINO TRANSPORT EQUATIONS2.1.
Boltzmann equation and neutrino-matter interactions
In nuclear astrophysics applications, neutrino transport can be modeled by the Boltzmann equation. In this paper, we considera non-relativistic Boltzmann equation c ∂ t f + T ( f ) = C ( f ) , (1)which governs the particle distribution function f ( x , p , t ) that describes the density of particles at position x = ( x , x , x ) ∈ R with momentum p = ( p , p , p ) ∈ R at time t ∈ R + , where c is the speed of light. Adopting curvilinear phase-spacecoordinates, the advection operator T takes the form (see, e.g., Endeve et al. (2015)) T ( f ) := 1 √ γ (cid:88) i =1 ∂ x i (cid:0) √ γ H i x f (cid:1) + 1 √ λ (cid:88) i =1 ∂ p i (cid:0) √ λ H i p f (cid:1) , (2)where H i x f and H i p f denote the position space flux and the momentum space flux in the corresponding i th direction, respectively,and γ , λ are the determinants of the position space and momentum space metric tensors, respectively. While the form of theadvection operator in Eq. (2) holds for more general (non-orthogonal) spacetime and momentum space bases (see, e.g., Cardallet al. 2013a), we restrict ourselves to orthogonal bases in this paper. The position space flux considered in this paper takes theform H i x f = ( p i / | p | ) f , which is proportional to the particle propagation direction. We will restrict ourselves to spherical polarmomentum coordinates, and in this case, H i p f can be obtained from Eqs. (A15) and (A16) in Endeve et al. (2015). The collisionoperator C models interactions between particles and a material background, and includes emission, absorption, elastic scatteringon nucleons and nuclei, neutrino-electron scattering, and thermal pair processes from electron-positron creation and annihilation.This is the neutrino opacity set described in Bruenn (1985) (cf. Table 1 therein).In this work, we consider the transport of electron neutrinos ( ν e ) and antineutrinos ( ¯ ν e ), which results in the coupled equations c ∂ t f + T ( f ) = C ( f, ¯ f ) , (3a) c ∂ t ¯ f + T ( ¯ f ) = ¯ C ( ¯ f , f ) , (3b)where the neutrino and antineutrino distribution functions are denoted with f and ¯ f , respectively, and the collision terms C and ¯ C both depend on f and ¯ f , since they include the thermal pair production and annihilation processes of neutrino-antineutrino pairs.Before giving a detail formulation of the collision terms, we first change to spherical polar momentum coordinates ( ε, ω ) and L AIU ET AL .decompose the neutrino three-momentum as p = ε (cid:96) ( ω ) , where ε ∈ R + denotes neutrino energy and the unit vector (cid:96) ( ω ) ∈ R only depends on the angular direction ω ∈ S relative to a local orthonormal basis. The momentum space volume element isthen decomposed into d p = dV ε d ω , where dV ε := ε dε is the spherical shell energy volume element and d ω is the momentumspace angular element. With this notation, the particle distribution f ( x , p , t ) can be written as f ( x , ω , ε, t ) . At each x and t , theneutrino collision operator then takes the form (Bruenn 1985) C ( f, ¯ f ) = (cid:0) − f ( ω , ε ) (cid:1) ˆ η ( ε ) − ˆ χ ( ε ) f ( ω , ε )+ ε c ( hc ) (cid:90) S R IS ( (cid:96) · (cid:96) (cid:48) , ε ) f ( ω (cid:48) , ε ) d ω (cid:48) − ε c ( hc ) f ( ω , ε ) (cid:90) S R IS ( (cid:96) · (cid:96) (cid:48) , ε ) d ω (cid:48) + 1 c ( hc ) (cid:0) − f ( ω , ε ) (cid:1) (cid:90) R + (cid:90) S R I N ( (cid:96) · (cid:96) (cid:48) , ε, ε (cid:48) ) f ( ω (cid:48) , ε (cid:48) ) d ω (cid:48) dV ε (cid:48) − c ( hc ) f ( ω , ε ) (cid:90) R + (cid:90) S R O UT ( (cid:96) · (cid:96) (cid:48) , ε, ε (cid:48) ) (cid:0) − f ( ω (cid:48) , ε (cid:48) ) (cid:1) d ω (cid:48) dV ε (cid:48) (4) + 1 c ( hc ) (cid:0) − f ( ω , ε ) (cid:1) (cid:90) R + (cid:90) S R P R ( (cid:96) · (cid:96) (cid:48) , ε, ε (cid:48) ) (cid:0) − ¯ f ( ω (cid:48) , ε (cid:48) ) (cid:1) d ω (cid:48) dV ε (cid:48) − c ( hc ) f ( ω , ε ) (cid:90) R + (cid:90) S R A N ( (cid:96) · (cid:96) (cid:48) , ε, ε (cid:48) ) ¯ f ( ω (cid:48) , ε (cid:48) ) d ω (cid:48) dV ε (cid:48) , with h the Planck constant, ˆ η the emissivity, ˆ χ the absorption opacity, R IS the elastic (isoenergetic) scattering kernel (due toscattering with nucleons and nuclei), R I N and R O UT the neutrino-electron scattering kernels, and R P R and R A N the thermal productionand annihilation kernels due to pair processes. The antineutrino collision operator ¯ C is defined analogously with ˆ η , ˆ χ , and theopacity kernels R OP , OP = IS , I N , O UT , P R , A N , replaced by their antineutrino counterparts.For thermal emission and absorption, we follow the approach in, e.g., Burrows et al. (2006), and rewrite (cid:0) − f (cid:1) ˆ η − ˆ χ f = χ (cid:0) f − f (cid:1) , (5)where the effective opacity and the equilibrium distribution are defined as χ := ˆ η + ˆ χ , and f := ˆ η ˆ η + ˆ χ , (6)respectively. Similarly, the antineutrino emission and absorption term is written as ¯ χ (cid:0) ¯ f − ¯ f (cid:1) , with ¯ χ and ¯ f defined analogously.The equilibrium distributions f and ¯ f are given by the Fermi-Dirac distribution, which is isotropic in angle ω . Specifically, f ( ε ) = 1 e ( ε − µ ) / k B T + 1 , and ¯ f ( ε ) = 1 e ( ε − ¯ µ ) / k B T + 1 , (7)where k B is the Boltzmann constant, T is the matter temperature, µ := µ e − + ( µ p − µ n ) the neutrino chemical potential,and ¯ µ := µ e + − ( µ p − µ n ) the antineutrino chemical potential. Here, µ e − (+) is the electron (positron) chemical potential, µ n the neutron chemical potential, and µ p the proton chemical potential, which are evaluated from an appropriate EoS. Since µ e + = − µ e − , we have ¯ µ = − µ .Following Bruenn (1985), the neutrino scattering and pair process kernels are approximated with L -term Legendre expansionsin the cosine of the scattering angle α := (cid:96) · (cid:96) (cid:48) , with expansion coefficients Φ depending on the neutrino energy. Specifically, R IS ( α, ε ) ≈ L (cid:88) l =0 Φ IS l ( ε ) P l ( α ) , and R OP ( α, ε, ε (cid:48) ) ≈ L (cid:88) l =0 Φ OP l ( ε, ε (cid:48) ) P l ( α ) , (8)where OP = I N , O UT , P R , A N , P l denotes the Legendre polynomial of degree l , and the expansion coefficients of degree l are given by Φ IS l ( ε ) = 1 C l (cid:90) − R IS ( α, ε ) P l ( α ) dα , and Φ OP l ( ε, ε (cid:48) ) = 1 C l (cid:90) − R OP ( α, ε, ε (cid:48) ) P l ( α ) dα . (9)Here C l := (cid:82) − P l ( α ) dα , l = 0 , . . . , L , are normalization constants. In this work, we use P ( α ) = 1 and P ( α ) = α , thus thenormalization constants are C = 2 and C = . The antineutrino opacity kernels are approximated analogously. Due to particleconservation and detailed balance (e.g., Cernohorsky (1994)), the neutrino-electron scattering and pair process kernels satisfy R I N ( α, ε, ε (cid:48) ) = R O UT ( α, ε (cid:48) , ε ) , R O UT ( α, ε, ε (cid:48) ) = R O UT ( α, ε (cid:48) , ε ) e ( ε − ε (cid:48) ) / k B T ,R A N ( α, ε, ε (cid:48) ) = ¯ R A N ( α, ε (cid:48) , ε ) , R P R ( α, ε, ε (cid:48) ) = R A N ( α, ε, ε (cid:48) ) e − ( ε + ε (cid:48) ) / k B T . (10) OLVERS FOR NEUTRINO - MATTER COUPLING L = 0 ). More realistictreatments would include linear corrections ( L = 1 ), while further corrections ( L > ) have been shown to result in only minordifferences for neutrino-electron scattering (e.g, Smit & Cernohorsky (1996)).2.2. Angular moment equations
The need for high spatial resolution and unconstrained spatial dimensionality, e.g., to capture fluid dynamics in our targetapplications, renders direct solutions of the Boltzmann equation too expensive. However, neutrino heating rates are sensitive tothe neutrino energy distribution, which demands retention of the energy dimension of momentum space. Therefore, to balancecomputational cost with physical fidelity, we settle for solving for a finite number of angular moments of the distribution functionby adopting a two-moment model. Two-moment models are widely used to model neutrino transport in core-collapse supernovae(e.g., O’Connor (2015); Just et al. (2015); Kuroda et al. (2016); Roberts et al. (2016); Skinner et al. (2019)). In the spectraltwo-moment model, we solve for the zeroth and first moments of the neutrino distribution function, while the second momentsare obtained from a closure procedure. These moments are defined respectively as (cid:8) J , H i , K ij (cid:9) ( x , ε, t ) = 14 π (cid:90) S f ( x , ω , ε, t ) (cid:8) , (cid:96) i ( ω ) , (cid:96) i ( ω ) (cid:96) j ( ω ) (cid:9) d ω , (11)with moments for antineutrinos ( ¯ J , ¯ H i , and ¯ K ij ) defined analogously. Then, the zeroth moment J ( ¯ J ) is the spectral numberdensity of neutrinos (antineutrinos), the first moment H i ( ¯ H i ) is the spectral number flux density of neutrinos (antineutrinos),while the second moment K ij ( ¯ K ij ) is proportional to the spectral pressure tensor of neutrinos (antineutrinos). Since the neutrinodistribution function is bounded between zero and one, it can be shown that the moments must satisfy the bounds (Larecki &Banach 2011) ≤ J ≤ and | H | ≤ (cid:0) − J (cid:1) J , (12)where H := ( H , H , H ) and | H | = √H i H i , with H i associated to H k via the spatial metric γ ik , i.e., H i = γ ik H k . (Fornotational convenience in this section we sometimes use Einstein’s summation convention where repeated Latin indices implysummation from to .) Bounds equivalent to those in Eq. (12) also hold for the antineutrino moments.Taking the zeroth and first moments of Eq. (3a) leads to c ∂ t J + 14 π (cid:90) S T ( f ) d ω = 14 π (cid:90) S C ( f, ¯ f ) d ω , (13a) c ∂ t H j + 14 π (cid:90) S T ( f ) (cid:96) j d ω = 14 π (cid:90) S C ( f, ¯ f ) (cid:96) j d ω , (13b)which, after plugging in the definitions of T and C in Eqs. (2) and (4), results in the moment equations for the number densityand number flux c ∂ t J + 1 √ γ (cid:88) i =1 ∂ x i (cid:0) √ γ H i (cid:1) = η T OT − χ T OT J , (14a) c ∂ t H j + 1 √ γ (cid:88) i =1 ∂ x i (cid:0) √ γ K ij (cid:1) = 12 K ik ∂ x j γ ik − ( χ T OT + σ IS ) H j , (14b)where we have used the facts that the position space flux p / | p | f = (cid:96) ( ω ) f and that contributions from momentum space fluxesvanish due to boundary conditions. (An analogous set of moment equations is derived for antineutrinos.) Eqs. (14) hold for gen-eral, time-independent curvilinear spatial coordinates encoded in the spatial metric γ ij , including the commonly used Cartesian,spherical polar, and cylindrical coordinates. The metric tensor is used to raise and lower indices on vectors and tensors; e.g., H j = γ jk H k , K ij = γ jk K ik . Remark 1
We note that Eqs. (14) can also be obtained from the non-relativistic (i.e., zero fluid velocity and no gravitationalfields) limit of the moment equations in Shibata et al. (2011) and Cardall et al. (2013b) (see also the number conservativetwo-moment model discussed by Mezzacappa et al. (2020); their Eqs. (123) and (125)). Moreover, when including relativisticeffects, one is, among other things, confronted with choosing appropriate momentum space coordinates. The most common L AIU ET AL . (and perhaps the most natural) choice is to use momentum space coordinates in the frame of reference of the inertial observerinstantaneously comoving with the fluid (i.e., the comoving frame), as opposed to the so-called laboratory frame (see, e.g.,Mihalas & Mihalas 1999). (In the absence of fluid motion and gravitational fields, which is assumed here, there is no distinctionbetween the comoving and laboratory frames.) Specifically, the choice of comoving frame momentum coordinates provides themost straightforward framework for describing neutrino-matter interactions, and this is the choice we intend to make whenincluding relativistic effects in the future. However, this choice complicates the advection operator associated with the momentequations, which then includes Doppler and/or gravitational frequency shift terms. The inclusion of these terms is beyond thescope of the present paper, but will be considered in a future study. In the right-hand side of Eqs. (14), the elastic scattering opacity is given by σ IS = πε c ( hc ) Φ IS ( ε ) , and we define the total emissivityand total opacity as η T OT ( J , ¯ J ) = χ J + η S C ( J ) + η T P ( ¯ J ) , and χ T OT ( J , ¯ J ) = χ + χ S C ( J ) + χ T P ( ¯ J ) , (15)where J = π (cid:82) S f d ω = f . Let d ˜ V ε := 4 πε dε denote the scaled spherical shell energy volume element, then the opacityterms in Eq. (15) are defined respectively as the scattering emissivity η S C ( J ) = 1 c ( hc ) (cid:90) R + Φ I N ( ε, ε (cid:48) ) J ( ε (cid:48) ) d ˜ V ε (cid:48) , (16)the scattering opacity χ S C ( J ) = 1 c ( hc ) (cid:90) R + (cid:2) Φ I N ( ε, ε (cid:48) ) J ( ε (cid:48) ) + Φ Out ( ε, ε (cid:48) ) (cid:0) − J ( ε (cid:48) ) (cid:1) (cid:3) d ˜ V ε (cid:48) , (17)the emissivity due to thermal pair processes η T P ( ¯ J ) = 1 c ( hc ) (cid:90) R + Φ P R ( ε, ε (cid:48) ) (cid:0) − ¯ J ( ε (cid:48) ) (cid:1) d ˜ V ε (cid:48) , (18)and the opacity due to thermal pair processes χ T P ( ¯ J ) = 1 c ( hc ) (cid:90) R + (cid:2) Φ P R ( ε, ε (cid:48) ) (cid:0) − ¯ J ( ε (cid:48) ) (cid:1) + Φ A N ( ε, ε (cid:48) ) ¯ J ( ε (cid:48) ) (cid:3) d ˜ V ε (cid:48) . (19)We make the following remarks on Eqs. (14): (i) the scattering and pair processes opacities depend on J and ¯ J , respectively,due to the Fermi blocking factors; (ii) since Φ OP ≥ and J , ¯ J are between zero and one, we have σ IS , η S C , χ S C , η T P , χ T P ≥ ; and(iii) the emissivities and opacities depend on the neutrino energy ε and local matter states (e.g., density ρ , temperature T , andelectron fraction Y e ).To close the two-moment model in Eqs. (14), we adopt an algebraic closure of the form (Levermore 1984) K ij = 12 (cid:2) (cid:0) − ψ (cid:1) δ ij + (cid:0) ψ − (cid:1) (cid:98) h i (cid:98) h j (cid:3) J , (20)where (cid:98) h i = H i / | H | are components of a unit vector parallel to H , and the Eddington factor ψ = ψ ( J , h ) with h = | H | / J . Weuse the maximum entropy Eddington factor of Cernohorsky & Bludman (1994) ψ ( J , h ) = 13 + 2 (1 − J ) (1 − J )3 Θ (cid:16) h − J (cid:17) , (21)where the closure polynomial is given by Θ( x ) = 15 (3 − x + 3 x ) x . (22)We point out that this closure is based on Fermi-Dirac statistics, and is suitable for designing numerical methods for the two-moment model satisfying the bounds in Eq. (12) (e.g., Chu et al. (2019)).For the antineutrino transport equation (3b), a two-moment model on the antineutrino number density ¯ J and number flux ¯ H can be derived following similar procedure as in Eqs. (13)–(19). We then close the resulting two-moment model using a closure OLVERS FOR NEUTRINO - MATTER COUPLING M := ( J , H , ¯ J , ¯ H ) andwrite the coupled two-moment models in operator form as c ∂ t M + 1 √ γ (cid:88) i =1 ∂ x i (cid:0) √ γ F i ( M ) (cid:1) = G ( M ) + C ( M ) , (23)where F , G , and C are the position flux operator, geometry source operator, and collision operator, respectively. In particular, C ( M ) := (cid:0) η T OT − χ T OT J , ( χ T OT + σ IS ) H , ¯ η T OT − ¯ χ T OT ¯ J , ( ¯ χ T OT + ¯ σ IS ) ¯ H (cid:1) . (24)2.3. Coupling to the matter equations
Neutrino-matter interactions mediate exchange of lepton number, momentum, and energy between matter and neutrinos. As afirst approximation, we assume that the fluid remains static and momentum exchange is ignored. Under these assumptions, thematter is described by the mass density ρ ( x ) (fixed in time), temperature T ( x , t ) , and electron fraction Y e ( x , t ) , and neutrino-matter interactions result in changes to the electron fraction c dY e dt = − (cid:18) m b ρ (cid:19)(cid:18) hc ) (cid:19) (cid:90) R + (cid:0) ( η T OT − χ T OT J ) − (¯ η T OT − ¯ χ T OT ¯ J ) (cid:1) d ˜ V ε , (25)and specific internal energy c d(cid:15)dt = − (cid:18) ρ (cid:19)(cid:18) hc ) (cid:19) (cid:90) R + (cid:0) ( η T OT − χ T OT J ) + (¯ η T OT − ¯ χ T OT ¯ J ) (cid:1) ε d ˜ V ε , (26)where m b is the average baryon mass. The specific internal energy (cid:15) ( x , t ) is defined such that ρ(cid:15) is the internal energy density.We note that, given any ρ and Y e , the specific internal energy (cid:15) is an one-to-one (injective) function of the temperature T , i.e.,one can map a given T to a unique (cid:15) , and vice versa; specifically, ( ∂(cid:15)/∂T ) ρ,Y e > .Together with the number density evolution equations (cf. Eq. (14a)) for neutrinos and antineutrinos, Eqs. (25) and (26) leadto conservation of lepton number (cid:90) D (cid:20) ρ m b dY e dt + 1( hc ) ddt (cid:90) R + (cid:0) J ( ε ) − ¯ J ( ε ) (cid:1) d ˜ V ε (cid:21) d x = − (cid:73) ∂D (cid:20) c ( hc ) (cid:90) R + ( H − ¯ H ) d ˜ V ε (cid:21) d x , (27)and energy (cid:90) D (cid:20) ρ d(cid:15)dt + 1( hc ) ddt (cid:90) R + (cid:0) J ( ε ) + ¯ J ( ε ) (cid:1) ε d ˜ V ε (cid:21) d x = − (cid:73) ∂D (cid:20) c ( hc ) (cid:90) R + ( H + ¯ H ) ε d ˜ V ε (cid:21) d x , (28)respectively. Note that ( ρ/ m b ) Y e = n e , is the electron density (technically the electron minus positron density), and thatneutrinos and electrons have lepton number , while antineutrinos and positrons have lepton number − . Eq. (27) implies thatthe total lepton number in the domain D (left-hand side of Eq. (27)) only changes due to fluxes through the domain boundary ∂D (right-hand side of Eq. (27)). A similar conservation statement holds for the total energy in D as given in Eq. (28). DG-IMEX SCHEME FOR SOLVING THE MOMENT EQUATIONS3.1.
Nodal Discontinuous Galerkin (DG) space and energy discretization
We apply the nodal DG discretization (see, e.g., Hesthaven & Warburton (2008) for an overview) to Eq. (23) in space x ∈ R and energy ε ∈ R + , in which a logically Cartesian mesh with coordinate aligned elements is considered. To derive the discretizedequation from the nodal DG method, we first divide the computational domain D = D x × D ε ⊂ R × R + into a disjoint union K of open elements K , where each element takes the form K = { ( x , ε ) : x i ∈ K i := ( x i L , x i H ) , i = 1 , , , ε ∈ K ε := ( ε L , ε H ) } (29)with ∆ x i = x i H − x i L and ∆ ε = ε H − ε L denoting the side-lengths of K . For i = 1 , , , the spatial surface elements in direction x i are denoted as ˜ K i = × j (cid:54) = i K j , with space variables ˜ x i := { x j ∈ x : j (cid:54) = i } on ˜ K i , where × is the Cartesian product operator.We use V K to denote the proper volume of the element V K = (cid:90) K dV, where dV = √ γ d x dV ε = √ γ d x ε dε . (30)0 L AIU ET AL .We let the approximation space for the DG method, V N x ,N ε , be constructed from the tensor product of one-dimensional poly-nomials of maximal degrees N x and N ε in space and energy, respectively. Note that functions in V N x ,N ε can be discontinuousacross element interfaces. The semi-discrete DG problem is to find M h ∈ V N x ,N ε (which approximates M in Eq. (23)) suchthat (cf. Cockburn & Shu (2001)) c ∂ t (cid:90) K M h v dV + (cid:88) i =1 (cid:90) K ε (cid:90) ˜ K i (cid:0) √ γ (cid:98) F i ( M h ) v (cid:12)(cid:12) x i H − √ γ (cid:98) F i ( M h ) v (cid:12)(cid:12) x i L (cid:1) d ˜ x i dV ε − (cid:88) i =1 (cid:90) K (cid:0) F i ( M h ) ∂v∂x i (cid:1) dV = (cid:90) K G ( M h ) v dV + (cid:90) K C ( M h ) v dV, (31)for all v ∈ V N x ,N ε and all K ∈ K . Here the numerical flux approximating the flux on the spatial surface element ˜ K i is denotedas (cid:98) F i ( M h ) . In this work, we consider the Lax-Friedrichs (LF) flux (cid:98) F i ( M h ) | x i = f LF ,i ( M h | x i, − , M h | x i, + ):= 12 (cid:0) F i ( M h | x i, − ) + F i ( M h | x i, + ) − α i ( M h | x i, + − M h | x i, − ) (cid:1) , (32)where M h | x i, ± = lim δ → M h | x i ± δ are the evaluations of M h at the immediate right/left of x i , which thus are functions of (˜ x i , ε, t ) . The parameter α i = || eig (cid:0) ∂ F i /∂ M (cid:1) || ∞ is the largest eigenvalue of the flux Jacobian. For massless neutrinos, whichpropagate at the speed of light, we can take α i = c (i.e., the global LF flux).In each element K , we approximate the conserved variables M by M ( x , ε, t ) ≈ M h ( x , ε, t ) = N ε (cid:88) n ε =1 N x (cid:88) n x = M n x ,n ε ( t ) (cid:96) n x ( x ) (cid:96) n ε ( ε ) , (33)where { (cid:96) n ε } N ε n ε =1 is chosen to be a collection of Lagrange polynomials in energy of degree N ε , n x := { n , n , n } is a multi-index that goes from = { , , } to N x = { N x , N x , N x } , and { (cid:96) n x } N x n x = is the collection of multidimensional polynomialsdefined as (cid:96) n x ( x ) := (cid:81) i =1 (cid:96) n i ( x i ) , with { (cid:96) n i } N x n i =1 one-dimensional Lagrange polynomials in x i of degree N x . It then followsfrom these definitions that { (cid:96) n x ( x ) (cid:96) n ε ( ε ) } N x ,N ε n x = ,n ε =1 forms a basis of V N x ,N ε on K . Motivated by the numerical experimentsreported in Bassi et al. (2013), we consider here the Lagrange polynomials with Gauss-Legendre interpolation points (instead ofGauss-Legendre-Lobatto points), e.g., (cid:96) n ε on interval K ε is defined as (cid:96) n ε ( ε ) = (cid:96) loc N ε ,n ε ( ξ ) := N ε (cid:89) j =1 j (cid:54) = n ε ξ − ξ εj ξ εn ε − ξ εj , ∀ ξ ∈ I := [ − / , / , n ε = 1 , . . . , N ε . (34)The local variable ξ ( ε ) := ε − ε C ∆ ε with the center ε C := ( ε L + ε H ) / , and interpolation points { ξ εj } N ε j =1 are given by the N ε -pointGauss-Legendre quadrature abscissas on I . On interval K i , (cid:96) n i is defined analogously as (cid:96) n i ( x i ) = (cid:96) loc N x ,n i ( ξ ) , with Gauss-Legendre interpolation points { ξ ij } N x j =1 on I . With this choice of basis, it follows that the expansion in Eq. (33) becomes a nodalrepresentation of M h , i.e., M h ( x n x , ε n ε , t ) = M n x ,n ε ( t ) , where x n x = { x n , x n , x n } and ε n ε are the global space andenergy variables corresponding to the local interpolation points { ξ ij } N x j =1 and { ξ εj } N ε j =1 on element K , respectively. Figure 2 showsan example of nodal DG elements in a reduced space x ∈ R and energy ε ∈ R + with the interpolations points in both local andglobal coordinates.We then follow the standard practice and approximate Eq. (31) by a semi-discrete system consisting of ( N x ) N ε equations,each on a nodal value M n x ,n ε , n x = , . . . , N x , n ε = 1 , . . . , N ε . To derive these equations, we approximate the integralsin Eq. (31) using the N x -points and N ε -points Gauss-Legendre quadrature rules in the space and energy, respectively, with theassociated weights { w n i } N x n i =1 and { w n ε } N ε n ε =1 normalized such that (cid:80) N x n i =1 w n i = 1 and (cid:80) N ε n ε =1 w n ε = 1 . Specifically, for v ( x , ε ) = (cid:96) n x ( x ) (cid:96) n ε ( ε ) , we have ∂ t (cid:90) K M h v dV ≈ w n x ,n ε ε n ε √ γ n x | K | ∂ t M n x ,n ε , (35) OLVERS FOR NEUTRINO - MATTER COUPLING x N x x x x x e x N e e K [ — - 1 — , ][ — — , ] (a) Local coordinate xx N 𝜀 𝜀 xx x S x S 𝜺 |K 𝜺 | |K x | N x e 𝜀 𝜀 e (b) Global coordinate Figure 2.
Example of nodal DG elements in a computational domain in R × R + – Figure 2a shows the interpolation points { ξ ij } N x j =1 and { ξ εj } N ε j =1 in a nodal DG element K in the local coordinate. Figure 2b shows the union K of all elements K and the space and energy nodes –the values of interpolation points in the global coordinate. Figure 2b also illustrates the collection of space nodes S x and the collection energynodes S ε , which are defined in Eqs. (40)–(41). and (cid:90) K (cid:0) G ( M h ) + C ( M h ) (cid:1) v dV ≈ w n x ,n ε ε n ε √ γ n x | K | (cid:0) G ( M h ) n x ,n ε + C ( M h ) n x ,n ε (cid:1) , (36)where | K | = ( (cid:81) i =1 ∆ x i )∆ ε , w n x ,n ε = ( (cid:81) i =1 w n i ) w n ε , and γ n x denotes the value of γ at x n x . For the surface integrals andthe ‘volume terms’, we denote the index ˜ n i := { n j ∈ n x : j (cid:54) = i } and obtain, e.g., (cid:90) K ε (cid:90) ˜ K i √ γ (cid:98) F i ( M h ) v | x i H d ˜ x i dV ε ≈ w ˜ n i ,n ε ε n ε (cid:0) √ γ ˜ n i (cid:98) F i ˜ n i ,n ε (cid:96) n i (cid:1)(cid:12)(cid:12) x i H | ˜ K i | | K ε | , (37)and (cid:90) K F i ( M h ) ∂v∂x i dV ≈ w ˜ n i ,n ε ε n ε N x (cid:88) k i =1 w k i (cid:0) √ γ ˜ n i F i ˜ n i ,n ε ∂(cid:96) k i ∂x i (cid:1)(cid:12)(cid:12) x iki | ˜ K i | | K ε | , (38)where | K ε | = ∆ ε , | ˜ K i | = (cid:81) j (cid:54) = i ∆ x j , w ˜ n i ,n ε = ( (cid:81) j (cid:54) = i w n j ) w n ε , γ ˜ n i is the evaluations of γ at ˜ x i ˜ n i , and (cid:98) F i ˜ n i ,n ε and F i ˜ n i ,n ε are the evaluations of (cid:98) F i ( M h ) and F i ( M h ) at (˜ x i ˜ n i , ε n ε ) , respectively. Here γ ˜ n i , (cid:98) F i ˜ n i ,n ε , and F i ˜ n i ,n ε are functions of x i .Plugging the terms in Eqs. (35)–(38) into Eq. (31) and dividing through by w n x ,n ε √ γ n x ε n ε | K | leads to the semi-discreteform of the moment equations when the test function is v = (cid:96) n x (cid:96) n ε : c ∂ t M n x ,n ε = − √ γ n x (cid:88) i =1 w n i ∆ x i (cid:0) √ γ ˜ n i (cid:98) F i ˜ n i ,n ε (cid:96) n i (cid:1)(cid:12)(cid:12) x i H x i L (39) + 1 √ γ n x (cid:88) i =1 w n i ∆ x i N x (cid:88) k i =1 w k i (cid:0) √ γ ˜ n i F i ˜ n i ,n ε ∂(cid:96) k i ∂x i (cid:1)(cid:12)(cid:12) x iki + G ( M h ) n x ,n ε + C ( M h ) n x ,n ε , for n x = , . . . , N x and n ε = 1 , . . . , N ε in all K ∈ K . Eq. (39) defines the spatial and energy discretization of the momentequations and provides the basis for implementation in thornado . Before discussing the time discretization, we further simplify2 L AIU ET AL .Eq. (39) utilizing the structure of the collision operator C . Here we introduce the collection of space and energy nodes from allelements, denoted as S := ∪ K ∈K { ( x n x , ε n ε ) ∈ K : n x = , . . . , N x , n ε = 1 , . . . , N ε } , (40)and we define the spatial component of S as S x and the energy component of S as S ε , i.e., x ∈ S x , ε ∈ S ε if and only if ( x , ε ) ∈ S . (41)A simplified illustration of S x and S ε is given in Figure 2b, in which x ∈ R . At each x k ∈ S x , the nodal value of M h on all ε ∈ S ε is then denoted as M k ( t ) := { M h ( x k , ε q , t ) : ε q ∈ S ε } . With these notations, we write Eq. (39) in the operator form inthe remainder of this paper as c ∂ t M k = F M ( M h ) k + G M ( M h ) k + C M ( M k ) , ∀ k such that x k ∈ S x , (42)where F M , G M , and C M denote respectively the discrete position space flux operator, the discrete geometry source, and the discretecollision operator. Here the collision term C M ( M k ) depends only on M k instead of the full discretized solution M h , since thephysical interactions modeled in the collision operator (see Eq. (4)) are independent of the position x , while coupled in the energydomain. Specifically, let ( J k , H k , ¯ J k , ¯ H k ) denote the discretized moments M k , then it follows from Eq. (24) that C M ( M k ) := ( η T OT − χ T OT J k , ( χ T OT + σ IS ) H k , ¯ η T OT − ¯ χ T OT ¯ J k , (¯ χ T OT + ¯ σ IS ) ¯ H k ) , (43)where σ IS and ¯ σ IS are the values of σ IS , ¯ σ IS on the energy nodes S ε , respectively, and η T OT , χ T OT , ¯ η T OT , and ¯ χ T OT are the discretecounterparts of η T OT , χ T OT , ¯ η T OT , and ¯ χ T OT . These discrete opacities are given by replacing the energy integrals in Eqs. (15)–(19)with the numerical integrals using S ε as quadrature points, for example, η S C ( J ) = 1 c ( hc ) (cid:90) R + Φ I N ( ε, ε (cid:48) ) J ( ε (cid:48) ) d ˜ V ε (cid:48) ≈ c ( hc ) | S ε | (cid:88) q =1 w εq Φ I N ( ε, ε q ) J q =: η S C ( J ) , (44)where w εq denotes the weight associated to energy node ε q , and | S ε | denotes the total number of energy nodes in S ε . Specifically,for ε q ∈ S ε in some energy element K ε , w εq := 4 πε q w n ε | K ε | , (45)where w n ε is the local Gauss-Legendre weight at ε q in K ε defined earlier in this section.Finally, we apply this same nodal DG discretization on the electron fraction and specific internal energy evolution equations(25)–(26). Augmenting Eq. (42) to the resulting semidiscrete equations then gives c ∂ t V k = F ( V h ) k + G ( V h ) k + C ( V k ) , ∀ k such that x k ∈ S x , (46)where V := ( Y e , (cid:15), M ) , V h := ( Y e, h , (cid:15) h , M h ) is the fully discretized version of V , and V k denotes the nodal value of V h at point x k ∈ S x . Here the operators F := ( , , F M ) , G := ( , , G M ) , and C is given by the discrete version of the right-hand sides inEqs. (25)–(26) and C M . Specifically, in C , the energy integrals in Eqs. (25)–(26) are evaluated using a quadrature with the energynodes ε q ∈ S ε as abscissas and weights defined in Eq. (45). Remark 2
In multidimensional CCSN simulations, the use of curvilinear coordinates, e.g., spherical-polar coordinates, suffersfrom an excessively stringent Courant-Friedrichs-Lewy (CFL) time step restriction due to singularities at the origin and the poles(see, e.g., M¨uller (2020, Section 2.1.1) and references therein). Although these singularities are not an issue in the sphericallysymmetric CCSN models considered here, they will become a problem when extending the nodal DG scheme to multidimensionalCCSN models. Several techniques have been developed to address this issue, such as mesh coarsening (Skinner et al. 2019),element averaging/merging (Asaithambi & Mahesh 2017; M¨uller et al. 2019), and spectral filtering (M¨uller et al. 2019). Forfuture multidimensional simulations, we will consider (i) adopting one of the aforementioned techniques or (ii) using Cartesiancoordinates in combination with adaptive mesh refinement.
OLVERS FOR NEUTRINO - MATTER COUPLING
Implicit-explicit time integration scheme
An implicit-explicit (IMEX) time integration scheme (Ascher et al. 1997; Pareschi & Russo 2005) is considered here forsolving the two-moment model in Eq. (23). When applied to transport equations with collision terms, IMEX schemes usuallyhandle the collision term with an implicit method, while applying an explicit method on the advection term (Hu et al. 2018) (seealso O’Connor (2015); Just et al. (2015); Kuroda et al. (2016); Skinner et al. (2019) for applications of IMEX-type schemes toneutrino transport). This approach relaxes the excessive time-step restriction from an explicit and stiff collision term, and avoidsthe spatially coupled nonlinear solves from an implicit advection term. While the class of IMEX schemes considered in thispaper is detailed in Chu et al. (2019), we include it in the following paragraph for completeness. We also stress that even thoughthe nonlinear solution strategies given in Section 4 are motivated from the implicit part of this class of IMEX schemes, they aregeneral enough to be used with any IMEX scheme that treats the collision term implicitly.To perform time integration, we discretize the time interval [ t , t f ] into N time steps t = t < t < · · · < t N = t f and denote V k ( t n ) as V n k , n = 1 , . . . , N . At each spatial node x k and time t n , the IMEX scheme integrates the semi-discrete equation (46)in time via V (0) k = V n k , (47a) V ( i ) k = i − (cid:88) j =0 c ij V ( ij ) k + a ii c ∆ t C ( V ( i ) k ) , i = 1 , . . . , s, (47b) V n +1 k = V ( s ) k , (47c)where s is the number of stages, the parameters c ij ≥ , (cid:80) i − j =0 c ij = 1 , a ii > , and V ( ij ) k is given by an explicit update V ( ij ) k = V ( j ) k + ˆ c ij c ∆ t (cid:0) F ( V ( j ) h ) k + G ( V ( j ) h ) k (cid:1) (48)with parameters ˆ c ij ≥ . Thus, in each time step, the s -stage IMEX scheme requires s evaluations of the discrete flux andgeometry operators F and G , and s inversions of the discrete collision operator C . Here the number of evaluations of F and G isidentical to the number of stages due to the fact that, by reusing values of F and G from earlier stages, each stage only requiresone additional evaluation of F and G . In general, the inversion of C is the dominant cost in the IMEX scheme. In the remainderof this section, we present the details of the nonlinear system arising from inverting C .At each stage i of the IMEX scheme, Eq. (47b) can be considered as the nonlinear system V (+) k = V ( ∗ ) k + τ C ( V (+) k ) , (49)where τ := a ii c ∆ t > denotes the effective time step, V ( ∗ ) k denotes the weighted sum of explicit updates V ( ij ) k , and V (+) k denotes the unknown nodal values to be solved.Let ( Y e, k , (cid:15) k ) denote the nodal value of ( Y e , (cid:15) ) at x k , and let ( J k , H k , ¯ J k , ¯ H k ) denote the discrete moment M k , which collectsvalues of moments ( J , H , ¯ J , ¯ H ) at x k on all energy nodes in S ε . Eq. (49) can then be considered as a nonlinear system on V k := ( Y e, k , (cid:15) k , J k , H k , ¯ J k , ¯ H k ) . For notational simplicity, we suppress all k subscripts when denoting the nodal values ofelectron fraction, specific internal energy, and moments in the remainder of this paper. It follows from the definition of C that, ateach x k , the resulting nonlinear system is given by Y (+) e = Y ( ∗ ) e − τ m b ρ (cid:80) | S ε | q =1 w (2) q (cid:0) ( η T OT − χ T OT J (+) q ) − (¯ η T OT − ¯ χ T OT ¯ J (+) q ) (cid:1) , (50a) (cid:15) (+) = (cid:15) ( ∗ ) − τ ρ (cid:80) | S ε | q =1 w (3) q (cid:0) ( η T OT − χ T OT J (+) q ) + (¯ η T OT − ¯ χ T OT ¯ J (+) q ) (cid:1) , (50b) J (+) = J ( ∗ ) + τ (cid:0) η T OT − χ T OT J (+) (cid:1) , (50c) ¯ J (+) = ¯ J ( ∗ ) + τ (cid:0) ¯ η T OT − ¯ χ T OT ¯ J (+) (cid:1) , (50d) H (+) = H ( ∗ ) − τ ( χ T OT + σ IS ) H (+) , (50e) ¯ H (+) = ¯ H ( ∗ ) − τ (¯ χ T OT + ¯ σ IS ) ¯ H (+) . (50f)In Eqs. (50a) and (50b), a quadrature is used to evaluate the energy integrals in Eqs. (25) and (26), as described when defining C in Eq. (46). Here ( J q , ¯ J q ) denotes the value of ( J , ¯ J ) at energy node ε q ∈ S ε . The physical constants are absorbed into theweights, i.e., w (2) q := 1( hc ) w εq , and w (3) q := 1( hc ) ε q w εq , (51)4 L AIU ET AL .with w εq defined in Eq. (45).We take a two-step approach to solve Eq. (50), which first solves the fully coupled nonlinear system in Eqs. (50a)–(50d), plugthe solution ( Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) ) into Eqs. (50e)–(50f) to compute ( χ T OT , σ IS , ¯ χ T OT , ¯ σ IS ) , and then update ( H (+) , ¯ H (+) ) . We notethat, once ( χ T OT , σ IS , ¯ χ T OT , ¯ σ IS ) are known, solving Eqs. (50e)–(50f) is straightforward. Thus, we focus on the solution procedure ofthe coupled system in Eqs. (50a)–(50d), where the opacities ( η T OT , χ T OT , ¯ η T OT , ¯ χ T OT ) are functions of ( Y e , (cid:15) , J , ¯ J ) . Specifically, whilethe opacities are written explicitly as functions of ( J , ¯ J ) in Eqs. (15)–(19), they also depend on the matter state ( Y e , (cid:15) ) throughthe opacity kernels Φ I N , Φ O UT , Φ P R , and Φ A N in Eqs. (16)–(19). In a fully implicit approach, these opacities need to be updated inthe solution procedure of Eqs. (50a)–(50d) in order to remain consistent with ( Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) ) .We proceed to investigate various approaches to solve the coupled nonlinear system in Eqs. (50a)–(50d). NONLINEAR SOLUTION STRATEGIESIn this section, we discuss two approaches for solving the system in Eqs. (50a)–(50d), which couples the evolution of the matterstates ( Y e , (cid:15) ) to the neutrino and antineutrino spectral distributions ( J , ¯ J ) . To start, we first rewrite Eqs. (50a)–(50d) as Y (+) e = Y ( ∗ ) e − m b ρ (cid:80) | S ε | q =1 w (2) q (cid:0) J (+) q − ¯ J (+) q (cid:1) + m b ρ (cid:80) | S ε | q =1 w (2) q (cid:0) J ( ∗ ) q − ¯ J ( ∗ ) q (cid:1) , (52a) (cid:15) (+) = (cid:15) ( ∗ ) − ρ (cid:80) | S ε | q =1 w (3) q (cid:0) J (+) q + ¯ J (+) q (cid:1) + ρ (cid:80) | S ε | q =1 w (3) q (cid:0) J ( ∗ ) q + ¯ J ( ∗ ) q (cid:1) , (52b) J (+) = J ( ∗ ) + τ η T OT (cid:0) Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) (cid:1) − τ χ T OT (cid:0) Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) (cid:1) J (+) , (52c) ¯ J (+) = ¯ J ( ∗ ) + τ ¯ η T OT (cid:0) Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) (cid:1) − τ ¯ χ T OT (cid:0) Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) (cid:1) ¯ J (+) , (52d)with unknowns Y (+) e ∈ R , (cid:15) (+) ∈ R , J (+) ∈ R | S ε | , and ¯ J (+) ∈ R | S ε | . Here Eqs. (52a)–(52b) are derived by substitutingEqs. (50c)–(50d) into the right-hand sides of Eqs. (50a)–(50b) to remove the explicit dependency on opacities, and Eqs. (52c)–(52d) are identical to Eqs. (50c)–(50d), with explicit expression of the dependency of opacities ( η T OT , χ T OT , ¯ η T OT , ¯ χ T OT ) on matterstates (through opacity kernels) and neutrino (antineutrino) distributions. Specifically, η T OT ( Y e , (cid:15) , J , ¯ J ) = χ J + η S C ( J ) + η T P (¯ J ) , and χ T OT ( Y e , (cid:15) , J , ¯ J ) = χ + χ S C ( J ) + χ T P (¯ J ) , (53)where the discrete opacities χ , χ OP , and η OP are computed as in Eq. (44), using opacity kernels Φ OP evaluated at ( Y e , (cid:15) ) , and theFermi-Dirac distribution J evaluated at chemical potential µ ( Y e , (cid:15) ) and matter temperature T ( Y e , (cid:15) ) .We propose two approaches for solving the nonlinear system in Eq. (52) – a coupled approach and a nested approach. The for-mer directly considers Eq. (52) as a fully coupled system, while the latter formulates Eq. (52) as a nested system with Eqs. (52a)–(52b) in the outer layer and Eqs. (52c)–(52d) in the inner. Opacity kernel evaluations, i.e., evaluating Φ OP at given ( Y e , (cid:15) ) , areneeded when solving Eqs. (52a)–(52b). Since the tabulated opacity kernels are used, evaluating Φ OP requires opacity table inter-polations, which are the dominant cost in solving Eq. (52). The nested approach aims to reduce the number of opacity kernelevaluations by giving a better prediction on ( J , ¯ J ) through the inner solver on Eqs. (52c)–(52d).We consider a fixed-point iteration method with Anderson acceleration and Newton’s method as the nonlinear system solvers inboth the coupled and nested approach. When solving the systems considered here, fixed-point methods are often more attractivethan Newton’s method because they (1) do not require the Jacobian matrix, which can be difficult to compute accurately withtabulated opacities; and (2) avoid inversion of dense linear systems. However, the rate of convergence can be slower for fixed-point methods than that of Newton-based methods. The performance of these two types of solvers on systems arising from eachapproach is compared in the numerical results reported in Section 5. In the following subsections, we state the coupled fixed-pointalgorithm (section 4.1), the coupled Newton’s method (section 4.2), the nested fixed-point algorithm (section 4.3), and the nestedNewton’s method (section 4.4). 4.1. Coupled fixed-point algorithm
To simplify the notation, we denote the matter states as u := ( Y (+) e , (cid:15) (+) ) , the discretized neutrino and antineutrino distribu-tions as U := ( J (+) , ¯ J (+) ) , and the collection of all the unknowns as U := ( u , U ) in the remainder of the paper. To formulatethe system in Eq. (52) as a fixed-point problem, we write it as U = G ( U ) := (cid:32) g ( U ) G ( u , U ) (cid:33) , (54)where g ( U ) := (cid:32) − m b ρ (cid:80) | S ε | q =1 w (2) q (cid:0) J (+) q − ¯ J (+) q (cid:1) + C Y e − ρ (cid:80) | S ε | q =1 w (3) q (cid:0) J (+) q + ¯ J (+) q (cid:1) + C (cid:15) (cid:33) (55) OLVERS FOR NEUTRINO - MATTER COUPLING C Y e = Y ( ∗ ) e + m b ρ (cid:80) | S ε | q =1 w (2) q (cid:0) J ( ∗ ) q − ¯ J ( ∗ ) q (cid:1) , C (cid:15) = (cid:15) ( ∗ ) + ρ (cid:80) | S ε | q =1 w (3) q ( J ( ∗ ) q + ¯ J ( ∗ ) q (cid:1) , and G ( u , U ) := (cid:32) (cid:0) J ( ∗ ) + τ η T OT ( u , U ) (cid:1) / (cid:0) τ χ T OT ( u , U ) (cid:1)(cid:0) ¯ J ( ∗ ) + τ ¯ η T OT ( u , U ) (cid:1) / (cid:0) τ ¯ χ T OT ( u , U ) (cid:1) (cid:33) . (56)Here G is an equivalent form of Eqs. (52c)–(52d), which ensures that G is a contraction map, i.e., the Lipschitz constant of G isstrictly less than one.The coupled fixed-point algorithm considers Eq. (54) as a fixed-point problem with unknowns U , e.g., applying Picard iterationon Eq. (54) leads to U [ k +1] = G ( U [ k ] ) , (57)where U [ k ] denotes the k -th iterate of unknowns U , starting from an initial guess U [0] . When G is a contraction mapping, Picarditeration guarantees that, as k → ∞ , the iterate U [ k ] converges to U = ( Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) ) , the solution to Eq. (52). Herethe opacities ( η , χ , ¯ η , ¯ χ ) in G are updated at each iteration k using U [ k ] and thus are consistent with the solution. While thePicard iteration guarantees convergence when G is a contraction, the convergence could be slow. To achieve faster convergence,we implement Anderson acceleration (Anderson 1965; Walker & Ni 2011) to solve Eq. (54). Anderson acceleration utilizesinformation from previous iterations to update the unknowns, which is expected to give faster convergence than Picard iteration,but at a cost of additional memory usage. Specifically, in iteration k , Anderson acceleration on the coupled problem first computesthe residual r [ k ] := G ( U [ k ] ) − U [ k ] , (58)then solves a least-squares problem α ∗ := argmin α ∈ R mk +1 (cid:8)(cid:13)(cid:13) (cid:80) m k i =0 α i r [ k − i ] (cid:13)(cid:13) : (cid:80) m k i =0 α i = 1 (cid:9) with m k := min { m, k } , and finallyupdates U [ k +1] = (cid:88) m k i =0 α ∗ i G ( U [ k − i ] ) . (59)Here the truncation parameter m ≥ is an integer that indicates the “memory” of Anderson acceleration, i.e., the maximumnumber of residuals kept in memory. When m = 0 , the solver reduces to Picard iteration. For m > , Anderson accelerationupdates U using a linear combination of the last m k iterates that leads to the minimum residual. In the numerical tests inSection 5, we use m = 2 , which we have found to significantly reduce the number of iterations when compared to Picard. For m = 2 , the additional memory required for Anderson acceleration is small since each implicit solve is local in space. In addition,the least-squares problem for α ∗ i is small, and can be written out explicitly or solved using LAPACK’s DGELS.4.2. Coupled Newton’s method
The other solver we considered for the nonlinear coupled system in Eq. (52) is Newton’s method, which formulates Eq. (52)as a root-finding problem (cid:32) f ( u , U ) F ( u , U ) (cid:33) =: F ( U ) = , (60)where f ( u , U ) := (cid:32) f Y e ( u , U ) f (cid:15) ( u , U ) (cid:33) := (cid:32) Y (+) e + m b ρ (cid:80) | S ε | q =1 w (2) q (cid:0) J (+) q − ¯ J (+) q (cid:1) − C Y e (cid:15) (+) + ρ (cid:80) | S ε | q =1 w (3) q (cid:0) J (+) q + ¯ J (+) q (cid:1) − C (cid:15) (cid:33) (61)with constants C Y e and C (cid:15) defined as in Eq. (55), and F ( u , U ) := (cid:32) F J ( u , U ) F ¯ J ( u , U ) (cid:33) := (cid:32) (cid:0) τ χ T OT ( u , U ) (cid:1) J (+) − (cid:0) J ( ∗ ) + τ η T OT ( u , U ) (cid:1)(cid:0) τ ¯ χ T OT ( u , U ) (cid:1) ¯ J (+) − (cid:0) ¯ J ( ∗ ) + τ ¯ η T OT ( u , U ) (cid:1) (cid:33) . (62)Applying Newton’s method to solve Eq. (60) leads to the following update of U in iteration k , U [ k +1] = U [ k ] + δ U [ k ] , (63)where the Newton step is given by δ U [ k ] = − J ( U [ k ] ) − F ( U [ k ] ) , with J ( U [ k ] ) := F (cid:48) ( U [ k ] ) . (64)6 L AIU ET AL .Eq. (64) shows that two key components are needed in Newton’s method – evaluating the Jacobian J and solving the linearsystem for the Newton step. For the coupled problem in Eq. (60), Jacobian evaluation at a given ˆ U requires computing ∂ f ∂ u , ∂ f ∂ U , ∂ F ∂ u , and ∂ F ∂ U at (ˆ u , ˆ U ) =: ˆ U . From Eq. (61), it is clear that evaluating ∂ f ∂ u and ∂ f ∂ U at (ˆ u , ˆ U ) is straightforward with minimal cost.However, it follows from Eq. (62) that evaluating ∂ F ∂ u and ∂ F ∂ U at (ˆ u , ˆ U ) requires gradients of the opacities ( η T OT , χ T OT , ¯ η T OT , ¯ χ T OT ) with respect to the matter state u and the neutrino and antineutrino number densities U , respectively. In particular, Eqs. (44) and(15)–(19) imply that, ∂ Φ OP ∂ u , which is the gradient of opacity kernels with respect to the matter state, is involved in the computationof ∂ F ∂ u (ˆ u , ˆ U ) . As discussed earlier, tabulated opacity kernels are considered in this paper. Thus, we can only obtain approximate ∂ Φ OP ∂ u from the tabulated quantities. Further, when the opacity kernels are not tabulated in terms of u , the gradient ∂ Φ OP ∂ u has tobe approximated using the chain rule. The detailed calculations of these kernel derivatives are given in Appendix A. Once theapproximate Jacobian is obtained, the linear solve in Eq. (64) is performed via LAPACK’s DGESV.4.3. Nested fixed-point algorithm
The next approach we consider is a nested algorithm, which formulates Eq. (54) as a nested fixed-point problem with twolayers u = g ( ˆ U ( u )) , (65a) ˆ U ( u ) = G ( u , ˆ U ( u )) , (65b)where the outer layer, Eq. (65a), is a fixed-point problem on the matter states u , and the inner, Eq. (65b), is on the distributions ˆ U for fixed matter states u . These two problems are nested in the sense that evaluating the right-hand side of Eq. (65a) at a given u requires solving Eq. (65b). For example, applying Picard iteration on both Eqs. (65a) and (65b) gives the following iterativescheme u [ k +1] = g ( ˆ U ( u [ k ] )) , (66a)where ˆ U ( u [ k ] ) = U [ k, ∗ ] , the limit point of the inner Picard iteration U [ k,(cid:96) +1] = G ( u [ k ] , U [ k,(cid:96) ] ) . (66b)In practice, we use Anderson acceleration with m = 2 , as described in Section 4.1, to accelerate both the outer and inner solvesseparately.The nested approach was considered in Laiu et al. (2020) for relaxing the nonlinear coupling between the electric field andelectron concentration when solving implicit systems for semiconductor models. Here the nested approach is motivated by thefact that in solving Eq. (54), the most costly part is evaluating the opacity kernels Φ at a given matter state u , which is performedin G whenever u is updated. Therefore, while the coupled approach seems simple and straightforward, the nested structure inEq. (65) justifies the additional complexity by reducing the number of updates (on u ) in Eq. (65a) via a more accurate distributionupdate given by solving Eq. (65b) at the current matter state. Note that the matter state u is fixed in the solution procedure of theinner problem in Eq. (65b), which does not require opacity kernel evaluations and results in much cheaper inner iterations.4.4. Nested Newton’s method
The nested Newton’s method formulates the inner layer of the nested system in Eq. (65) as a root-finding problem, resulting in u = g ( ˆ U ( u )) , (67a) F ( u , ˆ U ( u )) = , (67b)where the outer layer in Eq. (67a) is still a fixed-point problem on u , and the inner layer in Eq. (67b) is a root-finding problem on ˆ U for fixed u , with F defined in Eq. (62). Here Eq. (67a) is solved using Anderson acceleration, and, whenever the right-handside of Eq. (67a) is evaluated at some given u , the inner problem in Eq. (67b) is solved via Newton’s method to obtain ˆ U ( u ) that is used to evaluate g ( ˆ U ( u )) . This nested solver is identical to the nested fixed-point algorithm in Section 4.3, except thatthe inner problem is solved via Newton’s method instead of Anderson acceleration. Specifically, let u [ k ] be the k th iterate in theouter layer, then ˆ U ( u [ k ] ) = U [ k, ∗ ] , which is the limit point of the Newton iterate U [ k,(cid:96) +1] = U [ k,(cid:96) ] + δ U [ k,(cid:96) ] , with δ U [ k,(cid:96) ] = − (cid:2) ∂ F ∂ U ( u [ k ] , U [ k,(cid:96) ] ) (cid:3) − F ( u [ k ] , U [ k,(cid:96) ] ) . (68) OLVERS FOR NEUTRINO - MATTER COUPLING u is fixed in the inner iterations, Eqs. (44) and (15)–(19) imply that the Jacobian ∂ F ∂ U can be calculated withno additional opacity kernel evaluation, which justifies this nested approach. The reason we choose not to formulate the outerlayer as a root-finding problem and solve it with Newton’s method is to avoid the costly opacity kernel gradient approximationdiscussed in Section 4.2. NUMERICAL EXPERIMENTSThe four iterative solvers introduced in Section 4 are compared in this section. First, to investigate the iterative solvers inisolation, we report on results obtained on relaxation problems under conditions expected in CCSNe. Then we compare theiterative solvers in the context of the IMEX scheme in Eqs. (47)-(48) on proto-neutron star deleptonization problems using matterconditions from spherically symmetric CCSN simulations at various times after core-bounce.5.1.
Implementation details
In the numerical tests discussed in this section, the DG-IMEX scheme and the nonlinear solvers are implemented following thespecifics below, unless otherwise noted.• DG scheme – We consider problems with one spatial dimension (imposing spherical symmetry). For the relaxation prob-lems in Section 5.2, we solve the space-homogeneous problem in Eq. (49) for a single spatial element. For the proto-neutronstar deleptonization problems in Section 5.3, the spatial domain r ∈ [0 , km is divided into geometrically progress-ing elements with the first element of size ∆ x = 1 km and the last element of size ∆ x ≈ . km. In both tests, the energydomain covering ε ∈ [0 , MeV is divided into geometrically progressing elements, where the first element has ∆ ε = 4 MeV and the last element has ∆ ε ≈ MeV. The spatial and energy DG elements considered here are linear( k = 1 ).• IMEX scheme – In the proto-neutron star deleptonization problem in Section 5.3 we use the IMEX scheme in Eqs. (47)–(48) with two stages ( s = 2 ). When written in the so-called Shu-Osher form (as in Eqs. (47)–(48)), the coefficients aregiven by Chu et al. (2019) (cid:34) c c c (cid:35) = (cid:34) / / (cid:35) , (cid:34) ˆ c ˆ c ˆ c (cid:35) = (cid:34)
10 1 (cid:35) , and (cid:34) a a (cid:35) = (cid:34) / (cid:35) . (69)This scheme consists of two evaluations of the explicit part and two implicit solves. We use the realizability-enforcinglimiter in Chu et al. (2019) to enforce realizable moments (cf. Eq. (12)) after each stage.• Equation of state and neutrino opacity tables – In all the tests we use a tabulated version of the SFHo EoS (Steineret al. 2013). Thermodynamic (dependent) variables are tabulated as a function of mass density, temperature, and electronfraction ( ρ , T , and Y e ). The EoS table covers the ranges ρ ∈ [1 . × , . × ] g cm − , using N ρ = 185 points(logarithmically spaced to achieve about points per decade), T ∈ [1 . × , . × ] K, using N T = 81 points(logarithmically spaced to achieve about points per decade), and Y e ∈ [0 . , . , using N Y e = 30 points (linearlyspaced). The neutrino opacities are taken from Bruenn (1985), with all the input thermodynamic quantities computedwith the SFHo EoS. The absorption and scattering opacities ( χ and σ IS ) are tabulated in terms of the neutrino energy ε inaddition to ρ , T , and Y e (using the same resolution as the EoS table). The neutrino energy range covers ε ∈ [0 . , MeV,using N ε = 40 logarithmically spaced points. The neutrino-electron scattering and pair creation and annihilation kernelsare tabulated in terms of neutrino energy pairs ε and ε (cid:48) (using the same points as is used for χ and σ IS ), T (using thesame points as in the EoS table), and the degeneracy parameter η = µ e / ( k B T ) ∈ [1 × − , . × ] , using N η = 60 logarithmically spaced points. To evaluate dependent variables from the table, following, e.g., Mezzacappa & Messer(1999), we use bilinear interpolation (or the higher-dimensional equivalent), while derivatives with respect to any of theindependent variables are computed by taking the derivative of the interpolation formula. When interpolating the opacitykernels, we enforce the symmetries in Eq. (10).• Nonlinear solvers – The four nonlinear solvers are implemented following the description in Section 4, with one exceptionthat, in the implementation, the effective emission and absorption opacity χ in Eq. (15) is lagged. In other words, whensolving Eq. (52), χ is evaluated at the starting matter state ( Y ( ∗ ) e , (cid:15) ( ∗ ) ) and is not being updated in the solution procedure.We choose to lag χ in the nonlinear solvers to simplify the Jacobian calculation in Newton’s method. Since χ is usuallyvarying slowly in time, lagging χ has minimal impact on the solution accuracy. For the fixed-point solvers, χ can be8 L AIU ET AL .updated at each iteration at a minor additional cost. As mentioned in Section 4, we choose the truncation parameter inAnderson acceleration to be m = 2 , unless otherwise specified. In this case, the least-squares problem for determining α ∗ in Eq. (59) becomes an inversion of a × matrix and is thus solved analytically. A numerical justification of this choiceof m is given in Section 5.2. In the Newton’s method, the Jacobian matrix is constructed using the derivatives given inAppendix A.We also note that Jacobian-free Newton-Krylov (JFNK) methods, where the Newton step is computed by solving anapproximate Newton system with a Krylov solver, are not well suited for this problem (see, e.g., Knoll & Keyes (2004)for a comprehensive survey on these methods). The reason is that, in JFNK methods, one evaluation of the opacity kernelsis needed in every Krylov iteration to approximate the Jacobian matrix. Thus, JFNK methods require several opacityevaluations per Newton iteration, while the standard Newton’s method only requires one, which makes JFNK methodsmore expensive for solving these problem, where the opacity evaluation is a dominant computational cost.• Nonlinear solver initial guess – When solving the coupled nonlinear system inn Eq. (52), a natural choice of initial guess forthe unknowns ( Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) ) is ( Y ( ∗ ) e , (cid:15) ( ∗ ) , J ( ∗ ) , ¯ J ( ∗ ) ) , the weighted sum from explicit steps defined in Eq. (49).In this work, we add a “presolve” step that aims to provide a better starting point for the iterative solvers and speedupthe computation. Specifically, for given ( Y ( ∗ ) e , (cid:15) ( ∗ ) , J ( ∗ ) , ¯ J ( ∗ ) ) , the presolve step solves a subsystem of Eq. (52), which isobtained by setting the opacities η S C , η T P , χ S C , and χ T P in Eq. (53) to be zero; i.e., with only emission, absorption, and isoen-ergetic scattering. The solution of this simplified system then serves as the initial guess for ( Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) ) in thenonlinear solvers for system in Eq. (52). This presolve step is computationally inexpensive (no opacity table interpolationsare needed, since the emission and absorption opacity χ is lagged as discussed in the previous paragraph), while giving areasonable initial guess for the full system. We observe that, without the presolve step, the Coupled Newton’s method andthe Nested Newton’s method in Sections 4.2 and 4.4 could diverge if the initial guess ( Y ( ∗ ) e , (cid:15) ( ∗ ) , J ( ∗ ) , ¯ J ( ∗ ) ) is too far awayfrom the solution.• Nonlinear solver convergence criteria – We set the convergence criteria for the iterative solvers based on the relative residualof the system in Eq. (52) at the current iterate. Specifically, for the coupled solvers in Sections 4.1 and 4.2, the convergencecriteria are (cid:12)(cid:12) f Y e ( u [ k ] , U [ k ] ) (cid:12)(cid:12) ≤ tol (cid:12)(cid:12) Y [0] e (cid:12)(cid:12) , (cid:12)(cid:12) f (cid:15) ( u [ k ] , U [ k ] ) (cid:12)(cid:12) ≤ tol (cid:12)(cid:12) (cid:15) [0] (cid:12)(cid:12) , (70a) (cid:13)(cid:13) F J ( u [ k ] , U [ k ] ) (cid:13)(cid:13) ≤ tol (cid:13)(cid:13) J [0] (cid:13)(cid:13) , (cid:13)(cid:13) F ¯ J ( u [ k ] , U [ k ] ) (cid:13)(cid:13) ≤ tol (cid:13)(cid:13) ¯ J [0] (cid:13)(cid:13) , (70b)where tol > is a constant relative tolerance, and ( f Y e , f (cid:15) , F J , F ¯ J ) are the residuals as defined in Eqs. (61) and (62).As for the nested solvers in Sections 4.3 and 4.4, the outer layer (Eqs. (65a) and (67a)) uses the convergence criteria inEq. (70a), while the convergence criteria in the inner layer (Eqs. (65b) and (67b)) are given by (cid:13)(cid:13) F J ( u [ k ] , U [ k,(cid:96) ] ) (cid:13)(cid:13) ≤ tol (cid:13)(cid:13) J [ k, (cid:13)(cid:13) , (cid:13)(cid:13) F ¯ J ( u [ k ] , U [ k,(cid:96) ] ) (cid:13)(cid:13) ≤ tol (cid:13)(cid:13) ¯ J [ k, (cid:13)(cid:13) . (71)When all convergence criteria are satisfied, the solvers return the current iterate as the solution. In all numerical experi-ments, we choose the norm to be the discrete L norm on the energy domain, i.e., (cid:107) J (cid:107) := (cid:0) (cid:80) | S ε | q =1 w (2) q J q (cid:1) / , and we setthe relative tolerance to be tol = 10 − . Remark 3
In the coupled fixed-point and the two nested solvers, the basic version of Anderson acceleration outlined in Sec-tion 4.1 is implemented. It is known (see, e.g., Walker & Ni (2011) and references therein) that the iterations in Andersonacceleration may suffer from “stagnation” when the least-squares problem is ill-conditioned. Due to the choice of small trun-cation parameter (e.g., m = 2 ) and the contractive property of the collision operator, we do not encounter the stagnation issuein any of the numerical tests presented in this section. Stagnation may potentially become a practical concern when applyingAnderson acceleration to solve more complicated systems, e.g., fully-coupled neutrino radiation hydrodynamics. A common ap-proach to mitigate the stagnation issue is to control the condition number of the least-squares problems. This can be achieved by(i) solving a proper reformulation of the least-squares problems using QR factorization (Ni & Walker 2010), (ii) modifying thetruncation parameter m adaptively (Yang et al. 2009), and/or (iii) regularizing the least-squares problems (Scieur et al. 2016).Alternatively, one may consider the recently proposed globally convergent variant of Anderson acceleration (Zhang et al. 2020),in which ill-conditioning is handled via regularization and global convergence is guaranteed using safeguarding steps. OLVERS FOR NEUTRINO - MATTER COUPLING
Relaxation problem
The first class of test problems we consider here is the relaxation problem, where the neutrino and antineutrino transportequations (3a) and (3b) are solved with only the collision terms considered, i.e., the space-homogeneous case where the advectionterms T ( f ) and T ( ¯ f ) are zero. In these relaxation problems, the collision operator relaxes the distributions f and ¯ f to theequilibrium Fermi-Dirac distributions f and ¯ f in Eq. (7), respectively, as time evolves. Due to the lack of advection terms, thereis no spatial coupling. Thus, these problems can be solved independently in space, which makes them ideal test cases for thenonlinear collision system solvers considered in this paper. In this setup, the semidiscrete moment and matter equations reducefrom (46) to c ∂ t V k = C ( V k ) , ∀ k such that x k ∈ S x . (72)For the relaxation problems, we discretize Eq. (72) in time with the backward Euler method. At each time step, this timediscretization results in a coupled system that takes the form of Eq. (50) with the effective time step τ = ∆ t . To solve this system,we apply the nonlinear solvers in Section 4 on the subsystem in Eq. (52) (with τ = ∆ t ) to obtain ( Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) ) , whichare then used to compute ( H (+) , ¯ H (+) ) . We start with trivial initial states for H and ¯ H , which remain unchanged with time.We test the nonlinear solvers on the relaxation problem in Eq. (72) with two initial matter states that present problems withdifferent degrees of collisionality. The first state, which represents the high density, strongly collisional region inside a proto-neutron star, is sampled at radius r (1) = 9 . km from the center of a collapsed stellar core, and the second state, whichrepresents the lower density regions around the surface of a proto-neutron star with relatively weaker collisionality, is sampled atradius r (2) = 39 . km. We obtained the matter states at these two locations from a spherically symmetric core-collapse super-nova simulation ms after core bounce (Liebend¨orfer et al. 2005) (results from the VERTEX code using a 15 M (cid:12) progenitor).Specifically, the matter state at r (1) is given by ρ (1) = 1 . × g cm − , T (1) = 1 . × K, and Y (1) e = 0 . ; andthe matter state at r (2) is ρ (2) = 1 . × g cm − , T (2) = 8 . × K, and Y (2) e = 0 . . The inverse mean freepath associated with the neutrino opacities for these matter states are plotted versus neutrino energy in Figure 3. Consideringthe isoenergetic scattering opacity (which is equal for neutrinos and antineutrinos, and increases with neutrino energy as ε ),the mean free path λ IS = σ − IS varies from about km to about − km in the high collisional state (left panel) in the energyrange ε ∈ [1 , MeV. In the low collisional state, the scattering mean free path varies from λ IS ≈ km to λ IS ≈ − km.Correspondingly, the collision time τ IS = λ IS / c varies from about × − ms to × − ms in the high collisional state case,and from about ms to × − ms in the low collisional state case.We run the simulations from initial time t = 0 to final time t f = 0 . ms for the high collisional case, and from t = 0 tofinal time t f = 30 ms for the low collisional case. This is to guarantee that the distributions have relaxed to the equilibriumdistributions by the end of simulations. Each simulation is solved on a single spatial element (the relaxation problem is space-homogeneous) with geometrically progressing energy elements that divide the energy domain [0 , MeV, where the firstelement has ∆ ε = 4 MeV and the last element has ∆ ε ≈ MeV. In the simulations, the initial matter states take electronfraction Y e, = Y ( i ) e and specific internal energy (cid:15) calculated at ( ρ ( i ) , T ( i ) , Y ( i ) e ) from the EoS, with i = 1 , . The initialneutrino and antineutrino moments are given by J ( r ( i ) , ε ) = ¯ J ( r ( i ) , ε ) = 0 . × exp (cid:0) − ( ε − k B T ( i ) ) × (cid:1) , and H ( r ( i ) , ε ) = ¯ H ( r ( i ) , ε ) = 0 , (73)for i = 1 , , which were chosen such that the initial moments are away from the expected final equilibrium distributions, whilemaking sure the initial moments are realizable, i.e., satisfy Eq. (12).Figure 4 illustrates the initial and final equilibrium electron neutrino and electron antineutrino number densities for the tworelaxation problems. Here the initial number densities are given in Eq. (73), and the final number densities were generatedby solving Eq. (72) using the nodal DG and backward Euler discretization, together with the proposed nonlinear solvers. Wevalidated these results by comparing them to number densities from simulations on finer energy-temporal meshes, in which nonoticeable differences were observed. The number densities for the high collision case are shown in Figure 4a, and the ones forthe low collision case are shown in Figure 4b. For the high collision case, the temperature is . × K at t = 0 , andincreases by . to . × K at t f = 0 . ms, while the electron fraction drops by . from 0.2728 to 0.2347. Forthe low collision case, the temperature raises by . from . × K at t = 0 to . × K at t f = 30 ms, andmeanwhile the electron fraction increases by . from 0.1347 to 0.1376. From these plots, we observe that, in the high collisioncase, the neutrino number density J is at least two orders of magnitude higher than the antineutrino number density ¯ J at the finalequilibrium, thus they are further away from the initial densities than the ones in the low collision case. This is one of the reasons0 L AIU ET AL ..
The first class of test problems we consider here is the relaxation problem, where the neutrino and antineutrino transportequations (3a) and (3b) are solved with only the collision terms considered, i.e., the space-homogeneous case where the advectionterms T ( f ) and T ( ¯ f ) are zero. In these relaxation problems, the collision operator relaxes the distributions f and ¯ f to theequilibrium Fermi-Dirac distributions f and ¯ f in Eq. (7), respectively, as time evolves. Due to the lack of advection terms, thereis no spatial coupling. Thus, these problems can be solved independently in space, which makes them ideal test cases for thenonlinear collision system solvers considered in this paper. In this setup, the semidiscrete moment and matter equations reducefrom (46) to c ∂ t V k = C ( V k ) , ∀ k such that x k ∈ S x . (72)For the relaxation problems, we discretize Eq. (72) in time with the backward Euler method. At each time step, this timediscretization results in a coupled system that takes the form of Eq. (50) with the effective time step τ = ∆ t . To solve this system,we apply the nonlinear solvers in Section 4 on the subsystem in Eq. (52) (with τ = ∆ t ) to obtain ( Y (+) e , (cid:15) (+) , J (+) , ¯ J (+) ) , whichare then used to compute ( H (+) , ¯ H (+) ) . We start with trivial initial states for H and ¯ H , which remain unchanged with time.We test the nonlinear solvers on the relaxation problem in Eq. (72) with two initial matter states that present problems withdifferent degrees of collisionality. The first state, which represents the high density, strongly collisional region inside a proto-neutron star, is sampled at radius r (1) = 9 . km from the center of a collapsed stellar core, and the second state, whichrepresents the lower density regions around the surface of a proto-neutron star with relatively weaker collisionality, is sampled atradius r (2) = 39 . km. We obtained the matter states at these two locations from a spherically symmetric core-collapse super-nova simulation ms after core bounce (Liebend¨orfer et al. 2005) (results from the VERTEX code using a 15 M (cid:12) progenitor).Specifically, the matter state at r (1) is given by ρ (1) = 1 . × g cm − , T (1) = 1 . × K, and Y (1) e = 0 . ; andthe matter state at r (2) is ρ (2) = 1 . × g cm − , T (2) = 8 . × K, and Y (2) e = 0 . . The inverse mean freepath associated with the neutrino opacities for these matter states are plotted versus neutrino energy in Figure 3. Consideringthe isoenergetic scattering opacity (which is equal for neutrinos and antineutrinos, and increases with neutrino energy as ε ),the mean free path λ IS = σ − IS varies from about km to about − km in the high collisional state (left panel) in the energyrange ε ∈ [1 , MeV. In the low collisional state, the scattering mean free path varies from λ IS ≈ km to λ IS ≈ − km.Correspondingly, the collision time τ IS = λ IS / c varies from about × − ms to × − ms in the high collisional state case,and from about ms to × − ms in the low collisional state case.We run the simulations from initial time t = 0 to final time t f = 0 . ms for the high collisional case, and from t = 0 tofinal time t f = 30 ms for the low collisional case. This is to guarantee that the distributions have relaxed to the equilibriumdistributions by the end of simulations. Each simulation is solved on a single spatial element (the relaxation problem is space-homogeneous) with geometrically progressing energy elements that divide the energy domain [0 , MeV, where the firstelement has ∆ ε = 4 MeV and the last element has ∆ ε ≈ MeV. In the simulations, the initial matter states take electronfraction Y e, = Y ( i ) e and specific internal energy (cid:15) calculated at ( ρ ( i ) , T ( i ) , Y ( i ) e ) from the EoS, with i = 1 , . The initialneutrino and antineutrino moments are given by J ( r ( i ) , ε ) = ¯ J ( r ( i ) , ε ) = 0 . × exp (cid:0) − ( ε − k B T ( i ) ) × (cid:1) , and H ( r ( i ) , ε ) = ¯ H ( r ( i ) , ε ) = 0 , (73)for i = 1 , , which were chosen such that the initial moments are away from the expected final equilibrium distributions, whilemaking sure the initial moments are realizable, i.e., satisfy Eq. (12).Figure 4 illustrates the initial and final equilibrium electron neutrino and electron antineutrino number densities for the tworelaxation problems. Here the initial number densities are given in Eq. (73), and the final number densities were generatedby solving Eq. (72) using the nodal DG and backward Euler discretization, together with the proposed nonlinear solvers. Wevalidated these results by comparing them to number densities from simulations on finer energy-temporal meshes, in which nonoticeable differences were observed. The number densities for the high collision case are shown in Figure 4a, and the ones forthe low collision case are shown in Figure 4b. For the high collision case, the temperature is . × K at t = 0 , andincreases by . to . × K at t f = 0 . ms, while the electron fraction drops by . from 0.2728 to 0.2347. Forthe low collision case, the temperature raises by . from . × K at t = 0 to . × K at t f = 30 ms, andmeanwhile the electron fraction increases by . from 0.1347 to 0.1376. From these plots, we observe that, in the high collisioncase, the neutrino number density J is at least two orders of magnitude higher than the antineutrino number density ¯ J at the finalequilibrium, thus they are further away from the initial densities than the ones in the low collision case. This is one of the reasons0 L AIU ET AL .. -6 -4 -2 (a) ρ = 1 . × g cm − , T = 1 . × K, Y e = 0 . -8 -6 -4 -2 (b) ρ = 1 . × g cm − , T = 8 . × K, Y e = 0 . Figure 3.
Neutrino opacities for high (left panel) and moderate (right panel) mass density. In each panel, for neutrinos (solid lines) andantineutrinos (dashed lines), we plot, versus neutrino energy ε , the absorption opacity ( χ and ¯ χ ; red lines), the elastic (isoenergetic) scatteringopacity ( σ IS and ¯ σ IS ; blue lines), the neutrino-electron scattering opacity ( χ S C and ¯ χ S C ; green lines), and the opacity due to pair creation andannihilation ( χ T P and ¯ χ T P ; black lines). The neutrino-electron scattering and pair creation and annihilation opacities where computed with theneutrino and antineutrino number densities set to zero; cf. Eqs. (17) and (19). that the high collision rate problems are more challenging than the low collision rate problems in the earlier stage, as discussedin the following paragraphs.The conservation of lepton number and energy (see Eqs. (27) and (28), respectively) in the relaxation tests is shown in Figure 5.Since the relaxation problem is space-homogeneous, the right-hand sides of Eqs. (27) and (28) are both zero. It then follows fromthe nonlinear system formulation in Eq. (52) that the lepton number and energy are conserved if and only if Eqs. (52a) and (52b)are satisfied. Thus, the lepton number and energy are expected to be conserved up to the nonlinear solver tolerance at each pointon the space-time grid, which is confirmed in the results reported in Figure 5. In addition, we also observe from Figure 5 that,in both the high and low collision rate cases, tightening the nonlinear solver tolerance from − to − indeed improves theconservation results.Figure 6 shows iteration counts versus time for each nonlinear solver on the two relaxation problems, using various time stepsizes: ∆ t = 10 − , − , and − ms. These time step sizes are motivated by the fact that in the context of the IMEX scheme inEqs. (47)-(48), the maximum stable time step is ∆ t = C CFL × × − (cid:0) ∆ x/ (cid:1) ms, where C CFL (cid:46) is the CFL number. Fora spatial resolution ∆ x = O ( ) , our chosen time steps bracket what is typically used in core-collapse supernova simulations.The results for the high collision problem are shown in Figure 6a, while the ones for the low collision problem are shown inFigure 6b. In these figures, the top plot shows the “outer” iteration counts for each solver, while the bottom plot shows theaveraged “inner” iteration counts for the two nested solvers. Here the outer iteration counts represent the number of iterationsneeded for solving Eqs. (54), (60), (65a), and (67a) in the Coupled AA (Anderson acceleration), Coupled Newton, Nested AA,and Nested Newton solvers, respectively; the inner iteration counts are the number of iterations needed for solving Eqs. (65b) and OLVERS FOR NEUTRINO - MATTER COUPLING -6 -4 -2 (a) Number density J ( r (1) , ε, t ) at t = 0 and t = t f = 0 . ms. At r (1) ,with initial matter state ρ = 1 . × g cm − , T = 1 . × K, Y e = 0 . . -6 -4 -2 (b) Number density J ( r (2) , ε, t ) at t = 0 and t = t f = 30 ms. At r (2) ,with initial matter state ρ = 1 . × g cm − , T = 8 . × K, Y e = 0 . . Figure 4.
Initial and final number densities versus neutrino energy for relaxation problems with high (left column) and low (right column)collision rates. Each plot shows the initial number density given in Eq. (73) for both neutrinos and antineutrinos (black dashed line), the finalneutrino number density (blue solid line), and the final antineutrino number density (red solid line). -202 10 -8 -8 (a) Changes in lepton number and energy relative to their initial values inthe relaxation simulation at r (1) from t = 0 to t = t f = 0 . ms. Theinitial matter state is ρ = 1 . × g cm − , T = 1 . × K, Y e = 0 . . -101 10 -6 -7 (b) Changes in lepton number and energy relative to their initial values inthe relaxation simulation at r (2) from t = 0 to t = t f = 30 ms. Theinitial matter state is ρ = 1 . × g cm − , T = 8 . × K, Y e = 0 . . Figure 5.
Lepton number and energy conservation (see Eqs. (27) and (28), respectively) in relaxation problems with high (left column) andlow (right column) collision rate. Each plot shows the relative changes in lepton number and energy in the relaxation simulations with thenonlinear solver tolerance tol set to − (black solid line) and − (black dashed line). Tightening the tolerance from − to − leadsto approximately four orders of magnitude smaller relative changes in both lepton number and energy. (67b) in the Nested AA and Nested Newton solvers, respectively. Since the inner equations (65b) and (67b) are solved in everyouter iteration, the reported inner iteration counts in Figures 6a and 6b are averaged over the number of times that Eqs. (65b) and(67b) were solved, i.e., the total number of inner iterations taken in the solver is the product of the outer iteration count and theaveraged inner iteration count.2 L AIU ET AL .Before comparing the solvers, we first observe from the results that the implicit system in Eq. (52) in the high collision rate caseindeed requires more iterations to reach convergence than the same system does in the low collision rate case. Also, the overalliteration counts grows as ∆ t increases. These observations agree with our expectation, since increasing ∆ t effectively increasesthe collision rates, and higher collision rates result in stronger coupling of the number densities, which makes the system inEq. (52) harder to solve.For the nonlinear solver performance, we observe that the Coupled Newton solver requires more iterations than the CoupledAA solver on harder problems, e.g., problems with higher collision rates, larger time step, or at an earlier stage (further awayfrom equilibrium); on easier problems, the Coupled Newton solver converges in fewer iterations than the Coupled AA solver. Asexpected, the nested solvers indeed reduce the number of outer iterations when compared to the coupled solvers, presumably byproviding a more accurate update of the neutrino and antineutrino number densities from the inner iteration. The Nested AA andNested Newton solvers share nearly identical outer iteration counts. This is due to the fact that both nested solvers use Andersonacceleration in the outer layer (Eqs. (65a) and (67a)), which takes number densities ˆ U ( u ) from the inner layers (Eqs. (65b) and(67b)) of the two solvers. Since the problems in the inner layers are equivalent, the solutions ˆ U ( u ) are identical up to the residualtolerance. The Nested AA solver generally requires more (inner) iterations to converge than the Nested Newton solver does,especially on harder problems. We also note that the comparison of iteration counts here does not fully reflect the performanceof the solvers in terms of computational time. In particular, as discussed in Section 4.3, opacity kernel evaluations make outeriterations much more computationally expensive than the inner iterations. In addition, the iterations in Anderson accelerationare computationally cheaper than the ones in Newton’s method, since they do not require constructing and inverting the Jacobianmatrix. We defer the comparison of computational times to Section 5.3, where a more realistic test problem is considered.Next, we explore the effect of the truncation parameter m on the convergence of Anderson acceleration by comparing theiteration count for the Coupled AA solver with different values of m , on both the high and low collision rate relaxation problemsconsidered in the previous test. In this comparison, the time step size fixed to ∆ t = 10 − ms, and m varies from 0 to 4, where m = 0 resembles the simple Picard iteration. From the results reported in Figure 7, we observe that, on these problems, Andersonacceleration does converge faster as the value of m increases, while the marginal benefit becomes insignificant for m ≥ . Thisresult justifies our choice of m = 2 in the numerical tests throughout the paper, since higher values of m lead to larger memoryfootprints and more expensive least-squares solves in Anderson acceleration with minimal improvement in the iteration counts.We next investigate how early termination of the iterative solvers affects solution accuracy. We choose to test the Nested AAsolver with early termination of the outer loop, and compare the resulting solution to a converged reference solution at the finaltime. Here the solution process is terminated either when the convergence criteria (70a) are satisfied or when the outer iterationcount for solving Eq. (65a) reaches a preset maximum (MaxIter). We test the solver on the relaxation problem with low collisionrate and time step ∆ t = 10 − ms, which is much larger than the usual time step for a stable explicit scheme. The choice ofa large time step makes the effect of early termination more pronounced. Figure 8 reports the iteration counts for MaxIter =1 and MaxIter = 2 on the test problem, along with the electron neutrino and antineutrino (energy-integrated) number densities,temperatures, and electron fractions at the final time. These results are compared to a fully converged reference solution (MaxIter= 100), where the nested fixed-point solver converges well before the nominal maximal outer iteration is reached, as shown inFigure 8a. Here the “presolve” step discussed in Section 5.1 is turned off. We note that when setting MaxIter = 1, the NestedAA solver (w/o presolve) resembles an approach used in earlier works such as Just et al. (2015), where the radiation quantitiesare updated using opacities computed from the lagged matter states, i.e., matter states at time t n are used to update the radiationquantities at t n +1 . The results in Figure 8 show that when MaxIter = 1, while the relative differences in the earlier time steps arerather significant (from 30% to 0.5%), the solution still converges to an identical (up to the solver tolerance) equilibrium at thefinal time. However, when MaxIter = 2, it is clear that the solution converges to a different equilibrium. Indeed, it can be seenfrom Figure 8e that the lepton number and energy are not conserved in the solution with MaxIter = 2, while they are conservedup to the solver tolerance in the fully converged solution (MaxIter = 100). When MaxIter = 2, the changes in lepton numberand energy lead to a different equilibrium. As for the solution with MaxIter = 1, despite the fact that the solution process of thenonlinear system Eq. (52) is terminated early, the lepton number and energy are actually conserved in exact arithmetics. This isbecause, when MaxIter = 1, the early terminated nested iterates satisfy Eqs. (52a)–(52b) exactly, which enforces lepton numberand energy convergences and leads to the correct equilibrium. However, the early terminated solutions generally do not satisfyEqs. (52c)–(52d) and result in rather inaccurate solutions in the transient state. In Figure 9, we repeat the test but with the presolvestep turned on. Here we also include the option MaxIter = 0 into the comparison, in which the radiation and matter quantities areonly updated in the presolve step, i.e., with the NES and pair processes ignored. From Figure 9, it can be observed that even withMaxIter = 0, the presolve step gives fairly accurate solution at the final time, as the lepton number and energy are conserved upto the solver tolerance in the presolve step. When the maximum iteration is allowed to be higher, the presolve step improves the OLVERS FOR NEUTRINO - MATTER COUPLING -4 -3 -2 -1 -4 -3 -2 -1 -4 -3 -2 -1 (a) High collision rate, ρ = 1 . × g cm − , T = 1 . × K, Y e = 0 . -4 -3 -2 -1 -4 -3 -2 -1 -4 -3 -2 -1 (b) Low collision rate, ρ = 1 . × g cm − , T = 8 . × K, Y e = 0 . Figure 6.
Iteration counts of the nonlinear solvers on relaxation problems. The iteration counts of the Coupled AA, Coupled Newton, NestedAA, and Nested Newton solvers on relaxation problems in the high (left column) and low (right column) collision rate regimes with time stepsizes varying from − to − ms. AIU ET AL ..
Iteration counts of the nonlinear solvers on relaxation problems. The iteration counts of the Coupled AA, Coupled Newton, NestedAA, and Nested Newton solvers on relaxation problems in the high (left column) and low (right column) collision rate regimes with time stepsizes varying from − to − ms. AIU ET AL .. -4 -3 -2 -1 (a) High collision rate, ρ = 1 . × g cm − , T = 1 . × K, Y e = 0 . -4 -3 -2 -1 (b) Low collision rate, ρ = 1 . × g cm − , T = 8 . × K, Y e = 0 . Figure 7.
Iteration counts of the Coupled AA solver with truncation parameter m = 0 , . . . , on relaxation problems. solution accuracy, at least in the earlier stage. We observe that while the solution is still inaccurate when MaxIter = 2, it is muchcloser to the reference solution when the presolve is turned on.5.3. Deleptonization problem
We further investigate and compare the performance of the nonlinear solvers in a more realistic setting with a proto-neutronstar deleptonization problem, using initial matter profiles from spherically symmetric CCSN simulations. In this test we solvethe full moment equations presented in Section 2.2, using the DG phase-space discretization in Section 3.1 and the IMEX timeintegration scheme in Section 3.2.For this test, we adopt matter profiles from Liebend¨orfer et al. (2005). Specifically, we use profiles for mass density, tempera-ture, and electron fraction obtained with the
VERTEX code using a 15 M (cid:12) progenitor from Woosley & Weaver (1995) (model G15in Liebend¨orfer et al. (2005)). Other thermodynamics quantities (e.g., internal energy and electron, proton, and neutron chemicalpotentials) are obtained from the tabulated SFHo EoS (Steiner et al. 2013). To investigate sensitivity to conditions encounteredover an extended period covering the neutrino heating phase, we run the comparison on profiles taken at t = 50 , , , and ms after core bounce. The initial matter profiles are plotted in Figure 10.Since we do not have the radiation quantities from Liebend¨orfer et al. (2005), to initialize the radiation field, we adopt theanalytical distribution function from the homogeneous sphere test, f H S , (e.g., Smit et al. (1997)), which is a solution to the steadystate transport problem of radiation emanating from a sphere whose constant absorption opacity and emissivity inside a radius R are χ and χ f , respectively: f H S [ f , χ , R ]( r, µ ) = f (cid:0) − e − χ s [ R ]( r,µ ) (cid:1) , (74)where f is an isotropic equilibrium distribution, s [ R ]( r, µ ) = r µ + R g [ R ]( r, µ ) if r < R , µ ∈ [ − , +1] , R g [ R ]( r, µ ) if r ≥ R , µ ∈ [(1 − ( R /r ) ) / , +1] , otherwise , (75)and g [ R ]( r, µ ) = [1 − ( r/R ) (1 − µ )] / .To adopt the homogeneous sphere distribution to the current setting, we first estimate the energy-dependent neutrinosphereradius R ν ( ε ) (cid:90) ∞ R ν ( ε ) χ ( ε, r ) dr = 23 , (76) OLVERS FOR NEUTRINO - MATTER COUPLING -1 (a) Iteration Counts -1 -10 -5 -1 -10 -5 (b) Energy-integrated Number Densities -1 -10 -5 (c) Temperature -1 -10 -5 (d) Electron Fraction -15 -10 -5 -1 -15 -10 -5 (e) Relative Change in Lepton Number and Energy Figure 8.
Iteration count, energy-integrated number densities, temperature, electron fraction, and lepton number and energy conservation resultsfor the Nested AA solver on the relaxation problem with low collision rate ρ = 1 . × g cm − , T = 8 . × K, Y e = 0 . .The solver is applied without the presolve step and run to various maximum outer iterations (MaxIter), with time step ∆ t = 10 − ms. where χ is the absorption opacity. (Neutrinosphere radii for the various initial profiles used here are plotted versus neutrinoenergy in the lower right panel of Figure 10.) The neutrino distribution function is then set to f ( r, ε, µ ) = (cid:40) f H S [ f ( r, ε ) , χ ( r, ε ) , R ν ( ε )]( r, µ ) , r < R ν ( ε ) f H S [ f ( R ν ( ε )) , χ ( R ν ( ε )) , R ν ( ε )]( r, µ ) , r ≥ R ν ( ε ) , (77)where f ( r, ε ) is taken to be the Fermi-Dirac distribution in Eq. (7). Finally, the initial moments are computed as (cid:8) J , H (cid:9) ( r, ε ) = 12 (cid:90) − f ( r, ε, µ ) µ { , } dµ. (78)(The same procedure is adopted to initialize the antineutrinos.)The evolution of the electron fraction in the deleptonization problem is illustrated in Figure 11, where we plot the electronfraction versus mass density over 10 ms of evolution, starting from the ms profile in Figure 10 (dashed lines). Here theevolution of the electron fraction was generated from the DG-IMEX scheme with the proposed nonlinear solvers as described inSection 5.1. We validated the result by comparing it against solutions from simulations on finer spatial-energy-temporal meshes,in which no noticeable differences were observed. For higher densities ( ρ (cid:38) × g cm − ), neutrinos are effectively trapped,and the electron fraction remains largely unchanged. For lower densities, neutrinos (created by electron capture on protons) arenot trapped and escape the computational domain, which results in a lowering of the electron fraction (deleptonization). Figure 126 L AIU ET AL ..
Iteration count, energy-integrated number densities, temperature, electron fraction, and lepton number and energy conservation resultsfor the Nested AA solver on the relaxation problem with low collision rate ρ = 1 . × g cm − , T = 8 . × K, Y e = 0 . .The solver is applied without the presolve step and run to various maximum outer iterations (MaxIter), with time step ∆ t = 10 − ms. where χ is the absorption opacity. (Neutrinosphere radii for the various initial profiles used here are plotted versus neutrinoenergy in the lower right panel of Figure 10.) The neutrino distribution function is then set to f ( r, ε, µ ) = (cid:40) f H S [ f ( r, ε ) , χ ( r, ε ) , R ν ( ε )]( r, µ ) , r < R ν ( ε ) f H S [ f ( R ν ( ε )) , χ ( R ν ( ε )) , R ν ( ε )]( r, µ ) , r ≥ R ν ( ε ) , (77)where f ( r, ε ) is taken to be the Fermi-Dirac distribution in Eq. (7). Finally, the initial moments are computed as (cid:8) J , H (cid:9) ( r, ε ) = 12 (cid:90) − f ( r, ε, µ ) µ { , } dµ. (78)(The same procedure is adopted to initialize the antineutrinos.)The evolution of the electron fraction in the deleptonization problem is illustrated in Figure 11, where we plot the electronfraction versus mass density over 10 ms of evolution, starting from the ms profile in Figure 10 (dashed lines). Here theevolution of the electron fraction was generated from the DG-IMEX scheme with the proposed nonlinear solvers as described inSection 5.1. We validated the result by comparing it against solutions from simulations on finer spatial-energy-temporal meshes,in which no noticeable differences were observed. For higher densities ( ρ (cid:38) × g cm − ), neutrinos are effectively trapped,and the electron fraction remains largely unchanged. For lower densities, neutrinos (created by electron capture on protons) arenot trapped and escape the computational domain, which results in a lowering of the electron fraction (deleptonization). Figure 126 L AIU ET AL .. -1 (a) Iteration Counts -1 -10 -5 -1 -10 -5 (b) Energy-integrated Number Densities -1 -10 -5 (c) Temperature -1 -10 -5 (d) Electron Fraction -15 -10 -5 -1 -15 -10 -5 (e) Relative Change in Lepton Number and Energy Figure 9.
Iteration count, energy-integrated number densities, temperature, electron fraction, and lepton number and energy conservation resultsfor the Nested AA solver on the relaxation problem with low collision rate ρ = 1 . × g cm − , T = 8 . × K, Y e = 0 . .The solver is applied with the presolve step and run to various maximum outer iterations (MaxIter), with time step ∆ t = 10 − ms. shows conservations of lepton number and energy over ms of evolution in the deleptonization problem starting from the mspost-bounce matter profile. Here the lepton number and energy both comprise two parts: (i) the interior lepton number and energyin the computation domain, i.e., the time integrals of left-hand sides in Eqs. (27) and (28), respectively, and (ii) the accumulatedoutflow lepton number and energy at the boundary, which are respectively the negations of right-hand sides in Eqs. (27) and(28), integrated in time. The conservation results of the lepton number/energy, as well as the evolutions of the interior leptonnumber/energy and the accumulated outflow lepton number/energy, are illustrated in Figure 12a. The evolutions of individualcomponents of the interior lepton number/energy, including matter lepton number (electron number), neutrino lepton number,internal energy, and neutrino energy, are shown in Figure 12b.For each profile, the deleptonization problems are simulated from the profile time ( , , , and ms after core bounce)to 5 ms after the profile time, which are referred to as the initial time t = 0 and the final time t f = 5 ms in the remainder of thepaper. The IMEX time integration scheme discussed in Section 5.1 is used in the simulations, where the time step ∆ t determinedby the stability requirement for the explicit advection part. These problems are solved on the spatial domain r ∈ [0 , km andenergy domain [0 , MeV, which are divided into and geometrically progressing elements, respectively. Here the firstspatial element has ∆ x = 1 km, the last spatial element has ∆ x ≈ . km, the first energy element has ∆ ε = 4 MeV, and thelast energy element has ∆ ε ≈ MeV.Figure 13 shows the iteration counts of the nonlinear solvers on the deleptonization problem for various profiles. The resultsfor profiles taken at , , , and ms after core bounce are illustrated in Figures 13a, 13b, 13c, and 13d, respectively.As in Figure 6, the top plot in these figures shows the “outer” iteration counts of each solver, while the bottom plot shows the OLVERS FOR NEUTRINO - MATTER COUPLING (a) Mass Density (b) Temperature(c) Electron Fraction (d) Neutrinosphere Radius Figure 10.
Initial matter profiles used in the deleptonization problem, taken from Liebend¨orfer et al. (2005). Plotted versus radius are massdensity (upper left panel), temperature (upper right panel), and electron fraction (lower left panel) for various post-bounce times: ms (solid), ms (dashed), ms (dotted), and ms (dash-dot). The neutrinosphere radii, defined in Eq. (76), for the respective profiles — forelectron neutrinos (black) and electron antineutrinos (grey) — are also plotted versus neutrino energy (lower right panel). averaged “inner” iteration counts of the nested solvers, where the inner iteration counts are averaged over the number of timesthat the inner equations (65b) or (67b) were solved. In addition, here both the outer and inner iteration counts are averaged overtime (from t to t f ), and the averaged iteration counts are plotted against the mass density, which corresponds to spatial locationsfor each profile. (We have found that the number of iterations at a given location varies little from t to t f .)From Figure 13, we observe that the simulations with all four profiles give consistent results — the nested solvers requirefewer outer iterations than the coupled ones do, and the Nested AA solver requires more inner iteration to converge than theNested Newton solver does, especially for harder problems (at higher mass density). This observation also agrees with the resultsreported in Section 5.2 on the relaxation problems. Computational times for these simulations are reported in the top panel ofTable 1, where the tests t Tot as well as the detailed timingmeasurements in each simulation, such as the computational time spent on (i) solving the nonlinear system in Eq. (52) in theimplicit step ( t Im ), (ii) the opacity evaluations/interpolations when solving Eq. (52) ( t Op ), (iii) linear algebra operations, such as8 L AIU ET AL ..
Initial matter profiles used in the deleptonization problem, taken from Liebend¨orfer et al. (2005). Plotted versus radius are massdensity (upper left panel), temperature (upper right panel), and electron fraction (lower left panel) for various post-bounce times: ms (solid), ms (dashed), ms (dotted), and ms (dash-dot). The neutrinosphere radii, defined in Eq. (76), for the respective profiles — forelectron neutrinos (black) and electron antineutrinos (grey) — are also plotted versus neutrino energy (lower right panel). averaged “inner” iteration counts of the nested solvers, where the inner iteration counts are averaged over the number of timesthat the inner equations (65b) or (67b) were solved. In addition, here both the outer and inner iteration counts are averaged overtime (from t to t f ), and the averaged iteration counts are plotted against the mass density, which corresponds to spatial locationsfor each profile. (We have found that the number of iterations at a given location varies little from t to t f .)From Figure 13, we observe that the simulations with all four profiles give consistent results — the nested solvers requirefewer outer iterations than the coupled ones do, and the Nested AA solver requires more inner iteration to converge than theNested Newton solver does, especially for harder problems (at higher mass density). This observation also agrees with the resultsreported in Section 5.2 on the relaxation problems. Computational times for these simulations are reported in the top panel ofTable 1, where the tests t Tot as well as the detailed timingmeasurements in each simulation, such as the computational time spent on (i) solving the nonlinear system in Eq. (52) in theimplicit step ( t Im ), (ii) the opacity evaluations/interpolations when solving Eq. (52) ( t Op ), (iii) linear algebra operations, such as8 L AIU ET AL .. Figure 11.
Electron fraction versus mass density over ms of evolution from the initial state given by the ms post-bounce matter profilein Figure 10. The time evolution is indicated by the grayscale, which goes from lightest (initial state) to darkest (final state). the least-squares solve in Anderson acceleration or the assembly and inversion of Jacobian matrices in Newton’s method ( t LA ),(iv) the presolve step ( t Ps ), (v) the explicit update of the advection term ( t Ex ), and (vi) the positivity limiter ( t PL ) . The reportedtiming results in the top panel of Table 1 are all linearly scaled such that the highest reported computational time (16,085 secondsfrom test , which corresponds to the total computation time when the Coupled Newton solver is used tosolve the deleptonization problem with the 100 ms profile. For each simulation, we also record the solver configurations, suchas the value of the truncation parameter m in Anderson acceleration, the maximum allowed iteration (MaxIter), and whether thepresolve step is performed or not. We also note that these reported data are not mutually exclusive, e.g., t Tot ≈ t Im + t Ex + t PL ,and t Im ≈ t Op + t LA + t Ps . Figure 14 provides a column chart that visualizes the results with the 100 ms profile reported in the toppanel of Table 1. There is no qualitative difference between the results from problems with different profiles. The column chart inFigure 14 confirms that the majority of the computational time in these simulations is spent on opacity evaluation/interpolation,which is proportional to the outer iteration counts. It also shows that the nested solvers indeed speed up the computations bytaking inner iterations to reduce the number of outer iterations. Further, we observe that the lower computational cost on the linearalgebra operations ( t LA ) makes the Nested AA solver outperform the Nested Newton solver in terms of the total computationtime by about 8%, despite the higher inner iteration counts. These results indicate that the Nested AA solver leads to the leastcomputation time for the deleptonization problems among all the tested solvers.To further analyze the performance of the Nested AA solver, we experiment with different solver parameters on the deleptoniza-tion problem using the ms profile. Specifically, we tested the Nested AA solver with the Anderson acceleration truncationparameter set to m = 0 , , and for solving the outer system in Eq. (65a). This experiment is to verify the benefit of Andersonacceleration on solving the smaller outer system in Eq. (65a), which has only two unknowns. In addition, for each choice of m ,we initialize the solver with and without the “presolve” step introduced in Section 5.1, which helps us assess the effect of the“presolve” step in a more realistic setting. The resulting iteration counts are shown in Figure 15, which are averaged over the Here the positivity limiters are applied after each explicit and implicit update to enforce realizability of the moments. The computation time for the positivitylimiter ( t PL ) reported in Table 1 is relatively large (compared to t Ex ), and we expect that t PL could be further reduced by a more sophisticated implementation.In any case, it does not affect the observations we made in this paper. OLVERS FOR NEUTRINO - MATTER COUPLING (a) Lepton number and energy conservation (b) Evoluations of matter lepton number, neutrino lepton number, internalenergy, and neutrino energy. Figure 12.
Conservations in the deleptonization problem – The left panel shows the lepton number and energy (see Eqs. (27) and (28)) over ms of evolution starting from the ms post-bounce matter profile. Here the interior lepton number/energy in the computation domain areplotted in black dashed lines, and the accumulated outflow lepton number/energy are plotted in black dotted lines. Note that the accumulatedoutflows (black-dotted lines) are shifted up by the value of initial lepton number/energy for better illustration. The conserved lepton numberand energy are plotted in black solid line. The right panel shows the evolution of components of the interior lepton number and interior energy,including matter lepton number (electron number), neutrino lepton number, internal energy, and neutrino energy. time from t to t f as in Figure 13. The computation times are reported in the bottom panel of Table 1 (tests ∼ ( m = 0) to Anderson acceleration ( m > does reduce the outeriteration count, and thus improves the computation time by around 18%, especially for harder problems (at higher mass density).However, unlike the results for Coupled AA solver reported in Figure 7, increasing the truncation parameter from m = 1 to m = 2 does not lead to any observable reduction in the either iteration counts or computation time. This result is not unexpected,since Anderson acceleration is applied here to solve the outer system in Eq. (65a) with only two unknowns, while the resultsreported in Figure 7 are from solving the fully coupled system in Eq. (54). Another observation from Figure 15 is that there is noclear benefit in applying the presolve step on this problem. The presolve step does slightly reduce the iteration counts, however,the additional cost of the presolve step wipes out the gains from fewer iterations. We suspect that the diminished benefit of thepresolve step is due to the fact that in the deleptonization problem, the explicit time step restriction from the advection term limitsthe stiffness of the implicit problem and forces the implicit update to be small, which makes the effect of presolve insignificant.For deleptonization problems, larger implicit time steps could potentially be achieved by techniques such as subcycling the ex-plicit steps, or by adopting the general multirate framework proposed by Sandu (2019). We expect to see a greater impact of thepresolve step on these problems with larger time steps, as is observed in the relaxation tests.Finally, we investigate the effect of early termination for the Nested AA solver on the deleptonization problems. In these tests,the solver configurations are identical to the ones reported in Figures 8 and 9 for the relaxation problems, e.g., the outer loopof the Nested AA solver is terminated early by restricting the maximum number of outer iterations (MaxIter). Here we test thesolver on the deleptonization problem with the 100 ms initial matter profile for MaxIter = 1, 2, and 100, and the presolve stepdiscussed in Section 5.1 turned off. We then repeat the test for MaxIter = 0, 1, 2, and 100, with the presolve step turned on. Theresults with and without the presolve step are shown in Figures 16 and 17, respectively, where the iteration counts are reported,along with the electron neutrino and antineutrino (energy-integrated) number densities, temperatures, and electron fractions at thefinal time t = 5 ms. Here the fully converged solutions (MaxIter = 100) are considered as reference solutions, where the nestedsolver converges well before the nominal maximal outer iteration is reached, as shown in Figures 16a and 17a. As mentioned inthe earlier subsection, the Nested AA solver with MaxIter = 1 and the presolve step turned off is effectively lagging the opacitiesby computing them from the matter states at the previous time step, while updating the radiation quantities at the current timestep (see, e.g., Just et al. (2015)). From Figure 16, we observe that the early terminated solutions are in good agreement withthe reference solution. However, allowing two outer iterations does not give a better solution than the one with only a single0 L AIU ET AL ..
Conservations in the deleptonization problem – The left panel shows the lepton number and energy (see Eqs. (27) and (28)) over ms of evolution starting from the ms post-bounce matter profile. Here the interior lepton number/energy in the computation domain areplotted in black dashed lines, and the accumulated outflow lepton number/energy are plotted in black dotted lines. Note that the accumulatedoutflows (black-dotted lines) are shifted up by the value of initial lepton number/energy for better illustration. The conserved lepton numberand energy are plotted in black solid line. The right panel shows the evolution of components of the interior lepton number and interior energy,including matter lepton number (electron number), neutrino lepton number, internal energy, and neutrino energy. time from t to t f as in Figure 13. The computation times are reported in the bottom panel of Table 1 (tests ∼ ( m = 0) to Anderson acceleration ( m > does reduce the outeriteration count, and thus improves the computation time by around 18%, especially for harder problems (at higher mass density).However, unlike the results for Coupled AA solver reported in Figure 7, increasing the truncation parameter from m = 1 to m = 2 does not lead to any observable reduction in the either iteration counts or computation time. This result is not unexpected,since Anderson acceleration is applied here to solve the outer system in Eq. (65a) with only two unknowns, while the resultsreported in Figure 7 are from solving the fully coupled system in Eq. (54). Another observation from Figure 15 is that there is noclear benefit in applying the presolve step on this problem. The presolve step does slightly reduce the iteration counts, however,the additional cost of the presolve step wipes out the gains from fewer iterations. We suspect that the diminished benefit of thepresolve step is due to the fact that in the deleptonization problem, the explicit time step restriction from the advection term limitsthe stiffness of the implicit problem and forces the implicit update to be small, which makes the effect of presolve insignificant.For deleptonization problems, larger implicit time steps could potentially be achieved by techniques such as subcycling the ex-plicit steps, or by adopting the general multirate framework proposed by Sandu (2019). We expect to see a greater impact of thepresolve step on these problems with larger time steps, as is observed in the relaxation tests.Finally, we investigate the effect of early termination for the Nested AA solver on the deleptonization problems. In these tests,the solver configurations are identical to the ones reported in Figures 8 and 9 for the relaxation problems, e.g., the outer loopof the Nested AA solver is terminated early by restricting the maximum number of outer iterations (MaxIter). Here we test thesolver on the deleptonization problem with the 100 ms initial matter profile for MaxIter = 1, 2, and 100, and the presolve stepdiscussed in Section 5.1 turned off. We then repeat the test for MaxIter = 0, 1, 2, and 100, with the presolve step turned on. Theresults with and without the presolve step are shown in Figures 16 and 17, respectively, where the iteration counts are reported,along with the electron neutrino and antineutrino (energy-integrated) number densities, temperatures, and electron fractions at thefinal time t = 5 ms. Here the fully converged solutions (MaxIter = 100) are considered as reference solutions, where the nestedsolver converges well before the nominal maximal outer iteration is reached, as shown in Figures 16a and 17a. As mentioned inthe earlier subsection, the Nested AA solver with MaxIter = 1 and the presolve step turned off is effectively lagging the opacitiesby computing them from the matter states at the previous time step, while updating the radiation quantities at the current timestep (see, e.g., Just et al. (2015)). From Figure 16, we observe that the early terminated solutions are in good agreement withthe reference solution. However, allowing two outer iterations does not give a better solution than the one with only a single0 L AIU ET AL .. (a) Profile: 50 ms (b) Profile: 100 ms (c) Profile: 150 ms (d) Profile: 250 ms Figure 13.
Time-averaged iteration counts of the nonlinear solvers on deleptonization problems – The time-averaged iteration counts of theCoupled AA, Coupled Newton, Nested AA, and Nested Newton solvers on deleptonization problems with various profiles. outer iteration, which is possibly due to the issue on lepton number and energy conservation for early terminated solutions, asdiscussed in Section 5.2. Another potential reason is that Anderson acceleration does not guarantee monotone decreasing of theresidual (see Pollock & Rebholz (2019); Evans et al. (2020); Kelley (2018); Toth & Kelley (2015) for convergence analysis ofAA). The behavior of residuals from Anderson acceleration alternating increasing and decreasing has been observed in Pollock& Rebholz (2019), and a potential explanation is given from (Pollock & Rebholz 2019, Theorem 4.5). Similar results can beobserved in Figure 17, where the presolve step is in effect. Here the comparison includes the case that MaxIter = 0, where boththe radiation and matter quantities are solely updated in the presolve step, with the NES and pair processes omitted. The relativedifference in the solution with MaxIter = 0 is mostly around − to − , however, the difference could go up to for theenergy-integrated antineutrino number density in the high mass density region. We also note that the presolve step reduces thedifference in the solution with MaxIter = 2 by a few orders of magnitude in the low mass density region. The reason is that, in thelow density region, the nonlinear solve usually converges within two iterations with the presolve step. The results in Figures 16and 17 suggest that, for problems requiring few iterations, limiting MaxIter to one can give sufficiently accurate solutions, whilereducing the computation time by roughly a factor of two from the fully converged case, as shown in the results for tests ∼ OLVERS FOR NEUTRINO - MATTER COUPLING Figure 14.
Computation time breakdown for the nonlinear solvers on the deleptonization problem with 100 ms profile – The opacity inter-polation time ( t Op in Table 1) and the dense linear algebra operation time ( t LA in Table 1) are the blue and yellow sections of the columns,respectively. The green sections represent the time spent on all other computations. Figure 15.
Time-averaged iteration counts for the Nested AA algorithm with various solver configurations on the deleptonization problem with ms profile. Results for the Nested AA algorithm without/with presolve and with truncation parameter value m = 0 , , in the outer loopare reported. Here the iteration counts for m = 2 overlap with the corresponding reuslts for m = 1 , which implies that there is no benefitmoving from m = 1 to m = 2 when solving the outer loop problem. AIU ET AL ..
Time-averaged iteration counts for the Nested AA algorithm with various solver configurations on the deleptonization problem with ms profile. Results for the Nested AA algorithm without/with presolve and with truncation parameter value m = 0 , , in the outer loopare reported. Here the iteration counts for m = 2 overlap with the corresponding reuslts for m = 1 , which implies that there is no benefitmoving from m = 1 to m = 2 when solving the outer loop problem. AIU ET AL .. (a) Time-averaged Iteration Counts -10 -5 -10 -5 (b) Energy-integrated Number Densities at t f = 5 ms -10 -5 (c) Temperature at t f = 5 ms -10 -5 (d) Electron Fraction at t f = 5 ms Figure 16.
Time-averaged iteration counts and final-time ( t f = 5 ms) solutions for the Nested AA solver without the presolve step, running tovarious maximum outer iterations (MaxIter). The blue line (MaxIter = ) is considered as a reference. OLVERS FOR NEUTRINO - MATTER COUPLING (a) Time-averaged Iteration Counts -10 -5 -10 -5 (b) Energy-integrated Number Densities at t f = 5 ms -10 -5 (c) Temperature at t f = 5 ms -10 -5 (d) Electron Fraction at t f = 5 ms Figure 17.
Time-averaged iteration counts and final-time ( t f = 5 ms) solutions for the Nested AA solver with the presolve step, running tovarious maximum outer iterations (MaxIter). The blue line (MaxIter = ) is considered as a reference. AIU ET AL ..
Time-averaged iteration counts and final-time ( t f = 5 ms) solutions for the Nested AA solver with the presolve step, running tovarious maximum outer iterations (MaxIter). The blue line (MaxIter = ) is considered as a reference. AIU ET AL .. Table 1.
Overview of nonlinear solver timing results for the deleptonization problem.Top panel: The detailed computational time is reported for the results shown in Figure 13. Here the four iterative solvers are compared on fourinitial matter profiles with identical solver parameters. The results are linearly scaled so that the highest measurement (boxed, 16,085 seconds)is scaled to 100. t Tot t Im t Op t LA t Ps t Ex t PL m Presolve MaxIter1 Coupled-Newton 50 ms 98.8 96.8 78.6 11.7 2.7 0.8 1.0 – yes 1002 Coupled-Newton 100 ms t Tot t Im t Op t LA t Ps t Ex t PL m Presolve MaxIter14 Nested-AA 100 ms t Tot : total simulation time; t Im : implicit solution time; t Op : opacity interpolation time; t LA : dense linear algebra time; t Ps : initial presolve time; t Ex : explicit update time; t PL : positivity limiter time; m : Anderson acceleration truncation parameter;Presolve: whether the presolve step is performed in solver initialization;MaxIter: maximum allowed (outer) iteration for (nested) iterative solvers;Relations: t Tot ≈ t Im + t Ex + t PL , t Im ≈ t Op + t LA + t Ps . OLVERS FOR NEUTRINO - MATTER COUPLING SUMMARY AND DISCUSSIONWe have investigated several iterative solvers for nonlinear systems arising from the discretization of a non-relativistic two-moment model for neutrino transport with opacities from Bruenn (1985), coupled with static matter configurations. Specifically,we have incorporated the nonlinear solvers in a DG-IMEX scheme, as implemented in the toolkit for high-order neutrino-radiationhydrodynamics ( thornado ). Within the IMEX time integration scheme currently adopted in thornado , updating the neutrinotransport and matter equations requires solving a coupled nonlinear system on the radiation moments and matter states (internalenergy and electron fraction). We have considered two approaches to solve the nonlinear system — a coupled approach thatdirectly solves the fully coupled system, and a nested approach that formulates the nonlinear system as a nested system with theouter system governing the matter states and the inner system governing the neutrino number densities. The nested approach isintroduced to reduce the number of opacity evaluations/interpolations required in the solution procedure, and thus becomes moreefficient than the coupled approach. Two iterative solvers — the Anderson accelerated fixed point solver and Newton’s method— are implemented for both the coupled and nested approaches. We have tested the four solvers on relaxation problems withvarious collision rates and time steps, as well as on proto-neutron star deleptonization problems with post-bounce matter profilesfrom spherically symmetric CCSN simulations.Numerical results confirm that both nested solvers indeed require fewer iterations to converge (and thus less computationaltime) than the coupled solvers, due to the fewer number of opacity interpolations performed in the solution procedure. The nestedAnderson acceleration solver requires more inner iterations to converge, but, due to the low cost per iteration, less computationtime than the nested Newton’s method, which is a consequence of the heavier dense linear algebra operations in Newton’s method.In addition to the advantage in computation time, another benefit for using solvers based on Anderson acceleration over Newton’smethod is the simplicity of implementation, particularly for solving problems in which the derivatives are not readily available,such as the coupled nonlinear systems in CCSN simulations considered in this paper. For the test problems considered in thispaper, we also observe that forcing the nested Anderson acceleration solver to terminate after the first outer iteration could leadto reasonably accurate results. This observation confirms that, on these problems, solving the nonlinear coupled system in theimplicit step using lagged opacity kernels from the previous time step could give a sufficiently accurate solution, which has alsobe observed by others (e.g., Just et al. (2015)).Moving forward, we will continue to expand on the capabilities in thornado , and incorporate the more comprehensivephysics needed for realistic CCSN models. Next steps toward this goal include (i) incorporate special and general relativisticeffects into the neutrino transport model, (ii) include muon and tau neutrinos, (iii) update the opacity set to include, e.g., modernelectron capture rates, bremsstrahlung, and inelastic scattering on nucleons, (iv) couple the neutrino transport equations with fluidequations to self-consistently model neutrino-radiation hydrodynamics, (v) port the nonlinear solvers to modern hardware archi-tectures (e.g., GPUs), and further analyze implementation performance, and (vi) compare thornado to other well-establishedCCSN simulation codes, such as
AGILE - BOLTZTRAN (Liebend¨orfer et al. 2004), following the approaches in Just et al. (2015)and O’Connor et al. (2018). We are making progress in these directions, and plan to report on the results in future publications.ACKNOWLEDGMENTSThis manuscript has been authored, in part, by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US De-partment of Energy (DOE). The US government retains and the publisher, by accepting the article for publication, ac-knowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or repro-duce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will pro-vide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan( http://energy.gov/downloads/doe-public-access-plan ).This research was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Departmentof Energy Office of Science and the National Nuclear Security Administration. This research used resources of the Oak RidgeLeadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S.Department of Energy under Contract No. DE-AC05-00OR22725.Support for DOI 10.13139/OLCF/1735948 dataset is provided by the U.S. Department of Energy, project AST137 underContract DE-AC05-00OR22725. Project AST137 used resources of the Oak Ridge Leadership Computing Facility at Oak RidgeNational Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR227256 L
AIU ET AL .APPENDIX A. KERNEL DERIVATIVESThe NES and pair kernels in Eqs. (16), (17), (18), and (19) are tabulated in terms of temperature T and degeneracy parameter η = µ e /kT ; i.e., for given ε and ε (cid:48) they are functions of the form Φ = Φ(
T, η ) . (A1)However, for the coupled Newton’s method in Section 4.2, we need kernel derivatives with respect to internal energy (cid:15) andelectron fraction Y e . The electron chemical potential, along with all the other quantities given by the EoS, is tabulated in termsof ρ , T , and Y e . On the one hand, the variation of the kernel in Eq. (A1) is d Φ = (cid:16) ∂ Φ ∂T (cid:17) η dT + (cid:16) ∂ Φ ∂η (cid:17) T dη = 1 kT (cid:16) ∂ Φ ∂η (cid:17) T (cid:16) ∂µ e ∂ρ (cid:17) T,Y e dρ + (cid:110) (cid:16) ∂ Φ ∂T (cid:17) η + (cid:16) ∂ Φ ∂η (cid:17) T (cid:104) kT (cid:16) ∂µ e ∂T (cid:17) ρ,Y e − ηT (cid:105) (cid:111) dT + 1 kT (cid:16) ∂ Φ ∂η (cid:17) T (cid:16) ∂µ e ∂Y e (cid:17) ρ,T dY e , (A2)where we used dη = (cid:16) ∂η∂ρ (cid:17) T,Y e dρ + (cid:16) ∂η∂T (cid:17) ρ,Y e dT + (cid:16) ∂η∂Y e (cid:17) ρ,T dY e = 1 kT (cid:16) ∂µ e ∂ρ (cid:17) T,Y e dρ + (cid:104) kT (cid:16) ∂µ e ∂T (cid:17) ρ,Y e − ηT (cid:105) dT + 1 kT (cid:16) ∂µ e ∂Y e (cid:17) ρ,T dY e . (A3)On the other hand, by considering the kernel as a function of (cid:15) and Y e , the variation is d Φ = (cid:16) ∂ Φ ∂(cid:15) (cid:17) Y e d(cid:15) + (cid:16) ∂ Φ ∂Y e (cid:17) (cid:15) dY e = (cid:16) ∂ Φ ∂(cid:15) (cid:17) Y e (cid:16) ∂(cid:15)∂ρ (cid:17) T,Y e dρ + (cid:16) ∂ Φ ∂(cid:15) (cid:17) Y e (cid:16) ∂(cid:15)∂T (cid:17) ρ,Y e dT + (cid:104) (cid:16) ∂ Φ ∂(cid:15) (cid:17) Y e (cid:16) ∂(cid:15)∂Y e (cid:17) ρ,T + (cid:16) ∂ Φ ∂Y e (cid:17) (cid:15) (cid:105) dY e , (A4)where we used d(cid:15) = (cid:16) ∂(cid:15)∂ρ (cid:17) T,Y e dρ + (cid:16) ∂(cid:15)∂T (cid:17) ρ,Y e dT + (cid:16) ∂(cid:15)∂Y e (cid:17) ρ,T dY e . (A5)Comparing Eqs. (A2) and (A4), we have (cid:16) ∂ Φ ∂T (cid:17) η + (cid:104) kT (cid:16) ∂µ e ∂T (cid:17) ρ,Y e − ηT (cid:105) (cid:16) ∂ Φ ∂η (cid:17) T = (cid:16) ∂(cid:15)∂T (cid:17) ρ,Y e (cid:16) ∂ Φ ∂(cid:15) (cid:17) Y e , (A6) kT (cid:16) ∂µ e ∂Y e (cid:17) ρ,T (cid:16) ∂ Φ ∂η (cid:17) T = (cid:16) ∂(cid:15)∂Y e (cid:17) ρ,T (cid:16) ∂ Φ ∂(cid:15) (cid:17) Y e + (cid:16) ∂ Φ ∂Y e (cid:17) (cid:15) . (A7)Solving for ( ∂ Φ /∂(cid:15) ) Y e and ( ∂ Φ /∂Y e ) (cid:15) we obtain the derivatives we need for Newton’s method in terms of derivatives that canbe computed directly from the opacity and EoS tables (cid:16) ∂ Φ ∂(cid:15) (cid:17) Y e = (cid:104) kT (cid:16) ∂µ e ∂T (cid:17) ρ,Y e − ηT (cid:105) / (cid:16) ∂(cid:15)∂T (cid:17) ρ,Y e (cid:16) ∂ Φ ∂η (cid:17) T + (cid:16) ∂(cid:15)∂T (cid:17) − ρ,Y e (cid:16) ∂ Φ ∂T (cid:17) η , (A8) (cid:16) ∂ Φ ∂Y e (cid:17) (cid:15) = (cid:110) kT (cid:16) ∂µ e ∂Y e (cid:17) ρ,T − (cid:16) ∂(cid:15)∂Y e (cid:17) ρ,T (cid:104) kT (cid:16) ∂µ e ∂T (cid:17) ρ,Y e − ηT (cid:105) / (cid:16) ∂(cid:15)∂T (cid:17) ρ,Y e (cid:111) (cid:16) ∂ Φ ∂η (cid:17) T − (cid:110) (cid:16) ∂(cid:15)∂Y e (cid:17) ρ,T / (cid:16) ∂(cid:15)∂T (cid:17) ρ,Y e (cid:111) (cid:16) ∂ Φ ∂T (cid:17) η . (A9)As discussed in Section 5.1 (and following, e.g., Mezzacappa & Messer (1999)), tabulated quantities are evaluated using bilinearor trilinear interpolation, while derivatives are estimated by direct differentiation of the respective interpolation formulae. OLVERS FOR NEUTRINO - MATTER COUPLING
Adams, M. L. 2001, Nuclear science and engineering, 137, 298An, H., Jia, X., & Walker, H. F. 2017, Journal of ComputationalPhysics, 347, 1 , doi: https://doi.org/10.1016/j.jcp.2017.06.031Anderson, D. G. 1965, J. ACM, 12, 547,doi: 10.1145/321296.321305Asaithambi, R., & Mahesh, K. 2017, Journal of ComputationalPhysics, 341, 377 ,doi: https://doi.org/10.1016/j.jcp.2017.04.025Ascher, U., Ruuth, S., & Spiteri, R. 1997, Applied NumericalMathematics, 25, 151Bassi, F., Franchina, N., Ghidoni, A., & Rebay, S. 2013,International Journal for Numerical Methods in Fluids, 71, 1322Bethe, H. A., & Wilson, J. R. 1985, Astrophysical Journal, 295, 14Bruenn, S., Mezzacappa, A., Hix, W., et al. 2009, AIP Conf. Proc.,1111, 593, doi: 10.1063/1.3141615Bruenn, S. W. 1985, Astrophysical Journal Supplement Series, 58,771Bruenn, S. W., De Nisco, K. R., & Mezzacappa, A. 2001,Astrophysical Journal, 560, 326Bruenn, S. W., Blondin, J. M., Hix, W. R., et al. 2020, ApJS, 248,11Burrows, A. 2013, Reviews of Modern Physics, 85, 245Burrows, A., Radice, D., Vartanyan, D., et al. 2020, MNRAS, 491,2715Burrows, A., Reddy, S., & Thompson, T. A. 2006, Nuclear PhysicsA, 777, 356Burrows, A., Vartanyan, D., Dolence, J. C., Skinner, M. A., &Radice, D. 2018, Space Science Review, 214, 33Cardall, C. Y., Endeve, E., & Mezzacappa, A. 2013a, PhysicalReview D, 88, 023011—. 2013b, Physical Review D, 87, 103004Cernohorsky, J. 1994, Astrophysical Journal, 433, 247Cernohorsky, J., & Bludman, S. A. 1994, Astrophysical Journal,433, 250Chu, R., Endeve, E., Hauck, C. D., & Mezzacappa, A. 2019,Journal of Computational Physics, 389, 62 ,doi: https://doi.org/10.1016/j.jcp.2019.03.037Cockburn, B., Hou, S., & Shu, C.-W. 1990, Mathematics ofComputation, 54, 545Cockburn, B., Lin, S.-Y., & Shu, C.-W. 1989, Journal ofcomputational Physics, 84, 90Cockburn, B., & Shu, C.-W. 1989, Mathematics of computation,52, 411—. 1991, ESAIM: Mathematical Modelling and NumericalAnalysis, 25, 337—. 1998, Journal of Computational Physics, 141, 199Cockburn, B., & Shu, C.-W. 2001, Journal of ScientificComputing, 16, 173 Endeve, E., Hauck, C. D., Xing, Y., & Mezzacappa, A. 2015,Journal of Computational Physics, 287, 151Evans, C., Pollock, S., Rebholz, L. G., & Xiao, M. 2020, SIAMJournal on Numerical Analysis, 58, 788Hamilton, S., Berrill, M., Clarno, K., et al. 2016, Journal ofComputational Physics, 311, 241 ,doi: https://doi.org/10.1016/j.jcp.2016.02.012Hannestad, S., & Raffelt, G. 1998, Astrophysical Journal, 507, 339Hesthaven, J. S., & Warburton, T. 2008, Nodal discontinuousGalerkin methods: Algorithms, analysis and applications(Springer)Hix, W. R., Messer, O. E. B., Mezzacappa, A., et al. 2003, PhysicalReview Letters, 91, 201102Hix, W. R., Lentz, E. J., Endeve, E., et al. 2014, AIP Advances, 4,041013Hu, J., Shu, R., & Zhang, X. 2018, SIAM Journal on NumericalAnalysis, 56, 942Janka, H.-T. 2012, Annual Review of Nuclear and Particle Science,62, 407Juno, J., Hakim, A., TenBarge, J., Shi, E., & Dorland, W. 2018,Journal of Computational Physics, 353, 110Just, O., Bollig, R., Janka, H. T., et al. 2018, MNRAS, 481, 4786Just, O., Obergaulinger, M., & Janka, H.-T. 2015, MNRAS, 453,3386Kelley, C. 2018, Acta Numerica, 27Kl¨ockner, A., Warburton, T., Bridge, J., & Hesthaven, J. S. 2009,Journal of Computational Physics, 228, 7863Knoll, D. A., & Keyes, D. E. 2004, Journal of ComputationalPhysics, 193, 357Kuroda, T., Takiwaki, T., & Kotake, K. 2016, AstrophysicalJournal Supplement Series, 222, 20Laiu, M. P., Chen, Z., & Hauck, C. D. 2020, Journal ofComputational Physics, 417, 109567,doi: https://doi.org/10.1016/j.jcp.2020.109567Larecki, W., & Banach, Z. 2011, JQSRT, 112, 2486Larsen, E. W., & Morel, J. E. 1989, Journal of ComputationalPhysics, 83, 212Lentz, E. J., Mezzacappa, A., Messer, O. E. B., Hix, W. R., &Bruenn, S. W. 2012a, Astrophysical Journal, 760, 94Lentz, E. J., Mezzacappa, A., Messer, O. E. B., et al. 2012b,Astrophysical Journal, 747, 73Lentz, E. J., Bruenn, S. W., Hix, W. R., et al. 2015, AstrophysicalJournal Letters, 807, L31Levermore, C. D. 1984, JQSRT, 31, 149Liebend¨orfer, M., Messer, O. E. B., Mezzacappa, A., et al. 2004,Astrophysical Journal Supplement Series, 150, 263Liebend¨orfer, M., Rampp, M., Janka, H.-T., & Mezzacappa, A.2005, Astrophysical Journal, 620, 840
AIU ET AL ..