Free Access
Volume 586, February 2016
Article Number A82
Number of page(s) 13
Section Numerical methods and codes
Published online 28 January 2016

© ESO, 2016

1. Introduction

Astrochemistry deals with chemical reactions in astrophysical gases and, taking a broader view, embraces ionizing processes, dust-molecule interactions, and photo-chemical effects. Its significance in certain environments results from the fact that thermal and magnetic properties are determined by species abundances and the ionization state which, both, are regulated by chemical processes. Through diverse cooling and heating mechanisms entering the energy budget, and by controlling the coupling between magnetic field and matter, chemistry thus also influences the dynamical evolution of the system. The interstellar medium with its different prevailing phases of diffuse gas, molecular clouds, and dense gas cores is an environment where rich chemical processes take place. In the densest gas aggregates chemistry potentially manipulates fragmentation and the following gravitational collapse of pre-stellar cores to new stars. Further importance of interstellar chemistry comes from the role certain molecules play as tracers in observations such as the CO molecule for molecular gas detection or the NH3 molecule as an indicator of the gas temperature.

In order to bring magneto-gasdynamical simulations and observations closer, astrochemistry must find its way into numerical models. However, embedding chemistry in time-dependent, multidimensional computations is a very challenging task. A major difficulty arises from a possible large disparity between dynamical and reactive timescales in the flow. Since the overall complexity of chemo-magneto-gasdynamical simulations requires split integration techniques a sufficiently good coupling between dynamics and chemistry is necessary for accuracy and stability reasons. In practice, this forces a reduction in the size of the time-step, which can be used in advancing the solution since chemical kinetics may involve smaller timescales than the dynamical processes. In addition, the computational cost per numerical cell for solving the chemical rate equations quickly dominates the magneto-gasdynamical part for larger chemical networks. In consequence, this often prohibits the use of complex chemical networks in high-resolution simulations. What makes chemical kinetics so expensive is the need for sophisticated ordinary differential equation (ODE) solvers that can account for the stiffness inherent to the system. This comes about because individual reactions can proceed with very different speed. Cheap explicit schemes are ruled out because of their strict stability requirements which would severely limit the time-step according to the fastest reaction speed. A further complication arises because chemical kinetics and their associated thermal processes are often closely connected with regard to their time evolution. A treatment of both effects within a common numerical framework seems therefore reasonable.

The problem of incorporating chemistry and related thermal processes in astrophysical simulations has attracted substantial attention in the past. Early attempts to combine chemistry with time-dependent gasdynamics were undertaken by Lim et al. (1999) in 2D and by Abel et al. (2002) in 3D using the ENZO code (Bryan et al. 2014). In a work by Glover & Mac Low (2007), who adopt the ZEUS code (Stone & Norman 1992), a rudimentary hydrogen network was evolved consistently with the equations of magneto-gasdynamics in order to study molecular cloud formation. Nonequilibrium models describing the ionization of an atomic multispecies gas with cooling in optically thin plasma have also been implemented into magneto-gasdynamical codes like ASTROBEAR (Poludnenko et al. 2005) or PLUTO (Teşileanu et al. 2008). Moreover, there are powerful 1D time-dependent codes that integrate magneto-gasdynamics alongside more complex chemical networks used, e.g., in high-resolution studies of magnetic shocks (Lesaffre et al. 2004) or diffusion fronts (Lesaffre et al. 2007). Stand-alone astrochemistry packages (not coupled to dynamics) have been designed for special purposes like the ALCHEMIC code (Semenov et al. 2010) for the chemistry in proto-planetary discs, the NAHOON code (Wakelam et al. 2012) for the study of the chemical evolution in molecular clouds, or the MEUDONPDR code (Le Petit et al. 2006) that couples chemistry with radiation transfer suitable for radiation-dominated regions. Recently, a publicly available astrochemistry package called KROME has been developed (Grassi et al. 2014) based on the DLSODES solver from the Hindmarsh ODE tool (Hindmarsh 1983). This code aims at a more general-purpose chemistry framework for the astrophysics community, and claims to be easily includable in any fluid dynamics simulation code which provides a multispecies functionality. KROME comes with a variety of predefined chemical networks and is prepared for a treatment of chemistry-related thermal processes, photochemistry and dust evolution.

This paper describes the implementation of a chemical reaction network solver into the astrophysics code NIRVANA1 (Ziegler 2008, 2011, 2012). The chemistry module is constructed from scratch and adapted to the given code infrastructure. This proprietary coding effort allows full control over all aspects of implementation, and this was the main reason for not incorporating external software like KROME.

2. Chemical reaction network

A set of Nr chemical reactions of the form (1)between species is considered where, in each reaction r, a number of reactants XsR(r) are transformed into a number of products XsP(r) with R(r) and P(r) being subsets of { Xs }. The reaction speed is determined by the kinetic rate coefficient kr(T), which is assumed to be a function of temperature T. The α’s and β’s are stoichiometric coefficients. Chemical species can be atoms, molecules, positively or negatively charged ions, and free electrons. The total number of species is Ns + 1 with index s = Ns + 1 reserved for electrons. A species Xs is associated with a number density ns (expressed in m-3), which is a function of space and time, a molecular weight μs and a charge qs. At the current stage of development reactions including dust particles and photochemical processes are not considered. In principle, both effects could be treated within the framework presented here but would lead to a significant increase in complexity. It is naturally assumed that changes in the electron density are exclusively caused by changes in the ionic state of other species through ionization/recombination processes, electron attachment/detachment or charge transfer mechanisms. In this case, the electron density ne (=nNs + 1) need not be treated as an independent variable but is given by (2)expressing total charge conservation2. Here and in the following, an unspecified summation index s is assumed to run from 1 to Ns (i.e., excluding free electrons).

