this works really well now

zingale · Dec 13, 2023 · 20ba8c9 · 20ba8c9
1 parent 7660feb
commit 20ba8c9
Show file tree

Hide file tree

Showing 3 changed files with 24 additions and 10 deletions.
diff --git a/nse_tabular/nse_eos.H b/nse_tabular/nse_eos.H
@@ -20,7 +20,16 @@
 /// This function inverts this form of the EOS to find the T
 /// that satisfies the EOS and NSE.
 ///
-/// Note: T and abar come in as initial guesses and are updated
+/// The basic idea is that Abar and Zbar are both functions of
+/// rho, T, Ye through the NSE table, so we express the energy
+/// as:
+///
+///      e = e(rho, T, Abar(rho, T, Ye), Zbar(rho, T, Ye)
+///
+/// and NR on that.  Note that Zbar = Ye Abar, so we can group
+/// those derivative terms together.
+///
+/// T and abar come in as initial guesses and are updated
 /// on output
 ///
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
@@ -31,7 +40,7 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye,
     using namespace amrex;
     using namespace AuxZero;
 
-    const Real ttol{1.e-8_rt};
+    const Real ttol{1.e-6_rt};
     const int max_iter{100};
     const Real eps{1.e-8_rt};
 
@@ -67,18 +76,15 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye,
 
         Real f = eos_state.e - e_in;
 
-        // perturb T and call NSE again to get perturbed Abar
-        nse_state.T *= (1.0_rt + eps);
-        nse_interp(nse_state, skip_X_fill);
-
-        // estimate dAbar/dT
-        Real dabar_dT = (nse_state.abar - abar_old) / (eps * T);
+        Real dabar_dT = nse_interp_dT(T, rho, Ye, nse_table::abartab);
 
         // compute the correction to our guess
-        Real dT = -f / (eos_state.dedT + eos_state.dedA * dabar_dT);
+
+        Real dT = -f / (eos_state.dedT + eos_state.dedA * dabar_dT
+                                  + Ye * eos_state.dedZ * dabar_dT);
 
         // update the temperature
-        T = std::clamp(T + dT, 0.5 * T, 2.0 * T);
+        T = std::clamp(T + dT, 0.25 * T, 4.0 * T);
 
         // check convergence
 

diff --git a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out
@@ -65,4 +65,9 @@ dbea/dT = -6.886297725e-12
 now using derivative of the interpolant
 dAbar/dT = -1.070778122e-09
 dbea/dT = -6.886272826e-12
+
+EOS e consistency check: 1.395273834e+18 1.388449367e+18
+updated T: 6394612134
+change in abar: 55.60621897 50.27196857
+EOS e consistency check: 1.388449367e+18 1.388449367e+18
 AMReX (23.12-8-g43d71da32fa4) finalized
diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H
@@ -239,10 +239,13 @@ void nse_cell_c()
     eos_state.aux[iye] = ye;
     eos_state.aux[iabar] = abar_orig;
 
+    Real abar_start = eos_state.aux[iabar];
+
     nse_T_abar_from_e(eos_state.rho, eos_state.e, eos_state.aux[iye],
                       eos_state.T, eos_state.aux[iabar]);
 
     std::cout << "updated T: " << eos_state.T << std::endl;
+    std::cout << "change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl;
 
     // now check if we get back the correct e!