From 03b087e8e2088f1bba2125d39325975f3685fe22 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Thu, 18 Jul 2024 13:24:02 -0400
Subject: [PATCH] update the test_nse_net test  to also test the NSE EOS this
 follows the same test that is done in the test_nse_interp test.

---
 unit_test/test_nse_net/burn_cell.H | 55 ++++++++++++++++++++++++++++++
 1 file changed, 55 insertions(+)

diff --git a/unit_test/test_nse_net/burn_cell.H b/unit_test/test_nse_net/burn_cell.H
index 196684312..fc01490c5 100644
--- a/unit_test/test_nse_net/burn_cell.H
+++ b/unit_test/test_nse_net/burn_cell.H
@@ -66,5 +66,60 @@ void burn_cell_c()
       std::cout << "We're not in NSE. " << std::endl;
     }
 
+    // now test the EOS+NSE function.  The basic idea is the following:
+    //
+    // 1. find the e corresponding to the current NSE state (eos_input_rt)
+    // 2. perturb e and then call the NSE EOS to get the updated T
+    // 3. call the EOS with this new T (eos_input_rt) and updated NSE
+    //    composition and see if we recover the same perturbed e
+
+    eos_t eos_state;
+    eos_state.T = state.T;
+    eos_state.rho = state.rho;
+    for (int n = 0; n < NumSpec; ++n) {
+        eos_state.xn[n] = state.y[SFS+n] / state.rho;
+    }
+
+    // get the initial e and abar corresponding to this T
+
+    eos(eos_input_rt, eos_state);
+
+    amrex::Real abar_start = eos_state.abar;
+    amrex::Real e_start = eos_state.e;
+
+    // now perturb e and find the T/abar that are consistent with it
+    // and NSE
+
+    amrex::Real e_new = eos_state.e * 1.05;
+
+    eos_state.e = e_new;
+
+    amrex::Real T_new{eos_state.T};
+    amrex::Real abar_new{};
+
+    nse_T_abar_from_e(eos_state.rho, eos_state.e, ye,
+                      T_new, abar_new,
+                      state.mu_p, state.mu_n);
+
+    std::cout << "change in T: " << eos_state.T << " " << T_new <<  std::endl;
+    std::cout << "change in abar: " << abar_start << " " << abar_new <<  std::endl;
+
+    // now try calling the EOS with this T and composition and see if
+    // we get back the same energy.  we need to recompute the NSE
+    // state here since we didn't get the mass fractions from the EOS call
+
+    state.T = T_new;
+    auto nse_state = get_actual_nse_state(state, 1.e-10, true);
+
+    for (int n = 0; n < NumSpec; ++n) {
+        eos_state.xn[n] = nse_state.xn[n];
+    }
+
+    eos_state.T = T_new;
+
+    eos(eos_input_rt, eos_state);
+
+    std::cout << "recovered energy: " << eos_state.e << " " << e_new << std::endl;
+
 }
 #endif