diff --git a/Exec/science/massive_star/problem_initialize_state_data.H b/Exec/science/massive_star/problem_initialize_state_data.H index 0333d7da82..68fd2acfd8 100644 --- a/Exec/science/massive_star/problem_initialize_state_data.H +++ b/Exec/science/massive_star/problem_initialize_state_data.H @@ -7,6 +7,7 @@ #ifdef NSE_TABLE #include #include +#include #endif AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -128,8 +129,33 @@ void problem_initialize_state_data (int i, int j, int k, } if (problem::interpolate_pres == 1) { - eos(eos_input_rp, eos_state); - state(i,j,k,UTEMP) = eos_state.T; + + // we need to get T from P, rho, but consistent with NSE (if + // we are in NSE) + if (nse_check) { + // this will also change the composition, since the new T + // alters the NSE + nse_T_abar_from_p(eos_state.rho, eos_state.p, eos_state.aux[AuxZero::iye], + eos_state.T, eos_state.aux[AuxZero::iabar]); + state(i,j,k,UTEMP) = eos_state.T; + state(i,j,k,UFX+AuxZero::iabar) = eos_state.aux[AuxZero::iabar]; + + // now call the EOS with the new T and abar to get e + eos(eos_input_rt, eos_state); + + // finally, get the updated B/A + nse_table_t nse_state; + nse_state.T = state(i,j,k,UTEMP); + nse_state.rho = state(i,j,k,URHO); + nse_state.Ye = state(i,j,k,UFX+AuxZero::iye); + nse_interp(nse_state); + + state(i,j,k,UFX+AuxZero::ibea) = nse_state.bea; + + } else { + eos(eos_input_rp, eos_state); + state(i,j,k,UTEMP) = eos_state.T; + } } else { eos(eos_input_rt, eos_state); }