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

Algorithm improvements in the equilibrium solver.

Browse source
This commit is contained in:
Harry Moffat 2009-04-18 00:42:47 +00:00
parent 973d12f864
commit a5103167dd
16 changed files with 7157 additions and 8202 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cantera/python/setup.py.in
View file

@ -13,12 +13,16 @@ if @build_with_f2c@ <> 1:
linkargs = '@LCXX_FLAGS@'
numarray_incl = "@NUMARRAY_INC_DIR@"
numpy_incl = "@NUMPY_INC_DIR@"
incdirs=["../../build/include", "src", "../clib/src"]
if numarray_incl <> '':
incdirs.append(numarray_incl)
if numpy_incl <> '':
incdirs.append(numpy_incl)
bllibstr = "@BLAS_LAPACK_LIBS@"
bllibs = bllibstr.replace('-l',' ')
bllist = bllibs.split()

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

@ -450,6 +450,23 @@ namespace VCSnonideal {
vcs_SpeciesProperties * speciesProperty(const int kindex);
//! int indicating whether the phase exists or not
/*!
* returns the m_existence int for the phase
*
* - VCS_PHASE_EXIST_ZEROEDPHASE = -6: Set to not exist by fiat from a
* higher level.
* This is used in phase stability boundary calculations
* - VCS_PHASE_EXIST_NO = 0: Doesn't exist currently
* - VCS_PHASE_EXIST_MINORCONC = 1: Exists, but the concentration is
* so low that an alternate
* method is used to calculate the total phase concentrations.
* - VCS_PHASE_EXIST_YES = 2 : Does exist currently
* - VCS_PHASE_EXIST_ALWAYS = 3: Always exists because it contains
* inerts which can't exist in any other phase. Or,
* the phase exists always because it consists of a single
* species, which is identified with the voltage, i.e.,
* its an electron metal phase.
*/
int exists() const;
//! Set the existence flag in the object

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

@ -206,6 +206,16 @@ namespace VCSnonideal {
*/
#define VCS_SPECIES_ACTIVEBUTZERO -7
//! 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
*/
#define VCS_SPECIES_STOICHZERO -8
//@}
//! @name Phase Categories used during the iteration

2
Cantera/src/equil/vcs_internal.h
View file

@ -473,6 +473,8 @@ namespace VCSnonideal {
*/
void vcs_print_line(const char *str, int num);
const char *vcs_speciesType_string(int speciesStatus, int length = 100);
//! Print a string within a given space limit
/*!
* This routine limits the amount of the string that will be printed to a

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

@ -1,9 +1,5 @@
/* ======================================================================= */
/* -------------------------------------------------- */
/* | RCS Head Information on zuzax.pchem.sandia.gov | */
/* -------------------------------------------------- */
/* $RCSfile$ */
/* $Author$ */
/* $Date$ */
/* $Revision$ */
/* ======================================================================= */
@ -50,12 +46,15 @@ namespace VCSnonideal {
#endif
/*
* Loop through all of the species in the phase
* Loop through all of the species in the phase. We say the phase
* can be popped, if there is one species in the phase that can be
* popped.
*/
for (int k = 0; k < Vphase->nSpecies(); k++) {
int kspec = Vphase->spGlobalIndexVCS(k);
int irxn = kspec - m_numComponents;
if (irxn >= 0) {
int iPopPossible = true;
for (int j = 0; j < m_numComponents; ++j) {
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
@ -64,24 +63,18 @@ namespace VCSnonideal {
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) {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 3) {
plogf(" --- vcs_popPhasePosssible() Phase %d (%s) can't be popped\n",
iphasePop,
Vphase->PhaseName.c_str());
plogf(" --- Component %d (%s)will go negative\n",
j, m_speciesName[j].c_str());
}
#endif
return false;
iPopPossible = false;
}
}
}
}
}
if (iPopPossible == true) {
return true;
}
}
}
return true;
return false;
}
// Decision as to whether a phase pops back into existence
@ -604,7 +597,7 @@ namespace VCSnonideal {
irxn = kspec - m_numComponents;
double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
if (irxn >= 0) {
fracDelta_raw[k] = (sumFrac - fracDelta_old[k]) * b / (1.0 - b);
fracDelta_raw[k] = b;
}
}

25
Cantera/src/equil/vcs_rxnadj.cpp
View file

