research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Thermo] Make IonsFromNeutral constructible without XML

Browse source
This commit is contained in:
Ray Speth 2017-08-02 16:53:12 -04:00
parent 7fab7f3cd3
commit e3afaf5e61
6 changed files with 193 additions and 213 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

9
include/cantera/thermo/IonsFromNeutralVPSSTP.h
View file

@ -271,18 +271,14 @@ public:
//@}
virtual bool addSpecies(shared_ptr<Species> spec);
void setNeutralMoleculePhase(shared_ptr<ThermoPhase> neutral);
shared_ptr<ThermoPhase> getNeutralMoleculePhase();
virtual void initThermo();
virtual void initThermoXML(XML_Node& phaseNode, const std::string& id);
virtual void setParametersFromXML(const XML_Node& thermoNode);
private:
//! Initialize lengths of local variables after all species have
//! been identified.
void initLengths();
//! Update the activity coefficients
/*!
* This function will be called to update the internally stored natural
@ -354,9 +350,6 @@ protected:
//! Index of special species
size_t indexSpecialSpecies_;
//! Index of special species
size_t indexSecondSpecialSpecies_;
//! Formula Matrix for composition of neutral molecules
//! in terms of the molecules in this ThermoPhase
/*!

6
include/cantera/thermo/PDSS_IonsFromNeutral.h
View file

@ -90,6 +90,8 @@ public:
virtual bool useSTITbyPDSS() const { return true; }
void setNeutralSpeciesMultiplier(const std::string& species, double mult);
void setSpecialSpecies(bool special=true);
void setParametersFromXML(const XML_Node& speciesNode);
virtual void initThermo();
//@}
@ -100,7 +102,6 @@ protected:
std::map<std::string, double> neutralSpeciesMultipliers_;
public:
//! Number of neutral molecule species that make up the stoichiometric
//! vector for this species, in terms of calculating thermodynamic functions
size_t numMult_;
@ -120,9 +121,6 @@ public:
//! Vector of length equal to the number of species in the neutral molecule phase
mutable vector_fp tmpNM;
//! True if this species is the special species
int specialSpecies_;
};
}

311
src/thermo/IonsFromNeutralVPSSTP.cpp
View file

@ -30,7 +30,6 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() :
ionSolnType_(cIonSolnType_SINGLEANION),
numNeutralMoleculeSpecies_(0),
indexSpecialSpecies_(npos),
indexSecondSpecialSpecies_(npos),
neutralMoleculePhase_(0),
geThermo(0)
{
@ -40,8 +39,7 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const std::string& inputFile,
const std::string& id_) :
ionSolnType_(cIonSolnType_SINGLEANION),
numNeutralMoleculeSpecies_(0),
indexSpecialSpecies_(npos),
indexSecondSpecialSpecies_(npos)
indexSpecialSpecies_(npos)
{
initThermoFile(inputFile, id_);
}
@ -50,8 +48,7 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot,
const std::string& id_) :
ionSolnType_(cIonSolnType_SINGLEANION),
numNeutralMoleculeSpecies_(0),
indexSpecialSpecies_(npos),
indexSecondSpecialSpecies_(npos)
indexSpecialSpecies_(npos)
{
importPhase(phaseRoot, this);
}
@ -450,47 +447,6 @@ void IonsFromNeutralVPSSTP::compositionChanged()
// ------------ Partial Molar Properties of the Solution ------------
void IonsFromNeutralVPSSTP::initThermo()
{
initLengths();
GibbsExcessVPSSTP::initThermo();
}
void IonsFromNeutralVPSSTP::setNeutralMoleculePhase(shared_ptr<ThermoPhase> neutral)
{
neutralMoleculePhase_ = neutral;
geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_.get());
numNeutralMoleculeSpecies_ = neutralMoleculePhase_->nSpecies();
fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk);
NeutralMolecMoleFractions_.resize(numNeutralMoleculeSpecies_);
muNeutralMolecule_.resize(numNeutralMoleculeSpecies_);
lnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
dlnActCoeffdT_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
dlnActCoeffdlnX_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
dlnActCoeffdlnN_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
dlnActCoeffdlnN_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, numNeutralMoleculeSpecies_, 0.0);
y_.resize(numNeutralMoleculeSpecies_, 0.0);
dlnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
dX_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
}
shared_ptr<ThermoPhase> IonsFromNeutralVPSSTP::getNeutralMoleculePhase()
{
return neutralMoleculePhase_;
}
void IonsFromNeutralVPSSTP::initLengths()
{
moleFractions_.resize(m_kk);
fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk);
fm_invert_ionForNeutral.resize(m_kk);
cationList_.resize(m_kk);
anionList_.resize(m_kk);
passThroughList_.resize(m_kk);
moleFractionsTmp_.resize(m_kk);
m_work.resize(m_kk);
}
//! Return the factor overlap
/*!
* @param elnamesVN
@ -524,6 +480,129 @@ static double factorOverlap(const std::vector<std::string>& elnamesVN ,
return fMax;
}
void IonsFromNeutralVPSSTP::initThermo()
{
size_t nElementsN = neutralMoleculePhase_->nElements();
const std::vector<std::string>& elnamesVN = neutralMoleculePhase_->elementNames();
vector_fp elemVectorN(nElementsN);
size_t nElementsI = nElements();
const std::vector<std::string>& elnamesVI = elementNames();
vector_fp elemVectorI(nElementsI);
for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t m = 0; m < nElementsN; m++) {
elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m);
}
for (size_t m = 0; m < nElementsI; m++) {
elemVectorI[m] = nAtoms(indexSpecialSpecies_, m);
}
double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
elnamesVI ,elemVectorI, nElementsI);
if (fac > 0.0) {
for (size_t m = 0; m < nElementsN; m++) {
for (size_t mi = 0; mi < nElementsI; mi++) {
if (elnamesVN[m] == elnamesVI[mi]) {
elemVectorN[m] -= fac * elemVectorI[mi];
}
}
}
}
fm_neutralMolec_ions_[indexSpecialSpecies_ + jNeut * m_kk ] += fac;
for (size_t k = 0; k < m_kk; k++) {
for (size_t m = 0; m < nElementsI; m++) {
elemVectorI[m] = nAtoms(k, m);
}
fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
elnamesVI ,elemVectorI, nElementsI);
if (fac > 0.0) {
for (size_t m = 0; m < nElementsN; m++) {
for (size_t mi = 0; mi < nElementsI; mi++) {
if (elnamesVN[m] == elnamesVI[mi]) {
elemVectorN[m] -= fac * elemVectorI[mi];
}
}
}
bool notTaken = true;
for (size_t iNeut = 0; iNeut < jNeut; iNeut++) {
if (fm_invert_ionForNeutral[k] == iNeut) {
notTaken = false;
}
}
if (notTaken) {
fm_invert_ionForNeutral[k] = jNeut;
} else {
throw CanteraError("IonsFromNeutralVPSSTP::initThermo",
"Simple formula matrix generation failed, one cation is shared between two salts");
}
}
fm_neutralMolec_ions_[k + jNeut * m_kk] += fac;
}
// Ok check the work
for (size_t m = 0; m < nElementsN; m++) {
if (fabs(elemVectorN[m]) > 1.0E-13) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermo",
"Simple formula matrix generation failed");
}
}
}
GibbsExcessVPSSTP::initThermo();
}
void IonsFromNeutralVPSSTP::setNeutralMoleculePhase(shared_ptr<ThermoPhase> neutral)
{
neutralMoleculePhase_ = neutral;
geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_.get());
numNeutralMoleculeSpecies_ = neutralMoleculePhase_->nSpecies();
fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk);
NeutralMolecMoleFractions_.resize(numNeutralMoleculeSpecies_);
muNeutralMolecule_.resize(numNeutralMoleculeSpecies_);
lnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
dlnActCoeffdT_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
dlnActCoeffdlnX_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
dlnActCoeffdlnN_diag_NeutralMolecule_.resize(numNeutralMoleculeSpecies_);
dlnActCoeffdlnN_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, numNeutralMoleculeSpecies_, 0.0);
y_.resize(numNeutralMoleculeSpecies_, 0.0);
dlnActCoeff_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
dX_NeutralMolecule_.resize(numNeutralMoleculeSpecies_, 0.0);
}
shared_ptr<ThermoPhase> IonsFromNeutralVPSSTP::getNeutralMoleculePhase()
{
return neutralMoleculePhase_;
}
bool IonsFromNeutralVPSSTP::addSpecies(shared_ptr<Species> spec)
{
bool added = GibbsExcessVPSSTP::addSpecies(spec);
if (added) {
moleFractions_.push_back(0.0);
moleFractionsTmp_.push_back(0.0);
m_work.push_back(0.0);
fm_invert_ionForNeutral.push_back(npos);
fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk);
if (spec->charge > 0) {
cationList_.push_back(m_kk-1);
} else if (spec->charge < 0) {
anionList_.push_back(m_kk-1);
} else {
passThroughList_.push_back(m_kk-1);
}
if (spec->extra.hasKey("special_species")
&& spec->extra["special_species"].asBool()) {
indexSpecialSpecies_ = m_kk - 1;
}
}
return added;
}
void IonsFromNeutralVPSSTP::setParametersFromXML(const XML_Node& thermoNode)
{
GibbsExcessVPSSTP::setParametersFromXML(thermoNode);
@ -543,146 +622,6 @@ void IonsFromNeutralVPSSTP::setParametersFromXML(const XML_Node& thermoNode)
setNeutralMoleculePhase(shared_ptr<ThermoPhase>(newPhase(*neut_ptr)));
}
void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
if (id_.size() > 0 && phaseNode.id() != id_) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"phasenode and Id are incompatible");
}
// Find the Thermo XML node
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
// Make sure that the thermo model is IonsFromNeutralMolecule
if (!ba::iequals(thermoNode["model"], "ionsfromneutralmolecule")) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"model name isn't IonsFromNeutralMolecule: "
+ thermoNode["model"]);
}
cationList_.clear();
for (size_t k = 0; k < m_kk; k++) {
if (charge(k) > 0) {
cationList_.push_back(k);
}
}
anionList_.clear();
for (size_t k = 0; k < m_kk; k++) {
if (charge(k) < 0) {
anionList_.push_back(k);
}
}
passThroughList_.clear();
for (size_t k = 0; k < m_kk; k++) {
if (charge(k) == 0) {
passThroughList_.push_back(k);
}
}
indexSpecialSpecies_ = npos;
for (size_t k = 0; k < m_kk; k++) {
PDSS_IonsFromNeutral* speciesSS =
dynamic_cast<PDSS_IonsFromNeutral*>(providePDSS(k));
if (!speciesSS) {
throw CanteraError("initThermoXML", "Dynamic cast failed");
}
if (speciesSS->specialSpecies_ == 1) {
indexSpecialSpecies_ = k;
}
if (speciesSS->specialSpecies_ == 2) {
indexSecondSpecialSpecies_ = k;
}
}
size_t nElementsN = neutralMoleculePhase_->nElements();
const std::vector<std::string>& elnamesVN = neutralMoleculePhase_->elementNames();
vector_fp elemVectorN(nElementsN);
vector_fp elemVectorN_orig(nElementsN);
size_t nElementsI = nElements();
const std::vector<std::string>& elnamesVI = elementNames();
vector_fp elemVectorI(nElementsI);
vector_fp fm_tmp(m_kk);
for (size_t k = 0; k < m_kk; k++) {
fm_invert_ionForNeutral[k] = npos;
}
for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t m = 0; m < nElementsN; m++) {
elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m);
}
elemVectorN_orig = elemVectorN;
fm_tmp.assign(m_kk, 0.0);
for (size_t m = 0; m < nElementsI; m++) {
elemVectorI[m] = nAtoms(indexSpecialSpecies_, m);
}
double fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
elnamesVI ,elemVectorI, nElementsI);
if (fac > 0.0) {
for (size_t m = 0; m < nElementsN; m++) {
std::string mName = elnamesVN[m];
for (size_t mi = 0; mi < nElementsI; mi++) {
std::string eName = elnamesVI[mi];
if (mName == eName) {
elemVectorN[m] -= fac * elemVectorI[mi];
}
}
}
}
fm_neutralMolec_ions_[indexSpecialSpecies_ + jNeut * m_kk ] += fac;
for (size_t k = 0; k < m_kk; k++) {
for (size_t m = 0; m < nElementsI; m++) {
elemVectorI[m] = nAtoms(k, m);
}
fac = factorOverlap(elnamesVN, elemVectorN, nElementsN,
elnamesVI ,elemVectorI, nElementsI);
if (fac > 0.0) {
for (size_t m = 0; m < nElementsN; m++) {
std::string mName = elnamesVN[m];
for (size_t mi = 0; mi < nElementsI; mi++) {
std::string eName = elnamesVI[mi];
if (mName == eName) {
elemVectorN[m] -= fac * elemVectorI[mi];
}
}
}
bool notTaken = true;
for (size_t iNeut = 0; iNeut < jNeut; iNeut++) {
if (fm_invert_ionForNeutral[k] == iNeut) {
notTaken = false;
}
}
if (notTaken) {
fm_invert_ionForNeutral[k] = jNeut;
} else {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"Simple formula matrix generation failed, one cation is shared between two salts");
}
}
fm_neutralMolec_ions_[k + jNeut * m_kk] += fac;
}
// Ok check the work
for (size_t m = 0; m < nElementsN; m++) {
if (fabs(elemVectorN[m]) > 1.0E-13) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"Simple formula matrix generation failed");
}
}
}
// This includes the setStateFromXML calls
GibbsExcessVPSSTP::initThermoXML(phaseNode, id_);
}
void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const
{
size_t icat, jNeut;

30
src/thermo/PDSS_IonsFromNeutral.cpp
View file

@ -21,7 +21,6 @@ PDSS_IonsFromNeutral::PDSS_IonsFromNeutral()
: neutralMoleculePhase_(0)
, numMult_(0)
, add2RTln2_(true)
, specialSpecies_(0)
{
}
@ -30,6 +29,16 @@ void PDSS_IonsFromNeutral::setParent(VPStandardStateTP* phase, size_t k)
neutralMoleculePhase_ = dynamic_cast<IonsFromNeutralVPSSTP&>(*phase).getNeutralMoleculePhase();
}
void PDSS_IonsFromNeutral::setNeutralSpeciesMultiplier(const std::string& species, double mult)
{
neutralSpeciesMultipliers_[species] = mult;
numMult_++;
}
void PDSS_IonsFromNeutral::setSpecialSpecies(bool special) {
add2RTln2_ = !special;
}
void PDSS_IonsFromNeutral::setParametersFromXML(const XML_Node& speciesNode)
{
PDSS::setParametersFromXML(speciesNode);
@ -49,21 +58,12 @@ void PDSS_IonsFromNeutral::setParametersFromXML(const XML_Node& speciesNode)
"no Thermo::neutralSpeciesMultipliers Node for species " + speciesNode.name());
}
neutralSpeciesMultipliers_ = parseCompString(nsm->value());
numMult_ = neutralSpeciesMultipliers_.size();
for (auto& species_mult : parseCompString(nsm->value())) {
setNeutralSpeciesMultiplier(species_mult.first, species_mult.second);
}
specialSpecies_ = 0;
const XML_Node* ss = tn->findByName("specialSpecies");
if (ss) {
specialSpecies_ = 1;
}
const XML_Node* sss = tn->findByName("secondSpecialSpecies");
if (sss) {
specialSpecies_ = 2;
}
add2RTln2_ = true;
if (specialSpecies_ == 1) {
add2RTln2_ = false;
if (tn->findByName("specialSpecies")) {
setSpecialSpecies();
}
}

8
src/thermo/Species.cpp
View file

@ -66,6 +66,14 @@ shared_ptr<Species> newSpecies(const XML_Node& species_node)
s->extra["molar_volume"] = getFloat(*stdstate, "molarVolume", "toSI");
}
// Extra data possibly used by IonsFromNeutralVPSSTP
const XML_Node* thermo = species_node.findByName("thermo");
if (thermo && thermo->attrib("model") == "IonFromNeutral") {
if (thermo->hasChild("specialSpecies")) {
s->extra["special_species"] = true;
}
}
return s;
}

42
test/thermo/phaseConstructors.cpp
View file

@ -8,6 +8,7 @@
#include "cantera/thermo/WaterSSTP.h"
#include "cantera/thermo/RedlichKwongMFTP.h"
#include "cantera/thermo/IonsFromNeutralVPSSTP.h"
#include "cantera/thermo/PDSS_IonsFromNeutral.h"
#include "cantera/thermo/IdealSolnGasVPSS.h"
#include "cantera/thermo/IdealMolalSoln.h"
#include "cantera/thermo/DebyeHuckel.h"
@ -112,6 +113,47 @@ TEST(IonsFromNeutralConstructor, fromXML)
EXPECT_NEAR(mu[1], -2.88157298e+06, 1e-1);
}
TEST(IonsFromNeutralConstructor, fromScratch)
{
auto neutral = make_shared<MargulesVPSSTP>();
neutral->addUndefinedElements();
auto sKCl = make_shomate_species("KCl(L)", "K:1 Cl:1", kcl_shomate_coeffs);
neutral->addSpecies(sKCl);
std::unique_ptr<PDSS_ConstVol> ssKCl(new PDSS_ConstVol());
ssKCl->setMolarVolume(0.03757);
neutral->installPDSS(0, std::move(ssKCl));
neutral->initThermo();
IonsFromNeutralVPSSTP p;
p.addUndefinedElements();
p.setNeutralMoleculePhase(neutral);
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;
p.addSpecies(sKp);
p.addSpecies(sClm);
std::unique_ptr<PDSS_IonsFromNeutral> ssKp(new PDSS_IonsFromNeutral());
std::unique_ptr<PDSS_IonsFromNeutral> ssClm(new PDSS_IonsFromNeutral());
ssKp->setNeutralSpeciesMultiplier("KCl(L)", 1.2);
ssClm->setNeutralSpeciesMultiplier("KCl(L)", 1.5);
ssClm->setSpecialSpecies();
p.installPDSS(0, std::move(ssKp));
p.installPDSS(1, std::move(ssClm));
p.initThermo();
ASSERT_EQ((int) p.nSpecies(), 2);
p.setState_TPX(500, 2e5, "K+:0.1, Cl-:0.1");
vector_fp mu(p.nSpecies());
p.getChemPotentials(mu.data());
// Values for regression testing only -- same as XML test
EXPECT_NEAR(p.density(), 1984.3225978174073, 1e-6);
EXPECT_NEAR(p.enthalpy_mass(), -8035317241137.971, 1e-1);
EXPECT_NEAR(mu[0], -4.66404010e+08, 1e1);
EXPECT_NEAR(mu[1], -2.88157298e+06, 1e-1);
}
#ifndef HAS_NO_PYTHON // skip these tests if the Python converter is unavailable
class CtiConversionTest : public testing::Test
{