Skip to content

Commit

Permalink
move the true SDC VODE into the VODE header (#2567)
Browse files Browse the repository at this point in the history 
this puts all of the VODE stuff together
  • Loading branch information
@zingale
zingale authored Sep 21, 2023
1 parent a02f2dc commit fa367de
Show file tree
Hide file tree
Showing 2 changed files with 139 additions and 139 deletions.
140 changes: 1 addition & 139 deletions Source/sdc/Castro_sdc_util.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,20 +8,14 @@
#endif
#include <Castro_react_util.H>
#include <sdc_newton_solve.H>
#include <vode_dvode.H>
#ifndef NEW_NETWORK_IMPLEMENTATION
#include <linpack.H>
#endif
#include <vode_rhs_true_sdc.H>
#endif

// solvers
constexpr int NEWTON_SOLVE = 1;
constexpr int VODE_SOLVE = 2;
constexpr int HYBRID_SOLVE = 3;

constexpr int ldjac = NumSpec + 2;


AMREX_GPU_HOST_DEVICE AMREX_INLINE
void
normalize_species_sdc(const int i, const int j, const int k,
Expand All @@ -47,138 +41,6 @@ normalize_species_sdc(const int i, const int j, const int k,



AMREX_GPU_HOST_DEVICE AMREX_INLINE
void
sdc_vode_solve(const Real dt_m,
GpuArray<Real, NUM_STATE> const& U_old,
GpuArray<Real, NUM_STATE>& U_new,
GpuArray<Real, NUM_STATE> const& C,
const int sdc_iteration) {

// The purpose of this function is to solve the system the
// approximate system dU/dt = R + C using the VODE ODE
// integrator.
// The solution we get here will then be used as the
// initial guess to the Newton solve on the real system.

// 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)

#if (INTEGRATOR == 0)

// The tolerance we are solving to may depend on the
// iteration
auto relax_fac = std::pow(sdc_solver_relax_factor, sdc_order - sdc_iteration - 1);
auto tol_dens = sdc_solver_tol_dens * relax_fac;
auto tol_spec = sdc_solver_tol_spec * relax_fac;
auto tol_ener = sdc_solver_tol_ener * relax_fac;

// Update the momenta for this zone -- they don't react
for (int n = 0; n < 3; ++n) {
U_new[UMX+n] = U_old[UMX+n] + dt_m * C[UMX+n];
}

// Now only save the subset that participates in the
// nonlinear solve -- note: we include the old state
// in f_source

// load rpar

// If we are solving the system as an ODE, then we are
// solving
// dU/dt = R(U) + C
// so we simply pass in C
GpuArray<Real, NumSpec+2> C_react;

C_react[0] = C[URHO];
for (int n = 0; n < NumSpec; ++n) {
C_react[1+n] = C[UFS+n];
}
C_react[NumSpec+1] = C[UEINT];

dvode_t<NumSpec+2> dvode_state;

burn_t burn_state;

burn_state.success = true;

for (int n = 0; n < NumSpec+2; n++) {
burn_state.f_source[n] = C_react[n];
}

for (int n = 0; n < 3; n++) {
burn_state.mom[n] = U_new[UMX+n];
}

// we are always solving for rhoe with the VODE predict
burn_state.e = U_new[UEDEN];

burn_state.E_var = U_old[UEDEN];

// 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.

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

dvode_state.y(1) = U_old[URHO];
for (int n = 0; n < NumSpec; ++n) {
dvode_state.y(2+n) = U_old[UFS+n];
}
dvode_state.y(NumSpec+2) = 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;

// relative tolerances
dvode_state.rtol_dens = tol_dens;
dvode_state.rtol_spec = tol_spec;
dvode_state.rtol_enuc = tol_ener;

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

int istate = dvode(burn_state, dvode_state);

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

// update the full U_new
U_new[URHO] = dvode_state.y(1);
for (int n = 0; n < NumSpec; ++n) {
U_new[UFS+n] = dvode_state.y(2+n);
}
U_new[UEINT] = dvode_state.y(NumSpec+2);
auto v2 = 0.0_rt;
for (int n = 0; n < 3; ++n) {
v2 += U_new[UMX+n] * U_new[UMX+n];
}
U_new[UEDEN] = U_new[UEINT] + 0.5_rt * v2 / U_new[URHO];

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

#endif
}

AMREX_GPU_HOST_DEVICE AMREX_INLINE
void
sdc_solve(const Real dt_m,
Expand Down
138 changes: 138 additions & 0 deletions Source/sdc/vode_rhs_true_sdc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,11 @@
#include <Castro_react_util.H>
#include <sdc_react_util.H>

#include <vode_dvode.H>
#ifndef NEW_NETWORK_IMPLEMENTATION
#include <linpack.H>
#endif

// 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.
Expand Down Expand Up @@ -176,5 +181,138 @@ void jac (const Real time, burn_t& burn_state, DvodeT& vode_state, MatrixType& p
}

}

AMREX_GPU_HOST_DEVICE AMREX_INLINE
void
sdc_vode_solve(const Real dt_m,
GpuArray<Real, NUM_STATE> const& U_old,
GpuArray<Real, NUM_STATE>& U_new,
GpuArray<Real, NUM_STATE> const& C,
const int sdc_iteration) {

// The purpose of this function is to solve the system the
// approximate system dU/dt = R + C using the VODE ODE
// integrator.
// The solution we get here will then be used as the
// initial guess to the Newton solve on the real system.

// 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)

#if (INTEGRATOR == 0)

// The tolerance we are solving to may depend on the
// iteration
auto relax_fac = std::pow(sdc_solver_relax_factor, sdc_order - sdc_iteration - 1);
auto tol_dens = sdc_solver_tol_dens * relax_fac;
auto tol_spec = sdc_solver_tol_spec * relax_fac;
auto tol_ener = sdc_solver_tol_ener * relax_fac;

// Update the momenta for this zone -- they don't react
for (int n = 0; n < 3; ++n) {
U_new[UMX+n] = U_old[UMX+n] + dt_m * C[UMX+n];
}

// Now only save the subset that participates in the
// nonlinear solve -- note: we include the old state
// in f_source

// load rpar

// If we are solving the system as an ODE, then we are
// solving
// dU/dt = R(U) + C
// so we simply pass in C
GpuArray<Real, NumSpec+2> C_react;

C_react[0] = C[URHO];
for (int n = 0; n < NumSpec; ++n) {
C_react[1+n] = C[UFS+n];
}
C_react[NumSpec+1] = C[UEINT];

dvode_t<NumSpec+2> dvode_state;

burn_t burn_state;

burn_state.success = true;

for (int n = 0; n < NumSpec+2; n++) {
burn_state.f_source[n] = C_react[n];
}

for (int n = 0; n < 3; n++) {
burn_state.mom[n] = U_new[UMX+n];
}

// we are always solving for rhoe with the VODE predict
burn_state.e = U_new[UEDEN];

burn_state.E_var = U_old[UEDEN];

// 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.

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

dvode_state.y(1) = U_old[URHO];
for (int n = 0; n < NumSpec; ++n) {
dvode_state.y(2+n) = U_old[UFS+n];
}
dvode_state.y(NumSpec+2) = 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;

// relative tolerances
dvode_state.rtol_dens = tol_dens;
dvode_state.rtol_spec = tol_spec;
dvode_state.rtol_enuc = tol_ener;

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

int istate = dvode(burn_state, dvode_state);

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

// update the full U_new
U_new[URHO] = dvode_state.y(1);
for (int n = 0; n < NumSpec; ++n) {
U_new[UFS+n] = dvode_state.y(2+n);
}
U_new[UEINT] = dvode_state.y(NumSpec+2);
auto v2 = 0.0_rt;
for (int n = 0; n < 3; ++n) {
v2 += U_new[UMX+n] * U_new[UMX+n];
}
U_new[UEDEN] = U_new[UEINT] + 0.5_rt * v2 / U_new[URHO];

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

#endif
}

#endif

0 comments on commit fa367de

Please sign in to comment.