diff --git a/Source/sdc/Castro_sdc_util.H b/Source/sdc/Castro_sdc_util.H index 09581b06e8..164daff206 100644 --- a/Source/sdc/Castro_sdc_util.H +++ b/Source/sdc/Castro_sdc_util.H @@ -8,10 +8,7 @@ #endif #include #include -#include -#ifndef NEW_NETWORK_IMPLEMENTATION -#include -#endif +#include #endif // solvers @@ -19,9 +16,6 @@ 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, @@ -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 const& U_old, - GpuArray& U_new, - GpuArray 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 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 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, diff --git a/Source/sdc/vode_rhs_true_sdc.H b/Source/sdc/vode_rhs_true_sdc.H index bd12ae4e2c..91790cc4fc 100644 --- a/Source/sdc/vode_rhs_true_sdc.H +++ b/Source/sdc/vode_rhs_true_sdc.H @@ -14,6 +14,11 @@ #include #include +#include +#ifndef NEW_NETWORK_IMPLEMENTATION +#include +#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. @@ -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 const& U_old, + GpuArray& U_new, + GpuArray 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 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 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