@ -105,14 +105,27 @@ namespace VCSnonideal {
double trphmoles = tphmoles / m_totalMolNum;
if (trphmoles > VCS_DELETE_PHASE_CUTOFF) {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
#ifdef DEBUG_MODE
if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
m_deltaMolNumSpecies[kspec] = 0.0;
#ifdef DEBUG_MODE
sprintf(ANOTE,
"MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
#endif
} else {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
#ifdef DEBUG_MODE
sprintf(ANOTE,
"MultSpec: small species born again DG = %11.3E",
"MultSpec (%s): small species born again DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
#endif
}
} else {
#ifdef DEBUG_MODE
sprintf(ANOTE, "MultSpec: phase come alive DG = %11.3E",
sprintf(ANOTE, "MultSpec (%s): phase come alive DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
#endif
Vphase = m_VolPhaseList[iph];
@ -120,10 +133,12 @@ namespace VCSnonideal {
m_deltaMolNumSpecies[kspec] =
m_totalMolNum * 10.0 * VCS_DELETE_PHASE_CUTOFF / numSpPhase;
}
--(m_numRxnMinorZeroed);
} else {
#ifdef DEBUG_MODE
sprintf(ANOTE, "MultSpec: still dead DG = %11.3E", m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
#endif
m_deltaMolNumSpecies[kspec] = 0.0;
}

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

@ -1455,8 +1455,8 @@ public:
*/
int m_numRxnRdc;
//! Number of active species which are currently either zeroed out or
//! are minor species
//! Number of active species which are currently either treated as
//! minor species
int m_numRxnMinorZeroed;
//! Number of Phases in the problem
@ -1532,9 +1532,10 @@ public:
//! Setting for whether to do an initial estimate
/*!
* Initial estimate: 0 Do not estimate the solution at all. Use the supplied
* mole numbers as is.
* 1 Only do an estimate if the element abundances aren't satisfied.
* Initial estimate: 0 Do not estimate the solution at all. Use the
* supplied mole numbers as is.
* 1 Only do an estimate if the element abundances
* aren't satisfied.
* -1 Force an estimate of the soln. Throw out the input
* mole numbers.
*/
@ -1814,6 +1815,11 @@ public:
* of an equilibrium problem.
* The species is soon "birthed" or "deleted".
* - VCS_SPECIES_ACTIVEBUTZERO
* -8 -> 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
* - VCS_SPECIES_STOICHZERO
*
*/
std::vector<int> m_speciesStatus;

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

@ -418,44 +418,16 @@ namespace VCSnonideal {
}
if (iphasePop < 0) {
/*
* Figure out whether we will calculate new reaction step sizes
* for the major species.
* -> We won't if all species are minors (im), OR
* all major species have already converged
* Figure out the new reaction step sizes
* for the major species (do minor species in the future too)
*/
if (!(MajorSpeciesHaveConverged) && ! allMinorZeroedSpecies) {
soldel = vcs_RxnStepSizes();
/* - If SOLDEL is true then we encountered a reaction between */
/* - single-species-phase species, only, and have adjusted */
/* - the mole number vector, W(), directly. In this case, */
/* - we should immediately go back and recompute a new */
/* - component basis, if the species that was zeroed was */
/* - a component. SOLDEL is true when this is so. */
if (soldel > 0) {
/* - We have changed the base mole number amongst single- */
/* - species-phase species. However, we don't need to */
/* - recaculate their chemical potentials because they */
/* - are constant, anyway! */
if (soldel == 2) {
goto L_COMPONENT_CALC;
}
/* - We have not changed the actual DG values for */
/* - any species, even the one we deleted. Thus, */
/* - we don't need to start over. */
}
} else {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
if (allMinorZeroedSpecies) {
plogf(" --- vcs_RxnStepSizes not called because all"
"species are minors\n");
} else {
plogf(" --- vcs_RxnStepSizes not called because "
"all majors have converged\n");
}
}
#endif
soldel = vcs_RxnStepSizes();
if (soldel == 2) {
goto L_COMPONENT_CALC;
}
}
#ifdef DEBUG_MODE
else {
@ -552,7 +524,7 @@ namespace VCSnonideal {
/********************************************************************/
/********************** ZEROED OUT SPECIES **************************/
/********************************************************************/
bool resurrect = true;
bool resurrect = (m_deltaMolNumSpecies[kspec] > 0.0);
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 3) {
plogf(" --- %s currently zeroed (SpStatus=%-2d):",
@ -562,23 +534,22 @@ namespace VCSnonideal {
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]);
}
#endif
// HKM Alternative is to not allow ds[] = 0.0 phases
// to pop back into existence. For esthetics, I'm allowing this.
// so that dg < 0.0 phases with zero mole numbers become components.
// This is also better, because that component will be the first
// one to pop into existence if there is a minute quantity of the element.
// This could change in the future.
//if (dg[irxn] >= 0.0 || ds[kspec] <= 0.0) {
if (m_deltaGRxn_new[irxn] >= 0.0 ) {
if (m_deltaGRxn_new[irxn] >= 0.0 || !resurrect) {
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
m_deltaMolNumSpecies[kspec] = 0.0;
resurrect = false;
#ifdef DEBUG_MODE
sprintf(ANOTE, "Species stays zeroed: DG = %11.4E",
m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]);
if (m_deltaGRxn_new[irxn] < 0.0) {
sprintf(ANOTE, "Species stays zeroed even though dg neg:DG = %11.4E, ds zeroed ",
m_deltaGRxn_new[irxn]);
if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
sprintf(ANOTE, "Species stays zeroed even though dg neg due to "
"STOICH/PHASEPOP constraint: DG = %11.4E",
m_deltaGRxn_new[irxn]);
} else {
sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed",
m_deltaGRxn_new[irxn]);
}
}
#endif
} else {
@ -610,7 +581,6 @@ namespace VCSnonideal {
//Vphase->setExistence(1);
phaseResurrected = true;
}
--m_numRxnMinorZeroed;
if (phaseResurrected) {
#ifdef DEBUG_MODE
@ -1281,28 +1251,6 @@ namespace VCSnonideal {
for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
if (m_phaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) {
irxn = kspec - m_numComponents;
// Both of these conditions are false and should be discarded.
// I think a proper special case would be if a species in a small phase
// had a significant contribution to total element total of a single
// element constraint. That's it.
// if (kspec < m_numComponents) {
// if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF) {
//soldel = 0;
//break;
//}
//} else {
//for (k = 0; k < m_numComponents; k++) {
//if (m_stoichCoeffRxnMatrix[irxn][k] != 0.0) {
// if (m_molNumSpecies_old[kspec]/m_molNumSpecies_old[k] >
//VCS_DELETE_PHASE_CUTOFF) {
// soldel = 0;
// break;
// }
//}
//}
//}
}
}
if (soldel) {
@ -1559,8 +1507,8 @@ namespace VCSnonideal {
m_speciesName[kspec].c_str());
}
}
#endif
++m_numRxnMinorZeroed;
#endif
++m_numRxnMinorZeroed;
} else if (speciesType == VCS_SPECIES_MINOR) {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
@ -3676,9 +3624,7 @@ namespace VCSnonideal {
// ---------- Treat special cases first ---------------------
if (kspec < m_numComponents) {
return VCS_SPECIES_COMPONENT;
}
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
return VCS_SPECIES_INTERFACIALVOLTAGE;
}
@ -3691,20 +3637,39 @@ namespace VCSnonideal {
// ---------- Treat zeroed out species first ----------------
if (m_molNumSpecies_old[kspec] <= 0.0) {
if (m_deltaGRxn_old[irxn] >= 0.0) {
/*
* We are here when the species is or should be zeroed out
*/
if (m_SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDSS;
} else {
if (phaseExist >= VCS_PHASE_EXIST_YES) {
return VCS_SPECIES_ACTIVEBUTZERO;
} else if (phaseExist == VCS_PHASE_EXIST_ZEROEDPHASE) {
return VCS_SPECIES_ZEROEDPHASE;
} else {
return VCS_SPECIES_ZEROEDMS;
}
if (m_tPhaseMoles_old[iph] <= 0.0) {
if (!m_SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDMS;
}
}
/*
* see if the species has an element
* which is so low that species will always be zero
*
*/
for (int j = 0; j < m_numElemConstraints; ++j) {
int elType = m_elType[j];
if (elType == VCS_ELEM_TYPE_ABSPOS) {
double atomComp = m_formulaMatrix[j][kspec];
if (atomComp > 0.0) {
double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- %s can not be nonzero because"
" needed element %s is zero\n",
m_speciesName[kspec].c_str(), (m_elementName[j]).c_str());
}
#endif
if (m_SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDSS;
} else {
return VCS_SPECIES_STOICHZERO;
}
}
}
}
}
@ -3717,68 +3682,76 @@ namespace VCSnonideal {
* is also zeroed out. Then, don't pop the phase or the species back into
* existence.
*/
for (int j = 0; j < m_numComponents; ++j) {
double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
if (stoicC != 0.0) {
double negChangeComp = - stoicC;
if (negChangeComp > 0.0) {
if (m_molNumSpecies_old[j] < 1.0E-60) {
if (irxn >= 0) {
for (int j = 0; j < m_numComponents; ++j) {
double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
if (stoicC != 0.0) {
double negChangeComp = - stoicC;
if (negChangeComp > 0.0) {
if (m_molNumSpecies_old[j] < 1.0E-60) {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- %s would have popped back into existance but"
" needed component %s is zero\n",
m_speciesName[kspec].c_str(), m_speciesName[j].c_str());
}
if (m_debug_print_lvl >= 2) {
plogf(" --- %s is prevented from popping into existance because"
" a needed component to be consumed, %s, has a zero mole number\n",
m_speciesName[kspec].c_str(), m_speciesName[j].c_str());
}
#endif
if (m_SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDSS;
} else {
return VCS_SPECIES_ACTIVEBUTZERO;
if (m_SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDSS;
} else {
return VCS_SPECIES_STOICHZERO;
}
}
} else if (negChangeComp < 0.0) {
int jph = m_phaseID[j];
vcs_VolPhase *jVPhase = m_VolPhaseList[jph];
if (jVPhase->exists() <= 0) {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- %s is prevented from popping into existence because"
" a needed component %s is in a zeroed-phase that would be "
"popped into existence at the same time\n",
m_speciesName[kspec].c_str(), m_speciesName[j].c_str());
}
#endif
if (m_SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDSS;
} else {
return VCS_SPECIES_STOICHZERO;
}
}
}
}
}
}
/*
* The Gibbs free energy for this species is such that
* it will pop back into existence.
*
* -> Set it to a major species in anticipation.
* -> An exception to this is if the species has an element
* which is so low to cause problems.
*
* We need to have a PHASE_CUTOFF here. This algorithm is
* insufficient.
*/
for (int j = 0; j < m_numElemConstraints; ++j) {
int elType = m_elType[j];
if (elType == VCS_ELEM_TYPE_ABSPOS) {
double atomComp = m_formulaMatrix[j][kspec];
if (atomComp > 0.0) {
double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- %s would have popped back into existance but"
" needed element %s is zero\n",
m_speciesName[kspec].c_str(), (m_elementName[j]).c_str());
}
#endif
if (m_SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDSS;
} else {
return VCS_SPECIES_ACTIVEBUTZERO;
}
if (irxn >= 0) {
if (m_deltaGRxn_old[irxn] >= 0.0) {
/*
* We are here when the species is or should remain zeroed out
*/
if (m_SSPhase[kspec]) {
return VCS_SPECIES_ZEROEDSS;
} else {
if (phaseExist >= VCS_PHASE_EXIST_YES) {
return VCS_SPECIES_ACTIVEBUTZERO;
} else if (phaseExist == VCS_PHASE_EXIST_ZEROEDPHASE) {
return VCS_SPECIES_ZEROEDPHASE;
} else {
return VCS_SPECIES_ZEROEDMS;
}
}
}
}
/*
* If the current phase already exists, set the emerging species to a
* minor species.
* If the current phase already exists,
*/
if (m_tPhaseMoles_old[iph] > 0.0) {
return VCS_SPECIES_MINOR;
if (m_SSPhase[kspec]) {
return VCS_SPECIES_MAJOR;
} else {
return VCS_SPECIES_ACTIVEBUTZERO;
}
}
/*
@ -3790,7 +3763,13 @@ namespace VCSnonideal {
* fraction of the species in the phase, we could do
* better here.
*/
return VCS_SPECIES_MAJOR;
if (m_tPhaseMoles_old[iph] <= 0.0) {
if (m_SSPhase[kspec]) {
return VCS_SPECIES_MAJOR;
} else {
return VCS_SPECIES_ZEROEDMS;
}
}
}
// ---------- Treat species with non-zero mole numbers next ------------
@ -3822,6 +3801,9 @@ namespace VCSnonideal {
* phase and shares a non-zero stoichiometric coefficient, then
* the current species is a major species.
*/
if (irxn < 0) {
return VCS_SPECIES_MAJOR;
} else {
double szAdj = m_scSize[irxn] * std::sqrt((double)m_numRxnTot);
for (int k = 0; k < m_numComponents; ++k) {
if (!(m_SSPhase[k])) {
@ -3831,7 +3813,8 @@ namespace VCSnonideal {
}
}
}
}
}
}
return VCS_SPECIES_MINOR;
}
/*****************************************************************************/
@ -4579,13 +4562,19 @@ namespace VCSnonideal {
* These species will formally always have zero
* mole numbers in the solution vector.
* - VCS_SPECIES_ZEROEDPHASE
* -7 -> The species lies in a multicomponent phase which
* -7 -> The species lies in a multicomponent phase which
* currently does exist. Its concentration is currently
* identically zero, though the phase exists. Note, this
* is a temporary condition that exists at the start
* of an equilibrium problem.
* The species is soon "birthed" or "deleted".
* - VCS_SPECIES_ACTIVEBUTZERO
* -8 -> 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
* - VCS_SPECIES_STOICHZERO
*
*/
bool VCS_SOLVE::vcs_evaluate_speciesType() {
int kspec;
@ -4602,8 +4591,19 @@ namespace VCSnonideal {
#endif
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
m_speciesStatus[kspec] = vcs_species_type(kspec);
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 5) {
plogf(" --- %-16s: ", m_speciesName[kspec].c_str());
if (kspec < m_numComponents) {
plogf("(COMP) ");
} else {
plogf(" ");
}
plogf(" %10.3g ", m_molNumSpecies_old[kspec]);
const char *sString = vcs_speciesType_string(m_speciesStatus[kspec], 100);
plogf("%s\n", sString);
} else if (m_debug_print_lvl >= 2) {
if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
switch (m_speciesStatus[kspec]) {
case VCS_SPECIES_COMPONENT:
@ -4628,6 +4628,10 @@ namespace VCSnonideal {
plogf(" --- Zeroed Species in an active MS phase (tmp): %-s\n",
m_speciesName[kspec].c_str());
break;
case VCS_SPECIES_STOICHZERO:
plogf(" --- Zeroed Species in an active MS phase (Stoich Constraint): %-s\n",
m_speciesName[kspec].c_str());
break;
case VCS_SPECIES_INTERFACIALVOLTAGE:
plogf(" --- InterfaceVoltage Species: %-s\n", m_speciesName[kspec].c_str());
break;

73
Cantera/src/equil/vcs_util.cpp
View file

@ -413,8 +413,77 @@ namespace VCSnonideal {
plogendl();
}
/***************************************************************************/
/************************************************************************ **/
const char *vcs_speciesType_string(int speciesStatus, int length) {
char *sss;
switch (speciesStatus) {
case VCS_SPECIES_COMPONENT:
sss = "Component Species";
break;
case VCS_SPECIES_MAJOR:
sss ="Major Species";
break;
case VCS_SPECIES_MINOR:
sss ="Minor Species";
break;
case VCS_SPECIES_ZEROEDPHASE:
if (length < 48) {
sss = "Set Zeroed-Phase";
} else {
sss ="Purposely Zeroed-Phase Species (not in problem)";
}
break;
case VCS_SPECIES_ZEROEDMS:
if (length < 23) {
sss = "Zeroed-MS Phase";
} else {
sss ="Zeroed-MS Phase Species";
}
break;
case VCS_SPECIES_ZEROEDSS:
if (length < 23) {
sss = "Zeroed-SS Phase";
} else {
sss ="Zeroed-SS Phase Species";
}
break;
case VCS_SPECIES_DELETED:
if (length < 22) {
sss = "Deleted Species";
} else if (length < 40) {
sss = "Deleted-Small Species";
} else {
sss ="Deleted-Small Species in a MS phase";
}
break;
case VCS_SPECIES_ACTIVEBUTZERO:
if (length < 47) {
sss = "Tmp Zeroed in MS";
} else {
sss ="Zeroed Species in an active MS phase (tmp)";
}
break;
case VCS_SPECIES_STOICHZERO:
if (length < 56) {
sss = "Stoich Zeroed in MS";
} else {
sss ="Zeroed Species in an active MS phase (Stoich Constraint)";
}
break;
case VCS_SPECIES_INTERFACIALVOLTAGE:
if (length < 29) {
sss = "InterfaceVoltage";
} else {
sss ="InterfaceVoltage Species";
}
break;
default:
sss = "unknown species type";
}
return sss;
}
/************************************************************************ **/
void vcs_print_stringTrunc(const char *str, int space, int alignment)

10
config.h.in
View file

@ -103,9 +103,17 @@ typedef int ftnlen; // Fortran hidden string length type
#undef PYTHON_EXE
// If this is defined, the Cantera Python interface will use the
// Numeric package; otherwise, it will use numarray.
// Numeric package
#undef HAS_NUMERIC
// If this is defined, the Cantera Python interface will use the
// numarray package
#undef HAS_NUMARRAY
// If this is defined, the Cantera Python interface will use the
// numpy package
#undef HAS_NUMPY
// If this is defined, then python will not be assumed to be
// present to support conversions
#undef HAS_NO_PYTHON

10836
configure vendored
View file

File diff suppressed because one or more lines are too long

110
configure.in
View file

@ -1072,41 +1072,108 @@ fi
AC_SUBST(BUILD_PYTHON)
AC_DEFINE_UNQUOTED(PYTHON_EXE,"$PYTHON_CMD")
USE_NUMARRAY='y'
if test "$USE_NUMERIC" = "y"; then
USE_NUMARRAY='n'
USE_NUMPY='n'
AC_DEFINE(HAS_NUMERIC)
fi
if test -n "$NUMARRAY_INC_DIR" ; then
echo "setting NUMARRAY_INC_DIR to $NUMARRAY_INC_DIR"
else
NUMARRAY_INC_DIR=""
if test -n "$NUMARRAY_HOME"; then
dir5="$NUMARRAY_HOME/include/python2.5"
if test -d $dir5 ; then
NUMARRAY_INC_DIR=$dir5
if test "$USE_NUMPY" = "y"; then
USE_NUMARRAY='n'
AC_DEFINE(HAS_NUMPY)
fi
if test "$USE_NUMARRAY" = "y"; then
AC_DEFINE(HAS_NUMARRAY)
fi
if test "$USE_NUMARRAY" = "y"; then
if test -n "$NUMARRAY_INC_DIR" ; then
echo "setting NUMARRAY_INC_DIR to $NUMARRAY_INC_DIR"
else
dir4="$NUMARRAY_HOME/include/python2.4"
if test -d $dir4 ; then
NUMARRAY_INC_DIR=$dir4
NUMARRAY_INC_DIR=""
if test -n "$NUMARRAY_HOME"; then
dir5="$NUMARRAY_HOME/include/python2.5"
if test -d $dir5 ; then
NUMARRAY_INC_DIR=$dir5
else
dir3="$NUMARRAY_HOME/include/python2.3"
if test -d $dir3 ; then
NUMARRAY_INC_DIR=$dir3
dir4="$NUMARRAY_HOME/include/python2.4"
if test -d $dir4 ; then
NUMARRAY_INC_DIR=$dir4
else
dir1="$NUMARRAY_HOME/include/python"
if test -d $dir1 ; then
NUMARRAY_INC_DIR=$dir1
dir3="$NUMARRAY_HOME/include/python2.3"
if test -d $dir3 ; then
NUMARRAY_INC_DIR=$dir3
else
echo "WARNING: NUMARRAY include dir $dir1 does not exist."
NUMARRAY_INC_DIR=$dir1
dir1="$NUMARRAY_HOME/include/python"
if test -d $dir1 ; then
NUMARRAY_INC_DIR=$dir1
else
echo "WARNING: NUMARRAY include dir $dir1 does not exist."
NUMARRAY_INC_DIR=$dir1
fi
fi
fi
fi
echo "setting NUMARRAY_INC_DIR to $NUMARRAY_INC_DIR"
fi
fi
echo "setting NUMARRAY_INC_DIR to $NUMARRAY_INC_DIR"
fi
else
NUMARRAY_INC_DIR=""
NUMARRAY_HOME=""
fi
if test "$USE_NUMPY" = "y"; then
if test -n "$NUMPY_INC_DIR" ; then
echo "setting NUMPY_INC_DIR to $NUMPY_INC_DIR"
else
NUMPY_INC_DIR=""
if test -n "$NUMPY_HOME"; then
dir6="$NUMPY_HOME/include/python2.6"
if test -d $dir6 ; then
NUMPY_INC_DIR=$dir6
else
dir5="$NUMPY_HOME/include/python2.5"
if test -d $dir5 ; then
NUMPY_INC_DIR=$dir5
else
dir4="$NUMPY_HOME/include/python2.4"
if test -d $dir4 ; then
NUMPY_INC_DIR=$dir4
else
dir3="$NUMPY_HOME/include/python2.3"
if test -d $dir3 ; then
NUMPY_INC_DIR=$dir3
else
dir2="$NUMPY_HOME/include/python2.2"
if test -d $dir2 ; then
NUMPY_INC_DIR=$dir2
else
dir1="$NUMPY_HOME/include/python"
if test -d $dir1 ; then
NUMPY_INC_DIR=$dir1
else
echo "WARNING: NUMPY include dir $dir1 does not exist."
NUMPY_INC_DIR=$dir1
fi
fi
fi
fi
fi
fi
echo "setting NUMPY_INC_DIR to $NUMPY_INC_DIR"
fi
fi
else
NUMPY_INC_DIR=""
NUMPY_HOME=""
fi
AC_SUBST(NUMPY_INC_DIR)
AC_SUBST(NUMPY_HOME)
# this fails, at least on a Mac. By default, numarray include files
# are installed in the include directory in the Python framework. This
# does not require setting NUMARRAY_INC_DIR, so testing that it is
@ -1124,6 +1191,7 @@ fi
#fi
#fi
AC_SUBST(NUMARRAY_INC_DIR)
AC_SUBST(NUMARRAY_HOME)
AC_SUBST(CANTERA_PYTHON_HOME)

18
preconfig
View file

@ -72,15 +72,27 @@ PYTHON_PACKAGE=${PYTHON_PACKAGE:="minimal"}
# different name, set this to the full path to the Python interpreter.
PYTHON_CMD=${PYTHON_CMD:="default"}
# The Cantera Python interface can be built with either the numarray
# or Numeric packages. Set this to "y" to use Numeric, or anything
# else to use numarray. Using numarray is preferred.
# numeric, or numpy packages. Set this to "y" to use Numeric.
USE_NUMERIC=${USE_NUMERIC:="n"}
# The Cantera Python interface can be built with either the numarray
# or Numeric packages. Set this to "y" to use numpy, or 'n'
# to use numarray. set USE_NUMERIC to 'n' also if you want to use numarray.
# Using numpy is preferred, and is the supported option going forwards.
USE_NUMPY=${USE_NUMPY:="n"}
# If numpy was installed using the --home option, set this to the
# home directory for numpy. This will be needed for all numpy installations
# that don't put the include files into python's native include directory.
#NUMPY_HOME=${NUMPY_HOME:="$HOME/python_packages"}
# If numarray was installed using the --home option, set this to the
# home directory for numarray.
#NUMARRAY_HOME=${NUMARRAY_HOME:="$HOME/python_packages"}
# If you want to install the Cantera Python package somewhere other
# than the default 'site-packages' directory within the Python library
# directory, then set this to the desired directory. This is useful when
@ -543,7 +555,9 @@ export RPFONT
export FORTRAN_LIB_DIR
export CANTERA_INSTALL_DIR
export USE_NUMERIC
export USE_NUMPY
export NUMARRAY_HOME
export NUMPY_HOME
export CANTERA_PYTHON_HOME
export USE_SUNDIALS
export SUNDIALS_HOME

3869
test_problems/VCSnonideal/NaCl_equil/good_dout.txt
View file

File diff suppressed because it is too large Load diff

48
test_problems/VCSnonideal/NaCl_equil/good_out.txt
View file

@ -199,11 +199,11 @@ VCS CALCULATION METHOD
OH 3.7473450E-07 9.1461582E-08 -1.0769E+02 0
CO2 0.0000000E+00 0.0000000E+00 -5.0822E+02 0
Cl- 2.1230193E-01 9.1217525E-02 -7.2377E+01 KMolNum
H2O 9.7178675E-02 2.3718434E-02 -1.2400E+02 KMolNum
H2O 9.7178674E-02 2.3718434E-02 -1.2400E+02 KMolNum
H2 1.8736725E-07 4.5730791E-08 -3.2618E+01 KMolNum
H+ 1.8416226E-09 7.9127051E-10 -1.5174E+01 KMolNum
NaCl 3.9996395E-32 9.7619345E-33 -1.7452E+02 KMolNum
O2 8.9465966E-69 2.1835990E-69 -1.8277E+02 KMolNum
O2 8.9465967E-69 2.1835990E-69 -1.8277E+02 KMolNum
--------------------------------------------------------------------------------
@ -212,12 +212,12 @@ VCS CALCULATION METHOD
| Components| NaCl(S) N2 H2O(L) Na+ OH- OH CO2 | |
NonComponent | Moles | 4.79 4 1.9 0.212 1.84e-09 3.75e-07 0 | DG/RT Rxn |
-------------------------------------------------------------------------------------------------------------------
7 H+ | 1.84e-09 | 0.00 0.00 -1.00 0.00 1.00 0.00 0.00 | 9.65e-09 |
8 H2O | 0.0972 | 0.00 0.00 -1.00 0.00 0.00 0.00 0.00 | 2.34e-09 |
9 NaCl | 4e-32 | -1.00 0.00 0.00 0.00 0.00 0.00 0.00 | 5.73e-10 |
10 Cl- | 0.212 | -1.00 0.00 0.00 1.00 0.00 0.00 0.00 | -4.47e-09 |
11 H2 | 1.87e-07 | 0.00 0.00 -2.00 0.00 0.00 2.00 0.00 | 5.39e-09 |
12 O2 | 8.95e-69 | 0.00 0.00 2.00 0.00 0.00 -4.00 0.00 | 6.33e-08 |
7 H+ | 1.84e-09 | 0.00 0.00 -1.00 0.00 1.00 0.00 0.00 | 2.93e-09 |
8 H2O | 0.0972 | 0.00 0.00 -1.00 0.00 0.00 0.00 0.00 | 5.58e-10 |
9 NaCl | 4e-32 | -1.00 0.00 0.00 0.00 0.00 0.00 0.00 | 7.14e-10 |
10 Cl- | 0.212 | -1.00 0.00 0.00 1.00 0.00 0.00 0.00 | -5.8e-10 |
11 H2 | 1.87e-07 | 0.00 0.00 -2.00 0.00 0.00 2.00 0.00 | 1.21e-09 |
12 O2 | 8.95e-69 | 0.00 0.00 2.00 0.00 0.00 -4.00 0.00 | 7.89e-08 |
-------------------------------------------------------------------------------------------------------------------
@ -228,9 +228,9 @@ VCS CALCULATION METHOD
| Element | O H C N Na Cl cn_NaCl_el Fe E Si Ca | |
PhaseName |KMolTarget | 2 4 0 8 5 5 0 0 0 0 0 | Gibbs Total |
------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0 NaCl_electro | 2.327e+00 | 1.9 3.81 0 0 0.212 0.212 0 0 0 0 0 | -2.73006234858E+02 |
1 air | 4.097e+00 | 0.0972 0.194 0 8 4e-32 4e-32 0 0 0 0 0 | -1.04327434163E+02 |
2 NaCl(S) | 4.788e+00 | 0 0 0 0 4.79 4.79 0 0 0 0 0 | -8.35533676862E+02 |
0 NaCl_electro | 2.327e+00 | 1.9 3.81 0 0 0.212 0.212 0 0 0 0 0 | -2.73006234929E+02 |
1 air | 4.097e+00 | 0.0972 0.194 0 8 4e-32 4e-32 0 0 0 0 0 | -1.04327434136E+02 |
2 NaCl(S) | 4.788e+00 | 0 0 0 0 4.79 4.79 0 0 0 0 0 | -8.35533676818E+02 |
------------------------------------------------------------------------------------------------------------------------------------------------------------------------
TOTAL | 1.121e+01 | 2 4 0 8 5 5 0 0 0 0 0 | -1.21286734588E+03 |
------------------------------------------------------------------------------------------------------------------------------------------------------------------------
@ -257,32 +257,32 @@ Chemical Potentials of the Species: (dimensionless)
(RT = 2.47896e+06 J/kmol)
Name TKMoles StandStateChemPot ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)| (MolNum ChemPot)|
---------------------------------------------------------------------------------------------------------------------------------------------------
NaCl(S) 4.7876981E+00 -1.7451679E+02 0.0000000E+00 0.0000000E+00 | 0.0000000E+00 | -1.7452E+02 | | -8.355336769E+02 |
NaCl(S) 4.7876981E+00 -1.7451679E+02 0.0000000E+00 0.0000000E+00 | 0.0000000E+00 | -1.7452E+02 | | -8.355336768E+02 |
N2 4.0000000E+00 -2.3045225E+01 0.0000000E+00 -2.4004385E-02 | 0.0000000E+00 | -2.3069E+01 | | -9.227691608E+01 |
H2O(L) 1.9028209E+00 -1.2371551E+02 -8.6317425E-02 -2.0142493E-01 | 0.0000000E+00 | -1.2400E+02 | | -2.359559843E+02 |
Na+ 2.1230193E-01 -1.0397591E+02 2.1455966E-01 -2.3945082E+00 | 0.0000000E+00 | -1.0214E+02 |(4.01653E+00)| -2.168437492E+01 |
OH- 1.8416226E-09 -9.1483483E+01 -4.0506365E-01 -2.0957381E+01 |-0.0000000E+00 | -1.0883E+02 |(4.01653E+00)| -2.004226729E-07 |
OH 3.7473450E-07 -9.1487045E+01 0.0000000E+00 -1.6207347E+01 | 0.0000000E+00 | -1.0769E+02 | | -4.035680387E-05 |
H2O(L) 1.9028209E+00 -1.2371551E+02 -8.6317426E-02 -2.0142493E-01 | 0.0000000E+00 | -1.2400E+02 | | -2.359559843E+02 |
Na+ 2.1230193E-01 -1.0397591E+02 2.1455966E-01 -2.3945082E+00 | 0.0000000E+00 | -1.0214E+02 |(4.01653E+00)| -2.168437494E+01 |
OH- 1.8416226E-09 -9.1483483E+01 -4.0506365E-01 -2.0957381E+01 |-0.0000000E+00 | -1.0883E+02 |(4.01653E+00)| -2.004226720E-07 |
OH 3.7473450E-07 -9.1487045E+01 0.0000000E+00 -1.6207347E+01 | 0.0000000E+00 | -1.0769E+02 | | -4.035680380E-05 |
CO2 0.0000000E+00 -1.8445182E+02 0.0000000E+00 -3.2377221E+02 | 0.0000000E+00 | -5.0822E+02 | | -0.000000000E+00 |
Cl- 2.1230193E-01 -7.4214051E+01 2.1455966E-01 -2.3945082E+00 |-0.0000000E+00 | -7.2377E+01 |(4.01653E+00)| -1.536587544E+01 |
H2O 9.7178675E-02 -1.2026175E+02 0.0000000E+00 -3.7415027E+00 | 0.0000000E+00 | -1.2400E+02 | | -1.205047162E+01 |
H2 1.8736725E-07 -1.5717224E+01 0.0000000E+00 -1.6900494E+01 | 0.0000000E+00 | -3.2618E+01 | | -6.111492083E-06 |
H+ 1.8416226E-09 0.0000000E+00 1.7669883E+00 -2.0957381E+01 | 0.0000000E+00 | -1.5174E+01 |(4.01653E+00)| -2.794452014E-08 |
NaCl 3.9996395E-32 -1.0080997E+02 0.0000000E+00 -7.3706817E+01 | 0.0000000E+00 | -1.7452E+02 | | -6.980042318E-30 |
O2 8.9465966E-69 -2.4673669E+01 0.0000000E+00 -1.5809740E+02 | 0.0000000E+00 | -1.8277E+02 | | -1.635178995E-66 |
Cl- 2.1230193E-01 -7.4214051E+01 2.1455966E-01 -2.3945082E+00 |-0.0000000E+00 | -7.2377E+01 |(4.01653E+00)| -1.536587545E+01 |
H2O 9.7178674E-02 -1.2026175E+02 0.0000000E+00 -3.7415027E+00 | 0.0000000E+00 | -1.2400E+02 | | -1.205047159E+01 |
H2 1.8736725E-07 -1.5717224E+01 0.0000000E+00 -1.6900494E+01 | 0.0000000E+00 | -3.2618E+01 | | -6.111492072E-06 |
H+ 1.8416226E-09 0.0000000E+00 1.7669883E+00 -2.0957381E+01 | 0.0000000E+00 | -1.5174E+01 |(4.01653E+00)| -2.794452002E-08 |
NaCl 3.9996395E-32 -1.0080997E+02 0.0000000E+00 -7.3706817E+01 | 0.0000000E+00 | -1.7452E+02 | | -6.980042319E-30 |
O2 8.9465967E-69 -2.4673669E+01 0.0000000E+00 -1.5809740E+02 | 0.0000000E+00 | -1.8277E+02 | | -1.635179010E-66 |
-1.212867346E+03
---------------------------------------------------------------------------------------------------------------------------------------------------
Counters: Iterations Time (seconds)
vcs_basopt: 3 NA
vcs_TP: 44 NA
vcs_TP: 31 NA
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
TCounters: Num_Calls Total_Its Total_Time (seconds)
vcs_basopt: 3 3 NA
vcs_TP: 1 44 NA
vcs_TP: 1 31 NA
vcs_inest: 0 NA
vcs_TotalTime: NA

2
test_problems/VCSnonideal/NaCl_equil/vcs_equilibrate_blessed.csv
View file

@ -3,7 +3,7 @@ Temperature = 298.15 kelvin
Pressure = 1.0132e+05 Pascal
Total Volume = 100.41 m**3
Number Basis optimizations = 3
Number VCS iterations = 44
Number VCS iterations = 31
Name, Phase, PhaseMoles, Mole_Fract, Molalities, ActCoeff, Activity,ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume
, , (kmol), , , , , (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)
H2O(L), NaCl_electrolyte, 2.327e+00, 8.176e-01, 5.551e+01, 9.173e-01, 7.500e-01, -3.067e+02, -3.074e+02, 1.903e+00, 1.800e-02, 3.838e-02

1 --------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT -----------------------------
3 Pressure = 1.0132e+05 Pascal
4 Total Volume = 100.41 m**3
5 Number Basis optimizations = 3
6 Number VCS iterations = 44 Number VCS iterations = 31
7 Name, Phase, PhaseMoles, Mole_Fract, Molalities, ActCoeff, Activity,ChemPot_SS0, ChemPot, mole_num, PMVol, Phase_Volume
8 , , (kmol), , , , , (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)
9 H2O(L), NaCl_electrolyte, 2.327e+00, 8.176e-01, 5.551e+01, 9.173e-01, 7.500e-01, -3.067e+02, -3.074e+02, 1.903e+00, 1.800e-02, 3.838e-02