Line cooling and cosmic ray heating

[chongchonghe]

Description

This pull request implements line cooling and cosmic ray heating. The key changes include:

• Added line cooling to the single-group and multi-group radiation solvers when enable_dust_gas_thermal_coupling_model = true . Line cooling is enabled by defining DefineNetCoolingRate and DefineNetCoolingRateTempDerivative in the problem file.
• Implemented cosmic ray heating in all solvers, which is enabled by defining DefineCosmicRayHeatingRate in the problem file.
• Added tests for the new line cooling and cosmic ray heating functionality.
  • RadLineCooling: test line cooling + CR heating + PE heating for single-group radiation system.
  • RadLineCoolingMG: test line cooling + CR heating + PE heating for multi-group radiation system, for both well-coupled and decoupled gas-dust systems.

This is the last feature I have finished implementing a while ago and I want to merge into the development branch before we move to Microphysics. In the next PR, I'll try to structure the code in the format of actual_rhs(), actual_jac(), andbackward_euler_integrator() so that we can readily start to use Microphysics API.

Test problems

RadLineCooling

Initial conditions: uniform static gas with a temperature $T_0 = 1$ and $C_V = 1$; the radiation energy density is 0. The dust opacity is 0. The line cooling rate is $\Lambda_{\rm line} = C_{\Lambda} T$, where $C_{\Lambda} = 0.1$, and the CR heating rate is $\Gamma_{\rm cr} = 0.03$. As the gas cools down, its temperature follows the following ODE:

$$ C_V \frac{d T}{d t} = - C_{\Lambda} T + \Gamma_{\rm cr} $$

The solution is

$$\begin{gather} T(t) = C_{\Lambda}^{-1} e^{-\frac{C_{\Lambda }}{C_V} t} \left(\Gamma _{\text{cr}} e^{\frac{C_{\Lambda}}{C_V} t}+ T_0 C_{\Lambda }-\Gamma_{\text{cr}}\right) \\\ E_{\rm rad}(t) = -(C_V T(t) - C_V T_0 - \Gamma_{\rm cr} t) \end{gather}$$

In the end of the test, the numerical results of $T(t)$ and $E_{\rm rad}(t)$ are checked against the exact analytical solution.

RadLineCoolingMG

The initial condition is the same as that of RadLineCooling, except

1. There are four radiation bins. The first bin is the line emission and the last bin is photoelectric heating (FUV).
2. The initial radiation energy density is 1 in the last bin and 0 in all other bins. The opacity is 0 in all bins.
3. PE heating is enabled and $\Gamma_{\rm pe} = 0.02 E_{\rm FUV} = 0.02$.
   The ODE is

$$ C_V \frac{d T}{d t} = - C_{\Lambda} T + \Gamma_{\rm cr} + \Gamma_{\rm pe} $$

and the solution is the same as the previous one with two differences: (1) $\Gamma_{\rm cr}$ is replaced with $\Gamma_{\rm cr} + \Gamma_{\rm pe}$; (2) $E_{\rm rad}$ represents the energy density of the first bin.

This test is run twice, one with well-coupled gas-dust model and one with decoupled gas-dust model.

Related issues

Are there any GitHub issues that are fixed by this pull request? Add a link to them here.

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.

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[chongchonghe]

@markkrumholz This PR is ready for your review.

[markkrumholz]

@markkrumholz This PR is ready for your review.

Will get to it ASAP, but it may be a few days.

[markkrumholz]

I have left some comments on individual bits of code that you can look at, but there is also an overall design issue, which is the way you're using number density. The problem is that you don't specify the number density of what -- number of H nuclei, number of atomic H, number of H2, etc. -- and so this code becomes ambiguous and potentially incorrect as soon as we have chemistry. I think any time we need a number density, we need to be extremely careful to specify number density of which species, and we need to encapsulate the process of deriving the number density of that species from the density and mass scalars inside a function that can be swapped in or out depending on the chemistry network. Thus for example the current code that does things like num_density = rho / mean_molecular_weight could be encapsulated in a function called derive_num_density_H0 which derives the number density of neutral hydrogen from the density and mass scalar array. If the chemistry trait you have then selects the network to be none, then the current implementation is turned on, but if another chemistry network is selected then this handled differently, depending on the species in that network.

What we don't want to do, though, is hardwire in a solution that works only when chemistry is off, and then have it break in subtle and hard-to-detect ways if anyone tries to use it with chemistry turned on.

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[chongchonghe]

I defined a function, ComputeNumberDensityAtomicH(double rho, amrex::GpuArray<Real, nmscalars_> const &massScalars), to calculate atomic hydrogen number density. I use atomic hydrogen instead of neutral hydrogen because that's what's used in Bate & Keto. The cosmic ray heating might relate to neutral hydrogen density, but line cooling and photoelectric heating is related to atomic hydrogen (plus metal abundance). So, I'll use 'atomic hydrogen' number density to avoid confusion.

[markkrumholz]

I defined a function, ComputeNumberDensityAtomicH(double rho, amrex::GpuArray<Real, nmscalars_> const &massScalars), to calculate atomic hydrogen number density. I use atomic hydrogen instead of neutral hydrogen because that's what's used in Bate & Keto. The cosmic ray heating might relate to neutral hydrogen density, but line cooling and photoelectric heating is related to atomic hydrogen (plus metal abundance). So, I'll use 'atomic hydrogen' number density to avoid confusion.

I'm a bit confused about the distinction you're making between "atomic hydrogen" and "neutral hydrogen", since usually they are used as synonyms. In general in the astrochemistry literature the notation is that n_H refers to the number density of H nuclei regardless of their chemical state, while n_H+, n_H0, and n_H2 refer to the number density of ionised hydrogen, neutral atomic hydrogen, and hydrogen molecules, so that we have the relationship n_H = n_H+ + n_H0 + 2 n_H2. I think it would be good to follow this notational convention in the code, so that we could have routines like ComputeNumberDensityH0, ComputeNumberDensityH, ComputeNumberDensityHp (for H+), etc. This will make the code more readable and less susceptible to mistakes.

[chongchonghe]

I defined a function, ComputeNumberDensityAtomicH(double rho, amrex::GpuArray<Real, nmscalars_> const &massScalars), to calculate atomic hydrogen number density. I use atomic hydrogen instead of neutral hydrogen because that's what's used in Bate & Keto. The cosmic ray heating might relate to neutral hydrogen density, but line cooling and photoelectric heating is related to atomic hydrogen (plus metal abundance). So, I'll use 'atomic hydrogen' number density to avoid confusion.

I'm a bit confused about the distinction you're making between "atomic hydrogen" and "neutral hydrogen", since usually they are used as synonyms. In general in the astrochemistry literature the notation is that n_H refers to the number density of H nuclei regardless of their chemical state, while n_H+, n_H0, and n_H2 refer to the number density of ionised hydrogen, neutral atomic hydrogen, and hydrogen molecules, so that we have the relationship n_H = n_H+ + n_H0 + 2 n_H2. I think it would be good to follow this notational convention in the code, so that we could have routines like ComputeNumberDensityH0, ComputeNumberDensityH, ComputeNumberDensityHp (for H+), etc. This will make the code more readable and less susceptible to mistakes.

OK, I see. Then, what I meant is density of H nuclei, n_H. I've changed it to ComputeNumberDensityH.

[sonarqubecloud]

Quality Gate passed

Issues
20 New issues
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
4.9% Duplication on New Code

See analysis details on SonarCloud

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).