From 7fc2841f85f0eca39649a2a85594a24c7135066c Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 28 Oct 2010 01:33:54 +0000 Subject: [PATCH] Fixed an error in vcs_solve_TP. m_tPhaseMoles_new wasn't set before being used. --- Cantera/src/equil/vcs_defs.h | 20 ++- Cantera/src/equil/vcs_phaseStability.cpp | 159 ++++++++++++++++++++++- Cantera/src/equil/vcs_solve.h | 17 +++ Cantera/src/equil/vcs_solve_TP.cpp | 7 +- 4 files changed, 191 insertions(+), 12 deletions(-) diff --git a/Cantera/src/equil/vcs_defs.h b/Cantera/src/equil/vcs_defs.h index 2ee2b26ca..534d5ca54 100644 --- a/Cantera/src/equil/vcs_defs.h +++ b/Cantera/src/equil/vcs_defs.h @@ -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 diff --git a/Cantera/src/equil/vcs_phaseStability.cpp b/Cantera/src/equil/vcs_phaseStability.cpp index 17cbf9a2c..997bbf830 100644 --- a/Cantera/src/equil/vcs_phaseStability.cpp +++ b/Cantera/src/equil/vcs_phaseStability.cpp @@ -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 &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 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 > 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 &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 > zeroedPhaseLinkedZeroComponents(m_numPhases); + /* + * The logic below calculates zeroedPhaseLinkedZeroComponents + */ + for (iph = 0; iph < m_numPhases; iph++) { + std::vector &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 &iphList = zeroedPhaseLinkedZeroComponents[iph]; + std::vector popProblem(0); + popProblem.push_back(iph); + for (int i = 0; i < (int) iphList.size(); i++) { + j = iphList[i]; + std::vector &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 & phasePopPhaseIDs) { diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 2b43c55a6..6da5b2358 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -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 m_chargeSpecies; + std::vector > phasePopProblemLists_; + //! Vector of pointers to thermostructures which identify the model //! and parameters for evaluating the thermodynamic functions for that //! particular species. diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 556cd0b93..2772202f7 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -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