Update handling of critical properties

doc/sphinx/yaml/species.rst

@@ -26,6 +26,12 @@ The fields of a ``species`` entry are:
     :ref:`sec-yaml-species-eos`. If this field is absent and a model is
     required, the ``ideal-gas`` model is assumed.
 
+``critical-parameters``
+    Mapping containing parameters related to the critical state of a species. Used in
+    models that incorporate "real gas" effects, such as
+    :ref:`sec-yaml-eos-redlich-kwong`.
+    See :ref:`sec-yaml-species-crit-props`.
+
 ``transport``
     Mapping containing the species transport model specification and
     parameters. See :ref:`sec-yaml-species-transport`.
@@ -257,6 +263,19 @@ Species equation of state models
     - :ref:`molar-volume-temperature-polynomial <sec-yaml-eos-molar-volume-temperature-polynomial>`
     - :ref:`Redlich-Kwong <sec-yaml-eos-redlich-kwong>`
 
+.. _sec-yaml-species-crit-props:
+
+Species critical state parameters
+=================================
+
+``critical-temperature``
+    The critical temperature of the species.
+
+``critical-pressure``
+    The critical pressure of the species.
+
+``acentric-factor``
+    Pitzer's acentric factor :math:`omega`.
 
 .. _sec-yaml-eos-constant-volume:
 
@@ -460,7 +479,9 @@ The additional fields of a ``gas`` transport entry are:
     The rotational relaxation collision number at 298 K [-]. Default 0.0.
 
 ``acentric-factor``
-    Pitzer's acentric factor [-]. Default 0.0.
+    Pitzer's acentric factor [-]. Default 0.0. This value may also be specified as part
+    of the :ref:`critical-parameters <sec-yaml-species-crit-props>` field, in which case
+    the value provided there supersedes this one.
 
 ``dispersion-coefficient``
     The dispersion coefficient, normalized by :math:`e^2` [Å^5]. Default 0.0.

include/cantera/base/AnyMap.h

@@ -455,6 +455,9 @@ class AnyMap : public AnyBase
     //! Mainly for internal use in reading or writing from files.
     void setMetadata(const std::string& key, const AnyValue& value);
 
+    //! Copy metadata including input line/column from an existing AnyMap
+    void copyMetadata(const AnyMap& other);
+
     //! Propagate metadata to any child elements
     void propagateMetadata(shared_ptr<AnyMap>& file);
 

include/cantera/thermo/PengRobinson.h

@@ -151,20 +151,6 @@ class PengRobinson : public MixtureFugacityTP
     virtual bool addSpecies(shared_ptr<Species> spec);
     virtual void initThermo();
 
-    //! Retrieve \f$a\f$ and \f$b\f$ coefficients by looking up tabulated critical parameters
-    /*!
-    *  If `pureFluidParameters` are not provided for any species in the phase,
-    *  consult the critical properties tabulated in `data/inputs/critProperties.xml`.
-    *  If the species is found there, calculate pure fluid parameters \f$a_k\f$ and \f$b_k\f$ as:
-    *  \f[ a_k = 0.4278 R^2 T_c^2 / P_c \f]
-    *
-    *  and:
-    *  \f[ b_k = 0.08664 R T_c/ P_c \f]
-    *
-    *  @param iName    Name of the species
-    */
-    virtual vector_fp getCoeff(const std::string& iName);
-
     //! Set the pure fluid interaction parameters for a species
     /*!
      *  @param species   Name of the species

include/cantera/thermo/RedlichKwongMFTP.h

@@ -141,13 +141,16 @@ class RedlichKwongMFTP : public MixtureFugacityTP
     //! Retrieve a and b coefficients by looking up tabulated critical parameters
     /*!
      *  If pureFluidParameters are not provided for any species in the phase,
-     *  consult the critical properties tabulated in /thermo/critProperties.xml.
+     *  consult the critical properties tabulated in `critical-properties.yaml`.
      *  If the species is found there, calculate pure fluid parameters a_k and b_k as:
      *  \f[ a_k = 0.4278*R**2*T_c^2.5/P_c \f]
      *
      *  and:
      *  \f[ b_k = 0.08664*R*T_c/P_c \f]
      *
+     * @deprecated To be removed after Cantera 2.6. Use of critical-properties.yaml is
+     *     integrated into initThermo() for YAML input files.
+     *
      *  @param iName    Name of the species
      */
     virtual std::vector<double> getCoeff(const std::string& iName);

src/base/AnyMap.cpp

@@ -1490,6 +1490,28 @@ void AnyMap::setMetadata(const std::string& key, const AnyValue& value)
     propagateMetadata(m_metadata);
 }
 
