Generalized preconditioning and surface chemistry jacobian

include/cantera/kinetics/GasKinetics.h

@@ -77,6 +77,9 @@ class GasKinetics : public BulkKinetics
     virtual Eigen::SparseMatrix<double> fwdRatesOfProgress_ddX();
     virtual Eigen::SparseMatrix<double> revRatesOfProgress_ddX();
     virtual Eigen::SparseMatrix<double> netRatesOfProgress_ddX();
+    virtual Eigen::SparseMatrix<double> fwdRatesOfProgress_ddCi();
+    virtual Eigen::SparseMatrix<double> revRatesOfProgress_ddCi();
+    virtual Eigen::SparseMatrix<double> netRatesOfProgress_ddCi();
 
     //! Update temperature-dependent portions of reaction rates and falloff
     //! functions.
@@ -103,14 +106,14 @@ class GasKinetics : public BulkKinetics
     void processThirdBodies(double* rop);
 
     //! Multiply rate with inverse equilibrium constant
-    void processEquilibriumConstants(double* rop);
+    void applyEquilibriumConstants(double* rop);
 
     //! Multiply rate with scaled temperature derivatives of the inverse
     //! equilibrium constant
     /*!
      *  This (scaled) derivative is handled by a finite difference.
      */
-    void processEquilibriumConstants_ddT(double* drkcn);
+    void applyEquilibriumConstants_ddT(double* drkcn);
 
     //! Process temperature derivative
     //! @param in  rate expression used for the derivative calculation
@@ -130,11 +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<double> process_ddX(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<double> calculateCompositionDerivatives(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

include/cantera/kinetics/InterfaceKinetics.h

@@ -305,7 +305,40 @@ class InterfaceKinetics : public Kinetics
     */
     double interfaceCurrent(const size_t iphase);
 
+    virtual void setDerivativeSettings(const AnyMap& settings);
+    virtual void getDerivativeSettings(AnyMap& settings) const;
+    virtual Eigen::SparseMatrix<double> fwdRatesOfProgress_ddCi();
+    virtual Eigen::SparseMatrix<double> revRatesOfProgress_ddCi();
+    virtual Eigen::SparseMatrix<double> netRatesOfProgress_ddCi();
+
 protected:
+    //! @name Internal service methods
+    //!
+    //! @note These methods are for internal use, and seek to avoid code duplication
+    //! while evaluating terms used for rate constants, rates of progress, and
+    //! their derivatives.
+    //! @{
+
+
+    //! Multiply rate with inverse equilibrium constant
+    void applyEquilibriumConstants(double* rop);
+
+    //! Process mole fraction derivative
+    //! @param stoich  stoichiometry manager
+    //! @param in  rate expression used for the derivative calculation
+    //! @return a sparse matrix of derivative contributions for each reaction of
+    //! dimensions nTotalReactions by nTotalSpecies
+    Eigen::SparseMatrix<double> calculateCompositionDerivatives(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 string& name);
+
+    //! @}
+
     //! Temporary work vector of length m_kk
     vector_fp m_grt;
 
@@ -464,6 +497,22 @@ class InterfaceKinetics : public Kinetics
     //! Number of dimensions of reacting phase (2 for InterfaceKinetics, 1 for
     //! EdgeKinetics)
     size_t m_nDim = 2;
+
+    //! Buffers for partial rop results with length nReactions()
+    vector_fp m_rbuf0;
+    vector_fp m_rbuf1;
+
+    //! A flag used to neglect rate coefficient coverage dependence in derivative
+    //! formation
+    bool m_jac_skip_coverage_dependence = false;
+    //! A flag used to neglect electrochemical contributions in derivative formation
+    bool m_jac_skip_electrochemistry = false;
+    //! Relative tolerance used in developing numerical portions of specific derivatives
+    double m_jac_rtol_delta = 1e-8;
+    //! A flag stating if the object uses electrochemistry
+    bool m_has_electrochemistry = false;
+    //! A flag stating if the object has coverage dependent rates
+    bool m_has_coverage_dependence = false;
 };
 
 }

include/cantera/kinetics/InterfaceRate.h

@@ -116,7 +116,7 @@ class InterfaceRateBase
     //! Add a coverage dependency for species *sp*, with exponential dependence
     //! *a*, power-law exponent *m*, and activation energy dependence *e*,
     //! where *e* is in Kelvin, that is, energy divided by the molar gas constant.
