Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

reuse the simplified-SDC integrator for the ODE update #2584

Merged
merged 4 commits into from
Sep 29, 2023
Merged

reuse the simplified-SDC integrator for the ODE update #2584

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8d6afff
reuse the simplified-SDC integrator for the ODE update
zingale Sep 25, 2023
25248f9
fix compilation
zingale Sep 25, 2023
2127dfd
Merge branch 'development' into sdc_vode_cleanup
zingale Sep 25, 2023
7407195
Merge branch 'development' into sdc_vode_cleanup
zingale Sep 29, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
190 changes: 42 additions & 148 deletions Source/sdc/vode_rhs_true_sdc.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,97 +12,8 @@
#include <actual_rhs.H>
#endif
#include <Castro_react_util.H>
#include <sdc_react_util.H>

#include <vode_dvode.H>

// The f_rhs routine provides the right-hand-side for the DVODE solver.
// This is a generic interface that calls the specific RHS routine in the
// network you're actually using.

template <typename DvodeT>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void rhs(const Real time, burn_t& burn_state, DvodeT& vode_state, RArray1D& dUdt, const bool in_jacobian=false)
{

// note: dUdt is 1-based

GpuArray<Real, NUM_STATE> U_full;
GpuArray<Real, NUM_STATE> R_full;

// evaluate R

// update the density -- we'll need this for the EOS call in the
// react update

U_full[URHO] = burn_state.rho_orig + time * burn_state.ydot_a[SRHO];

// we are not solving the momentum equations
// and never need total energy, so we don't worry about filling those

for (int n = 0; n < NumSpec; n++) {
U_full[UFS+n] = vode_state.y(1+n);
}
U_full[UEINT] = vode_state.y(NumSpec+1);

// initialize the temperature -- a better value will be found when we do the EOS
// call in single_zone_react_source

U_full[UTEMP] = burn_state.T;

// create a temporary burn_t for this call. It is just going to
// get loaded up with U_full

burn_t burn_state_pass;
single_zone_react_source(U_full, R_full, burn_state_pass);

// update our temperature for next time

burn_state.T = burn_state_pass.T;

// now construct the RHS -- note this is 1-based

for (int n = 0; n < NumSpec; ++n) {
dUdt(n+1) = R_full[UFS+n] + burn_state.ydot_a[SFS+n];
}
dUdt(NumSpec+1) = R_full[UEINT] + burn_state.ydot_a[SEINT];

}


template<class MatrixType, typename DvodeT>
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void jac (const Real time, burn_t& burn_state, DvodeT& vode_state, MatrixType& pd)
{
// Call the specific network routine to get the Jacobian.

// note the Jacobian is 1-based at the moment

GpuArray<Real, NUM_STATE> U_full;
GpuArray<Real, NUM_STATE> R_full;

// we are not solving the momentum equations
// create a full state -- we need this for some interfaces

U_full[URHO] = burn_state.rho_orig + time * burn_state.ydot_a[URHO];

for (int n = 0; n < NumSpec; n++) {
U_full[UFS+n] = vode_state.y(1+n);
}
U_full[UEINT] = vode_state.y(NumSpec+1);

U_full[UTEMP] = burn_state.T;

burn_t burn_state_pass;

// this will initialize burn_state_pass and fill the temperature
// via the EOS

single_zone_react_source(U_full, R_full, burn_state_pass);

single_zone_jac(U_full, burn_state_pass, vode_state.H, pd);

}
#include <burner.H>


