Line cooling and cosmic ray heating #773

chongchonghe · 2024-10-11T05:44:44Z

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

There are four radiation bins. The first bin is the line emission and the last bin is photoelectric heating (FUV).
The initial radiation energy density is 1 in the last bin and 0 in all other bins. The opacity is 0 in all bins.
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.

for more information, see https://pre-commit.ci

src/problems/RadLineCooling/test_rad_line_cooling.cpp

src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp

src/radiation/radiation_system.hpp

chongchonghe · 2024-10-11T08:29:02Z

/azp run

azure-pipelines · 2024-10-11T08:29:16Z

Azure Pipelines successfully started running 2 pipeline(s).

for more information, see https://pre-commit.ci

chongchonghe · 2024-10-11T10:04:41Z

/azp run

azure-pipelines · 2024-10-11T10:04:53Z

Azure Pipelines successfully started running 2 pipeline(s).

sonarcloud · 2024-10-11T10:08:03Z

Quality Gate passed

Issues
21 New issues
0 Accepted issues

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

See analysis details on SonarCloud

src/problems/RadLineCooling/test_rad_line_cooling.cpp

markkrumholz · 2024-10-20T22:49:20Z

src/problems/RadLineCooling/test_rad_line_cooling.cpp

+	sim.evolve();
+
+	// read output variables
+	auto [_, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0);


Please either fix or NOLINT this.

src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp

markkrumholz · 2024-10-20T22:51:39Z

src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp

+	const bool is_coupled = sim.dustGasInteractionCoeff_ > 1.0;
+
+	// read output variables
+	auto [_, values] = fextract(sim.state_new_cc_[0], sim.Geom(0), 0, 0.0);


Again, please fix or NOLINT

chongchonghe · 2024-10-14T07:58:05Z

@markkrumholz This PR is ready for your review.

markkrumholz · 2024-10-14T21:54:36Z

@markkrumholz This PR is ready for your review.

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

markkrumholz · 2024-10-20T22:46:36Z

src/math/interpolate.cpp