The chemical network module of NIRVANA is designed to handle an arbitrary number of species with an arbitrary number of reactions involved allowing an arbitrary number of reactants and products. Technically, individual reactions are easily specified as a single line in a text file obeying a generic input format of the form ID S R [& S R [& ...]] > S P [& S P [& ...]], where ID (a running number) uniquely identifies the reaction, S are stoichiometric coefficients (the α’s and β’s), R stands for reactant, and P stands for product. For example, for the reaction H2+ e 2 H + e being number twelve in a network the format would be 12 1 H2 & 1 e- > 2 H & 1 e-. By parsing the “reactions file” the code checks for the mass and charge balance in each reaction, and automatically derives and saves the structure of the reaction term (4) below which enters the ODE solver for chemical kinetics. A by-hand definition of the ODE system is redundant. Since reaction rates kr can be complicated functions of temperature not necessarily fitting a prescribed generic expression, those rates are specified separately in order to permit the highest flexibility.

As a basis in this work, a simple network for a hydrogen–helium plasma was constructed to be used in some of the presented test problems. It consists of the 19 reactions cataloged in Table A.1 for the nine species H, H+, He, He+, He2+, H, H2, H, and free electrons. The network includes ionization/recombination of H, He and a channel for H2 formation via the H, H ions.

2.1. Advection-reaction system

In gasdynamical simulations the time evolution of species number densities is described by a coupled system of advection–reaction equations, (3)with reaction term (4)subject to constraint (2). The product on the right-hand side running over reactants in a reaction makes the system nonlinear. Equation (3) has to be supplemented, of course, by other fluid equations describing the evolution of the flow speed v and pressure p (respective energy density). Maxwell’s equations add in magneto-gasdynamical problems involving the evolution of a magnetic field B. Pressure, temperature and chemical composition of the gas are connected via the ideal equation of state

with ntot = ∑ sns + ne the total number density, kB the Boltzmann constant and mu the atomic mass unit. Taking constraint (2) into account, the gas density ϱ is linked to species densities by the relation (5)with and the electron weight μe = 5.49 × 10-4. The mean molecular weight and electron number fraction Y is defined by

Usually, in gasdynamical simulations a continuity equation for ϱ expressing mass conservation is evolved in time: (6)Multiplying Eq. (3) by and summing over species results in (7)Hence, consistency with the continuity equation requires (8)expressing mass conservation in each reaction.

2.2. Thermal processes

Chemical reactions in astrophysical gases are often accompanied by special thermal processes. For instance, collisional ionization of hydrogen leads to a cooling of the gas because the energy needed to pass the hydrogen atom from its neutral to ionized stage is extracted from the thermal energy. Radiative recombination of ionized hydrogen to its neutral stage also leads to cooling of the gas because a free electron is captured, which removes part of its kinetic energy. The formation of H2, on the other hand, is often an exothermic process that heats the gas. Macroscopically, cooling and heating processes are taken into account by a source term Se in the thermal energy equation (and total energy equation). Ignoring dissipative processes (9)holds for the thermal energy density e = p/ (γ−1) with adiabatic index γ assumed to be a constant value3. Since the full system of chemo-magneto-gasdynamical equations is solved in an operator-split fashion, from Eqs. (3) and (9) the thermo-reactive part in the solution procedure (ϱ, v, and B are assumed to be constant during this update) is governed by

Using the equation of state, the expression for the total number density and Eq. (10) one has

It follows for the evolution of temperature (12)which is used instead of (11) in the operator-split formalism presented below.

The source term Se, which generally depends on the species densities and temperature, is decomposed here into two parts according to (13)The first (canonical) term on the right-hand side has the same structural form as the source term Rs in the advection-reaction equation for ns with a thermal rate coefficient lr(T) instead of kinetic rate kr(T) associated with a reaction r. Many thermal processes fit this canonical form with lrkr in some cases. In order to offer a high degree of flexibility other processes not describable within the first term of the energy source are subsumed in the second term Le(n,T). Whenever possible, however, thermal terms should be treated within the canonical source term simply for reasons of computational efficiency: derivatives of products with respect to ns (needed to compute the Jacobian matrix) must already be evaluated for the source terms Rs and so can be reused.

A summary of the thermal processes that have been considered in combination with the chemical network in Table A.1 is presented in Table B.1. A more detailed description of individual thermal processes can also be found in Appendix B, which includes examples for Le-terms.

3. Numerical implementation

The following presentation concentrates on implementation details related to the introduction of reactive species. Other solution parts of the full system of chemo-magneto-gasdynamical equations are not described here. In this coverage numerical integration of the coupled system of advection-reaction Eq. (3) and temperature Eq. (12) must be carried out with respect to co-existing different timescales inherent to flow advection on the one hand and chemical kinetics on the other hand. Moreover, chemical reactions and related thermal processes may, on their own, proceed with disparate speeds introducing stiffness to the system and requiring robust and efficient solvers. Overall, the system is divided into an advection part described by the PDE (14)which is solved simultaneously with the advection of density, momentum equation, and energy equation with Se = 0, and a source step given by the ODE system (15)where . In order to advance the solution one time-step Δtn from tn to tn + 1 = tn + Δtn a second-order accurate Strang-type operator splitting is used. Denoting the advection and source step operators as and the update formally reads

