[Thermo] Allow instantiation of LatticeSolidPhase without XML

include/cantera/thermo/LatticeSolidPhase.h

@@ -406,6 +406,13 @@ public:
     virtual void getGibbs_ref(doublereal* g) const;
 
     virtual bool addSpecies(shared_ptr<Species> spec);
+
+    //! Add a lattice to this phase
+    void addLattice(shared_ptr<ThermoPhase> lattice);
+
+    //! Set the lattice stoichiometric coefficients, \f$ \theta_i \f$
+    void setLatticeStoichiometry(const compositionMap& comp);
+
     virtual void initThermo();
 
     virtual void setParametersFromXML(const XML_Node& eosdata);

src/thermo/LatticeSolidPhase.cpp

@@ -344,6 +344,23 @@ bool LatticeSolidPhase::addSpecies(shared_ptr<Species> spec)
     return added;
 }
 
+void LatticeSolidPhase::addLattice(shared_ptr<ThermoPhase> lattice)
+{
+    m_lattice.push_back(lattice);
+    if (theta_.size() == 0) {
+        theta_.push_back(1.0);
+    } else {
+        theta_.push_back(0.0);
+    }
+}
+
+void LatticeSolidPhase::setLatticeStoichiometry(const compositionMap& comp)
+{
+    for (size_t i = 0; i < m_lattice.size(); i++) {
+        theta_[i] = getValue(comp, m_lattice[i]->name(), 0.0);
+    }
+}
+
 void LatticeSolidPhase::_updateThermo() const
 {
     doublereal tnow = temperature();
@@ -380,29 +397,10 @@ void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata)
     eosdata._require("model","LatticeSolid");
     XML_Node& la = eosdata.child("LatticeArray");
     std::vector<XML_Node*> lattices = la.getChildren("phase");
-    for (size_t n = 0; n < lattices.size(); n++) {
-        m_lattice.emplace_back(newPhase(*lattices[n]));
-    }
-    std::vector<string> pnam;
-    std::vector<string> pval;
-    int np = getPairs(eosdata.child("LatticeStoichiometry"), pnam, pval);
-    theta_.resize(m_lattice.size());
-    for (int i = 0; i < np; i++) {
-        double val = fpValueCheck(pval[i]);
-        bool found = false;
-        for (size_t j = 0; j < m_lattice.size(); j++) {
-            ThermoPhase& tp = *m_lattice[j];
-            string idj = tp.id();
-            if (idj == pnam[i]) {
-                theta_[j] = val;
-                found = true;
-                break;
-            }
-        }
-        if (!found) {
-            throw CanteraError("LatticeSolidPhase::setParametersFromXML", "not found");
-        }
+    for (auto lattice : lattices) {
+        addLattice(shared_ptr<ThermoPhase>(newPhase(*lattice)));
     }
+    setLatticeStoichiometry(parseCompString(eosdata.child("LatticeStoichiometry").value()));
 }
 
 void LatticeSolidPhase::modifyOneHf298SS(const size_t k, const doublereal Hf298New)

test/thermo/phaseConstructors.cpp

@@ -12,7 +12,11 @@
 #include "cantera/thermo/IdealMolalSoln.h"
 #include "cantera/thermo/DebyeHuckel.h"
 #include "cantera/thermo/MargulesVPSSTP.h"
+#include "cantera/thermo/LatticePhase.h"
+#include "cantera/thermo/StoichSubstance.h"
+#include "cantera/thermo/LatticeSolidPhase.h"
 #include "cantera/thermo/NasaPoly2.h"
+#include "cantera/thermo/ConstCpPoly.h"
 #include "cantera/thermo/ShomatePoly.h"
 #include "cantera/thermo/IdealGasPhase.h"
 #include "cantera/thermo/Mu0Poly.h"
@@ -40,6 +44,14 @@ shared_ptr<Species> make_shomate_species(const std::string& name,
     return species;
 }
 
+shared_ptr<Species> make_shomate2_species(const std::string& name,
+     const std::string& composition, const double* shomate_coeffs)
+{
+    auto species = make_shared<Species>(name, parseCompString(composition));
+    species->thermo.reset(new ShomatePoly2(200, 3500, 101325, shomate_coeffs));
+    return species;
+}
+
 shared_ptr<Species> make_species(const std::string& name,
     const std::string& composition, double h298,
     double T1, double mu1, double T2, double mu2)
