improve convergence rate of the radiation Newton-Raphson solver

[chongchonghe]

Description

Changes:

1. Added a trick to make the Newton-Raphson iteration converge faster in some situations: when $\Delta T > \max(T_{\rm gas}, T_{\rm rad})$, set $T_{\rm gas}$ to $T_{\rm rad}$ and $R$ to 0. The motivation for doing this trick is explained below.
2. Added radiation.print_iteration_counts parameter to enable counting the mean and max number of Newton-Raphson iterations in each radiation step. When set to true, the log file would look like:

Coarse STEP 1 at t = 0 (0%) starts ...
time_subcycle = 0, total number of Newton-Raphson solvings is 128, (mean, max) number of Newton-Raphson iterations are 1, 1

Coarse STEP 2 at t = 0.0046875 (0.9375%) starts ...
time_subcycle = 0.0046875, total number of Newton-Raphson solvings is 128, (mean, max) number of Newton-Raphson iterations are 1, 1

3. Defined a new radiation problem, RadMarshakDust. Initially, the gas is set with a density and temperature of 1, and the initial radiation energy density is 0. Radiation with a temperature $T_r = 0.01$ streams in from the left boundary. Before the dust model was introduced, we assumed gas and dust were perfectly coupled, and the gas would be heated by radiation where it could penetrate, and in regions without radiation, the gas temperature would cool to 0.003159 (solution of $1 = a_R T^4 + T$). In this test, however, we turned on the dust model and set the dust-gas interaction coefficient to a very small number. As a result, the dust and gas are decoupled: the dust temperature will reach the radiation temperature and stay in equilibrium, while the gas remains at its initial temperature.

In a future PR, I will extend this test by incorporating a multigroup model. This will demonstrate how FUV radiation from the left boundary is absorbed by the dust, which then re-radiates in IR.

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](https://github.com/quokka-astro/quokka/blob/development/CONTRIBUTING.md).
• 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.

[BenWibking]

We should not do this. The auxiliary internal energy should never be used as a primary source of truth for the internal energy. It should only be used when sync'ing with the total energy in the hydro update.

[BenWibking]

The dominant source of error in the internal energy is due to truncation errors, not round off errors. These are not the same, and the distinction is very important.

[BenWibking]

A better way to solve this would be to add a call to SyncInternalEnergy after each timestep for radiation-only test problems. In real problems, this is done already during the hydro update.

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[chongchonghe]

The dominant source of error in the internal energy is due to truncation errors, not round off errors. These are not the same, and the distinction is very important.

I believe it is roundoff error in this case. What happens in this test is that the gas has an equilibrium temperature of 1 because I set the radiation temperature to 1. However, because the specific heat of radiation is very large, due to conservation of momentum, the gas momentum increases to 10^13 after the first step. Then, in the second time step, because the internal energy is derived from subtracting the total energy from the kinetic energy, the result is 0 due to roundoff error.

[chongchonghe]

We should not do this. The auxiliary internal energy should never be used as a primary source of truth for the internal energy. It should only be used when sync'ing with the total energy in the hydro update.

Why shouldn't we use the auxiliary internal energy?

[chongchonghe]

We should not do this. The auxiliary internal energy should never be used as a primary source of truth for the internal energy. It should only be used when sync'ing with the total energy in the hydro update.

The source of error in this test comes from the radiation subcycle it self and has nothing to do with hydro update. Here is the relavent code:

In the end of the first radiation subcycle:

consNew(i, j, k, gasInternalEnergy_index) = Egas_guess;
consNew(i, j, k, gasEnergy_index) = ComputeEgasFromEint(rho, x1GasMom1, x2GasMom1, x3GasMom1, Egas_guess);

In the beginning of the second radiation subcycle:

const double Egastot0 = consPrev(i, j, k, gasEnergy_index);
Egas0 = ComputeEintFromEgas(rho, x1GasMom0, x2GasMom0, x3GasMom0, Egastot0);

The internal energy is lost in this back-and-forth calculation.

I don't think SyncDualEnergy will solve this problem. Here is the relavent code in SyncDualEnergy:

		// Li et al. sync method
		// replace Eint with Eint_cons == (Etot - Ekin) if (Eint_cons / E) > eta
		if (Eint_cons > eta * Etot) {
			consVar[bx](i, j, k, internalEnergy_index) = Eint_cons;
		} else { // non-conservative sync
			consVar[bx](i, j, k, internalEnergy_index) = Eint_aux;
			consVar[bx](i, j, k, energy_index) = Eint_aux + Ekin;
		}

The problem is (Eint_aux + Ekin) - Ekin = 0.0 when Eint_aux << Ekin. So, I think we have to use the auxiliary internal energy in the beginning of radiation update.

[BenWibking]

The dominant source of error in the internal energy is due to truncation errors, not round off errors. These are not the same, and the distinction is very important.

I believe it is roundoff error in this case. What happens in this test is that the gas has an equilibrium temperature of 1 because I set the radiation temperature to 1. However, because the specific heat of radiation is very large, due to conservation of momentum, the gas momentum increases to 10^13 after the first step. Then, in the second time step, because the internal energy is derived from subtracting the total energy from the kinetic energy, the result is 0 due to roundoff error.

If it is roundoff error due to catastrophic cancellation, that means that the kinetic energy is $\sim 10^{15}$ times larger than the internal energy. The Mach number of that flow is unphysically enormous, and no code can deal with that (correctly). Not even molecular gas moving in a cosmological simulation has a Mach number that large (1000 km/s / 0.2 km/s ~ Mach 5000).

Unphysically enormous: a kinetic to internal ratio of $10^{15}$ implies a Mach number of order $10^{7.5}$. For molecular gas at 10 K, that implies a velocity of $\sim 20$ times the speed of light.

[chongchonghe]

OK, that totally makes sense. I will change the parameters in this test to make the Mach number more realistic.

However, I still don't understand why you can't use the auxiliary internal energy. The code I pasted above from SyncDualEnergy basically enforces that consVar[bx](i, j, k, internalEnergy_index) stores a value that is exactly what you would obtain from ComputeEintFromEgas, so, why not just use consVar[bx](i, j, k, internalEnergy_index) and save some calculation (and also avoid some roundoff error, although not significant)?

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[BenWibking]

Just two minor comments below.

[sonarqubecloud]

Quality Gate passed

Issues
7 New issues
0 Accepted issues

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

See analysis details on SonarCloud

[chongchonghe]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).