A robust numerical scheme for the dust-gas coupling model

[chongchonghe]

Description

In this PR, we solve the dust coupling to the gas and radiation in two different regimes: the well-coupled regime and the low-density, poorly-coupled regime.

Changes to the code

• Defined ComputeJacobianForGasAndDustDecoupled to compute the Jacobian for the gas-dust-radiation system in the regime where gas and dust are decoupled.
• Added a new parameter lambda_gd_times_dt to SolveMatterRadiationEnergyExchange. lambda_gd_times_dt is $\Lambda_{\rm gd} \Delta t$, which is used to update gas energy as well as constraining the sum of $R_g$.
• Added a new parameter dust_model to SolveMatterRadiationEnergyExchange. Unfortunately, the energy-exchange system becomes complicated and you can't deal with all three dust models by passing functions to SolveMatterRadiationEnergyExchange. dust_model takes three values, representing a gas-dust-radiation system in different regimes:
  • 0: dust and gas are perfectly coupled. In this case, T_d = T_gas, and we solve a system that involves T_gas and Erad. This is the model used in Wibking & Krumholz (2022), He, Wibking & Krumholz (2024a) and He, Wibking & Krumholz (2024b).
  • 1: dust and gas are well coupled. In this case, we solve a system that involves T_gas, T_d, and Erad, where T_d is used as an implicit, dependent variable.
  • 2: dust and gas are poorly coupled. In this case, we solve a system that involves only T_d and Erad.
• dust_model = 0 is enabled by setting enable_dust_gas_thermal_coupling_model_ = false. Otherwise, dust_model = 1 or dust_model = 2 is determined automatically in the code on a cell-by-cell basis. See the description below.

Alternative things to do

• Make threshold_for_dust_gas_coupling a runtime parameter. This is the parameter to control when to switch between dust model 1 and 2.

Physics

In the well-coupled regime, we solve the following equations:

$$ \begin{split} \frac{\partial E_{\rm gas}}{\partial t} = - \Lambda_{{\rm gd}} \\ \frac{\partial E_g}{\partial t} = \hat{c} \left[ \chi_B \frac{4 \pi}{c} B_g(T_d) - \chi_E E_g \right] \end{split} $$

where $\Lambda_{\mathrm{gd}}$ is the thermal interaction rate due to collisions between the gas and the dust, which we take the rate of Keto & Field (2005)

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

where $\Theta_{\mathrm{gd}} = 2.5 \times 10^{-34} \mathrm{erg \ cm^3 \ s^{-1} \ K^{-3/2}}$. This parameter can be specified in arbitrary units by the user. As in most literature, we assume that the dust is in LTE with the total radiation field (e.g. Bate & Keto 2015):

$$ c \sum_g \left[ \chi_B \frac{4 \pi}{c} B_g(T_d) - \chi_E E_g \right] = \Lambda_{{\rm gd}}. \quad{} (2) $$

We can then define a similar set of independent variables as in the IMEX paper, $(E_{\rm gas}, R_g), g = 1, 2, ..., N_g$, where

$$ R_g \equiv \hat{c} \Delta t \left[ \chi_B \frac{4 \pi}{c} B_g(T_d) - \chi_E E_g \right], \quad{} (3) $$

and solve for the following set of equations using Newton-Raphson iteration:

$$ \begin{split} 0 = F_G \equiv E_{\rm gas} - E_{\rm gas}^0 + \frac{c}{\hat{c}} \sum_g R_g,\\ 0 = F_{R,g} \equiv E_g - E_g^0 - R_g \end{split} $$

This is similar to the original gas-only scheme, except that the gas temperature $T$ is replaced with $T_d$.

We still need to calculate the dust temperature $T_d$. In the first step of the iteration, we solve for $T_d$ given $T$ and $E_g$ via combining Equation (1) and (2). Then, later in the iteration, we solve for $T_d$ via a simple relation:

$$ \sum_g R_g = \frac{\hat{c}}{c} \Delta t \ \Lambda_{\rm gd} = \frac{\hat{c}}{c} \Delta t \ \Theta_{{\rm gd}} ~n_{\mathrm{H}}^2 T^{1 / 2}\left(T-T_{{d}}\right) = N_{\rm gd} T^{1/2}(T - T_d) \quad{} (4) $$

where we have defined $N_{\rm gd} \equiv (\hat{c}/c) \Delta t \Theta_{\rm gd} n_H^2$. Apparently,

$$ T_d = T - \frac{\sum_g R_g}{N_{\rm gd} \sqrt{T}}. \quad{} (5) $$

Equation (4) also gives a simple formula for $d T_d / d T$, which is used in the computation of the Jacobian. Further, the Jacobian is no longer sparse. However, I did a trick to analytically simplify the Jacobian matrix via Gauss-Jordan elimination to make it sparse so that we can use the numerical scheme in the multigroup paper. More details in Overleaf document "quokka-thermochemistry".

