Add new components for radiation energy density

[chongchonghe]

To implement a photoionization network, we will need radiation energy density as a state components, call it state.rad. It has a RHS (ydot) that depends on T, state.xn, and state.rad, and it should evolve along with state.xn and T. However, it differs from state.xn in three different ways:

1. it does not have aion and zion;
2. it won't affect the EOS
3. it's not limited to 0 -- 1.

Is there an existing component that can work this way? Any suggestion on how I should go about adding this new component?

It seems to me that the auxiliary composition won't work, because it does not have a RHS and Jacobian, but I might be missing something.

[zingale]

can you write down here what ODEs you want to solve?

Right now we (usually) solve:

$$\frac{dX_k}{dt} = \dot{\omega}_k(\rho, T, X_j)$$

$$\frac{de}{dt} = \epsilon(\rho, T, X_j)$$

with $T$ found via the EOS, $T = T(\rho, e, X_k)$

[chongchonghe]

A general ODE for coupling radiation with thermochemistry would look like:

$$ \begin{aligned} &\frac{dE_g}{dt} = R_g(\rho, T, E_j, X_k) \\ &\frac{d X_i}{d t}=\dot{\omega}_i\left(\rho, T, E_j, X_k\right)\\ &\frac{d e}{d t}=\epsilon\left(\rho, T, E_j, X_k\right) \end{aligned} $$

and EOS: $T = T(\rho, e, X_k)$.

[markkrumholz]

@zingale, @BenWibking: @chongchonghe just had a long discussion, and we realized that this issue actually probing a much more fundamental design question, which is going to come up in many places: how should one deal with chemistry networks that have extra conservation laws? That is, right now microphysics just enforces a single conservation law, which is that the mass fraction of all species must add to unity. But for more general chemistry networks you may have additional stoichiometric constraints that should be obeyed. Most obvious is charge conservation, but for networks like ours that do not involve nuclear reactions you also have conservation of nuclear number, i.e., your reaction network isn't allowed to change the total number of atoms of hydrogen, carbon, etc.

The question is how to implement these extra constraints. One option is to track each species and implement the conservation laws on the RHS of the ODE -- so for example the way you would enforce conservation of H nuclei is by having an RHS that looks like

d(n_H0) / dt = (stuff)
d(n_H+) / dt = -(stuff).

The other option is only to follow the species that you actually have to, and then deduce the remaining abundances from conservation laws when they are needed; in this case you just have an RHS that looks like

d(n_H0) / dt = (stuff)

and any time you need to know the H+ number density you compute n_H+ = n_H - n_H0.

I can see advantages and disadvantages to both approaches. The first approach lets you use the existing microphysics infrastructure more or less as-is, since you can then just have things like free electrons be an "ion" that gets evolved like any other species -- just one that has Z = -1 and A = m_e / m_p. Everything just works out of the box. On the other hand, you will only obey your conservation laws to the accuracy of your ODE solver, and potentially significant violations of things like charge conservation build up over time, and you do a bunch of extra work evolving species that you really don't have to.

The second approach maintains conservation to machine precision and also saves you a bunch of work, but requires writing a bunch of custom code that enforces the conservation laws appropriate for each network. This is the approach I take in DESPOTIC: I make each network define both a list of species that it follows and an "extended" species list that contains both the species that the network explicitly evolves and the additional species that can be deduced from them by conservation laws. The network then also supplies a function to deduce the abundances of the extended species list from the primary list, so that you can have the abundances of all species when you need them.

So I think a fundamental design question we have to address, of which the EOS is just a special case, is whether we want to use strategy (1) or strategy (2).

[BenWibking]

I think strategy (2) has to be designed very carefully to avoid a catastrophic loss of numerical precision in the variables that are deduced from conservation laws. Since subtraction is necessary to compute, e.g., the electron number density, catastrophic loss of precision can happen in the cancellation. For instance, this is why I used the radiation-matter exchange method where both the radiation energy and matter energy are treated as independent variables.

Is this of immediate relevance for adding the additional radiation energy variables? I would think those would all need to be added regardless of whether all of the species are added as variables in the network.

[BenWibking]

My preferred strategy would be to follow what Microphysics does for mass fractions, as far as I understand: i) renormalize everything to satisfy the conservation law during integration sub-steps, and/or ii) then renormalize at the end of the integration.

The renormalization step is not trivial with multiple constraints. However, I think the best thing to do might be to formulate the constraints as a system of linear equations that we then solve for the least-squares solution. I think almost anything else will have (more) undesirable numerical error properties.

[BenWibking]