-    void addCoverageDependence(const std::string& sp, double a, double m, double e);
+    virtual void addCoverageDependence(const string& sp, double a, double m, double e);
 
     //! Boolean indicating whether rate uses exchange current density formulation
     bool exchangeCurrentDensityFormulation() {
@@ -433,6 +433,11 @@ class InterfaceRate : public RateType, public InterfaceRateBase
     virtual double activationEnergy() const override {
         return RateType::activationEnergy() + m_ecov * GasConstant;
     }
+
+    void addCoverageDependence(const string& sp, double a, double m, double e) override {
+        InterfaceRateBase::addCoverageDependence(sp, a, m, e);
+        RateType::setCompositionDependence(true);
+    }
 };
 
 using InterfaceArrheniusRate = InterfaceRate<ArrheniusRate, InterfaceData>;

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)
+    //!  - `_ddCi`: derivative with respect to species concentrations (a matrix)
     //!
     //! Settings for derivative evaluation are set by keyword/value pairs using
     //! the methods getDerivativeSettings() and setDerivativeSettings().
@@ -601,6 +602,16 @@ class Kinetics
     //!  - `rtol-delta` (double) ... relative tolerance used to perturb properties
     //!    when calculating numerical derivatives. The default value is 1e-8.
     //!
+    //! For InterfaceKinetics, the following keyword/value pairs are supported:
+    //!  - `skip-coverage-dependence` (boolean) ... if `false` (default), rate constant
+    //!   coverage dependence is not considered when evaluating derivatives.
+    //! - `skip-electrochemistry` (boolean) ... if `false` (default), electrical charge
+    //!   is not considered in evaluating the derivatives and these reactions are
+    //!   treated as normal surface reactions.
+    //!  - `rtol-delta` (double) ... relative tolerance used to perturb properties
+    //!    when calculating numerical derivatives. The default value is 1e-8.
+    //!
+    //!
     //! @warning  The calculation of derivatives is an experimental part of the
     //!      %Cantera API and may be changed or removed without notice.
     //! @{
@@ -716,6 +727,26 @@ 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$c_i\f$, all other \f$c_j\f$ are held
+     * constant.
+     *
+     * @warning  This method is an experimental part of the %Cantera API and
+     *      may be changed or removed without notice.
+     *
+     * @since  New in Cantera 3.0.
+     */
+    virtual Eigen::SparseMatrix<double> fwdRatesOfProgress_ddCi()
+    {
+        throw NotImplementedError("Kinetics::fwdRatesOfProgress_ddCi",
+            "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 +800,26 @@ 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$c_i\f$, all other \f$c_j\f$ are held
+     * constant.
+     *
+     * @warning  This method is an experimental part of the %Cantera API and
+     *      may be changed or removed without notice.
+     *
+     * @since  New in Cantera 3.0.
+     */
+    virtual Eigen::SparseMatrix<double> revRatesOfProgress_ddCi()
+    {
+        throw NotImplementedError("Kinetics::revRatesOfProgress_ddCi",
+            "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.
@@ -822,6 +873,26 @@ class Kinetics
             "Not implemented for kinetics type '{}'.", kineticsType());
     }
 
+    /**
+     * Calculate derivatives for net 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$c_i\f$, all other \f$c_j\f$ are held
+     * constant.
+     *
+     * @warning  This method is an experimental part of the %Cantera API and
+     *      may be changed or removed without notice.
+     *
+     * @since  New in Cantera 3.0.
+     */
+    virtual Eigen::SparseMatrix<double> netRatesOfProgress_ddCi()
+    {
+        throw NotImplementedError("Kinetics::netRatesOfProgress_ddCi",
+            "Not implemented for kinetics type '{}'.", kineticsType());
+    }
+
     /**
      * Calculate derivatives for species creation rates with respect to temperature
      * at constant pressure, molar concentration and mole fractions.
@@ -859,6 +930,22 @@ class Kinetics
      */
     Eigen::SparseMatrix<double> creationRates_ddX();
 
+    /**
+     * Calculate derivatives for species creation rates with respect to species
+     * 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$c_i\f$, all other \f$c_j\f$ are held
+     * constant.
+     *
+     * @warning  This method is an experimental part of the %Cantera API and
+     *      may be changed or removed without notice.
+     *
+     * @since  New in Cantera 3.0.
+     */
+    Eigen::SparseMatrix<double> creationRates_ddCi();
+
     /**
      * Calculate derivatives for species destruction rates with respect to temperature
      * at constant pressure, molar concentration and mole fractions.
@@ -896,6 +983,22 @@ class Kinetics
      */
     Eigen::SparseMatrix<double> destructionRates_ddX();
 
+    /**
+     * Calculate derivatives for species destruction rates with respect to species
+     * 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$c_i\f$, all other \f$c_j\f$ are held
+     * constant.
+     *
+     * @warning  This method is an experimental part of the %Cantera API and
+     *      may be changed or removed without notice.
+     *
+     * @since  New in Cantera 3.0.
+     */
+    Eigen::SparseMatrix<double> destructionRates_ddCi();
+
     /**
      * Calculate derivatives for species net production rates with respect to
      * temperature at constant pressure, molar concentration and mole fractions.
@@ -933,6 +1036,22 @@ class Kinetics
      */
     Eigen::SparseMatrix<double> netProductionRates_ddX();
 
+    /**
+     * Calculate derivatives for species net production rates with respect to species
+     * 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$c_i\f$, all other \f$c_j\f$ are held
+     * constant.
+     *
+     * @warning  This method is an experimental part of the %Cantera API and
+     *      may be changed or removed without notice.
+     *
+     * @since  New in Cantera 3.0.
+     */
+    Eigen::SparseMatrix<double> netProductionRates_ddCi();
+
     //! @}
     //! @name Reaction Mechanism Informational Query Routines
     //! @{

include/cantera/kinetics/ReactionRate.h

@@ -190,6 +190,18 @@ class ReactionRate
         return m_valid;
     }
 
+    //! Boolean indicating whether rate has compositional dependence
+    //! @since New in Cantera 3.0
+    bool compositionDependent() {
+        return m_composition_dependent_rate;
+    }
+
+    //! Set rate compositional dependence
+    //! @since New in Cantera 3.0
+    void setCompositionDependence(bool comp_dep) {
+        m_composition_dependent_rate = comp_dep;
+    }
+
 protected:
     //! Get parameters
     //! @param node  AnyMap containing rate information
@@ -209,6 +221,9 @@ class ReactionRate
     //! Flag indicating whether reaction rate is set up correctly
     bool m_valid = false;
 
+    //! Flag indicating composition dependent rate
+    bool m_composition_dependent_rate = false;
+
 private:
     //! Return an object that be used to evaluate the rate by converting general input
     //! such as temperature and pressure into the `DataType` struct that is particular

include/cantera/zeroD/IdealGasConstPressureMoleReactor.h

@@ -45,6 +45,8 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor
     //! between reactors, for example via inlets and outlets.
     virtual Eigen::SparseMatrix<double> jacobian();
 
+    virtual bool preconditionerSupported() const {return true;};
+
 protected:
     vector_fp m_hk; //!< Species molar enthalpies
 };

