From 92c5accae807788b49ee1b60bd77c17bc66fead9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 18 Oct 2023 12:49:10 -0400 Subject: [PATCH 1/4] 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 2/4] 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 3/4] 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 4/4] 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,