From 593f02b0d0c17838ba1dbc0d780a42c185e8d9c2 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sun, 31 Jul 2022 14:07:56 -0400 Subject: [PATCH] Simple sdc analytic jac memory (#955) this is similar to the change we did for the jacobian = 3 version in #917 --- integration/VODE/Make.package | 2 - .../VODE/vode_jac_util_simplified_sdc.H | 9 - .../VODE/vode_jac_util_simplified_sdc.cpp | 238 ------------------ integration/VODE/vode_rhs_simplified_sdc.H | 23 +- 4 files changed, 13 insertions(+), 259 deletions(-) delete mode 100644 integration/VODE/vode_jac_util_simplified_sdc.H delete mode 100644 integration/VODE/vode_jac_util_simplified_sdc.cpp diff --git a/integration/VODE/Make.package b/integration/VODE/Make.package index 37953dfcda..762e8feb02 100644 --- a/integration/VODE/Make.package +++ b/integration/VODE/Make.package @@ -2,8 +2,6 @@ ifeq ($(USE_SIMPLIFIED_SDC), TRUE) CEXE_headers += actual_integrator_simplified_sdc.H CEXE_headers += vode_rhs_simplified_sdc.H CEXE_headers += vode_type_simplified_sdc.H - CEXE_headers += vode_jac_util_simplified_sdc.H - CEXE_sources += vode_jac_util_simplified_sdc.cpp else ifneq ($(USE_TRUE_SDC), TRUE) CEXE_headers += vode_rhs.H diff --git a/integration/VODE/vode_jac_util_simplified_sdc.H b/integration/VODE/vode_jac_util_simplified_sdc.H deleted file mode 100644 index eb3fb82e85..0000000000 --- a/integration/VODE/vode_jac_util_simplified_sdc.H +++ /dev/null @@ -1,9 +0,0 @@ -#ifndef VODE_JAC_UTIL_SIMPLIFIED_SDC_H -#define VODE_JAC_UTIL_SIMPLIFIED_SDC_H - -AMREX_GPU_HOST_DEVICE -void jac_to_vode(const Real time, burn_t& state, - JacNetArray2D& jac_react, dvode_t& vode_state, - RArray2D& jac); - -#endif diff --git a/integration/VODE/vode_jac_util_simplified_sdc.cpp b/integration/VODE/vode_jac_util_simplified_sdc.cpp deleted file mode 100644 index af0402254c..0000000000 --- a/integration/VODE/vode_jac_util_simplified_sdc.cpp +++ /dev/null @@ -1,238 +0,0 @@ -#include -#include -#include -#include -#include - -AMREX_GPU_HOST_DEVICE -void jac_to_vode(const Real time, burn_t& state, - JacNetArray2D& jac_react, dvode_t& vode_state, - RArray2D& jac) -{ - - // this is only used with an analytic Jacobian. At the moment, we - // only support a dense Jacobian. - - - // we come in with jac_react being the Jacobian of the reacting - // system but we need to convert it to the SDC system. jac_react - // has the derivatives with respect to the native network - // variables, X, e. - -#if defined(SDC_EVOLVE_ENERGY) - - // The system we integrate has the form (rho X_k, rho E, rho e) - - // dR/dw has the form: - // - // SFS / d(rho X1dot)/dX1 d(rho X1dit)/dX2 ... d(rho X1dot)/dT \ - // | d(rho X2dot)/dX1 d(rho X2dot)/dX2 ... d(rho X2dot)/dT | - // SFS-1+nspec | ... | - // SEINT \ d(rho Edot)/dX1 d(rho Edot)/dX2 ... d(rho Edot)/dT / - // - // SFS SEINT - - - // keep in mind here that that we are using 1-based indexing but SFS, ... are 0-based. - - // start with the d(omegadot)/d* rows - - for (int m = 1; m <= NumSpec; m++) { - - // fill all elements of row m - - // first the species derivatives - - for (int n = 1; n <= NumSpec; n++) { - jac(SFS+m, SFS+n) = jac_react(m, n); - } - - // now the energy derivatives - - jac(SFS+m, SEINT+1) = jac_react(m, net_ienuc); - } - - // now the energy rows - - for (int n = 1; n <= NumSpec; n++) { - jac(SEINT+1, SFS+n) = jac_react(net_ienuc, n); - } - - jac(SEINT+1, SEINT+1) = jac_react(net_ienuc, net_ienuc); - -#elif defined(SDC_EVOLVE_ENTHALPY) - - constexpr int iwrho = 1; - constexpr int iwfs=2; - constexpr int iwe = iwfs+NumSpec; - constexpr int iwvar = 2+NumSpec; - - amrex::Array2D dRdw = {0.0_rt}; - amrex::Array2D dwdU = {0.0_rt}; - - constexpr Real eps = 1.e-8_rt; - - // this is 0-based to be consistent with SFS, ... - constexpr int SRHO_EXTRA = SVAR_EVOLVE; - - // jac_react 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. - - // also fill the ydot - YdotNetArray1D ydot; - vode_to_burn(time, vode_state, state); - actual_rhs(state, ydot); - - // 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 = 1; n <= NumSpec; n++) { - ydot(n) = ydot(n) * aion[n-1]; - } - - // now perturb density and call the RHS to compute the derivative wrt rho - // species rates come back in terms of molar fractions - burn_t state_pert = state; - state_pert.rho = state.rho * (1.0_rt + eps); - - YdotNetArray1D ydot_pert; - actual_rhs(state_pert, ydot_pert); - - // make the rates dX/dt and not dY/dt - for (int n = 1; n <= NumSpec; n++) { - ydot_pert(n) = ydot_pert(n) * aion[n-1]; - } - - - // Our R source has components for species and enthalpy only. But - // we will extend it here to include the mass density too to ensure - // that we have a square matrix in dU/dw that we can take the - // inverse of to use below. When we compute the final Jacobian, we will - // discard the density row. - - // Our jacobian, dR/dw has the form: - // - // SFS / 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 | - // SFS-1+nspec | ... | - // SENTH | d(rho h)/drho d(rho h)/dX1 d(rho h)/dX2 ... d(rho h)/dT | - // SRHO_EXTRA \ 0 0 0 0 / - - - // 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(SFS+m, iwrho) = ydot(m) + state.rho * (ydot_pert(m) - ydot(m))/(eps * state.rho); - } - - // d( d(rho h)/dt)/drho - dRdw(SENTH+1, iwrho) = ydot(net_ienuc) + - state.rho * (ydot_pert(net_ienuc) - ydot(net_ienuc))/(eps * state.rho); - - // d( d(rho)/dt)/drho - dRdw(SRHO_EXTRA+1, iwrho) = 0.0_rt; - - // fill the columns of dRdw corresponding to each derivative - // with respect to species mass fraction - for (int n = 1; n <= NumSpec; n++) { - for (int m = 1; m <= NumSpec; m++) { - // d( d(rho X_m)/dt)/dX_n - dRdw(SFS+m, iwfs-1+n) = state.rho * jac_react(m, n); - } - - // d( d(rho h)/dt)/dX_n - dRdw(SENTH+1, iwfs-1+n) = state.rho * jac_react(net_ienuc, n); - - // d( d(rho)/dt)/dX_n - dRdw(SRHO_EXTRA+1, iwfs-1+n) = 0.0_rt; - } - - // now fill the column corresponding to derivatives with respect to - // internal energy -- this column is iwe - - // d( d(rho X_m)/dt)/de - for (int m = 1; m <= NumSpec; m++) { - dRdw(SFS+m, iwe) = state.rho * jac_react(m, net_ienuc); - } - - // d( d(rho h)/dt)/de - dRdw(SENTH+1, iwe) = state.rho * jac_react(net_ienuc, net_ienuc); - - // d( d(rho)/dt)/de - dRdw(SRHO_EXTRA+1, iwe) = 0.0_rt; - - // that completes dRdw - - // construct dwdU. Here we take U = (rho X, rho h, rho)^T - - // density row (iwrho) - dwdU(iwrho, SRHO_EXTRA+1) = 1.0_rt; - - // species rows - for (int m = 1; m <= NumSpec; m++) { - dwdU(iwfs-1+m, SFS+m) = 1.0_rt / state.rho; - dwdU(iwfs-1+m, SRHO_EXTRA+1) = -state.xn[m-1] / state.rho; - } - - // h row - eos_t eos_state; - eos_state.rho = state.rho; - eos_state.T = 1.e4_rt; // initial guess - for (int n = 0; n < NumSpec; n++) { - eos_state.xn[n] = vode_state.y(SFS+1+n) / state.rho; - } -#ifdef NSE_THERMO - set_nse_aux_from_X(eos_state); -#endif - - eos_state.h = vode_state.y(SENTH+1) / state.rho; - - eos(eos_input_rh, eos_state); - - const eos_xderivs_t eos_xderivs = composition_derivatives(eos_state); - - // internal energy row - for (int n = 1; n <= NumSpec; n++) { - dwdU(iwe, SFS+n) = -eos_xderivs.dhdX[n-1] / (state.rho * eos_state.dedT); - } - dwdU(iwe, SENTH+1) = 1.0_rt / (state.rho * eos_state.dhdT); - Real X_dhdX_sum = 0.0; - for (int n = 0; n < NumSpec; n++) { - X_dhdX_sum += eos_state.xn[n] * eos_xderivs.dhdX[n]; - } - dwdU(iwe, SRHO_EXTRA+1) = - (X_dhdX_sum - state.rho * eos_state.dhdr - eos_state.h) / - (state.rho * eos_state.dhdT); - - - - // compute J = dR/dw dw/dU - - // J is SVAR_EVOLVE x SVAR_EVOLVE, which will call m x n - // - // J = dR/dw dw/dU - // - // dR/dw is SVAR_EVOLVE+1 x iwvar, which we call m x k - // dw/dU is iwvar x SVAR_EVOLVE+1, which we call k x n - // - - // we need to cut out the density (SRHO_EXTRA) row and column of - // the Jacobian, since that is not in our full SVAR_EVOLVE state - for (int n = 1; n <= SVAR_EVOLVE; n++) { - if (n == SRHO_EXTRA+1) continue; - for (int m = 1; m <= SVAR_EVOLVE; m++) { - if (m == SRHO_EXTRA+1) continue; - - jac(m, n) = 0.0_rt; - for (int k = 1; k <= iwvar; k++) { - jac(m, n) += dRdw(m, k) * dwdU(k, n); - } - } - } - -#endif - - -} diff --git a/integration/VODE/vode_rhs_simplified_sdc.H b/integration/VODE/vode_rhs_simplified_sdc.H index 827b5a987d..58cea60a9f 100644 --- a/integration/VODE/vode_rhs_simplified_sdc.H +++ b/integration/VODE/vode_rhs_simplified_sdc.H @@ -6,7 +6,6 @@ #include #include -#include #include #include #include @@ -143,36 +142,40 @@ void jac (burn_t& state, dvode_t& vode_state, MatrixType& pd) // Call the specific network routine to get the Jacobian. - JacNetArray2D jac_react; - - actual_jac(state, jac_react); + actual_jac(state, 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++) { - jac_react(n,m) = jac_react(n,m) * aion[n-1]; + pd(n,m) = pd(n,m) * aion[n-1]; } } for (int m = 1; m <= neqs; m++) { for (int n = 1; n <= NumSpec; n++) { - jac_react(m,n) = jac_react(m,n) * aion_inv[n-1]; + pd(m,n) = pd(m,n) * aion_inv[n-1]; } } // apply fudge factor: if (react_boost > 0.0_rt) { - jac_react.mul(react_boost); + pd.mul(react_boost); } - // convert to the system we are using + // The system we integrate has the form (rho X_k, rho e) - jac_to_vode(vode_state.tn, state, jac_react, vode_state, pd); + // pd is now of the form: + // + // SFS / d(rho X1dot)/dX1 d(rho X1dit)/dX2 ... 1/cv d(rho X1dot)/dT \ + // | d(rho X2dot)/dX1 d(rho X2dot)/dX2 ... 1/cv d(rho X2dot)/dT | + // SFS-1+nspec | ... | + // SEINT \ d(rho Edot)/dX1 d(rho Edot)/dX2 ... 1/cv d(rho Edot)/dT / + // + // SFS SEINT } - } #endif