From 92c5accae807788b49ee1b60bd77c17bc66fead9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 18 Oct 2023 12:49:10 -0400 Subject: [PATCH 01/13] add neutrino loses to NSE table with SDC --- integration/nse_update_simplified_sdc.H | 26 +++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 533eb5171b..26c74beee0 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -58,6 +58,21 @@ void sdc_nse_burn(BurnT& state, const Real dt) { nse_interp(T_in, state.rho, ye_in, abar_out, dq_out, dyedt, X); + // get the neutrino loss term + Real snu{0.0}; + Real dsnudt{0.0}; + Real dsnudd{0.0}; + Real dsnuda{0.0}; + Real dsnudz{0.0}; + + Real zbar = abar_out * ye_in; + + constexpr int do_derivatives = 0; + sneut5(T_in, state.rho, abar_out, zbar, + snu, dsnudt, dsnudd, dsnuda, dsnudz); + + Real snu_old = snu; + Real dyedt_old = dyedt; // density and momentum have no reactive sources @@ -82,7 +97,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real rho_aux_new[NumAux]; Real rhoe_new; - Real rho_enucdot = 0.0_rt; + Real rho_enucdot = -state.rho * snu; Real rho_half = 0.5_rt * (rho_old + state.y[SRHO]); @@ -130,7 +145,14 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real rho_dBEA = rho_bea_tilde - rho_bea_old; // this is MeV / nucleon * g / cm**3 // convert the energy to erg / cm**3 - rho_enucdot = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; + rho_enucdot = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; + + // now get the updated neutrino term + zbar = abar_out * eos_state.aux[iye]; + sneut5(T_new, eos_state.rho, abar_out, zbar, + snu, dsnudt, dsnudd, dsnuda, dsnudz); + + rho_enucdot -= 0.5_rt * rho_half * (snu_old + snu); // update the new state for the next pass From fbcef5fc91d3fa380394ea394f92490dfb49decc Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 22 Oct 2023 14:50:02 -0400 Subject: [PATCH 02/13] a lot of restructuring --- integration/nse_update_simplified_sdc.H | 131 ++++++++++++------------ 1 file changed, 64 insertions(+), 67 deletions(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 26c74beee0..a6b1c285ad 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -40,49 +40,59 @@ void sdc_nse_burn(BurnT& state, const Real dt) { state.n_rhs = 0; state.n_jac = 0; - // call the NSE table to get (dYe/dt)^n - Real abar_out; - Real dq_out; - Real dyedt; - Real X[NumSpec]; - + // we need the initial Ye to get the NSE state. We also + // want the initial rho since we may not be entering + // this routine already in NSE. This means that there will + // be an energy release from us instantaneously adjusting + // into NSE (the first call to nse_interp) + an energy + // release from the evolution of NSE over the timestep Real ye_in = state.y[SFX+iye] / state.rho; + Real rho_bea_old = state.y[SFX+ibea]; + + // density and momentum have no reactive sources + Real rho_old = state.y[SRHO]; + + 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]; // if we are doing drive_initial_convection, we want to use // the temperature that comes in through T_fixed Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T; - // get the current NSE state from the table - - nse_interp(T_in, state.rho, ye_in, - abar_out, dq_out, dyedt, X); + // get the neutrino loss term -- we want to use the state that we + // came in here with, so the original Abar and Zbar - // get the neutrino loss term Real snu{0.0}; Real dsnudt{0.0}; Real dsnudd{0.0}; Real dsnuda{0.0}; Real dsnudz{0.0}; - Real zbar = abar_out * ye_in; + Real abar = state.y[SFX+iabar] / rho_old; + Real zbar = abar * ye_in; constexpr int do_derivatives = 0; - sneut5(T_in, state.rho, abar_out, zbar, + sneut5(T_in, state.rho, abar, zbar, snu, dsnudt, dsnudd, dsnuda, dsnudz); Real snu_old = snu; - Real dyedt_old = dyedt; + // get the current NSE state from the table -- this will be used + // to compute the NSE evolution sources - // density and momentum have no reactive sources - Real rho_old = state.y[SRHO]; - Real rho_bea_old = state.y[SFX+ibea]; + Real abar_out; + Real dq_out; + Real dyedt; + Real X[NumSpec]; + + nse_interp(T_in, state.rho, ye_in, + abar_out, dq_out, dyedt, X); + + Real dyedt_old = dyedt; - 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]; // predict the U^{n+1,*} state with only estimates of the aux // reaction terms and advection to dt @@ -211,53 +221,42 @@ void sdc_nse_burn(BurnT& state, const Real dt) { state.n_rhs = 0; state.n_jac = 0; - // call the NSE table to get (dYe/dt)^n - Real abar_out; - Real dq_out; - Real dyedt; - Real X[NumSpec]; + // store the initial mass fractions -- we will need these + // to compute the energy release. - // TODO: are state.y[SFS:] and state.xn[:] synced? + Real X_old[NumSpec]; - // if we are doing drive_initial_convection, we want to use - // the temperature that comes in through T_fixed - - Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T; - - // We will get initial NSE prediction using the input X's - BurnT nse_state_in{state}; - nse_state_in.T = T_in; - - // solve for the NSE state directly -- we'll have it compute - // ye from the input Xs - auto nse_state = get_actual_nse_state(nse_state_in); - - // compute the new binding energy (B/A) and dyedt - dq_out = 0.0; for (int n = 0; n < NumSpec; ++n) { - dq_out += nse_state.xn[n] * network::bion(n+1) * aion_inv[n]; + X_old[n] = state.y[SFS+n] / state.y[SRHO]; } - dyedt = 0.0; // we can update this in the future by calling actual_rhs() - - Real dyedt_old = dyedt; - // density and momentum have no reactive sources Real rho_old = state.y[SRHO]; - // compute the original rho (B/A) of the composition - Real rho_bea_old = 0.0; - for (int n = 0; n < NumSpec; ++n) { - rho_bea_old += state.y[SFS+n] * network::bion(n+1) * aion_inv[n]; - } - 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]; + // if we are doing drive_initial_convection, we want to use + // the temperature that comes in through T_fixed + + Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T; + + // get the neutrino loss term -- we want to use the state that we + // came in here with, so the original Abar and Zbar + + + // 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. + + Real dyedt_old = 0.0; + + // predict the U^{n+1,*} state with only estimates of the X - // updates with advection to dt + // updates with advection to dt and the neutrino loss term in + // energy BurnT burn_state; burn_state.T = T_in; // initial guess @@ -309,27 +308,25 @@ void sdc_nse_burn(BurnT& state, const Real dt) { nse_state = get_actual_nse_state(burn_state); - // compute (B/A) and dyedt - dq_out = 0.0; + // compute the energy release. The mass fractions in nse_state.xn[] + // include the advective parts, so first we need to remove that. + + Real rhoX_tilde[NumSpec]; for (int n = 0; n < NumSpec; ++n) { - dq_out += nse_state.xn[n] * network::bion(n+1) * aion_inv[n]; + rhoX_tilde[n] = state.y[SRHO] * nse_state.xn[n] - dt * state.ydot_a[SFS+n]; } dyedt = 0.0_rt; // we can update this in the future by calling actual_rhs() - // compute the energy release -- we need to remove the - // advective part. To do this, we must first compute the - // advective update for B/A from those of the mass fractions - Real ydot_a_BA = 0.0_rt; + // 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) { - ydot_a_BA += network::bion(n+1) * aion_inv[n] * state.ydot_a[SFS+n]; + rho_enucdot += (rhoX_tilde[n] - rho_old * X_old[n]) * + network::mion(n+1) * aion_inv(n+1); } + rho_enucdot *= C::Legacy::enuc_conv2; - Real rho_bea_tilde = state.y[SRHO] * dq_out - dt * ydot_a_BA; - Real rho_dBEA = rho_bea_tilde - rho_bea_old; // this is MeV / nucleon * g / cm**3 - - // convert the energy to erg / cm**3 - rho_enucdot = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; + // now get the updated neutrino term // update the new state for the next pass -- this should // already implicitly have the advective portion included, From 6acf85b9b00d68b07cae4ef685fababd0091e6b0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 22 Oct 2023 15:07:53 -0400 Subject: [PATCH 03/13] fix compilation --- integration/nse_update_simplified_sdc.H | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index a6b1c285ad..0a424d5dba 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -279,6 +279,8 @@ void sdc_nse_burn(BurnT& state, const Real dt) { rhoX_new[n] = state.y[SFS+n] + dt * state.ydot_a[SFS+n] + dt * rhoX_source[n]; } + burn_t nse_state; + for (int iter = 0; iter < 3; iter++) { // update (rho e)^{n+1} based on the new energy generation rate @@ -316,18 +318,19 @@ void sdc_nse_burn(BurnT& state, const Real dt) { rhoX_tilde[n] = state.y[SRHO] * nse_state.xn[n] - dt * state.ydot_a[SFS+n]; } - dyedt = 0.0_rt; // we can update this in the future by calling actual_rhs() + Real dyedt = 0.0_rt; // we can update this in the future by calling actual_rhs() // 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+1); + network::mion(n+1) * aion_inv[n]; } rho_enucdot *= C::Legacy::enuc_conv2; // now get the updated neutrino term + // 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 @@ -343,6 +346,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // the new mass fractions are just those that come from the table // make sure they are normalized Real sum_X{0.0_rt}; + 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]; From 1e136b3c3512ac6392e7c5e00cc95320c71e9645 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 22 Oct 2023 15:21:07 -0400 Subject: [PATCH 04/13] update --- integration/nse_update_simplified_sdc.H | 41 ++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 4 deletions(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 0a424d5dba..adcf335f6f 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -75,7 +75,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real zbar = abar * ye_in; constexpr int do_derivatives = 0; - sneut5(T_in, state.rho, abar, zbar, + sneut5(T_in, rho_old, abar, zbar, snu, dsnudt, dsnudd, dsnuda, dsnudz); Real snu_old = snu; @@ -88,7 +88,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real dyedt; Real X[NumSpec]; - nse_interp(T_in, state.rho, ye_in, + nse_interp(T_in, rho_old, ye_in, abar_out, dq_out, dyedt, X); Real dyedt_old = dyedt; @@ -107,7 +107,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real rho_aux_new[NumAux]; Real rhoe_new; - Real rho_enucdot = -state.rho * snu; + Real rho_enucdot = -rho_old * snu; Real rho_half = 0.5_rt * (rho_old + state.y[SRHO]); @@ -245,6 +245,26 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // get the neutrino loss term -- we want to use the state that we // came in here with, so the original Abar and Zbar + Real snu{0.0}; + Real dsnudt{0.0}; + Real dsnudd{0.0}; + Real dsnuda{0.0}; + Real dsnudz{0.0}; + + Real abar{0.0}; + 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]; + } + abar = 1.0 / abar; + zbar *= abar; + + constexpr int do_derivatives = 0; + sneut5(T_in, rho_old, abar, zbar, + snu, dsnudt, dsnudd, dsnuda, dsnudz); + + Real snu_old = snu; // if our network could return the evolution of Ye due to the @@ -268,7 +288,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real rhoX_new[NumSpec]; Real rhoe_new; - Real rho_enucdot = 0.0_rt; + Real rho_enucdot = -rho_old * snu; Real rho_half = 0.5_rt * (rho_old + state.y[SRHO]); @@ -329,7 +349,20 @@ void sdc_nse_burn(BurnT& state, const Real dt) { rho_enucdot *= C::Legacy::enuc_conv2; // 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; + + constexpr int do_derivatives = 0; + sneut5(T_new, state.y[SRHO], abar, zbar, + snu, dsnudt, dsnudd, dsnuda, dsnudz); + 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, From 55bafe272a7605fe4cb7a8cbb45ace0818ca561a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 15 Nov 2023 17:31:38 -0500 Subject: [PATCH 05/13] start of a refactoring of SDC+NSE --- integration/nse_update_simplified_sdc.H | 87 +++++++++++++++++++++++++ 1 file changed, 87 insertions(+) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 5f3e9c5d6c..7faea0d073 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -26,6 +26,93 @@ using namespace amrex; #if defined(NSE_TABLE) + +/// +/// this acts as the RHS function for the system (rho e, rho aux) +/// on input, *_source are the reactive sources at time t0 and on output +/// they are the sources at time t0+dt +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void rhs(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux0, + const Real dt, const Real *ydot_a, + Real& rhoe_new, Real* rhoaux_new, + Real& rhoe_source, Real* rhoaux_source, const int T_fixed) { + + // integrate rho to the new time + + Real rho_new = rho0 + dt * ydot_a[SRHO]; + + // integrate rhoe and rhoaux using the input sources + + rhoe_new = rhoe + dt * ydot_a[SEINT] + dt * rhoe_source; + for (int n = 0; n < NumAux; n++) { + rhoaux_new[n] = rhoaux0[n] + dt * ydot_a[SFX+n] + dt * aux_source[n]; + } + + // find the temperature at t0 + dt + + eos_t eos_state; + eos_state.rho = rho_new; + eos_state.e = rhoe_new / rho_new; + eos_state.T = T0; + for (int n = 0; n < NumAux; n++) { + eos_state.aux[n] = rhoaux_new[n] / rho_new; + } + + Real T_new; + if (T_fixed > 0) { + T_new = T0; + } else { + eos(eos_input_re, eos_state); + T_new = eos_state.T; + } + + // compute the plasma neutrino losses at t0 + dt + + Real snu{0.0}; + Real dsnudt{0.0}; + Real dsnudd{0.0}; + Real dsnuda{0.0}; + Real dsnudz{0.0}; + + Real abar = eos_state.aux[iabar]; + Real zbar = abar * eos_state.aux[iye]; + + constexpr int do_derivatives = 0; + sneut5(T_new, rho_new, abar, zbar, + snu, dsnudt, dsnudd, dsnuda, dsnudz); + + // call the NSE table + + Real abar_out; + Real dq_out; + Real dyedt; + Real dabardt; + Real dbeadt; + Real enu; + Real X[NumSpec]; + + constexpr bool skip_X_fill = true; + + nse_interp(T_new, rho_new, eos_state,aux[iye], + abar_out, bea_out, dyedt, dabardt, dbeadt, enu, skip_X_fill); + + // construct the sources at time t0+dt + + // note that abar, (B/A) here have already seen advection + // implicitly + + Real rho_bea_tilde = rho_new * bea_out - dt * ydot_a[SFX+ibea]; + Real rho_dBEA = rho_bea_tilde - rhoaux0[ibea]; // this is MeV / nucleon * g / cm**3 + Real rhoe_source = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; + rhoe_source -= rho_new * (enu + snu); + + Real rhoaux_source[iabar] = 0.0; + Real rhoaux_source[iye] = rho_new * dyedt; + Real rhoaux_source[ibea] = rho_dBEA / dt +; +} + /// /// update the state due to NSE changes for the simplified-SDC case /// this version works with the tabulated NSE and requires AUX_THERMO From 4b1e230a7f25aa13482412b37b2d2a71d434340a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 16 Nov 2023 12:39:36 -0500 Subject: [PATCH 06/13] I think this is what I want --- integration/nse_update_simplified_sdc.H | 218 ++++++++++++------------ 1 file changed, 111 insertions(+), 107 deletions(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 7faea0d073..38b36dcc27 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -28,15 +28,16 @@ using namespace amrex; #if defined(NSE_TABLE) /// -/// this acts as the RHS function for the system (rho e, rho aux) +/// this acts as an explicit Euler step for the system (rho e, rho aux) /// 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 rhs(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux0, - const Real dt, const Real *ydot_a, - Real& rhoe_new, Real* rhoaux_new, - Real& rhoe_source, Real* rhoaux_source, const int T_fixed) { +void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux0, + const Real dt, const Real *ydot_a, + Real& rhoe_new, Real* rhoaux_new, Real* rhoX_new, + Real& rhoe_source, Real* rhoaux_source, const int T_fixed) { // integrate rho to the new time @@ -92,10 +93,15 @@ void rhs(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux0, Real enu; Real X[NumSpec]; - constexpr bool skip_X_fill = true; + nse_interp(T_new, rho_new, eos_state.aux[iye], + abar_out, bea_out, dyedt, dabardt, dbeadt, enu, rhoX_new, + skip_X_fill); - nse_interp(T_new, rho_new, eos_state,aux[iye], - abar_out, bea_out, dyedt, dabardt, dbeadt, enu, skip_X_fill); + if constexpr (!skip_X_fill) { + for (auto & X : rhoX_new) { + X *= rho_new; + } + } // construct the sources at time t0+dt @@ -134,8 +140,16 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // into NSE (the first call to nse_interp) + an energy // release from the evolution of NSE over the timestep + // we will write the energy update as: + // + // e^new = e^old + dbea_instantantous + dbea_over_dt + // + // where dbea_over_dt will be an integral over dt using + // a simple RK2 method + Real ye_in = state.y[SFX+iye] / state.rho; - Real rho_bea_old = state.y[SFX+ibea]; + Real rho_bea_orig = state.y[SFX+ibea]; + Real rhoe_orig = state.y[SEINT]; // density and momentum have no reactive sources Real rho_old = state.y[SRHO]; @@ -150,151 +164,141 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real T_in = state.T_fixed > 0.0_rt ? state.T_fixed : state.T; - // get the neutrino loss term -- we want to use the state that we - // came in here with, so the original Abar and Zbar + // get the instantaneous change in B/A, as well as the source + // terms for the new evolution. - Real snu{0.0}; - Real dsnudt{0.0}; - Real dsnudd{0.0}; - Real dsnuda{0.0}; - Real dsnudz{0.0}; - - Real abar = state.y[SFX+iabar] / rho_old; - Real zbar = abar * ye_in; - - constexpr int do_derivatives = 0; - sneut5(T_in, rho_old, abar, zbar, - snu, dsnudt, dsnudd, dsnuda, dsnudz); - - Real snu_old = snu; - - // get the current NSE state from the table -- this will be used - // to compute the NSE evolution sources - - Real abar_out; - Real dq_out; - Real dyedt; - Real dabardt; - Real dbeadt; - Real enu; - Real X[NumSpec]; + constexpr bool skip_X_fill{true}; nse_interp(T_in, state.rho, ye_in, - abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X); - - Real dyedt_old = dyedt; - - - // predict the U^{n+1,*} state with only estimates of the aux - // reaction terms and advection to dt + abar_out, bea_out, dyedt, dabardt, dbeadt, enu, X, skip_X_fill); - eos_re_t eos_state; - eos_state.T = T_in; // initial guess + // the starting state for integration uses the state that we just + // adjusted instantaneously into NSE - // initial aux_sources - Real aux_source[NumAux] = {0.0_rt}; - aux_source[iye] = rho_old * dyedt; + Real rhoaux_old[NumAux]; + rhoaux_old[iye] = rho_old * ye_in; + rhoaux_old[iabar] = rho_old * abar_out; + rhoaux_old[ibea] = rho_old * bea_out; - Real rho_aux_new[NumAux]; + Real rhoenuc = (rho_old * bea_out - rho_bea_orig) * C::MeV2eV * C::ev2erg * C::n_A; + Real rhoe_old = state.y[SEINT] + rhoenuc; - Real rhoe_new; - Real rho_enucdot = -rho_old * snu; - - Real rho_half = 0.5_rt * (rho_old + state.y[SRHO]); - // predict the new aux for the first iteration -- this is really - // just including the advection bits + // get the starting temperature - for (int n = 0; n < NumAux; n++) { - rho_aux_new[n] = state.y[SFX+n] + dt * state.ydot_a[SFX+n] + dt * aux_source[n]; + eos_t eos_state; + eos_state.rho = rho_old; + eos_state.e = rhoe_old / rho_old; + for (int n = 0; n < NumAux; ++n) { + eos_state.aux[n] = rhoaux_old[n] / rho_old; } + eos_state.T = T_in; - for (int iter = 0; iter < integrator_rp::nse_iters; iter++) { + Real T_old; + if (state.T_fixed > 0) { + T_old = state.T_fixed; + } else { + eos(eos_input_re, eos_state); + T_old = eos_state.T; + } - // 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; + // start of integration - // call the EOS to get the updated T* + // we want to integrate from t^n to t^{n+1/2} + // start by computing the sources - Real T_new; - eos_state.rho = state.y[SRHO]; - eos_state.e = rhoe_new / state.y[SRHO]; - for (int n = 0; n < NumAux; n++) { - eos_state.aux[n] = rho_aux_new[n] / state.y[SRHO]; - } + // the aux sources are just due to weak rates + Real rhoaux_source[NumAux]; + rhoaux_source[iye] = rho_old * dyedt; + rhoaux_source[ibea] = rho_old * dbeadt; + rhoaux_source[iabar] = 0.0; - if (state.T_fixed > 0) { - T_new = state.T_fixed; - } else { - eos(eos_input_re, eos_state); - T_new = eos_state.T; - } + // the energy source is due to plasma neutrinos + weak rate neutrinos - // call the NSE table using the * state to get the t^{n+1} - // source estimates - - nse_interp(T_new, eos_state.rho, eos_state.aux[iye], - abar_out, dq_out, dyedt, dabardt, dbeadt, enu, X); + Real snu{}; + Real dsnudt{}; + Real dsnudd{}; + Real dsnuda{}; + Real dsnudz{}; + constexpr int do_derivatives = 0; - // note that the abar, dq (B/A), and X here have all already - // seen advection implicitly + Real zbar = abar_out * ye_in; + sneut5(T_old, rho_old, abar_out, zbar, + snu, dsnudt, dsnudd, dsnuda, dsnudz); + Real rhoe_source = -rho_old * (snu + enu); - // compute the energy release -- we need to remove the advective part + Real rhoe_new; + Real rhoaux_new[NumAux]; + Real rhoe_source; + Real rhoaux_source[NumAux]; + Real rhoX_new[NumSpec]; - Real rho_bea_tilde = state.y[SRHO] * dq_out - dt * state.ydot_a[SFX+ibea]; - Real rho_dBEA = rho_bea_tilde - rho_bea_old; // this is MeV / nucleon * g / cm**3 + // evolve for dt/2 -- the sources will be overwritten with the + // sources at the end point of this integration - // convert the energy to erg / cm**3 - rho_enucdot = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; + { + constexpr bool skip_X_fill{true}; - // now get the updated neutrino term - zbar = abar_out * eos_state.aux[iye]; - sneut5(T_new, eos_state.rho, abar_out, zbar, - snu, dsnudt, dsnudd, dsnuda, dsnudz); + evolve(rho_old, T_old, rhoe_old, rhoaux_old, + 0.5_rt*dt, state.ydot_a, + rhoe_new, rhoaux_new, rhoX_new, + rhoe_source, rhoaux_source, state.T_fixed); - rho_enucdot -= 0.5_rt * rho_half * (snu_old + snu); + } - // update the new state for the next pass + // this gives us the sources at n+1/2. Now we do the full dt step, + // starting again from the beginning, but using the time-centered + // sources - aux_source[iye] = 0.5_rt * rho_half * (dyedt_old + dyedt); - rho_aux_new[iye] = state.y[SFX+iye] + dt * state.ydot_a[SFX+iye] + dt * aux_source[iye]; + { + constexpr bool skip_X_fill{false}; - rho_aux_new[iabar] = state.y[SRHO] * abar_out; - rho_aux_new[ibea] = state.y[SRHO] * dq_out; + evolve(rho_old, T_old, rhoe_old, rhoaux_old, + dt, state.ydot_a, + rhoe_new, rhoaux_new, rhoX_new, + rhoe_source, rhoaux_source, state.T_fixed); } + // done with the integration. This gave us a second-order + // accurate update to the rhoe and Ye, and the aux data and mass + // fractions are consistent with the NSE table at the end of the + // integration + // now update the aux quantities // the new mass fractions are just those that come from the table // make sure they are normalized - Real sum_X{0.0_rt}; - for (auto & xn : X) { - xn = amrex::max(small_x, amrex::min(1.0_rt, xn)); - sum_X += xn; + Real sum_rhoX{0.0_rt}; + for (auto & rx : rhoX_new) { + rx = amrex::max(rho_new*small_x, amrex::min(rho_new, rx)); + sum_rhoX += rx; } - for (auto & xn : X) { - xn /= sum_X; + for (auto & rx : rhoX_new) { + rx /= sum_rhoX; } for (int n = 0; n < NumSpec; ++n) { - state.y[SFS+n] = state.y[SRHO] * X[n]; + state.y[SFS+n] = rhoX_new[n]; } // aux data comes from the table (except Ye, which we predicted) - state.y[SFX+iye] = eos_state.rho * eos_state.aux[iye]; - state.y[SFX+iabar] = eos_state.rho * abar_out; - state.y[SFX+ibea] = eos_state.rho * dq_out; + for (int n = 0; n < NumAux; ++n) { + state.y[SFX+n] = rhoaux_new[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; + Real rho_enucdot = (rhoe_new - rhoe_orig) / dt - state.ydot_a[SEINT]; + + state.y[SEINT] = rhoe_new; state.y[SEDEN] += dt * state.ydot_a[SEDEN] + dt * rho_enucdot; } From c05dd017e66ce97cf65d70125c9dc314b9f8714f Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 16 Nov 2023 12:49:06 -0500 Subject: [PATCH 07/13] this compiles --- integration/nse_update_simplified_sdc.H | 31 +++++++++++++++---------- 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 38b36dcc27..47683b618c 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -45,9 +45,9 @@ void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux // integrate rhoe and rhoaux using the input sources - rhoe_new = rhoe + dt * ydot_a[SEINT] + dt * rhoe_source; + rhoe_new = rhoe0 + dt * ydot_a[SEINT] + dt * rhoe_source; for (int n = 0; n < NumAux; n++) { - rhoaux_new[n] = rhoaux0[n] + dt * ydot_a[SFX+n] + dt * aux_source[n]; + rhoaux_new[n] = rhoaux0[n] + dt * ydot_a[SFX+n] + dt * rhoaux_source[n]; } // find the temperature at t0 + dt @@ -86,7 +86,7 @@ void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux // call the NSE table Real abar_out; - Real dq_out; + Real bea_out; Real dyedt; Real dabardt; Real dbeadt; @@ -98,8 +98,8 @@ void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux skip_X_fill); if constexpr (!skip_X_fill) { - for (auto & X : rhoX_new) { - X *= rho_new; + for (int n = 0; n < NumSpec; ++n) { + rhoX_new[n] *= rho_new; } } @@ -110,12 +110,13 @@ void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux Real rho_bea_tilde = rho_new * bea_out - dt * ydot_a[SFX+ibea]; Real rho_dBEA = rho_bea_tilde - rhoaux0[ibea]; // this is MeV / nucleon * g / cm**3 - Real rhoe_source = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; + + rhoe_source = rho_dBEA * C::MeV2eV * C::ev2erg * C::n_A / dt; rhoe_source -= rho_new * (enu + snu); - Real rhoaux_source[iabar] = 0.0; - Real rhoaux_source[iye] = rho_new * dyedt; - Real rhoaux_source[ibea] = rho_dBEA / dt + rhoaux_source[iabar] = 0.0; + rhoaux_source[iye] = rho_new * dyedt; + rhoaux_source[ibea] = rho_dBEA / dt; ; } @@ -167,6 +168,14 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // get the instantaneous change in B/A, as well as the source // terms for the new evolution. + Real abar_out{}; + Real bea_out{}; + Real dyedt{}; + Real dabardt{}; + Real dbeadt{}; + Real enu{}; + Real X[NumSpec] = {}; + constexpr bool skip_X_fill{true}; nse_interp(T_in, state.rho, ye_in, @@ -230,8 +239,6 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real rhoe_new; Real rhoaux_new[NumAux]; - Real rhoe_source; - Real rhoaux_source[NumAux]; Real rhoX_new[NumSpec]; // evolve for dt/2 -- the sources will be overwritten with the @@ -272,7 +279,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // make sure they are normalized Real sum_rhoX{0.0_rt}; for (auto & rx : rhoX_new) { - rx = amrex::max(rho_new*small_x, amrex::min(rho_new, rx)); + rx = amrex::max(state.y[SRHO] * small_x, amrex::min(state.y[SRHO], rx)); sum_rhoX += rx; } From 3da1684269ca2f83f65cad6194fec6d7b954f174 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 17 Nov 2023 08:54:41 -0500 Subject: [PATCH 08/13] fix type error --- integration/nse_update_simplified_sdc.H | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 47683b618c..f9e3a6fa87 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -37,7 +37,7 @@ AMREX_GPU_HOST_DEVICE AMREX_INLINE void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux0, const Real dt, const Real *ydot_a, Real& rhoe_new, Real* rhoaux_new, Real* rhoX_new, - Real& rhoe_source, Real* rhoaux_source, const int T_fixed) { + Real& rhoe_source, Real* rhoaux_source, const Real T_fixed) { // integrate rho to the new time @@ -176,10 +176,12 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real enu{}; Real X[NumSpec] = {}; - constexpr bool skip_X_fill{true}; + { + constexpr bool skip_X_fill{true}; - nse_interp(T_in, state.rho, ye_in, - abar_out, bea_out, dyedt, dabardt, dbeadt, enu, X, skip_X_fill); + nse_interp(T_in, state.rho, ye_in, + abar_out, bea_out, dyedt, dabardt, dbeadt, enu, X, skip_X_fill); + } // the starting state for integration uses the state that we just // adjusted instantaneously into NSE @@ -251,7 +253,6 @@ void sdc_nse_burn(BurnT& state, const Real dt) { 0.5_rt*dt, state.ydot_a, rhoe_new, rhoaux_new, rhoX_new, rhoe_source, rhoaux_source, state.T_fixed); - } // this gives us the sources at n+1/2. Now we do the full dt step, @@ -265,7 +266,6 @@ void sdc_nse_burn(BurnT& state, const Real dt) { dt, state.ydot_a, rhoe_new, rhoaux_new, rhoX_new, rhoe_source, rhoaux_source, state.T_fixed); - } // done with the integration. This gave us a second-order From c27b37805dd5d07d40e41a3d1dc1cb65b542f480 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 17 Nov 2023 09:43:31 -0500 Subject: [PATCH 09/13] fix comp warning --- integration/nse_update_simplified_sdc.H | 1 - 1 file changed, 1 deletion(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index f9e3a6fa87..6ed3d8ecdb 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -91,7 +91,6 @@ void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux Real dabardt; Real dbeadt; Real enu; - Real X[NumSpec]; nse_interp(T_new, rho_new, eos_state.aux[iye], abar_out, bea_out, dyedt, dabardt, dbeadt, enu, rhoX_new, From f3fb3fc2f3e836d1d181059f6464cbb1ffa9c52a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 17 Nov 2023 10:46:45 -0500 Subject: [PATCH 10/13] a few minor fixes --- integration/nse_update_simplified_sdc.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 6ed3d8ecdb..4bfa41f651 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -62,7 +62,7 @@ void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux Real T_new; if (T_fixed > 0) { - T_new = T0; + T_new = T_fixed; } else { eos(eos_input_re, eos_state); T_new = eos_state.T; @@ -290,7 +290,7 @@ void sdc_nse_burn(BurnT& state, const Real dt) { state.y[SFS+n] = rhoX_new[n]; } - // aux data comes from the table (except Ye, which we predicted) + // aux data comes from the integration for (int n = 0; n < NumAux; ++n) { state.y[SFX+n] = rhoaux_new[n]; From 60644d720382160c2a45ff801b319025cc7c1fb4 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 18 Nov 2023 14:52:50 -0500 Subject: [PATCH 11/13] fix the update by preserving B/A --- integration/nse_update_simplified_sdc.H | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 4bfa41f651..48c77c689a 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -116,7 +116,13 @@ void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux rhoaux_source[iabar] = 0.0; rhoaux_source[iye] = rho_new * dyedt; rhoaux_source[ibea] = rho_dBEA / dt; -; + + // reset the new -- we want to ensure that if we use the + // updated T(rho, e, Y_e, Abar) that we get exactly the same B/A, + // so there is no instantaneous change in the binding energy + + rhoaux_new[ibea] = rho_new * bea_out; + } /// From c5b053292846c81b9de3ad78e5903b80eb6b306d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 19 Nov 2023 16:58:17 -0500 Subject: [PATCH 12/13] store abar too --- integration/nse_update_simplified_sdc.H | 1 + 1 file changed, 1 insertion(+) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 48c77c689a..f14ae6021d 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -122,6 +122,7 @@ void evolve(const Real rho0, const Real T0, const Real rhoe0, const Real *rhoaux // so there is no instantaneous change in the binding energy rhoaux_new[ibea] = rho_new * bea_out; + rhoaux_new[iabar] = rho_new * abar_out; } From a0d6bf34ced80d308618fa7ee295a9df315cd35a Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 20 Nov 2023 11:34:52 -0500 Subject: [PATCH 13/13] this does a better job at getting second order --- integration/nse_update_simplified_sdc.H | 50 +++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/integration/nse_update_simplified_sdc.H b/integration/nse_update_simplified_sdc.H index 48c77c689a..edd18c3de0 100644 --- a/integration/nse_update_simplified_sdc.H +++ b/integration/nse_update_simplified_sdc.H @@ -248,6 +248,24 @@ void sdc_nse_burn(BurnT& state, const Real dt) { Real rhoaux_new[NumAux]; Real rhoX_new[NumSpec]; + // evolve for eps dt -- the sources will be overwritten with the + // sources at the end point of this integration + + { + constexpr bool skip_X_fill{true}; + + evolve(rho_old, T_old, rhoe_old, rhoaux_old, + 0.05_rt*dt, state.ydot_a, + rhoe_new, rhoaux_new, rhoX_new, + rhoe_source, rhoaux_source, state.T_fixed); + } + + //std::cout << "starting integration, dt = " << dt << std::endl; + //std::cout << " rhoe: " << rhoe_old << " " << rhoe_source << std::endl; + //std::cout << " rhoYe: " << rhoaux_old[iye] << " " << rhoaux_source[iye] << std::endl; + //std::cout << " rhoAbar: " << rhoaux_old[iabar] << " " << rhoaux_source[iabar] << std::endl; + //std::cout << " rhoBA: " << rhoaux_old[ibea] << " " << rhoaux_source[ibea] << std::endl; + // evolve for dt/2 -- the sources will be overwritten with the // sources at the end point of this integration @@ -260,10 +278,34 @@ void sdc_nse_burn(BurnT& state, const Real dt) { rhoe_source, rhoaux_source, state.T_fixed); } + // evolve from dt/2 to dt/2 + eps -- the sources will be overwritten with the + // sources at the end point of this integration + + { + constexpr bool skip_X_fill{true}; + + Real rho_tmp = rho_old + 0.5*dt * state.ydot_a[SRHO]; + Real rhoe_tmp{rhoe_new}; + Real rhoaux_tmp[NumAux] = {rhoaux_new[0], rhoaux_new[1], rhoaux_new[2]}; + + evolve(rho_tmp, T_old, rhoe_tmp, rhoaux_tmp, + 0.05_rt*dt, state.ydot_a, + rhoe_new, rhoaux_new, rhoX_new, + rhoe_source, rhoaux_source, state.T_fixed); + } + // this gives us the sources at n+1/2. Now we do the full dt step, // starting again from the beginning, but using the time-centered // sources + + //std::cout << "second half of integration, dt = " << dt << std::endl; + //std::cout << " rhoe: " << rhoe_old << " " << rhoe_source << std::endl; + //std::cout << " rhoYe: " << rhoaux_old[iye] << " " << rhoaux_source[iye] << std::endl; + //std::cout << " rhoAbar: " << rhoaux_old[iabar] << " " << rhoaux_source[iabar] << std::endl; + //std::cout << " rhoBA: " << rhoaux_old[ibea] << " " << rhoaux_source[ibea] << std::endl; + + { constexpr bool skip_X_fill{false}; @@ -278,6 +320,14 @@ void sdc_nse_burn(BurnT& state, const Real dt) { // fractions are consistent with the NSE table at the end of the // integration + //std::cout << "done" << std::endl; + //std::cout << " rhoe: " << rhoe_new << std::endl; + //std::cout << " rhoYe: " << rhoaux_new[iye] << std::endl; + //std::cout << " rhoAbar: " << rhoaux_new[iabar] << std::endl; + //std::cout << " rhoBA: " << rhoaux_new[ibea] << std::endl; + + //amrex::Error("stop"); + // now update the aux quantities // the new mass fractions are just those that come from the table