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

Worked on making sure that vcs_deltag() adhered to the STATECALC concept.

Browse source
This commit is contained in:
Harry Moffat 2008-06-02 23:24:32 +00:00
parent beadead8b6
commit 2ddf9a409f
6 changed files with 152 additions and 101 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -82,11 +82,11 @@ namespace VCSnonideal {
* @{
*/
//! Cutoff relative mole number value,
//! Cutoff relative mole fraction value,
//! below which species are deleted from the equilibrium problem.
#ifndef VCS_DELETE_SPECIES_CUTOFF
#define VCS_DELETE_SPECIES_CUTOFF 1.0e-32
#ifndef VCS_RELDELETE_SPECIES_CUTOFF
#define VCS_RELDELETE_SPECIES_CUTOFF 1.0e-64
#endif
//! Cutoff relative mole number value,

46
Cantera/src/equil/vcs_inest.cpp
View file

@ -51,7 +51,7 @@ namespace VCSnonideal {
int nrxn = m_numRxnTot;
vcs_VolPhase *Vphase = 0;
double *molNum = VCS_DATA_PTR(m_molNumSpecies_old);
// double *molNum = VCS_DATA_PTR(m_molNumSpecies_old);
double TMolesMultiphase;
double *xtphMax = VCS_DATA_PTR(m_TmpPhase);
double *xtphMin = VCS_DATA_PTR(m_TmpPhase2);
@ -104,7 +104,7 @@ namespace VCSnonideal {
plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix);
for (kspec = 0; kspec < nspecies; ++kspec) {
plogf("%s ", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %15.5g %12.3g\n", molNum[kspec], -m_SSfeSpecies[kspec]);
plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]);
}
plogf("%s Element Abundance Agreement returned from linear "
"programming (vcs_inest initial guess):",
@ -116,7 +116,7 @@ namespace VCSnonideal {
if (m_elementActive[j]) {
double tmp = 0.0;
for (kspec = 0; kspec < nspecies; ++kspec) {
tmp += m_formulaMatrix[j][kspec] * molNum[kspec];
tmp += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec];
}
plogf("%s ", pprefix); plogf(" %-9.9s", (m_elementName[j]).c_str());
plogf(" %12.3g %12.3g\n", m_elemAbundancesGoal[j], tmp);
@ -137,18 +137,18 @@ namespace VCSnonideal {
iph = m_phaseID[kspec];
Vphase = m_VolPhaseList[iph];
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (molNum[kspec] <= 0.0) {
if (m_molNumSpecies_old[kspec] <= 0.0) {
/*
* HKM Should eventually include logic here for non SS phases
*/
if (!m_SSPhase[kspec]) {
molNum[kspec] = 1.0e-30;
m_molNumSpecies_old[kspec] = 1.0e-30;
}
}
} else {
molNum[kspec] = 0.0;
m_molNumSpecies_old[kspec] = 0.0;
}
if (molNum[kspec] > 0.0) {
if (m_molNumSpecies_old[kspec] > 0.0) {
if (Vphase->Existence == 0) {
Vphase->Existence = 1;
}
@ -179,7 +179,7 @@ namespace VCSnonideal {
}
for (kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
m_tPhaseMoles_new[m_phaseID[kspec]] += molNum[kspec];
m_tPhaseMoles_new[m_phaseID[kspec]] += m_molNumSpecies_old[kspec];
}
}
TMolesMultiphase = 0.0;
@ -188,36 +188,36 @@ namespace VCSnonideal {
TMolesMultiphase += m_tPhaseMoles_new[iph];
}
}
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNum, nspecies);
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_molNumSpecies_old), nspecies);
for (kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_MOLNUM) {
m_molNumSpecies_new[kspec] = 0.0;
}
}
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(m_SSfeSpecies),
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_SSfeSpecies),
nspecies);
for (kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
if (! m_SSPhase[kspec]) {
iph = m_phaseID[kspec];
m_feSpecies_curr[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]);
m_feSpecies_new[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]);
}
} else {
m_molNumSpecies_new[kspec] = 0.0;
}
}
vcs_deltag(0, true);
vcs_deltag(0, true, VCS_STATECALC_NEW);
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
for (kspec = 0; kspec < nspecies; ++kspec) {
plogf("%s", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str());
if (kspec < m_numComponents)
plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_curr[kspec],
plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_new[kspec],
m_SSfeSpecies[kspec]);
else
plogf("fe* = %15.5g ff = %15.5g dg* = %15.5g\n",
m_feSpecies_curr[kspec], m_SSfeSpecies[kspec], m_deltaGRxn_new[kspec-m_numComponents]);
m_feSpecies_new[kspec], m_SSfeSpecies[kspec], m_deltaGRxn_new[kspec-m_numComponents]);
}
}
#endif
@ -269,8 +269,8 @@ namespace VCSnonideal {
plogf("%sdirection (", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(") = %g", m_deltaMolNumSpecies[kspec]);
if (m_SSPhase[kspec]) {
if (molNum[kspec] > 0.0) {
plogf(" (ssPhase exists at w = %g moles)", molNum[kspec]);
if (m_molNumSpecies_old[kspec] > 0.0) {
plogf(" (ssPhase exists at w = %g moles)", m_molNumSpecies_old[kspec]);
} else {
plogf(" (ssPhase doesn't exist -> stability not checked)");
}
@ -304,14 +304,16 @@ namespace VCSnonideal {
do {
for (kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
molNum[kspec] = m_molNumSpecies_new[kspec] + par * m_deltaMolNumSpecies[kspec];
m_molNumSpecies_old[kspec] = m_molNumSpecies_new[kspec] + par * m_deltaMolNumSpecies[kspec];
} else {
m_deltaMolNumSpecies[kspec] = 0.0;
}
}
for (kspec = m_numComponents; kspec < nspecies; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (m_deltaMolNumSpecies[kspec] != 0.0) molNum[kspec] = m_deltaMolNumSpecies[kspec] * par;
if (m_deltaMolNumSpecies[kspec] != 0.0) {
m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par;
}
}
}
/*
@ -323,9 +325,9 @@ namespace VCSnonideal {
/* ******************************************* */
/* **** CONVERGENCE FORCING SECTION ********** */
/* ******************************************* */
vcs_dfe(molNum, VCS_STATECALC_OLD, 0, 0, nspecies);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, nspecies);
for (kspec = 0, s = 0.0; kspec < nspecies; ++kspec) {
s += m_deltaMolNumSpecies[kspec] * m_feSpecies_curr[kspec];
s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec];
}
if (s == 0.0) {
finished = TRUE; continue;
@ -377,7 +379,7 @@ namespace VCSnonideal {
plogf("%s SPECIES MOLE_NUMBER\n", pprefix);
for (kspec = 0; kspec < nspecies; ++kspec) {
plogf("%s ", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %g", molNum[kspec]);
plogf(" %g", m_molNumSpecies_old[kspec]);
plogendl();
}
}
@ -520,7 +522,7 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix,
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_new),
VCS_DATA_PTR(m_tPhaseMoles_old)));
plogendl();
}

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