The actual size of time-step is given by

where is the advection/dynamical time-step dictated by the usual CFL condition for numerical stability of the discretized gasdynamical equations and is a reaction time-step defined in more detail later in this section. A proper choice of the latter serves to control the coupling between advection and reaction. Generally, adopting as sole time-step may be insufficient because large splitting errors can lead to unphysical behavior or even destabilize numerics in situations where .

The advection of individual species is performed within the Finite-Volume (FV) framework of the NIRVANA code (Ziegler 2011). The numerical procedure is analogous to the advection of the gas density ϱ. In the FV scheme the update of volume-averaged number densities requires computing fluxes at each grid cell surface, which depend on a properly upwinded state . This is derived in accordance with other fluid variables from the approximate solution of a Riemann problem with left-sided/right-sided states at grid cell surfaces. These are obtained by a piecewise linear reconstruction procedure based on volume-averaged values in adjacent cells with slopes subject to the TVD limiter of van Leer (1977). However, since slope limiting is a nonlinear process and species profiles are reconstructed independently from each other, the consistency condition

is per se not fulfilled. If left uncorrelated, deviations between the total density and sum of partial densities ϱs can accumulate with time and may produce unphysical solutions. Hence, some modification in the advection algorithm is advisable in order to meet the conservation principle4. The simplest way to deal with the problem is to overwrite ϱL,R and to replace the corresponding total mass flux by the sum of partial mass fluxes. Although may be small in size, with this solution there is the danger of violating monotonicity and relaxing TVD constraints in the evolution of ϱ. This, in turn, may result in overshooting effects near local extrema or at discontinuities. An alternative, which is applied here, is profile flattening by reducing species density slopes, σs, while leaving the evolution of total density almost untouched. Slope modification should be carried out in a clever way to diminish overshooting, to maintain TVD properties where possible, and to retain second-order accuracy. Concretely, the individual species slopes in each spatial direction and in each numerical cell are modified according to


0 ≤ χs ≤ 1, with the slope sum σ = ∑ σs and with σϱ the slope in the ϱ-reconstruction. This choice of χs allows for an easy interpretation: First we note that flattening is selective, i.e., if σ> 0 (σ< 0) only components with positive (negative) slope are changed by the same factor, whereas species components with opposite slope are unaffected (χs = 1). Moreover, if sgn σ = sgn σϱ and | σ | > | σϱ | applying selective flattening recovers consistency in the sense that . If | σ | < | σϱ | no flattening is necessary and χs = 1. On the other hand if sgn σ ≠ sgn σϱ, the total density slope can be expected to be close to zero. Then, the slope is set σϱ ≡ 0 followed by species flattening as before, which also makes . We note that this does not generally mean setting , which would imply reduction to a first-order algorithm. The choice of χs may be thought of as a best choice introducing a minimal amount of dissipation where needed in order to let | σ | converge to | σϱ | from above. By contrast, convergence to | σϱ | from below is suppressed, which would require profile steepening. Finally, number density fluxes obtained on the basis of the modified reconstruction with slopes are rescaled to new values according to

which ensures conservation, i.e., where F(ϱ) is the total density mass flux. As a side note, in the similarly conservative multispecies FV advection scheme of Plewa & Müller (1999) the principle of selective flattening is applied to the upwinded cell interface values followed by a rescale of species fluxes. This is an equally valid approach applied to their third-order PPM scheme.

In summary, with updating ϱ mass is conserved up to machine round-off error accumulation as usual, whereas consistency of the gas species with ϱ in a FV sense is achieved, i.e., Eq. (5) – when volume-averaged over a grid cell – holds at tn + 1 when fulfilled at tn.

For the source step integration a stiff ODE solver is needed. There is a wide variety of such solvers available in the mathematical literature which vary in accuracy and computational complexity from the simple first-order implicit Euler scheme to very precise high-order Peer methods. I decided to implement the Rosenbrock-Wanner method named GRK4(3)T developed in Kaps & Rentrop (1979). It is a fourth-order generalized Runge-Kutta method belonging to the class of semi-implicit schemes. This choice of solver was guided by the following thoughts:

  • (1)

    Since the ODE solver is not envisaged as a stand-alone tool butis used operator-split together with second-order FVgasdynamics moderate accuracy requirements are sufficient,i.e., a fourth-order scheme like GRK4(3)T is convenient.

  • (2)

    Since the time-step delivered to the ODE solver module can be much larger than the smallest timescale in the reaction system, the solver should have an efficient internal step-size control to meet accuracy requirements. GRK4(3)T embeds a third-order solution approximation almost for free, which can be used in a standard error estimator with user-prescribed relative (RTOL) and absolute (ATOL) error tolerances (see, e.g., Strehmel et al. 2012). In short, if u(4) and u(3) denote the fourth-order and third-order approximations (see Table 1), an estimate for the normalized error is given by If ϵ ≤ 1 the integration step with internal step-size h is accepted, whereas for ϵ> 1 it is rejected and is repeated with smaller step-size. In both cases the choice of a new step-size is error-controlled and adjusted according to where p = 4 is the order of the scheme. All species components have the same value of ATOL whereas it is different, obviously, for the temperature component.

  • (3)

    The species advection step was designed consistently with the gas density advection to ensure mass conservation. In order to maintain this property for the combined solver, the ODE method must conserve total species mass during reactions, i.e., recalling (7), Our implementation of GRK4(3)T fulfills this condition to machine round-off5.

