diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index a43ef972d..61ef46370 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -145,6 +145,12 @@ public: virtual size_t componentIndex(const std::string& nm) const; protected: + //! Set reaction rate multipliers based on the sensitivity variables in + //! *params*. + virtual void applySensitivity(double* params); + //! Reset the reaction rate multipliers + virtual void resetSensitivity(double* params); + //! Pointer to the homogeneous Kinetics object that handles the reactions Kinetics* m_kin; diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 6d8f4d367..e1d0d5135 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -124,25 +124,7 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y, m_thermo->restoreState(m_state); Kinetics* kin; - size_t npar, ploc; - double mult; - - // process sensitivity parameters - if (params) { - - npar = m_pnum.size(); - for (size_t n = 0; n < npar; n++) { - mult = m_kin->multiplier(m_pnum[n]); - m_kin->setMultiplier(m_pnum[n], mult*params[n]); - } - ploc = npar; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_nsens_wall[m] > 0) { - m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc); - ploc += m_nsens_wall[m]; - } - } - } + applySensitivity(params); m_Q = 0.0; @@ -236,20 +218,7 @@ void ConstPressureReactor::evalEqs(doublereal time, doublereal* y, } // reset sensitivity parameters - if (params) { - npar = m_pnum.size(); - for (size_t n = 0; n < npar; n++) { - mult = m_kin->multiplier(m_pnum[n]); - m_kin->setMultiplier(m_pnum[n], mult/params[n]); - } - ploc = npar; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_nsens_wall[m] > 0) { - m_wall[m]->resetSensitivityParameters(m_lr[m]); - ploc += m_nsens_wall[m]; - } - } - } + resetSensitivity(params); } size_t ConstPressureReactor::componentIndex(const string& nm) const diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 59e7bbb59..e59739643 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -129,28 +129,9 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y, { size_t nk; m_thermo->restoreState(m_state); + applySensitivity(params); Kinetics* kin; - size_t npar, ploc; - double mult; - - // process sensitivity parameters - if (params) { - - npar = m_pnum.size(); - for (size_t n = 0; n < npar; n++) { - mult = m_kin->multiplier(m_pnum[n]); - m_kin->setMultiplier(m_pnum[n], mult*params[n]); - } - ploc = npar; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_nsens_wall[m] > 0) { - m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc); - ploc += m_nsens_wall[m]; - } - } - } - m_Q = 0.0; // compute wall terms @@ -245,21 +226,7 @@ void IdealGasConstPressureReactor::evalEqs(doublereal time, doublereal* y, ydot[1] = 0.0; } - // reset sensitivity parameters - if (params) { - npar = m_pnum.size(); - for (size_t n = 0; n < npar; n++) { - mult = m_kin->multiplier(m_pnum[n]); - m_kin->setMultiplier(m_pnum[n], mult/params[n]); - } - ploc = npar; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_nsens_wall[m] > 0) { - m_wall[m]->resetSensitivityParameters(m_lr[m]); - ploc += m_nsens_wall[m]; - } - } - } + resetSensitivity(params); } size_t IdealGasConstPressureReactor::componentIndex(const string& nm) const diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index 9d6e399f6..012798609 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -138,22 +138,7 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y, doublereal* ydot, doublereal* params) { m_thermo->restoreState(m_state); - - // process sensitivity parameters - if (params) { - size_t npar = m_pnum.size(); - for (size_t n = 0; n < npar; n++) { - double mult = m_kin->multiplier(m_pnum[n]); - m_kin->setMultiplier(m_pnum[n], mult*params[n]); - } - size_t ploc = npar; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_nsens_wall[m] > 0) { - m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc); - ploc += m_nsens_wall[m]; - } - } - } + applySensitivity(params); m_vdot = 0.0; m_Q = 0.0; @@ -261,21 +246,7 @@ void IdealGasReactor::evalEqs(doublereal time, doublereal* y, "ydot[" + int2str(i) + "] is not finite"); } - // reset sensitivity parameters - if (params) { - size_t npar = m_pnum.size(); - for (size_t n = 0; n < npar; n++) { - double mult = m_kin->multiplier(m_pnum[n]); - m_kin->setMultiplier(m_pnum[n], mult/params[n]); - } - size_t ploc = npar; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_nsens_wall[m] > 0) { - m_wall[m]->resetSensitivityParameters(m_lr[m]); - ploc += m_nsens_wall[m]; - } - } - } + resetSensitivity(params); } size_t IdealGasReactor::componentIndex(const string& nm) const diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index d56782ea2..73dbd5d97 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -180,22 +180,7 @@ void Reactor::evalEqs(doublereal time, doublereal* y, doublereal* ydot, doublereal* params) { m_thermo->restoreState(m_state); - - // process sensitivity parameters - if (params) { - size_t npar = m_pnum.size(); - for (size_t n = 0; n < npar; n++) { - double mult = m_kin->multiplier(m_pnum[n]); - m_kin->setMultiplier(m_pnum[n], mult*params[n]); - } - size_t ploc = npar; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_nsens_wall[m] > 0) { - m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc); - ploc += m_nsens_wall[m]; - } - } - } + applySensitivity(params); m_vdot = 0.0; m_Q = 0.0; @@ -306,21 +291,7 @@ void Reactor::evalEqs(doublereal time, doublereal* y, "ydot[" + int2str(i) + "] is not finite"); } - // reset sensitivity parameters - if (params) { - size_t npar = m_pnum.size(); - for (size_t n = 0; n < npar; n++) { - double mult = m_kin->multiplier(m_pnum[n]); - m_kin->setMultiplier(m_pnum[n], mult/params[n]); - } - size_t ploc = npar; - for (size_t m = 0; m < m_nwalls; m++) { - if (m_nsens_wall[m] > 0) { - m_wall[m]->resetSensitivityParameters(m_lr[m]); - ploc += m_nsens_wall[m]; - } - } - } + resetSensitivity(params); } void Reactor::addSensitivityReaction(size_t rxn) @@ -383,4 +354,43 @@ size_t Reactor::componentIndex(const string& nm) const return npos; } +void Reactor::applySensitivity(double* params) +{ + if (!params) { + return; + } + size_t npar = m_pnum.size(); + for (size_t n = 0; n < npar; n++) { + double mult = m_kin->multiplier(m_pnum[n]); + m_kin->setMultiplier(m_pnum[n], mult*params[n]); + } + size_t ploc = npar; + for (size_t m = 0; m < m_nwalls; m++) { + if (m_nsens_wall[m] > 0) { + m_wall[m]->setSensitivityParameters(m_lr[m], params + ploc); + ploc += m_nsens_wall[m]; + } + } + +} + +void Reactor::resetSensitivity(double* params) +{ + if (!params) { + return; + } + size_t npar = m_pnum.size(); + for (size_t n = 0; n < npar; n++) { + double mult = m_kin->multiplier(m_pnum[n]); + m_kin->setMultiplier(m_pnum[n], mult/params[n]); + } + size_t ploc = npar; + for (size_t m = 0; m < m_nwalls; m++) { + if (m_nsens_wall[m] > 0) { + m_wall[m]->resetSensitivityParameters(m_lr[m]); + ploc += m_nsens_wall[m]; + } + } +} + }