diff --git a/doc/sphinx/userguide/kinetics_transport.cpp b/doc/sphinx/userguide/kinetics_transport.cpp index 7798c86af18..c6da34a0b45 100644 --- a/doc/sphinx/userguide/kinetics_transport.cpp +++ b/doc/sphinx/userguide/kinetics_transport.cpp @@ -20,7 +20,7 @@ void simple_demo2() // Get the net reaction rates. vector 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"); diff --git a/include/cantera/cython/funcWrapper.h b/include/cantera/cython/funcWrapper.h index aba694e4203..d18919c02b8 100644 --- a/include/cantera/cython/funcWrapper.h +++ b/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) { diff --git a/include/cantera/cython/kinetics_utils.h b/include/cantera/cython/kinetics_utils.h index c990db42ac9..305d9137360 100644 --- a/include/cantera/cython/kinetics_utils.h +++ b/include/cantera/cython/kinetics_utils.h @@ -53,7 +53,9 @@ inline void sparseCscData(const Eigen::SparseMatrix& 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 data) \ + { object->FUNC_NAME(data); } KIN_1D(getFwdRatesOfProgress) KIN_1D(getRevRatesOfProgress) diff --git a/include/cantera/equil/ChemEquil.h b/include/cantera/equil/ChemEquil.h index cfd59d3ce2f..d1af1480af2 100644 --- a/include/cantera/equil/ChemEquil.h +++ b/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& elMoles, + int equilibrate(ThermoPhase& s, const char* XY, span 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& x, double t); + void setToEquilState(ThermoPhase& s, span x, double t); //! Estimate the initial mole numbers. This version borrows from the //! MultiPhaseEquil solver. - int setInitialMoles(ThermoPhase& s, vector& elMoleGoal, int loglevel = 0); + int setInitialMoles(ThermoPhase& s, span elMoleGoal, int loglevel = 0); //! Generate a starting estimate for the element potentials. - int estimateElementPotentials(ThermoPhase& s, vector& lambda, - vector& elMolesGoal, int loglevel = 0); + int estimateElementPotentials(ThermoPhase& s, span lambda, + span 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& lambda, vector& elMoles); + int estimateEP_Brinkley(ThermoPhase& s, span lambda, span 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& oldx, - double oldf, vector& grad, vector& step, vector& x, - double& f, vector& elmols, double xval, double yval); + int dampStep(ThermoPhase& s, span oldx, double oldf, span grad, + span step, span x, double& f, span elmols, + double xval, double yval); /** * Evaluates the residual vector F, of length #m_mm */ - void equilResidual(ThermoPhase& s, const vector& x, - const vector& elmtotal, vector& resid, + void equilResidual(ThermoPhase& s, span x, + span elmtotal, span resid, double xval, double yval, int loglevel = 0); - void equilJacobian(ThermoPhase& s, vector& x, - const vector& elmols, DenseMatrix& jac, + void equilJacobian(ThermoPhase& s, span x, + span elmols, DenseMatrix& jac, double xval, double yval, int loglevel = 0); - void adjustEloc(ThermoPhase& s, vector& elMolesGoal); + void adjustEloc(ThermoPhase& s, span 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& x, - const double& n_t, const vector& Xmol_i_calc, - vector& eMolesCalc, vector& n_i_calc, - double pressureConst); + double calcEmoles(ThermoPhase& s, span x, const double& n_t, + span Xmol_i_calc, span eMolesCalc, + span 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 diff --git a/include/cantera/equil/MultiPhase.h b/include/cantera/equil/MultiPhase.h index df3235e1d6f..baefa8b4f1a 100644 --- a/include/cantera/equil/MultiPhase.h +++ b/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& phases, const vector& phaseMoles); + void addPhases(span phases, span 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 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 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 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 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 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 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 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& orderVectorSpecies, - vector& orderVectorElements, - vector& formRxnMatrix); +size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, + span orderVectorSpecies, span orderVectorElements, + span 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& elementAbundances, - MultiPhase* mphase, - vector& orderVectorSpecies, - vector& orderVectorElements); +void ElemRearrange(size_t nComponents, span elementAbundances, + MultiPhase* mphase, span orderVectorSpecies, + span orderVectorElements); //! External int that is used to turn on debug printing for the //! BasisOptimize program. diff --git a/include/cantera/equil/MultiPhaseEquil.h b/include/cantera/equil/MultiPhaseEquil.h index 2962f3a8b38..aae3d29002a 100644 --- a/include/cantera/equil/MultiPhaseEquil.h +++ b/include/cantera/equil/MultiPhaseEquil.h @@ -51,8 +51,9 @@ class MultiPhaseEquil } } - void getStoichVector(size_t rxn, vector& nu) { - nu.resize(m_nsp, 0.0); + void getStoichVector(size_t rxn, span 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& order); + void getComponents(span 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& x); + void unsort(span x); - void step(double omega, vector& deltaN, int loglevel = 0); + void step(double omega, span deltaN, int loglevel = 0); //! Compute the change in extent of reaction for each reaction. - double computeReactionSteps(vector& dxi); + double computeReactionSteps(span dxi); void updateMixMoles(); diff --git a/include/cantera/equil/vcs_VolPhase.h b/include/cantera/equil/vcs_VolPhase.h index 97d32cc451f..fb906d4b901 100644 --- a/include/cantera/equil/vcs_VolPhase.h +++ b/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 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 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 molesSpeciesVCS, + span 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 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 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 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 xmol); public: //! Return a const reference to the mole fractions stored in the object. - const vector & moleFractions() const; + span 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 &creationGlobalRxnNumbers); + void setCreationMoleNumbers(span n_k, + span creationGlobalRxnNumbers); //! Return a const reference to the creationMoleNumbers stored in the object. /*! * @returns a const reference to the vector of creationMoleNumbers */ - const vector & creationMoleNumbers(vector &creationGlobalRxnNumbers) const; + span creationMoleNumbers(span creationGlobalRxnNumbers) const; //! Returns whether the phase is an ideal solution phase bool isIdealSoln() const; diff --git a/include/cantera/equil/vcs_internal.h b/include/cantera/equil/vcs_internal.h index 8c2e9f32ae7..0c91ad5a515 100644 --- a/include/cantera/equil/vcs_internal.h +++ b/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& vec); +double vcs_l2norm(span vec); //! Returns a const char string representing the type of the species given by //! the first argument diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index dde77ad375c..34c19643ce3 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/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 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 w, span fe, + span 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 ds, span delTPhMoles, + size_t kspec); //! Calculate the status of single species phases. void vcs_SSPhase(); diff --git a/include/cantera/kinetics/BulkKinetics.h b/include/cantera/kinetics/BulkKinetics.h index 87e8fd00a78..43d83af3c17 100644 --- a/include/cantera/kinetics/BulkKinetics.h +++ b/include/cantera/kinetics/BulkKinetics.h @@ -42,35 +42,35 @@ class BulkKinetics : public Kinetics //! @name Reaction rate constants, rates of progress, and thermodynamic properties //! @{ - void getFwdRateConstants(double* kfwd) override; - void getEquilibriumConstants(double* kc) override; - void getRevRateConstants(double* krev, bool doIrreversible=false) override; + void getFwdRateConstants(span kfwd) override; + void getEquilibriumConstants(span kc) override; + void getRevRateConstants(span krev, bool doIrreversible=false) override; - void getDeltaGibbs(double* deltaG) override; - void getDeltaEnthalpy(double* deltaH) override; - void getDeltaEntropy(double* deltaS) override; + void getDeltaGibbs(span deltaG) override; + void getDeltaEnthalpy(span deltaH) override; + void getDeltaEntropy(span deltaS) override; - void getDeltaSSGibbs(double* deltaG) override; - void getDeltaSSEnthalpy(double* deltaH) override; - void getDeltaSSEntropy(double* deltaS) override; + void getDeltaSSGibbs(span deltaG) override; + void getDeltaSSEnthalpy(span deltaH) override; + void getDeltaSSEntropy(span deltaS) override; //! @} //! @name Derivatives of rate constants and rates of progress //! @{ void getDerivativeSettings(AnyMap& settings) const override; void setDerivativeSettings(const AnyMap& settings) override; - void getFwdRateConstants_ddT(double* dkfwd) override; - void getFwdRatesOfProgress_ddT(double* drop) override; - void getRevRatesOfProgress_ddT(double* drop) override; - void getNetRatesOfProgress_ddT(double* drop) override; - void getFwdRateConstants_ddP(double* dkfwd) override; - void getFwdRatesOfProgress_ddP(double* drop) override; - void getRevRatesOfProgress_ddP(double* drop) override; - void getNetRatesOfProgress_ddP(double* drop) override; - void getFwdRateConstants_ddC(double* dkfwd) override; - void getFwdRatesOfProgress_ddC(double* drop) override; - void getRevRatesOfProgress_ddC(double* drop) override; - void getNetRatesOfProgress_ddC(double* drop) override; + void getFwdRateConstants_ddT(span dkfwd) override; + void getFwdRatesOfProgress_ddT(span drop) override; + void getRevRatesOfProgress_ddT(span drop) override; + void getNetRatesOfProgress_ddT(span drop) override; + void getFwdRateConstants_ddP(span dkfwd) override; + void getFwdRatesOfProgress_ddP(span drop) override; + void getRevRatesOfProgress_ddP(span drop) override; + void getNetRatesOfProgress_ddP(span drop) override; + void getFwdRateConstants_ddC(span dkfwd) override; + void getFwdRatesOfProgress_ddC(span drop) override; + void getRevRatesOfProgress_ddC(span drop) override; + void getNetRatesOfProgress_ddC(span drop) override; Eigen::SparseMatrix fwdRatesOfProgress_ddX() override; Eigen::SparseMatrix revRatesOfProgress_ddX() override; Eigen::SparseMatrix netRatesOfProgress_ddX() override; @@ -84,8 +84,8 @@ class BulkKinetics : public Kinetics void updateROP() override; - void getThirdBodyConcentrations(double* concm) override; - const vector& thirdBodyConcentrations() const override { + void getThirdBodyConcentrations(span concm) override; + span thirdBodyConcentrations() const override { return m_concm; } @@ -100,35 +100,35 @@ class BulkKinetics : public Kinetics //! @{ //! Multiply rate with third-body collider concentrations - void processThirdBodies(double* rop); + void processThirdBodies(span rop); //! Multiply rate with inverse equilibrium constant - void applyEquilibriumConstants(double* rop); + void applyEquilibriumConstants(span rop); //! Multiply rate with scaled temperature derivatives of the inverse //! equilibrium constant /*! * This (scaled) derivative is handled by a finite difference. */ - void applyEquilibriumConstants_ddT(double* drkcn); + void applyEquilibriumConstants_ddT(span drkcn); //! Process temperature derivative //! @param in rate expression used for the derivative calculation //! @param drop pointer to output buffer - void process_ddT(const vector& in, double* drop); + void process_ddT(const span in, span drop); //! Process pressure derivative //! @param in rate expression used for the derivative calculation //! @param drop pointer to output buffer - void process_ddP(const vector& in, double* drop); + void process_ddP(const span in, span drop); //! Process concentration (molar density) derivative //! @param stoich stoichiometry manager //! @param in rate expression used for the derivative calculation //! @param drop pointer to output buffer //! @param mass_action boolean indicating whether law of mass action applies - void process_ddC(StoichManagerN& stoich, const vector& in, - double* drop, bool mass_action=true); + void process_ddC(StoichManagerN& stoich, span in, + span drop, bool mass_action=true); //! Process derivatives //! @param stoich stoichiometry manager @@ -137,7 +137,7 @@ class BulkKinetics : public Kinetics //! @return a sparse matrix of derivative contributions for each reaction of //! dimensions nTotalReactions by nTotalSpecies Eigen::SparseMatrix calculateCompositionDerivatives( - StoichManagerN& stoich, const vector& in, bool ddX=true); + StoichManagerN& stoich, span 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/ElectronCollisionPlasmaRate.h b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h index fb0a2156eca..1cddd48df41 100644 --- a/include/cantera/kinetics/ElectronCollisionPlasmaRate.h +++ b/include/cantera/kinetics/ElectronCollisionPlasmaRate.h @@ -194,22 +194,22 @@ class ElectronCollisionPlasmaRate : public ReactionRate } //! The value of #m_energyLevels [eV] - const vector& energyLevels() const { + span energyLevels() const { return m_energyLevels; } //! The value of #m_crossSections [m2] - const vector& crossSections() const { + span crossSections() const { return m_crossSections; } //! The value of #m_crossSectionsInterpolated [m2] - const vector& crossSectionInterpolated() const { + span crossSectionInterpolated() const { return m_crossSectionsInterpolated; } //! Update the value of #m_crossSectionsInterpolated [m2] - void updateInterpolatedCrossSection(const vector&); + void updateInterpolatedCrossSection(span); private: //! The name of the kind of electron collision diff --git a/include/cantera/kinetics/Falloff.h b/include/cantera/kinetics/Falloff.h index 0ddf44c44fd..8fff7f218de 100644 --- a/include/cantera/kinetics/Falloff.h +++ b/include/cantera/kinetics/Falloff.h @@ -87,15 +87,17 @@ class FalloffRate : public ReactionRate * @param c Vector of coefficients of the parameterization. The number and * meaning of these coefficients is subclass-dependent. */ - virtual void setFalloffCoeffs(const vector& c); + virtual void setFalloffCoeffs(span c); /** * Retrieve coefficients of the falloff parameterization. * - * @param c Vector of coefficients of the parameterization. The number and - * meaning of these coefficients is subclass-dependent. + * @param c Vector of coefficients of the parameterization. The number (obtained + * with the nParameters() method) and meaning of these coefficients is + * subclass-dependent. */ - virtual void getFalloffCoeffs(vector& c) const; + virtual void getFalloffCoeffs(span c) const {}; + /** * Update the temperature-dependent portions of the falloff function, if @@ -104,7 +106,7 @@ class FalloffRate : public ReactionRate * @param T Temperature [K]. * @param work storage space for intermediate results. */ - virtual void updateTemp(double T, double* work) const {} + virtual void updateTemp(double T, span work) const {} /** * The falloff function. @@ -115,19 +117,19 @@ class FalloffRate : public ReactionRate * to updateTemp. * @returns the value of the falloff function @f$ F @f$ defined above */ - virtual double F(double pr, const double* work) const { + virtual double F(double pr, span work) const { return 1.0; } //! Evaluate falloff function at current conditions double evalF(double T, double conc3b) { - updateTemp(T, m_work.data()); + updateTemp(T, m_work); double logT = std::log(T); double recipT = 1. / T; m_rc_low = m_lowRate.evalRate(logT, recipT); m_rc_high = m_highRate.evalRate(logT, recipT); double pr = conc3b * m_rc_low / (m_rc_high + SmallNumber); - return F(pr, m_work.data()); + return F(pr, m_work); } const string type() const override { @@ -150,7 +152,7 @@ class FalloffRate : public ReactionRate //! Evaluate reaction rate //! @param shared_data data shared by all reactions of a given type double evalFromStruct(const FalloffData& shared_data) { - updateTemp(shared_data.temperature, m_work.data()); + updateTemp(shared_data.temperature, m_work); m_rc_low = m_lowRate.evalRate(shared_data.logT, shared_data.recipT); m_rc_high = m_highRate.evalRate(shared_data.logT, shared_data.recipT); double thirdBodyConcentration; @@ -164,12 +166,12 @@ class FalloffRate : public ReactionRate // Apply falloff function if (m_chemicallyActivated) { // 1 / (1 + Pr) * F - pr = F(pr, m_work.data()) / (1.0 + pr); + pr = F(pr, m_work) / (1.0 + pr); return pr * m_rc_low; } // Pr / (1 + Pr) * F - pr *= F(pr, m_work.data()) / (1.0 + pr); + pr *= F(pr, m_work) / (1.0 + pr); return pr * m_rc_high; } @@ -251,7 +253,7 @@ class LindemannRate final : public FalloffRate LindemannRate(const AnyMap& node, const UnitStack& rate_units={}); LindemannRate(const ArrheniusRate& low, const ArrheniusRate& high, - const vector& c); + span c); unique_ptr newMultiRate() const override{ return make_unique>(); @@ -315,7 +317,7 @@ class TroeRate final : public FalloffRate TroeRate(const AnyMap& node, const UnitStack& rate_units={}); TroeRate(const ArrheniusRate& low, const ArrheniusRate& high, - const vector& c); + span c); unique_ptr newMultiRate() const override { return make_unique>(); @@ -326,9 +328,9 @@ class TroeRate final : public FalloffRate * @param c Vector of three or four doubles: The doubles are the parameters, * a, T_3, T_1, and (optionally) T_2 of the Troe parameterization */ - void setFalloffCoeffs(const vector& c) override; + void setFalloffCoeffs(span c) override; - void getFalloffCoeffs(vector& c) const override; + void getFalloffCoeffs(span c) const override; //! Update the temperature parameters in the representation /*! @@ -336,9 +338,9 @@ class TroeRate final : public FalloffRate * @param work Vector of working space, length 1, representing the * temperature-dependent part of the parameterization. */ - void updateTemp(double T, double* work) const override; + void updateTemp(double T, span work) const override; - double F(double pr, const double* work) const override; + double F(double pr, span work) const override; const string subType() const override { return "Troe"; @@ -409,7 +411,7 @@ class SriRate final : public FalloffRate SriRate(const AnyMap& node, const UnitStack& rate_units={}); - SriRate(const ArrheniusRate& low, const ArrheniusRate& high, const vector& c) + SriRate(const ArrheniusRate& low, const ArrheniusRate& high, span c) : SriRate() { m_lowRate = low; @@ -427,9 +429,9 @@ class SriRate final : public FalloffRate * a, b, c, d (optional; default 1.0), and e (optional; default * 0.0) of the SRI parameterization */ - void setFalloffCoeffs(const vector& c) override; + void setFalloffCoeffs(span c) override; - void getFalloffCoeffs(vector& c) const override; + void getFalloffCoeffs(span c) const override; //! Update the temperature parameters in the representation /*! @@ -437,9 +439,9 @@ class SriRate final : public FalloffRate * @param work Vector of working space, length 2, representing the * temperature-dependent part of the parameterization. */ - void updateTemp(double T, double* work) const override; + void updateTemp(double T, span work) const override; - double F(double pr, const double* work) const override; + double F(double pr, span work) const override; const string subType() const override { return "SRI"; @@ -503,7 +505,7 @@ class TsangRate final : public FalloffRate TsangRate(const AnyMap& node, const UnitStack& rate_units={}); - TsangRate(const ArrheniusRate& low, const ArrheniusRate& high, const vector& c) + TsangRate(const ArrheniusRate& low, const ArrheniusRate& high, span c) : TsangRate() { m_lowRate = low; @@ -520,9 +522,9 @@ class TsangRate final : public FalloffRate * @param c Vector of one or two doubles: The doubles are the parameters, * a and (optionally) b of the Tsang F_cent parameterization */ - void setFalloffCoeffs(const vector& c) override; + void setFalloffCoeffs(span c) override; - void getFalloffCoeffs(vector& c) const override; + void getFalloffCoeffs(span c) const override; //! Update the temperature parameters in the representation /*! @@ -530,9 +532,9 @@ class TsangRate final : public FalloffRate * @param work Vector of working space, length 1, representing the * temperature-dependent part of the parameterization. */ - void updateTemp(double T, double* work) const override; + void updateTemp(double T, span work) const override; - double F(double pr, const double* work) const override; + double F(double pr, span work) const override; const string subType() const override { return "Tsang"; diff --git a/include/cantera/kinetics/Group.h b/include/cantera/kinetics/Group.h index a3f1bea2bb7..64a4abb6400 100644 --- a/include/cantera/kinetics/Group.h +++ b/include/cantera/kinetics/Group.h @@ -24,12 +24,16 @@ class Group Group(size_t n) : m_sign(0) { m_comp.resize(n,0); } - Group(const vector& elnumbers) : - m_comp(elnumbers), m_sign(0) { + Group(span elnumbers) + : m_comp(elnumbers.begin() + , elnumbers.end()) + , m_sign(0) + { validate(); } - Group(const vector& elnumbers) : - m_comp(elnumbers.size()), m_sign(0) { + Group(span elnumbers) : + m_comp(elnumbers.size()), m_sign(0) + { for (size_t i = 0; i < elnumbers.size(); i++) { m_comp[i] = int(elnumbers[i]); } @@ -125,7 +129,7 @@ class Group return m_comp[m]; } - std::ostream& fmt(std::ostream& s, const vector& esymbols) const; + std::ostream& fmt(std::ostream& s, span esymbols) const; private: vector m_comp; diff --git a/include/cantera/kinetics/InterfaceKinetics.h b/include/cantera/kinetics/InterfaceKinetics.h index 874857e5ada..14ebdee2026 100644 --- a/include/cantera/kinetics/InterfaceKinetics.h +++ b/include/cantera/kinetics/InterfaceKinetics.h @@ -81,24 +81,24 @@ class InterfaceKinetics : public Kinetics * where deltaG is the electrochemical potential difference between * products minus reactants. */ - void getEquilibriumConstants(double* kc) override; + void getEquilibriumConstants(span kc) override; - void getDeltaGibbs(double* deltaG) override; + void getDeltaGibbs(span deltaG) override; - void getDeltaElectrochemPotentials(double* deltaM) override; - void getDeltaEnthalpy(double* deltaH) override; - void getDeltaEntropy(double* deltaS) override; + void getDeltaElectrochemPotentials(span deltaM) override; + void getDeltaEnthalpy(span deltaH) override; + void getDeltaEntropy(span deltaS) override; - void getDeltaSSGibbs(double* deltaG) override; - void getDeltaSSEnthalpy(double* deltaH) override; - void getDeltaSSEntropy(double* deltaS) override; + void getDeltaSSGibbs(span deltaG) override; + void getDeltaSSEnthalpy(span deltaH) override; + void getDeltaSSEntropy(span deltaS) override; //! @} //! @name Reaction Mechanism Informational Query Routines //! @{ - void getFwdRateConstants(double* kfwd) override; - void getRevRateConstants(double* krev, bool doIrreversible=false) override; + void getFwdRateConstants(span kfwd) override; + void getRevRateConstants(span krev, bool doIrreversible=false) override; //! @} //! @name Reaction Mechanism Construction @@ -286,15 +286,15 @@ class InterfaceKinetics : public Kinetics //! Multiply rate with inverse equilibrium constant - void applyEquilibriumConstants(double* rop); + void applyEquilibriumConstants(span 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 calculateCompositionDerivatives(StoichManagerN& stoich, - const vector& in); + Eigen::SparseMatrix calculateCompositionDerivatives( + StoichManagerN& stoich, span in); //! Helper function ensuring that all rate derivatives can be calculated //! @param name method name used for error output diff --git a/include/cantera/kinetics/InterfaceRate.h b/include/cantera/kinetics/InterfaceRate.h index 85a3fe9022b..9cc801d1f7c 100644 --- a/include/cantera/kinetics/InterfaceRate.h +++ b/include/cantera/kinetics/InterfaceRate.h @@ -31,7 +31,7 @@ struct InterfaceData : public BlowersMaselData InterfaceData() = default; bool update(const ThermoPhase& bulk, const Kinetics& kin) override; void update(double T) override; - void update(double T, const vector& values) override; + void update(double T, span values) override; using BlowersMaselData::update; virtual void perturbTemperature(double deltaT); void resize(Kinetics& kin) override; @@ -103,7 +103,7 @@ class InterfaceRateBase //! *a*, power-law exponent *m*, and activation energy dependence *e*, //! where *e* is in Kelvin, that is, energy divided by the molar gas constant. virtual void addCoverageDependence(const string& sp, double a, double m, - const vector& e); + span e); //! Boolean indicating whether rate uses exchange current density formulation bool exchangeCurrentDensityFormulation() { @@ -117,7 +117,7 @@ class InterfaceRateBase //! Set association with an ordered list of all species associated with a given //! `Kinetics` object. - void setSpecies(const vector& species); + void setSpecies(span species); //! Update reaction rate parameters //! @param shared_data data shared by all reactions of a given type @@ -426,7 +426,7 @@ class InterfaceRate : public RateType, public InterfaceRateBase } void addCoverageDependence(const string& sp, double a, double m, - const vector& e) override + span e) override { InterfaceRateBase::addCoverageDependence(sp, a, m, e); RateType::setCompositionDependence(true); diff --git a/include/cantera/kinetics/Kinetics.h b/include/cantera/kinetics/Kinetics.h index 24a6f04923b..7791f4109e7 100644 --- a/include/cantera/kinetics/Kinetics.h +++ b/include/cantera/kinetics/Kinetics.h @@ -93,8 +93,8 @@ class AnyMap; //! phase 'a', another bulk phase 'b', and the surface phase 'a:b' at the a/b //! interface. Phase 'a' contains 12 species, phase 'b' contains 3, and at the //! interface there are 5 adsorbed species defined in phase 'a:b'. Then methods -//! like getNetProductionRates(double* net) will write and output array of -//! length 20, beginning at the location pointed to by 'net'. The first 12 +//! like getNetProductionRates(span net) will write and output array of +//! length 20, beginning at the start of the span. The first 12 //! values will be the net production rates for all 12 species of phase 'a' //! (even if some do not participate in the reactions), the next 3 will be for //! phase 'b', and finally the net production rates for the surface species will @@ -347,7 +347,7 @@ class Kinetics * @param fwdROP Output vector containing forward rates * of progress of the reactions. Length: nReactions(). */ - virtual void getFwdRatesOfProgress(double* fwdROP); + virtual void getFwdRatesOfProgress(span fwdROP); //! Return the Reverse rates of progress of the reactions /*! @@ -357,7 +357,7 @@ class Kinetics * @param revROP Output vector containing reverse rates * of progress of the reactions. Length: nReactions(). */ - virtual void getRevRatesOfProgress(double* revROP); + virtual void getRevRatesOfProgress(span revROP); /** * Net rates of progress. Return the net (forward - reverse) rates of @@ -366,7 +366,7 @@ class Kinetics * * @param netROP Output vector of the net ROP. Length: nReactions(). */ - virtual void getNetRatesOfProgress(double* netROP); + virtual void getNetRatesOfProgress(span netROP); //! Return a vector of Equilibrium constants. /*! @@ -381,7 +381,7 @@ class Kinetics * @param kc Output vector containing the equilibrium constants. * Length: nReactions(). */ - virtual void getEquilibriumConstants(double* kc) { + virtual void getEquilibriumConstants(span kc) { throw NotImplementedError("Kinetics::getEquilibriumConstants"); } @@ -399,7 +399,8 @@ class Kinetics * @param property Input vector of property value. Length: #m_kk. * @param deltaProperty Output vector of deltaRxn. Length: nReactions(). */ - virtual void getReactionDelta(const double* property, double* deltaProperty) const; + virtual void getReactionDelta(span property, + span deltaProperty) const; /** * Given an array of species properties 'g', return in array 'dg' the @@ -411,7 +412,7 @@ class Kinetics * primarily designed for use in calculating reverse rate coefficients * from thermochemistry for reversible reactions. */ - virtual void getRevReactionDelta(const double* g, double* dg) const; + virtual void getRevReactionDelta(span g, span dg) const; //! Return the vector of values for the reaction Gibbs free energy change. /*! @@ -423,7 +424,7 @@ class Kinetics * @param deltaG Output vector of deltaG's for reactions Length: * nReactions(). */ - virtual void getDeltaGibbs(double* deltaG) { + virtual void getDeltaGibbs(span deltaG) { throw NotImplementedError("Kinetics::getDeltaGibbs"); } @@ -438,7 +439,7 @@ class Kinetics * @param deltaM Output vector of deltaM's for reactions Length: * nReactions(). */ - virtual void getDeltaElectrochemPotentials(double* deltaM) { + virtual void getDeltaElectrochemPotentials(span deltaM) { throw NotImplementedError("Kinetics::getDeltaElectrochemPotentials"); } @@ -451,7 +452,7 @@ class Kinetics * @param deltaH Output vector of deltaH's for reactions Length: * nReactions(). */ - virtual void getDeltaEnthalpy(double* deltaH) { + virtual void getDeltaEnthalpy(span deltaH) { throw NotImplementedError("Kinetics::getDeltaEnthalpy"); } @@ -464,7 +465,7 @@ class Kinetics * @param deltaS Output vector of deltaS's for reactions Length: * nReactions(). */ - virtual void getDeltaEntropy(double* deltaS) { + virtual void getDeltaEntropy(span deltaS) { throw NotImplementedError("Kinetics::getDeltaEntropy"); } @@ -478,7 +479,7 @@ class Kinetics * @param deltaG Output vector of ss deltaG's for reactions Length: * nReactions(). */ - virtual void getDeltaSSGibbs(double* deltaG) { + virtual void getDeltaSSGibbs(span deltaG) { throw NotImplementedError("Kinetics::getDeltaSSGibbs"); } @@ -492,7 +493,7 @@ class Kinetics * @param deltaH Output vector of ss deltaH's for reactions Length: * nReactions(). */ - virtual void getDeltaSSEnthalpy(double* deltaH) { + virtual void getDeltaSSEnthalpy(span deltaH) { throw NotImplementedError("Kinetics::getDeltaSSEnthalpy"); } @@ -506,7 +507,7 @@ class Kinetics * @param deltaS Output vector of ss deltaS's for reactions Length: * nReactions(). */ - virtual void getDeltaSSEntropy(double* deltaS) { + virtual void getDeltaSSEntropy(span deltaS) { throw NotImplementedError("Kinetics::getDeltaSSEntropy"); } @@ -518,7 +519,7 @@ class Kinetics * @param concm Output vector of effective third-body concentrations. * Length: nReactions(). */ - virtual void getThirdBodyConcentrations(double* concm) { + virtual void getThirdBodyConcentrations(span concm) { throw NotImplementedError("Kinetics::getThirdBodyConcentrations", "Not applicable/implemented for Kinetics object of type '{}'", kineticsType()); @@ -528,10 +529,11 @@ class Kinetics * Provide direct access to current third-body concentration values. * @see getThirdBodyConcentrations. */ - virtual const vector& thirdBodyConcentrations() const { + virtual span thirdBodyConcentrations() const { throw NotImplementedError("Kinetics::thirdBodyConcentrations", "Not applicable/implemented for Kinetics object of type '{}'", kineticsType()); + return {}; } //! @} @@ -545,7 +547,7 @@ class Kinetics * * @param cdot Output vector of creation rates. Length: #m_kk. */ - virtual void getCreationRates(double* cdot); + virtual void getCreationRates(span cdot); /** * Species destruction rates [kmol/m^3/s or kmol/m^2/s]. Return the species @@ -554,7 +556,7 @@ class Kinetics * * @param ddot Output vector of destruction rates. Length: #m_kk. */ - virtual void getDestructionRates(double* ddot); + virtual void getDestructionRates(span ddot); /** * Species net production rates [kmol/m^3/s or kmol/m^2/s]. Return the @@ -564,7 +566,7 @@ class Kinetics * * @param wdot Output vector of net production rates. Length: #m_kk. */ - virtual void getNetProductionRates(double* wdot); + virtual void getNetProductionRates(span wdot); //! @} @@ -714,7 +716,7 @@ class Kinetics * at constant pressure, molar concentration and mole fractions * @param[out] dkfwd Output vector of derivatives. Length: nReactions(). */ - virtual void getFwdRateConstants_ddT(double* dkfwd) + virtual void getFwdRateConstants_ddT(span dkfwd) { throw NotImplementedError("Kinetics::getFwdRateConstants_ddT", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -725,7 +727,7 @@ class Kinetics * at constant temperature, molar concentration and mole fractions. * @param[out] dkfwd Output vector of derivatives. Length: nReactions(). */ - virtual void getFwdRateConstants_ddP(double* dkfwd) + virtual void getFwdRateConstants_ddP(span dkfwd) { throw NotImplementedError("Kinetics::getFwdRateConstants_ddP", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -739,7 +741,7 @@ class Kinetics * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - virtual void getFwdRateConstants_ddC(double* dkfwd) + virtual void getFwdRateConstants_ddC(span dkfwd) { throw NotImplementedError("Kinetics::getFwdRateConstants_ddC", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -750,7 +752,7 @@ class Kinetics * at constant pressure, molar concentration and mole fractions. * @param[out] drop Output vector of derivatives. Length: nReactions(). */ - virtual void getFwdRatesOfProgress_ddT(double* drop) + virtual void getFwdRatesOfProgress_ddT(span drop) { throw NotImplementedError("Kinetics::getFwdRatesOfProgress_ddT", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -761,7 +763,7 @@ class Kinetics * at constant temperature, molar concentration and mole fractions. * @param[out] drop Output vector of derivatives. Length: nReactions(). */ - virtual void getFwdRatesOfProgress_ddP(double* drop) + virtual void getFwdRatesOfProgress_ddP(span drop) { throw NotImplementedError("Kinetics::getFwdRatesOfProgress_ddP", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -775,7 +777,7 @@ class Kinetics * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - virtual void getFwdRatesOfProgress_ddC(double* drop) + virtual void getFwdRatesOfProgress_ddC(span drop) { throw NotImplementedError("Kinetics::getFwdRatesOfProgress_ddC", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -823,7 +825,7 @@ class Kinetics * at constant pressure, molar concentration and mole fractions. * @param[out] drop Output vector of derivatives. Length: nReactions(). */ - virtual void getRevRatesOfProgress_ddT(double* drop) + virtual void getRevRatesOfProgress_ddT(span drop) { throw NotImplementedError("Kinetics::getRevRatesOfProgress_ddT", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -834,7 +836,7 @@ class Kinetics * at constant temperature, molar concentration and mole fractions. * @param[out] drop Output vector of derivatives. Length: nReactions(). */ - virtual void getRevRatesOfProgress_ddP(double* drop) + virtual void getRevRatesOfProgress_ddP(span drop) { throw NotImplementedError("Kinetics::getRevRatesOfProgress_ddP", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -848,7 +850,7 @@ class Kinetics * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - virtual void getRevRatesOfProgress_ddC(double* drop) + virtual void getRevRatesOfProgress_ddC(span drop) { throw NotImplementedError("Kinetics::getRevRatesOfProgress_ddC", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -896,7 +898,7 @@ class Kinetics * at constant pressure, molar concentration and mole fractions. * @param[out] drop Output vector of derivatives. Length: nReactions(). */ - virtual void getNetRatesOfProgress_ddT(double* drop) + virtual void getNetRatesOfProgress_ddT(span drop) { throw NotImplementedError("Kinetics::getNetRatesOfProgress_ddT", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -907,7 +909,7 @@ class Kinetics * at constant temperature, molar concentration and mole fractions. * @param[out] drop Output vector of derivatives. Length: nReactions(). */ - virtual void getNetRatesOfProgress_ddP(double* drop) + virtual void getNetRatesOfProgress_ddP(span drop) { throw NotImplementedError("Kinetics::getNetRatesOfProgress_ddP", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -921,7 +923,7 @@ class Kinetics * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - virtual void getNetRatesOfProgress_ddC(double* drop) + virtual void getNetRatesOfProgress_ddC(span drop) { throw NotImplementedError("Kinetics::getNetRatesOfProgress_ddC", "Not implemented for kinetics type '{}'.", kineticsType()); @@ -969,14 +971,14 @@ class Kinetics * at constant pressure, molar concentration and mole fractions. * @param[out] dwdot Output vector of derivatives. Length: #m_kk. */ - void getCreationRates_ddT(double* dwdot); + void getCreationRates_ddT(span dwdot); /** * Calculate derivatives for species creation rates with respect to pressure * at constant temperature, molar concentration and mole fractions. * @param[out] dwdot Output vector of derivatives. Length: #m_kk. */ - void getCreationRates_ddP(double* dwdot); + void getCreationRates_ddP(span dwdot); /** * Calculate derivatives for species creation rates with respect to molar @@ -986,7 +988,7 @@ class Kinetics * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - void getCreationRates_ddC(double* dwdot); + void getCreationRates_ddC(span dwdot); /** * Calculate derivatives for species creation rates with respect to species @@ -1022,14 +1024,14 @@ class Kinetics * at constant pressure, molar concentration and mole fractions. * @param[out] dwdot Output vector of derivatives. Length: #m_kk. */ - void getDestructionRates_ddT(double* dwdot); + void getDestructionRates_ddT(span dwdot); /** * Calculate derivatives for species destruction rates with respect to pressure * at constant temperature, molar concentration and mole fractions. * @param[out] dwdot Output vector of derivatives. Length: #m_kk. */ - void getDestructionRates_ddP(double* dwdot); + void getDestructionRates_ddP(span dwdot); /** * Calculate derivatives for species destruction rates with respect to molar @@ -1039,7 +1041,7 @@ class Kinetics * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - void getDestructionRates_ddC(double* dwdot); + void getDestructionRates_ddC(span dwdot); /** * Calculate derivatives for species destruction rates with respect to species @@ -1075,14 +1077,14 @@ class Kinetics * temperature at constant pressure, molar concentration and mole fractions. * @param[out] dwdot Output vector of derivatives. Length: #m_kk. */ - void getNetProductionRates_ddT(double* dwdot); + void getNetProductionRates_ddT(span dwdot); /** * Calculate derivatives for species net production rates with respect to pressure * at constant temperature, molar concentration and mole fractions. * @param[out] dwdot Output vector of derivatives. Length: #m_kk. */ - void getNetProductionRates_ddP(double* dwdot); + void getNetProductionRates_ddP(span dwdot); /** * Calculate derivatives for species net production rates with respect to molar @@ -1092,7 +1094,7 @@ class Kinetics * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - void getNetProductionRates_ddC(double* dwdot); + void getNetProductionRates_ddC(span dwdot); /** * Calculate derivatives for species net production rates with respect to species @@ -1216,7 +1218,7 @@ class Kinetics * @param kfwd Output vector containing the forward reaction rate * constants. Length: nReactions(). */ - virtual void getFwdRateConstants(double* kfwd) { + virtual void getFwdRateConstants(span kfwd) { throw NotImplementedError("Kinetics::getFwdRateConstants"); } @@ -1236,8 +1238,7 @@ class Kinetics * @param doIrreversible boolean indicating whether irreversible reactions * should be included. */ - virtual void getRevRateConstants(double* krev, - bool doIrreversible = false) { + virtual void getRevRateConstants(span krev, bool doIrreversible = false) { throw NotImplementedError("Kinetics::getRevRateConstants"); } diff --git a/include/cantera/kinetics/MultiRate.h b/include/cantera/kinetics/MultiRate.h index e4d5f1e9bc6..dd2e916f49e 100644 --- a/include/cantera/kinetics/MultiRate.h +++ b/include/cantera/kinetics/MultiRate.h @@ -66,13 +66,13 @@ class MultiRate final : public MultiRateBase m_shared.invalidateCache(); } - void getRateConstants(double* kf) override { + void getRateConstants(span kf) override { for (auto& [iRxn, rate] : m_rxn_rates) { kf[iRxn] = rate.evalFromStruct(m_shared); } } - void modifyRateConstants(double* kf, double* kr) override { + void modifyRateConstants(span kf, span kr) override { if constexpr (has_modifyRateConstants::value) { for (auto& [iRxn, rate] : m_rxn_rates) { rate.modifyRateConstants(m_shared, kf[iRxn], kr[iRxn]); @@ -80,7 +80,8 @@ class MultiRate final : public MultiRateBase } } - void processRateConstants_ddT(double* rop, const double* kf, double deltaT) override + void processRateConstants_ddT(span rop, span kf, + double deltaT) override { if constexpr (has_ddT::value) { for (const auto& [iRxn, rate] : m_rxn_rates) { @@ -106,7 +107,8 @@ class MultiRate final : public MultiRateBase } } - void processRateConstants_ddP(double* rop, const double* kf, double deltaP) override + void processRateConstants_ddP(span rop, span kf, + double deltaP) override { if constexpr (has_ddP::value) { double dPinv = 1. / (m_shared.pressure * deltaP); @@ -130,8 +132,8 @@ class MultiRate final : public MultiRateBase } } - void processRateConstants_ddM(double* rop, const double* kf, double deltaM, - bool overwrite=true) override + void processRateConstants_ddM(span rop, span kf, + double deltaM, bool overwrite=true) override { if constexpr (has_ddM::value) { double dMinv = 1. / deltaM; @@ -172,7 +174,7 @@ class MultiRate final : public MultiRateBase _update(); } - void update(double T, const vector& extra) override { + void update(double T, span extra) override { m_shared.update(T, extra); _update(); } diff --git a/include/cantera/kinetics/MultiRateBase.h b/include/cantera/kinetics/MultiRateBase.h index c52d1361428..fb8c2239b76 100644 --- a/include/cantera/kinetics/MultiRateBase.h +++ b/include/cantera/kinetics/MultiRateBase.h @@ -49,7 +49,7 @@ class MultiRateBase //! Evaluate all rate constants handled by the evaluator //! @param kf array of rate constants - virtual void getRateConstants(double* kf) = 0; + virtual void getRateConstants(span kf) = 0; //! For certain reaction types that do not follow mass action kinetics (for example, //! Butler-Volmer), calculate modifications to the forward and reverse rate @@ -62,7 +62,7 @@ class MultiRateBase //! is updated by the product StoichManagerN to determine the reverse rates of //! progress. //! @since New in %Cantera 3.2 - virtual void modifyRateConstants(double* kf, double* kr) = 0; + virtual void modifyRateConstants(span kf, span kr) = 0; //! Evaluate all rate constant temperature derivatives handled by the evaluator; //! which are multiplied with the array of rate-of-progress variables. @@ -72,8 +72,7 @@ class MultiRateBase //! contains rop on input, and d(rop)/dT on output //! @param kf array of forward rate constants (numerical derivative only) //! @param deltaT relative temperature perturbation (numerical derivative only) - virtual void processRateConstants_ddT(double* rop, - const double* kf, + virtual void processRateConstants_ddT(span rop, span kf, double deltaT) = 0; //! Evaluate all rate constant pressure derivatives handled by the evaluator; @@ -82,8 +81,7 @@ class MultiRateBase //! contains rop on input, and d(rop)/dP on output //! @param kf array of forward rate constants //! @param deltaP relative pressure perturbation - virtual void processRateConstants_ddP(double* rop, - const double* kf, + virtual void processRateConstants_ddP(span rop, span kf, double deltaP) = 0; //! Evaluate all rate constant third-body derivatives handled by the evaluator; @@ -93,10 +91,8 @@ class MultiRateBase //! @param kf array of forward rate constants //! @param deltaM relative perturbation of third-body concentrations //! @param overwrite if `true`, rop entries not affected by M are set to zero - virtual void processRateConstants_ddM(double* rop, - const double* kf, - double deltaM, - bool overwrite=true) = 0; + virtual void processRateConstants_ddM(span rop, span kf, + double deltaM, bool overwrite=true) = 0; //! Update common reaction rate data based on temperature. //! Only used in conjunction with evalSingle and ReactionRate::eval @@ -115,7 +111,7 @@ class MultiRateBase //! @param extra extra vector parameter (depends on parameterization) //! @warning This method is an experimental part of the %Cantera API and //! may be changed or removed without notice. - virtual void update(double T, const vector& extra) = 0; + virtual void update(double T, span extra) = 0; //! Update data common to reaction rates of a specific type. //! This update mechanism is used by Kinetics reaction rate evaluators. diff --git a/include/cantera/kinetics/ReactionData.h b/include/cantera/kinetics/ReactionData.h index ddb6bc93318..b89800f33ad 100644 --- a/include/cantera/kinetics/ReactionData.h +++ b/include/cantera/kinetics/ReactionData.h @@ -61,7 +61,7 @@ struct ReactionData * @warning This method is an experimental part of the %Cantera API and * may be changed or removed without notice. */ - virtual void update(double T, const vector& extra) { + virtual void update(double T, span extra) { throw NotImplementedError("ReactionData::update", "ReactionData type does not use extra vector argument."); } diff --git a/include/cantera/kinetics/ReactionPath.h b/include/cantera/kinetics/ReactionPath.h index d2d3df5ae0c..d62b3c50aac 100644 --- a/include/cantera/kinetics/ReactionPath.h +++ b/include/cantera/kinetics/ReactionPath.h @@ -263,20 +263,20 @@ class ReactionPathDiagram void exclude(const string& aaname) { m_exclude.push_back(aaname); } - void include(vector& names) { + void include(span names) { for (size_t i = 0; i < names.size(); i++) { m_include.push_back(names[i]); } } - void exclude(vector& names) { + void exclude(span names) { for (size_t i = 0; i < names.size(); i++) { m_exclude.push_back(names[i]); } } - vector& included() { + const span included() const { return m_include; } - vector& excluded() { + const span excluded() const { return m_exclude; } vector species(); @@ -286,7 +286,7 @@ class ReactionPathDiagram /** * @todo Add documentation. */ - void findMajorPaths(double threshold, size_t lda, double* a); + void findMajorPaths(double threshold, size_t lda, span a); //! Set name of the font used. void setFont(const string& font) { diff --git a/include/cantera/kinetics/ReactionRate.h b/include/cantera/kinetics/ReactionRate.h index 487e16d9f83..9954b87fca4 100644 --- a/include/cantera/kinetics/ReactionRate.h +++ b/include/cantera/kinetics/ReactionRate.h @@ -200,7 +200,7 @@ class ReactionRate //! Kinetics reaction rate evaluators. //! @warning This method is an experimental part of the %Cantera API and //! may be changed or removed without notice. - double eval(double T, const vector& extra) { + double eval(double T, span extra) { _evaluator().update(T, extra); return _evaluator().evalSingle(*this); } diff --git a/include/cantera/kinetics/StoichManager.h b/include/cantera/kinetics/StoichManager.h index 0f930a3d6f8..1164839569f 100644 --- a/include/cantera/kinetics/StoichManager.h +++ b/include/cantera/kinetics/StoichManager.h @@ -136,23 +136,23 @@ class C1 m_ic0(ic0) { } - void incrementSpecies(const double* R, double* S) const { + void incrementSpecies(span R, span S) const { S[m_ic0] += R[m_rxn]; } - void decrementSpecies(const double* R, double* S) const { + void decrementSpecies(span R, span S) const { S[m_ic0] -= R[m_rxn]; } - void multiply(const double* S, double* R) const { + void multiply(span S, span R) const { R[m_rxn] *= S[m_ic0]; } - void incrementReaction(const double* S, double* R) const { + void incrementReaction(span S, span R) const { R[m_rxn] += S[m_ic0]; } - void decrementReaction(const double* S, double* R) const { + void decrementReaction(span S, span R) const { R[m_rxn] -= S[m_ic0]; } @@ -161,14 +161,14 @@ class C1 m_jc0 = indices.at({m_rxn, m_ic0}); } - void derivatives(const double* S, const double* R, vector& jac) const + void derivatives(span S, span R, + span jac) const { jac[m_jc0] += R[m_rxn]; // index (m_ic0, m_rxn) } - void scale(const double* R, double* out, double factor) const - { + void scale(span R, span out, double factor) const { out[m_rxn] = R[m_rxn] * factor; } @@ -192,17 +192,17 @@ class C2 C2(size_t rxn = 0, size_t ic0 = 0, size_t ic1 = 0) : m_rxn(rxn), m_ic0(ic0), m_ic1(ic1) {} - void incrementSpecies(const double* R, double* S) const { + void incrementSpecies(span R, span S) const { S[m_ic0] += R[m_rxn]; S[m_ic1] += R[m_rxn]; } - void decrementSpecies(const double* R, double* S) const { + void decrementSpecies(span R, span S) const { S[m_ic0] -= R[m_rxn]; S[m_ic1] -= R[m_rxn]; } - void multiply(const double* S, double* R) const { + void multiply(span S, span R) const { if (S[m_ic0] < 0 && S[m_ic1] < 0) { R[m_rxn] = 0; } else { @@ -210,11 +210,11 @@ class C2 } } - void incrementReaction(const double* S, double* R) const { + void incrementReaction(span S, span R) const { R[m_rxn] += S[m_ic0] + S[m_ic1]; } - void decrementReaction(const double* S, double* R) const { + void decrementReaction(span S, span R) const { R[m_rxn] -= (S[m_ic0] + S[m_ic1]); } @@ -224,7 +224,8 @@ class C2 m_jc1 = indices.at({m_rxn, m_ic1}); } - void derivatives(const double* S, const double* R, vector& jac) const + void derivatives(span S, span R, + span jac) const { if (S[m_ic1] > 0) { jac[m_jc0] += R[m_rxn] * S[m_ic1]; // index (m_ic0, m_rxn) @@ -234,7 +235,7 @@ class C2 } } - void scale(const double* R, double* out, double factor) const + void scale(span R, span out, double factor) const { out[m_rxn] = 2 * R[m_rxn] * factor; } @@ -260,19 +261,19 @@ class C3 C3(size_t rxn = 0, size_t ic0 = 0, size_t ic1 = 0, size_t ic2 = 0) : m_rxn(rxn), m_ic0(ic0), m_ic1(ic1), m_ic2(ic2) {} - void incrementSpecies(const double* R, double* S) const { + void incrementSpecies(span R, span S) const { S[m_ic0] += R[m_rxn]; S[m_ic1] += R[m_rxn]; S[m_ic2] += R[m_rxn]; } - void decrementSpecies(const double* R, double* S) const { + void decrementSpecies(span R, span S) const { S[m_ic0] -= R[m_rxn]; S[m_ic1] -= R[m_rxn]; S[m_ic2] -= R[m_rxn]; } - void multiply(const double* S, double* R) const { + void multiply(span S, span R) const { if ((S[m_ic0] < 0 && (S[m_ic1] < 0 || S[m_ic2] < 0)) || (S[m_ic1] < 0 && S[m_ic2] < 0)) { R[m_rxn] = 0; @@ -281,11 +282,11 @@ class C3 } } - void incrementReaction(const double* S, double* R) const { + void incrementReaction(span S, span R) const { R[m_rxn] += S[m_ic0] + S[m_ic1] + S[m_ic2]; } - void decrementReaction(const double* S, double* R) const { + void decrementReaction(span S, span R) const { R[m_rxn] -= (S[m_ic0] + S[m_ic1] + S[m_ic2]); } @@ -296,7 +297,8 @@ class C3 m_jc2 = indices.at({m_rxn, m_ic2}); } - void derivatives(const double* S, const double* R, vector& jac) const + void derivatives(span S, span R, + span jac) const { if (S[m_ic1] > 0 && S[m_ic2] > 0) { jac[m_jc0] += R[m_rxn] * S[m_ic1] * S[m_ic2];; // index (m_ic0, m_rxn) @@ -309,7 +311,7 @@ class C3 } } - void scale(const double* R, double* out, double factor) const + void scale(span R, span out, double factor) const { out[m_rxn] = 3 * R[m_rxn] * factor; } @@ -334,8 +336,8 @@ class C_AnyN public: C_AnyN() = default; - C_AnyN(size_t rxn, const vector& ic, - const vector& order_, const vector& stoich_) : + C_AnyN(size_t rxn, span ic, + span order_, span stoich_) : m_n(ic.size()), m_rxn(rxn) { m_ic.resize(m_n); @@ -349,7 +351,7 @@ class C_AnyN } } - void multiply(const double* input, double* output) const { + void multiply(span input, span output) const { for (size_t n = 0; n < m_n; n++) { double order = m_order[n]; if (order != 0.0) { @@ -363,27 +365,27 @@ class C_AnyN } } - void incrementSpecies(const double* input, double* output) const { + void incrementSpecies(span input, span output) const { double x = input[m_rxn]; for (size_t n = 0; n < m_n; n++) { output[m_ic[n]] += m_stoich[n]*x; } } - void decrementSpecies(const double* input, double* output) const { + void decrementSpecies(span input, span output) const { double x = input[m_rxn]; for (size_t n = 0; n < m_n; n++) { output[m_ic[n]] -= m_stoich[n]*x; } } - void incrementReaction(const double* input, double* output) const { + void incrementReaction(span input, span output) const { for (size_t n = 0; n < m_n; n++) { output[m_rxn] += m_stoich[n]*input[m_ic[n]]; } } - void decrementReaction(const double* input, double* output) const { + void decrementReaction(span input, span output) const { for (size_t n = 0; n < m_n; n++) { output[m_rxn] -= m_stoich[n]*input[m_ic[n]]; } @@ -401,7 +403,8 @@ class C_AnyN } } - void derivatives(const double* S, const double* R, vector& jac) const + void derivatives(span S, span R, + span jac) const { for (size_t i = 0; i < m_n; i++) { // calculate derivative @@ -426,7 +429,7 @@ class C_AnyN } } - void scale(const double* R, double* out, double factor) const + void scale(span R, span out, double factor) const { out[m_rxn] = m_sum_order * R[m_rxn] * factor; } @@ -531,7 +534,7 @@ inline static void _resizeCoeffs(InputIter begin, InputIter end, Indices& ix) template inline static void _derivatives(InputIter begin, InputIter end, - const Vec1& S, const Vec2& R, Vec3& jac) + const Vec1& S, const Vec2& R, Vec3& jac) { for (; begin != end; ++begin) { begin->derivatives(S, R, jac); @@ -643,13 +646,13 @@ class StoichManagerN * This function is the same as the add() function below. However, the order * of each species in the power list expression is set to one automatically. */ - void add(size_t rxn, const vector& k) { + void add(size_t rxn, span k) { vector order(k.size(), 1.0); vector stoich(k.size(), 1.0); add(rxn, k, order, stoich); } - void add(size_t rxn, const vector& k, const vector& order) { + void add(size_t rxn, span k, span order) { vector stoich(k.size(), 1.0); add(rxn, k, order, stoich); } @@ -670,8 +673,8 @@ class StoichManagerN * @param stoich This is used to handle fractional stoichiometric * coefficients on the product side of irreversible reactions. */ - void add(size_t rxn, const vector& k, const vector& order, - const vector& stoich) { + void add(size_t rxn, span k, span order, + span stoich) { if (order.size() != k.size()) { throw CanteraError( "StoichManagerN::add()", "size of order and species arrays differ"); @@ -719,35 +722,35 @@ class StoichManagerN m_ready = false; } - void multiply(const double* input, double* output) const { + void multiply(span input, span output) const { _multiply(m_c1_list.begin(), m_c1_list.end(), input, output); _multiply(m_c2_list.begin(), m_c2_list.end(), input, output); _multiply(m_c3_list.begin(), m_c3_list.end(), input, output); _multiply(m_cn_list.begin(), m_cn_list.end(), input, output); } - void incrementSpecies(const double* input, double* output) const { + void incrementSpecies(span input, span output) const { _incrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); _incrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); _incrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); _incrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); } - void decrementSpecies(const double* input, double* output) const { + void decrementSpecies(span input, span output) const { _decrementSpecies(m_c1_list.begin(), m_c1_list.end(), input, output); _decrementSpecies(m_c2_list.begin(), m_c2_list.end(), input, output); _decrementSpecies(m_c3_list.begin(), m_c3_list.end(), input, output); _decrementSpecies(m_cn_list.begin(), m_cn_list.end(), input, output); } - void incrementReactions(const double* input, double* output) const { + void incrementReactions(span input, span output) const { _incrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); _incrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); _incrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); _incrementReactions(m_cn_list.begin(), m_cn_list.end(), input, output); } - void decrementReactions(const double* input, double* output) const { + void decrementReactions(span input, span output) const { _decrementReactions(m_c1_list.begin(), m_c1_list.end(), input, output); _decrementReactions(m_c2_list.begin(), m_c2_list.end(), input, output); _decrementReactions(m_c3_list.begin(), m_c3_list.end(), input, output); @@ -777,7 +780,8 @@ class StoichManagerN * @param conc Species concentration. * @param rates Rates-of-progress. */ - Eigen::SparseMatrix derivatives(const double* conc, const double* rates) + Eigen::SparseMatrix derivatives(span conc, + span rates) { // calculate derivative entries using known sparse storage order std::fill(m_values.begin(), m_values.end(), 0.); @@ -792,7 +796,7 @@ class StoichManagerN } //! Scale input by reaction order and factor - void scale(const double* in, double* out, double factor) const + void scale(span in, span out, double factor) const { _scale(m_c1_list.begin(), m_c1_list.end(), in, out, factor); _scale(m_c2_list.begin(), m_c2_list.end(), in, out, factor); diff --git a/include/cantera/kinetics/ThirdBodyCalc.h b/include/cantera/kinetics/ThirdBodyCalc.h index c1e49904678..eb51111a12f 100644 --- a/include/cantera/kinetics/ThirdBodyCalc.h +++ b/include/cantera/kinetics/ThirdBodyCalc.h @@ -72,7 +72,7 @@ class ThirdBodyCalc } //! Update third-body concentrations in full vector - void update(const vector& conc, double ctot, double* concm) const { + void update(span conc, double ctot, span concm) const { for (size_t i = 0; i < m_reaction_index.size(); i++) { double sum = 0.0; for (size_t j = 0; j < m_species[i].size(); j++) { @@ -83,7 +83,7 @@ class ThirdBodyCalc } //! Multiply output with effective third-body concentration - void multiply(double* output, const double* concm) { + void multiply(span output, span concm) { for (size_t i = 0; i < m_mass_action_index.size(); i++) { size_t ix = m_reaction_index[m_mass_action_index[i]]; output[ix] *= concm[ix]; @@ -94,13 +94,13 @@ class ThirdBodyCalc /*! * @param product Product of law of mass action and rate terms. */ - Eigen::SparseMatrix derivatives(const double* product) { - Eigen::Map mapped(product, m_multipliers.rows()); + Eigen::SparseMatrix derivatives(span product) { + Eigen::Map mapped(product.data(), m_multipliers.rows()); return mapped.asDiagonal() * m_multipliers; } //! Scale entries involving third-body collider in law of mass action by factor - void scale(const double* in, double* out, double factor) const { + void scale(span in, span out, double factor) const { for (size_t i = 0; i < m_mass_action_index.size(); i++) { size_t ix = m_reaction_index[m_mass_action_index[i]]; out[ix] = factor * in[ix]; @@ -109,8 +109,8 @@ class ThirdBodyCalc //! Scale entries involving third-body collider in rate expression //! by third-body concentration and factor - void scaleM(const double* in, double* out, - const double* concm, double factor) const + void scaleM(span in, span out, + span concm, double factor) const { for (size_t i = 0; i < m_no_mass_action_index.size(); i++) { size_t ix = m_reaction_index[m_no_mass_action_index[i]]; diff --git a/include/cantera/numerics/eigen_dense.h b/include/cantera/numerics/eigen_dense.h index 4286b3a02f8..42175f5d0ee 100644 --- a/include/cantera/numerics/eigen_dense.h +++ b/include/cantera/numerics/eigen_dense.h @@ -15,6 +15,7 @@ #else #include "cantera/ext/Eigen/Dense" #endif +#include namespace Cantera { @@ -31,10 +32,28 @@ typedef Eigen::Map ConstMappedRowVector; //! @} -//! Convenience wrapper for accessing Eigen VectorXd data as a span +template +concept EigenDenseDouble = std::same_as; + +// Require a 1D type with contiguous storage +template +concept EigenVectorLike = EigenDenseDouble + && (Derived::IsVectorAtCompileTime == 1); + +//! Convenience wrapper for accessing Eigen vector/array/map data as a span //! @todo Remove once %Cantera requires Eigen 5.0 or newer. -inline span asSpan(Eigen::VectorXd& v) { - return span(v.data(), v.size()); +template +inline span asSpan(Eigen::DenseBase& v) +{ + return span(v.derived().data(), v.size()); +} + +//! Convenience wrapper for accessing Eigen vector/array/map data as a span +//! @todo Remove once %Cantera requires Eigen 5.0 or newer. +template +inline span asSpan(const Eigen::DenseBase& v) +{ + return span(v.derived().data(), v.size()); } } diff --git a/include/cantera/numerics/funcs.h b/include/cantera/numerics/funcs.h index 217f2779788..a5c88928703 100644 --- a/include/cantera/numerics/funcs.h +++ b/include/cantera/numerics/funcs.h @@ -35,7 +35,7 @@ namespace Cantera * @returns the value of of the interpolated function at x. * @ingroup mathUtils */ -double linearInterp(double x, const vector& xpts, const vector& fpts); +double linearInterp(double x, span xpts, span fpts); //! Numerical integration of a function using the trapezoidal rule. /*! diff --git a/include/cantera/oneD/Flow1D.h b/include/cantera/oneD/Flow1D.h index 3eb7f20452c..eacba7e1517 100644 --- a/include/cantera/oneD/Flow1D.h +++ b/include/cantera/oneD/Flow1D.h @@ -467,7 +467,7 @@ class Flow1D : public Domain1D m_wtm[j] = m_thermo->meanMolecularWeight(); m_cp[j] = m_thermo->cp_mass(); m_thermo->getPartialMolarEnthalpies(span(&m_hk(0, j), m_nsp)); - m_kin->getNetProductionRates(&m_wdot(0, j)); + m_kin->getNetProductionRates(span(&m_wdot(0, j), m_nsp)); } } diff --git a/interfaces/cython/cantera/kinetics.pxd b/interfaces/cython/cantera/kinetics.pxd index 300d7455218..924adbfee4c 100644 --- a/interfaces/cython/cantera/kinetics.pxd +++ b/interfaces/cython/cantera/kinetics.pxd @@ -84,59 +84,59 @@ cdef extern from "cantera/cython/kinetics_utils.h": cdef void CxxSparseCscData "sparseCscData" (CxxSparseMatrix, double*, int*, int*) except +translate_exception # Kinetics per-reaction properties - cdef void kin_getFwdRatesOfProgress(CxxKinetics*, double*) except +translate_exception - cdef void kin_getRevRatesOfProgress(CxxKinetics*, double*) except +translate_exception - cdef void kin_getNetRatesOfProgress(CxxKinetics*, double*) except +translate_exception + cdef void kin_getFwdRatesOfProgress(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getRevRatesOfProgress(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getNetRatesOfProgress(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getFwdRateConstants_ddT(CxxKinetics*, double*) except +translate_exception - cdef void kin_getFwdRateConstants_ddP(CxxKinetics*, double*) except +translate_exception - cdef void kin_getFwdRateConstants_ddC(CxxKinetics*, double*) except +translate_exception + cdef void kin_getFwdRateConstants_ddT(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getFwdRateConstants_ddP(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getFwdRateConstants_ddC(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getFwdRatesOfProgress_ddT(CxxKinetics*, double*) except +translate_exception - cdef void kin_getRevRatesOfProgress_ddT(CxxKinetics*, double*) except +translate_exception - cdef void kin_getNetRatesOfProgress_ddT(CxxKinetics*, double*) except +translate_exception + cdef void kin_getFwdRatesOfProgress_ddT(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getRevRatesOfProgress_ddT(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getNetRatesOfProgress_ddT(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getFwdRatesOfProgress_ddP(CxxKinetics*, double*) except +translate_exception - cdef void kin_getRevRatesOfProgress_ddP(CxxKinetics*, double*) except +translate_exception - cdef void kin_getNetRatesOfProgress_ddP(CxxKinetics*, double*) except +translate_exception + cdef void kin_getFwdRatesOfProgress_ddP(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getRevRatesOfProgress_ddP(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getNetRatesOfProgress_ddP(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getFwdRatesOfProgress_ddC(CxxKinetics*, double*) except +translate_exception - cdef void kin_getRevRatesOfProgress_ddC(CxxKinetics*, double*) except +translate_exception - cdef void kin_getNetRatesOfProgress_ddC(CxxKinetics*, double*) except +translate_exception + cdef void kin_getFwdRatesOfProgress_ddC(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getRevRatesOfProgress_ddC(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getNetRatesOfProgress_ddC(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getEquilibriumConstants(CxxKinetics*, double*) except +translate_exception - cdef void kin_getFwdRateConstants(CxxKinetics*, double*) except +translate_exception - cdef void kin_getRevRateConstants(CxxKinetics*, double*) except +translate_exception + cdef void kin_getEquilibriumConstants(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getFwdRateConstants(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getRevRateConstants(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getDeltaEnthalpy(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDeltaGibbs(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDeltaEntropy(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDeltaSSEnthalpy(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDeltaSSGibbs(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDeltaSSEntropy(CxxKinetics*, double*) except +translate_exception + cdef void kin_getDeltaEnthalpy(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDeltaGibbs(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDeltaEntropy(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDeltaSSEnthalpy(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDeltaSSGibbs(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDeltaSSEntropy(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getThirdBodyConcentrations(CxxKinetics*, double*) except +translate_exception + cdef void kin_getThirdBodyConcentrations(CxxKinetics*, span[double]) except +translate_exception # Kinetics per-species properties - cdef void kin_getCreationRates(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDestructionRates(CxxKinetics*, double*) except +translate_exception - cdef void kin_getNetProductionRates(CxxKinetics*, double*) except +translate_exception + cdef void kin_getCreationRates(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDestructionRates(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getNetProductionRates(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getCreationRates_ddT(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDestructionRates_ddT(CxxKinetics*, double*) except +translate_exception - cdef void kin_getNetProductionRates_ddT(CxxKinetics*, double*) except +translate_exception + cdef void kin_getCreationRates_ddT(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDestructionRates_ddT(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getNetProductionRates_ddT(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getCreationRates_ddP(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDestructionRates_ddP(CxxKinetics*, double*) except +translate_exception - cdef void kin_getNetProductionRates_ddP(CxxKinetics*, double*) except +translate_exception + cdef void kin_getCreationRates_ddP(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDestructionRates_ddP(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getNetProductionRates_ddP(CxxKinetics*, span[double]) except +translate_exception - cdef void kin_getCreationRates_ddC(CxxKinetics*, double*) except +translate_exception - cdef void kin_getDestructionRates_ddC(CxxKinetics*, double*) except +translate_exception - cdef void kin_getNetProductionRates_ddC(CxxKinetics*, double*) except +translate_exception + cdef void kin_getCreationRates_ddC(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getDestructionRates_ddC(CxxKinetics*, span[double]) except +translate_exception + cdef void kin_getNetProductionRates_ddC(CxxKinetics*, span[double]) except +translate_exception -ctypedef void (*kineticsMethod1d)(CxxKinetics*, double*) except +translate_exception +ctypedef void (*kineticsMethod1d)(CxxKinetics*, span[double]) except +translate_exception ctypedef CxxSparseMatrix (*kineticsMethodSparse)(CxxKinetics*) except +translate_exception cdef class Kinetics(_SolutionBase): diff --git a/interfaces/cython/cantera/kinetics.pyx b/interfaces/cython/cantera/kinetics.pyx index c70a1ceb966..b0436b1082f 100644 --- a/interfaces/cython/cantera/kinetics.pyx +++ b/interfaces/cython/cantera/kinetics.pyx @@ -8,6 +8,7 @@ import numpy as np from .reaction cimport * from ._utils cimport * +from .ctcxx cimport span from . import _utils # NOTE: These cdef functions cannot be members of Kinetics because they would @@ -17,7 +18,8 @@ cdef np.ndarray get_species_array(Kinetics kin, kineticsMethod1d method): cdef np.ndarray[np.double_t, ndim=1] data = np.empty(kin.n_total_species) if kin.n_total_species == 0: return data - method(kin.kinetics, &data[0]) + cdef span[double] view = span[double](&data[0], data.size) + method(kin.kinetics, view) # @todo: Fix _selected_species to work with interface kinetics if kin._selected_species.size: return data[kin._selected_species] @@ -28,7 +30,8 @@ cdef np.ndarray get_reaction_array(Kinetics kin, kineticsMethod1d method): cdef np.ndarray[np.double_t, ndim=1] data = np.empty(kin.n_reactions) if kin.n_reactions == 0: return data - method(kin.kinetics, &data[0]) + cdef span[double] view = span[double](&data[0], data.size) + method(kin.kinetics, view) return data cdef np.ndarray get_dense(CxxSparseMatrix& smat): diff --git a/interfaces/cython/cantera/mixture.pxd b/interfaces/cython/cantera/mixture.pxd index 05179182293..c4e3ec55168 100644 --- a/interfaces/cython/cantera/mixture.pxd +++ b/interfaces/cython/cantera/mixture.pxd @@ -27,7 +27,7 @@ cdef extern from "cantera/equil/MultiPhase.h" namespace "Cantera": double phaseMoles(size_t) except +translate_exception void setPhaseMoles(size_t, double) except +translate_exception - void setMoles(double*) except +translate_exception + void setMoles(span[double]) except +translate_exception void setMolesByName(string) except +translate_exception double speciesMoles(size_t) except +translate_exception diff --git a/interfaces/cython/cantera/mixture.pyx b/interfaces/cython/cantera/mixture.pyx index 67dba58375c..050e694a506 100644 --- a/interfaces/cython/cantera/mixture.pyx +++ b/interfaces/cython/cantera/mixture.pyx @@ -2,6 +2,7 @@ # at https://cantera.org/license.txt for license and copyright information. import numpy as np +cimport numpy as np from .solutionbase cimport * from .ctcxx cimport span @@ -267,12 +268,9 @@ cdef class Mixture: self.mix.setMolesByName(stringify(moles)) return - if len(moles) != self.n_species: - raise ValueError('mole array must be of length n_species') - cdef np.ndarray[np.double_t, ndim=1] data = \ np.ascontiguousarray(moles, dtype=np.double) - self.mix.setMoles(&data[0]) + self.mix.setMoles(span[double](&data[0], data.size)) def element_moles(self, e): """ diff --git a/interfaces/cython/cantera/reaction.pxd b/interfaces/cython/cantera/reaction.pxd index a08986eaf40..2f5e8446c5e 100644 --- a/interfaces/cython/cantera/reaction.pxd +++ b/interfaces/cython/cantera/reaction.pxd @@ -9,6 +9,8 @@ from .func1 cimport * from .delegator cimport CxxDelegator from libcpp.map cimport multimap +# TODO: Drop after requiring Cython 3.1.0 or newer +ctypedef const double const_double cdef extern from "cantera/kinetics/ReactionRateFactory.h" namespace "Cantera": cdef shared_ptr[CxxReactionRate] CxxNewReactionRate "newReactionRate" (CxxAnyMap&) except +translate_exception @@ -57,8 +59,8 @@ cdef extern from "cantera/kinetics/TwoTempPlasmaRate.h" namespace "Cantera": cdef extern from "cantera/kinetics/ElectronCollisionPlasmaRate.h" namespace "Cantera": cdef cppclass CxxElectronCollisionPlasmaRate "Cantera::ElectronCollisionPlasmaRate" (CxxReactionRate): CxxElectronCollisionPlasmaRate(CxxAnyMap) except +translate_exception - vector[double]& energyLevels() - vector[double]& crossSections() + span[const_double] energyLevels() + span[const_double] crossSections() cdef extern from "cantera/base/Array.h" namespace "Cantera": cdef cppclass CxxArray2D "Cantera::Array2D": @@ -131,8 +133,9 @@ cdef extern from "cantera/kinetics/Falloff.h" namespace "Cantera": void setLowRate(CxxArrheniusRate&) except +translate_exception CxxArrheniusRate& highRate() void setHighRate(CxxArrheniusRate&) except +translate_exception - void getFalloffCoeffs(vector[double]&) - void setFalloffCoeffs(vector[double]&) except +translate_exception + void getFalloffCoeffs(span[double]) + void setFalloffCoeffs(span[double]) except +translate_exception + size_t nParameters() double evalF(double, double) except +translate_exception cdef cppclass CxxLindemannRate "Cantera::LindemannRate" (CxxFalloffRate): diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index c722fed73d8..343aeca2794 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -9,6 +9,7 @@ from cython.operator cimport dereference as deref from .kinetics cimport Kinetics from ._utils cimport * from ._utils import CanteraError +from .ctcxx cimport span from .units cimport * from .delegator cimport * @@ -420,7 +421,12 @@ cdef class ElectronCollisionPlasmaRate(ReactionRate): of `cross_sections`. """ def __get__(self): - return np.fromiter(self.base.energyLevels(), np.double) + cdef span[const_double] levels = self.base.energyLevels() + cdef np.ndarray[np.double_t, ndim=1] data = np.empty(levels.size()) + cdef size_t i + for i in range(levels.size()): + data[i] = levels[i] + return data property cross_sections: """ @@ -428,7 +434,12 @@ cdef class ElectronCollisionPlasmaRate(ReactionRate): level of `energy_levels`. """ def __get__(self): - return np.fromiter(self.base.crossSections(), np.double) + cdef span[const_double] cross_sections = self.base.crossSections() + cdef np.ndarray[np.double_t, ndim=1] data = np.empty(cross_sections.size()) + cdef size_t i + for i in range(cross_sections.size()): + data[i] = cross_sections[i] + return data cdef class FalloffRate(ReactionRate): @@ -512,14 +523,21 @@ cdef class FalloffRate(ReactionRate): property falloff_coeffs: """ The array of coefficients used to define this falloff function. """ def __get__(self): - cdef vector[double] cxxdata - self.falloff.getFalloffCoeffs(cxxdata) - return np.fromiter(cxxdata, np.double) + cdef size_t n = self.falloff.nParameters() + cdef np.ndarray[np.double_t, ndim=1] data = np.empty(n) + if n: + self.falloff.getFalloffCoeffs(span[double](&data[0], n)) + return data + def __set__(self, data): cdef vector[double] cxxdata for c in data: cxxdata.push_back(c) - self.falloff.setFalloffCoeffs(cxxdata) + if cxxdata.size(): + self.falloff.setFalloffCoeffs( + span[double](&cxxdata[0], cxxdata.size())) + else: + self.falloff.setFalloffCoeffs(span[double]()) property allow_negative_pre_exponential_factor: """ diff --git a/samples/cxx/custom/custom.cpp b/samples/cxx/custom/custom.cpp index fbccd10daf6..099065c1a02 100644 --- a/samples/cxx/custom/custom.cpp +++ b/samples/cxx/custom/custom.cpp @@ -90,7 +90,7 @@ class ReactorODEs : public FuncEval { double rho = m_gas->density(); double cp = m_gas->cp_mass(); m_gas->getPartialMolarEnthalpies(m_hbar); - m_kinetics->getNetProductionRates(&m_wdot[0]); + m_kinetics->getNetProductionRates(m_wdot); /* -------------------------- ENERGY EQUATION ------------------------- */ // the rate of change of the system temperature is found using the energy diff --git a/samples/cxx/demo/demo.cpp b/samples/cxx/demo/demo.cpp index a721864725b..1285ca513df 100644 --- a/samples/cxx/demo/demo.cpp +++ b/samples/cxx/demo/demo.cpp @@ -81,9 +81,9 @@ void demoprog() // since the gas has been set to an equilibrium state, the forward // and reverse rates of progress should be equal for all // reversible reactions, and the net rates should be zero. - kin->getFwdRatesOfProgress(&qf[0]); - kin->getRevRatesOfProgress(&qr[0]); - kin->getNetRatesOfProgress(&q[0]); + kin->getFwdRatesOfProgress(qf); + kin->getRevRatesOfProgress(qr); + kin->getNetRatesOfProgress(q); writelog("\n\n"); for (int i = 0; i < irxns; i++) { diff --git a/samples/f77/demo_ftnlib.cpp b/samples/f77/demo_ftnlib.cpp index 72d481bb407..5f9745ee75d 100644 --- a/samples/f77/demo_ftnlib.cpp +++ b/samples/f77/demo_ftnlib.cpp @@ -289,32 +289,32 @@ extern "C" { void getnetproductionrates_(double* wdot) { - _kin->getNetProductionRates(wdot); + _kin->getNetProductionRates(span(wdot, _kin->nTotalSpecies())); } void getcreationrates_(double* cdot) { - _kin->getCreationRates(cdot); + _kin->getCreationRates(span(cdot, _kin->nTotalSpecies())); } void getdestructionrates_(double* ddot) { - _kin->getDestructionRates(ddot); + _kin->getDestructionRates(span(ddot, _kin->nTotalSpecies())); } void getnetratesofprogress_(double* q) { - _kin->getNetRatesOfProgress(q); + _kin->getNetRatesOfProgress(span(q, _kin->nReactions())); } void getfwdratesofprogress_(double* q) { - _kin->getFwdRatesOfProgress(q); + _kin->getFwdRatesOfProgress(span(q, _kin->nReactions())); } void getrevratesofprogress_(double* q) { - _kin->getRevRatesOfProgress(q); + _kin->getRevRatesOfProgress(span(q, _kin->nReactions())); } //-------------------- transport properties -------------------- diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp index 5101c70438b..2f0ee4c0da5 100644 --- a/src/equil/BasisOptimize.cpp +++ b/src/equil/BasisOptimize.cpp @@ -16,9 +16,8 @@ int BasisOptimize_print_lvl = 0; static const double USEDBEFORE = -1; size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, - vector& orderVectorSpecies, - vector& orderVectorElements, - vector& formRxnMatrix) + span orderVectorSpecies, span orderVectorElements, + span formRxnMatrix) { // Get the total number of elements defined in the multiphase object size_t ne = mphase->nElements(); @@ -26,17 +25,9 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // Get the total number of species in the multiphase object size_t nspecies = mphase->nSpecies(); - // Perhaps, initialize the element ordering - if (orderVectorElements.size() < ne) { - orderVectorElements.resize(ne); - iota(orderVectorElements.begin(), orderVectorElements.end(), 0); - } - - // Perhaps, initialize the species ordering - if (orderVectorSpecies.size() != nspecies) { - orderVectorSpecies.resize(nspecies); - iota(orderVectorSpecies.begin(), orderVectorSpecies.end(), 0); - } + checkArraySize("BasisOptimize: orderVectorElements", orderVectorElements.size(), ne); + checkArraySize("BasisOptimize: orderVectorSpecies", orderVectorSpecies.size(), nspecies); + checkArraySize("BasisOptimize: formRxnMatrix", formRxnMatrix.size(), nspecies * ne); if (BasisOptimize_print_lvl >= 1) { writelog(" "); @@ -80,15 +71,13 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, // Create an array of mole numbers vector molNum(nspecies,0.0); - mphase->getMoles(molNum.data()); + mphase->getMoles(molNum); // Other workspace DenseMatrix sm(ne, ne); vector ss(ne, 0.0); vector sa(ne, 0.0); - if (formRxnMatrix.size() < nspecies*ne) { - formRxnMatrix.resize(nspecies*ne, 0.0); - } + fill(formRxnMatrix.begin(), formRxnMatrix.end(), 0.0); // For debugging purposes keep an unmodified copy of the array. vector molNumBase = molNum; @@ -289,10 +278,9 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase, } // basopt() -void ElemRearrange(size_t nComponents, const vector& elementAbundances, - MultiPhase* mphase, - vector& orderVectorSpecies, - vector& orderVectorElements) +void ElemRearrange(size_t nComponents, span elementAbundances, + MultiPhase* mphase, span orderVectorSpecies, + span orderVectorElements) { size_t nelements = mphase->nElements(); // Get the total number of species in the multiphase object @@ -306,22 +294,10 @@ void ElemRearrange(size_t nComponents, const vector& elementAbundances, writelog(" --- and to rearrange the element ordering once\n"); } - // Perhaps, initialize the element ordering - if (orderVectorElements.size() < nelements) { - orderVectorElements.resize(nelements); - for (size_t j = 0; j < nelements; j++) { - orderVectorElements[j] = j; - } - } - - // Perhaps, initialize the species ordering. However, this is dangerous, as - // this ordering is assumed to yield the component species for the problem - if (orderVectorSpecies.size() != nspecies) { - orderVectorSpecies.resize(nspecies); - for (size_t k = 0; k < nspecies; k++) { - orderVectorSpecies[k] = k; - } - } + checkArraySize("ElemRearrange: orderVectorElements", + orderVectorElements.size(), nelements); + checkArraySize("ElemRearrange: orderVectorSpecies", + orderVectorSpecies.size(), nspecies); // If the elementAbundances aren't input, just create a fake one based on // summing the column of the stoich matrix. This will force elements with diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index da80b216415..7d10201bcfa 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -113,8 +113,7 @@ void ChemEquil::initialize(ThermoPhase& s) } } -void ChemEquil::setToEquilState(ThermoPhase& s, - const vector& lambda_RT, double t) +void ChemEquil::setToEquilState(ThermoPhase& s, span lambda_RT, double t) { // Construct the chemical potentials by summing element potentials fill(m_mu_RT.begin(), m_mu_RT.end(), 0.0); @@ -163,7 +162,7 @@ void ChemEquil::update(const ThermoPhase& s) } } -int ChemEquil::setInitialMoles(ThermoPhase& s, vector& elMoleGoal, int loglevel) +int ChemEquil::setInitialMoles(ThermoPhase& s, span elMoleGoal, int loglevel) { MultiPhase mp; mp.addPhase(&s, 1.0); @@ -198,8 +197,8 @@ int ChemEquil::setInitialMoles(ThermoPhase& s, vector& elMoleGoal, int l return 0; } -int ChemEquil::estimateElementPotentials(ThermoPhase& s, vector& lambda_RT, - vector& elMolesGoal, int loglevel) +int ChemEquil::estimateElementPotentials(ThermoPhase& s, span lambda_RT, + span elMolesGoal, int loglevel) { vector b(m_mm, -999.0); vector mu_RT(m_kk, 0.0); @@ -216,7 +215,7 @@ int ChemEquil::estimateElementPotentials(ThermoPhase& s, vector& lambda_ mp.addPhase(&s, 1.0); mp.init(); int usedZeroedSpecies = 0; - vector formRxnMatrix; + vector formRxnMatrix(mp.nSpecies() * mp.nElements()); m_nComponents = BasisOptimize(&usedZeroedSpecies, false, &mp, m_orderVectorSpecies, m_orderVectorElements, formRxnMatrix); @@ -299,7 +298,7 @@ int ChemEquil::equilibrate(ThermoPhase& s, const char* XY, int loglevel) } int ChemEquil::equilibrate(ThermoPhase& s, const char* XYstr, - vector& elMolesGoal, int loglevel) + span elMolesGoal, int loglevel) { int fail = 0; bool tempFixed = true; @@ -671,9 +670,9 @@ int ChemEquil::equilibrate(ThermoPhase& s, const char* XYstr, } -int ChemEquil::dampStep(ThermoPhase& mix, vector& oldx, double oldf, - vector& grad, vector& step, vector& x, - double& f, vector& elmols, double xval, double yval) +int ChemEquil::dampStep(ThermoPhase& mix, span oldx, double oldf, + span grad, span step, span x, + double& f, span elmols, double xval, double yval) { // Carry out a delta damping approach on the dimensionless element // potentials. @@ -710,8 +709,8 @@ int ChemEquil::dampStep(ThermoPhase& mix, vector& oldx, double oldf, return 1; } -void ChemEquil::equilResidual(ThermoPhase& s, const vector& x, - const vector& elmFracGoal, vector& resid, +void ChemEquil::equilResidual(ThermoPhase& s, span x, + span elmFracGoal, span resid, double xval, double yval, int loglevel) { setToEquilState(s, x, exp(x[m_mm])); @@ -756,9 +755,8 @@ void ChemEquil::equilResidual(ThermoPhase& s, const vector& x, } } -void ChemEquil::equilJacobian(ThermoPhase& s, vector& x, - const vector& elmols, DenseMatrix& jac, - double xval, double yval, int loglevel) +void ChemEquil::equilJacobian(ThermoPhase& s, span x, span elmols, + DenseMatrix& jac, double xval, double yval, int loglevel) { vector& r0 = m_jwork1; vector& r1 = m_jwork2; @@ -789,10 +787,9 @@ void ChemEquil::equilJacobian(ThermoPhase& s, vector& x, m_doResPerturb = false; } -double ChemEquil::calcEmoles(ThermoPhase& s, vector& x, const double& n_t, - const vector& Xmol_i_calc, - vector& eMolesCalc, vector& n_i_calc, - double pressureConst) +double ChemEquil::calcEmoles(ThermoPhase& s, span x, const double& n_t, + span Xmol_i_calc, span eMolesCalc, + span n_i_calc, double pressureConst) { double n_t_calc = 0.0; @@ -825,8 +822,7 @@ double ChemEquil::calcEmoles(ThermoPhase& s, vector& x, const double& n_ return n_t_calc; } -int ChemEquil::estimateEP_Brinkley(ThermoPhase& s, vector& x, - vector& elMoles) +int ChemEquil::estimateEP_Brinkley(ThermoPhase& s, span x, span elMoles) { // Before we do anything, we will save the state of the solution. Then, if // things go drastically wrong, we will restore the saved state. @@ -1294,7 +1290,7 @@ int ChemEquil::estimateEP_Brinkley(ThermoPhase& s, vector& x, } -void ChemEquil::adjustEloc(ThermoPhase& s, vector& elMolesGoal) +void ChemEquil::adjustEloc(ThermoPhase& s, span elMolesGoal) { if (m_eloc == npos) { return; diff --git a/src/equil/MultiPhase.cpp b/src/equil/MultiPhase.cpp index 34d7ec2772e..bbd746abd44 100644 --- a/src/equil/MultiPhase.cpp +++ b/src/equil/MultiPhase.cpp @@ -34,9 +34,9 @@ void MultiPhase::addPhases(MultiPhase& mix) } } -void MultiPhase::addPhases(vector& phases, - const vector& phaseMoles) +void MultiPhase::addPhases(span phases, span phaseMoles) { + checkArraySize("MultiPhase::addPhases", phaseMoles.size(), phases.size()); for (size_t n = 0; n < phases.size(); n++) { addPhase(phases[n], phaseMoles[n]); } @@ -245,7 +245,8 @@ void MultiPhase::getChemPotentials(span mu) const } } -void MultiPhase::getValidChemPotentials(double not_mu, double* mu, bool standard) const +void MultiPhase::getValidChemPotentials(double not_mu, span mu, + bool standard) const { updatePhases(); // iterate over the phases @@ -253,14 +254,13 @@ void MultiPhase::getValidChemPotentials(double not_mu, double* mu, bool standard for (size_t i = 0; i < nPhases(); i++) { if (tempOK(i) || m_phase[i]->nSpecies() > 1) { if (!standard) { - m_phase[i]->getChemPotentials( - span(mu + loc, m_phase[i]->nSpecies())); + m_phase[i]->getChemPotentials(mu.subspan(loc, m_phase[i]->nSpecies())); } else { m_phase[i]->getStandardChemPotentials( - span(mu + loc, m_phase[i]->nSpecies())); + mu.subspan(loc, m_phase[i]->nSpecies())); } } else { - fill(mu + loc, mu + loc + m_phase[i]->nSpecies(), not_mu); + fill(mu.begin() + loc, mu.begin() + loc + m_phase[i]->nSpecies(), not_mu); } loc += m_phase[i]->nSpecies(); } @@ -335,13 +335,13 @@ double MultiPhase::cp() const return sum; } -void MultiPhase::setPhaseMoleFractions(const size_t n, const double* const x) +void MultiPhase::setPhaseMoleFractions(const size_t n, span x) { if (!m_init) { init(); } ThermoPhase* p = m_phase[n]; - p->setState_TPX(m_temp, m_press, span(x, p->nSpecies())); + p->setState_TPX(m_temp, m_press, x.subspan(0, p->nSpecies())); size_t istart = m_spstart[n]; for (size_t k = 0; k < p->nSpecies(); k++) { m_moleFractions[istart+k] = x[k]; @@ -355,7 +355,7 @@ void MultiPhase::setMolesByName(const Composition& xMap) for (size_t k = 0; k < kk; k++) { moles[k] = std::max(getValue(xMap, speciesName(k), 0.0), 0.0); } - setMoles(moles.data()); + setMoles(moles); } void MultiPhase::setMolesByName(const string& x) @@ -365,11 +365,11 @@ void MultiPhase::setMolesByName(const string& x) setMolesByName(xx); } -void MultiPhase::getMoles(double* molNum) const +void MultiPhase::getMoles(span molNum) const { // First copy in the mole fractions - copy(m_moleFractions.begin(), m_moleFractions.end(), molNum); - double* dtmp = molNum; + copy(m_moleFractions.begin(), m_moleFractions.end(), molNum.begin()); + double* dtmp = molNum.data(); for (size_t ip = 0; ip < nPhases(); ip++) { double phasemoles = m_moles[ip]; ThermoPhase* p = m_phase[ip]; @@ -380,8 +380,9 @@ void MultiPhase::getMoles(double* molNum) const } } -void MultiPhase::setMoles(const double* n) +void MultiPhase::setMoles(span n) { + checkArraySize("MultiPhase::setMoles", n.size(), nSpecies()); if (!m_init) { init(); } @@ -398,7 +399,7 @@ void MultiPhase::setMoles(const double* n) m_moles[ip] = phasemoles; if (nsp > 1) { if (phasemoles > 0.0) { - p->setState_TPX(m_temp, m_press, span(n + loc, nsp)); + p->setState_TPX(m_temp, m_press, n.subspan(loc, nsp)); p->getMoleFractions(span(&m_moleFractions[loc], nsp)); } else { p->getMoleFractions(span(&m_moleFractions[loc], nsp)); @@ -413,10 +414,10 @@ void MultiPhase::setMoles(const double* n) void MultiPhase::addSpeciesMoles(const int indexS, const double addedMoles) { vector tmpMoles(m_nsp, 0.0); - getMoles(tmpMoles.data()); + getMoles(tmpMoles); tmpMoles[indexS] += addedMoles; tmpMoles[indexS] = std::max(tmpMoles[indexS], 0.0); - setMoles(tmpMoles.data()); + setMoles(tmpMoles); } void MultiPhase::setState_TP(const double T, const double Pres) @@ -429,14 +430,15 @@ void MultiPhase::setState_TP(const double T, const double Pres) updatePhases(); } -void MultiPhase::setState_TPMoles(const double T, const double Pres, const double* n) +void MultiPhase::setState_TPMoles(const double T, const double Pres, + span n) { m_temp = T; m_press = Pres; setMoles(n); } -void MultiPhase::getElemAbundances(double* elemAbundances) const +void MultiPhase::getElemAbundances(span elemAbundances) const { calcElemAbundances(); for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) { @@ -784,9 +786,10 @@ double MultiPhase::nAtoms(const size_t kGlob, const size_t mGlob) const return m_atoms(mGlob, kGlob); } -void MultiPhase::getMoleFractions(double* const x) const +void MultiPhase::getMoleFractions(span x) const { - std::copy(m_moleFractions.begin(), m_moleFractions.end(), x); + checkArraySize("MultiPhase::getMoleFractions", x.size(), nSpecies()); + std::copy(m_moleFractions.begin(), m_moleFractions.end(), x.begin()); } string MultiPhase::phaseName(size_t iph) const diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp index b40aec09a40..7821cc4899e 100644 --- a/src/equil/MultiPhaseEquil.cpp +++ b/src/equil/MultiPhaseEquil.cpp @@ -177,7 +177,7 @@ void MultiPhaseEquil::updateMixMoles() for (size_t k = 0; k < m_nsp; k++) { m_work3[m_species[k]] = m_moles[k]; } - m_mix->setMoles(m_work3.data()); + m_mix->setMoles(m_work3); } void MultiPhaseEquil::finish() @@ -186,13 +186,13 @@ void MultiPhaseEquil::finish() for (size_t k = 0; k < m_nsp; k++) { m_work3[m_species[k]] = (m_moles[k] > 0.0 ? m_moles[k] : 0.0); } - m_mix->setMoles(m_work3.data()); + m_mix->setMoles(m_work3); } int MultiPhaseEquil::setInitialMoles(int loglevel) { double not_mu = 1.0e12; - m_mix->getValidChemPotentials(not_mu, m_mu.data(), true); + m_mix->getValidChemPotentials(not_mu, m_mu, true); double dxi_min = 1.0e10; bool redo = true; int iter = 0; @@ -246,7 +246,7 @@ int MultiPhaseEquil::setInitialMoles(int loglevel) return 0; } -void MultiPhaseEquil::getComponents(const vector& order) +void MultiPhaseEquil::getComponents(span order) { // if the input species array has the wrong size, ignore it // and consider the species for components in declaration order. @@ -389,19 +389,24 @@ void MultiPhaseEquil::getComponents(const vector& order) } } -void MultiPhaseEquil::unsort(vector& x) +void MultiPhaseEquil::unsort(span x) { - m_work2 = x; + checkArraySize("MultiPhaseEquil::unsort", x.size(), m_nsp); + m_work2.assign(x.begin(), x.end()); for (size_t k = 0; k < m_nsp; k++) { x[m_order[k]] = m_work2[k]; } } -void MultiPhaseEquil::step(double omega, vector& deltaN, int loglevel) +void MultiPhaseEquil::step(double omega, span deltaN, int loglevel) { if (omega < 0.0) { throw CanteraError("MultiPhaseEquil::step","negative omega"); } + if (deltaN.size() != m_nsp) { + throw CanteraError("MultiPhaseEquil::step", + "Expected deltaN size {}, got {}", m_nsp, deltaN.size()); + } for (size_t ik = 0; ik < m_nel; ik++) { size_t k = m_order[ik]; @@ -488,7 +493,7 @@ double MultiPhaseEquil::stepComposition(int loglevel) // If it is positive, then we have overshot the minimum. In this case, // interpolate back. double not_mu = 1.0e12; - m_mix->getValidChemPotentials(not_mu, m_mu.data()); + m_mix->getValidChemPotentials(not_mu, m_mu); double grad1 = 0.0; for (size_t k = 0; k < m_nsp; k++) { grad1 += m_work[k] * m_mu[m_species[k]]; @@ -505,14 +510,14 @@ double MultiPhaseEquil::stepComposition(int loglevel) return omega; } -double MultiPhaseEquil::computeReactionSteps(vector& dxi) +double MultiPhaseEquil::computeReactionSteps(span dxi) { - vector nu; + checkArraySize("MultiPhaseEquil::computeReactionSteps", dxi.size(), nFree()); + vector nu(m_nsp, 0.0); double grad = 0.0; - dxi.resize(nFree()); computeN(); double not_mu = 1.0e12; - m_mix->getValidChemPotentials(not_mu, m_mu.data()); + m_mix->getValidChemPotentials(not_mu, m_mu); for (size_t j = 0; j < nFree(); j++) { // get stoichiometric vector diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index c7f1792f770..7156bd68c1f 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -145,8 +145,9 @@ double vcs_VolPhase::GStar_calc_one(size_t kspec) const return StarChemicalPotential[kspec]; } -void vcs_VolPhase::setMoleFractions(const double* const xmol) +void vcs_VolPhase::setMoleFractions(span xmol) { + checkArraySize("vcs_VolPhase::setMoleFractions", xmol.size(), m_numSpecies); double sum = -1.0; for (size_t k = 0; k < m_numSpecies; k++) { Xmol_[k] = xmol[k]; @@ -175,7 +176,7 @@ void vcs_VolPhase::_updateMoleFractionDependencies() } } -const vector & vcs_VolPhase::moleFractions() const +span vcs_VolPhase::moleFractions() const { return Xmol_; } @@ -186,9 +187,10 @@ double vcs_VolPhase::moleFraction(size_t k) const } void vcs_VolPhase::setMoleFractionsState(const double totalMoles, - const double* const moleFractions, - const int vcsStateStatus) + span moleFractions, const int vcsStateStatus) { + checkArraySize("vcs_VolPhase::setMoleFractionsState", moleFractions.size(), + m_numSpecies); if (totalMoles != 0.0) { // There are other ways to set the mole fractions when VCS_STATECALC // is set to a normal settting. @@ -228,15 +230,15 @@ void vcs_VolPhase::setMoleFractionsState(const double totalMoles, } void vcs_VolPhase::setMolesFromVCS(const int stateCalc, - const double* molesSpeciesVCS) + span molesSpeciesVCS) { v_totalMoles = 0.0; - if (molesSpeciesVCS == 0) { + if (molesSpeciesVCS.empty()) { if (stateCalc == VCS_STATECALC_OLD) { - molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_old[0]; + molesSpeciesVCS = span(m_owningSolverObject->m_molNumSpecies_old); } else if (stateCalc == VCS_STATECALC_NEW) { - molesSpeciesVCS = &m_owningSolverObject->m_molNumSpecies_new[0]; + molesSpeciesVCS = span(m_owningSolverObject->m_molNumSpecies_new); } else { throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here"); } @@ -294,8 +296,8 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc, } void vcs_VolPhase::setMolesFromVCSCheck(const int vcsStateStatus, - const double* molesSpeciesVCS, - const double* const TPhMoles) + span molesSpeciesVCS, + span TPhMoles) { setMolesFromVCS(vcsStateStatus, molesSpeciesVCS); @@ -319,8 +321,7 @@ void vcs_VolPhase::updateFromVCS_MoleNumbers(const int vcsStateStatus) } } -void vcs_VolPhase::sendToVCS_ActCoeff(const int vcsStateStatus, - double* const AC) +void vcs_VolPhase::sendToVCS_ActCoeff(const int vcsStateStatus, span AC) { updateFromVCS_MoleNumbers(vcsStateStatus); if (!m_UpToDate_AC) { @@ -332,7 +333,7 @@ void vcs_VolPhase::sendToVCS_ActCoeff(const int vcsStateStatus, } } -double vcs_VolPhase::sendToVCS_VolPM(double* const VolPM) const +double vcs_VolPhase::sendToVCS_VolPM(span VolPM) const { if (!m_UpToDate_VolPM) { _updateVolPM(); @@ -344,7 +345,7 @@ double vcs_VolPhase::sendToVCS_VolPM(double* const VolPM) const return m_totalVol; } -void vcs_VolPhase::sendToVCS_GStar(double* const gstar) const +void vcs_VolPhase::sendToVCS_GStar(span gstar) const { if (!m_UpToDate_GStar) { _updateGStar(); @@ -466,7 +467,7 @@ void vcs_VolPhase::_updateLnActCoeffJac() // Go get base values for the activity coefficients. Note this calls // setState_TPX() again; Just wanted to make sure that cantera is in sync // with VolPhase after this call. - setMoleFractions(&Xmol_Base[0]); + setMoleFractions(Xmol_Base); _updateMoleFractionDependencies(); _updateActCoeff(); } @@ -497,19 +498,23 @@ double vcs_VolPhase::molefraction(size_t k) const return Xmol_[k]; } -void vcs_VolPhase::setCreationMoleNumbers(const double* const n_k, - const vector &creationGlobalRxnNumbers) +void vcs_VolPhase::setCreationMoleNumbers(span n_k, + span creationGlobalRxnNumbers) { - creationMoleNumbers_.assign(n_k, n_k+m_numSpecies); + checkArraySize("vcs_VolPhase::setCreationMoleNumbers", n_k.size(), m_numSpecies); + checkArraySize("vcs_VolPhase::setCreationMoleNumbers", creationGlobalRxnNumbers.size(), m_numSpecies); + creationMoleNumbers_.assign(n_k.begin(), n_k.end()); for (size_t k = 0; k < m_numSpecies; k++) { creationGlobalRxnNumbers_[k] = creationGlobalRxnNumbers[k]; } } -const vector& vcs_VolPhase::creationMoleNumbers( - vector &creationGlobalRxnNumbers) const +span vcs_VolPhase::creationMoleNumbers( + span creationGlobalRxnNumbers) const { - creationGlobalRxnNumbers = creationGlobalRxnNumbers_; + checkArraySize("vcs_VolPhase::creationMoleNumbers", creationGlobalRxnNumbers.size(), m_numSpecies); + copy(creationGlobalRxnNumbers_.begin(), creationGlobalRxnNumbers_.end(), + creationGlobalRxnNumbers.begin()); return creationMoleNumbers_; } diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 8ec3b0155f6..357d5c8d063 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -177,7 +177,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : kT++; } - VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, &m_molNumSpecies_old[0]); + VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, m_molNumSpecies_old); } // Work arrays used by vcs_basopt @@ -511,7 +511,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) size_t kspec = Vphase->spGlobalIndexVCS(0); // Identify the formation reaction for that species size_t irxn = kspec - m_numComponents; - vector creationGlobalRxnNumbers; + vector creationGlobalRxnNumbers(Vphase->nSpecies(), npos); // Calculate the initial moles of the phase being born. // Here we set it to 10x of the value which would cause the phase to be @@ -569,7 +569,8 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) } else { vector fracDelta(Vphase->nSpecies()); vector X_est(Vphase->nSpecies()); - fracDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers); + auto storedDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers); + copy(storedDelta.begin(), storedDelta.end(), fracDelta.begin()); double sumFrac = 0.0; for (size_t k = 0; k < Vphase->nSpecies(); k++) { @@ -671,7 +672,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) // We update the matrix dlnActCoeffdmolNumber[][] at the top of the loop, // when necessary if (m_useActCoeffJac) { - vcs_CalcLnActCoeffJac(&m_molNumSpecies_old[0]); + vcs_CalcLnActCoeffJac(m_molNumSpecies_old); } // LOOP OVER THE FORMATION REACTIONS @@ -913,11 +914,13 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) vector feSpecies_Deficient = m_feSpecies_old; // get the activity coefficients - Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, &m_actCoeffSpecies_new[0]); + Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, m_actCoeffSpecies_new); // Get the stored estimate for the composition of the phase if // it gets created - vector fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers); + vector fracDelta_new(nsp, 0.0); + const auto& storedDelta = Vphase->creationMoleNumbers(creationGlobalRxnNumbers); + copy(storedDelta.begin(), storedDelta.end(), fracDelta_new.begin()); vector componentList; for (size_t k = 0; k < nsp; k++) { @@ -990,10 +993,10 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) // Feed the newly formed estimate of the mole fractions back into the // ThermoPhase object - Vphase->setMoleFractionsState(0.0, &X_est[0], VCS_STATECALC_PHASESTABILITY); + Vphase->setMoleFractionsState(0.0, X_est, VCS_STATECALC_PHASESTABILITY); // get the activity coefficients - Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, &m_actCoeffSpecies_new[0]); + Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, m_actCoeffSpecies_new); // First calculate altered chemical potentials for component species // belonging to this phase. @@ -1127,11 +1130,11 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) if (converged) { // Save the final optimized stated back into the VolPhase object for later use - Vphase->setMoleFractionsState(0.0, &X_est[0], VCS_STATECALC_PHASESTABILITY); + Vphase->setMoleFractionsState(0.0, X_est, VCS_STATECALC_PHASESTABILITY); // Save fracDelta for later use to initialize the problem better // @todo creationGlobalRxnNumbers needs to be calculated here and stored. - Vphase->setCreationMoleNumbers(&fracDelta_new[0], creationGlobalRxnNumbers); + Vphase->setCreationMoleNumbers(fracDelta_new, creationGlobalRxnNumbers); } } else { throw CanteraError("VCS_SOLVE::vcs_phaseStabilityTest", "not done yet"); @@ -1154,7 +1157,7 @@ int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) for (size_t iph = 0; iph < m_numPhases; iph++) { vcs_VolPhase* vph = m_VolPhaseList[iph].get(); vph->setState_TP(m_temperature, m_pressurePA); - vph->sendToVCS_GStar(&m_SSfeSpecies[0]); + vph->sendToVCS_GStar(m_SSfeSpecies); } for (size_t k = 0; k < m_nsp; k++) { m_SSfeSpecies[k] /= GasConstant * m_temperature; @@ -1472,7 +1475,7 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn) return s; } -void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS) +void VCS_SOLVE::vcs_CalcLnActCoeffJac(span moleSpeciesVCS) { for (auto& Vphase : m_VolPhaseList) { // We don't need to call single species phases; @@ -1519,8 +1522,8 @@ void VCS_SOLVE::vcs_report(int iconv) for (size_t iphase = 0; iphase < m_numPhases; iphase++) { vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get(); Vphase->setState_TP(m_temperature, m_pressurePA); - Vphase->setMolesFromVCS(VCS_STATECALC_OLD, m_molNumSpecies_old.data()); - double Volp = Vphase->sendToVCS_VolPM(m_PMVolumeSpecies.data()); + Vphase->setMolesFromVCS(VCS_STATECALC_OLD, m_molNumSpecies_old); + double Volp = Vphase->sendToVCS_VolPM(m_PMVolumeSpecies); totalVolume += Volp; } @@ -1681,8 +1684,7 @@ void VCS_SOLVE::vcs_report(int iconv) // GLOBAL SATISFACTION INFORMATION // Calculate the total dimensionless Gibbs Free Energy. - double g = vcs_Total_Gibbs(&m_molNumSpecies_old[0], &m_feSpecies_old[0], - &m_tPhaseMoles_old[0]); + double g = vcs_Total_Gibbs(m_molNumSpecies_old, m_feSpecies_old, m_tPhaseMoles_old); plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g); plogf("\nElemental Abundances (kmol): "); plogf(" Actual Target Type ElActive\n"); @@ -2268,8 +2270,9 @@ int VCS_SOLVE::vcs_setMolesLinProg() return 0; } -double VCS_SOLVE::vcs_Total_Gibbs(double* molesSp, double* chemPot, - double* tPhMoles) +double VCS_SOLVE::vcs_Total_Gibbs(span molesSp, + span chemPot, + span tPhMoles) { double g = 0.0; @@ -2691,8 +2694,7 @@ void VCS_SOLVE::vcs_inest() if (m_debug_print_lvl >= 2) { plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E\n", pprefix, - vcs_Total_Gibbs(&m_molNumSpecies_old[0], &m_feSpecies_new[0], - &m_tPhaseMoles_old[0])); + vcs_Total_Gibbs(m_molNumSpecies_old, m_feSpecies_new, m_tPhaseMoles_old)); } m_VCount->T_Calls_Inest++; diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 19859599f23..bd1a4eaa60e 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -25,8 +25,8 @@ enum stages {MAIN, EQUILIB_CHECK, ELEM_ABUND_CHECK, namespace Cantera { -void VCS_SOLVE::checkDelta1(double* const dsLocal, - double* const delTPhMoles, size_t kspec) +void VCS_SOLVE::checkDelta1(span dsLocal, + span delTPhMoles, size_t kspec) { vector dchange(m_numPhases, 0.0); for (size_t k = 0; k < kspec; k++) { @@ -767,8 +767,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } } - checkDelta1(&m_deltaMolNumSpecies[0], - &m_deltaPhaseMoles[0], kspec+1); + checkDelta1(m_deltaMolNumSpecies, m_deltaPhaseMoles, kspec+1); // Branch point for returning if (m_debug_print_lvl >= 2) { @@ -838,8 +837,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } else { par = 1.0; } - checkDelta1(&m_deltaMolNumSpecies[0], - &m_deltaPhaseMoles[0], m_nsp); + checkDelta1(m_deltaMolNumSpecies, m_deltaPhaseMoles, m_nsp); // Now adjust the wt[kspec]'s so that the reflect the decrease in the // overall length of m_deltaMolNumSpecies[kspec] just calculated. At the @@ -877,11 +875,11 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, // CONVERGENCE FORCER SECTION if (printDetails) { plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n", - vcs_Total_Gibbs(&m_molNumSpecies_old[0], &m_feSpecies_old[0], - &m_tPhaseMoles_old[0])); + vcs_Total_Gibbs(m_molNumSpecies_old, m_feSpecies_old, + m_tPhaseMoles_old)); plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E\n", - vcs_Total_Gibbs(&m_molNumSpecies_new[0], &m_feSpecies_new[0], - &m_tPhaseMoles_new[0])); + vcs_Total_Gibbs(m_molNumSpecies_new, m_feSpecies_new, + m_tPhaseMoles_new)); } bool forced = vcs_globStepDamp(); @@ -917,8 +915,8 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, plogf(" Total kmoles of liquid = %15.7E\n", 0.0); } plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n", - vcs_Total_Gibbs(&m_molNumSpecies_new[0], &m_feSpecies_new[0], - &m_tPhaseMoles_new[0])); + vcs_Total_Gibbs(m_molNumSpecies_new, m_feSpecies_new, + m_tPhaseMoles_new)); plogf(" -----------------------------------------------------\n"); } @@ -975,11 +973,9 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, plogf(" "); writeline('-', 103); plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n", - vcs_Total_Gibbs(&m_molNumSpecies_old[0], &m_feSpecies_old[0], - &m_tPhaseMoles_old[0])); + vcs_Total_Gibbs(m_molNumSpecies_old, m_feSpecies_old, m_tPhaseMoles_old)); plogf(" --- Total New Dimensionless Gibbs Free Energy = %20.13E\n", - vcs_Total_Gibbs(&m_molNumSpecies_new[0], &m_feSpecies_new[0], - &m_tPhaseMoles_new[0])); + vcs_Total_Gibbs(m_molNumSpecies_new, m_feSpecies_new, m_tPhaseMoles_new)); debuglog(" --- Troublesome solve\n", m_VCount->Its > 550); } @@ -1416,43 +1412,41 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete) } } -int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr) +int VCS_SOLVE::delta_species(const size_t kspec, double& delta) { size_t irxn = kspec - m_numComponents; int retn = 1; AssertThrowMsg(kspec >= m_numComponents, "VCS_SOLVE::delta_species", "delete_species() ERROR: called for a component {}", kspec); if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - // Attempt the given dx. If it doesn't work, try to see if a smaller one + // Attempt the given delta. If it doesn't work, try to see if a smaller one // would work, - double dx = *delta_ptr; double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn); for (size_t j = 0; j < m_numComponents; ++j) { if (m_molNumSpecies_old[j] > 0.0) { - double tmp = sc_irxn[j] * dx; + double tmp = sc_irxn[j] * delta; if (-tmp > m_molNumSpecies_old[j]) { retn = 0; - dx = std::min(dx, - m_molNumSpecies_old[j] / sc_irxn[j]); + delta = std::min(delta, - m_molNumSpecies_old[j] / sc_irxn[j]); } } // If the component has a zero concentration and is a reactant // in the formation reaction, then dx == 0.0, and we just return. if (m_molNumSpecies_old[j] <= 0.0 && sc_irxn[j] < 0.0) { - *delta_ptr = 0.0; + delta = 0.0; return 0; } } - // ok, we found a positive dx. implement it. - *delta_ptr = dx; - m_molNumSpecies_old[kspec] += dx; + // ok, we found a positive delta. implement it. + m_molNumSpecies_old[kspec] += delta; size_t iph = m_phaseID[kspec]; - m_tPhaseMoles_old[iph] += dx; + m_tPhaseMoles_old[iph] += delta; m_VolPhaseList[iph]->setMolesOutOfDate(VCS_STATECALC_OLD); for (size_t j = 0; j < m_numComponents; ++j) { - double tmp = sc_irxn[j] * dx; + double tmp = sc_irxn[j] * delta; if (tmp != 0.0) { iph = m_phaseID[j]; m_molNumSpecies_old[j] += tmp; @@ -1473,7 +1467,7 @@ int VCS_SOLVE::vcs_zero_species(const size_t kspec) if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double dx = -m_molNumSpecies_old[kspec]; if (dx != 0.0) { - retn = delta_species(kspec, &dx); + retn = delta_species(kspec, dx); if (!retn && m_debug_print_lvl >= 1) { plogf("vcs_zero_species: Couldn't zero the species %d, " "did delta of %g. orig conc of %g\n", @@ -1517,8 +1511,8 @@ int VCS_SOLVE::vcs_delete_species(const size_t kspec) } // Adjust the total moles in a phase downwards. - Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, &m_molNumSpecies_old[0], - &m_tPhaseMoles_old[0]); + Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, m_molNumSpecies_old, + m_tPhaseMoles_old); // Adjust the current number of active species and reactions counters --m_numRxnRdc; @@ -1556,7 +1550,7 @@ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec) // this adjusts m_molNumSpecies_old[] and m_tPhaseMoles_old[] // HKM -> make this a relative mole number! double dx = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10.; - delta_species(kspec, &dx); + delta_species(kspec, dx); m_speciesStatus[kspec] = VCS_SPECIES_MINOR; if (m_SSPhase[kspec]) { @@ -1565,9 +1559,8 @@ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec) } vcs_VolPhase* Vphase = m_VolPhaseList[iph].get(); - Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, - &m_molNumSpecies_old[0], - &m_tPhaseMoles_old[0]); + Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, m_molNumSpecies_old, + m_tPhaseMoles_old); // We may have popped a multispecies phase back into existence. If we did, // we have to check the other species in that phase. Take care of the @@ -1615,7 +1608,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) double dx = - m_molNumSpecies_old[kspec]; double dxTent = dx; - int retn = delta_species(kspec, &dxTent); + int retn = delta_species(kspec, dxTent); if (retn != 1) { successful = false; if (m_debug_print_lvl >= 2) { @@ -1676,7 +1669,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) } } if (dxPerm != 0.0) { - delta_species(kspec, &dxPerm); + delta_species(kspec, dxPerm); } } } @@ -1824,7 +1817,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() m_molNumSpecies_new[kspec] = VCS_DELETE_MINORSPECIES_CUTOFF * 1.0E-10; } if (!Vphase->m_singleSpecies) { - Vphase->sendToVCS_ActCoeff(VCS_STATECALC_NEW, &m_actCoeffSpecies_new[0]); + Vphase->sendToVCS_ActCoeff(VCS_STATECALC_NEW, m_actCoeffSpecies_new); } m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_new[kspec]) - m_lnMnaughtSpecies[kspec] + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]); @@ -1848,7 +1841,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() size_t iph = m_phaseID[kspec]; if (m_tPhaseMoles_old[iph] > 0.0) { double dx = m_molNumSpecies_new[kspec]; - size_t retn = delta_species(kspec, &dx); + size_t retn = delta_species(kspec, dx); if (retn == 0) { if (m_debug_print_lvl) { plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", @@ -1856,7 +1849,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() } if (dx > 1.0E-50) { dx = 1.0E-50; - retn = delta_species(kspec, &dx); + retn = delta_species(kspec, dx); if (retn == 0 && m_debug_print_lvl) { plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", m_speciesName[kspec], kspec, dx); @@ -2650,20 +2643,20 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const void VCS_SOLVE::vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop) { - double* tPhMoles_ptr=0; - double* actCoeff_ptr=0; - double* feSpecies=0; - double* molNum=0; + span tPhMoles; + span actCoeff; + span feSpecies; + span molNum; if (stateCalc == VCS_STATECALC_OLD) { - feSpecies = &m_feSpecies_old[0]; - tPhMoles_ptr = &m_tPhaseMoles_old[0]; - actCoeff_ptr = &m_actCoeffSpecies_old[0]; - molNum = &m_molNumSpecies_old[0]; + feSpecies = span(m_feSpecies_old); + tPhMoles = span(m_tPhaseMoles_old); + actCoeff = span(m_actCoeffSpecies_old); + molNum = span(m_molNumSpecies_old); } else if (stateCalc == VCS_STATECALC_NEW) { - feSpecies = &m_feSpecies_new[0]; - tPhMoles_ptr = &m_tPhaseMoles_new[0]; - actCoeff_ptr = &m_actCoeffSpecies_new[0]; - molNum = &m_molNumSpecies_new[0]; + feSpecies = span(m_feSpecies_new); + tPhMoles = span(m_tPhaseMoles_new); + actCoeff = span(m_actCoeffSpecies_new); + molNum = span(m_molNumSpecies_new); } else { throw CanteraError("VCS_SOLVE::vcs_dfe", "Subroutine vcs_dfe called with bad stateCalc value: {}", stateCalc); @@ -2700,15 +2693,15 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, } } for (size_t iph = 0; iph < m_numPhases; iph++) { - AssertThrowMsg(vcs_doubleEqual(tlogMoles[iph], tPhMoles_ptr[iph]), + AssertThrowMsg(vcs_doubleEqual(tlogMoles[iph], tPhMoles[iph]), "VCS_SOLVE::vcs_dfe", "phase Moles may be off, iph = {}, {} {}", - iph, tlogMoles[iph], tPhMoles_ptr[iph]); + iph, tlogMoles[iph], tPhMoles[iph]); } m_TmpPhase.assign(m_TmpPhase.size(), 0.0); for (size_t iph = 0; iph < m_numPhases; iph++) { - if (tPhMoles_ptr[iph] > 0.0) { - tlogMoles[iph] = log(tPhMoles_ptr[iph]); + if (tPhMoles[iph] > 0.0) { + tlogMoles[iph] = log(tPhMoles[iph]); } } @@ -2719,7 +2712,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get(); Vphase->updateFromVCS_MoleNumbers(stateCalc); if (!Vphase->m_singleSpecies) { - Vphase->sendToVCS_ActCoeff(stateCalc, &actCoeff_ptr[0]); + Vphase->sendToVCS_ActCoeff(stateCalc, actCoeff); } m_phasePhi[iphase] = Vphase->electricPotential(); } @@ -2746,14 +2739,14 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, double logTerm; if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) { - if (tPhMoles_ptr[iphase] > 0.0) { - logTerm = log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF) + if (tPhMoles[iphase] > 0.0) { + logTerm = log(actCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF) - tlogMoles[iphase]; } else { return m_SSfeSpecies[kspec] - m_lnMnaughtSpecies[kspec] + phiTerm; } } else { - logTerm = log(actCoeff_ptr[kspec] * molNum[kspec]) - tlogMoles[iphase]; + logTerm = log(actCoeff[kspec] * molNum[kspec]) - tlogMoles[iphase]; } return m_SSfeSpecies[kspec] + logTerm - m_lnMnaughtSpecies[kspec] + phiTerm; }; @@ -2858,13 +2851,11 @@ void VCS_SOLVE::vcs_updateVP(const int vcsState) { for (auto& Vphase : m_VolPhaseList) { if (vcsState == VCS_STATECALC_OLD) { - Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, - &m_molNumSpecies_old[0], - &m_tPhaseMoles_old[0]); + Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, m_molNumSpecies_old, + m_tPhaseMoles_old); } else if (vcsState == VCS_STATECALC_NEW) { - Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW, - &m_molNumSpecies_new[0], - &m_tPhaseMoles_new[0]); + Vphase->setMolesFromVCSCheck(VCS_STATECALC_NEW, m_molNumSpecies_new, + m_tPhaseMoles_new); } else { throw CanteraError("VCS_SOLVE::vcs_updateVP", "wrong stateCalc value: {}", vcsState); @@ -2946,20 +2937,17 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted, irxnl = m_numRxnTot; } - double* deltaGRxn; - double* feSpecies; - double* molNumSpecies; - double* actCoeffSpecies; + span deltaGRxn, feSpecies, molNumSpecies, actCoeffSpecies; if (vcsState == VCS_STATECALC_NEW) { - deltaGRxn = &m_deltaGRxn_new[0]; - feSpecies = &m_feSpecies_new[0]; - molNumSpecies = &m_molNumSpecies_new[0]; - actCoeffSpecies = &m_actCoeffSpecies_new[0]; + deltaGRxn = span(m_deltaGRxn_new); + feSpecies = span(m_feSpecies_new); + molNumSpecies = span(m_molNumSpecies_new); + actCoeffSpecies = span(m_actCoeffSpecies_new); } else if (vcsState == VCS_STATECALC_OLD) { - deltaGRxn = &m_deltaGRxn_old[0]; - feSpecies = &m_feSpecies_old[0]; - molNumSpecies = &m_molNumSpecies_old[0]; - actCoeffSpecies = &m_actCoeffSpecies_old[0]; + deltaGRxn = span(m_deltaGRxn_old); + feSpecies = span(m_feSpecies_old); + molNumSpecies = span(m_molNumSpecies_old); + actCoeffSpecies = span(m_actCoeffSpecies_old); } else { throw CanteraError("VCS_SOLVE::vcs_deltag", "bad vcsState"); } diff --git a/src/equil/vcs_util.cpp b/src/equil/vcs_util.cpp index 6ac343fc717..c10fcd5c68e 100644 --- a/src/equil/vcs_util.cpp +++ b/src/equil/vcs_util.cpp @@ -16,7 +16,7 @@ namespace Cantera { -double vcs_l2norm(const vector& vec) +double vcs_l2norm(span vec) { if (vec.empty()) { return 0.0; diff --git a/src/fortran/fct.cpp b/src/fortran/fct.cpp index d987bc0cbad..313bec4cb71 100644 --- a/src/fortran/fct.cpp +++ b/src/fortran/fct.cpp @@ -809,7 +809,7 @@ extern "C" { { try { Kinetics* k = _fkin(n); - k->getFwdRatesOfProgress(fwdROP); + k->getFwdRatesOfProgress(span(fwdROP, k->nReactions())); } catch (...) { return handleAllExceptions(-1, ERR); } @@ -820,7 +820,7 @@ extern "C" { { try { Kinetics* k = _fkin(n); - k->getRevRatesOfProgress(revROP); + k->getRevRatesOfProgress(span(revROP, k->nReactions())); } catch (...) { return handleAllExceptions(-1, ERR); } @@ -840,7 +840,7 @@ extern "C" { { try { Kinetics* k = _fkin(n); - k->getNetRatesOfProgress(netROP); + k->getNetRatesOfProgress(span(netROP, k->nReactions())); } catch (...) { return handleAllExceptions(-1, ERR); } @@ -851,7 +851,7 @@ extern "C" { { try { Kinetics* k = _fkin(n); - k->getCreationRates(cdot); + k->getCreationRates(span(cdot, k->nTotalSpecies())); } catch (...) { return handleAllExceptions(-1, ERR); } @@ -862,7 +862,7 @@ extern "C" { { try { Kinetics* k = _fkin(n); - k->getDestructionRates(ddot); + k->getDestructionRates(span(ddot, k->nTotalSpecies())); } catch (...) { return handleAllExceptions(-1, ERR); } @@ -873,7 +873,7 @@ extern "C" { { try { Kinetics* k = _fkin(n); - k->getNetProductionRates(wdot); + k->getNetProductionRates(span(wdot, k->nTotalSpecies())); } catch (...) { return handleAllExceptions(-1, ERR); } @@ -893,7 +893,7 @@ extern "C" { { try { Kinetics* k = _fkin(n); - k->getEquilibriumConstants(kc); + k->getEquilibriumConstants(span(kc, k->nReactions())); } catch (...) { return handleAllExceptions(-1, ERR); } diff --git a/src/kinetics/BulkKinetics.cpp b/src/kinetics/BulkKinetics.cpp index cdbbe5ae707..43c9b73f11e 100644 --- a/src/kinetics/BulkKinetics.cpp +++ b/src/kinetics/BulkKinetics.cpp @@ -4,6 +4,7 @@ #include "cantera/kinetics/BulkKinetics.h" #include "cantera/kinetics/Reaction.h" #include "cantera/thermo/ThermoPhase.h" +#include "cantera/numerics/eigen_dense.h" namespace Cantera { @@ -148,24 +149,26 @@ void BulkKinetics::invalidateCache() m_ROP_ok = false; } -void BulkKinetics::getFwdRateConstants(double* kfwd) +void BulkKinetics::getFwdRateConstants(span kfwd) { + checkArraySize("BulkKinetics::getFwdRateConstants", kfwd.size(), nReactions()); updateROP(); - copy(m_rfn.begin(), m_rfn.end(), kfwd); + copy(m_rfn.begin(), m_rfn.end(), kfwd.begin()); if (legacy_rate_constants_used()) { processThirdBodies(kfwd); } } -void BulkKinetics::getEquilibriumConstants(double* kc) +void BulkKinetics::getEquilibriumConstants(span kc) { + checkArraySize("BulkKinetics::getEquilibriumConstants", kc.size(), nReactions()); updateROP(); vector& delta_gibbs0 = m_rbuf0; fill(delta_gibbs0.begin(), delta_gibbs0.end(), 0.0); // compute Delta G^0 for all reactions - getReactionDelta(m_grt.data(), delta_gibbs0.data()); + getReactionDelta(m_grt, delta_gibbs0); double rrt = 1.0 / thermo().RT(); double logStandConc = log(thermo().standardConcentration()); @@ -174,14 +177,14 @@ void BulkKinetics::getEquilibriumConstants(double* kc) } } -void BulkKinetics::getRevRateConstants(double* krev, bool doIrreversible) +void BulkKinetics::getRevRateConstants(span krev, bool doIrreversible) { // go get the forward rate constants. -> note, we don't really care about // speed or redundancy in these informational routines. getFwdRateConstants(krev); if (doIrreversible) { - getEquilibriumConstants(m_rbuf0.data()); + getEquilibriumConstants(m_rbuf0); for (size_t i = 0; i < nReactions(); i++) { krev[i] /= m_rbuf0[i]; } @@ -193,31 +196,31 @@ void BulkKinetics::getRevRateConstants(double* krev, bool doIrreversible) } } -void BulkKinetics::getDeltaGibbs(double* deltaG) +void BulkKinetics::getDeltaGibbs(span deltaG) { // Get the chemical potentials for each species thermo().getChemPotentials(m_sbuf0); // Use the stoichiometric manager to find deltaG for each reaction. - getReactionDelta(m_sbuf0.data(), deltaG); + getReactionDelta(m_sbuf0, deltaG); } -void BulkKinetics::getDeltaEnthalpy(double* deltaH) +void BulkKinetics::getDeltaEnthalpy(span deltaH) { // Get the partial molar enthalpy for each species thermo().getPartialMolarEnthalpies(m_sbuf0); // Use the stoichiometric manager to find deltaH for each reaction. - getReactionDelta(m_sbuf0.data(), deltaH); + getReactionDelta(m_sbuf0, deltaH); } -void BulkKinetics::getDeltaEntropy(double* deltaS) +void BulkKinetics::getDeltaEntropy(span deltaS) { // Get the partial molar entropy for each species thermo().getPartialMolarEntropies(m_sbuf0); // Use the stoichiometric manager to find deltaS for each reaction. - getReactionDelta(m_sbuf0.data(), deltaS); + getReactionDelta(m_sbuf0, deltaS); } -void BulkKinetics::getDeltaSSGibbs(double* deltaG) +void BulkKinetics::getDeltaSSGibbs(span deltaG) { // Get the standard state chemical potentials of the species. This is the // array of chemical potentials at unit activity. We define these here as @@ -225,10 +228,10 @@ void BulkKinetics::getDeltaSSGibbs(double* deltaG) // pressure of the solution. thermo().getStandardChemPotentials(m_sbuf0); // Use the stoichiometric manager to find deltaG for each reaction. - getReactionDelta(m_sbuf0.data(), deltaG); + getReactionDelta(m_sbuf0, deltaG); } -void BulkKinetics::getDeltaSSEnthalpy(double* deltaH) +void BulkKinetics::getDeltaSSEnthalpy(span deltaH) { // Get the standard state enthalpies of the species. thermo().getEnthalpy_RT(m_sbuf0); @@ -236,10 +239,10 @@ void BulkKinetics::getDeltaSSEnthalpy(double* deltaH) m_sbuf0[k] *= thermo().RT(); } // Use the stoichiometric manager to find deltaH for each reaction. - getReactionDelta(m_sbuf0.data(), deltaH); + getReactionDelta(m_sbuf0, deltaH); } -void BulkKinetics::getDeltaSSEntropy(double* deltaS) +void BulkKinetics::getDeltaSSEntropy(span deltaS) { // Get the standard state entropy of the species. We define these here as // the entropies of the pure species at the temperature and pressure of the @@ -249,7 +252,7 @@ void BulkKinetics::getDeltaSSEntropy(double* deltaS) m_sbuf0[k] *= GasConstant; } // Use the stoichiometric manager to find deltaS for each reaction. - getReactionDelta(m_sbuf0.data(), deltaS); + getReactionDelta(m_sbuf0, deltaS); } void BulkKinetics::getDerivativeSettings(AnyMap& settings) const @@ -273,105 +276,105 @@ void BulkKinetics::setDerivativeSettings(const AnyMap& settings) } } -void BulkKinetics::getFwdRateConstants_ddT(double* dkfwd) +void BulkKinetics::getFwdRateConstants_ddT(span dkfwd) { assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddT"); updateROP(); process_ddT(m_rfn, dkfwd); } -void BulkKinetics::getFwdRatesOfProgress_ddT(double* drop) +void BulkKinetics::getFwdRatesOfProgress_ddT(span drop) { assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddT"); updateROP(); process_ddT(m_ropf, drop); } -void BulkKinetics::getRevRatesOfProgress_ddT(double* drop) +void BulkKinetics::getRevRatesOfProgress_ddT(span drop) { assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddT"); updateROP(); process_ddT(m_ropr, drop); - Eigen::Map dRevRop(drop, nReactions()); // reverse rop times scaled inverse equilibrium constant derivatives - Eigen::Map dRevRop2(m_rbuf2.data(), nReactions()); copy(m_ropr.begin(), m_ropr.end(), m_rbuf2.begin()); - applyEquilibriumConstants_ddT(dRevRop2.data()); + applyEquilibriumConstants_ddT(m_rbuf2); + Eigen::Map dRevRop(drop.data(), nReactions()); + Eigen::Map dRevRop2(m_rbuf2.data(), nReactions()); dRevRop += dRevRop2; } -void BulkKinetics::getNetRatesOfProgress_ddT(double* drop) +void BulkKinetics::getNetRatesOfProgress_ddT(span drop) { assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddT"); updateROP(); process_ddT(m_ropnet, drop); - Eigen::Map dNetRop(drop, nReactions()); // reverse rop times scaled inverse equilibrium constant derivatives - Eigen::Map dNetRop2(m_rbuf2.data(), nReactions()); copy(m_ropr.begin(), m_ropr.end(), m_rbuf2.begin()); - applyEquilibriumConstants_ddT(dNetRop2.data()); + applyEquilibriumConstants_ddT(m_rbuf2); + Eigen::Map dNetRop(drop.data(), nReactions()); + Eigen::Map dNetRop2(m_rbuf2.data(), nReactions()); dNetRop -= dNetRop2; } -void BulkKinetics::getFwdRateConstants_ddP(double* dkfwd) +void BulkKinetics::getFwdRateConstants_ddP(span dkfwd) { assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddP"); updateROP(); process_ddP(m_rfn, dkfwd); } -void BulkKinetics::getFwdRatesOfProgress_ddP(double* drop) +void BulkKinetics::getFwdRatesOfProgress_ddP(span drop) { assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddP"); updateROP(); process_ddP(m_ropf, drop); } -void BulkKinetics::getRevRatesOfProgress_ddP(double* drop) +void BulkKinetics::getRevRatesOfProgress_ddP(span drop) { assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddP"); updateROP(); process_ddP(m_ropr, drop); } -void BulkKinetics::getNetRatesOfProgress_ddP(double* drop) +void BulkKinetics::getNetRatesOfProgress_ddP(span drop) { assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddP"); updateROP(); process_ddP(m_ropnet, drop); } -void BulkKinetics::getFwdRateConstants_ddC(double* dkfwd) +void BulkKinetics::getFwdRateConstants_ddC(span dkfwd) { assertDerivativesValid("BulkKinetics::getFwdRateConstants_ddC"); updateROP(); process_ddC(m_reactantStoich, m_rfn, dkfwd, false); } -void BulkKinetics::getFwdRatesOfProgress_ddC(double* drop) +void BulkKinetics::getFwdRatesOfProgress_ddC(span drop) { assertDerivativesValid("BulkKinetics::getFwdRatesOfProgress_ddC"); updateROP(); process_ddC(m_reactantStoich, m_ropf, drop); } -void BulkKinetics::getRevRatesOfProgress_ddC(double* drop) +void BulkKinetics::getRevRatesOfProgress_ddC(span drop) { assertDerivativesValid("BulkKinetics::getRevRatesOfProgress_ddC"); updateROP(); return process_ddC(m_revProductStoich, m_ropr, drop); } -void BulkKinetics::getNetRatesOfProgress_ddC(double* drop) +void BulkKinetics::getNetRatesOfProgress_ddC(span drop) { assertDerivativesValid("BulkKinetics::getNetRatesOfProgress_ddC"); updateROP(); process_ddC(m_reactantStoich, m_ropf, drop); - Eigen::Map dNetRop(drop, nReactions()); + Eigen::Map dNetRop(drop.data(), nReactions()); - process_ddC(m_revProductStoich, m_ropr, m_rbuf2.data()); + process_ddC(m_revProductStoich, m_ropr, m_rbuf2); Eigen::Map dNetRop2(m_rbuf2.data(), nReactions()); dNetRop -= dNetRop2; } @@ -382,7 +385,7 @@ Eigen::SparseMatrix BulkKinetics::fwdRatesOfProgress_ddX() // forward reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); return calculateCompositionDerivatives(m_reactantStoich, rop_rates); } @@ -392,8 +395,8 @@ Eigen::SparseMatrix BulkKinetics::revRatesOfProgress_ddX() // reverse reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); - applyEquilibriumConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); + applyEquilibriumConstants(rop_rates); return calculateCompositionDerivatives(m_revProductStoich, rop_rates); } @@ -403,11 +406,11 @@ Eigen::SparseMatrix BulkKinetics::netRatesOfProgress_ddX() // forward reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates); // reverse reaction rate coefficients - applyEquilibriumConstants(rop_rates.data()); + applyEquilibriumConstants(rop_rates); return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates); } @@ -417,7 +420,7 @@ Eigen::SparseMatrix BulkKinetics::fwdRatesOfProgress_ddCi() // forward reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); return calculateCompositionDerivatives(m_reactantStoich, rop_rates, false); } @@ -427,8 +430,8 @@ Eigen::SparseMatrix BulkKinetics::revRatesOfProgress_ddCi() // reverse reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); - applyEquilibriumConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); + applyEquilibriumConstants(rop_rates); return calculateCompositionDerivatives(m_revProductStoich, rop_rates, false); } @@ -438,11 +441,11 @@ Eigen::SparseMatrix BulkKinetics::netRatesOfProgress_ddCi() // forward reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); auto jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates, false); // reverse reaction rate coefficients - applyEquilibriumConstants(rop_rates.data()); + applyEquilibriumConstants(rop_rates); return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates, false); } @@ -461,7 +464,7 @@ void BulkKinetics::updateROP() double logStandConc = log(thermo().standardConcentration()); // compute Delta G^0 for all reversible reactions - getRevReactionDelta(m_grt.data(), m_delta_gibbs0.data()); + getRevReactionDelta(m_grt, m_delta_gibbs0); double rrt = 1.0 / thermo().RT(); for (size_t i = 0; i < m_revindex.size(); i++) { @@ -482,7 +485,7 @@ void BulkKinetics::updateROP() double ctot = thermo().molarDensity(); // Third-body objects interacting with MultiRate evaluator - m_multi_concm.update(m_phys_conc, ctot, m_concm.data()); + m_multi_concm.update(m_phys_conc, ctot, m_concm); m_ROP_ok = false; } @@ -490,7 +493,7 @@ void BulkKinetics::updateROP() for (auto& rates : m_rateHandlers) { bool changed = rates->update(thermo(), *this); if (changed) { - rates->getRateConstants(m_kf0.data()); + rates->getRateConstants(m_kf0); m_ROP_ok = false; } } @@ -506,20 +509,20 @@ void BulkKinetics::updateROP() m_rfn[i] = m_kf0[i] * m_perturb[i]; } - copy(m_rfn.begin(), m_rfn.end(), m_ropf.data()); - processThirdBodies(m_ropf.data()); + copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin()); + processThirdBodies(m_ropf); copy(m_ropf.begin(), m_ropf.end(), m_ropr.begin()); // for reversible reactions, multiply ropr by concentration products - applyEquilibriumConstants(m_ropr.data()); + applyEquilibriumConstants(m_ropr); for (auto& rates : m_rateHandlers) { - rates->modifyRateConstants(m_ropf.data(), m_ropr.data()); + rates->modifyRateConstants(m_ropf, m_ropr); } // multiply ropf and ropr by concentration products - m_reactantStoich.multiply(m_act_conc.data(), m_ropf.data()); - m_revProductStoich.multiply(m_act_conc.data(), m_ropr.data()); + m_reactantStoich.multiply(m_act_conc, m_ropf); + m_revProductStoich.multiply(m_act_conc, m_ropr); for (size_t j = 0; j != nReactions(); ++j) { m_ropnet[j] = m_ropf[j] - m_ropr[j]; @@ -536,21 +539,22 @@ void BulkKinetics::updateROP() m_ROP_ok = true; } -void BulkKinetics::getThirdBodyConcentrations(double* concm) +void BulkKinetics::getThirdBodyConcentrations(span concm) { + checkArraySize("BulkKinetics::getThirdBodyConcentrations", concm.size(), nReactions()); updateROP(); - std::copy(m_concm.begin(), m_concm.end(), concm); + std::copy(m_concm.begin(), m_concm.end(), concm.begin()); } -void BulkKinetics::processThirdBodies(double* rop) +void BulkKinetics::processThirdBodies(span rop) { // reactions involving third body if (!m_concm.empty()) { - m_multi_concm.multiply(rop, m_concm.data()); + m_multi_concm.multiply(rop, m_concm); } } -void BulkKinetics::applyEquilibriumConstants(double* rop) +void BulkKinetics::applyEquilibriumConstants(span rop) { // For reverse rates computed from thermochemistry, multiply the forward // rate coefficients by the reciprocals of the equilibrium constants @@ -559,7 +563,7 @@ void BulkKinetics::applyEquilibriumConstants(double* rop) } } -void BulkKinetics::applyEquilibriumConstants_ddT(double* drkcn) +void BulkKinetics::applyEquilibriumConstants_ddT(span drkcn) { double T = thermo().temperature(); double P = thermo().pressure(); @@ -573,7 +577,7 @@ void BulkKinetics::applyEquilibriumConstants_ddT(double* drkcn) thermo().saveState(m_state); thermo().setState_TP(T * (1. + m_jac_rtol_delta), P); thermo().getStandardChemPotentials(grt); - getRevReactionDelta(grt.data(), delta_gibbs0.data()); + getRevReactionDelta(grt, delta_gibbs0); // apply scaling for derivative of inverse equilibrium constant double Tinv = 1. / T; @@ -594,34 +598,34 @@ void BulkKinetics::applyEquilibriumConstants_ddT(double* drkcn) thermo().restoreState(m_state); } -void BulkKinetics::process_ddT(const vector& in, double* drop) +void BulkKinetics::process_ddT(span in, span drop) { // apply temperature derivative - copy(in.begin(), in.end(), drop); + copy(in.begin(), in.end(), drop.begin()); for (auto& rates : m_rateHandlers) { - rates->processRateConstants_ddT(drop, m_rfn.data(), m_jac_rtol_delta); + rates->processRateConstants_ddT(drop, m_rfn, m_jac_rtol_delta); } } -void BulkKinetics::process_ddP(const vector& in, double* drop) +void BulkKinetics::process_ddP(span in, span drop) { // apply pressure derivative - copy(in.begin(), in.end(), drop); + copy(in.begin(), in.end(), drop.begin()); for (auto& rates : m_rateHandlers) { - rates->processRateConstants_ddP(drop, m_rfn.data(), m_jac_rtol_delta); + rates->processRateConstants_ddP(drop, m_rfn, m_jac_rtol_delta); } } -void BulkKinetics::process_ddC(StoichManagerN& stoich, const vector& in, - double* drop, bool mass_action) +void BulkKinetics::process_ddC(StoichManagerN& stoich, span in, + span drop, bool mass_action) { - Eigen::Map out(drop, nReactions()); + Eigen::Map out(drop.data(), nReactions()); out.setZero(); double ctot_inv = 1. / thermo().molarDensity(); // derivatives due to concentrations in law of mass action if (mass_action) { - stoich.scale(in.data(), out.data(), ctot_inv); + stoich.scale(in, span(out.data(), nReactions()), ctot_inv); } if (m_jac_skip_third_bodies || m_multi_concm.empty()) { return; @@ -631,24 +635,23 @@ void BulkKinetics::process_ddC(StoichManagerN& stoich, const vector& in, Eigen::Map outM(m_rbuf1.data(), nReactions()); if (mass_action) { outM.fill(0.); - m_multi_concm.scale(in.data(), outM.data(), ctot_inv); + m_multi_concm.scale(in, asSpan(outM), ctot_inv); out += outM; } // derivatives due to reaction rates depending on third-body colliders if (!m_jac_skip_falloff) { - m_multi_concm.scaleM(in.data(), outM.data(), m_concm.data(), ctot_inv); + m_multi_concm.scaleM(in, asSpan(outM), m_concm, ctot_inv); for (auto& rates : m_rateHandlers) { // processing step assigns zeros to entries not dependent on M - rates->processRateConstants_ddM( - outM.data(), m_rfn.data(), m_jac_rtol_delta); + rates->processRateConstants_ddM(asSpan(outM), m_rfn, m_jac_rtol_delta); } out += outM; } } Eigen::SparseMatrix BulkKinetics::calculateCompositionDerivatives( - StoichManagerN& stoich, const vector& in, bool ddX) + StoichManagerN& stoich, span in, bool ddX) { Eigen::SparseMatrix out; vector& scaled = m_rbuf1; @@ -665,27 +668,26 @@ Eigen::SparseMatrix BulkKinetics::calculateCompositionDerivatives( // derivatives handled by StoichManagerN copy(scaled.begin(), scaled.end(), outV.begin()); - processThirdBodies(outV.data()); - out = stoich.derivatives(m_act_conc.data(), outV.data()); + processThirdBodies(outV); + out = stoich.derivatives(m_act_conc, outV); if (m_jac_skip_third_bodies || m_multi_concm.empty()) { return out; } // derivatives due to law of mass action copy(scaled.begin(), scaled.end(), outV.begin()); - stoich.multiply(m_act_conc.data(), outV.data()); + stoich.multiply(m_act_conc, outV); // derivatives due to reaction rates depending on third-body colliders if (!m_jac_skip_falloff) { for (auto& rates : m_rateHandlers) { // processing step does not modify entries not dependent on M - rates->processRateConstants_ddM( - outV.data(), m_rfn.data(), m_jac_rtol_delta, false); + rates->processRateConstants_ddM(outV, m_rfn, m_jac_rtol_delta, false); } } // derivatives handled by ThirdBodyCalc - out += m_multi_concm.derivatives(outV.data()); + out += m_multi_concm.derivatives(outV); return out; } diff --git a/src/kinetics/ElectronCollisionPlasmaRate.cpp b/src/kinetics/ElectronCollisionPlasmaRate.cpp index 5ef0e3b71fb..ce4a1ecbcd6 100644 --- a/src/kinetics/ElectronCollisionPlasmaRate.cpp +++ b/src/kinetics/ElectronCollisionPlasmaRate.cpp @@ -77,7 +77,8 @@ void ElectronCollisionPlasmaRate::getParameters(AnyMap& node) const { } void ElectronCollisionPlasmaRate::updateInterpolatedCrossSection( - const vector& sharedLevels) { + span sharedLevels) +{ m_crossSectionsInterpolated.clear(); m_crossSectionsInterpolated.reserve(sharedLevels.size()); for (double level : sharedLevels) { diff --git a/src/kinetics/Falloff.cpp b/src/kinetics/Falloff.cpp index cbd163f8b60..ba99b5e61fc 100644 --- a/src/kinetics/Falloff.cpp +++ b/src/kinetics/Falloff.cpp @@ -48,7 +48,8 @@ bool FalloffData::update(const ThermoPhase& phase, const Kinetics& kin) if (rho_m != molar_density || mf != m_state_mf_number) { molar_density = rho_m; m_state_mf_number = mf; - conc_3b = kin.thirdBodyConcentrations(); + auto concm = kin.thirdBodyConcentrations(); + conc_3b.assign(concm.begin(), concm.end()); changed = true; } return changed; @@ -116,7 +117,7 @@ void FalloffRate::setHighRate(const ArrheniusRate& high) m_highRate = std::move(_high); } -void FalloffRate::setFalloffCoeffs(const vector& c) +void FalloffRate::setFalloffCoeffs(span c) { if (c.size() != 0) { throw InputFileError("FalloffRate::setFalloffCoeffs", m_input, @@ -126,11 +127,6 @@ void FalloffRate::setFalloffCoeffs(const vector& c) m_valid = true; } -void FalloffRate::getFalloffCoeffs(vector& c) const -{ - c.clear(); -} - void FalloffRate::setParameters(const AnyMap& node, const UnitStack& rate_units) { ReactionRate::setParameters(node, rate_units); @@ -213,7 +209,7 @@ LindemannRate::LindemannRate(const AnyMap& node, const UnitStack& rate_units) } LindemannRate::LindemannRate(const ArrheniusRate& low, const ArrheniusRate& high, - const vector& c) + span c) : LindemannRate() { m_lowRate = low; @@ -228,7 +224,7 @@ TroeRate::TroeRate(const AnyMap& node, const UnitStack& rate_units) } TroeRate::TroeRate(const ArrheniusRate& low, const ArrheniusRate& high, - const vector& c) + span c) : TroeRate() { m_lowRate = low; @@ -236,7 +232,7 @@ TroeRate::TroeRate(const ArrheniusRate& low, const ArrheniusRate& high, setFalloffCoeffs(c); } -void TroeRate::setFalloffCoeffs(const vector& c) +void TroeRate::setFalloffCoeffs(span c) { if (c.size() != 3 && c.size() != 4) { throw InputFileError("TroeRate::setFalloffCoeffs", m_input, @@ -273,35 +269,31 @@ void TroeRate::setFalloffCoeffs(const vector& c) m_valid = true; } -void TroeRate::getFalloffCoeffs(vector& c) const +void TroeRate::getFalloffCoeffs(span c) const { - if (std::abs(m_t2) < SmallNumber) { - c.resize(3); - } else { - c.resize(4, 0.); - c[3] = m_t2; - } + checkArraySize("TroeRate::getFalloffCoeffs", c.size(), nParameters()); c[0] = m_a; c[1] = 1.0 / m_rt3; c[2] = 1.0 / m_rt1; + c[3] = m_t2; } -void TroeRate::updateTemp(double T, double* work) const +void TroeRate::updateTemp(double T, span work) const { double Fcent = (1.0 - m_a) * exp(-T*m_rt3) + m_a * exp(-T*m_rt1); if (m_t2) { Fcent += exp(- m_t2 / T); } - *work = log10(std::max(Fcent, SmallNumber)); + work[0] = log10(std::max(Fcent, SmallNumber)); } -double TroeRate::F(double pr, const double* work) const +double TroeRate::F(double pr, span work) const { double lpr = log10(std::max(pr,SmallNumber)); - double cc = -0.4 - 0.67 * (*work); - double nn = 0.75 - 1.27 * (*work); + double cc = -0.4 - 0.67 * work[0]; + double nn = 0.75 - 1.27 * work[0]; double f1 = (lpr + cc)/ (nn - 0.14 * (lpr + cc)); - double lgf = (*work) / (1.0 + f1 * f1); + double lgf = work[0] / (1.0 + f1 * f1); return pow(10.0, lgf); } @@ -350,7 +342,7 @@ SriRate::SriRate(const AnyMap& node, const UnitStack& rate_units) setParameters(node, rate_units); } -void SriRate::setFalloffCoeffs(const vector& c) +void SriRate::setFalloffCoeffs(span c) { if (c.size() != 3 && c.size() != 5) { throw InputFileError("SriRate::setFalloffCoeffs", m_input, @@ -380,34 +372,30 @@ void SriRate::setFalloffCoeffs(const vector& c) m_valid = true; } -void SriRate::getFalloffCoeffs(vector& c) const +void SriRate::getFalloffCoeffs(span c) const { - if (m_e < SmallNumber && std::abs(m_e - 1.) < SmallNumber) { - c.resize(3); - } else { - c.resize(5, 0.); - c[3] = m_d; - c[4] = m_e; - } + checkArraySize("SriRate::getFalloffCoeffs", c.size(), nParameters()); c[0] = m_a; c[1] = m_b; c[2] = m_c; + c[3] = m_d; + c[4] = m_e; } -void SriRate::updateTemp(double T, double* work) const +void SriRate::updateTemp(double T, span work) const { - *work = m_a * exp(- m_b / T); + work[0] = m_a * exp(- m_b / T); if (m_c != 0.0) { - *work += exp(- T/m_c); + work[0] += exp(- T/m_c); } work[1] = m_d * pow(T,m_e); } -double SriRate::F(double pr, const double* work) const +double SriRate::F(double pr, span work) const { double lpr = log10(std::max(pr,SmallNumber)); double xx = 1.0/(1.0 + lpr*lpr); - return pow(*work, xx) * work[1]; + return pow(work[0], xx) * work[1]; } void SriRate::setParameters(const AnyMap& node, const UnitStack& rate_units) @@ -459,7 +447,7 @@ TsangRate::TsangRate(const AnyMap& node, const UnitStack& rate_units) setParameters(node, rate_units); } -void TsangRate::setFalloffCoeffs(const vector& c) +void TsangRate::setFalloffCoeffs(span c) { if (c.size() != 1 && c.size() != 2) { throw InputFileError("TsangRate::init", m_input, @@ -477,30 +465,26 @@ void TsangRate::setFalloffCoeffs(const vector& c) m_valid = true; } -void TsangRate::getFalloffCoeffs(vector& c) const +void TsangRate::getFalloffCoeffs(span c) const { - if (std::abs(m_b) < SmallNumber) { - c.resize(1); - } else { - c.resize(2, 0.); - c[1] = m_b; - } + checkArraySize("TsangRate::getFalloffCoeffs", c.size(), nParameters()); c[0] = m_a; + c[1] = m_b; } -void TsangRate::updateTemp(double T, double* work) const +void TsangRate::updateTemp(double T, span work) const { double Fcent = m_a + (m_b * T); - *work = log10(std::max(Fcent, SmallNumber)); + work[0] = log10(std::max(Fcent, SmallNumber)); } -double TsangRate::F(double pr, const double* work) const +double TsangRate::F(double pr, span work) const { //identical to TroeRate::F double lpr = log10(std::max(pr,SmallNumber)); - double cc = -0.4 - 0.67 * (*work); - double nn = 0.75 - 1.27 * (*work); + double cc = -0.4 - 0.67 * work[0]; + double nn = 0.75 - 1.27 * work[0]; double f1 = (lpr + cc)/ (nn - 0.14 * (lpr + cc)); - double lgf = (*work) / (1.0 + f1 * f1); + double lgf = work[0] / (1.0 + f1 * f1); return pow(10.0, lgf); } diff --git a/src/kinetics/Group.cpp b/src/kinetics/Group.cpp index fc4b6e8a083..ff84791a476 100644 --- a/src/kinetics/Group.cpp +++ b/src/kinetics/Group.cpp @@ -36,7 +36,7 @@ void Group::validate() } } -std::ostream& Group::fmt(std::ostream& s, const vector& esymbols) const +std::ostream& Group::fmt(std::ostream& s, span esymbols) const { s << "("; bool first = true; diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 68addd8e34e..6bc21dcc395 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -62,7 +62,7 @@ void InterfaceKinetics::_update_rates_T() } // Use the stoichiometric manager to find deltaH for each reaction. - getReactionDelta(m_grt.data(), m_dH.data()); + getReactionDelta(m_grt, m_dH); m_temp = T; m_ROP_ok = false; @@ -73,7 +73,7 @@ void InterfaceKinetics::_update_rates_T() for (auto& rates : m_rateHandlers) { bool changed = rates->update(thermo(0), *this); if (changed) { - rates->getRateConstants(m_rfn.data()); + rates->getRateConstants(m_rfn); m_ROP_ok = false; m_redo_rates = true; } @@ -129,7 +129,7 @@ void InterfaceKinetics::updateKc() double rrt = 1.0 / thermo(0).RT(); // compute Delta mu^0 for all reversible reactions - getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data()); + getRevReactionDelta(m_mu0_Kc, m_rkcn); for (size_t i = 0; i < m_revindex.size(); i++) { size_t irxn = m_revindex[i]; @@ -165,19 +165,20 @@ void InterfaceKinetics::updateMu0() } } -void InterfaceKinetics::getEquilibriumConstants(double* kc) +void InterfaceKinetics::getEquilibriumConstants(span kc) { updateMu0(); double rrt = 1.0 / thermo(0).RT(); - std::fill(kc, kc + nReactions(), 0.0); - getReactionDelta(m_mu0_Kc.data(), kc); + std::fill(kc.begin(), kc.end(), 0.0); + getReactionDelta(m_mu0_Kc, kc); for (size_t i = 0; i < nReactions(); i++) { kc[i] = exp(-kc[i]*rrt); } } -void InterfaceKinetics::getFwdRateConstants(double* kfwd) +void InterfaceKinetics::getFwdRateConstants(span kfwd) { + checkArraySize("InterfaceKinetics::getFwdRateConstants", kfwd.size(), nReactions()); updateROP(); for (size_t i = 0; i < nReactions(); i++) { // base rate coefficient multiplied by perturbation factor @@ -185,11 +186,11 @@ void InterfaceKinetics::getFwdRateConstants(double* kfwd) } } -void InterfaceKinetics::getRevRateConstants(double* krev, bool doIrreversible) +void InterfaceKinetics::getRevRateConstants(span krev, bool doIrreversible) { getFwdRateConstants(krev); if (doIrreversible) { - getEquilibriumConstants(m_ropnet.data()); + getEquilibriumConstants(m_ropnet); for (size_t i = 0; i < nReactions(); i++) { krev[i] /= m_ropnet[i]; } @@ -221,16 +222,16 @@ void InterfaceKinetics::updateROP() } for (auto& rates : m_rateHandlers) { - rates->modifyRateConstants(m_ropf.data(), m_ropr.data()); + rates->modifyRateConstants(m_ropf, m_ropr); } // multiply ropf by the activity concentration reaction orders to obtain // the forward rates of progress. - m_reactantStoich.multiply(m_actConc.data(), m_ropf.data()); + m_reactantStoich.multiply(m_actConc, m_ropf); // For reversible reactions, multiply ropr by the activity concentration // products - m_revProductStoich.multiply(m_actConc.data(), m_ropr.data()); + m_revProductStoich.multiply(m_actConc, m_ropr); for (size_t j = 0; j != nReactions(); ++j) { m_ropnet[j] = m_ropf[j] - m_ropr[j]; @@ -286,7 +287,7 @@ void InterfaceKinetics::updateROP() m_ROP_ok = true; } -void InterfaceKinetics::getDeltaGibbs(double* deltaG) +void InterfaceKinetics::getDeltaGibbs(span deltaG) { // Get the chemical potentials of the species in the all of the phases used // in the kinetics mechanism @@ -296,15 +297,13 @@ void InterfaceKinetics::getDeltaGibbs(double* deltaG) } // Use the stoichiometric manager to find deltaG for each reaction. - getReactionDelta(m_mu.data(), m_rbuf.data()); - if (deltaG != 0 && (m_rbuf.data() != deltaG)) { - for (size_t j = 0; j < nReactions(); ++j) { - deltaG[j] = m_rbuf[j]; - } + getReactionDelta(m_mu, m_rbuf); + if (deltaG.data() != m_rbuf.data()) { + copy(m_rbuf.begin(), m_rbuf.end(), deltaG.begin()); } } -void InterfaceKinetics::getDeltaElectrochemPotentials(double* deltaM) +void InterfaceKinetics::getDeltaElectrochemPotentials(span deltaM) { // Get the chemical potentials of the species for (size_t n = 0; n < nPhases(); n++) { @@ -313,10 +312,10 @@ void InterfaceKinetics::getDeltaElectrochemPotentials(double* deltaM) } // Use the stoichiometric manager to find deltaG for each reaction. - getReactionDelta(m_grt.data(), deltaM); + getReactionDelta(m_grt, deltaM); } -void InterfaceKinetics::getDeltaEnthalpy(double* deltaH) +void InterfaceKinetics::getDeltaEnthalpy(span deltaH) { // Get the partial molar enthalpy of all species for (size_t n = 0; n < nPhases(); n++) { @@ -325,10 +324,10 @@ void InterfaceKinetics::getDeltaEnthalpy(double* deltaH) } // Use the stoichiometric manager to find deltaH for each reaction. - getReactionDelta(m_grt.data(), deltaH); + getReactionDelta(m_grt, deltaH); } -void InterfaceKinetics::getDeltaEntropy(double* deltaS) +void InterfaceKinetics::getDeltaEntropy(span deltaS) { // Get the partial molar entropy of all species in all of the phases for (size_t n = 0; n < nPhases(); n++) { @@ -337,10 +336,10 @@ void InterfaceKinetics::getDeltaEntropy(double* deltaS) } // Use the stoichiometric manager to find deltaS for each reaction. - getReactionDelta(m_grt.data(), deltaS); + getReactionDelta(m_grt, deltaS); } -void InterfaceKinetics::getDeltaSSGibbs(double* deltaGSS) +void InterfaceKinetics::getDeltaSSGibbs(span deltaGSS) { // Get the standard state chemical potentials of the species. This is the // array of chemical potentials at unit activity We define these here as the @@ -352,10 +351,10 @@ void InterfaceKinetics::getDeltaSSGibbs(double* deltaGSS) } // Use the stoichiometric manager to find deltaG for each reaction. - getReactionDelta(m_mu0.data(), deltaGSS); + getReactionDelta(m_mu0, deltaGSS); } -void InterfaceKinetics::getDeltaSSEnthalpy(double* deltaH) +void InterfaceKinetics::getDeltaSSEnthalpy(span deltaH) { // Get the standard state enthalpies of the species. This is the array of // chemical potentials at unit activity We define these here as the @@ -370,10 +369,10 @@ void InterfaceKinetics::getDeltaSSEnthalpy(double* deltaH) } // Use the stoichiometric manager to find deltaH for each reaction. - getReactionDelta(m_grt.data(), deltaH); + getReactionDelta(m_grt, deltaH); } -void InterfaceKinetics::getDeltaSSEntropy(double* deltaS) +void InterfaceKinetics::getDeltaSSEntropy(span deltaS) { // Get the standard state entropy of the species. We define these here as // the entropies of the pure species at the temperature and pressure of the @@ -387,7 +386,7 @@ void InterfaceKinetics::getDeltaSSEntropy(double* deltaS) } // Use the stoichiometric manager to find deltaS for each reaction. - getReactionDelta(m_grt.data(), deltaS); + getReactionDelta(m_grt, deltaS); } bool InterfaceKinetics::addReaction(shared_ptr r_base, bool resize) @@ -633,7 +632,7 @@ double InterfaceKinetics::interfaceCurrent(const size_t iphase) double dotProduct = 0.0; thermo(iphase).getCharges(charges); - getNetProductionRates(netProdRates.data()); + getNetProductionRates(netProdRates); for (size_t k = 0; k < thermo(iphase).nSpecies(); k++) { @@ -649,7 +648,7 @@ Eigen::SparseMatrix InterfaceKinetics::fwdRatesOfProgress_ddCi() assertDerivativesValid("InterfaceKinetics::fwdRatesOfProgress_ddCi"); // forward reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); return calculateCompositionDerivatives(m_reactantStoich, rop_rates); } @@ -659,8 +658,8 @@ Eigen::SparseMatrix InterfaceKinetics::revRatesOfProgress_ddCi() assertDerivativesValid("InterfaceKinetics::revRatesOfProgress_ddCi"); // reverse reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); - applyEquilibriumConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); + applyEquilibriumConstants(rop_rates); return calculateCompositionDerivatives(m_revProductStoich, rop_rates); } @@ -670,12 +669,12 @@ Eigen::SparseMatrix InterfaceKinetics::netRatesOfProgress_ddCi() assertDerivativesValid("InterfaceKinetics::netRatesOfProgress_ddCi"); // forward reaction rate coefficients vector& rop_rates = m_rbuf0; - getFwdRateConstants(rop_rates.data()); + getFwdRateConstants(rop_rates); Eigen::SparseMatrix jac = calculateCompositionDerivatives(m_reactantStoich, rop_rates); // reverse reaction rate coefficients - applyEquilibriumConstants(rop_rates.data()); + applyEquilibriumConstants(rop_rates); return jac - calculateCompositionDerivatives(m_revProductStoich, rop_rates); } @@ -703,12 +702,12 @@ void InterfaceKinetics::getDerivativeSettings(AnyMap& settings) const } Eigen::SparseMatrix InterfaceKinetics::calculateCompositionDerivatives( - StoichManagerN& stoich, const vector& in) + StoichManagerN& stoich, span in) { vector& outV = m_rbuf1; // derivatives handled by StoichManagerN copy(in.begin(), in.end(), outV.begin()); - return stoich.derivatives(m_actConc.data(), outV.data()); + return stoich.derivatives(m_actConc, outV); } void InterfaceKinetics::assertDerivativesValid(const string& name) @@ -720,7 +719,7 @@ void InterfaceKinetics::assertDerivativesValid(const string& name) } } -void InterfaceKinetics::applyEquilibriumConstants(double* rop) +void InterfaceKinetics::applyEquilibriumConstants(span rop) { // For reverse rates computed from thermochemistry, multiply the forward // rate coefficients by the reciprocals of the equilibrium constants diff --git a/src/kinetics/InterfaceRate.cpp b/src/kinetics/InterfaceRate.cpp index 98c1d654f68..8a241621ea4 100644 --- a/src/kinetics/InterfaceRate.cpp +++ b/src/kinetics/InterfaceRate.cpp @@ -20,14 +20,14 @@ void InterfaceData::update(double T) "Missing state information: 'InterfaceData' requires species coverages."); } -void InterfaceData::update(double T, const vector& values) +void InterfaceData::update(double T, span values) { warn_user("InterfaceData::update", "This method does not update the site density."); ReactionData::update(T); sqrtT = sqrt(T); if (coverages.size() == 0) { - coverages = values; + coverages.assign(values.begin(), values.end()); logCoverages.resize(values.size()); } else if (values.size() == coverages.size()) { std::copy(values.begin(), values.end(), coverages.begin()); @@ -209,12 +209,12 @@ void InterfaceRateBase::getCoverageDependencies(AnyMap& dependencies) const } void InterfaceRateBase::addCoverageDependence(const string& sp, double a, double m, - const vector& e) + span e) { if (std::find(m_cov.begin(), m_cov.end(), sp) == m_cov.end()) { m_cov.push_back(sp); m_ac.push_back(a); - m_ec.push_back(e); + m_ec.emplace_back(e.begin(), e.end()); m_mc.push_back(m); m_indices.clear(); } else { @@ -223,7 +223,7 @@ void InterfaceRateBase::addCoverageDependence(const string& sp, double a, double } } -void InterfaceRateBase::setSpecies(const vector& species) +void InterfaceRateBase::setSpecies(span species) { m_indices.clear(); for (size_t k = 0; k < m_cov.size(); k++) { diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index f266d264846..6f098b0c5c4 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -372,111 +372,126 @@ double Kinetics::productStoichCoeff(size_t kSpec, size_t irxn) const return m_productStoich.stoichCoeffs().coeff(kSpec, irxn); } -void Kinetics::getFwdRatesOfProgress(double* fwdROP) +void Kinetics::getFwdRatesOfProgress(span fwdROP) { + checkArraySize("Kinetics::getFwdRatesOfProgress", fwdROP.size(), nReactions()); updateROP(); - std::copy(m_ropf.begin(), m_ropf.end(), fwdROP); + std::copy(m_ropf.begin(), m_ropf.end(), fwdROP.begin()); } -void Kinetics::getRevRatesOfProgress(double* revROP) +void Kinetics::getRevRatesOfProgress(span revROP) { + checkArraySize("Kinetics::getRevRatesOfProgress", revROP.size(), nReactions()); updateROP(); - std::copy(m_ropr.begin(), m_ropr.end(), revROP); + std::copy(m_ropr.begin(), m_ropr.end(), revROP.begin()); } -void Kinetics::getNetRatesOfProgress(double* netROP) +void Kinetics::getNetRatesOfProgress(span netROP) { + checkArraySize("Kinetics::getNetRatesOfProgress", netROP.size(), nReactions()); updateROP(); - std::copy(m_ropnet.begin(), m_ropnet.end(), netROP); + std::copy(m_ropnet.begin(), m_ropnet.end(), netROP.begin()); } -void Kinetics::getReactionDelta(const double* prop, double* deltaProp) const +void Kinetics::getReactionDelta(span prop, + span deltaProp) const { - fill(deltaProp, deltaProp + nReactions(), 0.0); + checkArraySize("Kinetics::getReactionDelta", prop.size(), nTotalSpecies()); + checkArraySize("Kinetics::getReactionDelta", deltaProp.size(), nReactions()); + fill(deltaProp.begin(), deltaProp.end(), 0.0); // products add m_productStoich.incrementReactions(prop, deltaProp); // reactants subtract m_reactantStoich.decrementReactions(prop, deltaProp); } -void Kinetics::getRevReactionDelta(const double* prop, double* deltaProp) const +void Kinetics::getRevReactionDelta(span prop, + span deltaProp) const { - fill(deltaProp, deltaProp + nReactions(), 0.0); + checkArraySize("Kinetics::getRevReactionDelta", prop.size(), nTotalSpecies()); + checkArraySize("Kinetics::getRevReactionDelta", deltaProp.size(), nReactions()); + fill(deltaProp.begin(), deltaProp.end(), 0.0); // products add m_revProductStoich.incrementReactions(prop, deltaProp); // reactants subtract m_reactantStoich.decrementReactions(prop, deltaProp); } -void Kinetics::getCreationRates(double* cdot) +void Kinetics::getCreationRates(span cdot) { + checkArraySize("Kinetics::getCreationRates", cdot.size(), nTotalSpecies()); updateROP(); // zero out the output array - fill(cdot, cdot + m_kk, 0.0); + fill(cdot.begin(), cdot.end(), 0.0); // the forward direction creates product species - m_productStoich.incrementSpecies(m_ropf.data(), cdot); + m_productStoich.incrementSpecies(m_ropf, cdot); // the reverse direction creates reactant species - m_reactantStoich.incrementSpecies(m_ropr.data(), cdot); + m_reactantStoich.incrementSpecies(m_ropr, cdot); } -void Kinetics::getDestructionRates(double* ddot) +void Kinetics::getDestructionRates(span ddot) { + checkArraySize("Kinetics::getDestructionRates", ddot.size(), nTotalSpecies()); updateROP(); - fill(ddot, ddot + m_kk, 0.0); + fill(ddot.begin(), ddot.end(), 0.0); // the reverse direction destroys products in reversible reactions - m_revProductStoich.incrementSpecies(m_ropr.data(), ddot); + m_revProductStoich.incrementSpecies(m_ropr, ddot); // the forward direction destroys reactants - m_reactantStoich.incrementSpecies(m_ropf.data(), ddot); + m_reactantStoich.incrementSpecies(m_ropf, ddot); } -void Kinetics::getNetProductionRates(double* net) +void Kinetics::getNetProductionRates(span net) { + checkArraySize("Kinetics::getNetProductionRates", net.size(), nTotalSpecies()); updateROP(); - fill(net, net + m_kk, 0.0); + fill(net.begin(), net.end(), 0.0); // products are created for positive net rate of progress - m_productStoich.incrementSpecies(m_ropnet.data(), net); + m_productStoich.incrementSpecies(m_ropnet, net); // reactants are destroyed for positive net rate of progress - m_reactantStoich.decrementSpecies(m_ropnet.data(), net); + m_reactantStoich.decrementSpecies(m_ropnet, net); } -void Kinetics::getCreationRates_ddT(double* dwdot) +void Kinetics::getCreationRates_ddT(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getCreationRates_ddT", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); // the forward direction creates product species - getFwdRatesOfProgress_ddT(buf.data()); + getFwdRatesOfProgress_ddT(m_rbuf); out = m_productStoich.stoichCoeffs() * buf; // the reverse direction creates reactant species - getRevRatesOfProgress_ddT(buf.data()); + getRevRatesOfProgress_ddT(m_rbuf); out += m_reactantStoich.stoichCoeffs() * buf; } -void Kinetics::getCreationRates_ddP(double* dwdot) +void Kinetics::getCreationRates_ddP(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getCreationRates_ddP", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); // the forward direction creates product species - getFwdRatesOfProgress_ddP(buf.data()); + getFwdRatesOfProgress_ddP(m_rbuf); out = m_productStoich.stoichCoeffs() * buf; // the reverse direction creates reactant species - getRevRatesOfProgress_ddP(buf.data()); + getRevRatesOfProgress_ddP(m_rbuf); out += m_reactantStoich.stoichCoeffs() * buf; } -void Kinetics::getCreationRates_ddC(double* dwdot) +void Kinetics::getCreationRates_ddC(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getCreationRates_ddC", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); // the forward direction creates product species - getFwdRatesOfProgress_ddC(buf.data()); + getFwdRatesOfProgress_ddC(m_rbuf); out = m_productStoich.stoichCoeffs() * buf; // the reverse direction creates reactant species - getRevRatesOfProgress_ddC(buf.data()); + getRevRatesOfProgress_ddC(m_rbuf); out += m_reactantStoich.stoichCoeffs() * buf; } @@ -500,39 +515,42 @@ Eigen::SparseMatrix Kinetics::creationRates_ddCi() return jac; } -void Kinetics::getDestructionRates_ddT(double* dwdot) +void Kinetics::getDestructionRates_ddT(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getDestructionRates_ddT", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); // the reverse direction destroys products in reversible reactions - getRevRatesOfProgress_ddT(buf.data()); + getRevRatesOfProgress_ddT(m_rbuf); out = m_revProductStoich.stoichCoeffs() * buf; // the forward direction destroys reactants - getFwdRatesOfProgress_ddT(buf.data()); + getFwdRatesOfProgress_ddT(m_rbuf); out += m_reactantStoich.stoichCoeffs() * buf; } -void Kinetics::getDestructionRates_ddP(double* dwdot) +void Kinetics::getDestructionRates_ddP(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getDestructionRates_ddP", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); // the reverse direction destroys products in reversible reactions - getRevRatesOfProgress_ddP(buf.data()); + getRevRatesOfProgress_ddP(m_rbuf); out = m_revProductStoich.stoichCoeffs() * buf; // the forward direction destroys reactants - getFwdRatesOfProgress_ddP(buf.data()); + getFwdRatesOfProgress_ddP(m_rbuf); out += m_reactantStoich.stoichCoeffs() * buf; } -void Kinetics::getDestructionRates_ddC(double* dwdot) +void Kinetics::getDestructionRates_ddC(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getDestructionRates_ddC", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); // the reverse direction destroys products in reversible reactions - getRevRatesOfProgress_ddC(buf.data()); + getRevRatesOfProgress_ddC(m_rbuf); out = m_revProductStoich.stoichCoeffs() * buf; // the forward direction destroys reactants - getFwdRatesOfProgress_ddC(buf.data()); + getFwdRatesOfProgress_ddC(m_rbuf); out += m_reactantStoich.stoichCoeffs() * buf; } @@ -556,27 +574,30 @@ Eigen::SparseMatrix Kinetics::destructionRates_ddCi() return jac; } -void Kinetics::getNetProductionRates_ddT(double* dwdot) +void Kinetics::getNetProductionRates_ddT(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getNetProductionRates_ddT", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); - getNetRatesOfProgress_ddT(buf.data()); + getNetRatesOfProgress_ddT(m_rbuf); out = m_stoichMatrix * buf; } -void Kinetics::getNetProductionRates_ddP(double* dwdot) +void Kinetics::getNetProductionRates_ddP(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getNetProductionRates_ddP", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); - getNetRatesOfProgress_ddP(buf.data()); + getNetRatesOfProgress_ddP(m_rbuf); out = m_stoichMatrix * buf; } -void Kinetics::getNetProductionRates_ddC(double* dwdot) +void Kinetics::getNetProductionRates_ddC(span dwdot) { - Eigen::Map out(dwdot, m_kk); + checkArraySize("Kinetics::getNetProductionRates_ddC", dwdot.size(), nTotalSpecies()); + Eigen::Map out(dwdot.data(), m_kk); Eigen::Map buf(m_rbuf.data(), nReactions()); - getNetRatesOfProgress_ddC(buf.data()); + getNetRatesOfProgress_ddC(m_rbuf); out = m_stoichMatrix * buf; } diff --git a/src/kinetics/ReactionPath.cpp b/src/kinetics/ReactionPath.cpp index 139629ae04b..291e97613f1 100644 --- a/src/kinetics/ReactionPath.cpp +++ b/src/kinetics/ReactionPath.cpp @@ -142,7 +142,7 @@ void ReactionPathDiagram::add(shared_ptr d) add(*d.get()); } -void ReactionPathDiagram::findMajorPaths(double athreshold, size_t lda, double* a) +void ReactionPathDiagram::findMajorPaths(double athreshold, size_t lda, span a) { double netmax = 0.0; for (size_t n = 0; n < nNodes(); n++) { @@ -745,12 +745,12 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, return -1; } - s.getFwdRatesOfProgress(m_ropf.data()); - s.getRevRatesOfProgress(m_ropr.data()); + s.getFwdRatesOfProgress(m_ropf); + s.getRevRatesOfProgress(m_ropr); // species explicitly included or excluded - vector& in_nodes = r.included(); - vector& out_nodes = r.excluded(); + auto in_nodes = r.included(); + auto out_nodes = r.excluded(); vector status(s.nTotalSpecies(), 0); for (size_t ni = 0; ni < in_nodes.size(); ni++) { diff --git a/src/numerics/funcs.cpp b/src/numerics/funcs.cpp index c8cf46b7c03..201ccb8c49a 100644 --- a/src/numerics/funcs.cpp +++ b/src/numerics/funcs.cpp @@ -10,7 +10,7 @@ namespace Cantera { -double linearInterp(double x, const vector& xpts, const vector& fpts) +double linearInterp(double x, span xpts, span fpts) { AssertThrowMsg(!xpts.empty(), "linearInterp", "x data empty"); AssertThrowMsg(!fpts.empty(), "linearInterp", "f(x) data empty"); diff --git a/src/oneD/Boundary1D.cpp b/src/oneD/Boundary1D.cpp index eb7c19354e5..6ef07bf4509 100644 --- a/src/oneD/Boundary1D.cpp +++ b/src/oneD/Boundary1D.cpp @@ -718,7 +718,7 @@ void ReactingSurf1D::eval(size_t jg, double* xg, double* rg, m_flow_right->setGas(xg + rightloc, 0); } - m_kin->getNetProductionRates(m_work.data()); + m_kin->getNetProductionRates(m_work); double rs0 = 1.0/m_sphase->siteDensity(); if (m_enabled) { diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index ee5ddc78e8d..afb1e5e7cf1 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/src/thermo/EEDFTwoTermApproximation.cpp @@ -485,8 +485,8 @@ void EEDFTwoTermApproximation::calculateTotalCrossSection() m_totalCrossSectionCenter.assign(m_points, 0.0); m_totalCrossSectionEdge.assign(m_points + 1, 0.0); for (size_t k = 0; k < m_phase->nCollisions(); k++) { - auto& x = m_phase->collisionRate(k)->energyLevels(); - auto& y = m_phase->collisionRate(k)->crossSections(); + auto x = m_phase->collisionRate(k)->energyLevels(); + auto y = m_phase->collisionRate(k)->crossSections(); for (size_t i = 0; i < m_points; i++) { m_totalCrossSectionCenter[i] += m_X_targets[m_klocTargets[k]] * @@ -504,8 +504,8 @@ void EEDFTwoTermApproximation::calculateTotalElasticCrossSection() m_sigmaElastic.clear(); m_sigmaElastic.resize(m_points, 0.0); for (size_t k : m_phase->kElastic()) { - auto& x = m_phase->collisionRate(k)->energyLevels(); - auto& y = m_phase->collisionRate(k)->crossSections(); + auto x = m_phase->collisionRate(k)->energyLevels(); + auto y = m_phase->collisionRate(k)->crossSections(); // Note: // moleFraction(m_kTargets[k]) <=> m_X_targets[m_klocTargets[k]] double mass_ratio = ElectronMass / (m_phase->molecularWeight(m_kTargets[k]) / Avogadro); @@ -528,8 +528,8 @@ void EEDFTwoTermApproximation::setGridCache() m_i.resize(m_phase->nCollisions()); for (size_t k = 0; k < m_phase->nCollisions(); k++) { auto& collision = m_phase->collisionRate(k); - auto& x = collision->energyLevels(); - auto& y = collision->crossSections(); + auto x = collision->energyLevels(); + auto y = collision->crossSections(); vector eps1(m_points + 1); int shiftFactor = (collision->kind() == "ionization") ? 2 : 1; @@ -580,7 +580,8 @@ void EEDFTwoTermApproximation::setGridCache() } // construct sigma_offset - auto x_offset = collision->energyLevels(); + vector x_offset(collision->energyLevels().begin(), + collision->energyLevels().end()); for (auto& element : x_offset) { element -= collision->threshold(); } diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index a8857f12875..354448dbfa4 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -479,8 +479,10 @@ void PlasmaPhase::addCollision(shared_ptr collision) m_kInelastic.push_back(i); } - m_energyLevels.push_back(rate.energyLevels()); - m_crossSections.push_back(rate.crossSections()); + auto levels = rate.energyLevels(); + m_energyLevels.emplace_back(levels.begin(), levels.end()); + auto sections = rate.crossSections(); + m_crossSections.emplace_back(sections.begin(), sections.end()); m_eedfSolver->setGridCache(); } diff --git a/src/zeroD/ConstPressureMoleReactor.cpp b/src/zeroD/ConstPressureMoleReactor.cpp index 5f145d148a6..4afcbdd4a5c 100644 --- a/src/zeroD/ConstPressureMoleReactor.cpp +++ b/src/zeroD/ConstPressureMoleReactor.cpp @@ -64,7 +64,7 @@ void ConstPressureMoleReactor::eval(double time, double* LHS, double* RHS) auto imw = m_thermo->inverseMolecularWeights(); if (m_chem) { - m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } // external heat transfer diff --git a/src/zeroD/ConstPressureReactor.cpp b/src/zeroD/ConstPressureReactor.cpp index 1f1960a2cbd..39a78694e7d 100644 --- a/src/zeroD/ConstPressureReactor.cpp +++ b/src/zeroD/ConstPressureReactor.cpp @@ -72,7 +72,7 @@ void ConstPressureReactor::eval(double time, double* LHS, double* RHS) dmdt += mdot_surf; if (m_chem) { - m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } for (size_t k = 0; k < m_nsp; k++) { diff --git a/src/zeroD/FlowReactor.cpp b/src/zeroD/FlowReactor.cpp index 978c11b620c..c116f706736 100644 --- a/src/zeroD/FlowReactor.cpp +++ b/src/zeroD/FlowReactor.cpp @@ -55,7 +55,7 @@ void FlowReactor::getStateDae(double* y, double* ydot) y[3] = m_thermo->temperature(); if (m_chem) { - m_kin->getNetProductionRates(m_wdot.data()); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } // set the initial coverages @@ -179,7 +179,7 @@ void FlowReactor::evalDae(double time, double* y, double* ydot, double* residual m_thermo->getPartialMolarEnthalpies(m_hk); // get net production if (m_chem) { - m_kin->getNetProductionRates(m_wdot.data()); + m_kin->getNetProductionRates(m_wdot); } // set dphi/dz variables diff --git a/src/zeroD/IdealGasConstPressureMoleReactor.cpp b/src/zeroD/IdealGasConstPressureMoleReactor.cpp index a60e22fc4f4..187e3d15f52 100644 --- a/src/zeroD/IdealGasConstPressureMoleReactor.cpp +++ b/src/zeroD/IdealGasConstPressureMoleReactor.cpp @@ -62,7 +62,7 @@ void IdealGasConstPressureMoleReactor::eval(double time, double* LHS, double* RH auto imw = m_thermo->inverseMolecularWeights(); if (m_chem) { - m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } // external heat transfer @@ -161,7 +161,7 @@ void IdealGasConstPressureMoleReactor::getJacobianElements( // gas phase m_thermo->getPartialMolarCp(asSpan(specificHeat)); m_thermo->getPartialMolarEnthalpies(asSpan(enthalpy)); - m_kin->getNetProductionRates(netProductionRates.data()); + m_kin->getNetProductionRates(asSpan(netProductionRates)); // scale production rates by the volume for gas species for (size_t i = 0; i < m_nsp; i++) { netProductionRates[i] *= m_vol; diff --git a/src/zeroD/IdealGasConstPressureReactor.cpp b/src/zeroD/IdealGasConstPressureReactor.cpp index 56024b3cff6..92cf49fb303 100644 --- a/src/zeroD/IdealGasConstPressureReactor.cpp +++ b/src/zeroD/IdealGasConstPressureReactor.cpp @@ -68,7 +68,7 @@ void IdealGasConstPressureReactor::eval(double time, double* LHS, double* RHS) m_thermo->getPartialMolarEnthalpies(m_hk); if (m_chem) { - m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } // external heat transfer diff --git a/src/zeroD/IdealGasMoleReactor.cpp b/src/zeroD/IdealGasMoleReactor.cpp index 900c5794c2c..2218726a6b7 100644 --- a/src/zeroD/IdealGasMoleReactor.cpp +++ b/src/zeroD/IdealGasMoleReactor.cpp @@ -122,7 +122,7 @@ void IdealGasMoleReactor::eval(double time, double* LHS, double* RHS) auto imw = m_thermo->inverseMolecularWeights(); if (m_chem) { - m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } // external heat transfer and compression work @@ -228,7 +228,7 @@ void IdealGasMoleReactor::getJacobianElements(vector>& tr Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(m_nsp); // getting species data m_thermo->getPartialMolarIntEnergies(asSpan(internal_energy)); - m_kin->getNetProductionRates(netProductionRates.data()); + m_kin->getNetProductionRates(asSpan(netProductionRates)); m_thermo->getPartialMolarCp(asSpan(specificHeat)); // convert Cp to Cv for ideal gas as Cp - Cv = R for (size_t i = 0; i < m_nsp; i++) { diff --git a/src/zeroD/IdealGasReactor.cpp b/src/zeroD/IdealGasReactor.cpp index b3dbfb6604a..0ee976df789 100644 --- a/src/zeroD/IdealGasReactor.cpp +++ b/src/zeroD/IdealGasReactor.cpp @@ -66,7 +66,7 @@ void IdealGasReactor::eval(double time, double* LHS, double* RHS) auto Y = m_thermo->massFractions(); if (m_chem) { - m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } double mdot_surf = dot(m_sdot.begin(), m_sdot.end(), mw.begin()); diff --git a/src/zeroD/MoleReactor.cpp b/src/zeroD/MoleReactor.cpp index f6e5bb65321..b02986e9ea8 100644 --- a/src/zeroD/MoleReactor.cpp +++ b/src/zeroD/MoleReactor.cpp @@ -116,7 +116,7 @@ void MoleReactor::eval(double time, double* LHS, double* RHS) RHS[1] = m_vdot; if (m_chem) { - m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } // Energy equation. diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 0c87201cfeb..c648914cdae 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -138,7 +138,7 @@ void Reactor::eval(double time, double* LHS, double* RHS) RHS[1] = m_vdot; if (m_chem) { - m_kin->getNetProductionRates(&m_wdot[0]); // "omega dot" + m_kin->getNetProductionRates(m_wdot); // "omega dot" } for (size_t k = 0; k < m_nsp; k++) { diff --git a/src/zeroD/ReactorSurface.cpp b/src/zeroD/ReactorSurface.cpp index f78dc87b8d4..92638cb6d7e 100644 --- a/src/zeroD/ReactorSurface.cpp +++ b/src/zeroD/ReactorSurface.cpp @@ -104,7 +104,7 @@ void ReactorSurface::updateState(double* y) { m_surf->setCoveragesNoNorm(span(y, m_nsp)); m_thermo->setState_TP(m_reactors[0]->temperature(), m_reactors[0]->pressure()); - m_kinetics->getNetProductionRates(m_sdot.data()); + m_kinetics->getNetProductionRates(m_sdot); } void ReactorSurface::eval(double t, double* LHS, double* RHS) @@ -248,7 +248,7 @@ void MoleReactorSurface::updateState(double* y) } m_surf->setCoveragesNoNorm(m_cov_tmp); m_thermo->setState_TP(m_reactors[0]->temperature(), m_reactors[0]->pressure()); - m_kinetics->getNetProductionRates(m_sdot.data()); + m_kinetics->getNetProductionRates(m_sdot); } void MoleReactorSurface::eval(double t, double* LHS, double* RHS) diff --git a/test/clib/test_ctkin.cpp b/test/clib/test_ctkin.cpp index 3371d4e7d57..416ddafbe65 100644 --- a/test/clib/test_ctkin.cpp +++ b/test/clib/test_ctkin.cpp @@ -43,7 +43,7 @@ TEST(ctkin, kinetics) vector c_ropf(nr); kin_getFwdRatesOfProgress(kin, 325, c_ropf.data()); vector cpp_ropf(nr); - kinetics->getFwdRatesOfProgress(cpp_ropf.data()); + kinetics->getFwdRatesOfProgress(cpp_ropf); for (size_t n = 0; n < nr; n++) { ASSERT_NEAR(cpp_ropf[n], c_ropf[n], 1e-6); diff --git a/test/equil/equil_gas.cpp b/test/equil/equil_gas.cpp index d64da80a31a..2662570278d 100644 --- a/test/equil/equil_gas.cpp +++ b/test/equil/equil_gas.cpp @@ -9,6 +9,8 @@ #include "cantera/base/global.h" #include "cantera/base/utilities.h" +#include + using namespace Cantera; bool double_close(double expected, double actual, double tol) @@ -93,11 +95,13 @@ TEST_F(OverconstrainedEquil, BasisOptimize) mphase.addPhase(gas.get(), 10.0); mphase.init(); int usedZeroedSpecies = 0; - vector orderVectorSpecies; - vector orderVectorElements; + vector orderVectorSpecies(mphase.nSpecies()); + vector orderVectorElements(mphase.nElements()); + std::iota(orderVectorSpecies.begin(), orderVectorSpecies.end(), 0); + std::iota(orderVectorElements.begin(), orderVectorElements.end(), 0); bool doFormMatrix = true; - vector formRxnMatrix; + vector formRxnMatrix(mphase.nSpecies() * mphase.nElements()); size_t nc = BasisOptimize(&usedZeroedSpecies, doFormMatrix, &mphase, orderVectorSpecies, orderVectorElements, @@ -112,11 +116,13 @@ TEST_F(OverconstrainedEquil, DISABLED_BasisOptimize2) mphase.addPhase(gas.get(), 10.0); mphase.init(); int usedZeroedSpecies = 0; - vector orderVectorSpecies; - vector orderVectorElements; + vector orderVectorSpecies(mphase.nSpecies()); + vector orderVectorElements(mphase.nElements()); + std::iota(orderVectorSpecies.begin(), orderVectorSpecies.end(), 0); + std::iota(orderVectorElements.begin(), orderVectorElements.end(), 0); bool doFormMatrix = true; - vector formRxnMatrix; + vector formRxnMatrix(mphase.nSpecies() * mphase.nElements()); size_t nc = BasisOptimize(&usedZeroedSpecies, doFormMatrix, &mphase, orderVectorSpecies, orderVectorElements, diff --git a/test/general/test_composite.cpp b/test/general/test_composite.cpp index aa2093b834e..71bfa1d830d 100644 --- a/test/general/test_composite.cpp +++ b/test/general/test_composite.cpp @@ -17,8 +17,8 @@ TEST(Solution, clone_idealgas) EXPECT_EQ(soln->kinetics()->nReactions(), dup->kinetics()->nReactions()); EXPECT_EQ(soln->thermo()->name(), dup->thermo()->name()); EXPECT_NEAR(soln->thermo()->enthalpy_mass(), dup->thermo()->enthalpy_mass(), 1e-5); - soln->kinetics()->getFwdRateConstants(kf1.data()); - dup->kinetics()->getFwdRateConstants(kf2.data()); + soln->kinetics()->getFwdRateConstants(kf1); + dup->kinetics()->getFwdRateConstants(kf2); for (size_t i = 0; i < nrxn; i++) { EXPECT_NEAR(kf1[i], kf2[i], 1e-12*kf1[i]); } @@ -26,7 +26,7 @@ TEST(Solution, clone_idealgas) dup->transport()->thermalConductivity(), 1e-11); dup->thermo()->setState_TP(600, OneAtm); EXPECT_DOUBLE_EQ(soln->thermo()->temperature(), 400.0); - dup->kinetics()->getFwdRateConstants(kf2.data()); + dup->kinetics()->getFwdRateConstants(kf2); EXPECT_GT(kf2[2], kf1[2]); } @@ -41,14 +41,14 @@ TEST(Solution, clone_surf_complete) EXPECT_DOUBLE_EQ(soln2->adjacent("gas")->thermo()->pressure(), 3*OneAtm); vector ropf1(nrxn), ropf2(nrxn); - soln1->kinetics()->getFwdRatesOfProgress(ropf1.data()); - soln2->kinetics()->getFwdRatesOfProgress(ropf2.data()); + soln1->kinetics()->getFwdRatesOfProgress(ropf1); + soln2->kinetics()->getFwdRatesOfProgress(ropf2); for (size_t i = 0; i < nrxn; i++) { EXPECT_NEAR(ropf1[i], ropf2[i], 1e-12*ropf1[i]); soln2->kinetics()->setMultiplier(i, 0.0); } // Original mixture should still have same (non-zero) rates - soln1->kinetics()->getFwdRatesOfProgress(ropf1.data()); + soln1->kinetics()->getFwdRatesOfProgress(ropf1); for (size_t i = 0; i < nrxn; i++) { EXPECT_NEAR(ropf1[i], ropf2[i], 1e-12*ropf1[i]); } @@ -70,14 +70,14 @@ TEST(Solution, clone_surf_manual) EXPECT_DOUBLE_EQ(soln2->adjacent("gas")->thermo()->pressure(), 3*OneAtm); vector ropf1(nrxn), ropf2(nrxn); - soln1->kinetics()->getFwdRatesOfProgress(ropf1.data()); - soln2->kinetics()->getFwdRatesOfProgress(ropf2.data()); + soln1->kinetics()->getFwdRatesOfProgress(ropf1); + soln2->kinetics()->getFwdRatesOfProgress(ropf2); for (size_t i = 0; i < nrxn; i++) { EXPECT_NEAR(ropf1[i], ropf2[i], 1e-12*ropf1[i]); soln2->kinetics()->setMultiplier(i, 0.0); } // Original mixture should still have same (non-zero) rates - soln1->kinetics()->getFwdRatesOfProgress(ropf1.data()); + soln1->kinetics()->getFwdRatesOfProgress(ropf1); for (size_t i = 0; i < nrxn; i++) { EXPECT_NEAR(ropf1[i], ropf2[i], 1e-12*ropf1[i]); } diff --git a/test/general/test_serialization.cpp b/test/general/test_serialization.cpp index 58de4baa32d..97c4393a2ee 100644 --- a/test/general/test_serialization.cpp +++ b/test/general/test_serialization.cpp @@ -183,8 +183,8 @@ TEST(YamlWriter, reactions) ASSERT_EQ(kin1->nReactions(), kin2->nReactions()); vector kf1(kin1->nReactions()), kf2(kin1->nReactions()); - kin1->getFwdRateConstants(kf1.data()); - kin2->getFwdRateConstants(kf2.data()); + kin1->getFwdRateConstants(kf1); + kin2->getFwdRateConstants(kf2); for (size_t i = 0; i < kin1->nReactions(); i++) { EXPECT_NEAR(kf1[i], kf2[i], 1e-13 * kf1[i]) << "for reaction i = " << i; } @@ -212,8 +212,8 @@ TEST(YamlWriter, reaction_units_from_Yaml) ASSERT_EQ(kin1->nReactions(), kin2->nReactions()); vector kf1(kin1->nReactions()), kf2(kin1->nReactions()); - kin1->getFwdRateConstants(kf1.data()); - kin2->getFwdRateConstants(kf2.data()); + kin1->getFwdRateConstants(kf1); + kin2->getFwdRateConstants(kf2); for (size_t i = 0; i < kin1->nReactions(); i++) { EXPECT_NEAR(kf1[i], kf2[i], 1e-13 * kf1[i]) << "for reaction i = " << i; } @@ -239,8 +239,8 @@ TEST(YamlWriter, chebyshev_units_from_Yaml) ASSERT_EQ(kin1->nReactions(), kin2->nReactions()); vector kf1(kin1->nReactions()), kf2(kin1->nReactions()); - kin1->getFwdRateConstants(kf1.data()); - kin2->getFwdRateConstants(kf2.data()); + kin1->getFwdRateConstants(kf1); + kin2->getFwdRateConstants(kf2); for (size_t i = 0; i < kin1->nReactions(); i++) { EXPECT_NEAR(kf1[i], kf2[i], 1e-13 * kf1[i]) << "for reaction i = " << i; } @@ -311,16 +311,16 @@ TEST(YamlWriter, Interface) ASSERT_EQ(kin1->nReactions(), kin2->nReactions()); vector kf1(kin1->nReactions()), kf2(kin1->nReactions()); - kin1->getFwdRateConstants(kf1.data()); - kin2->getFwdRateConstants(kf2.data()); + kin1->getFwdRateConstants(kf1); + kin2->getFwdRateConstants(kf2); for (size_t i = 0; i < kin1->nReactions(); i++) { EXPECT_NEAR(kf1[i], kf2[i], 1e-13 * kf1[i]) << "for reaction i = " << i; } vector wdot1(kin1->nTotalSpecies()); vector wdot2(kin2->nTotalSpecies()); - kin1->getNetProductionRates(wdot1.data()); - kin2->getNetProductionRates(wdot2.data()); + kin1->getNetProductionRates(wdot1); + kin2->getNetProductionRates(wdot2); for (size_t i = 0; i < kin1->nTotalSpecies(); i++) { EXPECT_NEAR(wdot1[i], wdot2[i], 1e-13 * fabs(wdot1[i])) << "for species i = " << i; } @@ -348,16 +348,16 @@ TEST(YamlWriter, sofc) ASSERT_EQ(tpb_kin1->nReactions(), tpb_kin2->nReactions()); vector kf1(tpb_kin1->nReactions()), kf2(tpb_kin1->nReactions()); - tpb_kin1->getFwdRateConstants(kf1.data()); - tpb_kin2->getFwdRateConstants(kf2.data()); + tpb_kin1->getFwdRateConstants(kf1); + tpb_kin2->getFwdRateConstants(kf2); for (size_t i = 0; i < tpb_kin1->nReactions(); i++) { EXPECT_NEAR(kf1[i], kf2[i], 1e-13 * kf1[i]) << "for tpb reaction i = " << i; } vector wdot1(tpb_kin1->nTotalSpecies()); vector wdot2(tpb_kin2->nTotalSpecies()); - tpb_kin1->getNetProductionRates(wdot1.data()); - tpb_kin2->getNetProductionRates(wdot2.data()); + tpb_kin1->getNetProductionRates(wdot1); + tpb_kin2->getNetProductionRates(wdot2); for (size_t i = 0; i < tpb_kin1->nTotalSpecies(); i++) { EXPECT_NEAR(wdot1[i], wdot2[i], 1e-13 * fabs(wdot1[i])) << "for species i = " << i; } @@ -365,16 +365,16 @@ TEST(YamlWriter, sofc) ASSERT_EQ(ox_kin1->nReactions(), ox_kin2->nReactions()); kf1.resize(ox_kin1->nReactions()); kf2.resize(ox_kin1->nReactions()); - ox_kin1->getFwdRateConstants(kf1.data()); - ox_kin2->getFwdRateConstants(kf2.data()); + ox_kin1->getFwdRateConstants(kf1); + ox_kin2->getFwdRateConstants(kf2); for (size_t i = 0; i < ox_kin1->nReactions(); i++) { EXPECT_NEAR(kf1[i], kf2[i], 1e-13 * kf1[i]) << "for ox reaction i = " << i; } wdot1.resize(ox_kin1->nTotalSpecies()); wdot2.resize(ox_kin2->nTotalSpecies()); - ox_kin1->getNetProductionRates(wdot1.data()); - ox_kin2->getNetProductionRates(wdot2.data()); + ox_kin1->getNetProductionRates(wdot1); + ox_kin2->getNetProductionRates(wdot2); for (size_t i = 0; i < ox_kin1->nTotalSpecies(); i++) { EXPECT_NEAR(wdot1[i], wdot2[i], 1e-13 * fabs(wdot1[i])) << "for ox species i = " << i; } diff --git a/test/kinetics/kineticsFromScratch.cpp b/test/kinetics/kineticsFromScratch.cpp index 1a29824993a..b4d120099db 100644 --- a/test/kinetics/kineticsFromScratch.cpp +++ b/test/kinetics/kineticsFromScratch.cpp @@ -49,12 +49,12 @@ class KineticsFromScratch : public testing::Test vector k(1), k_ref(kin_ref->nReactions()); - kin.getFwdRateConstants(&k[0]); - kin_ref->getFwdRateConstants(&k_ref[0]); + kin.getFwdRateConstants(k); + kin_ref->getFwdRateConstants(k_ref); EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); - kin.getRevRateConstants(&k[0]); - kin_ref->getRevRateConstants(&k_ref[0]); + kin.getRevRateConstants(k); + kin_ref->getRevRateConstants(k_ref); EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); } }; @@ -637,12 +637,12 @@ class InterfaceKineticsFromScratch : public testing::Test vector k(1), k_ref(kin_ref->nReactions()); - kin.getFwdRateConstants(&k[0]); - kin_ref->getFwdRateConstants(&k_ref[0]); + kin.getFwdRateConstants(k); + kin_ref->getFwdRateConstants(k_ref); EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); - kin.getRevRateConstants(&k[0]); - kin_ref->getRevRateConstants(&k_ref[0]); + kin.getRevRateConstants(k); + kin_ref->getRevRateConstants(k_ref); EXPECT_DOUBLE_EQ(k_ref[iRef], k[0]); } }; @@ -730,32 +730,32 @@ class KineticsAddSpecies : public testing::Test vector k(kin.nReactions()), k_ref(kin_ref->nReactions()); vector w(kin.nTotalSpecies()), w_ref(kin_ref->nTotalSpecies()); - kin.getFwdRateConstants(k.data()); - kin_ref->getFwdRateConstants(k_ref.data()); + kin.getFwdRateConstants(k); + kin_ref->getFwdRateConstants(k_ref); for (size_t i = 0; i < kin.nReactions(); i++) { EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; } - kin.getFwdRatesOfProgress(k.data()); - kin_ref->getFwdRatesOfProgress(k_ref.data()); + kin.getFwdRatesOfProgress(k); + kin_ref->getFwdRatesOfProgress(k_ref); for (size_t i = 0; i < kin.nReactions(); i++) { EXPECT_DOUBLE_EQ(k_ref[i], k[i]) << "i = " << i << "; N = " << N; } - kin.getRevRateConstants(k.data()); - kin_ref->getRevRateConstants(k_ref.data()); + kin.getRevRateConstants(k); + kin_ref->getRevRateConstants(k_ref); for (size_t i = 0; i < kin.nReactions(); i++) { EXPECT_NEAR(k_ref[i], k[i], k_ref[i]*1e-12) << "i = " << i << "; N = " << N; } - kin.getRevRatesOfProgress(k.data()); - kin_ref->getRevRatesOfProgress(k_ref.data()); + kin.getRevRatesOfProgress(k); + kin_ref->getRevRatesOfProgress(k_ref); for (size_t i = 0; i < kin.nReactions(); i++) { EXPECT_NEAR(k_ref[i], k[i], k_ref[i]*1e-12) << "i = " << i << "; N = " << N; } - kin.getCreationRates(w.data()); - kin_ref->getCreationRates(w_ref.data()); + kin.getCreationRates(w); + kin_ref->getCreationRates(w_ref); for (size_t i = 0; i < kin.nTotalSpecies(); i++) { size_t iref = pp_ref->speciesIndex(p->speciesName(i)); EXPECT_NEAR(w_ref[iref], w[i], w_ref[iref]*1e-12) << "sp = " << p->speciesName(i) << "; N = " << N; diff --git a/test/kinetics/kineticsFromYaml.cpp b/test/kinetics/kineticsFromYaml.cpp index e35343d3783..237612652e4 100644 --- a/test/kinetics/kineticsFromYaml.cpp +++ b/test/kinetics/kineticsFromYaml.cpp @@ -335,11 +335,13 @@ TEST(Reaction, FalloffFromYaml2) EXPECT_DOUBLE_EQ(rate->highRate().preExponentialFactor(), 6e11); EXPECT_DOUBLE_EQ(rate->lowRate().preExponentialFactor(), 1.04e20); EXPECT_DOUBLE_EQ(rate->lowRate().activationEnergy(), 1600); - vector params; + vector params(rate->nParameters()); + ASSERT_EQ(params.size(), (size_t) 4); rate->getFalloffCoeffs(params); - ASSERT_EQ(params.size(), (size_t) 3); EXPECT_DOUBLE_EQ(params[0], 0.562); EXPECT_DOUBLE_EQ(params[1], 91.0); + EXPECT_DOUBLE_EQ(params[2], 5836.0); + EXPECT_DOUBLE_EQ(params[3], 0.0); EXPECT_EQ(R->input["source"].asString(), "somewhere"); } @@ -597,7 +599,7 @@ TEST(Reaction, PythonExtensibleRate) auto rate = R->rate(); EXPECT_EQ(rate->type(), "square-rate"); vector kf(sol->kinetics()->nReactions()); - sol->kinetics()->getFwdRateConstants(kf.data()); + sol->kinetics()->getFwdRateConstants(kf); EXPECT_DOUBLE_EQ(kf[0], 3.14 * 300 * 300); } @@ -710,8 +712,8 @@ TEST(Kinetics, ElectrochemFromYaml) auto kin = soln->kinetics(); soln->adjacent("graphite")->thermo()->setElectricPotential(0.4); vector ropf(kin->nReactions()), ropr(kin->nReactions()); - kin->getFwdRatesOfProgress(ropf.data()); - kin->getRevRatesOfProgress(ropr.data()); + kin->getFwdRatesOfProgress(ropf); + kin->getRevRatesOfProgress(ropr); EXPECT_NEAR(ropf[0], 0.279762338, 1e-8); EXPECT_NEAR(ropr[0], 0.045559670, 1e-8); @@ -835,10 +837,10 @@ class ReactionToYaml : public testing::Test vector kf(kin->nReactions()), kr(kin->nReactions()); vector ropf(kin->nReactions()), ropr(kin->nReactions()); - kin->getFwdRateConstants(kf.data()); - kin->getRevRateConstants(kr.data()); - kin->getFwdRatesOfProgress(ropf.data()); - kin->getRevRatesOfProgress(ropr.data()); + kin->getFwdRateConstants(kf); + kin->getRevRateConstants(kr); + kin->getFwdRatesOfProgress(ropf); + kin->getRevRatesOfProgress(ropr); EXPECT_NEAR(kf[iOld], kf[iNew], 1e-13 * kf[iOld]); EXPECT_NEAR(kr[iOld], kr[iNew], 1e-13 * kr[iOld]); EXPECT_NEAR(ropf[iOld], ropf[iNew], 1e-13 * ropf[iOld]); diff --git a/test/kinetics/pdep.cpp b/test/kinetics/pdep.cpp index 1876736391f..e04dd62056d 100644 --- a/test/kinetics/pdep.cpp +++ b/test/kinetics/pdep.cpp @@ -54,7 +54,7 @@ TEST_F(PdepTest, PlogLowPressure) // Test that P-log reactions have the right low-pressure limit set_TP(500.0, 1e-7); vector kf(7); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); // Pre-exponential factor decreases by 10^3 for second-order reaction // when converting from cm + mol to m + kmol @@ -74,7 +74,7 @@ TEST_F(PdepTest, PlogHighPressure) // Test that P-log reactions have the right high-pressure limit set_TP(500.0, 1e10); vector kf(7); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); // Pre-exponential factor decreases by 10^3 for second-order reaction // when converting from cm + mol to m + kmol @@ -91,7 +91,7 @@ TEST_F(PdepTest, PlogDuplicatePressures) set_TP(500.0, 1e10); vector kf(7); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); double kf1 = k(1.3700e+14, -0.79, 17603.0) + k(1.2800e+03, 1.71, 9774.0); double kf2 = k(-7.4100e+27, -5.54, 12108.0) + k(1.9000e+12, -0.29, 8306.0); @@ -105,7 +105,7 @@ TEST_F(PdepTest, PlogCornerCases) // is exactly of the specified interpolation values set_TP(500.0, 101325); vector kf(7); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); double kf0 = k(4.910800e+28, -4.8507, 24772.8); double kf1 = k(1.2600e+17, -1.83, 15003.0) + k(1.2300e+01, 2.68, 6335.0); @@ -120,7 +120,7 @@ TEST_F(PdepTest, PlogIntermediatePressure1) { set_TP(1100.0, 20*101325); vector ropf(7); - soln_->kinetics()->getFwdRatesOfProgress(&ropf[0]); + soln_->kinetics()->getFwdRatesOfProgress(ropf); // Expected rates computed using Chemkin // ROP increases by 10**3 when converting from mol/cm3 to kmol/m3 @@ -134,7 +134,7 @@ TEST_F(PdepTest, PlogIntermediatePressure2) { set_TP(1100.0, 0.5*101325); vector ropf(7); - soln_->kinetics()->getFwdRatesOfProgress(&ropf[0]); + soln_->kinetics()->getFwdRatesOfProgress(ropf); EXPECT_NEAR(5.244649e+02, ropf[0], 5e-2); EXPECT_NEAR(2.252537e+02, ropf[1], 2e-2); @@ -146,7 +146,7 @@ TEST_F(PdepTest, PlogIntermediatePressure3) { set_TP(800.0, 70*101325); vector ropf(7); - soln_->kinetics()->getFwdRatesOfProgress(&ropf[0]); + soln_->kinetics()->getFwdRatesOfProgress(ropf); EXPECT_NEAR(2.274501e+04, ropf[0], 1e+1); EXPECT_NEAR(2.307191e+05, ropf[1], 1e+2); @@ -160,7 +160,7 @@ TEST_F(PdepTest, ChebyshevIntermediate1) vector kf(7); set_TP(1100.0, 20 * 101325); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); // Expected rates computed using RMG-py EXPECT_NEAR(3.130698657e+06, kf[4], 1e-1); EXPECT_NEAR(1.187949573e+00, kf[5], 1e-7); @@ -177,7 +177,7 @@ TEST_F(PdepTest, ChebyshevIntermediate2) vector kf(7); set_TP(400.0, 0.1 * 101325); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); // Expected rates computed using RMG-py EXPECT_NEAR(1.713599902e+05, kf[4], 1e-3); EXPECT_NEAR(9.581780687e-24, kf[5], 1e-31); @@ -189,7 +189,7 @@ TEST_F(PdepTest, ChebyshevIntermediateROP) set_TP(1100.0, 30 * 101325); vector ropf(7); // Expected rates computed using Chemkin - soln_->kinetics()->getFwdRatesOfProgress(&ropf[0]); + soln_->kinetics()->getFwdRatesOfProgress(ropf); EXPECT_NEAR(4.552930e+03, ropf[4], 1e-1); EXPECT_NEAR(4.877390e-02, ropf[5], 1e-5); } @@ -200,22 +200,22 @@ TEST_F(PdepTest, ChebyshevEdgeCases) // Minimum P set_TP(500.0, 1000.0); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); EXPECT_NEAR(1.225785655e+06, kf[4], 1e-2); // Maximum P set_TP(500.0, 1.0e7); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); EXPECT_NEAR(1.580981157e+03, kf[4], 1e-5); // Minimum T set_TP(300.0, 101325); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); EXPECT_NEAR(5.405987017e+03, kf[4], 1e-5); // Maximum T set_TP(2000.0, 101325); - soln_->kinetics()->getFwdRateConstants(&kf[0]); + soln_->kinetics()->getFwdRateConstants(kf); EXPECT_NEAR(3.354054351e+07, kf[4], 1e-1); } diff --git a/test/kinetics/rates.cpp b/test/kinetics/rates.cpp index 7fe39cfce1c..e3ddbd44784 100644 --- a/test/kinetics/rates.cpp +++ b/test/kinetics/rates.cpp @@ -43,8 +43,8 @@ TEST_F(FracCoeffTest, RateConstants) { vector kf(kin->nReactions(), 0.0); vector kr(kin->nReactions(), 0.0); - kin->getFwdRateConstants(&kf[0]); - kin->getRevRateConstants(&kr[0]); + kin->getFwdRateConstants(kf); + kin->getRevRateConstants(kr); // sum of reaction orders is 1.0; kf has units of 1/s EXPECT_DOUBLE_EQ(1e13, kf[0]); @@ -64,8 +64,8 @@ TEST_F(FracCoeffTest, RatesOfProgress) vector conc(therm->nSpecies(), 0.0); vector ropf(kin->nReactions(), 0.0); therm->getConcentrations(conc); - kin->getFwdRateConstants(&kf[0]); - kin->getFwdRatesOfProgress(&ropf[0]); + kin->getFwdRateConstants(kf); + kin->getFwdRatesOfProgress(ropf); EXPECT_DOUBLE_EQ(conc[kH2O]*kf[0], ropf[0]); EXPECT_DOUBLE_EQ(pow(conc[kH2], 0.8)*conc[kO2]*pow(conc[kOH],2)*kf[1], @@ -77,9 +77,9 @@ TEST_F(FracCoeffTest, CreationDestructionRates) vector ropf(kin->nReactions(), 0.0); vector cdot(therm->nSpecies(), 0.0); vector ddot(therm->nSpecies(), 0.0); - kin->getFwdRatesOfProgress(&ropf[0]); - kin->getCreationRates(&cdot[0]); - kin->getDestructionRates(&ddot[0]); + kin->getFwdRatesOfProgress(ropf); + kin->getCreationRates(cdot); + kin->getDestructionRates(ddot); EXPECT_DOUBLE_EQ(ropf[0], ddot[kH2O]); EXPECT_DOUBLE_EQ(1.4*ropf[0], cdot[kH]); @@ -100,7 +100,7 @@ TEST_F(FracCoeffTest, EquilibriumConstants) vector Kc(kin->nReactions(), 0.0); vector mu0(therm->nSpecies(), 0.0); - kin->getEquilibriumConstants(&Kc[0]); + kin->getEquilibriumConstants(Kc); therm->getGibbs_ref(mu0); // at pRef double deltaG0_0 = 1.4 * mu0[kH] + 0.6 * mu0[kOH] + 0.2 * mu0[kO2] - mu0[kH2O]; @@ -131,7 +131,7 @@ class NegativePreexponentialFactor : public testing::Test ASSERT_EQ(12, (int) nSpec); ASSERT_EQ(12, (int) nRxn); vector wdot(nSpec); - soln->kinetics()->getNetProductionRates(&wdot[0]); + soln->kinetics()->getNetProductionRates(wdot); for (size_t i = 0; i < nSpec; i++) { EXPECT_NEAR(wdot_ref[i], wdot[i], std::abs(wdot_ref[i])*2e-5 + 1e-9); } @@ -141,8 +141,8 @@ class NegativePreexponentialFactor : public testing::Test vector ropf(nRxn); vector ropr(nRxn); - soln->kinetics()->getFwdRatesOfProgress(&ropf[0]); - soln->kinetics()->getRevRatesOfProgress(&ropr[0]); + soln->kinetics()->getFwdRatesOfProgress(ropf); + soln->kinetics()->getRevRatesOfProgress(ropr); for (size_t i = 0; i < nRxn; i++) { EXPECT_NEAR(ropf_ref[i], ropf[i], std::abs(ropf_ref[i])*2e-5 + 1e-9); EXPECT_NEAR(ropr_ref[i], ropr[i], std::abs(ropr_ref[i])*2e-5 + 1e-9); @@ -167,7 +167,7 @@ TEST(InterfaceReaction, CoverageDependency) { iface->thermo()->setState_TP(T, 101325); iface->thermo()->setCoveragesByName("PT(S):0.7, H(S):0.3"); vector kf(iface->kinetics()->nReactions()); - iface->kinetics()->getFwdRateConstants(&kf[0]); + iface->kinetics()->getFwdRateConstants(kf); EXPECT_NEAR(kf[0], 4.4579e7 * pow(T, 0.5), 1e-14*kf[0]); // Energies in XML file are converted from J/mol to J/kmol EXPECT_NEAR(kf[1], 3.7e20 * exp(-(67.4e6-6e6*0.3)/(GasConstant*T)), 1e-14*kf[1]); @@ -186,7 +186,7 @@ TEST(LinearBurkeRate, RateCombinations) auto kf = [&](size_t iRxn, double T, double P, const string& X="P1:1.0") { thermo.setState_TPX(T, P, X); vector rates(kin.nReactions()); - kin.getFwdRateConstants(rates.data()); + kin.getFwdRateConstants(rates); return rates[iRxn]; }; diff --git a/test/python/test_mixture.py b/test/python/test_mixture.py index d690278ef62..ad6c0214aa8 100644 --- a/test/python/test_mixture.py +++ b/test/python/test_mixture.py @@ -137,10 +137,10 @@ def test_species_moles(self, mix, phase1): assert mix.species_moles[2] == approx(S[2]) assert mix.phase_moles(0) == approx(sum(S[:phase1.n_species])) - with pytest.raises(ValueError): + with pytest.raises(ValueError, match="too small"): mix.species_moles = (1, 2, 3) - with pytest.raises(TypeError): + with pytest.raises(ValueError, match="too small"): mix.species_moles = 9 def test_element_moles(self, mix): diff --git a/test_problems/diamondSurf/runDiamond.cpp b/test_problems/diamondSurf/runDiamond.cpp index 59ca4a9f047..1213b123f6f 100644 --- a/test_problems/diamondSurf/runDiamond.cpp +++ b/test_problems/diamondSurf/runDiamond.cpp @@ -12,8 +12,6 @@ using namespace Cantera; int main(int argc, char** argv) { - int i, k; - try { CanteraError::setStackTraceDepth(20); auto iface = newInterface("diamond.yaml", "diamond_100"); @@ -62,14 +60,11 @@ int main(int argc, char** argv) AssertThrow(p == p, "main"); AssertThrowMsg(p == p, "main", "are you kidding"); - double src[20]; - for (i = 0; i < 20; i++) { - src[i] = 0.0; - } + vector src(ikin->nTotalSpecies()); ikin->getNetProductionRates(src); double sum = 0.0; double naH = 0.0; - for (k = 0; k < 13; k++) { + for (size_t k = 0; k < 13; k++) { if (k < 4) { naH = gas->nAtoms(k, 0); } else if (k == 4) { @@ -91,7 +86,7 @@ int main(int argc, char** argv) diamond100->getMoleFractions(x); cout << "Coverages:" << endl; - for (k = 0; k < 8; k++) { + for (size_t k = 0; k < 8; k++) { cout << k << " " << diamond100->speciesName(k) << " " << x[k] << endl; diff --git a/test_problems/stoichSolidKinetics/stoichSolidKinetics.cpp b/test_problems/stoichSolidKinetics/stoichSolidKinetics.cpp index 75159c66a33..25e52e37df2 100644 --- a/test_problems/stoichSolidKinetics/stoichSolidKinetics.cpp +++ b/test_problems/stoichSolidKinetics/stoichSolidKinetics.cpp @@ -20,19 +20,19 @@ void printValue(const string& label, double value) void printRates(InterfaceKinetics& iKin) { vector work(iKin.nReactions(), 0.0); - iKin.getNetRatesOfProgress(&work[0]); + iKin.getNetRatesOfProgress(work); printValue("ROP_net = ", work[0]); - iKin.getFwdRatesOfProgress(&work[0]); + iKin.getFwdRatesOfProgress(work); printValue("ROP_forward = ", work[0]); - iKin.getRevRatesOfProgress(&work[0]); + iKin.getRevRatesOfProgress(work); printValue("ROP_reverse = ", work[0]); - iKin.getFwdRateConstants(&work[0]); + iKin.getFwdRateConstants(work); printValue(" kfwd = ", work[0]); - iKin.getRevRateConstants(&work[0]); + iKin.getRevRateConstants(work); printValue(" krev = ", work[0]); } @@ -70,10 +70,10 @@ void testProblem() << " CaCO3(s) = CO2(g) + CaO(s) \n" << endl; for (int ktrials = 0; ktrials < 2; ktrials++) { - iKin.getDeltaSSGibbs(&work[0]); + iKin.getDeltaSSGibbs(work); printValue(" deltaGSS = ", work[0]); - iKin.getDeltaGibbs(&work[0]); + iKin.getDeltaGibbs(work); printValue(" deltaG = ", work[0]); gasTP->getChemPotentials(work); diff --git a/test_problems/surfSolverTest/surfaceSolver.cpp b/test_problems/surfSolverTest/surfaceSolver.cpp index 710c9705c9a..e6dab35886b 100644 --- a/test_problems/surfSolverTest/surfaceSolver.cpp +++ b/test_problems/surfSolverTest/surfaceSolver.cpp @@ -16,7 +16,8 @@ using namespace Cantera; #define MSSIZE 200 -void printGas(ostream& oooo, ThermoPhase* gasTP, InterfaceKinetics* iKin_ptr, double* src) +void printGas(ostream& oooo, ThermoPhase* gasTP, InterfaceKinetics* iKin_ptr, + span src) { vector x(gasTP->nSpecies()); vector C(gasTP->nSpecies()); @@ -48,8 +49,8 @@ void printGas(ostream& oooo, ThermoPhase* gasTP, InterfaceKinetics* iKin_ptr, do oooo << endl; } -void printBulk(ostream& oooo, - ThermoPhase* bulkPhaseTP, InterfaceKinetics* iKin_ptr, double* src) +void printBulk(ostream& oooo, ThermoPhase* bulkPhaseTP, InterfaceKinetics* iKin_ptr, + span src) { vector x(bulkPhaseTP->nSpecies()); vector C(bulkPhaseTP->nSpecies()); @@ -94,8 +95,8 @@ void printBulk(ostream& oooo, oooo << endl; } -void printSurf(ostream& oooo, - ThermoPhase* surfPhaseTP, InterfaceKinetics* iKin_ptr, double* src) +void printSurf(ostream& oooo, ThermoPhase* surfPhaseTP, InterfaceKinetics* iKin_ptr, + span src) { vector x(surfPhaseTP->nSpecies()); string surfParticlePhaseName = surfPhaseTP->name(); @@ -165,16 +166,16 @@ int main(int argc, char** argv) /* * Download the source terms for the rate equations */ - double src[MSSIZE]; + vector src(iKin_ptr->nTotalSpecies()); iKin_ptr->getNetProductionRates(src); printGas(cout, gasTP, iKin_ptr, src); printBulk(cout, bulkPhaseTP, iKin_ptr, src); - printSurf(cout, surfPhaseTP, iKin_ptr, src) ; + printSurf(cout, surfPhaseTP, iKin_ptr, src); printGas(ofile, gasTP, iKin_ptr, src); printBulk(ofile, bulkPhaseTP, iKin_ptr, src); - printSurf(ofile, surfPhaseTP, iKin_ptr, src) ; + printSurf(ofile, surfPhaseTP, iKin_ptr, src); /*****************************************************************************/ /* Now Tweak the inputs and do a quick calculation */ /****************************************************************************/