From 6434e1edca526ceeb845cf3108a506045cfb3532 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Wed, 27 Mar 2013 00:35:58 +0000 Subject: [PATCH] Tightened the tolerances on the phaseStability calculation routine. Took out an unnecessary duplication of work. --- include/cantera/equil/vcs_VolPhase.h | 2 + src/equil/vcs_VolPhase.cpp | 16 +++--- src/equil/vcs_phaseStability.cpp | 69 +++++++++++++++++--------- src/equil/vcs_solve_phaseStability.cpp | 9 ++-- 4 files changed, 62 insertions(+), 34 deletions(-) diff --git a/include/cantera/equil/vcs_VolPhase.h b/include/cantera/equil/vcs_VolPhase.h index 1da3837a7..f05e44c41 100644 --- a/include/cantera/equil/vcs_VolPhase.h +++ b/include/cantera/equil/vcs_VolPhase.h @@ -861,6 +861,8 @@ private: * the value of the variable is the global vcs reaction number. Note, * that the global reaction number will go out of order when the species positions * are swapped. So, this number has to be recalculated. + * + * Length = number of species in phase */ std::vector creationGlobalRxnNumbers_; diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 59c0d427e..ed78a17b3 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -303,7 +303,11 @@ void vcs_VolPhase::resize(const size_t phaseNum, const size_t nspecies, for (size_t i = 0; i < nspecies; i++) { Xmol_[i] = 1.0/nspecies; creationMoleNumbers_[i] = 1.0/nspecies; - creationGlobalRxnNumbers_[i] = IndSpecies[i] - m_numElemConstraints; + if (IndSpecies[i] - m_numElemConstraints >= 0) { + creationGlobalRxnNumbers_[i] = IndSpecies[i] - m_numElemConstraints; + } else { + creationGlobalRxnNumbers_[i] = npos; + } } SS0ChemicalPotential.resize(nspecies, -1.0); @@ -344,9 +348,8 @@ void vcs_VolPhase::elemResize(const size_t numElemConstraints) m_numElemConstraints = numElemConstraints; } /***************************************************************************/ - -//! Evaluate activity coefficients -/*! +// Evaluate activity coefficients +/* * We carry out a calculation whenever UpTODate_AC is false. Specifically * whenever a phase goes zero, we do not carry out calculations on it. * @@ -1129,8 +1132,9 @@ void vcs_VolPhase::setCreationMoleNumbers(const double* const n_k, const std::vector &creationGlobalRxnNumbers) { vcs_dcopy(VCS_DATA_PTR(creationMoleNumbers_), n_k, m_numSpecies); - creationGlobalRxnNumbers_ = creationGlobalRxnNumbers; - + for (size_t k = 0; k < m_numSpecies; k++) { + creationGlobalRxnNumbers_[k] = creationGlobalRxnNumbers[k]; + } } /***************************************************************************/ diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index 36e3d7530..d6d58cc06 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -626,7 +626,16 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) return 0; } - +//====================================================================================================================== +// Main program to test whether a deleted phase should be brought +// back into existence +/* + * + * @param iph Phase id of the deleted phase + * + * So far this algorithm seems to be stable. I haven't run across any instance where it hasn't + * converged. However, it's probably only a matter of time + */ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) { @@ -635,18 +644,23 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) */ size_t kspec, irxn, k, i, kc, kc_spec; vcs_VolPhase* Vphase = m_VolPhaseList[iph]; + const size_t nsp = Vphase->nSpecies(); doublereal deltaGRxn; + int minNumberIterations = 3; + if (nsp <= 1) { + minNumberIterations = 1; + } // We will do a full newton calculation later, but for now, ... bool doSuccessiveSubstitution = true; double funcPhaseStability; - vector X_est(Vphase->nSpecies(), 0.0); - vector delFrac(Vphase->nSpecies(), 0.0); - vector E_phi(Vphase->nSpecies(), 0.0); - vector fracDelta_new(Vphase->nSpecies(), 0.0); - vector fracDelta_old(Vphase->nSpecies(), 0.0); - vector fracDelta_raw(Vphase->nSpecies(), 0.0); - vector creationGlobalRxnNumbers(Vphase->nSpecies(), npos); + vector X_est(nsp, 0.0); + vector delFrac(nsp, 0.0); + vector E_phi(nsp, 0.0); + vector fracDelta_new(nsp, 0.0); + vector fracDelta_old(nsp, 0.0); + vector fracDelta_raw(nsp, 0.0); + vector creationGlobalRxnNumbers(nsp, npos); vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); vector m_feSpecies_Deficient(m_numComponents, 0.0); @@ -668,7 +682,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) std::vector componentList; - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { kspec = Vphase->spGlobalIndexVCS(k); if (kspec < m_numComponents) { componentList.push_back(k); @@ -696,7 +710,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } #endif - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { if (fracDelta_new[k] < 1.0E-13) { fracDelta_new[k] = 1.0E-13; } @@ -715,7 +729,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) kc = componentList[i]; kc_spec = Vphase->spGlobalIndexVCS(kc); fracDelta_old[kc] = 0.0; - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { kspec = Vphase->spGlobalIndexVCS(k); if (kspec >= m_numComponents) { irxn = kspec - m_numComponents; @@ -726,7 +740,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) // Now, calculate the predicted mole fractions, X_est[k] double sumFrac = 0.0; - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { sumFrac += fracDelta_old[k]; } // Necessary because this can be identically zero. -> we need to fix this algorithm! @@ -734,7 +748,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) sumFrac = 1.0; } double sum_Xcomp = 0.0; - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { X_est[k] = fracDelta_old[k] / sumFrac; kc_spec = Vphase->spGlobalIndexVCS(k); if (kc_spec < m_numComponents) { @@ -796,7 +810,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) */ sum = 0.0; funcPhaseStability = sum_Xcomp - 1.0; - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { kspec = Vphase->spGlobalIndexVCS(k); if (kspec >= m_numComponents) { irxn = kspec - m_numComponents; @@ -818,7 +832,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) /* * Calculate the raw estimate of the new fracs */ - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { kspec = Vphase->spGlobalIndexVCS(k); double b = E_phi[k] / sum * (1.0 - sum_Xcomp); if (kspec >= m_numComponents) { @@ -834,7 +848,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) kc = componentList[i]; kc_spec = Vphase->spGlobalIndexVCS(kc); fracDelta_raw[kc] = 0.0; - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { kspec = Vphase->spGlobalIndexVCS(k); if (kspec >= m_numComponents) { irxn = kspec - m_numComponents; @@ -849,14 +863,14 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) * Now possibly dampen the estimate. */ doublereal sumADel = 0.0; - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { delFrac[k] = fracDelta_raw[k] - fracDelta_old[k]; sumADel += fabs(delFrac[k]); } normUpdate = vcs_l2norm(delFrac); dirProd = 0.0; - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { dirProd += fracDelta_old[k] * delFrac[k]; } bool crossedSign = false; @@ -883,7 +897,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } } - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) { damp = std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]), 1.0E-8/fabs(delFrac[k])); @@ -903,7 +917,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) damp = 0.000001; } - for (k = 0; k < Vphase->nSpecies(); k++) { + for (k = 0; k < nsp; k++) { fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]); } @@ -914,15 +928,24 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } #endif - if (normUpdate < 1.0E-5) { + if (normUpdate < 1.0E-5 * damp) { converged = true; + if (its < minNumberIterations) { + converged = false; + } } } if (converged) { - Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est), - VCS_STATECALC_PHASESTABILITY); + /* + * Save the final optimized stated back into the VolPhase object for later use + */ + Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est), VCS_STATECALC_PHASESTABILITY); + /* + * Save fracDelta for later use to initialize the problem better + * @TODO creationGlobalRxnNumbers needs to be calculated here and storred. + */ Vphase->setCreationMoleNumbers(VCS_DATA_PTR(fracDelta_new), creationGlobalRxnNumbers); } diff --git a/src/equil/vcs_solve_phaseStability.cpp b/src/equil/vcs_solve_phaseStability.cpp index a2871e77f..143a7a7a3 100644 --- a/src/equil/vcs_solve_phaseStability.cpp +++ b/src/equil/vcs_solve_phaseStability.cpp @@ -201,7 +201,7 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, { double test = -1.0E-10; bool usedZeroedSpecies; - std::vector phasePopPhaseIDs(0); + // std::vector phasePopPhaseIDs(0); int iStab = 0; std::vector sm(m_numElemConstraints*m_numElemConstraints, 0.0); @@ -212,8 +212,7 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, std::vector wx(m_numElemConstraints, 0.0); - vcs_basopt(false, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), - VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), + vcs_basopt(false, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test, &usedZeroedSpecies); vcs_evaluate_speciesType(); @@ -227,8 +226,8 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, vcs_printDeltaG(VCS_STATECALC_OLD); } vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); - phasePopPhaseIDs.clear(); - vcs_popPhaseID(phasePopPhaseIDs); + // phasePopPhaseIDs.clear(); + // vcs_popPhaseID(phasePopPhaseIDs); funcVal = vcs_phaseStabilityTest(iph); if (funcVal > 0.0) { iStab = 1;