@@ -108,6 +108,9 @@ void interpolate_arrays(double *x, double *y, int len, double *arr_x, double *ar
 	/* Note: arr_x must be sorted in ascending order,
 		and arr_len must be >= 3. */

+	// check if x is within the range of arr_x


You should change this to one of the AMReX assert macros so that it functions properly on GPU or GPU, and in parallel. See https://amrex-codes.github.io/amrex/doxygen/AMReX__BLassert_8H.html

markkrumholz · 2024-10-20T22:48:29Z

src/problems/RadLineCooling/test_rad_line_cooling.cpp

+
+const double max_time = 10.0;
+
+template <> struct SimulationData<PulseProblem> {


Maybe call this something other than PulseProblem, since it appears to be unrelated to the actual radiation pulse problem described in our papers?

markkrumholz · 2024-10-20T22:50:23Z

src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp

+const double CR_heating_rate = 0.03;
+const double PE_rate = 0.02;
+
+template <> struct SimulationData<PulseProblem> {


Same comment about the name of this test as with the non-multigroup version -- call it something more usefully descriptive.

markkrumholz · 2024-10-20T22:51:26Z

src/problems/RadLineCoolingMG/test_rad_line_cooling_MG.cpp

+
+auto problem_main() -> int
+{
+	// This problem is a *linear* radiation diffusion problem, i.e.


I think this comment does not apply to this problem, and is left over from the problem you cloned. Please correct or delete.

markkrumholz · 2024-10-20T22:55:27Z

src/radiation/radiation_system.hpp

@@ -138,6 +139,10 @@ template <typename problem_t> struct FluxUpdateResult {
 	amrex::GpuArray<amrex::GpuArray<amrex::Real, Physics_Traits<problem_t>::nGroups>, 3> Frad; // radiation flux
 };

+// template <typename problem_t> struct SingleVariableState {


Do we need to keep this commented-out code, or can it be deleted?

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 added 30 commits October 6, 2024 15:46

remove unused headers in radhydro_bb.cpp

0612941

defined the rad_line_cooling test

2a2005a

declare functions

2483092

Merge branch 'development' into chong/line_cooling_and_CR_heating

2689a53

mod result.F0 calculation

7048b3d

update line_cooling.cpp

c693bce

tmp

a5eb36a

update DefineNetCoolingRate

d97d66e

cleanup

953364d

add empty functions to ComputeJacobianForGasAndDust

785f0ef

finish ComputeJacobianForGasAndDust; passing tests

7d0edfd

add line cooling to single group radiation; passing tests

c81f49f

PE line cooling: step 1

32ca0d0

PE line cooling: step 2, mod ComputeJacobianForGasAndDustWithPE

6008c1c

step3: add cooling to ComputeJacobianForGasAndDustWithPE

45687b0

fix missing dt problem; passing tests

39f6534

fix Jgg

e5ecb9f

update cooling_MG.cpp to include both cooling and PE; passing

a368347

put dt inside cooling

944a114

Cooling in gas_dust_decoupled; passing test

26aef50

Cooling in gad_dust_PE_decoupled; passing test

e29431f

line cooling to radiation for both decoupled case, working

4e7afa1

cleanup

84a4477

pass num_den to DefineNetCoolingRate

2fb9884

use IntegratorOneVariable to solve for T_d; passing tests on CPU and GPU

f61a9b2

finish ComputeDustTemperatureBateKeto; passing tests

591222a

fix spelling

1d4c617

implicit Egas update in decoupled GasDustRad system

fcba84c

implicit Egas update in decoupled GasDustRadPE system

be3c703

update cooling_MG test to make it plot T(t)

0fc7dfe

chongchonghe added 6 commits October 11, 2024 00:52

add RadLineCoolingMG_coupled test; passing

a5c8451

add CR heating to single group

2b88743

rename cooling_t0 to cooling_tend

7120177

add line cooling to radiation group in single_group

221f193

fix missing num_den

de072df

cleanup

a8df3b9

chongchonghe requested review from BenWibking and markkrumholz as code owners October 11, 2024 05:44

dosubot bot added the size:XXL This PR changes 1000+ lines, ignoring generated files. label Oct 11, 2024

pre-commit-ci bot and others added 4 commits October 11, 2024 05:44

[pre-commit.ci] auto fixes from pre-commit.com hooks

4c5af29

for more information, see https://pre-commit.ci

clean unused variable

93c72fe

fix error "may be used uninitialized "

5167920

remove unused variables

cdf0c7e

github-advanced-security bot found potential problems Oct 11, 2024

View reviewed changes

chongchonghe and others added 2 commits October 11, 2024 21:03

fix github-advanced-security warnings

bda155c

[pre-commit.ci] auto fixes from pre-commit.com hooks

3b58d97

for more information, see https://pre-commit.ci

github-advanced-security bot found potential problems Oct 11, 2024

View reviewed changes

markkrumholz reviewed Oct 20, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Line cooling and cosmic ray heating #773

Line cooling and cosmic ray heating #773

chongchonghe commented Oct 11, 2024 •

edited

Loading

chongchonghe commented Oct 11, 2024

azure-pipelines bot commented Oct 11, 2024

chongchonghe commented Oct 11, 2024

azure-pipelines bot commented Oct 11, 2024

sonarcloud bot commented Oct 11, 2024

markkrumholz Oct 20, 2024

markkrumholz Oct 20, 2024

chongchonghe commented Oct 14, 2024

markkrumholz commented Oct 14, 2024

markkrumholz Oct 20, 2024

markkrumholz Oct 20, 2024

markkrumholz Oct 20, 2024

markkrumholz Oct 20, 2024

markkrumholz Oct 20, 2024

markkrumholz left a comment


		const double max_time = 10.0;

		template <> struct SimulationData<PulseProblem> {

Line cooling and cosmic ray heating #773

Are you sure you want to change the base?

Line cooling and cosmic ray heating #773

Conversation

chongchonghe commented Oct 11, 2024 • edited Loading

Description

Test problems

RadLineCooling

RadLineCoolingMG

Related issues

Checklist

chongchonghe commented Oct 11, 2024

azure-pipelines bot commented Oct 11, 2024

chongchonghe commented Oct 11, 2024

azure-pipelines bot commented Oct 11, 2024

sonarcloud bot commented Oct 11, 2024

Quality Gate passed

markkrumholz Oct 20, 2024

Choose a reason for hiding this comment

markkrumholz Oct 20, 2024

Choose a reason for hiding this comment

chongchonghe commented Oct 14, 2024

markkrumholz commented Oct 14, 2024

markkrumholz Oct 20, 2024

Choose a reason for hiding this comment

markkrumholz Oct 20, 2024

Choose a reason for hiding this comment

markkrumholz Oct 20, 2024

Choose a reason for hiding this comment

markkrumholz Oct 20, 2024

Choose a reason for hiding this comment

markkrumholz Oct 20, 2024

Choose a reason for hiding this comment

markkrumholz left a comment

Choose a reason for hiding this comment

chongchonghe commented Oct 11, 2024 •

edited

Loading