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

reorder some loops over Jacobians #1639

Merged
merged 3 commits into from
Jul 31, 2024
Merged

reorder some loops over Jacobians #1639

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1af7a78
reorder some loops over Jacobians
zingale Jul 31, 2024
b72771f
fix compilation
zingale Jul 31, 2024
6b8b917
fix compilation
zingale Jul 31, 2024
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
28 changes: 11 additions & 17 deletions integration/integrator_rhs_sdc.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -105,15 +105,8 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd)


if (state.T <= EOSData::mintemp || state.T >= integrator_rp::MAX_TEMP) {

for (int j = 1; j <= INT_NEQS; ++j) {
for (int i = 1; i <= INT_NEQS; ++i) {
pd(i,j) = 0.0_rt;
}
}

pd.zero();
return;

}

// Call the specific network routine to get the Jacobian.
Expand All @@ -122,15 +115,16 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd)

// 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++) {
pd(n,m) = pd(n,m) * aion[n-1];

for (int jcol = 1; jcol <= neqs; ++jcol) {
for (int irow = 1; irow <= NumSpec; irow++) {
pd(irow, jcol) = pd(irow, jcol) * aion[irow-1];
}
}

for (int m = 1; m <= neqs; m++) {
for (int n = 1; n <= NumSpec; n++) {
pd(m,n) = pd(m,n) * aion_inv[n-1];
for (int jcol = 1; jcol <= NumSpec; ++jcol) {
for (int irow = 1; irow <= neqs; ++irow) {
pd(irow, jcol) = pd(irow, jcol) * aion_inv[jcol-1];
}
}

Expand Down Expand Up @@ -170,9 +164,9 @@ void jac (const amrex::Real time, BurnT& state, T& int_state, MatrixType& pd)

eos_xderivs_t eos_xderivs = composition_derivatives(eos_state);

for (int m = 1; m <= neqs; m++) {
for (int n = 1; n <= NumSpec; n++) {
pd(m, n) -= eos_xderivs.dedX[n-1] * pd(m, net_ienuc);
for (int jcol = 1; jcol <= NumSpec;++jcol) {
for (int irow = 1; irow <= neqs; ++irow) {
pd(irow, jcol) -= eos_xderivs.dedX[jcol-1] * pd(irow, net_ienuc);
}
}

Expand Down
9 changes: 1 addition & 8 deletions integration/integrator_rhs_strang.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,15 +119,8 @@ void jac ([[maybe_unused]] const amrex::Real time, BurnT& state, T& int_state, M
// bounds. Otherwise set the Jacobian to zero and return.

if (state.T <= EOSData::mintemp || state.T >= integrator_rp::MAX_TEMP) {

for (int j = 1; j <= INT_NEQS; ++j) {
for (int i = 1; i <= INT_NEQS; ++i) {
pd(i,j) = 0.0_rt;
}
}

pd.zero();
return;

}

// Call the specific network routine to get the Jacobian.
Expand Down
13 changes: 3 additions & 10 deletions integration/utils/numerical_jacobian.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,15 +136,8 @@ void numerical_jac(BurnT& state, const jac_info_t& jac_info, JacNetArray2D& jac)
state_delp.T += dy;

if (state_delp.T <= EOSData::mintemp || state_delp.T >= integrator_rp::MAX_TEMP) {

for (int i = 1; i <= int_neqs; i++) {
for (int j = 1; j <= int_neqs; j++) {
jac(i,j) = 0.0_rt;
}
}

jac.zero();
return;

}


Expand Down Expand Up @@ -195,8 +188,8 @@ void numerical_jac(BurnT& state, const jac_info_t& jac_info, JacNetArray2D& jac)
// now correct the species derivatives
// this constructs dy/dX_k |_e = dy/dX_k |_T - e_{X_k} |_T dy/dT / c_v

for (int m = 1; m <= int_neqs; m++) {
for (int n = 1; n <= NumSpec; n++) {
for (int n = 1; n <= NumSpec; n++) {
for (int m = 1; m <= int_neqs; m++) {
jac(m, n) -= eos_xderivs.dedX[n-1] * jac(m, net_ienuc);
}
}
Expand Down
Loading