From 02a2d50bd2c4036cadba893d2e76c447fd14d6c6 Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Thu, 30 Nov 2023 13:03:46 -0500
Subject: [PATCH] switch the massive star to use (rho, P) to map onto the
 domain (#2658)

this should help with HSE:
---
 Exec/science/massive_star/_prob_params        |  3 +++
 .../problem_initialize_state_data.H           | 26 ++++++++++++-------
 2 files changed, 20 insertions(+), 9 deletions(-)

diff --git a/Exec/science/massive_star/_prob_params b/Exec/science/massive_star/_prob_params
index f763030baa..888bf80000 100644
--- a/Exec/science/massive_star/_prob_params
+++ b/Exec/science/massive_star/_prob_params
@@ -1,7 +1,10 @@
 model_name      character       ""        y
 
+interpolate_pres   int          0         y
+
 perturb_model   int             0         y
 velpert_amplitude  real         1.e5      y
 velpert_scale   real            1.e7      y
 
 isi28           int            -1
+
diff --git a/Exec/science/massive_star/problem_initialize_state_data.H b/Exec/science/massive_star/problem_initialize_state_data.H
index 51b4a9e38d..5850ba23b4 100644
--- a/Exec/science/massive_star/problem_initialize_state_data.H
+++ b/Exec/science/massive_star/problem_initialize_state_data.H
@@ -33,6 +33,7 @@ void problem_initialize_state_data (int i, int j, int k,
 
     state(i,j,k,URHO) = interpolate(dist, model::idens);
     state(i,j,k,UTEMP) = interpolate(dist, model::itemp);
+    Real pres = interpolate(dist, model::ipres);
     for (int n = 0; n < NumSpec; n++) {
         state(i,j,k,UFS+n) = amrex::max(interpolate(dist, model::ispec+n), small_x);
     }
@@ -116,22 +117,29 @@ void problem_initialize_state_data (int i, int j, int k,
         }
     }
 
-    // now call the EOS -- reload the eos_state since composition may
-    // have changed
+    // now call the EOS -- we need to use an eos_t here since we are
+    // going to take pressure as an input
 
-    burn_state.rho = state(i,j,k,URHO);
-    burn_state.T = state(i,j,k,UTEMP);
+    eos_t eos_state;
+    eos_state.rho = state(i,j,k,URHO);
+    eos_state.T = state(i,j,k,UTEMP);
+    eos_state.p = pres;
     for (int n = 0; n < NumSpec; n++) {
-        burn_state.xn[n] = state(i,j,k,UFS+n);
+        eos_state.xn[n] = state(i,j,k,UFS+n);
     }
     for (int n = 0; n < NumAux; n++) {
-        burn_state.aux[n] = state(i,j,k,UFX+n);
+        eos_state.aux[n] = state(i,j,k,UFX+n);
     }
 
-    eos(eos_input_rt, burn_state);
+    if (problem::interpolate_pres == 1) {
+        eos(eos_input_rp, eos_state);
+        state(i,j,k,UTEMP) = eos_state.T;
+    } else {
+        eos(eos_input_rt, eos_state);
+    }
 
-    state(i,j,k,UEINT) = state(i,j,k,URHO) * burn_state.e;
-    state(i,j,k,UEDEN) = state(i,j,k,URHO) * burn_state.e;
+    state(i,j,k,UEINT) = state(i,j,k,URHO) * eos_state.e;
+    state(i,j,k,UEDEN) = state(i,j,k,URHO) * eos_state.e;
 
     // add density scaling to composition