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

Working towards a more robust treatment of phase stability.

Browse source
This commit is contained in:
Harry Moffat 2008-08-17 21:19:04 +00:00
parent 3793f5dfac
commit e72d7a7cf6
11 changed files with 1084 additions and 804 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Cantera/src/equil/Makefile.in
View file

@ -70,7 +70,7 @@ VCSNONIDEAL_OBJ = vcs_solve_TP.o vcs_VolPhase.o vcs_solve.o vcs_prob.o \
vcs_inest.o vcs_rearrange.o \
vcs_root1d.o vcs_rxnadj.o \
vcs_SpeciesProperties.o vcs_equilibrate.o \
vcs_prep.o vcs_species_thermo.o vcs_Gibbs.o \
vcs_prep.o vcs_species_thermo.o vcs_Gibbs.o vcs_phaseStability.o \
$(DALT_OBJ)
VCSNONIDEAL_H = vcs_internal.h vcs_VolPhase.h vcs_solve.h vcs_prob.h \

72
Cantera/src/equil/vcs_VolPhase.cpp
View file

@ -45,7 +45,7 @@ namespace VCSnonideal {
m_numSpecies(0),
m_totalMolesInert(0.0),
m_isIdealSoln(false),
m_existence(0),
m_existence(VCS_PHASE_EXIST_NO),
m_MFStartIndex(0),
IndSpecies(0),
m_useCanteraCalls(false),
@ -512,6 +512,54 @@ namespace VCSnonideal {
}
/***************************************************************************/
// Set the moles and/or mole fractions within the phase
/*
*
*
*/
void vcs_VolPhase::setMoleFractionsState(double totalMoles,
const double * moleFractions,
const int vcsStateStatus) {
if (totalMoles != 0.0) {
if (vcsStateStatus != VCS_STATECALC_TMP) {
printf("vcs_VolPhase::setMolesFractionsState: inappropriate usage\n");
std::exit(-1);
}
m_UpToDate = false;
m_vcsStateStatus = VCS_STATECALC_TMP;
if (m_existence == -VCS_PHASE_EXIST_ZEROEDPHASE ) {
printf("vcs_VolPhase::setMolesFractionsState: inappropriate usage\n");
std::exit(-1);
}
m_existence = VCS_PHASE_EXIST_YES;
} else {
m_UpToDate = true;
m_vcsStateStatus = vcsStateStatus;
if (m_existence > VCS_PHASE_EXIST_NO ) {
m_existence = VCS_PHASE_EXIST_NO;
}
}
double sum = 0.0;
for (int k = 0; k < m_numSpecies; k++) {
Xmol[k] = moleFractions[k];
sum += moleFractions[k];
}
if (sum == 0.0) {
printf("vcs_VolPhase::setMolesFractionsState: inappropriate usage\n");
std::exit(-1);
}
if (sum != 1.0) {
for (int k = 0; k < m_numSpecies; k++) {
Xmol[k] /= sum;
}
}
_updateMoleFractionDependencies();
}
/***************************************************************************/
// Set the moles within the phase
/*
* This function takes as input the mole numbers in vcs format, and
@ -582,7 +630,7 @@ namespace VCSnonideal {
Xmol[k] = tmp / v_totalMoles;
}
}
m_existence = 1;
m_existence = VCS_PHASE_EXIST_YES;
} else {
// This is where we will start to store a better approximation
// for the mole fractions, when the phase doesn't exist.
@ -590,7 +638,7 @@ namespace VCSnonideal {
for (int k = 0; k < m_numSpecies; k++) {
Xmol[k] = 1.0 / m_numSpecies;
}
m_existence = 0;
m_existence = VCS_PHASE_EXIST_NO;
}
/*
* Update the electric potential if it is a solution variable
@ -606,12 +654,12 @@ namespace VCSnonideal {
double phi = molesSpeciesVCS[kglob];
setElectricPotential(phi);
if (m_numSpecies == 1) {
m_existence = 1;
m_existence = VCS_PHASE_EXIST_YES;
}
}
_updateMoleFractionDependencies();
if (m_totalMolesInert > 0.0) {
m_existence = 2;
m_existence = VCS_PHASE_EXIST_ALWAYS;
}
/*
* Set flags indicating we are up to date with the VCS state vector.
@ -1099,7 +1147,7 @@ namespace VCSnonideal {
void vcs_VolPhase::setTotalMoles(const double totalMols) {
v_totalMoles = totalMols;
if (m_totalMolesInert > 0.0) {
m_existence = 2;
m_existence = VCS_PHASE_EXIST_ALWAYS;
#ifdef DEBUG_MODE
if (totalMols < m_totalMolesInert) {
printf(" vcs_VolPhase::setTotalMoles:: ERROR totalMoles "
@ -1110,9 +1158,9 @@ namespace VCSnonideal {
#endif
} else {
if (totalMols > 0.0) {
m_existence = 1;
m_existence = VCS_PHASE_EXIST_YES;
} else {
m_existence = 0;
m_existence = VCS_PHASE_EXIST_NO;
}
}
}
@ -1226,7 +1274,7 @@ namespace VCSnonideal {
// Set the existence flag in the object
void vcs_VolPhase::setExistence(const int existence) {
if (existence == 0) {
if (existence == VCS_PHASE_EXIST_NO || existence == VCS_PHASE_EXIST_ZEROEDPHASE) {
if (v_totalMoles != 0.0) {
#ifdef DEBUG_MODE
plogf("vcs_VolPhase::setExistence setting false existence for phase with moles");
@ -1295,12 +1343,12 @@ namespace VCSnonideal {
m_totalMolesInert = tMolesInert;
}
if (m_totalMolesInert > 0.0) {
m_existence = 2;
m_existence = VCS_PHASE_EXIST_ALWAYS;
} else {
if (v_totalMoles > 0.0) {
m_existence = 1;
m_existence = VCS_PHASE_EXIST_YES;
} else {
m_existence = 0;
m_existence = VCS_PHASE_EXIST_NO;
}
}
}

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

@ -157,6 +157,9 @@ namespace VCSnonideal {
* @param kspec species number
*/
double AC_calc_one(int kspec) const;
void setMoleFractionsState(const double molNum, const double * const moleFracVec,
const int vcsStateStatus);
//! Set the moles within the phase
/*!
@ -391,6 +394,7 @@ namespace VCSnonideal {
*/
void setMolesCurrent(int stateCalc);
private:
//! Set the mole fractions from a conventional mole fraction vector
/*!
*
@ -399,6 +403,7 @@ namespace VCSnonideal {
*/
void setMoleFractions(const double * const xmol);
public:
//! Return a const reference to the mole fractions
const std::vector<double> & moleFractions() const;
@ -768,9 +773,13 @@ namespace VCSnonideal {
//! Current state of existence:
/*!
* 0 : Doesn't exist currently
* 1 : Does exist currently
* 2 : Always exists because it contains
* 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
*/
@ -885,6 +894,7 @@ namespace VCSnonideal {
* valid values are
* VCS_STATECALC_OLD
* VCS_STATECALC_NEW
* VCS_STATECALC_TMP
*/
int m_vcsStateStatus;

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

@ -143,19 +143,18 @@ namespace VCSnonideal {
*/
#define VCS_SPECIES_MINOR 0
//! Species lies in a multicomponent phase that is zeroed atm
//! Species lies in a multicomponent phase, with a small phase concentration
/*!
* The species lies in a multicomponent phase that is currently
* deleted.
* The species lies in a multicomponent phase that exists.
* It concentration is currently very low, necessitating a
* different method of calculation.
*/
#define VCS_SPECIES_ZEROEDPHASE -1
#define VCS_SPECIES_SMALLMS -1
//! Species lies in a multicomponent phase, with concentration zero
/*!
* The species lies in a multicomponent phase that exists.
* It concentration is currently zero, even though it may
* or may not actually have a low mole fraction in the phase
* this situation occurs when phases pop back into life.
* The species lies in a multicomponent phase which currently doesn't exist.
* It concentration is currently zero.
*/
#define VCS_SPECIES_ZEROEDMS -2
@ -166,11 +165,12 @@ namespace VCSnonideal {
*/
#define VCS_SPECIES_ZEROEDSS -3
//! Species has such a small mole fraction it is deleted.
//! 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 aded back in at the end of the calculation.
* It will be added back in at the end of the calculation.
*/
#define VCS_SPECIES_DELETED -4
@ -182,8 +182,54 @@ namespace VCSnonideal {
*/
#define VCS_SPECIES_INTERFACIALVOLTAGE -5
//! Species lies in a multicomponent phase that is zeroed atm
/*!
* The species lies in a multicomponent phase that is currently
* deleted 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.
*/
#define VCS_SPECIES_ZEROEDPHASE -6
//! Species lies in a multicomponent phase that is active, but its concentration is zero
/*!
* The species lies in a multicomponent phase which currently does exist.
* It concentration is currently 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".
*/
#define VCS_SPECIES_ACTIVEBUTZERO -7
//@}
//! @name Phase Categories used during the iteration
/*!
* These defines are valid values for the phase existence flag
*/
//@{
//! Always exists because it contains inerts which can't exist in any other phase
#define VCS_PHASE_EXIST_ALWAYS 3
//! Phase is a normal phase that currently exists
#define VCS_PHASE_EXIST_YES 2
//! Phase is a normal phase that exists in a small concentration
/*!
* Concentration is so small that it must be calculated using an alternate
* method
*/
#define VCS_PHASE_EXIST_MINORCONC 1
//! Phase doesn't currently exist in the mixture
#define VCS_PHASE_EXIST_NO 0
//! Phase currently is zeroed due to a programmatic issue
/*!
* We zero phases because we want to follow phase stability boundaries.
*/
#define VCS_PHASE_EXIST_ZEROEDPHASE -6
//@}
/*!
* @name Units for the chemical potential data and pressure variables

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

@ -428,6 +428,8 @@ namespace VCSnonideal {
const std::vector<int> vec_from, const int len = -1);
#endif
double vcs_l2norm(const std::vector<double> vec);
//! Finds the location of the maximum component in a double vector
/*!
* @param x pointer to a vector of doubles

116
Cantera/src/equil/vcs_phaseStability.cpp Normal file
View file

@ -0,0 +1,116 @@
/* ======================================================================= */
/* -------------------------------------------------- */
/* | RCS Head Information on zuzax.pchem.sandia.gov | */
/* -------------------------------------------------- */
/* $RCSfile$ */
/* $Author$ */
/* $Date$ */
/* $Revision$ */
/* ======================================================================= */
#include "vcs_solve.h"
#include "vcs_internal.h"
#include "vcs_species_thermo.h"
#include "vcs_VolPhase.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <vector>
using namespace std;
namespace VCSnonideal {
int VCS_SOLVE::vcs_phaseStabilityTest(const int iph) {
/*
* We will use the _new state calc here
*/
int kspec, irxn, k;
vcs_VolPhase *Vphase = m_VolPhaseList[iph];
double deltaGRxn;
// We will do a full newton calculation later, but for now, ...
bool doSuccessiveSubstitution = true;
int res = 0;
vector<double> X_est(Vphase->nSpecies(), 0.0);
vector<double> X_est_old(Vphase->nSpecies(), 0.0);
vector<double> delX(Vphase->nSpecies(), 0.0);
vector<double> E_phi(Vphase->nSpecies(), 0.0);
double damp = 1.0;
double normUpdate = 1.0;
double normUpdateOld = 1.0;
// get the activity coefficients
Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new));
if (doSuccessiveSubstitution) {
for (int its = 0; its < 20; its++) {
normUpdateOld = normUpdate;
for (k = 0; k < Vphase->nSpecies(); k++) {
X_est_old[k] = X_est[k];
}
double poly = -1.0;
for (k = 0; k < Vphase->nSpecies(); k++) {
kspec = Vphase->spGlobalIndexVCS(k);
irxn = kspec - m_numComponents;
deltaGRxn = m_deltaGRxn_old[irxn];
// We may need to look at deltaGRxn for components!
if (irxn >= 0) {
if (deltaGRxn > 50.0) deltaGRxn = 50.0;
if (deltaGRxn < -50.0) deltaGRxn = -50.0;
E_phi[k] = exp(-deltaGRxn)/m_actCoeffSpecies_new[kspec];
poly += E_phi[k];
}
}
double sum = poly + 1.0;
for (k = 0; k < Vphase->nSpecies(); k++) {
delX[k] = E_phi[k]/sum - X_est_old[k];
}
normUpdate = vcs_l2norm(delX);
// Figure out the damping coefficient
double ratio = normUpdate / normUpdateOld;
if (ratio < 0.4) {
damp = 1.0;
} else if (ratio > 1.0) {
damp = 0.03;
} else {
damp = 0.1;
}
for (k = 0; k < Vphase->nSpecies(); k++) {
X_est[k] = X_est_old[k] + damp * delX[k];
}
for (k = 0; k < Vphase->nSpecies(); k++) {
kspec = Vphase->spGlobalIndexVCS(k);
m_molNumSpecies_new[kspec] = X_est[k];
}
Vphase->setMolesFromVCS(VCS_STATECALC_NEW);
}
} else {
printf("not done yet\n");
exit(-1);
}
return res;
}
}

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

@ -21,6 +21,338 @@
namespace VCSnonideal {
// Calculates formation reaction step sizes.
/*
* This is equation 6.4-16, p. 143 in Smith and Missen.
*
* Output
* -------
* m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers
* to the irxn'th species
* formation reaction. This adjustment is for species
* irxn + M, where M is the number of components.
*
* Special branching occurs sometimes. This causes the component basis
* to be reevaluated
*
* @return Returns an int representing the status of the step
* - 0 : normal return
* - 1 : A single species phase species has been zeroed out
* in this routine. The species is a noncomponent
* - 2 : Same as one but, the zeroed species is a component.
*/
int VCS_SOLVE::vcs_RxnStepSizes() {
int j, irxn, kspec, soldel = 0, iph;
double s, xx, dss;
int k = 0;
vcs_VolPhase *Vphase = 0;
double *dnPhase_irxn;
#ifdef DEBUG_MODE
char ANOTE[128];
if (m_debug_print_lvl >= 2) {
plogf(" "); for (j = 0; j < 82; j++) plogf("-"); plogf("\n");
plogf(" --- Subroutine vcs_RxnStepSizes called - Details:\n");
plogf(" "); for (j = 0; j < 82; j++) plogf("-"); plogf("\n");
plogf(" --- Species KMoles Rxn_Adjustment DeltaG"
" | Comment\n");
}
#endif
/*
* We update the matrix dlnActCoeffdmolNumber[][] at the
* top of the loop, when necessary
*/
if (m_useActCoeffJac) {
vcs_CalcLnActCoeffJac(VCS_DATA_PTR(m_molNumSpecies_old));
}
/************************************************************************
******** LOOP OVER THE FORMATION REACTIONS *****************************
************************************************************************/
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Normal Calc");
#endif
kspec = m_indexRxnToSpecies[irxn];
if (m_rxnStatus[irxn] == VCS_SPECIES_ZEROEDPHASE) {
m_deltaMolNumSpecies[kspec] = 0.0;
#ifdef DEBUG_MODE
sprintf(ANOTE, "ZeroedPhase: Phase is artificially zeroed");
#endif
} else if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
dnPhase_irxn = m_deltaMolNumPhase[irxn];
if (m_molNumSpecies_old[kspec] == 0.0 && (! m_SSPhase[kspec])) {
/********************************************************************/
/******* MULTISPECIES PHASE WITH total moles equal to zero *********/
/*******************************************************************/
/*
* If dg[irxn] is negative, then the multispecies phase should
* come alive again. Add a small positive step size to
* make it come alive.
*/
if (m_deltaGRxn_new[irxn] < -1.0e-4) {
/*
* First decide if this species is part of a multiphase that
* is nontrivial in size.
*/
iph = m_phaseID[kspec];
double tphmoles = m_tPhaseMoles_old[iph];
double trphmoles = tphmoles / m_totalMolNum;
if (trphmoles > VCS_DELETE_PHASE_CUTOFF) {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
#ifdef DEBUG_MODE
sprintf(ANOTE,
"MultSpec: small species born again DG = %11.3E",
m_deltaGRxn_new[irxn]);
#endif
} else {
#ifdef DEBUG_MODE
sprintf(ANOTE, "MultSpec: phase come alive DG = %11.3E",
m_deltaGRxn_new[irxn]);
#endif
Vphase = m_VolPhaseList[iph];
int numSpPhase = Vphase->nSpecies();
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]);
#endif
m_deltaMolNumSpecies[kspec] = 0.0;
}
} else {
/********************************************************************/
/************************* REGULAR PROCESSING ************/
/********************************************************************/
/*
* First take care of cases where we want to bail out
*
*
* Don't bother if superconvergence has already been achieved
* in this mode.
*/
if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]);
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
#endif
continue;
}
/*
* Don't calculate for minor or nonexistent species if
* their values are to be decreasing anyway.
*/
if ((m_rxnStatus[irxn] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E",
m_rxnStatus[irxn], m_deltaGRxn_new[irxn]);
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
#endif
continue;
}
/*
* Start of the regular processing
*/
if (m_SSPhase[kspec]) {
s = 0.0;
} else {
s = 1.0 / m_molNumSpecies_old[kspec] ;
}
for (j = 0; j < m_numComponents; ++j) {
if (!m_SSPhase[j]) {
if (m_molNumSpecies_old[j] > 0.0) {
s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
}
}
}
for (j = 0; j < m_numPhases; j++) {
Vphase = m_VolPhaseList[j];
if (! Vphase->m_singleSpecies) {
if (m_tPhaseMoles_old[j] > 0.0)
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
}
}
if (s != 0.0) {
/*
* Take into account of the
* derivatives of the activity coefficients with respect to the
* mole numbers, even in our diagonal approximation.
*/
if (m_useActCoeffJac) {
double s_old = s;
s = vcs_Hessian_diag_adj(irxn, s_old);
#ifdef DEBUG_MODE
if (s_old != s) {
sprintf(ANOTE, "Normal calc: diag adjusted from %g "
"to %g due to act coeff", s_old, s);
}
#endif
}
m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
// New section to do damping of the m_deltaMolNumSpecies[]
/*
*
*/
for (j = 0; j < m_numComponents; ++j) {
double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
if (stoicC != 0.0) {
double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
if (negChangeComp > m_molNumSpecies_old[j]) {
if (m_molNumSpecies_old[j] > 0.0) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %d (%10s) going neg", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
#endif
m_deltaMolNumSpecies[kspec] = - m_molNumSpecies_old[j] / stoicC;
} else {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %d (%10s) zero", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
#endif
m_deltaMolNumSpecies[kspec] = 0.0;
}
}
}
}
// Implement a damping term that limits m_deltaMolNumSpecies to the size of the mole number
if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to %s going negative", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
#endif
m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
}
} else {
/* ************************************************************ */
/* **** REACTION IS ENTIRELY AMONGST SINGLE SPECIES PHASES **** */
/* **** DELETE ONE OF THE PHASES AND RECOMPUTE BASIS ********* */
/* ************************************************************ */
/*
* Either the species L will disappear or one of the
* component single species phases will disappear. The sign
* of DG(I) will indicate which way the reaction will go.
* Then, we need to follow the reaction to see which species
* will zero out first.
* -> The species to be zeroed out will be "k".
*/
if (m_deltaGRxn_new[irxn] > 0.0) {
dss = m_molNumSpecies_old[kspec];
k = kspec;
for (j = 0; j < m_numComponents; ++j) {
if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) {
xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
}
}
}
dss = -dss;
} else {
dss = 1.0e10;
for (j = 0; j < m_numComponents; ++j) {
if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) {
xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
}
}
}
}
/*
* Here we adjust the mole fractions
* according to DSS and the stoichiometric array
* to take into account that we are eliminating
* the kth species. DSS contains the amount
* of moles of the kth species that needs to be
* added back into the component species.
*/
if (dss != 0.0) {
m_molNumSpecies_old[kspec] += dss;
m_tPhaseMoles_old[m_phaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j];
m_tPhaseMoles_old[m_phaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
}
m_molNumSpecies_old[k] = 0.0;
iph = m_phaseID[k];
m_tPhaseMoles_old[iph] = 0.0;
Vphase = m_VolPhaseList[iph];
Vphase->setTotalMoles(0.0);
if (k == kspec) {
m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDSS;
if (m_SSPhase[kspec] != 1) {
printf("we shouldn't be here!\n");
exit(-1);
}
}
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_RxnStepSizes Special section to delete %s",
m_speciesName[k].c_str());
plogendl();
}
#endif
/*
* We need to immediately recompute the
* component basis, because we just zeroed
* it out.
*/
soldel = 1;
if (k != kspec) {
soldel = 2;
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- Immediate return to get new basis - Restart iteration\n");
plogendl();
}
#endif
return soldel;
}
}
}
} /* End of regular processing */
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E %12.4E | %s\n",
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
#endif
} /* End of loop over m_speciesUnknownType */
} /* End of loop over non-component stoichiometric formation reactions */
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" "); vcs_print_line("-", 82);
}
#endif
return soldel;
}
/*****************************************************************************/
//! Calculates reaction adjustments using a full Hessian approximation
/*!
* Calculates reaction adjustments. This does what equation 6.4-16, p. 143

9
Cantera/src/equil/vcs_solve.cpp
View file

@ -948,7 +948,7 @@ namespace VCSnonideal {
pub->w[i] = 0.0;
plogf("voltage species = %g\n", m_molNumSpecies_old[k1]);
}
pub->mf[i] = m_molNumSpecies_new[k1];
//pub->mf[i] = m_molNumSpecies_new[k1];
pub->m_gibbsSpecies[i] = m_feSpecies_old[k1];
pub->VolPM[i] = m_PMVolumeSpecies[k1];
}
@ -964,10 +964,13 @@ namespace VCSnonideal {
pubPhase->setTotalMoles(vPhase->TotalMoles());
pubPhase->setElectricPotential(vPhase->electricPotential());
double sumMoles = pubPhase->totalMolesInert();
pubPhase->setMoleFractions(VCS_DATA_PTR(vPhase->moleFractions()));
pubPhase->setMoleFractionsState(vPhase->TotalMoles(),
VCS_DATA_PTR(vPhase->moleFractions()),
VCS_STATECALC_TMP);
const std::vector<double> & mfVector = pubPhase->moleFractions();
for (int k = 0; k < pubPhase->nSpecies(); k++) {
kT = pubPhase->spGlobalIndexVCS(k);
pub->mf[kT] = mfVector[k];
if (pubPhase->phiVarIndex() == k) {
k1 = vPhase->spGlobalIndexVCS(k);
double tmp = m_molNumSpecies_old[k1];

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

@ -528,8 +528,10 @@ public:
* minor noncomponent and zeroed species only
*
* @param doDeleted Do deleted species
* @param stateCalc Calculate deltaG corresponding to either old or new
* @param vcsState Calculate deltaG corresponding to either old or new
* free energies
* @param alterZeroedPhases boolean indicating whether we should
* add in a special section for zeroed phases.
*
* Note we special case one important issue.
* If the component has zero moles, then we do not
@ -538,7 +540,30 @@ public:
* 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);
void vcs_deltag(const int l, const bool doDeleted, const int vcsState,
const bool alterZeroedPhases = true);
//! Calculate deltag of formation for all species in a single phase.
/*!
* Calculate deltag of formation for all species in a single
* phase. It is assumed that the fe[] is up to date for all species.
* Howevever, if the phase is currently zereoed out, a subproblem
* is calculated to solve for AC[i] and pseudo-X[i] for that
* phase.
*
* @param iphase phase index of the phase to be calculated
* @param doDeleted boolean indicating whether to do deleted
* species or not
* @param stateCalc integer describing which set of free energies
* to use and where to stick the results.
* @param alterZeroedPhases boolean indicating whether we should
* add in a special section for zeroed phases.
*
* NOTE: this is currently not used used anywhere.
* It may be in the future?
*/
void vcs_deltag_Phase(const int iphase, const bool doDeleted,
const int stateCalc, const bool alterZeroedPhases = true);
//! Swaps the indecises for all of the global data for two species, k1
//! and k2.
@ -558,26 +583,6 @@ public:
void vcs_switch_pos(const int ifunc, const int k1, const int k2);
//! Calculate deltag of formation for all species in a single phase.
/*!
* Calculate deltag of formation for all species in a single
* phase. It is assumed that the fe[] is up to date for all species.
* Howevever, if the phase is currently zereoed out, a subproblem
* is calculated to solve for AC[i] and pseudo-X[i] for that
* phase.
*
* @param iphase phase index of the phase to be calculated
* @param doDeleted boolean indicating whether to do deleted
* species or not
* @param stateCalc integer describing which set of free energies
* to use and where to stick the results.
*
* NOTE: this is currently not used used anywhere.
* It may be in the future?
*/
void vcs_deltag_Phase(const int iphase, const bool doDeleted,
const int stateCalc);
//! Birth guess returns the number of moles of a species
//! that is coming back to life.
/*!
@ -603,6 +608,10 @@ public:
*/
double vcs_birthGuess(const int kspec);
int vcs_phaseStabilityTest(const int iph);
//! Solve an equilibrium problem at a particular fixed temperature
//! and pressure
/*!

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

File diff suppressed because it is too large Load diff

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

@ -110,6 +110,19 @@ namespace VCSnonideal {
}
#endif
double vcs_l2norm(const std::vector<double> vec) {
int len = vec.size();
if (len == 0) {
return 0.0;
}
double sum = 0.0;
std::vector<double>::const_iterator pos;
for (pos = vec.begin(); pos != vec.end(); ++pos) {
sum += (*pos) * (*pos);
}
return std::sqrt(sum/len);
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/