From ef278fc8853ba07c7de63005e10fbaf3ce043bab Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 22 Jan 2009 22:27:07 +0000 Subject: [PATCH] Significant Update to the Equilibrium Solver Added a capability to solve for the mathematical condition for the birth of a multispecies phase. Added the logic into the algorithm. --- Cantera/src/equil/vcs_VolPhase.cpp | 1 + Cantera/src/equil/vcs_defs.h | 2 +- Cantera/src/equil/vcs_phaseStability.cpp | 695 +++++++++++++++++-- Cantera/src/equil/vcs_rxnadj.cpp | 10 +- Cantera/src/equil/vcs_solve.h | 43 +- Cantera/src/equil/vcs_solve_TP.cpp | 818 ++++++++++++----------- Cantera/src/thermo/HMWSoln.cpp | 48 +- Cantera/src/thermo/HMWSoln.h | 3 + 8 files changed, 1159 insertions(+), 461 deletions(-) diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 62bcb1811..143d0ecc1 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -1107,6 +1107,7 @@ namespace VCSnonideal { resize(VP_ID, nsp, nelem, PhaseName.c_str()); } TP_ptr->getMoleFractions(VCS_DATA_PTR(Xmol)); + fractionCreationDelta_ = Xmol; _updateMoleFractionDependencies(); /* diff --git a/Cantera/src/equil/vcs_defs.h b/Cantera/src/equil/vcs_defs.h index 0df5f2cf2..0f2d9e672 100644 --- a/Cantera/src/equil/vcs_defs.h +++ b/Cantera/src/equil/vcs_defs.h @@ -104,7 +104,7 @@ namespace VCSnonideal { //! Cutoff relative moles below which a phase is deleted //! from the equilibrium problem. #ifndef VCS_DELETE_PHASE_CUTOFF -#define VCS_DELETE_PHASE_CUTOFF 1.0e-11 +#define VCS_DELETE_PHASE_CUTOFF 1.0e-12 #endif //@} diff --git a/Cantera/src/equil/vcs_phaseStability.cpp b/Cantera/src/equil/vcs_phaseStability.cpp index 1f87da720..1a1dfa89c 100644 --- a/Cantera/src/equil/vcs_phaseStability.cpp +++ b/Cantera/src/equil/vcs_phaseStability.cpp @@ -22,90 +22,683 @@ using namespace std; namespace VCSnonideal { + + + // Utility function that evaluates whether a phase can be popped + // into existence + /* + * @param iphasePop id of the phase, which is currently zeroed, + * + * @return Returns true if the phase can come into existence + * and false otherwise. + */ + bool VCS_SOLVE::vcs_popPhasePossible(const int iphasePop) const { + + vcs_VolPhase *Vphase = m_VolPhaseList[iphasePop]; + +#ifdef DEBUG_MODE + int existence = Vphase->exists(); + if (existence > 0) { + printf("ERROR vcs_popPhasePossible called for a phase that exists!"); + exit(-1); + } +#endif + + /* + * section to do damping of the m_deltaMolNumSpecies[] + */ + for (int k = 0; k < Vphase->nSpecies(); k++) { + int kspec = Vphase->spGlobalIndexVCS(k); + int irxn = kspec - m_numComponents; + if (irxn >= 0) { + for (int j = 0; j < m_numComponents; ++j) { + if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { + double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; + if (stoicC != 0.0) { + double negChangeComp = - stoicC * 1.0; + if (negChangeComp > 0.0) { + // TODO: We may have to come up with a tolerance here + if (m_molNumSpecies_old[j] <= 1.0E-300) { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 3) { + plogf(" --- vcs_popPhasePosssible() Phase %d (%s) can't be popped\n", iphasePop, + Vphase->PhaseName.c_str()); + plogf(" --- Component %d (%s)will go negative\n", j, m_speciesName[j].c_str()); + } +#endif + return false; + } + } + } + } + } + } + } + return true; + } + + // Decision as to whether a phase pops back into existence + /* + * @return returns the phase id of the phase that pops back into + * existence. Returns -1 if there are no phases + */ + int VCS_SOLVE::vcs_popPhaseID() { + int iphasePop = -1; + int iph; + int irxn, kspec; + doublereal FephaseMax = -1.0E30; + doublereal Fephase = -1.0E30; + vcs_VolPhase *Vphase = 0; + - int VCS_SOLVE::vcs_phaseStabilityTest(const int iph) { +#ifdef DEBUG_MODE + char anote[128]; + if (m_debug_print_lvl >= 2) { + plogf(" --- vcs_popPhaseID() called\n"); + plogf(" --- Phase Status F_e MoleNum\n"); + plogf(" --------------------------------------------------------------\n"); + } +#endif + for (iph = 0; iph < m_numPhases; iph++) { + Vphase = m_VolPhaseList[iph]; + int existence = Vphase->exists(); + strcpy(anote, ""); + if (existence > 0) { + +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- %18s %5d NA %11.3e\n", + Vphase->PhaseName.c_str(), + existence, + m_tPhaseMoles_old[iph]); + } +#endif + } else { + if (Vphase->m_singleSpecies) { + /*********************************************************************** + * + * Single Phase Stability Resolution + * + ***********************************************************************/ + kspec = Vphase->spGlobalIndexVCS(0); + irxn = kspec - m_numComponents; + doublereal deltaGRxn = m_deltaGRxn_old[irxn]; + Fephase = exp(-deltaGRxn) - 1.0; + if (Fephase > 0.0) { +#ifdef DEBUG_MODE + strcpy(anote," (ready to be birthed)"); +#endif + if (Fephase > FephaseMax) { + iphasePop = iph; + FephaseMax = Fephase; +#ifdef DEBUG_MODE + strcpy(anote," (chosen to be birthed)"); +#endif + } + } + +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- %18s %5d NA %11.3g %s\n", + Vphase->PhaseName.c_str(), + existence, Fephase, + m_tPhaseMoles_old[iph], anote); + } +#endif + + } else { + /*********************************************************************** + * + * MultiSpecies Phase Stability Resolution + * + ***********************************************************************/ + if (vcs_popPhasePossible(iph)) { + Fephase = vcs_phaseStabilityTest(iph); + if (Fephase > 0.0) { + if (Fephase > FephaseMax) { + iphasePop = iph; + FephaseMax = Fephase; + } + } else { + if (Fephase > FephaseMax) { + FephaseMax = Fephase; + } + } +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- %18s %5d %11.3g %11.3g\n", + Vphase->PhaseName.c_str(), + existence, Fephase, + m_tPhaseMoles_old[iph]); + } +#endif + } else { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- %18s %5d blocked %11.3g\n", + Vphase->PhaseName.c_str(), + existence, m_tPhaseMoles_old[iph]); + } +#endif + } + } + } + } + + +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --------------------------------------------------------------\n"); + } +#endif + return iphasePop; + } + + // Calculates the deltas of the reactions due to phases popping + // into existence + /* + * @param iphasePop Phase id of the phase that will come into existence + * + * Output + * ------- + * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers + * to the irxn'th species + * formation reaction. This adjustment + * is for species + * irxn + M, where M is the number + * of components. + * + * @return Returns an int representing the status of the step + * - 0 : normal return + * - 1 : A single species phase species has been zeroed out + * in this routine. The species is a noncomponent + * - 2 : Same as one but, the zeroed species is a component. + * - 3 : Nothing was done because the phase couldn't be birthed + * because a needed component is zero. + */ + int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const int iphasePop) { + vcs_VolPhase *Vphase = m_VolPhaseList[iphasePop]; + // Identify the first species in the phase + int kspec = Vphase->spGlobalIndexVCS(0); + // Identify the formation reaction for that species + int irxn = kspec - m_numComponents; + + doublereal s; + int j, k; + // Calculate the initial moles of the phase being born. + // Here we set it to 10x of the value which would cause the phase to be + // zeroed out within the algorithm. We may later adjust the value. + doublereal tPhaseMoles = 10. * m_totalMolNum * VCS_DELETE_PHASE_CUTOFF; + + +#ifdef DEBUG_MODE + int existence = Vphase->exists(); + if (existence > 0) { + printf("ERROR vcs_popPhaseRxnStepSizes called for a phase that exists!"); + exit(-1); + } + char anote[256]; + if (m_debug_print_lvl >= 2) { + plogf(" --- vcs_popPhaseRxnStepSizes() called to pop phase %s %d into existence\n", + Vphase->PhaseName.c_str(), iphasePop); + } +#endif + // Section for a single-species phase + // + if (Vphase->m_singleSpecies) { + s = 0.0; + double *dnPhase_irxn = m_deltaMolNumPhase[irxn]; + 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]; + } + } + } + for (j = 0; j < m_numPhases; j++) { + Vphase = m_VolPhaseList[j]; + if (! Vphase->m_singleSpecies) { + if (m_tPhaseMoles_old[j] > 0.0) + s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j]; + } + } + if (s != 0.0) { + double s_old = s; + s = vcs_Hessian_diag_adj(irxn, s_old); +#ifdef DEBUG_MODE + if (s_old != s) { + sprintf(anote, "Normal calc: diag adjusted from %g " + "to %g due to act coeff", s_old, s); + } +#endif + m_deltaMolNumSpecies[kspec] = -m_deltaGRxn_new[irxn] / s; + } else { + // Ok, s is equal to zero. We can not apply a sophisticated theory + // to birth the phase. Just pick a small delta and go with it. + m_deltaMolNumSpecies[kspec] = tPhaseMoles; + } + + /* + * section to do damping of the m_deltaMolNumSpecies[] + */ + for (j = 0; j < m_numComponents; ++j) { + double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; + if (stoicC != 0.0) { + if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { + double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec]; + if (negChangeComp > m_molNumSpecies_old[j]) { + if (m_molNumSpecies_old[j] > 0.0) { +#ifdef DEBUG_MODE + sprintf(anote, "Delta damped from %g " + "to %g due to component %d (%10s) going neg", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str()); +#endif + m_deltaMolNumSpecies[kspec] = - 0.5 * m_molNumSpecies_old[j] / stoicC; + } else { +#ifdef DEBUG_MODE + sprintf(anote, "Delta damped from %g " + "to %g due to component %d (%10s) zero", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str()); +#endif + m_deltaMolNumSpecies[kspec] = 0.0; + } + } + } + } + } + // Implement a damping term that limits m_deltaMolNumSpecies to the size of the mole number + if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) { +#ifdef DEBUG_MODE + sprintf(anote, "Delta damped from %g " + "to %g due to %s going negative", m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str()); +#endif + m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec]; + } + + + } else { + vector fracDelta(Vphase->nSpecies()); + vector X_est(Vphase->nSpecies()); + fracDelta = Vphase->fractionCreationDeltas(); + + double sumFrac = 0.0; + for (k = 0; k < Vphase->nSpecies(); k++) { + sumFrac += fracDelta[k]; + } + for (k = 0; k < Vphase->nSpecies(); k++) { + X_est[k] = fracDelta[k] / sumFrac; + } + + doublereal deltaMolNumPhase = tPhaseMoles; + doublereal damp = 1.0; + m_deltaGRxn_tmp = m_molNumSpecies_old; + double * molNumSpecies_tmp = DATA_PTR(m_deltaGRxn_tmp); + + + for (k = 0; k < Vphase->nSpecies(); k++) { + kspec = Vphase->spGlobalIndexVCS(k); + double delmol = deltaMolNumPhase * X_est[k]; + irxn = kspec - m_numComponents; + if (kspec >= m_numComponents) { + for (j = 0; j < m_numComponents; ++j) { + double stoicC = m_stoichCoeffRxnMatrix[irxn][j]; + if (stoicC != 0.0) { + if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { + molNumSpecies_tmp[j] += stoicC * delmol; + } + } + } + } + } + + doublereal ratioComp = 0.0; + for (j = 0; j < m_numComponents; ++j) { + double deltaJ = m_molNumSpecies_old[j] - molNumSpecies_tmp[j]; + if (molNumSpecies_tmp[j] < 0.0) { + ratioComp = 1.0; + if (deltaJ > 0.0) { + double delta0 = m_molNumSpecies_old[j]; + double dampj = delta0 / deltaJ * 0.9; + if (dampj < damp) { + damp = dampj; + } + } + } else { + if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { + int jph = m_phaseID[j]; + if ((jph != iphasePop) && (!m_SSPhase[j])) { + double fdeltaJ = fabs(deltaJ); + if ( m_molNumSpecies_old[j] > 0.0) { + ratioComp = MAX(ratioComp, fdeltaJ/ m_molNumSpecies_old[j]); + } + } + } + } + } + + // We may have greatly underestimated the deltaMoles for the phase pop + // Here we create a damp > 1 to account for this possibility. + // We adjust upwards to make sure that a component in an existing multispecies + // phase is modified by a factor of 1/1000. + if (ratioComp > 1.0E-30) { + if (ratioComp < 0.001) { + damp = 0.001 / ratioComp; + } + } + + + if (damp <= 1.0E-6) { + return 3; + } + + for (k = 0; k < Vphase->nSpecies(); k++) { + kspec = Vphase->spGlobalIndexVCS(k); + if (kspec < m_numComponents) { + m_speciesStatus[kspec] = VCS_SPECIES_COMPONENT; + } else { + m_deltaMolNumSpecies[kspec] = deltaMolNumPhase * X_est[k] * damp; + if (X_est[k] > 1.0E-3) { + m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; + } else { + m_speciesStatus[kspec] = VCS_SPECIES_MINOR; + } + } + } + + } + + return 0; + } + + + // + + double VCS_SOLVE::vcs_phaseStabilityTest(const int iph) { /* * We will use the _new state calc here */ - int kspec, irxn, k; + int kspec, irxn, k, i, kc, kc_spec; vcs_VolPhase *Vphase = m_VolPhaseList[iph]; - double deltaGRxn; + doublereal deltaGRxn; + // We will do a full newton calculation later, but for now, ... bool doSuccessiveSubstitution = true; - int res = 0; + 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 X_est(Vphase->nSpecies(), 0.0); - vector X_est_old(Vphase->nSpecies(), 0.0); - vector delX(Vphase->nSpecies(), 0.0); - vector E_phi(Vphase->nSpecies(), 0.0); - double damp = 1.0; - double normUpdate = 1.0; - double normUpdateOld = 1.0; + vector m_feSpecies_Deficient(m_numComponents, 0.0); + doublereal damp = 1.0; + doublereal dampOld = 1.0; + doublereal normUpdate = 1.0; + doublereal normUpdateOld = 1.0; + doublereal sum = 0.0; + doublereal dirProd = 0.0; + doublereal dirProdOld = 0.0; // get the activity coefficients Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new)); - - + + // Get the storred estimate for the composition of the phase if + // it gets created + fracDelta_new = Vphase->fractionCreationDeltas(); + +#ifdef DEBUG_MODE + if (m_temperature < 380.) { + // fracDelta_new[0] = 0.8; + // fracDelta_new[1] = 0.1; + // fracDelta_new[2] = 1.0E-8; + // fracDelta_new[3] = 0.1; + //fracDelta_new[4] = 1.0E-8; + } + if (m_temperature < 390.) { + printf("we are here\n"); + } +#endif + + bool oneIsComponent = false; + std::vector componentList; + + for (k = 0; k < Vphase->nSpecies(); k++) { + kspec = Vphase->spGlobalIndexVCS(k); + if (kspec < m_numComponents) { + oneIsComponent = true; + componentList.push_back(k); + } + } + + for (k = 0; k < m_numComponents; k++) { + m_feSpecies_Deficient[k] = m_feSpecies_old[k]; + } + normUpdate = 0.1 * vcs_l2norm(fracDelta_new); + damp = 1.0E-2; + if (doSuccessiveSubstitution) { - for (int its = 0; its < 20; its++) { - - + bool converged = false; + for (int its = 0; its < 200 && (!converged); its++) { + + dampOld = damp; normUpdateOld = normUpdate; - for (k = 0; k < Vphase->nSpecies(); k++) { - X_est_old[k] = X_est[k]; + fracDelta_old = fracDelta_new; + dirProdOld = dirProd; + + // Given a set of fracDelta's, we calculate the fracDelta's + // for the component species, if any + for (i = 0; i < (int) componentList.size(); i++) { + kc = componentList[i]; + kc_spec = Vphase->spGlobalIndexVCS(kc); + fracDelta_old[kc] = 0.0; + for (k = 0; k < Vphase->nSpecies(); k++) { + kspec = Vphase->spGlobalIndexVCS(k); + irxn = kspec - m_numComponents; + if (irxn >= 0) { + fracDelta_old[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_old[k]; + } + } } - double poly = -1.0; + // Now, calculate the predicted mole fractions, X_est[k] + double sumFrac = 0.0; for (k = 0; k < Vphase->nSpecies(); k++) { - kspec = Vphase->spGlobalIndexVCS(k); - irxn = kspec - m_numComponents; - deltaGRxn = m_deltaGRxn_old[irxn]; - // We may need to look at deltaGRxn for components! - if (irxn >= 0) { - if (deltaGRxn > 50.0) deltaGRxn = 50.0; - if (deltaGRxn < -50.0) deltaGRxn = -50.0; - E_phi[k] = exp(-deltaGRxn)/m_actCoeffSpecies_new[kspec]; - poly += E_phi[k]; + sumFrac += fracDelta_old[k]; + } + double sum_Xcomp = 0.0; + for (k = 0; k < Vphase->nSpecies(); k++) { + X_est[k] = fracDelta_old[k] / sumFrac; + kc_spec = Vphase->spGlobalIndexVCS(k); + if (kc_spec < m_numComponents) { + sum_Xcomp += X_est[k]; } } - double sum = poly + 1.0; - - for (k = 0; k < Vphase->nSpecies(); k++) { - delX[k] = E_phi[k]/sum - X_est_old[k]; - } - normUpdate = vcs_l2norm(delX); - - // Figure out the damping coefficient - double ratio = normUpdate / normUpdateOld; - if (ratio < 0.4) { - damp = 1.0; - } else if (ratio > 1.0) { - damp = 0.03; - } else { - damp = 0.1; - } - + + /* + * Feed the newly formed estimate of the mole fractions back into the + * ThermoPhase object + */ + Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est), VCS_STATECALC_PHASESTABILITY); + + /* + * get the activity coefficients + */ + Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, VCS_DATA_PTR(m_actCoeffSpecies_new)); + + /* + * first Calculate altered chemical potentials for component species + * belonging to this phase. + */ + for (i = 0; i < (int) componentList.size(); i++) { + kc = componentList[i]; + kc_spec = Vphase->spGlobalIndexVCS(kc); + if ( X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) { + m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + + log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]); + } else { + m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec] + + log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF); + } + } + + for (i = 0; i < (int) componentList.size(); i++) { + kc = componentList[i]; + kc_spec = Vphase->spGlobalIndexVCS(kc); + + for (k = 0; k < Vphase->nSpecies(); k++) { + kspec = Vphase->spGlobalIndexVCS(k); + irxn = kspec - m_numComponents; + if (irxn >= 0) { + if (i == 0) { + m_deltaGRxn_Deficient[irxn] = m_deltaGRxn_old[irxn]; + } + double *dtmp_ptr = m_stoichCoeffRxnMatrix[irxn]; + if (dtmp_ptr[kc_spec] != 0.0) { + m_deltaGRxn_Deficient[irxn] += + dtmp_ptr[kc_spec] * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]); + } + } + + } + } + + /* + * Calculate the E_phi's + */ + sum = 0.0; + funcPhaseStability = sum_Xcomp - 1.0; for (k = 0; k < Vphase->nSpecies(); k++) { - X_est[k] = X_est_old[k] + damp * delX[k]; + kspec = Vphase->spGlobalIndexVCS(k); + irxn = kspec - m_numComponents; + if (irxn >= 0) { + deltaGRxn = m_deltaGRxn_Deficient[irxn]; + if (deltaGRxn > 50.0) deltaGRxn = 50.0; + if (deltaGRxn < -50.0) deltaGRxn = -50.0; + E_phi[k] = std::exp(-deltaGRxn) / m_actCoeffSpecies_new[kspec]; + sum += E_phi[k]; + funcPhaseStability += E_phi[k]; + } else { + E_phi[k] = 0.0; + } + } + + /* + * Calculate the raw estimate of the new fracs + */ + for (k = 0; k < Vphase->nSpecies(); k++) { + kspec = Vphase->spGlobalIndexVCS(k); + irxn = kspec - m_numComponents; + double b = E_phi[k] / sum * (1.0 - sum_Xcomp); + if (irxn >= 0) { + fracDelta_raw[k] = (sumFrac - fracDelta_old[k]) * b / (1.0 - b); + } + } + + + // Given a set of fracDelta's, we calculate the fracDelta's + // for the component species, if any + for (i = 0; i < (int) componentList.size(); i++) { + kc = componentList[i]; + kc_spec = Vphase->spGlobalIndexVCS(kc); + fracDelta_raw[kc] = 0.0; + for (k = 0; k < Vphase->nSpecies(); k++) { + kspec = Vphase->spGlobalIndexVCS(k); + irxn = kspec - m_numComponents; + if (irxn >= 0) { + fracDelta_raw[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_raw[k]; + } + } + } + + + + /* + * Now possibly dampen the estimate. + */ + doublereal sumADel = 0.0; + for (k = 0; k < Vphase->nSpecies(); 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++) { + dirProd += fracDelta_old[k] * delFrac[k]; + } + bool crossedSign = false; + if (dirProd * dirProdOld < 0.0) { + crossedSign = true; + } + + + damp = 0.5; + if (dampOld < 0.25) { + damp = 2.0 * dampOld; + } + if (crossedSign) { + if (normUpdate *1.5 > normUpdateOld) { + damp = 0.5 * dampOld; + } else if (normUpdate *2.0 > normUpdateOld) { + damp = 0.8 * dampOld; + } + } else { + if (normUpdate > normUpdateOld * 2.0) { + damp = 0.6 * dampOld; + } else if (normUpdate > normUpdateOld * 1.2) { + damp = 0.9 * dampOld; + } } for (k = 0; k < Vphase->nSpecies(); k++) { - kspec = Vphase->spGlobalIndexVCS(k); - m_molNumSpecies_new[kspec] = X_est[k]; + if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) { + damp = MAX(0.3*fabs(fracDelta_old[k]) / fabs( delFrac[k]), 1.0E-8/fabs( delFrac[k])); + } + if (delFrac[k] < 0.0) { + if (2.0 * damp * (-delFrac[k]) > fracDelta_old[k]) { + damp = fracDelta_old[k] / (2.0 * (-delFrac[k])); + } + } + if (delFrac[k] > 0.0) { + if (2.0 * damp * delFrac[k] > fracDelta_old[k]) { + damp = fracDelta_old[k] / (2.0 * delFrac[k]); + } + } + } + if (damp < 0.000001) { + damp = 0.000001; + } + + for (k = 0; k < Vphase->nSpecies(); k++) { + fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]); + } + + + + + if (normUpdate < 1.0E-5) { + converged = true; } - Vphase->setMolesFromVCS(VCS_STATECALC_NEW); - } + + if (converged) { + Vphase->setMoleFractionsState(0.0, VCS_DATA_PTR(X_est), VCS_STATECALC_PHASESTABILITY); + Vphase->setFractionCreationDeltas( VCS_DATA_PTR(fracDelta_new)); + } + + } else { printf("not done yet\n"); exit(-1); } - return res; + return funcPhaseStability; } } diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index 80e5e8301..6bd036a82 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -27,10 +27,12 @@ namespace VCSnonideal { * * Output * ------- - * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers + * m_deltaMolNumSpecies[kspec] : reaction adjustments, where irxn refers * to the irxn'th species - * formation reaction. This adjustment is for species - * irxn + M, where M is the number of components. + * formation reaction. This adjustment + * is for species + * irxn + M, where M is the number + * of components. * * Special branching occurs sometimes. This causes the component basis * to be reevaluated @@ -127,7 +129,7 @@ namespace VCSnonideal { } } else { /********************************************************************/ - /************************* REGULAR PROCESSING ************/ + /************************* REGULAR PROCESSING ***********************/ /********************************************************************/ /* * First take care of cases where we want to bail out diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 498b5010e..4b28f146a 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -498,7 +498,37 @@ public: * */ void vcs_updateVP(const int stateCalc); - + + + //! Utility function that evaluates whether a phase can be popped + //! into existence + /*! + * @param iphasePop id of the phase, which is currently zeroed, + * + * @return Returns true if the phase can come into existence + * and false otherwise. + */ + bool vcs_popPhasePossible(const int iphasePop) const; + + //! Decision as to whether a phase pops back into existence + /*! + * @return returns the phase id of the phase that pops back into + * existence. Returns -1 if there are no phases + */ + int vcs_popPhaseID(); + + //! Calculates the deltas of the reactions due to phases popping + //! into existence + /*! + * @param iphasePop Phase id of the phase that will come into existence + * + * @return Returns an int representing the status of the step + * - 0 : normal return + * - 1 : A single species phase species has been zeroed out + * in this routine. The species is a noncomponent + * - 2 : Same as one but, the zeroed species is a component. + */ + int vcs_popPhaseRxnStepSizes(const int iphasePop); //! Calculates formation reaction step sizes. /*! @@ -633,8 +663,13 @@ public: double vcs_birthGuess(const int kspec); - - int vcs_phaseStabilityTest(const int iph); + //! Main program to test whether a deleted phase should be brought + //! back into existence + /*! + * + * @param iph Phase id of the deleted phase + */ + double vcs_phaseStabilityTest(const int iph); //! Solve an equilibrium problem at a particular fixed temperature //! and pressure @@ -1441,7 +1476,7 @@ public: * stoichiometric coefficient of one is assumed for the * species K in this mechanism. * - * NOTE: K = IRXN + NC + * NOTE: kspec = Irxn + m_numComponents * * sc[irxn][j] : * j refers to the component number, and irxn diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 7b5c2a64e..36d727eb9 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -447,45 +447,74 @@ namespace VCSnonideal { vcs_dzero(VCS_DATA_PTR(m_deltaMolNumSpecies), m_numSpeciesTot); /* - * Figure out whether we will calculate new reaction step sizes - * for the major species. - * -> We won't if all species are minors (im), OR - * all major species have already converged + * First step is a major branch in the algorithm. + * We first determine if a phase pops into existence. */ - if (!(MajorSpeciesHaveConverged) && ! allMinorZeroedSpecies) { - soldel = vcs_RxnStepSizes(); - /* - If SOLDEL is true then we encountered a reaction between */ - /* - single-species-phase species, only, and have adjusted */ - /* - the mole number vector, W(), directly. In this case, */ - /* - we should immediately go back and recompute a new */ - /* - component basis, if the species that was zeroed was */ - /* - a component. SOLDEL is true when this is so. */ - if (soldel > 0) { - /* - We have changed the base mole number amongst single- */ - /* - species-phase species. However, we don't need to */ - /* - recaculate their chemical potentials because they */ - /* - are constant, anyway! */ - if (soldel == 2) { - goto L_COMPONENT_CALC; - } - /* - We have not changed the actual DG values for */ - /* - any species, even the one we deleted. Thus, */ - /* - we don't need to start over. */ - } - } else { + int iphasePop = vcs_popPhaseID(); + /* + * + */ + soldel = -1; + if (iphasePop >= 0) { + soldel = vcs_popPhaseRxnStepSizes(iphasePop); + if (soldel == 3) { + iphasePop = -1; #ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - if (allMinorZeroedSpecies) { - plogf(" --- vcs_RxnStepSizes not called because all" - "species are minors\n"); - } else { - plogf(" --- vcs_RxnStepSizes not called because " - "all majors have converged\n"); + if (m_debug_print_lvl >= 2) { + plogf(" --- vcs_popPhaseRxnStepSizes() was called but stoich prevented phase %d popping\n"); } - } #endif + } } - + if (iphasePop < 0) { + /* + * Figure out whether we will calculate new reaction step sizes + * for the major species. + * -> We won't if all species are minors (im), OR + * all major species have already converged + */ + if (!(MajorSpeciesHaveConverged) && ! allMinorZeroedSpecies) { + soldel = vcs_RxnStepSizes(); + /* - If SOLDEL is true then we encountered a reaction between */ + /* - single-species-phase species, only, and have adjusted */ + /* - the mole number vector, W(), directly. In this case, */ + /* - we should immediately go back and recompute a new */ + /* - component basis, if the species that was zeroed was */ + /* - a component. SOLDEL is true when this is so. */ + if (soldel > 0) { + /* - We have changed the base mole number amongst single- */ + /* - species-phase species. However, we don't need to */ + /* - recaculate their chemical potentials because they */ + /* - are constant, anyway! */ + if (soldel == 2) { + goto L_COMPONENT_CALC; + } + /* - We have not changed the actual DG values for */ + /* - any species, even the one we deleted. Thus, */ + /* - we don't need to start over. */ + } + } else { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + if (allMinorZeroedSpecies) { + plogf(" --- vcs_RxnStepSizes not called because all" + "species are minors\n"); + } else { + plogf(" --- vcs_RxnStepSizes not called because " + "all majors have converged\n"); + } + } +#endif + } + } + #ifdef DEBUG_MODE + else { + if (m_debug_print_lvl >= 2) { + plogf(" --- vcs_RxnStepSizes not called because alternative" + "phase creation delta was used instead\n"); + } + } +#endif lec = FALSE; doPhaseDeleteIph = -1; doPhaseDeleteKspec = -1; @@ -544,389 +573,403 @@ namespace VCSnonideal { #ifdef DEBUG_MODE ANOTE[0] = '\0'; #endif - - if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) { - /********************************************************************/ - /************************ VOLTAGE SPECIES ***************************/ - /********************************************************************/ -#ifdef DEBUG_MODE - dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE); -#else - dx = vcs_minor_alt_calc(kspec, irxn, &soldel); -#endif - m_deltaMolNumSpecies[kspec] = dx; - } - else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) { - /********************************************************************/ - /********************** ZEROED OUT SPECIES **************************/ - /********************************************************************/ - bool resurrect = true; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 3) { - plogf(" --- %s currently zeroed (SpStatus=%-2d):", - m_speciesName[kspec].c_str(), m_speciesStatus[kspec]); - plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n", - irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec], - m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]); - } -#endif - // HKM Alternative is to not allow ds[] = 0.0 phases - // to pop back into existence. For esthetics, I'm allowing this. - // so that dg < 0.0 phases with zero mole numbers become components. - // This is also better, because that component will be the first - // one to pop into existence if there is a minute quantity of the element. - // This could change in the future. - //if (dg[irxn] >= 0.0 || ds[kspec] <= 0.0) { - if (m_deltaGRxn_new[irxn] >= 0.0 ) { - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; - m_deltaMolNumSpecies[kspec] = 0.0; - resurrect = false; -#ifdef DEBUG_MODE - sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", - m_deltaGRxn_new[irxn]); - if (m_deltaGRxn_new[irxn] < 0.0) { - sprintf(ANOTE, "Species stays zeroed even though dg neg:DG = %11.4E, ds zeroed ", - m_deltaGRxn_new[irxn]); - } -#endif + if (iphasePop >= 0) { + if (iph == iphasePop) { + dx = m_deltaMolNumSpecies[kspec]; + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; +#ifdef DEBUG_MODE + sprintf(ANOTE, "Phase pop"); +#endif } else { - for (int j = 0; j < m_numElemConstraints; ++j) { - int elType = m_elType[j]; - if (elType == VCS_ELEM_TYPE_ABSPOS) { - double atomComp = m_formulaMatrix[j][kspec]; - if (atomComp > 0.0) { - double maxPermissible = m_elemAbundancesGoal[j] / atomComp; - if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { -#ifdef DEBUG_MODE - sprintf(ANOTE, "Species stays zeroed even though dG " - "neg, because of %s elemAbund", - m_elementName[j].c_str()); -#endif - resurrect = false; - break; - } - } - } - } - } - /* - * Resurrect the species - */ - if (resurrect) { - bool phaseResurrected = false; - if (Vphase->exists() == VCS_PHASE_EXIST_NO) { - //Vphase->setExistence(1); - phaseResurrected = true; - } - --m_numRxnMinorZeroed; - - if (phaseResurrected) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf(" --- Zeroed species changed to major: "); - plogf("%-12s\n", m_speciesName[kspec].c_str()); - } -#endif - m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; - MajorSpeciesHaveConverged = false; - allMinorZeroedSpecies = false; - } else { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf(" --- Zeroed species changed to minor: "); - plogf("%-12s\n", m_speciesName[kspec].c_str()); - } -#endif - m_speciesStatus[kspec] = VCS_SPECIES_MINOR; - } - if (m_deltaMolNumSpecies[kspec] > 0.0) { - dx = m_deltaMolNumSpecies[kspec] * 0.01; - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; - } else { - m_molNumSpecies_new[kspec] = m_totalMolNum * VCS_DELETE_PHASE_CUTOFF * 10.; - dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec]; - } - m_deltaMolNumSpecies[kspec] = dx; -#ifdef DEBUG_MODE - sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]); -#endif - } else { - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; - m_deltaMolNumSpecies[kspec] = 0.0; dx = 0.0; - } - } else if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) { - /********************************************************************/ - /***************************** MINOR SPECIES ************************/ - /********************************************************************/ - /* - * Unless ITI isn't equal to zero we zero out changes - * to minor species. - */ - if (iti != 0) { m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; - m_deltaMolNumSpecies[kspec] = 0.0; - dx = 0.0; -#ifdef DEBUG_MODE - sprintf(ANOTE,"minor species not considered"); - if (m_debug_print_lvl >= 2) { - plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str()); - plogf("%3d%11.4E%11.4E%11.4E | %s", - m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); - plogendl(); - } -#endif - continue; - } - /* - * Minor species alternative calculation - * --------------------------------------- - * This is based upon the following approximation: - * The mole fraction changes due to these reactions don't affect - * the mole numbers of the component species. Therefore the - * following approximation is valid for an ideal solution - * 0 = DG(I) + log(WT(I)/W(I)) - * (DG contains the contribution from FF(I) + log(W(I)/TL) ) - * Thus, - * WT(I) = W(I) EXP(-DG(I)) - * If soldel is true on return, then we branch to the section - * that deletes a species from the current set of active species. - */ -#ifdef DEBUG_MODE - dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE); -#else - dx = vcs_minor_alt_calc(kspec, irxn, &soldel); -#endif - m_deltaMolNumSpecies[kspec] = dx; - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; - - if (soldel) { - /*******************************************************************/ - /***** DELETE MINOR SPECIES LESS THAN VCS_DELETE_SPECIES_CUTOFF */ - /***** MOLE NUMBER */ - /*******************************************************************/ -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf(" --- Delete minor species in multispec phase: %-12s", - m_speciesName[kspec].c_str()); - plogendl(); - } -#endif - m_deltaMolNumSpecies[kspec] = 0.0; - /* - * Delete species, kspec. The alternate return is for the case - * where all species become deleted. Then, we need to - * branch to the code where we reevaluate the deletion - * of all species. - */ - lnospec = vcs_delete_species(kspec); - if (lnospec) goto L_RECHECK_DELETED; - /* - * Go back to consider the next species in the list. - * Note, however, that the next species in the list is now - * in slot l. In deleting the previous species L, We have - * exchanged slot MR with slot l, and then have - * decremented MR. - * Therefore, we will decrement the species counter, here. - */ - --irxn; -#ifdef DEBUG_MODE - goto L_MAIN_LOOP_END_NO_PRINT; -#else - goto L_MAIN_LOOP_END; -#endif } } else { - /********************************************************************/ - /*********************** MAJOR SPECIES ******************************/ - /********************************************************************/ -#ifdef DEBUG_MODE - sprintf(ANOTE, "Normal Major Calc"); + + + if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) { + /********************************************************************/ + /************************ VOLTAGE SPECIES ***************************/ + /********************************************************************/ +#ifdef DEBUG_MODE + dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE); +#else + dx = vcs_minor_alt_calc(kspec, irxn, &soldel); #endif - /* - * Check for superconvergence of the formation reaction. Do - * nothing if it is superconverged. Skip to the end of the - * irxn loop if it is superconverged. - */ - if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) { - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; - m_deltaMolNumSpecies[kspec] = 0.0; - dx = 0.0; + m_deltaMolNumSpecies[kspec] = dx; + } + else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) { + /********************************************************************/ + /********************** ZEROED OUT SPECIES **************************/ + /********************************************************************/ + bool resurrect = true; #ifdef DEBUG_MODE - sprintf(ANOTE, "major species is converged"); - if (m_debug_print_lvl >= 2) { - plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str()); - plogf("%3d%11.4E%11.4E%11.4E | %s", - m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], - m_deltaMolNumSpecies[kspec], ANOTE); - plogendl(); + if (m_debug_print_lvl >= 3) { + plogf(" --- %s currently zeroed (SpStatus=%-2d):", + m_speciesName[kspec].c_str(), m_speciesStatus[kspec]); + plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n", + irxn, m_deltaGRxn_new[irxn], m_molNumSpecies_new[kspec], + m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec]); } #endif - continue; - } - /* - * Set the initial step size, dx, equal to the value produced - * by the routine, vcs_RxnStepSize(). - * - * Note the multiplition logic is to make sure that - * dg[] didn't change sign due to w[] changing in the - * middle of the iteration. (it can if a single species - * phase goes out of existence). - */ - if ((m_deltaGRxn_new[irxn] * m_deltaMolNumSpecies[kspec]) <= 0.0) { - dx = m_deltaMolNumSpecies[kspec]; - } else { - dx = 0.0; - m_deltaMolNumSpecies[kspec] = 0.0; + // HKM Alternative is to not allow ds[] = 0.0 phases + // to pop back into existence. For esthetics, I'm allowing this. + // so that dg < 0.0 phases with zero mole numbers become components. + // This is also better, because that component will be the first + // one to pop into existence if there is a minute quantity of the element. + // This could change in the future. + //if (dg[irxn] >= 0.0 || ds[kspec] <= 0.0) { + if (m_deltaGRxn_new[irxn] >= 0.0 ) { + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; + m_deltaMolNumSpecies[kspec] = 0.0; + resurrect = false; #ifdef DEBUG_MODE - sprintf(ANOTE, "dx set to 0, DG flipped sign due to " - "changed initial point"); + sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", + m_deltaGRxn_new[irxn]); + if (m_deltaGRxn_new[irxn] < 0.0) { + sprintf(ANOTE, "Species stays zeroed even though dg neg:DG = %11.4E, ds zeroed ", + m_deltaGRxn_new[irxn]); + } #endif - } - /* - * Form a tentative value of the new species moles - */ - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; - - /* - * Check for non-positive mole fraction of major species. - * If we find one, we branch to a section below. Then, - * depending upon the outcome, we branch to sections below, - * or we restart the entire iteration. - */ - if (m_molNumSpecies_new[kspec] <= 0.0) { -#ifdef DEBUG_MODE - sprintf(ANOTE, "initial nonpos kmoles= %11.3E", - m_molNumSpecies_new[kspec]); -#endif - /* ************************************************* */ - /* *** NON-POSITIVE MOLES OF MAJOR SPECIES ********* */ - /* ************************************************* */ - /* - * We are here when a tentative value of a mole fraction - * created by a tentative value of M_DELTAMOLNUMSPECIES(*) is negative. - * We branch from here depending upon whether this - * species is in a single species phase or in - * a multispecies phase. - */ - if (! (m_SSPhase[kspec])) { - /* - * Section for multispecies phases: - * - Cut reaction adjustment for positive kmoles of - * major species in multispecies phases. - * Decrease its concentration by a factor of 10. - */ - dx = -0.9 * m_molNumSpecies_old[kspec]; - m_deltaMolNumSpecies[kspec] = dx; - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; } else { - /* - * Section for single species phases: - * Calculate a dx that will wipe out the - * moles in the phase. - */ - dx = -m_molNumSpecies_old[kspec]; - /* - * Calculate an update that doesn't create a negative mole - * number for a component species. Actually, restrict this - * a little more so that the component values can only be - * reduced by two 99%, - */ - for (j = 0; j < m_numComponents; ++j) { - if (sc_irxn[j] != 0.0) { - wx[j] = m_molNumSpecies_old[j] + sc_irxn[j] * dx; - if (wx[j] <= m_molNumSpecies_old[j] * 0.01 - 1.0E-150) { - dx = MAX(dx, m_molNumSpecies_old[j] * -0.99 / sc_irxn[j]); + for (int j = 0; j < m_numElemConstraints; ++j) { + int elType = m_elType[j]; + if (elType == VCS_ELEM_TYPE_ABSPOS) { + double atomComp = m_formulaMatrix[j][kspec]; + if (atomComp > 0.0) { + double maxPermissible = m_elemAbundancesGoal[j] / atomComp; + if (maxPermissible < VCS_DELETE_MINORSPECIES_CUTOFF) { +#ifdef DEBUG_MODE + sprintf(ANOTE, "Species stays zeroed even though dG " + "neg, because of %s elemAbund", + m_elementName[j].c_str()); +#endif + resurrect = false; + break; + } } - } else { - wx[j] = m_molNumSpecies_old[j]; } } - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; - if (m_molNumSpecies_new[kspec] > 0.0) { - m_deltaMolNumSpecies[kspec] = dx; -#ifdef DEBUG_MODE - sprintf(ANOTE, - "zeroing SS phase created a neg component species " - "-> reducing step size instead"); -#endif - } else { - /* - * We are going to zero the single species phase. - * Set the existence flag - */ - iph = m_phaseID[kspec]; - Vphase = m_VolPhaseList[iph]; - //Vphase->setExistence(0); -#ifdef DEBUG_MODE - sprintf(ANOTE, "zeroing out SS phase: "); -#endif - /* - * Change the base mole numbers for the iteration. - * We need to do this here, because we have decided - * to eliminate the phase in this special section - * outside the main loop. - */ - m_molNumSpecies_new[kspec] = 0.0; - doPhaseDeleteIph = iph; - doPhaseDeleteKspec = kspec; + } + /* + * Resurrect the species + */ + if (resurrect) { + bool phaseResurrected = false; + if (Vphase->exists() == VCS_PHASE_EXIST_NO) { + //Vphase->setExistence(1); + phaseResurrected = true; + } + --m_numRxnMinorZeroed; + if (phaseResurrected) { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - if (m_speciesStatus[kspec] >= 0) { - plogf(" --- SS species changed to zeroedss: "); - plogf("%-12s", m_speciesName[kspec].c_str()); - plogendl(); - } + plogf(" --- Zeroed species changed to major: "); + plogf("%-12s\n", m_speciesName[kspec].c_str()); } #endif - m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; - ++m_numRxnMinorZeroed; - allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc); - - for (int kk = 0; kk < m_numSpeciesTot; kk++) { - m_deltaMolNumSpecies[kk] = 0.0; - m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk]; + m_speciesStatus[kspec] = VCS_SPECIES_MAJOR; + MajorSpeciesHaveConverged = false; + allMinorZeroedSpecies = false; + } else { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- Zeroed species changed to minor: "); + plogf("%-12s\n", m_speciesName[kspec].c_str()); } - m_deltaMolNumSpecies[kspec] = dx; - m_molNumSpecies_new[kspec] = 0.0; - - for (k = 0; k < m_numComponents; ++k) { - m_deltaMolNumSpecies[k] = 0.0; - } - for (iph = 0; iph < m_numPhases; iph++) { - m_deltaPhaseMoles[iph] = 0.0; - } - +#endif + m_speciesStatus[kspec] = VCS_SPECIES_MINOR; } + if (m_deltaMolNumSpecies[kspec] > 0.0) { + dx = m_deltaMolNumSpecies[kspec] * 0.01; + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; + } else { + m_molNumSpecies_new[kspec] = m_totalMolNum * VCS_DELETE_PHASE_CUTOFF * 10.; + dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec]; + } + m_deltaMolNumSpecies[kspec] = dx; +#ifdef DEBUG_MODE + sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", m_deltaGRxn_new[irxn]); +#endif + } else { + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; + m_deltaMolNumSpecies[kspec] = 0.0; + dx = 0.0; } + } else if (m_speciesStatus[kspec] == VCS_SPECIES_MINOR) { + /********************************************************************/ + /***************************** MINOR SPECIES ************************/ + /********************************************************************/ + /* + * Unless ITI isn't equal to zero we zero out changes + * to minor species. + */ + if (iti != 0) { + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; + m_deltaMolNumSpecies[kspec] = 0.0; + dx = 0.0; +#ifdef DEBUG_MODE + sprintf(ANOTE,"minor species not considered"); + if (m_debug_print_lvl >= 2) { + plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str()); + plogf("%3d%11.4E%11.4E%11.4E | %s", + m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); + plogendl(); + } +#endif + continue; + } + /* + * Minor species alternative calculation + * --------------------------------------- + * This is based upon the following approximation: + * The mole fraction changes due to these reactions don't affect + * the mole numbers of the component species. Therefore the + * following approximation is valid for an ideal solution + * 0 = DG(I) + log(WT(I)/W(I)) + * (DG contains the contribution from FF(I) + log(W(I)/TL) ) + * Thus, + * WT(I) = W(I) EXP(-DG(I)) + * If soldel is true on return, then we branch to the section + * that deletes a species from the current set of active species. + */ +#ifdef DEBUG_MODE + dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE); +#else + dx = vcs_minor_alt_calc(kspec, irxn, &soldel); +#endif + m_deltaMolNumSpecies[kspec] = dx; + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; + + if (soldel) { + /*******************************************************************/ + /***** DELETE MINOR SPECIES LESS THAN VCS_DELETE_SPECIES_CUTOFF */ + /***** MOLE NUMBER */ + /*******************************************************************/ +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- Delete minor species in multispec phase: %-12s", + m_speciesName[kspec].c_str()); + plogendl(); + } +#endif + m_deltaMolNumSpecies[kspec] = 0.0; + /* + * Delete species, kspec. The alternate return is for the case + * where all species become deleted. Then, we need to + * branch to the code where we reevaluate the deletion + * of all species. + */ + lnospec = vcs_delete_species(kspec); + if (lnospec) goto L_RECHECK_DELETED; + /* + * Go back to consider the next species in the list. + * Note, however, that the next species in the list is now + * in slot l. In deleting the previous species L, We have + * exchanged slot MR with slot l, and then have + * decremented MR. + * Therefore, we will decrement the species counter, here. + */ + --irxn; +#ifdef DEBUG_MODE + goto L_MAIN_LOOP_END_NO_PRINT; +#else + goto L_MAIN_LOOP_END; +#endif + } + } else { + /********************************************************************/ + /*********************** MAJOR SPECIES ******************************/ + /********************************************************************/ +#ifdef DEBUG_MODE + sprintf(ANOTE, "Normal Major Calc"); +#endif + /* + * Check for superconvergence of the formation reaction. Do + * nothing if it is superconverged. Skip to the end of the + * irxn loop if it is superconverged. + */ + if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) { + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; + m_deltaMolNumSpecies[kspec] = 0.0; + dx = 0.0; +#ifdef DEBUG_MODE + sprintf(ANOTE, "major species is converged"); + if (m_debug_print_lvl >= 2) { + plogf(" --- "); plogf("%-12s", m_speciesName[kspec].c_str()); + plogf("%3d%11.4E%11.4E%11.4E | %s", + m_speciesStatus[kspec], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], + m_deltaMolNumSpecies[kspec], ANOTE); + plogendl(); + } +#endif + continue; + } + /* + * Set the initial step size, dx, equal to the value produced + * by the routine, vcs_RxnStepSize(). + * + * Note the multiplition logic is to make sure that + * dg[] didn't change sign due to w[] changing in the + * middle of the iteration. (it can if a single species + * phase goes out of existence). + */ + if ((m_deltaGRxn_new[irxn] * m_deltaMolNumSpecies[kspec]) <= 0.0) { + dx = m_deltaMolNumSpecies[kspec]; + } else { + dx = 0.0; + m_deltaMolNumSpecies[kspec] = 0.0; +#ifdef DEBUG_MODE + sprintf(ANOTE, "dx set to 0, DG flipped sign due to " + "changed initial point"); +#endif + } + /* + * Form a tentative value of the new species moles + */ + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; + + /* + * Check for non-positive mole fraction of major species. + * If we find one, we branch to a section below. Then, + * depending upon the outcome, we branch to sections below, + * or we restart the entire iteration. + */ + if (m_molNumSpecies_new[kspec] <= 0.0) { +#ifdef DEBUG_MODE + sprintf(ANOTE, "initial nonpos kmoles= %11.3E", + m_molNumSpecies_new[kspec]); +#endif + /* ************************************************* */ + /* *** NON-POSITIVE MOLES OF MAJOR SPECIES ********* */ + /* ************************************************* */ + /* + * We are here when a tentative value of a mole fraction + * created by a tentative value of M_DELTAMOLNUMSPECIES(*) is negative. + * We branch from here depending upon whether this + * species is in a single species phase or in + * a multispecies phase. + */ + if (! (m_SSPhase[kspec])) { + /* + * Section for multispecies phases: + * - Cut reaction adjustment for positive kmoles of + * major species in multispecies phases. + * Decrease its concentration by a factor of 10. + */ + dx = -0.9 * m_molNumSpecies_old[kspec]; + m_deltaMolNumSpecies[kspec] = dx; + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; + } else { + /* + * Section for single species phases: + * Calculate a dx that will wipe out the + * moles in the phase. + */ + dx = -m_molNumSpecies_old[kspec]; + /* + * Calculate an update that doesn't create a negative mole + * number for a component species. Actually, restrict this + * a little more so that the component values can only be + * reduced by two 99%, + */ + for (j = 0; j < m_numComponents; ++j) { + if (sc_irxn[j] != 0.0) { + wx[j] = m_molNumSpecies_old[j] + sc_irxn[j] * dx; + if (wx[j] <= m_molNumSpecies_old[j] * 0.01 - 1.0E-150) { + dx = MAX(dx, m_molNumSpecies_old[j] * -0.99 / sc_irxn[j]); + } + } else { + wx[j] = m_molNumSpecies_old[j]; + } + } + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; + if (m_molNumSpecies_new[kspec] > 0.0) { + m_deltaMolNumSpecies[kspec] = dx; +#ifdef DEBUG_MODE + sprintf(ANOTE, + "zeroing SS phase created a neg component species " + "-> reducing step size instead"); +#endif + } else { + /* + * We are going to zero the single species phase. + * Set the existence flag + */ + iph = m_phaseID[kspec]; + Vphase = m_VolPhaseList[iph]; + //Vphase->setExistence(0); +#ifdef DEBUG_MODE + sprintf(ANOTE, "zeroing out SS phase: "); +#endif + /* + * Change the base mole numbers for the iteration. + * We need to do this here, because we have decided + * to eliminate the phase in this special section + * outside the main loop. + */ + m_molNumSpecies_new[kspec] = 0.0; + doPhaseDeleteIph = iph; + doPhaseDeleteKspec = kspec; + +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + if (m_speciesStatus[kspec] >= 0) { + plogf(" --- SS species changed to zeroedss: "); + plogf("%-12s", m_speciesName[kspec].c_str()); + plogendl(); + } + } +#endif + m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; + ++m_numRxnMinorZeroed; + allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc); + + for (int kk = 0; kk < m_numSpeciesTot; kk++) { + m_deltaMolNumSpecies[kk] = 0.0; + m_molNumSpecies_new[kk] = m_molNumSpecies_old[kk]; + } + m_deltaMolNumSpecies[kspec] = dx; + m_molNumSpecies_new[kspec] = 0.0; + + for (k = 0; k < m_numComponents; ++k) { + m_deltaMolNumSpecies[k] = 0.0; + } + for (iph = 0; iph < m_numPhases; iph++) { + m_deltaPhaseMoles[iph] = 0.0; + } + + } + } + } #ifdef VCS_LINE_SEARCH - /*********************************************************************/ - /*** LINE SEARCH ALGORITHM FOR MAJOR SPECIES IN NON-IDEAL PHASES *****/ - /*********************************************************************/ - /* - * Skip the line search if we are birthing a species - */ - if ((dx != 0.0) && - (m_molNumSpecies_old[kspec] > 0.0) && - (doPhaseDeleteIph == -1) && - (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { - double dx_old = dx; + /*********************************************************************/ + /*** LINE SEARCH ALGORITHM FOR MAJOR SPECIES IN NON-IDEAL PHASES *****/ + /*********************************************************************/ + /* + * Skip the line search if we are birthing a species + */ + if ((dx != 0.0) && + (m_molNumSpecies_old[kspec] > 0.0) && + (doPhaseDeleteIph == -1) && + (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { + double dx_old = dx; #ifdef DEBUG_MODE - dx = vcs_line_search(irxn, dx_old, ANOTE); + dx = vcs_line_search(irxn, dx_old, ANOTE); #else - dx = vcs_line_search(irxn, dx_old); + dx = vcs_line_search(irxn, dx_old); #endif - vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); - } - m_deltaMolNumSpecies[kspec] = dx; + vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); + } + m_deltaMolNumSpecies[kspec] = dx; #endif - - } /* End of Loop on ic[irxn] -> the type of species */ + }/* End of Loop on ic[irxn] -> the type of species */ + } /***********************************************************************/ /****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/ /***********************************************************************/ @@ -4748,7 +4791,8 @@ namespace VCSnonideal { * This should be implemented. */ int k; - if (alterZeroedPhases) { + //alterZeroedPhases = false; + if (alterZeroedPhases && false) { for (iph = 0; iph < m_numPhases; iph++) { lneed = FALSE; vcs_VolPhase *Vphase = m_VolPhaseList[iph]; diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp index a8b797fed..0432ff2e1 100644 --- a/Cantera/src/thermo/HMWSoln.cpp +++ b/Cantera/src/thermo/HMWSoln.cpp @@ -75,8 +75,10 @@ namespace Cantera { MC_apCut_(0.0), MC_bpCut_(0.0), MC_cpCut_(0.0), - CROP_ln_gamma_o_min(-25.0), - CROP_ln_gamma_k_max(23.0), + CROP_ln_gamma_o_min(-10.0), + CROP_ln_gamma_o_max(3.0), + CROP_ln_gamma_k_min(-5.0), + CROP_ln_gamma_k_max(15.0), m_debugCalc(0) { for (int i = 0; i < 17; i++) { @@ -130,8 +132,10 @@ namespace Cantera { MC_apCut_(0.0), MC_bpCut_(0.0), MC_cpCut_(0.0), - CROP_ln_gamma_o_min(-25.0), - CROP_ln_gamma_k_max(23.0), + CROP_ln_gamma_o_min(-10.0), + CROP_ln_gamma_o_max(3.0), + CROP_ln_gamma_k_min(-5.0), + CROP_ln_gamma_k_max(15.0), m_debugCalc(0) { for (int i = 0; i < 17; i++) { @@ -179,8 +183,10 @@ namespace Cantera { MC_apCut_(0.0), MC_bpCut_(0.0), MC_cpCut_(0.0), - CROP_ln_gamma_o_min(-25.0), - CROP_ln_gamma_k_max(23.0), + CROP_ln_gamma_o_min(-10.0), + CROP_ln_gamma_o_max(3.0), + CROP_ln_gamma_k_min(-5.0), + CROP_ln_gamma_k_max(15.0), m_debugCalc(0) { for (int i = 0; i < 17; i++) { @@ -234,8 +240,10 @@ namespace Cantera { MC_apCut_(0.0), MC_bpCut_(0.0), MC_cpCut_(0.0), - CROP_ln_gamma_o_min(-25.0), - CROP_ln_gamma_k_max(23.0), + CROP_ln_gamma_o_min(-10.0), + CROP_ln_gamma_o_max(3.0), + CROP_ln_gamma_k_min(-5.0), + CROP_ln_gamma_k_max(15.0), m_debugCalc(0) { /* @@ -389,6 +397,8 @@ namespace Cantera { MC_bpCut_ = b.MC_bpCut_; MC_cpCut_ = b.MC_cpCut_; CROP_ln_gamma_o_min = b.CROP_ln_gamma_o_min; + CROP_ln_gamma_o_max = b.CROP_ln_gamma_o_max; + CROP_ln_gamma_k_min = b.CROP_ln_gamma_k_min; CROP_ln_gamma_k_max = b.CROP_ln_gamma_k_max; m_CounterIJ = b.m_CounterIJ; m_molalitiesCropped = b.m_molalitiesCropped; @@ -463,8 +473,10 @@ namespace Cantera { MC_apCut_(0.0), MC_bpCut_(0.0), MC_cpCut_(0.0), - CROP_ln_gamma_o_min(-25.0), - CROP_ln_gamma_k_max(23.0), + CROP_ln_gamma_o_min(-10.0), + CROP_ln_gamma_o_max(3.0), + CROP_ln_gamma_k_min(-5.0), + CROP_ln_gamma_k_max(15.0), m_debugCalc(0) { if (testProb != 1) { @@ -1812,14 +1824,22 @@ namespace Cantera { for (int k = 1; k < m_kk; k++) { m_lnActCoeffMolal_Unscaled[k] += IMS_lnActCoeffMolal_[k]; - if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max + lnxs)) { - m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max + lnxs; + if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max)) { + m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max; + } + if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) { + // -1.0 and -1.5 caused multiple solutions + m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs; } } m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0); - if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) { - m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min - lnxs; + if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) { + m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min; + } + if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) { + // -0.5 caused multiple solutions + m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max; } /* diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h index 97d2015ef..13818fa9c 100644 --- a/Cantera/src/thermo/HMWSoln.h +++ b/Cantera/src/thermo/HMWSoln.h @@ -3194,9 +3194,12 @@ namespace Cantera { doublereal MC_cpCut_; doublereal CROP_ln_gamma_o_min; + doublereal CROP_ln_gamma_o_max; + doublereal CROP_ln_gamma_k_min; doublereal CROP_ln_gamma_k_max; + //! Local error routine /*! * @param msg print out a message and error exit