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

[Equil] Eliminate SpeciesThermo and VPhaseList from VCS_SOLVE

Browse source 
This means that the VCS_SPECIES_THERMO and vcs_VolPhase classes no longer need
to be able to be copied.
This commit is contained in:
Ray Speth 2017-08-18 23:29:07 -04:00
parent 4e53c893cf
commit 7eb939dc5f
17 changed files with 86 additions and 295 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
include/cantera/equil/vcs_VolPhase.h
View file

@ -83,10 +83,8 @@ class vcs_VolPhase
public:
vcs_VolPhase(VCS_SOLVE* owningSolverObject = 0);
vcs_VolPhase(const vcs_VolPhase& b);
vcs_VolPhase& operator=(const vcs_VolPhase& b);
vcs_VolPhase(const vcs_VolPhase& b) = delete;
vcs_VolPhase& operator=(const vcs_VolPhase& b) = delete;
~vcs_VolPhase();
//! The resize() function fills in all of the initial information if it
@ -529,10 +527,6 @@ private:
private:
//! Backtrack value of VCS_SOLVE *
/*!
* Note the default for this is 0. That's a valid value too, since VCS_PROB
* also uses vcs_VolPhase objects.
*/
VCS_SOLVE* m_owningSolverObject;
public:

12
include/cantera/equil/vcs_solve.h
View file