@@ -50,6 +62,16 @@ shared_ptr<Species> make_species(const std::string& name,
     return species;
 }
 
+shared_ptr<Species> make_const_cp_species(const std::string& name,
+    const std::string& composition, double T0, double h0, double s0, double cp)
+{
+    auto species = make_shared<Species>(name, parseCompString(composition));
+    double coeffs[] = {T0, h0, s0, cp};
+    species->thermo.reset(new ConstCpPoly(200, 3500, 101325, coeffs));
+    return species;
+}
+
+
 class FixedChemPotSstpConstructorTest : public testing::Test
 {
 };
@@ -412,4 +434,49 @@ TEST(MargulesVPSSTP, fromScratch)
     EXPECT_NEAR(p.cp_mole(), 67478.48085733457, 1e-8);
 }
 
+TEST(LatticeSolidPhase, fromScratch)
+{
+    auto base = make_shared<StoichSubstance>();
+    base->addUndefinedElements();
+    base->setName("Li7Si3(S)");
+    base->setDensity(1390.0);
+    auto sLi7Si3 = make_shomate2_species("Li7Si3(S)", "Li:7 Si:3", li7si3_shomate_coeffs);
+    base->addSpecies(sLi7Si3);
+    base->initThermo();
+
+    auto interstital = make_shared<LatticePhase>();
+    interstital->addUndefinedElements();
+    interstital->setName("Li7Si3_Interstitial");
+    auto sLii = make_const_cp_species("Li(i)", "Li:1", 298.15, 0, 2e4, 2e4);
+    auto sVac = make_const_cp_species("V(i)", "", 298.15, 8.98e4, 0, 0);
+    sLii->extra["molar_volume"] = 0.2;
+    interstital->setSiteDensity(10.46344);
+    interstital->addSpecies(sLii);
+    interstital->addSpecies(sVac);
+    interstital->initThermo();
+    interstital->setMoleFractionsByName("Li(i):0.01 V(i):0.99");
+
+    LatticeSolidPhase p;
+    p.addUndefinedElements();
+    p.addLattice(base);
+    p.addLattice(interstital);
+    p.setLatticeStoichiometry(parseCompString("Li7Si3(S):1.0 Li7Si3_Interstitial:1.0"));
+    p.initThermo();
+    p.setState_TP(725, 10 * OneAtm);
+
+    // Regression test based on modified version of Li7Si3_ls.xml
+    EXPECT_NEAR(p.enthalpy_mass(), -2077821.9295456698, 1e-6);
+    double mu_ref[] = {-4.62717474e+08, -4.64248485e+07, 1.16370186e+05};
+    double vol_ref[] = {0.09557086, 0.2, 0.09557086};
+    vector_fp mu(p.nSpecies());
+    vector_fp vol(p.nSpecies());
+    p.getChemPotentials(mu.data());
+    p.getPartialMolarVolumes(vol.data());
+
+    for (size_t k = 0; k < p.nSpecies(); k++) {
+        EXPECT_NEAR(mu[k], mu_ref[k], 1e-7*fabs(mu_ref[k]));
+        EXPECT_NEAR(vol[k], vol_ref[k], 1e-7);
+    }
+}
+
 } // namespace Cantera

test/thermo/thermo_data.h

@@ -47,6 +47,11 @@ const double co_shomate_coeffs[] = {
     35.15070, 1.300095, -0.205921, 0.013550, -3.282780, -127.8375, 231.7120};
 const double co_comp[] = {0.0, 1.0, 1.0};
 
+const double li7si3_shomate_coeffs[] = {
+    700.0, 295.73961, -6.753295, -44.538551, 29.738846, -1.022387, -348.88919,
+    554.35647, 250.51429, 51.125155, -28.341244, 6.242135, 1.346861, -328.46578,
+    498.84106};
+
 // single-region Shomate coefficients
 const double kcl_shomate_coeffs[] = {
     73.59698, 0.0, 0.0, 0.0, 0.0, -443.7341, 175.7209};