diff --git a/include/cantera/numerics/eigen_dense.h b/include/cantera/numerics/eigen_dense.h index 42175f5d0ee..efd2b87b73c 100644 --- a/include/cantera/numerics/eigen_dense.h +++ b/include/cantera/numerics/eigen_dense.h @@ -56,6 +56,29 @@ inline span asSpan(const Eigen::DenseBase& v) return span(v.derived().data(), v.size()); } +//! Convenience wrapper for accessing std::vector as an Eigen VectorXd +inline MappedVector asVectorXd(vector& v) +{ + return MappedVector(v.data(), static_cast(v.size())); +} + +//! Convenience wrapper for accessing std::vector as an Eigen VectorXd +inline ConstMappedVector asVectorXd(const vector& v) +{ + return ConstMappedVector(v.data(), static_cast(v.size())); } +//! Convenience wrapper for accessing std::span as an Eigen VectorXd +inline MappedVector asVectorXd(span v) +{ + return MappedVector(v.data(), static_cast(v.size())); +} + +//! Convenience wrapper for accessing std::span as an Eigen VectorXd +inline ConstMappedVector asVectorXd(span v) +{ + return ConstMappedVector(v.data(), static_cast(v.size())); +} + +} #endif diff --git a/include/cantera/thermo/RedlichKwongMFTP.h b/include/cantera/thermo/RedlichKwongMFTP.h index 72b82c5b0d6..f7221d284b3 100644 --- a/include/cantera/thermo/RedlichKwongMFTP.h +++ b/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 a_vec_Curr_; - vector 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>> 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 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 m_partialMolarVolumes; - mutable vector 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 dpdni_; + mutable Eigen::ArrayXd m_dpdni; private: //! Omega constant for a -> value of a in terms of critical properties diff --git a/src/thermo/EEDFTwoTermApproximation.cpp b/src/thermo/EEDFTwoTermApproximation.cpp index 07393be211d..7fa38c916e5 100644 --- a/src/thermo/EEDFTwoTermApproximation.cpp +++ b/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() diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 354448dbfa4..26f83aeb469 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/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 energyLevels(m_nPoints); - MappedVector(energyLevels.data(), m_nPoints) = m_electronEnergyLevels; for (shared_ptr collision : m_collisions) { const auto& rate = boost::polymorphic_pointer_downcast (collision->rate()); - rate->updateInterpolatedCrossSection(energyLevels); + rate->updateInterpolatedCrossSection(asSpan(m_electronEnergyLevels)); } } } @@ -222,10 +220,8 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(span levels, { m_distributionType = "discretized"; m_nPoints = levels.size(); - m_electronEnergyLevels = - Eigen::Map(levels.data(), m_nPoints); - m_electronEnergyDist = - Eigen::Map(dist.data(), m_nPoints); + m_electronEnergyLevels = asVectorXd(levels); + m_electronEnergyDist = asVectorXd(dist); checkElectronEnergyLevels(); if (m_do_normalizeElectronEnergyDist) { normalizeElectronEnergyDistribution(); diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index 2e8472adc8c..9b363111591 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -33,9 +33,8 @@ void RedlichKwongMFTP::setSpeciesCoeffs(const string& species, m_formTempParam = 1; // expression is temperature-dependent } - size_t counter = k + m_kk * k; - a_coeff_vec(0, counter) = a0; - a_coeff_vec(1, counter) = a1; + m_a0(k, k) = a0; + m_a1(k, k) = a1; // standard mixing rule for cross-species interaction term for (size_t j = 0; j < m_kk; j++) { @@ -43,23 +42,19 @@ void RedlichKwongMFTP::setSpeciesCoeffs(const string& species, continue; } - // a_coeff_vec(0) is initialized to NaN to mark uninitialized species - if (isnan(a_coeff_vec(0, j + m_kk * j))) { + // m_a0 is initialized to NaN to mark uninitialized species + if (isnan(m_a0(j, j))) { // The diagonal element of the jth species has not yet been defined. continue; - } else if (isnan(a_coeff_vec(0, j + m_kk * k))) { + } else if (isnan(m_a0(j, k))) { // Only use the mixing rules if the off-diagonal element has not already been defined by a // user-specified crossFluidParameters entry: - double a0kj = sqrt(a_coeff_vec(0, j + m_kk * j) * a0); - double a1kj = sqrt(a_coeff_vec(1, j + m_kk * j) * a1); - a_coeff_vec(0, j + m_kk * k) = a0kj; - a_coeff_vec(1, j + m_kk * k) = a1kj; - a_coeff_vec(0, k + m_kk * j) = a0kj; - a_coeff_vec(1, k + m_kk * j) = a1kj; + m_a0(j, k) = m_a0(k, j) = sqrt(m_a0(j, j) * a0); + m_a1(j, k) = m_a1(k, j) = sqrt(m_a1(j, j) * a1); } } - a_coeff_vec.getRow(0, a_vec_Curr_); - b_vec_Curr_[k] = b; + m_a = m_a0; + m_b[k] = b; } void RedlichKwongMFTP::setBinaryCoeffs(const string& species_i, const string& species_j, @@ -74,11 +69,9 @@ void RedlichKwongMFTP::setBinaryCoeffs(const string& species_i, const string& sp m_binaryParameters[species_i][species_j] = {a0, a1}; m_binaryParameters[species_j][species_i] = {a0, a1}; - size_t counter1 = ki + m_kk * kj; - size_t counter2 = kj + m_kk * ki; - a_coeff_vec(0, counter1) = a_coeff_vec(0, counter2) = a0; - a_coeff_vec(1, counter1) = a_coeff_vec(1, counter2) = a1; - a_vec_Curr_[counter1] = a_vec_Curr_[counter2] = a0; + m_a0(ki, kj) = m_a0(kj, ki) = a0; + m_a1(ki, kj) = m_a1(kj, ki) = a1; + m_a(ki, kj) = m_a(kj, ki) = a0; } // ------------Molar Thermodynamic Properties ------------------------- @@ -87,36 +80,36 @@ double RedlichKwongMFTP::cp_mole() const { _updateReferenceStateThermo(); double cpref = GasConstant * mean_X(m_cp0_R); - if (m_b_current == 0.0) { + if (m_bMix == 0.0) { return cpref; } - double TKelvin = temperature(); - double sqt = sqrt(TKelvin); + double T = temperature(); + double sqt = sqrt(T); double mv = molarVolume(); - double vpb = mv + m_b_current; + double vpb = mv + m_bMix; pressureDerivatives(); double dadt = da_dt(); - double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0; - double dHdT_V = (cpref + mv * dpdT_ - GasConstant - 1.0 / (2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv) * fac - +1.0/(m_b_current * sqt) * log(vpb/mv) * (-0.5 * dadt)); - return dHdT_V - (mv + TKelvin * dpdT_ / dpdV_) * dpdT_; + double fac = T * dadt - 3.0 * m_aMix / 2.0; + double dHdT_V = cpref - 0.5 * dadt * log(vpb/mv) / (m_bMix * sqt) + + mv * dpdT_ - GasConstant - log(vpb/mv) * fac / (2.0 * m_bMix * T * sqt); + return dHdT_V - (mv + T * dpdT_ / dpdV_) * dpdT_; } double RedlichKwongMFTP::cv_mole() const { _updateReferenceStateThermo(); double cvref = GasConstant * (mean_X(m_cp0_R) - 1.0); - if (m_b_current == 0.0) { + if (m_bMix == 0.0) { return cvref; } - double TKelvin = temperature(); - double sqt = sqrt(TKelvin); + double T = temperature(); + double sqt = sqrt(T); double mv = molarVolume(); - double vpb = mv + m_b_current; + double vpb = mv + m_bMix; double dadt = da_dt(); - double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0; - return (cvref - 1.0/(2.0 * m_b_current * TKelvin * sqt) * log(vpb/mv)*fac - +1.0/(m_b_current * sqt) * log(vpb/mv)*(-0.5*dadt)); + double fac = T * dadt - 3.0 * m_aMix / 2.0; + return (cvref - 1.0/(2.0 * m_bMix * T * sqt) * log(vpb/mv)*fac + +1.0/(m_bMix * sqt) * log(vpb/mv)*(-0.5*dadt)); } double RedlichKwongMFTP::pressure() const @@ -126,7 +119,8 @@ double RedlichKwongMFTP::pressure() const // Get a copy of the private variables stored in the State object double T = temperature(); double molarV = meanMolecularWeight() / density(); - double pp = GasConstant * T/(molarV - m_b_current) - m_a_current/(sqrt(T) * molarV * (molarV + m_b_current)); + double pp = GasConstant * T/(molarV - m_bMix) + - m_aMix/(sqrt(T) * molarV * (molarV + m_bMix)); return pp; } @@ -142,35 +136,19 @@ void RedlichKwongMFTP::getActivityCoefficients(span ac) const for (size_t k = 0; k < m_kk; k++) { ac[k] = 1.0; } - if (m_b_current == 0.0) { + if (m_bMix == 0.0) { return; } double mv = molarVolume(); double sqt = sqrt(temperature()); - double vpb = mv + m_b_current; - double vmb = mv - m_b_current; - - for (size_t k = 0; k < m_kk; k++) { - m_pp[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk*i; - m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; - } - } - double pres = pressure(); - - for (size_t k = 0; k < m_kk; k++) { - ac[k] = (- RT() * log(pres * mv / RT()) - + RT() * log(mv / vmb) - + RT() * b_vec_Curr_[k] / vmb - - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv) - + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) - - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb) - ); - } - for (size_t k = 0; k < m_kk; k++) { - ac[k] = exp(ac[k]/RT()); - } + double vpb = mv + m_bMix; + double vmb = mv - m_bMix; + double logv = log(vpb / mv); + Eigen::ArrayXd g = RT() * (log(RT() / (vmb * pressure())) + m_b / vmb) + - 2.0 * m_Ak / (m_bMix * sqt) * logv + + m_aMix * m_b / (m_bMix * m_bMix * sqt) * logv + - m_aMix * m_b / (m_bMix * sqt * vpb); + MappedVector(ac.data(), m_kk) = (g / RT()).exp(); } // ---- Partial Molar Properties of the Solution ----------------- @@ -178,37 +156,23 @@ void RedlichKwongMFTP::getActivityCoefficients(span ac) const void RedlichKwongMFTP::getChemPotentials(span mu) const { getStandardChemPotentials(mu); - for (size_t k = 0; k < m_kk; k++) { - double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] += RT() * log(xx); - } - if (m_b_current == 0.0) { + Eigen::Map muVec(mu.data(), m_kk); + auto x = asVectorXd(moleFractions_); + muVec += RT() * x.array().max(SmallNumber).log(); + if (m_bMix == 0.0) { return; } double mv = molarVolume(); double sqt = sqrt(temperature()); - double vpb = mv + m_b_current; - double vmb = mv - m_b_current; - - for (size_t k = 0; k < m_kk; k++) { - m_pp[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk*i; - m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; - } - } - double pres = pressure(); - - for (size_t k = 0; k < m_kk; k++) { - mu[k] += (- RT() * log(pres * mv / RT()) - + RT() * log(mv / vmb) - + RT() * b_vec_Curr_[k] / vmb - - 2.0 * m_pp[k] / (m_b_current * sqt) * log(vpb/mv) - + m_a_current * b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) - - m_a_current / (m_b_current * sqt) * (b_vec_Curr_[k]/vpb) - ); - } + double vpb = mv + m_bMix; + double vmb = mv - m_bMix; + double logv = log(vpb / mv); + muVec += RT() * log(RT() / (vmb * pressure())) + + RT() * m_b / vmb + - 2.0 * m_Ak / (m_bMix * sqt) * logv + + m_aMix * m_b / (m_bMix * m_bMix * sqt) * logv + - m_aMix * m_b / (m_bMix * sqt * vpb); } void RedlichKwongMFTP::getPartialMolarEnthalpies(span hbar) const @@ -216,101 +180,64 @@ void RedlichKwongMFTP::getPartialMolarEnthalpies(span hbar) const // First we get the reference state contributions getEnthalpy_RT_ref(hbar); scale(hbar.begin(), hbar.end(), hbar.begin(), RT()); - if (m_b_current == 0.0) { + if (m_bMix == 0.0) { return; } - // We calculate dpdni_ - double TKelvin = temperature(); + // We calculate m_dpdni + double T = temperature(); double mv = molarVolume(); - double sqt = sqrt(TKelvin); - double vpb = mv + m_b_current; - double vmb = mv - m_b_current; - for (size_t k = 0; k < m_kk; k++) { - m_pp[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk*i; - m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; - } - } - for (size_t k = 0; k < m_kk; k++) { - dpdni_[k] = RT()/vmb + RT() * b_vec_Curr_[k] / (vmb * vmb) - 2.0 * m_pp[k] / (sqt * mv * vpb) - + m_a_current * b_vec_Curr_[k]/(sqt * mv * vpb * vpb); - } + double sqt = sqrt(T); + double vpb = mv + m_bMix; + double vmb = mv - m_bMix; + m_dpdni = RT() / vmb * (1.0 + m_b / vmb) + - 2.0 * m_Ak / (sqt * mv * vpb) + + m_aMix * m_b / (sqt * mv * vpb * vpb); double dadt = da_dt(); - double fac = TKelvin * dadt - 3.0 * m_a_current / 2.0; - - for (size_t k = 0; k < m_kk; k++) { - m_workS[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk*i; - m_workS[k] += 2.0 * moleFractions_[i] * TKelvin * a_coeff_vec(1,counter) - 3.0 * moleFractions_[i] * a_vec_Curr_[counter]; - } - } + double fac = T * dadt - 3.0 * m_aMix / 2.0; + Eigen::ArrayXd Sk = 2.0 * T * m_dAkdT - 3.0 * m_Ak; pressureDerivatives(); - double fac2 = mv + TKelvin * dpdT_ / dpdV_; - for (size_t k = 0; k < m_kk; k++) { - double hE_v = (mv * dpdni_[k] - RT() - b_vec_Curr_[k]/ (m_b_current * m_b_current * sqt) * log(vpb/mv)*fac - + 1.0 / (m_b_current * sqt) * log(vpb/mv) * m_workS[k] - + b_vec_Curr_[k] / vpb / (m_b_current * sqt) * fac); - hbar[k] = hbar[k] + hE_v; - hbar[k] -= fac2 * dpdni_[k]; - } + double fac2 = mv + T * dpdT_ / dpdV_; + double logv = log(vpb / mv); + Eigen::ArrayXd hE_v = mv * m_dpdni + - RT() + - m_b / (m_bMix * m_bMix * sqt) * logv * fac + + (1.0 / (m_bMix * sqt) * logv) * Sk + + m_b / (vpb * m_bMix * sqt) * fac; + Eigen::Map(hbar.data(), m_kk) += hE_v - fac2 * m_dpdni; } void RedlichKwongMFTP::getPartialMolarEntropies(span sbar) const { getEntropy_R_ref(sbar); - scale(sbar.begin(), sbar.end(), sbar.begin(), GasConstant); - double TKelvin = temperature(); - double sqt = sqrt(TKelvin); + double T = temperature(); + double sqt = sqrt(T); double mv = molarVolume(); double pres = pressure(); - double logPres = log(pres / refPressure()); - for (size_t k = 0; k < m_kk; k++) { - double xx = std::max(SmallNumber, moleFraction(k)); - sbar[k] -= GasConstant * (log(xx) + logPres); - } - if (m_b_current == 0.0) { + MappedVector sbarVec(sbar.data(), m_kk); + sbarVec -= asVectorXd(moleFractions_).array().max(SmallNumber).log().matrix(); + sbarVec.array() -= log(pres / refPressure()); + sbarVec *= GasConstant; + if (m_bMix == 0.0) { return; } - for (size_t k = 0; k < m_kk; k++) { - m_pp[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk*i; - m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; - } - } - for (size_t k = 0; k < m_kk; k++) { - m_workS[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk*i; - m_workS[k] += moleFractions_[i] * a_coeff_vec(1,counter); - } - } - - double dadt = da_dt(); - double fac = dadt - m_a_current / (2.0 * TKelvin); - double vmb = mv - m_b_current; - double vpb = mv + m_b_current; - for (size_t k = 0; k < m_kk; k++) { - sbar[k] += (GasConstant * log(pres * mv / RT()) - - GasConstant - - GasConstant * log(mv/vmb) - - GasConstant * b_vec_Curr_[k]/vmb - - m_pp[k]/(m_b_current * TKelvin * sqt) * log(vpb/mv) - + 2.0 * m_workS[k]/(m_b_current * sqt) * log(vpb/mv) - - b_vec_Curr_[k] / (m_b_current * m_b_current * sqt) * log(vpb/mv) * fac - + 1.0 / (m_b_current * sqt) * b_vec_Curr_[k] / vpb * fac); - } + double fac = da_dt() - m_aMix / (2.0 * T); + double vmb = mv - m_bMix; + double vpb = mv + m_bMix; + double logv = log(vpb / mv); + Eigen::ArrayXd sdep = GasConstant * (log(RT() / (pres * vmb)) + 1) + + GasConstant * m_b / vmb + + m_Ak / (m_bMix * T * sqt) * logv + - 2.0 * m_dAkdT / (m_bMix * sqt) * logv + + m_b / (m_bMix * m_bMix * sqt) * logv * fac + - m_b / (m_bMix * sqt * vpb) * fac; + sbarVec -= sdep.matrix(); pressureDerivatives(); getPartialMolarVolumes(m_partialMolarVolumes); - for (size_t k = 0; k < m_kk; k++) { - sbar[k] += m_partialMolarVolumes[k] * dpdT_; - } + sbarVec += dpdT_ * asVectorXd(m_partialMolarVolumes); } void RedlichKwongMFTP::getPartialMolarIntEnergies(span ubar) const @@ -329,104 +256,67 @@ void RedlichKwongMFTP::getPartialMolarIntEnergies_TV(span utilde) const checkArraySize("RedlichKwongMFTP::getPartialMolarIntEnergies_TV", utilde.size(), m_kk); _updateReferenceStateThermo(); - double T = temperature(); - double sqt = sqrt(T); - double mv = molarVolume(); - double b = m_b_current; - double a = m_a_current; - double dadt = da_dt(); + Eigen::Map utildeVec(utilde.data(), m_kk); + utildeVec = RT() * (asVectorXd(m_h0_RT).array() - 1.0); - for (size_t k = 0; k < m_kk; k++) { - utilde[k] = RT() * (m_h0_RT[k] - 1.0); - } - - if (fabs(b) < SmallNumber) { + if (fabs(m_bMix) < SmallNumber) { return; } - // A_k = sum_i x_i a_ki and A'_k = dA_k/dT = sum_i x_i a_ki,1 - for (size_t k = 0; k < m_kk; k++) { - m_pp[k] = 0.0; - m_dAkdT[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk * i; - m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; - m_dAkdT[k] += moleFractions_[i] * a_coeff_vec(1, counter); - } - } - - double vpb = mv + b; + double T = temperature(); + double mv = molarVolume(); + double vpb = mv + m_bMix; double logv = log(vpb / mv); - double F = T * dadt - 1.5 * a; - double C = -logv / b + 1.0 / vpb; - double pref = 1.0 / (b * sqt); + double F = T * da_dt() - 1.5 * m_aMix; + double C = -logv / m_bMix + 1.0 / vpb; + double pref = 1.0 / (m_bMix * sqrt(T)); - for (size_t k = 0; k < m_kk; k++) { - double Sk = 2.0 * T * m_dAkdT[k] - 3.0 * m_pp[k]; - double ures = pref * (logv * Sk + b_vec_Curr_[k] * F * C); - utilde[k] += ures; - } + Eigen::ArrayXd Sk = 2.0 * T * m_dAkdT - 3.0 * m_Ak; + utildeVec += pref * (logv * Sk + m_b * F * C); } void RedlichKwongMFTP::getPartialMolarCp(span cpbar) const { checkArraySize("RedlichKwongMFTP::getPartialMolarCp", cpbar.size(), m_kk); _updateReferenceStateThermo(); - - double T = temperature(); - double sqt = sqrt(T); - double v = molarVolume(); - double b = m_b_current; - double a = m_a_current; - double dadt = da_dt(); - double d2adt2 = 0.0; // current RK model uses a(T)=a0+a1*T - for (size_t k = 0; k < m_kk; k++) { cpbar[k] = GasConstant * m_cp0_R[k]; } - if (fabs(b) < SmallNumber) { + if (fabs(m_bMix) < SmallNumber) { return; } + double T = temperature(); + double sqt = sqrt(T); + double v = molarVolume(); + double dadt = da_dt(); + double d2adt2 = 0.0; // current RK model uses a(T)=a0+a1*T + // Mixture pressure derivatives pressureDerivatives(); - double pT = dpdT_; - double pV = dpdV_; - double dvdT_P = -pT / pV; - - double vpb = v + b; - double vmb = v - b; + double dvdT_P = -dpdT_ / dpdV_; + double vpb = v + m_bMix; + double vmb = v - m_bMix; // P_TT, P_TV and P_VV for v''(T)|P from EOS: // v'' = -(P_TT + 2 P_TV v' + P_VV v'^2) / P_V - double pTT = (-d2adt2 + dadt / T - 3.0 * a / (4.0 * T * T)) / (sqt * v * vpb); + double pTT = (-d2adt2 + dadt / T - 3.0 * m_aMix / (4.0 * T * T)) / (sqt * v * vpb); double pTV = -GasConstant / (vmb * vmb) - + (dadt - a / (2.0 * T)) * (2.0 * v + b) / (sqt * v * v * vpb * vpb); + + (dadt - m_aMix / (2.0 * T)) * (2.0 * v + m_bMix) / (sqt * v * v * vpb * vpb); double pVV = 2.0 * GasConstant * T / (vmb * vmb * vmb) - - 2.0 * a * (3.0 * v * v + 3.0 * b * v + b * b) + - 2.0 * m_aMix * (3.0 * v * v + 3.0 * m_bMix * v + m_bMix * m_bMix) / (sqt * v * v * v * vpb * vpb * vpb); - double d2vdT2_P = -(pTT + 2.0 * pTV * dvdT_P + pVV * dvdT_P * dvdT_P) / pV; - - // Species-dependent A_k and dA_k/dT - for (size_t k = 0; k < m_kk; k++) { - m_pp[k] = 0.0; - m_dAkdT[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk * i; - m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; - m_dAkdT[k] += moleFractions_[i] * a_coeff_vec(1, counter); - } - } + double d2vdT2_P = -(pTT + 2.0 * pTV * dvdT_P + pVV * dvdT_P * dvdT_P) / dpdV_; - double F = T * dadt - 1.5 * a; + double F = T * dadt - 1.5 * m_aMix; double dF_dT = -0.5 * dadt + T * d2adt2; double logvpv = log(vpb / v); - double termLPrime = -b * dvdT_P / (v * vpb); + double termLPrime = -m_bMix * dvdT_P / (v * vpb); for (size_t k = 0; k < m_kk; k++) { - double bk = b_vec_Curr_[k]; - double Ak = m_pp[k]; + double bk = m_b[k]; + double Ak = m_Ak[k]; double dAk = m_dAkdT[k]; double Sk = 2.0 * T * dAk - 3.0 * Ak; double dSk_dT = -dAk; // for linear a(T) @@ -434,7 +324,7 @@ void RedlichKwongMFTP::getPartialMolarCp(span cpbar) const // dp/dn_k at constant T,V,n_j double dpdni = RT() / vmb + RT() * bk / (vmb * vmb) - 2.0 * Ak / (sqt * v * vpb) - + a * bk / (sqt * v * vpb * vpb); + + m_aMix * bk / (sqt * v * vpb * vpb); // d/dT|P [dp/dn_k] double d_dpdni_dT = GasConstant / vmb @@ -443,17 +333,17 @@ void RedlichKwongMFTP::getPartialMolarCp(span cpbar) const - 2.0 * GasConstant * T * bk * dvdT_P / (vmb * vmb * vmb) - 2.0 * dAk / (sqt * v * vpb) + Ak / (T * sqt * v * vpb) - + 2.0 * Ak * (2.0 * v + b) * dvdT_P / (sqt * v * v * vpb * vpb) - + bk * dadt / (sqt * v * vpb * vpb) - - bk * a / (2.0 * T * sqt * v * vpb * vpb) - - bk * a * (3.0 * v + b) * dvdT_P / (sqt * v * v * vpb * vpb * vpb); + + 2.0 * Ak * (2.0 * v + m_bMix) * dvdT_P / (sqt * v * v * vpb * vpb) + + bk * da_dt() / (sqt * v * vpb * vpb) + - bk * m_aMix / (2.0 * T * sqt * v * vpb * vpb) + - bk * m_aMix * (3.0 * v + m_bMix) * dvdT_P / (sqt * v * v * vpb * vpb * vpb); // d/dT|P of departure terms used in getPartialMolarEnthalpies - double dterm1 = -bk / (b * b * sqt) + double dterm1 = -bk / (m_bMix * m_bMix * sqt) * (termLPrime * F + logvpv * dF_dT - logvpv * F / (2.0 * T)); - double dterm2 = 1.0 / (b * sqt) + double dterm2 = 1.0 / (m_bMix * sqt) * (termLPrime * Sk + logvpv * dSk_dT - logvpv * Sk / (2.0 * T)); - double dterm3 = bk / (b * sqt) + double dterm3 = bk / (m_bMix * sqt) * (dF_dT / vpb - F * dvdT_P / (vpb * vpb) - F / (2.0 * T * vpb)); // cp_k = d hbar_k / dT at constant P and composition @@ -469,100 +359,74 @@ void RedlichKwongMFTP::getPartialMolarCv_TV(span cvtilde) const checkArraySize("RedlichKwongMFTP::getPartialMolarCv_TV", cvtilde.size(), m_kk); _updateReferenceStateThermo(); - double T = temperature(); - double sqt = sqrt(T); - double mv = molarVolume(); - double b = m_b_current; - double a = m_a_current; - double dadt = da_dt(); + Eigen::Map cvtildeVec(cvtilde.data(), m_kk); + cvtildeVec = GasConstant * (asVectorXd(m_cp0_R).array() - 1.0); - for (size_t k = 0; k < m_kk; k++) { - cvtilde[k] = GasConstant * (m_cp0_R[k] - 1.0); - } - - if (fabs(b) < SmallNumber) { + if (fabs(m_bMix) < SmallNumber) { return; } - // A_k = sum_i x_i a_ki and A'_k = dA_k/dT = sum_i x_i a_ki,1 - for (size_t k = 0; k < m_kk; k++) { - m_pp[k] = 0.0; - m_dAkdT[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk * i; - m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; - m_dAkdT[k] += moleFractions_[i] * a_coeff_vec(1, counter); - } - } + double T = temperature(); + double sqt = sqrt(T); + double mv = molarVolume(); + double dadt = da_dt(); // Residual contribution for // \tilde{u}_k = (\partial U / \partial n_k)_{T,V,n_j} // with linear a(T): a_ij = a_ij,0 + a_ij,1 T - double vpb = mv + b; + double vpb = mv + m_bMix; double logv = log(vpb / mv); - double F = T * dadt - 1.5 * a; - double C = -logv / b + 1.0 / vpb; - double pref = 1.0 / (b * sqt); + double F = T * dadt - 1.5 * m_aMix; + double C = -logv / m_bMix + 1.0 / vpb; + double pref = 1.0 / (m_bMix * sqt); - for (size_t k = 0; k < m_kk; k++) { - double Sk = 2.0 * T * m_dAkdT[k] - 3.0 * m_pp[k]; - double ures = pref * (logv * Sk + b_vec_Curr_[k] * F * C); - double dresdT = pref * (-logv * m_dAkdT[k] - 0.5 * b_vec_Curr_[k] * dadt * C) - - 0.5 * ures / T; - cvtilde[k] += dresdT; - } + Eigen::ArrayXd Sk = 2.0 * T * m_dAkdT - 3.0 * m_Ak; + Eigen::ArrayXd ures = pref * (logv * Sk + m_b * F * C); + Eigen::ArrayXd dresdT = - pref * (logv * m_dAkdT + 0.5 * m_b * dadt * C) + - 0.5 * ures / T; + cvtildeVec += dresdT; } void RedlichKwongMFTP::getPartialMolarVolumes(span vbar) const { checkArraySize("RedlichKwongMFTP::getPartialMolarVolumes", vbar.size(), m_kk); - for (size_t k = 0; k < m_kk; k++) { - m_pp[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk*i; - m_pp[k] += moleFractions_[i] * a_vec_Curr_[counter]; - } - } - for (size_t k = 0; k < m_kk; k++) { - m_workS[k] = 0.0; - for (size_t i = 0; i < m_kk; i++) { - size_t counter = k + m_kk*i; - m_workS[k] += moleFractions_[i] * a_coeff_vec(1,counter); - } - } - double sqt = sqrt(temperature()); double mv = molarVolume(); - double vmb = mv - m_b_current; - double vpb = mv + m_b_current; - for (size_t k = 0; k < m_kk; k++) { - double num = (RT() + RT() * m_b_current/ vmb + RT() * b_vec_Curr_[k] / vmb - + RT() * m_b_current * b_vec_Curr_[k] /(vmb * vmb) - - 2.0 * m_pp[k] / (sqt * vpb) - + m_a_current * b_vec_Curr_[k] / (sqt * vpb * vpb) - ); - double denom = (pressure() + RT() * m_b_current/(vmb * vmb) - m_a_current / (sqt * vpb * vpb) - ); - vbar[k] = num / denom; - } + double vmb = mv - m_bMix; + double vpb = mv + m_bMix; + Eigen::ArrayXd num = RT() * (1.0 + (m_bMix + m_b) / vmb + m_bMix * m_b / (vmb * vmb)) + - 2.0 * m_Ak / (sqt * vpb) + + m_aMix * m_b / (sqt * vpb * vpb); + double denom = pressure() + RT() * m_bMix / (vmb * vmb) - m_aMix / (sqt * vpb * vpb); + MappedVector(vbar.data(), m_kk) = num / denom; } bool RedlichKwongMFTP::addSpecies(shared_ptr spec) { bool added = MixtureFugacityTP::addSpecies(spec); if (added) { - a_vec_Curr_.resize(m_kk * m_kk, 0.0); - - // Initialize a_vec and b_vec to NaN, to screen for species with - // pureFluidParameters which are undefined in the input file: - b_vec_Curr_.push_back(NAN); - a_coeff_vec.resize(2, m_kk * m_kk, NAN); - - m_pp.push_back(0.0); - m_dAkdT.push_back(0.0); + size_t k = m_kk - 1; + + m_a.conservativeResize(m_kk, m_kk); + m_a.row(k).setZero(); + m_a.col(k).setZero(); + + // Initialize new entries to NaN to screen for species with undefined + // pure-fluid parameters in the input file. + m_b.conservativeResize(m_kk); + m_b[k] = NAN; + m_a0.conservativeResize(m_kk, m_kk); + m_a0.row(k).setConstant(NAN); + m_a0.col(k).setConstant(NAN); + m_a1.conservativeResize(m_kk, m_kk); + m_a1.row(k).setConstant(NAN); + m_a1.col(k).setConstant(NAN); + + m_Ak.conservativeResizeLike(Eigen::VectorXd::Zero(m_kk)); + m_dAkdT.conservativeResizeLike(Eigen::VectorXd::Zero(m_kk)); m_coeffSource.push_back(CoeffSource::EoS); m_partialMolarVolumes.push_back(0.0); - dpdni_.push_back(0.0); + m_dpdni.conservativeResizeLike(Eigen::VectorXd::Zero(m_kk)); } return added; } @@ -576,7 +440,7 @@ void RedlichKwongMFTP::initThermo() for (auto& [name, species] : m_species) { auto& data = species->input; size_t k = speciesIndex(name, true); - if (!isnan(a_coeff_vec(0, k + m_kk * k))) { + if (!isnan(m_a0(k, k)) && !isnan(m_b[k])) { continue; } bool foundCoeffs = false; @@ -675,21 +539,20 @@ void RedlichKwongMFTP::getSpeciesParameters(const string& name, auto& eosNode = speciesNode["equation-of-state"].getMapWhere( "model", "Redlich-Kwong", true); - size_t counter = k + m_kk * k; - if (a_coeff_vec(1, counter) != 0.0) { + if (m_a1(k, k) != 0.0) { vector coeffs(2); - coeffs[0].setQuantity(a_coeff_vec(0, counter), "Pa*m^6/kmol^2*K^0.5"); - coeffs[1].setQuantity(a_coeff_vec(1, counter), "Pa*m^6/kmol^2/K^0.5"); + coeffs[0].setQuantity(m_a0(k, k), "Pa*m^6/kmol^2*K^0.5"); + coeffs[1].setQuantity(m_a1(k, k), "Pa*m^6/kmol^2/K^0.5"); eosNode["a"] = std::move(coeffs); } else { - eosNode["a"].setQuantity(a_coeff_vec(0, counter), + eosNode["a"].setQuantity(m_a0(k, k), "Pa*m^6/kmol^2*K^0.5"); } - eosNode["b"].setQuantity(b_vec_Curr_[k], "m^3/kmol"); + eosNode["b"].setQuantity(m_b[k], "m^3/kmol"); } else if (m_coeffSource[k] == CoeffSource::CritProps) { auto& critProps = speciesNode["critical-parameters"]; - double a = a_coeff_vec(0, k + m_kk * k); - double b = b_vec_Curr_[k]; + double a = m_a0(k, k); + double b = m_b[k]; double Tc = pow(a * omega_b / (b * omega_a * GasConstant), 2.0/3.0); double Pc = omega_b * GasConstant * Tc / b; critProps["critical-temperature"].setQuantity(Tc, "K"); @@ -716,53 +579,44 @@ void RedlichKwongMFTP::getSpeciesParameters(const string& name, double RedlichKwongMFTP::sresid() const { - if (m_b_current == 0.0) { + if (m_bMix == 0.0) { return 0.0; } // note this agrees with tpx - double rho = density(); - double mmw = meanMolecularWeight(); - double molarV = mmw / rho; - double hh = m_b_current / molarV; + double hh = m_bMix / molarVolume(); double zz = z(); - double dadt = da_dt(); double T = temperature(); - double sqT = sqrt(T); - double fac = dadt - m_a_current / (2.0 * T); - double sresid_mol_R = log(zz*(1.0 - hh)) + log(1.0 + hh) * fac / (sqT * GasConstant * m_b_current); + double fac = da_dt() - m_aMix / (2.0 * T); + double sresid_mol_R = log(zz*(1.0 - hh)) + + log(1.0 + hh) * fac / (sqrt(T) * GasConstant * m_bMix); return GasConstant * sresid_mol_R; } double RedlichKwongMFTP::hresid() const { - if (m_b_current == 0.0) { + if (m_bMix == 0.0) { return 0.0; } // note this agrees with tpx - double rho = density(); - double mmw = meanMolecularWeight(); - double molarV = mmw / rho; - double hh = m_b_current / molarV; + double hh = m_bMix / molarVolume(); double zz = z(); - double dadt = da_dt(); double T = temperature(); - double sqT = sqrt(T); - double fac = T * dadt - 3.0 *m_a_current / (2.0); - return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqT * m_b_current); + double fac = T * da_dt() - 3.0 *m_aMix / (2.0); + return GasConstant * T * (zz - 1.0) + fac * log(1.0 + hh) / (sqrt(T) * m_bMix); } -double RedlichKwongMFTP::liquidVolEst(double TKelvin, double& presGuess) const +double RedlichKwongMFTP::liquidVolEst(double T, double& presGuess) const { - double v = m_b_current * 1.1; + double v = m_bMix * 1.1; double atmp; double btmp; - calculateAB(TKelvin, atmp, btmp); - double pres = std::max(psatEst(TKelvin), presGuess); + calculateAB(T, atmp, btmp); + double pres = std::max(psatEst(T), presGuess); double Vroot[3]; bool foundLiq = false; int m = 0; while (m < 100 && !foundLiq) { - int nsol = solveCubic(TKelvin, pres, atmp, btmp, Vroot); + int nsol = solveCubic(T, pres, atmp, btmp, Vroot); if (nsol == 1 || nsol == 2) { double pc = critPressure(); if (pres > pc) { @@ -783,29 +637,29 @@ double RedlichKwongMFTP::liquidVolEst(double TKelvin, double& presGuess) const return v; } -double RedlichKwongMFTP::densityCalc(double TKelvin, double presPa, int phaseRequested, +double RedlichKwongMFTP::densityCalc(double T, double presPa, int phaseRequested, double rhoguess) { // It's necessary to set the temperature so that m_a_current is set correctly. - setTemperature(TKelvin); + setTemperature(T); double tcrit = critTemperature(); double mmw = meanMolecularWeight(); if (rhoguess == -1.0) { if (phaseRequested >= FLUID_LIQUID_0) { - double lqvol = liquidVolEst(TKelvin, presPa); + double lqvol = liquidVolEst(T, presPa); rhoguess = mmw / lqvol; } } else { // Assume the Gas phase initial guess, if nothing is specified to the routine - rhoguess = presPa * mmw / (GasConstant * TKelvin); + rhoguess = presPa * mmw / (GasConstant * T); } double volguess = mmw / rhoguess; - NSolns_ = solveCubic(TKelvin, presPa, m_a_current, m_b_current, Vroot_); + NSolns_ = solveCubic(T, presPa, m_aMix, m_bMix, Vroot_); // Ensure root ordering used for branch selection only contains physical roots. - double vmin = std::max(0.0, m_b_current * (1.0 + 1e-12)); + double vmin = std::max(0.0, m_bMix * (1.0 + 1e-12)); vector physicalRoots; for (double root : Vroot_) { if (std::isfinite(root) && root > vmin) { @@ -858,7 +712,7 @@ double RedlichKwongMFTP::densityCalc(double TKelvin, double presPa, int phaseReq } else if (NSolns_ == -1) { if (phaseRequested >= FLUID_LIQUID_0 || phaseRequested == FLUID_UNDEFINED || phaseRequested == FLUID_SUPERCRIT) { molarVolLast = Vroot_[0]; - } else if (TKelvin > tcrit) { + } else if (T > tcrit) { molarVolLast = Vroot_[0]; } else { return -2.0; @@ -870,17 +724,16 @@ double RedlichKwongMFTP::densityCalc(double TKelvin, double presPa, int phaseReq return mmw / molarVolLast; } -double RedlichKwongMFTP::dpdVCalc(double TKelvin, double molarVol, double& presCalc) const +double RedlichKwongMFTP::dpdVCalc(double T, double molarVol, double& presCalc) const { - double sqt = sqrt(TKelvin); - presCalc = GasConstant * TKelvin / (molarVol - m_b_current) - - m_a_current / (sqt * molarVol * (molarVol + m_b_current)); - - double vpb = molarVol + m_b_current; - double vmb = molarVol - m_b_current; - double dpdv = (- GasConstant * TKelvin / (vmb * vmb) - + m_a_current * (2 * molarVol + m_b_current) / (sqt * molarVol * molarVol * vpb * vpb)); - return dpdv; + double sqt = sqrt(T); + presCalc = GasConstant * T / (molarVol - m_bMix) + - m_aMix / (sqt * molarVol * (molarVol + m_bMix)); + + double vpb = molarVol + m_bMix; + double vmb = molarVol - m_bMix; + return - GasConstant * T / (vmb * vmb) + + m_aMix * (2 * molarVol + m_bMix) / (sqt * molarVol * molarVol * vpb * vpb); } double RedlichKwongMFTP::isothermalCompressibility() const @@ -909,44 +762,32 @@ double RedlichKwongMFTP::soundSpeed() const void RedlichKwongMFTP::pressureDerivatives() const { - double TKelvin = temperature(); + double T = temperature(); double mv = molarVolume(); double pres; - dpdV_ = dpdVCalc(TKelvin, mv, pres); - double sqt = sqrt(TKelvin); - double vpb = mv + m_b_current; - double vmb = mv - m_b_current; - double dadt = da_dt(); - double fac = dadt - m_a_current/(2.0 * TKelvin); - dpdT_ = (GasConstant / vmb - fac / (sqt * mv * vpb)); + dpdV_ = dpdVCalc(T, mv, pres); + double fac = da_dt() - m_aMix/(2.0 * T); + dpdT_ = (GasConstant / (mv - m_bMix) - fac / (sqrt(T) * mv * (mv + m_bMix))); } void RedlichKwongMFTP::updateMixingExpressions() { - double temp = temperature(); + double T = temperature(); + auto x = asVectorXd(moleFractions_); if (m_formTempParam == 1) { - for (size_t i = 0; i < m_kk; i++) { - for (size_t j = 0; j < m_kk; j++) { - size_t counter = i * m_kk + j; - a_vec_Curr_[counter] = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp; - } - } + m_a = m_a0 + T * m_a1; + m_dAkdT = m_a1 * x; } - m_b_current = 0.0; - m_a_current = 0.0; - for (size_t i = 0; i < m_kk; i++) { - m_b_current += moleFractions_[i] * b_vec_Curr_[i]; - for (size_t j = 0; j < m_kk; j++) { - m_a_current += a_vec_Curr_[i * m_kk + j] * moleFractions_[i] * moleFractions_[j]; - } - } - if (isnan(m_b_current)) { + m_bMix = m_b.matrix().dot(x); + m_Ak = m_a * x; + m_aMix = x.dot(m_Ak.matrix()); + if (isnan(m_bMix)) { // One or more species do not have specified coefficients. fmt::memory_buffer b; for (size_t k = 0; k < m_kk; k++) { - if (isnan(b_vec_Curr_[k])) { + if (isnan(m_b[k])) { if (b.size() > 0) { fmt_append(b, ", {}", speciesName(k)); } else { @@ -959,62 +800,36 @@ void RedlichKwongMFTP::updateMixingExpressions() } } -void RedlichKwongMFTP::calculateAB(double temp, double& aCalc, double& bCalc) const +void RedlichKwongMFTP::calculateAB(double T, double& aCalc, double& bCalc) const { - bCalc = 0.0; - aCalc = 0.0; + auto x = asVectorXd(moleFractions_); + bCalc = m_b.matrix().dot(x); if (m_formTempParam == 1) { - for (size_t i = 0; i < m_kk; i++) { - bCalc += moleFractions_[i] * b_vec_Curr_[i]; - for (size_t j = 0; j < m_kk; j++) { - size_t counter = i * m_kk + j; - double a_vec_Curr = a_coeff_vec(0,counter) + a_coeff_vec(1,counter) * temp; - aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j]; - } - } + aCalc = x.dot((m_a0 + T * m_a1) * x); } else { - for (size_t i = 0; i < m_kk; i++) { - bCalc += moleFractions_[i] * b_vec_Curr_[i]; - for (size_t j = 0; j < m_kk; j++) { - size_t counter = i * m_kk + j; - double a_vec_Curr = a_coeff_vec(0,counter); - aCalc += a_vec_Curr * moleFractions_[i] * moleFractions_[j]; - } - } + aCalc = x.dot(m_a0 * x); } } double RedlichKwongMFTP::da_dt() const { - double dadT = 0.0; - if (m_formTempParam == 1) { - for (size_t i = 0; i < m_kk; i++) { - for (size_t j = 0; j < m_kk; j++) { - size_t counter = i * m_kk + j; - dadT+= a_coeff_vec(1,counter) * moleFractions_[i] * moleFractions_[j]; - } - } + if (m_formTempParam == 0) { + return 0.0; } - return dadT; + auto x = asVectorXd(moleFractions_); + return x.dot(m_a1 * x); } void RedlichKwongMFTP::calcCriticalConditions(double& pc, double& tc, double& vc) const { - double a0 = 0.0; - double aT = 0.0; - for (size_t i = 0; i < m_kk; i++) { - for (size_t j = 0; j