In the decoupled regime where $N_{\rm gd}$ is small, however, this scheme failed. This is obvious from Equation (5): even a small value of $\sum_g R_g$ will result in a large, negative $T_d$. This is unavoidable in the Newton-Raphson step; the valid domain for $\sum_g R_g$ is tiny: $[0, N_{\rm gd} T^{3/2})$. Here I'm considering the case where the radiation temperature is well below the gas temperature, so $T \gg T_r \approx T_d$, therefore $\Lambda_{\rm gd} \ge 0$.

After extensive exploration, I constructed the following scheme for the decoupled regime. First, we define the decoupled regime as when $\max(N_{\rm gd} T^{3/2}, N_{\rm gd} T_d^{3/2}) < 10^{-6} \max(E_{\rm gas}, \sum_g R_g)$. Since the gas temperature won't change much, it can be updated via a simple forward Euler step of the following equation:

$$ C_V \frac{d T}{dt} = - \Lambda_{\rm gd} = - \Theta_{{\rm gd}} ~n_{\mathrm{H}}^2 \ T^{1 / 2}\left(T-T_{{d}}\right) \quad{} (6) $$

One can also choose to solve this analytically. Thus we have obtained $\Lambda_{\rm gd} \Delta t$.

Then, we want to deal with energy exchange between radiation groups, to capture dust absorbing UV radiation and re-emitting in IR. At first, I tried to update $E_g$ with fixed $T_d$, by solving

$$ E_g - E_g^0 = c \Delta t \ \chi_g (B_g(T_d) - E_g). $$

This can be solved with simple algebra. However, I have shown that the resulting $\sum_g R_g$ (the sum of the RHS over groups) does not necessarily equal to the $\Lambda_{\rm gd} \Delta t$ used previously to update gas temperature. This results in the failure of energy conservation. To enforce $\sum_g R_g = \Lambda_{\rm gd} \Delta t$, we have to include $T_d$ as an independent variable; otherwise, the system is over-determined.

Then, here is the new scheme that I have shown to be very robust. After updating $T$ via Equation (6), we solve the updated $T_d$ and $E_g$ via the following set of equations:

$$ \begin{split} \sum_g R_g = \Lambda_{\rm gd} \Delta t, \\ E_g - E_g^0 = R_g. \end{split} $$

We solve this equation via Newton-Raphson iteration on the base variable $(T_d, R_g)$.

$$ \begin{split} 0 = F_0 = \sum_g R_g - \Lambda_{\rm gd} \Delta t \\ 0 = F_g = E_g - E_g^0 - R_g \end{split} $$

and the Jacobian is

$$ \begin{split} J_{00} = 0 \\ J_{0g} = 1 \\ J_{g0} = \frac{\partial E_g}{\partial T_d} = \frac{\chi_{B,g}}{\chi_{E,g}} \frac{\partial B_g(T_d)}{\partial T_d} \\ J_{gg} = \frac{\partial E_g}{\partial R_g} - 1 = - \frac{\chi_{B,g}}{\chi_{E,g}} \frac{1}{\hat{c} \chi_{\rm B,g} \Delta t} - 1 \end{split} $$

This set of equations is very similar to the previous ones, so we can use the same routine with some switches to reset different Jacobian for different regimes. After this system is converged, the whole iteration is done and energy is guaranteed to be conserved (to the degree of how well the Newton-Raphson iteration converges).

Tests

RadMarshakDust is a Marshak wave problem where FUV radiation is streaming from the left boundary, gets absorbed by dust which then radiates in IR. The IR opacity is set to a large number so that diffusion is extremely slow. I compare the (semi-)exact analytic solution with the numerical calculation and find perfect agreement. Figure shown below:

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]

@markkrumholz You can start to review this PR. I'll try to merge this in once #731 is merged.

[BenWibking]

@chongchonghe @markkrumholz Before adding an very elaborate and baroque scheme that uses different numerical methods in different stiffness regimes, I think it would be worth investigating adding backtracking to the Newton-Raphson solver would help: #11

[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).

[BenWibking]

The code looks fine.

Two questions:

1. How did you determine 1e-6 as the threshold value to switch between the two dust solver? Is there a way to set this threshold based on an error tolerance for the dust temperature directly?
2. In your description of the test problem, you say this: "The IR opacity is set to a large number so that there is zero diffusion." I don't understand how there is zero diffusion with a large opacity. Can you explain?

[chongchonghe]

The code looks fine.

Two questions:

1. How did you determine 1e-6 as the threshold value to switch between the two dust solver? Is there a way to set this threshold based on an error tolerance for the dust temperature directly?

This threshold is determined empirically through testing the algorithm for a range of gas density and temperature. If we I set a lower number, at moderate density the algorithm falls into the well-coupled regime and fails to converge. We also want to set this number as low as possible to avoid truncations error due to the choice of the decoupled assumption. I agree that we need to justify the choice of this number by estimating the relative error. Because this transition happens in the intermediate regime, it's hard to design a test with an exact solution. I plan to do a simple analytic analysis.

Bate & Keto (2015) did a similar trick by dealing the coupled and poorly coupled regime separately. They described their choice of transition condition (above Eq. 57) without justification. However, their equations is much simpler because they use single-group radiation so they don't have to worry about redistribution of radiation in frequency space by dust. In single-group RHD, our "well-coupled regime" scheme works for all regimes.

