Convert kinetics & equilibrium code to use std::span

doc/sphinx/userguide/kinetics_transport.cpp

@@ -20,7 +20,7 @@ void simple_demo2()
 
     // Get the net reaction rates.
     vector<double> wdot(kin->nReactions());
-    kin->getNetRatesOfProgress(wdot.data());
+    kin->getNetRatesOfProgress(wdot);
 
     writelog("Net reaction rates for reactions involving CO2\n");
     size_t kCO2 = gas->speciesIndex("CO2");

include/cantera/cython/funcWrapper.h

@@ -219,6 +219,8 @@ inline int translate_exception()
         PyErr_SetString(PyExc_IndexError, exn.what());
     } catch (const Cantera::NotImplementedError& exn) {
         PyErr_SetString(PyExc_NotImplementedError, exn.what());
+    } catch (const Cantera::ArraySizeError& exn) {
+        PyErr_SetString(PyExc_ValueError, exn.what());
     } catch (const Cantera::CanteraError& exn) {
         PyErr_SetString(pyCanteraError, exn.what());
     } catch (const std::exception& exn) {

include/cantera/cython/kinetics_utils.h

@@ -53,7 +53,9 @@ inline void sparseCscData(const Eigen::SparseMatrix<double>& mat,
     }
 }
 
-#define KIN_1D(FUNC_NAME) ARRAY_FUNC(kin, Kinetics, FUNC_NAME)
+#define KIN_1D(FUNC_NAME) \
+    inline void kin_ ## FUNC_NAME(Cantera::Kinetics* object, std::span<double> data) \
+    { object->FUNC_NAME(data); }
 
 KIN_1D(getFwdRatesOfProgress)
 KIN_1D(getRevRatesOfProgress)

include/cantera/equil/ChemEquil.h

@@ -117,7 +117,7 @@ class ChemEquil
      *     Unsuccessful returns are indicated by a return value of -1 for lack
      *     of convergence or -3 for a singular Jacobian.
      */
