diff --git a/integration/nse_update_sdc.H b/integration/nse_update_sdc.H index bfc24ec230..ed4913bb27 100644 --- a/integration/nse_update_sdc.H +++ b/integration/nse_update_sdc.H @@ -21,6 +21,7 @@ #endif #ifdef NSE_NET #include +#include #endif using namespace amrex::literals; @@ -301,207 +302,311 @@ void sdc_nse_burn(BurnT& state, const amrex::Real dt) { } #else +/// Self-Consistent NSE VERSION + /// -/// update the state due to NSE changes for the simplified-SDC case -/// this version works with the self-consistent NSE +/// this acts as an explicit Euler step for the system (rho e, rho) +/// on input, *_source are the reactive sources at time t0 and on output +/// they are the sources at time t0+dt /// -template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void sdc_nse_burn(BurnT& state, const amrex::Real dt) { +void nse_derivs(const amrex::Real rho0, const amrex::Real rhoe0, + const amrex::Real rhoYe0, const amrex::Real T_old, + const amrex::Real dt, amrex::Real &mu_p, amrex::Real &mu_n, + amrex::Real &drhoyedt_weak, const amrex::Real drhoyedt_a, + amrex::Array1D &ydot_weak, + const amrex::Real *ydot_a, amrex::Real& drhoedt, + const amrex::Real T_fixed) { - state.success = true; - state.n_rhs = 0; - state.n_jac = 0; + // start with the current state and compute T - // store the initial mass fractions -- we will need these - // to compute the energy release. + amrex::Real T0; + amrex::Real Ye0 = rhoYe0 / rho0; + amrex::Real abar{0.0_rt}; - amrex::Real X_old[NumSpec]; + // Get initial temperature and abar for - for (int n = 0; n < NumSpec; ++n) { - X_old[n] = state.y[SFS+n] / state.y[SRHO]; + if (T_fixed > 0) { + T0 = T_fixed; + } else { + // initial guess for eos + T0 = T_old; + amrex::Real e0 = rhoe0 / rho0; + nse_T_abar_from_e(rho0, e0, Ye0, T0, abar, mu_p, mu_n); } - // density and momentum have no reactive sources - amrex::Real rho_old = state.y[SRHO]; + // initialize burn_state for NSE state - state.y[SRHO] += dt * state.ydot_a[SRHO]; - state.y[SMX] += dt * state.ydot_a[SMX]; - state.y[SMY] += dt * state.ydot_a[SMY]; - state.y[SMZ] += dt * state.ydot_a[SMZ]; + burn_t burn_state; - // if we are doing drive_initial_convection, we want to use - // the temperature that comes in through T_fixed + burn_state.T = T0; + burn_state.rho = rho0; + burn_state.y_e = Ye0; + burn_state.mu_p = mu_p; + burn_state.mu_n = mu_n; - amrex::Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T; + // get NSE state at t0 - // get the neutrino loss term -- we want to use the state that we - // came in here with, so the original Abar and Zbar - amrex::Real snu{0.0}; - amrex::Real dsnudt{0.0}; - amrex::Real dsnudd{0.0}; - amrex::Real dsnuda{0.0}; - amrex::Real dsnudz{0.0}; + auto nse_state0 = get_actual_nse_state(burn_state, 1.0e-10_rt, true); - amrex::Real abar{0.0}; - amrex::Real zbar{0.0}; - for (int n = 0; n < NumSpec; ++n) { - abar += X_old[n] * aion_inv[n]; - zbar += X_old[n] * zion[n] * aion_inv[n]; + // Still need to get abar for T_fixed > 0 + + if (T_fixed > 0) { + for (int n = 0; n < NumSpec; ++n) { + abar += nse_state0.xn[n] * aion_inv[n]; + } + abar = 1.0_rt / abar; } - abar = 1.0 / abar; - zbar *= abar; + + // compute the plasma neutrino losses + // here abar should already be in-sync with NSE + + amrex::Real snu{0.0_rt}; + amrex::Real dsnudt{0.0_rt}; + amrex::Real dsnudd{0.0_rt}; + amrex::Real dsnuda{0.0_rt}; + amrex::Real dsnudz{0.0_rt}; + + amrex::Real zbar = abar * Ye0; #ifdef NEUTRINOS { constexpr int do_derivatives = 0; - sneut5(T_in, rho_old, abar, zbar, - snu, dsnudt, dsnudd, dsnuda, dsnudz); + sneut5(T0, rho0, abar, zbar, + snu, dsnudt, dsnudd, + dsnuda, dsnudz); } #endif - amrex::Real snu_old = snu; + // Get molar fractions, ydots and neutrino loss due to weak rates - // if our network could return the evolution of Ye due to the - // weak interactions, we would evaluate the NSE state here and - // get dYe/dt. + amrex::Array1D Y; + amrex::Real e_nu {0.0_rt}; - //amrex::Real dyedt_old = 0.0; + for (int n = 1; n <= NumSpec; ++n) { + Y(n) = nse_state0.xn[n-1] * aion_inv[n-1]; + } + // Fill in ydot with only weak rates contributing - // predict the U^{n+1,*} state with only estimates of the X - // updates with advection to dt and the neutrino loss term in - // energy + get_ydot_weak(nse_state0, ydot_weak, e_nu, Y); - BurnT burn_state; - burn_state.T = T_in; // initial guess - burn_state.mu_p = state.mu_p; - burn_state.mu_n = state.mu_n; + /// + /// construct initial sources at t0 + /// - amrex::Real rhoX_source[NumSpec] = {0.0_rt}; + // Find drhoyedt_weak from weak reaction - amrex::Real rhoX_new[NumSpec]; + drhoyedt_weak = 0.0_rt; + for (int n = 0; n < NumSpec; ++n) { + drhoyedt_weak += rho0 * zion[n] * ydot_weak(n+1); + } - amrex::Real rhoe_new; - amrex::Real rho_enucdot = -rho_old * snu; + // find rhoe_source - amrex::Real rho_half = 0.5_rt * (rho_old + state.y[SRHO]); + amrex::Real rhoe_source{0.0_rt}; + for (int n = 1; n <= NumSpec; ++n) { + rhoe_source += rho0 * ydot_weak(n) * network::mion(n); + } + rhoe_source *= C::Legacy::enuc_conv2; - // predict the new aux for the first iteration -- this is really - // just including the advection bits + // include plasma neutrino terms - for (int n = 0; n < NumSpec; n++) { - rhoX_new[n] = state.y[SFS+n] + dt * state.ydot_a[SFS+n] + dt * rhoX_source[n]; + rhoe_source -= rho0 * snu; + + // include neutrino loss terms from weak rates + // Note that e_nu here is already negative so we add here. + + if (integrator_rp::nse_include_enu_weak == 1) { + rhoe_source += rho0 * e_nu; } - burn_t nse_state; + // + // evolve for eps * dt; + // - for (int iter = 0; iter < integrator_rp::nse_iters; iter++) { + amrex::Real tau = integrator_rp::nse_deriv_dt_factor * dt; - // update (rho e)^{n+1} based on the new energy generation rate - rhoe_new = state.y[SEINT] + dt * state.ydot_a[SEINT] + dt * rho_enucdot; + amrex::Real rho1 = rho0 + tau * ydot_a[SRHO]; + amrex::Real rhoe1 = rhoe0 + tau * (ydot_a[SEINT] + rhoe_source); + amrex::Real rhoYe1 = rhoYe0 + tau * (drhoyedt_weak + drhoyedt_a); - // call the EOS to get the updated T* + // compute the temperature at t0 + tau - amrex::Real T_new; - burn_state.rho = state.y[SRHO]; - burn_state.e = rhoe_new / state.y[SRHO]; - for (int n = 0; n < NumSpec; n++) { - burn_state.xn[n] = rhoX_new[n] / state.y[SRHO]; - burn_state.y[SFS+n] = rhoX_new[n]; - } + amrex::Real T1; + amrex::Real Ye1 = rhoYe1 / rho1; - if (state.T_fixed > 0) { - T_new = state.T_fixed; - } else { - eos(eos_input_re, burn_state); - T_new = burn_state.T; - } - burn_state.T = T_new; + if (T_fixed > 0) { + T1 = T_fixed; + } else { + // initial guess for eos + T1 = T0; + amrex::Real e1 = rhoe1 / rho1; + amrex::Real abar1; + nse_T_abar_from_e(rho1, e1, Ye1, T1, abar1, mu_p, mu_n); + } - // solve for the NSE state for this network we need to call - // get_nse_state with a burn_t. We will have it compute - // Ye from the input X's + // update burn_state at t0 + tau - nse_state = get_actual_nse_state(burn_state); + burn_state.T = T1; + burn_state.rho = rho1; + burn_state.y_e = Ye1; - // compute the energy release. The mass fractions in nse_state.xn[] - // include the advective parts, so first we need to remove that. + // call NSE at t0 + tau with - amrex::Real rhoX_tilde[NumSpec]; - for (int n = 0; n < NumSpec; ++n) { - rhoX_tilde[n] = state.y[SRHO] * nse_state.xn[n] - dt * state.ydot_a[SFS+n]; - } + // Note this new NSE state has advection contribution - //amrex::Real dyedt = 0.0_rt; // we can update this in the future by calling actual_rhs() + auto nse_state1 = get_actual_nse_state(burn_state, 1.0e-10_rt, true); - // we want to compute (rho eps) = - N_A c^2 sum{m_i (rhoX_tilde - rhoX_old) / A_i} - rho_enucdot = 0.0; - for (int n = 0; n < NumSpec; ++n) { - rho_enucdot += (rhoX_tilde[n] - rho_old * X_old[n]) * - network::mion(n+1) * aion_inv[n]; - } - rho_enucdot *= C::Legacy::enuc_conv2; + // construct the finite-difference approximation to the derivatives + // subtract off advection contribution to get pure drhoedt from reaction - // now get the updated neutrino term - abar = 0.0; - zbar = 0.0; - for (int n = 0; n < NumSpec; ++n) { - abar += nse_state.xn[n] * aion_inv[n]; - zbar += nse_state.xn[n] * zion[n] * aion_inv[n]; - } - abar = 1.0 / abar; - zbar *= abar; + drhoedt = 0.0_rt; + for (int n = 0; n < NumSpec; ++n) { + drhoedt += (nse_state1.y[SFS+n] - tau * ydot_a[SFS+n] - + nse_state0.y[SFS+n]) * aion_inv[n] * network::mion(n+1); + } + drhoedt *= C::Legacy::enuc_conv2 / tau; -#ifdef NEUTRINOS - constexpr int do_derivatives = 0; - sneut5(T_new, state.y[SRHO], abar, zbar, - snu, dsnudt, dsnudd, dsnuda, dsnudz); -#endif + // include neutrino terms again + + drhoedt -= rho1 * snu; + if (integrator_rp::nse_include_enu_weak == 1) { + drhoedt += rho1 * e_nu; + } - rho_enucdot -= 0.5_rt * rho_half * (snu_old + snu); +} - // update the new state for the next pass -- this should - // already implicitly have the advective portion included, - // since it was there when we called the NSE state - for (int n = 0; n < NumSpec; n++) { - rhoX_new[n] = state.y[SRHO] * nse_state.xn[n]; - } +/// +/// update the state due to NSE changes for the simplified-SDC case +/// this version works with the self-consistent NSE +/// + +template +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void sdc_nse_burn(BurnT& state, const amrex::Real dt) { + + state.success = true; + state.n_rhs = 0; + state.n_jac = 0; + + // store the initial state + + amrex::Real rho_old = state.y[SRHO]; + amrex::Real rhoe_old = state.y[SEINT]; + amrex::Real rhoye_old = rho_old * state.y_e; + amrex::Real T_old = state.T; + amrex::Real mu_p = state.mu_p; + amrex::Real mu_n = state.mu_n; + // density and momentum have no reactive sources + + state.y[SRHO] += dt * state.ydot_a[SRHO]; + state.y[SMX] += dt * state.ydot_a[SMX]; + state.y[SMY] += dt * state.ydot_a[SMY]; + state.y[SMZ] += dt * state.ydot_a[SMZ]; + + // do an RK2 integration + + // get the derivatives, rho_dyedt, drhoedt, ydot_weak at t = t^n + + amrex::Array1D ydot_weak; + amrex::Real drhoedt{0.0_rt}; + amrex::Real drhoyedt_weak{0.0_rt}; + amrex::Real drhoyedt_a{0.0_rt}; + + // find the advection contribution: drhoyedt_a and drhoabardt_a + + for (int n = 0; n < NumSpec; ++n) { + drhoyedt_a += state.ydot_a[SFS+n] * zion[n] * aion_inv[n]; } - // now update the aux quantities + nse_derivs(rho_old, rhoe_old, + rhoye_old, T_old, + dt, mu_p, mu_n, + drhoyedt_weak, drhoyedt_a, + ydot_weak, state.ydot_a, + drhoedt, state.T_fixed); - // the new mass fractions are just those that come from the table + // evolve to the midpoint in time + + amrex::Real rho_tmp = rho_old + 0.5_rt * dt * state.ydot_a[SRHO]; + amrex::Real rhoe_tmp = rhoe_old + 0.5_rt * dt * (state.ydot_a[SEINT] + drhoedt); + amrex::Real rhoye_tmp = rhoye_old + 0.5_rt * dt * (drhoyedt_weak + drhoyedt_a); + + // compute the derivatives at the midpoint in time + + nse_derivs(rho_tmp, rhoe_tmp, + rhoye_tmp, T_old, + dt, mu_p, mu_n, + drhoyedt_weak, drhoyedt_a, + ydot_weak, state.ydot_a, + drhoedt, state.T_fixed); + + // evolve to the new time + + amrex::Real rho_new = rho_old + dt * state.ydot_a[SRHO]; + amrex::Real rhoe_new = rhoe_old + dt * (state.ydot_a[SEINT] + drhoedt); + amrex::Real rhoye_new = rhoye_old + dt * (drhoyedt_weak + drhoyedt_a); + + // Update the temperature with new internal energy + + amrex::Real T_new; + amrex::Real Ye_new = rhoye_new / rho_new; + + if (state.T_fixed > 0) { + T_new = state.T_fixed; + } else { + // initial eos guess + T_new = T_old; + amrex::Real e_new = rhoe_new / rho_new; + amrex::Real abar_new; + nse_T_abar_from_e(rho_new, e_new, Ye_new, T_new, abar_new, mu_p, mu_n); + } + + // With updated temp, rho, and Ye get final NSE state + + burn_t burn_state; + burn_state.T = T_new; + burn_state.rho = rho_new; + burn_state.y_e = Ye_new; + burn_state.mu_p = mu_p; + burn_state.mu_n = mu_n; + + auto nse_state = get_actual_nse_state(burn_state, 1.0e-10_rt, true); + + // store the new state // make sure they are normalized + amrex::Real sum_X{0.0_rt}; - amrex::Real X[NumSpec] = {0.0_rt}; - for (int n = 0; n < NumSpec; ++n) { - X[n] = amrex::max(small_x, amrex::min(1.0_rt, nse_state.xn[n])); - sum_X += X[n]; + for (auto & xn : nse_state.xn) { + xn = std::clamp(xn, small_x, 1.0_rt); + sum_X += xn; } - for (int n = 0; n < NumSpec; ++n) { - X[n] /= sum_X; + for (auto & xn : nse_state.xn) { + xn /= sum_X; } for (int n = 0; n < NumSpec; ++n) { - state.y[SFS+n] = state.y[SRHO] * X[n]; + state.y[SFS+n] = rho_new * nse_state.xn[n]; } - // density and momenta have already been updated - // update the total and internal energy now + // get the energy release from the change in rhoe over the timestep + // excluding any advection, and use that to update the total energy - state.y[SEINT] += dt * state.ydot_a[SEINT] + dt * rho_enucdot; + amrex::Real rho_enucdot = (rhoe_new - rhoe_old) / dt - state.ydot_a[SEINT]; + + state.y[SEINT] = rhoe_new; state.y[SEDEN] += dt * state.ydot_a[SEDEN] + dt * rho_enucdot; // store the chemical potentials + state.mu_p = nse_state.mu_p; state.mu_n = nse_state.mu_n; - } #endif diff --git a/networks/CNO_extras/actual_rhs.H b/networks/CNO_extras/actual_rhs.H index a9eaa3f3c0..b74804639b 100644 --- a/networks/CNO_extras/actual_rhs.H +++ b/networks/CNO_extras/actual_rhs.H @@ -693,11 +693,87 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(C13) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(N14) = 0.0_rt; + + ydot_nuc(N15) = 0.0_rt; + + ydot_nuc(O14) = 0.0_rt; + + ydot_nuc(O15) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(O17) = 0.0_rt; + + ydot_nuc(O18) = 0.0_rt; + + ydot_nuc(F17) = 0.0_rt; + + ydot_nuc(F18) = 0.0_rt; + + ydot_nuc(F19) = 0.0_rt; + + ydot_nuc(Ne18) = 0.0_rt; + + ydot_nuc(Ne19) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Mg22) = 0.0_rt; + + ydot_nuc(Mg24) = 0.0_rt; + + ydot_nuc(Fe56) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, @@ -858,7 +934,7 @@ void rhs_nuc(const burn_t& state, (screened_rates(k_He4_Ne20_to_Mg24)*Y(He4)*Y(Ne20)*state.rho + -screened_rates(k_Mg24_to_He4_Ne20)*Y(Mg24)) + (screened_rates(k_C12_O16_to_He4_Mg24)*Y(C12)*Y(O16)*state.rho + -screened_rates(k_He4_Mg24_to_C12_O16)*Y(He4)*Y(Mg24)*state.rho); - ydot_nuc(Fe56) = 0.0; + ydot_nuc(Fe56) = 0.0_rt; } diff --git a/networks/ECSN/actual_rhs.H b/networks/ECSN/actual_rhs.H index 251dd5610e..0b404497ef 100644 --- a/networks/ECSN/actual_rhs.H +++ b/networks/ECSN/actual_rhs.H @@ -275,7 +275,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; tabular_evaluate(j_F20_O20_meta, j_F20_O20_rhoy, j_F20_O20_temp, j_F20_O20_data, rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); @@ -312,6 +312,87 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + tabular_evaluate(j_F20_O20_meta, j_F20_O20_rhoy, j_F20_O20_temp, j_F20_O20_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_F20_to_O20) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Ne20_F20_meta, j_Ne20_F20_rhoy, j_Ne20_F20_temp, j_Ne20_F20_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Ne20_to_F20) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne20) * (edot_nu + edot_gamma); + + tabular_evaluate(j_O20_F20_meta, j_O20_F20_rhoy, j_O20_F20_temp, j_O20_F20_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_O20_to_F20) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(O20) * (edot_nu + edot_gamma); + + tabular_evaluate(j_F20_Ne20_meta, j_F20_Ne20_rhoy, j_F20_Ne20_temp, j_F20_Ne20_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_F20_to_Ne20) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(F20) * (edot_nu + edot_gamma); + + auto screened_rates = rate_eval.screened_rates; + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(O20) = + (-screened_rates(k_O20_to_F20)*Y(O20) + screened_rates(k_F20_to_O20)*Y(F20)); + + ydot_nuc(F20) = + (screened_rates(k_O20_to_F20)*Y(O20) + -screened_rates(k_F20_to_O20)*Y(F20)) + + (-screened_rates(k_F20_to_Ne20)*Y(F20) + screened_rates(k_Ne20_to_F20)*Y(Ne20)); + + ydot_nuc(Ne20) = + (screened_rates(k_F20_to_Ne20)*Y(F20) + -screened_rates(k_Ne20_to_F20)*Y(Ne20)); + + ydot_nuc(Mg24) = 0.0_rt; + + ydot_nuc(Al27) = 0.0_rt; + + ydot_nuc(Si28) = 0.0_rt; + + ydot_nuc(P31) = 0.0_rt; + + ydot_nuc(S32) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/He-C-Fe-group/actual_rhs.H b/networks/He-C-Fe-group/actual_rhs.H index 12c101e7ca..8736460ae3 100644 --- a/networks/He-C-Fe-group/actual_rhs.H +++ b/networks/He-C-Fe-group/actual_rhs.H @@ -1276,7 +1276,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; tabular_evaluate(j_Co55_Fe55_meta, j_Co55_Fe55_rhoy, j_Co55_Fe55_temp, j_Co55_Fe55_data, rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); @@ -1377,6 +1377,187 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + tabular_evaluate(j_Co55_Fe55_meta, j_Co55_Fe55_rhoy, j_Co55_Fe55_temp, j_Co55_Fe55_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Co55_to_Fe55) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Co55) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Co56_Fe56_meta, j_Co56_Fe56_rhoy, j_Co56_Fe56_temp, j_Co56_Fe56_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Co56_to_Fe56) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Co56) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Co56_Ni56_meta, j_Co56_Ni56_rhoy, j_Co56_Ni56_temp, j_Co56_Ni56_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Co56_to_Ni56) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Co56) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Co57_Ni57_meta, j_Co57_Ni57_rhoy, j_Co57_Ni57_temp, j_Co57_Ni57_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Co57_to_Ni57) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Co57) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Fe55_Co55_meta, j_Fe55_Co55_rhoy, j_Fe55_Co55_temp, j_Fe55_Co55_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Fe55_to_Co55) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe55) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Fe55_Mn55_meta, j_Fe55_Mn55_rhoy, j_Fe55_Mn55_temp, j_Fe55_Mn55_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Fe55_to_Mn55) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe55) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Fe56_Co56_meta, j_Fe56_Co56_rhoy, j_Fe56_Co56_temp, j_Fe56_Co56_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Fe56_to_Co56) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Fe56) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Mn55_Fe55_meta, j_Mn55_Fe55_rhoy, j_Mn55_Fe55_temp, j_Mn55_Fe55_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Mn55_to_Fe55) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Mn55) * (edot_nu + edot_gamma); + + tabular_evaluate(j_n_p_meta, j_n_p_rhoy, j_n_p_temp, j_n_p_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_n_to_p) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(N) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Ni56_Co56_meta, j_Ni56_Co56_rhoy, j_Ni56_Co56_temp, j_Ni56_Co56_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Ni56_to_Co56) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Ni56) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Ni57_Co57_meta, j_Ni57_Co57_rhoy, j_Ni57_Co57_temp, j_Ni57_Co57_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Ni57_to_Co57) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Ni57) * (edot_nu + edot_gamma); + + tabular_evaluate(j_p_n_meta, j_p_n_rhoy, j_p_n_temp, j_p_n_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_p_to_n) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma); + + auto screened_rates = rate_eval.screened_rates; + + ydot_nuc(N) = + -screened_rates(k_n_to_p)*Y(N) + + screened_rates(k_p_to_n)*Y(H1); + + ydot_nuc(H1) = + screened_rates(k_n_to_p)*Y(N) + + -screened_rates(k_p_to_n)*Y(H1); + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(N14) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(F18) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Ne21) = 0.0_rt; + + ydot_nuc(Na22) = 0.0_rt; + + ydot_nuc(Na23) = 0.0_rt; + + ydot_nuc(Mg24) = 0.0_rt; + + ydot_nuc(Al27) = 0.0_rt; + + ydot_nuc(Si28) = 0.0_rt; + + ydot_nuc(P31) = 0.0_rt; + + ydot_nuc(S32) = 0.0_rt; + + ydot_nuc(Ar36) = 0.0_rt; + + ydot_nuc(Ca40) = 0.0_rt; + + ydot_nuc(Ti44) = 0.0_rt; + + ydot_nuc(Cr48) = 0.0_rt; + + ydot_nuc(Mn51) = 0.0_rt; + + ydot_nuc(Mn55) = + (screened_rates(k_Fe55_to_Mn55)*Y(Fe55) + -screened_rates(k_Mn55_to_Fe55)*Y(Mn55)); + + ydot_nuc(Fe52) = 0.0_rt; + + ydot_nuc(Fe53) = 0.0_rt; + + ydot_nuc(Fe54) = 0.0_rt; + + ydot_nuc(Fe55) = + (screened_rates(k_Co55_to_Fe55)*Y(Co55) + -screened_rates(k_Fe55_to_Co55)*Y(Fe55)) + + (-screened_rates(k_Fe55_to_Mn55)*Y(Fe55) + screened_rates(k_Mn55_to_Fe55)*Y(Mn55)); + + ydot_nuc(Fe56) = + (screened_rates(k_Co56_to_Fe56)*Y(Co56) + -screened_rates(k_Fe56_to_Co56)*Y(Fe56)); + + ydot_nuc(Co55) = + (-screened_rates(k_Co55_to_Fe55)*Y(Co55) + screened_rates(k_Fe55_to_Co55)*Y(Fe55)); + + ydot_nuc(Co56) = + (-screened_rates(k_Co56_to_Fe56)*Y(Co56) + screened_rates(k_Fe56_to_Co56)*Y(Fe56)) + + (screened_rates(k_Ni56_to_Co56)*Y(Ni56) + -screened_rates(k_Co56_to_Ni56)*Y(Co56)); + + ydot_nuc(Co57) = + (screened_rates(k_Ni57_to_Co57)*Y(Ni57) + -screened_rates(k_Co57_to_Ni57)*Y(Co57)); + + ydot_nuc(Ni56) = + (-screened_rates(k_Ni56_to_Co56)*Y(Ni56) + screened_rates(k_Co56_to_Ni56)*Y(Co56)); + + ydot_nuc(Ni57) = + (-screened_rates(k_Ni57_to_Co57)*Y(Ni57) + screened_rates(k_Co57_to_Ni57)*Y(Co57)); + + ydot_nuc(Ni58) = 0.0_rt; + + ydot_nuc(Cu59) = 0.0_rt; + + ydot_nuc(Zn60) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/ase/actual_rhs.H b/networks/ase/actual_rhs.H index f5925ec203..4d08b90d65 100644 --- a/networks/ase/actual_rhs.H +++ b/networks/ase/actual_rhs.H @@ -894,11 +894,91 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(N) = 0.0_rt; + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(N14) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(F18) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Ne21) = 0.0_rt; + + ydot_nuc(Na22) = 0.0_rt; + + ydot_nuc(Na23) = 0.0_rt; + + ydot_nuc(Mg24) = 0.0_rt; + + ydot_nuc(Al27) = 0.0_rt; + + ydot_nuc(Si28) = 0.0_rt; + + ydot_nuc(P31) = 0.0_rt; + + ydot_nuc(S32) = 0.0_rt; + + ydot_nuc(Ar36) = 0.0_rt; + + ydot_nuc(Ca40) = 0.0_rt; + + ydot_nuc(Ti44) = 0.0_rt; + + ydot_nuc(Cr48) = 0.0_rt; + + ydot_nuc(Fe52) = 0.0_rt; + + ydot_nuc(Ni56) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H index ad623a6c06..d96207d283 100644 --- a/networks/ignition_reaclib/C-burn-simple/actual_rhs.H +++ b/networks/ignition_reaclib/C-burn-simple/actual_rhs.H @@ -140,11 +140,61 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(N) = 0.0_rt; + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Na23) = 0.0_rt; + + ydot_nuc(Mg23) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/ignition_reaclib/URCA-medium/actual_rhs.H b/networks/ignition_reaclib/URCA-medium/actual_rhs.H index 1b9f728b47..a6d35692e1 100644 --- a/networks/ignition_reaclib/URCA-medium/actual_rhs.H +++ b/networks/ignition_reaclib/URCA-medium/actual_rhs.H @@ -300,7 +300,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; tabular_evaluate(j_Na23_Ne23_meta, j_Na23_Ne23_rhoy, j_Na23_Ne23_temp, j_Na23_Ne23_data, rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); @@ -345,6 +345,96 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + tabular_evaluate(j_Na23_Ne23_meta, j_Na23_Ne23_rhoy, j_Na23_Ne23_temp, j_Na23_Ne23_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Na23_to_Ne23) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Na23) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Ne23_Na23_meta, j_Ne23_Na23_rhoy, j_Ne23_Na23_temp, j_Ne23_Na23_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Ne23_to_Na23) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne23) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Mg23_Na23_meta, j_Mg23_Na23_rhoy, j_Mg23_Na23_temp, j_Mg23_Na23_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Mg23_to_Na23) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Mg23) * (edot_nu + edot_gamma); + + tabular_evaluate(j_n_p_meta, j_n_p_rhoy, j_n_p_temp, j_n_p_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_n_to_p) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(N) * (edot_nu + edot_gamma); + + tabular_evaluate(j_p_n_meta, j_p_n_rhoy, j_p_n_temp, j_p_n_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_p_to_n) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma); + + auto screened_rates = rate_eval.screened_rates; + + ydot_nuc(N) = + -screened_rates(k_n_to_p)*Y(N) + + screened_rates(k_p_to_n)*Y(H1); + + ydot_nuc(H1) = + screened_rates(k_n_to_p)*Y(N) + + -screened_rates(k_p_to_n)*Y(H1); + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Ne23) = + (-screened_rates(k_Ne23_to_Na23)*Y(Ne23) + screened_rates(k_Na23_to_Ne23)*Y(Na23)); + + ydot_nuc(Na23) = + (screened_rates(k_Ne23_to_Na23)*Y(Ne23) + -screened_rates(k_Na23_to_Ne23)*Y(Na23)) + + screened_rates(k_Mg23_to_Na23)*Y(Mg23); + + ydot_nuc(Mg23) = + -screened_rates(k_Mg23_to_Na23)*Y(Mg23); + + ydot_nuc(Mg24) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/ignition_reaclib/URCA-simple/actual_rhs.H b/networks/ignition_reaclib/URCA-simple/actual_rhs.H index 627b2a97f8..ff15e513fb 100644 --- a/networks/ignition_reaclib/URCA-simple/actual_rhs.H +++ b/networks/ignition_reaclib/URCA-simple/actual_rhs.H @@ -140,7 +140,7 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; tabular_evaluate(j_Na23_Ne23_meta, j_Na23_Ne23_rhoy, j_Na23_Ne23_temp, j_Na23_Ne23_data, rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); @@ -177,6 +177,85 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + tabular_evaluate(j_Na23_Ne23_meta, j_Na23_Ne23_rhoy, j_Na23_Ne23_temp, j_Na23_Ne23_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Na23_to_Ne23) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Na23) * (edot_nu + edot_gamma); + + tabular_evaluate(j_Ne23_Na23_meta, j_Ne23_Na23_rhoy, j_Ne23_Na23_temp, j_Ne23_Na23_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_Ne23_to_Na23) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(Ne23) * (edot_nu + edot_gamma); + + tabular_evaluate(j_n_p_meta, j_n_p_rhoy, j_n_p_temp, j_n_p_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_n_to_p) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(N) * (edot_nu + edot_gamma); + + tabular_evaluate(j_p_n_meta, j_p_n_rhoy, j_p_n_temp, j_p_n_data, + rhoy, state.T, rate, drate_dt, edot_nu, edot_gamma); + rate_eval.screened_rates(k_p_to_n) = rate; + rate_eval.enuc_weak += C::Legacy::n_A * Y(H1) * (edot_nu + edot_gamma); + + auto screened_rates = rate_eval.screened_rates; + + ydot_nuc(N) = + -screened_rates(k_n_to_p)*Y(N) + + screened_rates(k_p_to_n)*Y(H1); + + ydot_nuc(H1) = + screened_rates(k_n_to_p)*Y(N) + + -screened_rates(k_p_to_n)*Y(H1); + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Ne23) = + (-screened_rates(k_Ne23_to_Na23)*Y(Ne23) + screened_rates(k_Na23_to_Ne23)*Y(Na23)); + + ydot_nuc(Na23) = + (screened_rates(k_Ne23_to_Na23)*Y(Ne23) + -screened_rates(k_Na23_to_Ne23)*Y(Na23)); + + ydot_nuc(Mg23) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/nova/actual_rhs.H b/networks/nova/actual_rhs.H index 0c2dd612d3..8f3e7c1c60 100644 --- a/networks/nova/actual_rhs.H +++ b/networks/nova/actual_rhs.H @@ -362,11 +362,71 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(C13) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(N14) = 0.0_rt; + + ydot_nuc(N15) = 0.0_rt; + + ydot_nuc(O14) = 0.0_rt; + + ydot_nuc(O15) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(O17) = 0.0_rt; + + ydot_nuc(F17) = 0.0_rt; + + ydot_nuc(F18) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/nova2/actual_rhs.H b/networks/nova2/actual_rhs.H index 1219595277..6cb140fd89 100644 --- a/networks/nova2/actual_rhs.H +++ b/networks/nova2/actual_rhs.H @@ -559,11 +559,79 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(H2) = 0.0_rt; + + ydot_nuc(He3) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(Be7) = 0.0_rt; + + ydot_nuc(B8) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(C13) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(N14) = 0.0_rt; + + ydot_nuc(N15) = 0.0_rt; + + ydot_nuc(O14) = 0.0_rt; + + ydot_nuc(O15) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(O17) = 0.0_rt; + + ydot_nuc(F17) = 0.0_rt; + + ydot_nuc(F18) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/partition_test/actual_rhs.H b/networks/partition_test/actual_rhs.H index ba170ca871..72296fab0e 100644 --- a/networks/partition_test/actual_rhs.H +++ b/networks/partition_test/actual_rhs.H @@ -140,11 +140,55 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(Fe52) = 0.0_rt; + + ydot_nuc(Co55) = 0.0_rt; + + ydot_nuc(Ni56) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/sn160/actual_rhs.H b/networks/sn160/actual_rhs.H index 9d14330199..3ff15b8115 100644 --- a/networks/sn160/actual_rhs.H +++ b/networks/sn160/actual_rhs.H @@ -8815,11 +8815,365 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(N) = 0.0_rt; + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(H2) = 0.0_rt; + + ydot_nuc(He3) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(Li6) = 0.0_rt; + + ydot_nuc(Li7) = 0.0_rt; + + ydot_nuc(Be7) = 0.0_rt; + + ydot_nuc(Be9) = 0.0_rt; + + ydot_nuc(B8) = 0.0_rt; + + ydot_nuc(B10) = 0.0_rt; + + ydot_nuc(B11) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(C13) = 0.0_rt; + + ydot_nuc(C14) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(N14) = 0.0_rt; + + ydot_nuc(N15) = 0.0_rt; + + ydot_nuc(O14) = 0.0_rt; + + ydot_nuc(O15) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(O17) = 0.0_rt; + + ydot_nuc(O18) = 0.0_rt; + + ydot_nuc(F17) = 0.0_rt; + + ydot_nuc(F18) = 0.0_rt; + + ydot_nuc(F19) = 0.0_rt; + + ydot_nuc(Ne18) = 0.0_rt; + + ydot_nuc(Ne19) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Ne21) = 0.0_rt; + + ydot_nuc(Ne22) = 0.0_rt; + + ydot_nuc(Na21) = 0.0_rt; + + ydot_nuc(Na22) = 0.0_rt; + + ydot_nuc(Na23) = 0.0_rt; + + ydot_nuc(Mg23) = 0.0_rt; + + ydot_nuc(Mg24) = 0.0_rt; + + ydot_nuc(Mg25) = 0.0_rt; + + ydot_nuc(Mg26) = 0.0_rt; + + ydot_nuc(Al25) = 0.0_rt; + + ydot_nuc(Al26) = 0.0_rt; + + ydot_nuc(Al27) = 0.0_rt; + + ydot_nuc(Si28) = 0.0_rt; + + ydot_nuc(Si29) = 0.0_rt; + + ydot_nuc(Si30) = 0.0_rt; + + ydot_nuc(Si31) = 0.0_rt; + + ydot_nuc(Si32) = 0.0_rt; + + ydot_nuc(P29) = 0.0_rt; + + ydot_nuc(P30) = 0.0_rt; + + ydot_nuc(P31) = 0.0_rt; + + ydot_nuc(P32) = 0.0_rt; + + ydot_nuc(P33) = 0.0_rt; + + ydot_nuc(S32) = 0.0_rt; + + ydot_nuc(S33) = 0.0_rt; + + ydot_nuc(S34) = 0.0_rt; + + ydot_nuc(S35) = 0.0_rt; + + ydot_nuc(S36) = 0.0_rt; + + ydot_nuc(Cl33) = 0.0_rt; + + ydot_nuc(Cl34) = 0.0_rt; + + ydot_nuc(Cl35) = 0.0_rt; + + ydot_nuc(Cl36) = 0.0_rt; + + ydot_nuc(Cl37) = 0.0_rt; + + ydot_nuc(Ar36) = 0.0_rt; + + ydot_nuc(Ar37) = 0.0_rt; + + ydot_nuc(Ar38) = 0.0_rt; + + ydot_nuc(Ar39) = 0.0_rt; + + ydot_nuc(Ar40) = 0.0_rt; + + ydot_nuc(K37) = 0.0_rt; + + ydot_nuc(K38) = 0.0_rt; + + ydot_nuc(K39) = 0.0_rt; + + ydot_nuc(K40) = 0.0_rt; + + ydot_nuc(K41) = 0.0_rt; + + ydot_nuc(Ca40) = 0.0_rt; + + ydot_nuc(Ca41) = 0.0_rt; + + ydot_nuc(Ca42) = 0.0_rt; + + ydot_nuc(Ca43) = 0.0_rt; + + ydot_nuc(Ca44) = 0.0_rt; + + ydot_nuc(Ca45) = 0.0_rt; + + ydot_nuc(Ca46) = 0.0_rt; + + ydot_nuc(Ca47) = 0.0_rt; + + ydot_nuc(Ca48) = 0.0_rt; + + ydot_nuc(Sc43) = 0.0_rt; + + ydot_nuc(Sc44) = 0.0_rt; + + ydot_nuc(Sc45) = 0.0_rt; + + ydot_nuc(Sc46) = 0.0_rt; + + ydot_nuc(Sc47) = 0.0_rt; + + ydot_nuc(Sc48) = 0.0_rt; + + ydot_nuc(Sc49) = 0.0_rt; + + ydot_nuc(Ti44) = 0.0_rt; + + ydot_nuc(Ti45) = 0.0_rt; + + ydot_nuc(Ti46) = 0.0_rt; + + ydot_nuc(Ti47) = 0.0_rt; + + ydot_nuc(Ti48) = 0.0_rt; + + ydot_nuc(Ti49) = 0.0_rt; + + ydot_nuc(Ti50) = 0.0_rt; + + ydot_nuc(Ti51) = 0.0_rt; + + ydot_nuc(V46) = 0.0_rt; + + ydot_nuc(V47) = 0.0_rt; + + ydot_nuc(V48) = 0.0_rt; + + ydot_nuc(V49) = 0.0_rt; + + ydot_nuc(V50) = 0.0_rt; + + ydot_nuc(V51) = 0.0_rt; + + ydot_nuc(V52) = 0.0_rt; + + ydot_nuc(Cr48) = 0.0_rt; + + ydot_nuc(Cr49) = 0.0_rt; + + ydot_nuc(Cr50) = 0.0_rt; + + ydot_nuc(Cr51) = 0.0_rt; + + ydot_nuc(Cr52) = 0.0_rt; + + ydot_nuc(Cr53) = 0.0_rt; + + ydot_nuc(Cr54) = 0.0_rt; + + ydot_nuc(Mn50) = 0.0_rt; + + ydot_nuc(Mn51) = 0.0_rt; + + ydot_nuc(Mn52) = 0.0_rt; + + ydot_nuc(Mn53) = 0.0_rt; + + ydot_nuc(Mn54) = 0.0_rt; + + ydot_nuc(Mn55) = 0.0_rt; + + ydot_nuc(Fe52) = 0.0_rt; + + ydot_nuc(Fe53) = 0.0_rt; + + ydot_nuc(Fe54) = 0.0_rt; + + ydot_nuc(Fe55) = 0.0_rt; + + ydot_nuc(Fe56) = 0.0_rt; + + ydot_nuc(Fe57) = 0.0_rt; + + ydot_nuc(Fe58) = 0.0_rt; + + ydot_nuc(Co53) = 0.0_rt; + + ydot_nuc(Co54) = 0.0_rt; + + ydot_nuc(Co55) = 0.0_rt; + + ydot_nuc(Co56) = 0.0_rt; + + ydot_nuc(Co57) = 0.0_rt; + + ydot_nuc(Co58) = 0.0_rt; + + ydot_nuc(Co59) = 0.0_rt; + + ydot_nuc(Ni56) = 0.0_rt; + + ydot_nuc(Ni57) = 0.0_rt; + + ydot_nuc(Ni58) = 0.0_rt; + + ydot_nuc(Ni59) = 0.0_rt; + + ydot_nuc(Ni60) = 0.0_rt; + + ydot_nuc(Ni61) = 0.0_rt; + + ydot_nuc(Ni62) = 0.0_rt; + + ydot_nuc(Ni63) = 0.0_rt; + + ydot_nuc(Ni64) = 0.0_rt; + + ydot_nuc(Cu57) = 0.0_rt; + + ydot_nuc(Cu58) = 0.0_rt; + + ydot_nuc(Cu59) = 0.0_rt; + + ydot_nuc(Cu60) = 0.0_rt; + + ydot_nuc(Cu61) = 0.0_rt; + + ydot_nuc(Cu62) = 0.0_rt; + + ydot_nuc(Cu63) = 0.0_rt; + + ydot_nuc(Cu64) = 0.0_rt; + + ydot_nuc(Cu65) = 0.0_rt; + + ydot_nuc(Zn59) = 0.0_rt; + + ydot_nuc(Zn60) = 0.0_rt; + + ydot_nuc(Zn61) = 0.0_rt; + + ydot_nuc(Zn62) = 0.0_rt; + + ydot_nuc(Zn63) = 0.0_rt; + + ydot_nuc(Zn64) = 0.0_rt; + + ydot_nuc(Zn65) = 0.0_rt; + + ydot_nuc(Zn66) = 0.0_rt; + + ydot_nuc(Ga62) = 0.0_rt; + + ydot_nuc(Ga63) = 0.0_rt; + + ydot_nuc(Ga64) = 0.0_rt; + + ydot_nuc(Ge63) = 0.0_rt; + + ydot_nuc(Ge64) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/subch_base/actual_rhs.H b/networks/subch_base/actual_rhs.H index 6b083c0284..c1fae4ff7e 100644 --- a/networks/subch_base/actual_rhs.H +++ b/networks/subch_base/actual_rhs.H @@ -804,11 +804,81 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Na23) = 0.0_rt; + + ydot_nuc(Mg24) = 0.0_rt; + + ydot_nuc(Al27) = 0.0_rt; + + ydot_nuc(Si28) = 0.0_rt; + + ydot_nuc(P31) = 0.0_rt; + + ydot_nuc(S32) = 0.0_rt; + + ydot_nuc(Ar36) = 0.0_rt; + + ydot_nuc(Ca40) = 0.0_rt; + + ydot_nuc(Ti44) = 0.0_rt; + + ydot_nuc(Cr48) = 0.0_rt; + + ydot_nuc(Fe52) = 0.0_rt; + + ydot_nuc(Ni56) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/networks/subch_simple/actual_rhs.H b/networks/subch_simple/actual_rhs.H index 5a72edee3c..9e86ce9ddf 100644 --- a/networks/subch_simple/actual_rhs.H +++ b/networks/subch_simple/actual_rhs.H @@ -1038,11 +1038,89 @@ void evaluate_rates(const burn_t& state, T& rate_eval) { [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; - rate_eval.enuc_weak = 0.0; + rate_eval.enuc_weak = 0.0_rt; } +#ifdef NSE_NET +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void get_ydot_weak(const burn_t& state, + amrex::Array1D& ydot_nuc, + amrex::Real& enuc_weak, + [[maybe_unused]] const amrex::Array1D& Y) { + /// + /// Calculate Ydots contribute only from weak reactions. + /// This is used to calculate dyedt and energy generation from + /// weak reactions for self-consistent NSE + /// + + + // initialize ydot_nuc to 0 + + for (int i = 1; i <= neqs; ++i) { + ydot_nuc(i) = 0.0_rt; + } + + rate_t rate_eval; + + [[maybe_unused]] amrex::Real rate, drate_dt, edot_nu, edot_gamma; + [[maybe_unused]] amrex::Real rhoy = state.rho * state.y_e; + + rate_eval.enuc_weak = 0.0_rt; + + // Calculate tabular rates and get ydot_weak + + + ydot_nuc(H1) = 0.0_rt; + + ydot_nuc(He4) = 0.0_rt; + + ydot_nuc(C12) = 0.0_rt; + + ydot_nuc(N13) = 0.0_rt; + + ydot_nuc(N14) = 0.0_rt; + + ydot_nuc(O16) = 0.0_rt; + + ydot_nuc(F18) = 0.0_rt; + + ydot_nuc(Ne20) = 0.0_rt; + + ydot_nuc(Ne21) = 0.0_rt; + + ydot_nuc(Na22) = 0.0_rt; + + ydot_nuc(Na23) = 0.0_rt; + + ydot_nuc(Mg24) = 0.0_rt; + + ydot_nuc(Al27) = 0.0_rt; + + ydot_nuc(Si28) = 0.0_rt; + + ydot_nuc(P31) = 0.0_rt; + + ydot_nuc(S32) = 0.0_rt; + + ydot_nuc(Ar36) = 0.0_rt; + + ydot_nuc(Ca40) = 0.0_rt; + + ydot_nuc(Ti44) = 0.0_rt; + + ydot_nuc(Cr48) = 0.0_rt; + + ydot_nuc(Fe52) = 0.0_rt; + + ydot_nuc(Ni56) = 0.0_rt; + + enuc_weak = rate_eval.enuc_weak; +} +#endif + + AMREX_GPU_HOST_DEVICE AMREX_INLINE void rhs_nuc(const burn_t& state, amrex::Array1D& ydot_nuc, diff --git a/nse_solver/Make.package b/nse_solver/Make.package index 66fe590429..f5344d27ff 100644 --- a/nse_solver/Make.package +++ b/nse_solver/Make.package @@ -1,4 +1,5 @@ ifeq ($(USE_NSE_NET), TRUE) CEXE_headers += nse_solver.H CEXE_headers += nse_check.H + CEXE_headers += nse_eos.H endif diff --git a/nse_solver/nse_eos.H b/nse_solver/nse_eos.H new file mode 100644 index 0000000000..69d4fccc44 --- /dev/null +++ b/nse_solver/nse_eos.H @@ -0,0 +1,174 @@ +#ifndef NSE_EOS_H +#define NSE_EOS_H + +#include + +#include + +#include +#include + + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real nse_abar(const amrex::Real T, const amrex::Real rho, + const amrex::Real Ye, amrex::Real &mu_p, + amrex::Real &mu_n) { + /// + /// This function calculates abar from NSE using + /// Temp, rho, and Ye + /// + + burn_t burn_state; + burn_state.rho = rho; + burn_state.y_e = Ye; + burn_state.T = T; + burn_state.mu_p = mu_p; + burn_state.mu_n = mu_n; + + auto nse_state = get_actual_nse_state(burn_state, 1.0e-10_rt, true); + + amrex::Real abar{0.0_rt}; + for (int n = 0; n < NumSpec; ++n) { + abar += nse_state.xn[n] * aion_inv[n]; + } + + // update mu_p and mu_n + + mu_p = burn_state.mu_p; + mu_n = burn_state.mu_n; + + return abar; +} + + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::Real nse_dabar_dT(const amrex::Real T, const amrex::Real rho, + const amrex::Real Ye, amrex::Real &mu_p, + amrex::Real &mu_n) { + /// + /// This function constructs dabar_dT + /// This should be 2nd order accurate + /// + + // deviation in temperature + + const amrex::Real dT = 1.0e-6_rt * T; + + // Calculate derivative using five-point stencil method + + amrex::Real dabar_dT = (-nse_abar(T + 2.0_rt*dT, rho, Ye, mu_p, mu_n) + + 8.0_rt * nse_abar(T+dT, rho, Ye, mu_p, mu_n) - + 8.0_rt * nse_abar(T-dT, rho, Ye, mu_p, mu_n) + + nse_abar(T - 2.0_rt*dT, rho, Ye, mu_p, mu_n) + ) / (12.0_rt * dT); + + // Calculate derivative using central differencing + + // amrex::Real dabar_dT = 0.5_rt * (nse_abar(T + dT, rho, Ye, mu_p, mu_n) - + // nse_abar(T - dT, rho, Ye, mu_p, mu_n)) / dT; + + return dabar_dT; +} + +/// +/// This function inverts this form of the EOS to find the T +/// that satisfies the EOS and NSE given an input e and rho. +/// +/// if we are in NSE, then the entire thermodynamic state is just +/// a function of rho, T, Ye. We can write the energy as: +/// +/// e = e(rho, T, Y_e, Abar(rho, T, Ye)) +/// +/// where we note that Abar is a function of those same inputs. +/// +/// The basic idea is that Abar and Zbar are both functions of +/// rho, T, Ye through NSE calculations, 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 +void +nse_T_abar_from_e(const amrex::Real rho, const amrex::Real e_in, + const amrex::Real Ye, + amrex::Real &T, amrex::Real &abar, + amrex::Real &mu_p, amrex::Real &mu_n) { + + constexpr Real ttol{1.e-8_rt}; + constexpr int max_iter{100}; + + bool converged{false}; + + int iter{0}; + + // initialize burn_state + burn_t burn_state; + burn_state.rho = rho; + burn_state.y_e = Ye; + burn_state.mu_p = mu_p; + burn_state.mu_n = mu_n; + + while (not converged && iter < max_iter) { + + // update Temperature + + burn_state.T = T; + + auto nse_state = get_actual_nse_state(burn_state, 1.0e-10_rt, true); + + // call the EOS with the initial guess for T + // as well as density and NSE mass fractions + + eos_extra_t eos_state; + for (int n = 0; n < NumSpec; ++n) { + eos_state.xn[n] = nse_state.xn[n]; + } + + // Call the EOS to get internal energy + // abar, zbar, etc are calculated automatically in EOS + + eos_state.rho = rho; + eos_state.T = T; + eos(eos_input_rt, eos_state); + + // f is the quantity we want to zero + + amrex::Real f = eos_state.e - e_in; + + amrex::Real dabar_dT = nse_dabar_dT(T, rho, Ye, + nse_state.mu_p, nse_state.mu_n); + + // compute the correction to our guess + + Real dT = -f / (eos_state.dedT + eos_state.dedA * dabar_dT + + Ye * eos_state.dedZ * dabar_dT); + + // update the temperature and abar + + T = std::clamp(T + dT, 0.25 * T, 4.0 * T); + abar = eos_state.abar; + + // update mu_p and mu_n + + mu_p = nse_state.mu_p; + mu_n = nse_state.mu_n; + + // check convergence + + if (std::abs(dT) < ttol * T) { + converged = true; + } + iter++; + + } + +} + +#endif