Vectorize Redlich-Kwong property calculations

include/cantera/numerics/eigen_dense.h

@@ -56,6 +56,29 @@ inline span<const double> asSpan(const Eigen::DenseBase<Derived>& v)
     return span<const double>(v.derived().data(), v.size());
 }
 
+//! Convenience wrapper for accessing std::vector as an Eigen VectorXd
+inline MappedVector asVectorXd(vector<double>& v)
+{
+    return MappedVector(v.data(), static_cast<Eigen::Index>(v.size()));
+}
+
+//! Convenience wrapper for accessing std::vector as an Eigen VectorXd
+inline ConstMappedVector asVectorXd(const vector<double>& v)
+{
+    return ConstMappedVector(v.data(), static_cast<Eigen::Index>(v.size()));
 }
 
+//! Convenience wrapper for accessing std::span as an Eigen VectorXd
+inline MappedVector asVectorXd(span<double> v)
+{
+    return MappedVector(v.data(), static_cast<Eigen::Index>(v.size()));
+}
+
+//! Convenience wrapper for accessing std::span as an Eigen VectorXd
+inline ConstMappedVector asVectorXd(span<const double> v)
+{
+    return ConstMappedVector(v.data(), static_cast<Eigen::Index>(v.size()));
+}
+
+}
 #endif

include/cantera/thermo/RedlichKwongMFTP.h

@@ -7,7 +7,7 @@
 #define CT_REDLICHKWONGMFTP_H
 
 #include "MixtureFugacityTP.h"
-#include "cantera/base/Array.h"
+#include "cantera/numerics/eigen_dense.h"
 
 namespace Cantera
 {
@@ -350,9 +350,9 @@ class RedlichKwongMFTP : public MixtureFugacityTP
     double hresid() const override;
 
 public:
-    double liquidVolEst(double TKelvin, double& pres) const override;
+    double liquidVolEst(double T, double& pres) const override;
     double densityCalc(double T, double pressure, int phase, double rhoguess) override;
-    double dpdVCalc(double TKelvin, double molarVol, double& presCalc) const override;
+    double dpdVCalc(double T, double molarVol, double& presCalc) const override;
 
     double isothermalCompressibility() const override;
     double thermalExpansionCoeff() const override;
@@ -378,11 +378,11 @@ class RedlichKwongMFTP : public MixtureFugacityTP
      * This function doesn't change the internal state of the object, so it is a
      * const function.  It does use the stored mole fractions in the object.
      *
-     * @param temp  Temperature (TKelvin)
+     * @param T  Temperature
      * @param aCalc (output)  Returns the a value
      * @param bCalc (output)  Returns the b value.
      */
-    void calculateAB(double temp, double& aCalc, double& bCalc) const;
+    void calculateAB(double T, double& aCalc, double& bCalc) const;
 
     // Special functions not inherited from MixtureFugacityTP
 
@@ -403,20 +403,30 @@ class RedlichKwongMFTP : public MixtureFugacityTP
 
     //! Value of b in the equation of state
     /*!
-     *  m_b is a function of the temperature and the mole fraction.
+     *  `m_bMix` is a function of the temperature and the mole fraction.
      */
-    double m_b_current = 0.0;
+    double m_bMix = 0.0;
 
     //! Value of a in the equation of state
     /*!
-     *  a_b is a function of the temperature and the mole fraction.
+     *  `m_aMix` is a function of the temperature and the mole fraction.
      */
-    double m_a_current = 0.0;
+    double m_aMix = 0.0;
 
-    vector<double> a_vec_Curr_;
-    vector<double> b_vec_Curr_;
+    //! Current temperature-dependent value of @f$ a_{jk} @f$ for each species pair.
+    //! Size #m_kk by #m_kk.
+    Eigen::MatrixXd m_a;
 
-    Array2D a_coeff_vec;
+    //! Coefficients @f$ b_k @f$ for each species. Size #m_kk.
+    Eigen::ArrayXd m_b;
+
+    //! Constant term in the expression for @f$ a_{jk} @f$ for each species pair.
+    //! Size #m_kk by #m_kk.
+    Eigen::MatrixXd m_a0;
+
+    //! Temperature coefficient in the expression for @f$ a_{jk} @f$.
+    //! Size #m_kk by #m_kk.
+    Eigen::MatrixXd m_a1;
 
     //! Explicitly-specified binary interaction parameters
     map<string, map<string, pair<double, double>>> m_binaryParameters;
@@ -429,13 +439,14 @@ class RedlichKwongMFTP : public MixtureFugacityTP
 
     double Vroot_[3] = {0.0, 0.0, 0.0};
 
-    //! Temporary storage - length = m_kk.
-    mutable vector<double> m_pp;
+    //! @f$ A_k = \sum_i X_i a_{ki} @f$. Length #m_kk.
+    mutable Eigen::ArrayXd m_Ak;
 
     // Partial molar volumes of the species
     mutable vector<double> m_partialMolarVolumes;
 
-    mutable vector<double> m_dAkdT; //!< Temporary storage for dA_k/dT; length #m_kk.
+    //! @f$ A'_k = dA_k/dT = \sum_i X_i a_{ki,1} @f$. Length #m_kk.
+    mutable Eigen::ArrayXd m_dAkdT;
 
     //! The derivative of the pressure wrt the volume
     /*!
@@ -456,7 +467,7 @@ class RedlichKwongMFTP : public MixtureFugacityTP
      *  Calculated at the current conditions. Total volume, temperature and
      *  other mole number kept constant
      */
-    mutable vector<double> dpdni_;
+    mutable Eigen::ArrayXd m_dpdni;
 
 private:
     //! Omega constant for a -> value of a in terms of critical properties

src/thermo/EEDFTwoTermApproximation.cpp

@@ -407,9 +407,7 @@ double EEDFTwoTermApproximation::electronMobility(const Eigen::VectorXd& f0)
         }
     }
     double nDensity = m_phase->molarDensity() * Avogadro;