@ -1026,14 +1026,6 @@ public:
//! length number of species.
vector_fp mf;
//! Array of phase structures
std::vector<vcs_VolPhase*> VPhaseList;
//! Vector of pointers to thermo structures which identify the model and
//! parameters for evaluating the thermodynamic functions for that
//! particular species
std::vector<VCS_SPECIES_THERMO*> SpeciesThermo;
//! Print level for print routines
int m_printLvl;
@ -1446,7 +1438,7 @@ public:
vector_int m_elementActive;
//! Array of Phase Structures. Length = number of phases.
std::vector<vcs_VolPhase*> m_VolPhaseList;
std::vector<std::unique_ptr<vcs_VolPhase>> m_VolPhaseList;
//! This specifies the current state of units for the Gibbs free energy
//! properties in the program.
@ -1529,7 +1521,7 @@ public:
/*!
* SpeciesThermo[k] pointer to the thermo information for the kth species
*/
std::vector<VCS_SPECIES_THERMO*> m_speciesThermoList;
std::vector<std::unique_ptr<VCS_SPECIES_THERMO>> m_speciesThermoList;
//! Choice of Hessians
/*!

6
include/cantera/equil/vcs_species_thermo.h
View file

@ -81,11 +81,7 @@ public:
//! parameter that is used in the VCS_SSVOL_CONSTANT model.
double SSStar_Vol0;
VCS_SPECIES_THERMO(size_t indexPhase, size_t indexSpeciesPhase);
virtual ~VCS_SPECIES_THERMO() {}
//! Duplication function for derived classes.
virtual VCS_SPECIES_THERMO* duplMyselfAsVCS_SPECIES_THERMO();
VCS_SPECIES_THERMO();
};
}

4
src/equil/vcs_Gibbs.cpp
View file

@ -17,7 +17,7 @@ double VCS_SOLVE::vcs_Total_Gibbs(double* molesSp, double* chemPot,
double g = 0.0;
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
if ((TPhInertMoles[iph] > 0.0) && (tPhMoles[iph] > 0.0)) {
g += TPhInertMoles[iph] *
log(TPhInertMoles[iph] / tPhMoles[iph]);
@ -51,7 +51,7 @@ double VCS_SOLVE::vcs_GibbsPhase(size_t iphase, const double* const w,
if (TPhInertMoles[iphase] > 0.0) {
phaseMols += TPhInertMoles[iphase];
g += TPhInertMoles[iphase] * log(TPhInertMoles[iphase] / phaseMols);
vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
if (Vphase->m_gasPhase) {
g += TPhInertMoles[iphase] * log(m_pressurePA/1.01325E5);
}

6
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -484,7 +484,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
plogf(" %15.3e %15.3e ", m_vsolve.w[i], m_vsolve.mf[i]);
if (m_vsolve.w[i] <= 0.0) {
size_t iph = m_vsolve.m_phaseID[i];
vcs_VolPhase* VPhase = m_vsolve.VPhaseList[iph];
vcs_VolPhase* VPhase = m_vsolve.m_VolPhaseList[iph].get();
if (VPhase->nSpecies() > 1) {
plogf(" -1.000e+300\n");
} else {
@ -529,7 +529,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
size_t nSpecies = tref.nSpecies();
VolPM.resize(nSpecies, 0.0);
tref.getPartialMolarVolumes(&VolPM[0]);
vcs_VolPhase* volP = m_vsolve.VPhaseList[iphase];
vcs_VolPhase* volP = m_vsolve.m_VolPhaseList[iphase].get();
double TMolesPhase = volP->totalMoles();
double VolPhaseVolumes = 0.0;
@ -552,7 +552,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
size_t istart = m_mix->speciesIndex(0, iphase);
ThermoPhase& tref = m_mix->phase(iphase);
string phaseName = tref.name();
vcs_VolPhase* volP = m_vsolve.VPhaseList[iphase];
vcs_VolPhase* volP = m_vsolve.m_VolPhaseList[iphase].get();
double TMolesPhase = volP->totalMoles();
size_t nSpecies = tref.nSpecies();
activity.resize(nSpecies, 0.0);

2
src/equil/vcs_TP.cpp
View file

@ -50,7 +50,7 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg
int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
{
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* vph = m_VolPhaseList[iph];
vcs_VolPhase* vph = m_VolPhaseList[iph].get();
vph->setState_TP(m_temperature, m_pressurePA);
vph->sendToVCS_GStar(&m_SSfeSpecies[0]);
}

138
src/equil/vcs_VolPhase.cpp
View file

@ -58,126 +58,6 @@ vcs_VolPhase::~vcs_VolPhase()
}
}
vcs_VolPhase::vcs_VolPhase(const vcs_VolPhase& b) :
m_owningSolverObject(b.m_owningSolverObject),
VP_ID_(b.VP_ID_),
m_singleSpecies(b.m_singleSpecies),
m_gasPhase(b.m_gasPhase),
m_eqnState(b.m_eqnState),
ChargeNeutralityElement(b.ChargeNeutralityElement),
p_activityConvention(b.p_activityConvention),
m_numElemConstraints(b.m_numElemConstraints),
m_numSpecies(b.m_numSpecies),
m_totalMolesInert(b.m_totalMolesInert),
m_isIdealSoln(b.m_isIdealSoln),
m_existence(b.m_existence),
m_MFStartIndex(b.m_MFStartIndex),
TP_ptr(b.TP_ptr),
v_totalMoles(b.v_totalMoles),
creationMoleNumbers_(0),
creationGlobalRxnNumbers_(0),
m_phiVarIndex(npos),
m_totalVol(b.m_totalVol),
m_vcsStateStatus(VCS_STATECALC_OLD),
m_phi(b.m_phi),
m_UpToDate(false),
m_UpToDate_AC(false),
m_UpToDate_VolStar(false),
m_UpToDate_VolPM(false),
m_UpToDate_GStar(false),
m_UpToDate_G0(false),
Temp_(b.Temp_),
Pres_(b.Pres_)
{
//! Objects that are owned by this object are deep copied here, except for
//! the ThermoPhase object. The assignment operator does most of the work.
*this = b;
}
vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b)
{
if (&b != this) {
size_t old_num = m_numSpecies;
// Note: we comment this out for the assignment operator
// specifically, because it isn't true for the assignment
// operator but is true for a copy constructor
// m_owningSolverObject = b.m_owningSolverObject;
VP_ID_ = b.VP_ID_;
m_singleSpecies = b.m_singleSpecies;
m_gasPhase = b.m_gasPhase;
m_eqnState = b.m_eqnState;
ChargeNeutralityElement = b.ChargeNeutralityElement;
p_activityConvention= b.p_activityConvention;
m_numSpecies = b.m_numSpecies;
m_numElemConstraints = b.m_numElemConstraints;
m_elementNames.resize(b.m_numElemConstraints);
for (size_t e = 0; e < b.m_numElemConstraints; e++) {
m_elementNames[e] = b.m_elementNames[e];
}
m_elementActive = b.m_elementActive;
m_elementType = b.m_elementType;
m_formulaMatrix = b.m_formulaMatrix;
m_speciesUnknownType = b.m_speciesUnknownType;
m_elemGlobalIndex = b.m_elemGlobalIndex;
PhaseName = b.PhaseName;
m_totalMolesInert = b.m_totalMolesInert;
m_isIdealSoln = b.m_isIdealSoln;
m_existence = b.m_existence;
m_MFStartIndex = b.m_MFStartIndex;
// Do a shallow copy because we haven' figured this out.
IndSpecies = b.IndSpecies;
for (size_t k = 0; k < old_num; k++) {
if (ListSpeciesPtr[k]) {
delete ListSpeciesPtr[k];
ListSpeciesPtr[k] = 0;
}
}
ListSpeciesPtr.resize(m_numSpecies, 0);
for (size_t k = 0; k < m_numSpecies; k++) {
ListSpeciesPtr[k] =
new vcs_SpeciesProperties(*(b.ListSpeciesPtr[k]));
}
// Do a shallow copy of the ThermoPhase object pointer. We don't
// duplicate the object.
//
// Um, there is no reason we couldn't do a
// duplicateMyselfAsThermoPhase() call here. This will have to be looked
// into.
TP_ptr = b.TP_ptr;
v_totalMoles = b.v_totalMoles;
Xmol_ = b.Xmol_;
creationMoleNumbers_ = b.creationMoleNumbers_;
creationGlobalRxnNumbers_ = b.creationGlobalRxnNumbers_;
m_phiVarIndex = b.m_phiVarIndex;
m_totalVol = b.m_totalVol;
SS0ChemicalPotential = b.SS0ChemicalPotential;
StarChemicalPotential = b.StarChemicalPotential;
StarMolarVol = b.StarMolarVol;
PartialMolarVol = b.PartialMolarVol;
ActCoeff = b.ActCoeff;
np_dLnActCoeffdMolNumber = b.np_dLnActCoeffdMolNumber;
m_vcsStateStatus = b.m_vcsStateStatus;
m_phi = b.m_phi;
m_UpToDate = false;
m_UpToDate_AC = false;
m_UpToDate_VolStar = false;
m_UpToDate_VolPM = false;
m_UpToDate_GStar = false;
m_UpToDate_G0 = false;
Temp_ = b.Temp_;
Pres_ = b.Pres_;
setState_TP(Temp_, Pres_);
_updateMoleFractionDependencies();
}
return *this;
}
void vcs_VolPhase::resize(const size_t phaseNum, const size_t nspecies,
const size_t numElem, const char* const phaseName,
const double molesInert)
@ -416,15 +296,15 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc,
throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
}
} else if (m_owningSolverObject) {
if (stateCalc == VCS_STATECALC_OLD) {
if (molesSpeciesVCS != &m_owningSolverObject->m_molNumSpecies_old[0]) {
throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
}
} else if (stateCalc == VCS_STATECALC_NEW) {
if (molesSpeciesVCS != &m_owningSolverObject->m_molNumSpecies_new[0]) {
throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
}
}
// if (stateCalc == VCS_STATECALC_OLD) {
// if (molesSpeciesVCS != &m_owningSolverObject->m_molNumSpecies_old[0]) {
// throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
// }
// } else if (stateCalc == VCS_STATECALC_NEW) {
// if (molesSpeciesVCS != &m_owningSolverObject->m_molNumSpecies_new[0]) {
// throw CanteraError("vcs_VolPhase::setMolesFromVCS", "shouldn't be here");
// }
// }
}
for (size_t k = 0; k < m_numSpecies; k++) {

2
src/equil/vcs_elem_rearrange.cpp
View file

@ -143,7 +143,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
// Change the element Global Index list in each vcs_VolPhase object
// to reflect the switch in the element positions.
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* volPhase = m_VolPhaseList[iph];
vcs_VolPhase* volPhase = m_VolPhaseList[iph].get();
for (size_t e = 0; e < volPhase->nElemConstraints(); e++) {
if (volPhase->elemGlobalIndex(e) == ipos) {
volPhase->setElemGlobalIndex(e, jpos);

4
src/equil/vcs_nondim.cpp
View file

@ -78,7 +78,7 @@ void VCS_SOLVE::vcs_nondim_TP()
for (size_t iph = 0; iph < m_numPhases; iph++) {
TPhInertMoles[iph] *= (1.0 / m_totalMoleScale);
if (TPhInertMoles[iph] != 0.0) {
vcs_VolPhase* vphase = m_VolPhaseList[iph];
vcs_VolPhase* vphase = m_VolPhaseList[iph].get();
vphase->setTotalMolesInert(TPhInertMoles[iph]);
}
}
@ -119,7 +119,7 @@ void VCS_SOLVE::vcs_redim_TP()
for (size_t iph = 0; iph < m_numPhases; iph++) {
TPhInertMoles[iph] *= m_totalMoleScale;
if (TPhInertMoles[iph] != 0.0) {
vcs_VolPhase* vphase = m_VolPhaseList[iph];
vcs_VolPhase* vphase = m_VolPhaseList[iph].get();
vphase->setTotalMolesInert(TPhInertMoles[iph]);
}
}

10
src/equil/vcs_phaseStability.cpp
View file

@ -19,7 +19,7 @@ namespace Cantera
bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const
{
vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop];
vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop].get();
AssertThrowMsg(!Vphase->exists(), "VCS_SOLVE::vcs_popPhasePossible",
"called for a phase that exists!");
@ -119,7 +119,7 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
plogf(" --------------------------------------------------------------------------\n");
}
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
int existence = Vphase->exists();
strcpy(anote, "");
if (existence > 0) {
@ -199,7 +199,7 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
{
vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop];
vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop].get();
// Identify the first species in the phase
size_t kspec = Vphase->spGlobalIndexVCS(0);
// Identify the formation reaction for that species
@ -226,7 +226,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
}
}
for (size_t j = 0; j < m_numPhases; j++) {
Vphase = m_VolPhaseList[j];
Vphase = m_VolPhaseList[j].get();
if (! Vphase->m_singleSpecies && m_tPhaseMoles_old[j] > 0.0) {
s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j];
}
@ -349,7 +349,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
{
// We will use the _new state calc here
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
const size_t nsp = Vphase->nSpecies();
int minNumberIterations = 3;
if (nsp <= 1) {

6
src/equil/vcs_prep.cpp
View file

@ -23,7 +23,7 @@ void VCS_SOLVE::vcs_SSPhase()
// earmarked as a multispecies phase. Treat that species as a single-species
// phase
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
Vphase->m_singleSpecies = false;
if (TPhInertMoles[iph] > 0.0) {
Vphase->setExistence(2);
@ -38,7 +38,7 @@ void VCS_SOLVE::vcs_SSPhase()
// indicating whether a species is in a single species phase or not.
for (size_t kspec = 0; kspec < m_nsp; kspec++) {
size_t iph = m_phaseID[kspec];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
if (Vphase->m_singleSpecies) {
m_SSPhase[kspec] = true;
} else {
@ -72,7 +72,7 @@ int VCS_SOLVE::vcs_prep(int printLvl)
for (size_t kspec = 0; kspec < m_nsp; ++kspec) {
size_t pID = m_phaseID[kspec];
size_t spPhIndex = m_speciesLocalPhaseIndex[kspec];
vcs_VolPhase* vPhase = m_VolPhaseList[pID];
vcs_VolPhase* vPhase = m_VolPhaseList[pID].get();
vcs_SpeciesProperties* spProp = vPhase->speciesProperty(spPhIndex);
double sz = 0.0;
size_t eSize = spProp->FormulaMatrixCol.size();

6
src/equil/vcs_prob.cpp
View file

@ -44,7 +44,7 @@ void VCS_SOLVE::prob_report(int print_lvl)
plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_Moles Species_Type\n");
for (size_t i = 0; i < m_nsp; i++) {
vcs_VolPhase* Vphase = VPhaseList[m_phaseID[i]];
vcs_VolPhase* Vphase = m_VolPhaseList[m_phaseID[i]].get();
plogf("%16s %5d %16s", m_mix->speciesName(i), m_phaseID[i],
Vphase->PhaseName);
if (m_doEstimateEquil >= 0) {
@ -70,7 +70,7 @@ void VCS_SOLVE::prob_report(int print_lvl)
plogf(" TMolesInert TKmoles\n");
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* Vphase = VPhaseList[iphase];
vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
plogf("%16s %5d %5d %8d ", Vphase->PhaseName,
Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase);
plogf("%16s %8d %16e ", Vphase->eos_name(),
@ -95,7 +95,7 @@ void VCS_SOLVE::prob_report(int print_lvl)
plogf(" Species (phase) "
" SS0ChemPot StarChemPot\n");
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* Vphase = VPhaseList[iphase];
vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
Vphase->setState_TP(m_temperature, m_pressurePA);
for (size_t kindex = 0; kindex < Vphase->nSpecies(); kindex++) {
size_t kglob = Vphase->spGlobalIndexVCS(kindex);

2
src/equil/vcs_report.cpp
View file

@ -196,7 +196,7 @@ int VCS_SOLVE::vcs_report(int iconv)
writeline('-', m_nelem*10 + 58);
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
plogf(" %3d ", iphase);
vcs_VolPhase* VPhase = m_VolPhaseList[iphase];
vcs_VolPhase* VPhase = m_VolPhaseList[iphase].get();
plogf("%-12.12s |",VPhase->PhaseName);
plogf("%10.3e |", m_tPhaseMoles_old[iphase]*molScale);
totalMoles += m_tPhaseMoles_old[iphase];

6
src/equil/vcs_rxnadj.cpp
View file

@ -63,7 +63,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
size_t iph = m_phaseID[kspec];
double tphmoles = m_tPhaseMoles_old[iph];
double trphmoles = tphmoles / m_totalMolNum;
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
if (Vphase->exists() && (trphmoles > VCS_DELETE_PHASE_CUTOFF)) {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
@ -134,7 +134,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
}
}
for (size_t j = 0; j < m_numPhases; j++) {
vcs_VolPhase* Vphase = m_VolPhaseList[j];
vcs_VolPhase* Vphase = m_VolPhaseList[j].get();
if (!Vphase->m_singleSpecies && m_tPhaseMoles_old[j] > 0.0) {
s -= pow(m_deltaMolNumPhase(j,irxn), 2) / m_tPhaseMoles_old[j];
}
@ -346,7 +346,7 @@ void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS)
{
// Loop over all of the phases in the problem
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
// We don't need to call single species phases;
if (!Vphase->m_singleSpecies && !Vphase->isIdealSoln()) {

124
src/equil/vcs_solve.cpp
View file

@ -52,13 +52,9 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_gibbsSpecies.resize(m_nsp, 0.0);
w.resize(m_nsp, 0.0);
mf.resize(m_nsp, 0.0);
SpeciesThermo.resize(m_nsp,0);
m_speciesThermoList.resize(m_nsp);
for (size_t kspec = 0; kspec < m_nsp; kspec++) {
SpeciesThermo[kspec] = new VCS_SPECIES_THERMO(0, 0);
}
VPhaseList.resize(m_numPhases, 0);
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
VPhaseList[iphase] = new vcs_VolPhase();
m_speciesThermoList[kspec].reset(new VCS_SPECIES_THERMO());
}
string ser = "VCS_SOLVE: ERROR:\n\t";
@ -124,12 +120,11 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_actCoeffSpecies_old.resize(m_nsp, 1.0);
m_wtSpecies.resize(m_nsp, 0.0);
m_chargeSpecies.resize(m_nsp, 0.0);
m_speciesThermoList.resize(m_nsp, (VCS_SPECIES_THERMO*)0);
// Phase Info
m_VolPhaseList.resize(m_numPhases, 0);
m_VolPhaseList.resize(m_numPhases);
for (size_t iph = 0; iph < m_numPhases; iph++) {
m_VolPhaseList[iph] = new vcs_VolPhase(this);
m_VolPhaseList[iph].reset(new vcs_VolPhase(this));
}
// For Future expansion
@ -173,7 +168,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
// ->NumSpecies = number of species in the phase
// ->TMolesInert = Inerts in the phase = 0.0 for cantera
// ->PhaseName = Name of the phase
vcs_VolPhase* VolPhase = VPhaseList[iphase];
vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0);
VolPhase->m_gasPhase = gasPhase;
@ -263,7 +258,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
addOnePhaseSpecies(VolPhase, k, kT);
// Get a pointer to the thermo object
ts_ptr = SpeciesThermo[kT];
ts_ptr = m_speciesThermoList[kT].get();
// Fill in the vcs_SpeciesProperty structure
vcs_SpeciesProperties* sProp = VolPhase->speciesProperty(k);
@ -384,7 +379,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
for (size_t i = 0; i < m_nsp; i++) {
size_t iphase = m_phaseID[i];
vcs_VolPhase* VolPhase = VPhaseList[iphase];
vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase,
VolPhase->PhaseName.c_str());
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
@ -401,7 +396,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
plogf(" TMolesInert Tmoles(kmol)\n");
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* VolPhase = VPhaseList[iphase];
vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID_, VolPhase->m_singleSpecies,
VolPhase->m_gasPhase, VolPhase->eos_name(),
@ -417,23 +412,12 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
plogf("\n");
}
// Copy the VCS_SPECIES_THERMO structures
for (size_t kspec = 0; kspec < m_nsp; kspec++) {
delete m_speciesThermoList[kspec];
VCS_SPECIES_THERMO* spf = SpeciesThermo[kspec];
m_speciesThermoList[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO();
if (m_speciesThermoList[kspec] == NULL) {
throw CanteraError("VCS_SOLVE::VCS_SOLVE",
" duplMyselfAsVCS_SPECIES_THERMO returned an error!");
}
}
// Copy the equilibrium mole number estimate
m_molNumSpecies_old = w;
// TPhInertMoles[] -> must be copied over here
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = VPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
TPhInertMoles[iph] = Vphase->totalMolesInert();
}
@ -464,7 +448,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
numPhSp[iph]++;
}
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = VPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
if (numPhSp[iph] != Vphase->nSpecies()) {
throw CanteraError("VCS_SOLVE::VCS_SOLVE",
"Number of species in phase {}, {}, doesn't match ({} != {}) [vphase = {}]",
@ -496,24 +480,9 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_speciesName[i] = m_mix->speciesName(i);
}
// Copy over all of the phase information. Use the object's assignment
// operator
for (size_t iph = 0; iph < m_numPhases; iph++) {
*m_VolPhaseList[iph] = *VPhaseList[iph];
// Fix up the species thermo pointer in the vcs_SpeciesThermo object. It
// should point to the species thermo pointer in the private data space.
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
for (size_t k = 0; k < Vphase->nSpecies(); k++) {
vcs_SpeciesProperties* sProp = Vphase->speciesProperty(k);
size_t kT = Vphase->spGlobalIndexVCS(k);
sProp->SpeciesThermo = m_speciesThermoList[kT];
}
}
// Specify the Activity Convention information
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
m_phaseActConvention[iph] = Vphase->p_activityConvention;
if (Vphase->p_activityConvention != 0) {
// We assume here that species 0 is the solvent. The solvent isn't
@ -541,25 +510,6 @@ VCS_SOLVE::~VCS_SOLVE()
void VCS_SOLVE::vcs_delete_memory()
{
for (size_t j = 0; j < m_VolPhaseList.size(); j++) {
delete m_VolPhaseList[j];
m_VolPhaseList[j] = 0;
}
for (size_t j = 0; j < m_speciesThermoList.size(); j++) {
delete m_speciesThermoList[j];
m_speciesThermoList[j] = 0;
}
for (size_t i = 0; i < SpeciesThermo.size(); i++) {
delete SpeciesThermo[i];
SpeciesThermo[i] = 0;
}
for (size_t iph = 0; iph < VPhaseList.size(); iph++) {
delete VPhaseList[iph];
VPhaseList[iph] = 0;
}
delete m_VCount;
m_VCount = 0;
@ -637,7 +587,7 @@ void VCS_SOLVE::vcs_prob_specifyFully()
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
ThermoPhase* tPhase = &m_mix->phase(iphase);
vcs_VolPhase* volPhase = VPhaseList[iphase];
vcs_VolPhase* volPhase = m_VolPhaseList[iphase].get();
volPhase->setState_TP(m_temperature, m_pressurePA);
vector_fp muPhase(tPhase->nSpecies(),0.0);
@ -673,7 +623,7 @@ void VCS_SOLVE::vcs_prob_specifyFully()
plogf(" Initial_Estimated_kMols\n");
for (size_t i = 0; i < m_nsp; i++) {
size_t iphase = m_phaseID[i];
vcs_VolPhase* VolPhase = m_VolPhaseList[iphase];
vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
plogf("%16s %5d %16s", m_speciesName[i].c_str(), iphase,
VolPhase->PhaseName.c_str());
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
@ -690,7 +640,7 @@ void VCS_SOLVE::vcs_prob_specifyFully()
plogf(" TMolesInert Tmoles(kmol)\n");
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* VolPhase = VPhaseList[iphase];
vcs_VolPhase* VolPhase = m_VolPhaseList[iphase].get();
plogf("%16s %5d %5d %8d %16s %8d %16e ", VolPhase->PhaseName.c_str(),
VolPhase->VP_ID_, VolPhase->m_singleSpecies,
VolPhase->m_gasPhase, VolPhase->eos_name(),
@ -720,7 +670,6 @@ void VCS_SOLVE::vcs_prob_specifyFully()
int VCS_SOLVE::vcs_prob_update()
{
size_t k1 = 0;
vcs_tmoles();
m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA,
&m_molNumSpecies_old[0], &m_PMVolumeSpecies[0]);
@ -728,6 +677,7 @@ int VCS_SOLVE::vcs_prob_update()
for (size_t i = 0; i < m_nsp; ++i) {
// Find the index of I in the index vector, m_speciesIndexVector[]. Call
// it K1 and continue.
size_t k1 = 0;
for (size_t j = 0; j < m_nsp; ++j) {
k1 = j;
if (m_speciesMapIndex[j] == i) {
@ -744,44 +694,30 @@ int VCS_SOLVE::vcs_prob_update()
m_gibbsSpecies[i] = m_feSpecies_old[k1];
}
size_t kT = 0;
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* pubPhase = VPhaseList[iph];
vcs_VolPhase* vPhase = m_VolPhaseList[iph];
pubPhase->setTotalMolesInert(vPhase->totalMolesInert());
pubPhase->setTotalMoles(vPhase->totalMoles());
pubPhase->setElectricPotential(vPhase->electricPotential());
double sumMoles = pubPhase->totalMolesInert();
pubPhase->setMoleFractionsState(vPhase->totalMoles(),
&vPhase->moleFractions()[0],
VCS_STATECALC_TMP);
const vector_fp & mfVector = pubPhase->moleFractions();
for (size_t k = 0; k < pubPhase->nSpecies(); k++) {
kT = pubPhase->spGlobalIndexVCS(k);
mf[kT] = mfVector[k];
if (pubPhase->phiVarIndex() == k) {
k1 = vPhase->spGlobalIndexVCS(k);
double tmp = m_molNumSpecies_old[k1];
if (! vcs_doubleEqual(pubPhase->electricPotential() , tmp)) {
vcs_VolPhase* vPhase = m_VolPhaseList[iph].get();
double sumMoles = vPhase->totalMolesInert();
const vector_fp & mfVector = vPhase->moleFractions();
for (size_t k = 0; k < vPhase->nSpecies(); k++) {
size_t kT = vPhase->spGlobalIndexVCS(k);
size_t kOrig = m_speciesMapIndex[kT];
mf[kOrig] = mfVector[k];
if (vPhase->phiVarIndex() == k) {
double tmp = m_molNumSpecies_old[vPhase->spGlobalIndexVCS(k)];
if (!vcs_doubleEqual(vPhase->electricPotential(), tmp)) {
throw CanteraError("VCS_SOLVE::vcs_prob_update",
"We have an inconsistency in voltage, {} {}",
pubPhase->electricPotential(), tmp);
vPhase->electricPotential(), tmp);
}
}
if (! vcs_doubleEqual(mf[kT], vPhase->molefraction(k))) {
throw CanteraError("VCS_SOLVE::vcs_prob_update",
"We have an inconsistency in mole fraction, {} {}",
mf[kT], vPhase->molefraction(k));
}
if (pubPhase->speciesUnknownType(k) != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
sumMoles += w[kT];
if (vPhase->speciesUnknownType(k) != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
sumMoles += w[kOrig];
}
}
if (! vcs_doubleEqual(sumMoles, vPhase->totalMoles())) {
throw CanteraError("VCS_SOLVE::vcs_prob_update",
"We have an inconsistency in total moles, {} {}",
sumMoles, pubPhase->totalMoles());
sumMoles, vPhase->totalMoles());
}
}
return VCS_SUCCESS;
@ -810,7 +746,7 @@ double VCS_SOLVE::vcs_VolTotal(const double tkelvin, const double pres,
{
double VolTot = 0.0;
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
Vphase->setState_TP(tkelvin, pres);
Vphase->setMolesFromVCS(VCS_STATECALC_OLD, w);
double Volp = Vphase->sendToVCS_VolPM(volPM);

30
src/equil/vcs_solve_TP.cpp
View file

@ -86,7 +86,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
plogf("%5d PHASES\n", m_numPhases);
plogf(" PRESSURE%22.8g %3s\n", m_pressurePA, "Pa ");
plogf(" TEMPERATURE%19.3f K\n", m_temperature);
vcs_VolPhase* Vphase = m_VolPhaseList[0];
vcs_VolPhase* Vphase = m_VolPhaseList[0].get();
if (Vphase->nSpecies() > 0) {
plogf(" PHASE1 INERTS%17.3f\n", TPhInertMoles[0]);
}
@ -471,7 +471,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
size_t kspec = m_indexRxnToSpecies[irxn];
double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
size_t iph = m_phaseID[kspec];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
ANOTE[0] = '\0';
double dx;
@ -715,7 +715,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
// We are going to zero the single species phase.
// Set the existence flag
iph = m_phaseID[kspec];
Vphase = m_VolPhaseList[iph];
Vphase = m_VolPhaseList[iph].get();
sprintf(ANOTE, "zeroing out SS phase: ");
// Change the base mole numbers for the iteration.
@ -980,7 +980,7 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
plogf(" --- ");
writeline('-', 50);
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
plogf(" --- %18s = %15.7E\n", Vphase->PhaseName, m_tPhaseMoles_new[iph]);
}
plogf(" ");
@ -1512,7 +1512,7 @@ int VCS_SOLVE::vcs_delete_species(const size_t kspec)
{
const size_t klast = m_numSpeciesRdc - 1;
const size_t iph = m_phaseID[kspec];
vcs_VolPhase* const Vphase = m_VolPhaseList[iph];
vcs_VolPhase* const Vphase = m_VolPhaseList[iph].get();
const size_t irxn = kspec - m_numComponents;
// Zero the concentration of the species.
@ -1588,7 +1588,7 @@ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec)
--m_numRxnMinorZeroed;
}
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
&m_molNumSpecies_old[0],
&m_tPhaseMoles_old[0]);
@ -1622,7 +1622,7 @@ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec)
bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
{
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
bool successful = true;
// set the phase existence flag to dead
@ -1843,7 +1843,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted()
for (int cits = 0; cits < 3; cits++) {
for (size_t kspec = m_numSpeciesRdc; kspec < m_nsp; ++kspec) {
size_t iph = m_phaseID[kspec];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
if (m_molNumSpecies_new[kspec] == 0.0) {
m_molNumSpecies_new[kspec] = VCS_DELETE_MINORSPECIES_CUTOFF * 1.0E-10;
}
@ -2799,7 +2799,7 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc,
// we also trigger an update check for each VolPhase to see if its mole
// numbers are current with vcs
for (size_t iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
Vphase->updateFromVCS_MoleNumbers(stateCalc);
if (!Vphase->m_singleSpecies) {
Vphase->sendToVCS_ActCoeff(stateCalc, &actCoeff_ptr[0]);
@ -2972,7 +2972,7 @@ double VCS_SOLVE::vcs_tmoles()
double sum = 0.0;
for (size_t i = 0; i < m_numPhases; i++) {
sum += m_tPhaseMoles_old[i];
vcs_VolPhase* Vphase = m_VolPhaseList[i];
vcs_VolPhase* Vphase = m_VolPhaseList[i].get();
if (m_tPhaseMoles_old[i] == 0.0) {
Vphase->setTotalMoles(0.0);
} else {
@ -3007,7 +3007,7 @@ void VCS_SOLVE::check_tmoles() const
void VCS_SOLVE::vcs_updateVP(const int vcsState)
{
for (size_t i = 0; i < m_numPhases; i++) {
vcs_VolPhase* Vphase = m_VolPhaseList[i];
vcs_VolPhase* Vphase = m_VolPhaseList[i].get();
if (vcsState == VCS_STATECALC_OLD) {
Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD,
&m_molNumSpecies_old[0],
@ -3226,7 +3226,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted,
if (alterZeroedPhases && false) {
for (size_t iph = 0; iph < m_numPhases; iph++) {
bool lneed = false;
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph].get();
if (! Vphase->m_singleSpecies) {
double sum = 0.0;
for (size_t k = 0; k < Vphase->nSpecies(); k++) {
@ -3333,7 +3333,7 @@ void VCS_SOLVE::vcs_printDeltaG(const int stateCalc)
}
double mfValue = 1.0;
size_t iphase = m_phaseID[kspec];
const vcs_VolPhase* Vphase = m_VolPhaseList[iphase];
const vcs_VolPhase* Vphase = m_VolPhaseList[iphase].get();
if ((m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDMS) ||
(m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) ||
(m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDSS)) {
@ -3408,8 +3408,8 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k
}
// Handle the index pointer in the phase structures first
vcs_VolPhase* pv1 = m_VolPhaseList[m_phaseID[k1]];
vcs_VolPhase* pv2 = m_VolPhaseList[m_phaseID[k2]];
vcs_VolPhase* pv1 = m_VolPhaseList[m_phaseID[k1]].get();
vcs_VolPhase* pv2 = m_VolPhaseList[m_phaseID[k2]].get();
size_t kp1 = m_speciesLocalPhaseIndex[k1];
size_t kp2 = m_speciesLocalPhaseIndex[k2];
AssertThrowMsg(pv1->spGlobalIndexVCS(kp1) == k1, "VCS_SOLVE::vcs_switch_pos",

13
src/equil/vcs_species_thermo.cpp
View file

@ -16,10 +16,9 @@
using namespace std;
namespace Cantera
{
VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(size_t indexPhase,
size_t indexSpeciesPhase) :
IndexPhase(indexPhase),
IndexSpeciesPhase(indexSpeciesPhase),
VCS_SPECIES_THERMO::VCS_SPECIES_THERMO() :
IndexPhase(0),
IndexSpeciesPhase(0),
OwningPhase(0),
SS0_Model(VCS_SS0_CONSTANT),
SS0_feSave(0.0),
@ -33,12 +32,6 @@ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(size_t indexPhase,
SSStar_Vol_Model(VCS_SSVOL_IDEALGAS),
SSStar_Vol0(-1.0)
{
SS0_Pref = 1.01325E5;
}
VCS_SPECIES_THERMO* VCS_SPECIES_THERMO::duplMyselfAsVCS_SPECIES_THERMO()
{
return new VCS_SPECIES_THERMO(*this);
}
}