Skip to content

Commit

Permalink
move SDC-specific routiens into Source/sdc (#2558)
Browse files Browse the repository at this point in the history
  • Loading branch information
@zingale
zingale authored Sep 19, 2023
1 parent 575c8c9 commit 07132a8
Show file tree
Hide file tree
Showing 5 changed files with 215 additions and 208 deletions.
208 changes: 0 additions & 208 deletions Source/reactions/Castro_react_util.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real, NUM_STATE> const& state,
GpuArray<Real, NUM_STATE>& 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<Real, 1, neqs> 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 Real> const& state,
Array4<Real> const& R,
burn_t& burn_state) {
GpuArray<Real, NUM_STATE> state_arr;
GpuArray<Real, NUM_STATE> 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<Real, NUM_STATE> const& state,
burn_t& burn_state, const Real dt,
Array2D<Real, 0, NumSpec+1, 0, NumSpec+1>& 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<Real, 1, neqs> ydot;
Array1D<Real, 1, neqs> 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
1 change: 1 addition & 0 deletions Source/sdc/Castro_sdc_util.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -806,6 +806,7 @@ instantaneous_react(const int i, const int j, const int k,
}
}
}

#endif

#endif
1 change: 1 addition & 0 deletions Source/sdc/Make.package
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading

0 comments on commit 07132a8

Please sign in to comment.