-    auto f = ConstMappedVector(y.data(), y.size());
-    auto x = ConstMappedVector(m_gridEdge.data(), m_gridEdge.size());
-    return -1./3. * m_gamma * simpson(f, x) / nDensity;
+    return -1./3. * m_gamma * simpson(asVectorXd(y), asVectorXd(m_gridEdge)) / nDensity;
 }
 
 void EEDFTwoTermApproximation::initSpeciesIndexCrossSections()

src/thermo/PlasmaPhase.cpp

@@ -36,7 +36,7 @@ PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_)
     size_t nGridCells = 301;
     m_nPoints = nGridCells + 1;
     m_eedfSolver->setLinearGrid(kTe_max, nGridCells);
-    m_electronEnergyLevels = MappedVector(m_eedfSolver->getGridEdge().data(), m_nPoints);
+    m_electronEnergyLevels = asVectorXd(m_eedfSolver->getGridEdge());
     m_electronEnergyDist.setZero(m_nPoints);
 }
 
@@ -179,12 +179,10 @@ void PlasmaPhase::electronEnergyLevelChanged()
     m_levelNum++;
     // Cross sections are interpolated on the energy levels
     if (m_collisions.size() > 0) {
-        vector<double> energyLevels(m_nPoints);
-        MappedVector(energyLevels.data(), m_nPoints) = m_electronEnergyLevels;
         for (shared_ptr<Reaction> collision : m_collisions) {
             const auto& rate = boost::polymorphic_pointer_downcast
                 <ElectronCollisionPlasmaRate>(collision->rate());
-            rate->updateInterpolatedCrossSection(energyLevels);
+            rate->updateInterpolatedCrossSection(asSpan(m_electronEnergyLevels));
         }
     }
 }
@@ -222,10 +220,8 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(span<const double> levels,
 {
     m_distributionType = "discretized";
     m_nPoints = levels.size();
-    m_electronEnergyLevels =
-        Eigen::Map<const Eigen::ArrayXd>(levels.data(), m_nPoints);
-    m_electronEnergyDist =
-        Eigen::Map<const Eigen::ArrayXd>(dist.data(), m_nPoints);
+    m_electronEnergyLevels = asVectorXd(levels);
+    m_electronEnergyDist = asVectorXd(dist);
     checkElectronEnergyLevels();
     if (m_do_normalizeElectronEnergyDist) {
         normalizeElectronEnergyDistribution();