diff --git a/Source/reactions/Castro_react_util.H b/Source/reactions/Castro_react_util.H index bca29ac27c..f64bfc6836 100644 --- a/Source/reactions/Castro_react_util.H +++ b/Source/reactions/Castro_react_util.H @@ -53,212 +53,4 @@ okay_to_burn_type(burn_t const& state) { return true; } -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -single_zone_react_source(GpuArray const& state, - GpuArray& R, - burn_t& burn_state) { - - // note: we want this burn_state to come out of here so we can - // reuse it elsewhere - - auto rhoInv = 1.0_rt / state[URHO]; - - burn_state.rho = state[URHO]; - burn_state.T = state[UTEMP]; - burn_state.e = state[UEINT] * rhoInv; - - for (int n = 0; n < NumSpec; ++n) { - burn_state.xn[n] = amrex::max(amrex::min(state[UFS+n] * rhoInv, 1.0_rt), small_x); - } - -#if NAUX_NET > 0 - for (int n = 0; n < NumAux; ++n) { - burn_state.aux[n] = state[UFX+n] * rhoInv; - } -#endif - - eos_t eos_state; - - // Ensure that the temperature going in is consistent with the internal energy. - burn_to_eos(burn_state, eos_state); - eos(eos_input_re, eos_state); - eos_to_burn(eos_state, burn_state); - - // eos_get_small_temp(&small_temp); - burn_state.T = amrex::min(MAX_TEMP, amrex::max(burn_state.T, small_temp)); - - Array1D ydot; - - actual_rhs(burn_state, ydot); - - // store the instantaneous R - for (int n = 0; n < NUM_STATE; ++n) { - R[n] = 0.0_rt; - } - - // species rates come back in terms of molar fractions - for (int n = 0; n < NumSpec; ++n) { - R[UFS+n] = state[URHO] * aion[n] * ydot(1+n); - } - - R[UEDEN] = state[URHO] * ydot(net_ienuc); - R[UEINT] = state[URHO] * ydot(net_ienuc); -} - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -single_zone_react_source(const int i, const int j, const int k, - Array4 const& state, - Array4 const& R, - burn_t& burn_state) { - GpuArray state_arr; - GpuArray R_arr; - - for (int n = 0; n < NUM_STATE; ++n) { - state_arr[n] = state(i,j,k,n); - R_arr[n] = R(i,j,k,n); - } - - single_zone_react_source(state_arr, R_arr, burn_state); - - for (int n = 0; n < NUM_STATE; ++n) { - R(i,j,k,n) = R_arr[n]; - } -} - -AMREX_GPU_HOST_DEVICE AMREX_INLINE -void -single_zone_jac(GpuArray const& state, - burn_t& burn_state, const Real dt, - Array2D& dRdw) { - -#ifdef SIMPLIFIED_SDC -#ifndef AMREX_USE_GPU - amrex::Error("we shouldn't be here with the simplified SDC method (USE_SIMPLIFIED_SDC=TRUE)"); -#endif -#else - - Array1D ydot; - Array1D ydot_pert; - JacNetArray2D Jac; - - const auto eps = 1.e-8_rt; - - actual_rhs(burn_state, ydot); - - // Jac has the derivatives with respect to the native - // network variables, X, T. e. It does not have derivatives with - // respect to density, so we'll have to compute those ourselves. - - if (sdc_newton_use_analytic_jac == 0) { - // note the numerical Jacobian will be returned in terms of X - jac_info_t jac_info; - jac_info.h = dt; - numerical_jac(burn_state, jac_info, Jac); - } else { - actual_jac(burn_state, Jac); - - // The Jacobian from the nets is in terms of dYdot/dY, but we want - // it was dXdot/dX, so convert here. - for (int n = 1; n <= NumSpec; n++) { - for (int m = 1; m <= neqs; m++) { - // d(X_n)/dq_m - Jac(n,m) = Jac(n,m) * aion[n-1]; - } - } - - for (int m = 1; m <= neqs; m++) { - for (int n = 1; n <= NumSpec; n++) { - // d RHS_m / dX_n - Jac(m,n) = Jac(m,n) * aion_inv[n-1]; - } - } - } -#endif - - // at this point, our Jacobian should be entirely in terms of X, - // not Y. Let's now fix the rhs terms themselves to be in terms of - // dX/dt and not dY/dt. - for (int n = 0; n < NumSpec; ++n) { - ydot(n+1) *= aion[n]; - } - - // Our jacobian, dR/dw has the form: - // - // / 0 0 0 0 \ - // | d(rho X1dot)/drho d(rho X1dot)/dX1 d(rho X1dit)/dX2 ... d(rho X1dot)/dT | - // | d(rho X2dot)/drho d(rho X2dot)/dX1 d(rho X2dot)/dX2 ... d(rho X2dot)/dT | - // | ... | - // \ d(rho Edot)/drho d(rho Edot)/dX1 d(rho Edot)/dX2 ... d(rho Edot)/dT / - - for (int n = 0; n < NumSpec+2; ++n) { - for (int m = 0; m < NumSpec+2; ++m) { - dRdw(n,m) = 0.0_rt; - } - } - - // now perturb density and call the RHS to compute the derivative wrt rho - // species rates come back in terms of molar fractions - burn_t burn_state_pert = burn_state; - - burn_state_pert.rho = burn_state.rho * (1.0_rt + eps); - burn_state_pert.T = burn_state.T; - burn_state_pert.e = burn_state.e; - for (int n = 0; n < NumSpec; ++n) { - burn_state_pert.xn[n] = burn_state.xn[n]; - } -#if NAUX_NET > 0 - for (int n = 0; n < NumAux; ++n) { - burn_state_pert.aux[n] = burn_state.aux[n]; - } -#endif - - actual_rhs(burn_state_pert, ydot_pert); - - // make the rates dX/dt and not dY/dt - for (int n = 0; n < NumSpec; ++n) { - ydot_pert(n+1) *= aion[n]; - } - - // fill the column of dRdw corresponding to the derivative - // with respect to rho - for (int m = 1; m <= NumSpec; ++m) { - // d( d(rho X_m)/dt)/drho - dRdw(m, iwrho) = ydot(m) + - state[URHO] * (ydot_pert(m) - ydot(m))/(eps * burn_state.rho); - } - - // d( d(rho E)/dt)/drho - dRdw(NumSpec+1, iwrho) = ydot(net_ienuc) + - state[URHO] * (ydot_pert(net_ienuc) - ydot(net_ienuc))/(eps * burn_state.rho); - - // fill the columns of dRdw corresponding to each derivative - // with respect to species mass fraction - for (int n = 1; n <= NumSpec; ++n) { - dRdw(0, iwfs-1+n) = 0.0_rt; // density source - - for (int m = 1; m <= NumSpec; ++m) { - // d( d(rho X_m)/dt)/dX_n - dRdw(m, iwfs-1+n) = state[URHO] * Jac(m, n); - } - - // d( d(rho E)/dt)/dX_n - dRdw(NumSpec+1, iwfs-1+n) = state[URHO] * Jac(net_ienuc, n); - - } - - // now fill the column corresponding to derivatives with respect to - // energy -- this column is iwe - dRdw(0, iwe) = 0.0_rt; - - // d( d(rho X_m)/dt)/de - for (int m = 1; m <= NumSpec; ++m) { - dRdw(m, iwe) = state[URHO] * Jac(m, net_ienuc); - } - - // d( d(rho E)/dt)/de - dRdw(NumSpec+1, iwe) = state[URHO] * Jac(net_ienuc, net_ienuc); -} - #endif diff --git a/Source/sdc/Castro_sdc_util.H b/Source/sdc/Castro_sdc_util.H index ef632056b5..d223d759a7 100644 --- a/Source/sdc/Castro_sdc_util.H +++ b/Source/sdc/Castro_sdc_util.H @@ -806,6 +806,7 @@ instantaneous_react(const int i, const int j, const int k, } } } + #endif #endif diff --git a/Source/sdc/Make.package b/Source/sdc/Make.package index 78cbb732ac..709b42db9d 100644 --- a/Source/sdc/Make.package +++ b/Source/sdc/Make.package @@ -7,5 +7,6 @@ ifneq ($(USE_GPU), TRUE) CEXE_headers += Castro_sdc_util.H ifeq ($(USE_REACT), TRUE) CEXE_headers += vode_rhs_true_sdc.H + CEXE_headers += sdc_react_util.H endif endif diff --git a/Source/sdc/sdc_react_util.H b/Source/sdc/sdc_react_util.H new file mode 100644 index 0000000000..3357612879 --- /dev/null +++ b/Source/sdc/sdc_react_util.H @@ -0,0 +1,212 @@ +#ifndef SDC_REACT_UTIL_H +#define SDC_REACT_UTIL_H + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void +single_zone_react_source(GpuArray const& state, + GpuArray& R, + burn_t& burn_state) { + + // note: we want this burn_state to come out of here so we can + // reuse it elsewhere + + auto rhoInv = 1.0_rt / state[URHO]; + + burn_state.rho = state[URHO]; + burn_state.T = state[UTEMP]; + burn_state.e = state[UEINT] * rhoInv; + + for (int n = 0; n < NumSpec; ++n) { + burn_state.xn[n] = amrex::max(amrex::min(state[UFS+n] * rhoInv, 1.0_rt), small_x); + } + +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; ++n) { + burn_state.aux[n] = state[UFX+n] * rhoInv; + } +#endif + + eos_t eos_state; + + // Ensure that the temperature going in is consistent with the internal energy. + burn_to_eos(burn_state, eos_state); + eos(eos_input_re, eos_state); + eos_to_burn(eos_state, burn_state); + + // eos_get_small_temp(&small_temp); + burn_state.T = amrex::min(MAX_TEMP, amrex::max(burn_state.T, small_temp)); + + Array1D ydot; + + actual_rhs(burn_state, ydot); + + // store the instantaneous R + for (int n = 0; n < NUM_STATE; ++n) { + R[n] = 0.0_rt; + } + + // species rates come back in terms of molar fractions + for (int n = 0; n < NumSpec; ++n) { + R[UFS+n] = state[URHO] * aion[n] * ydot(1+n); + } + + R[UEDEN] = state[URHO] * ydot(net_ienuc); + R[UEINT] = state[URHO] * ydot(net_ienuc); +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void +single_zone_react_source(const int i, const int j, const int k, + Array4 const& state, + Array4 const& R, + burn_t& burn_state) { + GpuArray state_arr; + GpuArray R_arr; + + for (int n = 0; n < NUM_STATE; ++n) { + state_arr[n] = state(i,j,k,n); + R_arr[n] = R(i,j,k,n); + } + + single_zone_react_source(state_arr, R_arr, burn_state); + + for (int n = 0; n < NUM_STATE; ++n) { + R(i,j,k,n) = R_arr[n]; + } +} + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void +single_zone_jac(GpuArray const& state, + burn_t& burn_state, const Real dt, + Array2D& dRdw) { + +#ifdef SIMPLIFIED_SDC +#ifndef AMREX_USE_GPU + amrex::Error("we shouldn't be here with the simplified SDC method (USE_SIMPLIFIED_SDC=TRUE)"); +#endif +#else + + Array1D ydot; + Array1D ydot_pert; + JacNetArray2D Jac; + + const auto eps = 1.e-8_rt; + + actual_rhs(burn_state, ydot); + + // Jac has the derivatives with respect to the native + // network variables, X, T. e. It does not have derivatives with + // respect to density, so we'll have to compute those ourselves. + + if (sdc_newton_use_analytic_jac == 0) { + // note the numerical Jacobian will be returned in terms of X + jac_info_t jac_info; + jac_info.h = dt; + numerical_jac(burn_state, jac_info, Jac); + } else { + actual_jac(burn_state, Jac); + + // The Jacobian from the nets is in terms of dYdot/dY, but we want + // it was dXdot/dX, so convert here. + for (int n = 1; n <= NumSpec; n++) { + for (int m = 1; m <= neqs; m++) { + // d(X_n)/dq_m + Jac(n,m) = Jac(n,m) * aion[n-1]; + } + } + + for (int m = 1; m <= neqs; m++) { + for (int n = 1; n <= NumSpec; n++) { + // d RHS_m / dX_n + Jac(m,n) = Jac(m,n) * aion_inv[n-1]; + } + } + } +#endif + + // at this point, our Jacobian should be entirely in terms of X, + // not Y. Let's now fix the rhs terms themselves to be in terms of + // dX/dt and not dY/dt. + for (int n = 0; n < NumSpec; ++n) { + ydot(n+1) *= aion[n]; + } + + // Our jacobian, dR/dw has the form: + // + // / 0 0 0 0 \ + // | d(rho X1dot)/drho d(rho X1dot)/dX1 d(rho X1dit)/dX2 ... d(rho X1dot)/dT | + // | d(rho X2dot)/drho d(rho X2dot)/dX1 d(rho X2dot)/dX2 ... d(rho X2dot)/dT | + // | ... | + // \ d(rho Edot)/drho d(rho Edot)/dX1 d(rho Edot)/dX2 ... d(rho Edot)/dT / + + for (int n = 0; n < NumSpec+2; ++n) { + for (int m = 0; m < NumSpec+2; ++m) { + dRdw(n,m) = 0.0_rt; + } + } + + // now perturb density and call the RHS to compute the derivative wrt rho + // species rates come back in terms of molar fractions + burn_t burn_state_pert = burn_state; + + burn_state_pert.rho = burn_state.rho * (1.0_rt + eps); + burn_state_pert.T = burn_state.T; + burn_state_pert.e = burn_state.e; + for (int n = 0; n < NumSpec; ++n) { + burn_state_pert.xn[n] = burn_state.xn[n]; + } +#if NAUX_NET > 0 + for (int n = 0; n < NumAux; ++n) { + burn_state_pert.aux[n] = burn_state.aux[n]; + } +#endif + + actual_rhs(burn_state_pert, ydot_pert); + + // make the rates dX/dt and not dY/dt + for (int n = 0; n < NumSpec; ++n) { + ydot_pert(n+1) *= aion[n]; + } + + // fill the column of dRdw corresponding to the derivative + // with respect to rho + for (int m = 1; m <= NumSpec; ++m) { + // d( d(rho X_m)/dt)/drho + dRdw(m, iwrho) = ydot(m) + + state[URHO] * (ydot_pert(m) - ydot(m))/(eps * burn_state.rho); + } + + // d( d(rho E)/dt)/drho + dRdw(NumSpec+1, iwrho) = ydot(net_ienuc) + + state[URHO] * (ydot_pert(net_ienuc) - ydot(net_ienuc))/(eps * burn_state.rho); + + // fill the columns of dRdw corresponding to each derivative + // with respect to species mass fraction + for (int n = 1; n <= NumSpec; ++n) { + dRdw(0, iwfs-1+n) = 0.0_rt; // density source + + for (int m = 1; m <= NumSpec; ++m) { + // d( d(rho X_m)/dt)/dX_n + dRdw(m, iwfs-1+n) = state[URHO] * Jac(m, n); + } + + // d( d(rho E)/dt)/dX_n + dRdw(NumSpec+1, iwfs-1+n) = state[URHO] * Jac(net_ienuc, n); + + } + + // now fill the column corresponding to derivatives with respect to + // energy -- this column is iwe + dRdw(0, iwe) = 0.0_rt; + + // d( d(rho X_m)/dt)/de + for (int m = 1; m <= NumSpec; ++m) { + dRdw(m, iwe) = state[URHO] * Jac(m, net_ienuc); + } + + // d( d(rho E)/dt)/de + dRdw(NumSpec+1, iwe) = state[URHO] * Jac(net_ienuc, net_ienuc); +} + +#endif diff --git a/Source/sdc/vode_rhs_true_sdc.H b/Source/sdc/vode_rhs_true_sdc.H index a01b11b9be..bd12ae4e2c 100644 --- a/Source/sdc/vode_rhs_true_sdc.H +++ b/Source/sdc/vode_rhs_true_sdc.H @@ -12,6 +12,7 @@ #include #endif #include +#include // 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