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

Fixed an error in vcs_solve_TP. m_tPhaseMoles_new wasn't set

Browse source 
before being used.
This commit is contained in:
Harry Moffat 2010-10-28 01:33:54 +00:00
parent c7ce446b32
commit 7fc2841f85
4 changed files with 191 additions and 12 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

20
Cantera/src/equil/vcs_defs.h
View file

@ -146,7 +146,16 @@ namespace VCSnonideal {
* These defines are valid values for spStatus()
*/
//@{
//! Species is a component
//! Species is a component which can never be nonzero because of a
//! stoichiometric constraint
/*!
* An example of this would be a species that contains Ni. But,
* the amount of Ni elements is exactly zero.
*/
#define VCS_SPECIES_COMPONENT_STOICHZERO 3
//! Species is a component which can be nonzero
#define VCS_SPECIES_COMPONENT 2
//! Species is a major species
@ -223,10 +232,11 @@ namespace VCSnonideal {
//! Species lies in a multicomponent phase that is active,
//! but species concentration is zero due to stoich constraint
/*!
* The species lies in a multicomponent phase which
* currently does exist. Its concentration is currently
* identically zero, though the phase exists. This is
* a permament condition due to stoich constraints
* The species lies in a multicomponent phase which currently does exist. Its concentration is currently
* identically zero, though the phase exists. This is a permament condition due to stoich constraints.
*
* An example of this would be a species that contains Ni. But,
* the amount of Ni elements in the current problem statement is exactly zero.
*/
#define VCS_SPECIES_STOICHZERO -8

159
Cantera/src/equil/vcs_phaseStability.cpp
View file

@ -30,8 +30,7 @@ using namespace std;
namespace VCSnonideal {
//====================================================================================================================
// Utility function that evaluates whether a phase can be popped
// into existence
// Utility function that evaluates whether a phase can be popped into existence
/*
* A phase can be popped iff the stoichiometric coefficients for the
* component species, whose concentrations will be lowered during the
@ -141,10 +140,164 @@ namespace VCSnonideal {
}
return false;
}
//====================================================================================================================
int inList(const std::vector<int> &list, int val)
{
for (int i = 0; i < (int) list.size(); i++) {
if (val == list[i]) {
return i;
}
}
return -1;
}
//====================================================================================================================
// Determine the list of problems that need to be checked to see if there are any phases pops
/*
* This routine evaluates and fills in the following quantities
* phasePopProblemLists_
*
* Need to work in species that are zeroed by element constraints
*
* @return Returns the number of problems that must be checked.
*/
int VCS_SOLVE::vcs_phasePopDeterminePossibleList() {
int nfound = 0;
int irxn, kspec;
vcs_VolPhase *Vphase = 0;
int iph, j, k;
int nsp;
double stoicC;
double molComp;
std::vector<int> linkedPhases;
phasePopProblemLists_.clear();
/*
* This is a vector over each component.
* For zeroed components it lists the phases, which are currently zeroed,
* which have a species with a positive stoichiometric value wrt the component.
* Therefore, we could pop the component species and pop that phase at the same time
* if we considered no other factors than keeping the component mole number positve.
*
* It does not count species with positive stoichiometric values if that species
* already has a positive mole number. The phase is already popped.
*/
std::vector< std::vector<int> > zeroedComponentLinkedPhasePops(m_numComponents);
/*
* The logic below calculates zeroedComponentLinkedPhasePops
*/
for (j = 0; j < m_numComponents; j++) {
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
molComp = m_molNumSpecies_old[j];
if (molComp <= 0.0) {
std::vector<int> &jList = zeroedComponentLinkedPhasePops[j];
iph = m_phaseID[j];
jList.push_back(iph);
for (irxn = 0; irxn < m_numRxnTot; irxn++) {
kspec = irxn + m_numComponents;
iph = m_phaseID[kspec];
Vphase = m_VolPhaseList[iph];
int existence = Vphase->exists();
if (existence < 0) {
stoicC = m_stoichCoeffRxnMatrix[irxn][j];
if (stoicC > 0.0) {
if (inList(jList, iph) != -1) {
jList.push_back(iph);
}
}
}
}
}
}
}
/*
* This is a vector over each zeroed phase
* For zeroed phases, it lists the components, which are currently zereoed,
* which have a species with a negative stoichiometric value wrt one or more species in the phase.
* Cut out components which have a pos stoichiometric value with another species in the phase.
*/
std::vector< std::vector<int> > zeroedPhaseLinkedZeroComponents(m_numPhases);
/*
* The logic below calculates zeroedPhaseLinkedZeroComponents
*/
for (iph = 0; iph < m_numPhases; iph++) {
std::vector<int> &iphList = zeroedPhaseLinkedZeroComponents[iph];
iphList.clear();
Vphase = m_VolPhaseList[iph];
int existence = Vphase->exists();
if (existence < 0) {
linkedPhases.clear();
nsp = Vphase->nSpecies();
for (k = 0; k < nsp; k++) {
kspec = Vphase->spGlobalIndexVCS(k);
irxn = kspec - m_numComponents;
for (j = 0; j < m_numComponents; j++) {
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
molComp = m_molNumSpecies_old[j];
if (molComp <= 0.0) {
stoicC = m_stoichCoeffRxnMatrix[irxn][j];
if (stoicC < 0.0) {
bool foundPos = false;
for (int kk = 0; kk < nsp; kk++) {
int kkspec = Vphase->spGlobalIndexVCS(kk);
int iirxn = kkspec - m_numComponents;
if (iirxn >= 0) {
if (m_stoichCoeffRxnMatrix[iirxn][j] > 0.0) {
foundPos = true;
}
}
}
if (!foundPos) {
if (inList(iphList, j) != -1) {
iphList.push_back(j);
}
}
}
}
}
}
}
}
}
/*
* Now fill in the phasePopProblemLists_ list.
*
*/
for (iph = 0; iph < m_numPhases; iph++) {
Vphase = m_VolPhaseList[iph];
int existence = Vphase->exists();
if (existence < 0) {
std::vector<int> &iphList = zeroedPhaseLinkedZeroComponents[iph];
std::vector<int> popProblem(0);
popProblem.push_back(iph);
for (int i = 0; i < (int) iphList.size(); i++) {
j = iphList[i];
std::vector<int> &jList = zeroedComponentLinkedPhasePops[j];
for (int jjl = 0; jjl < (int) jList.size(); jjl++) {
int jph = jList[jjl];
if (inList(popProblem, jph) != -1) {
popProblem.push_back(jph);
}
}
}
phasePopProblemLists_.push_back(popProblem);
}
}
return nfound;
}
//====================================================================================================================
// Decision as to whether a phase pops back into existence
/*
* @return returns the phase id of the phase that pops back into
* @return returns the phase id of the phases that pops back into
* existence. Returns -1 if there are no phases
*/
int VCS_SOLVE::vcs_popPhaseID(std::vector<int> & phasePopPhaseIDs) {

17
Cantera/src/equil/vcs_solve.h
View file

@ -512,6 +512,21 @@ public:
*/
bool vcs_popPhasePossible(const int iphasePop) const;
//! Determine the list of problems that need to be checked to see if there are any phases pops
/*!
* This routine evaluates and fills in the following quantities
* phasePopProblemLists_
*
* @return Returns the number of problems that must be checked.
*/
int vcs_phasePopDeterminePossibleList();
//! Decision as to whether a phase pops back into existence
/*!
* @param phasePopPhaseIDs Vector containing the phase ids of the phases
@ -1960,6 +1975,8 @@ public:
*/
std::vector<double> m_chargeSpecies;
std::vector<std::vector<int> > phasePopProblemLists_;
//! Vector of pointers to thermostructures which identify the model
//! and parameters for evaluating the thermodynamic functions for that
//! particular species.

7
Cantera/src/equil/vcs_solve_TP.cpp
View file

@ -387,12 +387,11 @@ namespace VCSnonideal {
/*
* Copy the old solution into the new solution as an initial guess
*/
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new),
VCS_DATA_PTR(m_feSpecies_old), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_new),
VCS_DATA_PTR(m_actCoeffSpecies_old), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_feSpecies_old), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_new), VCS_DATA_PTR(m_actCoeffSpecies_old), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_new), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc);
vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc);
vcs_dcopy(VCS_DATA_PTR(m_tPhaseMoles_new), VCS_DATA_PTR(m_tPhaseMoles_old), m_numPhases);
/*
* Zero out the entire vector of updates. We sometimes would