For interested readers, GRK4(3)T is sumarized in Table 1. It needs computation of the Jacobi matrix J = S/u and the solution of a linear equation system of dimension Ns + 1 (number of species plus temperature) for four different right-hand side vectors corresponding to a four-stage procedure. At this stage of development, standard LU factorization is applied which is an (Ns + 1)3 operation and so becomes rapidly expensive for larger astrochemical networks. Costs may be reduced by utilizing a sparse structure of J inherent to astrochemistry. For several hundreds or thousands of species LU decomposition may still be much too expensive and, as a cheaper alternative, a Krylov-type solution approach could be applied. For the network in Table A.1 with Ns = 8 and Nr = 19, however, computation of J is the most expensive part and not the LU factorization. For efficiency reasons the Jacobian should be defined analytically, but the implementation also offers the possibility of an approximate evaluation by finite differencing if derivatives are too cumbersome.

In flow regions where the physical conditions are such that the ODE system is non-stiff, using an explicit method may be more efficient6. Therefore, a hybrid integration strategy was designed based on a measure of stiffness. The ODE system is advanced in time according to

where the norm7 of the Jacobi matrix is defined by

The explicit solver DOPRI5(4) is the well-known fifth-order accurate Runge-Kutta method of Dormand & Prince (1980). The value of two is roughly the stability limit of DOPRI5(4).

In order to better control the coupling between advection and reaction an a posteriori reaction time-step is introduced given by

It means that the follow-up reaction time-step is limited by a maximum allowed relative change in species number densities (ϵn) due to chemical kinetics and relative change in temperature (ϵT) due to thermal processes computed from the last reaction source update in the Strang-splitting procedure. The coupling parameters ϵn and ϵT are user-prescribed quantities. Naturally, the values ϵn< 1,ϵT< 1 should be chosen with smaller values giving a better coupling at the expense of possibly reducing the overall time-step and, hence, increasing computational costs.

Table 1

Four-stage GRK4(3)T method.

4. Validation and tests

4.1. Advection scheme validation

In order to validate the multispecies advection scheme proposed in this work an extended interacting blast wave problem was studied using the setup of Plewa & Müller (1999). The initial state is given by

on the interval x = [0,1] discretized by 400 equidistant cells and adopting reflecting boundary conditions on both ends. Additionally, three nonreactive species components are introduced that have initial mass fractions of as illustrated in Fig. 1 (left plot). An ideal equation of state with adiabatic index γ = 1.4 is used.

Figures 1 (right plot) and 2 illustrate the results at the final time t = 0.038. Figure 2 shows the solution for the hydrodynamic variables ϱ, p, and v. The spatial distribution of mass fractions in Fig. 1 agree well with those in Plewa & Müller (1999) (cf. Fig. 5) although their third-order and less dissipative PPM algorithm produces somewhat higher-amplitude variations in the region . Moreover, no artificial discontinuity is created in the components f2 and f3 at x ≈ 0.6, which is visible in their solution obtained with the preferred sCMA method. As expected, the deviation from Eq. (5) over the integration time remains well bounded. Expressed in numbers for the mass fractions measured on the grid for the maximum and spatially averaged error, and clearly demonstrating the conservative nature of the multispecies advection scheme.

thumbnail Fig. 1

Spatial distribution of initial (left) and final (right) mass fractions (f1 – solid line; f2 – dashed line; f3 – dash-dotted line) in the interacting blast wave problem.

thumbnail Fig. 2

Spatial distribution of density (top), pressure (middle), and velocity (bottom) at t = 0.038 for the multispecies interacting blast wave simulation.

4.2. Stiff solver validation

The accuracy of the GRK4(3)T solver and its correct implementation into the NIRVANA code was checked by performing several stiff problems. Here, the well-known problem named OREGO is presented. The OREGO problem is described in Hairer & Wanner (1996) and was originally developed by Field et al. (1972) as a simplified mathematical model for the oscillating behavior in the Belousov–Zhabotinskii chemical reaction (Gray 2002).

The problem is written in the form for concentrations n1,n2,n3 with rate coefficients Initial concentrations, obtained concentrations with GRK4(3)T at the final integration time t = 360, and high-precision reference values taken from the literature are given in Table 2. The numerical solution is illustrated in Fig. 3 demonstrating its oscillating character and the occurrence of very sharp variations over short time intervals. The maximum relative error obtained with GRK4(3)T is on the order of 5 × 10-4 in the variable n3 adopting numerical parameters RTOL= 10-4 and ATOL= 10-14.

Table 2

Concentrations in the OREGO problem.

thumbnail Fig. 3

Numerical solution of the OREGO problem (with logarithmic scaling for n2).

4.3. Advective OREGO

This challenging test problem combines the previous OREGO problem with an underlying advection in a 1D setup. Choosing a constant velocity of v = 1 the OREGO equations in Sect. 4.2 are solved together with Eq. (14) on a domain x = [0,400] with initial conditions n0(x,t = 0) = nOREGO(t = 0) from Table 2. Outflow conditions are adopted on the right side of the domain and inflow conditions on the left side of the domain where concentrations are fixed at their initial values. The system is evolved up to time t = 360 corresponding to the integration time in the pure OREGO problem. The exact solution at the final time is given by n(x,t = 360) = nOREGO(x/v) for x ≤ 360 and a constant continuation for x> 360, i.e., the spatial profiles of concentrations look identical to the temporal variation in the pure OREGO problem. We note that since the OREGO equations represent a nonconservative system (ss ≠ 0), this problem cannot serve as a consistency check of the multispecies advection scheme because the sum of concentrations does not obey a continuity equation.

