Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
A robust numerical scheme for the dust-gas coupling model (#733)
### 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:_
- [x] I have added a description (see above).
- [x] I have added a link to any related issues see (see above).
- [x] I have read the [Contributing
Guide](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
- [x] I have added tests for any new physics that this PR adds to the
code.
- [x] I have tested this PR on my local computer and all tests pass.
- [x] I have manually triggered the GPU tests with the magic comment
`/azp run`.
- [x] I have requested a reviewer for this PR.
---------
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>- Loading branch information
1 parent
4a7057f
commit 2d273a2