From ea5d8387353c82f47d4cd6e86d50b4d12e4c24c4 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 23 Jul 2010 17:51:03 +0000 Subject: [PATCH] Added a bit more logic into the phaseStability calc. More to come on this. --- Cantera/src/equil/equil.h | 3 +- Cantera/src/equil/vcs_MultiPhaseEquil.h | 18 +++---- Cantera/src/equil/vcs_defs.h | 18 ++++++- Cantera/src/equil/vcs_phaseStability.cpp | 69 +++++++++++++++++++++--- Cantera/src/equil/vcs_solve.h | 25 +++++---- 5 files changed, 102 insertions(+), 31 deletions(-) diff --git a/Cantera/src/equil/equil.h b/Cantera/src/equil/equil.h index 937b70a4a..c8d4498a8 100644 --- a/Cantera/src/equil/equil.h +++ b/Cantera/src/equil/equil.h @@ -19,6 +19,7 @@ //#include "ChemEquil.h" #include "MultiPhase.h" +#include "vcs_defs.h" namespace Cantera { @@ -72,7 +73,7 @@ namespace Cantera { * @ingroup equil */ int equilibrate(thermo_t& s, const char* XY, - int solver = -1, doublereal rtol = 1.0e-9, int maxsteps = 5000, + int solver = -1, doublereal rtol = 1.0e-9, int maxsteps = VCS_MAXSTEPS, int maxiter = 100, int loglevel = -99); //! Equilibrate a MultiPhase object diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.h b/Cantera/src/equil/vcs_MultiPhaseEquil.h index 3b35973a4..9abff4bb5 100644 --- a/Cantera/src/equil/vcs_MultiPhaseEquil.h +++ b/Cantera/src/equil/vcs_MultiPhaseEquil.h @@ -15,7 +15,7 @@ #include "ct_defs.h" #include "MultiPhase.h" - +#include "vcs_defs.h" namespace Cantera { @@ -81,7 +81,7 @@ namespace Cantera { int vcs_equilibrate(thermo_t& s, const char* XY, int estimateEquil = 0, int printLvl = 0, int solver = -1, doublereal rtol = 1.0e-9, - int maxsteps = 5000, + int maxsteps = VCS_MAXSTEPS, int maxiter = 100, int loglevel = -99); @@ -142,7 +142,7 @@ namespace Cantera { int vcs_equilibrate(MultiPhase& s, const char* XY, int estimateEquil = 0, int printLvl = 0, int solver = 2, - doublereal rtol = 1.0e-9, int maxsteps = 5000, + doublereal rtol = 1.0e-9, int maxsteps = VCS_MAXSTEPS, int maxiter = 100, int loglevel = -99); //! Set a multi-phase chemical solution to chemical equilibrium. @@ -202,7 +202,7 @@ namespace Cantera { int vcs_equilibrate_1(MultiPhase& s, int ixy, int estimateEquil = 0, int printLvl = 0, int solver = 2, - doublereal rtol = 1.0e-9, int maxsteps = 5000, + doublereal rtol = 1.0e-9, int maxsteps = VCS_MAXSTEPS, int maxiter = 100, int loglevel = -99); } @@ -355,7 +355,7 @@ namespace VCSnonideal { */ int equilibrate(int XY, int estimateEquil = 0, int printLvl= 0, doublereal err = 1.0e-6, - int maxsteps = 5000, int loglevel=-99); + int maxsteps = VCS_MAXSTEPS, int loglevel=-99); //! Equilibrate the solution using the current element abundances //! storred in the MultiPhase object using constant T and P @@ -385,7 +385,7 @@ namespace VCSnonideal { */ int equilibrate_TP(int estimateEquil = 0, int printLvl= 0, doublereal err = 1.0e-6, - int maxsteps = 5000, int loglevel=-99); + int maxsteps = VCS_MAXSTEPS, int loglevel=-99); //! Equilibrate the solution using the current element abundances //! storred in the MultiPhase object using either constant H and P @@ -439,7 +439,7 @@ namespace VCSnonideal { int equilibrate_HP(doublereal Htarget, int XY, double Tlow, double Thigh, int estimateEquil = 0, int printLvl = 0, doublereal err = 1.0E-6, - int maxsteps = 5000, int loglevel=-99); + int maxsteps = VCS_MAXSTEPS, int loglevel=-99); //! Equilibrate the solution using the current element abundances //! storred in the MultiPhase object using constant S and P. @@ -490,7 +490,7 @@ namespace VCSnonideal { int equilibrate_SP(doublereal Starget, double Tlow, double Thigh, int estimateEquil = 0, int printLvl = 0, doublereal err = 1.0E-6, - int maxsteps = 5000, int loglevel=-99); + int maxsteps = VCS_MAXSTEPS, int loglevel=-99); //! Equilibrate the solution using the current element abundances @@ -538,7 +538,7 @@ namespace VCSnonideal { int equilibrate_TV(int XY, doublereal xtarget, int estimateEquil = 0, int printLvl = 0, doublereal err = 1.0E-6, - int maxsteps = 5000, int loglevel = -99); + int maxsteps = VCS_MAXSTEPS, int loglevel = -99); //! Report the equilibrium answer in a comma separated table format /*! diff --git a/Cantera/src/equil/vcs_defs.h b/Cantera/src/equil/vcs_defs.h index 056859ff6..2ee2b26ca 100644 --- a/Cantera/src/equil/vcs_defs.h +++ b/Cantera/src/equil/vcs_defs.h @@ -78,7 +78,8 @@ namespace VCSnonideal { /*! * @name Sizes of Phases and Cutoff Mole Numbers - * + * + * All size parameters are listed here * @{ */ @@ -103,14 +104,27 @@ namespace VCSnonideal { //! Cutoff relative moles below which a phase is deleted //! from the equilibrium problem. #ifndef VCS_DELETE_PHASE_CUTOFF -#define VCS_DELETE_PHASE_CUTOFF 1.0e-12 +#define VCS_DELETE_PHASE_CUTOFF 1.0e-13 #endif + //! Relative mole number of species in a phase that is created + //! We want this to be comfortably larger than the VCS_DELETE_PHASE_CUTOFF value + //! so that the phase can have a chance to survive. +#ifndef VCS_POP_PHASE_MOLENUM +#define VCS_POP_PHASE_MOLENUM 1.0e-11 +#endif + + //! Cutoff moles below which a phase or species which //! comprises the bulk of an element's total concentration //! is deleted. #ifndef VCS_DELETE_ELEMENTABS_CUTOFF #define VCS_DELETE_ELEMENTABS_CUTOFF 1.0e-280 +#endif + + //! Maximum steps in the inner loop +#ifndef VCS_MAXSTEPS +#define VCS_MAXSTEPS 50000 #endif //@} diff --git a/Cantera/src/equil/vcs_phaseStability.cpp b/Cantera/src/equil/vcs_phaseStability.cpp index a0a560abf..e92dd7f54 100644 --- a/Cantera/src/equil/vcs_phaseStability.cpp +++ b/Cantera/src/equil/vcs_phaseStability.cpp @@ -20,13 +20,17 @@ using namespace std; namespace VCSnonideal { - + //==================================================================================================================== // 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 * process, are positive by at least a small degree. + * + * If one of the phase species is a zeroed component, then the phase can + * be popped if the component increases in mole number as the phase moles + * are increased. * * @param iphasePop id of the phase, which is currently zeroed, * @@ -52,6 +56,13 @@ namespace VCSnonideal { */ for (int k = 0; k < Vphase->nSpecies(); k++) { int kspec = Vphase->spGlobalIndexVCS(k); +#ifdef DEBUG_MODE + if (m_molNumSpecies_old[kspec] > 0.0) { + printf("ERROR vcs_popPhasePossible we shouldn't be here %d %g > 0.0", + kspec, m_molNumSpecies_old[kspec]); + exit(-1); + } +#endif int irxn = kspec - m_numComponents; if (irxn >= 0) { int iPopPossible = true; @@ -62,7 +73,7 @@ namespace VCSnonideal { double negChangeComp = - stoicC * 1.0; if (negChangeComp > 0.0) { // TODO: We may have to come up with a tolerance here - if (m_molNumSpecies_old[j] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.1) { + if (m_molNumSpecies_old[j] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { iPopPossible = false; } } @@ -72,11 +83,56 @@ namespace VCSnonideal { if (iPopPossible == true) { return true; } + } else { + /* + * We are here when the species in the phase is a component. Its mole number is zero. + * We loop through the regular reaction looking for a reaction that can pop the + * component. + */ + printf("WE are here at new logic - CHECK\n"); + for (int jrxn = 0; jrxn < m_numRxnRdc; jrxn++) { + bool foundJrxn = false; + // First, if the component is a product of the reaction + if (m_stoichCoeffRxnMatrix[jrxn][kspec] > 0.0) { + foundJrxn = true; + for (int kcomp = 0; kcomp < m_numComponents; kcomp++) { + if (m_stoichCoeffRxnMatrix[jrxn][kcomp] < 0.0) { + if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { + foundJrxn = false; + } + } + } + if (foundJrxn) { + printf("We have found a component phase pop! CHECK1 \n"); + return true; + } + } + // Second we are here if the component is a reactant in the reaction, and the reaction goes backwards. + else if (m_stoichCoeffRxnMatrix[jrxn][kspec] < 0.0) { + foundJrxn = true; + int jspec = jrxn + m_numComponents; + if (m_molNumSpecies_old[jspec] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { + foundJrxn = false; + continue; + } + for (int kcomp = 0; kcomp < m_numComponents; kcomp++) { + if (m_stoichCoeffRxnMatrix[jrxn][kcomp] > 0.0) { + if (m_molNumSpecies_old[kcomp] <= VCS_DELETE_ELEMENTABS_CUTOFF*0.5) { + foundJrxn = false; + } + } + } + if (foundJrxn) { + printf("We have found a component phase pop! CHECK2 \n"); + return true; + } + } + } } } return false; } - + //==================================================================================================================== // Decision as to whether a phase pops back into existence /* * @return returns the phase id of the phase that pops back into @@ -195,7 +251,7 @@ namespace VCSnonideal { #endif return iphasePop; } - + //==================================================================================================================== // Calculates the deltas of the reactions due to phases popping // into existence /* @@ -415,7 +471,7 @@ namespace VCSnonideal { // - + //==================================================================================================================== double VCS_SOLVE::vcs_phaseStabilityTest(const int iph) { /* @@ -719,6 +775,7 @@ namespace VCSnonideal { #endif return funcPhaseStability; } - + //==================================================================================================================== } +//====================================================================================================================== diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 64acf7ee0..053e681c8 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -200,7 +200,7 @@ public: * m_stoichCoeffRxnMatrix[irxn][jcomp] * Stoichiometric coefficient matrix for the reaction mechanism * expressed in Reduced Canonical Form. - * j refers to the component number, and irxn + * jcomp refers to the component number, and irxn * refers to the irxn_th non-component species. * * m_deltaMolNumPhase[irxn] @@ -1438,7 +1438,7 @@ public: //! Number of components calculated for the problem int m_numComponents; - //! Total number of non-component species in the problem + //! Total number of non-component species in the problem int m_numRxnTot; //! Current number of species in the problems @@ -1448,7 +1448,7 @@ public: */ int m_numSpeciesRdc; - //! Current number of non-component species in the problem + //! Current number of non-component species in the problem /*! * Species can be deleted if they aren't * stable under the current conditions @@ -1470,20 +1470,19 @@ public: */ DoubleStarStar m_formulaMatrix; - //! Stoichiometric coefficient matrix for the reaction mechanism - //! expressed in Reduced Canonical Form. + //! Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form. /*! * This is the stoichiometric coefficient matrix for the - * reaction which forms species K from the component species. A - * stoichiometric coefficient of one is assumed for the - * species K in this mechanism. + * reaction which forms species kspec from the component species. A + * stoichiometric coefficient of one is assumed for the species kspec in this mechanism. * - * NOTE: kspec = Irxn + m_numComponents + * NOTE: kspec = irxn + m_numComponents * - * sc[irxn][j] : - * j refers to the component number, and irxn - * refers to the irxn_th non-component species. - * + * m_stoichCoeffRxnMatrix[irxn][j] : + * j refers to the component number, and irxn refers to the irxn_th non-component species. + * The stoichiometric coefficents multilpled by the Formula coefficients of the + * component species add up to the negative value of the number of elements in + * the species kspec. * * length = [nspecies0][nelements0] */