From f741f96b3c7eefe4b4ad7daadbaacf2e672ff238 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 14 May 2008 16:50:34 +0000 Subject: [PATCH] More variable and documentation checks. Added the VCS_STATECALC_ defines. Will try to clean up the evaluation of phase properties using this concept. --- Cantera/src/equil/vcs_defs.h | 16 + Cantera/src/equil/vcs_inest.cpp | 2 +- Cantera/src/equil/vcs_report.cpp | 2 +- Cantera/src/equil/vcs_rxnadj.cpp | 902 +++++++++++++++-------------- Cantera/src/equil/vcs_solve.cpp | 10 +- Cantera/src/equil/vcs_solve.h | 122 +++- Cantera/src/equil/vcs_solve_TP.cpp | 220 +++---- 7 files changed, 706 insertions(+), 568 deletions(-) diff --git a/Cantera/src/equil/vcs_defs.h b/Cantera/src/equil/vcs_defs.h index 42046a326..362ed2588 100644 --- a/Cantera/src/equil/vcs_defs.h +++ b/Cantera/src/equil/vcs_defs.h @@ -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 //@} diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index 779cf2c68..9d97d6b20 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -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]; } diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 5c2d7a618..10b51b90c 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -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; diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index 82dacc10e..ca87712f6 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -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 -#include -#include #include "vcs_solve.h" #include "vcs_internal.h" #include "vcs_VolPhase.h" +#include +#include +#include + 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; -} -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ + /*****************************************************************************/ } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 44d288f14..c34712c46 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -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; diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 8d2857b51..4ee4d4ba0 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -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 SpecLnMnaught; - //! Activity Coefficients for Species + //! Molar-based Activity Coefficients for Species /*! * * Length = number of species */ std::vector ActCoeff; - //! Activity Coefficients for Species + //! Molar-based Activity Coefficients for Species /*! - * + * Molar based activity coeffients. * Length = number of species */ std::vector 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 WtSpecies; + std::vector m_wtSpecies; //! Charge of each species /*! * Length = number of species */ - std::vector Charge; + std::vector m_chargeSpecies; //! Vector of pointers to thermostructures which identify the model //! and parameters for evaluating the thermodynamic functions for that diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 0aad9559d..dba6eb3b6 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -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);