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

Added a bit more logic into the phaseStability calc.

Browse source 
More to come on this.
This commit is contained in:
Harry Moffat 2010-07-23 17:51:03 +00:00
parent 96be262030
commit ea5d838735
5 changed files with 102 additions and 31 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
Cantera/src/equil/equil.h
View file

@ -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

18
Cantera/src/equil/vcs_MultiPhaseEquil.h
View file

@ -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
/*!

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

@ -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
//@}

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

@ -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;
}
//====================================================================================================================
}
//======================================================================================================================

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

@ -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]
*/