add a pressure interface

AMReX-Astro · Dec 18, 2023 · da880f0 · da880f0
1 parent 20ba8c9
commit da880f0
Show file tree

Hide file tree

Showing 3 changed files with 208 additions and 44 deletions.
diff --git a/nse_tabular/nse_eos.H b/nse_tabular/nse_eos.H
@@ -42,7 +42,6 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye,
 
     const Real ttol{1.e-6_rt};
     const int max_iter{100};
-    const Real eps{1.e-8_rt};
 
     // we need the full EOS type, since we need de/dA
     //eos_extra_t eos_state;
@@ -100,4 +99,94 @@ nse_T_abar_from_e(const Real rho, const Real e_in, const Real Ye,
 
 }
 
+
+///
+/// if we are in NSE, then the entire thermodynamic state is just
+/// a function of rho, T, Ye.  We can write the pressure as:
+///
+///    p = [(rho, T, Y_e, Abar(rho, T, Ye))
+///
+/// where we note that Abar is a function of those same inputs.
+/// This function inverts this form of the EOS to find the T
+/// that satisfies the EOS and NSE.
+///
+/// The basic idea is that Abar and Zbar are both functions of
+/// rho, T, Ye through the NSE table, so we express the pressure
+/// as:
+///
+///      p = p(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
+void
+nse_T_abar_from_p(const Real rho, const Real p_in, const Real Ye,
+                  Real& T, Real& abar) {
+
+    using namespace amrex;
+    using namespace AuxZero;
+
+    const Real ttol{1.e-6_rt};
+    const int max_iter{100};
+
+    // we need the full EOS type, since we need de/dA
+    //eos_extra_t eos_state;
+
+    bool converged{false};
+
+    int iter{};
+
+    nse_table_t nse_state;
+
+    while (not converged && iter < max_iter) {
+
+        // call NSE table to get Abar
+        nse_state.T = T;
+        nse_state.rho = rho;
+        nse_state.Ye = Ye;
+
+        constexpr bool skip_X_fill{true};
+        nse_interp(nse_state, skip_X_fill);
+        Real abar_old = nse_state.abar;
+
+        // call the EOS with the initial guess for T
+        eos_extra_t eos_state;
+        eos_state.rho = rho;
+        eos_state.T = T;
+        eos_state.aux[iye] = Ye;
+        eos_state.aux[iabar] = abar_old;
+        eos(eos_input_rt, eos_state);
+
+        // f is the quantity we want to zero
+
+        Real f = eos_state.p - p_in;
+
+        Real dabar_dT = nse_interp_dT(T, rho, Ye, nse_table::abartab);
+
+        // compute the correction to our guess
+
+        Real dT = -f / (eos_state.dpdT + eos_state.dpdA * dabar_dT
+                                  + Ye * eos_state.dpdZ * dabar_dT);
+
+        // update the temperature
+        T = std::clamp(T + dT, 0.25 * T, 4.0 * T);
+
+        // check convergence
+
+        if (std::abs(dT) < ttol * T) {
+            converged = true;
+        }
+        iter++;
+    }
+
+    // T is set to the last T
+    // we just need to save abar for output
+    abar = nse_state.abar;
+
+}
+
 #endif
diff --git a/unit_test/test_nse_interp/ci-benchmarks/aprox19.out b/unit_test/test_nse_interp/ci-benchmarks/aprox19.out
@@ -1,5 +1,5 @@
-Initializing AMReX (23.12-8-g43d71da32fa4)...
-AMReX (23.12-8-g43d71da32fa4) initialized
+Initializing AMReX (23.11-5-gd36463103dae)...
+AMReX (23.11-5-gd36463103dae) initialized
 starting the single zone burn...
 reading the NSE table (C++) ...
 rho, T, Ye = 1230000000 5180000000 0.472
@@ -66,8 +66,13 @@ 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
+EOS e consistency check (old method): 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
+EOS e consistency check (new method): 1.388449367e+18 1.388449367e+18
+
+EOS p consistency check (old method): 6.622132093e+26 6.57785215e+26
+updated T: 6466848772
+change in abar: 55.60621897 49.50670173
+EOS p consistency check (new method): 6.57785215e+26 6.57785215e+26
+AMReX (23.11-5-gd36463103dae) finalized
diff --git a/unit_test/test_nse_interp/nse_cell.H b/unit_test/test_nse_interp/nse_cell.H
@@ -193,64 +193,134 @@ void nse_cell_c()
     //
     // T', rho, Ye, Abar' -> e' ???
 
