From 7a99ca115f5d6c0d4b4b412c13f3ec3888adc09a Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 17 May 2008 23:52:13 +0000 Subject: [PATCH] Name changes on variables, only --- Cantera/src/equil/vcs_report.cpp | 6 +- Cantera/src/equil/vcs_rxnadj.cpp | 8 +- Cantera/src/equil/vcs_solve.cpp | 46 +++---- Cantera/src/equil/vcs_solve.h | 34 +++-- Cantera/src/equil/vcs_solve_TP.cpp | 153 ++++++++++++++++------- Cantera/src/equil/vcs_species_thermo.cpp | 4 +- 6 files changed, 165 insertions(+), 86 deletions(-) diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 189e1acfe..31ed89cfc 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -106,12 +106,12 @@ namespace VCSnonideal { * Calculate some quantities that may need updating */ vcs_tmoles(); - Vol = vcs_VolTotal(m_temperature, m_pressurePA, - VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM)); + m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA, + VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_PMVolumeSpecies)); plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature); plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA); - plogf("\t\tVolume = %15.5g m**3\n", Vol); + plogf("\t\tVolume = %15.5g m**3\n", m_totalVol); /* * -------- TABLE OF SPECIES IN DECREASING MOLE NUMBERS -------------- diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index 43247dfe5..dd9538f75 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -285,7 +285,7 @@ namespace VCSnonideal { /* * First the diagonal term of the Jacobian */ - s = dLnActCoeffdMolNum[kspec][kspec]; + s = m_dLnActCoeffdMolNum[kspec][kspec]; /* * Next, the other terms. Note this only a loop over the components * So, it's not too expensive to calculate. @@ -294,11 +294,11 @@ namespace VCSnonideal { if (!m_SSPhase[l]) { for (k = 0; k < m_numComponents; ++k) { if (m_phaseID[k] == m_phaseID[l]) { - s += sc_irxn[k] * sc_irxn[l] * dLnActCoeffdMolNum[k][l]; + s += sc_irxn[k] * sc_irxn[l] * m_dLnActCoeffdMolNum[k][l]; } } if (kph == m_phaseID[l]) { - s += sc_irxn[l] * (dLnActCoeffdMolNum[kspec][l] + dLnActCoeffdMolNum[l][kspec]); + s += sc_irxn[l] * (m_dLnActCoeffdMolNum[kspec][l] + m_dLnActCoeffdMolNum[l][kspec]); } } } @@ -334,7 +334,7 @@ namespace VCSnonideal { * -> This scatter calculation is carried out in the * vcs_VolPhase object. */ - Vphase->sendToVCSLnActCoeffJac(dLnActCoeffdMolNum.baseDataAddr()); + Vphase->sendToVCSLnActCoeffJac(m_dLnActCoeffdMolNum.baseDataAddr()); } } } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 2cd89e624..80f6bed2b 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -50,8 +50,8 @@ namespace VCSnonideal { m_tolmaj2(0.0), m_tolmin2(0.0), UnitsState(VCS_DIMENSIONAL_G), - UseActCoeffJac(0), - Vol(0.0), + m_useActCoeffJac(0), + m_totalVol(0.0), Faraday_dim(1.602e-19 * 6.022136736e26), m_VCount(0), vcs_debug_print_lvl(0), @@ -186,7 +186,7 @@ namespace VCSnonideal { CurrPhAC.resize(nphase0, 0); m_wtSpecies.resize(nspecies0, 0.0); m_chargeSpecies.resize(nspecies0, 0.0); - SpeciesThermo.resize(nspecies0, (VCS_SPECIES_THERMO *)0); + m_speciesThermoList.resize(nspecies0, (VCS_SPECIES_THERMO *)0); /* * Malloc Phase Info @@ -199,12 +199,12 @@ namespace VCSnonideal { /* * For Future expansion */ - UseActCoeffJac = true; - if (UseActCoeffJac ) { - dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0); + m_useActCoeffJac = true; + if (m_useActCoeffJac) { + m_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0); } - VolPM.resize(nspecies0, 0.0); + m_PMVolumeSpecies.resize(nspecies0, 0.0); /* * Malloc space for counters kept within vcs @@ -241,8 +241,8 @@ namespace VCSnonideal { } for (j = 0; j < nspecies; j++) { - delete SpeciesThermo[j]; - SpeciesThermo[j] = 0; + delete m_speciesThermoList[j]; + m_speciesThermoList[j] = 0; } delete m_VCount; m_VCount = 0; @@ -516,12 +516,12 @@ namespace VCSnonideal { * */ for (kspec = 0; kspec < nspecies; kspec++) { - if (SpeciesThermo[kspec] != NULL) { - delete SpeciesThermo[kspec]; + if (m_speciesThermoList[kspec] != NULL) { + delete m_speciesThermoList[kspec]; } VCS_SPECIES_THERMO *spf = pub->SpeciesThermo[kspec]; - SpeciesThermo[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO(); - if (SpeciesThermo[kspec] == NULL) { + m_speciesThermoList[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO(); + if (m_speciesThermoList[kspec] == NULL) { plogf(" duplMyselfAsVCS_SPECIES_THERMO returned an error!\n"); return VCS_PUB_BAD; } @@ -705,7 +705,7 @@ namespace VCSnonideal { for (int k = 0; k < Vphase->NVolSpecies; k++) { vcs_SpeciesProperties *sProp = Vphase->ListSpeciesPtr[k]; int kT = Vphase->IndSpecies[k]; - sProp->SpeciesThermo = SpeciesThermo[kT]; + sProp->SpeciesThermo = m_speciesThermoList[kT]; } } @@ -747,9 +747,9 @@ namespace VCSnonideal { /* * Copy the volume info */ - Vol = pub->Vol; - if (VolPM.size() != 0) { - vcs_dcopy(VCS_DATA_PTR(VolPM), VCS_DATA_PTR(pub->VolPM), nspecies); + m_totalVol = pub->Vol; + if (m_PMVolumeSpecies.size() != 0) { + vcs_dcopy(VCS_DATA_PTR(m_PMVolumeSpecies), VCS_DATA_PTR(pub->VolPM), nspecies); } /* @@ -767,7 +767,7 @@ namespace VCSnonideal { * * It's assumed we are solving the same problem. * - * @param pub Pointer to VCS_PROB that will be used to + * @param pub Pointer to VCS_PROdB that will be used to * initialize the current equilibrium problem */ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB *pub) { @@ -781,7 +781,7 @@ namespace VCSnonideal { m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat; m_doEstimateEquil = pub->iest; - Vol = pub->Vol; + m_totalVol = pub->Vol; m_tolmaj = pub->tolmaj; m_tolmin = pub->tolmin; @@ -906,8 +906,8 @@ namespace VCSnonideal { int k1 = 0; vcs_tmoles(); - Vol = vcs_VolTotal(m_temperature, m_pressurePA, - VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM)); + m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA, + VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_PMVolumeSpecies)); for (i = 0; i < m_numSpeciesTot; ++i) { /* @@ -930,12 +930,12 @@ namespace VCSnonideal { } pub->mf[i] = m_molNumSpecies_new[k1]; pub->m_gibbsSpecies[i] = m_feSpecies_curr[k1]; - pub->VolPM[i] = VolPM[k1]; + pub->VolPM[i] = m_PMVolumeSpecies[k1]; } pub->T = m_temperature; pub->PresPA = m_pressurePA; - pub->Vol = Vol; + pub->Vol = m_totalVol; int kT = 0; for (int iph = 0; iph < pub->NPhase; iph++) { vcs_VolPhase *pubPhase = pub->VPhaseList[iph]; diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 66c60f1ec..c37a93737 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -724,7 +724,19 @@ private: int delete_species(int kspec); void delete_multiphase(int iph); int delta_species(int kspec, double *delta_ptr); - void add_deleted(void); + + //! Provide an estimate for the deleted species in phases that + //! are not zeroed out + /*! + * Try to add back in all deleted species. An estimate of the kmol numbers + * are obtained and the species is added back into the equation system, + * into the old state vector. + * + * This routine is called at the end of the calculation, just before + * returning to the user. + */ + int vcs_add_all_deleted(); + int recheck_deleted(void); //! Alternative treatment for the update of a minor species @@ -1311,10 +1323,10 @@ public: * needed. It is not swapped wrt species * (unused atm) */ - DoubleStarStar dLnActCoeffdMolNum; + DoubleStarStar m_dLnActCoeffdMolNum; - //! This boolean indicates whether the activity coefficients for a phase - //! are current. + //! This boolean indicates whether the activity coefficients for a phase + //! are current. std::vector CurrPhAC; //! Molecular weight of each species @@ -1338,7 +1350,7 @@ public: /*! * SpeciesThermo[k] pointer to the thermo information for the kth species */ - std::vector SpeciesThermo; + std::vector m_speciesThermoList; //! Choice of Hessians /*! @@ -1346,16 +1358,20 @@ public: * Hessian based on Jacobian of the ln(ActCoeff) with respect to mole * numbers */ - int UseActCoeffJac; + int m_useActCoeffJac; - double Vol; /* Vol = Volume (m^3) */ + //! Total volume of all phases + /*! + * units are m^3 + */ + double m_totalVol; - //! Partialm molar volumes of the species + //! Partial molar volumes of the species /*! * units = mks (m^3/kmol) -determined by m_VCS_UnitsFormat * Length = number of species */ - std::vector VolPM; + std::vector m_PMVolumeSpecies; //! dimensionless value of Faraday's constant /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 35e1d4a9f..e65b486e4 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -1867,7 +1867,19 @@ namespace VCSnonideal { * Add back deleted species in non-zeroed phases. Estimate their * mole numbers. */ - add_deleted(); + npb = vcs_add_all_deleted(); + if (npb > 0) { + MajorSpeciesHaveConverged = true; + iti = 0; +#ifdef DEBUG_MODE + if (vcs_debug_print_lvl >= 1) { + plogf(" --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!"); + plogendl(); + } +#endif + goto L_MAINLOOP_ALL_SPECIES; + } + /* * Make sure the volume phase objects hold the same state and * information as the vcs object. This also update the Cantera objects @@ -2083,9 +2095,6 @@ namespace VCSnonideal { } return dx; } - - /*****************************************************************************/ - /*****************************************************************************/ /*****************************************************************************/ int VCS_SOLVE::delta_species(int kspec, double *delta_ptr) @@ -2284,9 +2293,10 @@ namespace VCSnonideal { * have to signal the calling code */ return (m_numRxnRdc == 0); - } /* delete_species() ********************************************************/ + } + /***************************************************************************/ - /**************************************************************************** + /* * * reinsert_deleted(): * @@ -2560,33 +2570,32 @@ namespace VCSnonideal { } } return npb; - } /* recheck_deleted() *******************************************************/ + } + /***********************************************************************************/ - - void VCS_SOLVE::add_deleted(void) - - /************************************************************************* - * - * Provide an estimate for the deleted species in phases that - * are not zeroed out - * - *************************************************************************/ - { + // Provide an estimate for the deleted species in phases that + // are not zeroed out + /* + * Try to add back in all deleted species. An estimate of the kmol numbers + * are obtained and the species is added back into the equation system, + * into the old state vector. + */ + int VCS_SOLVE::vcs_add_all_deleted() { int iph, kspec, retn; - if (m_numSpeciesRdc == m_numSpeciesTot) return; + if (m_numSpeciesRdc == m_numSpeciesTot) return 0; /* * Use the standard chemical potentials for the chemical potentials * of deleted species. Then, calculate Delta G for - * for formation reactions - * - * HKM Note: We need to update this step for nonunity activity - * coefficients. - * The formula will be fe = ff + RT * ln(actCoeff) - * where the activity coefficient is evaluated at - * ~ infinite dilution. + * for formation reactions. + * We are relying here on a old saved value of m_actCoeffSpecies_old[kspec] + * being sufficiently good. Note, we will recalculate everything at the + * end of the routine. */ for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { - m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; + iph = m_phaseID[kspec]; + m_feSpecies_curr[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) + - SpecLnMnaught[kspec] + + m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iph]); } /* * Recalculate the DeltaG's of the formation reactions for the @@ -2594,20 +2603,73 @@ namespace VCSnonideal { */ vcs_deltag(0, true); - for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) { kspec = m_indexRxnToSpecies[irxn]; iph = m_phaseID[kspec]; if (m_tPhaseMoles_old[iph] > 0.0) { - double maxDG = MIN(m_deltaGRxn_new[irxn], 300); + double maxDG = MIN(m_deltaGRxn_new[irxn], 300.0); + double dx = m_tPhaseMoles_old[iph] * exp(- maxDG); retn = delta_species(kspec, &dx); + if (retn == 0) { +#ifdef DEBUG_MODE + if (vcs_debug_print_lvl) { + plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", + m_speciesName[kspec].c_str(), kspec, dx); + } +#endif + if (dx > 1.0E-50) { + dx = 1.0E-50; + retn = delta_species(kspec, &dx); +#ifdef DEBUG_MODE + if (retn == 0) { + if (vcs_debug_print_lvl) { + plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n", + m_speciesName[kspec].c_str(), kspec, dx); + } + } +#endif + } + } +#ifdef DEBUG_MODE + if (vcs_debug_print_lvl >= 2) { + if (retn != 0) { + plogf(" --- add_deleted(): species %s added back in with mol number %g", + m_speciesName[kspec].c_str(), dx); + plogendl(); + } else { + plogf(" --- add_deleted(): species %s failed to be added back in"); + plogendl(); + } + } +#endif } } vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); vcs_deltag(0, true); + + retn = 0; + for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) { + kspec = m_indexRxnToSpecies[irxn]; + iph = m_phaseID[kspec]; + if (m_tPhaseMoles_old[iph] > 0.0) { + if (fabs(m_deltaGRxn_old[irxn]) > m_tolmin) { + retn++; +#ifdef DEBUG_MODE + if (vcs_debug_print_lvl >= 2) { + plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g", + m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec], + m_deltaGRxn_old[irxn]); + plogendl(); + } +#endif + } + } + } + return retn; } + /***********************************************************************************/ /* globalStepDamp * @@ -2806,7 +2868,7 @@ namespace VCSnonideal { * We update the matrix dlnActCoeffdmolNumber[][] at the * top of the loop, when necessary */ - if (UseActCoeffJac) { + if (m_useActCoeffJac) { vcs_CalcLnActCoeffJac(VCS_DATA_PTR(m_molNumSpecies_old)); } /************************************************************************ @@ -2855,7 +2917,8 @@ namespace VCSnonideal { #endif Vphase = m_VolPhaseList[iph]; int numSpPhase = Vphase->NVolSpecies; - m_deltaMolNumSpecies[kspec] = m_totalMolNum * 10.0 * VCS_DELETE_PHASE_CUTOFF / numSpPhase; + m_deltaMolNumSpecies[kspec] = + m_totalMolNum * 10.0 * VCS_DELETE_PHASE_CUTOFF / numSpPhase; } --(m_numRxnMinorZeroed); } else { @@ -2881,7 +2944,8 @@ namespace VCSnonideal { if (vcs_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], m_deltaGRxn_new[irxn], ANOTE); + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); } #endif continue; @@ -2897,7 +2961,8 @@ namespace VCSnonideal { if (vcs_debug_print_lvl >= 2) { plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); plogf(" %12.4E %12.4E %12.4E | %s\n", - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], m_deltaGRxn_new[irxn], ANOTE); + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec], + m_deltaGRxn_new[irxn], ANOTE); } #endif continue; @@ -2913,7 +2978,7 @@ namespace VCSnonideal { for (j = 0; j < m_numComponents; ++j) { if (!m_SSPhase[j]) { if (m_molNumSpecies_old[j] > 0.0) { - s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j]; + s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j]; } } } @@ -2930,7 +2995,7 @@ namespace VCSnonideal { * derivatives of the activity coefficients with respect to the * mole numbers, even in our diagonal approximation. */ - if (UseActCoeffJac) { + if (m_useActCoeffJac) { double s_old = s; s = vcs_Hessian_diag_adj(irxn, s_old); #ifdef DEBUG_MODE @@ -3074,7 +3139,7 @@ namespace VCSnonideal { } /*****************************************************************************/ - /************************************************************************** + /* * * vcs_deltag: * @@ -3086,7 +3151,7 @@ namespace VCSnonideal { * species I in this reaction. * * INPUT - * L = < 0 : Calculate reactions corresponding to + * L < 0 : Calculate reactions corresponding to * major noncomponent and zeroed species only * L = 0 : Do all noncomponent reactions, i, between * 0 <= i < irxnl @@ -3185,7 +3250,7 @@ namespace VCSnonideal { } } /* ************************************************* */ - /* **** MULTISPECIES PHASES WITH ZERO MOLES************ */ + /* **** MULTISPECIES PHASES WITH ZERO MOLES ******** */ /* ************************************************* */ /* * Massage the free energies for species with zero mole fractions @@ -3273,9 +3338,7 @@ namespace VCSnonideal { checkFinite(m_deltaGRxn_new[irxn]); } #endif - } /* vcs_deltag() ************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ + } /*****************************************************************************/ int VCS_SOLVE::vcs_basopt(int ifirst, double aw[], double sa[], double sm[], @@ -4450,7 +4513,7 @@ namespace VCSnonideal { m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } } else { - st_ptr = SpeciesThermo[kspec]; + st_ptr = m_speciesThermoList[kspec]; m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(actCoeff_ptr[kspec] * z[kspec]) - tlogMoles[m_phaseID[kspec]] - SpecLnMnaught[kspec]; } @@ -4684,14 +4747,14 @@ namespace VCSnonideal { SWAP(m_actCoeffSpecies_old[k1], m_actCoeffSpecies_old[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); + SWAP(m_speciesThermoList[k1], m_speciesThermoList[k2], st_tmp); + SWAP(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2], t1); for (j = 0; j < m_numElemConstraints; ++j) { SWAP(m_formulaMatrix[j][k1], m_formulaMatrix[j][k2], t1); } - if (UseActCoeffJac) { - vcs_switch2D(dLnActCoeffdMolNum.baseDataAddr(), k1, k2); + if (m_useActCoeffJac) { + vcs_switch2D(m_dLnActCoeffdMolNum.baseDataAddr(), k1, k2); } /* diff --git a/Cantera/src/equil/vcs_species_thermo.cpp b/Cantera/src/equil/vcs_species_thermo.cpp index 653df81f7..0258dd270 100644 --- a/Cantera/src/equil/vcs_species_thermo.cpp +++ b/Cantera/src/equil/vcs_species_thermo.cpp @@ -389,7 +389,7 @@ double VCS_SOLVE::vcs_Gxs_phase_calc(vcs_VolPhase *Vphase, double *mf_PO) if (Vphase->Activity_Coeff_Model != VCS_AC_CONSTANT) { for (kspec = 0; kspec < Vphase->NVolSpecies; kspec++) { kglob = Vphase->IndSpecies[kspec]; - ts_ptr = SpeciesThermo[kglob]; + ts_ptr = m_speciesThermoList[kglob]; ac = ts_ptr->eval_ac(kspec); Gxs += mf_PO[kspec] * log(ac); } @@ -434,7 +434,7 @@ double VCS_SOLVE::vcs_Gxs_calc(int iphase) for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (m_phaseID[kspec] == iphase) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - ts_ptr = SpeciesThermo[kspec]; + ts_ptr = m_speciesThermoList[kspec]; ac = ts_ptr->eval_ac(kspec); Gxs += m_molNumSpecies_old[kspec]/totmol * log(ac); } else {