2. In your description of the test problem, you say this: "The IR opacity is set to a large number so that there is zero diffusion." I don't understand how there is zero diffusion with a large opacity. Can you explain?

The diffusion coefficient of material with opacity $\chi$ is $c/(3\chi)$, which comes from the diffusion equation $\partial E / \partial t = c / (3 \chi) \nabla^2 E$. So, with a large opacity, the amount of diffusion is zero. In this limit, the amount of IR radiation at any cell is easily calculated via integration of IR emission rate at that cell over time.

[BenWibking]

This threshold is determined empirically through testing the algorithm for a range of gas density and temperature. If we I set a lower number, at moderate density the algorithm falls into the well-coupled regime and fails to converge. We also want to set this number as low as possible to avoid truncations error due to the choice of the decoupled assumption. I agree that we need to justify the choice of this number by estimating the relative error. Because this transition happens in the intermediate regime, it's hard to design a test with an exact solution. I plan to do a simple analytic analysis.

Ok, that would be good. It should be possible to estimate from the next-to-leading order term.

The diffusion coefficient of material with opacity χ is c / ( 3 χ ) , which comes from the diffusion equation ∂ E / ∂ t = c / ( 3 χ ) ∇ 2 E . So, with a large opacity, the amount of diffusion is zero. In this limit, the amount of IR radiation at any cell is easily calculated via integration of IR emission rate at that cell over time.

I would call this a regime where the diffusion is arbitrarily slow, rather than zero. To me, "zero diffusion" makes it sound like it is streaming, rather than just very slowly diffusion.

[chongchonghe]

Here's an estimation of the 'error' of the decoupled approximation. It's not an error per se. Every numerical scheme to solve PDE has an error. All I did here is to use forward Euler as opposed to backward Euler for gas temperature update, and the relative change in the temperature is at most 10^-6, determined by that threshold parameter, so it's very stable. Here's why.

In the decoupled regime, we use the following equation to update the gas temperature via the Forward Euler scheme:

$$ C_V d T / d t = \Lambda_{\rm gd}^0 \sqrt{T} (T - T_d) $$

and

$$ C_V (T - T^0) = \Lambda_{\rm gd}^0 \Delta t \sqrt{T} (T - T_d) < \Lambda_{\rm gd}^0 \Delta t T^{3/2} < \epsilon C_V T^0 $$

The parameter $\epsilon = 10^{-6}$ was described in my previous writing, right before Equation (6). Note that $N_{\rm gd} \equiv \Lambda_{\rm gd}^0 \Delta t$. From here, we can estimate the relative error of this scheme on the gas temperature update, compared to the "exact solution". Say, we start from $T = 1$, and we take the maximum $N_{\rm gd} = 10^{-6}$, then the updated $T = 1 + 10^{-6}$. If we were to use an implicit update, the updated temperature would be about $T = 1 + 10^{-6}(1 + 10^{-6}) = 1 + 10^{-6} + 10^{-12}$. Therefore, although the error relative to the change of the temperature can be up to $\epsilon = 10^{-6}$, the error relative to the value of the temperature is at most $\epsilon^2 = 10^{-12}$. That should be totally tolerable for a numerical scheme. After all, although implicit ODE solver is stable, we can't say much about it's accuracy.

[chongchonghe]

The relative error on gas temperature introduced by the photoionization scheme from the RAMSES paper will be 5 orders of magnitude bigger than this (their ten-percent rule), and we will probably have to do something similar. So, be prepared for that.

[sonarqubecloud]

Quality Gate passed

Issues
11 New issues
0 Accepted issues

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

See analysis details on SonarCloud

[BenWibking]

The relative error on gas temperature introduced by the photoionization scheme from the RAMSES paper will be 5 orders of magnitude bigger than this (their ten-percent rule), and we will probably have to do something similar. So, be prepared for that.

If we need to do adaptive error control, I think the radiation-matter exchange solver should just be moved into Microphysics. I will create a Slack channel to discuss this.

[BenWibking]

The parameter ϵ = 10 − 6 was described in my previous writing, right before Equation (6). Note that N g d ≡ Λ g d 0 Δ t . From here, we can estimate the relative error of this scheme on the gas temperature update, compared to the "exact solution". Say, we start from T = 1 , and we take the maximum N g d = 10 − 6 , then the updated T = 1 + 10 − 6 . If we were to use an implicit update, the updated temperature would be about T = 1 + 10 − 6 ( 1 + 10 − 6 ) = 1 + 10 − 6 + 10 − 12 . Therefore, although the error relative to the change of the temperature can be up to ϵ = 10 − 6 , the error relative to the value of the temperature is at most ϵ 2 = 10 − 12 . That should be totally tolerable for a numerical scheme. After all, although implicit ODE solver is stable, we can't say much about it's accuracy.

I think the right thing to do here is adaptively control the relative error based on explicit computation of the of the higher-order terms in an expansion of the solution in the dust-gas coupling constant. But that can be done in a future PR.

[BenWibking]

Ok, I've edited the PR description to say that diffusion is "extremely slow" instead of "zero". I will approve the PR, assuming tests pass.

[BenWibking]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).