update self-consistent NSE SDC update with Tabular NSE

[zhichen3]

address issue #1543
This depends on pynucastro/pynucastro#739

[zhichen3]

Detonation with this pr:

but it was a lot slower. Took ~400s runtime.

For comparison, I could run the same problem about ~100s with NSE before, and its about the same without NSE

[zhichen3]

Detonation with this pr also gets statements like:
burn entered NSE during integration (after 117 steps), zone = (227, 0, 0) recovering burn failure in NSE, zone = (227, 0, 0)
throughout all zones behind the shock.

Probably due to higher energy release behind the shock compared to the version currently on development. This causes a steeper change in composition and puts it out of NSE at the beginning of the next step?

Previously, we would only get these for zones near the shock front and for other zones behind the shock its simply in NSE. It appears that it made it drastically slower.

[zhichen3]

I realized that the current version on development is likely missing a factor of 1/dt in the definition of enuc_dot, but it seems to produce consistent results when I don't use NSE, which is weird.

[zhichen3]

Update with some progress:
With current version, I have:

 density                         3.575417e+19   1.942      9.303788e+18   1.980       2.35768e+18    1.601      7.773908e+17
 xmom                            1.804038e+28   2.065      4.311366e+27   1.983      1.090598e+27    1.902      2.917623e+26
 ymom                            1.803993e+28   2.065      4.311497e+27   1.983      1.090619e+27    1.902      2.917184e+26
 rho_E                           5.395202e+37   1.953      1.393596e+37   2.006      3.470112e+36    1.864      9.535919e+35
 rho_e                           5.323479e+37   1.953      1.375263e+37   2.006      3.424472e+36    1.864      9.404656e+35
 Temp                            5.891681e+20   1.996      1.477084e+20   2.013      3.658885e+19    1.169      1.627243e+19
 rho_He4                         2.839664e+18   1.852      7.863418e+17   1.921      2.076376e+17    1.393      7.903603e+16
 rho_Cr48                        1.238844e+18   2.007      3.081691e+17   1.740      9.226988e+16    0.532      6.379937e+16
 rho_Fe52                         1.43196e+19   1.892      3.857925e+18   1.750      1.146709e+18    0.820      6.496755e+17
 rho_Ni56                        3.878334e+20   2.038      9.443133e+19   2.077       2.23866e+19    2.139      5.083817e+18
 rho_enuc                         5.75844e+38   1.087      2.711178e+38   1.501      9.576835e+37    0.598      6.329164e+37

Not sure why convergence looks so bad when I included 512 resolution run.

And detonation looks like:

which is really fuzzy. Also encountered burn failures so I can't even run til the end.

[zhichen3]

Current state of this pr:

Its takes ~400s to run detonation with NSE and ~250s without NSE. (just run with USE_NSE_NET=TRUE and USE_SIMPLIFIED_SDC=TRUE with inputs-det-x.nse_net)

Profile for detonation, without NSE:

With NSE:

[zhichen3]

previously only ran up to 128 res because its slow, but I decided to run 256 res to check the convergence, and for whatever reason, its not good...

 density                          3.57458e+19   1.954      9.227014e+18   1.783      2.681491e+18
 xmom                            1.802359e+28   2.060      4.323031e+27   2.001       1.07969e+27
 ymom                            1.802364e+28   2.060      4.322971e+27   2.001      1.079702e+27
 rho_E                           5.404931e+37   1.964      1.385427e+37   2.017      3.423357e+36
 rho_e                           5.333054e+37   1.964      1.367084e+37   2.018      3.375406e+36
 Temp                            5.856433e+20   1.990       1.47429e+20   1.623       4.78547e+19
 rho_He4                         2.828525e+18   1.844      7.876808e+17   1.470      2.843232e+17
 rho_Cr48                        1.222316e+18   2.013      3.027591e+17   0.713       1.84743e+17
 rho_Fe52                        1.426529e+19   1.904      3.812699e+18   1.250      1.602976e+18
 rho_Ni56                        3.876178e+20   2.037       9.44657e+19   2.083      2.229842e+19
 rho_enuc                        5.431598e+38   1.071      2.585564e+38   0.988      1.303703e+38

By decreasing the tolerance for temperature convergence from 1.e-6 to 1.e-7, I get:

 density                         3.576225e+19   1.947      9.278317e+18   2.028      2.274303e+18    1.999      5.690442e+17
 xmom                             1.80319e+28   2.045      4.370978e+27   2.013      1.083208e+27    2.011       2.68692e+26
 ymom                            1.803184e+28   2.045      4.371014e+27   2.013       1.08321e+27    2.011      2.686647e+26
 rho_E                           5.406496e+37   1.954      1.395561e+37   2.025      3.429176e+36    2.012      8.499385e+35
 rho_e                           5.334729e+37   1.954      1.377077e+37   2.026      3.381836e+36    2.012      8.384062e+35
 Temp                             5.85942e+20   2.000      1.464391e+20   2.109      3.395133e+19    2.024       8.34525e+18
 rho_He4                         2.828242e+18   1.860       7.78868e+17   1.991      1.958958e+17    2.083      4.623912e+16
 rho_Cr48                        1.225026e+18   2.039      2.981315e+17   2.145      6.738324e+16    1.655      2.139777e+16
 rho_Fe52                        1.429159e+19   1.901      3.826192e+18   1.997      9.587731e+17    1.754      2.843321e+17
 rho_Ni56                        3.876275e+20   2.038       9.44107e+19   2.084      2.226414e+19    2.208      4.820328e+18
 rho_enuc                        5.419401e+38   1.082      2.559245e+38   1.708      7.832314e+37    0.817      4.446305e+37

With ttol=1.e-8, I get

 density                         3.576018e+19   1.948      9.268791e+18   2.033      2.265095e+18    2.042      5.498526e+17
 xmom                            1.803014e+28   2.046      4.365948e+27   2.002      1.090134e+27    1.988       2.74821e+26
 ymom                            1.803014e+28   2.046      4.365947e+27   2.002      1.090141e+27    1.988      2.748503e+26
 rho_E                           5.406225e+37   1.955      1.394656e+37   2.021      3.436645e+36    1.999        8.5966e+35
 rho_e                           5.334454e+37   1.955      1.376181e+37   2.022      3.389202e+36    1.999      8.479188e+35
 Temp                            5.859102e+20   2.000      1.464346e+20   2.131      3.344017e+19    2.278      6.896348e+18
 rho_He4                         2.828401e+18   1.859      7.794693e+17   2.014      1.930269e+17    2.227      4.121872e+16
 rho_Cr48                        1.224792e+18   2.041       2.97577e+17   2.214      6.413942e+16    2.242       1.35617e+16
 rho_Fe52                         1.42888e+19   1.903      3.820818e+18   1.999      9.559052e+17    2.076      2.267084e+17
 rho_Ni56                        3.876303e+20   2.038      9.440292e+19   2.086      2.223161e+19    2.228       4.74419e+18
 rho_enuc                        5.418569e+38   1.083       2.55734e+38   1.692      7.914355e+37    0.929      4.156051e+37

After using ttol=1.e-8, I would get abs(H)=HMIN failure for detonation, then I tried switching to 5-point stencil for constructing the derivative dabardT, then it goes away.

Then the new detonation profile is:

[zingale]

I did a comparison to the tabular case and this seems to follow the same flow.