-    // first get the abar consistent with our inputs and find e
+    {
 
-    nse_state.T = temperature;
-    nse_state.rho = density;
-    nse_state.Ye = ye;
+        // first get the abar consistent with our inputs and find e
 
-    nse_interp(nse_state);
+        nse_state.T = temperature;
+        nse_state.rho = density;
+        nse_state.Ye = ye;
 
-    Real abar_orig = nse_state.abar;
+        nse_interp(nse_state);
 
-    eos_t eos_state;
+        Real abar_orig = nse_state.abar;
 
-    eos_state.T = temperature;
-    eos_state.rho = density;
-    eos_state.aux[iye] = nse_state.Ye;
-    eos_state.aux[iabar] = nse_state.abar;
+        eos_t eos_state;
 
-    eos(eos_input_rt, eos_state);
+        eos_state.T = temperature;
+        eos_state.rho = density;
+        eos_state.aux[iye] = nse_state.Ye;
+        eos_state.aux[iabar] = nse_state.abar;
 
-    Real e_orig = eos_state.e;
+        eos(eos_input_rt, eos_state);
 
-    // now perturn e and find T, then redo NSE to get abar and finally
-    // see if we get back our new e
+        Real e_orig = eos_state.e;
 
-    Real e_new = eos_state.e * 1.05;
+        // now perturn e and find T, then redo NSE to get abar and finally
+        // see if we get back our new e
 
-    eos_state.e = e_new;
-    eos(eos_input_re, eos_state);
+        Real e_new = eos_state.e * 1.05;
 
-    nse_state.T = eos_state.T;
-    nse_interp(nse_state);
+        eos_state.e = e_new;
+        eos(eos_input_re, eos_state);
+
+        nse_state.T = eos_state.T;
+        nse_interp(nse_state);
+
+        eos_state.aux[iabar] = nse_state.abar;
+        eos(eos_input_rt, eos_state);
+
+        std::cout << "EOS e consistency check (old method): " << eos_state.e << " " << e_new << std::endl;
 
-    eos_state.aux[iabar] = nse_state.abar;
-    eos(eos_input_rt, eos_state);
+        // attempt 2:
+        // now we try the new interface.  This effectively does:
+        // e', rho, Ye -> Abar', T'
 
-    std::cout << "EOS e consistency check: " << eos_state.e << " " << e_new << std::endl;
+        eos_state.T = temperature;
+        eos_state.e = e_new;
+        eos_state.rho = density;
+        eos_state.aux[iye] = ye;
+        eos_state.aux[iabar] = abar_orig;
 
-    // now we try the new interface.  This effectively does:
-    // e', rho, Ye -> Abar', T'
+        Real abar_start = eos_state.aux[iabar];
 
-    eos_state.T = temperature;
-    eos_state.e = e_new;
-    eos_state.rho = density;
-    eos_state.aux[iye] = ye;
-    eos_state.aux[iabar] = abar_orig;
+        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!
+
+        eos(eos_input_rt, eos_state);
+
+        std::cout << "EOS e consistency check (new method): " << eos_state.e << " " << e_new << std::endl;
+        std::cout << std::endl;
+    }
+
+
+    {
 
-    Real abar_start = eos_state.aux[iabar];
+        // now redo it for pressure
 
-    nse_T_abar_from_e(eos_state.rho, eos_state.e, eos_state.aux[iye],
-                      eos_state.T, eos_state.aux[iabar]);
+        nse_state.T = temperature;
+        nse_state.rho = density;
+        nse_state.Ye = ye;
 
-    std::cout << "updated T: " << eos_state.T << std::endl;
-    std::cout << "change in abar: " << abar_start << " " << eos_state.aux[iabar] << std::endl;
+        nse_interp(nse_state);
 
-    // now check if we get back the correct e!
+        Real abar_orig = nse_state.abar;
 
-    eos(eos_input_rt, eos_state);
+        eos_t eos_state;
 
-    std::cout << "EOS e consistency check: " << eos_state.e << " " << e_new << std::endl;
+        eos_state.T = temperature;
+        eos_state.rho = density;
+        eos_state.aux[iye] = nse_state.Ye;
+        eos_state.aux[iabar] = nse_state.abar;
+
+        eos(eos_input_rt, eos_state);
+
+        Real p_orig = eos_state.p;
+
+        // now perturn p and find T, then redo NSE to get abar and finally
+        // see if we get back our new p
+
+        Real p_new = eos_state.p * 1.05;
+
+        eos_state.p = p_new;
+        eos(eos_input_rp, eos_state);
+
+        nse_state.T = eos_state.T;
+        nse_interp(nse_state);
+
+        eos_state.aux[iabar] = nse_state.abar;
+        eos(eos_input_rt, eos_state);
+
+        std::cout << "EOS p consistency check (old method): " << eos_state.p << " " << p_new << std::endl;
+
+        // attempt 2:
+        // now we try the new interface.  This effectively does:
+        // p', rho, Ye -> Abar', T'
+
+        eos_state.T = temperature;
+        eos_state.p = p_new;
+        eos_state.rho = density;
+        eos_state.aux[iye] = ye;
+        eos_state.aux[iabar] = abar_orig;
+
+        Real abar_start = eos_state.aux[iabar];
+
+        nse_T_abar_from_p(eos_state.rho, eos_state.p, 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 p!
+
+        eos(eos_input_rt, eos_state);
+
+        std::cout << "EOS p consistency check (new method): " << eos_state.p << " " << p_new << std::endl;
+    }
 
 }