+void AnyMap::copyMetadata(const AnyMap& other)
+{
+    m_line = other.m_line;
+    m_column = other.m_column;
+    if (!other.m_metadata) {
+        return;
+    }
+
+    if (m_metadata) {
+        // Fork the metadata tree at this point to avoid affecting parent nodes
+        m_metadata = make_shared<AnyMap>(*m_metadata);
+    } else {
+        m_metadata = make_shared<AnyMap>();
+    }
+
+    for (const auto& item : *other.m_metadata) {
+        (*m_metadata)[item.first] = item.second;
+    }
+
+    propagateMetadata(m_metadata);
+}
+
 bool AnyMap::getBool(const std::string& key, bool default_) const
 {
     return (hasKey(key)) ? m_data.at(key).asBool() : default_;

src/thermo/PengRobinson.cpp

@@ -345,124 +345,86 @@ bool PengRobinson::addSpecies(shared_ptr<Species> spec)
     return added;
 }
 
-vector<double> PengRobinson::getCoeff(const std::string& iName)
-{
-    vector_fp spCoeff{ NAN, NAN, NAN };
-
-    // Get number of species in the database
-    // open xml file critProperties.xml
-    XML_Node* doc = get_XML_File("critProperties.xml");
-    size_t nDatabase = doc->nChildren();
-
-    // Loop through all species in the database and attempt to match supplied
-    // species to each. If present, calculate pureFluidParameters a_k and b_k
-    // based on crit properties T_c and P_c:
-    for (size_t isp = 0; isp < nDatabase; isp++) {
-        XML_Node& acNodeDoc = doc->child(isp);
-        std::string iNameLower = toLowerCopy(iName);
-        std::string dbName = toLowerCopy(acNodeDoc.attrib("name"));
-
-        // Attempt to match provided species iName to current database species
-        //  dbName:
-        if (iNameLower == dbName) {
-            // Read from database and calculate a and b coefficients
-            double vParams;
-            double T_crit = 0.0, P_crit = 0.0, w_ac = 0.0;
-
-            if (acNodeDoc.hasChild("Tc")) {
-                vParams = 0.0;
-                XML_Node& xmlChildCoeff = acNodeDoc.child("Tc");
-                if (xmlChildCoeff.hasAttrib("value")) {
-                    std::string critTemp = xmlChildCoeff.attrib("value");
-                    vParams = strSItoDbl(critTemp);
-                }
-                if (vParams <= 0.0) { //Assuming that Pc and Tc are non zero.
-                    throw CanteraError("PengRobinson::getCoeff",
-                        "Critical Temperature must be positive");
-                }
-                T_crit = vParams;
-            }
-            if (acNodeDoc.hasChild("Pc")) {
-                vParams = 0.0;
-                XML_Node& xmlChildCoeff = acNodeDoc.child("Pc");
-                if (xmlChildCoeff.hasAttrib("value")) {
-                    std::string critPressure = xmlChildCoeff.attrib("value");
-                    vParams = strSItoDbl(critPressure);
-                }
-                if (vParams <= 0.0) { //Assuming that Pc and Tc are non zero.
-                    throw CanteraError("PengRobinson::getCoeff",
-                        "Critical Pressure must be positive");
-                }
-                P_crit = vParams;
-            }
-            if (acNodeDoc.hasChild("omega")) {
-                vParams = 0.0;
-                XML_Node& xmlChildCoeff = acNodeDoc.child("omega");
-                if (xmlChildCoeff.hasChild("value")) {
-                    std::string acentric_factor = xmlChildCoeff.attrib("value");
-                    vParams = strSItoDbl(acentric_factor);
-                }
-                w_ac = vParams;
-            }
-
-            spCoeff[0] = omega_a * (GasConstant * GasConstant) * (T_crit * T_crit) / P_crit; //coeff a
-            spCoeff[1] = omega_b * GasConstant * T_crit / P_crit; // coeff b
-            spCoeff[2] = w_ac; // acentric factor
-            break;
-        }
-    }
-    // If the species is not present in the database, throw an error
-    if(isnan(spCoeff[0]))
-    {
-        throw CanteraError("PengRobinson::getCoeff",
-            "Species '{}' is not present in the database", iName);
-    }
-    return spCoeff;
-}
-
 void PengRobinson::initThermo()
 {
+    // Contents of 'critical-properties.yaml', loaded later if needed
+    AnyMap critPropsDb;
+    std::unordered_map<std::string, AnyMap*> dbSpecies;
+
     for (auto& item : m_species) {
-        // Read a and b coefficients and acentric factor w_ac from species input
-        // information, specified in a YAML input file.
-        if (item.second->input.hasKey("equation-of-state")) {
-            auto eos = item.second->input["equation-of-state"].getMapWhere(
+        auto& data = item.second->input;
+        size_t k = speciesIndex(item.first);
+        if (m_a_coeffs(k, k) != 0.0) {
+            continue;
+        }
+        bool foundCoeffs = false;
+
+        if (data.hasKey("equation-of-state") &&
+            data["equation-of-state"].hasMapWhere("model", "Peng-Robinson"))
+        {
+            // Read a and b coefficients and acentric factor w_ac from species input
+            // information, specified in a YAML input file.
+            auto eos = data["equation-of-state"].getMapWhere(
                 "model", "Peng-Robinson");
-            double a0 = 0;
-            if (eos["a"].isScalar()) {
-                a0 = eos.convert("a", "Pa*m^6/kmol^2");
+            if (eos.hasKey("a") && eos.hasKey("b") && eos.hasKey("acentric-factor")) {
+                double a0 = eos.convert("a", "Pa*m^6/kmol^2");
+                double b = eos.convert("b", "m^3/kmol");
+                // unitless acentric factor:
+                double w = eos["acentric-factor"].asDouble();
+                setSpeciesCoeffs(item.first, a0, b, w);
+                foundCoeffs = true;
             }
-            double b = eos.convert("b", "m^3/kmol");
-            // unitless acentric factor:
-            double w = eos["acentric-factor"].asDouble();
 
-            setSpeciesCoeffs(item.first, a0, b, w);
             if (eos.hasKey("binary-a")) {
                 AnyMap& binary_a = eos["binary-a"].as<AnyMap>();
                 const UnitSystem& units = binary_a.units();
                 for (auto& item2 : binary_a) {
-                    double a0 = 0;
-                    if (item2.second.isScalar()) {
-                        a0 = units.convert(item2.second, "Pa*m^6/kmol^2");
-                    }
+                    double a0 = units.convert(item2.second, "Pa*m^6/kmol^2");
                     setBinaryCoeffs(item.first, item2.first, a0);
                 }
             }
+            if (foundCoeffs) {
+                continue;
+            }
+        }
+
+        // Coefficients have not been populated from model-specific input
+        double Tc = NAN, Pc = NAN, omega_ac = NAN;
+        if (data.hasKey("critical-parameters")) {
+            // Use critical state information stored in the species entry to
+            // calculate a, b, and the acentric factor.
+            auto& critProps = data["critical-parameters"].as<AnyMap>();
+            Tc = critProps.convert("critical-temperature", "K");
+            Pc = critProps.convert("critical-pressure", "Pa");
+            omega_ac = critProps["acentric-factor"].asDouble();
         } else {
-            // Check if a and b are already populated for this species (only the
-            // diagonal elements of a). If not, then search 'critProperties.xml'
-            // to find critical temperature and pressure to calculate a and b.
-            size_t k = speciesIndex(item.first);
-            if (m_a_coeffs(k, k) == 0.0) {
-                vector<double> coeffs = getCoeff(item.first);
-
-                // Check if species was found in the database of critical
-                // properties, and assign the results
-                if (!isnan(coeffs[0])) {
-                    setSpeciesCoeffs(item.first, coeffs[0], coeffs[1], coeffs[2]);
-                }
+            // Search 'crit-properties.yaml' to find Tc and Pc. Load data if needed.
+            if (critPropsDb.empty()) {
+                critPropsDb = AnyMap::fromYamlFile("critical-properties.yaml");
+                dbSpecies = critPropsDb["species"].asMap("name");
+            }
+
+            // All names in critical-properties.yaml are upper case
+            auto ucName = boost::algorithm::to_upper_copy(item.first);
+            if (dbSpecies.count(ucName)) {
+                auto& spec = *dbSpecies.at(ucName);
+                auto& critProps = spec["critical-parameters"].as<AnyMap>();
+                Tc = critProps.convert("critical-temperature", "K");
+                Pc = critProps.convert("critical-pressure", "Pa");
+                omega_ac = critProps["acentric-factor"].asDouble();
             }
         }
+
+        // Check if critical properties were found in either location
+        if (!isnan(Tc)) {
+            double a = omega_a * std::pow(GasConstant * Tc, 2) / Pc;
+            double b = omega_b * GasConstant * Tc / Pc;
+            setSpeciesCoeffs(item.first, a, b, omega_ac);
+        } else {
+            throw InputFileError("PengRobinson::initThermo", data,
+            "No Peng-Robinson model parameters or critical properties found for "
+            "species '{}'", item.first);
+        }
     }
 }