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

move the true SDC VODE into the VODE header #2567

Merged
merged 1 commit into from
Sep 21, 2023
Merged

move the true SDC VODE into the VODE header #2567

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
140 changes: 1 addition & 139 deletions Source/sdc/Castro_sdc_util.H
View file
Open in desktop
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
Open in desktop
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