Dust temperature coupled to gas and radiation

[chongchonghe]

Description

Add support for dust temperature coupled to radiation and gas, following Bate & Keto (2015).

The dust is assumed to be in thermal equilibrium all the time:

$$ \Lambda_{{\rm gd},g} = - c \left[ \chi_{0,g} E_g - \chi_{0,g} \frac{4 \pi}{c} B_g(T_d) \right] $$

where $\sum_g \Lambda_{{\rm gd},g} = \Lambda_{{\rm gd}}$ is the dust-gas interaction rate:

$$ \Lambda_{\mathrm{gd}} = \Lambda_{{\mathrm{gd}}}^0 n_{\mathrm{H}}^2 T^{1 / 2}\left(T-T_{{d}}\right) $$

and $\Lambda_{\mathrm{gd}}^0 = 2.5 \times 10^{-34} \ {\rm erg \ cm^3 \ s^{-1} \ K^{-3/2}}$.

In the matter-radiation coupling step, we ignore the opacity of the gas and the coupling equations become:

$$ \begin{align} \frac{\partial E_{\rm gas}}{\partial t} &= - \sum_g \Lambda_{{\rm gd},g} \\ \frac{\partial E_g}{\partial t} &= \frac{\hat{c}}{c} \Lambda_{{\rm gd},g} \end{align} $$

At high density, the dust and gas are thermally coupled to each other, so $T_d = T_{\rm gas}$, and the original recipe is retrieved. In the low-density regime, however, the dust temperature is significantly lower than the gas temperature, therefore the thermal emission of the gas is very low and the gas relies on line emission to cool.

This PR only implements dust temperature and I will implement line cooling, photoelectric heating and cosmic ray heating in the next PR.

More changes:

• Added a new parameter enable_dust in RadSystem_Traits to control enabling dust.
• Added two tests, RadDust and RadDustMG, to test the new dust physics. The simulation results are compared to the reference results obtained from solving the above equations using an ODE integrator.

Related issues

None.

Checklist

Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an x inside the square brackets [ ] in the Markdown source below:

• I have added a description (see above).
• I have added a link to any related issues see (see above).
• I have read the Contributing Guide.
• I have added tests for any new physics that this PR adds to the code.
• I have tested this PR on my local computer and all tests pass.
• I have manually triggered the GPU tests with the magic comment /azp run.
• I have requested a reviewer for this PR.

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[markkrumholz]

Looks good to me other than the issue of using a more descriptive enable that I raised in the comments. One more marker to lay down, which doesn't have to be in this PR but we will probably want at some point, is to switch the dust-gas coupling coefficient from a constant to a function, the same way opacity is a function. This will enable use to do more sophisticated problems where the dust-gas coupling coefficient is not constant in time and space, which for example we might want in a simulation where we are following dust physics or metallicity evolution. Constant is fine for now, though.

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[chongchonghe]

Should be able to pass all tests. I will merge it when all tests pass, once you approve it.

[BenWibking]

We can't handle errors this way on the GPU. This needs to be fixed.

[chongchonghe]

OK. I will solve it along with #643

[BenWibking]

It may make sense to move the dust physics to a header file in a separate src/dust/ subdirectory.

[markkrumholz]

This isn't really dust physics per se, it's more a treatment of radiation in the case when you don't assume perfect thermodynamic equilibrium between the dust and gas. We're going to add more non-LTE processes soon, including photoelectric heating and line cooling, and I guess it might make sense to put all of those into a separate module for non-LTE ISM physics.

[chongchonghe]

Right, as Mark said, it's a simple, imperfect thermal interaction network between dust, radiation and gas, and the only change is a few terms in the Jacobian, so it makes sense to keep it inside the radiation system.

[BenWibking]

I agree with Mark. I think this makes sense to put in a separate microphysics source file for ISM physics, since it doesn't have to do with solving the radiation moment equations generically. It's fine for now. But it will get very cluttered if you add all of the possible non-LTE microphysics anyone would ever want to use with Quokka in radiation_system.hpp.

[markkrumholz]

I am sympathetic with your desire to keep the code clean, but we are going to need to think a bit about how exactly to do that. The issue is that all of these different non-LTE terms are going to modify the RHS and the Jacobian in the radiation sub-system, and it's not just a matter of adding new terms to the existing LTE forms. Thus the calculations that go into solving the RHD system are inevitably going to fundamentally depend on what microphysical model one has for the matter interacting with the radiation. So far we have been assuming that those dependences can be characterized just by specifying the opacities, which is true in LTE, but not in general.

Given this, here is one idea for how we might organize this: we pull the calculation of the Jacobian out as a function, which lives in a separate matter-radiation interaction file, something like radiation_microphysics.hpp. This file is responsible for supplying two functions needed by the radiation implicit solver: the Jacobian and the convergence test (which is equivalent to the RHS of the system). We can then put different implementations of that microphysics, corresponding to different choices of Jacobian and RHS, into that file, separated by constexpr's. In this view our current default implementation, where we assume that the matter is in LTE and its properties can be described just opacity functions, is one choice of radiation microphysics. The model @chongchonghe just implemented, where dust is assumed to be in LTE at a single temperature but that temperature is not the same as the gas temperature, but there are no non-LTE emission processes, is a second microphysics model. Adding optically thin line emission processes represents yet a third model. Depending on the level of complexity these different models could all be put in a single file, or could go in separate files using template overrides to let users select the model and compile time.

My inspiration for this is how we're handling chemistry, where we segregate the parts that are in common -- have an ODE solver that you use to integrate the chemical network forward -- from the parts that are unique to each network, and which are selected by a physics trait. Same idea here: we keep the parts that are in common to all microphysics modules -- the Netwon-Raphson iteration -- in one place, and then put implementations of different micorphysics modules separately.

Thoughts on this idea?

[BenWibking]

I think something like that makes sense. I agree that having it choose the Jacobian and rhs at compile time is the right thing to do. I had in mind different Jacobians and rhs that could be selected via templates, but I'm not sure of the exact implementation mechanics.

Is there a reason to write a new convergence test for each non-LTE model?

[BenWibking]

Actually, I can definitely think of cases where you might want to do that.

However, it's very easy to screw up the implementation of the convergence test (I did this many times in the early versions of the code), so I think it might be better to impose more structure on this. In principle, it should be possible to programmatically specify which variables that you want to check for convergence without re-writing the convergence test itself. This might require some macro or template metaprogramming tricks.

[markkrumholz]

Agree that automating convergence testing, with the possibility of overriding, is probably the right thing to do. As for how to handle this programmatically, I’d suggest that we have a physics trait with an enum for radiation-matter coupling model, the same way we do now for chemistry networks. The original quokka implementation can be called LTE, for matter in LTE, and the one @chongchonghe just implemented can be called dust_gas_coupling. Once we add line cooling and FUV heating we then have a complete physics model for the atomic ISM, and we can call it atomic_ISM. Each of these would live in a different file, with template overrides to select the version to use at compile time.

[BenWibking]

Sounds good to me 👍