research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
cantera/src/thermo/IonsFromNeutralVPSSTP.cpp

1500 lines
47 KiB
C++
Raw Blame History

/**
* @file IonsFromNeutralVPSSTP.cpp
* Definitions for the object which treats ionic liquids as made of ions as species
* even though the thermodynamics is obtained from the neutral molecule representation.
* (see \ref thermoprops
* and class \link Cantera::IonsFromNeutralVPSSTP IonsFromNeutralVPSSTP\endlink).
*
* Header file for a derived class of ThermoPhase that handles
* variable pressure standard state methods for calculating
* thermodynamic properties that are further based upon expressions
* for the excess gibbs free energy expressed as a function of
* the mole fractions.
*/
/*
* Copyright (2009) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
#include "cantera/thermo/IonsFromNeutralVPSSTP.h"
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/thermo/PDSS_IonsFromNeutral.h"
#include "cantera/thermo/mix_defs.h"
#include "cantera/base/stringUtils.h"
#include <cmath>
#include <iomanip>
#include <fstream>
using namespace std;
namespace Cantera
{
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP() :
GibbsExcessVPSSTP(),
ionSolnType_(cIonSolnType_SINGLEANION),
numNeutralMoleculeSpecies_(0),
indexSpecialSpecies_(npos),
indexSecondSpecialSpecies_(npos),
numCationSpecies_(0),
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralMoleculePhase_(0),
geThermo(0),
IOwnNThermoPhase_(true),
moleFractionsTmp_(0),
muNeutralMolecule_(0),
lnActCoeff_NeutralMolecule_(0)
{
}
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const std::string& inputFile,
const std::string& id_,
ThermoPhase* neutralPhase) :
GibbsExcessVPSSTP(),
ionSolnType_(cIonSolnType_SINGLEANION),
numNeutralMoleculeSpecies_(0),
indexSpecialSpecies_(npos),
indexSecondSpecialSpecies_(npos),
numCationSpecies_(0),
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralMoleculePhase_(neutralPhase),
IOwnNThermoPhase_(true),
moleFractionsTmp_(0),
muNeutralMolecule_(0),
lnActCoeff_NeutralMolecule_(0)
{
if (neutralPhase) {
IOwnNThermoPhase_ = false;
}
constructPhaseFile(inputFile, id_);
geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
//y.resize(numNeutralMoleculeSpecies_,0.0);
//size_t numNeutMolSpec = geThermo->nSpecies();
//dlnActCoeff_NeutralMolecule.resize(numNeutMolSpec);
//dX_NeutralMolecule.resize(numNeutMolSpec);
}
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot,
const std::string& id_, ThermoPhase* neutralPhase) :
GibbsExcessVPSSTP(),
ionSolnType_(cIonSolnType_SINGLEANION),
numNeutralMoleculeSpecies_(0),
indexSpecialSpecies_(npos),
indexSecondSpecialSpecies_(npos),
numCationSpecies_(0),
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralMoleculePhase_(neutralPhase),
IOwnNThermoPhase_(true),
moleFractionsTmp_(0),
muNeutralMolecule_(0),
lnActCoeff_NeutralMolecule_(0)
{
if (neutralPhase) {
IOwnNThermoPhase_ = false;
}
constructPhaseXML(phaseRoot, id_);
geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
y_.resize(numNeutralMoleculeSpecies_,0.0);
size_t numNeutMolSpec = geThermo->nSpecies();
dlnActCoeff_NeutralMolecule_.resize(numNeutMolSpec);
dX_NeutralMolecule_.resize(numNeutMolSpec);
}
IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const IonsFromNeutralVPSSTP& b) :
GibbsExcessVPSSTP(),
ionSolnType_(cIonSolnType_SINGLEANION),
numNeutralMoleculeSpecies_(0),
indexSpecialSpecies_(npos),
indexSecondSpecialSpecies_(npos),
numCationSpecies_(0),
numAnionSpecies_(0),
numPassThroughSpecies_(0),
neutralMoleculePhase_(0),
geThermo(0),
IOwnNThermoPhase_(true),
moleFractionsTmp_(0),
muNeutralMolecule_(0),
lnActCoeff_NeutralMolecule_(0)
{
IonsFromNeutralVPSSTP::operator=(b);
}
IonsFromNeutralVPSSTP& IonsFromNeutralVPSSTP::
operator=(const IonsFromNeutralVPSSTP& b)
{
if (&b == this) {
return *this;
}
/*
* If we own the underlying neutral molecule phase, then we do a deep
* copy. If not, we do a shallow copy. We get a valid pointer for
* neutralMoleculePhase_ first, because we need it to assign the pointers
* within the PDSS_IonsFromNeutral object. which is done in the
* GibbsExcessVPSSTP::operator=(b) step.
*/
if (IOwnNThermoPhase_) {
if (b.neutralMoleculePhase_) {
if (neutralMoleculePhase_) {
delete neutralMoleculePhase_;
}
neutralMoleculePhase_ = (b.neutralMoleculePhase_)->duplMyselfAsThermoPhase();
} else {
neutralMoleculePhase_ = 0;
}
} else {
neutralMoleculePhase_ = b.neutralMoleculePhase_;
}
geThermo = dynamic_cast<GibbsExcessVPSSTP*>(neutralMoleculePhase_);
GibbsExcessVPSSTP::operator=(b);
ionSolnType_ = b.ionSolnType_;
numNeutralMoleculeSpecies_ = b.numNeutralMoleculeSpecies_;
indexSpecialSpecies_ = b.indexSpecialSpecies_;
indexSecondSpecialSpecies_ = b.indexSecondSpecialSpecies_;
fm_neutralMolec_ions_ = b.fm_neutralMolec_ions_;
fm_invert_ionForNeutral = b.fm_invert_ionForNeutral;
NeutralMolecMoleFractions_ = b.NeutralMolecMoleFractions_;
cationList_ = b.cationList_;
numCationSpecies_ = b.numCationSpecies_;
anionList_ = b.anionList_;
numAnionSpecies_ = b.numAnionSpecies_;
passThroughList_ = b.passThroughList_;
numPassThroughSpecies_ = b.numPassThroughSpecies_;
y_ = b.y_;
dlnActCoeff_NeutralMolecule_ = b.dlnActCoeff_NeutralMolecule_;
dX_NeutralMolecule_ = b.dX_NeutralMolecule_;
IOwnNThermoPhase_ = b.IOwnNThermoPhase_;
moleFractionsTmp_ = b.moleFractionsTmp_;
muNeutralMolecule_ = b.muNeutralMolecule_;
// gammaNeutralMolecule_ = b.gammaNeutralMolecule_;
lnActCoeff_NeutralMolecule_ = b.lnActCoeff_NeutralMolecule_;
dlnActCoeffdT_NeutralMolecule_ = b.dlnActCoeffdT_NeutralMolecule_;
dlnActCoeffdlnX_diag_NeutralMolecule_ = b.dlnActCoeffdlnX_diag_NeutralMolecule_;
dlnActCoeffdlnN_diag_NeutralMolecule_ = b.dlnActCoeffdlnN_diag_NeutralMolecule_;
dlnActCoeffdlnN_NeutralMolecule_ = b.dlnActCoeffdlnN_NeutralMolecule_;
return *this;
}
IonsFromNeutralVPSSTP::~IonsFromNeutralVPSSTP()
{
if (IOwnNThermoPhase_) {
delete neutralMoleculePhase_;
neutralMoleculePhase_ = 0;
}
}
ThermoPhase*
IonsFromNeutralVPSSTP::duplMyselfAsThermoPhase() const
{
return new IonsFromNeutralVPSSTP(*this);
}
void IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::string id_)
{
if (inputFile.size() == 0) {
throw CanteraError("MargulesVPSSTP:constructPhaseFile",
"input file is null");
}
string path = findInputFile(inputFile);
std::ifstream fin(path.c_str());
if (!fin) {
throw CanteraError("MargulesVPSSTP:constructPhaseFile","could not open "
+path+" for reading.");
}
/*
* The phase object automatically constructs an XML object.
* Use this object to store information.
*/
XML_Node& phaseNode_XML = xml();
XML_Node* fxml = new XML_Node();
fxml->build(fin);
XML_Node* fxml_phase = findXMLPhase(fxml, id_);
if (!fxml_phase) {
throw CanteraError("MargulesVPSSTP:constructPhaseFile",
"ERROR: Can not find phase named " +
id_ + " in file named " + inputFile);
}
fxml_phase->copy(&phaseNode_XML);
constructPhaseXML(*fxml_phase, id_);
delete fxml;
}
void IonsFromNeutralVPSSTP::constructPhaseXML(XML_Node& phaseNode, std::string id_)
{
string stemp;
if (id_.size() > 0) {
string idp = phaseNode.id();
if (idp != id_) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"phasenode and Id are incompatible");
}
}
/*
* Find the thermo XML node
*/
if (!phaseNode.hasChild("thermo")) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"no thermo XML node");
}
XML_Node& thermoNode = phaseNode.child("thermo");
/*
* Make sure that the thermo model is IonsFromNeutralMolecule
*/
stemp = thermoNode.attrib("model");
string formString = lowercase(stemp);
if (formString != "ionsfromneutralmolecule") {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"model name isn't IonsFromNeutralMolecule: " + formString);
}
/*
* Find the Neutral Molecule Phase
*/
if (!thermoNode.hasChild("neutralMoleculePhase")) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"no neutralMoleculePhase XML node");
}
XML_Node& neutralMoleculeNode = thermoNode.child("neutralMoleculePhase");
string nsource = neutralMoleculeNode["datasrc"];
XML_Node* neut_ptr = get_XML_Node(nsource, 0);
if (!neut_ptr) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"neut_ptr = 0");
}
/*
* Create the neutralMolecule ThermoPhase if we haven't already
*/
if (!neutralMoleculePhase_) {
neutralMoleculePhase_ = newPhase(*neut_ptr);
}
/*
* Call the Cantera importPhase() function. This will import
* all of the species into the phase. This will also handle
* all of the solvent and solute standard states
*/
bool m_ok = importPhase(phaseNode, this);
if (!m_ok) {
throw CanteraError("IonsFromNeutralVPSSTP::constructPhaseXML",
"importPhase failed ");
}
}
/*
* -------------- Utilities -------------------------------
*/
int IonsFromNeutralVPSSTP::eosType() const
{
return cIonsFromNeutral;
}
/*
* ------------ Molar Thermodynamic Properties ----------------------
*/
doublereal IonsFromNeutralVPSSTP::enthalpy_mole() const
{
getPartialMolarEnthalpies(DATA_PTR(m_pp));
return mean_X(DATA_PTR(m_pp));
}
doublereal IonsFromNeutralVPSSTP::intEnergy_mole() const
{
double hh = enthalpy_mole();
double pres = pressure();
double molarV = 1.0/molarDensity();
double uu = hh - pres * molarV;
return uu;
}
doublereal IonsFromNeutralVPSSTP::entropy_mole() const
{
getPartialMolarEntropies(DATA_PTR(m_pp));
return mean_X(DATA_PTR(m_pp));
}
doublereal IonsFromNeutralVPSSTP::gibbs_mole() const
{
getChemPotentials(DATA_PTR(m_pp));
return mean_X(DATA_PTR(m_pp));
}
doublereal IonsFromNeutralVPSSTP::cp_mole() const
{
getPartialMolarCp(DATA_PTR(m_pp));
double val = mean_X(DATA_PTR(m_pp));
return val;
}
doublereal IonsFromNeutralVPSSTP::cv_mole() const
{
// Need to revisit this, as it is wrong
getPartialMolarCp(DATA_PTR(m_pp));
return mean_X(DATA_PTR(m_pp));
//err("not implemented");
//return 0.0;
}
/*
* - Activities, Standard States, Activity Concentrations -----------
*/
void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs,
vector_fp& charges, std::vector<size_t>& neutMolIndex) const
{
coeffs = fm_neutralMolec_ions_;
charges = m_speciesCharge;
neutMolIndex = fm_invert_ionForNeutral;
}
void IonsFromNeutralVPSSTP::getActivityCoefficients(doublereal* ac) const
{
// This stuff has moved to the setState routines
// calcNeutralMoleculeMoleFractions();
// neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_));
// neutralMoleculePhase_->getStandardChemPotentials(DATA_PTR(muNeutralMolecule_));
/*
* Update the activity coefficients
*/
s_update_lnActCoeff();
/*
* take the exp of the internally stored coefficients.
*/
for (size_t k = 0; k < m_kk; k++) {
ac[k] = exp(lnActCoeff_Scaled_[k]);
}
}
/*
* --------- Partial Molar Properties of the Solution -------------
*/
void IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const
{
size_t icat, jNeut;
doublereal xx, fact2;
/*
* Transfer the mole fractions to the slave neutral molecule
* phase
* Note we may move this in the future.
*/
//calcNeutralMoleculeMoleFractions();
//neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_));
/*
* Get the standard chemical potentials of netural molecules
*/
neutralMoleculePhase_->getStandardChemPotentials(DATA_PTR(muNeutralMolecule_));
doublereal RT_ = GasConstant * temperature();
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
neutralMoleculePhase_->getChemPotentials(mu);
break;
case cIonSolnType_SINGLEANION:
// neutralMoleculePhase_->getActivityCoefficients(DATA_PTR(gammaNeutralMolecule_));
neutralMoleculePhase_->getLnActivityCoefficients(DATA_PTR(lnActCoeff_NeutralMolecule_));
fact2 = 2.0 * RT_ * log(2.0);
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
xx = std::max(SmallNumber, moleFractions_[icat]);
mu[icat] = muNeutralMolecule_[jNeut] + fact2 + RT_ * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
}
// Do the anion list
icat = anionList_[0];
jNeut = fm_invert_ionForNeutral[icat];
xx = std::max(SmallNumber, moleFractions_[icat]);
mu[icat] = RT_ * log(xx);
// Do the list of neutral molecules
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
icat = passThroughList_[k];
jNeut = fm_invert_ionForNeutral[icat];
xx = std::max(SmallNumber, moleFractions_[icat]);
mu[icat] = muNeutralMolecule_[jNeut] + RT_ * (lnActCoeff_NeutralMolecule_[jNeut] + log(xx));
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("eosType", "Unknown type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("eosType", "Unknown type");
break;
default:
throw CanteraError("eosType", "Unknown type");
break;
}
}
void IonsFromNeutralVPSSTP::getPartialMolarEnthalpies(doublereal* hbar) const
{
/*
* Get the nondimensional standard state enthalpies
*/
getEnthalpy_RT(hbar);
/*
* dimensionalize it.
*/
double T = temperature();
double RT = GasConstant * T;
for (size_t k = 0; k < m_kk; k++) {
hbar[k] *= RT;
}
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
s_update_lnActCoeff();
s_update_dlnActCoeffdT();
double RTT = RT * T;
for (size_t k = 0; k < m_kk; k++) {
hbar[k] -= RTT * dlnActCoeffdT_Scaled_[k];
}
}
void IonsFromNeutralVPSSTP::getPartialMolarEntropies(doublereal* sbar) const
{
double xx;
/*
* Get the nondimensional standard state entropies
*/
getEntropy_R(sbar);
double T = temperature();
/*
* Update the activity coefficients, This also update the
* internally stored molalities.
*/
s_update_lnActCoeff();
s_update_dlnActCoeffdT();
for (size_t k = 0; k < m_kk; k++) {
xx = std::max(moleFractions_[k], SmallNumber);
sbar[k] += - lnActCoeff_Scaled_[k] -log(xx) - T * dlnActCoeffdT_Scaled_[k];
}
/*
* dimensionalize it.
*/
for (size_t k = 0; k < m_kk; k++) {
sbar[k] *= GasConstant;
}
}
void IonsFromNeutralVPSSTP::getdlnActCoeffdlnX_diag(doublereal* dlnActCoeffdlnX_diag) const
{
s_update_lnActCoeff();
s_update_dlnActCoeff_dlnX_diag();
for (size_t k = 0; k < m_kk; k++) {
dlnActCoeffdlnX_diag[k] = dlnActCoeffdlnX_diag_[k];
}
}
void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN_diag(doublereal* dlnActCoeffdlnN_diag) const
{
s_update_lnActCoeff();
s_update_dlnActCoeff_dlnN_diag();
for (size_t k = 0; k < m_kk; k++) {
dlnActCoeffdlnN_diag[k] = dlnActCoeffdlnN_diag_[k];
}
}
void IonsFromNeutralVPSSTP::getdlnActCoeffdlnN(const size_t ld, doublereal* dlnActCoeffdlnN)
{
s_update_lnActCoeff();
s_update_dlnActCoeff_dlnN();
double* data = & dlnActCoeffdlnN_(0,0);
for (size_t k = 0; k < m_kk; k++) {
for (size_t m = 0; m < m_kk; m++) {
dlnActCoeffdlnN[ld * k + m] = data[m_kk * k + m];
}
}
}
void IonsFromNeutralVPSSTP::setTemperature(const doublereal temp)
{
double p = pressure();
IonsFromNeutralVPSSTP::setState_TP(temp, p);
}
void IonsFromNeutralVPSSTP::setPressure(doublereal p)
{
double t = temperature();
IonsFromNeutralVPSSTP::setState_TP(t, p);
}
void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p)
{
/*
* This is a two phase process. First, we calculate the standard states
* within the neutral molecule phase.
*/
neutralMoleculePhase_->setState_TP(t, p);
VPStandardStateTP::setState_TP(t,p);
/*
* Calculate the partial molar volumes, and then the density of the fluid
*/
//calcDensity();
double dd = neutralMoleculePhase_->density();
Phase::setDensity(dd);
}
void IonsFromNeutralVPSSTP::calcIonMoleFractions(doublereal* const mf) const
{
doublereal fmij;
/*
* Download the neutral mole fraction vector into the
* vector, NeutralMolecMoleFractions_[]
*/
neutralMoleculePhase_->getMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
// Zero the mole fractions
for (size_t k = 0; k < m_kk; k++) {
mf[k] = 0.0;
}
/*
* Use the formula matrix to calculate the relative mole numbers.
*/
for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t k = 0; k < m_kk; k++) {
fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
mf[k] += fmij * NeutralMolecMoleFractions_[jNeut];
}
}
/*
* Normalize the new mole fractions
*/
doublereal sum = 0.0;
for (size_t k = 0; k < m_kk; k++) {
sum += mf[k];
}
for (size_t k = 0; k < m_kk; k++) {
mf[k] /= sum;
}
}
void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const
{
size_t icat, jNeut;
doublereal fmij;
doublereal sum = 0.0;
//! Zero the vector we are trying to find.
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
NeutralMolecMoleFractions_[k] = 0.0;
}
#ifdef DEBUG_MODE
sum = -1.0;
for (size_t k = 0; k < m_kk; k++) {
sum += moleFractions_[k];
}
if (fabs(sum) > 1.0E-11) {
throw CanteraError("IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions",
"molefracts don't sum to one: " + fp2str(sum));
}
#endif
// bool fmSimple = true;
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
for (size_t k = 0; k < m_kk; k++) {
NeutralMolecMoleFractions_[k] = moleFractions_[k];
}
break;
case cIonSolnType_SINGLEANION:
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
NeutralMolecMoleFractions_[k] = 0.0;
}
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
if (jNeut != npos) {
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
AssertTrace(fmij != 0.0);
NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
}
}
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
icat = passThroughList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
NeutralMolecMoleFractions_[jNeut] += moleFractions_[icat] / fmij;
}
#ifdef DEBUG_MODE
for (size_t k = 0; k < m_kk; k++) {
moleFractionsTmp_[k] = moleFractions_[k];
}
for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t k = 0; k < m_kk; k++) {
fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
moleFractionsTmp_[k] -= fmij * NeutralMolecMoleFractions_[jNeut];
}
}
for (size_t k = 0; k < m_kk; k++) {
if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
//! Check to see if we have in fact found the inverse.
if (anionList_[0] != k) {
throw CanteraError("", "neutral molecule calc error");
} else {
//! For the single anion case, we will allow some slippage
if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
throw CanteraError("", "neutral molecule calc error - anion");
}
}
}
}
#endif
// Normalize the Neutral Molecule mole fractions
sum = 0.0;
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
sum += NeutralMolecMoleFractions_[k];
}
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
NeutralMolecMoleFractions_[k] /= sum;
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("eosType", "Unknown type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("eosType", "Unknown type");
break;
default:
throw CanteraError("eosType", "Unknown type");
break;
}
}
void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const dx, doublereal* const dy) const
{
doublereal fmij;
doublereal sumy, sumdy;
//check sum dx = 0
//! Zero the vector we are trying to find.
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
y_[k] = 0.0;
dy[k] = 0.0;
}
// bool fmSimple = true;
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
for (size_t k = 0; k < m_kk; k++) {
dy[k] = dx[k];
}
break;
case cIonSolnType_SINGLEANION:
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
size_t icat = cationList_[k];
size_t jNeut = fm_invert_ionForNeutral[icat];
if (jNeut != npos) {
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
AssertTrace(fmij != 0.0);
const doublereal temp = 1.0/fmij;
dy[jNeut] += dx[icat] * temp;
y_[jNeut] += moleFractions_[icat] * temp;
}
}
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
size_t icat = passThroughList_[k];
size_t jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[ icat + jNeut * m_kk];
const doublereal temp = 1.0/fmij;
dy[jNeut] += dx[icat] * temp;
y_[jNeut] += moleFractions_[icat] * temp;
}
#ifdef DEBUG_MODE_NOT
//check dy sum to zero
for (size_t k = 0; k < m_kk; k++) {
moleFractionsTmp_[k] = dx[k];
}
for (jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
for (size_t k = 0; k < m_kk; k++) {
fmij = fm_neutralMolec_ions_[k + jNeut * m_kk];
moleFractionsTmp_[k] -= fmij * dy[jNeut];
}
}
for (size_t k = 0; k < m_kk; k++) {
if (fabs(moleFractionsTmp_[k]) > 1.0E-13) {
//! Check to see if we have in fact found the inverse.
if (anionList_[0] != k) {
throw CanteraError("", "neutral molecule calc error");
} else {
//! For the single anion case, we will allow some slippage
if (fabs(moleFractionsTmp_[k]) > 1.0E-5) {
throw CanteraError("", "neutral molecule calc error - anion");
}
}
}
}
#endif
// Normalize the Neutral Molecule mole fractions
sumy = 0.0;
sumdy = 0.0;
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
sumy += y_[k];
sumdy += dy[k];
}
sumy = 1.0 / sumy;
for (size_t k = 0; k < numNeutralMoleculeSpecies_; k++) {
dy[k] = dy[k] * sumy - y_[k]*sumdy*sumy*sumy;
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("eosType", "Unknown type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("eosType", "Unknown type");
break;
default:
throw CanteraError("eosType", "Unknown type");
break;
}
}
void IonsFromNeutralVPSSTP::setMassFractions(const doublereal* const y)
{
GibbsExcessVPSSTP::setMassFractions(y);
calcNeutralMoleculeMoleFractions();
neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
}
void IonsFromNeutralVPSSTP::setMassFractions_NoNorm(const doublereal* const y)
{
GibbsExcessVPSSTP::setMassFractions_NoNorm(y);
calcNeutralMoleculeMoleFractions();
neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
}
void IonsFromNeutralVPSSTP::setMoleFractions(const doublereal* const x)
{
GibbsExcessVPSSTP::setMoleFractions(x);
calcNeutralMoleculeMoleFractions();
neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
}
void IonsFromNeutralVPSSTP::setMoleFractions_NoNorm(const doublereal* const x)
{
GibbsExcessVPSSTP::setMoleFractions_NoNorm(x);
calcNeutralMoleculeMoleFractions();
neutralMoleculePhase_->setMoleFractions_NoNorm(DATA_PTR(NeutralMolecMoleFractions_));
}
void IonsFromNeutralVPSSTP::setConcentrations(const doublereal* const c)
{
GibbsExcessVPSSTP::setConcentrations(c);
calcNeutralMoleculeMoleFractions();
neutralMoleculePhase_->setMoleFractions(DATA_PTR(NeutralMolecMoleFractions_));
}
/*
* ------------ Partial Molar Properties of the Solution ------------
*/
doublereal IonsFromNeutralVPSSTP::err(const std::string& msg) const
{
throw CanteraError("IonsFromNeutralVPSSTP","Base class method "
+msg+" called. Equation of state type: "+int2str(eosType()));
return 0;
}
void IonsFromNeutralVPSSTP::initThermo()
{
initLengths();
GibbsExcessVPSSTP::initThermo();
}
void IonsFromNeutralVPSSTP::initLengths()
{
m_kk = nSpecies();
numNeutralMoleculeSpecies_ = neutralMoleculePhase_->nSpecies();
moleFractions_.resize(m_kk);
fm_neutralMolec_ions_.resize(numNeutralMoleculeSpecies_ * m_kk);
fm_invert_ionForNeutral.resize(m_kk);
NeutralMolecMoleFractions_.resize(numNeutralMoleculeSpecies_);
cationList_.resize(m_kk);
anionList_.resize(m_kk);
passThroughList_.resize(m_kk);
moleFractionsTmp_.resize(m_kk);
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);
}
//! Return the factor overlap
/*!
* @param elnamesVN
* @param elemVectorN
* @param nElementsN
* @param elnamesVI
* @param elemVectorI
* @param nElementsI
*/
static double factorOverlap(const std::vector<std::string>& elnamesVN ,
const std::vector<double>& elemVectorN,
const size_t nElementsN,
const std::vector<std::string>& elnamesVI ,
const std::vector<double>& elemVectorI,
const size_t nElementsI)
{
double fMax = 1.0E100;
for (size_t mi = 0; mi < nElementsI; mi++) {
if (elnamesVI[mi] != "E") {
if (elemVectorI[mi] > 1.0E-13) {
double eiNum = elemVectorI[mi];
for (size_t mn = 0; mn < nElementsN; mn++) {
if (elnamesVI[mi] == elnamesVN[mn]) {
if (elemVectorN[mn] <= 1.0E-13) {
return 0.0;
}
fMax = std::min(fMax, elemVectorN[mn]/eiNum);
}
}
}
}
}
return fMax;
}
void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{
string stemp;
if (id_.size() > 0) {
string idp = phaseNode.id();
if (idp != 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
*/
stemp = thermoNode.attrib("model");
string formString = lowercase(stemp);
if (formString != "ionsfromneutralmolecule") {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"model name isn't IonsFromNeutralMolecule: " + formString);
}
/*
* Find the Neutral Molecule Phase
*/
if (!thermoNode.hasChild("neutralMoleculePhase")) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"no neutralMoleculePhase XML node");
}
XML_Node& neutralMoleculeNode = thermoNode.child("neutralMoleculePhase");
string nsource = neutralMoleculeNode["datasrc"];
XML_Node* neut_ptr = get_XML_Node(nsource, 0);
if (!neut_ptr) {
throw CanteraError("IonsFromNeutralVPSSTP::initThermoXML",
"neut_ptr = 0");
}
/*
* Create the neutralMolecule ThermoPhase if we haven't already
*/
if (!neutralMoleculePhase_) {
neutralMoleculePhase_ = newPhase(*neut_ptr);
}
/*
* variables that need to be populated
*
* cationList_
* numCationSpecies_;
*/
numCationSpecies_ = 0;
cationList_.clear();
for (size_t k = 0; k < m_kk; k++) {
if (charge(k) > 0) {
cationList_.push_back(k);
numCationSpecies_++;
}
}
numAnionSpecies_ = 0;
anionList_.clear();
for (size_t k = 0; k < m_kk; k++) {
if (charge(k) < 0) {
anionList_.push_back(k);
numAnionSpecies_++;
}
}
numPassThroughSpecies_= 0;
passThroughList_.clear();
for (size_t k = 0; k < m_kk; k++) {
if (charge(k) == 0) {
passThroughList_.push_back(k);
numPassThroughSpecies_++;
}
}
PDSS_IonsFromNeutral* speciesSS = 0;
indexSpecialSpecies_ = npos;
for (size_t k = 0; k < m_kk; k++) {
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();
std::vector<double> elemVectorN(nElementsN);
std::vector<double> elemVectorN_orig(nElementsN);
size_t nElementsI = nElements();
const std::vector<std::string>& elnamesVI = elementNames();
std::vector<double> elemVectorI(nElementsI);
vector<doublereal> fm_tmp(m_kk);
for (size_t k = 0; k < m_kk; k++) {
fm_invert_ionForNeutral[k] = npos;
}
/* for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) {
fm_invert_ionForNeutral[jNeut] = -1;
}*/
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_);
/*
* There is one extra step here. We assure ourselves that we
* have charge conservation.
*/
}
void IonsFromNeutralVPSSTP::s_update_lnActCoeff() const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficiens of the neutral molecules
*/
neutralMoleculePhase_->getLnActivityCoefficients(DATA_PTR(lnActCoeff_NeutralMolecule_));
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
break;
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
lnActCoeff_Scaled_[icat] = lnActCoeff_NeutralMolecule_[jNeut] / fmij;
}
// Do the anion list
icat = anionList_[0];
jNeut = fm_invert_ionForNeutral[icat];
lnActCoeff_Scaled_[icat]= 0.0;
// Do the list of neutral molecules
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
icat = passThroughList_[k];
jNeut = fm_invert_ionForNeutral[icat];
lnActCoeff_Scaled_[icat] = lnActCoeff_NeutralMolecule_[jNeut];
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
break;
default:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff", "Unimplemented type");
break;
}
}
void IonsFromNeutralVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* const dXds,
doublereal* dlnActCoeffds) const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficients of the neutral molecules
*/
if (!geThermo) {
for (size_t k = 0; k < m_kk; k++) {
dlnActCoeffds[k] = dXds[k] / moleFractions_[k];
}
return;
}
// static vector_fp dlnActCoeff_NeutralMolecule(numNeutMolSpec);
// static vector_fp dX_NeutralMolecule(numNeutMolSpec);
getNeutralMoleculeMoleGrads(DATA_PTR(dXds),DATA_PTR(dX_NeutralMolecule_));
// All mole fractions returned to normal
geThermo->getdlnActCoeffds(dTds, DATA_PTR(dX_NeutralMolecule_), DATA_PTR(dlnActCoeff_NeutralMolecule_));
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
break;
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule_[jNeut]/fmij;
}
// Do the anion list
icat = anionList_[0];
jNeut = fm_invert_ionForNeutral[icat];
dlnActCoeffds[icat]= 0.0;
// Do the list of neutral molecules
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
icat = passThroughList_[k];
jNeut = fm_invert_ionForNeutral[icat];
dlnActCoeffds[icat] = dlnActCoeff_NeutralMolecule_[jNeut];
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
break;
default:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffds", "Unimplemented type");
break;
}
}
void IonsFromNeutralVPSSTP::s_update_dlnActCoeffdT() const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficients of the neutral molecules
*/
if (!geThermo) {
dlnActCoeffdT_Scaled_.assign(m_kk, 0.0);
return;
}
geThermo->getdlnActCoeffdT(DATA_PTR(dlnActCoeffdT_NeutralMolecule_));
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
break;
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
dlnActCoeffdT_Scaled_[icat] = dlnActCoeffdT_NeutralMolecule_[jNeut]/fmij;
}
// Do the anion list
icat = anionList_[0];
jNeut = fm_invert_ionForNeutral[icat];
dlnActCoeffdT_Scaled_[icat]= 0.0;
// Do the list of neutral molecules
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
icat = passThroughList_[k];
jNeut = fm_invert_ionForNeutral[icat];
dlnActCoeffdT_Scaled_[icat] = dlnActCoeffdT_NeutralMolecule_[jNeut];
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
break;
default:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeffdT", "Unimplemented type");
break;
}
}
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnX_diag() const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficients of the neutral molecules
*/
if (!geThermo) {
dlnActCoeffdlnX_diag_.assign(m_kk, 0.0);
return;
}
geThermo->getdlnActCoeffdlnX_diag(DATA_PTR(dlnActCoeffdlnX_diag_NeutralMolecule_));
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
break;
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
dlnActCoeffdlnX_diag_[icat] = dlnActCoeffdlnX_diag_NeutralMolecule_[jNeut]/fmij;
}
// Do the anion list
icat = anionList_[0];
jNeut = fm_invert_ionForNeutral[icat];
dlnActCoeffdlnX_diag_[icat]= 0.0;
// Do the list of neutral molecules
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
icat = passThroughList_[k];
jNeut = fm_invert_ionForNeutral[icat];
dlnActCoeffdlnX_diag_[icat] = dlnActCoeffdlnX_diag_NeutralMolecule_[jNeut];
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
break;
default:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnX_diag()", "Unimplemented type");
break;
}
}
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN_diag() const
{
size_t icat, jNeut;
doublereal fmij;
/*
* Get the activity coefficients of the neutral molecules
*/
if (!geThermo) {
dlnActCoeffdlnN_diag_.assign(m_kk, 0.0);
return;
}
geThermo->getdlnActCoeffdlnN_diag(DATA_PTR(dlnActCoeffdlnN_diag_NeutralMolecule_));
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
break;
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
//! Get the id for the next cation
icat = cationList_[k];
jNeut = fm_invert_ionForNeutral[icat];
fmij = fm_neutralMolec_ions_[icat + jNeut * m_kk];
dlnActCoeffdlnN_diag_[icat] = dlnActCoeffdlnN_diag_NeutralMolecule_[jNeut]/fmij;
}
// Do the anion list
icat = anionList_[0];
jNeut = fm_invert_ionForNeutral[icat];
dlnActCoeffdlnN_diag_[icat]= 0.0;
// Do the list of neutral molecules
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
icat = passThroughList_[k];
jNeut = fm_invert_ionForNeutral[icat];
dlnActCoeffdlnN_diag_[icat] = dlnActCoeffdlnN_diag_NeutralMolecule_[jNeut];
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
break;
default:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN_diag()", "Unimplemented type");
break;
}
}
void IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN() const
{
size_t kcat = 0, kNeut = 0, mcat = 0, mNeut = 0;
doublereal fmij = 0.0, mfmij;
dlnActCoeffdlnN_.zero();
/*
* Get the activity coefficients of the neutral molecules
*/
if (!geThermo) {
throw CanteraError("IonsFromNeutralVPSSTP::s_update_dlnActCoeff_dlnN()", "dynamic cast failed");
}
size_t nsp_ge = geThermo->nSpecies();
geThermo->getdlnActCoeffdlnN(nsp_ge, &(dlnActCoeffdlnN_NeutralMolecule_(0,0)));
switch (ionSolnType_) {
case cIonSolnType_PASSTHROUGH:
break;
case cIonSolnType_SINGLEANION:
// Do the cation list
for (size_t k = 0; k < cationList_.size(); k++) {
for (size_t m = 0; m < cationList_.size(); m++) {
kcat = cationList_[k];
kNeut = fm_invert_ionForNeutral[kcat];
fmij = fm_neutralMolec_ions_[kcat + kNeut * m_kk];
dlnActCoeffdlnN_diag_[kcat] = dlnActCoeffdlnN_diag_NeutralMolecule_[kNeut]/fmij;
mcat = cationList_[m];
mNeut = fm_invert_ionForNeutral[mcat];
mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk];
dlnActCoeffdlnN_(kcat,mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut) * mfmij / fmij;
}
for (size_t m = 0; m < numPassThroughSpecies_; m++) {
mcat = passThroughList_[m];
mNeut = fm_invert_ionForNeutral[mcat];
dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut) / fmij;
}
}
// Do the anion list -> anion activity coefficient is one
kcat = anionList_[0];
kNeut = fm_invert_ionForNeutral[kcat];
for (size_t k = 0; k < m_kk; k++) {
dlnActCoeffdlnN_(kcat, k) = 0.0;
dlnActCoeffdlnN_(k, kcat) = 0.0;
}
// Do the list of neutral molecules
for (size_t k = 0; k < numPassThroughSpecies_; k++) {
kcat = passThroughList_[k];
kNeut = fm_invert_ionForNeutral[kcat];
dlnActCoeffdlnN_diag_[kcat] = dlnActCoeffdlnN_diag_NeutralMolecule_[kNeut];
for (size_t m = 0; m < m_kk; m++) {
mcat = passThroughList_[m];
mNeut = fm_invert_ionForNeutral[mcat];
dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut, mNeut);
}
for (size_t m = 0; m < cationList_.size(); m++) {
mcat = cationList_[m];
mNeut = fm_invert_ionForNeutral[mcat];
mfmij = fm_neutralMolec_ions_[mcat + mNeut * m_kk];
dlnActCoeffdlnN_(kcat, mcat) = dlnActCoeffdlnN_NeutralMolecule_(kNeut,mNeut);
}
}
break;
case cIonSolnType_SINGLECATION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
break;
case cIonSolnType_MULTICATIONANION:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
break;
default:
throw CanteraError("IonsFromNeutralVPSSTP::s_update_lnActCoeff_dlnN", "Unimplemented type");
break;
}
}
}
Reference in a new issue View git blame Copy permalink