diff --git a/Cantera/src/equil/vcs_Gibbs.cpp b/Cantera/src/equil/vcs_Gibbs.cpp index 18c8dd29b..521f15e41 100644 --- a/Cantera/src/equil/vcs_Gibbs.cpp +++ b/Cantera/src/equil/vcs_Gibbs.cpp @@ -19,12 +19,12 @@ namespace VCSnonideal { -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ -double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot, - double *tPhMoles) + double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot, + double *tPhMoles) /************************************************************************* * @@ -36,14 +36,14 @@ double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot, * Note, for this algorithm this function should be MONOTONICALLY * DECREASING. *************************************************************************/ -{ + { double g = 0.0; for (int iph = 0; iph < NPhase; iph++) { vcs_VolPhase *Vphase = VPhaseList[iph]; if ((TPhInertMoles[iph] > 0.0) && (tPhMoles[iph] > 0.0)) { g += TPhInertMoles[iph] * - log(TPhInertMoles[iph] / tPhMoles[iph]); + log(TPhInertMoles[iph] / tPhMoles[iph]); if (Vphase->GasPhase) { g += TPhInertMoles[iph] * log(Pres); } @@ -55,48 +55,40 @@ double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot, } return g; -} + } -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -double VCS_SOLVE::vcs_GibbsPhase(int iphase, double *w, double *fe) - - /************************************************************************* - * - * vcs_Total_Gibbs: - * - * Calculate the total dimensionless Gibbs free energy - * -> Inert species are handled as if they had a standard free - * energy of zero. - * Note, for this algorithm this function should be MONOTONICALLY - * DECREASING. - *************************************************************************/ -{ + // Calculate the total dimensionless Gibbs free energy of a single phase + /* + * -> Inert species are handled as if they had a standard free + * energy of zero and if they obeyed ideal solution/gas theory + * + * @param iphase ID of the phase + * @param w Species mole number vector + * @param fe vector of partial molar free energies of the species. + */ + double VCS_SOLVE::vcs_GibbsPhase(int iphase, const double * const w, + const double * const fe) { double g = 0.0; - - vcs_VolPhase *Vphase = VPhaseList[iphase]; - if ((TPhInertMoles[iphase] > 0.0) && (TPhMoles[iphase] > 0.0)) { - g += TPhInertMoles[iphase] * - log(TPhInertMoles[iphase] / TPhMoles[iphase]); + double phaseMols = 0.0; + for (int kspec = 0; kspec < m_numSpeciesRdc; ++kspec) { + if (PhaseID[kspec] == iphase) { + g += w[kspec] * fe[kspec]; + phaseMols += w[kspec]; + } + } + + if (TPhInertMoles[iphase] > 0.0) { + phaseMols += TPhInertMoles[iphase]; + g += TPhInertMoles[iphase] * log(TPhInertMoles[iphase] / phaseMols); + vcs_VolPhase *Vphase = VPhaseList[iphase]; if (Vphase->GasPhase == iphase) { g += TPhInertMoles[iphase] * log(Pres); } } - - for (int kspec = 0; kspec < m_numSpeciesRdc; ++kspec) { - if (PhaseID[kspec] == iphase) { - g += w[kspec] * fe[kspec]; - } - } - + return g; -} + } -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ } diff --git a/Cantera/src/equil/vcs_elem.cpp b/Cantera/src/equil/vcs_elem.cpp index 275dbbfbe..bba1ec12e 100644 --- a/Cantera/src/equil/vcs_elem.cpp +++ b/Cantera/src/equil/vcs_elem.cpp @@ -1,5 +1,8 @@ /** * @file vcs_elem.cpp + * This file contains the algorithm for checking the satisfaction of the + * element abundances constraints and the algorithm for fixing violations + * of the element abundances constraints. */ /* * $Id$ @@ -11,592 +14,584 @@ namespace VCSnonideal { -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ -void VCS_SOLVE::vcs_elab(void) - - /************************************************************************* - * - * vcs_elab: - * - * Computes the elemental abundances vector, m_elemAbundances[], and stores it - * back into the global structure - *************************************************************************/ -{ - for (int j = 0; j < m_numElemConstraints; ++j) { - m_elemAbundances[j] = 0.0; - for (int i = 0; i < m_numSpeciesTot; ++i) { - if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - m_elemAbundances[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i]; + //! Computes the current elemental abundances vector + /*! + * Computes the elemental abundances vector, m_elemAbundances[], and stores it + * back into the global structure + */ + void VCS_SOLVE::vcs_elab() { + for (int j = 0; j < m_numElemConstraints; ++j) { + m_elemAbundances[j] = 0.0; + for (int i = 0; i < m_numSpeciesTot; ++i) { + if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + m_elemAbundances[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i]; + } } } } -} -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ -/* - * - * vcs_elabcheck: - * - * This function checks to see if the element abundances are in - * compliance. If they are, then TRUE is returned. If not, - * FALSE is returned. Note the number of constraints checked is - * usually equal to the number of components in the problem. This - * routine can check satisfaction of all of the constraints in the - * problem, which is equal to ne. However, the solver can't fix - * breakage of constraints above nc, because that nc is the - * range space by definition. Satisfaction of extra constraints would - * have had to occur in the problem specification. - * - * The constraints should be broken up into 2 sections. If - * a constraint involves a formula matrix with positive and - * negative signs, and eaSet = 0.0, then you can't expect that the - * sum will be zero. There may be roundoff that inhibits this. - * However, if the formula matrix is all of one sign, then - * this requires that all species with nonzero entries in the - * formula matrix be identically zero. We put this into - * the logic below. - * - * Input - * ------- - * ibound = 1 : Checks constraints up to the number of elements - * 0 : Checks constraints up to the number of components. - * - */ -int VCS_SOLVE::vcs_elabcheck(int ibound) { - int i; - int top = m_numComponents; - double eval, scale; - int numNonZero; - bool multisign = false; - if (ibound) { - top = m_numElemConstraints; - } + + /* - * Require 12 digits of accuracy on non-zero constraints. + * + * vcs_elabcheck: + * + * This function checks to see if the element abundances are in + * compliance. If they are, then TRUE is returned. If not, + * FALSE is returned. Note the number of constraints checked is + * usually equal to the number of components in the problem. This + * routine can check satisfaction of all of the constraints in the + * problem, which is equal to ne. However, the solver can't fix + * breakage of constraints above nc, because that nc is the + * range space by definition. Satisfaction of extra constraints would + * have had to occur in the problem specification. + * + * The constraints should be broken up into 2 sections. If + * a constraint involves a formula matrix with positive and + * negative signs, and eaSet = 0.0, then you can't expect that the + * sum will be zero. There may be roundoff that inhibits this. + * However, if the formula matrix is all of one sign, then + * this requires that all species with nonzero entries in the + * formula matrix be identically zero. We put this into + * the logic below. + * + * Input + * ------- + * ibound = 1 : Checks constraints up to the number of elements + * 0 : Checks constraints up to the number of components. + * */ - for (i = 0; i < top; ++i) { - if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > (fabs(m_elemAbundancesGoal[i]) * 1.0e-12)) { - /* - * This logic is for charge neutrality condition - */ - if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) { - AssertThrowVCS(m_elemAbundancesGoal[i] == 0.0, "vcs_elabcheck"); - } - if (m_elemAbundancesGoal[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) { - scale = VCS_DELETE_MINORSPECIES_CUTOFF; + int VCS_SOLVE::vcs_elabcheck(int ibound) { + int i; + int top = m_numComponents; + double eval, scale; + int numNonZero; + bool multisign = false; + if (ibound) { + top = m_numElemConstraints; + } + /* + * Require 12 digits of accuracy on non-zero constraints. + */ + for (i = 0; i < top; ++i) { + if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > (fabs(m_elemAbundancesGoal[i]) * 1.0e-12)) { /* - * Find out if the constraint is a multisign constraint. - * If it is, then we have to worry about roundoff error - * in the addition of terms. We are limited to 13 - * digits of finite arithmetic accuracy. + * This logic is for charge neutrality condition */ - numNonZero = 0; - multisign = false; - for (int kspec = 0; kspec < m_numSpeciesTot; kspec++) { - eval = FormulaMatrix[i][kspec]; + if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) { + AssertThrowVCS(m_elemAbundancesGoal[i] == 0.0, "vcs_elabcheck"); + } + if (m_elemAbundancesGoal[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) { + scale = VCS_DELETE_MINORSPECIES_CUTOFF; + /* + * Find out if the constraint is a multisign constraint. + * If it is, then we have to worry about roundoff error + * in the addition of terms. We are limited to 13 + * digits of finite arithmetic accuracy. + */ + numNonZero = 0; + multisign = false; + for (int kspec = 0; kspec < m_numSpeciesTot; kspec++) { + eval = FormulaMatrix[i][kspec]; + if (eval < 0.0) { + multisign = true; + } + if (eval != 0.0) { + scale = MAX(scale, fabs(eval * m_molNumSpecies_old[kspec])); + numNonZero++; + } + } + if (multisign) { + if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) { + return FALSE; + } + } else { + if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) { + return FALSE; + } + } + } else { + /* + * For normal element balances, we require absolute compliance + * even for rediculously small numbers. + */ + if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) { + return FALSE; + } else { + return FALSE; + } + } + } + } + return TRUE; + } /* vcs_elabcheck() *********************************************************/ + + /*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ + + void VCS_SOLVE::vcs_elabPhase(int iphase, double * const elemAbundPhase) + + /************************************************************************* + * + * vcs_elabPhase: + * + * Computes the elemental abundances vector for a single phase, + * elemAbundPhase[], and returns it through the argument list. + * The mole numbers of species are taken from the current value + * in m_molNumSpecies_old[]. + *************************************************************************/ + { + int i, j; + for (j = 0; j < m_numElemConstraints; ++j) { + elemAbundPhase[j] = 0.0; + for (i = 0; i < m_numSpeciesTot; ++i) { + if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + if (PhaseID[i] == iphase) { + elemAbundPhase[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i]; + } + } + } + } + } + + /*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ + + int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) + + /************************************************************************** + * + * vcs_elcorr: + * + * This subroutine corrects for element abundances. At the end of the + * surbroutine, the total moles in all phases are recalculated again, + * because we have changed the number of moles in this routine. + * + * Input + * -> temporary work vectors: + * aa[ne*ne] + * x[ne] + * + * Return Values: + * 0 = Nothing of significance happened, + * Element abundances were and still are good. + * 1 = The solution changed significantly; + * The element abundances are now good. + * 2 = The solution changed significantly, + * The element abundances are still bad. + * 3 = The solution changed significantly, + * The element abundances are still bad and a component + * species got zeroed out. + * + * Internal data to be worked on:: + * + * ga Current element abundances + * m_elemAbundancesGoal Required elemental abundances + * m_molNumSpecies_old Current mole number of species. + * FormulaMatrix[][] Formular matrix of the species + * ne Number of elements + * nc Number of components. + * + * NOTES: + * This routine is turning out to be very problematic. There are + * lots of special cases and problems with zeroing out species. + * + * Still need to check out when we do loops over nc vs. ne. + * + *************************************************************************/ + { + int i, j, retn = 0, kspec, goodSpec, its; + double xx, par, saveDir, dir; + +#ifdef DEBUG_MODE + double l2before = 0.0, l2after = 0.0; + std::vector ga_save(m_numElemConstraints, 0.0); + vcs_dcopy(VCS_DATA_PTR(ga_save), VCS_DATA_PTR(m_elemAbundances), m_numElemConstraints); + if (vcs_debug_print_lvl >= 2) { + plogf(" --- vcsc_elcorr: Element abundances correction routine"); + if (m_numElemConstraints != m_numComponents) { + plogf(" (m_numComponents != m_numElemConstraints)"); + } + plogf("\n"); + } + + for (i = 0; i < m_numElemConstraints; ++i) { + x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; + } + l2before = 0.0; + for (i = 0; i < m_numElemConstraints; ++i) { + l2before += x[i] * x[i]; + } + l2before = sqrt(l2before/m_numElemConstraints); +#endif + + /* + * Special section to take out single species, single component, + * moles. These are species which have non-zero entries in the + * formula matrix, and no other species have zero values either. + * + */ + int numNonZero = 0; + bool changed = false; + bool multisign = false; + for (i = 0; i < m_numElemConstraints; ++i) { + numNonZero = 0; + multisign = false; + for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { + if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + double eval = FormulaMatrix[i][kspec]; if (eval < 0.0) { multisign = true; } if (eval != 0.0) { - scale = MAX(scale, fabs(eval * m_molNumSpecies_old[kspec])); numNonZero++; } } - if (multisign) { - if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) { - return FALSE; - } - } else { - if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) { - return FALSE; - } - } - } else { - /* - * For normal element balances, we require absolute compliance - * even for rediculously small numbers. - */ - if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) { - return FALSE; - } else { - return FALSE; - } } - } - } - return TRUE; -} /* vcs_elabcheck() *********************************************************/ - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -void VCS_SOLVE::vcs_elabPhase(int iphase, double * const elemAbundPhase) - - /************************************************************************* - * - * vcs_elabPhase: - * - * Computes the elemental abundances vector for a single phase, - * elemAbundPhase[], and returns it through the argument list. - * The mole numbers of species are taken from the current value - * in m_molNumSpecies_old[]. - *************************************************************************/ -{ - int i, j; - for (j = 0; j < m_numElemConstraints; ++j) { - elemAbundPhase[j] = 0.0; - for (i = 0; i < m_numSpeciesTot; ++i) { - if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - if (PhaseID[i] == iphase) { - elemAbundPhase[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i]; - } - } - } - } -} - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) - - /************************************************************************** - * - * vcs_elcorr: - * - * This subroutine corrects for element abundances. At the end of the - * surbroutine, the total moles in all phases are recalculated again, - * because we have changed the number of moles in this routine. - * - * Input - * -> temporary work vectors: - * aa[ne*ne] - * x[ne] - * - * Return Values: - * 0 = Nothing of significance happened, - * Element abundances were and still are good. - * 1 = The solution changed significantly; - * The element abundances are now good. - * 2 = The solution changed significantly, - * The element abundances are still bad. - * 3 = The solution changed significantly, - * The element abundances are still bad and a component - * species got zeroed out. - * - * Internal data to be worked on:: - * - * ga Current element abundances - * m_elemAbundancesGoal Required elemental abundances - * m_molNumSpecies_old Current mole number of species. - * FormulaMatrix[][] Formular matrix of the species - * ne Number of elements - * nc Number of components. - * - * NOTES: - * This routine is turning out to be very problematic. There are - * lots of special cases and problems with zeroing out species. - * - * Still need to check out when we do loops over nc vs. ne. - * - *************************************************************************/ -{ - int i, j, retn = 0, kspec, goodSpec, its; - double xx, par, saveDir, dir; - -#ifdef DEBUG_MODE - double l2before = 0.0, l2after = 0.0; - std::vector ga_save(m_numElemConstraints, 0.0); - vcs_dcopy(VCS_DATA_PTR(ga_save), VCS_DATA_PTR(m_elemAbundances), m_numElemConstraints); - if (vcs_debug_print_lvl >= 2) { - plogf(" --- vcsc_elcorr: Element abundances correction routine"); - if (m_numElemConstraints != m_numComponents) { - plogf(" (m_numComponents != m_numElemConstraints)"); - } - plogf("\n"); - } - - for (i = 0; i < m_numElemConstraints; ++i) { - x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; - } - l2before = 0.0; - for (i = 0; i < m_numElemConstraints; ++i) { - l2before += x[i] * x[i]; - } - l2before = sqrt(l2before/m_numElemConstraints); -#endif - - /* - * Special section to take out single species, single component, - * moles. These are species which have non-zero entries in the - * formula matrix, and no other species have zero values either. - * - */ - int numNonZero = 0; - bool changed = false; - bool multisign = false; - for (i = 0; i < m_numElemConstraints; ++i) { - numNonZero = 0; - multisign = false; - for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { - if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double eval = FormulaMatrix[i][kspec]; - if (eval < 0.0) { - multisign = true; - } - if (eval != 0.0) { - numNonZero++; - } - } - } - if (!multisign) { - if (numNonZero < 2) { - for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { - if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double eval = FormulaMatrix[i][kspec]; - if (eval > 0.0) { - m_molNumSpecies_old[kspec] = m_elemAbundancesGoal[i] / eval; - changed = true; - } - } - } - } else { - int numCompNonZero = 0; - int compID = -1; - for (kspec = 0; kspec < m_numComponents; kspec++) { - if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double eval = FormulaMatrix[i][kspec]; - if (eval > 0.0) { - compID = kspec; - numCompNonZero++; - } - } - } - if (numCompNonZero == 1) { - double diff = m_elemAbundancesGoal[i]; - for (kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) { + if (!multisign) { + if (numNonZero < 2) { + for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { double eval = FormulaMatrix[i][kspec]; - diff -= eval * m_molNumSpecies_old[kspec]; + if (eval > 0.0) { + m_molNumSpecies_old[kspec] = m_elemAbundancesGoal[i] / eval; + changed = true; + } + } + } + } else { + int numCompNonZero = 0; + int compID = -1; + for (kspec = 0; kspec < m_numComponents; kspec++) { + if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + double eval = FormulaMatrix[i][kspec]; + if (eval > 0.0) { + compID = kspec; + numCompNonZero++; + } + } + } + if (numCompNonZero == 1) { + double diff = m_elemAbundancesGoal[i]; + for (kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) { + if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + double eval = FormulaMatrix[i][kspec]; + diff -= eval * m_molNumSpecies_old[kspec]; + } + m_molNumSpecies_old[compID] = MAX(0.0,diff/FormulaMatrix[i][compID]); + changed = true; } - m_molNumSpecies_old[compID] = MAX(0.0,diff/FormulaMatrix[i][compID]); - changed = true; } } } } - } - if (changed) { - vcs_elab(); - } + if (changed) { + vcs_elab(); + } - /* - * Section to check for maximum bounds errors on all species - * due to elements. - * This may only be tried on element types which are VCS_ELEM_TYPE_ABSPOS. - * This is because no other species may have a negative number of these. - * - * Note, also we can do this over ne, the number of elements, not just - * the number of components. - */ - changed = false; - for (i = 0; i < m_numElemConstraints; ++i) { - int elType = m_elType[i]; - if (elType == VCS_ELEM_TYPE_ABSPOS) { - for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { - if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - double atomComp = FormulaMatrix[i][kspec]; - if (atomComp > 0.0) { - double maxPermissible = m_elemAbundancesGoal[i] / atomComp; - if (m_molNumSpecies_old[kspec] > maxPermissible) { + /* + * Section to check for maximum bounds errors on all species + * due to elements. + * This may only be tried on element types which are VCS_ELEM_TYPE_ABSPOS. + * This is because no other species may have a negative number of these. + * + * Note, also we can do this over ne, the number of elements, not just + * the number of components. + */ + changed = false; + for (i = 0; i < m_numElemConstraints; ++i) { + int elType = m_elType[i]; + if (elType == VCS_ELEM_TYPE_ABSPOS) { + for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { + if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + double atomComp = FormulaMatrix[i][kspec]; + if (atomComp > 0.0) { + double maxPermissible = m_elemAbundancesGoal[i] / atomComp; + if (m_molNumSpecies_old[kspec] > maxPermissible) { #ifdef DEBUG_MODE - if (vcs_debug_print_lvl >= 3) { - plogf(" --- vcs_elcorr: Reduced species %s from %g to %g due to %s max bounds constraint\n", - SpName[kspec].c_str(), m_molNumSpecies_old[kspec], maxPermissible, ElName[i].c_str()); - } -#endif - m_molNumSpecies_old[kspec] = maxPermissible; - changed = true; - if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF) { - m_molNumSpecies_old[kspec] = 0.0; - if (SSPhase[kspec]) { - spStatus[kspec] = VCS_SPECIES_ZEROEDSS; - } else { - spStatus[kspec] = VCS_SPECIES_ZEROEDMS; - } -#ifdef DEBUG_MODE - if (vcs_debug_print_lvl >= 2) { - plogf(" --- vcs_elcorr: Zeroed species %s and changed status to %d due to max bounds constraint\n", - SpName[kspec].c_str(), spStatus[kspec]); + if (vcs_debug_print_lvl >= 3) { + plogf(" --- vcs_elcorr: Reduced species %s from %g to %g due to %s max bounds constraint\n", + SpName[kspec].c_str(), m_molNumSpecies_old[kspec], maxPermissible, ElName[i].c_str()); } #endif + m_molNumSpecies_old[kspec] = maxPermissible; + changed = true; + if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF) { + m_molNumSpecies_old[kspec] = 0.0; + if (SSPhase[kspec]) { + spStatus[kspec] = VCS_SPECIES_ZEROEDSS; + } else { + spStatus[kspec] = VCS_SPECIES_ZEROEDMS; + } +#ifdef DEBUG_MODE + if (vcs_debug_print_lvl >= 2) { + plogf(" --- vcs_elcorr: Zeroed species %s and changed status to %d due to max bounds constraint\n", + SpName[kspec].c_str(), spStatus[kspec]); + } +#endif + } } } } } } } - } - // Recalculate the element abundances if something has changed. - if (changed) { - vcs_elab(); - } + // Recalculate the element abundances if something has changed. + if (changed) { + vcs_elab(); + } - /* - * Ok, do the general case. Linear algebra problem is - * of length nc, not ne, as there may be degenerate rows when - * nc .ne. ne. - */ - for (i = 0; i < m_numComponents; ++i) { - x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; - if (fabs(x[i]) > 1.0E-13) retn = 1; - for (j = 0; j < m_numComponents; ++j) { - aa[j + i*m_numElemConstraints] = FormulaMatrix[j][i]; - } - } - i = vcsUtil_mlequ(aa, m_numElemConstraints, m_numComponents, x, 1); - if (i == 1) { - plogf("vcs_elcorr ERROR: mlequ returned error condition\n"); - return VCS_FAILED_CONVERGENCE; - } - /* - * Now apply the new direction without creating negative species. - */ - par = 0.5; - for (i = 0; i < m_numComponents; ++i) { - if (m_molNumSpecies_old[i] > 0.0) { - xx = -x[i] / m_molNumSpecies_old[i]; - if (par < xx) par = xx; - } - } - if (par > 100.0) { - par = 100.0; - } - par = 1.0 / par; - if (par < 1.0 && par > 0.0) { - retn = 2; - par *= 0.9999; - for (i = 0; i < m_numComponents; ++i) { - double tmp = m_molNumSpecies_old[i] + par * x[i]; - if (tmp > 0.0) { - m_molNumSpecies_old[i] = tmp; - } else { - if (SSPhase[i]) { - m_molNumSpecies_old[i] = 0.0; - } else { - m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001; - } - } - } - } else { - for (i = 0; i < m_numComponents; ++i) { - double tmp = m_molNumSpecies_old[i] + x[i]; - if (tmp > 0.0) { - m_molNumSpecies_old[i] = tmp; - } else { - if (SSPhase[i]) { - m_molNumSpecies_old[i] = 0.0; - } else { - m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001; - } - } - } - } - - /* - * We have changed the element abundances. Calculate them again - */ - vcs_elab(); - /* - * We have changed the total moles in each phase. Calculate them again - */ - vcs_tmoles(); - - /* - * Try some ad hoc procedures for fixing the problem - */ - if (retn >= 2) { /* - * First find a species whose adjustment is a win-win - * situation. + * Ok, do the general case. Linear algebra problem is + * of length nc, not ne, as there may be degenerate rows when + * nc .ne. ne. */ - for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { - if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - continue; + for (i = 0; i < m_numComponents; ++i) { + x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i]; + if (fabs(x[i]) > 1.0E-13) retn = 1; + for (j = 0; j < m_numComponents; ++j) { + aa[j + i*m_numElemConstraints] = FormulaMatrix[j][i]; } - saveDir = 0.0; - goodSpec = TRUE; + } + i = vcsUtil_mlequ(aa, m_numElemConstraints, m_numComponents, x, 1); + if (i == 1) { + plogf("vcs_elcorr ERROR: mlequ returned error condition\n"); + return VCS_FAILED_CONVERGENCE; + } + /* + * Now apply the new direction without creating negative species. + */ + par = 0.5; + for (i = 0; i < m_numComponents; ++i) { + if (m_molNumSpecies_old[i] > 0.0) { + xx = -x[i] / m_molNumSpecies_old[i]; + if (par < xx) par = xx; + } + } + if (par > 100.0) { + par = 100.0; + } + par = 1.0 / par; + if (par < 1.0 && par > 0.0) { + retn = 2; + par *= 0.9999; for (i = 0; i < m_numComponents; ++i) { - dir = FormulaMatrix[i][kspec] * (m_elemAbundancesGoal[i] - m_elemAbundances[i]); - if (fabs(dir) > 1.0E-10) { - if (dir > 0.0) { - if (saveDir < 0.0) { - goodSpec = FALSE; - break; - } - } else { - if (saveDir > 0.0) { - goodSpec = FALSE; - break; - } - } - saveDir = dir; + double tmp = m_molNumSpecies_old[i] + par * x[i]; + if (tmp > 0.0) { + m_molNumSpecies_old[i] = tmp; } else { - if (FormulaMatrix[i][kspec] != 0.) { - goodSpec = FALSE; - break; + if (SSPhase[i]) { + m_molNumSpecies_old[i] = 0.0; + } else { + m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001; } } } - if (goodSpec) { - its = 0; - xx = 0.0; - for (i = 0; i < m_numComponents; ++i) { - if (FormulaMatrix[i][kspec] != 0.0) { - xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / FormulaMatrix[i][kspec]; - its++; - } - } - if (its > 0) xx /= its; - m_molNumSpecies_old[kspec] += xx; - m_molNumSpecies_old[kspec] = MAX(m_molNumSpecies_old[kspec], 1.0E-10); - /* - * If we are dealing with a deleted species, then - * we need to reinsert it into the active list. - */ - if (kspec >= m_numSpeciesRdc) { - vcs_reinsert_deleted(kspec); - m_molNumSpecies_old[m_numSpeciesRdc - 1] = xx; - vcs_elab(); - goto L_CLEANUP; - } - vcs_elab(); - } - } - } - if (vcs_elabcheck(0)) { - retn = 1; - goto L_CLEANUP; - } - - for (i = 0; i < m_numElemConstraints; ++i) { - if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY || - (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) { - for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { - if (m_elemAbundances[i] > 0.0) { - if (FormulaMatrix[i][kspec] < 0.0) { - m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / FormulaMatrix[i][kspec] ; - if (m_molNumSpecies_old[kspec] < 0.0) { - m_molNumSpecies_old[kspec] = 0.0; - } - vcs_elab(); - break; - } - } - if (m_elemAbundances[i] < 0.0) { - if (FormulaMatrix[i][kspec] > 0.0) { - m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / FormulaMatrix[i][kspec]; - if (m_molNumSpecies_old[kspec] < 0.0) { - m_molNumSpecies_old[kspec] = 0.0; - } - vcs_elab(); - break; + } else { + for (i = 0; i < m_numComponents; ++i) { + double tmp = m_molNumSpecies_old[i] + x[i]; + if (tmp > 0.0) { + m_molNumSpecies_old[i] = tmp; + } else { + if (SSPhase[i]) { + m_molNumSpecies_old[i] = 0.0; + } else { + m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001; } } } } - } - if (vcs_elabcheck(1)) { - retn = 1; - goto L_CLEANUP; - } + + /* + * We have changed the element abundances. Calculate them again + */ + vcs_elab(); + /* + * We have changed the total moles in each phase. Calculate them again + */ + vcs_tmoles(); + + /* + * Try some ad hoc procedures for fixing the problem + */ + if (retn >= 2) { + /* + * First find a species whose adjustment is a win-win + * situation. + */ + for (kspec = 0; kspec < m_numSpeciesTot; kspec++) { + if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + continue; + } + saveDir = 0.0; + goodSpec = TRUE; + for (i = 0; i < m_numComponents; ++i) { + dir = FormulaMatrix[i][kspec] * (m_elemAbundancesGoal[i] - m_elemAbundances[i]); + if (fabs(dir) > 1.0E-10) { + if (dir > 0.0) { + if (saveDir < 0.0) { + goodSpec = FALSE; + break; + } + } else { + if (saveDir > 0.0) { + goodSpec = FALSE; + break; + } + } + saveDir = dir; + } else { + if (FormulaMatrix[i][kspec] != 0.) { + goodSpec = FALSE; + break; + } + } + } + if (goodSpec) { + its = 0; + xx = 0.0; + for (i = 0; i < m_numComponents; ++i) { + if (FormulaMatrix[i][kspec] != 0.0) { + xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / FormulaMatrix[i][kspec]; + its++; + } + } + if (its > 0) xx /= its; + m_molNumSpecies_old[kspec] += xx; + m_molNumSpecies_old[kspec] = MAX(m_molNumSpecies_old[kspec], 1.0E-10); + /* + * If we are dealing with a deleted species, then + * we need to reinsert it into the active list. + */ + if (kspec >= m_numSpeciesRdc) { + vcs_reinsert_deleted(kspec); + m_molNumSpecies_old[m_numSpeciesRdc - 1] = xx; + vcs_elab(); + goto L_CLEANUP; + } + vcs_elab(); + } + } + } + if (vcs_elabcheck(0)) { + retn = 1; + goto L_CLEANUP; + } + + for (i = 0; i < m_numElemConstraints; ++i) { + if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY || + (m_elType[i] == VCS_ELEM_TYPE_ABSPOS && m_elemAbundancesGoal[i] == 0.0)) { + for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { + if (m_elemAbundances[i] > 0.0) { + if (FormulaMatrix[i][kspec] < 0.0) { + m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / FormulaMatrix[i][kspec] ; + if (m_molNumSpecies_old[kspec] < 0.0) { + m_molNumSpecies_old[kspec] = 0.0; + } + vcs_elab(); + break; + } + } + if (m_elemAbundances[i] < 0.0) { + if (FormulaMatrix[i][kspec] > 0.0) { + m_molNumSpecies_old[kspec] -= m_elemAbundances[i] / FormulaMatrix[i][kspec]; + if (m_molNumSpecies_old[kspec] < 0.0) { + m_molNumSpecies_old[kspec] = 0.0; + } + vcs_elab(); + break; + } + } + } + } + } + if (vcs_elabcheck(1)) { + retn = 1; + goto L_CLEANUP; + } - /* - * For electron charges element types, we try positive deltas - * in the species concentrations to match the desired - * electron charge exactly. - */ - for (i = 0; i < m_numElemConstraints; ++i) { - double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i]; - if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) { - bool useZeroed = true; - for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { - if (dev < 0.0) { - if (FormulaMatrix[i][kspec] < 0.0) { - if (m_molNumSpecies_old[kspec] > 0.0) { - useZeroed = false; - } - } - } else { - if (FormulaMatrix[i][kspec] > 0.0) { - if (m_molNumSpecies_old[kspec] > 0.0) { - useZeroed = false; - } - } - } - } - for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { - if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) { + /* + * For electron charges element types, we try positive deltas + * in the species concentrations to match the desired + * electron charge exactly. + */ + for (i = 0; i < m_numElemConstraints; ++i) { + double dev = m_elemAbundancesGoal[i] - m_elemAbundances[i]; + if (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE && (fabs(dev) > 1.0E-300)) { + bool useZeroed = true; + for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (dev < 0.0) { if (FormulaMatrix[i][kspec] < 0.0) { - double delta = dev / FormulaMatrix[i][kspec] ; - m_molNumSpecies_old[kspec] += delta; - if (m_molNumSpecies_old[kspec] < 0.0) { - m_molNumSpecies_old[kspec] = 0.0; + if (m_molNumSpecies_old[kspec] > 0.0) { + useZeroed = false; + } + } + } else { + if (FormulaMatrix[i][kspec] > 0.0) { + if (m_molNumSpecies_old[kspec] > 0.0) { + useZeroed = false; } - vcs_elab(); - break; } } - if (dev > 0.0) { - if (FormulaMatrix[i][kspec] > 0.0) { - double delta = dev / FormulaMatrix[i][kspec] ; - m_molNumSpecies_old[kspec] += delta; - if (m_molNumSpecies_old[kspec] < 0.0) { - m_molNumSpecies_old[kspec] = 0.0; + } + for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { + if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) { + if (dev < 0.0) { + if (FormulaMatrix[i][kspec] < 0.0) { + double delta = dev / FormulaMatrix[i][kspec] ; + m_molNumSpecies_old[kspec] += delta; + if (m_molNumSpecies_old[kspec] < 0.0) { + m_molNumSpecies_old[kspec] = 0.0; + } + vcs_elab(); + break; + } + } + if (dev > 0.0) { + if (FormulaMatrix[i][kspec] > 0.0) { + double delta = dev / FormulaMatrix[i][kspec] ; + m_molNumSpecies_old[kspec] += delta; + if (m_molNumSpecies_old[kspec] < 0.0) { + m_molNumSpecies_old[kspec] = 0.0; + } + vcs_elab(); + break; } - vcs_elab(); - break; } } } } } - } - if (vcs_elabcheck(1)) { - retn = 1; - goto L_CLEANUP; - } - - L_CLEANUP: ; - vcs_tmoles(); -#ifdef DEBUG_MODE - l2after = 0.0; - for (i = 0; i < m_numElemConstraints; ++i) { - l2after += SQUARE(m_elemAbundances[i] - m_elemAbundancesGoal[i]); - } - l2after = sqrt(l2after/m_numElemConstraints); - if (vcs_debug_print_lvl >= 2) { - plogf(" --- Elem_Abund: Correct Initial " - " Final\n"); - for (i = 0; i < m_numElemConstraints; ++i) { - plogf(" --- "); plogf("%-2.2s", ElName[i].c_str()); - plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]); + if (vcs_elabcheck(1)) { + retn = 1; + goto L_CLEANUP; + } + + L_CLEANUP: ; + vcs_tmoles(); +#ifdef DEBUG_MODE + l2after = 0.0; + for (i = 0; i < m_numElemConstraints; ++i) { + l2after += SQUARE(m_elemAbundances[i] - m_elemAbundancesGoal[i]); + } + l2after = sqrt(l2after/m_numElemConstraints); + if (vcs_debug_print_lvl >= 2) { + plogf(" --- Elem_Abund: Correct Initial " + " Final\n"); + for (i = 0; i < m_numElemConstraints; ++i) { + plogf(" --- "); plogf("%-2.2s", ElName[i].c_str()); + plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]); + } + plogf(" --- Diff_Norm: %20.12E %20.12E\n", + l2before, l2after); } - plogf(" --- Diff_Norm: %20.12E %20.12E\n", - l2before, l2after); - } #endif - return retn; -} /* vcs_elcorr() ************************************************************/ + return retn; + } } diff --git a/Cantera/src/equil/vcs_elem_rearrange.cpp b/Cantera/src/equil/vcs_elem_rearrange.cpp index a11c25ffa..c3fac5c53 100644 --- a/Cantera/src/equil/vcs_elem_rearrange.cpp +++ b/Cantera/src/equil/vcs_elem_rearrange.cpp @@ -7,7 +7,6 @@ * $Date$ * $Revision$ */ - /* * Copywrite (2005) Sandia Corporation. Under the terms of * Contract DE-AC04-94AL85000 with Sandia Corporation, the diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index f5fa65f15..eeb78dfab 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -1,6 +1,6 @@ /** * @file vcs_inest.cpp - * Methods for obtaining a good initial guess + * Implementation methods for obtaining a good initial guess */ /* $Author$ * $Date$ @@ -13,36 +13,35 @@ * 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 "clockWC.h" +#include +#include +#include + namespace VCSnonideal { static char pprefix[20] = " --- vcs_inest: "; - /*****************************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ - - void VCS_SOLVE::inest(double *aw, double *sa, double *sm, - double *ss, double test) - - /************************************************************************** - * - * inest: - * - * Estimates equilibrium compositions. - * Algorithm covered in a section of Smith and Missen's Book. - * - * Linear programming module is based on using dbolm. - ***************************************************************************/ - { + + // Estimate equilibrium compositions + /* + * Estimates equilibrium compositions. + * Algorithm covered in a section of Smith and Missen's Book. + * + * Linear programming module is based on using dbolm. + * + * @param aw aw[i[ Mole fraction work space (ne in length) + * @param sa sa[j] = Gramm-Schmidt orthog work space (ne in length) + * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length) + * @param ss ss[j] = Gramm-Schmidt orthog work space (ne in length) + * @param test This is a small negative number. + */ + void VCS_SOLVE::inest(double * const aw, double * const sa, double * const sm, + double * const ss, double test) { int conv, k, lt, ikl, kspec, iph, irxn; double s; double s1 = 0.0; @@ -169,24 +168,24 @@ namespace VCSnonideal { /* **** CHEMICAL POTENTIALS OF BASIS ****************** */ /* ***************************************************************** */ /* - * Calculate TMoles and TPhMoles[] + * Calculate TMoles and m_tPhaseMoles_old[] */ vcs_tmoles(); /* - * TPhMoles1[] will consist of just the component moles + * m_tPhaseMoles_new[] will consist of just the component moles */ for (iph = 0; iph < NPhase; iph++) { - TPhMoles1[iph] = TPhInertMoles[iph] + 1.0E-20; + m_tPhaseMoles_new[iph] = TPhInertMoles[iph] + 1.0E-20; } for (kspec = 0; kspec < m_numComponents; ++kspec) { if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { - TPhMoles1[PhaseID[kspec]] += molNum[kspec]; + m_tPhaseMoles_new[PhaseID[kspec]] += molNum[kspec]; } } TMolesMultiphase = 0.0; for (iph = 0; iph < NPhase; iph++) { if (! VPhaseList[iph]->SingleSpecies) { - TMolesMultiphase += TPhMoles1[iph]; + TMolesMultiphase += m_tPhaseMoles_new[iph]; } } vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNum, nspecies); @@ -202,7 +201,7 @@ namespace VCSnonideal { if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { if (! SSPhase[kspec]) { iph = PhaseID[kspec]; - m_feSpecies_curr[kspec] += log(m_molNumSpecies_new[kspec] / TPhMoles[iph]); + m_feSpecies_curr[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]); } } else { m_molNumSpecies_new[kspec] = 0.0; @@ -227,8 +226,8 @@ namespace VCSnonideal { /* ********************************************************** */ vcs_dzero(VCS_DATA_PTR(DelTPhMoles), NPhase); for (iph = 0; iph < NPhase; iph++) { - xtphMax[iph] = log(TPhMoles1[iph] * 1.0E32); - xtphMin[iph] = log(TPhMoles1[iph] * 1.0E-32); + xtphMax[iph] = log(m_tPhaseMoles_new[iph] * 1.0E32); + xtphMin[iph] = log(m_tPhaseMoles_new[iph] * 1.0E-32); } for (irxn = 0; irxn < nrxn; ++irxn) { kspec = ir[irxn]; @@ -251,7 +250,7 @@ namespace VCSnonideal { * phase. * It cut diamond4.vin iterations down from 62 to 14. */ - m_deltaMolNumSpecies[kspec] = 0.5 * (TPhMoles1[iph] + TMolesMultiphase) + m_deltaMolNumSpecies[kspec] = 0.5 * (m_tPhaseMoles_new[iph] + TMolesMultiphase) * exp(-m_deltaGRxn_new[irxn]); for (k = 0; k < m_numComponents; ++k) { @@ -317,7 +316,7 @@ namespace VCSnonideal { } /* * We have a new w[] estimate, go get the - * TMoles and TPhMoles[] values + * TMoles and m_tPhaseMoles_old[] values */ vcs_tmoles(); if (lt > 0) goto finished; @@ -505,7 +504,7 @@ namespace VCSnonideal { if (vcs_debug_print_lvl >= 2) { plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix, vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), - VCS_DATA_PTR(TPhMoles))); + VCS_DATA_PTR(m_tPhaseMoles_old))); plogendl(); } #endif @@ -517,7 +516,7 @@ namespace VCSnonideal { m_VCount->T_Time_inest += tsecond; (m_VCount->T_Calls_Inest)++; return retn; - }/**** vcs_inest() ***********************************************************/ + } } diff --git a/Cantera/src/equil/vcs_prep.cpp b/Cantera/src/equil/vcs_prep.cpp index 012b2381e..7c226ce78 100644 --- a/Cantera/src/equil/vcs_prep.cpp +++ b/Cantera/src/equil/vcs_prep.cpp @@ -290,7 +290,7 @@ int VCS_SOLVE::vcs_prep(void) { vcs_dzero(&(DnPhase[0][0]), m_numSpeciesTot*NPhase); vcs_izero(&(PhaseParticipation[0][0]), m_numSpeciesTot*NPhase); vcs_dzero(VCS_DATA_PTR(DelTPhMoles), NPhase); - vcs_dzero(VCS_DATA_PTR(TPhMoles1), NPhase); + vcs_dzero(VCS_DATA_PTR(m_tPhaseMoles_new), NPhase); /* * Calculate the total number of moles in all phases. */ diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 69ff0afaa..fe07ebd62 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -153,7 +153,7 @@ int VCS_SOLVE::vcs_report(int iconv) (VPhaseList[i])->PhaseName.c_str()); } plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i], - TPhInertMoles[i] / TPhMoles[i], 0.0); + TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0); } } if (m_numSpeciesRdc != nspecies) { @@ -242,10 +242,10 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" %3d ", iphase); vcs_VolPhase *VPhase = VPhaseList[iphase]; plogf("%-12.12s |",VPhase->PhaseName.c_str()); - plogf("%10.3e |", TPhMoles[iphase]); - totalMoles += TPhMoles[iphase]; - if (TPhMoles[iphase] != VPhase->TotalMoles()) { - if (! vcs_doubleEqual(TPhMoles[iphase], VPhase->TotalMoles())) { + plogf("%10.3e |", m_tPhaseMoles_old[iphase]); + totalMoles += m_tPhaseMoles_old[iphase]; + if (m_tPhaseMoles_old[iphase] != VPhase->TotalMoles()) { + if (! vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->TotalMoles())) { plogf("We have a problem\n"); exit(-1); } @@ -281,7 +281,7 @@ int VCS_SOLVE::vcs_report(int iconv) */ g = vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), - VCS_DATA_PTR(TPhMoles)); + VCS_DATA_PTR(m_tPhaseMoles_old)); plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g); if (inertYes) plogf("\t\t(Inert species have standard free energy of zero)\n"); @@ -315,7 +315,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" %14.7E ", m_molNumSpecies_old[l]); plogf("%14.7E ", m_SSfeSpecies[l]); plogf("%14.7E ", log(ActCoeff[l])); - double tpmoles = TPhMoles[pid]; + double tpmoles = m_tPhaseMoles_old[pid]; double phi = phasePhi[pid]; double eContrib = phi * Charge[l] * Faraday_dim; double lx = 0.0; diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index f27e1b9e6..dc51d3491 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -142,8 +142,8 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void) } for (j = 0; j < NPhase; j++) { if (! (VPhaseList[j])->SingleSpecies) { - if (TPhMoles[j] > 0.0) - s -= SQUARE(dnPhase_irxn[j]) / TPhMoles[j]; + if (m_tPhaseMoles_old[j] > 0.0) + s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j]; } } if (s != 0.0) { @@ -195,13 +195,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void) */ if (dss != 0.0) { m_molNumSpecies_old[kspec] += dss; - TPhMoles[PhaseID[kspec]] += dss; + m_tPhaseMoles_old[PhaseID[kspec]] += dss; for (j = 0; j < m_numComponents; ++j) { m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j]; - TPhMoles[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j]; + m_tPhaseMoles_old[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j]; } m_molNumSpecies_old[k] = 0.0; - TPhMoles[PhaseID[k]] = 0.0; + m_tPhaseMoles_old[PhaseID[k]] = 0.0; #ifdef DEBUG_MODE plogf(" --- vcs_st2 Special section to delete "); plogf("%-12.12s", SpName[k].c_str()); diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 6e6ed42c3..4541f7292 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -131,8 +131,8 @@ namespace VCSnonideal { m_elemAbundances.resize(nelements, 0.0); m_elemAbundancesGoal.resize(nelements, 0.0); - TPhMoles.resize(nphase0, 0.0); - TPhMoles1.resize(nphase0, 0.0); + m_tPhaseMoles_old.resize(nphase0, 0.0); + m_tPhaseMoles_new.resize(nphase0, 0.0); DelTPhMoles.resize(nphase0, 0.0); TmpPhase.resize(nphase0, 0.0); TmpPhase2.resize(nphase0, 0.0); diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index f6c65671e..3e867b639 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -270,7 +270,13 @@ public: void vcs_redim_TP(void); void vcs_printChemPotUnits(int unitsFormat); - void vcs_elab(void); + //! Computes the current elemental abundances vector + /*! + * Computes the elemental abundances vector, m_elemAbundances[], and stores it + * back into the global structure + */ + void vcs_elab(); + int vcs_elabcheck(int ibound); void vcs_elabPhase(int iphase, double * const elemAbundPhase); int vcs_elcorr(double aa[], double x[]); @@ -294,7 +300,19 @@ public: #endif double vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles); - double vcs_GibbsPhase(int iphase, double *w, double *fe); + + //! Calculate the total dimensionless Gibbs free energy of a single phase + /*! + * -> Inert species are handled as if they had a standard free + * energy of zero and if they obeyed ideal solution/gas theory + * + * @param iphase ID of the phase + * @param w Species mole number vector for all species + * @param fe vector of partial molar free energies of all of the + * species + */ + double vcs_GibbsPhase(int iphase, const double * const w, + const double * const fe); double vcs_Gxs_phase_calc(vcs_VolPhase *Vphase, double *mf_PO); double vcs_Gxs_calc(int iphase); @@ -358,8 +376,25 @@ private: void prneav(void); void checkDelta1(double * const ds, double * const delTPhMoles, int kspec); #endif - void inest(double *aw, double *sa, double *sm, - double *ss, double test); + + //! Estimate equilibrium compositions + /*! + * Estimates equilibrium compositions. + * Algorithm covered in a section of Smith and Missen's Book. + * + * Linear programming module is based on using dbolm. + * + * @param aw aw[i[ Mole fraction work space (ne in length) + * @param sa sa[j] = Gramm-Schmidt orthog work space (ne in length) + * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length) + * @param ss ss[j] = Gramm-Schmidt orthog work space (ne in length) + * @param test This is a small negative number. + */ + void inest(double * const aw, double * const sa, double * const sm, + double * const ss, double test); + + + void vcs_SSPhase(void); double deltaG_Recalc_Rxn(int irxn, const double *const molNum, double * const ac, double * const mu_i); @@ -612,7 +647,7 @@ public: * * Length = number of phases */ - std::vector TPhMoles; + std::vector m_tPhaseMoles_old; //! total gmols of species in each phase in the tentative soln vector /*! @@ -621,7 +656,7 @@ public: * * Length = number of phases */ - std::vector TPhMoles1; + std::vector m_tPhaseMoles_new; //! Temporary vector of length NPhase std::vector TmpPhase; diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index b0fa33433..33c5b812e 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -814,7 +814,7 @@ namespace VCSnonideal { */ if (m_molNumSpecies_new[kspec] < 0.005 * TMoles) { iph = PhaseID[kspec]; - if (m_molNumSpecies_new[kspec] < (TPhMoles[iph] * 0.01)) { + if (m_molNumSpecies_new[kspec] < (m_tPhaseMoles_old[iph] * 0.01)) { #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { plogf(" --- Major species changed to minor: "); @@ -887,9 +887,9 @@ namespace VCSnonideal { */ dnPhase_irxn = DnPhase[irxn]; for (int iphase = 0; iphase < NPhase; iphase++) { - TPhMoles[iphase] += dnPhase_irxn[iphase] * dx; + m_tPhaseMoles_old[iphase] += dnPhase_irxn[iphase] * dx; } - TPhMoles[iph] = 0.0; + m_tPhaseMoles_old[iph] = 0.0; vcs_updateVP(0); /* * Recalcuate the chemical potentials, FE(), and the @@ -1086,7 +1086,7 @@ namespace VCSnonideal { * Calculate the tentative total mole numbers for each phase */ for (iph = 0; iph < NPhase; iph++) { - TPhMoles1[iph] = TPhMoles[iph] + DelTPhMoles[iph]; + m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + DelTPhMoles[iph]; } /* * Calculate the new chemical potentials using the tentative @@ -1119,10 +1119,10 @@ namespace VCSnonideal { if (printDetails) { plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n", vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old), - VCS_DATA_PTR(TPhMoles))); + VCS_DATA_PTR(m_tPhaseMoles_old))); plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E", vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr), - VCS_DATA_PTR(TPhMoles1))); + VCS_DATA_PTR(m_tPhaseMoles_new))); plogendl(); } @@ -1153,15 +1153,15 @@ namespace VCSnonideal { plogf("Norms of Delta G():%14.6E%14.6E\n", l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)), l2normdg(VCS_DATA_PTR(m_deltaGRxn_new))); - plogf(" Total moles of gas = %15.7E\n", TPhMoles[0]); + plogf(" Total moles of gas = %15.7E\n", m_tPhaseMoles_old[0]); if ((NPhase > 1) && (! (VPhaseList[1])->SingleSpecies)) { - plogf(" Total moles of liquid = %15.7E\n", TPhMoles[1]); + plogf(" Total moles of liquid = %15.7E\n", m_tPhaseMoles_old[1]); } else { plogf(" Total moles of liquid = %15.7E\n", 0.0); } plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n", vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr), - VCS_DATA_PTR(TPhMoles1))); + VCS_DATA_PTR(m_tPhaseMoles_new))); plogf(" -----------------------------------------------------"); plogendl(); } @@ -1211,15 +1211,15 @@ namespace VCSnonideal { plogf(" --- "); vcs_print_line("-", 50); for (iph = 0; iph < NPhase; iph++) { Vphase = VPhaseList[iph]; - plogf(" --- %18s = %15.7E\n", Vphase->PhaseName.c_str(), TPhMoles1[iph]); + plogf(" --- %18s = %15.7E\n", Vphase->PhaseName.c_str(), m_tPhaseMoles_new[iph]); } plogf(" "); vcs_print_line("-", 103); plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n", vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old), - VCS_DATA_PTR(TPhMoles))); + VCS_DATA_PTR(m_tPhaseMoles_old))); plogf(" --- Total New Dimensionless Gibbs Free Energy = %20.13E", vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr), - VCS_DATA_PTR(TPhMoles1))); + VCS_DATA_PTR(m_tPhaseMoles_new))); plogendl(); if (m_VCount->Its > 550) { plogf(" --- Troublesome solve"); @@ -1243,7 +1243,7 @@ namespace VCSnonideal { * we have already done this inside the FORCED * loop. */ - vcs_dcopy(VCS_DATA_PTR(TPhMoles), VCS_DATA_PTR(TPhMoles1), NPhase); + vcs_dcopy(VCS_DATA_PTR(m_tPhaseMoles_old), VCS_DATA_PTR(m_tPhaseMoles_new), NPhase); vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_molNumSpecies_new), m_numSpeciesRdc); vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc); vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc); @@ -1275,8 +1275,8 @@ namespace VCSnonideal { for (iph = 0; iph < NPhase; iph++) { Vphase = VPhaseList[iph]; if (!(Vphase->SingleSpecies)) { - if (TPhMoles[iph] != 0.0 && - TPhMoles[iph]/TMoles <= VCS_DELETE_PHASE_CUTOFF) { + if (m_tPhaseMoles_old[iph] != 0.0 && + m_tPhaseMoles_old[iph]/TMoles <= VCS_DELETE_PHASE_CUTOFF) { soldel = 1; for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (PhaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) { @@ -1922,8 +1922,8 @@ namespace VCSnonideal { m_molNumSpecies_new[kspec] = 1.0; } else { iph = PhaseID[kspec]; - if (TPhMoles[iph] != 0.0) { - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] / TPhMoles[iph]; + if (m_tPhaseMoles_old[iph] != 0.0) { + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] / m_tPhaseMoles_old[iph]; } else { /* * For MultiSpecies phases that are zeroed out, @@ -2166,12 +2166,12 @@ namespace VCSnonideal { *delta_ptr = dx; m_molNumSpecies_old[kspec] += dx; int iph = PhaseID[kspec]; - TPhMoles[iph] += dx; + m_tPhaseMoles_old[iph] += dx; for (j = 0; j < m_numComponents; ++j) { iph = PhaseID[j]; tmp = sc_irxn[j] * dx; m_molNumSpecies_old[j] += tmp; - TPhMoles[iph] += tmp; + m_tPhaseMoles_old[iph] += tmp; if (m_molNumSpecies_old[j] < 0.0) { m_molNumSpecies_old[j] = 0.0; } @@ -2193,7 +2193,7 @@ namespace VCSnonideal { * Zero out the concentration of a species. Make sure to conserve * elements and keep track of the total moles in all phases. * w[] - * TPhMoles[] + * m_tPhaseMoles_old[] * * return: * 1: succeeded @@ -2247,7 +2247,7 @@ namespace VCSnonideal { int irxn = kspec - m_numComponents; /* This is the noncomponent rxn index */ /* * Zero the concentration of the species. - * -> This zeroes w[kspec] and modifies TPhMoles[] + * -> This zeroes w[kspec] and modifies m_tPhaseMoles_old[] */ int retn = zero_species(kspec); if (! retn) { @@ -2275,7 +2275,7 @@ namespace VCSnonideal { /* * Adjust the total moles in a phase downwards. */ - Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles)); + Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old)); /* * Adjust the current number of active species and reactions counters @@ -2340,7 +2340,7 @@ namespace VCSnonideal { #endif /* * Set the species back to minor species status - * this adjusts m_molNumSpecies_old[] and TPhMoles[] + * this adjusts m_molNumSpecies_old[] and m_tPhaseMoles_old[] * HKM -> make this a relative mole number! */ dx = VCS_DELETE_SPECIES_CUTOFF * 10.; @@ -2353,7 +2353,7 @@ namespace VCSnonideal { } int iph = PhaseID[kspec]; vcs_VolPhase *Vphase = VPhaseList[iph]; - Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles)); + Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old)); /* * We may have popped a multispecies phase back * into existence. If we did, we have to check @@ -2421,8 +2421,8 @@ namespace VCSnonideal { /* * Zero out the total moles counters for the phase */ - TPhMoles[iph] = 0.0; - TPhMoles1[iph] = 0.0; + m_tPhaseMoles_old[iph] = 0.0; + m_tPhaseMoles_new[iph] = 0.0; DelTPhMoles[iph] = 0.0; /* @@ -2499,7 +2499,7 @@ namespace VCSnonideal { /* * Upload the state to the VP object */ - Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles), iph); + Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old), iph); } /* delete_multiphase() *****************************************************/ @@ -2536,8 +2536,8 @@ namespace VCSnonideal { vcs_deltag(0, true); for (iph = 0; iph < NPhase; iph++) { - if (TPhMoles[iph] > 0.0) - xtcutoff[iph] = log (TPhMoles[iph] / VCS_DELETE_SPECIES_CUTOFF); + if (m_tPhaseMoles_old[iph] > 0.0) + xtcutoff[iph] = log (m_tPhaseMoles_old[iph] / VCS_DELETE_SPECIES_CUTOFF); else xtcutoff[iph] = 0.0; } @@ -2546,7 +2546,7 @@ namespace VCSnonideal { * We are checking the equation: * * sum_u = sum_j_comp [ sigma_i_j * u_j ] - * = u_i_O + log((AC_i * W_i)/TPhMoles) + * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old) * * by first evaluating: * @@ -2571,14 +2571,14 @@ namespace VCSnonideal { for (irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) { kspec = ir[irxn]; iph = PhaseID[kspec]; - if (TPhMoles[iph] == 0.0) { + if (m_tPhaseMoles_old[iph] == 0.0) { if (m_deltaGRxn_new[irxn] < 0.0) { vcs_reinsert_deleted(kspec); npb++; } else { m_molNumSpecies_old[kspec] = 0.0; } - } else if (TPhMoles[iph] > 0.0) { + } else if (m_tPhaseMoles_old[iph] > 0.0) { if (m_deltaGRxn_new[irxn] < xtcutoff[iph]) { vcs_reinsert_deleted(kspec); npb++; @@ -2624,9 +2624,9 @@ namespace VCSnonideal { for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) { kspec = ir[irxn]; iph = PhaseID[kspec]; - if (TPhMoles[iph] > 0.0) { + if (m_tPhaseMoles_old[iph] > 0.0) { double maxDG = MIN(m_deltaGRxn_new[irxn], 300); - double dx = TPhMoles[iph] * exp(- maxDG); + double dx = m_tPhaseMoles_old[iph] * exp(- maxDG); retn = delta_species(kspec, &dx); } } @@ -2750,7 +2750,7 @@ namespace VCSnonideal { m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + al * m_deltaMolNumSpecies[kspec]; } for (iph = 0; iph < NPhase; iph++) { - TPhMoles1[iph] = TPhMoles[iph] + al * DelTPhMoles[iph]; + m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + al * DelTPhMoles[iph]; } vcs_updateVP(1); @@ -2865,7 +2865,7 @@ namespace VCSnonideal { * is nontrivial in size. */ iph = PhaseID[kspec]; - double tphmoles = TPhMoles[iph]; + double tphmoles = m_tPhaseMoles_old[iph]; double trphmoles = tphmoles / TMoles; if (trphmoles > VCS_DELETE_PHASE_CUTOFF) { m_deltaMolNumSpecies[kspec] = TMoles * VCS_SMALL_MULTIPHASE_SPECIES; @@ -2946,8 +2946,8 @@ namespace VCSnonideal { for (j = 0; j < NPhase; j++) { Vphase = VPhaseList[j]; if (! Vphase->SingleSpecies) { - if (TPhMoles[j] > 0.0) - s -= SQUARE(dnPhase_irxn[j]) / TPhMoles[j]; + if (m_tPhaseMoles_old[j] > 0.0) + s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j]; } } if (s != 0.0) { @@ -3053,16 +3053,16 @@ namespace VCSnonideal { */ if (dss != 0.0) { m_molNumSpecies_old[kspec] += dss; - TPhMoles[PhaseID[kspec]] += dss; + m_tPhaseMoles_old[PhaseID[kspec]] += dss; for (j = 0; j < m_numComponents; ++j) { m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j]; - TPhMoles[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j]; + m_tPhaseMoles_old[PhaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j]; } m_molNumSpecies_old[k] = 0.0; iph = PhaseID[k]; Vphase = VPhaseList[iph]; Vphase->Existence = 0; - TPhMoles[iph] = 0.0; + m_tPhaseMoles_old[iph] = 0.0; #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { plogf(" --- vcs_RxnStepSizes Special section to delete %s\n", @@ -3901,7 +3901,7 @@ namespace VCSnonideal { if (SSPhase[kspec]) { return VCS_SPECIES_ZEROEDSS; } else { - if (TPhMoles[iph] == 0.0) return VCS_SPECIES_ZEROEDPHASE; + if (m_tPhaseMoles_old[iph] == 0.0) return VCS_SPECIES_ZEROEDPHASE; else return VCS_SPECIES_ZEROEDMS; } } @@ -3973,7 +3973,7 @@ namespace VCSnonideal { * Check to see whether the current species is a major component * of its phase. If it is, it is a major component */ - if (m_molNumSpecies_old[kspec] > (TPhMoles[iph] * 0.1)) return VCS_SPECIES_MAJOR; + if (m_molNumSpecies_old[kspec] > (m_tPhaseMoles_old[iph] * 0.1)) return VCS_SPECIES_MAJOR; /* * Main check in the loop: * Check to see if there is a component with a mole number that is @@ -4243,9 +4243,9 @@ namespace VCSnonideal { } #endif if (kk <= 0) { - tPhMoles_ptr = VCS_DATA_PTR(TPhMoles); + tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old); } else { - tPhMoles_ptr = VCS_DATA_PTR(TPhMoles1); + tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new); } tlogMoles = VCS_DATA_PTR(TmpPhase); /* @@ -4535,20 +4535,20 @@ namespace VCSnonideal { double sum; vcs_VolPhase *Vphase; for (i = 0; i < NPhase; i++) { - TPhMoles[i] = TPhInertMoles[i]; + m_tPhaseMoles_old[i] = TPhInertMoles[i]; } for (i = 0; i < m_numSpeciesTot; i++) { if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { - TPhMoles[PhaseID[i]] += m_molNumSpecies_old[i]; + m_tPhaseMoles_old[PhaseID[i]] += m_molNumSpecies_old[i]; } } sum = 0.0; for (i = 0; i < NPhase; i++) { - sum += TPhMoles[i]; + sum += m_tPhaseMoles_old[i]; Vphase = VPhaseList[i]; // Took out because we aren't updating mole fractions in Vphase - // Vphase->TMoles = TPhMoles[i]; - if (TPhMoles[i] == 0.0) { + // Vphase->TMoles = m_tPhaseMoles_old[i]; + if (m_tPhaseMoles_old[i] == 0.0) { Vphase->Existence = 0; } else { if (TPhInertMoles[i] > 0.0) { @@ -4581,10 +4581,10 @@ namespace VCSnonideal { Vphase = VPhaseList[i]; if (place == 0) { Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), - VCS_DATA_PTR(TPhMoles), i); + VCS_DATA_PTR(m_tPhaseMoles_old), i); } else if (place == 1) { Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_new), - VCS_DATA_PTR(TPhMoles1), i); + VCS_DATA_PTR(m_tPhaseMoles_new), i); } else { plogf("we shouldn't be here\n"); exit(-1); diff --git a/Cantera/src/equil/vcs_species_thermo.cpp b/Cantera/src/equil/vcs_species_thermo.cpp index e142fb1a5..e1e463bc1 100644 --- a/Cantera/src/equil/vcs_species_thermo.cpp +++ b/Cantera/src/equil/vcs_species_thermo.cpp @@ -1,16 +1,17 @@ +/** + * @file vcs_species_thermo.cpp + * Implementation for the VCS_SPECIES_THERMO object. + */ /* * $Id$ */ - /* * Copywrite (2005) Sandia Corporation. Under the terms of * 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_species_thermo.h" @@ -20,14 +21,14 @@ #include "vcs_Exception.h" #include "vcs_internal.h" +#include +#include +#include + using namespace std; namespace VCSnonideal { -/***************************************************************************** - * - * constructor(): - */ VCS_SPECIES_THERMO::VCS_SPECIES_THERMO(int indexPhase, int indexSpeciesPhase) : @@ -436,7 +437,7 @@ double VCS_SOLVE::vcs_Gxs_calc(int iphase) { int kspec; double Gxs = 0.0, ac; - double totmol = TPhMoles[iphase]; + double totmol = m_tPhaseMoles_old[iphase]; vcs_VolPhase *Vphase = VPhaseList[iphase]; VCS_SPECIES_THERMO *ts_ptr;