The obtained numerical solution within the second-order operator-split framework with second-order advection scheme and GRK4(3)T solver is shown in Fig. 4 for two different grids with 400 and 1200 grid cells. A significant dependence on the numerical resolution is observed due to the presence of very sharp, thin structures advected with the flow. Numerical dissipation tends to broaden features and therefore to lower peak values. The largest of these errors are visible at the minimum peaks in variable n2 in the lower resolved run. In this case insufficient resolution also causes a dispersion-like error noticeable as a moderate leftward shift of the right-hand feature. This dispersion error results from the interaction of numerical dissipation errors with the update of rate equations in virtue of the operator-split approximation. The more highly resolved case produces results much closer to the reference solution. The outcome of both low- and high-resolution runs is not sensitive to the choice of coupling parameter ϵn in the expression for Δtreac because there is no back-reaction of the species evolution on the flow.

thumbnail Fig. 4

Numerical solution of the OREGO problem with underlying advection for different spatial resolution (solid: 1200 grid cells; dashed: 400 grid cells).

4.4. Ionization equilibrium

The relaxation of a hydrogen-helium plasma with species H, H+, He, He+, He2+ in ionization equilibrium is considered. Assuming a fixed temperature of 105 K and initial number densities far from equilibrium for this temperature, namely,

measured in units m-3, the ionization/recombination network consisting of reactions 1–6 in Table A.1 is solved. The evolution of the network is followed up to a maximum time of tmax = 106 yr which is more than two orders of magnitude larger then the relaxation time for the ionization of He+.

In equilibrium (t → ∞), where H = He = H+ = He+ = He2 + = 0, the number densities are given by

with rate coefficients ki as specified in Table A.1. Equilibrium values together with initial values and obtained numerical values are summarized in Table 3. The numerical solution of the relaxation process is shown in Fig. 5. Ionization of H takes place on a timescale of 50 yr followed by ionization of He on a timescale of 100 yr. Finally, built-up He+ is transformed into He2+ in a few thousand years. At tmax, equilibrium balance is reached with relative deviation less than 10-8 from exact equilibrium values.

Table 3

Number densities in the relaxation problem.

thumbnail Fig. 5

Numerical solution of the ionization/recombination relaxation problem. Number densities of all species are shown as a function of time.

4.5. One-dimensional Riemann problem

As a first test in combining gasdynamics with chemistry, a 1D shock tube problem – a modified Sod problem – is considered. The chemical network as cataloged in Table A.1 is adopted accompanied by all the thermal processes summarized in Table B.1. The chemistry model is that used in Grassi et al. (2014) who performed the same Riemann problem to check their implementation of KROME in conjunction with the astrophysics codes RAMSES (Teyssier 2002) and FLASH (Fryxell et al. 2000).

The initial conditions are set up on a computational domain x = [0,1 pc] discretized by 1000 equidistant grid cells. Values on the regions left and right of the interface at x = 0.5 pc are given by

where is the initial mass fraction of the different species X (see Table 4).

Table 4

Initial mass fractions of species in the Riemann problem.

The outcome of this test is illustrated in Figs. 6 and 7. Figure 6 shows gasdynamical quantities and Fig. 7 shows the spatial variation of the mass fraction of species H2, H, and e. In order to allow a qualitative comparison8 with the results of Grassi et al. (2014) the simulation was stopped at canonical time t = 2.45 × 104 yr (see their Figs. 19–22). Differences from the pure gasdynamical case without chemistry are mainly due to the chemical evolution of H2 and H2-related cooling mechanisms which are the dominant thermal processes. It leads to a significant temperature reduction in regions of relatively high H2 mass fraction. Finally, the overall characteristics of chemo-thermal evolution obtained with NIRVANA agree well with the behavior observed in Grassi et al. (2014)9.

thumbnail Fig. 6

