diff --git a/include/cantera/thermo/Phase.h b/include/cantera/thermo/Phase.h index 656311c8e64..878d1088633 100644 --- a/include/cantera/thermo/Phase.h +++ b/include/cantera/thermo/Phase.h @@ -377,7 +377,7 @@ class Phase //! @param t Temperature in kelvin //! @param rho Density (kg/m^3) //! @since New in %Cantera 3.0. - void setState_TD(double t, double rho); + virtual void setState_TD(double t, double rho); //! @} end group set thermo state diff --git a/include/cantera/thermo/PlasmaPhase.h b/include/cantera/thermo/PlasmaPhase.h index 61bb85c9338..a050e9dc982 100644 --- a/include/cantera/thermo/PlasmaPhase.h +++ b/include/cantera/thermo/PlasmaPhase.h @@ -85,9 +85,9 @@ class PlasmaPhase: public IdealGasPhase public: //! Construct and initialize a PlasmaPhase object //! directly from an input file. The constructor initializes the electron - //! energy distribution to be Druyvesteyn distribution (m_x = 2.0). The initial - //! electron energy grid is set to a linear space which starts at 0.01 eV and ends - //! at 1 eV with 1000 points. + //! energy distribution to be a Maxwellian distribution (m_isotropicShapeFactor = 1.0). + // The initial electron energy grid is set to a linear space which starts + // at 0.01 eV and ends at 1 eV with 1000 points. /*! * @param inputFile Name of the input file containing the phase definition * to set up the object. If blank, an empty phase will be @@ -105,19 +105,50 @@ class PlasmaPhase: public IdealGasPhase void initThermo() override; + // ================================================================= // + // ================================================================= // + //! @name Overridden from IdealGasPhase or ThermoPhase + //! @{ + bool addSpecies(shared_ptr spec) override; + virtual void setSolution(std::weak_ptr soln) override; + void getParameters(AnyMap& phaseNode) const override; + void setParameters(const AnyMap& phaseNode, + const AnyMap& rootNode=AnyMap()) override; + //! @} + // ================================================================= // + // ================================================================= // + //! @name Electron Species Information + //! @{ + + //! Electron Species Index + size_t electronSpeciesIndex() const { + return m_electronSpeciesIndex; + } + + //! Electron species name + string electronSpeciesName() const { + return speciesName(m_electronSpeciesIndex); + } + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Electron Energy Distribution Functions + //! @{ + //! Set electron energy levels. - //! @param levels The vector of electron energy levels (eV). + //! @param levels The vector of electron energy levels [eV]. //! Length: #m_nPoints. void setElectronEnergyLevels(span levels); //! Get electron energy levels. - //! @param levels The vector of electron energy levels (eV). Length: #m_nPoints + //! @param levels The vector of electron energy levels [eV]. Length: #m_nPoints void getElectronEnergyLevels(span levels) const { Eigen::Map(levels.data(), levels.size()) = m_electronEnergyLevels; } //! Set discretized electron energy distribution. - //! @param levels The vector of electron energy levels (eV). + //! @param levels The vector of electron energy levels [eV]. //! Length: #m_nPoints. //! @param distrb The vector of electron energy distribution. //! Length: #m_nPoints. @@ -137,25 +168,35 @@ class PlasmaPhase: public IdealGasPhase //! @param x The shape factor void setIsotropicShapeFactor(double x); - //! The shape factor of isotropic electron energy distribution + //! The shape factor of isotropic electron energy distribution. double isotropicShapeFactor() const { return m_isotropicShapeFactor; } - //! Set the internally stored electron temperature of the phase (K). - //! @param Te Electron temperature in Kelvin + //! Electron Temperature [K]. + //! @return The electron temperature of the phase. + double electronTemperature() const override { + return m_electronTemp; + } + //! Set the internally stored electron temperature of the phase [K]. + //! @param Te Electron temperature in Kelvin. void setElectronTemperature(double Te) override; - //! Set mean electron energy [eV]. This method also sets electron temperature - //! accordingly. + //! Mean electron energy [eV]. + double meanElectronEnergy() const { + return 3.0 / 2.0 * electronTemperature() * Boltzmann / ElectronCharge; + } + + //! Set mean electron energy [eV]. + //! This method also sets electron temperature accordingly. void setMeanElectronEnergy(double energy); - //! Get electron energy distribution type + //! Get electron energy distribution type. string electronEnergyDistributionType() const { return m_distributionType; } - //! Set electron energy distribution type + //! Set electron energy distribution type. void setElectronEnergyDistributionType(const string& type); //! Numerical quadrature method. Method: #m_quadratureMethod @@ -169,12 +210,7 @@ class PlasmaPhase: public IdealGasPhase m_quadratureMethod = method; } - //! Mean electron energy [eV] - double meanElectronEnergy() const { - return 3.0 / 2.0 * electronTemperature() * Boltzmann / ElectronCharge; - } - - //! Set flag of automatically normalize electron energy distribution + //! Set flag of automatically normalize electron energy distribution. //! Flag: #m_do_normalizeElectronEnergyDist void enableNormalizeElectronEnergyDist(bool enable) { m_do_normalizeElectronEnergyDist = enable; @@ -186,37 +222,12 @@ class PlasmaPhase: public IdealGasPhase return m_do_normalizeElectronEnergyDist; } - bool addSpecies(shared_ptr spec) override; - - //! Electron Temperature (K) - //! @return The electron temperature of the phase - double electronTemperature() const override { - return m_electronTemp; - } - - //! Return the Gas Constant multiplied by the current electron temperature - /*! - * The units are Joules kmol-1 - */ - double RTe() const { - return electronTemperature() * GasConstant; - } - - /** - * Electron pressure. Units: Pa. - * @f[P = n_{k_e} R T_e @f] - */ - virtual double electronPressure() const { - return GasConstant * concentration(m_electronSpeciesIndex) * - electronTemperature(); - } - - //! Number of electron levels + //! Number of electron levels. size_t nElectronEnergyLevels() const { return m_nPoints; } - //! Number of electron collision cross sections + //! Number of electron collision cross sections. size_t nCollisions() const { return m_collisions.size(); } @@ -234,20 +245,131 @@ class PlasmaPhase: public IdealGasPhase return m_collisionRates[i]; } - //! Electron Species Index - size_t electronSpeciesIndex() const { - return m_electronSpeciesIndex; + //! Update the electron energy distribution. + void updateElectronEnergyDistribution(); + + + //! Return the distribution number #m_distNum. + int distributionNumber() const { + return m_distNum; + } + + //! Return the electron energy level number #m_levelNum. + int levelNumber() const { + return m_levelNum; + } + + //! Get the indicies for inelastic electron collisions. + //! @since New in %Cantera 3.2. + const vector& kInelastic() const { + return m_kInelastic; } + //! Get the indices for elastic electron collisions. + //! @since New in %Cantera 3.2. + const vector& kElastic() const { + return m_kElastic; + } + + //! Return the target of a specific process. + //! @since New in %Cantera 3.2. + size_t targetIndex(size_t i) const { + return m_targetSpeciesIndices[i]; + } + + //! Get the frequency of the applied electric field [Hz]. + //! @since New in %Cantera 3.2. + double electricFieldFrequency() const { + return m_electricFieldFrequency; + } + + //! Get the applied electric field strength [V/m]. + double electricField() const { + return m_electricField; + } + + //! Set the absolute electric field strength [V/m]. + void setElectricField(double E) { + m_electricField = E; + } + + //! Calculate the degree of ionization + //double ionDegree() const { + // double ne = concentration(m_electronSpeciesIndex); // [kmol/m³] + // double n_total = molarDensity(); // [kmol/m³] + // return ne / n_total; + //} + + //! Get the reduced electric field strength [V·m²]. + double reducedElectricField() const { + return m_electricField / (molarDensity() * Avogadro); + } + + //! Set reduced electric field given in [V·m²]. + void setReducedElectricField(double EN) { + m_electricField = EN * molarDensity() * Avogadro; // [V/m] + } + + /** + * The electron mobility (m²/V/s) + * @f[ + * \mu = \nu_d / E, + * @f] + * where @f$ \nu_d @f$ is the drift velocity (m²/s), and @f$ E @f$ is the electric + * field strength (V/m). + */ + double electronMobility() const; + + /** + * The elastic power loss [J/s/m³] + * @f[ + * P_k = N_A N_A C_e e \sum_k C_k K_k, + * @f] + * where @f$ C_k @f$ and @f$ C_e @f$ are the concentration (kmol/m³) of the + * target species and electrons, respectively. @f$ K_k @f$ is the elastic + * electron energy loss coefficient (eV-m³/s). + */ + double elasticPowerLoss(); + + /** + * The joule heating power (W/m³) + * @f[ + * q_J = \sigma * E^2, + * @f] + * where @f$ \sigma @f$ is the conductivity (S/m), defined by: + * @f[ + * \sigma = e * n_e * \mu_e + * @f] + * and @f$ E @f$ is the electric field strength (V/m). + */ + double jouleHeatingPower() const; + + void beginEquilibrate() override; + + void endEquilibrate() override; + + double intrinsicHeating() override; + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Molar Thermodynamic Properties of the Solution + //! @{ + //! Return the Molar enthalpy. Units: J/kmol. /*! - * For an ideal gas mixture with additional electron, + * For an ideal gas mixture with electrons at a different temperature + * than the heavy species, the molar enthalpy is calculated as: * @f[ - * \hat h(T) = \sum_{k \neq k_e} X_k \hat h^0_k(T) + X_{k_e} \hat h^0_{k_e}(T_e), + * \begin{align} + * \hat h(T, T_e) &= \sum_{k} X_k \hat h^0_k(T_k) \\ + * &= \sum_{k \neq k_e} X_k \hat h^0_k(T_g) + X_{k_e} \hat h^0_{k_e}(T_e), + * \end{align} * @f] - * and is a function only of temperature. The standard-state pure-species - * enthalpies @f$ \hat h^0_k(T) @f$ are computed by the species - * thermodynamic property manager. + * where heavy-species properties are evaluated at @f$T@f$, and electron + * properties at @f$T_e@f$. + * The standard-state pure-species enthalpies @f$ \hat h^0_k(T_k) @f$ are + * computed by the species thermodynamic property manager. * * @see MultiSpeciesThermo */ @@ -296,139 +418,329 @@ class PlasmaPhase: public IdealGasPhase */ double intEnergy_mole() const override; - void getEntropy_R(span sr) const override; + // double cp_mole() const override; // Already defined in IdealGasPhase + // double cp_mass() const; // Already defined in ThermoPhase + // double cv_mole() const; // Already defined in IdealGasPhase - void getGibbs_RT(span grt) const override; - void getGibbs_ref(span g) const override; + //! @} + // ================================================================= // + // ================================================================= // + //! @name Mechanical Equation of State + //! @{ - void getStandardVolumes_ref(span vol) const override; + //! Return the mean temperature of the plasma phase. Units: K. + /*! + * In a plasma phase, the electron temperature can differ from the + * heavy-species (gas) temperature. Therefore, the mean temperature is + * defined as a mole-fraction-weighted average of the electron and + * heavy-species temperatures: + * @f[ + * \overline{T} = \sum_{k \neq k_e} X_k T_g + X_{k_e} T_e + * = (1 - X_{k_e}) T_g + X_{k_e} T_e + * = T_g + X_{k_e} (T_e - T_g) + * @f] + * See the `pressure()` method for usage of the mean temperature in + * calculating the pressure of the plasma phase. + */ + double meanTemperature() const; - void getChemPotentials(span mu) const override; + //! Return the pressure of the plasma phase. Units: Pa. + /*! + * The pressure of the plasma phase is calculated using the mean temperature, + * which is a mole-fraction-weighted average of the electron and heavy-species + * temperatures. + * @f[ + * P = \sum_k n_k k_B T_k + * = \sum_{k \neq e} n_k k_B T_g + n_{e} k_B T_e + * = (n_{total} - n_{e}) k_B T_g + n_{e} k_B T_e + * = n_{total} (1 - X_e) k_B T_g + n_{total} X_e k_B T_e + * = n_{total} k_B \overline{T} + * @f] + * where @f$ \overline{T} @f$ is the mean temperature of the plasma phase, + * defined in the `meanTemperature()` method. + * Here, @f$ n_k @f$ is the number density of species @f$ k @f$ [in 1/m³], + * @f$ n_{total} @f$ is the total number density of the phase [in 1/m³], + * @f$ X_e @f$ is the mole fraction of electrons, @f$ T_g @f$ is the + * heavy-species (gas) temperature [K], @f$ T_e @f$ is the electron temperature [K], + * and @f$ k_B @f$ is the Boltzmann constant [J/K]. + * + * The number density times Boltzmann constant can be expressed as + * @f$ n_{total} k_B = C_{total} R @f$, where @f$ C_{total} @f$ is the + * total molar concentration of the phase [in kmol/m³] and + * @f$ R @f$ is the gas constant [J/kmol/K], so that: + * @f[ + * P = C_{total} R \overline{T}. + * @f] + */ + double pressure() const override; - void getStandardChemPotentials(span muStar) const override; + //! Set the pressure at constant temperature and composition. Units: Pa. + /*! + * This method is implemented by setting the mass density to + * @f[ + * \rho = \frac{P \overline{W}}{\hat R \overline{T} }. + * @f] + * + * @param p Pressure [Pa]. + */ + void setPressure(double p) override { + setDensity(p * meanMolecularWeight() / (GasConstant * meanTemperature())); + } - void getPartialMolarEnthalpies(span hbar) const override; + //! Return the electron pressure. Units: Pa. + /* + * The partial pressure of electrons is calculated using the electron temperature: + * @f[ + * P_e = n_e k_B T_e = C_e R T_e, + * @f] + * where @f$ C_e @f$ is the molar concentration of electrons [in kmol/m³], + * @f$ R @f$ is the gas constant [J/kmol/K], and @f$ T_e @f$ is the electron temperature [K]. + */ + virtual double electronPressure() const { + return GasConstant * concentration(m_electronSpeciesIndex) * + electronTemperature(); + } - void getPartialMolarEntropies(span sbar) const override; + //! Return the Gas Constant multiplied by the current electron temperature [J/kmol]. + double RTe() const { + return electronTemperature() * GasConstant; + } - void getPartialMolarIntEnergies(span ubar) const override; + double thermalExpansionCoeff() const override { + // Which temperature to use here? + throw NotImplementedError("PlasmaPhase::thermalExpansionCoeff"); + } - void getParameters(AnyMap& phaseNode) const override; + double soundSpeed() const override { + // Which temperature to use here? + throw NotImplementedError("PlasmaPhase::soundSpeed"); + } - void setParameters(const AnyMap& phaseNode, - const AnyMap& rootNode=AnyMap()) override; + //! @} + // ================================================================= // + // ================================================================= // + //! @name Chemical Potentials and Activities + //! @{ - //! Update the electron energy distribution. - void updateElectronEnergyDistribution(); + //! Returns the standard concentration @f$ C^0_k @f$, which is used to + //! normalize the generalized concentration. Units: m³/kmol. + /*! + * This is defined as the concentration by which the generalized + * concentration is normalized to produce the activity. Since the activity + * for an ideal gas mixture is simply the mole fraction, for an ideal gas + * @f$ C^0_k = P/\hat R T @f$. + * For a multi-temperature system, this translates to: + * @f[ + * \begin{align} + * C^0_k &= \frac{C_k}{a_k} + * &= \frac{X_k \frac{P}{R \overline{T}}}{X_k} + * &= \frac{P}{R \overline{T}} + * \end{align} + * @f] + * where @f$C_k@f$ is the molar concentration of species @f$k@f$ [in kmol/m³], + * @f$ a_k @f$ is the activity of species @f$ k @f$, which is equal to the + * mole fraction @f$ X_k @f$. + * + * @param k Parameter indicating the species. The default + * is to assume this refers to species 0. + */ + double standardConcentration(size_t k=0) const override; - //! Electron species name - string electronSpeciesName() const { - return speciesName(m_electronSpeciesIndex); - } + //! @} + // ================================================================= // + // ================================================================= // + //! @name Partial Molar Properties of the Solution + //! @{ - //! Return the distribution Number #m_distNum - int distributionNumber() const { - return m_distNum; - } + //! Return the chemical potentials of the species in the solution. Units: J/kmol. + /*! + * The chemical potential of species @f$ k @f$ is calculated as: + * @f[ + * \begin{align} + * \mu_k(T_k, X_k, P) + * &= \mu_k^*(T_k, P) + R T_k \ln(X_k) \\ + * &= g_k^0(T_k) + R T_k \ln\left(\frac{P}{P^0}\right) + R T_k \ln(X_k) \\ + * &= h_k^0(T_k) - T_k s_k^0(T_k) + R T_k \ln\left(\frac{P}{P^0}\right) + R T_k \ln(X_k) \\ + * \end{align} + * @f] + */ + void getChemPotentials(span mu) const override; - //! Return the electron energy level Number #m_levelNum - int levelNumber() const { - return m_levelNum; - } + //! Return the partial molar enthalpies of the species in the solution. Units: J/kmol. + /*! + * The partial molar enthalpy of species @f$ k @f$ is @f[h_k^0(T_k)@f] + * where heavy-species properties are evaluated at the gas temperature, + * and electron properties are evaluated at the electron temperature. + * Since @f[h_k^0(T_k)@f] is computed from by @f[m_h0_RT(T_k) * R * T@f], + * we need to update the calculation for electron so that + * @f[h_k^0(T_k) = m_h0_RT(T_k) * R * T_k@f] + */ + void getPartialMolarEnthalpies(span hbar) const override; - //! Get the indicies for inelastic electron collisions - //! @since New in %Cantera 3.2. - const vector& kInelastic() const { - return m_kInelastic; - } + //! Return the partial molar entropies of the species in the solution. Units: J/kmol/K. + /*! + * The partial molar enthalpy of species @f$ k @f$ is @f[s_k^0(T_k) + log(P/P^0)@f] + * where heavy-species properties are evaluated at the gas temperature, + * and electron properties are evaluated at the electron temperature. + * Since @f[s_k^0(T_k)@f] is computed from by @f[m_s0_R(T_k) * R + log(P/P^0)@f], + * it does not depend on temperature, so we don't need to override the function. + * @f[h_k^0(T_k) = m_h0_RT(T_k) * R * T_k@f] + */ + // void getPartialMolarEntropies(span sbar) const override; - //! Get the indices for elastic electron collisions - //! @since New in %Cantera 3.2. - const vector& kElastic() const { - return m_kElastic; - } + //! Return the partial molar internal energies of the species in the solution. Units: J/kmol. + /*! + * The partial molar internal energy of species @f$ k @f$ is calculated as: + * @f[ + * u_k^0 = h_k^0 - R T_k, + * @f] + * where @f$ h_k^0 @f$ is the partial molar enthalpy of species @f$ k @f$, and + * @f$ T_k @f$ is the temperature at which the partial molar enthalpy is evaluated. + * For heavy species, the partial molar enthalpy is evaluated at the gas temperature, + * while for electrons, it is evaluated at the electron temperature. + */ + void getPartialMolarIntEnergies(span ubar) const override; - //! target of a specific process - //! @since New in %Cantera 3.2. - size_t targetIndex(size_t i) const { - return m_targetSpeciesIndices[i]; - } + //! Return the partial molar heat capacities of the species in the solution. Units: J/kmol/K. + /*! + * Since the computation of the partial molar enthalpy does not depend on the temperature, + * there is no need to override the method. + */ + // void getPartialMolarCp(span cpbar) const override; - //! Get the frequency of the applied electric field [Hz] - //! @since New in %Cantera 3.2. - double electricFieldFrequency() const { - return m_electricFieldFrequency; - } + //! Return the partial molar volumes of the species in the solution. Units: m³/kmol. + /*! + * For a multitemperature system, + * @f[ + * v_k = \frac{R T_k}{P}, + * @f] + * where @f$ T_k @f$ is the temperature at which the partial molar enthalpy is evaluated. + */ + void getPartialMolarVolumes(span vbar) const override; - //! Get the applied electric field strength [V/m] - double electricField() const { - return m_electricField; - } + //! @} + // ================================================================= // + // ================================================================= // + //! @name Properties of the Standard State of the Species in the Solution + //! @{ - //! Set the absolute electric field strength [V/m] - void setElectricField(double E) { - m_electricField = E; - } + //! Return the standard chemical potentials of the species. Units: J/kmol. + /*! + * For heavy species, this is identical to the IdealGasPhase + * implementation. For electrons, the standard chemical potential + * is evaluated at the electron temperature: + * @f[ + * \mu^*_{e}(T_e) = g_k^0(T_e) + RT_e \ln \left(\frac{P}{P^0}\right). + * @f] + */ + void getStandardChemPotentials(span muStar) const override; - //! Calculate the degree of ionization - //double ionDegree() const { - // double ne = concentration(m_electronSpeciesIndex); // [kmol/m³] - // double n_total = molarDensity(); // [kmol/m³] - // return ne / n_total; - //} + // Below are 5 methods already defined in IdealGasPhase, that do not need + // to be overridden since `updateThermo` already updates the standard-state + // properties for electrons at the electron temperature. + // And for entropy (and thus Gibbs free energy), for all species, the total + // pressure is used in the logarithmic term (not the partial pressure). + // { + // void getEnthalpy_RT(span hrt) const override; + // void getEntropy_R(span sr) const override; + // void getGibbs_RT(span grt) const override; + // void getIntEnergy_RT(span urt) const override; + // void getCp_R(span cpr) const override; + // } + + //! Return the standard molar volumes of the species. Units: m³/kmol. + /*! + * For a multitemperature system, + * @f[ + * v_k = \frac{R T_k}{P}, + * @f] + * where @f$ T_k @f$ is the temperature at which the standard molar volume is evaluated. + */ + void getStandardVolumes(span vol) const override; + + //! @} + // ================================================================= // + // ================================================================= // + //! @name Thermodynamic Values for the Species Reference States + //! @{ + + // Below are 5 methods already defined in IdealGasPhase, that do not need + // to be overridden, since they do not depend on the temperature. + // { + // void getEnthalpy_RT_ref(span hrt) const override; + // void getGibbs_RT_ref(span grt) const override; + // void getEntropy_R_ref(span er) const override; + // void getIntEnergy_RT_ref(span urt) const override; + // void getCp_R_ref(span cprt) const override; + // } - //! Get the reduced electric field strength [V·m²] - double reducedElectricField() const { - return m_electricField / (molarDensity() * Avogadro); - } + void getGibbs_ref(span g) const override; + void getStandardVolumes_ref(span vol) const override; - //! Set reduced electric field given in [V·m²] - void setReducedElectricField(double EN) { - m_electricField = EN * molarDensity() * Avogadro; // [V/m] - } + //! @} + // ================================================================= // + // ================================================================= // + //! @name Setting the State + //! + //! For a plasma phase, setting the state requires specifying both + //! the heavy-species (gas) temperature and the electron temperature. + //! @{ + + //! Set the state using an AnyMap containing any combination of properties + //! supported by the thermodynamic model + /*! + * Accepted keys are: + * * `X` (mole fractions) + * * `Y` (mass fractions) + * * `T` or `Tg` or `gas-temperature` [K] + * * `Te` or `electron-temperature` [K] + * * `P` or `pressure` [Pa] + * * `D` or `density` [kg/m^3] + * + * Composition can be specified as either an AnyMap of species names to + * values or as a composition string. All other values can be given as + * floating point values in Cantera's default units, or as strings with the + * units specified, which will be converted using the Units class. + */ + void setState(const AnyMap& state) override; - virtual void setSolution(std::weak_ptr soln) override; + //! Set the gas and electron temperature [K] and pressure [Pa]. + /*! + * @param t Temperature [K] + * @param p Pressure [Pa] + */ + void setState_TP(double t, double p) override; - /** - * The elastic power loss (J/s/m³) - * @f[ - * P_k = N_A N_A C_e e \sum_k C_k K_k, - * @f] - * where @f$ C_k @f$ and @f$ C_e @f$ are the concentration (kmol/m³) of the - * target species and electrons, respectively. @f$ K_k @f$ is the elastic - * electron energy loss coefficient (eV-m³/s). + //! Set the gas temperature [K], electron temperature [K], and pressure [Pa]. + /*! + * @param Tg Gas (heavy-species) temperature [K] + * @param Te Electron temperature [K] + * @param p Pressure [Pa] */ - double elasticPowerLoss(); + virtual void setState_TgTeP(double Tg, double Te, double p); - /** - * The electron mobility (m²/V/s) - * @f[ - * \mu = \nu_d / E, - * @f] - * where @f$ \nu_d @f$ is the drift velocity (m²/s), and @f$ E @f$ is the electric - * field strength (V/m). + //! Set the gas and electron temperature [K] and mass density [kg/m^3]. + /*! + * @param t Temperature [K] + * @param rho Mass density [kg/m^3] */ - double electronMobility() const; + void setState_TD(double t, double rho) override; - /** - * The joule heating power (W/m³) - * @f[ - * q_J = \sigma * E^2, - * @f] - * where @f$ \sigma @f$ is the conductivity (S/m), defined by: - * @f[ - * \sigma = e * n_e * \mu_e - * @f] - * and @f$ E @f$ is the electric field strength (V/m). + //! Set the gas temperature [K], electron temperature [K], and mass density [kg/m^3]. + /*! + * @param Tg Gas (heavy-species) temperature [K] + * @param Te Electron temperature [K] + * @param rho Mass density [kg/m^3] */ - double jouleHeatingPower() const; + virtual void setState_TgTeD(double Tg, double Te, double rho); + + //! @} + + - void beginEquilibrate() override; - void endEquilibrate() override; - double intrinsicHeating() override; protected: void updateThermo() const override; @@ -492,31 +804,31 @@ class PlasmaPhase: public IdealGasPhase //! Length: #m_nPoints Eigen::ArrayXd m_electronEnergyDist; - //! Index of electron species + //! Index of electron species. size_t m_electronSpeciesIndex = npos; - //! Electron temperature [K] + //! Electron temperature [K]. double m_electronTemp; - //! Electron energy distribution type + //! Electron energy distribution type. Can be "isotropic", "discretized" or "Boltzmann-two-term". string m_distributionType = "isotropic"; - //! Numerical quadrature method for electron energy distribution + //! Numerical quadrature method for electron energy distribution. string m_quadratureMethod = "simpson"; - //! Flag of normalizing electron energy distribution + //! Flag of normalizing electron energy distribution. bool m_do_normalizeElectronEnergyDist = true; - //! Indices of inelastic collisions in m_crossSections + //! Indices of inelastic collisions in m_crossSections. vector m_kInelastic; - //! Indices of elastic collisions in m_crossSections + //! Indices of elastic collisions in m_crossSections. vector m_kElastic; - //! electric field [V/m] + //! electric field [V/m]. double m_electricField = 0.0; - //! electric field freq [Hz] + //! electric field freq [Hz]. double m_electricFieldFrequency = 0.0; //! Cross section data. m_crossSections[i][j], where i is the specific process, diff --git a/interfaces/cython/cantera/composite.pyi b/interfaces/cython/cantera/composite.pyi index 646fddc254f..24644ca0871 100644 --- a/interfaces/cython/cantera/composite.pyi +++ b/interfaces/cython/cantera/composite.pyi @@ -743,6 +743,8 @@ class Quantity: def electron_species_name(self) -> str: ... @property def elastic_power_loss(self) -> float: ... + @property + def mean_temperature(self) -> float: ... class SolutionArray(SolutionArrayBase, Generic[_P]): _phase: _P | None diff --git a/interfaces/cython/cantera/thermo.pxd b/interfaces/cython/cantera/thermo.pxd index c80bed26316..ce1e5cd8e49 100644 --- a/interfaces/cython/cantera/thermo.pxd +++ b/interfaces/cython/cantera/thermo.pxd @@ -213,6 +213,7 @@ cdef extern from "cantera/thermo/PlasmaPhase.h": double electricField() void updateElectronEnergyDistribution() except +translate_exception double elasticPowerLoss() except +translate_exception + double meanTemperature() cdef extern from "cantera/cython/thermo_utils.h": diff --git a/interfaces/cython/cantera/thermo.pyi b/interfaces/cython/cantera/thermo.pyi index 80a66e547d9..322b8e86514 100644 --- a/interfaces/cython/cantera/thermo.pyi +++ b/interfaces/cython/cantera/thermo.pyi @@ -521,6 +521,8 @@ class ThermoPhase(_SolutionBase): def electron_species_name(self) -> str: ... @property def elastic_power_loss(self) -> float: ... + @property + def mean_temperature(self) -> float: ... class InterfacePhase(ThermoPhase): @property diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index 021f80ca359..9c00e6ba6f5 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -1,2324 +1,2333 @@ -# This file is part of Cantera. See License.txt in the top-level directory or -# at https://cantera.org/license.txt for license and copyright information. - -import warnings -import weakref as _weakref -import numbers as _numbers -import numpy as np -cimport numpy as np - -from .speciesthermo cimport * -from .ctcxx cimport span -from .kinetics cimport CxxKinetics -from .transport cimport * -from ._utils cimport * -from ._utils import CanteraError -from .units cimport * - -cdef enum ThermoBasisType: - mass_basis = 0 - molar_basis = 1 - -ctypedef CxxPlasmaPhase* CxxPlasmaPhasePtr -ctypedef CxxSurfPhase* CxxSurfPhasePtr - -class ThermoModelMethodError(Exception): - """Exception raised for an invalid method used by a thermo model - - :param thermo_model: - The thermo model used by class `ThermoPhase` - - """ - - def __init__(self, thermo_model): - self.thermo_model = thermo_model - super().__init__(f"This method is invalid for {self.thermo_model}") - - -cdef class Species: - """ - A class which stores data about a single chemical species that may be - needed to add it to a `Solution` or `Interface` object (and to the - underlying `ThermoPhase` and `Transport` objects). - - :param name: - A string giving the name of the species, such as ``'CH4'`` - :param composition: - The elemental composition of the species, given either as a dict or a - composition string, such as ``{'C':1, 'H':4}`` or ``'C:1, H:4'``. - :param charge: - The electrical charge, in units of the elementary charge. Default 0.0. - :param size: - The effective size [m] of the species. Default 1.0. - :param init: - Used internally when wrapping :ct:`Species` objects returned from C++ - - Example: creating an ideal gas phase with a single species:: - - ch4 = ct.Species('CH4', 'C:1, H:4') - ch4.thermo = ct.ConstantCp(300, 1000, 101325, - (300, -7.453347e7, 1.865912e5, 3.576053e4)) - tran = ct.GasTransportData() - tran.set_customary_units('nonlinear', 3.75, 141.40, 0.0, 2.60, 13.00) - ch4.transport = tran - gas = ct.Solution(thermo='ideal-gas', species=[ch4]) - - The static methods `list_from_file` and `list_from_yaml` can be used to create - `Species` objects from existing definitions in the YAML format. Either of the - following will produce a list of 53 `Species` objects containing the species defined - in the GRI 3.0 mechanism:: - - S = ct.Species.list_from_file("gri30.yaml") - - import pathlib - S = ct.Species.list_from_yaml( - pathlib.Path('path/to/gri30.yaml').read_text(), - section='species') - - """ - def __cinit__(self, *args, init=True, **kwargs): - if init: - self._species.reset(new CxxSpecies()) - self.species = self._species.get() - - def __init__(self, name=None, composition=None, charge=None, size=None, - *args, init=True, **kwargs): - if not init: - return - - if name is not None: - self.species.name = stringify(name) - - if composition is not None: - self.species.composition = comp_map(composition) - - if charge is not None: - self.species.charge = charge - - if size is not None: - self.species.size = size - - cdef _assign(self, shared_ptr[CxxSpecies] other): - self._species = other - self.species = self._species.get() - - @staticmethod - def from_yaml(text): - """ - Create a `Species` object from its YAML string representation. - """ - cxx_species = CxxNewSpecies(AnyMapFromYamlString(stringify(text))) - species = Species(init=False) - species._assign(cxx_species) - return species - - @staticmethod - def from_dict(data): - """ - Create a `Species` object from a dictionary corresponding to its YAML - representation. - - :param data: - A dictionary corresponding to the YAML representation. - """ - cdef CxxAnyMap any_map = py_to_anymap(data) - cxx_species = CxxNewSpecies(any_map) - species = Species(init=False) - species._assign(cxx_species) - return species - - @staticmethod - def list_from_file(filename, section="species"): - """ - Create a list of Species objects from all of the species defined in the section - ``section`` of a YAML file. Directories on Cantera's input file path will be - searched for the specified file. - """ - root = AnyMapFromYamlFile(stringify(filename)) - species = [] - for a in CxxGetSpecies(root[stringify(section)]): - b = Species(init=False) - b._assign(a) - species.append(b) - return species - - @staticmethod - def list_from_yaml(text, section=None): - """ - Create a list of Species objects from all the species defined in a YAML - string. If ``text`` is a YAML mapping, the ``section`` name of the list - to be read must be specified. If ``text`` is a YAML list, no ``section`` - name should be supplied. - """ - root = AnyMapFromYamlString(stringify(text)) - - # ``items`` is the pseudo-key used to access a list when it is at the - # top level of a YAML document - cxx_species = CxxGetSpecies(root[stringify(section or "items")]) - species = [] - for a in cxx_species: - b = Species(init=False) - b._assign(a) - species.append(b) - return species - - property name: - """ The name of the species. """ - def __get__(self): - return pystr(self.species.name) - - property composition: - """ - A dict containing the elemental composition of the species. Keys are - element names; values are the corresponding atomicities. - """ - def __get__(self): - return comp_map_to_dict(self.species.composition) - - property charge: - """ - The electrical charge on the species, in units of the elementary charge. - """ - def __get__(self): - return self.species.charge - - property size: - """ The effective size [m] of the species. """ - def __get__(self): - return self.species.size - - property molecular_weight: - """The molecular weight [amu] of the species. - - .. versionadded:: 3.0 - """ - def __get__(self): - return self.species.molecularWeight() - - property thermo: - """ - Get/Set the species reference-state thermodynamic data, as an instance - of class `SpeciesThermo`. - """ - def __get__(self): - if self.species.thermo.get() != NULL: - return wrapSpeciesThermo(self.species.thermo) - else: - return None - - def __set__(self, SpeciesThermo spthermo): - self.species.thermo = spthermo._spthermo - - property transport: - """ - Get/Set the species transport parameters, as an instance of class - `GasTransportData`. - """ - def __get__(self): - if self.species.transport.get() != NULL: - data = GasTransportData(init=False) - data._assign(self.species.transport) - return data - else: - return None - def __set__(self, GasTransportData tran): - self.species.transport = tran._data - - property input_data: - """ - Get input data defining this Species, along with any user-specified data - provided with its input (YAML) definition. - """ - def __get__(self): - cdef CxxThermoPhase* phase = self._phase.thermo if self._phase else NULL - return anymap_to_py(self.species.parameters(phase)) - - def update_user_data(self, data): - """ - Add the contents of the provided `dict` as additional fields when generating - YAML phase definition files with `Solution.write_yaml` or in the data returned - by `input_data`. Existing keys with matching names are overwritten. - """ - self.species.input.update(py_to_anymap(data), False) - - def clear_user_data(self): - """ - Clear all saved input data, so that the data given by `input_data` or - `Solution.write_yaml` will only include values generated by Cantera based on - the current object state. - """ - self.species.input.clear() - - def __repr__(self): - return ''.format(self.name) - - -cdef class ThermoPhase(_SolutionBase): - """ - A phase with an equation of state. - - Class `ThermoPhase` may be used to represent the intensive thermodynamic - state of a phase of matter, which might be a gas, liquid, or solid. - - Class `ThermoPhase` is not usually instantiated directly. It is used - as a base class for classes `Solution` and `Interface`. - """ - - def __init__(self, *args, **kwargs): - super().__init__(*args, **kwargs) - if 'source' not in kwargs: - self.thermo_basis = mass_basis - # In composite objects, the ThermoPhase constructor needs to be called first - # to prevent instantiation of stand-alone 'Kinetics' or 'Transport' objects. - # The following is used as a sentinel. After initialization, the _references - # object is used to track whether the ThermoPhase is being used by other - # objects that require the number of species to remain constant and do not - # have C++ implementations (e.g. Quantity objects). - self._references = _weakref.WeakKeyDictionary() - # validate plasma phase - self._enable_plasma = False - if dynamic_cast[CxxPlasmaPhasePtr](self.thermo): - self._enable_plasma = True - self.plasma = self.thermo - - property thermo_model: - """ - Return thermodynamic model describing phase. - """ - def __get__(self): - return pystr(self.thermo.type()) - - property phase_of_matter: - """ - Get the thermodynamic phase (gas, liquid, etc.) at the current conditions. - """ - def __get__(self): - return pystr(self.thermo.phaseOfMatter()) - - def report(self, show_thermo=True, float threshold=1e-14): - """ - Generate a report describing the thermodynamic state of this phase. To - print the report to the terminal, simply call the phase object. The - following two statements are equivalent:: - - >>> phase() - >>> print(phase.report()) - - :param show_thermo: - A Boolean argument specifying whether to show phase thermodynamic - information in the output. - :param threshold: - The threshold used to clip data in the output. Values below the threshold - are not displayed. - """ - return pystr(self.thermo.report(bool(show_thermo), threshold)) - - def __call__(self, *args, **kwargs): - print(self.report(*args, **kwargs)) - - property is_pure: - """ - Returns true if the phase represents a pure (fixed composition) substance - """ - def __get__(self): - return self.thermo.isPure() - - property has_phase_transition: - """ - Returns true if the phase represents a substance with phase transitions - """ - def __get__(self): - return self.thermo.hasPhaseTransition() - - property is_compressible: - """ - Returns true if the density of the phase is an independent variable defining - the thermodynamic state of a substance - """ - def __get__(self): - return self.thermo.isCompressible() - - @property - def _native_mode(self): - """ Return string acronym representing native state """ - return pystr(self.thermo.nativeMode()) - - property _native_state: - """ - Default properties defining a state - """ - def __get__(self): - cdef pair[string, size_t] item - native = {pystr(item.first): item.second for item in self.thermo.nativeState()} - return tuple([i for i, j in sorted(native.items(), key=lambda kv: kv[1])]) - - property _full_states: - """ - Sets of parameters which set the full thermodynamic state - """ - def __get__(self): - states = self.thermo.fullStates() - states = [pystr(s) for s in states] - return {frozenset(k): k for k in states} - - property _partial_states: - """ - Sets of parameters which set a valid partial thermodynamic state - """ - def __get__(self): - states = self.thermo.partialStates() - states = [pystr(s) for s in states] - return {frozenset(k): k for k in states} - - property basis: - """ - Determines whether intensive thermodynamic properties are treated on a - ``mass`` (per kg) or ``molar`` (per kmol) basis. This affects the values - returned by the properties `h`, `u`, `s`, `g`, `v`, `density`, `cv`, - and `cp`, as well as the values used with the state-setting properties - such as `HPX` and `UV`. - """ - def __get__(self): - if self.thermo_basis == mass_basis: - return 'mass' - else: - return 'molar' - - def __set__(self, value): - if value == 'mass': - self.thermo_basis = mass_basis - elif value == 'molar': - self.thermo_basis = molar_basis - else: - raise ValueError("Valid choices are 'mass' or 'molar'." - " Got {!r}.".format(value)) - - cdef double _mass_factor(self): - """ Conversion factor from current basis to kg """ - if self.thermo_basis == molar_basis: - return self.thermo.meanMolecularWeight() - else: - return 1.0 - - cdef double _mole_factor(self): - """ Conversion factor from current basis to moles """ - if self.thermo_basis == mass_basis: - return 1.0/self.thermo.meanMolecularWeight() - else: - return 1.0 - - def equilibrate(self, XY, solver='auto', double rtol=1e-9, - int max_steps=1000, int max_iter=100, int estimate_equil=0, - int log_level=0): - """ - Set to a state of chemical equilibrium holding property pair - ``XY`` constant. - - :param XY: - A two-letter string, which must be one of the set:: - - ['TP','TV','HP','SP','SV','UV'] - - :param solver: - Specifies the equilibrium solver to use. May be one of the following: - - * ``'element_potential'`` - a fast solver using the element potential - method - * ``'gibbs'`` - a slower but more robust Gibbs minimization solver - * ``'vcs'`` - the VCS non-ideal equilibrium solver - * ``'auto'`` - The element potential solver will be tried first, then - if it fails the Gibbs solver will be tried. - :param rtol: - The relative error tolerance. - :param max_steps: - The maximum number of steps in composition to take to find a converged - solution. - :param max_iter: - For the Gibbs minimization solver, this specifies the number of - outer iterations on T or P when some property pair other - than TP is specified. - :param estimate_equil: - Integer indicating whether the solver should estimate its own - initial condition. If 0, the initial mole fraction vector in the - `ThermoPhase` object is used as the initial condition. If 1, the - initial mole fraction vector is used if the element abundances are - satisfied. If -1, the initial mole fraction vector is thrown out, - and an estimate is formulated. - :param log_level: - Set to a value greater than 0 to write diagnostic output. - """ - self.thermo.equilibrate(stringify(XY.upper()), stringify(solver), rtol, - max_steps, max_iter, estimate_equil, log_level) - - ####### Composition, species, and elements ######## - - property n_elements: - """Number of elements.""" - def __get__(self): - return self.thermo.nElements() - - cpdef int element_index(self, element) except *: - """ - The index of element ``element``, which may be specified as a string or - an integer. In the latter case, the index is checked for validity and - returned. If no such element is present, an exception is thrown. - """ - if isinstance(element, (str, bytes)): - return self.thermo.elementIndex(stringify(element), True) - if isinstance(element, (int, float)): - return self.thermo.checkElementIndex(element) - - raise TypeError("'element' must be a string or a number. " - f"Got {element!r}.") - - def element_name(self, m): - """Name of the element with index ``m``.""" - return pystr(self.thermo.elementName(m)) - - property element_names: - """A list of all the element names.""" - def __get__(self): - return [self.element_name(m) for m in range(self.n_elements)] - - def atomic_weight(self, m): - """Atomic weight [kg/kmol] of element ``m``""" - return self.thermo.atomicWeight(self.element_index(m)) - - property atomic_weights: - """Array of atomic weight [kg/kmol] for each element in the mixture.""" - def __get__(self): - return np.array([self.thermo.atomicWeight(m) for m in range(self.n_elements)]) - - property n_species: - """Number of species.""" - def __get__(self): - return self.thermo.nSpecies() - - property n_selected_species: - """ - Number of species selected for output (by slicing of Solution object) - """ - def __get__(self): - return self._selected_species.size or self.n_species - - def species_name(self, k): - """Name of the species with index ``k``.""" - return pystr(self.thermo.speciesName(k)) - - property species_names: - """A list of all the species names.""" - def __get__(self): - if self._selected_species.size: - indices = self._selected_species - else: - indices = range(self.n_species) - return [self.species_name(k) for k in indices] - - cpdef int species_index(self, species) except *: - """ - The index of species ``species``, which may be specified as a string or - an integer. In the latter case, the index is checked for validity and - returned. If no such species is present, an exception is thrown. - """ - if isinstance(species, (str, bytes)): - return self.thermo.speciesIndex(stringify(species), True) - if isinstance(species, (int, float)): - return self.thermo.checkSpeciesIndex(species) - raise TypeError("'species' must be a string or a number." - f" Got {species!r}.") - - property case_sensitive_species_names: - """Enforce case-sensitivity for look up of species names""" - def __get__(self): - return self.thermo.caseSensitiveSpecies() - def __set__(self, val): - self.thermo.setCaseSensitiveSpecies(bool(val)) - - def species(self, k=None): - """ - Return the `Species` object for species ``k``, where ``k`` is either the - species index or the species name. If ``k`` is not specified, a list of - all species objects is returned. Changes to this object do not affect - the `ThermoPhase` or `Solution` object until the `modify_species` - function is called. - """ - if k is None: - return [self.species(i) for i in range(self.n_species)] - - s = Species(init=False) - s._phase = self - if isinstance(k, (str, bytes)): - s._assign(self.thermo.species(stringify(k))) - elif isinstance(k, (int, float)): - s._assign(self.thermo.species(k)) - else: - raise TypeError("Argument must be a string or a number." - " Got {!r}.".format(k)) - return s - - def modify_species(self, k, Species species): - """ - Modify the thermodynamic data associated with a species. The species name, - elemental composition, and type of thermo parameterization must be unchanged. - """ - self.thermo.modifySpecies(k, species._species) - if self.kinetics: - self.kinetics.invalidateCache() - - def add_species(self, Species species): - """ - Add a new species to this phase. Missing elements will be added - automatically. - """ - if self._references: - raise CanteraError('Cannot add species to ThermoPhase object because it' - ' is being used by a Quantity object.') - self.thermo.addUndefinedElements() - self.thermo.addSpecies(species._species) - species._phase = self - self.thermo.initThermo() - if self.kinetics: - self.kinetics.invalidateCache() - - def add_species_alias(self, name, alias): - """ - Add the alternate species name ``alias`` for an original species ``name``. - """ - self.thermo.addSpeciesAlias(stringify(name), stringify(alias)) - - def find_isomers(self, comp): - """ - Find species/isomers matching a composition specified by ``comp``. - """ - - if isinstance(comp, dict): - iso = self.thermo.findIsomers(comp_map(comp)) - elif isinstance(comp, (str, bytes)): - iso = self.thermo.findIsomers(stringify(comp)) - else: - raise CanteraError('Invalid composition') - - return [pystr(b) for b in iso] - - def n_atoms(self, species, element): - """ - Number of atoms of element ``element`` in species ``species``. The element - and species may be specified by name or by index. - - >>> phase.n_atoms('CH4','H') - 4 - """ - return self.thermo.nAtoms(self.species_index(species), - self.element_index(element)) - - cdef np.ndarray _getArray1(self, thermoMethod1d method): - cdef np.ndarray[np.double_t, ndim=1] data = np.empty(self.n_species) - method(self.thermo, span[double](&data[0], data.size)) - if self._selected_species.size: - return data[self._selected_species] - else: - return data - - cdef void _setArray1(self, thermoMethod1d method, values) except *: - cdef np.ndarray[np.double_t, ndim=1] data - cdef span[double] view - - values = np.squeeze(values) - if values.ndim == 0: - values = values[np.newaxis] # corner case for single-species phases - - if len(values) == self.n_species: - data = np.ascontiguousarray(values, dtype=np.double) - elif len(values) == len(self._selected_species) != 0: - data = np.zeros(self.n_species, dtype=np.double) - for i,k in enumerate(self._selected_species): - data[k] = values[i] - else: - msg = "Got {}. Expected {}".format(len(values), self.n_species) - if len(self._selected_species): - msg += ' or {}'.format(len(self._selected_species)) - raise ValueError('Array has incorrect length. ' + msg + '.') - method(self.thermo, span[double](&data[0], data.size)) - - property molecular_weights: - """Array of species molecular weights (molar masses) [kg/kmol].""" - def __get__(self): - return self._getArray1(thermo_getMolecularWeights) - - property charges: - """Array of species charges [elem. charge].""" - def __get__(self): - return self._getArray1(thermo_getCharges) - - property mean_molecular_weight: - """The mean molecular weight (molar mass) [kg/kmol].""" - def __get__(self): - return self.thermo.meanMolecularWeight() - - property Y: - """ - Get/Set the species mass fractions. Can be set as an array, as a dictionary, - or as a string. Always returns an array:: - - >>> phase.Y = [0.1, 0, 0, 0.4, 0, 0, 0, 0, 0.5] - >>> phase.Y = {'H2':0.1, 'O2':0.4, 'AR':0.5} - >>> phase.Y = 'H2:0.1, O2:0.4, AR:0.5' - >>> phase.Y - array([0.1, 0, 0, 0.4, 0, 0, 0, 0, 0.5]) - """ - def __get__(self): - return self._getArray1(thermo_getMassFractions) - def __set__(self, Y): - if isinstance(Y, (str, bytes)): - self.thermo.setMassFractionsByName(stringify(Y)) - elif isinstance(Y, dict): - self.thermo.setMassFractionsByName(comp_map(Y)) - else: - self._setArray1(thermo_setMassFractions, Y) - - property X: - """ - Get/Set the species mole fractions. Can be set as an array, as a dictionary, - or as a string. Always returns an array:: - - >>> phase.X = [0.1, 0, 0, 0.4, 0, 0, 0, 0, 0.5] - >>> phase.X = {'H2':0.1, 'O2':0.4, 'AR':0.5} - >>> phase.X = 'H2:0.1, O2:0.4, AR:0.5' - >>> phase.X - array([0.1, 0, 0, 0.4, 0, 0, 0, 0, 0.5]) - """ - def __get__(self): - return self._getArray1(thermo_getMoleFractions) - def __set__(self, X): - if isinstance(X, (str, bytes)): - self.thermo.setMoleFractionsByName(stringify(X)) - elif isinstance(X, dict): - self.thermo.setMoleFractionsByName(comp_map(X)) - else: - self._setArray1(thermo_setMoleFractions, X) - - property concentrations: - """ - Get/Set the species concentrations. Units are kmol/m³ for bulk phases, kmol/m² - for surface phases, and kmol/m for edge phases. - """ - def __get__(self): - return self._getArray1(thermo_getConcentrations) - def __set__(self, C): - self._setArray1(thermo_setConcentrations, C) - - def __composition_to_array(self, comp, basis): - """take a mixture composition in mole or mass fraction as string, - dict or array and return array (for internal use)""" - if (isinstance(comp, str) and ':' not in comp - and comp in self.species_names): - comp += ':1.0' - - original_state = self.state - - if basis == 'mole': - self.TPX = None, None, comp - arr = np.copy(self.X) - elif basis == 'mass': - self.TPY = None, None, comp - arr = np.copy(self.Y) - else: - raise ValueError("basis must either be 'mass' or mole'.") - - self.state = original_state - return arr - - def set_equivalence_ratio(self, phi, fuel, oxidizer, basis="mole", *, diluent=None, - fraction=None): - """ - Set the composition to a mixture of ``fuel`` and ``oxidizer`` at the - specified equivalence ratio ``phi``, holding temperature and pressure - constant. Considers the oxidation of C to CO2, H to H2O and S to SO2. - Other elements are assumed not to participate in oxidation (that is, - N ends up as N2). The ``basis`` determines the fuel and oxidizer - compositions: ``basis='mole'`` means mole fractions (default), - ``basis='mass'`` means mass fractions. The fuel/oxidizer mixture can be - be diluted by a ``diluent`` based on a mixing ``fraction``. The amount of - diluent is quantified as a fraction of fuel, oxidizer or the fuel/oxidizer - mixture. For more information, see the :doc:`Python example - ` :: - - >>> gas.set_equivalence_ratio(0.5, 'CH4', 'O2:1.0, N2:3.76', basis='mole') - >>> gas.mass_fraction_dict() - {'CH4': 0.02837633052851, 'N2': 0.7452356312613, 'O2': 0.22638803821018} - >>> gas.set_equivalence_ratio(1.2, 'NH3:0.8,CO:0.2', 'O2:1', basis='mole') - >>> gas.mass_fraction_dict() - {'CO': 0.14784006249290, 'NH3': 0.35956645545401, 'O2': 0.49259348205308} - - :param phi: - Equivalence ratio - :param fuel: - Fuel species name or mole/mass fractions as string, array, or dict. - :param oxidizer: - Oxidizer species name or mole/mass fractions as a string, array, or dict. - :param basis: - Determines if ``fuel`` and ``oxidizer`` are given in mole - fractions (``basis='mole'``) or mass fractions (``basis='mass'``). - :param diluent: - Optional parameter. Required if dilution is used. Specifies the composition - of the diluent in mole/mass fractions as a string, array or dict. - :param fraction: - Optional parameter. Dilutes the fuel/oxidizer mixture with the diluent - according to ``fraction``. Fraction can refer to the fraction of diluent in - the mixture (for example ``fraction="diluent:0.7"`` will create a mixture - with 30 % fuel/oxidizer and 70 % diluent), the fraction of fuel in the - mixture (for example ``fraction="fuel:0.1"`` means that the mixture contains - 10 % fuel. The amount of oxidizer is determined from the equivalence ratio - and the remaining mixture is the diluent) or fraction of oxidizer in the - mixture (for example ``fraction="oxidizer:0.1"``). The fraction itself is - interpreted as mole or mass fraction based on ``basis``. The diluent is not - considered in the computation of the equivalence ratio. Default is no - dilution or ``fraction=None``. May be given as string or dictionary (for - example ``fraction={"fuel":0.7}``). - """ - cdef np.ndarray[np.double_t, ndim=1] fuel_comp = np.ascontiguousarray( - self.__composition_to_array(fuel, basis), dtype=np.double) - cdef np.ndarray[np.double_t, ndim=1] ox_comp = np.ascontiguousarray( - self.__composition_to_array(oxidizer, basis), dtype=np.double) - - cdef span[double] fuel_view = span[double](&fuel_comp[0], fuel_comp.size) - cdef span[double] ox_view = span[double](&ox_comp[0], ox_comp.size) - self.thermo.setEquivalenceRatio(phi, fuel_view, ox_view, - ThermoBasis.mass if basis == "mass" - else ThermoBasis.molar) - - if (fraction is None) != (diluent is None): - raise ValueError("If dilution is used, both 'fraction' and 'diluent' " - "parameters are required.") - - if fraction is None: - return - - if isinstance(fraction, str): - fraction_dict = comp_map_to_dict(parseCompString(stringify(fraction))) - elif isinstance(fraction, dict): - fraction_dict = fraction - else: - raise ValueError("The fraction argument must be given as string or " - "dictionary.") - - if len(fraction_dict) != 1: - raise ValueError("Invalid format for the fraction. Must be provided for " - "example as fraction='fuel:0.1'") - - fraction_type = list(fraction_dict.keys())[0] - fraction_value = float(list(fraction_dict.values())[0]) - - if fraction_value < 0 or fraction_value > 1: - raise ValueError("The fraction must be between 0 and 1") - - if fraction_type not in ["fuel", "oxidizer", "diluent"]: - raise ValueError("The fraction must specify 'fuel', 'oxidizer' or " - "'diluent'") - - cdef np.ndarray[np.double_t, ndim=1] diluent_comp = np.ascontiguousarray( - self.__composition_to_array(diluent, basis), dtype=np.double) - - # this function changes the composition and fixes temperature and pressure - T, P = self.T, self.P - # if 'fraction' is specified for diluent, just scale the mass or mole fractions - # of the fuel/oxidizer mixture accordingly - if fraction_type == "diluent": - if basis == "mole": - X_fuelox = self.X - self.X = diluent_comp - self.X = (1.0 - fraction_value) * X_fuelox + fraction_value * self.X - else: - Y_fuelox = self.Y - self.Y = diluent_comp - self.Y = (1.0 - fraction_value) * Y_fuelox + fraction_value * self.Y - self.TP = T, P - return - - # get the mixture fraction before scaling / diluent addition - Z_fuel = self.mixture_fraction(fuel, oxidizer, basis) - - if Z_fuel == 0.0 and fraction_type == "fuel": - raise ValueError("No fuel in the fuel/oxidizer mixture") - - if Z_fuel == 1.0 and fraction_type == "oxidizer": - raise ValueError("No oxidizer in the fuel/oxidizer mixture") - - if basis == "mass": # for mass basis, it is straight forward - if fraction_type == "fuel": - Z = Z_fuel - else: # oxidizer - Z = 1.0 - Z_fuel - if fraction_value > Z: - raise ValueError(f"The {fraction_type} fraction after dilution cannot " - "be higher than {fraction_type} fraction in the " - "original mixture.") - Y_mix = self.Y - self.Y = diluent_comp - factor = fraction_value / Z - self.Y = factor*Y_mix + (1.0 - factor) * self.Y - else: - # convert mass based mixture fraction to molar one, Z = kg fuel / kg mixture - X_mix = self.X - M_mix = self.mean_molecular_weight - self.X = fuel_comp - M_fuel = self.mean_molecular_weight - Z_fuel_mole = Z_fuel * M_mix / M_fuel # mol fuel / mol mixture - if fraction_type == "fuel": - Z = Z_fuel_mole - else: # oxidizer - Z = 1.0 - Z_fuel_mole - if fraction_value > Z: - raise ValueError(f"The {fraction_type} fuel or oxidizer fraction after " - "dilution cannot be higher than {fraction_type} " - "fraction in the original mixture.") - self.X = diluent_comp - factor = fraction_value / Z - self.X = factor * X_mix + (1.0 - factor) * self.X - self.TP = T, P - - def set_mixture_fraction(self, mixture_fraction, fuel, oxidizer, basis='mole'): - """ - Set the composition to a mixture of ``fuel`` and ``oxidizer`` at the - specified mixture fraction ``mixture_fraction`` (kg fuel / kg mixture), holding - temperature and pressure constant. Considers the oxidation of C to CO2, - H to H2O and S to SO2. Other elements are assumed not to participate in - oxidation (that is, N ends up as N2). The ``basis`` determines the composition - of fuel and oxidizer: ``basis='mole'`` (default) means mole fractions, - ``basis='mass'`` means mass fractions. For more information, see the - :doc:`Python example ` :: - - >>> gas.set_mixture_fraction(0.5, 'CH4', 'O2:1.0, N2:3.76') - >>> gas.mass_fraction_dict() - {'CH4': 0.5, 'N2': 0.38350014242997776, 'O2': 0.11649985757002226} - >>> gas.set_mixture_fraction(0.5, {'NH3':0.8, 'CO':0.2}, 'O2:1.0') - >>> gas.mass_fraction_dict() - {'CO': 0.145682068778996, 'NH3': 0.354317931221004, 'O2': 0.5} - - :param mixture_fraction: - Mixture fraction (kg fuel / kg mixture) - :param fuel: - Fuel species name or mole/mass fractions as string, array, or dict. - :param oxidizer: - Oxidizer species name or mole/mass fractions as a string, array, or - dict. - :param basis: determines if ``fuel`` and ``oxidizer`` are given in mole - fractions (``basis='mole'``) or mass fractions (``basis='mass'``) - """ - cdef np.ndarray[np.double_t, ndim=1] f = \ - np.ascontiguousarray(self.__composition_to_array(fuel, basis), dtype=np.double) - cdef np.ndarray[np.double_t, ndim=1] o = \ - np.ascontiguousarray(self.__composition_to_array(oxidizer, basis), dtype=np.double) - - cdef span[double] fuel_view = span[double](&f[0], f.size) - cdef span[double] ox_view = span[double](&o[0], o.size) - self.thermo.setMixtureFraction(mixture_fraction, fuel_view, ox_view, - ThermoBasis.mass if basis == 'mass' else ThermoBasis.molar) - - def equivalence_ratio(self, fuel=None, oxidizer=None, basis="mole", - include_species=None): - r""" - Get the equivalence ratio :math:`\phi` of the current mixture, which is a - conserved quantity. Considers the oxidation of C to CO2, H to H2O - and S to SO2. Other elements are assumed not to participate in oxidation - (that is, N ends up as N2). If fuel and oxidizer are not specified, the - equivalence ratio is computed from the available oxygen and the - required oxygen for complete oxidation: - - .. math:: \phi = \frac{Z_{\mathrm{mole},C} + Z_{\mathrm{mole},S} - + \frac{1}{4}Z_{\mathrm{mole},H}} {\frac{1}{2}Z_{\mathrm{mole},O}} - - where :math:`Z_{\mathrm{mole},e}` is the elemental mole fraction of element - :math:`e`. If the fuel and oxidizer compositions are specified, :math:`\phi` is - computed from: - - .. math:: \phi = \frac{Z}{1-Z}\frac{1-Z_{\mathrm{st}}}{Z_{\mathrm{st}}} - - where :math:`Z` is the Bilger mixture fraction and :math:`Z_{\mathrm{st}}` - the Bilger mixture fraction at stoichiometric conditions. - The ``basis`` determines the composition of fuel and oxidizer: - ``basis='mole'`` (default) means mole fractions, ``basis='mass'`` means - mass fractions. Note that this definition takes all species into account. - In case certain species like inert diluents should be ignored, a - list of species can be provided with ``include_species``. This means that - only these species are considered for the computation of the equivalence - ratio. For more information, see the - :doc:`Python example ` :: - - >>> gas.set_equivalence_ratio(0.5, fuel='CH3:0.5, CH3OH:.5, N2:0.125', oxidizer='O2:0.21, N2:0.79, NO:0.01') - >>> gas.equivalence_ratio(fuel='CH3:0.5, CH3OH:.5, N2:0.125', oxidizer='O2:0.21, N2:0.79, NO:0.01') - 0.5 - - :param fuel: - Fuel species name or mole/mass fractions as string, array, or dict. - :param oxidizer: - Oxidizer species name or mole/mass fractions as a string, array, or dict. - :param basis: - Determines if ``fuel`` and ``oxidizer`` are given in mole fractions - (``basis="mole"``) or mass fractions (``basis="mass"``) - :param include_species: - List of species names (optional). Only these species are considered for the - computation of the equivalence ratio. By default, all species are considered - """ - if include_species is not None: - # remove unwanted species temporarily - Y = np.zeros(self.n_species) - indices = [self.species_index(s) for s in include_species] - for k in indices: - Y[k] = self.Y[k] - T_orig, P_orig, Y_orig = self.T, self.P, self.Y - self.Y = Y - - if fuel is None and oxidizer is None: - phi = self.thermo.equivalenceRatio() - if include_species is not None: - self.TPY = T_orig, P_orig, Y_orig - return phi - - cdef np.ndarray[np.double_t, ndim=1] f = np.ascontiguousarray( - self.__composition_to_array(fuel, basis), dtype=np.double) - cdef np.ndarray[np.double_t, ndim=1] o = np.ascontiguousarray( - self.__composition_to_array(oxidizer, basis), dtype=np.double) - - cdef span[double] fuel_view = span[double](&f[0], f.size) - cdef span[double] ox_view = span[double](&o[0], o.size) - phi = self.thermo.equivalenceRatio(fuel_view, ox_view, - ThermoBasis.mass if basis == "mass" else ThermoBasis.molar) - if include_species is not None: - self.TPY = T_orig, P_orig, Y_orig - return phi - - def mixture_fraction(self, fuel, oxidizer, basis='mole', element="Bilger"): - r""" - Get the mixture fraction of the current mixture in - (kg fuel / (kg oxidizer + kg fuel)). This is a quantity that is conserved after - oxidation. Considers the oxidation of C to CO2, H to H2O and S to SO2. Other - elements are assumed not to participate in oxidation (that is, N ends up as N2). - The ``basis`` determines the composition of fuel and oxidizer: - ``basis="mole"`` (default) means mole fractions, ``basis="mass"`` means mass - fractions. The mixture fraction can be computed from a single element (for - example, carbon with ``element="C"``) - - .. math:: Z_m = \frac{Z_{\mathrm{mass},m}-Z_{\mathrm{mass},m,\mathrm{ox}}} - {Z_{\mathrm{mass},\mathrm{fuel}}-Z_{\mathrm{mass},m,\mathrm{ox}}} - - where :math:`Z_{\mathrm{mass},m}` is the elemental mass fraction of - element :math:`m` in the mixture, and :math:`Z_{\mathrm{mass},m,\mathrm{ox}}` - and :math:`Z_{\mathrm{mass},\mathrm{fuel}}` are the elemental mass fractions of - the oxidizer and fuel, or from the Bilger mixture fraction - (``element="Bilger"``), which considers the elements C, S, H and O - :cite:p:`bilger1979`. The Bilger mixture fraction is computed by default: - - .. math:: Z_m = Z_{\mathrm{Bilger}} = \frac{\beta-\beta_{\mathrm{ox}}} - {\beta_{\mathrm{fuel}}-\beta_{\mathrm{ox}}} - - with - - .. math:: \beta = 2\frac{Z_C}{M_C}+2\frac{Z_S}{M_S}+\frac{1}{2}\frac{Z_H}{M_H} - - \frac{Z_O}{M_O} - - and :math:`M_m` the atomic weight of element :math:`m`. - For more information, see the - :doc:`Python example ` :: - - >>> gas.set_mixture_fraction(0.5, 'CH3:0.5, CH3OH:0.5, N2:0.125', 'O2:0.21, N2:0.79, NO:0.01') - >>> gas.mixture_fraction('CH3:0.5, CH3OH:0.5, N2:0.125', 'O2:0.21, N2:0.79, NO:.01') - 0.5 - - :param fuel: - Fuel species name or mole/mass fractions as string, array, or dict. - :param oxidizer: - Oxidizer species name or mole/mass fractions as a string, array, or - dict. - :param basis: - Determines if ``fuel`` and ``oxidizer`` are given in mole - fractions (``basis='mole'``) or mass fractions (``basis='mass'``) - :param element: - Computes the mixture fraction from the specified elemental - mass fraction (given by element name or element index) or as - the Bilger mixture fraction (default) - """ - cdef np.ndarray[np.double_t, ndim=1] f = \ - np.ascontiguousarray(self.__composition_to_array(fuel, basis), dtype=np.double) - cdef np.ndarray[np.double_t, ndim=1] o = \ - np.ascontiguousarray(self.__composition_to_array(oxidizer, basis), dtype=np.double) - - if isinstance(element, (str, bytes)): - e_name = element - else: - e_name = self.element_name(self.element_index(element)) - - cdef span[double] fuel_view = span[double](&f[0], f.size) - cdef span[double] ox_view = span[double](&o[0], o.size) - return self.thermo.mixtureFraction( - fuel_view, ox_view, - ThermoBasis.mass if basis=='mass' else ThermoBasis.molar, - stringify(e_name)) - - def stoich_air_fuel_ratio(self, fuel, oxidizer, basis='mole'): - """ - Get the stoichiometric air to fuel ratio (kg oxidizer / kg fuel). Considers the - oxidation of C to CO2, H to H2O and S to SO2. Other elements are assumed - not to participate in oxidation (that is, N ends up as N2). - The ``basis`` determines the composition of fuel and oxidizer: ``basis='mole'`` (default) - means mole fractions, ``basis='mass'`` means mass fractions:: - - >>> gas.set_mixture_fraction(0.5, 'CH3:0.5, CH3OH:.5, N2:0.125', 'O2:0.21, N2:0.79, NO:0.01') - >>> gas.stoich_air_fuel_ratio('CH3:0.5, CH3OH:.5, N2:0.125', 'O2:0.21, N2:0.79, NO:0.01') - 8.148040722239438 - - :param fuel: - Fuel species name or mole/mass fractions as string, array, or dict. - :param oxidizer: - Oxidizer species name or mole/mass fractions as a string, array, or - dict. - :param basis: - Determines if ``fuel`` and ``oxidizer`` are given in mole - fractions (``basis='mole'``) or mass fractions (``basis='mass'``) - - """ - cdef np.ndarray[np.double_t, ndim=1] f = \ - np.ascontiguousarray(self.__composition_to_array(fuel, basis), dtype=np.double) - cdef np.ndarray[np.double_t, ndim=1] o = \ - np.ascontiguousarray(self.__composition_to_array(oxidizer, basis), dtype=np.double) - - cdef span[double] fuel_view = span[double](&f[0], f.size) - cdef span[double] ox_view = span[double](&o[0], o.size) - return self.thermo.stoichAirFuelRatio(fuel_view, ox_view, - ThermoBasis.mass if basis=='mass' else ThermoBasis.molar) - - def elemental_mass_fraction(self, m): - r""" - Get the elemental mass fraction :math:`Z_{\mathrm{mass},m}` of element - :math:`m` as defined by: - - .. math:: Z_{\mathrm{mass},m} = \sum_k \frac{a_{m,k} M_m}{M_k} Y_k - - with :math:`a_{m,k}` being the number of atoms of element :math:`m` in - species :math:`k`, :math:`M_m` the atomic weight of element :math:`m`, - :math:`M_k` the molecular weight of species :math:`k`, and :math:`Y_k` - the mass fraction of species :math:`k`:: - - >>> phase.elemental_mass_fraction('H') - 1.0 - - :param m: - Base element, may be specified by name or by index. - """ - return self.thermo.elementalMassFraction(self.element_index(m)) - - def elemental_mole_fraction(self, m): - r""" - Get the elemental mole fraction :math:`Z_{\mathrm{mole},m}` of element - :math:`m` (the number of atoms of element m divided by the total number - of atoms) as defined by: - - .. math:: Z_{\mathrm{mole},m} = \frac{\sum_k a_{m,k} X_k} - {\sum_k \sum_j a_{j,k} X_k} - - with :math:`a_{m,k}` being the number of atoms of element :math:`m` in - species :math:`k`, :math:`\sum_j` being a sum over all elements, and - :math:`X_k` being the mole fraction of species :math:`k`:: - - >>> phase.elemental_mole_fraction('H') - 1.0 - - :param m: - Base element, may be specified by name or by index. - """ - return self.thermo.elementalMoleFraction(self.element_index(m)) - - def set_unnormalized_mass_fractions(self, Y): - """ - Set the mass fractions without normalizing to force ``sum(Y) == 1.0``. - Useful primarily when calculating derivatives with respect to ``Y[k]`` by - finite difference. - """ - cdef np.ndarray[np.double_t, ndim=1] data - if len(Y) == self.n_species: - data = np.ascontiguousarray(Y, dtype=np.double) - else: - raise ValueError("Array has incorrect length." - " Got {}, expected {}.".format(len(Y), self.n_species)) - self.thermo.setMassFractions_NoNorm(span[double](&data[0], data.size)) - - def set_unnormalized_mole_fractions(self, X): - """ - Set the mole fractions without normalizing to force ``sum(X) == 1.0``. - Useful primarily when calculating derivatives with respect to ``X[k]`` - by finite difference. - """ - cdef np.ndarray[np.double_t, ndim=1] data - if len(X) == self.n_species: - data = np.ascontiguousarray(X, dtype=np.double) - else: - raise ValueError("Array has incorrect length." - " Got {}, expected {}.".format(len(X), self.n_species)) - self.thermo.setMoleFractions_NoNorm(span[double](&data[0], data.size)) - - def mass_fraction_dict(self, double threshold=0.0): - """ - Return a dictionary giving the mass fraction for each species by name where the - mass fraction is greater than ``threshold``. - """ - cdef pair[string,double] item - Y = self.thermo.getMassFractionsByName(threshold) - return {pystr(item.first):item.second for item in Y} - - def mole_fraction_dict(self, double threshold=0.0): - """ - Return a dictionary giving the mole fraction for each species by name where the - mole fraction is greater than ``threshold``. - """ - cdef pair[string,double] item - X = self.thermo.getMoleFractionsByName(threshold) - return {pystr(item.first):item.second for item in X} - - - ######## Read-only thermodynamic properties ######## - - property P: - """Pressure [Pa].""" - def __get__(self): - return self.thermo.pressure() - - property T: - """Temperature [K].""" - def __get__(self): - return self.thermo.temperature() - - property density: - """Density [kg/m³ or kmol/m³] depending on `basis`.""" - def __get__(self): - return self.thermo.density() / self._mass_factor() - - property density_mass: - """(Mass) density [kg/m³].""" - def __get__(self): - return self.thermo.density() - - property density_mole: - """Molar density [kmol/m³].""" - def __get__(self): - return self.thermo.molarDensity() - - property v: - """Specific volume [m³/kg or m³/kmol] depending on `basis`.""" - def __get__(self): - return self._mass_factor() / self.thermo.density() - - property volume_mass: - """Specific volume [m³/kg].""" - def __get__(self): - return 1.0 / self.thermo.density() - - property volume_mole: - """Molar volume [m³/kmol].""" - def __get__(self): - return self.thermo.molarVolume() - - property u: - """Internal energy in [J/kg or J/kmol].""" - def __get__(self): - return self.thermo.intEnergy_mole() * self._mole_factor() - - property int_energy_mole: - """Molar internal energy [J/kmol].""" - def __get__(self): - return self.thermo.intEnergy_mole() - - property int_energy_mass: - """Specific internal energy [J/kg].""" - def __get__(self): - return self.thermo.intEnergy_mass() - - property h: - """Enthalpy [J/kg or J/kmol] depending on `basis`.""" - def __get__(self): - return self.thermo.enthalpy_mole() * self._mole_factor() - - property enthalpy_mole: - """Molar enthalpy [J/kmol].""" - def __get__(self): - return self.thermo.enthalpy_mole() - - property enthalpy_mass: - """Specific enthalpy [J/kg].""" - def __get__(self): - return self.thermo.enthalpy_mass() - - property s: - """Entropy [J/kg/K or J/kmol/K] depending on `basis`.""" - def __get__(self): - return self.thermo.entropy_mole() * self._mole_factor() - - property entropy_mole: - """Molar entropy [J/kmol/K].""" - def __get__(self): - return self.thermo.entropy_mole() - - property entropy_mass: - """Specific entropy [J/kg/K].""" - def __get__(self): - return self.thermo.entropy_mass() - - property g: - """Gibbs free energy [J/kg or J/kmol] depending on `basis`.""" - def __get__(self): - return self.thermo.gibbs_mole() * self._mole_factor() - - property gibbs_mole: - """Molar Gibbs free energy [J/kmol].""" - def __get__(self): - return self.thermo.gibbs_mole() - - property gibbs_mass: - """Specific Gibbs free energy [J/kg].""" - def __get__(self): - return self.thermo.gibbs_mass() - - property cv: - """ - Heat capacity at constant volume and composition - [J/kg/K or J/kmol/K depending on `basis`]. - """ - def __get__(self): - return self.thermo.cv_mole() * self._mole_factor() - - property cv_mole: - """Molar heat capacity at constant volume and composition [J/kmol/K].""" - def __get__(self): - return self.thermo.cv_mole() - - property cv_mass: - """Specific heat capacity at constant volume and composition [J/kg/K].""" - def __get__(self): - return self.thermo.cv_mass() - - property cp: - """ - Heat capacity at constant pressure and composition. - [J/kg/K or J/kmol/K depending on `basis`] - """ - def __get__(self): - return self.thermo.cp_mole() * self._mole_factor() - - property cp_mole: - """Molar heat capacity at constant pressure and composition [J/kmol/K].""" - def __get__(self): - return self.thermo.cp_mole() - - property cp_mass: - """Specific heat capacity at constant pressure and composition [J/kg/K].""" - def __get__(self): - return self.thermo.cp_mass() - - property critical_temperature: - """Critical temperature [K].""" - def __get__(self): - return self.thermo.critTemperature() - - property critical_pressure: - """Critical pressure [Pa].""" - def __get__(self): - return self.thermo.critPressure() - - property critical_density: - """Critical density [kg/m³ or kmol/m³] depending on `basis`.""" - def __get__(self): - return self.thermo.critDensity() / self._mass_factor() - - property P_sat: - """Saturation pressure [Pa] at the current temperature.""" - def __get__(self): - return self.thermo.satPressure(self.T) - - property T_sat: - """Saturation temperature [K] at the current pressure.""" - def __get__(self): - return self.thermo.satTemperature(self.P) - - property auxiliary_data: - """ - Intermediate or model-specific parameters used by particular - derived classes. - """ - def __get__(self): - cdef CxxAnyMap values = self.thermo.getAuxiliaryData() - return anymap_to_py(values) - - ######## Methods to get/set the complete thermodynamic state ######## - - property state_size: - """ - Return size of vector defining internal state of the phase. - """ - def __get__(self): - return self.thermo.stateSize() - - property state: - """ - Get/Set the full thermodynamic state as a single array, arranged as - [temperature, density, mass fractions] for most phases. Useful mainly - in cases where it is desired to store many states in a multidimensional - array. - """ - def __get__(self): - cdef np.ndarray[np.double_t, ndim=1] state = np.empty(self.state_size) - self.thermo.saveState(span[double](&state[0], len(state))) - return state - - def __set__(self, state): - cdef np.ndarray[np.double_t, ndim=1] cstate = np.asarray(state) - self.thermo.restoreState(span[double](&cstate[0], len(cstate))) - - property TD: - """Get/Set temperature [K] and density [kg/m³ or kmol/m³].""" - def __get__(self): - return self.T, self.density - def __set__(self, values): - assert len(values) == 2, 'incorrect number of values' - T = values[0] if values[0] is not None else self.T - D = values[1] if values[1] is not None else self.density - self.thermo.setState_TD(T, D * self._mass_factor()) - - property TDX: - """ - Get/Set temperature [K], density [kg/m³ or kmol/m³], and mole fractions. - """ - def __get__(self): - return self.T, self.density, self.X - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - T = values[0] if values[0] is not None else self.T - D = values[1] if values[1] is not None else self.density - self.X = values[2] - self.thermo.setState_TD(T, D * self._mass_factor()) - - property TDY: - """ - Get/Set temperature [K] and density [kg/m³ or kmol/m³], and mass fractions. - """ - def __get__(self): - return self.T, self.density, self.Y - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - T = values[0] if values[0] is not None else self.T - D = values[1] if values[1] is not None else self.density - self.Y = values[2] - self.thermo.setState_TD(T, D * self._mass_factor()) - - property TP: - """Get/Set temperature [K] and pressure [Pa].""" - def __get__(self): - return self.T, self.P - def __set__(self, values): - assert len(values) == 2, 'incorrect number of values' - T = values[0] if values[0] is not None else self.T - P = values[1] if values[1] is not None else self.P - self.thermo.setState_TP(T, P) - - property TPX: - """Get/Set temperature [K], pressure [Pa], and mole fractions.""" - def __get__(self): - return self.T, self.P, self.X - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - T = values[0] if values[0] is not None else self.T - P = values[1] if values[1] is not None else self.P - self.X = values[2] - self.thermo.setState_TP(T, P) - - property TPY: - """Get/Set temperature [K], pressure [Pa], and mass fractions.""" - def __get__(self): - return self.T, self.P, self.Y - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - T = values[0] if values[0] is not None else self.T - P = values[1] if values[1] is not None else self.P - self.Y = values[2] - self.thermo.setState_TP(T, P) - - property UV: - """ - Get/Set internal energy [J/kg or J/kmol] and specific volume [m³/kg or m³/kmol]. - """ - def __get__(self): - return self.u, self.v - def __set__(self, values): - assert len(values) == 2, 'incorrect number of values' - U = values[0] if values[0] is not None else self.u - V = values[1] if values[1] is not None else self.v - self.thermo.setState_UV(U / self._mass_factor(), - V / self._mass_factor()) - - property UVX: - """ - Get/Set internal energy [J/kg or J/kmol], specific volume [m³/kg or m³/kmol], - and mole fractions. - """ - def __get__(self): - return self.u, self.v, self.X - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - U = values[0] if values[0] is not None else self.u - V = values[1] if values[1] is not None else self.v - self.X = values[2] - self.thermo.setState_UV(U / self._mass_factor(), - V / self._mass_factor()) - - property UVY: - """ - Get/Set internal energy [J/kg or J/kmol], specific volume [m³/kg or m³/kmol], - and mass fractions. - """ - def __get__(self): - return self.u, self.v, self.Y - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - U = values[0] if values[0] is not None else self.u - V = values[1] if values[1] is not None else self.v - self.Y = values[2] - self.thermo.setState_UV(U / self._mass_factor(), - V / self._mass_factor()) - - property DP: - """Get/Set density [kg/m³] and pressure [Pa].""" - def __get__(self): - return self.density, self.P - def __set__(self, values): - assert len(values) == 2, 'incorrect number of values' - D = values[0] if values[0] is not None else self.density - P = values[1] if values[1] is not None else self.P - self.thermo.setState_DP(D*self._mass_factor(), P) - - property DPX: - """Get/Set density [kg/m³], pressure [Pa], and mole fractions.""" - def __get__(self): - return self.density, self.P, self.X - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - D = values[0] if values[0] is not None else self.density - P = values[1] if values[1] is not None else self.P - self.X = values[2] - self.thermo.setState_DP(D*self._mass_factor(), P) - - property DPY: - """Get/Set density [kg/m³], pressure [Pa], and mass fractions.""" - def __get__(self): - return self.density, self.P, self.Y - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - D = values[0] if values[0] is not None else self.density - P = values[1] if values[1] is not None else self.P - self.Y = values[2] - self.thermo.setState_DP(D*self._mass_factor(), P) - - property HP: - """Get/Set enthalpy [J/kg or J/kmol] and pressure [Pa].""" - def __get__(self): - return self.h, self.P - def __set__(self, values): - assert len(values) == 2, 'incorrect number of values' - H = values[0] if values[0] is not None else self.h - P = values[1] if values[1] is not None else self.P - self.thermo.setState_HP(H / self._mass_factor(), P) - - property HPX: - """Get/Set enthalpy [J/kg or J/kmol], pressure [Pa] and mole fractions.""" - def __get__(self): - return self.h, self.P, self.X - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - H = values[0] if values[0] is not None else self.h - P = values[1] if values[1] is not None else self.P - self.X = values[2] - self.thermo.setState_HP(H / self._mass_factor(), P) - - property HPY: - """Get/Set enthalpy [J/kg or J/kmol], pressure [Pa] and mass fractions.""" - def __get__(self): - return self.h, self.P, self.Y - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - H = values[0] if values[0] is not None else self.h - P = values[1] if values[1] is not None else self.P - self.Y = values[2] - self.thermo.setState_HP(H / self._mass_factor(), P) - - property SP: - """Get/Set entropy [J/kg/K or J/kmol/K] and pressure [Pa].""" - def __get__(self): - return self.s, self.P - def __set__(self, values): - assert len(values) == 2, 'incorrect number of values' - S = values[0] if values[0] is not None else self.s - P = values[1] if values[1] is not None else self.P - self.thermo.setState_SP(S / self._mass_factor(), P) - - property SPX: - """Get/Set entropy [J/kg/K or J/kmol/K], pressure [Pa], and mole fractions.""" - def __get__(self): - return self.s, self.P, self.X - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - S = values[0] if values[0] is not None else self.s - P = values[1] if values[1] is not None else self.P - self.X = values[2] - self.thermo.setState_SP(S / self._mass_factor(), P) - - property SPY: - """Get/Set entropy [J/kg/K or J/kmol/K], pressure [Pa], and mass fractions.""" - def __get__(self): - return self.s, self.P, self.Y - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - S = values[0] if values[0] is not None else self.s - P = values[1] if values[1] is not None else self.P - self.Y = values[2] - self.thermo.setState_SP(S / self._mass_factor(), P) - - property SV: - """ - Get/Set entropy [J/kg/K or J/kmol/K] and specific volume [m³/kg or m³/kmol]. - """ - def __get__(self): - return self.s, self.v - def __set__(self, values): - assert len(values) == 2, 'incorrect number of values' - S = values[0] if values[0] is not None else self.s - V = values[1] if values[1] is not None else self.v - self.thermo.setState_SV(S / self._mass_factor(), - V / self._mass_factor()) - - property SVX: - """ - Get/Set entropy [J/kg/K or J/kmol/K], specific volume [m³/kg or m³/kmol], and - mole fractions. - """ - def __get__(self): - return self.s, self.v, self.X - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - S = values[0] if values[0] is not None else self.s - V = values[1] if values[1] is not None else self.v - self.X = values[2] - self.thermo.setState_SV(S / self._mass_factor(), - V / self._mass_factor()) - - property SVY: - """ - Get/Set entropy [J/kg/K or J/kmol/K], specific volume [m³/kg or m³/kmol], and - mass fractions. - """ - def __get__(self): - return self.s, self.v, self.Y - def __set__(self, values): - assert len(values) == 3, 'incorrect number of values' - S = values[0] if values[0] is not None else self.s - V = values[1] if values[1] is not None else self.v - self.Y = values[2] - self.thermo.setState_SV(S / self._mass_factor(), - V / self._mass_factor()) - - # partial molar / non-dimensional properties - property partial_molar_enthalpies: - """Array of species partial molar enthalpies [J/kmol].""" - def __get__(self): - return self._getArray1(thermo_getPartialMolarEnthalpies) - - property partial_molar_entropies: - """Array of species partial molar entropies [J/kmol/K].""" - def __get__(self): - return self._getArray1(thermo_getPartialMolarEntropies) - - property partial_molar_int_energies: - """Array of species partial molar internal energies [J/kmol].""" - def __get__(self): - return self._getArray1(thermo_getPartialMolarIntEnergies) - - property partial_molar_int_energies_TV: - r""" - Array of species `\tilde{u}_k = (\partial U/\partial n_k)_{T,V,n_{j\ne k}}` - values [J/kmol]. - - .. versionadded:: 4.0 - """ - def __get__(self): - return self._getArray1(thermo_getPartialMolarIntEnergies_TV) - - property chemical_potentials: - """Array of species chemical potentials [J/kmol].""" - def __get__(self): - return self._getArray1(thermo_getChemPotentials) - - property electrochemical_potentials: - """Array of species electrochemical potentials [J/kmol].""" - def __get__(self): - return self._getArray1(thermo_getElectrochemPotentials) - - property partial_molar_cp: - """ - Array of species partial molar specific heat capacities at constant - pressure [J/kmol/K]. - """ - def __get__(self): - return self._getArray1(thermo_getPartialMolarCp) - - property partial_molar_cv_TV: - r""" - Array of species `\tilde{c}_{v,k} = (\partial \tilde{u}_k/\partial T)_{V,n}` - values [J/kmol/K]. - - .. versionadded:: 4.0 - """ - def __get__(self): - return self._getArray1(thermo_getPartialMolarCv_TV) - - property partial_molar_volumes: - """Array of species partial molar volumes [m³/kmol].""" - def __get__(self): - return self._getArray1(thermo_getPartialMolarVolumes) - - property standard_enthalpies_RT: - """ - Array of nondimensional species standard-state enthalpies at the - current temperature and pressure. - """ - def __get__(self): - return self._getArray1(thermo_getEnthalpy_RT) - - property standard_entropies_R: - """ - Array of nondimensional species standard-state entropies at the - current temperature and pressure. - """ - def __get__(self): - return self._getArray1(thermo_getEntropy_R) - - property standard_int_energies_RT: - """ - Array of nondimensional species standard-state internal energies at the - current temperature and pressure. - """ - def __get__(self): - return self._getArray1(thermo_getIntEnergy_RT) - - property standard_gibbs_RT: - """ - Array of nondimensional species standard-state Gibbs free energies at - the current temperature and pressure. - """ - def __get__(self): - return self._getArray1(thermo_getGibbs_RT) - - property standard_cp_R: - """ - Array of nondimensional species standard-state specific heat capacities - at constant pressure at the current temperature and pressure. - """ - def __get__(self): - return self._getArray1(thermo_getCp_R) - - property activities: - """ - Array of nondimensional activities. Returns either molar or molal - activities depending on the convention of the thermodynamic model. - """ - def __get__(self): - return self._getArray1(thermo_getActivities) - - property activity_coefficients: - """ - Array of nondimensional, molar activity coefficients. - """ - def __get__(self): - return self._getArray1(thermo_getActivityCoefficients) - - ######## Miscellaneous properties ######## - property isothermal_compressibility: - """Isothermal compressibility [1/Pa].""" - def __get__(self): - return self.thermo.isothermalCompressibility() - - property thermal_expansion_coeff: - """Thermal expansion coefficient [1/K].""" - def __get__(self): - return self.thermo.thermalExpansionCoeff() - - property internal_pressure: - """ - Internal pressure [Pa]. - - .. versionadded:: 4.0 - """ - def __get__(self): - return self.thermo.internalPressure() - - property sound_speed: - """Speed of sound [m/s].""" - def __get__(self): - return self.thermo.soundSpeed() - - property min_temp: - """ - Minimum temperature for which the thermodynamic data for the phase are - valid. - """ - def __get__(self): - return self.thermo.minTemp() - - property max_temp: - """ - Maximum temperature for which the thermodynamic data for the phase are - valid. - """ - def __get__(self): - return self.thermo.maxTemp() - - property reference_pressure: - """Reference state pressure [Pa].""" - def __get__(self): - return self.thermo.refPressure() - - property electric_potential: - """Get/Set the electric potential [V] for this phase.""" - def __get__(self): - return self.thermo.electricPotential() - def __set__(self, double value): - self.thermo.setElectricPotential(value) - - property standard_concentration_units: - """Get standard concentration units for this phase.""" - def __get__(self): - cdef CxxUnits units = self.thermo.standardConcentrationUnits() - return Units.copy(units) - - # methods for plasma - property Te: - """Get/Set electron Temperature [K].""" - def __get__(self): - return self.thermo.electronTemperature() - - def __set__(self, value): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - self.plasma.setElectronTemperature(value) - - property Pe: - """Get electron Pressure [Pa].""" - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return self.plasma.electronPressure() - - property reduced_electric_field: - """ - Get/Set the reduced electric field (E/N) [V·m²]. - - .. versionadded:: 3.2 - """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return self.plasma.reducedElectricField() - - def __set__(self, value): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - self.plasma.setReducedElectricField(value) - - property electric_field: - """ - Get/Set the electric field (E) [V/m]. - - This is the absolute electric field strength. It is related to the reduced - electric field (E/N) through the number density of neutrals. - - .. versionadded:: 3.2 - """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return self.plasma.electricField() - - def __set__(self, value): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - self.plasma.setElectricField(value) - - def set_discretized_electron_energy_distribution(self, levels, distribution): - """ - Set electron energy distribution. When this method is used, electron - temperature is calculated from the distribution. - - :param levels: - vector of electron energy levels [eV] - :param distribution: - vector of distribution - """ - if not self._enable_plasma: - raise TypeError('This method is invalid for ' - f'thermo model: {self.thermo_model}.') - # check length - if (len(levels) != len(distribution)): - raise ValueError('Length of levels and distribution are different') - - cdef np.ndarray[np.double_t, ndim=1] data_levels = \ - np.ascontiguousarray(levels, dtype=np.double) - cdef np.ndarray[np.double_t, ndim=1] data_dist = \ - np.ascontiguousarray(distribution, dtype=np.double) - self.plasma.setDiscretizedElectronEnergyDist( - span[double](&data_levels[0], data_levels.size), - span[double](&data_dist[0], data_dist.size)) - - def update_electron_energy_distribution(self): - """ - Update the electron energy distribution function to account for changes in - composition, temperature, pressure, or electric field strength. - - .. versionadded:: 3.2 - """ - if not self._enable_plasma: - raise TypeError('This method is invalid for ' - f'thermo model: {self.thermo_model}.') - self.plasma.updateElectronEnergyDistribution() - - property n_electron_energy_levels: - """ Number of electron energy levels """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return self.plasma.nElectronEnergyLevels() - - property electron_energy_levels: - """ Electron energy levels [eV]""" - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - cdef np.ndarray[np.double_t, ndim=1] data = np.empty( - self.n_electron_energy_levels) - self.plasma.getElectronEnergyLevels(span[double](&data[0], data.size)) - return data - def __set__(self, levels): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - cdef np.ndarray[np.double_t, ndim=1] data = \ - np.ascontiguousarray(levels, dtype=np.double) - self.plasma.setElectronEnergyLevels(span[double](&data[0], data.size)) - - property electron_energy_distribution: - """ Electron energy distribution """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - cdef np.ndarray[np.double_t, ndim=1] data = np.empty( - self.n_electron_energy_levels) - self.plasma.getElectronEnergyDistribution(span[double](&data[0], data.size)) - return data - - property isotropic_shape_factor: - """ Shape factor of isotropic-velocity distribution for electron energy """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return self.plasma.isotropicShapeFactor() - def __set__(self, x): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - self.plasma.setIsotropicShapeFactor(x) - - property electron_energy_distribution_type: - """ Electron energy distribution type """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return pystr(self.plasma.electronEnergyDistributionType()) - def __set__(self, distribution_type): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - self.plasma.setElectronEnergyDistributionType(stringify(distribution_type)) - - property mean_electron_energy: - """ Mean electron energy [eV] """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return self.plasma.meanElectronEnergy() - def __set__(self, double energy): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - self.plasma.setMeanElectronEnergy(energy) - - property quadrature_method: - """ Quadrature method """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return pystr(self.plasma.quadratureMethod()) - def __set__(self, method): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - self.plasma.setQuadratureMethod(stringify(method)) - - property normalize_electron_energy_distribution_enabled: - """ Automatically normalize electron energy distribution """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return self.plasma.normalizeElectronEnergyDistEnabled() - def __set__(self, enable): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - self.plasma.enableNormalizeElectronEnergyDist(enable) - - property electron_species_name: - """ Electron species name """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return pystr(self.plasma.electronSpeciesName()) - - property elastic_power_loss: - """ - Elastic power loss (W/m³) - - .. versionadded:: 3.2 - """ - def __get__(self): - if not self._enable_plasma: - raise ThermoModelMethodError(self.thermo_model) - return self.plasma.elasticPowerLoss() - -cdef class InterfacePhase(ThermoPhase): - """ A class representing a surface, edge phase """ - def __cinit__(self, *args, **kwargs): - if not kwargs.get("init", True): - return - if not dynamic_cast[CxxSurfPhasePtr](self.thermo): - raise TypeError('Underlying ThermoPhase object is of the wrong type.') - self.surf = (self.thermo) - - property adjacent: - """ - A dictionary containing higher-dimensional phases adjacent to this interface, - for example bulk phases adjacent to a surface. - """ - def __get__(self): - return self._adjacent - - property site_density: - """ - Get/Set the site density. [kmol/m²] for surface phases; [kmol/m] for - edge phases. - """ - def __get__(self): - return self.surf.siteDensity() - def __set__(self, double value): - self.surf.setSiteDensity(value) - - property coverages: - """Get/Set the fraction of sites covered by each species.""" - def __get__(self): - cdef np.ndarray[np.double_t, ndim=1] data = np.empty(self.n_species) - self.surf.getCoverages(span[double](&data[0], data.size)) - if self._selected_species.size: - return data[self._selected_species] - else: - return data - - def __set__(self, theta): - if isinstance(theta, (dict, str, bytes)): - self.surf.setCoveragesByName(comp_map(theta)) - return - - if len(theta) != self.n_species: - raise ValueError("Array has incorrect length." - " Got {}, expected {}".format(len(theta), self.n_species)) - cdef np.ndarray[np.double_t, ndim=1] data = \ - np.ascontiguousarray(theta, dtype=np.double) - self.surf.setCoverages(span[double](&data[0], data.size)) - - def set_unnormalized_coverages(self, cov): - """ - Set the surface coverages without normalizing to force sum(cov) == 1.0. - Useful primarily when calculating derivatives with respect to cov[k] by - finite difference. - """ - cdef np.ndarray[np.double_t, ndim=1] data - if len(cov) == self.n_species: - data = np.ascontiguousarray(cov, dtype=np.double) - else: - raise ValueError("Array has incorrect length." - " Got {}, expected {}.".format(len(cov), self.n_species)) - self.surf.setCoveragesNoNorm(span[double](&data[0], data.size)) - - -cdef class PureFluid(ThermoPhase): - """ - A pure substance that can be a gas, a liquid, a mixed gas-liquid fluid, - or a fluid beyond its critical point. - """ - - property Q: - """ - Get/Set vapor fraction (quality). Can be set only when in the two-phase - region. - """ - def __get__(self): - return self.thermo.vaporFraction() - def __set__(self, Q): - if (self.P >= self.critical_pressure or - abs(self.P-self.P_sat)/self.P > 1e-4): - raise ValueError('Cannot set vapor quality outside the ' - 'two-phase region') - self.thermo.setState_Psat(self.P, Q) - - property TQ: - """Get/Set the temperature [K] and vapor fraction of a two-phase state.""" - def __get__(self): - return self.T, self.Q - def __set__(self, values): - T = values[0] if values[0] is not None else self.T - Q = values[1] if values[1] is not None else self.Q - self.thermo.setState_Tsat(T, Q) - - property PQ: - """Get/Set the pressure [Pa] and vapor fraction of a two-phase state.""" - def __get__(self): - return self.P, self.Q - def __set__(self, values): - P = values[0] if values[0] is not None else self.P - Q = values[1] if values[1] is not None else self.Q - self.thermo.setState_Psat(P, Q) - - property ST: - """Get/Set the entropy [J/kg/K] and temperature [K] of a PureFluid.""" - def __get__(self): - return self.s, self.T - def __set__(self, values): - S = values[0] if values[0] is not None else self.s - T = values[1] if values[1] is not None else self.T - self.thermo.setState_ST(S / self._mass_factor(), T) - - property TV: - """ - Get/Set the temperature [K] and specific volume [m³/kg] of a PureFluid. - """ - def __get__(self): - return self.T, self.v - def __set__(self, values): - T = values[0] if values[0] is not None else self.T - V = values[1] if values[1] is not None else self.v - self.thermo.setState_TV(T, V / self._mass_factor()) - - property PV: - """ - Get/Set the pressure [Pa] and specific volume [m³/kg] of a PureFluid. - """ - def __get__(self): - return self.P, self.v - def __set__(self, values): - P = values[0] if values[0] is not None else self.P - V = values[1] if values[1] is not None else self.v - self.thermo.setState_PV(P, V / self._mass_factor()) - - property UP: - """ - Get/Set the specific internal energy [J/kg] and the - pressure [Pa] of a PureFluid. - """ - def __get__(self): - return self.u, self.P - def __set__(self, values): - U = values[0] if values[0] is not None else self.u - P = values[1] if values[1] is not None else self.P - self.thermo.setState_UP(U / self._mass_factor(), P) - - property VH: - """ - Get/Set the specific volume [m³/kg] and the specific - enthalpy [J/kg] of a PureFluid. - """ - def __get__(self): - return self.v, self.h - def __set__(self, values): - V = values[0] if values[0] is not None else self.v - H = values[1] if values[1] is not None else self.h - self.thermo.setState_VH(V/self._mass_factor(), H/self._mass_factor()) - - property TH: - """ - Get/Set the temperature [K] and the specific enthalpy [J/kg] - of a PureFluid. - """ - def __get__(self): - return self.T, self.h - def __set__(self, values): - T = values[0] if values[0] is not None else self.T - H = values[1] if values[1] is not None else self.h - self.thermo.setState_TH(T, H / self._mass_factor()) - - property SH: - """ - Get/Set the specific entropy [J/kg/K] and the specific - enthalpy [J/kg] of a PureFluid. - """ - def __get__(self): - return self.s, self.h - def __set__(self, values): - S = values[0] if values[0] is not None else self.s - H = values[1] if values[1] is not None else self.h - self.thermo.setState_SH(S/self._mass_factor(), H/self._mass_factor()) - - property TDQ: - """ - Get the temperature [K], density [kg/m³ or kmol/m³], and vapor fraction. - """ - def __get__(self): - return self.T, self.density, self.Q - - property TPQ: - """ - Get/Set the temperature [K], pressure [Pa], and vapor fraction of a - PureFluid. - - An Exception is raised if the thermodynamic state is not consistent. - """ - def __get__(self): - return self.T, self.P, self.Q - def __set__(self, values): - T = values[0] if values[0] is not None else self.T - P = values[1] if values[1] is not None else self.P - Q = values[2] if values[2] is not None else self.Q - self.thermo.setState_TPQ(T, P, Q) - - property UVQ: - """ - Get the internal energy [J/kg or J/kmol], specific volume - [m³/kg or m³/kmol], and vapor fraction. - """ - def __get__(self): - return self.u, self.v, self.Q - - property DPQ: - """Get the density [kg/m³], pressure [Pa], and vapor fraction.""" - def __get__(self): - return self.density, self.P, self.Q - - property HPQ: - """ - Get the enthalpy [J/kg or J/kmol], pressure [Pa] and vapor fraction. - """ - def __get__(self): - return self.h, self.P, self.Q - - property SPQ: - """ - Get the entropy [J/kg/K or J/kmol/K], pressure [Pa], and vapor fraction. - """ - def __get__(self): - return self.s, self.P, self.Q - - property SVQ: - """ - Get the entropy [J/kg/K or J/kmol/K], specific volume [m³/kg or - m³/kmol], and vapor fraction. - """ - def __get__(self): - return self.s, self.v, self.Q - - -class Element: - """ - An element or a named isotope defined in Cantera. - - Class `Element` gets data for the elements and isotopes defined in - :ct:`Elements.cpp`. This class can be used in two ways. The - first way is to get information about all of the elements stored in - Cantera. The three attributes `num_elements_defined`, - `element_symbols`, and `element_names` can be accessed by:: - - >>> ct.Element.num_elements_defined - >>> ct.Element.element_symbols - >>> ct.Element.element_names - - Otherwise, if the class `Element` is called with an argument, it - stores the data about that particular element. For example:: - - >>> ar_sym = ct.Element('Ar') - >>> ar_name = ct.Element('argon') - >>> ar_num = ct.Element(18) - - would all create instances with the information for argon. The - available argument options to create an instance of the `Element` - class with the element information are the `name`, `symbol`, and - `atomic_number`. Once an instance of the class is made, the `name`, - `atomic_number`, `symbol`, and atomic `weight` can be accessed as - attributes of the instance of the `Element` class. - - >>> ar_sym.name - 'argon' - >>> ar_sym.weight - 39.948 - >>> ar_sym.atomic_number - 18 - >>> ar_sym.symbol - 'Ar' - - The elements available are listed below, in the `element_symbols` - and `element_names` attribute documentation. - """ - - #: The number of named elements (not isotopes) defined in Cantera - num_elements_defined = numElementsDefined() - - #: A list of the symbols of all the elements (not isotopes) defined in Cantera - element_symbols = tuple(pystr(s) for s in elementSymbols()) - - #: A list of the names of all the elements (not isotopes) defined in Cantera - element_names = tuple(pystr(n) for n in elementNames()) - - def __init__(self, arg): - if isinstance(arg, (str, bytes)): - try: - # Assume the argument is the element symbol and try to get the name - self._name = pystr(getElementName(stringify(arg))) - except CanteraError: - # If getting the name failed, the argument must be the name - self._symbol = pystr(getElementSymbol(stringify(arg))) - self._name = arg.lower() - else: - self._symbol = arg - - self._atomic_number = getAtomicNumber(stringify(arg)) - self._weight = getElementWeight(stringify(arg)) - elif isinstance(arg, int): - self._atomic_number = arg - self._name = pystr(getElementName(arg)) - self._symbol = pystr(getElementSymbol(arg)) - self._weight = getElementWeight(arg) - else: - raise TypeError('The input argument to Element must be a string ' - 'or an integer') - - @property - def name(self): - """The name of the element or isotope.""" - return self._name - - @property - def atomic_number(self): - """The atomic number of the element or isotope.""" - return self._atomic_number - - @property - def symbol(self): - """The symbol of the element or isotope.""" - return self._symbol - - @property - def weight(self): - """The atomic weight of the element or isotope.""" - return self._weight +# This file is part of Cantera. See License.txt in the top-level directory or +# at https://cantera.org/license.txt for license and copyright information. + +import warnings +import weakref as _weakref +import numbers as _numbers +import numpy as np +cimport numpy as np + +from .speciesthermo cimport * +from .ctcxx cimport span +from .kinetics cimport CxxKinetics +from .transport cimport * +from ._utils cimport * +from ._utils import CanteraError +from .units cimport * + +cdef enum ThermoBasisType: + mass_basis = 0 + molar_basis = 1 + +ctypedef CxxPlasmaPhase* CxxPlasmaPhasePtr +ctypedef CxxSurfPhase* CxxSurfPhasePtr + +class ThermoModelMethodError(Exception): + """Exception raised for an invalid method used by a thermo model + + :param thermo_model: + The thermo model used by class `ThermoPhase` + + """ + + def __init__(self, thermo_model): + self.thermo_model = thermo_model + super().__init__(f"This method is invalid for {self.thermo_model}") + + +cdef class Species: + """ + A class which stores data about a single chemical species that may be + needed to add it to a `Solution` or `Interface` object (and to the + underlying `ThermoPhase` and `Transport` objects). + + :param name: + A string giving the name of the species, such as ``'CH4'`` + :param composition: + The elemental composition of the species, given either as a dict or a + composition string, such as ``{'C':1, 'H':4}`` or ``'C:1, H:4'``. + :param charge: + The electrical charge, in units of the elementary charge. Default 0.0. + :param size: + The effective size [m] of the species. Default 1.0. + :param init: + Used internally when wrapping :ct:`Species` objects returned from C++ + + Example: creating an ideal gas phase with a single species:: + + ch4 = ct.Species('CH4', 'C:1, H:4') + ch4.thermo = ct.ConstantCp(300, 1000, 101325, + (300, -7.453347e7, 1.865912e5, 3.576053e4)) + tran = ct.GasTransportData() + tran.set_customary_units('nonlinear', 3.75, 141.40, 0.0, 2.60, 13.00) + ch4.transport = tran + gas = ct.Solution(thermo='ideal-gas', species=[ch4]) + + The static methods `list_from_file` and `list_from_yaml` can be used to create + `Species` objects from existing definitions in the YAML format. Either of the + following will produce a list of 53 `Species` objects containing the species defined + in the GRI 3.0 mechanism:: + + S = ct.Species.list_from_file("gri30.yaml") + + import pathlib + S = ct.Species.list_from_yaml( + pathlib.Path('path/to/gri30.yaml').read_text(), + section='species') + + """ + def __cinit__(self, *args, init=True, **kwargs): + if init: + self._species.reset(new CxxSpecies()) + self.species = self._species.get() + + def __init__(self, name=None, composition=None, charge=None, size=None, + *args, init=True, **kwargs): + if not init: + return + + if name is not None: + self.species.name = stringify(name) + + if composition is not None: + self.species.composition = comp_map(composition) + + if charge is not None: + self.species.charge = charge + + if size is not None: + self.species.size = size + + cdef _assign(self, shared_ptr[CxxSpecies] other): + self._species = other + self.species = self._species.get() + + @staticmethod + def from_yaml(text): + """ + Create a `Species` object from its YAML string representation. + """ + cxx_species = CxxNewSpecies(AnyMapFromYamlString(stringify(text))) + species = Species(init=False) + species._assign(cxx_species) + return species + + @staticmethod + def from_dict(data): + """ + Create a `Species` object from a dictionary corresponding to its YAML + representation. + + :param data: + A dictionary corresponding to the YAML representation. + """ + cdef CxxAnyMap any_map = py_to_anymap(data) + cxx_species = CxxNewSpecies(any_map) + species = Species(init=False) + species._assign(cxx_species) + return species + + @staticmethod + def list_from_file(filename, section="species"): + """ + Create a list of Species objects from all of the species defined in the section + ``section`` of a YAML file. Directories on Cantera's input file path will be + searched for the specified file. + """ + root = AnyMapFromYamlFile(stringify(filename)) + species = [] + for a in CxxGetSpecies(root[stringify(section)]): + b = Species(init=False) + b._assign(a) + species.append(b) + return species + + @staticmethod + def list_from_yaml(text, section=None): + """ + Create a list of Species objects from all the species defined in a YAML + string. If ``text`` is a YAML mapping, the ``section`` name of the list + to be read must be specified. If ``text`` is a YAML list, no ``section`` + name should be supplied. + """ + root = AnyMapFromYamlString(stringify(text)) + + # ``items`` is the pseudo-key used to access a list when it is at the + # top level of a YAML document + cxx_species = CxxGetSpecies(root[stringify(section or "items")]) + species = [] + for a in cxx_species: + b = Species(init=False) + b._assign(a) + species.append(b) + return species + + property name: + """ The name of the species. """ + def __get__(self): + return pystr(self.species.name) + + property composition: + """ + A dict containing the elemental composition of the species. Keys are + element names; values are the corresponding atomicities. + """ + def __get__(self): + return comp_map_to_dict(self.species.composition) + + property charge: + """ + The electrical charge on the species, in units of the elementary charge. + """ + def __get__(self): + return self.species.charge + + property size: + """ The effective size [m] of the species. """ + def __get__(self): + return self.species.size + + property molecular_weight: + """The molecular weight [amu] of the species. + + .. versionadded:: 3.0 + """ + def __get__(self): + return self.species.molecularWeight() + + property thermo: + """ + Get/Set the species reference-state thermodynamic data, as an instance + of class `SpeciesThermo`. + """ + def __get__(self): + if self.species.thermo.get() != NULL: + return wrapSpeciesThermo(self.species.thermo) + else: + return None + + def __set__(self, SpeciesThermo spthermo): + self.species.thermo = spthermo._spthermo + + property transport: + """ + Get/Set the species transport parameters, as an instance of class + `GasTransportData`. + """ + def __get__(self): + if self.species.transport.get() != NULL: + data = GasTransportData(init=False) + data._assign(self.species.transport) + return data + else: + return None + def __set__(self, GasTransportData tran): + self.species.transport = tran._data + + property input_data: + """ + Get input data defining this Species, along with any user-specified data + provided with its input (YAML) definition. + """ + def __get__(self): + cdef CxxThermoPhase* phase = self._phase.thermo if self._phase else NULL + return anymap_to_py(self.species.parameters(phase)) + + def update_user_data(self, data): + """ + Add the contents of the provided `dict` as additional fields when generating + YAML phase definition files with `Solution.write_yaml` or in the data returned + by `input_data`. Existing keys with matching names are overwritten. + """ + self.species.input.update(py_to_anymap(data), False) + + def clear_user_data(self): + """ + Clear all saved input data, so that the data given by `input_data` or + `Solution.write_yaml` will only include values generated by Cantera based on + the current object state. + """ + self.species.input.clear() + + def __repr__(self): + return ''.format(self.name) + + +cdef class ThermoPhase(_SolutionBase): + """ + A phase with an equation of state. + + Class `ThermoPhase` may be used to represent the intensive thermodynamic + state of a phase of matter, which might be a gas, liquid, or solid. + + Class `ThermoPhase` is not usually instantiated directly. It is used + as a base class for classes `Solution` and `Interface`. + """ + + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + if 'source' not in kwargs: + self.thermo_basis = mass_basis + # In composite objects, the ThermoPhase constructor needs to be called first + # to prevent instantiation of stand-alone 'Kinetics' or 'Transport' objects. + # The following is used as a sentinel. After initialization, the _references + # object is used to track whether the ThermoPhase is being used by other + # objects that require the number of species to remain constant and do not + # have C++ implementations (e.g. Quantity objects). + self._references = _weakref.WeakKeyDictionary() + # validate plasma phase + self._enable_plasma = False + if dynamic_cast[CxxPlasmaPhasePtr](self.thermo): + self._enable_plasma = True + self.plasma = self.thermo + + property thermo_model: + """ + Return thermodynamic model describing phase. + """ + def __get__(self): + return pystr(self.thermo.type()) + + property phase_of_matter: + """ + Get the thermodynamic phase (gas, liquid, etc.) at the current conditions. + """ + def __get__(self): + return pystr(self.thermo.phaseOfMatter()) + + def report(self, show_thermo=True, float threshold=1e-14): + """ + Generate a report describing the thermodynamic state of this phase. To + print the report to the terminal, simply call the phase object. The + following two statements are equivalent:: + + >>> phase() + >>> print(phase.report()) + + :param show_thermo: + A Boolean argument specifying whether to show phase thermodynamic + information in the output. + :param threshold: + The threshold used to clip data in the output. Values below the threshold + are not displayed. + """ + return pystr(self.thermo.report(bool(show_thermo), threshold)) + + def __call__(self, *args, **kwargs): + print(self.report(*args, **kwargs)) + + property is_pure: + """ + Returns true if the phase represents a pure (fixed composition) substance + """ + def __get__(self): + return self.thermo.isPure() + + property has_phase_transition: + """ + Returns true if the phase represents a substance with phase transitions + """ + def __get__(self): + return self.thermo.hasPhaseTransition() + + property is_compressible: + """ + Returns true if the density of the phase is an independent variable defining + the thermodynamic state of a substance + """ + def __get__(self): + return self.thermo.isCompressible() + + @property + def _native_mode(self): + """ Return string acronym representing native state """ + return pystr(self.thermo.nativeMode()) + + property _native_state: + """ + Default properties defining a state + """ + def __get__(self): + cdef pair[string, size_t] item + native = {pystr(item.first): item.second for item in self.thermo.nativeState()} + return tuple([i for i, j in sorted(native.items(), key=lambda kv: kv[1])]) + + property _full_states: + """ + Sets of parameters which set the full thermodynamic state + """ + def __get__(self): + states = self.thermo.fullStates() + states = [pystr(s) for s in states] + return {frozenset(k): k for k in states} + + property _partial_states: + """ + Sets of parameters which set a valid partial thermodynamic state + """ + def __get__(self): + states = self.thermo.partialStates() + states = [pystr(s) for s in states] + return {frozenset(k): k for k in states} + + property basis: + """ + Determines whether intensive thermodynamic properties are treated on a + ``mass`` (per kg) or ``molar`` (per kmol) basis. This affects the values + returned by the properties `h`, `u`, `s`, `g`, `v`, `density`, `cv`, + and `cp`, as well as the values used with the state-setting properties + such as `HPX` and `UV`. + """ + def __get__(self): + if self.thermo_basis == mass_basis: + return 'mass' + else: + return 'molar' + + def __set__(self, value): + if value == 'mass': + self.thermo_basis = mass_basis + elif value == 'molar': + self.thermo_basis = molar_basis + else: + raise ValueError("Valid choices are 'mass' or 'molar'." + " Got {!r}.".format(value)) + + cdef double _mass_factor(self): + """ Conversion factor from current basis to kg """ + if self.thermo_basis == molar_basis: + return self.thermo.meanMolecularWeight() + else: + return 1.0 + + cdef double _mole_factor(self): + """ Conversion factor from current basis to moles """ + if self.thermo_basis == mass_basis: + return 1.0/self.thermo.meanMolecularWeight() + else: + return 1.0 + + def equilibrate(self, XY, solver='auto', double rtol=1e-9, + int max_steps=1000, int max_iter=100, int estimate_equil=0, + int log_level=0): + """ + Set to a state of chemical equilibrium holding property pair + ``XY`` constant. + + :param XY: + A two-letter string, which must be one of the set:: + + ['TP','TV','HP','SP','SV','UV'] + + :param solver: + Specifies the equilibrium solver to use. May be one of the following: + + * ``'element_potential'`` - a fast solver using the element potential + method + * ``'gibbs'`` - a slower but more robust Gibbs minimization solver + * ``'vcs'`` - the VCS non-ideal equilibrium solver + * ``'auto'`` - The element potential solver will be tried first, then + if it fails the Gibbs solver will be tried. + :param rtol: + The relative error tolerance. + :param max_steps: + The maximum number of steps in composition to take to find a converged + solution. + :param max_iter: + For the Gibbs minimization solver, this specifies the number of + outer iterations on T or P when some property pair other + than TP is specified. + :param estimate_equil: + Integer indicating whether the solver should estimate its own + initial condition. If 0, the initial mole fraction vector in the + `ThermoPhase` object is used as the initial condition. If 1, the + initial mole fraction vector is used if the element abundances are + satisfied. If -1, the initial mole fraction vector is thrown out, + and an estimate is formulated. + :param log_level: + Set to a value greater than 0 to write diagnostic output. + """ + self.thermo.equilibrate(stringify(XY.upper()), stringify(solver), rtol, + max_steps, max_iter, estimate_equil, log_level) + + ####### Composition, species, and elements ######## + + property n_elements: + """Number of elements.""" + def __get__(self): + return self.thermo.nElements() + + cpdef int element_index(self, element) except *: + """ + The index of element ``element``, which may be specified as a string or + an integer. In the latter case, the index is checked for validity and + returned. If no such element is present, an exception is thrown. + """ + if isinstance(element, (str, bytes)): + return self.thermo.elementIndex(stringify(element), True) + if isinstance(element, (int, float)): + return self.thermo.checkElementIndex(element) + + raise TypeError("'element' must be a string or a number. " + f"Got {element!r}.") + + def element_name(self, m): + """Name of the element with index ``m``.""" + return pystr(self.thermo.elementName(m)) + + property element_names: + """A list of all the element names.""" + def __get__(self): + return [self.element_name(m) for m in range(self.n_elements)] + + def atomic_weight(self, m): + """Atomic weight [kg/kmol] of element ``m``""" + return self.thermo.atomicWeight(self.element_index(m)) + + property atomic_weights: + """Array of atomic weight [kg/kmol] for each element in the mixture.""" + def __get__(self): + return np.array([self.thermo.atomicWeight(m) for m in range(self.n_elements)]) + + property n_species: + """Number of species.""" + def __get__(self): + return self.thermo.nSpecies() + + property n_selected_species: + """ + Number of species selected for output (by slicing of Solution object) + """ + def __get__(self): + return self._selected_species.size or self.n_species + + def species_name(self, k): + """Name of the species with index ``k``.""" + return pystr(self.thermo.speciesName(k)) + + property species_names: + """A list of all the species names.""" + def __get__(self): + if self._selected_species.size: + indices = self._selected_species + else: + indices = range(self.n_species) + return [self.species_name(k) for k in indices] + + cpdef int species_index(self, species) except *: + """ + The index of species ``species``, which may be specified as a string or + an integer. In the latter case, the index is checked for validity and + returned. If no such species is present, an exception is thrown. + """ + if isinstance(species, (str, bytes)): + return self.thermo.speciesIndex(stringify(species), True) + if isinstance(species, (int, float)): + return self.thermo.checkSpeciesIndex(species) + raise TypeError("'species' must be a string or a number." + f" Got {species!r}.") + + property case_sensitive_species_names: + """Enforce case-sensitivity for look up of species names""" + def __get__(self): + return self.thermo.caseSensitiveSpecies() + def __set__(self, val): + self.thermo.setCaseSensitiveSpecies(bool(val)) + + def species(self, k=None): + """ + Return the `Species` object for species ``k``, where ``k`` is either the + species index or the species name. If ``k`` is not specified, a list of + all species objects is returned. Changes to this object do not affect + the `ThermoPhase` or `Solution` object until the `modify_species` + function is called. + """ + if k is None: + return [self.species(i) for i in range(self.n_species)] + + s = Species(init=False) + s._phase = self + if isinstance(k, (str, bytes)): + s._assign(self.thermo.species(stringify(k))) + elif isinstance(k, (int, float)): + s._assign(self.thermo.species(k)) + else: + raise TypeError("Argument must be a string or a number." + " Got {!r}.".format(k)) + return s + + def modify_species(self, k, Species species): + """ + Modify the thermodynamic data associated with a species. The species name, + elemental composition, and type of thermo parameterization must be unchanged. + """ + self.thermo.modifySpecies(k, species._species) + if self.kinetics: + self.kinetics.invalidateCache() + + def add_species(self, Species species): + """ + Add a new species to this phase. Missing elements will be added + automatically. + """ + if self._references: + raise CanteraError('Cannot add species to ThermoPhase object because it' + ' is being used by a Quantity object.') + self.thermo.addUndefinedElements() + self.thermo.addSpecies(species._species) + species._phase = self + self.thermo.initThermo() + if self.kinetics: + self.kinetics.invalidateCache() + + def add_species_alias(self, name, alias): + """ + Add the alternate species name ``alias`` for an original species ``name``. + """ + self.thermo.addSpeciesAlias(stringify(name), stringify(alias)) + + def find_isomers(self, comp): + """ + Find species/isomers matching a composition specified by ``comp``. + """ + + if isinstance(comp, dict): + iso = self.thermo.findIsomers(comp_map(comp)) + elif isinstance(comp, (str, bytes)): + iso = self.thermo.findIsomers(stringify(comp)) + else: + raise CanteraError('Invalid composition') + + return [pystr(b) for b in iso] + + def n_atoms(self, species, element): + """ + Number of atoms of element ``element`` in species ``species``. The element + and species may be specified by name or by index. + + >>> phase.n_atoms('CH4','H') + 4 + """ + return self.thermo.nAtoms(self.species_index(species), + self.element_index(element)) + + cdef np.ndarray _getArray1(self, thermoMethod1d method): + cdef np.ndarray[np.double_t, ndim=1] data = np.empty(self.n_species) + method(self.thermo, span[double](&data[0], data.size)) + if self._selected_species.size: + return data[self._selected_species] + else: + return data + + cdef void _setArray1(self, thermoMethod1d method, values) except *: + cdef np.ndarray[np.double_t, ndim=1] data + cdef span[double] view + + values = np.squeeze(values) + if values.ndim == 0: + values = values[np.newaxis] # corner case for single-species phases + + if len(values) == self.n_species: + data = np.ascontiguousarray(values, dtype=np.double) + elif len(values) == len(self._selected_species) != 0: + data = np.zeros(self.n_species, dtype=np.double) + for i,k in enumerate(self._selected_species): + data[k] = values[i] + else: + msg = "Got {}. Expected {}".format(len(values), self.n_species) + if len(self._selected_species): + msg += ' or {}'.format(len(self._selected_species)) + raise ValueError('Array has incorrect length. ' + msg + '.') + method(self.thermo, span[double](&data[0], data.size)) + + property molecular_weights: + """Array of species molecular weights (molar masses) [kg/kmol].""" + def __get__(self): + return self._getArray1(thermo_getMolecularWeights) + + property charges: + """Array of species charges [elem. charge].""" + def __get__(self): + return self._getArray1(thermo_getCharges) + + property mean_molecular_weight: + """The mean molecular weight (molar mass) [kg/kmol].""" + def __get__(self): + return self.thermo.meanMolecularWeight() + + property Y: + """ + Get/Set the species mass fractions. Can be set as an array, as a dictionary, + or as a string. Always returns an array:: + + >>> phase.Y = [0.1, 0, 0, 0.4, 0, 0, 0, 0, 0.5] + >>> phase.Y = {'H2':0.1, 'O2':0.4, 'AR':0.5} + >>> phase.Y = 'H2:0.1, O2:0.4, AR:0.5' + >>> phase.Y + array([0.1, 0, 0, 0.4, 0, 0, 0, 0, 0.5]) + """ + def __get__(self): + return self._getArray1(thermo_getMassFractions) + def __set__(self, Y): + if isinstance(Y, (str, bytes)): + self.thermo.setMassFractionsByName(stringify(Y)) + elif isinstance(Y, dict): + self.thermo.setMassFractionsByName(comp_map(Y)) + else: + self._setArray1(thermo_setMassFractions, Y) + + property X: + """ + Get/Set the species mole fractions. Can be set as an array, as a dictionary, + or as a string. Always returns an array:: + + >>> phase.X = [0.1, 0, 0, 0.4, 0, 0, 0, 0, 0.5] + >>> phase.X = {'H2':0.1, 'O2':0.4, 'AR':0.5} + >>> phase.X = 'H2:0.1, O2:0.4, AR:0.5' + >>> phase.X + array([0.1, 0, 0, 0.4, 0, 0, 0, 0, 0.5]) + """ + def __get__(self): + return self._getArray1(thermo_getMoleFractions) + def __set__(self, X): + if isinstance(X, (str, bytes)): + self.thermo.setMoleFractionsByName(stringify(X)) + elif isinstance(X, dict): + self.thermo.setMoleFractionsByName(comp_map(X)) + else: + self._setArray1(thermo_setMoleFractions, X) + + property concentrations: + """ + Get/Set the species concentrations. Units are kmol/m³ for bulk phases, kmol/m² + for surface phases, and kmol/m for edge phases. + """ + def __get__(self): + return self._getArray1(thermo_getConcentrations) + def __set__(self, C): + self._setArray1(thermo_setConcentrations, C) + + def __composition_to_array(self, comp, basis): + """take a mixture composition in mole or mass fraction as string, + dict or array and return array (for internal use)""" + if (isinstance(comp, str) and ':' not in comp + and comp in self.species_names): + comp += ':1.0' + + original_state = self.state + + if basis == 'mole': + self.TPX = None, None, comp + arr = np.copy(self.X) + elif basis == 'mass': + self.TPY = None, None, comp + arr = np.copy(self.Y) + else: + raise ValueError("basis must either be 'mass' or mole'.") + + self.state = original_state + return arr + + def set_equivalence_ratio(self, phi, fuel, oxidizer, basis="mole", *, diluent=None, + fraction=None): + """ + Set the composition to a mixture of ``fuel`` and ``oxidizer`` at the + specified equivalence ratio ``phi``, holding temperature and pressure + constant. Considers the oxidation of C to CO2, H to H2O and S to SO2. + Other elements are assumed not to participate in oxidation (that is, + N ends up as N2). The ``basis`` determines the fuel and oxidizer + compositions: ``basis='mole'`` means mole fractions (default), + ``basis='mass'`` means mass fractions. The fuel/oxidizer mixture can be + be diluted by a ``diluent`` based on a mixing ``fraction``. The amount of + diluent is quantified as a fraction of fuel, oxidizer or the fuel/oxidizer + mixture. For more information, see the :doc:`Python example + ` :: + + >>> gas.set_equivalence_ratio(0.5, 'CH4', 'O2:1.0, N2:3.76', basis='mole') + >>> gas.mass_fraction_dict() + {'CH4': 0.02837633052851, 'N2': 0.7452356312613, 'O2': 0.22638803821018} + >>> gas.set_equivalence_ratio(1.2, 'NH3:0.8,CO:0.2', 'O2:1', basis='mole') + >>> gas.mass_fraction_dict() + {'CO': 0.14784006249290, 'NH3': 0.35956645545401, 'O2': 0.49259348205308} + + :param phi: + Equivalence ratio + :param fuel: + Fuel species name or mole/mass fractions as string, array, or dict. + :param oxidizer: + Oxidizer species name or mole/mass fractions as a string, array, or dict. + :param basis: + Determines if ``fuel`` and ``oxidizer`` are given in mole + fractions (``basis='mole'``) or mass fractions (``basis='mass'``). + :param diluent: + Optional parameter. Required if dilution is used. Specifies the composition + of the diluent in mole/mass fractions as a string, array or dict. + :param fraction: + Optional parameter. Dilutes the fuel/oxidizer mixture with the diluent + according to ``fraction``. Fraction can refer to the fraction of diluent in + the mixture (for example ``fraction="diluent:0.7"`` will create a mixture + with 30 % fuel/oxidizer and 70 % diluent), the fraction of fuel in the + mixture (for example ``fraction="fuel:0.1"`` means that the mixture contains + 10 % fuel. The amount of oxidizer is determined from the equivalence ratio + and the remaining mixture is the diluent) or fraction of oxidizer in the + mixture (for example ``fraction="oxidizer:0.1"``). The fraction itself is + interpreted as mole or mass fraction based on ``basis``. The diluent is not + considered in the computation of the equivalence ratio. Default is no + dilution or ``fraction=None``. May be given as string or dictionary (for + example ``fraction={"fuel":0.7}``). + """ + cdef np.ndarray[np.double_t, ndim=1] fuel_comp = np.ascontiguousarray( + self.__composition_to_array(fuel, basis), dtype=np.double) + cdef np.ndarray[np.double_t, ndim=1] ox_comp = np.ascontiguousarray( + self.__composition_to_array(oxidizer, basis), dtype=np.double) + + cdef span[double] fuel_view = span[double](&fuel_comp[0], fuel_comp.size) + cdef span[double] ox_view = span[double](&ox_comp[0], ox_comp.size) + self.thermo.setEquivalenceRatio(phi, fuel_view, ox_view, + ThermoBasis.mass if basis == "mass" + else ThermoBasis.molar) + + if (fraction is None) != (diluent is None): + raise ValueError("If dilution is used, both 'fraction' and 'diluent' " + "parameters are required.") + + if fraction is None: + return + + if isinstance(fraction, str): + fraction_dict = comp_map_to_dict(parseCompString(stringify(fraction))) + elif isinstance(fraction, dict): + fraction_dict = fraction + else: + raise ValueError("The fraction argument must be given as string or " + "dictionary.") + + if len(fraction_dict) != 1: + raise ValueError("Invalid format for the fraction. Must be provided for " + "example as fraction='fuel:0.1'") + + fraction_type = list(fraction_dict.keys())[0] + fraction_value = float(list(fraction_dict.values())[0]) + + if fraction_value < 0 or fraction_value > 1: + raise ValueError("The fraction must be between 0 and 1") + + if fraction_type not in ["fuel", "oxidizer", "diluent"]: + raise ValueError("The fraction must specify 'fuel', 'oxidizer' or " + "'diluent'") + + cdef np.ndarray[np.double_t, ndim=1] diluent_comp = np.ascontiguousarray( + self.__composition_to_array(diluent, basis), dtype=np.double) + + # this function changes the composition and fixes temperature and pressure + T, P = self.T, self.P + # if 'fraction' is specified for diluent, just scale the mass or mole fractions + # of the fuel/oxidizer mixture accordingly + if fraction_type == "diluent": + if basis == "mole": + X_fuelox = self.X + self.X = diluent_comp + self.X = (1.0 - fraction_value) * X_fuelox + fraction_value * self.X + else: + Y_fuelox = self.Y + self.Y = diluent_comp + self.Y = (1.0 - fraction_value) * Y_fuelox + fraction_value * self.Y + self.TP = T, P + return + + # get the mixture fraction before scaling / diluent addition + Z_fuel = self.mixture_fraction(fuel, oxidizer, basis) + + if Z_fuel == 0.0 and fraction_type == "fuel": + raise ValueError("No fuel in the fuel/oxidizer mixture") + + if Z_fuel == 1.0 and fraction_type == "oxidizer": + raise ValueError("No oxidizer in the fuel/oxidizer mixture") + + if basis == "mass": # for mass basis, it is straight forward + if fraction_type == "fuel": + Z = Z_fuel + else: # oxidizer + Z = 1.0 - Z_fuel + if fraction_value > Z: + raise ValueError(f"The {fraction_type} fraction after dilution cannot " + "be higher than {fraction_type} fraction in the " + "original mixture.") + Y_mix = self.Y + self.Y = diluent_comp + factor = fraction_value / Z + self.Y = factor*Y_mix + (1.0 - factor) * self.Y + else: + # convert mass based mixture fraction to molar one, Z = kg fuel / kg mixture + X_mix = self.X + M_mix = self.mean_molecular_weight + self.X = fuel_comp + M_fuel = self.mean_molecular_weight + Z_fuel_mole = Z_fuel * M_mix / M_fuel # mol fuel / mol mixture + if fraction_type == "fuel": + Z = Z_fuel_mole + else: # oxidizer + Z = 1.0 - Z_fuel_mole + if fraction_value > Z: + raise ValueError(f"The {fraction_type} fuel or oxidizer fraction after " + "dilution cannot be higher than {fraction_type} " + "fraction in the original mixture.") + self.X = diluent_comp + factor = fraction_value / Z + self.X = factor * X_mix + (1.0 - factor) * self.X + self.TP = T, P + + def set_mixture_fraction(self, mixture_fraction, fuel, oxidizer, basis='mole'): + """ + Set the composition to a mixture of ``fuel`` and ``oxidizer`` at the + specified mixture fraction ``mixture_fraction`` (kg fuel / kg mixture), holding + temperature and pressure constant. Considers the oxidation of C to CO2, + H to H2O and S to SO2. Other elements are assumed not to participate in + oxidation (that is, N ends up as N2). The ``basis`` determines the composition + of fuel and oxidizer: ``basis='mole'`` (default) means mole fractions, + ``basis='mass'`` means mass fractions. For more information, see the + :doc:`Python example ` :: + + >>> gas.set_mixture_fraction(0.5, 'CH4', 'O2:1.0, N2:3.76') + >>> gas.mass_fraction_dict() + {'CH4': 0.5, 'N2': 0.38350014242997776, 'O2': 0.11649985757002226} + >>> gas.set_mixture_fraction(0.5, {'NH3':0.8, 'CO':0.2}, 'O2:1.0') + >>> gas.mass_fraction_dict() + {'CO': 0.145682068778996, 'NH3': 0.354317931221004, 'O2': 0.5} + + :param mixture_fraction: + Mixture fraction (kg fuel / kg mixture) + :param fuel: + Fuel species name or mole/mass fractions as string, array, or dict. + :param oxidizer: + Oxidizer species name or mole/mass fractions as a string, array, or + dict. + :param basis: determines if ``fuel`` and ``oxidizer`` are given in mole + fractions (``basis='mole'``) or mass fractions (``basis='mass'``) + """ + cdef np.ndarray[np.double_t, ndim=1] f = \ + np.ascontiguousarray(self.__composition_to_array(fuel, basis), dtype=np.double) + cdef np.ndarray[np.double_t, ndim=1] o = \ + np.ascontiguousarray(self.__composition_to_array(oxidizer, basis), dtype=np.double) + + cdef span[double] fuel_view = span[double](&f[0], f.size) + cdef span[double] ox_view = span[double](&o[0], o.size) + self.thermo.setMixtureFraction(mixture_fraction, fuel_view, ox_view, + ThermoBasis.mass if basis == 'mass' else ThermoBasis.molar) + + def equivalence_ratio(self, fuel=None, oxidizer=None, basis="mole", + include_species=None): + r""" + Get the equivalence ratio :math:`\phi` of the current mixture, which is a + conserved quantity. Considers the oxidation of C to CO2, H to H2O + and S to SO2. Other elements are assumed not to participate in oxidation + (that is, N ends up as N2). If fuel and oxidizer are not specified, the + equivalence ratio is computed from the available oxygen and the + required oxygen for complete oxidation: + + .. math:: \phi = \frac{Z_{\mathrm{mole},C} + Z_{\mathrm{mole},S} + + \frac{1}{4}Z_{\mathrm{mole},H}} {\frac{1}{2}Z_{\mathrm{mole},O}} + + where :math:`Z_{\mathrm{mole},e}` is the elemental mole fraction of element + :math:`e`. If the fuel and oxidizer compositions are specified, :math:`\phi` is + computed from: + + .. math:: \phi = \frac{Z}{1-Z}\frac{1-Z_{\mathrm{st}}}{Z_{\mathrm{st}}} + + where :math:`Z` is the Bilger mixture fraction and :math:`Z_{\mathrm{st}}` + the Bilger mixture fraction at stoichiometric conditions. + The ``basis`` determines the composition of fuel and oxidizer: + ``basis='mole'`` (default) means mole fractions, ``basis='mass'`` means + mass fractions. Note that this definition takes all species into account. + In case certain species like inert diluents should be ignored, a + list of species can be provided with ``include_species``. This means that + only these species are considered for the computation of the equivalence + ratio. For more information, see the + :doc:`Python example ` :: + + >>> gas.set_equivalence_ratio(0.5, fuel='CH3:0.5, CH3OH:.5, N2:0.125', oxidizer='O2:0.21, N2:0.79, NO:0.01') + >>> gas.equivalence_ratio(fuel='CH3:0.5, CH3OH:.5, N2:0.125', oxidizer='O2:0.21, N2:0.79, NO:0.01') + 0.5 + + :param fuel: + Fuel species name or mole/mass fractions as string, array, or dict. + :param oxidizer: + Oxidizer species name or mole/mass fractions as a string, array, or dict. + :param basis: + Determines if ``fuel`` and ``oxidizer`` are given in mole fractions + (``basis="mole"``) or mass fractions (``basis="mass"``) + :param include_species: + List of species names (optional). Only these species are considered for the + computation of the equivalence ratio. By default, all species are considered + """ + if include_species is not None: + # remove unwanted species temporarily + Y = np.zeros(self.n_species) + indices = [self.species_index(s) for s in include_species] + for k in indices: + Y[k] = self.Y[k] + T_orig, P_orig, Y_orig = self.T, self.P, self.Y + self.Y = Y + + if fuel is None and oxidizer is None: + phi = self.thermo.equivalenceRatio() + if include_species is not None: + self.TPY = T_orig, P_orig, Y_orig + return phi + + cdef np.ndarray[np.double_t, ndim=1] f = np.ascontiguousarray( + self.__composition_to_array(fuel, basis), dtype=np.double) + cdef np.ndarray[np.double_t, ndim=1] o = np.ascontiguousarray( + self.__composition_to_array(oxidizer, basis), dtype=np.double) + + cdef span[double] fuel_view = span[double](&f[0], f.size) + cdef span[double] ox_view = span[double](&o[0], o.size) + phi = self.thermo.equivalenceRatio(fuel_view, ox_view, + ThermoBasis.mass if basis == "mass" else ThermoBasis.molar) + if include_species is not None: + self.TPY = T_orig, P_orig, Y_orig + return phi + + def mixture_fraction(self, fuel, oxidizer, basis='mole', element="Bilger"): + r""" + Get the mixture fraction of the current mixture in + (kg fuel / (kg oxidizer + kg fuel)). This is a quantity that is conserved after + oxidation. Considers the oxidation of C to CO2, H to H2O and S to SO2. Other + elements are assumed not to participate in oxidation (that is, N ends up as N2). + The ``basis`` determines the composition of fuel and oxidizer: + ``basis="mole"`` (default) means mole fractions, ``basis="mass"`` means mass + fractions. The mixture fraction can be computed from a single element (for + example, carbon with ``element="C"``) + + .. math:: Z_m = \frac{Z_{\mathrm{mass},m}-Z_{\mathrm{mass},m,\mathrm{ox}}} + {Z_{\mathrm{mass},\mathrm{fuel}}-Z_{\mathrm{mass},m,\mathrm{ox}}} + + where :math:`Z_{\mathrm{mass},m}` is the elemental mass fraction of + element :math:`m` in the mixture, and :math:`Z_{\mathrm{mass},m,\mathrm{ox}}` + and :math:`Z_{\mathrm{mass},\mathrm{fuel}}` are the elemental mass fractions of + the oxidizer and fuel, or from the Bilger mixture fraction + (``element="Bilger"``), which considers the elements C, S, H and O + :cite:p:`bilger1979`. The Bilger mixture fraction is computed by default: + + .. math:: Z_m = Z_{\mathrm{Bilger}} = \frac{\beta-\beta_{\mathrm{ox}}} + {\beta_{\mathrm{fuel}}-\beta_{\mathrm{ox}}} + + with + + .. math:: \beta = 2\frac{Z_C}{M_C}+2\frac{Z_S}{M_S}+\frac{1}{2}\frac{Z_H}{M_H} + - \frac{Z_O}{M_O} + + and :math:`M_m` the atomic weight of element :math:`m`. + For more information, see the + :doc:`Python example ` :: + + >>> gas.set_mixture_fraction(0.5, 'CH3:0.5, CH3OH:0.5, N2:0.125', 'O2:0.21, N2:0.79, NO:0.01') + >>> gas.mixture_fraction('CH3:0.5, CH3OH:0.5, N2:0.125', 'O2:0.21, N2:0.79, NO:.01') + 0.5 + + :param fuel: + Fuel species name or mole/mass fractions as string, array, or dict. + :param oxidizer: + Oxidizer species name or mole/mass fractions as a string, array, or + dict. + :param basis: + Determines if ``fuel`` and ``oxidizer`` are given in mole + fractions (``basis='mole'``) or mass fractions (``basis='mass'``) + :param element: + Computes the mixture fraction from the specified elemental + mass fraction (given by element name or element index) or as + the Bilger mixture fraction (default) + """ + cdef np.ndarray[np.double_t, ndim=1] f = \ + np.ascontiguousarray(self.__composition_to_array(fuel, basis), dtype=np.double) + cdef np.ndarray[np.double_t, ndim=1] o = \ + np.ascontiguousarray(self.__composition_to_array(oxidizer, basis), dtype=np.double) + + if isinstance(element, (str, bytes)): + e_name = element + else: + e_name = self.element_name(self.element_index(element)) + + cdef span[double] fuel_view = span[double](&f[0], f.size) + cdef span[double] ox_view = span[double](&o[0], o.size) + return self.thermo.mixtureFraction( + fuel_view, ox_view, + ThermoBasis.mass if basis=='mass' else ThermoBasis.molar, + stringify(e_name)) + + def stoich_air_fuel_ratio(self, fuel, oxidizer, basis='mole'): + """ + Get the stoichiometric air to fuel ratio (kg oxidizer / kg fuel). Considers the + oxidation of C to CO2, H to H2O and S to SO2. Other elements are assumed + not to participate in oxidation (that is, N ends up as N2). + The ``basis`` determines the composition of fuel and oxidizer: ``basis='mole'`` (default) + means mole fractions, ``basis='mass'`` means mass fractions:: + + >>> gas.set_mixture_fraction(0.5, 'CH3:0.5, CH3OH:.5, N2:0.125', 'O2:0.21, N2:0.79, NO:0.01') + >>> gas.stoich_air_fuel_ratio('CH3:0.5, CH3OH:.5, N2:0.125', 'O2:0.21, N2:0.79, NO:0.01') + 8.148040722239438 + + :param fuel: + Fuel species name or mole/mass fractions as string, array, or dict. + :param oxidizer: + Oxidizer species name or mole/mass fractions as a string, array, or + dict. + :param basis: + Determines if ``fuel`` and ``oxidizer`` are given in mole + fractions (``basis='mole'``) or mass fractions (``basis='mass'``) + + """ + cdef np.ndarray[np.double_t, ndim=1] f = \ + np.ascontiguousarray(self.__composition_to_array(fuel, basis), dtype=np.double) + cdef np.ndarray[np.double_t, ndim=1] o = \ + np.ascontiguousarray(self.__composition_to_array(oxidizer, basis), dtype=np.double) + + cdef span[double] fuel_view = span[double](&f[0], f.size) + cdef span[double] ox_view = span[double](&o[0], o.size) + return self.thermo.stoichAirFuelRatio(fuel_view, ox_view, + ThermoBasis.mass if basis=='mass' else ThermoBasis.molar) + + def elemental_mass_fraction(self, m): + r""" + Get the elemental mass fraction :math:`Z_{\mathrm{mass},m}` of element + :math:`m` as defined by: + + .. math:: Z_{\mathrm{mass},m} = \sum_k \frac{a_{m,k} M_m}{M_k} Y_k + + with :math:`a_{m,k}` being the number of atoms of element :math:`m` in + species :math:`k`, :math:`M_m` the atomic weight of element :math:`m`, + :math:`M_k` the molecular weight of species :math:`k`, and :math:`Y_k` + the mass fraction of species :math:`k`:: + + >>> phase.elemental_mass_fraction('H') + 1.0 + + :param m: + Base element, may be specified by name or by index. + """ + return self.thermo.elementalMassFraction(self.element_index(m)) + + def elemental_mole_fraction(self, m): + r""" + Get the elemental mole fraction :math:`Z_{\mathrm{mole},m}` of element + :math:`m` (the number of atoms of element m divided by the total number + of atoms) as defined by: + + .. math:: Z_{\mathrm{mole},m} = \frac{\sum_k a_{m,k} X_k} + {\sum_k \sum_j a_{j,k} X_k} + + with :math:`a_{m,k}` being the number of atoms of element :math:`m` in + species :math:`k`, :math:`\sum_j` being a sum over all elements, and + :math:`X_k` being the mole fraction of species :math:`k`:: + + >>> phase.elemental_mole_fraction('H') + 1.0 + + :param m: + Base element, may be specified by name or by index. + """ + return self.thermo.elementalMoleFraction(self.element_index(m)) + + def set_unnormalized_mass_fractions(self, Y): + """ + Set the mass fractions without normalizing to force ``sum(Y) == 1.0``. + Useful primarily when calculating derivatives with respect to ``Y[k]`` by + finite difference. + """ + cdef np.ndarray[np.double_t, ndim=1] data + if len(Y) == self.n_species: + data = np.ascontiguousarray(Y, dtype=np.double) + else: + raise ValueError("Array has incorrect length." + " Got {}, expected {}.".format(len(Y), self.n_species)) + self.thermo.setMassFractions_NoNorm(span[double](&data[0], data.size)) + + def set_unnormalized_mole_fractions(self, X): + """ + Set the mole fractions without normalizing to force ``sum(X) == 1.0``. + Useful primarily when calculating derivatives with respect to ``X[k]`` + by finite difference. + """ + cdef np.ndarray[np.double_t, ndim=1] data + if len(X) == self.n_species: + data = np.ascontiguousarray(X, dtype=np.double) + else: + raise ValueError("Array has incorrect length." + " Got {}, expected {}.".format(len(X), self.n_species)) + self.thermo.setMoleFractions_NoNorm(span[double](&data[0], data.size)) + + def mass_fraction_dict(self, double threshold=0.0): + """ + Return a dictionary giving the mass fraction for each species by name where the + mass fraction is greater than ``threshold``. + """ + cdef pair[string,double] item + Y = self.thermo.getMassFractionsByName(threshold) + return {pystr(item.first):item.second for item in Y} + + def mole_fraction_dict(self, double threshold=0.0): + """ + Return a dictionary giving the mole fraction for each species by name where the + mole fraction is greater than ``threshold``. + """ + cdef pair[string,double] item + X = self.thermo.getMoleFractionsByName(threshold) + return {pystr(item.first):item.second for item in X} + + + ######## Read-only thermodynamic properties ######## + + property P: + """Pressure [Pa].""" + def __get__(self): + return self.thermo.pressure() + + property T: + """Temperature [K].""" + def __get__(self): + return self.thermo.temperature() + + property density: + """Density [kg/m³ or kmol/m³] depending on `basis`.""" + def __get__(self): + return self.thermo.density() / self._mass_factor() + + property density_mass: + """(Mass) density [kg/m³].""" + def __get__(self): + return self.thermo.density() + + property density_mole: + """Molar density [kmol/m³].""" + def __get__(self): + return self.thermo.molarDensity() + + property v: + """Specific volume [m³/kg or m³/kmol] depending on `basis`.""" + def __get__(self): + return self._mass_factor() / self.thermo.density() + + property volume_mass: + """Specific volume [m³/kg].""" + def __get__(self): + return 1.0 / self.thermo.density() + + property volume_mole: + """Molar volume [m³/kmol].""" + def __get__(self): + return self.thermo.molarVolume() + + property u: + """Internal energy in [J/kg or J/kmol].""" + def __get__(self): + return self.thermo.intEnergy_mole() * self._mole_factor() + + property int_energy_mole: + """Molar internal energy [J/kmol].""" + def __get__(self): + return self.thermo.intEnergy_mole() + + property int_energy_mass: + """Specific internal energy [J/kg].""" + def __get__(self): + return self.thermo.intEnergy_mass() + + property h: + """Enthalpy [J/kg or J/kmol] depending on `basis`.""" + def __get__(self): + return self.thermo.enthalpy_mole() * self._mole_factor() + + property enthalpy_mole: + """Molar enthalpy [J/kmol].""" + def __get__(self): + return self.thermo.enthalpy_mole() + + property enthalpy_mass: + """Specific enthalpy [J/kg].""" + def __get__(self): + return self.thermo.enthalpy_mass() + + property s: + """Entropy [J/kg/K or J/kmol/K] depending on `basis`.""" + def __get__(self): + return self.thermo.entropy_mole() * self._mole_factor() + + property entropy_mole: + """Molar entropy [J/kmol/K].""" + def __get__(self): + return self.thermo.entropy_mole() + + property entropy_mass: + """Specific entropy [J/kg/K].""" + def __get__(self): + return self.thermo.entropy_mass() + + property g: + """Gibbs free energy [J/kg or J/kmol] depending on `basis`.""" + def __get__(self): + return self.thermo.gibbs_mole() * self._mole_factor() + + property gibbs_mole: + """Molar Gibbs free energy [J/kmol].""" + def __get__(self): + return self.thermo.gibbs_mole() + + property gibbs_mass: + """Specific Gibbs free energy [J/kg].""" + def __get__(self): + return self.thermo.gibbs_mass() + + property cv: + """ + Heat capacity at constant volume and composition + [J/kg/K or J/kmol/K depending on `basis`]. + """ + def __get__(self): + return self.thermo.cv_mole() * self._mole_factor() + + property cv_mole: + """Molar heat capacity at constant volume and composition [J/kmol/K].""" + def __get__(self): + return self.thermo.cv_mole() + + property cv_mass: + """Specific heat capacity at constant volume and composition [J/kg/K].""" + def __get__(self): + return self.thermo.cv_mass() + + property cp: + """ + Heat capacity at constant pressure and composition. + [J/kg/K or J/kmol/K depending on `basis`] + """ + def __get__(self): + return self.thermo.cp_mole() * self._mole_factor() + + property cp_mole: + """Molar heat capacity at constant pressure and composition [J/kmol/K].""" + def __get__(self): + return self.thermo.cp_mole() + + property cp_mass: + """Specific heat capacity at constant pressure and composition [J/kg/K].""" + def __get__(self): + return self.thermo.cp_mass() + + property critical_temperature: + """Critical temperature [K].""" + def __get__(self): + return self.thermo.critTemperature() + + property critical_pressure: + """Critical pressure [Pa].""" + def __get__(self): + return self.thermo.critPressure() + + property critical_density: + """Critical density [kg/m³ or kmol/m³] depending on `basis`.""" + def __get__(self): + return self.thermo.critDensity() / self._mass_factor() + + property P_sat: + """Saturation pressure [Pa] at the current temperature.""" + def __get__(self): + return self.thermo.satPressure(self.T) + + property T_sat: + """Saturation temperature [K] at the current pressure.""" + def __get__(self): + return self.thermo.satTemperature(self.P) + + property auxiliary_data: + """ + Intermediate or model-specific parameters used by particular + derived classes. + """ + def __get__(self): + cdef CxxAnyMap values = self.thermo.getAuxiliaryData() + return anymap_to_py(values) + + ######## Methods to get/set the complete thermodynamic state ######## + + property state_size: + """ + Return size of vector defining internal state of the phase. + """ + def __get__(self): + return self.thermo.stateSize() + + property state: + """ + Get/Set the full thermodynamic state as a single array, arranged as + [temperature, density, mass fractions] for most phases. Useful mainly + in cases where it is desired to store many states in a multidimensional + array. + """ + def __get__(self): + cdef np.ndarray[np.double_t, ndim=1] state = np.empty(self.state_size) + self.thermo.saveState(span[double](&state[0], len(state))) + return state + + def __set__(self, state): + cdef np.ndarray[np.double_t, ndim=1] cstate = np.asarray(state) + self.thermo.restoreState(span[double](&cstate[0], len(cstate))) + + property TD: + """Get/Set temperature [K] and density [kg/m³ or kmol/m³].""" + def __get__(self): + return self.T, self.density + def __set__(self, values): + assert len(values) == 2, 'incorrect number of values' + T = values[0] if values[0] is not None else self.T + D = values[1] if values[1] is not None else self.density + self.thermo.setState_TD(T, D * self._mass_factor()) + + property TDX: + """ + Get/Set temperature [K], density [kg/m³ or kmol/m³], and mole fractions. + """ + def __get__(self): + return self.T, self.density, self.X + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + T = values[0] if values[0] is not None else self.T + D = values[1] if values[1] is not None else self.density + self.X = values[2] + self.thermo.setState_TD(T, D * self._mass_factor()) + + property TDY: + """ + Get/Set temperature [K] and density [kg/m³ or kmol/m³], and mass fractions. + """ + def __get__(self): + return self.T, self.density, self.Y + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + T = values[0] if values[0] is not None else self.T + D = values[1] if values[1] is not None else self.density + self.Y = values[2] + self.thermo.setState_TD(T, D * self._mass_factor()) + + property TP: + """Get/Set temperature [K] and pressure [Pa].""" + def __get__(self): + return self.T, self.P + def __set__(self, values): + assert len(values) == 2, 'incorrect number of values' + T = values[0] if values[0] is not None else self.T + P = values[1] if values[1] is not None else self.P + self.thermo.setState_TP(T, P) + + property TPX: + """Get/Set temperature [K], pressure [Pa], and mole fractions.""" + def __get__(self): + return self.T, self.P, self.X + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + T = values[0] if values[0] is not None else self.T + P = values[1] if values[1] is not None else self.P + self.X = values[2] + self.thermo.setState_TP(T, P) + + property TPY: + """Get/Set temperature [K], pressure [Pa], and mass fractions.""" + def __get__(self): + return self.T, self.P, self.Y + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + T = values[0] if values[0] is not None else self.T + P = values[1] if values[1] is not None else self.P + self.Y = values[2] + self.thermo.setState_TP(T, P) + + property UV: + """ + Get/Set internal energy [J/kg or J/kmol] and specific volume [m³/kg or m³/kmol]. + """ + def __get__(self): + return self.u, self.v + def __set__(self, values): + assert len(values) == 2, 'incorrect number of values' + U = values[0] if values[0] is not None else self.u + V = values[1] if values[1] is not None else self.v + self.thermo.setState_UV(U / self._mass_factor(), + V / self._mass_factor()) + + property UVX: + """ + Get/Set internal energy [J/kg or J/kmol], specific volume [m³/kg or m³/kmol], + and mole fractions. + """ + def __get__(self): + return self.u, self.v, self.X + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + U = values[0] if values[0] is not None else self.u + V = values[1] if values[1] is not None else self.v + self.X = values[2] + self.thermo.setState_UV(U / self._mass_factor(), + V / self._mass_factor()) + + property UVY: + """ + Get/Set internal energy [J/kg or J/kmol], specific volume [m³/kg or m³/kmol], + and mass fractions. + """ + def __get__(self): + return self.u, self.v, self.Y + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + U = values[0] if values[0] is not None else self.u + V = values[1] if values[1] is not None else self.v + self.Y = values[2] + self.thermo.setState_UV(U / self._mass_factor(), + V / self._mass_factor()) + + property DP: + """Get/Set density [kg/m³] and pressure [Pa].""" + def __get__(self): + return self.density, self.P + def __set__(self, values): + assert len(values) == 2, 'incorrect number of values' + D = values[0] if values[0] is not None else self.density + P = values[1] if values[1] is not None else self.P + self.thermo.setState_DP(D*self._mass_factor(), P) + + property DPX: + """Get/Set density [kg/m³], pressure [Pa], and mole fractions.""" + def __get__(self): + return self.density, self.P, self.X + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + D = values[0] if values[0] is not None else self.density + P = values[1] if values[1] is not None else self.P + self.X = values[2] + self.thermo.setState_DP(D*self._mass_factor(), P) + + property DPY: + """Get/Set density [kg/m³], pressure [Pa], and mass fractions.""" + def __get__(self): + return self.density, self.P, self.Y + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + D = values[0] if values[0] is not None else self.density + P = values[1] if values[1] is not None else self.P + self.Y = values[2] + self.thermo.setState_DP(D*self._mass_factor(), P) + + property HP: + """Get/Set enthalpy [J/kg or J/kmol] and pressure [Pa].""" + def __get__(self): + return self.h, self.P + def __set__(self, values): + assert len(values) == 2, 'incorrect number of values' + H = values[0] if values[0] is not None else self.h + P = values[1] if values[1] is not None else self.P + self.thermo.setState_HP(H / self._mass_factor(), P) + + property HPX: + """Get/Set enthalpy [J/kg or J/kmol], pressure [Pa] and mole fractions.""" + def __get__(self): + return self.h, self.P, self.X + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + H = values[0] if values[0] is not None else self.h + P = values[1] if values[1] is not None else self.P + self.X = values[2] + self.thermo.setState_HP(H / self._mass_factor(), P) + + property HPY: + """Get/Set enthalpy [J/kg or J/kmol], pressure [Pa] and mass fractions.""" + def __get__(self): + return self.h, self.P, self.Y + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + H = values[0] if values[0] is not None else self.h + P = values[1] if values[1] is not None else self.P + self.Y = values[2] + self.thermo.setState_HP(H / self._mass_factor(), P) + + property SP: + """Get/Set entropy [J/kg/K or J/kmol/K] and pressure [Pa].""" + def __get__(self): + return self.s, self.P + def __set__(self, values): + assert len(values) == 2, 'incorrect number of values' + S = values[0] if values[0] is not None else self.s + P = values[1] if values[1] is not None else self.P + self.thermo.setState_SP(S / self._mass_factor(), P) + + property SPX: + """Get/Set entropy [J/kg/K or J/kmol/K], pressure [Pa], and mole fractions.""" + def __get__(self): + return self.s, self.P, self.X + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + S = values[0] if values[0] is not None else self.s + P = values[1] if values[1] is not None else self.P + self.X = values[2] + self.thermo.setState_SP(S / self._mass_factor(), P) + + property SPY: + """Get/Set entropy [J/kg/K or J/kmol/K], pressure [Pa], and mass fractions.""" + def __get__(self): + return self.s, self.P, self.Y + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + S = values[0] if values[0] is not None else self.s + P = values[1] if values[1] is not None else self.P + self.Y = values[2] + self.thermo.setState_SP(S / self._mass_factor(), P) + + property SV: + """ + Get/Set entropy [J/kg/K or J/kmol/K] and specific volume [m³/kg or m³/kmol]. + """ + def __get__(self): + return self.s, self.v + def __set__(self, values): + assert len(values) == 2, 'incorrect number of values' + S = values[0] if values[0] is not None else self.s + V = values[1] if values[1] is not None else self.v + self.thermo.setState_SV(S / self._mass_factor(), + V / self._mass_factor()) + + property SVX: + """ + Get/Set entropy [J/kg/K or J/kmol/K], specific volume [m³/kg or m³/kmol], and + mole fractions. + """ + def __get__(self): + return self.s, self.v, self.X + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + S = values[0] if values[0] is not None else self.s + V = values[1] if values[1] is not None else self.v + self.X = values[2] + self.thermo.setState_SV(S / self._mass_factor(), + V / self._mass_factor()) + + property SVY: + """ + Get/Set entropy [J/kg/K or J/kmol/K], specific volume [m³/kg or m³/kmol], and + mass fractions. + """ + def __get__(self): + return self.s, self.v, self.Y + def __set__(self, values): + assert len(values) == 3, 'incorrect number of values' + S = values[0] if values[0] is not None else self.s + V = values[1] if values[1] is not None else self.v + self.Y = values[2] + self.thermo.setState_SV(S / self._mass_factor(), + V / self._mass_factor()) + + # partial molar / non-dimensional properties + property partial_molar_enthalpies: + """Array of species partial molar enthalpies [J/kmol].""" + def __get__(self): + return self._getArray1(thermo_getPartialMolarEnthalpies) + + property partial_molar_entropies: + """Array of species partial molar entropies [J/kmol/K].""" + def __get__(self): + return self._getArray1(thermo_getPartialMolarEntropies) + + property partial_molar_int_energies: + """Array of species partial molar internal energies [J/kmol].""" + def __get__(self): + return self._getArray1(thermo_getPartialMolarIntEnergies) + + property partial_molar_int_energies_TV: + r""" + Array of species `\tilde{u}_k = (\partial U/\partial n_k)_{T,V,n_{j\ne k}}` + values [J/kmol]. + + .. versionadded:: 4.0 + """ + def __get__(self): + return self._getArray1(thermo_getPartialMolarIntEnergies_TV) + + property chemical_potentials: + """Array of species chemical potentials [J/kmol].""" + def __get__(self): + return self._getArray1(thermo_getChemPotentials) + + property electrochemical_potentials: + """Array of species electrochemical potentials [J/kmol].""" + def __get__(self): + return self._getArray1(thermo_getElectrochemPotentials) + + property partial_molar_cp: + """ + Array of species partial molar specific heat capacities at constant + pressure [J/kmol/K]. + """ + def __get__(self): + return self._getArray1(thermo_getPartialMolarCp) + + property partial_molar_cv_TV: + r""" + Array of species `\tilde{c}_{v,k} = (\partial \tilde{u}_k/\partial T)_{V,n}` + values [J/kmol/K]. + + .. versionadded:: 4.0 + """ + def __get__(self): + return self._getArray1(thermo_getPartialMolarCv_TV) + + property partial_molar_volumes: + """Array of species partial molar volumes [m³/kmol].""" + def __get__(self): + return self._getArray1(thermo_getPartialMolarVolumes) + + property standard_enthalpies_RT: + """ + Array of nondimensional species standard-state enthalpies at the + current temperature and pressure. + """ + def __get__(self): + return self._getArray1(thermo_getEnthalpy_RT) + + property standard_entropies_R: + """ + Array of nondimensional species standard-state entropies at the + current temperature and pressure. + """ + def __get__(self): + return self._getArray1(thermo_getEntropy_R) + + property standard_int_energies_RT: + """ + Array of nondimensional species standard-state internal energies at the + current temperature and pressure. + """ + def __get__(self): + return self._getArray1(thermo_getIntEnergy_RT) + + property standard_gibbs_RT: + """ + Array of nondimensional species standard-state Gibbs free energies at + the current temperature and pressure. + """ + def __get__(self): + return self._getArray1(thermo_getGibbs_RT) + + property standard_cp_R: + """ + Array of nondimensional species standard-state specific heat capacities + at constant pressure at the current temperature and pressure. + """ + def __get__(self): + return self._getArray1(thermo_getCp_R) + + property activities: + """ + Array of nondimensional activities. Returns either molar or molal + activities depending on the convention of the thermodynamic model. + """ + def __get__(self): + return self._getArray1(thermo_getActivities) + + property activity_coefficients: + """ + Array of nondimensional, molar activity coefficients. + """ + def __get__(self): + return self._getArray1(thermo_getActivityCoefficients) + + ######## Miscellaneous properties ######## + property isothermal_compressibility: + """Isothermal compressibility [1/Pa].""" + def __get__(self): + return self.thermo.isothermalCompressibility() + + property thermal_expansion_coeff: + """Thermal expansion coefficient [1/K].""" + def __get__(self): + return self.thermo.thermalExpansionCoeff() + + property internal_pressure: + """ + Internal pressure [Pa]. + + .. versionadded:: 4.0 + """ + def __get__(self): + return self.thermo.internalPressure() + + property sound_speed: + """Speed of sound [m/s].""" + def __get__(self): + return self.thermo.soundSpeed() + + property min_temp: + """ + Minimum temperature for which the thermodynamic data for the phase are + valid. + """ + def __get__(self): + return self.thermo.minTemp() + + property max_temp: + """ + Maximum temperature for which the thermodynamic data for the phase are + valid. + """ + def __get__(self): + return self.thermo.maxTemp() + + property reference_pressure: + """Reference state pressure [Pa].""" + def __get__(self): + return self.thermo.refPressure() + + property electric_potential: + """Get/Set the electric potential [V] for this phase.""" + def __get__(self): + return self.thermo.electricPotential() + def __set__(self, double value): + self.thermo.setElectricPotential(value) + + property standard_concentration_units: + """Get standard concentration units for this phase.""" + def __get__(self): + cdef CxxUnits units = self.thermo.standardConcentrationUnits() + return Units.copy(units) + + # methods for plasma + property Te: + """Get/Set electron Temperature [K].""" + def __get__(self): + return self.thermo.electronTemperature() + + def __set__(self, value): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setElectronTemperature(value) + + property Pe: + """Get electron Pressure [Pa].""" + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.electronPressure() + + property reduced_electric_field: + """ + Get/Set the reduced electric field (E/N) [V·m²]. + + .. versionadded:: 3.2 + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.reducedElectricField() + + def __set__(self, value): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setReducedElectricField(value) + + property electric_field: + """ + Get/Set the electric field (E) [V/m]. + + This is the absolute electric field strength. It is related to the reduced + electric field (E/N) through the number density of neutrals. + + .. versionadded:: 3.2 + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.electricField() + + def __set__(self, value): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setElectricField(value) + + def set_discretized_electron_energy_distribution(self, levels, distribution): + """ + Set electron energy distribution. When this method is used, electron + temperature is calculated from the distribution. + + :param levels: + vector of electron energy levels [eV] + :param distribution: + vector of distribution + """ + if not self._enable_plasma: + raise TypeError('This method is invalid for ' + f'thermo model: {self.thermo_model}.') + # check length + if (len(levels) != len(distribution)): + raise ValueError('Length of levels and distribution are different') + + cdef np.ndarray[np.double_t, ndim=1] data_levels = \ + np.ascontiguousarray(levels, dtype=np.double) + cdef np.ndarray[np.double_t, ndim=1] data_dist = \ + np.ascontiguousarray(distribution, dtype=np.double) + self.plasma.setDiscretizedElectronEnergyDist( + span[double](&data_levels[0], data_levels.size), + span[double](&data_dist[0], data_dist.size)) + + property mean_temperature: + """Get the mean temperature of the plasma [K]. + This is a mole-weighted average of the electron and heavy species temperatures. + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.meanTemperature() + + def update_electron_energy_distribution(self): + """ + Update the electron energy distribution function to account for changes in + composition, temperature, pressure, or electric field strength. + + .. versionadded:: 3.2 + """ + if not self._enable_plasma: + raise TypeError('This method is invalid for ' + f'thermo model: {self.thermo_model}.') + self.plasma.updateElectronEnergyDistribution() + + property n_electron_energy_levels: + """ Number of electron energy levels """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.nElectronEnergyLevels() + + property electron_energy_levels: + """ Electron energy levels [eV]""" + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + cdef np.ndarray[np.double_t, ndim=1] data = np.empty( + self.n_electron_energy_levels) + self.plasma.getElectronEnergyLevels(span[double](&data[0], data.size)) + return data + def __set__(self, levels): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + cdef np.ndarray[np.double_t, ndim=1] data = \ + np.ascontiguousarray(levels, dtype=np.double) + self.plasma.setElectronEnergyLevels(span[double](&data[0], data.size)) + + property electron_energy_distribution: + """ Electron energy distribution """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + cdef np.ndarray[np.double_t, ndim=1] data = np.empty( + self.n_electron_energy_levels) + self.plasma.getElectronEnergyDistribution(span[double](&data[0], data.size)) + return data + + property isotropic_shape_factor: + """ Shape factor of isotropic-velocity distribution for electron energy """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.isotropicShapeFactor() + def __set__(self, x): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setIsotropicShapeFactor(x) + + property electron_energy_distribution_type: + """ Electron energy distribution type """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return pystr(self.plasma.electronEnergyDistributionType()) + def __set__(self, distribution_type): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setElectronEnergyDistributionType(stringify(distribution_type)) + + property mean_electron_energy: + """ Mean electron energy [eV] """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.meanElectronEnergy() + def __set__(self, double energy): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setMeanElectronEnergy(energy) + + property quadrature_method: + """ Quadrature method """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return pystr(self.plasma.quadratureMethod()) + def __set__(self, method): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.setQuadratureMethod(stringify(method)) + + property normalize_electron_energy_distribution_enabled: + """ Automatically normalize electron energy distribution """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.normalizeElectronEnergyDistEnabled() + def __set__(self, enable): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + self.plasma.enableNormalizeElectronEnergyDist(enable) + + property electron_species_name: + """ Electron species name """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return pystr(self.plasma.electronSpeciesName()) + + property elastic_power_loss: + """ + Elastic power loss (W/m³) + + .. versionadded:: 3.2 + """ + def __get__(self): + if not self._enable_plasma: + raise ThermoModelMethodError(self.thermo_model) + return self.plasma.elasticPowerLoss() + +cdef class InterfacePhase(ThermoPhase): + """ A class representing a surface, edge phase """ + def __cinit__(self, *args, **kwargs): + if not kwargs.get("init", True): + return + if not dynamic_cast[CxxSurfPhasePtr](self.thermo): + raise TypeError('Underlying ThermoPhase object is of the wrong type.') + self.surf = (self.thermo) + + property adjacent: + """ + A dictionary containing higher-dimensional phases adjacent to this interface, + for example bulk phases adjacent to a surface. + """ + def __get__(self): + return self._adjacent + + property site_density: + """ + Get/Set the site density. [kmol/m²] for surface phases; [kmol/m] for + edge phases. + """ + def __get__(self): + return self.surf.siteDensity() + def __set__(self, double value): + self.surf.setSiteDensity(value) + + property coverages: + """Get/Set the fraction of sites covered by each species.""" + def __get__(self): + cdef np.ndarray[np.double_t, ndim=1] data = np.empty(self.n_species) + self.surf.getCoverages(span[double](&data[0], data.size)) + if self._selected_species.size: + return data[self._selected_species] + else: + return data + + def __set__(self, theta): + if isinstance(theta, (dict, str, bytes)): + self.surf.setCoveragesByName(comp_map(theta)) + return + + if len(theta) != self.n_species: + raise ValueError("Array has incorrect length." + " Got {}, expected {}".format(len(theta), self.n_species)) + cdef np.ndarray[np.double_t, ndim=1] data = \ + np.ascontiguousarray(theta, dtype=np.double) + self.surf.setCoverages(span[double](&data[0], data.size)) + + def set_unnormalized_coverages(self, cov): + """ + Set the surface coverages without normalizing to force sum(cov) == 1.0. + Useful primarily when calculating derivatives with respect to cov[k] by + finite difference. + """ + cdef np.ndarray[np.double_t, ndim=1] data + if len(cov) == self.n_species: + data = np.ascontiguousarray(cov, dtype=np.double) + else: + raise ValueError("Array has incorrect length." + " Got {}, expected {}.".format(len(cov), self.n_species)) + self.surf.setCoveragesNoNorm(span[double](&data[0], data.size)) + + +cdef class PureFluid(ThermoPhase): + """ + A pure substance that can be a gas, a liquid, a mixed gas-liquid fluid, + or a fluid beyond its critical point. + """ + + property Q: + """ + Get/Set vapor fraction (quality). Can be set only when in the two-phase + region. + """ + def __get__(self): + return self.thermo.vaporFraction() + def __set__(self, Q): + if (self.P >= self.critical_pressure or + abs(self.P-self.P_sat)/self.P > 1e-4): + raise ValueError('Cannot set vapor quality outside the ' + 'two-phase region') + self.thermo.setState_Psat(self.P, Q) + + property TQ: + """Get/Set the temperature [K] and vapor fraction of a two-phase state.""" + def __get__(self): + return self.T, self.Q + def __set__(self, values): + T = values[0] if values[0] is not None else self.T + Q = values[1] if values[1] is not None else self.Q + self.thermo.setState_Tsat(T, Q) + + property PQ: + """Get/Set the pressure [Pa] and vapor fraction of a two-phase state.""" + def __get__(self): + return self.P, self.Q + def __set__(self, values): + P = values[0] if values[0] is not None else self.P + Q = values[1] if values[1] is not None else self.Q + self.thermo.setState_Psat(P, Q) + + property ST: + """Get/Set the entropy [J/kg/K] and temperature [K] of a PureFluid.""" + def __get__(self): + return self.s, self.T + def __set__(self, values): + S = values[0] if values[0] is not None else self.s + T = values[1] if values[1] is not None else self.T + self.thermo.setState_ST(S / self._mass_factor(), T) + + property TV: + """ + Get/Set the temperature [K] and specific volume [m³/kg] of a PureFluid. + """ + def __get__(self): + return self.T, self.v + def __set__(self, values): + T = values[0] if values[0] is not None else self.T + V = values[1] if values[1] is not None else self.v + self.thermo.setState_TV(T, V / self._mass_factor()) + + property PV: + """ + Get/Set the pressure [Pa] and specific volume [m³/kg] of a PureFluid. + """ + def __get__(self): + return self.P, self.v + def __set__(self, values): + P = values[0] if values[0] is not None else self.P + V = values[1] if values[1] is not None else self.v + self.thermo.setState_PV(P, V / self._mass_factor()) + + property UP: + """ + Get/Set the specific internal energy [J/kg] and the + pressure [Pa] of a PureFluid. + """ + def __get__(self): + return self.u, self.P + def __set__(self, values): + U = values[0] if values[0] is not None else self.u + P = values[1] if values[1] is not None else self.P + self.thermo.setState_UP(U / self._mass_factor(), P) + + property VH: + """ + Get/Set the specific volume [m³/kg] and the specific + enthalpy [J/kg] of a PureFluid. + """ + def __get__(self): + return self.v, self.h + def __set__(self, values): + V = values[0] if values[0] is not None else self.v + H = values[1] if values[1] is not None else self.h + self.thermo.setState_VH(V/self._mass_factor(), H/self._mass_factor()) + + property TH: + """ + Get/Set the temperature [K] and the specific enthalpy [J/kg] + of a PureFluid. + """ + def __get__(self): + return self.T, self.h + def __set__(self, values): + T = values[0] if values[0] is not None else self.T + H = values[1] if values[1] is not None else self.h + self.thermo.setState_TH(T, H / self._mass_factor()) + + property SH: + """ + Get/Set the specific entropy [J/kg/K] and the specific + enthalpy [J/kg] of a PureFluid. + """ + def __get__(self): + return self.s, self.h + def __set__(self, values): + S = values[0] if values[0] is not None else self.s + H = values[1] if values[1] is not None else self.h + self.thermo.setState_SH(S/self._mass_factor(), H/self._mass_factor()) + + property TDQ: + """ + Get the temperature [K], density [kg/m³ or kmol/m³], and vapor fraction. + """ + def __get__(self): + return self.T, self.density, self.Q + + property TPQ: + """ + Get/Set the temperature [K], pressure [Pa], and vapor fraction of a + PureFluid. + + An Exception is raised if the thermodynamic state is not consistent. + """ + def __get__(self): + return self.T, self.P, self.Q + def __set__(self, values): + T = values[0] if values[0] is not None else self.T + P = values[1] if values[1] is not None else self.P + Q = values[2] if values[2] is not None else self.Q + self.thermo.setState_TPQ(T, P, Q) + + property UVQ: + """ + Get the internal energy [J/kg or J/kmol], specific volume + [m³/kg or m³/kmol], and vapor fraction. + """ + def __get__(self): + return self.u, self.v, self.Q + + property DPQ: + """Get the density [kg/m³], pressure [Pa], and vapor fraction.""" + def __get__(self): + return self.density, self.P, self.Q + + property HPQ: + """ + Get the enthalpy [J/kg or J/kmol], pressure [Pa] and vapor fraction. + """ + def __get__(self): + return self.h, self.P, self.Q + + property SPQ: + """ + Get the entropy [J/kg/K or J/kmol/K], pressure [Pa], and vapor fraction. + """ + def __get__(self): + return self.s, self.P, self.Q + + property SVQ: + """ + Get the entropy [J/kg/K or J/kmol/K], specific volume [m³/kg or + m³/kmol], and vapor fraction. + """ + def __get__(self): + return self.s, self.v, self.Q + + +class Element: + """ + An element or a named isotope defined in Cantera. + + Class `Element` gets data for the elements and isotopes defined in + :ct:`Elements.cpp`. This class can be used in two ways. The + first way is to get information about all of the elements stored in + Cantera. The three attributes `num_elements_defined`, + `element_symbols`, and `element_names` can be accessed by:: + + >>> ct.Element.num_elements_defined + >>> ct.Element.element_symbols + >>> ct.Element.element_names + + Otherwise, if the class `Element` is called with an argument, it + stores the data about that particular element. For example:: + + >>> ar_sym = ct.Element('Ar') + >>> ar_name = ct.Element('argon') + >>> ar_num = ct.Element(18) + + would all create instances with the information for argon. The + available argument options to create an instance of the `Element` + class with the element information are the `name`, `symbol`, and + `atomic_number`. Once an instance of the class is made, the `name`, + `atomic_number`, `symbol`, and atomic `weight` can be accessed as + attributes of the instance of the `Element` class. + + >>> ar_sym.name + 'argon' + >>> ar_sym.weight + 39.948 + >>> ar_sym.atomic_number + 18 + >>> ar_sym.symbol + 'Ar' + + The elements available are listed below, in the `element_symbols` + and `element_names` attribute documentation. + """ + + #: The number of named elements (not isotopes) defined in Cantera + num_elements_defined = numElementsDefined() + + #: A list of the symbols of all the elements (not isotopes) defined in Cantera + element_symbols = tuple(pystr(s) for s in elementSymbols()) + + #: A list of the names of all the elements (not isotopes) defined in Cantera + element_names = tuple(pystr(n) for n in elementNames()) + + def __init__(self, arg): + if isinstance(arg, (str, bytes)): + try: + # Assume the argument is the element symbol and try to get the name + self._name = pystr(getElementName(stringify(arg))) + except CanteraError: + # If getting the name failed, the argument must be the name + self._symbol = pystr(getElementSymbol(stringify(arg))) + self._name = arg.lower() + else: + self._symbol = arg + + self._atomic_number = getAtomicNumber(stringify(arg)) + self._weight = getElementWeight(stringify(arg)) + elif isinstance(arg, int): + self._atomic_number = arg + self._name = pystr(getElementName(arg)) + self._symbol = pystr(getElementSymbol(arg)) + self._weight = getElementWeight(arg) + else: + raise TypeError('The input argument to Element must be a string ' + 'or an integer') + + @property + def name(self): + """The name of the element or isotope.""" + return self._name + + @property + def atomic_number(self): + """The atomic number of the element or isotope.""" + return self._atomic_number + + @property + def symbol(self): + """The symbol of the element or isotope.""" + return self._symbol + + @property + def weight(self): + """The atomic weight of the element or isotope.""" + return self._weight diff --git a/src/thermo/PlasmaPhase.cpp b/src/thermo/PlasmaPhase.cpp index 354448dbfa4..5902a42f349 100644 --- a/src/thermo/PlasmaPhase.cpp +++ b/src/thermo/PlasmaPhase.cpp @@ -25,13 +25,13 @@ PlasmaPhase::PlasmaPhase(const string& inputFile, const string& id_) { initThermoFile(inputFile, id_); - // initial electron temperature + // Initial electron temperature. m_electronTemp = temperature(); - // Initialize the Boltzmann Solver + // Initialize the Boltzmann Solver. m_eedfSolver = make_unique(this); - // Set Energy Grid (Hardcoded Defaults for Now) + // Set Energy Grid (Hardcoded Defaults for Now). double kTe_max = 60; size_t nGridCells = 301; m_nPoints = nGridCells + 1; @@ -47,11 +47,172 @@ PlasmaPhase::~PlasmaPhase() soln->kinetics()->removeReactionAddedCallback(this); } for (size_t k = 0; k < nCollisions(); k++) { - // remove callback + // Remove callback. m_collisions[k]->removeSetRateCallback(this); } } +void PlasmaPhase::initThermo() +{ + IdealGasPhase::initThermo(); + + // Check if there is an electron species in the phase. + if (m_electronSpeciesIndex == npos) { + throw CanteraError("PlasmaPhase::initThermo", + "No electron species found."); + } +} + +void PlasmaPhase::updateThermo() const +{ + // Update the heavy species thermodynamic properties + // before updating the electron species properties. + IdealGasPhase::updateThermo(); + static const int cacheId = m_cache.getId(); + CachedScalar cached = m_cache.getScalar(cacheId); + double tempNow = temperature(); + double electronTempNow = electronTemperature(); + size_t k = m_electronSpeciesIndex; + // If the electron temperature has changed since the last time these + // properties were computed, recompute them. + if (cached.state1 != tempNow || cached.state2 != electronTempNow) { + // Evaluate the electron species thermodynamic properties + // at the electron temperature. + m_spthermo.update_single(k, electronTemperature(), + m_cp0_R[k], m_h0_RT[k], m_s0_R[k]); + cached.state1 = tempNow; + cached.state2 = electronTempNow; + + // Update the electron Gibbs functions, with the electron temperature. + m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; + } +} + +// ================================================================= // +// Overridden from IdealGasPhase or ThermoPhase // +// ================================================================= // + +bool PlasmaPhase::addSpecies(shared_ptr spec) +{ + bool added = IdealGasPhase::addSpecies(spec); + size_t k = m_kk - 1; + + if ((spec->name == "e" || spec->name == "Electron") || + (spec->composition.find("E") != spec->composition.end() && + spec->composition.size() == 1 && + spec->composition["E"] == 1)) { + if (m_electronSpeciesIndex == npos) { + m_electronSpeciesIndex = k; + } else { + throw CanteraError("PlasmaPhase::addSpecies", + "Cannot add species, {}. " + "Only one electron species is allowed.", spec->name); + } + } + return added; +} + +void PlasmaPhase::setSolution(std::weak_ptr soln) { + ThermoPhase::setSolution(soln); + // Register callback function to be executed + // when the thermo or kinetics object changed. + if (shared_ptr soln = m_soln.lock()) { + soln->registerChangedCallback(this, [&]() { + setCollisions(); + }); + } +} + +void PlasmaPhase::getParameters(AnyMap& phaseNode) const +{ + IdealGasPhase::getParameters(phaseNode); + AnyMap eedf; + eedf["type"] = m_distributionType; + vector levels(m_nPoints); + Eigen::Map(levels.data(), m_nPoints) = m_electronEnergyLevels; + eedf["energy-levels"] = levels; + if (m_distributionType == "isotropic") { + eedf["shape-factor"] = m_isotropicShapeFactor; + eedf["mean-electron-energy"].setQuantity(meanElectronEnergy(), "eV"); + } else if (m_distributionType == "discretized") { + vector dist(m_nPoints); + Eigen::Map(dist.data(), m_nPoints) = m_electronEnergyDist; + eedf["distribution"] = dist; + eedf["normalize"] = m_do_normalizeElectronEnergyDist; + } + phaseNode["electron-energy-distribution"] = std::move(eedf); +} + +void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) +{ + IdealGasPhase::setParameters(phaseNode, rootNode); + if (phaseNode.hasKey("electron-energy-distribution")) { + const AnyMap eedf = phaseNode["electron-energy-distribution"].as(); + m_distributionType = eedf["type"].asString(); + if (m_distributionType == "isotropic") { + if (eedf.hasKey("shape-factor")) { + setIsotropicShapeFactor(eedf["shape-factor"].asDouble()); + } else { + throw CanteraError("PlasmaPhase::setParameters", + "isotropic type requires shape-factor key."); + } + if (eedf.hasKey("mean-electron-energy")) { + double energy = eedf.convert("mean-electron-energy", "eV"); + setMeanElectronEnergy(energy); + } else { + throw CanteraError("PlasmaPhase::setParameters", + "isotropic type requires electron-temperature key."); + } + if (eedf.hasKey("energy-levels")) { + auto levels = eedf["energy-levels"].asVector(); + setElectronEnergyLevels(levels); + } + setIsotropicElectronEnergyDistribution(); + } else if (m_distributionType == "discretized") { + if (!eedf.hasKey("energy-levels")) { + throw CanteraError("PlasmaPhase::setParameters", + "Cannot find key energy-levels."); + } + if (!eedf.hasKey("distribution")) { + throw CanteraError("PlasmaPhase::setParameters", + "Cannot find key distribution."); + } + if (eedf.hasKey("normalize")) { + enableNormalizeElectronEnergyDist(eedf["normalize"].asBool()); + } + auto levels = eedf["energy-levels"].asVector(); + auto distribution = eedf["distribution"].asVector(levels.size()); + setDiscretizedElectronEnergyDist(levels, distribution); + } + } + + if (rootNode.hasKey("electron-collisions")) { + for (const auto& item : rootNode["electron-collisions"].asVector()) { + auto rate = make_shared(item); + Composition reactants, products; + reactants[item["target"].asString()] = 1; + reactants[electronSpeciesName()] = 1; + if (item.hasKey("product")) { + products[item["product"].asString()] = 1; + } else { + products[item["target"].asString()] = 1; + } + products[electronSpeciesName()] = 1; + if (rate->kind() == "ionization") { + products[electronSpeciesName()] += 1; + } else if (rate->kind() == "attachment") { + products[electronSpeciesName()] -= 1; + } + auto R = make_shared(reactants, products, rate); + addCollision(R); + } + } +} + +// ================================================================= // +// Electron Energy Distribution Functions // +// ================================================================= // + void PlasmaPhase::updateElectronEnergyDistribution() { if (m_distributionType == "discretized") { @@ -128,6 +289,11 @@ void PlasmaPhase::setIsotropicElectronEnergyDistribution() } void PlasmaPhase::setElectronTemperature(const double Te) { + // If the electron temperature is negative, throw an error. + if (Te < 0.0) { + throw CanteraError("PlasmaPhase::setElectronTemperature", + "Electron temperature cannot be negative."); + } m_electronTemp = Te; updateElectronEnergyDistribution(); } @@ -177,7 +343,7 @@ void PlasmaPhase::electronEnergyDistributionChanged() void PlasmaPhase::electronEnergyLevelChanged() { m_levelNum++; - // Cross sections are interpolated on the energy levels + // Cross sections are interpolated on the energy levels. if (m_collisions.size() > 0) { vector energyLevels(m_nPoints); MappedVector(energyLevels.data(), m_nPoints) = m_electronEnergyLevels; @@ -239,12 +405,12 @@ void PlasmaPhase::setDiscretizedElectronEnergyDist(span levels, void PlasmaPhase::updateElectronTemperatureFromEnergyDist() { - // calculate mean electron energy and electron temperature + // Calculate mean electron energy and electron temperature. Eigen::ArrayXd eps52 = m_electronEnergyLevels.pow(5./2.); double epsilon_m = 2.0 / 5.0 * numericalQuadrature(m_quadratureMethod, m_electronEnergyDist, eps52); if (epsilon_m < 0.0 && m_quadratureMethod == "simpson") { - // try trapezoidal method + // Try trapezoidal method. epsilon_m = 2.0 / 5.0 * numericalQuadrature( "trapezoidal", m_electronEnergyDist, eps52); } @@ -262,133 +428,6 @@ void PlasmaPhase::setIsotropicShapeFactor(double x) { updateElectronEnergyDistribution(); } -void PlasmaPhase::getParameters(AnyMap& phaseNode) const -{ - IdealGasPhase::getParameters(phaseNode); - AnyMap eedf; - eedf["type"] = m_distributionType; - vector levels(m_nPoints); - Eigen::Map(levels.data(), m_nPoints) = m_electronEnergyLevels; - eedf["energy-levels"] = levels; - if (m_distributionType == "isotropic") { - eedf["shape-factor"] = m_isotropicShapeFactor; - eedf["mean-electron-energy"].setQuantity(meanElectronEnergy(), "eV"); - } else if (m_distributionType == "discretized") { - vector dist(m_nPoints); - Eigen::Map(dist.data(), m_nPoints) = m_electronEnergyDist; - eedf["distribution"] = dist; - eedf["normalize"] = m_do_normalizeElectronEnergyDist; - } - phaseNode["electron-energy-distribution"] = std::move(eedf); -} - -void PlasmaPhase::setParameters(const AnyMap& phaseNode, const AnyMap& rootNode) -{ - IdealGasPhase::setParameters(phaseNode, rootNode); - if (phaseNode.hasKey("electron-energy-distribution")) { - const AnyMap eedf = phaseNode["electron-energy-distribution"].as(); - m_distributionType = eedf["type"].asString(); - if (m_distributionType == "isotropic") { - if (eedf.hasKey("shape-factor")) { - setIsotropicShapeFactor(eedf["shape-factor"].asDouble()); - } else { - throw CanteraError("PlasmaPhase::setParameters", - "isotropic type requires shape-factor key."); - } - if (eedf.hasKey("mean-electron-energy")) { - double energy = eedf.convert("mean-electron-energy", "eV"); - setMeanElectronEnergy(energy); - } else { - throw CanteraError("PlasmaPhase::setParameters", - "isotropic type requires electron-temperature key."); - } - if (eedf.hasKey("energy-levels")) { - auto levels = eedf["energy-levels"].asVector(); - setElectronEnergyLevels(levels); - } - setIsotropicElectronEnergyDistribution(); - } else if (m_distributionType == "discretized") { - if (!eedf.hasKey("energy-levels")) { - throw CanteraError("PlasmaPhase::setParameters", - "Cannot find key energy-levels."); - } - if (!eedf.hasKey("distribution")) { - throw CanteraError("PlasmaPhase::setParameters", - "Cannot find key distribution."); - } - if (eedf.hasKey("normalize")) { - enableNormalizeElectronEnergyDist(eedf["normalize"].asBool()); - } - auto levels = eedf["energy-levels"].asVector(); - auto distribution = eedf["distribution"].asVector(levels.size()); - setDiscretizedElectronEnergyDist(levels, distribution); - } - } - - if (rootNode.hasKey("electron-collisions")) { - for (const auto& item : rootNode["electron-collisions"].asVector()) { - auto rate = make_shared(item); - Composition reactants, products; - reactants[item["target"].asString()] = 1; - reactants[electronSpeciesName()] = 1; - if (item.hasKey("product")) { - products[item["product"].asString()] = 1; - } else { - products[item["target"].asString()] = 1; - } - products[electronSpeciesName()] = 1; - if (rate->kind() == "ionization") { - products[electronSpeciesName()] += 1; - } else if (rate->kind() == "attachment") { - products[electronSpeciesName()] -= 1; - } - auto R = make_shared(reactants, products, rate); - addCollision(R); - } - } -} - -bool PlasmaPhase::addSpecies(shared_ptr spec) -{ - bool added = IdealGasPhase::addSpecies(spec); - size_t k = m_kk - 1; - - if ((spec->name == "e" || spec->name == "Electron") || - (spec->composition.find("E") != spec->composition.end() && - spec->composition.size() == 1 && - spec->composition["E"] == 1)) { - if (m_electronSpeciesIndex == npos) { - m_electronSpeciesIndex = k; - } else { - throw CanteraError("PlasmaPhase::addSpecies", - "Cannot add species, {}. " - "Only one electron species is allowed.", spec->name); - } - } - return added; -} - -void PlasmaPhase::initThermo() -{ - IdealGasPhase::initThermo(); - - // Check electron species - if (m_electronSpeciesIndex == npos) { - throw CanteraError("PlasmaPhase::initThermo", - "No electron species found."); - } -} - -void PlasmaPhase::setSolution(std::weak_ptr soln) { - ThermoPhase::setSolution(soln); - // register callback function to be executed - // when the thermo or kinetics object changed - if (shared_ptr soln = m_soln.lock()) { - soln->registerChangedCallback(this, [&]() { - setCollisions(); - }); - } -} void PlasmaPhase::setCollisions() { @@ -398,8 +437,8 @@ void PlasmaPhase::setCollisions() return; } - // add collision from the initial list of reactions. Only add reactions we - // haven't seen before + // Add collision from the initial list of reactions. + // Only add reactions we haven't seen before. set existing; for (auto& R : m_collisions) { existing.insert(R.get()); @@ -413,8 +452,8 @@ void PlasmaPhase::setCollisions() addCollision(R); } - // register callback when reaction is added later - // Modifying collision reactions is not supported + // Register callback when reaction is added later. + // Modifying collision reactions is not supported. kin->registerReactionAddedCallback(this, [this, kin]() { size_t i = kin->nReactions() - 1; if (kin->reaction(i)->type() == "electron-collision-plasma") { @@ -428,18 +467,18 @@ void PlasmaPhase::addCollision(shared_ptr collision) { size_t i = nCollisions(); - // setup callback to signal updating the cross-section-related - // parameters + // Setup callback to signal updating the cross-section-related + // parameters. collision->registerSetRateCallback(this, [this, i, collision]() { m_interp_cs_ready[i] = false; m_collisionRates[i] = std::dynamic_pointer_cast(collision->rate()); }); - // Identify target species for electron-collision reactions + // Identify target species for electron-collision reactions. string target; for (const auto& [name, _] : collision->reactants) { - // Reactants are expected to be electrons and the target species + // Reactants are expected to be electrons and the target species. if (name != electronSpeciesName()) { m_targetSpeciesIndices.emplace_back(speciesIndex(name, true)); target = name; @@ -456,11 +495,11 @@ void PlasmaPhase::addCollision(shared_ptr collision) std::dynamic_pointer_cast(collision->rate())); m_interp_cs_ready.emplace_back(false); - // resize parameters + // Resize parameters. m_elasticElectronEnergyLossCoefficients.resize(nCollisions()); updateInterpolatedCrossSection(i); - // Set up data used by Boltzmann solver + // Set up data used by Boltzmann solver. auto& rate = *m_collisionRates.back(); string kind = m_collisionRates.back()->kind(); @@ -501,12 +540,12 @@ bool PlasmaPhase::updateInterpolatedCrossSection(size_t i) void PlasmaPhase::updateElectronEnergyDistDifference() { m_electronEnergyDistDiff.resize(nElectronEnergyLevels()); - // Forward difference for the first point + // Forward difference for the first point. m_electronEnergyDistDiff[0] = (m_electronEnergyDist[1] - m_electronEnergyDist[0]) / (m_electronEnergyLevels[1] - m_electronEnergyLevels[0]); - // Central difference for the middle points + // Central difference for the middle points. for (size_t i = 1; i < m_nPoints - 1; i++) { double h1 = m_electronEnergyLevels[i+1] - m_electronEnergyLevels[i]; double h0 = m_electronEnergyLevels[i] - m_electronEnergyLevels[i-1]; @@ -516,7 +555,7 @@ void PlasmaPhase::updateElectronEnergyDistDifference() (h1 * h0) / (h1 + h0); } - // Backward difference for the last point + // Backward difference for the last point. m_electronEnergyDistDiff[m_nPoints-1] = (m_electronEnergyDist[m_nPoints-1] - m_electronEnergyDist[m_nPoints-2]) / @@ -526,11 +565,11 @@ void PlasmaPhase::updateElectronEnergyDistDifference() void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() { - // cache of cross section plus distribution plus energy-level number + // Cache of cross section plus distribution plus energy-level number. static const int cacheId = m_cache.getId(); CachedScalar last_stateNum = m_cache.getScalar(cacheId); - // combine the distribution and energy level number + // Combine the distribution and energy level number. int stateNum = m_distNum + m_levelNum; vector interpChanged(m_collisions.size()); @@ -539,15 +578,15 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficients() } if (last_stateNum.validate(temperature(), stateNum)) { - // check each cross section, and only update coefficients that - // the interpolated cross sections change + // Check each cross section, and only update coefficients that + // the interpolated cross sections change. for (size_t i = 0; i < m_collisions.size(); i++) { if (interpChanged[i]) { updateElasticElectronEnergyLossCoefficient(i); } } } else { - // update every coefficient if distribution, temperature, + // Update every coefficient if distribution, temperature, // or energy levels change. for (size_t i = 0; i < m_collisions.size(); i++) { updateElasticElectronEnergyLossCoefficient(i); @@ -566,10 +605,10 @@ void PlasmaPhase::updateElasticElectronEnergyLossCoefficient(size_t i) m_collisionRates[i]->crossSectionInterpolated().size() ); - // Mass ratio calculation + // Mass ratio calculation. double mass_ratio = ElectronMass / molecularWeight(k) * Avogadro; - // Calculate the rate using Simpson's rule or trapezoidal rule + // Calculate the rate using Simpson's rule or trapezoidal rule. Eigen::ArrayXd f0_plus = m_electronEnergyDist + Boltzmann * temperature() / ElectronCharge * m_electronEnergyDistDiff; m_elasticElectronEnergyLossCoefficients[i] = 2.0 * mass_ratio * gamma * @@ -605,26 +644,57 @@ double PlasmaPhase::elasticPowerLoss() return q_elastic; } -void PlasmaPhase::updateThermo() const + +double PlasmaPhase::electronMobility() const { - IdealGasPhase::updateThermo(); - static const int cacheId = m_cache.getId(); - CachedScalar cached = m_cache.getScalar(cacheId); - double tempNow = temperature(); - double electronTempNow = electronTemperature(); - size_t k = m_electronSpeciesIndex; - // If the electron temperature has changed since the last time these - // properties were computed, recompute them. - if (cached.state1 != tempNow || cached.state2 != electronTempNow) { - m_spthermo.update_single(k, electronTemperature(), - m_cp0_R[k], m_h0_RT[k], m_s0_R[k]); - cached.state1 = tempNow; - cached.state2 = electronTempNow; + // Only implemented when using the Boltzmann two-term EEDF + if (m_distributionType == "Boltzmann-two-term") { + return m_eedfSolver->getElectronMobility(); + } else { + throw NotImplementedError("PlasmaPhase::electronMobility", + "Electron mobility is only available for 'Boltzmann-two-term' " + "electron energy distributions."); } - // update the species Gibbs functions - m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k]; } + +double PlasmaPhase::jouleHeatingPower() const +{ + // sigma = e * n_e * mu_e [S/m]; q_J = sigma * E^2 [W/m^3] + const double mu_e = electronMobility(); // m^2 / (V·s) + if (mu_e <= 0.0) { + return 0.0; + } + const double ne = concentration(m_electronSpeciesIndex) * Avogadro; // m^-3 + if (ne <= 0.0) { + return 0.0; + } + const double E = electricField(); // V/m + if (E <= 0.0) { + return 0.0; + } + const double sigma = ElectronCharge * ne * mu_e; // S/m + return sigma * E * E; // W/m^3 +} + +double PlasmaPhase::intrinsicHeating() +{ + // Joule heating: sigma * E^2 [W/m^3] + const double qJ = jouleHeatingPower(); + checkFinite(qJ); + + // Elastic + inelastic recoil power loss [W/m^3] + double qElastic = elasticPowerLoss(); + checkFinite(qElastic); + + return qJ + qElastic; +} + + +// ================================================================= // +// Molar Thermodynamic Properties of the Solution // +// ================================================================= // + double PlasmaPhase::enthalpy_mole() const { double value = IdealGasPhase::enthalpy_mole(); value += GasConstant * (electronTemperature() - temperature()) * @@ -633,17 +703,6 @@ double PlasmaPhase::enthalpy_mole() const { return value; } -double PlasmaPhase::intEnergy_mole() const -{ - m_work.resize(m_kk); - getPartialMolarIntEnergies(m_work); - double u = 0.0; - for (size_t k = 0; k < m_kk; ++k) { - u += moleFraction(k) * m_work[k]; - } - return u; -} - double PlasmaPhase::entropy_mole() const { m_work.resize(m_kk); @@ -666,16 +725,55 @@ double PlasmaPhase::gibbs_mole() const return g; } -void PlasmaPhase::getGibbs_ref(span g) const +double PlasmaPhase::intEnergy_mole() const { - IdealGasPhase::getGibbs_ref(g); - g[m_electronSpeciesIndex] *= electronTemperature() / temperature(); + m_work.resize(m_kk); + getPartialMolarIntEnergies(m_work); + double u = 0.0; + for (size_t k = 0; k < m_kk; ++k) { + u += moleFraction(k) * m_work[k]; + } + return u; } -void PlasmaPhase::getStandardVolumes_ref(span vol) const +// ================================================================= // +// Mechanical Equation of State // +// ================================================================= // + +double PlasmaPhase::meanTemperature() const { - IdealGasPhase::getStandardVolumes_ref(vol); - vol[m_electronSpeciesIndex] *= electronTemperature() / temperature(); + double T_g = temperature(); + double T_e = electronTemperature(); + double X_e = moleFraction(m_electronSpeciesIndex); + return T_g + X_e * (T_e - T_g); +} + +double PlasmaPhase::pressure() const { + return GasConstant * meanTemperature() * density() / meanMolecularWeight(); +} + +// ================================================================= // +// Chemical Potentials and Activities // +// ================================================================= // + + +double PlasmaPhase::standardConcentration(size_t k) const +{ + return pressure() / (GasConstant * meanTemperature()); +} + + + +// ================================================================= // +// Partial Molar Properties of the Solution // +// ================================================================= // + +void PlasmaPhase::getChemPotentials(span mu) const +{ + IdealGasPhase::getChemPotentials(mu); + size_t k = m_electronSpeciesIndex; + double xx = std::max(SmallNumber, moleFraction(k)); + mu[k] += (RTe() - RT()) * log(xx); } void PlasmaPhase::getPartialMolarEnthalpies(span hbar) const @@ -684,13 +782,13 @@ void PlasmaPhase::getPartialMolarEnthalpies(span hbar) const hbar[m_electronSpeciesIndex] *= electronTemperature() / temperature(); } -void PlasmaPhase::getPartialMolarEntropies(span sbar) const -{ - IdealGasPhase::getPartialMolarEntropies(sbar); - double logp = log(pressure()); - double logpe = log(electronPressure()); - sbar[m_electronSpeciesIndex] += GasConstant * (logp - logpe); -} +// void PlasmaPhase::getPartialMolarEntropies(span sbar) const +// { +// IdealGasPhase::getPartialMolarEntropies(sbar); +// double logp = log(pressure()); +// double logpe = log(electronPressure()); +// sbar[m_electronSpeciesIndex] += GasConstant * (logp - logpe); +// } void PlasmaPhase::getPartialMolarIntEnergies(span ubar) const { @@ -699,100 +797,238 @@ void PlasmaPhase::getPartialMolarIntEnergies(span ubar) const for (size_t k = 0; k < m_kk; k++) { ubar[k] = RT() * (_h[k] - 1.0); } + // Redefine it for the electron species. size_t k = m_electronSpeciesIndex; ubar[k] = RTe() * (_h[k] - 1.0); } -void PlasmaPhase::getChemPotentials(span mu) const +void PlasmaPhase::getPartialMolarVolumes(span vbar) const { - IdealGasPhase::getChemPotentials(mu); - size_t k = m_electronSpeciesIndex; - double xx = std::max(SmallNumber, moleFraction(k)); - mu[k] += (RTe() - RT()) * log(xx); + double vol = RT() / pressure(); + for (size_t k = 0; k < m_kk; k++) { + vbar[k] = vol; + } + vbar[m_electronSpeciesIndex] = RTe() / pressure(); } +// ================================================================= // +// Properties of the Standard State of the Species in the Solution // +// ================================================================= // + void PlasmaPhase::getStandardChemPotentials(span muStar) const { - IdealGasPhase::getStandardChemPotentials(muStar); + // After calling PlasmaPhase::getGibbs_ref, muStar = g_k(T_k). + // g_k is evaluated at T for heavy species and at Te for electrons. + getGibbs_ref(muStar); + + // Then, we need to add R*T_k*ln(P/Pref) to g_k. + // .. For heavy species, mu_star = g_k(T) + R*T*ln(P/Pref) + double tmp = log(pressure() / refPressure()) * RT(); + for (size_t k = 0; k < m_kk; k++) { + muStar[k] += tmp; + } + // .. For electrons, mu_star = g_k(Te) + R*T_e*ln(P/Pref) size_t k = m_electronSpeciesIndex; muStar[k] -= log(pressure() / refPressure()) * RT(); - muStar[k] += log(electronPressure() / refPressure()) * RTe(); + muStar[k] += log(pressure() / refPressure()) * RTe(); } -void PlasmaPhase::getEntropy_R(span sr) const +void PlasmaPhase::getStandardVolumes(span vol) const { - checkArraySize("PlasmaPhase::getEntropy_R", sr.size(), m_kk); - auto _s = entropy_R_ref(); - copy(_s.begin(), _s.end(), sr.begin()); - double tmp = log(pressure() / refPressure()); + double tmp = RT() / pressure(); for (size_t k = 0; k < m_kk; k++) { - if (k != m_electronSpeciesIndex) { - sr[k] -= tmp; + vol[k] = tmp; + } + vol[m_electronSpeciesIndex] = RTe() / pressure(); +} + + +// ================================================================= // +// Thermodynamic Values for the Species Reference States // +// ================================================================= // + +void PlasmaPhase::getGibbs_ref(span g) const +{ + IdealGasPhase::getGibbs_ref(g); + g[m_electronSpeciesIndex] *= electronTemperature() / temperature(); +} + +void PlasmaPhase::getStandardVolumes_ref(span vol) const +{ + IdealGasPhase::getStandardVolumes_ref(vol); + vol[m_electronSpeciesIndex] *= electronTemperature() / temperature(); +} + + +// ================================================================= // +// Setting the State // +// ================================================================= // + + +void PlasmaPhase::setState(const AnyMap& input_state) +{ + AnyMap state = input_state; + + // Remap allowable synonyms. + if (state.hasKey("mass-fractions")) { + state["Y"] = state["mass-fractions"]; + state.erase("mass-fractions"); + } + if (state.hasKey("mole-fractions")) { + state["X"] = state["mole-fractions"]; + state.erase("mole-fractions"); + } + if (state.hasKey("gas-temperature")) { + state["T"] = state["gas-temperature"]; + } + if (state.hasKey("Tg")) { + state["T"] = state["Tg"]; + } + if (state.hasKey("electron-temperature")) { + state["Te"] = state["electron-temperature"]; + } + if (state.hasKey("pressure")) { + state["P"] = state["pressure"]; + } + if (state.hasKey("density")) { + state["D"] = state["density"]; + } + // Set composition. + if (state.hasKey("X")) { + if (state["X"].is()) { + setMoleFractionsByName(state["X"].asString()); + } else { + setMoleFractionsByName(state["X"].asMap()); + } + state.erase("X"); + } else if (state.hasKey("Y")) { + if (state["Y"].is()) { + setMassFractionsByName(state["Y"].asString()); } else { - sr[k] -= log(electronPressure() / refPressure()); + setMassFractionsByName(state["Y"].asMap()); } + state.erase("Y"); + } + + // Set thermodynamic state using whichever property set is found + if (state.hasKey("T") && state.hasKey("Te") && state.hasKey("P")) { + setState_TgTeP( + state.convert("T", "K"), + state.convert("Te", "K"), + state.convert("P", "Pa") + ); + } else if (state.hasKey("T") && state.hasKey("Te") && state.hasKey("D")) { + setState_TgTeD( + state.convert("T", "K"), + state.convert("Te", "K"), + state.convert("D", "kg/m^3") + ); + } else if (state.hasKey("T") && state.hasKey("P")) { + setState_TP( + state.convert("T", "K"), + state.convert("P", "Pa") + ); + } else if (state.hasKey("T") && state.hasKey("D")) { + setState_TD( + state.convert("T", "K"), + state.convert("D", "kg/m^3") + ); + } else { + throw CanteraError("PlasmaPhase::setState", + "'state' did not specify a recognized set of properties.\n" + "Keys provided were: {}", input_state.keys_str()); } } -void PlasmaPhase::getGibbs_RT(span grt) const +void PlasmaPhase::setState_TP(double t, double p) { - checkArraySize("PlasmaPhase::getGibbs_RT", grt.size(), m_kk); - auto gibbsrt = gibbs_RT_ref(); - copy(gibbsrt.begin(), gibbsrt.end(), grt.begin()); - double tmp = log(pressure() / refPressure()); - for (size_t k = 0; k < m_kk; k++) { - if (k != m_electronSpeciesIndex) { - grt[k] += tmp; + vector state(partialStateSize()); + savePartialState(state); + try { + setTemperature(t); + if (m_distributionType == "isotropic") { + setElectronTemperature(t); } else { - grt[k] += log(electronPressure() / refPressure()); + // Warning for other distribution types. + warn_user("PlasmaPhase::setState_TP", + "The electron temperature is not set equal to the gas temperature " + "when using '{}' electron energy distribution. " + "This may not be intended.", m_distributionType); } + setPressure(p); + } catch (std::exception&) { + restorePartialState(state); + throw; } } -double PlasmaPhase::electronMobility() const +void PlasmaPhase::setState_TgTeP(double Tg, double Te, double p) { - // Only implemented when using the Boltzmann two-term EEDF - if (m_distributionType == "Boltzmann-two-term") { - return m_eedfSolver->getElectronMobility(); - } else { - throw NotImplementedError("PlasmaPhase::electronMobility", - "Electron mobility is only available for 'Boltzmann-two-term' " - "electron energy distributions."); + vector state(partialStateSize()); + savePartialState(state); + try { + setTemperature(Tg); + if (m_distributionType == "isotropic") { + setElectronTemperature(Te); + } else { + // Warning for other distribution types. + warn_user("PlasmaPhase::setState_TgTeP", + "The electron temperature cannot be set " + "when using '{}' electron energy distribution. " + "This may not be intended.", m_distributionType); + } + setPressure(p); + } catch (std::exception&) { + restorePartialState(state); + throw; } } - -double PlasmaPhase::jouleHeatingPower() const +void PlasmaPhase::setState_TD(double t, double rho) { - // sigma = e * n_e * mu_e [S/m]; q_J = sigma * E^2 [W/m^3] - const double mu_e = electronMobility(); // m^2 / (V·s) - if (mu_e <= 0.0) { - return 0.0; - } - const double ne = concentration(m_electronSpeciesIndex) * Avogadro; // m^-3 - if (ne <= 0.0) { - return 0.0; - } - const double E = electricField(); // V/m - if (E <= 0.0) { - return 0.0; + vector state(partialStateSize()); + savePartialState(state); + try { + setTemperature(t); + if (m_distributionType == "isotropic") { + setElectronTemperature(t); + } else { + // Warning for other distribution types. + warn_user("PlasmaPhase::setState_TD", + "The electron temperature is not set equal to the gas temperature " + "when using '{}' electron energy distribution. " + "This may not be intended.", m_distributionType); + } + setDensity(rho); + } catch (std::exception&) { + restorePartialState(state); + throw; } - const double sigma = ElectronCharge * ne * mu_e; // S/m - return sigma * E * E; // W/m^3 } -double PlasmaPhase::intrinsicHeating() +void PlasmaPhase::setState_TgTeD(double Tg, double Te, double rho) { - // Joule heating: sigma * E^2 [W/m^3] - const double qJ = jouleHeatingPower(); - checkFinite(qJ); + vector state(partialStateSize()); + savePartialState(state); + try { + setTemperature(Tg); + if (m_distributionType == "isotropic") { + setElectronTemperature(Te); + } else { + // Warning for other distribution types. + warn_user("PlasmaPhase::setState_TgTeD", + "The electron temperature cannot be set " + "when using '{}' electron energy distribution. " + "This may not be intended.", m_distributionType); + } + setDensity(rho); + } catch (std::exception&) { + restorePartialState(state); + throw; + } +} - // Elastic + inelastic recoil power loss [W/m^3] - double qElastic = elasticPowerLoss(); - checkFinite(qElastic); - return qJ + qElastic; -} } diff --git a/test/data/consistency-cases.yaml b/test/data/consistency-cases.yaml index 151bcb4631c..845cd925da1 100644 --- a/test/data/consistency-cases.yaml +++ b/test/data/consistency-cases.yaml @@ -145,23 +145,29 @@ nitrogen-purefluid: plasma: setup: file: oxygen-plasma.yaml - phase: discretized-electron-energy-plasma + phase: isotropic-electron-energy-plasma known-failures: - g_eq_h_minus_Ts/.+: Test does not account for distinct electron temperature - u_eq_sum_uk_Xk/.+: Test does not account for distinct electron temperature - s_eq_sum_sk_Xk/.+: Test does not account for distinct electron temperature - cp_eq_dhdT/(0|2): Test does not account for distinct electron temperature - cv_eq_dudT/(0|2): Test does not account for distinct electron temperature - c_eq_sqrt_dP_drho_const_s/1: Test does not account for distinct electron temperature - cp_eq_.+/1: Test does not account for distinct electron temperature - cv_eq_.+/1: Test does not account for distinct electron temperature - c._eq_dsdT_const_._times_T: Test does not account for distinct electron temperature - cpk0_eq_dhk0dT: Test does not account for distinct electron temperature - cpRef_eq_dhRefdT: Test does not account for distinct electron temperature + c_eq_sqrt_dP_drho_const_s/.+: For now, `setState_SV` is not supported in PlasmaPhase. + cp_eq_dhdT/[345]: The differential methods are not yet implemented for two-temperature systems. + cv_eq_dudT/[345]: The differential methods are not yet implemented for two-temperature systems. + cp_eq_dsdT_const_p_times_T/[345]: + The differential methods are not yet implemented for two-temperature systems. + cv_eq_dsdT_const_v_times_T/[345]: + The differential methods are not yet implemented for two-temperature systems. + dSdv_const_T_eq_dPdT_const_V/[345]: + The differential methods are not yet implemented for two-temperature systems. + betaT_eq_minus_dmv_dP_const_T_div_mv/[345]: + The differential methods are not yet implemented for two-temperature systems. states: + # In the test cases below, it is implicitly assumed that the electron temperature + # equals the gas temperature. - {T: 300, P: 1 atm, X: {O2: 1.0, O2-: 1e-5, E: 1e-5}} - {T: 300, P: 1 atm, X: {E: 1.0}} - {T: 3500, P: 10 atm, X: {O2: 1.0, O2-: 2e-5, E: 2e-5}} + # In the test case below, set Te != T + - {T: 300, P: 1 atm, X: {O2: 1.0, O2-: 1e-5, E: 1e-5}, Te: 10000} + - {T: 300, P: 1 atm, X: {E: 1.0}, Te: 10000} + - {T: 3500, P: 10 atm, X: {O2: 1.0, O2-: 2e-5, E: 2e-5}, Te: 10000} debye-huckel-dilute: setup: diff --git a/test/data/thermo_plasma_simple.yaml b/test/data/thermo_plasma_simple.yaml new file mode 100644 index 00000000000..a77f046d15c --- /dev/null +++ b/test/data/thermo_plasma_simple.yaml @@ -0,0 +1,49 @@ +description: |- + This is a simple test of the plasma thermodynamic model, with a single plasma phase and a single gas phase. + The plasma phase has an isotropic electron energy distribution, and the gas phase is an ideal gas. + +phases: + - name: plasma + thermo: plasma + elements: [O, E] + species: [O, O+, e-] + kinetics: gas + state: { T: 300.0, P: 1 atm } + electron-energy-distribution: + type: isotropic + shape-factor: 1.0 + energy-levels: [0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0] + mean-electron-energy: 1.0 eV + reactions: none + +species: + - name: O + composition: { O: 1 } + thermo: + model: NASA9 + temperature-ranges: [200.0, 1000.0, 6000.0, 2.0e+04] + data: + - [0.0, 0.0, 3, 0.0, 0.0, 0.0, 0.0, 3e+04, 8.5] + - [0.0, 0.0, 3, 0.0, 0.0, 0.0, 0.0, 3e+04, 8.5] + - [0.0, 0.0, 3, 0.0, 0.0, 0.0, 0.0, 3e+04, 8.5] + note: False, but for testing purposes. + - name: O+ + composition: { O: 1, E: -1 } + thermo: + model: NASA9 + temperature-ranges: [298.15, 1000.0, 6000.0, 2.0e+04] + data: + - [0.0, 0.0, 8, 0.0, 0.0, 0.0, 0.0, 2e+05, 10.0] + - [0.0, 0.0, 8, 0.0, 0.0, 0.0, 0.0, 2e+05, 10.0] + - [0.0, 0.0, 8, 0.0, 0.0, 0.0, 0.0, 2e+05, 10.0] + note: False, but for testing purposes. + - name: e- + composition: { E: 1 } + thermo: + model: NASA9 + temperature-ranges: [298.15, 1000.0, 6000.0, 2.0e+04] + data: + - [0.0, 0.0, 2.5, 0.0, 0.0, 0.0, 0.0, -750, -11.7] + - [0.0, 0.0, 2.5, 0.0, 0.0, 0.0, 0.0, -750, -11.7] + - [0.0, 0.0, 2.5, 0.0, 0.0, 0.0, 0.0, -750, -11.7] + note: False, but for testing purposes. diff --git a/test/python/conftest.py b/test/python/conftest.py index e9a9f69cd4c..770c9de496a 100644 --- a/test/python/conftest.py +++ b/test/python/conftest.py @@ -10,8 +10,8 @@ pytest.register_assert_rewrite("pint.testing") pytest.register_assert_rewrite("python.utilities") -TEST_DATA_PATH = Path(__file__).parents[1] / "data" -CANTERA_DATA_PATH = Path(cantera.__file__).parent / "data" +TEST_DATA_PATH = (Path(__file__).parents[1] / "data").resolve() +CANTERA_DATA_PATH = (Path(cantera.__file__).parent / "data").resolve() def pytest_addoption(parser): diff --git a/test/python/test_thermo.py b/test/python/test_thermo.py index 7806d6bc9d5..cb1833eb9ab 100644 --- a/test/python/test_thermo.py +++ b/test/python/test_thermo.py @@ -1282,15 +1282,18 @@ def test_add_multiple_electron_species(self, phase): def test_elastic_power_loss_low_T(self, phase): phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 assert phase.elastic_power_loss == approx(6846332332) def test_elastic_power_loss_high_T(self, phase): # when T is as high as Te the energy loss rate becomes small phase.TPX = 4000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 assert phase.elastic_power_loss == approx(2865540) def test_elastic_power_loss_replace_rate(self, phase): phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 rate = ct.ReactionRate.from_dict(self.collision_data) phase.reaction(1).rate = rate assert phase.elastic_power_loss == approx(11765800095) @@ -1299,11 +1302,13 @@ def test_elastic_power_loss_add_reaction(self, phase): phase2 = ct.Solution(thermo="plasma", kinetics="bulk", species=phase.species(), reactions=[]) phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 phase.add_reaction(ct.Reaction.from_dict(self.collision_data, phase)) assert phase.elastic_power_loss == approx(18612132428) def test_elastic_power_loss_change_levels(self, phase): phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 phase.electron_energy_levels = np.linspace(0,10,101) assert phase.elastic_power_loss == approx(113058853) @@ -1326,6 +1331,7 @@ def test_elastic_power_loss_change_mean_electron_energy(self, phase): def test_elastic_power_loss_change_shape_factor(self, phase): phase.TPX = 1000, ct.one_atm, "O2:1, E:1e-5" + phase.mean_electron_energy = 1.0 phase.isotropic_shape_factor = 1.1 assert phase.elastic_power_loss == approx(7408711810) @@ -1364,6 +1370,78 @@ def test_eedf_solver(self): l2_norm = np.linalg.norm(interp - reference_eedf) assert l2_norm < 1e-3 + def _make_simple_plasma_phase(self, Te, set_te_first=False): + phase = ct.Solution("thermo_plasma_simple.yaml", "plasma", transport_model=None) + phase.isotropic_shape_factor = 1.0 + + T = 300.0 # K + P = ct.one_atm + composition = "O: 0.8, O+: 0.1, e-: 0.1" + + if set_te_first: + phase.Te = Te + phase.TPX = T, P, composition + if not set_te_first: + phase.Te = Te + + return phase + + @staticmethod + def _expected_cp_mole(X): + cp_o = 3.0 * ct.gas_constant + cp_o_plus = 8.0 * ct.gas_constant + cp_e = 2.5 * ct.gas_constant + return X[0] * cp_o + X[1] * cp_o_plus + X[2] * cp_e + + @staticmethod + def _expected_h_mole(X, T, Te): + h_o = 3.0 * ct.gas_constant * T + 3e4 * ct.gas_constant + h_o_plus = 8.0 * ct.gas_constant * T + 2e5 * ct.gas_constant + h_e = 2.5 * ct.gas_constant * Te - 750.0 * ct.gas_constant + return X[0] * h_o + X[1] * h_o_plus + X[2] * h_e + + @pytest.mark.parametrize( + "Te,set_te_first", + [ + (300.0, True), + (300.0, False), + (30_000.0, True), + (30_000.0, False), + ], + ) + def test_rho_two_temperature_modes(self, Te, set_te_first): + phase = self._make_simple_plasma_phase(Te, set_te_first=set_te_first) + T_ref = phase.mean_temperature + expected_rho = phase.P * phase.mean_molecular_weight / (ct.gas_constant * T_ref) + assert phase.density == approx(expected_rho) + + @pytest.mark.parametrize( + "Te,set_te_first", + [ + (300.0, True), + (300.0, False), + (30_000.0, True), + (30_000.0, False), + ], + ) + def test_cp_mole_two_temperature_modes(self, Te, set_te_first): + phase = self._make_simple_plasma_phase(Te, set_te_first=set_te_first) + assert phase.cp_mole == approx(self._expected_cp_mole(phase.X)) + + @pytest.mark.parametrize( + "Te,set_te_first", + [ + (300.0, True), + (300.0, False), + (30_000.0, True), + (30_000.0, False), + ], + ) + def test_enthalpy_mole_two_temperature_modes(self, Te, set_te_first): + phase = self._make_simple_plasma_phase(Te, set_te_first=set_te_first) + expected_h = self._expected_h_mole(phase.X, phase.T, phase.Te) + assert phase.enthalpy_mole == approx(expected_h) + class TestImport: """ diff --git a/test/thermo_consistency/consistency.cpp b/test/thermo_consistency/consistency.cpp index 5a863a3fcc2..7508f6cba94 100644 --- a/test/thermo_consistency/consistency.cpp +++ b/test/thermo_consistency/consistency.cpp @@ -1,956 +1,1005 @@ -// This file is part of Cantera. See License.txt in the top-level directory or -// at https://cantera.org/license.txt for license and copyright information. - -#include "gtest/gtest.h" -#include "cantera/thermo/ThermoPhase.h" -#include "cantera/thermo/PlasmaPhase.h" -#include "cantera/thermo/ThermoFactory.h" -#include "cantera/thermo/MolalityVPSSTP.h" -#include "cantera/base/Solution.h" -#include "cantera/base/utilities.h" -#include -#include - -// This is a set of tests that check all ThermoPhase models implemented in Cantera -// for thermodynamic identities such as g = h - T*s and c_p = dh/dT at constant P. - -// Below the definition of the TestConsistency test fixture class, the individual -// consistency tests are implemented in the functions declared using the TEST_P macro. -// Each of these tests starts with the 'phase' object created and set to the -// thermodynamic state that should be used for the test. - -// Following the individual test definitions, a "test suite" is created for each phase -// model using the INSTANTIATE_TEST_SUITE_P macro. The string passed to the getSetup() -// and getStates() functions in the test suite instantiation corresponds to a top-level -// section of the file test/data/consistency-cases.yaml. Each of these sections defines -// a phase definition and a sequence of thermodynamic states to be used for each of the -// individual test functions. Generally, there is one test suite for each ThermoPhase -// model. However, there are a few phase models that have multiple test suites to allow -// testing the effects of parameters that can only be varied within the phase -// definition. - -using namespace std; -using namespace Cantera; - -// Helper functions to reduce code duplication in test suite instantiations -vector getStates(const string& name) { - static AnyMap cases = AnyMap::fromYamlFile("consistency-cases.yaml"); - return cases[name]["states"].asVector(); -} - -AnyMap getSetup(const string& name) { - static AnyMap cases = AnyMap::fromYamlFile("consistency-cases.yaml"); - return cases[name]["setup"].as(); -} - -namespace Cantera { - -// For more informative output about failing test cases -std::ostream& operator<<(std::ostream& s, const AnyMap& m) -{ - if (m.hasKey("phase")) { - s << fmt::format("file: {}, phase: {}", - m["file"].asString(), m["phase"].asString()); - } else if (m.hasKey("file")) { - s << fmt::format("file: {}", m["file"].asString()); - } else { - AnyMap out = m; - out.setFlowStyle(); - s << "state: " - << boost::algorithm::replace_all_copy(out.toYamlString(), "\n", " "); - } - return s; -} - -} // end namespace Cantera - -class TestConsistency : public testing::TestWithParam> -{ -public: - TestConsistency() { - auto param = GetParam(); - setup = get<0>(param); - AnyMap state = get<1>(param); - - if (setup.getBool("ignore-deprecations", false)) { - suppress_deprecation_warnings(); - } - - // For efficiency, cache the instantiated phase object rather than recreating - // it for every single test case. - pair key = {setup["file"].asString(), setup.getString("phase", "")}; - if (cache.count(key) == 0) { - cache[key] = newThermo(key.first, key.second); - } - atol = setup.getDouble("atol", 1e-5); - rtol_fd = setup.getDouble("rtol_fd", 1e-6); - atol_v = setup.getDouble("atol_v", 1e-11); - atol_c = setup.getDouble("atol_c", 1e-14); - atol_e = setup.getDouble("atol_e", 1e-18); - - phase = cache[key]; - phase->setState(state); - nsp = phase->nSpecies(); - p = phase->pressure(); - T = phase->temperature(); - Te = phase->electronTemperature(); - RTe = Te * GasConstant; - RT = T * GasConstant; - if (phase->type() == "plasma") { - ke = dynamic_cast(*phase).electronSpeciesIndex(); - } else { - ke = npos; - } - } - - ~TestConsistency() { - make_deprecation_warnings_fatal(); - } - - void SetUp() override { - // See if we should skip this test specific test case - if (setup.hasKey("known-failures")) { - auto current = testing::UnitTest::GetInstance()->current_test_info()->name(); - for (auto& [pattern, reason] : setup["known-failures"].asMap()) { - if (regex_search(current, regex(pattern))) { - GTEST_SKIP() << reason; - break; - } - } - } - } - - static map, shared_ptr> cache; - - AnyMap setup; - shared_ptr phase; - size_t nsp; - double T, p, RT, RTe, Te; - size_t ke; - double atol; // absolute tolerance for molar energy comparisons - double atol_v; // absolute tolerance for molar volume comparisons - double atol_c; // absolute tolerance for compressibility comparison - double atol_e; // absolute tolerance for expansion comparison - double rtol_fd; // relative tolerance for finite difference comparisons -}; - -map, shared_ptr> TestConsistency::cache = {}; - -// --------------- Definitions for individual consistency tests --------------- - -TEST_P(TestConsistency, h_eq_u_plus_Pv) { - double h, u, v; - try { - h = phase->enthalpy_mole(); - u = phase->intEnergy_mole(); - v = phase->molarVolume(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - if (phase->type() == "plasma" && ke != npos) { - // Two-temperature identity for PlasmaPhase: - // h = u + p*v + X_e * R * (Te - T) - vector X(nsp); - phase->getMoleFractions(X); - double Xe = X[ke]; - EXPECT_NEAR(h, u + p * v + Xe * (RTe - RT), atol); - } else { - // Standard single-temperature identity - EXPECT_NEAR(h, u + p * v, atol); - } -} - -TEST_P(TestConsistency, g_eq_h_minus_Ts) { - double g, h, s; - try { - g = phase->gibbs_mole(); - h = phase->enthalpy_mole(); - s = phase->entropy_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - EXPECT_NEAR(g, h - T * s, atol); -} - -TEST_P(TestConsistency, hk_eq_uk_plus_P_vk) -{ - vector hk(nsp), uk(nsp), vk(nsp); - try { - phase->getPartialMolarEnthalpies(hk); - phase->getPartialMolarIntEnergies(uk); - phase->getPartialMolarVolumes(vk); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - if (k != ke) { - EXPECT_NEAR(hk[k], uk[k] + p * vk[k], atol) << "k = " << k; - } // not applicable for electron - } -} - -TEST_P(TestConsistency, utilde_eq_uk_minus_piT_vk) -{ - vector utilde(nsp), uk(nsp), vk(nsp); - double piT; - try { - phase->getPartialMolarIntEnergies_TV(utilde); - phase->getPartialMolarIntEnergies(uk); - phase->getPartialMolarVolumes(vk); - piT = phase->internalPressure(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - if (!std::isfinite(piT)) { - GTEST_SKIP() << "internalPressure is not finite at this state"; - } - for (size_t k = 0; k < nsp; k++) { - if (k != ke) { - double expected = uk[k] - piT * vk[k]; - double tol = max({atol, rtol_fd * abs(utilde[k]), rtol_fd * abs(expected)}); - EXPECT_NEAR(utilde[k], expected, tol) << "k = " << k; - } - } -} - -TEST_P(TestConsistency, piT_eq_dudv_const_T) -{ - if (!phase->isCompressible()) { - GTEST_SKIP() << "Undefined for incompressible phase"; - } - - double piT; - double T0 = phase->temperature(); - double v0 = 1.0 / phase->density(); // specific volume [m^3/kg] - // Use a moderate perturbation to avoid cancellation in u(v) differences. - double dv = 1e-4 * v0; - double v1 = v0 - dv; - double v2 = v0 + dv; - - try { - piT = phase->internalPressure(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - - phase->setState_TD(T0, 1.0 / v1); - double u1 = phase->intEnergy_mass(); - phase->setState_TD(T0, 1.0 / v2); - double u2 = phase->intEnergy_mass(); - - double piT_fd = (u2 - u1) / (v2 - v1); - double tol = max({rtol_fd * abs(piT), rtol_fd * abs(piT_fd), atol}); - EXPECT_NEAR(piT_fd, piT, tol); -} - -TEST_P(TestConsistency, gk_eq_hk_minus_T_sk) -{ - vector gk(nsp), hk(nsp), sk(nsp); - try { - phase->getChemPotentials(gk); - phase->getPartialMolarEnthalpies(hk); - phase->getPartialMolarEntropies(sk); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - if (k != ke) { - EXPECT_NEAR(gk[k], hk[k] - T * sk[k], atol) << "k = " << k; - } // not applicable for electron - } -} - -TEST_P(TestConsistency, h_eq_sum_hk_Xk) -{ - vector hk(nsp); - try { - phase->getPartialMolarEnthalpies(hk); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - EXPECT_NEAR(phase->enthalpy_mole(), phase->mean_X(hk), atol); -} - -TEST_P(TestConsistency, u_eq_sum_uk_Xk) -{ - vector uk(nsp); - double u; - try { - phase->getPartialMolarIntEnergies(uk); - u = phase->intEnergy_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - EXPECT_NEAR(u, phase->mean_X(uk), atol); -} - -TEST_P(TestConsistency, g_eq_sum_gk_Xk) -{ - vector gk(nsp); - double g; - try { - phase->getChemPotentials(gk); - g = phase->gibbs_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - EXPECT_NEAR(g, phase->mean_X(gk), atol); -} - -TEST_P(TestConsistency, s_eq_sum_sk_Xk) -{ - vector sk(nsp); - double s; - try { - phase->getPartialMolarEntropies(sk); - s = phase->entropy_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - EXPECT_NEAR(s, phase->mean_X(sk), atol); -} - -TEST_P(TestConsistency, v_eq_sum_vk_Xk) -{ - vector vk(nsp); - try { - phase->getPartialMolarVolumes(vk); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - EXPECT_NEAR(phase->molarVolume(), phase->mean_X(vk), atol_v); -} - -TEST_P(TestConsistency, cp_eq_sum_cpk_Xk) -{ - vector cpk(nsp); - double cp; - try { - phase->getPartialMolarCp(cpk); - cp = phase->cp_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - EXPECT_NEAR(cp, phase->mean_X(cpk), atol); -} - -TEST_P(TestConsistency, cp_eq_dhdT) -{ - double h1, cp1; - try { - h1 = phase->enthalpy_mole(); - cp1 = phase->cp_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - double T1 = phase->temperature(); - double dT = 1e-5 * phase->temperature(); - phase->setState_TP(T1 + dT, phase->pressure()); - double h2 = phase->enthalpy_mole(); - double cp2 = phase->cp_mole(); - double cp_mid = 0.5 * (cp1 + cp2); - double cp_fd = (h2 - h1) / dT; - EXPECT_NEAR(cp_fd, cp_mid, max({rtol_fd * cp_mid, rtol_fd * cp_fd, atol})); -} - -TEST_P(TestConsistency, cv_eq_dudT) -{ - double u1, cv1; - try { - u1 = phase->intEnergy_mole(); - cv1 = phase->cv_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - double T1 = phase->temperature(); - double dT = 1e-5 * phase->temperature(); - if (phase->isCompressible()) { - phase->setState_TD(T1 + dT, phase->density()); - } else { - phase->setTemperature(T1 + dT); - } - double u2 = phase->intEnergy_mole(); - double cv2 = phase->cv_mole(); - double cv_mid = 0.5 * (cv1 + cv2); - double cv_fd = (u2 - u1) / dT; - EXPECT_NEAR(cv_fd, cv_mid, max({rtol_fd * cv_mid, rtol_fd * cv_fd, atol})); -} - -TEST_P(TestConsistency, cp_eq_dsdT_const_p_times_T) -{ - double s1, cp1; - try { - s1 = phase->entropy_mole(); - cp1 = phase->cp_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - double T1 = phase->temperature(); - double dT = 1e-4 * phase->temperature(); - phase->setState_TP(T1 + dT, phase->pressure()); - double s2 = phase->entropy_mole(); - double cp2 = phase->cp_mole(); - double cp_mid = 0.5 * (cp1 + cp2); - double cp_fd = (T1 + dT/2) * (s2 - s1) / dT; - EXPECT_NEAR(cp_fd, cp_mid, max({rtol_fd * cp_mid, rtol_fd * cp_fd, atol})); -} - -TEST_P(TestConsistency, cv_eq_dsdT_const_v_times_T) -{ - double s1, cv1; - try { - s1 = phase->entropy_mole(); - cv1 = phase->cv_mole(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - double T1 = phase->temperature(); - double dT = 1e-4 * phase->temperature(); - if (phase->isCompressible()) { - phase->setState_TD(T1 + dT, phase->density()); - } else { - phase->setTemperature(T1 + dT); - } - double s2 = phase->entropy_mole(); - double cv2 = phase->cv_mole(); - double cv_mid = 0.5 * (cv1 + cv2); - double cv_fd = (T1 + dT/2) * (s2 - s1) / dT; - EXPECT_NEAR(cv_fd, cv_mid, max({rtol_fd * cv_mid, rtol_fd * cv_fd, atol})); -} - -TEST_P(TestConsistency, dsdP_const_T_eq_minus_dV_dT_const_P) -{ - double P0 = phase->pressure(); - double T0 = phase->temperature(); - double s1, v1; - double dP = 1e-4 * P0; - double dT = 1e-4 * T0; - try { - phase->setState_TP(T0, P0 - dP); - s1 = phase->entropy_mole(); - phase->setState_TP(T0 - dT, P0); - v1 = phase->molarVolume(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - phase->setState_TP(T0, P0 + dP); - double s2 = phase->entropy_mole(); - double dsdP = (s2 - s1) / (2 * dP); - - phase->setState_TP(T0 + dT, P0); - double v2 = phase->molarVolume(); - double dvdT = (v2 - v1) / (2 * dT); - double tol = rtol_fd * std::max(std::abs(dsdP), std::abs(dvdT)); - EXPECT_NEAR(dsdP, -dvdT, tol); -} - -TEST_P(TestConsistency, dSdv_const_T_eq_dPdT_const_V) { - if (phase->isCompressible()) { - double s1; - try { - s1 = phase->entropy_mass(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - double v1 = 1 / phase->density(); - double P1 = phase->pressure(); - double v2 = v1 * (1 + 1e-7); - phase->setState_TD(T, 1 / v2); - double s2 = phase->entropy_mass(); - double dsdv = (s2 - s1) / (v2 - v1); - - double T2 = T * (1 + 1e-7); - phase->setState_TD(T2, 1 / v1); - double P2 = phase->pressure(); - double dPdT = (P2 - P1) / (T2 - T); - double tol = rtol_fd * std::max(std::abs(dPdT), std::abs(dsdv)); - EXPECT_NEAR(dsdv, dPdT, tol); - } else { - GTEST_SKIP() << "Undefined for incompressible phase"; - } -} - -TEST_P(TestConsistency, betaT_eq_minus_dmv_dP_const_T_div_mv) -{ - double betaT1; - try { - betaT1 = phase->isothermalCompressibility(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - - double T = phase->temperature(); - double P1 = phase->pressure(); - double mv1 = phase->molarVolume(); - - double P2 = P1 * (1 + 1e-6); - phase->setState_TP(T, P2); - double betaT2 = phase->isothermalCompressibility(); - double mv2 = phase->molarVolume(); - - double betaT_mid = 0.5 * (betaT1 + betaT2); - double mv_mid = 0.5 * (mv1 + mv2); - double betaT_fd = -1 / mv_mid * (mv2 - mv1) / (P2 - P1); - - EXPECT_NEAR(betaT_fd, betaT_mid, - max({rtol_fd * betaT_mid, rtol_fd * betaT_fd, atol_c})); -} - -TEST_P(TestConsistency, alphaV_eq_dmv_dT_const_P_div_mv) -{ - double alphaV1; - try { - alphaV1 = phase->thermalExpansionCoeff(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - - double P = phase->pressure(); - double T1 = phase->temperature(); - double mv1 = phase->molarVolume(); - - double T2 = T1 * (1 + 1e-6); - phase->setState_TP(T2, P); - double alphaV2 = phase->thermalExpansionCoeff(); - double mv2 = phase->molarVolume(); - - double alphaV_mid = 0.5 * (alphaV1 + alphaV2); - double mv_mid = 0.5 * (mv1 + mv2); - double alphaV_fd = 1 / mv_mid * (mv2 - mv1) / (T2 - T1); - - EXPECT_NEAR(alphaV_fd, alphaV_mid, - max({rtol_fd * alphaV_mid, rtol_fd * alphaV_fd, atol_e})); -} - -TEST_P(TestConsistency, c_eq_sqrt_dP_drho_const_s) -{ - double c1; - try { - c1 = phase->soundSpeed(); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - - double P1 = phase->pressure(); - double rho1 = phase->density(); - - double rho2 = rho1 * (1 + 1e-6); - phase->setState_SV(phase->entropy_mass(), 1 / rho2); - double c2 = phase->soundSpeed(); - double P2 = phase->pressure(); - - double c_mid = 0.5 * (c1 + c2); - double c_fd = sqrt((P2 - P1) / (rho2 - rho1)); - - EXPECT_NEAR(c_fd, c_mid, max({rtol_fd * c_mid, rtol_fd * c_fd, atol})); -} - -// ---------- Tests for consistency of standard state properties --------------- - -TEST_P(TestConsistency, hk0_eq_uk0_plus_p_vk0) -{ - vector h0(nsp), u0(nsp), v0(nsp); - try { - phase->getEnthalpy_RT(h0); - phase->getIntEnergy_RT(u0); - phase->getStandardVolumes(v0); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - EXPECT_NEAR(h0[k] * RT, u0[k] * RT + p * v0[k], atol) << "k = " << k; - } -} - -TEST_P(TestConsistency, gk0_eq_hk0_minus_T_sk0) -{ - vector g0(nsp), h0(nsp), s0(nsp); - try { - phase->getEnthalpy_RT(h0); - phase->getGibbs_RT(g0); - phase->getEntropy_R(s0); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - EXPECT_NEAR(g0[k] * RT , - h0[k] * RT - T * s0[k] * GasConstant, atol) << "k = " << k; - } -} - -TEST_P(TestConsistency, cpk0_eq_dhk0dT) -{ - vector h1(nsp), h2(nsp), cp1(nsp), cp2(nsp); - try { - phase->getEnthalpy_RT(h1); - phase->getCp_R(cp1); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - double T1 = phase->temperature(); - double dT = 1e-5 * phase->temperature(); - phase->setState_TP(T1 + dT, phase->pressure()); - phase->getEnthalpy_RT(h2); - phase->getCp_R(cp2); - for (size_t k = 0; k < nsp; k++) { - double cp_mid = 0.5 * (cp1[k] + cp2[k]) * GasConstant; - double cp_fd = (h2[k] * (T1 + dT) - h1[k] * T1) / dT * GasConstant; - double tol = max({rtol_fd * std::abs(cp_mid), rtol_fd * std::abs(cp_fd), atol}); - EXPECT_NEAR(cp_fd, cp_mid, tol) << "k = " << k; - } -} - -TEST_P(TestConsistency, standard_gibbs_nondim) -{ - vector g0_RT(nsp), mu0(nsp); - try { - phase->getGibbs_RT(g0_RT); - phase->getStandardChemPotentials(mu0); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - if (k != ke) { - EXPECT_NEAR(g0_RT[k] * RT , mu0[k], atol) << "k = " << k; - } else { - EXPECT_NEAR(g0_RT[k] * RTe , mu0[k], atol) << "k = " << k; - } - } -} - -TEST_P(TestConsistency, chem_potentials_to_activities) { - vector mu0(nsp), mu(nsp), a(nsp); - try { - phase->getChemPotentials(mu); - phase->getStandardChemPotentials(mu0); - phase->getActivities(a); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - if (k != ke) { - double a_from_mu = exp((mu[k] - mu0[k]) / RT); - double scale = std::max(std::abs(a[k]), std::abs(a_from_mu)); - EXPECT_NEAR(a_from_mu, a[k], 1e-9 * scale + 1e-14) << "k = " << k; - } else { - double a_from_mu = exp((mu[k] - mu0[k]) / RTe); - double scale = std::max(std::abs(a[k]), std::abs(a_from_mu)); - EXPECT_NEAR(a_from_mu, a[k], 1e-9 * scale + 1e-14) << "k = " << k; - } - } -} - -TEST_P(TestConsistency, activity_coeffs) { - vector a(nsp), gamma(nsp), X(nsp); - try { - phase->getActivities(a); - phase->getActivityCoefficients(gamma); - phase->getMoleFractions(X); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - double scale = std::max(std::abs(a[k]), std::abs(gamma[k] * X[k])); - EXPECT_NEAR(a[k], gamma[k] * X[k], 1e-10 * scale + 1e-20) << "k = " << k; - } -} - -TEST_P(TestConsistency, activity_concentrations) { - vector a(nsp), Cact(nsp); - try { - phase->getActivities(a); - phase->getActivityConcentrations(Cact); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - double C0k = phase->standardConcentration(k); - EXPECT_NEAR(a[k], Cact[k] / C0k, 1e-9 * std::abs(a[k])) << "k = " << k; - } -} - -TEST_P(TestConsistency, log_standard_concentrations) { - for (size_t k = 0; k < nsp; k++) { - double c0 = phase->standardConcentration(k); - EXPECT_NEAR(c0, exp(phase->logStandardConc(k)), 1e-10 * c0) << "k = " << k; - } -} - -TEST_P(TestConsistency, log_activity_coeffs) { - vector gamma(nsp), log_gamma(nsp); - try { - phase->getActivityCoefficients(gamma); - phase->getLnActivityCoefficients(log_gamma); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - EXPECT_NEAR(gamma[k], exp(log_gamma[k]), 1e-10 * gamma[k]) << "k = " << k; - } -} - -// -------------------- Tests for reference state properties ------------------- - -TEST_P(TestConsistency, hRef_eq_uRef_plus_P_vRef) -{ - vector hRef(nsp), uRef(nsp), vRef(nsp); - try { - phase->getEnthalpy_RT_ref(hRef); - phase->getIntEnergy_RT_ref(uRef); - phase->getStandardVolumes_ref(vRef); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - if (k != ke) { - EXPECT_NEAR(hRef[k] * RT, uRef[k] * RT + OneAtm * vRef[k], atol) << "k = " << k; - } else { - EXPECT_NEAR(hRef[k] * RTe, uRef[k] * RTe + OneAtm * vRef[k], atol) << "k = " << k; - } - } -} - -TEST_P(TestConsistency, gRef_eq_hRef_minus_T_sRef) -{ - vector hRef(nsp), gRef_RT(nsp), gRef(nsp), sRef(nsp); - try { - phase->getEnthalpy_RT_ref(hRef); - phase->getGibbs_ref(gRef); - phase->getGibbs_RT_ref(gRef_RT); - phase->getEntropy_R_ref(sRef); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - for (size_t k = 0; k < nsp; k++) { - if (k != ke) { - EXPECT_NEAR(gRef[k], gRef_RT[k] * RT, atol) << "k = " << k; - EXPECT_NEAR(gRef[k], hRef[k] * RT - T * sRef[k] * GasConstant, - atol) << "k = " << k; - } else { - EXPECT_NEAR(gRef[k], gRef_RT[k] * RTe, atol) << "k = " << k; - EXPECT_NEAR(gRef[k], hRef[k] * RTe - Te * sRef[k] * GasConstant, - atol) << "k = " << k; - } - } -} - -TEST_P(TestConsistency, cpRef_eq_dhRefdT) -{ - vector h1(nsp), h2(nsp), cp1(nsp), cp2(nsp); - try { - phase->getEnthalpy_RT_ref(h1); - phase->getCp_R_ref(cp1); - } catch (NotImplementedError& err) { - GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; - } - double T1 = phase->temperature(); - double dT = 1e-5 * phase->temperature(); - phase->setState_TP(T1 + dT, phase->pressure()); - phase->getEnthalpy_RT_ref(h2); - phase->getCp_R_ref(cp2); - for (size_t k = 0; k < nsp; k++) { - double cp_mid = 0.5 * (cp1[k] + cp2[k]) * GasConstant; - double cp_fd = (h2[k] * (T1 + dT) - h1[k] * T1) / dT * GasConstant; - double tol = max({rtol_fd * std::abs(cp_mid), rtol_fd * std::abs(cp_fd), atol}); - EXPECT_NEAR(cp_fd, cp_mid, tol) << "k = " << k; - } -} - -// ------- Instantiation of test suites defined in consistency-cases.yaml ------ - -INSTANTIATE_TEST_SUITE_P(IdealGas, TestConsistency, - testing::Combine( - testing::Values(getSetup("ideal-gas-h2o2")), - testing::ValuesIn(getStates("ideal-gas-h2o2"))) -); - -INSTANTIATE_TEST_SUITE_P(RedlichKwong, TestConsistency, - testing::Combine( - testing::Values(getSetup("redlich-kwong")), - testing::ValuesIn(getStates("redlich-kwong"))) -); - -INSTANTIATE_TEST_SUITE_P(PengRobinson, TestConsistency, - testing::Combine( - testing::Values(getSetup("peng-robinson")), - testing::ValuesIn(getStates("peng-robinson"))) -); - -INSTANTIATE_TEST_SUITE_P(IdealMolalSolution, TestConsistency, - testing::Combine( - testing::Values(getSetup("ideal-molal-solution")), - testing::ValuesIn(getStates("ideal-molal-solution"))) -); - -INSTANTIATE_TEST_SUITE_P(IdealSolidSolnPhase1, TestConsistency, - testing::Combine( - testing::Values(getSetup("ideal-condensed-1")), - testing::ValuesIn(getStates("ideal-condensed-1"))) -); - -INSTANTIATE_TEST_SUITE_P(IdealSolidSolnPhase2, TestConsistency, - testing::Combine( - testing::Values(getSetup("ideal-condensed-2")), - testing::ValuesIn(getStates("ideal-condensed-2"))) -); - -INSTANTIATE_TEST_SUITE_P(BinarySolutionTabulated, TestConsistency, - testing::Combine( - testing::Values(getSetup("binary-solution-tabulated")), - testing::ValuesIn(getStates("binary-solution-tabulated"))) -); - -INSTANTIATE_TEST_SUITE_P(ElectronCloud, TestConsistency, - testing::Combine( - testing::Values(getSetup("electron-cloud")), - testing::ValuesIn(getStates("electron-cloud"))) -); - -INSTANTIATE_TEST_SUITE_P(NitrogenPureFluid, TestConsistency, - testing::Combine( - testing::Values(getSetup("nitrogen-purefluid")), - testing::ValuesIn(getStates("nitrogen-purefluid"))) -); - -INSTANTIATE_TEST_SUITE_P(PlasmaPhase, TestConsistency, - testing::Combine( - testing::Values(getSetup("plasma")), - testing::ValuesIn(getStates("plasma"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckelDilute, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-dilute")), - testing::ValuesIn(getStates("debye-huckel-dilute"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckelDilute_IAPWS, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-dilute-IAPWS")), - testing::ValuesIn(getStates("debye-huckel-dilute-IAPWS"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckel_b_dot_ak, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-B-dot-ak")), - testing::ValuesIn(getStates("debye-huckel-B-dot-ak"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckel_b_dot_ak_IAPWS, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-B-dot-ak-IAPWS")), - testing::ValuesIn(getStates("debye-huckel-B-dot-ak-IAPWS"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckel_b_dot_a, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-B-dot-a")), - testing::ValuesIn(getStates("debye-huckel-B-dot-a"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckel_b_dot_a_IAPWS, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-B-dot-a-IAPWS")), - testing::ValuesIn(getStates("debye-huckel-B-dot-a-IAPWS"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckel_pitzer_beta_ij, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-pitzer-beta_ij")), - testing::ValuesIn(getStates("debye-huckel-pitzer-beta_ij"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckel_pitzer_beta_ij_IAPWS, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-pitzer-beta_ij-IAPWS")), - testing::ValuesIn(getStates("debye-huckel-pitzer-beta_ij-IAPWS"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckel_beta_ij, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-beta_ij")), - testing::ValuesIn(getStates("debye-huckel-beta_ij"))) -); - -INSTANTIATE_TEST_SUITE_P(DebyeHuckel_beta_ij_IAPWS, TestConsistency, - testing::Combine( - testing::Values(getSetup("debye-huckel-beta_ij-IAPWS")), - testing::ValuesIn(getStates("debye-huckel-beta_ij-IAPWS"))) -); - -INSTANTIATE_TEST_SUITE_P(Margules, TestConsistency, - testing::Combine( - testing::Values(getSetup("margules")), - testing::ValuesIn(getStates("margules"))) -); - -INSTANTIATE_TEST_SUITE_P(FixedStoichiometry, TestConsistency, - testing::Combine( - testing::Values(getSetup("fixed-stoichiometry")), - testing::ValuesIn(getStates("fixed-stoichiometry"))) -); - -INSTANTIATE_TEST_SUITE_P(IdealSurface, TestConsistency, - testing::Combine( - testing::Values(getSetup("ideal-surface")), - testing::ValuesIn(getStates("ideal-surface"))) -); - -INSTANTIATE_TEST_SUITE_P(IdealEdge, TestConsistency, - testing::Combine( - testing::Values(getSetup("ideal-edge")), - testing::ValuesIn(getStates("ideal-edge"))) -); - -INSTANTIATE_TEST_SUITE_P(CoverageDependentSurface, TestConsistency, - testing::Combine( - testing::Values(getSetup("coverage-dependent-surface")), - testing::ValuesIn(getStates("coverage-dependent-surface"))) -); - -INSTANTIATE_TEST_SUITE_P(LiquidWaterIapws95, TestConsistency, - testing::Combine( - testing::Values(getSetup("liquid-water-IAPWS95")), - testing::ValuesIn(getStates("liquid-water-IAPWS95"))) -); - -INSTANTIATE_TEST_SUITE_P(IdealSolnVPSS_simple, TestConsistency, - testing::Combine( - testing::Values(getSetup("ideal-solution-VPSS-simple")), - testing::ValuesIn(getStates("ideal-solution-VPSS-simple"))) -); - -INSTANTIATE_TEST_SUITE_P(IdealSolnVPSS_HKFT, TestConsistency, - testing::Combine( - testing::Values(getSetup("ideal-solution-VPSS-HKFT")), - testing::ValuesIn(getStates("ideal-solution-VPSS-HKFT"))) -); - -INSTANTIATE_TEST_SUITE_P(RedlichKister_LiC6, TestConsistency, - testing::Combine( - testing::Values(getSetup("Redlich-Kister-LiC6")), - testing::ValuesIn(getStates("Redlich-Kister-LiC6"))) -); - -INSTANTIATE_TEST_SUITE_P(RedlichKister_complex, TestConsistency, - testing::Combine( - testing::Values(getSetup("Redlich-Kister-complex")), - testing::ValuesIn(getStates("Redlich-Kister-complex"))) -); - -INSTANTIATE_TEST_SUITE_P(HMWSoln, TestConsistency, - testing::Combine( - testing::Values(getSetup("HMW-electrolyte")), - testing::ValuesIn(getStates("HMW-electrolyte"))) -); - -int main(int argc, char** argv) -{ - printf("Running main() from consistency.cpp\n"); - testing::InitGoogleTest(&argc, argv); - Cantera::make_deprecation_warnings_fatal(); - Cantera::CanteraError::setStackTraceDepth(0); - int result = RUN_ALL_TESTS(); - Cantera::appdelete(); - return result; -} +// This file is part of Cantera. See License.txt in the top-level directory or +// at https://cantera.org/license.txt for license and copyright information. + +#include "gtest/gtest.h" +#include "cantera/thermo/ThermoPhase.h" +#include "cantera/thermo/PlasmaPhase.h" +#include "cantera/thermo/ThermoFactory.h" +#include "cantera/thermo/MolalityVPSSTP.h" +#include "cantera/base/Solution.h" +#include "cantera/base/utilities.h" +#include +#include + +// This is a set of tests that check all ThermoPhase models implemented in Cantera +// for thermodynamic identities such as g = h - T*s and c_p = dh/dT at constant P. + +// Below the definition of the TestConsistency test fixture class, the individual +// consistency tests are implemented in the functions declared using the TEST_P macro. +// Each of these tests starts with the 'phase' object created and set to the +// thermodynamic state that should be used for the test. + +// Following the individual test definitions, a "test suite" is created for each phase +// model using the INSTANTIATE_TEST_SUITE_P macro. The string passed to the getSetup() +// and getStates() functions in the test suite instantiation corresponds to a top-level +// section of the file test/data/consistency-cases.yaml. Each of these sections defines +// a phase definition and a sequence of thermodynamic states to be used for each of the +// individual test functions. Generally, there is one test suite for each ThermoPhase +// model. However, there are a few phase models that have multiple test suites to allow +// testing the effects of parameters that can only be varied within the phase +// definition. + +using namespace std; +using namespace Cantera; + +// Helper functions to reduce code duplication in test suite instantiations +vector getStates(const string& name) { + static AnyMap cases = AnyMap::fromYamlFile("consistency-cases.yaml"); + return cases[name]["states"].asVector(); +} + +AnyMap getSetup(const string& name) { + static AnyMap cases = AnyMap::fromYamlFile("consistency-cases.yaml"); + return cases[name]["setup"].as(); +} + +namespace Cantera { + +// For more informative output about failing test cases +std::ostream& operator<<(std::ostream& s, const AnyMap& m) +{ + if (m.hasKey("phase")) { + s << fmt::format("file: {}, phase: {}", + m["file"].asString(), m["phase"].asString()); + } else if (m.hasKey("file")) { + s << fmt::format("file: {}", m["file"].asString()); + } else { + AnyMap out = m; + out.setFlowStyle(); + s << "state: " + << boost::algorithm::replace_all_copy(out.toYamlString(), "\n", " "); + } + return s; +} + +} // end namespace Cantera + +class TestConsistency : public testing::TestWithParam> +{ +public: + TestConsistency() { + auto param = GetParam(); + setup = get<0>(param); + AnyMap state = get<1>(param); + + if (setup.getBool("ignore-deprecations", false)) { + suppress_deprecation_warnings(); + } + + // For efficiency, cache the instantiated phase object rather than recreating + // it for every single test case. + pair key = {setup["file"].asString(), setup.getString("phase", "")}; + if (cache.count(key) == 0) { + cache[key] = newThermo(key.first, key.second); + } + atol = setup.getDouble("atol", 1e-5); + rtol_fd = setup.getDouble("rtol_fd", 1e-6); + atol_v = setup.getDouble("atol_v", 1e-11); + atol_c = setup.getDouble("atol_c", 1e-14); + atol_e = setup.getDouble("atol_e", 1e-18); + + phase = cache[key]; + phase->setState(state); + nsp = phase->nSpecies(); + p = phase->pressure(); + T = phase->temperature(); + Te = phase->electronTemperature(); + + if (phase->type() == "plasma") { + ke = dynamic_cast(*phase).electronSpeciesIndex(); + } else { + ke = npos; + } + RTe = Te * GasConstant; + RT = T * GasConstant; + } + + ~TestConsistency() { + make_deprecation_warnings_fatal(); + } + + void SetUp() override { + // See if we should skip this test specific test case + if (setup.hasKey("known-failures")) { + auto current = testing::UnitTest::GetInstance()->current_test_info()->name(); + for (auto& [pattern, reason] : setup["known-failures"].asMap()) { + if (regex_search(current, regex(pattern))) { + GTEST_SKIP() << reason; + break; + } + } + } + } + + static map, shared_ptr> cache; + + AnyMap setup; + shared_ptr phase; + size_t nsp; + double T, p, RT, RTe, Te; + size_t ke; + double atol; // absolute tolerance for molar energy comparisons + double atol_v; // absolute tolerance for molar volume comparisons + double atol_c; // absolute tolerance for compressibility comparison + double atol_e; // absolute tolerance for expansion comparison + double rtol_fd; // relative tolerance for finite difference comparisons +}; + +map, shared_ptr> TestConsistency::cache = {}; + +// --------------- Definitions for individual consistency tests --------------- + +TEST_P(TestConsistency, h_eq_u_plus_Pv) { + double h, u, v; + try { + h = phase->enthalpy_mole(); + u = phase->intEnergy_mole(); + v = phase->molarVolume(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + EXPECT_NEAR(h, u + p * v, atol); +} + +TEST_P(TestConsistency, g_eq_h_minus_Ts) { + double g, h, s; + try { + g = phase->gibbs_mole(); + h = phase->enthalpy_mole(); + s = phase->entropy_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + if (phase->type() == "plasma") { + // Whenever a temperature can be defined, the following relation holds: + // g_k = h_k - T_k * s_k + // When multiplying by mole fraction and summing over all species, we get: + // sum_k X_k * g_k = sum_k X_k * h_k - sum_k X_k * T_k * s_k + // The left side is simply g. The first term on the right side is h. + // The second term on the right side can be separated into electron and + // non-electron contributions (for a 2 temperature plasma model): + // sum_k X_k * T_k * s_k = T * sum_{k!=e} X_k * s_k + Te * X_e * s_e + // = T * (s - X_e * s_e) + Te * X_e * s_e + // = T * s + X_e * s_e * (Te - T) + // Rearranging gives: + // g = h - T * s - X_e * s_e * (Te - T) + + // Get partial molar entropy of electron species. + vector sk(nsp); + phase->getPartialMolarEntropies(sk); + double se = sk[ke]; + // Get mole fraction of electron species. + double Xe = phase->moleFraction(ke); + + EXPECT_NEAR(g, h - T * s - Xe * se * (Te - T), atol); + } + else { + EXPECT_NEAR(g, h - T * s, atol); + } +} + +TEST_P(TestConsistency, hk_eq_uk_plus_P_vk) +{ + vector hk(nsp), uk(nsp), vk(nsp); + try { + phase->getPartialMolarEnthalpies(hk); + phase->getPartialMolarIntEnergies(uk); + phase->getPartialMolarVolumes(vk); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + EXPECT_NEAR(hk[k], uk[k] + p * vk[k], atol) << "k = " << k; + } +} + +TEST_P(TestConsistency, utilde_eq_uk_minus_piT_vk) +{ + vector utilde(nsp), uk(nsp), vk(nsp); + double piT; + try { + phase->getPartialMolarIntEnergies_TV(utilde); + phase->getPartialMolarIntEnergies(uk); + phase->getPartialMolarVolumes(vk); + piT = phase->internalPressure(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + if (!std::isfinite(piT)) { + GTEST_SKIP() << "internalPressure is not finite at this state"; + } + for (size_t k = 0; k < nsp; k++) { + if (k != ke) { + double expected = uk[k] - piT * vk[k]; + double tol = max({atol, rtol_fd * abs(utilde[k]), rtol_fd * abs(expected)}); + EXPECT_NEAR(utilde[k], expected, tol) << "k = " << k; + } + } +} + +TEST_P(TestConsistency, piT_eq_dudv_const_T) +{ + if (!phase->isCompressible()) { + GTEST_SKIP() << "Undefined for incompressible phase"; + } + + double piT; + double T0 = phase->temperature(); + double v0 = 1.0 / phase->density(); // specific volume [m^3/kg] + // Use a moderate perturbation to avoid cancellation in u(v) differences. + double dv = 1e-4 * v0; + double v1 = v0 - dv; + double v2 = v0 + dv; + + try { + piT = phase->internalPressure(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + + phase->setState_TD(T0, 1.0 / v1); + double u1 = phase->intEnergy_mass(); + phase->setState_TD(T0, 1.0 / v2); + double u2 = phase->intEnergy_mass(); + + double piT_fd = (u2 - u1) / (v2 - v1); + double tol = max({rtol_fd * abs(piT), rtol_fd * abs(piT_fd), atol}); + EXPECT_NEAR(piT_fd, piT, tol); +} + +TEST_P(TestConsistency, gk_eq_hk_minus_T_sk) +{ + vector gk(nsp), hk(nsp), sk(nsp); + try { + phase->getChemPotentials(gk); + phase->getPartialMolarEnthalpies(hk); + phase->getPartialMolarEntropies(sk); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + if (k != ke) { + EXPECT_NEAR(gk[k], hk[k] - T * sk[k], atol) << "k = " << k; + } + else { + EXPECT_NEAR(gk[k], hk[k] - Te * sk[k], atol) << "k = " << k; + } + } +} + +TEST_P(TestConsistency, h_eq_sum_hk_Xk) +{ + vector hk(nsp); + try { + phase->getPartialMolarEnthalpies(hk); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + EXPECT_NEAR(phase->enthalpy_mole(), phase->mean_X(hk), atol); +} + +TEST_P(TestConsistency, u_eq_sum_uk_Xk) +{ + vector uk(nsp); + double u; + try { + phase->getPartialMolarIntEnergies(uk); + u = phase->intEnergy_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + EXPECT_NEAR(u, phase->mean_X(uk), atol); +} + +TEST_P(TestConsistency, g_eq_sum_gk_Xk) +{ + vector gk(nsp); + double g; + try { + phase->getChemPotentials(gk); + g = phase->gibbs_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + EXPECT_NEAR(g, phase->mean_X(gk), atol); +} + +TEST_P(TestConsistency, s_eq_sum_sk_Xk) +{ + vector sk(nsp); + double s; + try { + phase->getPartialMolarEntropies(sk); + s = phase->entropy_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + EXPECT_NEAR(s, phase->mean_X(sk), atol); +} + +TEST_P(TestConsistency, v_eq_sum_vk_Xk) +{ + vector vk(nsp); + try { + phase->getPartialMolarVolumes(vk); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + EXPECT_NEAR(phase->molarVolume(), phase->mean_X(vk), atol_v); +} + +TEST_P(TestConsistency, cp_eq_sum_cpk_Xk) +{ + vector cpk(nsp); + double cp; + try { + phase->getPartialMolarCp(cpk); + cp = phase->cp_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + EXPECT_NEAR(cp, phase->mean_X(cpk), atol); +} + +TEST_P(TestConsistency, cp_eq_dhdT) +{ + double h1, cp1; + try { + h1 = phase->enthalpy_mole(); + cp1 = phase->cp_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + double T1 = phase->temperature(); + double dT = 1e-5 * phase->temperature(); + phase->setState_TP(T1 + dT, phase->pressure()); + double h2 = phase->enthalpy_mole(); + double cp2 = phase->cp_mole(); + double cp_mid = 0.5 * (cp1 + cp2); + double cp_fd = (h2 - h1) / dT; + EXPECT_NEAR(cp_fd, cp_mid, max({rtol_fd * cp_mid, rtol_fd * cp_fd, atol})); +} + +TEST_P(TestConsistency, cv_eq_dudT) +{ + double u1, cv1; + try { + u1 = phase->intEnergy_mole(); + cv1 = phase->cv_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + double T1 = phase->temperature(); + double dT = 1e-5 * phase->temperature(); + if (phase->isCompressible()) { + phase->setState_TD(T1 + dT, phase->density()); + } else { + phase->setTemperature(T1 + dT); + } + double u2 = phase->intEnergy_mole(); + double cv2 = phase->cv_mole(); + double cv_mid = 0.5 * (cv1 + cv2); + double cv_fd = (u2 - u1) / dT; + EXPECT_NEAR(cv_fd, cv_mid, max({rtol_fd * cv_mid, rtol_fd * cv_fd, atol})); +} + +TEST_P(TestConsistency, cp_eq_dsdT_const_p_times_T) +{ + double s1, cp1; + try { + s1 = phase->entropy_mole(); + cp1 = phase->cp_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + double T1 = phase->temperature(); + double dT = 1e-4 * phase->temperature(); + phase->setState_TP(T1 + dT, phase->pressure()); + double s2 = phase->entropy_mole(); + double cp2 = phase->cp_mole(); + double cp_mid = 0.5 * (cp1 + cp2); + double cp_fd = (T1 + dT/2) * (s2 - s1) / dT; + EXPECT_NEAR(cp_fd, cp_mid, max({rtol_fd * cp_mid, rtol_fd * cp_fd, atol})); +} + +TEST_P(TestConsistency, cv_eq_dsdT_const_v_times_T) +{ + double s1, cv1; + try { + s1 = phase->entropy_mole(); + cv1 = phase->cv_mole(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + double T1 = phase->temperature(); + double dT = 1e-4 * phase->temperature(); + if (phase->isCompressible()) { + phase->setState_TD(T1 + dT, phase->density()); + } else { + phase->setTemperature(T1 + dT); + } + double s2 = phase->entropy_mole(); + double cv2 = phase->cv_mole(); + double cv_mid = 0.5 * (cv1 + cv2); + double cv_fd = (T1 + dT/2) * (s2 - s1) / dT; + EXPECT_NEAR(cv_fd, cv_mid, max({rtol_fd * cv_mid, rtol_fd * cv_fd, atol})); +} + +TEST_P(TestConsistency, dsdP_const_T_eq_minus_dV_dT_const_P) +{ + double P0 = phase->pressure(); + double T0 = phase->temperature(); + double s1, v1; + double dP = 1e-4 * P0; + double dT = 1e-4 * T0; + try { + phase->setState_TP(T0, P0 - dP); + s1 = phase->entropy_mole(); + phase->setState_TP(T0 - dT, P0); + v1 = phase->molarVolume(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + phase->setState_TP(T0, P0 + dP); + double s2 = phase->entropy_mole(); + double dsdP = (s2 - s1) / (2 * dP); + + phase->setState_TP(T0 + dT, P0); + double v2 = phase->molarVolume(); + double dvdT = (v2 - v1) / (2 * dT); + double tol = rtol_fd * std::max(std::abs(dsdP), std::abs(dvdT)); + EXPECT_NEAR(dsdP, -dvdT, tol); +} + +TEST_P(TestConsistency, dSdv_const_T_eq_dPdT_const_V) { + if (phase->isCompressible()) { + double s1; + try { + s1 = phase->entropy_mass(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + double v1 = 1 / phase->density(); + double P1 = phase->pressure(); + double v2 = v1 * (1 + 1e-7); + phase->setState_TD(T, 1 / v2); + double s2 = phase->entropy_mass(); + double dsdv = (s2 - s1) / (v2 - v1); + + double T2 = T * (1 + 1e-7); + phase->setState_TD(T2, 1 / v1); + double P2 = phase->pressure(); + double dPdT = (P2 - P1) / (T2 - T); + double tol = rtol_fd * std::max(std::abs(dPdT), std::abs(dsdv)); + EXPECT_NEAR(dsdv, dPdT, tol); + } else { + GTEST_SKIP() << "Undefined for incompressible phase"; + } +} + +TEST_P(TestConsistency, betaT_eq_minus_dmv_dP_const_T_div_mv) +{ + double betaT1; + try { + betaT1 = phase->isothermalCompressibility(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + + double T = phase->temperature(); + double P1 = phase->pressure(); + double mv1 = phase->molarVolume(); + + double P2 = P1 * (1 + 1e-6); + phase->setState_TP(T, P2); + double betaT2 = phase->isothermalCompressibility(); + double mv2 = phase->molarVolume(); + + double betaT_mid = 0.5 * (betaT1 + betaT2); + double mv_mid = 0.5 * (mv1 + mv2); + double betaT_fd = -1 / mv_mid * (mv2 - mv1) / (P2 - P1); + + EXPECT_NEAR(betaT_fd, betaT_mid, + max({rtol_fd * betaT_mid, rtol_fd * betaT_fd, atol_c})); +} + +TEST_P(TestConsistency, alphaV_eq_dmv_dT_const_P_div_mv) +{ + double alphaV1; + try { + alphaV1 = phase->thermalExpansionCoeff(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + + double P = phase->pressure(); + double T1 = phase->temperature(); + double mv1 = phase->molarVolume(); + + double T2 = T1 * (1 + 1e-6); + phase->setState_TP(T2, P); + double alphaV2 = phase->thermalExpansionCoeff(); + double mv2 = phase->molarVolume(); + + double alphaV_mid = 0.5 * (alphaV1 + alphaV2); + double mv_mid = 0.5 * (mv1 + mv2); + double alphaV_fd = 1 / mv_mid * (mv2 - mv1) / (T2 - T1); + + EXPECT_NEAR(alphaV_fd, alphaV_mid, + max({rtol_fd * alphaV_mid, rtol_fd * alphaV_fd, atol_e})); +} + +TEST_P(TestConsistency, c_eq_sqrt_dP_drho_const_s) +{ + double c1; + try { + c1 = phase->soundSpeed(); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + + double P1 = phase->pressure(); + double rho1 = phase->density(); + + double rho2 = rho1 * (1 + 1e-6); + phase->setState_SV(phase->entropy_mass(), 1 / rho2); + double c2 = phase->soundSpeed(); + double P2 = phase->pressure(); + + double c_mid = 0.5 * (c1 + c2); + double c_fd = sqrt((P2 - P1) / (rho2 - rho1)); + + EXPECT_NEAR(c_fd, c_mid, max({rtol_fd * c_mid, rtol_fd * c_fd, atol})); +} + +// ---------- Tests for consistency of standard state properties --------------- + +TEST_P(TestConsistency, hk0_eq_uk0_plus_p_vk0) +{ + vector h0(nsp), u0(nsp), v0(nsp); + try { + phase->getEnthalpy_RT(h0); + phase->getIntEnergy_RT(u0); + phase->getStandardVolumes(v0); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + if (k != ke) { + EXPECT_NEAR(h0[k] * RT, u0[k] * RT + p * v0[k], atol) << "k = " << k; + } else { + EXPECT_NEAR(h0[k] * RTe, u0[k] * RTe + p * v0[k], atol) << "k = " << k; + } + } +} + +TEST_P(TestConsistency, gk0_eq_hk0_minus_T_sk0) +{ + vector g0(nsp), h0(nsp), s0(nsp); + try { + phase->getEnthalpy_RT(h0); + phase->getGibbs_RT(g0); + phase->getEntropy_R(s0); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + EXPECT_NEAR(g0[k] * RT , + h0[k] * RT - T * s0[k] * GasConstant, atol) << "k = " << k; + } +} + +TEST_P(TestConsistency, cpk0_eq_dhk0dT) +{ + vector h1(nsp), h2(nsp), cp1(nsp), cp2(nsp); + try { + phase->getEnthalpy_RT(h1); + phase->getCp_R(cp1); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + // Perturb temperature. + double T1 = phase->temperature(); + double dT = 1e-5 * phase->temperature(); + // For plasma phases, also perturb the electron temperature. + double Te = phase->electronTemperature(); + double dTe = 1e-5 * Te; + + phase->setState_TP(T1 + dT, phase->pressure()); + + if (phase->type() == "plasma") { + phase->setElectronTemperature(Te + dTe); + } + phase->getEnthalpy_RT(h2); + phase->getCp_R(cp2); + + for (size_t k = 0; k < nsp; k++) { + // Determine effective temperature and perturbation for species k. + double T_eff = (k == ke) ? Te : T1; + double dT_eff = (k == ke) ? dTe : dT; + + double cp_mid = 0.5 * (cp1[k] + cp2[k]) * GasConstant; + double cp_fd = (h2[k] * (T_eff + dT_eff) - h1[k] * T_eff) / dT_eff * GasConstant; + double tol = max({rtol_fd * std::abs(cp_mid), rtol_fd * std::abs(cp_fd), atol}); + EXPECT_NEAR(cp_fd, cp_mid, tol) << "k = " << k; + } +} + +TEST_P(TestConsistency, standard_gibbs_nondim) +{ + vector g0_RT(nsp), mu0(nsp); + try { + phase->getGibbs_RT(g0_RT); + phase->getStandardChemPotentials(mu0); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + if (k != ke) { + EXPECT_NEAR(g0_RT[k] * RT , mu0[k], atol) << "k = " << k; + } else { + EXPECT_NEAR(g0_RT[k] * RTe , mu0[k], atol) << "k = " << k; + } + } +} + +TEST_P(TestConsistency, chem_potentials_to_activities) { + vector mu0(nsp), mu(nsp), a(nsp); + try { + phase->getChemPotentials(mu); + phase->getStandardChemPotentials(mu0); + phase->getActivities(a); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + if (k != ke) { + double a_from_mu = exp((mu[k] - mu0[k]) / RT); + double scale = std::max(std::abs(a[k]), std::abs(a_from_mu)); + EXPECT_NEAR(a_from_mu, a[k], 1e-9 * scale + 1e-14) << "k = " << k; + } else { + double a_from_mu = exp((mu[k] - mu0[k]) / RTe); + double scale = std::max(std::abs(a[k]), std::abs(a_from_mu)); + EXPECT_NEAR(a_from_mu, a[k], 1e-9 * scale + 1e-14) << "k = " << k; + } + } +} + +TEST_P(TestConsistency, activity_coeffs) { + vector a(nsp), gamma(nsp), X(nsp); + try { + phase->getActivities(a); + phase->getActivityCoefficients(gamma); + phase->getMoleFractions(X); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + double scale = std::max(std::abs(a[k]), std::abs(gamma[k] * X[k])); + EXPECT_NEAR(a[k], gamma[k] * X[k], 1e-10 * scale + 1e-20) << "k = " << k; + } +} + +TEST_P(TestConsistency, activity_concentrations) { + vector a(nsp), Cact(nsp); + try { + phase->getActivities(a); + phase->getActivityConcentrations(Cact); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + double C0k = phase->standardConcentration(k); + EXPECT_NEAR(a[k], Cact[k] / C0k, 1e-9 * std::abs(a[k])) << "k = " << k; + } +} + +TEST_P(TestConsistency, log_standard_concentrations) { + for (size_t k = 0; k < nsp; k++) { + double c0 = phase->standardConcentration(k); + EXPECT_NEAR(c0, exp(phase->logStandardConc(k)), 1e-10 * c0) << "k = " << k; + } +} + +TEST_P(TestConsistency, log_activity_coeffs) { + vector gamma(nsp), log_gamma(nsp); + try { + phase->getActivityCoefficients(gamma); + phase->getLnActivityCoefficients(log_gamma); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + EXPECT_NEAR(gamma[k], exp(log_gamma[k]), 1e-10 * gamma[k]) << "k = " << k; + } +} + +// -------------------- Tests for reference state properties ------------------- + +TEST_P(TestConsistency, hRef_eq_uRef_plus_P_vRef) +{ + vector hRef(nsp), uRef(nsp), vRef(nsp); + try { + phase->getEnthalpy_RT_ref(hRef); + phase->getIntEnergy_RT_ref(uRef); + phase->getStandardVolumes_ref(vRef); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + if (k != ke) { + EXPECT_NEAR(hRef[k] * RT, uRef[k] * RT + OneAtm * vRef[k], atol) << "k = " << k; + } else { + EXPECT_NEAR(hRef[k] * RTe, uRef[k] * RTe + OneAtm * vRef[k], atol) << "k = " << k; + } + } +} + +TEST_P(TestConsistency, gRef_eq_hRef_minus_T_sRef) +{ + vector hRef(nsp), gRef_RT(nsp), gRef(nsp), sRef(nsp); + try { + phase->getEnthalpy_RT_ref(hRef); + phase->getGibbs_ref(gRef); + phase->getGibbs_RT_ref(gRef_RT); + phase->getEntropy_R_ref(sRef); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + for (size_t k = 0; k < nsp; k++) { + if (k != ke) { + EXPECT_NEAR(gRef[k], gRef_RT[k] * RT, atol) << "k = " << k; + EXPECT_NEAR(gRef[k], hRef[k] * RT - T * sRef[k] * GasConstant, + atol) << "k = " << k; + } else { + EXPECT_NEAR(gRef[k], gRef_RT[k] * RTe, atol) << "k = " << k; + EXPECT_NEAR(gRef[k], hRef[k] * RTe - Te * sRef[k] * GasConstant, + atol) << "k = " << k; + } + } +} + +TEST_P(TestConsistency, cpRef_eq_dhRefdT) +{ + vector h1(nsp), h2(nsp), cp1(nsp), cp2(nsp); + try { + phase->getEnthalpy_RT_ref(h1); + phase->getCp_R_ref(cp1); + } catch (NotImplementedError& err) { + GTEST_SKIP() << err.getMethod() << " threw NotImplementedError"; + } + // Perturb temperature. + double T1 = phase->temperature(); + double dT = 1e-5 * phase->temperature(); + // For plasma phases, also perturb the electron temperature. + double Te = phase->electronTemperature(); + double dTe = 1e-5 * Te; + + phase->setState_TP(T1 + dT, phase->pressure()); + + if (phase->type() == "plasma") { + phase->setElectronTemperature(Te + dTe); + } + phase->getEnthalpy_RT_ref(h2); + phase->getCp_R_ref(cp2); + + for (size_t k = 0; k < nsp; k++) { + // Determine effective temperature and perturbation for species k. + double T_eff = (k == ke) ? Te : T1; + double dT_eff = (k == ke) ? dTe : dT; + + double cp_mid = 0.5 * (cp1[k] + cp2[k]) * GasConstant; + double cp_fd = (h2[k] * (T_eff + dT_eff) - h1[k] * T_eff) / dT_eff * GasConstant; + double tol = max({rtol_fd * std::abs(cp_mid), rtol_fd * std::abs(cp_fd), atol}); + EXPECT_NEAR(cp_fd, cp_mid, tol) << "k = " << k << "(ke = " << ke << ")"; + } +} + +// ------- Instantiation of test suites defined in consistency-cases.yaml ------ + +INSTANTIATE_TEST_SUITE_P(IdealGas, TestConsistency, + testing::Combine( + testing::Values(getSetup("ideal-gas-h2o2")), + testing::ValuesIn(getStates("ideal-gas-h2o2"))) +); + +INSTANTIATE_TEST_SUITE_P(RedlichKwong, TestConsistency, + testing::Combine( + testing::Values(getSetup("redlich-kwong")), + testing::ValuesIn(getStates("redlich-kwong"))) +); + +INSTANTIATE_TEST_SUITE_P(PengRobinson, TestConsistency, + testing::Combine( + testing::Values(getSetup("peng-robinson")), + testing::ValuesIn(getStates("peng-robinson"))) +); + +INSTANTIATE_TEST_SUITE_P(IdealMolalSolution, TestConsistency, + testing::Combine( + testing::Values(getSetup("ideal-molal-solution")), + testing::ValuesIn(getStates("ideal-molal-solution"))) +); + +INSTANTIATE_TEST_SUITE_P(IdealSolidSolnPhase1, TestConsistency, + testing::Combine( + testing::Values(getSetup("ideal-condensed-1")), + testing::ValuesIn(getStates("ideal-condensed-1"))) +); + +INSTANTIATE_TEST_SUITE_P(IdealSolidSolnPhase2, TestConsistency, + testing::Combine( + testing::Values(getSetup("ideal-condensed-2")), + testing::ValuesIn(getStates("ideal-condensed-2"))) +); + +INSTANTIATE_TEST_SUITE_P(BinarySolutionTabulated, TestConsistency, + testing::Combine( + testing::Values(getSetup("binary-solution-tabulated")), + testing::ValuesIn(getStates("binary-solution-tabulated"))) +); + +INSTANTIATE_TEST_SUITE_P(ElectronCloud, TestConsistency, + testing::Combine( + testing::Values(getSetup("electron-cloud")), + testing::ValuesIn(getStates("electron-cloud"))) +); + +INSTANTIATE_TEST_SUITE_P(NitrogenPureFluid, TestConsistency, + testing::Combine( + testing::Values(getSetup("nitrogen-purefluid")), + testing::ValuesIn(getStates("nitrogen-purefluid"))) +); + +INSTANTIATE_TEST_SUITE_P(PlasmaPhase, TestConsistency, + testing::Combine( + testing::Values(getSetup("plasma")), + testing::ValuesIn(getStates("plasma"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckelDilute, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-dilute")), + testing::ValuesIn(getStates("debye-huckel-dilute"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckelDilute_IAPWS, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-dilute-IAPWS")), + testing::ValuesIn(getStates("debye-huckel-dilute-IAPWS"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckel_b_dot_ak, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-B-dot-ak")), + testing::ValuesIn(getStates("debye-huckel-B-dot-ak"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckel_b_dot_ak_IAPWS, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-B-dot-ak-IAPWS")), + testing::ValuesIn(getStates("debye-huckel-B-dot-ak-IAPWS"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckel_b_dot_a, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-B-dot-a")), + testing::ValuesIn(getStates("debye-huckel-B-dot-a"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckel_b_dot_a_IAPWS, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-B-dot-a-IAPWS")), + testing::ValuesIn(getStates("debye-huckel-B-dot-a-IAPWS"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckel_pitzer_beta_ij, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-pitzer-beta_ij")), + testing::ValuesIn(getStates("debye-huckel-pitzer-beta_ij"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckel_pitzer_beta_ij_IAPWS, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-pitzer-beta_ij-IAPWS")), + testing::ValuesIn(getStates("debye-huckel-pitzer-beta_ij-IAPWS"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckel_beta_ij, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-beta_ij")), + testing::ValuesIn(getStates("debye-huckel-beta_ij"))) +); + +INSTANTIATE_TEST_SUITE_P(DebyeHuckel_beta_ij_IAPWS, TestConsistency, + testing::Combine( + testing::Values(getSetup("debye-huckel-beta_ij-IAPWS")), + testing::ValuesIn(getStates("debye-huckel-beta_ij-IAPWS"))) +); + +INSTANTIATE_TEST_SUITE_P(Margules, TestConsistency, + testing::Combine( + testing::Values(getSetup("margules")), + testing::ValuesIn(getStates("margules"))) +); + +INSTANTIATE_TEST_SUITE_P(FixedStoichiometry, TestConsistency, + testing::Combine( + testing::Values(getSetup("fixed-stoichiometry")), + testing::ValuesIn(getStates("fixed-stoichiometry"))) +); + +INSTANTIATE_TEST_SUITE_P(IdealSurface, TestConsistency, + testing::Combine( + testing::Values(getSetup("ideal-surface")), + testing::ValuesIn(getStates("ideal-surface"))) +); + +INSTANTIATE_TEST_SUITE_P(IdealEdge, TestConsistency, + testing::Combine( + testing::Values(getSetup("ideal-edge")), + testing::ValuesIn(getStates("ideal-edge"))) +); + +INSTANTIATE_TEST_SUITE_P(CoverageDependentSurface, TestConsistency, + testing::Combine( + testing::Values(getSetup("coverage-dependent-surface")), + testing::ValuesIn(getStates("coverage-dependent-surface"))) +); + +INSTANTIATE_TEST_SUITE_P(LiquidWaterIapws95, TestConsistency, + testing::Combine( + testing::Values(getSetup("liquid-water-IAPWS95")), + testing::ValuesIn(getStates("liquid-water-IAPWS95"))) +); + +INSTANTIATE_TEST_SUITE_P(IdealSolnVPSS_simple, TestConsistency, + testing::Combine( + testing::Values(getSetup("ideal-solution-VPSS-simple")), + testing::ValuesIn(getStates("ideal-solution-VPSS-simple"))) +); + +INSTANTIATE_TEST_SUITE_P(IdealSolnVPSS_HKFT, TestConsistency, + testing::Combine( + testing::Values(getSetup("ideal-solution-VPSS-HKFT")), + testing::ValuesIn(getStates("ideal-solution-VPSS-HKFT"))) +); + +INSTANTIATE_TEST_SUITE_P(RedlichKister_LiC6, TestConsistency, + testing::Combine( + testing::Values(getSetup("Redlich-Kister-LiC6")), + testing::ValuesIn(getStates("Redlich-Kister-LiC6"))) +); + +INSTANTIATE_TEST_SUITE_P(RedlichKister_complex, TestConsistency, + testing::Combine( + testing::Values(getSetup("Redlich-Kister-complex")), + testing::ValuesIn(getStates("Redlich-Kister-complex"))) +); + +INSTANTIATE_TEST_SUITE_P(HMWSoln, TestConsistency, + testing::Combine( + testing::Values(getSetup("HMW-electrolyte")), + testing::ValuesIn(getStates("HMW-electrolyte"))) +); + +int main(int argc, char** argv) +{ + printf("Running main() from consistency.cpp\n"); + testing::InitGoogleTest(&argc, argv); + Cantera::make_deprecation_warnings_fatal(); + Cantera::CanteraError::setStackTraceDepth(0); + int result = RUN_ALL_TESTS(); + Cantera::appdelete(); + return result; +}