@ -325,7 +325,7 @@ namespace VCSnonideal {
* @param vec_to vector of doubles
* @param length length of the vector to zero.
*/
inline void vcs_dzero(double *vec_to, int length) {
inline void vcs_dzero(double * const vec_to, const int length) {
(void) memset((void *) vec_to, 0, length * sizeof(double));
}
@ -334,7 +334,7 @@ namespace VCSnonideal {
* @param vec_to vector of ints
* @param length length of the vector to zero.
*/
inline void vcs_izero(int *vec_to, int length) {
inline void vcs_izero(int * const vec_to, const int length) {
(void) memset((void *) vec_to, 0, length * sizeof(int));
}
@ -345,7 +345,8 @@ namespace VCSnonideal {
* @param vec_from Vector to copy from
* @param length Number of doubles to copy.
*/
inline void vcs_dcopy(double *vec_to, const double *vec_from, int length) {
inline void vcs_dcopy(double * const vec_to,
const double * const vec_from, const int length) {
(void) memcpy((void *) vec_to, (const void *) vec_from,
(length) * sizeof(double));
}
@ -358,7 +359,8 @@ namespace VCSnonideal {
* @param vec_from Vector to copy from
* @param length Number of int to copy.
*/
inline void vcs_icopy(int *vec_to, const int *vec_from, int length) {
inline void vcs_icopy(int * const vec_to,
const int * const vec_from, const int length) {
(void) memcpy((void *) vec_to, (const void *) vec_from,
(length) * sizeof(int));
}
@ -368,7 +370,7 @@ namespace VCSnonideal {
* @param vec_to vector of doubles
* @param length length of the vector to zero.
*/
inline void vcs_vdzero(std::vector<double> &vec_to, int length) {
inline void vcs_vdzero(std::vector<double> &vec_to, const int length) {
(void) memset((void *)VCS_DATA_PTR(vec_to), 0, (length) * sizeof(double));
}
@ -377,7 +379,7 @@ namespace VCSnonideal {
* @param vec_to vector of ints
* @param length length of the vector to zero.
*/
inline void vcs_vizero(std::vector<int> &vec_to, int length) {
inline void vcs_vizero(std::vector<int> &vec_to, const int length) {
(void) memset((void *)VCS_DATA_PTR(vec_to), 0, (length) * sizeof(int));
}
@ -409,21 +411,21 @@ namespace VCSnonideal {
* @param length Number of integers to copy.
*/
inline void vcs_vicopy(std::vector<int> & vec_to,
const std::vector<int> & vec_from, int length) {
const std::vector<int> & vec_from, const int length) {
(void) memcpy((void *)&(vec_to[0]), (const void *) &(vec_from[0]),
(length) * sizeof(int));
}
#else
extern void vcs_dzero(double *, int);
extern void vcs_izero(int *, int);
extern void vcs_dcopy(double *, double *, int);
extern void vcs_icopy(int *, int *, int);
extern void vcs_vdzero(std::vector<double> &vvv, int len = -1);
extern void vcs_vizero(std::vector<double> &vvv, int len = -1);
extern void vcs_dzero(double * const, const int);
extern void vcs_izero(int * const , const int);
extern void vcs_dcopy(double * const, const double * const, const int);
extern void vcs_icopy(int * const, const int * const, const int);
extern void vcs_vdzero(std::vector<double> &vvv, const int len = -1);
extern void vcs_vizero(std::vector<double> &vvv, const int len = -1);
void vcs_vdcopy(std::vector<double> &vec_to,
const std::vector<double> vec_from, int len = -1);
const std::vector<double> vec_from, const int len = -1);
void vcs_vicopy(std::vector<int> &vec_to,
const std::vector<int> vec_from, int len = -1);
const std::vector<int> vec_from, const int len = -1);
#endif
//! Finds the location of the maximum component in a double vector

3
Cantera/src/equil/vcs_report.cpp
View file

@ -329,7 +329,8 @@ namespace VCSnonideal {
lx = 0.0;
} else {
if (tpmoles > 0.0 && m_molNumSpecies_old[l] > 0.0) {
lx = log(m_molNumSpecies_old[l]) - log(tpmoles);
double tmp = MAX(VCS_DELETE_MINORSPECIES_CUTOFF, m_molNumSpecies_old[l]);
lx = log(tmp) - log(tpmoles);
} else {
lx = m_feSpecies_old[l] - m_SSfeSpecies[l]
- log(m_actCoeffSpecies_old[l]) + m_lnMnaughtSpecies[l];

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

@ -504,7 +504,37 @@ public:
*/
void vcs_tmoles();
void vcs_deltag(int l, bool doDeleted);
//! This subroutine calculates reaction free energy changes for
//! all noncomponent formation reactions.
/*!
* Formation reactions are
* reactions which create each noncomponent species from the component
* species. m_stoichCoeffRxnMatrix[irxn][jcomp] are the stoichiometric
* coefficients for these reactions. A stoichiometric coefficient of
* one is assumed for species irxn in this reaction.
*
* INPUT
* @param l
* L < 0 : Calculate reactions corresponding to
* major noncomponent and zeroed species only
* L = 0 : Do all noncomponent reactions, i, between
* 0 <= i < irxnl
* L > 0 : Calculate reactions corresponding to
* minor noncomponent and zeroed species only
*
* @param doDeleted Do deleted species
* @param stateCalc Calculate deltaG corresponding to either old or new
* free energies
*
* Note we special case one important issue.
* If the component has zero moles, then we do not
* allow deltaG < 0.0 for formation reactions which
* would lead to the loss of more of that same component.
* This dG < 0.0 condition feeds back into the algorithm in several
* places, and leads to a infinite loop in at least one case.
*/
void vcs_deltag(const int l, const bool doDeleted, const int stateCalc);
//! Swaps the indecises for all of the global data for two species, k1
//! and k2.

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

@ -267,7 +267,12 @@ namespace VCSnonideal {
"negative MF, setting it small",
m_speciesName[i].c_str());
plogendl();
m_molNumSpecies_old[i] = VCS_DELETE_SPECIES_CUTOFF;
iph = m_phaseID[i];
double tmp = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10;
if ( VCS_DELETE_MINORSPECIES_CUTOFF*10. > tmp) {
tmp = VCS_DELETE_MINORSPECIES_CUTOFF*10.;
}
m_molNumSpecies_old[i] = tmp;
}
}
@ -390,7 +395,7 @@ namespace VCSnonideal {
#endif
// Update the phase objects with the contents of the soln vector
vcs_updateVP(0);
vcs_deltag(0, false);
vcs_deltag(0, false, VCS_STATECALC_OLD);
iti = 0;
goto L_MAINLOOP_ALL_SPECIES;
@ -416,7 +421,7 @@ namespace VCSnonideal {
*/
if (uptodate_minors == FALSE) {
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
vcs_deltag(1, false);
vcs_deltag(1, false, VCS_STATECALC_NEW);
}
uptodate_minors = TRUE;
} else {
@ -435,10 +440,13 @@ namespace VCSnonideal {
}
}
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_old), VCS_DATA_PTR(m_actCoeffSpecies_new), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc);
/*
* Copy the old solution into the new solution as an initial guess
*/
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(m_feSpecies_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);
/* Go find a new reaction adjustment ->
* i.e., change in extent of reaction for each reaction.
@ -1094,7 +1102,7 @@ namespace VCSnonideal {
* Evaluate DeltaG for all components if ITI=0, and for
* major components only if ITI NE 0
*/
vcs_deltag(0, false);
vcs_deltag(0, false, VCS_STATECALC_NEW);
/* *************************************************************** */
/* **** CONVERGENCE FORCER SECTION ******************************* */
@ -1231,6 +1239,8 @@ namespace VCSnonideal {
vcs_dcopy(VCS_DATA_PTR(m_tPhaseMoles_old), VCS_DATA_PTR(m_tPhaseMoles_new), m_numPhases);
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_molNumSpecies_new),
m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_old),
VCS_DATA_PTR(m_actCoeffSpecies_new), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc);
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc);
@ -1268,7 +1278,7 @@ namespace VCSnonideal {
if (m_phaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) {
irxn = kspec - m_numComponents;
if (kspec < m_numComponents) {
if (m_molNumSpecies_old[kspec] > VCS_DELETE_SPECIES_CUTOFF) {
if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF) {
soldel = 0;
break;
}
@ -1311,7 +1321,7 @@ namespace VCSnonideal {
&usedZeroedSpecies);
if (retn != VCS_SUCCESS) return retn;
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, true);
vcs_deltag(0, true, VCS_STATECALC_OLD);
uptodate_minors = TRUE;
if (conv) {
/*
@ -1344,7 +1354,7 @@ namespace VCSnonideal {
#endif
vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx));
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, true);
vcs_deltag(0, true, VCS_STATECALC_OLD);
uptodate_minors = TRUE;
}
#ifdef DEBUG_MODE
@ -1579,7 +1589,7 @@ namespace VCSnonideal {
*/
if (iti != 0) {
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, kspec, kspec+1);
vcs_deltag(0, false);
vcs_deltag(0, false, VCS_STATECALC_OLD);
}
}
}
@ -1660,7 +1670,7 @@ namespace VCSnonideal {
*/
if (iti != 0) {
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
vcs_deltag(1, false);
vcs_deltag(1, false, VCS_STATECALC_OLD);
uptodate_minors = TRUE;
}
#ifdef DEBUG_MODE
@ -1770,7 +1780,7 @@ namespace VCSnonideal {
* Go back to evaluate the total moles of gas and liquid.
*/
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
vcs_deltag(0, false, VCS_STATECALC_OLD);
/*
*
*/
@ -1829,7 +1839,7 @@ namespace VCSnonideal {
}
}
// Calculate delta g's
vcs_deltag(0, false);
vcs_deltag(0, false, VCS_STATECALC_OLD);
// Go back to equilibrium check as a prep to eventually checking out
goto L_EQUILIB_CHECK;
@ -1860,7 +1870,7 @@ namespace VCSnonideal {
*/
MajorSpeciesHaveConverged = true;
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
vcs_deltag(0, false, VCS_STATECALC_OLD);
iti = 0;
goto L_MAINLOOP_ALL_SPECIES;
/*************************************************************************/
@ -1879,7 +1889,7 @@ namespace VCSnonideal {
*/
MajorSpeciesHaveConverged = true;
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
vcs_deltag(0, false, VCS_STATECALC_OLD);
iti = 0;
goto L_MAINLOOP_ALL_SPECIES;
}
@ -2336,7 +2346,7 @@ namespace VCSnonideal {
*/
void VCS_SOLVE::vcs_reinsert_deleted(int kspec) {
int i, k, irxn = kspec - m_numComponents;
int *phaseID = VCS_DATA_PTR(m_phaseID);
int iph = m_phaseID[kspec];
double dx;
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
@ -2348,7 +2358,7 @@ namespace VCSnonideal {
* this adjusts m_molNumSpecies_old[] and m_tPhaseMoles_old[]
* HKM -> make this a relative mole number!
*/
dx = VCS_DELETE_SPECIES_CUTOFF * 10.;
dx = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10.;
delta_species(kspec, &dx);
m_rxnStatus[irxn] = VCS_SPECIES_MINOR;
@ -2356,7 +2366,7 @@ namespace VCSnonideal {
m_rxnStatus[irxn] = VCS_SPECIES_MAJOR;
--(m_numRxnMinorZeroed);
}
int iph = m_phaseID[kspec];
vcs_VolPhase *Vphase = m_VolPhaseList[iph];
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(m_tPhaseMoles_old));
@ -2373,7 +2383,7 @@ namespace VCSnonideal {
if (Vphase->Existence == 0) {
Vphase->Existence = 1;
for (k = 0; k < m_numSpeciesTot; k++) {
if (phaseID[k] == iph) {
if (m_phaseID[k] == iph) {
i = k - m_numComponents;
if (m_rxnStatus[i] == VCS_SPECIES_ZEROEDPHASE)
m_rxnStatus[i] = VCS_SPECIES_ZEROEDMS;
@ -2547,7 +2557,7 @@ namespace VCSnonideal {
*/
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
iph = m_phaseID[kspec];
m_feSpecies_curr[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
- m_lnMnaughtSpecies[kspec]
+ m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
}
@ -2556,11 +2566,11 @@ namespace VCSnonideal {
* Recalculate the DeltaG's of the formation reactions for the
* deleted species in the mechanism
*/
vcs_deltag(0, true);
vcs_deltag(0, true, VCS_STATECALC_NEW);
for (iph = 0; iph < m_numPhases; iph++) {
if (m_tPhaseMoles_old[iph] > 0.0)
xtcutoff[iph] = log (m_tPhaseMoles_old[iph] / VCS_DELETE_SPECIES_CUTOFF);
xtcutoff[iph] = log (1.0 / VCS_RELDELETE_SPECIES_CUTOFF);
else
xtcutoff[iph] = 0.0;
}
@ -2706,7 +2716,7 @@ namespace VCSnonideal {
*/
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
iph = m_phaseID[kspec];
m_feSpecies_curr[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
- m_lnMnaughtSpecies[kspec]
+ m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]);
}
@ -2714,13 +2724,13 @@ namespace VCSnonideal {
* Recalculate the DeltaG's of the formation reactions for the
* deleted species in the mechanism
*/
vcs_deltag(0, true);
vcs_deltag(0, true, VCS_STATECALC_NEW);
for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] > 0.0) {
double maxDG = MIN(m_deltaGRxn_new[irxn], 300.0);
double maxDG = MIN(m_deltaGRxn_new[irxn], 690.0);
double dx = m_tPhaseMoles_old[iph] * exp(- maxDG);
retn = delta_species(kspec, &dx);
@ -2762,7 +2772,7 @@ namespace VCSnonideal {
}
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot);
vcs_deltag(0, true);
vcs_deltag(0, true, VCS_STATECALC_OLD);
retn = 0;
for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
@ -2770,15 +2780,18 @@ namespace VCSnonideal {
iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] > 0.0) {
if (fabs(m_deltaGRxn_old[irxn]) > m_tolmin) {
retn++;
if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) > VCS_DELETE_MINORSPECIES_CUTOFF) ||
(m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) {
retn++;
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g",
m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec],
m_deltaGRxn_old[irxn]);
plogendl();
}
if (m_debug_print_lvl >= 2) {
plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g",
m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec],
m_deltaGRxn_old[irxn]);
plogendl();
}
#endif
}
}
}
}
@ -2926,7 +2939,7 @@ namespace VCSnonideal {
* Evaluate DeltaG for all components if ITI=0, and for
* major components only if ITI NE 0
*/
vcs_deltag(0, false);
vcs_deltag(0, false, VCS_STATECALC_NEW);
dptr = VCS_DATA_PTR(m_deltaGRxn_new);
s2 = 0.0;
@ -3274,34 +3287,38 @@ namespace VCSnonideal {
}
/*****************************************************************************/
// This subroutine calculates reaction free energy changes for
// all noncomponent formation reactions.
/*
* Formation reactions are
* reactions which create each noncomponent species from the component
* species. m_stoichCoeffRxnMatrix[irxn][jcomp] are the stoichiometric
* coefficients for these reactions. A stoichiometric coefficient of
* one is assumed for species irxn in this reaction.
*
* vcs_deltag:
*
* This subroutine calculates reaction free energy changes for
* all noncomponent formation reactions. Formation reactions are
* reactions which create each noncomponent species from the component
* species. M_STOICHCOEFFRXNMATRIX(J,I) are the stoichiometric coefficients for these
* reactions. A stoichiometric coefficient of one is assumed for
* species I in this reaction.
*
* INPUT
* INPUT
* @param l
* L < 0 : Calculate reactions corresponding to
* major noncomponent and zeroed species only
* L = 0 : Do all noncomponent reactions, i, between
* 0 <= i < irxnl
* L > 0 : Calculate reactions corresponding to
* minor noncomponent and zeroed species only
* irxnl : used with L = 0 to indicate upper limit.
*
* Note we special case one important issue.
* If the component has zero moles, then we do not
* allow deltaG < 0.0 for formation reactions which
* would lead to the loss of more of the component.
* This dG < 0.0 feeds back into the algorithm in several
* places, and leads to a infinite loop in at least one case.
* @param doDeleted Do deleted species
* @param stateCalc Calculate deltaG corresponding to either old or new
* free energies
*
* Note we special case one important issue.
* If the component has zero moles, then we do not
* allow deltaG < 0.0 for formation reactions which
* would lead to the loss of more of that same component.
* This dG < 0.0 condition feeds back into the algorithm in several
* places, and leads to a infinite loop in at least one case.
*/
void VCS_SOLVE::vcs_deltag(int l, bool doDeleted) {
void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted,
const int vcsState) {
int iph;
int lneed, irxn, kspec;
double *dtmp_ptr;
@ -3310,7 +3327,6 @@ namespace VCSnonideal {
if (doDeleted) {
irxnl = m_numRxnTot;
}
int vcsState = VCS_STATECALC_NEW;
double * deltaGRxn;
double * feSpecies;
@ -3318,10 +3334,10 @@ namespace VCSnonideal {
double * actCoeffSpecies;
if (vcsState == VCS_STATECALC_NEW) {
deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new);
feSpecies = VCS_DATA_PTR(m_feSpecies_curr);
feSpecies = VCS_DATA_PTR(m_feSpecies_new);
molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_new);
actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_new);
} else if (vcsState == VCS_STATECALC_NEW) {
} else if (vcsState == VCS_STATECALC_OLD) {
deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old);
feSpecies = VCS_DATA_PTR(m_feSpecies_old);
molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old);
@ -5292,7 +5308,7 @@ namespace VCSnonideal {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
return dx;
}
double w_kspec = VCS_DELETE_SPECIES_CUTOFF;
double w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF;
#ifdef DEBUG_MODE
// Check to make sure that species is zero in the solution vector
// If it isn't, we don't know what's happening
@ -5324,7 +5340,7 @@ namespace VCSnonideal {
* There is no real way to estimate the moles. So
* we set it to a small number.
*/
dx = VCS_DELETE_SPECIES_CUTOFF * 100.;
dx = 1.0E-30;
}
/*