Gas density (top), temperature (middle), and velocity (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry. The dashed line corresponds to the solution of a standard Sod problem without chemistry.

thumbnail Fig. 7

Mass fraction of species H2 (top), H (middle), and e (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry.

4.6. One-dimensional propagating shock

As a second example combining magneto-gasdynamics with chemistry a propagating shock problem is presented (see Teşileanu et al. 2008 and Massaglia et al. 2005). The background medium is a stratified gas with constant temperature T0 = 1000 K, zero velocity v0 = 0, and constant (perpendicular) magnetic field By0 = 10-8 T. The background density profile is given by

where ϱ00 = 1011mu m-3 and x0 = 1.8 × 1012 m. The computational domain ranges from x = 0 to x = 4 × 1013 m assuming outflow boundary conditions on both sides and adopting an effective spatial resolution10 of δx = 1.25 × 109 m. The background medium is superimposed on a localized adiabatic perturbation that steepens into a single shock propagating down the stratification. In order to eliminate the reverse shock, the disturbance is constructed such that the non-magnetic Riemann invariant

is fulfilled with γ assumed 5/3. At t = 0 this results in a structure (background plus perturbation) given by with computed for the background medium. The velocity perturbation v is prescribed as

where the peak value , σ = 2 × 1011 m and x1 = 1012 m. We note that exterior to the interval (x1,x1 + 2σ) the configuration reduces to the isothermal background.

The gas is composed of species X with initial mass fractions as summarized in Table 5. The chemical network in Table A.1 with thermal processes in Table B.1 is employed. Three different simulations are performed. The first includes the full chemistry model and all thermal processes, i.e., energy losses by ionization/recombination and line emission of H and He as well as H2-related cooling/heating mechanisms. In the second simulation H2 is treated as a passive scalar, H2 chemistry is switched off, and there is no H2 cooling/heating. It serves to isolate the effects of H2 chemistry which are expected to be important in virtue of the relatively high initial H2 mass fraction of and chosen temperature T0. This simulation is, at least qualitatively, comparable with the SNEq model in Teşileanu et al. (2008) despite the absence of metal cooling here. The third simulation is adiabatic with all species passively advected and no chemical reactions at all. In the simulations with a chemical reaction network included, the reaction time–step parameters are chosen ϵn = ϵT = 0.1, which adjusts it to any rate that is on the order of the dynamical time–step.

The results are shown in Fig. 8 for selected quantities at a time t = 13 yr where the simulations have been stopped. Evidently, the evolution of the shock and the structure of the post-shock region is crucially affected by the chemistry and cooling processes when compared with the adiabatic case. The shock speed is substantially reduced and the temperature drops faster behind the shock. Shock heating raises the temperature well above 104 K, which partially ionizes the region with extent 3 × 1012 m right behind the shock front. Also, the magnetic field reaches its maximum in this region. The influence of H2 chemistry is to provide further cooling in regions where the number density of molecular hydrogen is sufficiently high. This is all of the region up to x ≈ 7 × 1012 m and around the peak in molecular hydrogen density at x ≈ 1.3 × 1013 m where a thin shell has formed in the close vicinity of the left edge of ionized region. The shell lies within a broader region where the H2 fraction otherwise was substantially reduced by collisional dissociation in passing the shock front – up to three orders of magnitude compared to the initial mass fraction. This feature is also clearly seen in the temperature that has dropped to 100 K there. The pre-shock region is at most marginally affected by cooling because of low effectiveness. Overall, the striking similarity in the dynamical evolution and spatial structure of the simulation without H2 chemistry/cooling with the SNEq model in Teşileanu et al. (2008) further strengthens confidence in the validity of the presented algorithmic approach.

Table 5

Initial mass fractions of species in the propagating shock problem.

thumbnail Fig. 8

From top to bottom: gas density, electron number fraction, H2 mass fraction, temperature, velocity and magnetic field at t = 13 yr (solid line: full chemistry; dashed line: no H2 chemistry; dotted line: no chemistry, adiabatic). The thin solid line is the initial profile.

5. Conclusions

A chemical reaction network solver was successfully implemented into the astrophysics code NIRVANA. For that purpose a numerical framework was devised that accounts for the accurate treatment of the set of advection-reaction equations for chemical species. These equations add to the already complex equations of magneto-gasdynamics. In particular, the ODE system built from the rate equations representing chemical kinetics and the temporal temperature change due to chemistry-connected thermal processes is solved coherently with a high-order stiff integrator of linear-implicit type. Another feature of the implementation is that it evolves chemical species consistently with the total gas density. This nice property is achieved by modifications in the advection part of the FV method and by a proper choice of ODE solver for the reaction part. It means that the mass conservation principle, which is usually fulfilled in a FV discretization of the gasdynamics equations, remains valid here for the case of reactive flows.

A series of tests have provided confidence in future applications of the code involving chemistry and complex thermal processes. However, computational costs for multidimensional simulations are a limiting factor in practice and, this time, prohibit the use of chemical networks that are too large. There is otherwise some room for optimizations of the employed ODE solver, for example by utilizing the sparse matrix structure to accelerate LU decomposition and/or by applying a Krylov-type approximation in order to reduce costs further. Such possible improvements are the subject of future work.

Complementary to the numerical framework an easy-to-use interface was developed which allows straightforward upgrades of the H-He chemical network that was constructed as a standard network in this work. Adding new species and reactions is performed just by editing a specifically formatted text file. By parsing that file the ODE system is generated automatically, i.e, the ODE system need not be given explicitly. There is, in principle, no restriction to the number of species, the number of reactions, and the number of reactants in a reaction. On the other hand, kinetic and thermal rate coefficients for each chemical reaction as well as special thermal processes have to be defined in separate modules. This allows the highest degree of flexibility regarding their dependence on temperature and, possibly, gas composition.

All in all, with the numerical and technical framework installed in this work the NIRVANA code has become a very powerful tool for time-dependent chemo-magneto-gasdynamical simulations. In combination with the code’s self-gravity solver, its efficient solver for dissipation processes and with its adaptive mesh refinement capability, challenging astrophysical problems seem attackable in the future limited more by computational costs than by physics.


There may be more algebraic constraints between species inherent to (1) resulting, e.g., from the elemental conservation principle or chemical equilibrium assumptions for some species. The presented numerical solver does not yet account for such general situations.


More generally, the adiabatic index depends on the gas composition and temperature.


The case of an adaptive grid, which requires some additional effort, is not discussed here.


Because it is a special linear-implicit type of solver, this feature stems from the observation that a cascade of linear systems for intermediate states is solved exactly with a matrix based on the Jacobian, and where in each stage the condition is fulfilled for the intermediate state.


This is not necessarily true if the LU factorization is inexpensive and J is to be computed anyway to check stiffness.


This is an incomplete maximum norm composed of the maximum norm of the chemical kinetics subsystem and an inverse temperature timescale due to thermal processes (JNs + 1,Ns + 1-element). Use of the full maximum norm is not suitable here because matrix elements do not all have the same physical dimension.


An elaborate comparison is not possible because some details of initial conditions are missing in their description of the problem.


In Grassi et al. (2014) their temperature Eq. (61) misses the term tot/ntot which is, however, very small in this test problem.


In order to speed up the simulation an adaptive mesh with a 4000 equidistant-cell base grid and three additional levels of refinement was used (32 000 effective grid). This is a slightly better resolution than in Teşileanu et al. (2008), but not as good as in Massaglia et al. (2005).


I would like to thank the referee, Pierre Lesaffre, for valuable comments useful in improving the paper.


  1. Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New Astron., 2, 181 [Google Scholar]
  2. Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  3. Black, J. H. 1981, MNRAS, 197, 553 [NASA ADS] [Google Scholar]
  4. Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19 [NASA ADS] [CrossRef] [Google Scholar]
  5. Cen, R. 1992, ApJS, 78, 341 [NASA ADS] [CrossRef] [Google Scholar]
  6. Dalgarno, A., & Lepp, S. 1987, Astrochemistry (Dordrecht: Reidel) [Google Scholar]
  7. Donahue, M., & Shull, J. M. 1991, ApJ, 383, 511 [NASA ADS] [CrossRef] [Google Scholar]
  8. Dormand, J. R., & Prince, P. J. 1980, J. Comput. Appl. Math., 6, 19 [CrossRef] [MathSciNet] [Google Scholar]
  9. Dove, J. E., & Mandy, M. E. 1986, ApJ, 311, L93 [NASA ADS] [CrossRef] [Google Scholar]
  10. Field, R. J., Körös, E., & Noyes, R. M. 1972, J. Am. Soc., 94, 8649 [CrossRef] [Google Scholar]
  11. Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
  12. Glover, S. C. O., & Abel, T. 2008, MNRAS, 388, 1627 [NASA ADS] [CrossRef] [Google Scholar]
  13. Glover, S. C. O., & Mac Low, M.-M. 2007, ApJ, 659, 1317 [NASA ADS] [CrossRef] [Google Scholar]
  14. Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  15. Gray, C. 2002, Rose-Hulman Undergraduate Mathematics Journal, 3, 1 [Google Scholar]
  16. Hairer, E., & Wanner, G. 1996, Solving Ordinary Differential Equations II. Stiff and Differential-algebraic Problems (Springer Verlag) [Google Scholar]
  17. Hindmarsch, A. C. 1983, IMACS Trans. Sci. Comput., 1, 55 [Google Scholar]
  18. Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555 [NASA ADS] [CrossRef] [Google Scholar]
  19. Janev, R. K., Langer, W. D., & Evans, K. 1987, Elementary Processes in Hydrogen–Helium Plasmas (Berlin: Springer Verlag) [Google Scholar]
  20. Kaps, P., & Rentrop, P. 1978, Numer. Math., 33, 55 [CrossRef] [Google Scholar]
  21. Karpas, Z., Anicich, V., & Huntress, W. T. 1979, J. Chem. Phys., 70, 2877 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lesaffre, P., Chièze, J.-P., Cabrit, S., & Pineau des Forêts, G. 2004, A&A, 427, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Lesaffre, P., Gerin, M., & Hennebelle, P. 2007, A&A, 469, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Le Petit, F., Nehme-, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Lim, A. J., Rawlings, J. M. C., & Williams, D. A. 1999, MNRAS, 308, 1126 [NASA ADS] [CrossRef] [Google Scholar]
  26. Massaglia, S., Mignone, A., & Bodo, G. 2005, A&A, 442, 549 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Omukai, K. 2000, ApJ, 534, 809 [NASA ADS] [CrossRef] [Google Scholar]
  28. Plewa, T., & Müller, E. 1999, A&A, 342, 179 [NASA ADS] [Google Scholar]
  29. Poludnenko, A., Varniere, P., Frank, A., & Mitran, S. 2005, Springer’s Lecture Notes in Computational Science and Engineering, 41 [Google Scholar]
  30. Poulaert, G., Brouillard, F., Claeys, W., McGowan, J. W., & Van Wassenhove, G. 1978, J. Phys. B, 11, L671 [NASA ADS] [CrossRef] [Google Scholar]
  31. Ripamonti, E., & Abel, T. 2004, MNRAS, 348, 1019 [NASA ADS] [CrossRef] [Google Scholar]
  32. Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Shapiro, P. R., & Kang, H. 1987, ApJ, 318, 32 [NASA ADS] [CrossRef] [Google Scholar]
  34. Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791 [Google Scholar]
  35. Strehmel, K., Weiner, R., & Podhaisky, H. 2012, Numerik gewöhnliche Differentialgleichungen (Wiesbaden: Vieweg+Teubner Verlag) [Google Scholar]
  36. Teşileanu, O., Mignone, A., & Massaglia, S. 2008, A&A, 488, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. van Leer, B. 1977, J. Comput. Phys., 23, 276 [Google Scholar]
  39. Wakelam, V., Herbst, E., Loison, J.-C., et al. 2012, ApJS, 199, 21 [NASA ADS] [CrossRef] [Google Scholar]
  40. Ziegler, U. 2008, Comput. Phys. Common., 179, 227 [Google Scholar]
  41. Ziegler, U. 2011, J. Comput. Phys., 230, 1035 [Google Scholar]
  42. Ziegler, U. 2012, SIAM J. Sci. Comput., 34, C102 [CrossRef] [Google Scholar]

Appendix A: Chemical network

Table A.1 catalogs the chemical network of a hydrogen-helium plasma used in this work. The temperature-dependent rate coefficients for collisional ionization, radiative, and di-electronic recombination processes for H and He are collected from Janev et al. (1987) and Cen (1992). Other rates are mostly adopted from Abel et al. (1997) where the same network was used in a cosmological context to reflect a primordial gas composition.

Table A.1

Hydrogen-helium network.

Appendix B: Cooling and heating processes

All thermal processes that have been considered in conjunction with the chemical network in Table A.1 are summarized in Table B.1.

Appendix B.1: Hydrogen and helium cooling processes

Following Cen (1992) a standard description for cooling processes in a hydrogen-helium plasma is adopted. It includes ionization cooling of H, He, and He+ by collisions with electrons; radiative recombination cooling of H+, He+, and He2+; di-electronic recombination cooling of He+; collisional ionization cooling of the metastable He(23S); and cooling due to collisional excitation of H (all atomic levels n), He (n = 2,3,4 triplets), and He+ (n = 2). Moreover, Bremsstrahlung for the ions H+, He+, and He2+ according to Black (1981) is taken into account. We note that the effective cooling rates for He(23S) and He (n = 2,3,4 triplets) cooling are , hence, they cannot be treated as canonical source terms in (13).

Appendix B.2: Molecular hydrogen cooling and formation heating

Table B.1

Thermal processes associated with the chemical network in Table A.1.

Cooling due to collisional excitation of ro-vibrational levels of molecular hydrogen is considered. The net energy loss per volume per time is written in the form

where ΛH2 denotes the cooling function for an optically thin environment, and τ(n) is an opacity correction term once the gas enters the optically thick regime. Here, according to Ripamonti & Abel (2004)

The cooling function ΛH2 measured in units [J/s] is given by

which, in the high-density limit (ΛH2,CE ≫ ΛH2,LTE), reduces to the LTE cooling function following the work of Hollenbach & McKee (1979). In the low-density limit, ΛH2,CE ≪ ΛH2,LTE, the cooling function of Glover & Abel (2008) is adopted. It represents a sum of contributions from collisional excitation of H2with collision partners Xk ={H, H+, He, H2,e}:

The logarithms of individual cooling rates are polynomial functions in log T, i.e., log ΛH2,k = ∑ lcl(log T)l, with fitting coefficients cl taken from Table 8 in Glover & Abel (2008). In this fitting approach no distinction is made between excitation of ortho–H2and para–H2; instead, an ortho–para ratio of 3:1 is assumed.

Collisional dissociation of molecular hydrogen described by reactions 12 and 13 in Table A.1 is an endothermic process that leads to gas cooling. According to Omukai (2000) the energy removed from the gas per dissociated molecule is 4.48eV for both reactions. The thermal rate coefficients for these reactions are in canonical form and are given by

with k12,13 the kinetic rate coefficients.

In contrast, molecular hydrogen formation via associative detachment (reaction 8) and charge transfer (reaction 10) are exothermic processes that partly heat the gas. According to Omukai (2000), the binding energy released per produced molecule is 3.53 eV (1.83 eV) for reaction 8 (10). The thermal rate coefficients are given by


according to Hollenbach & McKee (1979), approximates the fraction of liberated binding energy going into heating the gas. The factor fcr arises from the assumption that all of the binding energy is first put into rotational/vibrational excitation of H2which is then partly converted into heat via collisional de-excitation rather than radiated away. However, these heat sources only become important at high densities (n ≳ 1014 m-3, say) at which fcr approaches 1.

All Tables

Table 1

Four-stage GRK4(3)T method.

Table 2

Concentrations in the OREGO problem.

Table 3

Number densities in the relaxation problem.

Table 4

Initial mass fractions of species in the Riemann problem.

Table 5

Initial mass fractions of species in the propagating shock problem.

Table A.1

Hydrogen-helium network.

Table B.1

Thermal processes associated with the chemical network in Table A.1.

All Figures

thumbnail Fig. 1

Spatial distribution of initial (left) and final (right) mass fractions (f1 – solid line; f2 – dashed line; f3 – dash-dotted line) in the interacting blast wave problem.

In the text
thumbnail Fig. 2

Spatial distribution of density (top), pressure (middle), and velocity (bottom) at t = 0.038 for the multispecies interacting blast wave simulation.

In the text
thumbnail Fig. 3

Numerical solution of the OREGO problem (with logarithmic scaling for n2).

In the text
thumbnail Fig. 4

Numerical solution of the OREGO problem with underlying advection for different spatial resolution (solid: 1200 grid cells; dashed: 400 grid cells).

In the text
thumbnail Fig. 5

Numerical solution of the ionization/recombination relaxation problem. Number densities of all species are shown as a function of time.

In the text
thumbnail Fig. 6

Gas density (top), temperature (middle), and velocity (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry. The dashed line corresponds to the solution of a standard Sod problem without chemistry.

In the text
thumbnail Fig. 7

Mass fraction of species H2 (top), H (middle), and e (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry.

In the text
thumbnail Fig. 8

From top to bottom: gas density, electron number fraction, H2 mass fraction, temperature, velocity and magnetic field at t = 13 yr (solid line: full chemistry; dashed line: no H2 chemistry; dotted line: no chemistry, adiabatic). The thin solid line is the initial profile.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.