[Input] Create IonsFromNeutral objects from YAML definitions

include/cantera/thermo/ThermoFactory.h

@@ -127,7 +127,7 @@ ThermoPhase* newPhase(XML_Node& phase);
  *     which will be used as the default location from which to read species
  *     definitions.
  */
-unique_ptr<ThermoPhase> newPhase(const AnyMap& phaseNode,
+unique_ptr<ThermoPhase> newPhase(AnyMap& phaseNode,
                                  const AnyMap& rootNode=AnyMap());
 
 //! Create and Initialize a ThermoPhase object from an XML input file.
@@ -207,7 +207,7 @@ void importPhase(XML_Node& phase, ThermoPhase* th);
  *     which will be used as the default location from which to read species
  *     definitions.
  */
-void setupPhase(ThermoPhase& phase, const AnyMap& phaseNode,
+void setupPhase(ThermoPhase& phase, AnyMap& phaseNode,
                 const AnyMap& rootNode=AnyMap());
 
 //! Add the elements given in an XML_Node tree to the specified phase

src/thermo/IonsFromNeutralVPSSTP.cpp

@@ -482,6 +482,13 @@ static double factorOverlap(const std::vector<std::string>& elnamesVN ,
 
 void IonsFromNeutralVPSSTP::initThermo()
 {
+    if (m_input.hasKey("neutral-phase") && m_input.hasKey("__file__")) {
+        string neutralName = m_input["neutral-phase"].asString();
+        AnyMap infile = AnyMap::fromYamlFile(m_input["__file__"].asString());
+        auto phaseNodes = infile["phases"].asMap("name");
+        setNeutralMoleculePhase(newPhase(*phaseNodes.at(neutralName), infile));
+    }
+
     size_t nElementsN = neutralMoleculePhase_->nElements();
     const std::vector<std::string>& elnamesVN = neutralMoleculePhase_->elementNames();
     vector_fp elemVectorN(nElementsN);
@@ -570,6 +577,9 @@ void IonsFromNeutralVPSSTP::setNeutralMoleculePhase(shared_ptr<ThermoPhase> neut
     y_.resize(numNeutralMoleculeSpecies_, 0.0);
     dlnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
     dX_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
+    for (size_t k = 0; k < nSpecies(); k++) {
+        providePDSS(k)->setParent(this, k);
+    }
 }
 
 shared_ptr<ThermoPhase> IonsFromNeutralVPSSTP::getNeutralMoleculePhase()
@@ -595,9 +605,11 @@ bool IonsFromNeutralVPSSTP::addSpecies(shared_ptr<Species> spec)
             passThroughList_.push_back(m_kk-1);
         }
 
-        if (spec->extra.hasKey("special_species")
-            && spec->extra["special_species"].asBool()) {
-            indexSpecialSpecies_ = m_kk - 1;
+        if (spec->input.hasKey("equation-of-state")) {
+            auto& ss = spec->input["equation-of-state"].as<AnyMap>();
+            if (ss.getBool("special-species", false)) {
+                indexSpecialSpecies_ = m_kk - 1;
+            }
         }
     }
     return added;

src/thermo/PDSSFactory.cpp

@@ -28,6 +28,7 @@ PDSSFactory::PDSSFactory()
     m_synonyms["water-IAPWS95"] = "water";
     reg("ions-from-neutral", []() { return new PDSS_IonsFromNeutral(); });
     m_synonyms["IonFromNeutral"] = "ions-from-neutral";
+    m_synonyms["ions-from-neutral-molecule"] = "ions-from-neutral";
     reg("temperature_polynomial", []() { return new PDSS_SSVol(); });
     m_synonyms["temperature-polynomial"] = "temperature_polynomial";
     m_synonyms["density_temperature_polynomial"] = "temperature_polynomial";

src/thermo/PDSS_IonsFromNeutral.cpp

@@ -70,6 +70,15 @@ void PDSS_IonsFromNeutral::setParametersFromXML(const XML_Node& speciesNode)
 void PDSS_IonsFromNeutral::initThermo()
 {
     PDSS::initThermo();
+    if (m_input.getBool("special-species", false)) {
+        setSpecialSpecies();
+    }
+    if (m_input.hasKey("multipliers")) {
+        for (const auto& item : m_input["multipliers"].asMap<double>()) {
+            setNeutralSpeciesMultiplier(item.first, item.second);
+        }
+    }
+
     m_p0 = neutralMoleculePhase_->refPressure();
     m_minTemp = neutralMoleculePhase_->minTemp();
     m_maxTemp = neutralMoleculePhase_->maxTemp();

src/thermo/Species.cpp

@@ -75,7 +75,7 @@ shared_ptr<Species> newSpecies(const XML_Node& species_node)
     const XML_Node* thermo = species_node.findByName("thermo");
     if (thermo && thermo->attrib("model") == "IonFromNeutral") {
         if (thermo->hasChild("specialSpecies")) {
-            s->extra["special_species"] = true;
+            s->input["equation-of-state"]["special-species"] = true;
         }
     }
 

src/thermo/ThermoFactory.cpp

@@ -70,6 +70,7 @@ ThermoFactory::ThermoFactory()
     m_synonyms["IdealGasVPSS"] = "IdealSolnVPSS";
     reg("Margules", []() { return new MargulesVPSSTP(); });
     reg("IonsFromNeutralMolecule", []() { return new IonsFromNeutralVPSSTP(); });
+    m_synonyms["ions-from-neutral-molecule"] = "IonsFromNeutralMolecule";
     reg("FixedChemPot", []() { return new FixedChemPotSSTP(); });
     m_synonyms["fixed-chemical-potential"] = "FixedChemPot";
     reg("Redlich-Kister", []() { return new RedlichKisterVPSSTP(); });
@@ -94,7 +95,7 @@ ThermoPhase* newPhase(XML_Node& xmlphase)
     return t.release();
 }
 
-unique_ptr<ThermoPhase> newPhase(const AnyMap& phaseNode, const AnyMap& rootNode)
+unique_ptr<ThermoPhase> newPhase(AnyMap& phaseNode, const AnyMap& rootNode)
 {
     unique_ptr<ThermoPhase> t(newThermoPhase(phaseNode["thermo"].asString()));
     setupPhase(*t, phaseNode, rootNode);
@@ -404,10 +405,12 @@ void addSpecies(ThermoPhase& thermo, const AnyValue& names, const AnyValue& spec
     }
 }
 
-void setupPhase(ThermoPhase& thermo, const AnyMap& phaseNode,
-                const AnyMap& rootNode)
+void setupPhase(ThermoPhase& thermo, AnyMap& phaseNode, const AnyMap& rootNode)
 {
     thermo.setName(phaseNode["name"].asString());
+    if (rootNode.hasKey("__file__")) {
+        phaseNode["__file__"] = rootNode["__file__"];
+    }
 
     // Add elements
     if (phaseNode.hasKey("elements")) {
@@ -500,7 +503,7 @@ void setupPhase(ThermoPhase& thermo, const AnyMap& phaseNode,
                 "Could not parse species declaration of type '{}'",
                 phaseNode["species"].type_str());
         }
-    } else {
+    } else if (rootNode.hasKey("species")) {
         // By default, add all species from the 'species' section
         addSpecies(thermo, AnyValue("all"), rootNode["species"]);
     }

test/data/mock_ion.xml

@@ -46,7 +46,6 @@
             KCl(L):1.2
         </neutralSpeciesMultipliers>
       </thermo>
-      <density units="g/cm3"> 0.0 </density>
     </species>
 
     <species name="Cl-">
@@ -59,7 +58,6 @@
             KCl(L):1.5
         </neutralSpeciesMultipliers>
       </thermo>
-      <density units="g/cm3"> 0.0 </density>
     </species>
 
     <species name="KCl(L)">

test/data/thermo-models.yaml

@@ -82,6 +82,16 @@ phases:
     P: 1 atm
     M: {Na+: 3.0, Cl-: 3.0, H+: 1.0499E-8, OH-: 1.3765E-6, NaCl(aq): 0.98492}
 
+- name: ions-from-neutral-molecule
+  thermo: ions-from-neutral-molecule
+  neutral-phase: KCl-neutral
+  species: [{ions-from-neutral-species: all}]
+  state: {T: 500, P: 2 bar, X: {K+: 0.1, Cl-: 0.1}}
+
+- name: KCl-neutral
+  species: [KCl(l)]
+  thermo: Margules
+
 
 species:
 - name: NaCl(s)
@@ -247,3 +257,16 @@ dh-electrolyte-species:
     molar-volume: 1.3
   electrolyte-species-type: weak-acid-associated
   weak-acid-charge: -1.0
+
+ions-from-neutral-species:
+- name: K+
+  composition: {K: 1, E: -1}
+  equation-of-state:
+    model: ions-from-neutral-molecule
+    multipliers: {KCl(l): 1.2}
+- name: Cl-
+  composition: {Cl: 1, E: 1}
+  equation-of-state:
+    special-species: true
+    model: ions-from-neutral-molecule
+    multipliers: {KCl(l): 1.5}

test/thermo/phaseConstructors.cpp

@@ -131,7 +131,8 @@ TEST(IonsFromNeutralConstructor, fromScratch)
 
     auto sKp = make_shared<Species>("K+", parseCompString("K:1"), 1);
     auto sClm = make_shared<Species>("Cl-", parseCompString("Cl:1"), -1);
-    sClm->extra["special_species"] = true;
+    sClm->input["equation-of-state"]["special-species"] = true;
+    sClm->input["equation-of-state"]["model"] = "ions-from-neutral-molecule";
     p.addSpecies(sKp);
     p.addSpecies(sClm);
     std::unique_ptr<PDSS_IonsFromNeutral> ssKp(new PDSS_IonsFromNeutral());

test/thermo/thermoFromYaml.cpp

@@ -191,3 +191,20 @@ TEST(ThermoFromYaml, DebyeHuckel_beta_ij)
         EXPECT_NEAR(mu_ss[k], mu_ss_ref[k], 1e3);
     }
 }
+
+TEST(ThermoFromYaml, IonsFromNeutral)
+{
+    AnyMap infile = AnyMap::fromYamlFile("thermo-models.yaml");
+    auto phaseNodes = infile["phases"].asMap("name");
+    auto thermo = newPhase(*phaseNodes.at("ions-from-neutral-molecule"), infile);
+
+    ASSERT_EQ((int) thermo->nSpecies(), 2);
+    vector_fp mu(thermo->nSpecies());
+    thermo->getChemPotentials(mu.data());
+
+    // Values for regression testing only -- same as "fromScratch" test
+    EXPECT_NEAR(thermo->density(), 1984.3225978174073, 1e-6);
+    EXPECT_NEAR(thermo->enthalpy_mass(), -14737778.668383721, 1e-6);
+    EXPECT_NEAR(mu[0], -4.66404010e+08, 1e1);
+    EXPECT_NEAR(mu[1], -2.88157298e+06, 1e-1);
+}