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

More variable and documentation checks.

Browse source 
Added the VCS_STATECALC_ defines. Will try to clean up the evaluation
of phase properties using this concept.
This commit is contained in:
Harry Moffat 2008-05-14 16:50:34 +00:00
parent 6f07cb109e
commit f741f96b3c
7 changed files with 706 additions and 568 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -271,7 +271,23 @@ namespace VCSnonideal {
* is ddefined by the interface voltage.
*/
#define VCS_SPECIES_TYPE_INTERFACIALVOLTAGE -5
//@}
/*!
* @name Types of State Calculations within VCS
* These values determine where the
* results are storred within the VCS_SOLVE
* object.
* @{
*/
//! State Calculation based on the old or base mole numbers
#define VCS_STATECALC_OLD 0
//! State Calculation based on the new or tentative mole numbers
#define VCS_STATECALC_NEW 1
//! State Calculation based on a temporary set of mole numbers
#define VCS_STATECALC_TMP 2
//@}

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

@ -323,7 +323,7 @@ namespace VCSnonideal {
/* ******************************************* */
/* **** CONVERGENCE FORCING SECTION ********** */
/* ******************************************* */
vcs_dfe(molNum, 0, 0, 0, nspecies);
vcs_dfe(molNum, VCS_STATECALC_OLD, 0, 0, nspecies);
for (kspec = 0, s = 0.0; kspec < nspecies; ++kspec) {
s += m_deltaMolNumSpecies[kspec] * m_feSpecies_curr[kspec];
}

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

@ -318,7 +318,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("%14.7E ", log(ActCoeff[l]));
double tpmoles = m_tPhaseMoles_old[pid];
double phi = m_phasePhi[pid];
double eContrib = phi * Charge[l] * Faraday_dim;
double eContrib = phi * m_chargeSpecies[l] * Faraday_dim;
double lx = 0.0;
if (m_speciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
lx = 0.0;

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

@ -1,6 +1,6 @@
/**
* @file vcs_rxnadj.cpp
* routines for carrying out various line adjustments
* Routines for carrying out various adjustments to the reaction steps
*/
/*
* $Id$
@ -10,533 +10,535 @@
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "vcs_solve.h"
#include "vcs_internal.h"
#include "vcs_VolPhase.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
namespace VCSnonideal {
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
int VCS_SOLVE::vcs_rxn_adj_cg(void)
/**************************************************************************
*
* vcs_rxn_adj_cg:
*
* Calculates reaction adjustments. This does what equation 6.4-16, p. 143
//! Calculates reaction adjustments using a full Hessian approximation
/*!
* Calculates reaction adjustments. This does what equation 6.4-16, p. 143
* in Smith and Missen is suppose to do. However, a full matrix is
* formed and then solved via a conjugate gradient algorithm. No
* formed and then solved via a conjugate gradient algorithm. No
* preconditioning is done.
*
* If special branching is warranted, then the program bails out.
*
* Output
* -------
* Output
* -------
* DS(I) : reaction adjustment, where I refers to the Ith species
* Special branching occurs sometimes. This causes the component basis
* to be reevaluated
* Special branching occurs sometimes. This causes the component basis
* to be reevaluated
* return = 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.
* in this routine. The species is a noncomponent
* 2 : Same as one but, the zeroed species is a component.
*
* Special attention is taken to flag cases where the direction of the
* update is contrary to the steepest descent rule. This is an important
* attribute of the regular vcs algorithm. We don't want to violate this
***************************************************************************/
{
int irxn, j;
int k = 0;
int kspec, soldel = 0;
double s, xx, dss;
double *dnPhase_irxn;
* attribute of the regular vcs algorithm. We don't want to violate this.
*
* NOTE: currently this routine is not used.
*/
int VCS_SOLVE::vcs_rxn_adj_cg() {
int irxn, j;
int k = 0;
int kspec, soldel = 0;
double s, xx, dss;
double *dnPhase_irxn;
#ifdef DEBUG_MODE
char ANOTE[128];
plogf(" "); for (j = 0; j < 77; j++) plogf("-");
plogf("\n --- Subroutine rxn_adj_cg() called\n");
plogf(" --- Species Moles Rxn_Adjustment | Comment\n");
char ANOTE[128];
plogf(" "); for (j = 0; j < 77; j++) plogf("-");
plogf("\n --- Subroutine rxn_adj_cg() called\n");
plogf(" --- Species Moles Rxn_Adjustment | Comment\n");
#endif
/*
* Precalculation loop -> we calculate quantities based on
* loops over the number of species.
* We also evaluate whether the matrix is appropriate for
* this algorithm. If not, we bail out.
*/
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
/*
* Precalculation loop -> we calculate quantities based on
* loops over the number of species.
* We also evaluate whether the matrix is appropriate for
* this algorithm. If not, we bail out.
*/
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Normal Calc");
sprintf(ANOTE,"Normal Calc");
#endif
kspec = ir[irxn];
dnPhase_irxn = m_deltaMolNumPhase[irxn];
kspec = ir[irxn];
dnPhase_irxn = m_deltaMolNumPhase[irxn];
if (m_molNumSpecies_old[kspec] == 0.0 && (! SSPhase[kspec])) {
/* *******************************************************************/
/* **** MULTISPECIES PHASE WITH total moles equal to zero ************/
/* *******************************************************************/
/*
* HKM -> the statment below presupposes units in m_deltaGRxn_new[]. It probably
* should be replaced with something more relativistic
*/
if (m_deltaGRxn_new[irxn] < -1.0e-4) {
#ifdef DEBUG_MODE
(void) sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]);
#endif
m_deltaMolNumSpecies[kspec] = 1.0e-10;
spStatus[irxn] = VCS_SPECIES_MAJOR;
--(m_numRxnMinorZeroed);
} else {
#ifdef DEBUG_MODE
(void) 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: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]);
plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], ANOTE);
#endif
continue;
}
/*
* Don't calculate for minor or nonexistent species if
* their values are to be decreasing anyway.
*/
if (spStatus[irxn] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E\n",
spStatus[irxn], m_deltaGRxn_new[irxn]);
plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
#endif
continue;
}
/*
* Start of the regular processing
*/
if (SSPhase[kspec]) s = 0.0;
else s = 1.0 / m_molNumSpecies_old[kspec];
for (j = 0; j < m_numComponents; ++j) {
if (! SSPhase[j]) s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
}
for (j = 0; j < m_numPhases; j++) {
if (! (VPhaseList[j])->SingleSpecies) {
if (m_tPhaseMoles_old[j] > 0.0)
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
}
}
if (s != 0.0) {
m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
} else {
/* ************************************************************ */
/* **** REACTION IS ENTIRELY AMONGST SINGLE SPECIES PHASES **** */
/* **** DELETE ONE SOLID 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.
if (m_molNumSpecies_old[kspec] == 0.0 && (! SSPhase[kspec])) {
/* *******************************************************************/
/* **** MULTISPECIES PHASE WITH total moles equal to zero ************/
/* *******************************************************************/
/*
* HKM -> the statment below presupposes units in m_deltaGRxn_new[]. It probably
* should be replaced with something more relativistic
*/
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;
if (m_deltaGRxn_new[irxn] < -1.0e-4) {
#ifdef DEBUG_MODE
(void) sprintf(ANOTE, "MultSpec: come alive DG = %11.3E", m_deltaGRxn_new[irxn]);
#endif
m_deltaMolNumSpecies[kspec] = 1.0e-10;
spStatus[irxn] = VCS_SPECIES_MAJOR;
--(m_numRxnMinorZeroed);
} 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;
}
}
}
#ifdef DEBUG_MODE
(void) 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: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]);
plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], ANOTE);
#endif
continue;
}
/*
* 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.
* Don't calculate for minor or nonexistent species if
* their values are to be decreasing anyway.
*/
if (dss != 0.0) {
m_molNumSpecies_old[kspec] += dss;
m_tPhaseMoles_old[PhaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j];
m_tPhaseMoles_old[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
}
m_molNumSpecies_old[k] = 0.0;
m_tPhaseMoles_old[PhaseID[k]] = 0.0;
if (spStatus[irxn] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) {
#ifdef DEBUG_MODE
plogf(" --- vcs_st2 Special section to delete ");
plogf("%-12.12s", m_speciesName[k].c_str());
plogf("\n --- Immediate return - Restart iteration\n");
#endif
/*
* We need to immediately recompute the
* component basis, because we just zeroed
* it out.
*/
if (k != kspec) soldel = 2;
else soldel = 1;
return soldel;
sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E\n",
spStatus[irxn], m_deltaGRxn_new[irxn]);
plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
#endif
continue;
}
}
} /* End of regular processing */
/*
* Start of the regular processing
*/
if (SSPhase[kspec]) s = 0.0;
else s = 1.0 / m_molNumSpecies_old[kspec];
for (j = 0; j < m_numComponents; ++j) {
if (! SSPhase[j]) s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
}
for (j = 0; j < m_numPhases; j++) {
if (! (VPhaseList[j])->SingleSpecies) {
if (m_tPhaseMoles_old[j] > 0.0)
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
}
}
if (s != 0.0) {
m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s;
} else {
/* ************************************************************ */
/* **** REACTION IS ENTIRELY AMONGST SINGLE SPECIES PHASES **** */
/* **** DELETE ONE SOLID 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.
*/
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[PhaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j];
m_tPhaseMoles_old[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
}
m_molNumSpecies_old[k] = 0.0;
m_tPhaseMoles_old[PhaseID[k]] = 0.0;
#ifdef DEBUG_MODE
plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
#endif
} /* End of loop over non-component stoichiometric formation reactions */
/*
*
* When we form the Hessian we must be careful to ensure that it
* is a symmetric positive definate matrix, still. This means zeroing
* out columns when we zero out rows as well.
* -> I suggest writing a small program to make sure of this
* property.
*/
#ifdef DEBUG_MODE
plogf(" "); for (j = 0; j < 77; j++) plogf("-"); plogf("\n");
plogf(" --- vcs_st2 Special section to delete ");
plogf("%-12.12s", m_speciesName[k].c_str());
plogf("\n --- Immediate return - Restart iteration\n");
#endif
return soldel;
} /* vcs_rxn_adj_cg() ********************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
double VCS_SOLVE::vcs_Hessian_diag_adj(int irxn, double hessianDiag_Ideal)
/*
* We need to immediately recompute the
* component basis, because we just zeroed
* it out.
*/
if (k != kspec) soldel = 2;
else soldel = 1;
return soldel;
}
}
} /* End of regular processing */
#ifdef DEBUG_MODE
plogf(" --- "); plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
#endif
} /* End of loop over non-component stoichiometric formation reactions */
/**************************************************************************
*
* vcs_actCoeff_diag_adj(irxn):
*
* Calculates the diagonal contribution to the Hessian due to
* the dependence of the activity coefficients on the mole numbers.
*
/*
*
* When we form the Hessian we must be careful to ensure that it
* is a symmetric positive definate matrix, still. This means zeroing
* out columns when we zero out rows as well.
* -> I suggest writing a small program to make sure of this
* property.
*/
#ifdef DEBUG_MODE
plogf(" "); for (j = 0; j < 77; j++) plogf("-"); plogf("\n");
#endif
return soldel;
}
/*****************************************************************************/
// Calculates the diagonal contribution to the Hessian due to
// the dependence of the activity coefficients on the mole numbers.
/*
* (See framemaker notes, Eqn. 20 - VCS Equations document)
*
* We allow the diagonal to be increased positively to any degree.
* We allow the diagonal to be decreased to 1/3 of the ideal solution
* value, but no more -> it must remain positive.
**************************************************************************/
{
double diag = hessianDiag_Ideal;
double hessActCoef = vcs_Hessian_actCoeff_diag(irxn);
if (hessianDiag_Ideal <= 0.0) {
plogf("We shouldn't be here\n");
exit(-1);
*
* NOTE: currently this routine is not used
*/
double VCS_SOLVE::vcs_Hessian_diag_adj(int irxn, double hessianDiag_Ideal) {
double diag = hessianDiag_Ideal;
double hessActCoef = vcs_Hessian_actCoeff_diag(irxn);
if (hessianDiag_Ideal <= 0.0) {
plogf("We shouldn't be here\n");
exit(-1);
}
if (hessActCoef >= 0.0) {
diag += hessActCoef;
} else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) {
diag += hessActCoef;
} else {
diag -= 0.6666 * hessianDiag_Ideal;
}
return diag;
}
if (hessActCoef >= 0.0) {
diag += hessActCoef;
} else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) {
diag += hessActCoef;
} else {
diag -= 0.6666 * hessianDiag_Ideal;
}
return diag;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
double VCS_SOLVE::vcs_Hessian_actCoeff_diag(int irxn)
/**************************************************************************
*
* vcs_Hessian_actCoeff_diag(irxn):
*
* Calculates the diagonal contribution to the Hessian due to
* the dependence of the activity coefficients on the mole numbers.
//! Calculates the diagonal contribution to the Hessian due to
//! the dependence of the activity coefficients on the mole numbers.
/*!
* (See framemaker notes, Eqn. 20 - VCS Equations document)
**************************************************************************/
{
int kspec, k, l, kph;
double s;
double *sc_irxn;
kspec = ir[irxn];
kph = PhaseID[kspec];
sc_irxn = m_stoichCoeffRxnMatrix[irxn];
/*
* First the diagonal term of the Jacobian
*
* NOTE: currently this routine is not used
*/
s = dLnActCoeffdMolNum[kspec][kspec];
/*
* Next, the other terms. Note this only a loop over the components
* So, it's not too expensive to calculate.
*/
for (l = 0; l < m_numComponents; l++) {
if (!SSPhase[l]) {
for (k = 0; k < m_numComponents; ++k) {
if (PhaseID[k] == PhaseID[l]) {
s += sc_irxn[k] * sc_irxn[l] * dLnActCoeffdMolNum[k][l];
double VCS_SOLVE::vcs_Hessian_actCoeff_diag(int irxn)
{
int kspec, k, l, kph;
double s;
double *sc_irxn;
kspec = ir[irxn];
kph = PhaseID[kspec];
sc_irxn = m_stoichCoeffRxnMatrix[irxn];
/*
* First the diagonal term of the Jacobian
*/
s = dLnActCoeffdMolNum[kspec][kspec];
/*
* Next, the other terms. Note this only a loop over the components
* So, it's not too expensive to calculate.
*/
for (l = 0; l < m_numComponents; l++) {
if (!SSPhase[l]) {
for (k = 0; k < m_numComponents; ++k) {
if (PhaseID[k] == PhaseID[l]) {
s += sc_irxn[k] * sc_irxn[l] * dLnActCoeffdMolNum[k][l];
}
}
if (kph == PhaseID[l]) {
s += sc_irxn[l] * (dLnActCoeffdMolNum[kspec][l] + dLnActCoeffdMolNum[l][kspec]);
}
}
if (kph == PhaseID[l]) {
s += sc_irxn[l] * (dLnActCoeffdMolNum[kspec][l] + dLnActCoeffdMolNum[l][kspec]);
}
return s;
}
/*****************************************************************************/
//! Recalculate all of the activity coefficients in all of the phases
//! based on input mole numbers
/*!
*
* @param moleSpeciesVCS kmol of species to be used in the update.
*
* NOTE: This routine needs to be regulated.
*/
void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double * const moleSpeciesVCS) {
/*
* Loop over all of the phases in the problem
*/
for (int iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase *Vphase = VPhaseList[iphase];
/*
* We don't need to call single species phases;
*/
if (!Vphase->SingleSpecies) {
/*
* update the Ln Act Coeff jacobian entries with respect to the
* mole number of species in the phase
*/
Vphase->updateLnActCoeffJac(moleSpeciesVCS);
/*
* Download the resulting calculation into the full vector
* -> This scatter calculation is carried out in the
* vcs_VolPhase object.
*/
Vphase->sendToVCSLnActCoeffJac(dLnActCoeffdMolNum.baseDataAddr());
}
}
}
return s;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double * const moleSpeciesVCS)
/*************************************************************************
*
*
*
*
*************************************************************************/
{
/*
* Loop over all of the phases in the problem
*/
for (int iphase = 0; iphase < m_numPhases; iphase++) {
vcs_VolPhase *Vphase = VPhaseList[iphase];
/*
* We don't need to call single species phases;
*/
if (!Vphase->SingleSpecies) {
/*
* update the Ln Act Coeff jacobian entries with respect to the
* mole number of species in the phase
*/
Vphase->updateLnActCoeffJac(moleSpeciesVCS);
/*
* Download the resulting calculation into the full matrix
* -> This scatter calculation is carried out in the
* volume object.
*/
Vphase->sendToVCSLnActCoeffJac(dLnActCoeffdMolNum.baseDataAddr());
}
}
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
double VCS_SOLVE::deltaG_Recalc_Rxn(int irxn, const double *const molNum,
double * const ac, double * const mu_i)
/*************************************************************************
*
* deltaG_Recalc_Rxn
//! This function recalculates the deltaG for reaction, irxn
/*!
* This function recalculates the deltaG for reaction irxn,
* given the mole numbers in molNum. It uses the temporary
* space mu_i, to hold the chemical potentials
*************************************************************************/
{
int kspec = irxn + m_numComponents;
int *pp_ptr = m_phaseParticipation[irxn];
for (int iphase = 0; iphase < m_numPhases; iphase++) {
if (pp_ptr[iphase]) {
vcs_chemPotPhase(iphase, molNum, ac, mu_i);
* space mu_i, to hold the recalculated chemical potentials.
* It only recalculates the chemical potentials for species in phases
* which participate in the irxn reaction.
*
* This function is used by the vcs_line_search algorithm() and
* should not be used widely due to the unknown state it leaves the
* system.
*
* Input
* ------------
* @param irxn Reaction number
* @param molNum Current mole numbers of species to be used as
* input to the calculation (units = kmol)
* (length = totalNuMSpecies)
*
* Output
* ------------
* @param ac output Activity coefficients (length = totalNumSpecies)
* Note this is only partially formed. Only species in
* phases that participate in the reaction will be updated
* @param mu_i diemsionless chemical potentials (length - totalNumSpecies
* Note this is only partially formed. Only species in
* phases that participate in the reaction will be updated
*
* @return Returns the dimensionless deltaG of the reaction
*
* Note, this is a dangerous routine that leaves the underlying objects in
* an unknown state.
*/
double VCS_SOLVE::deltaG_Recalc_Rxn(const int irxn, const double *const molNum,
double * const ac, double * const mu_i) {
int kspec = irxn + m_numComponents;
int *pp_ptr = m_phaseParticipation[irxn];
for (int iphase = 0; iphase < m_numPhases; iphase++) {
if (pp_ptr[iphase]) {
vcs_chemPotPhase(iphase, molNum, ac, mu_i);
}
}
double deltaG = mu_i[kspec];
double *sc_irxn = m_stoichCoeffRxnMatrix[irxn];
for (int k = 0; k < m_numComponents; k++) {
deltaG += sc_irxn[k] * mu_i[k];
}
return deltaG;
}
double deltaG = mu_i[kspec];
double *sc_irxn = m_stoichCoeffRxnMatrix[irxn];
for (int k = 0; k < m_numComponents; k++) {
deltaG += sc_irxn[k] * mu_i[k];
}
return deltaG;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
#ifdef DEBUG_MODE
double VCS_SOLVE::vcs_line_search(int irxn, double dx_orig, char *ANOTE)
// A line search algorithm is carried out on one reaction
/*
* In this routine we carry out a rough line search algorithm
* to make sure that the m_deltaGRxn_new doesn't switch signs prematurely.
*
* @param irxn Reaction number
* @param dx_orig Original step length
*
* @param ANOTE Output character string stating the conclusions of the
* line search
*
* @return Returns the optimized step length found by the search
*/
double VCS_SOLVE::vcs_line_search(const int irxn, const double dx_orig,
char * const ANOTE)
#else
double VCS_SOLVE::vcs_line_search(int irxn, double dx_orig)
double VCS_SOLVE::vcs_line_search(const int irxn, cost double dx_orig)
#endif
/*************************************************************************
*
* In this routine we carry out a rough line search algorithm
* to make sure that the m_deltaGRxn_new doesn't switch signs prematurely.
*
*
*************************************************************************/
{
int its = 0;
int k;
int kspec = ir[irxn];
const int MAXITS = 10;
double dx = dx_orig;
double *sc_irxn = m_stoichCoeffRxnMatrix[irxn];
double *molNumBase = VCS_DATA_PTR(m_molNumSpecies_old);
double *acBase = VCS_DATA_PTR(ActCoeff0);
double *ac = VCS_DATA_PTR(ActCoeff);
double molSum = 0.0;
double slope;
/*
* Calculate the deltaG value at the dx = 0.0 point
*/
double deltaGOrig = deltaG_Recalc_Rxn(irxn, molNumBase, acBase,
VCS_DATA_PTR(m_feSpecies_old));
double forig = fabs(deltaGOrig) + 1.0E-15;
if (deltaGOrig > 0.0) {
if (dx_orig > 0.0) {
dx = 0.0;
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
//plogf(" --- %s :Warning possible error dx>0 dg > 0\n", SpName[kspec]);
}
sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx>0 dg > 0");
#endif
return dx;
}
} else if (deltaGOrig < 0.0) {
if (dx_orig < 0.0) {
dx = 0.0;
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
//plogf(" --- %s :Warning possible error dx<0 dg < 0\n", SpName[kspec]);
}
sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx<0 dg < 0");
#endif
return dx;
}
} else if (deltaGOrig == 0.0) {
return 0.0;
}
if (dx_orig == 0.0) return 0.0;
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNumBase, m_numSpeciesRdc);
molSum = molNumBase[kspec];
m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx_orig;
for (k = 0; k < m_numComponents; k++) {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx_orig;
molSum += molNumBase[k];
}
double deltaG1 = deltaG_Recalc_Rxn(irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
/*
* If deltaG hasn't switched signs when going the full distance
* then we are heading in the appropriate direction, and
* we should accept the current full step size
*/
if (deltaG1 * deltaGOrig > 0.0) {
dx = dx_orig;
goto finalize;
}
/*
* If we have decreased somewhat, the deltaG return after finding
* a better estimate for the line search.
*/
if (fabs(deltaG1) < 0.8*forig) {
if (deltaG1 * deltaGOrig < 0.0) {
slope = (deltaG1 - deltaGOrig) / dx_orig;
dx = -deltaGOrig / slope;
} else {
dx = dx_orig;
}
goto finalize;
}
dx = dx_orig;
for (its = 0; its < MAXITS; its++) {
{
int its = 0;
int k;
int kspec = ir[irxn];
const int MAXITS = 10;
double dx = dx_orig;
double *sc_irxn = m_stoichCoeffRxnMatrix[irxn];
double *molNumBase = VCS_DATA_PTR(m_molNumSpecies_old);
double *acBase = VCS_DATA_PTR(ActCoeff0);
double *ac = VCS_DATA_PTR(ActCoeff);
double molSum = 0.0;
double slope;
/*
* Calculate the approximation to the total Gibbs free energy at
* the dx *= 0.5 point
* Calculate the deltaG value at the dx = 0.0 point
*/
dx *= 0.5;
m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx;
for (k = 0; k < m_numComponents; k++) {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx;
double deltaGOrig = deltaG_Recalc_Rxn(irxn, molNumBase, acBase,
VCS_DATA_PTR(m_feSpecies_old));
double forig = fabs(deltaGOrig) + 1.0E-15;
if (deltaGOrig > 0.0) {
if (dx_orig > 0.0) {
dx = 0.0;
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
//plogf(" --- %s :Warning possible error dx>0 dg > 0\n", SpName[kspec]);
}
sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx>0 dg > 0");
#endif
return dx;
}
} else if (deltaGOrig < 0.0) {
if (dx_orig < 0.0) {
dx = 0.0;
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
//plogf(" --- %s :Warning possible error dx<0 dg < 0\n", SpName[kspec]);
}
sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx<0 dg < 0");
#endif
return dx;
}
} else if (deltaGOrig == 0.0) {
return 0.0;
}
double deltaG = deltaG_Recalc_Rxn(irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
if (dx_orig == 0.0) return 0.0;
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNumBase, m_numSpeciesRdc);
molSum = molNumBase[kspec];
m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx_orig;
for (k = 0; k < m_numComponents; k++) {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx_orig;
molSum += molNumBase[k];
}
double deltaG1 = deltaG_Recalc_Rxn(irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
/*
* If deltaG hasn't switched signs when going the full distance
* then we are heading in the appropriate direction, and
* we should accept the current step
* we should accept the current full step size
*/
if (deltaG * deltaGOrig > 0.0) {
if (deltaG1 * deltaGOrig > 0.0) {
dx = dx_orig;
goto finalize;
}
/*
* If we have decreased somewhat, the deltaG return after finding
* a better estimate for the line search.
*/
if (fabs(deltaG) / forig < (1.0 - 0.1 * dx / dx_orig)) {
if (deltaG * deltaGOrig < 0.0) {
slope = (deltaG - deltaGOrig) / dx;
if (fabs(deltaG1) < 0.8*forig) {
if (deltaG1 * deltaGOrig < 0.0) {
slope = (deltaG1 - deltaGOrig) / dx_orig;
dx = -deltaGOrig / slope;
} else {
dx = dx_orig;
}
goto finalize;
}
}
dx = dx_orig;
for (its = 0; its < MAXITS; its++) {
/*
* Calculate the approximation to the total Gibbs free energy at
* the dx *= 0.5 point
*/
dx *= 0.5;
m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx;
for (k = 0; k < m_numComponents; k++) {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx;
}
double deltaG = deltaG_Recalc_Rxn(irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
/*
* If deltaG hasn't switched signs when going the full distance
* then we are heading in the appropriate direction, and
* we should accept the current step
*/
if (deltaG * deltaGOrig > 0.0) {
goto finalize;
}
/*
* If we have decreased somewhat, the deltaG return after finding
* a better estimate for the line search.
*/
if (fabs(deltaG) / forig < (1.0 - 0.1 * dx / dx_orig)) {
if (deltaG * deltaGOrig < 0.0) {
slope = (deltaG - deltaGOrig) / dx;
dx = -deltaGOrig / slope;
}
goto finalize;
}
}
finalize:
if (its >= MAXITS) {
finalize:
if (its >= MAXITS) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Rxn reduced to zero step size from %g to %g (MAXITS)",
dx_orig, dx);
sprintf(ANOTE,"Rxn reduced to zero step size from %g to %g (MAXITS)",
dx_orig, dx);
return dx;
#endif
}
#ifdef DEBUG_MODE
if (dx != dx_orig) {
sprintf(ANOTE,"Line Search reduced step size from %g to %g",
dx_orig, dx);
}
#endif
return dx;
#endif
}
#ifdef DEBUG_MODE
if (dx != dx_orig) {
sprintf(ANOTE,"Line Search reduced step size from %g to %g",
dx_orig, dx);
}
#endif
return dx;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
}

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

@ -184,8 +184,8 @@ namespace VCSnonideal {
ActCoeff.resize(nspecies0, 1.0);
ActCoeff0.resize(nspecies0, 1.0);
CurrPhAC.resize(nphase0, 0);
WtSpecies.resize(nspecies0, 0.0);
Charge.resize(nspecies0, 0.0);
m_wtSpecies.resize(nspecies0, 0.0);
m_chargeSpecies.resize(nspecies0, 0.0);
SpeciesThermo.resize(nspecies0, (VCS_SPECIES_THERMO *)0);
/*
@ -504,12 +504,12 @@ namespace VCSnonideal {
/*
* Copy over the species molecular weights
*/
vcs_vdcopy(WtSpecies, pub->WtSpecies, nspecies);
vcs_vdcopy(m_wtSpecies, pub->WtSpecies, nspecies);
/*
* Copy over the charges
*/
vcs_vdcopy(Charge, pub->Charge, nspecies);
vcs_vdcopy(m_chargeSpecies, pub->Charge, nspecies);
/*
* Malloc and Copy the VCS_SPECIES_THERMO structures
@ -726,7 +726,7 @@ namespace VCSnonideal {
* loop below starts at 1, not 0.
*/
int iSolvent = Vphase->IndSpecies[0];
double mnaught = WtSpecies[iSolvent] / 1000.;
double mnaught = m_wtSpecies[iSolvent] / 1000.;
for (int k = 1; k < Vphase->NVolSpecies; k++) {
int kspec = Vphase->IndSpecies[k];
SpecActConvention[kspec] = Vphase->ActivityConvention;

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

@ -270,7 +270,16 @@ public:
*/
void vcs_dfe(double const * const z, int kk, int ll, int lbot, int ltop);
void vcs_updateVP(int place);
//! This routine uploads the state of the system into all of the
//! vcs_VolumePhase objects in the current problem.
/*!
* @param vcsState Determines where to get the mole numbers from.
* - VCS_STATECALC_OLD -> from m_molNumSpecies_old
* - VCS_STATECALC_NEW -> from m_molNumSpecies_new
*
*/
void vcs_updateVP(const int vcsState);
int vcs_RxnStepSizes(void);
//! Calculates the total number of moles of species in all phases.
@ -459,14 +468,74 @@ public:
*/
void vcs_switch_elem_pos(int ipos, int jpos);
int vcs_rxn_adj_cg(void);
double vcs_Hessian_diag_adj(int, double);
//! Calculates reaction adjustments using a full Hessian approximation
/*!
* Calculates reaction adjustments. This does what equation 6.4-16, p. 143
* in Smith and Missen is suppose to do. However, a full matrix is
* formed and then solved via a conjugate gradient algorithm. No
* preconditioning is done.
*
* If special branching is warranted, then the program bails out.
*
* Output
* -------
* DS(I) : reaction adjustment, where I refers to the Ith species
* Special branching occurs sometimes. This causes the component basis
* to be reevaluated
* return = 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.
*
* Special attention is taken to flag cases where the direction of the
* update is contrary to the steepest descent rule. This is an important
* attribute of the regular vcs algorithm. We don't want to violate this.
*
* NOTE: currently this routine is not used.
*/
int vcs_rxn_adj_cg(void);
//! Calculates the diagonal contribution to the Hessian due to
//! the dependence of the activity coefficients on the mole numbers.
/*!
* (See framemaker notes, Eqn. 20 - VCS Equations document)
*
* We allow the diagonal to be increased positively to any degree.
* We allow the diagonal to be decreased to 1/3 of the ideal solution
* value, but no more -> it must remain positive.
*
* NOTE: currently this routine is not used
*/
double vcs_Hessian_diag_adj(int irxn, double hessianDiag_Ideal);
//! Calculates the diagonal contribution to the Hessian due to
//! the dependence of the activity coefficients on the mole numbers.
/*!
* (See framemaker notes, Eqn. 20 - VCS Equations document)
*
* NOTE: currently this routine is not used
*/
double vcs_Hessian_actCoeff_diag(int irxn);
void vcs_CalcLnActCoeffJac(const double * const moleSpeciesVCS);
#ifdef DEBUG_MODE
double vcs_line_search(int irxn, double dx_orig, char *ANOTE);
//! A line search algorithm is carried out on one reaction
/*!
* In this routine we carry out a rough line search algorithm
* to make sure that the m_deltaGRxn_new doesn't switch signs prematurely.
*
* @param irxn Reaction number
* @param dx_orig Original step length
*
* @param ANOTE Output character string stating the conclusions of the
* line search
*
*/
double vcs_line_search(const int irxn, const double dx_orig,
char * const ANOTE);
#else
double vcs_line_search(int irxn, double dx_orig);
double vcs_line_search(const int irxn, const double dx_orig);
#endif
@ -698,7 +767,34 @@ private:
void vcs_SSPhase(void);
double deltaG_Recalc_Rxn(int irxn, const double *const molNum,
//! This function recalculates the deltaG for reaction, irxn
/*!
* This function recalculates the deltaG for reaction irxn,
* given the mole numbers in molNum. It uses the temporary
* space mu_i, to hold the recalculated chemical potentials.
* It only recalculates the chemical potentials for species in phases
* which participate in the irxn reaction.
*
* Input
* ------------
* @param irxn Reaction number
* @param molNum Current mole numbers of species to be used as
* input to the calculation (units = kmol)
* (length = totalNuMSpecies)
*
* Output
* ------------
* @param ac output Activity coefficients (length = totalNumSpecies)
* Note this is only partially formed. Only species in
* phases that participate in the reaction will be updated
* @param mu_i diemsionless chemical potentials (length - totalNumSpecies
* Note this is only partially formed. Only species in
* phases that participate in the reaction will be updated
*
* @return Returns the dimensionless deltaG of the reaction
*/
double deltaG_Recalc_Rxn(const int irxn, const double *const molNum,
double * const ac, double * const mu_i);
void delete_memory();
@ -947,7 +1043,7 @@ public:
* -> Don't use this except for scaling
* purposes
*/
double m_totalMolNum;
double m_totalMolNum;
//! Total kmols of species in each phase
/*!
@ -1155,16 +1251,16 @@ public:
*/
std::vector<double> SpecLnMnaught;
//! Activity Coefficients for Species
//! Molar-based Activity Coefficients for Species
/*!
*
* Length = number of species
*/
std::vector<double> ActCoeff;
//! Activity Coefficients for Species
//! Molar-based Activity Coefficients for Species
/*!
*
* Molar based activity coeffients.
* Length = number of species
*/
std::vector<double> ActCoeff0;
@ -1188,14 +1284,16 @@ public:
/*!
* units = kg/kmol
* length = number of species
*
* note: this is a candidate for removal. I don't think we use it.
*/
std::vector<double> WtSpecies;
std::vector<double> m_wtSpecies;
//! Charge of each species
/*!
* Length = number of species
*/
std::vector<double> Charge;
std::vector<double> m_chargeSpecies;
//! Vector of pointers to thermostructures which identify the model
//! and parameters for evaluating the thermodynamic functions for that

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

@ -270,16 +270,17 @@ namespace VCSnonideal {
}
}
/* *********************************************** */
/* **** EVALUATE TOTAL MOLES, GAS AND LIQUID ***** */
/* *********************************************** */
/* - Evaluate the total moles of gas and liquid */
/* - These quantities are storred in the global variables */
/*
* Evaluate the total moles of species in the problem
*/
vcs_tmoles();
/* ******************************************* */
/* **** EVALUATE ALL CHEMICAL POTENTIALS ***** */
/* ******************************************* */
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
/* ***************************************************************************** */
/* **** EVALUATE ALL CHEMICAL POTENTIALS AT THE OLD (CURRENT) MOLE NUMBERS ***** */
/* ***************************************************************************** */
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
/*
* HKM -> If there was a machine estimate, we used to branch
* to the code segment which determined whether we needed a
@ -340,7 +341,7 @@ namespace VCSnonideal {
}
#endif
vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx));
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
}
#ifdef DEBUG_MODE
else {
@ -355,6 +356,7 @@ namespace VCSnonideal {
vcs_deltag(0, false);
iti = 0;
goto L_MAINLOOP_ALL_SPECIES;
/* ********************************************************* */
/* **** SET INITIAL VALUES FOR ITERATION ******************* */
/* **** EVALUATE REACTION ADJUSTMENTS ******************* */
@ -376,7 +378,7 @@ namespace VCSnonideal {
* We have already evaluated the major non-components
*/
if (uptodate_minors == FALSE) {
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 1, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
vcs_deltag(1, false);
}
uptodate_minors = TRUE;
@ -454,27 +456,39 @@ namespace VCSnonideal {
* Zero out the net change in moles of multispecies phases
*/
vcs_dzero(VCS_DATA_PTR(m_deltaPhaseMoles), m_numPhases);
/* **************************************************************** */
/* ***************** MAIN LOOP IN CALCULATION ******************** */
/* **************************************************************** */
/*
* Check on too many iterations.
* If we have too many iterations, Clean up and exit code even though we haven't
* converged. -> we have run out of iterations!
*/
if (m_VCount->Its > maxit) {
solveFail = -1;
goto L_RETURN_BLOCK;
}
/* ********************************************************************** */
/* ***************** MAIN LOOP IN CALCULATION *************************** */
/* ***************** LOOP OVER IRXN TO DETERMINE STEP SIZE ************** */
/* ********************************************************************** */
/*
* Loop through all of the reactions, irxn, pertaining to the
* formation reaction for species kspec in canonical form.
*
* At the end of this loop, we will have a new estimate for the
* mole numbers wt[kspec] for all species consistent with an extent
* of reaction, ds[kspec] for all noncomponent species formation
* mole numbers for all species consistent with an extent
* of reaction for all noncomponent species formation
* reactions. We will have also ensured that all predicted
* non-component mole numbers are greater than zero.
*/
if (m_VCount->Its > maxit) {
solveFail = -1;
/*
* Clean up and exit code even though we haven't
* converged. -> we have run out of iterations!
*/
goto L_RETURN_BLOCK;
}
*
* Old_Solution New_Solution Description
* -----------------------------------------------------------------------------
* m_molNumSpecies_old[kspec] m_molNumSpecies_new[kspec] Species Mole Numbers
* m_deltaMolNumSpecies[kspec] Delta in the Species Mole Numbers
*
*
*
*/
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf(" --- Main Loop Treatment of each non-component species ");
@ -493,10 +507,11 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
ANOTE[0] = '\0';
#endif
/********************************************************************/
/********************** VOLTAGE SPECIES **************************/
/********************************************************************/
if (spStatus[irxn] == VCS_SPECIES_INTERFACIALVOLTAGE) {
/********************************************************************/
/************************ VOLTAGE SPECIES ***************************/
/********************************************************************/
#ifdef DEBUG_MODE
dx = minor_alt_calc(kspec, irxn, &soldel, ANOTE);
#else
@ -505,7 +520,6 @@ namespace VCSnonideal {
m_deltaMolNumSpecies[kspec] = dx;
}
else if (spStatus[irxn] < VCS_SPECIES_MINOR) {
/********************************************************************/
/********************** ZEROED OUT SPECIES **************************/
/********************************************************************/
@ -726,6 +740,7 @@ namespace VCSnonideal {
* Form a tentative value of the new species moles
*/
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
/*
* Check for non-positive mole fraction of major species.
* If we find one, we branch to a section below. Then,
@ -851,7 +866,7 @@ namespace VCSnonideal {
* set of reactions being considered. The set of reactions
* is determined by the value of iti.
*/
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, iti, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, iti, 0, m_numSpeciesRdc);
vcs_deltag(iti, false);
/*
* Redefine the starting conditions for noncomponents
@ -902,6 +917,7 @@ namespace VCSnonideal {
m_deltaMolNumSpecies[kspec] = dx;
} /* End of Loop on ic[irxn] -> the type of species */
/***********************************************************************/
/****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/
/***********************************************************************/
@ -926,12 +942,12 @@ namespace VCSnonideal {
* Calculate the tentative change in the total number of
* moles in all of the phases
*/
dnPhase_irxn = m_deltaMolNumPhase[irxn];
for (iph = 0; iph < m_numPhases; iph++) {
m_deltaPhaseMoles[iph] += dx * dnPhase_irxn[iph];
}
}
#ifdef DEBUG_MODE
checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), VCS_DATA_PTR(m_deltaPhaseMoles), kspec+1);
#endif
@ -952,8 +968,9 @@ namespace VCSnonideal {
}
L_MAIN_LOOP_END_NO_PRINT: ;
#endif
/**************** END OF MAIN LOOP OVER FORMATION REACTIONS ************/
}
} /**************** END OF MAIN LOOP OVER FORMATION REACTIONS ************/
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
for (k = 0; k < m_numComponents; k++) {
@ -966,6 +983,7 @@ namespace VCSnonideal {
plogendl();
}
#endif
/*************************************************************************/
/*********** LIMIT REDUCTION OF BASIS SPECIES TO 99% *********************/
/*************************************************************************/
@ -1050,8 +1068,8 @@ namespace VCSnonideal {
* we have only updated a subset of the W().
*/
vcs_updateVP(1);
//vcs_dfe(VCS_DATA_PTR(wt), 1, iti, 0, m_numSpeciesTot);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_new), 1, 0, 0, m_numSpeciesTot);
//vcs_dfe(VCS_DATA_PTR(wt), VCS_STATECALC_NEW, iti, 0, m_numSpeciesTot);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_new), VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot);
/*
* Evaluate DeltaG for all components if ITI=0, and for
* major components only if ITI NE 0
@ -1063,7 +1081,7 @@ namespace VCSnonideal {
// Actually always need to calculate this
// or else nonprintouts get different results and sometimes
// fail in the line search algorithm -> Why is this?
//vcs_dfe(VCS_DATA_PTR(wt), 1, 1, 0, m_numSpeciesRdc);
//vcs_dfe(VCS_DATA_PTR(wt), VCS_STATECALC_NEW, 1, 0, m_numSpeciesRdc);
//if (iti != 0) {
// vcs_deltag(1, false);
//}
@ -1279,7 +1297,7 @@ namespace VCSnonideal {
VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test,
&usedZeroedSpecies);
if (retn != VCS_SUCCESS) return retn;
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, true);
uptodate_minors = TRUE;
if (conv) {
@ -1312,7 +1330,7 @@ namespace VCSnonideal {
}
#endif
vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx));
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, true);
uptodate_minors = TRUE;
}
@ -1537,7 +1555,7 @@ namespace VCSnonideal {
* For this special case, we must reevaluate thermo functions
*/
if (iti != 0) {
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, kspec, kspec+1);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, kspec, kspec+1);
vcs_deltag(0, false);
}
}
@ -1618,7 +1636,7 @@ namespace VCSnonideal {
* for minor species, if needed.
*/
if (iti != 0) {
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 1, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
vcs_deltag(1, false);
uptodate_minors = TRUE;
}
@ -1728,7 +1746,7 @@ namespace VCSnonideal {
/*
* Go back to evaluate the total moles of gas and liquid.
*/
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
/*
*
@ -1818,7 +1836,7 @@ namespace VCSnonideal {
* for minor species and go back to do a full iteration
*/
MajorSpeciesHaveConverged = true;
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 1, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
iti = 0;
goto L_MAINLOOP_ALL_SPECIES;
@ -1837,7 +1855,7 @@ namespace VCSnonideal {
* for minor species and go back to do a full iteration
*/
MajorSpeciesHaveConverged = true;
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 1, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
iti = 0;
goto L_MAINLOOP_ALL_SPECIES;
@ -2587,7 +2605,7 @@ namespace VCSnonideal {
}
}
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesTot);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot);
vcs_deltag(0, true);
}
@ -2722,8 +2740,8 @@ namespace VCSnonideal {
* only step is being carried out, then we don't need to
* update the minor noncomponents.
*/
// vcs_dfe(dptr, 1, iti, 0, m_numSpeciesRdc);
vcs_dfe(dptr, 1, 0, 0, m_numSpeciesRdc);
// vcs_dfe(dptr, VCS_STATECALC_NEW, iti, 0, m_numSpeciesRdc);
vcs_dfe(dptr, VCS_STATECALC_NEW, 0, 0, m_numSpeciesRdc);
/*
* Evaluate DeltaG for all components if ITI=0, and for
* major components only if ITI NE 0
@ -4058,26 +4076,25 @@ namespace VCSnonideal {
plogf("We have an inconsistency!\n");
exit(-1);
}
if (Charge[kspec] != -1.0) {
if (m_chargeSpecies[kspec] != -1.0) {
plogf("We have an unexpected situation!\n");
exit(-1);
}
#endif
mu_i[kspec] = m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_phi;
mu_i[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
} else {
if (SSPhase[kspec]) {
mu_i[kspec] = m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_phi;
mu_i[kspec] = m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
} else if (molNum[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
mu_i[kspec] = m_SSfeSpecies[kspec] + log(ac[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles - SpecLnMnaught[kspec] + Charge[kspec] * Faraday_phi;
- tlogMoles - SpecLnMnaught[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
} else {
mu_i[kspec] = m_SSfeSpecies[kspec] + log(ac[kspec] * molNum[kspec])
- tlogMoles - SpecLnMnaught[kspec] + Charge[kspec] * Faraday_phi;
- tlogMoles - SpecLnMnaught[kspec] + m_chargeSpecies[kspec] * Faraday_phi;
}
}
}
}
/*****************************************************************************/
// Calculalte the dimensionless chemical potentials of all species or
@ -4105,7 +4122,7 @@ namespace VCSnonideal {
* Ideal Mixtures:
*
* m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph])
* + Charge[I] * Faraday_dim * m_phasePhi[iphase];
* + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
*
* ( This is equivalent to the adding the log of the
* mole fraction onto the standard chemical
@ -4116,7 +4133,7 @@ namespace VCSnonideal {
*
* m_feSpecies(I) = m_SSfeSpecies(I)
* + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph])
* + Charge[I] * Faraday_dim * m_phasePhi[iphase];
* + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
*
* ( This is equivalent to the adding the log of the
* mole fraction multiplied by the activity coefficient
@ -4130,7 +4147,7 @@ namespace VCSnonideal {
* m_feSpecies(I) = m_SSfeSpecies(I)
* + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph])
* - ln(Mnaught * m_units)
* + Charge[I] * Faraday_dim * m_phasePhi[iphase];
* + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
*
* note: m_SSfeSpecies(I) is the molality based standard state.
* However, ActCoeff[I] is the molar based activity coefficient
@ -4142,7 +4159,7 @@ namespace VCSnonideal {
*
* m_feSpecies(I) = m_SSfeSpecies(I)
* + ln(ActCoeff_M[I] * m(I))
* + Charge[I] * Faraday_dim * m_phasePhi[iphase];
* + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
* where m[I] is the molality of the ith solute
*
* m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units)
@ -4157,7 +4174,7 @@ namespace VCSnonideal {
*
* The chemical potential is calculated as:
*
* m_feSpecies(I)(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF))
* m_feSpecies(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF))
*
* Handling of "Species" Representing Interfacial Voltages
* ---------------------------------------------------------
@ -4191,11 +4208,11 @@ namespace VCSnonideal {
* -> This can either be the current solution vector WT()
* or the actual solution vector W()
*
* @param kk Determines whether z is old or new or tentative:
* 1: Use the tentative values for the total number of
* moles in the phases, i.e., use TG1 instead of TG etc.
* 0: Use the base values of the total number of
* moles in each system.
* @param kk Determines whether z is old or new or tmp:
* VCS_STATECALC_NEW: Use the tentative values for the total number of
* moles in the phases, i.e., use TG1 instead of TG etc.
* VCS_STATECALC_OLD: Use the base values of the total number of
* moles in each system.
*
* Also needed:
* ff : standard state chemical potentials. These are the
@ -4211,6 +4228,14 @@ namespace VCSnonideal {
vcs_VolPhase *Vphase;
VCS_SPECIES_THERMO *st_ptr;
#ifdef DEBUG_MODE
if (kk != VCS_STATECALC_OLD && kk != VCS_STATECALC_NEW) {
plogf(" --- Subroutine vcs_dfe called with bad kk value: %d", kk);
plogendl();
exit(-1);
}
#endif
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
if (ll == 0) {
@ -4225,11 +4250,11 @@ namespace VCSnonideal {
} else {
plogf(" --- Subroutine vcs_dfe called for components and majors");
}
if (kk == 1) plogf(" using tentative solution\n");
else plogf("\n");
if (kk == VCS_STATECALC_NEW) plogf(" using tentative solution");
plogendl();
}
#endif
if (kk <= 0) {
if (kk <= VCS_STATECALC_OLD) {
tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
} else {
tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
@ -4311,13 +4336,13 @@ namespace VCSnonideal {
plogf("We have an inconsistency!\n");
exit(-1);
}
if (Charge[kspec] != -1.0) {
if (m_chargeSpecies[kspec] != -1.0) {
plogf("We have an unexpected situation!\n");
exit(-1);
}
#endif
m_feSpecies_curr[kspec] =
m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * m_phasePhi[iphase];
m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase];
} else {
if (SSPhase[kspec]) {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
@ -4328,14 +4353,14 @@ namespace VCSnonideal {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]
+ log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ Charge[kspec] * Faraday_dim * m_phasePhi[iphase];
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase];
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ Charge[kspec] * Faraday_dim * m_phasePhi[iphase];
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase];
}
}
}
@ -4354,13 +4379,13 @@ namespace VCSnonideal {
plogf("We have an inconsistency!\n");
exit(-1);
}
if (Charge[kspec] != -1.0) {
if (m_chargeSpecies[kspec] != -1.0) {
plogf("We have an unexpected situation!\n");
exit(-1);
}
#endif
m_feSpecies_curr[kspec] =
m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * m_phasePhi[iphase];
m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase];
} else {
if (SSPhase[kspec]) {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
@ -4371,14 +4396,14 @@ namespace VCSnonideal {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]
+ log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ Charge[kspec] * Faraday_dim * m_phasePhi[iphase]; ;
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase]; ;
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ Charge[kspec] * Faraday_dim * m_phasePhi[iphase];
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase];
}
}
}
@ -4398,13 +4423,13 @@ namespace VCSnonideal {
plogf("We have an inconsistency!\n");
exit(-1);
}
if (Charge[kspec] != -1.0) {
if (m_chargeSpecies[kspec] != -1.0) {
plogf("We have an unexpected situation!\n");
exit(-1);
}
#endif
m_feSpecies_curr[kspec] =
m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * m_phasePhi[iphase]; ;
m_SSfeSpecies[kspec] + m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase]; ;
} else {
if (SSPhase[kspec]) {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
@ -4434,7 +4459,6 @@ namespace VCSnonideal {
}
#endif
}
/*****************************************************************************/
#ifdef DEBUG_MODE
@ -4472,7 +4496,6 @@ namespace VCSnonideal {
plogendl();
}
#endif
/*****************************************************************************/
// Calculate the norm of a deltaGibbs free energy vector
@ -4495,7 +4518,6 @@ namespace VCSnonideal {
}
return (std::sqrt(tmp / m_numRxnRdc));
}
/*****************************************************************************/
// Calculates the total number of moles of species in all phases.
@ -4532,35 +4554,35 @@ namespace VCSnonideal {
}
}
m_totalMolNum = sum;
}
}
/*****************************************************************************/
void VCS_SOLVE::vcs_updateVP (int place)
/*************************************************************************
* vcs_updateVP()
*
* This routine uploads the state of the system into all of the
* VolumePhase objects in the current problem.
* place
* 0 -> from m_molNumSpecies_old
* 1 -> from wt
*************************************************************************/
{
// This routine uploads the state of the system into all of the
// vcs_VolPhase objects in the current problem.
/*
* @param vcsState Determines where to get the mole numbers from.
* - VCS_STATECALC_OLD -> from m_molNumSpecies_old
* - VCS_STATECALC_NEW -> from m_molNumSpecies_new
*
*/
void VCS_SOLVE::vcs_updateVP(const int vcsState) {
vcs_VolPhase *Vphase;
for (int i = 0; i < m_numPhases; i++) {
Vphase = VPhaseList[i];
if (place == 0) {
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old),
if (vcsState == VCS_STATECALC_OLD) {
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(m_tPhaseMoles_old), i);
} else if (place == 1) {
} else if (vcsState == VCS_STATECALC_NEW) {
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_new),
VCS_DATA_PTR(m_tPhaseMoles_new), i);
} else {
plogf("we shouldn't be here\n");
}
#ifdef DEBUG_MODE
else {
plogf("we shouldn't be here");
plogendl();
exit(-1);
}
#endif
}
}
@ -4659,8 +4681,8 @@ namespace VCSnonideal {
SWAP(SpecLnMnaught[k1], SpecLnMnaught[k2], t1);
SWAP(ActCoeff[k1], ActCoeff[k2], t1);
SWAP(ActCoeff0[k1], ActCoeff0[k2], t1);
SWAP(WtSpecies[k1], WtSpecies[k2], t1);
SWAP(Charge[k1], Charge[k2], t1);
SWAP(m_wtSpecies[k1], m_wtSpecies[k2], t1);
SWAP(m_chargeSpecies[k1], m_chargeSpecies[k2], t1);
SWAP(SpeciesThermo[k1], SpeciesThermo[k2], st_tmp);
SWAP(VolPM[k1], VolPM[k2], t1);