-    int equilibrate(ThermoPhase& s, const char* XY, vector<double>& elMoles,
+    int equilibrate(ThermoPhase& s, const char* XY, span<double> elMoles,
                     int loglevel = 0);
 
     /**
@@ -155,16 +155,15 @@ class ChemEquil
      * @f[ \lambda_m/RT @f].
      * @param t temperature in K.
      */
-    void setToEquilState(ThermoPhase& s,
-                         const vector<double>& x, double t);
+    void setToEquilState(ThermoPhase& s, span<const double> x, double t);
 
     //! Estimate the initial mole numbers. This version borrows from the
     //! MultiPhaseEquil solver.
-    int setInitialMoles(ThermoPhase& s, vector<double>& elMoleGoal, int loglevel = 0);
+    int setInitialMoles(ThermoPhase& s, span<double> elMoleGoal, int loglevel = 0);
 
     //! Generate a starting estimate for the element potentials.
-    int estimateElementPotentials(ThermoPhase& s, vector<double>& lambda,
-                                  vector<double>& elMolesGoal, int loglevel = 0);
+    int estimateElementPotentials(ThermoPhase& s, span<double> lambda,
+                                  span<double> elMolesGoal, int loglevel = 0);
 
     /**
      * Do a calculation of the element potentials using the Brinkley method,
@@ -199,7 +198,7 @@ class ChemEquil
      *
      * NOTE: update for activity coefficients.
      */
-    int estimateEP_Brinkley(ThermoPhase& s, vector<double>& lambda, vector<double>& elMoles);
+    int estimateEP_Brinkley(ThermoPhase& s, span<double> lambda, span<double> elMoles);
 
     //! Find an acceptable step size and take it.
     /*!
@@ -214,22 +213,22 @@ class ChemEquil
      * limited to a factor of 2 jump in the values from this method. Near
      * convergence, the delta damping gets out of the way.
      */
-    int dampStep(ThermoPhase& s, vector<double>& oldx,
-                 double oldf, vector<double>& grad, vector<double>& step, vector<double>& x,
-                 double& f, vector<double>& elmols, double xval, double yval);
+    int dampStep(ThermoPhase& s, span<double> oldx, double oldf, span<double> grad,
+                 span<double> step, span<double> x, double& f, span<double> elmols,
+                 double xval, double yval);
 
     /**
      *  Evaluates the residual vector F, of length #m_mm
      */
-    void equilResidual(ThermoPhase& s, const vector<double>& x,
-                       const vector<double>& elmtotal, vector<double>& resid,
+    void equilResidual(ThermoPhase& s, span<const double> x,
+                       span<const double> elmtotal, span<double> resid,
                        double xval, double yval, int loglevel = 0);
 
-    void equilJacobian(ThermoPhase& s, vector<double>& x,
-                       const vector<double>& elmols, DenseMatrix& jac,
+    void equilJacobian(ThermoPhase& s, span<double> x,
+                       span<const double> elmols, DenseMatrix& jac,
                        double xval, double yval, int loglevel = 0);
 
-    void adjustEloc(ThermoPhase& s, vector<double>& elMolesGoal);
+    void adjustEloc(ThermoPhase& s, span<double> elMolesGoal);
 
     //! Update internally stored state information.
     void update(const ThermoPhase& s);
@@ -243,10 +242,9 @@ class ChemEquil
      * @param[in] Xmol_i_calc Mole fractions of the species
      * @param[in] pressureConst Pressure
      */
-    double calcEmoles(ThermoPhase& s, vector<double>& x,
-                      const double& n_t, const vector<double>& Xmol_i_calc,
-                      vector<double>& eMolesCalc, vector<double>& n_i_calc,
-                      double pressureConst);
+    double calcEmoles(ThermoPhase& s, span<double> x, const double& n_t,
+                      span<const double> Xmol_i_calc, span<double> eMolesCalc,
+                      span<double> n_i_calc, double pressureConst);
 
     size_t m_mm; //!< number of elements in the phase
     size_t m_kk; //!< number of species in the phase

include/cantera/equil/MultiPhase.h

@@ -79,7 +79,7 @@ class MultiPhase
      *   @param phases Vector of pointers to phases
      *   @param phaseMoles Vector of mole numbers in each phase (kmol)
      */
-    void addPhases(vector<ThermoPhase*>& phases, const vector<double>& phaseMoles);
+    void addPhases(span<ThermoPhase*> phases, span<const double> phaseMoles);
 
     //! Add all phases present in 'mix' to this mixture.
     /*!
@@ -174,7 +174,7 @@ class MultiPhase
      *
      *    @param x  vector of mole fractions. Length = number of global species.
      */
-    void getMoleFractions(double* const x) const;
+    void getMoleFractions(span<double> x) const;
 
     //! Process phases and build atomic composition array.
     /*!
@@ -342,7 +342,7 @@ class MultiPhase
      *                  potentials are returned instead of the composition-
      *                  dependent chemical potentials.
      */
-    void getValidChemPotentials(double not_mu, double* mu,
+    void getValidChemPotentials(double not_mu, span<double> mu,
                                 bool standard = false) const;
 
     //! Temperature [K].
@@ -403,7 +403,7 @@ class MultiPhase
      * @param Moles  Vector of mole numbers of all the species in all the phases
      *               (kmol)
      */
-    void setState_TPMoles(const double T, const double Pres, const double* Moles);
+    void setState_TPMoles(const double T, const double Pres, span<const double> Moles);
 
     //! Pressure [Pa].
     double pressure() const {
@@ -475,7 +475,7 @@ class MultiPhase
      * @param n    index of the phase
      * @param x    Vector of input mole fractions.
      */
-    void setPhaseMoleFractions(const size_t n, const double* const x);
+    void setPhaseMoleFractions(const size_t n, span<const double> x);
 
     //! Set the number of moles of species in the mixture
     /*!
@@ -500,7 +500,7 @@ class MultiPhase
      * @param[out] molNum Vector of doubles of length nSpecies() containing the
      *               global mole numbers (kmol).
      */
-    void getMoles(double* molNum) const;
+    void getMoles(span<double> molNum) const;
 
     //! Sets all of the global species mole numbers
     /*!
@@ -510,7 +510,7 @@ class MultiPhase
      * @param n    Vector of doubles of length nSpecies() containing the global
      *             mole numbers (kmol).
      */
-    void setMoles(const double* n);
+    void setMoles(span<const double> n);
 
     //! Adds moles of a certain species to the mixture
     /*!
@@ -525,7 +525,7 @@ class MultiPhase
      * Length = number of elements in the MultiPhase object.
      * Index is the global element index. Units is in kmol.
      */
-    void getElemAbundances(double* elemAbundances) const;
+    void getElemAbundances(span<double> elemAbundances) const;
 
     //! Return true if the phase @e p has valid thermo data for the current
     //! temperature.
@@ -723,10 +723,9 @@ std::ostream& operator<<(std::ostream& s, MultiPhase& x);
  *
  * @ingroup equilGroup
  */
-size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
-                     MultiPhase* mphase, vector<size_t>& orderVectorSpecies,
-                     vector<size_t>& orderVectorElements,
-                     vector<double>& formRxnMatrix);
+size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
+                     span<size_t> orderVectorSpecies, span<size_t> orderVectorElements,
+                     span<double> formRxnMatrix);
 
 //! Handles the potential rearrangement of the constraint equations
 //! represented by the Formula Matrix.
@@ -766,10 +765,9 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
  *
  * @ingroup equilGroup
  */
-void ElemRearrange(size_t nComponents, const vector<double>& elementAbundances,
-                   MultiPhase* mphase,
-                   vector<size_t>& orderVectorSpecies,
-                   vector<size_t>& orderVectorElements);
+void ElemRearrange(size_t nComponents, span<const double> elementAbundances,
+                   MultiPhase* mphase, span<size_t> orderVectorSpecies,
+                   span<size_t> orderVectorElements);
 
 //! External int that is used to turn on debug printing for the
 //! BasisOptimize program.

include/cantera/equil/MultiPhaseEquil.h

@@ -51,8 +51,9 @@ class MultiPhaseEquil
         }
     }
 
-    void getStoichVector(size_t rxn, vector<double>& nu) {
-        nu.resize(m_nsp, 0.0);
+    void getStoichVector(size_t rxn, span<double> nu) {
+        checkArraySize("MultiPhaseEquil::getStoichVector", nu.size(), m_nsp);
+        fill(nu.begin(), nu.end(), 0.0);
         if (rxn > nFree()) {
             return;
         }
@@ -104,7 +105,7 @@ class MultiPhaseEquil
     //!     have length = nSpecies(), then the species will be considered as
     //!     candidates to be components in declaration order, beginning with
     //!     the first phase added.
-    void getComponents(const vector<size_t>& order);
+    void getComponents(span<const size_t> order);
 
     //! Estimate the initial mole numbers. This is done by running each
     //! reaction as far forward or backward as possible, subject to the
@@ -124,12 +125,12 @@ class MultiPhaseEquil
 
     //! Re-arrange a vector of species properties in sorted form
     //! (components first) into unsorted, sequential form.
-    void unsort(vector<double>& x);
+    void unsort(span<double> x);
 
-    void step(double omega, vector<double>& deltaN, int loglevel = 0);
+    void step(double omega, span<double> deltaN, int loglevel = 0);
 
     //! Compute the change in extent of reaction for each reaction.
-    double computeReactionSteps(vector<double>& dxi);
+    double computeReactionSteps(span<double> dxi);
 
     void updateMixMoles();
 

include/cantera/equil/vcs_VolPhase.h

@@ -72,7 +72,7 @@ class vcs_VolPhase
      * @param moleFracVec      Vector of input mole fractions
      * @param vcsStateStatus   Status flag for this update
      */
-    void setMoleFractionsState(const double molNum, const double* const moleFracVec,
+    void setMoleFractionsState(const double molNum, span<const double> moleFracVec,
                                const int vcsStateStatus);
 
     //! Set the moles within the phase
@@ -88,7 +88,7 @@ class vcs_VolPhase
      *     to gather the species into the local contiguous vector format.
      */
     void setMolesFromVCS(const int stateCalc,
-                         const double* molesSpeciesVCS = 0);
+                         span<const double> molesSpeciesVCS = {});
 
     //! Set the moles within the phase
     /*!
@@ -108,8 +108,8 @@ class vcs_VolPhase
      * @param TPhMoles  VCS's array containing the number of moles in each phase.
      */
     void setMolesFromVCSCheck(const int vcsStateStatus,
-                              const double* molesSpeciesVCS,
-                              const double* const TPhMoles);
+                              span<const double> molesSpeciesVCS,
+                              span<const double> TPhMoles);
 
     //! Update the moles within the phase, if necessary
     /*!
@@ -134,7 +134,7 @@ class vcs_VolPhase
      *     of the phases in a VCS problem. Only the entries for the current
      *     phase are filled in.
      */
-    void sendToVCS_ActCoeff(const int stateCalc, double* const AC);
+    void sendToVCS_ActCoeff(const int stateCalc, span<double> AC);
 
     //! set the electric potential of the phase
     /*!
@@ -176,7 +176,7 @@ class vcs_VolPhase
      *     all of the phases in a VCS problem. Only the entries for the current
      *     phase are filled in.
      */
-    double sendToVCS_VolPM(double* const VolPM) const;
+    double sendToVCS_VolPM(span<double> VolPM) const;
 
     //! Fill in the standard state Gibbs free energy vector for VCS
     /*!
@@ -188,7 +188,7 @@ class vcs_VolPhase
      *     species in all of the phases in a VCS problem. Only the entries for the
      *     current phase are filled in.
      */
-    void sendToVCS_GStar(double* const gstar) const;
+    void sendToVCS_GStar(span<double> gstar) const;
 
     //! Sets the temperature and pressure in this object and underlying
     //! ThermoPhase objects
@@ -250,11 +250,11 @@ class vcs_VolPhase
      * @param xmol Value of the mole fractions for the species in the phase.
      *             These are contiguous.
      */
-    void setMoleFractions(const double* const xmol);
+    void setMoleFractions(span<const double> xmol);
 
 public:
     //! Return a const reference to the mole fractions stored in the object.
-    const vector<double> & moleFractions() const;
+    span<const double> moleFractions() const;
 
     double moleFraction(size_t klocal) const;
 
@@ -263,13 +263,14 @@ class vcs_VolPhase
      * @param n_k Pointer to a vector of n_k's
      * @param creationGlobalRxnNumbers  Vector of global creation reaction numbers
      */
-    void setCreationMoleNumbers(const double* const n_k, const vector<size_t> &creationGlobalRxnNumbers);
+    void setCreationMoleNumbers(span<const double> n_k,
+                                span<const size_t> creationGlobalRxnNumbers);
 
     //! Return a const reference to the creationMoleNumbers stored in the object.
     /*!
      * @returns a const reference to the vector of creationMoleNumbers
      */
-    const vector<double> & creationMoleNumbers(vector<size_t> &creationGlobalRxnNumbers) const;
+    span<const double> creationMoleNumbers(span<size_t> creationGlobalRxnNumbers) const;
 
     //! Returns whether the phase is an ideal solution phase
     bool isIdealSoln() const;

include/cantera/equil/vcs_internal.h

@@ -51,7 +51,7 @@ class VCS_COUNTERS
  * @param vec vector of doubles
  * @returns   the l2 norm of the vector
  */
-double vcs_l2norm(const vector<double>& vec);
+double vcs_l2norm(span<const double> vec);
 
 //! Returns a const char string representing the type of the species given by
 //! the first argument

include/cantera/equil/vcs_solve.h

@@ -493,7 +493,7 @@ class VCS_SOLVE
      *
      * NOTE: This routine needs to be regulated.
      */
-    void vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS);
+    void vcs_CalcLnActCoeffJac(span<const double> moleSpeciesVCS);
 
     //! Print out a report on the state of the equilibrium problem to
     //! standard output.
@@ -589,7 +589,8 @@ class VCS_SOLVE
     /*!
      * Note, for this algorithm this function should be monotonically decreasing.
      */
-    double vcs_Total_Gibbs(double* w, double* fe, double* tPhMoles);
+    double vcs_Total_Gibbs(span<const double> w, span<const double> fe,
+                           span<const double> tPhMoles);
 
     //! Fully specify the problem to be solved
     void vcs_prob_specifyFully();
@@ -646,14 +647,13 @@ class VCS_SOLVE
      * Make sure to conserve elements and keep track of the total kmoles in
      * all phases.
      *
-     * @param kspec The species index
-     * @param delta_ptr   pointer to the delta for the species. This may
-     *                     change during the calculation
+     * @param kspec  The species index
+     * @param delta  The delta for the species. This may change during the calculation.
      * @return
      *     1: succeeded without change of dx
      *     0: Had to adjust dx, perhaps to zero, in order to do the delta.
      */
-    int delta_species(const size_t kspec, double* const delta_ptr);
+    int delta_species(const size_t kspec, double& delta);
 
     //! Provide an estimate for the deleted species in phases that are not
     //! zeroed out
@@ -751,7 +751,8 @@ class VCS_SOLVE
      */
     double l2normdg(double dg[]) const;
 
-    void checkDelta1(double* const ds, double* const delTPhMoles, size_t kspec);
+    void checkDelta1(span<const double> ds, span<const double> delTPhMoles,
+                     size_t kspec);
 
     //! Calculate the status of single species phases.
     void vcs_SSPhase();