diff --git a/include/cantera/kinetics/GasKinetics.h b/include/cantera/kinetics/GasKinetics.h index 34e8a31d93..2f28474c83 100644 --- a/include/cantera/kinetics/GasKinetics.h +++ b/include/cantera/kinetics/GasKinetics.h @@ -77,7 +77,9 @@ class GasKinetics : public BulkKinetics virtual Eigen::SparseMatrix fwdRatesOfProgress_ddX(); virtual Eigen::SparseMatrix revRatesOfProgress_ddX(); virtual Eigen::SparseMatrix netRatesOfProgress_ddX(); - virtual Eigen::SparseMatrix netRatesOfProgress_ddC(); + virtual Eigen::SparseMatrix fwdRatesOfProgress_ddN(); + virtual Eigen::SparseMatrix revRatesOfProgress_ddN(); + virtual Eigen::SparseMatrix netRatesOfProgress_ddN(); //! Update temperature-dependent portions of reaction rates and falloff //! functions. @@ -131,17 +133,14 @@ class GasKinetics : public BulkKinetics void process_ddC(StoichManagerN& stoich, const vector_fp& in, double* drop, bool mass_action=true); - //! Process mole fraction derivative + //! Process derivatives //! @param stoich stoichiometry manager //! @param in rate expression used for the derivative calculation - Eigen::SparseMatrix process_ddX(StoichManagerN& stoich, - const vector_fp& in); - - //! Process mole fraction derivative - //! @param stoich stoichiometry manager - //! @param in rate expression used for the derivative calculation - Eigen::SparseMatrix process_ddC(StoichManagerN& stoich, - const vector_fp& in); + //! @param ddX true: w.r.t mole fractions false: w.r.t species concentrations + //! @return a sparse matrix of derivative contributions for each reaction of + //! dimensions nTotalReactions by nTotalSpecies + Eigen::SparseMatrix process_derivatives(StoichManagerN& stoich, + const vector_fp& in, bool ddX=true); //! Helper function ensuring that all rate derivatives can be calculated //! @param name method name used for error output diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 7627825879..a902a6ecd0 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -305,7 +305,16 @@ class InterfaceKinetics : public Kinetics */ double interfaceCurrent(const size_t iphase); - virtual Eigen::SparseMatrix netRatesOfProgress_ddC(); + /** + * Set/modify derivative settings. For the interface kinetics object this + * includes `"skip-cov-dep"`, `"skip-electrochem"`, and `"rtol-delta"`. + * + * @param settings AnyMap containing settings determining derivative evaluation. + */ + virtual void setDerivativeSettings(const AnyMap& settings); + virtual Eigen::SparseMatrix fwdRatesOfProgress_ddN(); + virtual Eigen::SparseMatrix revRatesOfProgress_ddN(); + virtual Eigen::SparseMatrix netRatesOfProgress_ddN(); protected: //! @name Internal service methods @@ -324,9 +333,18 @@ class InterfaceKinetics : public Kinetics //! Process mole fraction derivative //! @param stoich stoichiometry manager //! @param in rate expression used for the derivative calculation - Eigen::SparseMatrix process_ddC(StoichManagerN& stoich, + //! @return a sparse matrix of derivative contributions for each reaction of + //! dimensions nTotalReactions by nTotalSpecies + Eigen::SparseMatrix process_derivatives(StoichManagerN& stoich, const vector_fp& in); + //! Helper function ensuring that all rate derivatives can be calculated + //! @param name method name used for error output + //! @throw CanteraError if coverage dependence or electrochemical reactions are + //! included + void assertDerivativesValid(const std::string& name); + //! @} + //! Temporary work vector of length m_kk vector_fp m_grt; @@ -490,6 +508,11 @@ class InterfaceKinetics : public Kinetics vector_fp m_rbuf0; vector_fp m_rbuf1; vector_fp m_rbuf2; + + // Derivative settings initialized to their default values + bool m_jac_skip_cov_dependance = true; + bool m_jac_skip_electrochem = true; + double m_jac_rtol_delta = 1e-8; }; } diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index b6ed1abb66..0c7471097f 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -589,6 +589,7 @@ class Kinetics //! - `_ddP`: derivative with respect to pressure (a vector) //! - `_ddC`: derivative with respect to molar concentration (a vector) //! - `_ddX`: derivative with respect to species mole fractions (a matrix) + //! - `_ddN`: derivative with respect to species concentrations (a matrix) //! //! Settings for derivative evaluation are set by keyword/value pairs using //! the methods getDerivativeSettings() and setDerivativeSettings(). @@ -716,6 +717,24 @@ class Kinetics "Not implemented for kinetics type '{}'.", kineticsType()); } + /** + * Calculate derivatives for forward rates-of-progress with respect to species + * concentration at constant temperature, pressure and remaining species + * concentrations. + * + * The method returns a matrix with nReactions rows and nTotalSpecies columns. + * For a derivative with respect to \f$n_i\f$, all other \f$n_j\f$ are held + * constant. + * + * @warning This method is an experimental part of the %Cantera API and + * may be changed or removed without notice. + */ + virtual Eigen::SparseMatrix fwdRatesOfProgress_ddN() + { + throw NotImplementedError("Kinetics::fwdRatesOfProgress_ddN", + "Not implemented for kinetics type '{}'.", kineticsType()); + } + /** * Calculate derivatives for reverse rates-of-progress with respect to temperature * at constant pressure, molar concentration and mole fractions. @@ -769,6 +788,24 @@ class Kinetics "Not implemented for kinetics type '{}'.", kineticsType()); } + /** + * Calculate derivatives for forward rates-of-progress with respect to species + * concentration at constant temperature, pressure and remaining species + * concentrations. + * + * The method returns a matrix with nReactions rows and nTotalSpecies columns. + * For a derivative with respect to \f$n_i\f$, all other \f$n_j\f$ are held + * constant. + * + * @warning This method is an experimental part of the %Cantera API and + * may be changed or removed without notice. + */ + virtual Eigen::SparseMatrix revRatesOfProgress_ddN() + { + throw NotImplementedError("Kinetics::revRatesOfProgress_ddN", + "Not implemented for kinetics type '{}'.", kineticsType()); + } + /** * Calculate derivatives for net rates-of-progress with respect to temperature * at constant pressure, molar concentration and mole fractions. @@ -824,7 +861,8 @@ class Kinetics /** * Calculate derivatives for net rates-of-progress with respect to species - * concentration at constant temperature and pressure. + * concentration at constant temperature, pressure, and remaining species + * concentrations. * * The method returns a matrix with nReactions rows and nTotalSpecies columns. * For a derivative with respect to \f$[n_i]\f$, all other \f$[n_i]\f$ are held @@ -833,9 +871,9 @@ class Kinetics * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - virtual Eigen::SparseMatrix netRatesOfProgress_ddC() + virtual Eigen::SparseMatrix netRatesOfProgress_ddN() { - throw NotImplementedError("Kinetics::netRatesOfProgress_ddC", + throw NotImplementedError("Kinetics::netRatesOfProgress_ddN", "Not implemented for kinetics type '{}'.", kineticsType()); } @@ -952,15 +990,17 @@ class Kinetics /** * Calculate derivatives for species net production rates with respect to species - * concentration at constant temperature and pressure. + * concentration at constant temperature, pressure, and concentration of all other + * species. * * The method returns a matrix with nTotalSpecies rows and nTotalSpecies columns. - * For a derivative with respect to \f$[n_i]\f$, all other \f$[n_i]\f$ are held constant. + * For a derivative with respect to \f$[n_i]\f$, all other \f$[n_i]\f$ are held + * constant. * * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - Eigen::SparseMatrix netProductionRates_ddC(); + Eigen::SparseMatrix netProductionRates_ddN(); //! @} //! @name Reaction Mechanism Informational Query Routines diff --git a/include/cantera/zeroD/MoleReactor.h b/include/cantera/zeroD/MoleReactor.h index e2812ea09c..4dbff7c115 100644 --- a/include/cantera/zeroD/MoleReactor.h +++ b/include/cantera/zeroD/MoleReactor.h @@ -36,10 +36,10 @@ class MoleReactor : public Reactor std::string componentName(size_t k); - //! Add the surface chemistry Jacobian values to m_jac_trips - virtual void addSurfJacobian(); - protected: + //! Add the surface chemistry dn/dnj Jacobian values to m_jac_trips + virtual void addSurfJacobian(double cp, bool pressure=false); + //! Get moles of the system from mass fractions stored by thermo object //! @param y vector for moles to be put into virtual void getMoles(double* y); diff --git a/include/cantera/zeroD/ReactorNet.h b/include/cantera/zeroD/ReactorNet.h index 779d8da28d..01eeeb9608 100644 --- a/include/cantera/zeroD/ReactorNet.h +++ b/include/cantera/zeroD/ReactorNet.h @@ -302,9 +302,8 @@ class ReactorNet : public FuncEval virtual Eigen::SparseMatrix jacobian(); protected: - //! Check if surfaces and preconditioning are included, if so throw an error because - //! they are currently not supported. - virtual void checkPreconditionerSupported(){}; + //! check that preconditioning is supported. + virtual void checkPreconditionerSupported() {}; //! Update the preconditioner based on the already computed jacobian values virtual void updatePreconditioner(double gamma); diff --git a/interfaces/cython/cantera/kinetics.pxd b/interfaces/cython/cantera/kinetics.pxd index fd205547dd..bde1783a64 100644 --- a/interfaces/cython/cantera/kinetics.pxd +++ b/interfaces/cython/cantera/kinetics.pxd @@ -60,6 +60,10 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera": CxxSparseMatrix revRatesOfProgress_ddX() except +translate_exception CxxSparseMatrix netRatesOfProgress_ddX() except +translate_exception + CxxSparseMatrix fwdRatesOfProgress_ddN() except +translate_exception + CxxSparseMatrix revRatesOfProgress_ddN() except +translate_exception + CxxSparseMatrix netRatesOfProgress_ddN() except +translate_exception + CxxSparseMatrix creationRates_ddX() except +translate_exception CxxSparseMatrix destructionRates_ddX() except +translate_exception CxxSparseMatrix netProductionRates_ddX() except +translate_exception diff --git a/interfaces/cython/cantera/kinetics.pyx b/interfaces/cython/cantera/kinetics.pyx index 0669af45bf..af1c8e13ca 100644 --- a/interfaces/cython/cantera/kinetics.pyx +++ b/interfaces/cython/cantera/kinetics.pyx @@ -517,6 +517,23 @@ cdef class Kinetics(_SolutionBase): return get_from_sparse(self.kinetics.fwdRatesOfProgress_ddX(), self.n_reactions, self.n_total_species) + property forward_rates_of_progress_ddN: + """ + Calculate derivatives for forward rates-of-progress with respect to species + concentrations at constant temperature, pressure and remaining species + concentrations. + For sparse output, set ``ct.use_sparse(True)``. + + Note that for derivatives with respect to :math:`[n_i]`, all other :math:`[n_j]` + are held constant. + + **Warning:** this property is an experimental part of the Cantera API and + may be changed or removed without notice. + """ + def __get__(self): + return get_from_sparse(self.kinetics.fwdRatesOfProgress_ddN(), + self.n_reactions, self.n_total_species) + property reverse_rates_of_progress_ddT: """ Calculate derivatives for reverse rates-of-progress with respect to temperature @@ -560,6 +577,23 @@ cdef class Kinetics(_SolutionBase): return get_from_sparse(self.kinetics.revRatesOfProgress_ddX(), self.n_reactions, self.n_total_species) + property reverse_rates_of_progress_ddN: + """ + Calculate derivatives for reverse rates-of-progress with respect to species + concentrations at constant temperature, pressure and remaining species + concentrations. + For sparse output, set ``ct.use_sparse(True)``. + + Note that for derivatives with respect to :math:`[n_i]`, all other :math:`[n_j]` + are held constant. + + **Warning:** this property is an experimental part of the Cantera API and + may be changed or removed without notice. + """ + def __get__(self): + return get_from_sparse(self.kinetics.revRatesOfProgress_ddN(), + self.n_reactions, self.n_total_species) + property net_rates_of_progress_ddT: """ Calculate derivatives for net rates-of-progress with respect to temperature @@ -603,6 +637,23 @@ cdef class Kinetics(_SolutionBase): return get_from_sparse(self.kinetics.netRatesOfProgress_ddX(), self.n_reactions, self.n_total_species) + property net_rates_of_progress_ddN: + """ + Calculate derivatives for net rates-of-progress with respect to species + concentrations at constant temperature, pressure and remaining species + concentrations. + For sparse output, set ``ct.use_sparse(True)``. + + Note that for derivatives with respect to :math:`[n_i]`, all other :math:`[n_j]` + are held constant. + + **Warning:** this property is an experimental part of the Cantera API and + may be changed or removed without notice. + """ + def __get__(self): + return get_from_sparse(self.kinetics.netRatesOfProgress_ddN(), + self.n_reactions, self.n_total_species) + property creation_rates_ddT: """ Calculate derivatives of species creation rates with respect to temperature diff --git a/interfaces/cython/cantera/preconditioners.pxd b/interfaces/cython/cantera/preconditioners.pxd index 1f2a59c3f6..ebe80f6a02 100644 --- a/interfaces/cython/cantera/preconditioners.pxd +++ b/interfaces/cython/cantera/preconditioners.pxd @@ -16,7 +16,7 @@ cdef extern from "cantera/numerics/PreconditionerBase.h" namespace "Cantera": cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera": cdef cppclass CxxAdaptivePreconditioner "Cantera::AdaptivePreconditioner" \ (CxxPreconditionerBase): - CxxAdaptivePreconditioner() except + + CxxAdaptivePreconditioner() except +translate_exception void setThreshold(double threshold) double threshold() void setIlutFillFactor(int fillfactor) @@ -24,7 +24,7 @@ cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera" void setIlutDropTol(double droptol) double ilutDropTol() void printPreconditioner() - CxxSparseMatrix matrix() except + + CxxSparseMatrix matrix() except +translate_exception cdef extern from "cantera/numerics/PreconditionerFactory.h" namespace "Cantera": cdef shared_ptr[CxxPreconditionerBase] newPreconditioner(string) except\ diff --git a/interfaces/cython/cantera/reactor.pyx b/interfaces/cython/cantera/reactor.pyx index 05c434a144..cc46e1dc62 100644 --- a/interfaces/cython/cantera/reactor.pyx +++ b/interfaces/cython/cantera/reactor.pyx @@ -1605,7 +1605,7 @@ cdef class ReactorNet: property jacobian: """ - Get the Jacobian or an approximation thereof + Get the system Jacobian or an approximation thereof. **Warning**: Depending on the particular implementation, this may return an approximate Jacobian intended only for use in forming a preconditioner for diff --git a/src/kinetics/GasKinetics.cpp b/src/kinetics/GasKinetics.cpp index df5e08b2eb..6d2506d72b 100644 --- a/src/kinetics/GasKinetics.cpp +++ b/src/kinetics/GasKinetics.cpp @@ -403,17 +403,20 @@ void GasKinetics::getNetRatesOfProgress_ddC(double* drop) dNetRop -= dNetRop2; } -Eigen::SparseMatrix GasKinetics::process_ddX( - StoichManagerN& stoich, const vector_fp& in) +Eigen::SparseMatrix GasKinetics::process_derivatives( + StoichManagerN& stoich, const vector_fp& in, bool ddX) { Eigen::SparseMatrix out; vector_fp& scaled = m_rbuf1; vector_fp& outV = m_rbuf2; // convert from concentration to mole fraction output - double ctot = thermo().molarDensity(); - for (size_t i = 0; i < nReactions(); ++i) { - scaled[i] = ctot * in[i]; + copy(in.begin(), in.end(), scaled.begin()); + if (ddX) { + double ctot = thermo().molarDensity(); + for (size_t i = 0; i < nReactions(); ++i) { + scaled[i] *= ctot; + } } // derivatives handled by StoichManagerN @@ -443,39 +446,6 @@ Eigen::SparseMatrix GasKinetics::process_ddX( return out; } -Eigen::SparseMatrix GasKinetics::process_ddC( - StoichManagerN& stoich, const vector_fp& in) -{ - Eigen::SparseMatrix out; - vector_fp& outV = m_rbuf1; - - // derivatives handled by StoichManagerN - copy(in.begin(), in.end(), outV.begin()); - processThirdBodies(outV.data()); - out = stoich.derivatives(m_act_conc.data(), outV.data()); - if (m_jac_skip_third_bodies || m_multi_concm.empty()) { - return out; - } - - // derivatives due to law of mass action - copy(in.begin(), in.end(), outV.begin()); - stoich.multiply(m_act_conc.data(), outV.data()); - - // derivatives due to reaction rates depending on third-body colliders - if (!m_jac_skip_falloff) { - for (auto& rates : m_bulk_rates) { - // processing step does not modify entries not dependent on M - rates->processRateConstants_ddM( - outV.data(), m_rfn.data(), m_jac_rtol_delta, false); - } - } - - // derivatives handled by ThirdBodyCalc - out += m_multi_concm.derivatives(outV.data()); - - return out; -} - Eigen::SparseMatrix GasKinetics::fwdRatesOfProgress_ddX() { assertDerivativesValid("GasKinetics::fwdRatesOfProgress_ddX"); @@ -483,7 +453,7 @@ Eigen::SparseMatrix GasKinetics::fwdRatesOfProgress_ddX() // forward reaction rate coefficients vector_fp& rop_rates = m_rbuf0; processFwdRateCoefficients(rop_rates.data()); - return process_ddX(m_reactantStoich, rop_rates); + return process_derivatives(m_reactantStoich, rop_rates); } Eigen::SparseMatrix GasKinetics::revRatesOfProgress_ddX() @@ -494,7 +464,7 @@ Eigen::SparseMatrix GasKinetics::revRatesOfProgress_ddX() vector_fp& rop_rates = m_rbuf0; processFwdRateCoefficients(rop_rates.data()); processEquilibriumConstants(rop_rates.data()); - return process_ddX(m_revProductStoich, rop_rates); + return process_derivatives(m_revProductStoich, rop_rates); } Eigen::SparseMatrix GasKinetics::netRatesOfProgress_ddX() @@ -504,26 +474,46 @@ Eigen::SparseMatrix GasKinetics::netRatesOfProgress_ddX() // forward reaction rate coefficients vector_fp& rop_rates = m_rbuf0; processFwdRateCoefficients(rop_rates.data()); - Eigen::SparseMatrix jac = process_ddX(m_reactantStoich, rop_rates); + Eigen::SparseMatrix jac = process_derivatives(m_reactantStoich, rop_rates); // reverse reaction rate coefficients processEquilibriumConstants(rop_rates.data()); - return jac - process_ddX(m_revProductStoich, rop_rates); + return jac - process_derivatives(m_revProductStoich, rop_rates); } -Eigen::SparseMatrix GasKinetics::netRatesOfProgress_ddC() +Eigen::SparseMatrix GasKinetics::fwdRatesOfProgress_ddN() { - assertDerivativesValid("GasKinetics::netRatesOfProgress_ddC"); + assertDerivativesValid("GasKinetics::fwdRatesOfProgress_ddN"); // forward reaction rate coefficients vector_fp& rop_rates = m_rbuf0; processFwdRateCoefficients(rop_rates.data()); - Eigen::SparseMatrix jac = process_ddC(m_reactantStoich, rop_rates); + return process_derivatives(m_reactantStoich, rop_rates, false); +} + +Eigen::SparseMatrix GasKinetics::revRatesOfProgress_ddN() +{ + assertDerivativesValid("GasKinetics::revRatesOfProgress_ddN"); + + // reverse reaction rate coefficients + vector_fp& rop_rates = m_rbuf0; + processFwdRateCoefficients(rop_rates.data()); + processEquilibriumConstants(rop_rates.data()); + return process_derivatives(m_revProductStoich, rop_rates, false); +} +Eigen::SparseMatrix GasKinetics::netRatesOfProgress_ddN() +{ + assertDerivativesValid("GasKinetics::netRatesOfProgress_ddN"); + + // forward reaction rate coefficients + vector_fp& rop_rates = m_rbuf0; + processFwdRateCoefficients(rop_rates.data()); + Eigen::SparseMatrix jac = process_derivatives(m_reactantStoich, rop_rates, false); // reverse reaction rate coefficients processEquilibriumConstants(rop_rates.data()); - return jac - process_ddC(m_revProductStoich, rop_rates); + return jac - process_derivatives(m_revProductStoich, rop_rates, false); } bool GasKinetics::addReaction(shared_ptr r, bool resize) diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index db495b4e31..b129b5e8e9 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -615,17 +615,55 @@ double InterfaceKinetics::interfaceCurrent(const size_t iphase) return dotProduct * Faraday; } -Eigen::SparseMatrix InterfaceKinetics::netRatesOfProgress_ddC() +Eigen::SparseMatrix InterfaceKinetics::fwdRatesOfProgress_ddN() { + // check derivatives are valid + assertDerivativesValid("InterfaceKinetics::fwdRatesOfProgress_ddN"); + // forward reaction rate coefficients + vector_fp& rop_rates = m_rbuf0; + processFwdRateCoefficients(rop_rates.data()); + return process_derivatives(m_reactantStoich, rop_rates); +} + +Eigen::SparseMatrix InterfaceKinetics::revRatesOfProgress_ddN() +{ + // check derivatives are valid + assertDerivativesValid("InterfaceKinetics::revRatesOfProgress_ddN"); + // reverse reaction rate coefficients + vector_fp& rop_rates = m_rbuf0; + processFwdRateCoefficients(rop_rates.data()); + processEquilibriumConstants(rop_rates.data()); + return process_derivatives(m_revProductStoich, rop_rates); +} +Eigen::SparseMatrix InterfaceKinetics::netRatesOfProgress_ddN() +{ + // check derivatives are valid + assertDerivativesValid("InterfaceKinetics::netRatesOfProgress_ddN"); // forward reaction rate coefficients vector_fp& rop_rates = m_rbuf0; processFwdRateCoefficients(rop_rates.data()); - Eigen::SparseMatrix jac = process_ddC(m_reactantStoich, rop_rates); + Eigen::SparseMatrix jac = process_derivatives(m_reactantStoich, rop_rates); // reverse reaction rate coefficients processEquilibriumConstants(rop_rates.data()); - return jac - process_ddC(m_revProductStoich, rop_rates); + return jac - process_derivatives(m_revProductStoich, rop_rates); +} + +void InterfaceKinetics::setDerivativeSettings(const AnyMap& settings) +{ + bool force = settings.empty(); + if (force || settings.hasKey("skip-cov-dep")) { + m_jac_skip_cov_dependance = settings.getBool("skip-cov-dep", + true); + } + if (force || settings.hasKey("skip-electrochem")) { + m_jac_skip_electrochem = settings.getBool("skip-electrochem", + true); + } + if (force || settings.hasKey("rtol-delta")) { + m_jac_rtol_delta = settings.getDouble("rtol-delta", 1e-8); + } } void InterfaceKinetics::processFwdRateCoefficients(double* ropf) @@ -644,7 +682,7 @@ void InterfaceKinetics::processFwdRateCoefficients(double* ropf) } } -Eigen::SparseMatrix InterfaceKinetics::process_ddC( +Eigen::SparseMatrix InterfaceKinetics::process_derivatives( StoichManagerN& stoich, const vector_fp& in) { Eigen::SparseMatrix out; @@ -655,6 +693,16 @@ Eigen::SparseMatrix InterfaceKinetics::process_ddC( return out; } +void InterfaceKinetics::assertDerivativesValid(const std::string& name) +{ + if (!m_jac_skip_cov_dependance) { + throw NotImplementedError(name, "Coverage-dependent reactions not supported."); + } + else if (!m_jac_skip_electrochem) { + throw NotImplementedError(name, "Electrochemical reactions not supported."); + } +} + void InterfaceKinetics::processEquilibriumConstants(double* rop) { // For reverse rates computed from thermochemistry, multiply the forward diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index 942586421b..0b0ec27797 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -566,9 +566,9 @@ Eigen::SparseMatrix Kinetics::netProductionRates_ddX() return m_stoichMatrix * netRatesOfProgress_ddX(); } -Eigen::SparseMatrix Kinetics::netProductionRates_ddC() +Eigen::SparseMatrix Kinetics::netProductionRates_ddN() { - return m_stoichMatrix * netRatesOfProgress_ddC(); + return m_stoichMatrix * netRatesOfProgress_ddN(); } void Kinetics::addThermo(shared_ptr thermo) diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index 612514ae4f..caf0bb0897 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -131,20 +131,21 @@ Eigen::SparseMatrix IdealGasConstPressureMoleReactor::jacobian() // volume / moles * rates portion of equation Eigen::VectorXd netProductionRates(m_nsp); m_kin->getNetProductionRates(netProductionRates.data()); // "omega dot" - Eigen::SparseMatrix dwdC = m_kin->netProductionRates_ddC(); + Eigen::SparseMatrix dnk_dnj = m_kin->netProductionRates_ddN(); double molarVolume = m_thermo->molarVolume(); // Calculate ROP derivatives, excluding the term // molarVolume * (wdot(j) - sum_k(X_k * dwdot_j/dX_k)) // which is small and would completely destroy the sparsity of the Jacobian - for (int k = 0; k < dwdC.outerSize(); k++) { - for (Eigen::SparseMatrix::InnerIterator it(dwdC, k); it; ++it) { + for (int k = 0; k < dnk_dnj.outerSize(); k++) { + for (Eigen::SparseMatrix::InnerIterator it(dnk_dnj, k); it; ++it) { + it.valueRef() = it.value() + netProductionRates[it.row()] * molarVolume; m_jac_trips.emplace_back(static_cast(it.row() + m_sidx), - static_cast(it.col() + m_sidx), - it.value() + netProductionRates[it.row()] * molarVolume); + static_cast(it.col() + m_sidx), it.value()); } } // Temperature Derivatives + double NCp = 0.0; if (m_energy) { // getting perturbed state for finite difference double deltaTemp = m_thermo->temperature() @@ -172,26 +173,33 @@ Eigen::SparseMatrix IdealGasConstPressureMoleReactor::jacobian() (ydotPerturbed - ydotCurrent) / deltaTemp); } // d T_dot/dnj - Eigen::VectorXd specificHeat(m_nsp); Eigen::VectorXd enthalpy(m_nsp); - Eigen::VectorXd dwdot_dC(m_nsp); - m_thermo->getPartialMolarCp(specificHeat.data()); m_thermo->getPartialMolarEnthalpies(enthalpy.data()); - double qdot = enthalpy.dot(netProductionRates); - double cp_mole = m_thermo->cp_mole(); - double total_moles = m_vol / molarVolume; - double denom = (total_moles * total_moles * cp_mole * cp_mole); - // determine derivatives - // spans columns - Eigen::VectorXd hk_dnkdnj_sum = enthalpy.transpose() * dwdC; - for (int j = 0; j < m_nsp; j++) { + // get specific heat for whole system + Eigen::VectorXd specificHeat(m_nv - m_sidx); + // gas phase + m_thermo->getPartialMolarCp(specificHeat.data()); + // surface phases + size_t shift = m_nsp; + for (auto S : m_surfaces) { + S->thermo()->getPartialMolarCp(specificHeat.data() + shift); + shift += S->thermo()->nSpecies(); + } + double qdot = m_vol * enthalpy.dot(netProductionRates); + // find denominator ahead of time + for (size_t i = 0; i < m_nv - m_sidx; i++) { + NCp += yCurrent[i + m_sidx] * specificHeat[i]; + } + double denom = 1 / (NCp * NCp); + Eigen::VectorXd hk_dnkdnj_sum = dnk_dnj.transpose() * enthalpy; + // Add derivatives to jac by spanning columns + for (size_t j = 0; j < m_nsp; j++) { m_jac_trips.emplace_back(0, static_cast(j + m_sidx), - ((specificHeat[j] - cp_mole) * m_vol * qdot - - total_moles * cp_mole * hk_dnkdnj_sum[j]) / denom); + (specificHeat[j] * qdot - NCp * hk_dnkdnj_sum[j]) * denom); } } // add surface jacobian to system - addSurfJacobian(); + addSurfJacobian(NCp, true); // convert triplets to sparse matrix Eigen::SparseMatrix jac(m_nv, m_nv); jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index 7e96ff96f7..2d41c943e8 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -164,11 +164,10 @@ Eigen::SparseMatrix IdealGasMoleReactor::jacobian() // Determine Species Derivatives // get ROP derivatives, excluding the term molarVolume * sum_k(X_k * dwdot_j/dX_j) // which is small and would completely destroy the sparsity of the Jacobian - Eigen::SparseMatrix dwdC = - m_kin->netProductionRates_ddC() / m_vol; + Eigen::SparseMatrix dnk_dnj = m_kin->netProductionRates_ddN(); // add to preconditioner - for (int k=0; k::InnerIterator it(dwdC, k); it; ++it) { + for (int k=0; k::InnerIterator it(dnk_dnj, k); it; ++it) { m_jac_trips.emplace_back( static_cast(it.row() + m_sidx), static_cast(it.col() + m_sidx), @@ -176,6 +175,7 @@ Eigen::SparseMatrix IdealGasMoleReactor::jacobian() } } // Temperature Derivatives + double NCv = 0.0; if (m_energy) { // getting perturbed state for finite difference double deltaTemp = m_thermo->temperature() @@ -202,36 +202,43 @@ Eigen::SparseMatrix IdealGasMoleReactor::jacobian() m_jac_trips.emplace_back(static_cast(j), 0, (ydotPerturbed - ydotCurrent) / deltaTemp); } - // find derivatives d T_dot/dNj - Eigen::VectorXd specificHeat(m_nsp); + // d T_dot/dnj Eigen::VectorXd netProductionRates(m_nsp); Eigen::VectorXd internal_energy(m_nsp); // getting species data - m_thermo->getPartialMolarCp(specificHeat.data()); m_thermo->getPartialMolarIntEnergies(internal_energy.data()); m_kin->getNetProductionRates(netProductionRates.data()); - // scale net production rates by volume to get molar rate - netProductionRates *= m_vol; + // get specific heat for whole system + Eigen::VectorXd specificHeat(m_nv - m_sidx); + // gas phase + m_thermo->getPartialMolarCp(specificHeat.data()); // convert Cp to Cv for ideal gas as Cp - Cv = R - specificHeat = specificHeat.array() - GasConstant; - // dwdot_dC - Eigen::VectorXd dwdot_dC(m_nsp); - m_kin->getNetProductionRates_ddC(dwdot_dC.data()); - // finding a sums inside the derivative - double totalCv = m_mass * m_thermo->cv_mass(); - double uk_dwdC_sum = internal_energy.transpose() * dwdot_dC; - double uk_nk_sum = internal_energy.transpose() * netProductionRates; - Eigen::VectorXd uk_dnkdnj_sums = dwdC.transpose() * internal_energy; - // finding derivatives - // spans columns + for (size_t i = 0; i < m_nsp; i++) { + specificHeat[i] -= GasConstant; + } + // surface phases, Cp = Cv for surfaces + size_t shift = m_nsp; + for (auto S : m_surfaces) { + S->thermo()->getPartialMolarCp(specificHeat.data() + shift); + shift += S->thermo()->nSpecies(); + } + // scale net production rates by volume to get molar rate + double qdot = m_vol * internal_energy.dot(netProductionRates); + // find denominator ahead of time + for (size_t i = 0; i < m_nv - m_sidx; i++) { + NCv += yCurrent[i + m_sidx] * specificHeat[i]; + } + // Make denominator beforehand + double denom = 1 / (NCv * NCv); + Eigen::VectorXd uk_dnkdnj_sums = dnk_dnj.transpose() * internal_energy; + // Add derivatives to jac by spanning columns for (size_t j = 0; j < m_nsp; j++) { - // set appropriate column of preconditioner m_jac_trips.emplace_back(0, static_cast(j + m_sidx), - (uk_dwdC_sum - uk_dnkdnj_sums[j] + specificHeat[j] * uk_nk_sum / totalCv) / totalCv); + (specificHeat[j] * qdot - NCv * uk_dnkdnj_sums[j]) * denom); } } // add surface jacobian to system - addSurfJacobian(); + addSurfJacobian(NCv, false); // convert triplets to sparse matrix Eigen::SparseMatrix jac(m_nv, m_nv); jac.setFromTriplets(m_jac_trips.begin(), m_jac_trips.end()); diff --git a/src/zeroD/MoleReactor.cpp b/src/zeroD/MoleReactor.cpp index 7c4c909dcb..d5ed18a8cf 100644 --- a/src/zeroD/MoleReactor.cpp +++ b/src/zeroD/MoleReactor.cpp @@ -84,35 +84,70 @@ void MoleReactor::evalSurfaces(double* LHS, double* RHS, double* sdot) } } -void MoleReactor::addSurfJacobian() +void MoleReactor::addSurfJacobian(double NCp, bool pressure) { - // Translate surface netProductionRates_ddC into mole values for reactor - // Jacobian + // For all surfaces find necessary species derivatives + std::string gas_phase = m_thermo->name(); for (auto& S : m_surfaces) { auto curr_kin = S->kinetics(); - // get surface jacobian - Eigen::SparseMatrix surfJac = curr_kin->netProductionRates_ddC(); + auto curr_thermo = S->thermo(); + double V_A = m_vol / S->area(); + // get surface jacobian in concentration units + Eigen::SparseMatrix surfJac = curr_kin->netProductionRates_ddN(); // Add elements to jacobian triplets for (int k=0; k::InnerIterator it(surfJac, k); it; ++it) { // Get species row and column inside of reactor size_t row = componentIndex(curr_kin->kineticsSpeciesName(it.row())); size_t col = componentIndex(curr_kin->kineticsSpeciesName(it.col())); - double scalar = 1; if (row != npos && col != npos) { - if (row > m_nsp && col > m_nsp) { - scalar = 1 / S->thermo()->size(it.row()); + // derivatives w.r.t gas phase species + if (col < m_nsp) { + it.valueRef() = V_A * it.value(); } - // surf w.r.t gas phase scalar - else if (row > m_nsp) { - scalar = S->area() / (S->thermo()->size(it.row()) * m_vol); - } - // gas phase w.r.t surf scalar - else if (col > m_nsp) { - scalar = m_vol / S->area(); - } - // add derivative to the jacobian - m_jac_trips.emplace_back(row, col, scalar * it.value()); + // add to appropriate row and col + m_jac_trips.emplace_back(row, col, it.value()); + } + } + } + // Add thermo jacobian for surfaces + if (m_energy) { + // d T_dot/dnj + // constants used for indexing + size_t nspec = curr_kin->nTotalSpecies(); + size_t surf_spec = curr_thermo->nSpecies(); + size_t ns = curr_kin->phaseIndex(gas_phase); + size_t gasloc = curr_kin->kineticsSpeciesIndex(0, ns); + size_t gasbound = gasloc + m_nsp; + // create needed vectors + Eigen::VectorXd specificHeat(nspec); + Eigen::VectorXd enthalpy(nspec); + Eigen::VectorXd netProductionRates(nspec); + // get needed surface data + curr_thermo->getPartialMolarCp(specificHeat.data()); + curr_thermo->getPartialMolarEnthalpies(enthalpy.data()); + curr_kin->getNetProductionRates(netProductionRates.data()); + // get needed gas phase data + m_thermo->getPartialMolarCp(specificHeat.data() + gasloc); + if (pressure) { + m_thermo->getPartialMolarEnthalpies(enthalpy.data() + gasloc); + } else { + m_thermo->getPartialMolarIntEnergies(enthalpy.data() + gasloc); + for (size_t i = gasloc; i < gasbound; i++) { + specificHeat[i] -= GasConstant; + } + } + // calculate need parameters + double qdot = S->area() * enthalpy.dot(netProductionRates); + // Make denominators beforehand + double denom = 1 / (NCp * NCp); + Eigen::VectorXd hk_dnkdnj_sum = surfJac.transpose() * enthalpy; + // Add derivatives to jac by spanning columns + for (size_t j = 0; j < nspec; j++) { + if (j < surf_spec || (j >= gasloc && j < gasbound)) { + size_t col = componentIndex(curr_kin->kineticsSpeciesName(j)); + m_jac_trips.emplace_back(0, static_cast(col), + (specificHeat[j] * qdot - NCp * hk_dnkdnj_sum[j]) * denom); } } } diff --git a/test/python/test_jacobian.py b/test/python/test_jacobian.py index 8d17d546d4..122fd509bb 100644 --- a/test/python/test_jacobian.py +++ b/test/python/test_jacobian.py @@ -48,7 +48,7 @@ def test_input(self): self.assertEqual(self.equation, self.rxn.equation) self.assertEqual(self.rate_type, self.rxn.rate.type) - def rop_ddX(self, spc_ix, mode, const_t=True, rtol_deltac=1e-5, atol_deltac=1e-20): + def rop_derivs(self, spc_ix, mode, const_t=True, rtol_deltac=1e-5, atol_deltac=1e-20, ddX=True): # numerical derivative for rates-of-progress with respect to mole fractions def calc(): if mode == "forward": @@ -77,7 +77,10 @@ def calc(): self.gas.TPX = tnew, self.gas.P, conc / ctot1 drop = (calc() - rop0) / dconc self.gas.TPX = self.tpx - return drop * self.gas.density_mole + if ddX: + return drop * self.gas.density_mole + else: + return drop def test_forward_rop_ddX(self): # check derivatives of forward rates of progress with respect to mole fractions @@ -96,7 +99,7 @@ def test_forward_rop_ddX(self): self.assertNear(rop[self.rxn_idx], drop[self.rxn_idx, spc_ix] * self.gas.X[spc_ix] / order) - drop_num = self.rop_ddX(spc_ix, mode="forward") + drop_num = self.rop_derivs(spc_ix, mode="forward") self.assertArrayNear(dropm[:, spc_ix] + dropp * self.gas.P, drop_num, self.rtol) if isinstance(self.rxn.rate, ct.FalloffRate): @@ -122,7 +125,7 @@ def test_reverse_rop_ddX(self): self.assertNear(rop[self.rxn_idx], drop[self.rxn_idx, spc_ix] * self.gas.X[spc_ix] / order) - drop_num = self.rop_ddX(spc_ix, mode="reverse") + drop_num = self.rop_derivs(spc_ix, mode="reverse") self.assertArrayNear(dropm[:, spc_ix] + dropp * self.gas.P, drop_num, self.rtol) if not self.rxn.reversible or isinstance(self.rxn.rate, ct.FalloffRate): @@ -141,7 +144,7 @@ def test_net_rop_ddX(self): dropp = self.gas.net_rates_of_progress_ddP for spc_ix in self.rix + self.pix: - drop_num = self.rop_ddX(spc_ix, mode="net") + drop_num = self.rop_derivs(spc_ix, mode="net") ix = drop[:, spc_ix] != 0 drop_ = drop[:, spc_ix] + dropp * self.gas.P self.assertArrayNear(drop_[ix], drop_num[ix], self.rtol) @@ -155,6 +158,76 @@ def test_net_rop_ddX(self): drop[self.rxn_idx, spc_ix] = 0 self.assertFalse(drop.any()) + def test_forward_rop_ddN(self): + # check derivatives of forward rates of progress with respect to species + # concentrations against analytic result + dropm = self.gas.forward_rates_of_progress_ddN + dropp = self.gas.forward_rates_of_progress_ddP + + self.gas.derivative_settings = {"skip-third-bodies": True} + drop = self.gas.forward_rates_of_progress_ddN + rop = self.gas.forward_rates_of_progress + for spc_ix in self.rix: + if self.orders is None: + order = self.r_stoich[spc_ix, self.rxn_idx] + else: + order = self.orders[self.gas.species_names[spc_ix]] + self.assertNear(rop[self.rxn_idx], + drop[self.rxn_idx, spc_ix] * self.gas.concentrations[spc_ix] / order) + + drop_num = self.rop_derivs(spc_ix, mode="forward", ddX=False) + self.assertArrayNear(dropm[:, spc_ix] + dropp * self.gas.P, drop_num, 1e-3) + + if isinstance(self.rxn.rate, ct.FalloffRate): + return + + # ensure all zeros are in the correct spots + for spc_ix in set(self.rix + self.ix3b): + self.assertTrue(dropm[self.rxn_idx, spc_ix]) # non-zero + dropm[self.rxn_idx, spc_ix] = 0 + self.assertFalse(dropm.any()) + + def test_reverse_rop_ddN(self): + # check derivatives of reverse rates of progress with respect to species + # concentrations against analytic result + dropm = self.gas.reverse_rates_of_progress_ddN + dropp = self.gas.reverse_rates_of_progress_ddP + + self.gas.derivative_settings = {"skip-third-bodies": True} + drop = self.gas.reverse_rates_of_progress_ddN + rop = self.gas.reverse_rates_of_progress + for spc_ix in self.pix: + order = self.p_stoich[spc_ix, self.rxn_idx] + self.assertNear(rop[self.rxn_idx], + drop[self.rxn_idx, spc_ix] * self.gas.concentrations[spc_ix] / order) + + drop_num = self.rop_derivs(spc_ix, mode="reverse", ddX=False) + self.assertArrayNear(dropm[:, spc_ix] + dropp * self.gas.P, drop_num, 1e-3) + + if not self.rxn.reversible or isinstance(self.rxn.rate, ct.FalloffRate): + return + + # ensure all zeros are in the correct spots + for spc_ix in set(self.pix + self.ix3b): + self.assertTrue(dropm[self.rxn_idx, spc_ix]) # non-zero + dropm[self.rxn_idx, spc_ix] = 0 + self.assertFalse(dropm.any()) + + def test_net_rop_ddN(self): + # check derivatives of net rates of progress with respect to species + # concentrations against numeric result + drop = self.gas.net_rates_of_progress_ddN + dropp = self.gas.net_rates_of_progress_ddP + + for spc_ix in self.rix + self.pix: + drop_num = self.rop_derivs(spc_ix, mode="net", ddX=False) + ix = drop[:, spc_ix] != 0 + drop_ = drop[:, spc_ix] + dropp * self.gas.P + self.assertArrayNear(drop_[ix], drop_num[ix], 1e-4) + + if not self.rxn.reversible or isinstance(self.rxn.rate, ct.FalloffRate): + return + def rop_ddT(self, mode=None, const_p=False, rtol=1e-6): # numerical derivative for rates-of-progress at constant pressure def calc(): @@ -539,7 +612,7 @@ def setUp(self): self.gas.TPX = self.tpx self.gas.derivative_settings = {} # reset - def rop_ddX(self, mode, rtol_deltac=1e-9, atol_deltac=1e-20): + def rop_derivs(self, mode, rtol_deltac=1e-9, atol_deltac=1e-20, ddX=True): # numerical derivative for rates-of-progress with respect to mole fractions def calc(): if mode == "forward": @@ -566,13 +639,16 @@ def calc(): drop[:, spc_ix] = (calc() - rop0) / dconc self.gas.TPX = self.tpx - return drop * self.gas.density_mole + if ddX: + return drop * self.gas.density_mole + else: + return drop def test_forward_rop_ddX(self): # check forward rop against numerical derivative with respect to mole fractions drop = self.gas.forward_rates_of_progress_ddX dropp = self.gas.forward_rates_of_progress_ddP - drop_num = self.rop_ddX(mode="forward") + drop_num = self.rop_derivs(mode="forward") stoich = self.gas.reactant_stoich_coeffs for i in range(self.gas.n_reactions): try: @@ -590,7 +666,7 @@ def test_reverse_rop_ddX(self): # check reverse rop against numerical derivative with respect to mole fractions drop = self.gas.reverse_rates_of_progress_ddX dropp = self.gas.reverse_rates_of_progress_ddP - drop_num = self.rop_ddX(mode="reverse") + drop_num = self.rop_derivs(mode="reverse") stoich = self.gas.product_stoich_coeffs for i in range(self.gas.n_reactions): try: @@ -608,7 +684,65 @@ def test_net_rop_ddX(self): # check net rop against numerical derivative with respect to mole fractions drop = self.gas.net_rates_of_progress_ddX dropp = self.gas.net_rates_of_progress_ddP - drop_num = self.rop_ddX(mode="net") + drop_num = self.rop_derivs(mode="net") + stoich = self.gas.product_stoich_coeffs - self.gas.reactant_stoich_coeffs + for i in range(self.gas.n_reactions): + try: + # test entries that are not spurious + ix = np.abs((stoich[:, i] != 0) * drop[i, :]) > 1e-6 + drop_ = drop[i, ix] + dropp[i] * self.gas.P + self.assertArrayNear(drop_, drop_num[i, ix], self.rtol) + except AssertionError as err: + if self.gas.reaction(i).reversible: + print(i, self.gas.reaction(i).rate.type) + print(self.gas.reaction(i)) + print(np.vstack([drop[i, ix], drop_num[i, ix]]).T) + raise err + + def test_forward_rop_ddN(self): + # check forward rop against numerical derivative with respect to species + # concentrations + drop = self.gas.forward_rates_of_progress_ddN + dropp = self.gas.forward_rates_of_progress_ddP + drop_num = self.rop_derivs(mode="forward", ddX=False) + stoich = self.gas.reactant_stoich_coeffs + for i in range(self.gas.n_reactions): + try: + # test entries that are not spurious + ix = np.abs((stoich[:, i] != 0) * drop[i, :]) > 1e-6 + drop_ = drop[i, ix] + dropp[i] * self.gas.P + self.assertArrayNear(drop_, drop_num[i, ix], self.rtol) + except AssertionError as err: + print(i, self.gas.reaction(i).rate.type) + print(self.gas.reaction(i)) + print(np.vstack([drop[i, ix], drop_num[i, ix]]).T) + raise err + + def test_reverse_rop_ddN(self): + # check reverse rop against numerical derivative with respect to species + # concentrations + drop = self.gas.reverse_rates_of_progress_ddN + dropp = self.gas.reverse_rates_of_progress_ddP + drop_num = self.rop_derivs(mode="reverse", ddX=False) + stoich = self.gas.product_stoich_coeffs + for i in range(self.gas.n_reactions): + try: + # test entries that are not spurious + ix = np.abs((stoich[:, i] != 0) * drop[i, :]) > 1e-6 + drop_ = drop[i, ix] + dropp[i] * self.gas.P + self.assertArrayNear(drop_, drop_num[i, ix], self.rtol) + except AssertionError as err: + print(i, self.gas.reaction(i).rate.type) + print(self.gas.reaction(i)) + print(np.vstack([drop[i, ix], drop_num[i, ix]]).T) + raise err + + def test_net_rop_ddN(self): + # check net rop against numerical derivative with respect to species + # concentrations + drop = self.gas.net_rates_of_progress_ddN + dropp = self.gas.net_rates_of_progress_ddP + drop_num = self.rop_derivs(mode="net", ddX=False) stoich = self.gas.product_stoich_coeffs - self.gas.reactant_stoich_coeffs for i in range(self.gas.n_reactions): try: @@ -702,3 +836,275 @@ def setUpClass(cls): cls.gas.TPX = 300, 2 * ct.one_atm, "H2:1, O2:3, AR:0.4" cls.gas.equilibrate("HP") super().setUpClass() + + +class SurfaceRateExpressionTests: + # Generic test class to check derivatives evaluated for a single reaction within + # a reaction mechanism for surfaces + + rxn_idx = None # index of reaction to be tested + phase = None + rtol = 1e-5 + orders = None + ix3b = [] # three-body indices + equation = None + rate_type = None + + @classmethod + def setUpClass(cls): + # all species indexs + all_species = cls.surf.species() + cls.gas.species() + cls.sidxs = {spec.name:i for i, spec in enumerate(all_species)} + + # kinetics objects + cls.r_stoich = cls.surf.reactant_stoich_coeffs + cls.p_stoich = cls.surf.product_stoich_coeffs + cls.rxn = cls.surf.reactions()[cls.rxn_idx] + cls.rix = [cls.sidxs[k] for k in cls.rxn.reactants.keys()] + cls.pix = [cls.sidxs[k] for k in cls.rxn.products.keys()] + + def setUp(self): + # gas phase + self.gas.TPX = self.gas_tpx + self.gas.derivative_settings = {} # reset defaults + + # surface phase + self.surf.TPX = self.surf_tpx + self.surf.set_multiplier(0.) + self.surf.set_multiplier(1., self.rxn_idx) + self.surf.derivative_settings = {"skip-cov-dep": True, "skip-electrochem": True} + + # check stoichiometric coefficient output + for k, v in self.rxn.reactants.items(): + ix = self.sidxs[k] + self.assertEqual(self.r_stoich[ix, self.rxn_idx], v) + for k, v in self.rxn.products.items(): + ix = self.sidxs[k] + self.assertEqual(self.p_stoich[ix, self.rxn_idx], v) + + def test_input(self): + # ensure that correct equation is referenced + self.assertEqual(self.equation, self.rxn.equation) + self.assertEqual(self.rate_type, self.rxn.rate.type) + + def test_forward_rop_ddN(self): + # check derivatives of forward rates of progress with respect to species + # concentrations against analytic result + drop = self.surf.forward_rates_of_progress_ddN + rop = self.surf.forward_rates_of_progress + concentrations = np.concatenate((self.surf.concentrations, self.gas.concentrations)) + specs = self.surf.species_names + self.gas.species_names + for spc_ix in self.rix: + if self.orders is None: + order = self.r_stoich[spc_ix, self.rxn_idx] + else: + order = self.orders[specs[spc_ix]] + self.assertNear(rop[self.rxn_idx], + drop[self.rxn_idx, spc_ix] * concentrations[spc_ix] / order) + + def test_reverse_rop_ddN(self): + # check derivatives of forward rates of progress with respect to species + # concentrations against analytic result + drop = self.surf.reverse_rates_of_progress_ddN + rop = self.surf.reverse_rates_of_progress + concentrations = np.concatenate((self.surf.concentrations, self.gas.concentrations)) + specs = self.surf.species_names + self.gas.species_names + for spc_ix in self.pix: + if self.orders is None: + order = self.p_stoich[spc_ix, self.rxn_idx] + else: + order = self.orders[specs[spc_ix]] + self.assertNear(rop[self.rxn_idx], + drop[self.rxn_idx, spc_ix] * concentrations[spc_ix] / order) + + def test_net_rop_ddN(self): + # check derivatives of net rates of progress with respect to species + # concentrations against analytic + rop = self.surf.net_rates_of_progress + drop = self.surf.net_rates_of_progress_ddN + concentrations = np.concatenate((self.surf.concentrations, self.gas.concentrations)) + + drop *= concentrations + for spc_ix in self.rix + self.pix: + ix = np.abs(drop[:, spc_ix]) > 1 + drop_ = drop[:, spc_ix] + self.assertArrayNear(drop_[ix], rop[ix], 1e-3) + +class PlatinumHydrogen(SurfaceRateExpressionTests): + + @classmethod + def setUpClass(cls): + phase_defs = """ + units: {length: cm, quantity: mol, activation-energy: J/mol} + phases: + - name: gas + thermo: ideal-gas + species: + - gri30.yaml/species: [H2, H2O, H2O2, O2] + kinetics: gas + reactions: + - gri30.yaml/reactions: declared-species + skip-undeclared-third-bodies: true + - name: Pt_surf + thermo: ideal-surface + species: + - ptcombust.yaml/species: [PT(S), H(S), H2O(S), OH(S), O(S)] + kinetics: surface + reactions: [ptcombust.yaml/reactions: declared-species] + site-density: 3e-09 + """ + # create phase objects + cls.gas = ct.Solution(yaml=phase_defs, name="gas") + cls.surf = ct.Interface(yaml=phase_defs, name="Pt_surf", adjacent=[cls.gas]) + cls.gas.TPX = 800, 2*ct.one_atm, "H2:1.5, O2:1.0, H2O2:0.75, H2O:0.3" + cls.surf.TPX = 800, 2*ct.one_atm , "PT(S):4.0, H(S):0.5, H2O(S):0.1, OH(S):0.2, O(S):0.8" + cls.gas_tpx = cls.gas.TPX + cls.surf_tpx = cls.surf.TPX + super().setUpClass() + +class SurfInterfaceArrhenius(PlatinumHydrogen, utilities.CanteraTest): + rxn_idx = 7 + equation = "H(S) + O(S) <=> OH(S) + PT(S)" + rate_type = "interface-Arrhenius" + +class SurfGasFwdStickingArrhenius(PlatinumHydrogen, utilities.CanteraTest): + rxn_idx = 5 + equation = "H2O + PT(S) => H2O(S)" + rate_type = "sticking-Arrhenius" + +class SurfGasInterfaceArrhenius(PlatinumHydrogen, utilities.CanteraTest): + rxn_idx = 0 + equation = "H2 + 2 PT(S) => 2 H(S)" + rate_type = "interface-Arrhenius" + orders = {"PT(S)": 1, "H2": 1, "H(S)": 2} + +class GasSurfInterfaceArrhenius(PlatinumHydrogen, utilities.CanteraTest): + rxn_idx = 6 + equation = "H2O(S) => H2O + PT(S)" + rate_type = "interface-Arrhenius" + +class SurfaceFullTests: + # Generic test class to check derivatives evaluated for an entire reaction mechanisms + rtol = 1e-4 + + @classmethod + def setUpClass(cls): + # all species indexs + all_species = cls.surf.species() + cls.gas.species() + cls.sidxs = {spec.name:i for i, spec in enumerate(all_species)} + + def setUp(self): + # gas phase + self.gas.TPX = self.gas_tpx + self.gas.derivative_settings = {} # reset defaults + # surface phase + self.surf.TPX = self.surf_tpx + self.surf.derivative_settings = {"skip-cov-dep": True, "skip-electrochem": True} + + # closure to get concentrations vector + def get_concentrations(self): + return np.concatenate((self.surf.concentrations, self.gas.concentrations)) + + def test_forward_rop_ddN(self): + # check forward rop against anayltical result + drop = self.surf.forward_rates_of_progress_ddN + rop = self.surf.forward_rates_of_progress + conc = self.get_concentrations() + # multiply derivatives with concentrations + drop = np.matmul(drop, conc) + # get total reactant reaction orders + total_orders = [] + for rxn in self.surf.reactions(): + orders = rxn.orders + curr_order = 0 + for k, v in rxn.reactants.items(): + if k in orders: + curr_order += orders[k] + else: + curr_order += v + total_orders.append(curr_order) + total_orders = np.array(total_orders) + # divide derivatives by told orders which should provide rops + drop /= total_orders + # compare the two + self.assertArrayNear(drop, rop, self.rtol) + + def test_reverse_rop_ddN(self): + # check forward rop against anayltical result + drop = self.surf.reverse_rates_of_progress_ddN + rop = self.surf.reverse_rates_of_progress + conc = self.get_concentrations() + # multiply derivatives with concentrations + drop = np.matmul(drop, conc) + # get total reactant reaction orders + total_orders = [] + for rxn in self.surf.reactions(): + orders = rxn.orders + curr_order = 0 + for k, v in rxn.products.items(): + if k in orders: + curr_order += orders[k] + else: + curr_order += v + total_orders.append(curr_order) + total_orders = np.array(total_orders) + # divide derivatives by told orders which should provide rops + drop /= total_orders + # compare the two + self.assertArrayNear(drop, rop, self.rtol) + + def test_net_rop_ddN(self): + # check forward rop against anayltical result + drop = self.surf.net_rates_of_progress_ddN + rop = self.surf.net_rates_of_progress + conc = self.get_concentrations() + # multiply derivatives with concentrations + drop = np.matmul(drop, conc) + # get total reaction orders + total_orders = [] + for rxn in self.surf.reactions(): + orders = rxn.orders + curr_order = 0 + for k, v in rxn.reactants.items(): + if k in orders: + curr_order += orders[k] + else: + curr_order += v + total_orders.append(curr_order) + total_orders = np.array(total_orders) + # divide derivatives by told orders which should provide rops + drop /= total_orders + # compare the two + self.assertArrayNear(drop, rop, self.rtol) + +class FullPlatinumHydrogen(SurfaceFullTests, utilities.CanteraTest): + + @classmethod + def setUpClass(cls): + phase_defs = """ + units: {length: cm, quantity: mol, activation-energy: J/mol} + phases: + - name: gas + thermo: ideal-gas + species: + - gri30.yaml/species: [H2, H2O, H2O2, O2] + kinetics: gas + reactions: + - gri30.yaml/reactions: declared-species + skip-undeclared-third-bodies: true + - name: Pt_surf + thermo: ideal-surface + species: + - ptcombust.yaml/species: [PT(S), H(S), H2O(S), OH(S), O(S)] + kinetics: surface + reactions: [ptcombust.yaml/reactions: declared-species] + site-density: 3e-09 + """ + # create phase objects + cls.gas = ct.Solution(yaml=phase_defs, name="gas") + cls.surf = ct.Interface(yaml=phase_defs, name="Pt_surf", adjacent=[cls.gas]) + cls.gas.TPX = 800, 2*ct.one_atm, "H2:1.5, O2:1.0, H2O2:0.75, H2O:0.3" + cls.surf.TPX = 800, 2*ct.one_atm , "PT(S):1.0, H(S):0.5, H2O(S):0.1, OH(S):0.2, O(S):0.8" + cls.gas_tpx = cls.gas.TPX + cls.surf_tpx = cls.surf.TPX + super().setUpClass() diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index d2ab3184aa..1342e25841 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -675,7 +675,9 @@ def test_array_properties(self): # Skipped because they are 2D (conversion not implemented) "binary_diff_coeffs", "creation_rates_ddX", "destruction_rates_ddX", "forward_rates_of_progress_ddX", "net_production_rates_ddX", - "net_rates_of_progress_ddX", "reverse_rates_of_progress_ddX" + "net_rates_of_progress_ddX", "reverse_rates_of_progress_ddX", + "net_rates_of_progress_ddN", "forward_rates_of_progress_ddN", + "reverse_rates_of_progress_ddN" } for attr in dir(self.gas): diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 363ec3e65e..b28a722c85 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1944,7 +1944,8 @@ def test_array_properties_exist(self): skip = { # Skipped because they are complicated (conversion not implemented) "forward_rates_of_progress_ddX", "net_rates_of_progress_ddX", - "reverse_rates_of_progress_ddX", "state" + "reverse_rates_of_progress_ddX", "state", "net_rates_of_progress_ddN", + "forward_rates_of_progress_ddN", "reverse_rates_of_progress_ddN" } skip.update(ct.SolutionArray._passthrough)