It turns out that this renomalization approach is in fact exactly how Mathematica solves ODE systems with invariants (or "index 0" DAEs). Their projection method (applied after each integration step) works even for nonlinear constraints: https://reference.wolfram.com/language/tutorial/NDSolveProjection.html

[markkrumholz]

I certainly agree that you can get precision loss when the quantities you're trying to back out are very small compared to the ones you're following, but in general for chemistry networks this is not a problem. As I said, DESPOTIC successfully uses this strategy.

I also agree that we need to add radiation quantities to the networks regardless. The specific place where this came up is in how you treat the modified version of the EOS routines once you have added radiation quantities, where you use different strategies depending on whether you want to follow my options (1) or (2).

[BenWibking]

I still prefer option (1), because i) it makes the code simpler, and ii) produces (very nearly) uniform relative errors on all of the variables.

[chongchonghe]

Aligned with option 2, I got an idea that solves the problem of radiation quantities and EOS. We could use state.xn for both mass fractions/number densities $X_i$ and radiation quantities $E_g$. Then, we would create a custom EOS for the photoionization network, where I define the number of ionizing species N_ion = 0, 1, or 3 as a compiler flag, and that determines the mean molecular weight in EOS. For N_ion = 0, there is no species and the EOS is the gamma law with mu=1. For N_ion = 1, it's the Hydrogen ionizing network, xn[0] is the H^+ mass fraction, and the EOS is a gamma law with mu^-1 = 1 + xn[0], C_v = 1.5 * mu^-1.

[chongchonghe]

If we want to keep track of all species (option 1), we should evolve the following species: H0, H+, He0, He+, He++, e-, and electron would have Z = 1 and A = m_e/m_u. Then,

$$ \frac{1}{\mu} = \sum_k \frac{Z_k X_k}{A_k}. $$

where $X_k$ is mass fraction. Then, we can use the available EOS in Microphysics. For this approach, we need a separate variable for radiation, say state.rad, which evolves like state.xn but doesn't go into EOS. Although, a trick can be done by setting $Z=0$ for radiation quantities and putting them into state.xn. I don't know if there is any potential issue with this approach, or if there is a better place to hold radiation quantities @zingale

[BenWibking]

I don't see the downside in creating a separate array for radiation quantities. It should be almost trivial to modify the integrator to loop over an array-like wrapper object whose indices enumerate all the components of both arrays. For all existing networks, it reduces to a zero-cost wrapper object for the existing state.xn[] array.

Also, I don't see any downside in option 1. I see several downsides in option 2.

[BenWibking]

@zingale Do you have any thoughts on Chongchong's two proposals above on where to put the E_g variables in Microphysics?

As I understand it, either:

• We add a state.rad[] array and perform the requisite modifications to the integrator to evolve both state.xn[] and state.rad[], or
• We add the radiation $E_g$ variables at the end of the state.xn[] array (in a somewhat hack-y way)

We do need to figure it out soon in order to continue making progress on the code.

[chongchonghe]

Aligned with option 2, I got an idea that solves the problem of radiation quantities and EOS. We could use state.xn for both mass fractions/number densities X i  and radiation quantities E g . Then, we would create a custom EOS for the photoionization network, where I define the number of ionizing species N_ion = 0, 1, or 3 as a compiler flag, and that determines the mean molecular weight in EOS. For N_ion = 0, there is no species and the EOS is the gamma law with mu=1. For N_ion = 1, it's the Hydrogen ionizing network, xn[0] is the H^+ mass fraction, and the EOS is a gamma law with mu^-1 = 1 + xn[0], C_v = 1.5 * mu^-1.

I hold back on this idea, because as we include more and more species in the radiaiton-chemistry network, it will become very dirty. Although, a hack to set $Z=0$ for radiation components and put them into state.xn might still work.

[zingale]

I think the option 1 method of explicitly tracking every species is most consistent with the current design, so that would be my preference.

For the radiation energy, you might be able to piggy back onto the SDC stuff, because there we pass the entire conserved state into the solver (via the burn_t) and, for instance, update density during the burn evolution. So you might be able to do something similar for radiation energy.

[BenWibking]

For the radiation energy, you might be able to piggy back onto the SDC stuff, because there we pass the entire conserved state into the solver (via the burn_t) and, for instance, update density during the burn evolution. So you might be able to do something similar for radiation energy.

That could be very useful. We currently have both radiation and hydro variables as part of the conserved state. In the future, this would also let us save the advective fluxes and do the hydro update using the Microphysics integrator itself, following https://doi.org/10.1006/jcph.1996.0165, without having to rewrite our hydro integrator to do SDC.