AMREX_GPU_HOST_DEVICE AMREX_INLINE
Expand All @@ -122,13 +33,13 @@ sdc_vode_solve(const Real dt_m,
// We will do the implicit update of only the terms that
// have reactive sources
//
// 0 : rho
// 1:NumSpec : species
// NumSpec+1 : (rho E) or (rho e)
// 1:NumSpec : (rho X)
// NumSpec+1 : (rho e)

#if (INTEGRATOR == 0)

// Update the momenta for this zone -- they don't react
// Update the density and momenta for this zone -- they don't react

U_new[URHO] = U_old[URHO] + dt_m * C[URHO];

for (int n = 0; n < 3; ++n) {
U_new[UMX+n] = U_old[UMX+n] + dt_m * C[UMX+n];
Expand All @@ -138,15 +49,32 @@ sdc_vode_solve(const Real dt_m,
// nonlinear solve -- note: we include the old state
// in f_source

// load rpar
burn_t burn_state;

burn_state.success = true;

dvode_t<NumSpec+1> dvode_state;
// store the state

burn_t burn_state;
burn_state.y[SRHO] = U_old[URHO];
burn_state.y[SMX] = U_old[UMX];
burn_state.y[SMY] = U_old[UMY];
burn_state.y[SMZ] = U_old[UMZ];
burn_state.y[SEDEN] = U_old[UEDEN];
burn_state.y[SEINT] = U_old[UEINT];
for (int n = 0; n < NumSpec; n++) {
burn_state.y[SFS+n] = U_old[UFS+n];
}
#if NAUX_NET > 0
for (int n = 0; n < NumAux; n++) {
burn_state.y[SFX+n] = U_old[UFX+n];
}
#endif

burn_state.rho_orig = U_old[URHO];
burn_state.rhoe_orig = U_old[UEINT];
// we need an initial T guess for the EOS
burn_state.T = U_old[UTEMP];

burn_state.T_fixed = -1.e30_rt;
burn_state.rho = burn_state.y[SRHO];

// If we are solving the system as an ODE, then we are
// solving
Expand All @@ -171,65 +99,31 @@ sdc_vode_solve(const Real dt_m,
}
#endif

burn_state.success = true;

// temperature will be used as an initial guess in the EOS
burn_state.T = U_old[UTEMP];


// store the subset for the nonlinear solve. We only
// consider (rho e), not (rho E). This is because at
// present we do not have a method of updating the
// velocities during the multistep integration.
burn_state.sdc_iter = sdc_iteration;
burn_state.num_sdc_iters = sdc_order + sdc_extra;

// Also note that the dvode_state is 1-based, but we'll
// access it as 0-based in our implementation of the
// RHS routine
burner(burn_state, dt_m);

for (int n = 0; n < NumSpec; ++n) {
dvode_state.y(1+n) = U_old[UFS+n];
if (! burn_state.success) {
amrex::Error("Error: integration failed");
}
dvode_state.y(NumSpec+1) = U_old[UEINT];

// // set the maximum number of steps allowed
// dvode_state.MXSTEP = 25000;

dvode_state.t = 0.0_rt;
dvode_state.tout = dt_m;
dvode_state.HMXI = 1.0_rt / ode_max_dt;

dvode_state.jacobian_type = static_cast<short>(integrator_rp::jacobian);
// update the full U_new -- we already did density and momentum

// relative tolerances
dvode_state.rtol_spec = integrator_rp::rtol_spec;
dvode_state.rtol_enuc = integrator_rp::rtol_enuc;

// absolute tolerances
dvode_state.atol_spec = integrator_rp::atol_spec * U_old[URHO]; // this way, atol is the minimum x
dvode_state.atol_enuc = integrator_rp::atol_enuc * U_old[UEINT];

int istate = dvode(burn_state, dvode_state);

if (istate < 0) {
Abort("Vode terminated poorly");
}

// update the full U_new
for (int n = 0; n < NumSpec; ++n) {
U_new[UFS+n] = dvode_state.y(1+n);
U_new[UEDEN] = burn_state.y[SEDEN];
U_new[UEINT] = burn_state.y[SEINT];
for (int n = 0; n < NumSpec; n++) {
U_new[UFS+n] = burn_state.y[SFS+n];
}
U_new[UEINT] = dvode_state.y(NumSpec+1);

auto v2 = 0.0_rt;
for (int n = 0; n < 3; ++n) {
v2 += U_new[UMX+n] * U_new[UMX+n];
#if NAUX_NET > 0
for (int n = 0; n < NumAux; n++) {
U_new[UFX+n] = burn_state.y[SFX+n];
}
U_new[UEDEN] = U_new[UEINT] + 0.5_rt * v2 / U_new[URHO];
#endif

// keep our temperature guess
U_new[UTEMP] = U_old[UTEMP];

#endif
}

#endif
Expand Down