include/cantera/zeroD/IdealGasMoleReactor.h

@@ -45,6 +45,8 @@ class IdealGasMoleReactor : public MoleReactor
     //! between reactors, for example via inlets and outlets.
     virtual Eigen::SparseMatrix<double> jacobian();
 
+    virtual bool preconditionerSupported() const {return true;};
+
 protected:
     vector_fp m_uk; //!< Species molar internal energies
 };

include/cantera/zeroD/MoleReactor.h

@@ -37,6 +37,11 @@ class MoleReactor : public Reactor
     std::string componentName(size_t k);
 
 protected:
+    //! For each surface in the reactor, update vector of triplets with all relevant
+    //! surface jacobian derivatives of species with respect to species
+    //! which are appropriately offset to align with the reactor's state vector.
+    virtual void addSurfaceJacobian(vector<Eigen::Triplet<double>> &triplets);
+
     //! 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);

include/cantera/zeroD/Reactor.h

@@ -186,6 +186,15 @@ class Reactor : public ReactorBase
     //! Reset the reaction rate multipliers
     virtual void resetSensitivity(double* params);
 
+    //! Return a false if preconditioning is not supported or true otherwise.
+    //!
+    //! @warning  This method is an experimental part of the %Cantera
+    //! API and may be changed or removed without notice.
+    //!
+    //! @since New in Cantera 3.0
+    //!
+    virtual bool preconditionerSupported() const {return false;};
+
 protected:
     //! Return the index in the solution vector for this reactor of the species
     //! named *nm*, in either the homogeneous phase or a surface phase, relative

include/cantera/zeroD/ReactorNet.h

@@ -284,8 +284,7 @@ class ReactorNet : public FuncEval
     virtual void setDerivativeSettings(AnyMap& settings);
 
 protected:
-    //! Check if surfaces and preconditioning are included, if so throw an error because
-    //! they are currently not supported.
+    //! Check that preconditioning is supported by all reactors in the network
     virtual void checkPreconditionerSupported();
 
     //! Update the preconditioner based on the already computed jacobian values

interfaces/cython/cantera/composite.py

@@ -523,7 +523,8 @@ class SolutionArray(SolutionArrayBase):
 
     _n_species2 = [
         'multi_diff_coeffs', 'binary_diff_coeffs', 'creation_rates_ddX',
-        'destruction_rates_ddX', 'net_production_rates_ddX'
+        'destruction_rates_ddX', 'net_production_rates_ddX', 'creation_rates_ddCi',
+        'destruction_rates_ddCi', 'net_production_rates_ddCi'
     ]
 
     _n_reactions = [