update nse_sdc_burn for NSE_NET

[zhichen3]

Addresses the following:

1. Fix the issue that abar wasn't calculated correctly. Needed to take the reciprocal in the end.
2. Switch to use central difference rather than 5-stencil to compute dabardT. Convergence test shows that central difference is sufficient for 2nd order.
3. Construct the reactive source term for rho e directly, instead of evolving for tau and do a finite difference to get the derivative.

[zhichen3]

After fixing the abar calculation, CI test shows that we fully recover the internal energy.

[zhichen3]

First running the usual nse_test.

With this PR but don't change how we calculate drhoedt, i.e. fixed abar calculation and used central difference for dabardT, but still evolve for tau and find drhoedt

And convergence rate:

| field    |   \epsilon[32 \rightarrow 64] |  rate |  \epsilon[64 \rightarrow 128] |  rate | \epsilon[128 \rightarrow 256] |  rate | \epsilon[256 \rightarrow 512] |
|----------+--------------+-------+--------------+-------+--------------+-------+--------------|
| density  | 3.576264e+19 | 1.948 | 9.266531e+18 | 2.033 | 2.264561e+18 | 2.044 | 5.493195e+17 |
| xmom     | 1.802962e+28 | 2.046 |  4.36657e+27 | 2.004 | 1.088621e+27 | 1.993 | 2.735492e+26 |
| ymom     | 1.802963e+28 | 2.046 | 4.366637e+27 | 2.004 | 1.088603e+27 | 1.993 | 2.735202e+26 |
| rho_E    | 5.406592e+37 | 1.955 | 1.394247e+37 | 2.020 | 3.436507e+36 | 2.005 | 8.562214e+35 |
| rho_e    |  5.33494e+37 | 1.955 | 1.375767e+37 | 2.021 | 3.389099e+36 | 2.005 | 8.444978e+35 |
| Temp     | 5.858881e+20 | 2.001 |  1.46416e+20 | 2.131 | 3.341744e+19 | 2.272 | 6.917425e+18 |
| rho_He4  | 2.829101e+18 | 1.860 | 7.794294e+17 | 2.015 | 1.927953e+17 | 2.222 | 4.131369e+16 |
| rho_Cr48 | 1.222589e+18 | 2.042 | 2.968651e+17 | 2.213 | 6.401146e+16 | 2.223 | 1.371048e+16 |
| rho_Fe52 | 1.433159e+19 | 1.903 | 3.832956e+18 | 1.998 | 9.594055e+17 | 2.095 | 2.245903e+17 |
| rho_Ni56 | 3.875922e+20 | 2.038 |  9.43913e+19 | 2.086 | 2.222705e+19 | 2.226 |  4.75007e+18 |
| rho_enuc | 5.422511e+38 | 1.093 | 2.541502e+38 | 1.655 | 8.069996e+37 | 0.944 | 4.195144e+37 |

Now with this PR, i.e. just use Ydot_weak to find drhoedt:

| field    |   \epsilon[32 \rightarrow 64] |  rate |  \epsilon[64 \rightarrow 128] |  rate | \epsilon[128 \rightarrow 256] |  rate | \epsilon[256 \rightarrow 512] |
|----------+--------------+-------+--------------+-------+--------------+-------+--------------|
| density  | 3.561917e+19 | 1.938 | 9.293718e+18 | 2.026 |  2.28117e+18 | 2.046 | 5.525218e+17 |
| xmom     | 1.833654e+28 | 2.042 | 4.452669e+27 | 2.001 | 1.112583e+27 | 1.996 | 2.789012e+26 |
| ymom     | 1.833654e+28 | 2.042 | 4.452669e+27 | 2.001 | 1.112583e+27 | 1.996 | 2.789012e+26 |
| rho_E    | 5.519616e+37 | 1.964 | 1.414586e+37 | 2.023 |  3.48144e+36 | 2.009 | 8.646525e+35 |
| rho_e    | 5.441977e+37 | 1.963 | 1.395544e+37 | 2.023 | 3.433881e+36 | 2.010 | 8.526743e+35 |
| Temp     | 5.663637e+20 | 1.974 | 1.442111e+20 | 2.119 | 3.320305e+19 | 2.280 | 6.835962e+18 |
| rho_He4  | 2.701604e+18 | 1.836 | 7.567474e+17 | 2.000 | 1.891694e+17 | 2.223 | 4.052204e+16 |
| rho_Cr48 | 1.135922e+18 | 1.952 | 2.935108e+17 | 2.198 | 6.397405e+16 | 2.217 |  1.37644e+16 |
| rho_Fe52 | 1.477479e+19 | 1.901 | 3.956677e+18 | 2.018 | 9.769028e+17 | 2.114 | 2.257332e+17 |
| rho_Ni56 | 3.855638e+20 | 2.035 | 9.410037e+19 | 2.084 | 2.219803e+19 | 2.228 | 4.737176e+18 |
| rho_enuc | 2.957712e+33 | 1.737 | 8.870316e+32 | 1.679 | 2.770562e+32 | 1.434 | 1.025225e+32 |

[zhichen3]

Now with detonation:

Running with direct integration:

Running with this PR but without modifying how we calculate drhoedt:

Completely running with this PR:

With this PR we see that the temperature profile matches with direct integration exactly. We only see energy generation from zones that are not in NSE, i.e. just regular burning. And behind the shock is in NSE.
Since we calculate drhoedt directly from Ydot_weak, e_nuc contribution is identically zero for network without weak rates, and only contributing factors are from advection terms and neutrinos. I think it is a reasonable thing since if we're in NSE, then there is no energy generation from strong reactions?

With this PR, I'm able to run detonation at a similar runtime (both ~65s with 4 processors) compared to direct integration now. Before it was running substantially slower.

[zingale]

this is really cool. Do you understand why we don';t need to do the same finite-difference estimate of the update as we do with the tabular NSE here?

[zhichen3]

I'm not sure. One thing that's not being tested here is I don't have any weak rates in the network, so the energy generation from reaction is exactly zero. I should probably make a network and test that.

[zhichen3]

updated the docs