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

Working towards fixing an algorithm bug.

Browse source
This commit is contained in:
Harry Moffat 2009-04-17 00:40:48 +00:00
parent eaf88a6b70
commit d560ccd175
4 changed files with 103 additions and 21 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
Cantera/src/equil/MultiPhase.cpp
View file

@ -230,6 +230,18 @@ namespace Cantera {
return sum;
}
int MultiPhase::speciesIndex(std::string speciesName, std::string phaseName) {
int p = phaseIndex(phaseName);
if (p < 0) {
throw CanteraError("MultiPhase::speciesIndex", "phase not found: " + phaseName);
}
int k = m_phase[p]->speciesIndex(speciesName);
if (k < 0) {
throw CanteraError("MultiPhase::speciesIndex", "species not found: " + speciesName);
}
return m_spstart[p] + k;
}
/// Net charge of one phase (Coulombs). The net charge is computed as
/// \f[ Q_p = N_p \sum_k F z_k X_k \f]
/// where the sum runs only over species in phase \a p.
@ -459,6 +471,16 @@ namespace Cantera {
}
}
void MultiPhase::addSpeciesMoles(const int indexS, const doublereal addedMoles) {
vector_fp tmpMoles(m_nsp, 0.0);
getMoles(DATA_PTR(tmpMoles));
tmpMoles[indexS] += addedMoles;
if (tmpMoles[indexS] < 0.0) {
tmpMoles[indexS] = 0.0;
}
setMoles(DATA_PTR(tmpMoles));
}
void MultiPhase::setState_TP(const doublereal T, const doublereal Pres) {
if (!m_init) init();
m_temp = T;

22
Cantera/src/equil/MultiPhase.h
View file

@ -231,6 +231,21 @@ namespace Cantera {
return m_spstart[p] + k;
}
//! Return the global index of the species belonging to phase name \c phaseName
//! with species name \c speciesName
/*!
* Returns the index of the global species
*
* @param speciesName Species Name
* @param phaseName Phase Name
*
* @return returns the global index
*
* If the species or phase name is not recognized, this routine throws
* a CanteraError.
*/
int speciesIndex(std::string speciesName, std::string phaseName);
/// Minimum temperature for which all solution phases have
/// valid thermo data. Stoichiometric phases are not
/// considered, since they may have thermo data only valid for
@ -478,6 +493,13 @@ namespace Cantera {
*/
void setMoles(const doublereal* n);
//! Adds moles of a certain species to the mixture
/*!
*
*/
void addSpeciesMoles(const int indexS, const doublereal addedMoles);
//! Retrieves a vector of element abundances
/*!
* @param elemAbundances Vector of element abundances

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

@ -1278,8 +1278,8 @@ private:
//! energy by making sure the slope of the following functional stays
//! negative:
/*!
* The slope of the following functional is equivalent to the slope of the total
* Gibbs free energy of the system:
* The slope of the following functional is equivalent to the slope
* of the total Gibbs free energy of the system:
*
* d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] )
*
@ -1768,17 +1768,51 @@ public:
* to is
* kspec = irxn + m_numComponents
*
* kspec : 1 -> Major species VCS_SPECIES_MAJOR
* kspec : 2 -> Component species VCS_SPECIES_COMPONENT
* -> deprecated, want to assign -2 to some
* component species. We can already determine
* whether the species is a component from
* its position in the species vector.
* 1 -> Major species VCS_SPECIES_MAJOR
* 0 -> Minor species VCS_SPECIES_MINOR
* -1 -> Mole number is zero
* in inactive phase VCS_SPECIES_ZEROEDPHASE
* -2 -> Deleted species in an
* active multicom phase VCS_SPECIES_ZEROEDMS
* -3 -> Mole number is zero
* in a stoich phase - VCS_SPECIES_ZEREODSS
* -4 -> Species is deleted
* -1 -> The species lies in a multicomponent phase
* that exists. Its concentration is currently
* very low, necessitating a different method
* of calculation.
* - VCS_SPECIES_ZEROEDPHASE
* -2 -> The species lies in a multicomponent phase
* which currently doesn't exist.
* Its concentration is currently zero.
* - VCS_SPECIES_ZEROEDMS
* -3 -> Species lies in a single-species phase which
* is currently zereod out.
* - VCS_SPECIES_ZEREODSS
* -4 -> Species has such a small mole fraction it is
* deleted even though its phase may possibly exist.
* The species is believed to have such a small
* mole fraction that it best to throw the
* calculation of it out. It will be added back in
* at the end of the calculation.
* - VCS_SPECIES_DELETED
* -> Length equal to number of species
* -5 -> Species refers to an electron in the metal
* The unknown is equal to the interfacial voltage
* drop across the interface on the SHE (standard
* hydroogen electrode) scale (volts).
* - VCS_SPECIES_INTERFACIALVOLTAGE
* -6 -> Species lies in a multicomponent phase that
* is zeroed atm and will stay deleted due to a
* choice from a higher level.
* 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
* 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
*
*/
std::vector<int> m_speciesStatus;

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

@ -292,24 +292,25 @@ namespace VCSnonideal {
MajorSpeciesHaveConverged = false;
/*************************************************************************/
/************** EVALUATE INITIAL MAJOR-MINOR VECTOR **********************/
/************** EVALUATE INITIAL SPECIES STATUS VECTOR *******************/
/*************************************************************************/
m_numRxnMinorZeroed = 0;
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- MAJOR-MINOR decision is reavaluated: All species are minor except for:\n");
plogf(" --- Species Status decision is reavaluated: All species are minor except for:\n");
} else if (m_debug_print_lvl >= 5) {
plogf(" --- MAJOR-MINOR decision is reavaluated");
plogf(" --- Species Status decision is reavaluated");
plogendl();
}
#endif
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
kspec = m_indexRxnToSpecies[irxn];
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
m_speciesStatus[kspec] = vcs_species_type(kspec);
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
switch (m_speciesStatus[kspec]) {
case VCS_SPECIES_COMPONENT:
break;
case VCS_SPECIES_MAJOR:
plogf(" --- Major Species : %-s\n", m_speciesName[kspec].c_str());
break;
@ -327,7 +328,8 @@ namespace VCSnonideal {
plogf(" --- Deleted-Small Species : %-s\n", m_speciesName[kspec].c_str());
break;
case VCS_SPECIES_ACTIVEBUTZERO:
plogf(" --- Zeroed Species in an active MS phase (tmp): %-s\n", m_speciesName[kspec].c_str());
plogf(" --- Zeroed Species in an active MS phase (tmp): %-s\n",
m_speciesName[kspec].c_str());
break;
case VCS_SPECIES_INTERFACIALVOLTAGE:
plogf(" --- InterfaceVoltage Species: %-s\n", m_speciesName[kspec].c_str());
@ -340,8 +342,10 @@ namespace VCSnonideal {
}
}
#endif
if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
++m_numRxnMinorZeroed;
if (kspec >= m_numComponents) {
if (m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) {
++m_numRxnMinorZeroed;
}
}
}
#ifdef DEBUG_MODE
@ -1619,7 +1623,7 @@ namespace VCSnonideal {
if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
if (m_speciesStatus[kspec] == VCS_SPECIES_MAJOR) {
plogf(" --- Noncomponent turned from major to minor: ");
} else if (m_speciesStatus[kspec] == VCS_SPECIES_COMPONENT) {
} else if (kspec < m_numComponents) {
plogf(" --- Component turned into a minor species: ");
} else {
plogf(" --- Zeroed Species turned into a "
@ -1636,7 +1640,7 @@ namespace VCSnonideal {
if (m_debug_print_lvl >= 2) {
if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) {
plogf(" --- Noncomponent turned from minor to major: ");
} else if (m_speciesStatus[kspec] == VCS_SPECIES_COMPONENT) {
} else if (kspec < m_numComponents) {
plogf(" --- Component turned into a major: ");
} else {
plogf(" --- Noncomponent turned from zeroed to major: ");