From cf13b318f44a4fa0d0b34a706c1985064200dfc7 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sun, 8 May 2016 17:59:34 -0400 Subject: [PATCH] [Reactor] Implement enthalpy of formation sensitivity analysis --- include/cantera/numerics/FuncEval.h | 3 ++ include/cantera/zeroD/Reactor.h | 4 ++ include/cantera/zeroD/ReactorBase.h | 6 +++ include/cantera/zeroD/ReactorNet.h | 22 ++++++---- interfaces/cython/cantera/_cantera.pxd | 1 + interfaces/cython/cantera/reactor.pyx | 8 ++++ .../cython/cantera/test/test_reactor.py | 42 +++++++++++-------- src/numerics/CVodesIntegrator.cpp | 13 +++++- src/zeroD/Reactor.cpp | 38 ++++++++++++++--- src/zeroD/ReactorNet.cpp | 18 +++++--- src/zeroD/Wall.cpp | 14 ++++--- 11 files changed, 126 insertions(+), 43 deletions(-) diff --git a/include/cantera/numerics/FuncEval.h b/include/cantera/numerics/FuncEval.h index a63bbc79e..bf4994c05 100644 --- a/include/cantera/numerics/FuncEval.h +++ b/include/cantera/numerics/FuncEval.h @@ -65,6 +65,9 @@ public: //! This is the array which is perturbed and passed back as the fourth //! argument to eval(). vector_fp m_sens_params; + + //! Scaling factors for each sensitivity parameter + vector_fp m_paramScales; }; } diff --git a/include/cantera/zeroD/Reactor.h b/include/cantera/zeroD/Reactor.h index 95c1d9a1c..cc2a43708 100644 --- a/include/cantera/zeroD/Reactor.h +++ b/include/cantera/zeroD/Reactor.h @@ -126,6 +126,10 @@ public: //! (in the homogeneous phase). virtual void addSensitivityReaction(size_t rxn); + //! Add a sensitivity parameter associated with the enthalpy formation of + //! species *k* (in the homogeneous phase) + virtual void addSensitivitySpeciesEnthalpy(size_t k); + //! Return the index in the solution vector for this reactor of the //! component named *nm*. Possible values for *nm* are "mass", "volume", //! "int_energy", the name of a homogeneous phase species, or the name of a diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 1543dc113..56ec9b1d6 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -21,11 +21,17 @@ const int ConstPressureReactorType = 4; const int IdealGasReactorType = 5; const int IdealGasConstPressureReactorType = 6; +enum class SensParameterType { + reaction, + enthalpy +}; + struct SensitivityParameter { size_t local; //!< local parameter index size_t global; //!< global parameter index double value; //!< nominal value of the parameter + SensParameterType type; //!< type of sensitivity parameter }; /** diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index b8de2251b..7dcb1bfab 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -118,15 +118,19 @@ public: //! Return the sensitivity of the *k*-th solution component with respect to //! the *p*-th sensitivity parameter. /*! - * The normalized sensitivity coefficient \f$ S_{ki} \f$ of solution - * variable \f$ y_k \f$ with respect to sensitivity parameter \f$ p_i \f$ - * is defined as: + * The sensitivity coefficient \f$ S_{ki} \f$ of solution variable \f$ y_k + * \f$ with respect to sensitivity parameter \f$ p_i \f$ is defined as: * - * \f[ S_{ki} = \frac{p_i}{y_k} \frac{\partial y_k}{\partial p_i} \f] + * \f[ S_{ki} = \frac{1}{y_k} \frac{\partial y_k}{\partial p_i} \f] * * For reaction sensitivities, the parameter is a multiplier on the forward * rate constant (and implicitly on the reverse rate constant for - * reversible reactions). + * reversible reactions) which has a nominal value of 1.0, and the + * sensitivity is nondimensional. + * + * For species enthalpy sensitivities, the parameter is a perturbation to + * the molar enthalpy of formation, such that the dimensions of the + * sensitivity are kmol/J. */ double sensitivity(size_t k, size_t p); @@ -171,11 +175,15 @@ public: size_t globalComponentIndex(const std::string& component, size_t reactor=0); //! Used by Reactor and Wall objects to register the addition of - //! sensitivity reactions so that the ReactorNet can keep track of the + //! sensitivity parameters so that the ReactorNet can keep track of the //! order in which sensitivity parameters are added. + //! @param value The nominal value of the parameter + //! @param scale A scaling factor to be applied to the sensitivity + //! coefficient //! @returns the index of this parameter in the vector of sensitivity //! parameters (global across all reactors) - size_t registerSensitivityReaction(const std::string& name); + size_t registerSensitivityParameter(const std::string& name, double value, + double scale); //! The name of the p-th sensitivity parameter added to this ReactorNet. const std::string& sensitivityParameterName(size_t p) { diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 455f4220b..ae2528fb9 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -491,6 +491,7 @@ cdef extern from "cantera/zeroD/Reactor.h": void getState(double*) void addSensitivityReaction(size_t) except + + void addSensitivitySpeciesEnthalpy(size_t) except + size_t nSensParams() diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index ecdeec3d8..e02ad9e99 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -221,6 +221,14 @@ cdef class Reactor(ReactorBase): """ self.reactor.addSensitivityReaction(m) + def add_sensitivity_species_enthalpy(self, k): + """ + Specifies that the sensitivity of the state variables with respect to + species *k* should be computed. The reactor must be part of a network + first. + """ + self.reactor.addSensitivitySpeciesEnthalpy(self.thermo.species_index(k)) + def component_index(self, name): """ Returns the index of the component named *name* in the system. This diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py index 61a827225..84d04ebfb 100644 --- a/interfaces/cython/cantera/test/test_reactor.py +++ b/interfaces/cython/cantera/test/test_reactor.py @@ -1011,12 +1011,18 @@ class TestReactorSensitivities(utilities.CanteraTest): # Single reactor, changing the order in which parameters are added gas = ct.Solution('h2o2.xml') - def setup(): + def setup(params): net = ct.ReactorNet() gas.TPX = 900, 101325, 'H2:0.1, OH:1e-7, O2:0.1, AR:1e-5' r = reactorClass(gas) net.add_reactor(r) + + for kind, p in params: + if kind == 'r': + r.add_sensitivity_reaction(p) + elif kind == 's': + r.add_sensitivity_species_enthalpy(p) return r, net def integrate(r, net): @@ -1024,26 +1030,28 @@ class TestReactorSensitivities(utilities.CanteraTest): net.step() return net.sensitivities() - r1,net1 = setup() - params1 = [2,10,18,19] - for p in params1: - r1.add_sensitivity_reaction(p) + def check_names(reactor, net, params): + for i,(kind,p) in enumerate(params): + rname, comp = net.sensitivity_parameter_name(i).split(': ') + self.assertEqual(reactor.name, rname) + if kind == 'r': + self.assertEqual(gas.reaction_equation(p), comp) + elif kind == 's': + self.assertEqual(p + ' enthalpy', comp) + + params1 = [('r', 2), ('r', 10), ('r', 18), ('r', 19), ('s', 'O2'), + ('s', 'OH'), ('s', 'H2O2')] + r1,net1 = setup(params1) S1 = integrate(r1, net1) + check_names(r1, net1, params1) - pname = lambda r,i: '%s: %s' % (r.name, gas.reaction_equation(i)) - for i,p in enumerate(params1): - self.assertEqual(pname(r1,p), net1.sensitivity_parameter_name(i)) - - r2,net2 = setup() - params2 = [19,10,2,18] - for p in params2: - r2.add_sensitivity_reaction(p) + params2 = [('r', 19), ('s', 'H2O2'), ('s', 'OH'), ('r', 10), + ('s', 'O2'), ('r', 2), ('r', 18)] + r2,net2 = setup(params2) S2 = integrate(r2, net2) + check_names(r2, net2, params2) - for i,p in enumerate(params2): - self.assertEqual(pname(r2,p), net2.sensitivity_parameter_name(i)) - - for i,j in enumerate((2,1,3,0)): + for i,j in enumerate((5,3,6,0,4,2,1)): self.assertArrayNear(S1[:,i], S2[:,j]) def test_parameter_order1a(self): diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index c42b2f0c1..66cf616f8 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -235,7 +235,12 @@ void CVodesIntegrator::sensInit(double t0, FuncEval& func) if (flag != CV_SUCCESS) { throw CanteraError("CVodesIntegrator::sensInit", "Error in CVodeSensInit"); } - vector_fp atol(m_np, m_abstolsens); + vector_fp atol(m_np); + for (size_t n = 0; n < m_np; n++) { + // This scaling factor is tuned so that reaction and species enthalpy + // sensitivities can be computed simultaneously with the same abstol. + atol[n] = m_abstolsens / func.m_paramScales[n]; + } flag = CVodeSensSStolerances(m_cvode_mem, m_reltolsens, atol.data()); } @@ -312,7 +317,11 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func) if (func.nparams() > 0) { sensInit(t0, func); flag = CVodeSetSensParams(m_cvode_mem, func.m_sens_params.data(), - NULL, NULL); + func.m_paramScales.data(), NULL); + if (flag != CV_SUCCESS) { + throw CanteraError("CVodesIntegrator::initialize", + "CVodeSetSensParams failed."); + } } applyOptions(); } diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index a4fd242ed..49db701a6 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -320,8 +320,25 @@ void Reactor::addSensitivityReaction(size_t rxn) "Reaction number out of range ({})", rxn); } - size_t p = network().registerSensitivityReaction(name()+": "+m_kin->reactionString(rxn)); - m_sensParams.emplace_back(SensitivityParameter{rxn, p, 1.0}); + size_t p = network().registerSensitivityParameter( + name()+": "+m_kin->reactionString(rxn), 1.0, 1.0); + m_sensParams.emplace_back( + SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction}); +} + +void Reactor::addSensitivitySpeciesEnthalpy(size_t k) +{ + if (k >= m_thermo->nSpecies()) { + throw CanteraError("Reactor::addSensitivitySpeciesEnthalpy", + "Species index out of range ({})", k); + } + + size_t p = network().registerSensitivityParameter( + name() + ": " + m_thermo->speciesName(k) + " enthalpy", + 0.0, GasConstant * 298.15); + m_sensParams.emplace_back( + SensitivityParameter{k, p, m_thermo->Hf298SS(k), + SensParameterType::enthalpy}); } size_t Reactor::speciesIndex(const string& nm) const @@ -387,12 +404,17 @@ void Reactor::applySensitivity(double* params) return; } for (auto& p : m_sensParams) { - p.value = m_kin->multiplier(p.local); - m_kin->setMultiplier(p.local, p.value*params[p.global]); + if (p.type == SensParameterType::reaction) { + p.value = m_kin->multiplier(p.local); + m_kin->setMultiplier(p.local, p.value*params[p.global]); + } else if (p.type == SensParameterType::enthalpy) { + m_thermo->modifyOneHf298SS(p.local, p.value + params[p.global]); + } } for (size_t m = 0; m < m_wall.size(); m++) { m_wall[m]->setSensitivityParameters(params); } + m_thermo->invalidateCache(); m_kin->invalidateCache(); } @@ -402,11 +424,17 @@ void Reactor::resetSensitivity(double* params) return; } for (auto& p : m_sensParams) { - m_kin->setMultiplier(p.local, p.value); + if (p.type == SensParameterType::reaction) { + m_kin->setMultiplier(p.local, p.value); + } else if (p.type == SensParameterType::enthalpy) { + m_thermo->resetHf298(p.local); + } } for (size_t m = 0; m < m_wall.size(); m++) { m_wall[m]->resetSensitivityParameters(); } + m_thermo->invalidateCache(); + m_kin->invalidateCache(); } } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index bdc14e20f..adc790f66 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -13,7 +13,7 @@ namespace Cantera ReactorNet::ReactorNet() : m_integ(0), m_time(0.0), m_init(false), m_integrator_init(false), m_nv(0), m_rtol(1.0e-9), m_rtolsens(1.0e-4), - m_atols(1.0e-15), m_atolsens(1.0e-4), + m_atols(1.0e-15), m_atolsens(1.0e-6), m_maxstep(0.0), m_maxErrTestFails(0), m_verbose(false) { @@ -179,7 +179,11 @@ double ReactorNet::sensitivity(size_t k, size_t p) throw IndexError("ReactorNet::sensitivity", "m_sens_params", p, m_sens_params.size()-1); } - return m_integ->sensitivity(k, p) / m_integ->solution(k); + double denom = m_integ->solution(k); + if (denom == 0.0) { + denom = SmallNumber; + } + return m_integ->sensitivity(k, p) / denom; } void ReactorNet::evalJacobian(doublereal t, doublereal* y, @@ -238,15 +242,17 @@ size_t ReactorNet::globalComponentIndex(const string& component, size_t reactor) return m_start[reactor] + m_reactors[reactor]->componentIndex(component); } -size_t ReactorNet::registerSensitivityReaction(const std::string& name) +size_t ReactorNet::registerSensitivityParameter( + const std::string& name, double value, double scale) { if (m_integrator_init) { - throw CanteraError("ReactorNet::registerSensitivityReaction", - "Sensitivity reactions cannot be added after the" + throw CanteraError("ReactorNet::registerSensitivityParameter", + "Sensitivity parameters cannot be added after the" "integrator has been initialized."); } m_paramNames.push_back(name); - m_sens_params.push_back(1.0); + m_sens_params.push_back(value); + m_paramScales.push_back(scale); return m_sens_params.size() - 1; } diff --git a/src/zeroD/Wall.cpp b/src/zeroD/Wall.cpp index 234a16b23..16d7ccc9c 100644 --- a/src/zeroD/Wall.cpp +++ b/src/zeroD/Wall.cpp @@ -140,13 +140,15 @@ void Wall::addSensitivityReaction(int leftright, size_t rxn) "Reaction number out of range ({})", rxn); } if (leftright == 0) { - size_t p = m_left->network().registerSensitivityReaction( - m_chem[0]->reactionString(rxn)); - m_pleft.emplace_back(SensitivityParameter{rxn, p, 1.0}); + size_t p = m_left->network().registerSensitivityParameter( + m_chem[0]->reactionString(rxn), 1.0, 1.0); + m_pleft.emplace_back( + SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction}); } else { - size_t p = m_right->network().registerSensitivityReaction( - m_chem[1]->reactionString(rxn)); - m_pright.emplace_back(SensitivityParameter{rxn, p, 1.0}); + size_t p = m_right->network().registerSensitivityParameter( + m_chem[1]->reactionString(rxn), 1.0, 1.0); + m_pright.emplace_back( + SensitivityParameter{rxn, p, 1.0, SensParameterType::reaction}); } }