From d2358f18182fb987df10871d3fc93751b7085949 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Mon, 18 Aug 2008 23:17:18 +0000 Subject: [PATCH] On the road to a better scheme for handling the deletion of multispecies phases. --- Cantera/src/equil/vcs_VolPhase.cpp | 6 +- Cantera/src/equil/vcs_elem.cpp | 78 +++++----- Cantera/src/equil/vcs_rxnadj.cpp | 6 +- Cantera/src/equil/vcs_solve.cpp | 20 ++- Cantera/src/equil/vcs_solve.h | 7 +- Cantera/src/equil/vcs_solve_TP.cpp | 225 +++++++++++++++++++---------- 6 files changed, 214 insertions(+), 128 deletions(-) diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index 90b97781c..2caa8d08c 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -635,9 +635,9 @@ namespace VCSnonideal { // This is where we will start to store a better approximation // for the mole fractions, when the phase doesn't exist. // This is currently unimplemented. - for (int k = 0; k < m_numSpecies; k++) { - Xmol[k] = 1.0 / m_numSpecies; - } + //for (int k = 0; k < m_numSpecies; k++) { + // Xmol[k] = 1.0 / m_numSpecies; + //} m_existence = VCS_PHASE_EXIST_NO; } /* diff --git a/Cantera/src/equil/vcs_elem.cpp b/Cantera/src/equil/vcs_elem.cpp index be6bee774..5e2118f8a 100644 --- a/Cantera/src/equil/vcs_elem.cpp +++ b/Cantera/src/equil/vcs_elem.cpp @@ -74,52 +74,54 @@ namespace VCSnonideal { * 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)) { - /* - * 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; + if (m_elementActive[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 = m_formulaMatrix[i][kspec]; - if (eval < 0.0) { - multisign = true; - } - if (eval != 0.0) { - scale = MAX(scale, fabs(eval * m_molNumSpecies_old[kspec])); - numNonZero++; - } + if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) { + AssertThrowVCS(m_elemAbundancesGoal[i] == 0.0, "vcs_elabcheck"); } - if (multisign) { - if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) { - return FALSE; + 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 = m_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 { - if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) { + /* + * 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; } } - } 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; - } } } } diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index e439f18d7..7772d82dd 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -303,8 +303,8 @@ namespace VCSnonideal { if (k == kspec) { m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDSS; if (m_SSPhase[kspec] != 1) { - printf("we shouldn't be here!\n"); - exit(-1); + printf("vcs_RxnStepSizes:: we shouldn't be here!\n"); + std::exit(-1); } } #ifdef DEBUG_MODE @@ -586,7 +586,7 @@ namespace VCSnonideal { double diag = hessianDiag_Ideal; double hessActCoef = vcs_Hessian_actCoeff_diag(irxn); if (hessianDiag_Ideal <= 0.0) { - plogf("We shouldn't be here\n"); + plogf("vcs_Hessian_diag_adj::We shouldn't be here\n"); exit(-1); } if (hessActCoef >= 0.0) { diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 0eb3436f0..d9950ee48 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -62,6 +62,7 @@ namespace VCSnonideal { } // Initialize the sizes within the VCS_SOLVE object + /* * This resizes all of the internal arrays within the object. This routine * operates in two modes. If all of the parameters are the same as it @@ -691,7 +692,7 @@ namespace VCSnonideal { m_elementActive.resize(nelements, 1); /* - * Copy over the element names + * Copy over the element names and types */ for (i = 0; i < nelements; i++) { m_elementName[i] = pub->ElName[i]; @@ -705,6 +706,23 @@ namespace VCSnonideal { } } } + + for (i = 0; i < nelements; i++) { + if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) { + if (m_elemAbundancesGoal[i] != 0.0) { + if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) { + plogf("Charge neutrality condition %s is signicantly nonzero, %g. Giving up\n", + m_elementName[i].c_str(), m_elemAbundancesGoal[i]); + std::exit(-1); + } else { + plogf("Charge neutrality condition %s not zero, %g. Setting it zero\n", + m_elementName[i].c_str(), m_elemAbundancesGoal[i]); + m_elemAbundancesGoal[i] = 0.0; + } + + } + } + } /* * Copy over the species names diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 520575a69..0e42ca58c 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -507,7 +507,9 @@ public: * Reconciles Phase existence flags with total moles in each phase. */ double vcs_tmoles(); - +#ifdef DEBUG_MODE + void check_tmoles() const; +#endif //! This subroutine calculates reaction free energy changes for //! all noncomponent formation reactions. @@ -1415,7 +1417,6 @@ public: * * length = [nspecies0][nelements0] */ - //DoubleStarStar sc; DoubleStarStar m_stoichCoeffRxnMatrix; //! Absolute size of the stoichiometric coefficients @@ -1718,7 +1719,7 @@ public: //! Species string name for the kth species /*! - * SpName[k] = Species string name for the kth species + * Species string name for the kth species */ std::vector m_speciesName; diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index a15f385ac..f89aaa86b 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -404,7 +404,7 @@ namespace VCSnonideal { if (uptodate_minors == FALSE) { vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); - vcs_deltag(1, false, VCS_STATECALC_NEW); + vcs_deltag(1, false, VCS_STATECALC_OLD); } uptodate_minors = TRUE; } else { @@ -422,7 +422,15 @@ namespace VCSnonideal { plogf(" (only major species)\n"); } } - + /* + * Calculate the total moles in each phase -> old solution + * -> Needed for numerical stability when phases disappear. + * -> the phase moles tend to drift off without this step. + */ +#ifdef DEBUG_MODE + check_tmoles(); +#endif + vcs_tmoles(); /* * Copy the old solution into the new solution as an initial guess */ @@ -431,6 +439,7 @@ namespace VCSnonideal { vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_new), VCS_DATA_PTR(m_actCoeffSpecies_old), m_numSpeciesRdc); vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_new), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); + /* Go find a new reaction adjustment -> * i.e., change in extent of reaction for each reaction. @@ -1278,25 +1287,32 @@ namespace VCSnonideal { if (m_tPhaseMoles_old[iph] != 0.0 && m_tPhaseMoles_old[iph]/m_totalMolNum <= VCS_DELETE_PHASE_CUTOFF) { soldel = 1; + for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) { if (m_phaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) { irxn = kspec - m_numComponents; - if (kspec < m_numComponents) { - if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF) { - soldel = 0; - break; - } - } else { - for (k = 0; k < m_numComponents; k++) { - if (m_stoichCoeffRxnMatrix[irxn][k] != 0.0) { - if (m_molNumSpecies_old[kspec]/m_molNumSpecies_old[k] > - VCS_DELETE_PHASE_CUTOFF) { - soldel = 0; - break; - } - } - } - } + + // Both of these conditions are false and should be discarded. + // I think a proper special case would be if a species in a small phase + // had a significant contribution to total element total of a single + // element constraint. That's it. + + // if (kspec < m_numComponents) { + // if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF) { + //soldel = 0; + //break; + //} + //} else { + //for (k = 0; k < m_numComponents; k++) { + //if (m_stoichCoeffRxnMatrix[irxn][k] != 0.0) { + // if (m_molNumSpecies_old[kspec]/m_molNumSpecies_old[k] > + //VCS_DELETE_PHASE_CUTOFF) { + // soldel = 0; + // break; + // } + //} + //} + //} } } if (soldel) { @@ -1323,23 +1339,19 @@ namespace VCSnonideal { retn = vcs_basopt(FALSE, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test, &usedZeroedSpecies); - if (retn != VCS_SUCCESS) return retn; + if (retn != VCS_SUCCESS) { +#ifdef DEBUG_MODE + plogf(" --- BASOPT returned with an error condition\n"); +#endif + std::exit(-1); + } vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); + vcs_deltag(0, true, VCS_STATECALC_OLD); - uptodate_minors = TRUE; - if (conv) { - /* - * HKM -> I don't understand why the code would just give - * up here in some cases. - * This should probably be taken out - */ - plogf(" DELETION OF MULTISPECIES PHASE. "); - plogf("Convergence to number of positive n(i) less than C.\n"); - plogf("Check results to follow carefully. \n"); - plogendl(); - goto L_RETURN_BLOCK; - } + iti = 0; + goto L_MAINLOOP_ALL_SPECIES ; + } /*************************************************************************/ /***************** CHECK FOR ELEMENT ABUNDANCE****************************/ @@ -2418,7 +2430,7 @@ namespace VCSnonideal { /* * When they are deleted, all of their species become active * species, even though their mole numbers are set to zero. - * The routine does not make the decision to eliminate multiphases. + * This routine does not make the decision to eliminate multiphases. * * Note, species in phases with zero mole numbers are still * considered active. Whether the phase pops back into @@ -2441,23 +2453,17 @@ namespace VCSnonideal { plogf(" --- delete_multiphase %d, %s\n", iph, Vphase->PhaseName.c_str()); } #endif - /* - * Zero out the total moles counters for the phase - */ - m_tPhaseMoles_old[iph] = 0.0; - m_tPhaseMoles_new[iph] = 0.0; - m_deltaPhaseMoles[iph] = 0.0; + /* - * Loop over all of the active species in the phase. + * Loop over all of the active noncomponent species in the phase. */ - for (kspec = 0; kspec < m_numSpeciesRdc; ++kspec) { + for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) { if (m_phaseID[kspec] == iph) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { irxn = kspec - m_numComponents; - if (irxn >= 0) { /* * calculate an extent of rxn, dx, that zeroes out the species. */ @@ -2493,40 +2499,66 @@ namespace VCSnonideal { * zeroed phase */ m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDMS; - /* - * Changed the component mole numbers to account for the - * final extent of reaction. Make sure to keep component - * mole numbers above zero - * - */ - // for (j = 0; j < m_numComponents; ++j) { - // m_molNumSpecies_old[j] += m_stoichCoeffRxnMatrix[irxn][j] * dx; - //if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_MOLNUM) { - //if (m_molNumSpecies_old[j] < 0.0) { - // m_molNumSpecies_old[j] = 0.0; - //} - //} - //} - } - else { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf(" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n", - kspec, m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec]); - } -#endif - if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF * VCS_DELETE_PHASE_CUTOFF ) { - plogf(" --- delete_multiphase unknown situation error exit"); - plogendl(); - std::exit(-1); - } else { - m_molNumSpecies_old[kspec] = 0.0; - m_molNumSpecies_new[kspec] = 0.0; - } - } } } } + int jcomp; + double deltaLarge, dj, dxWant, dxPerm, dxPerm2; + for (int kcomp = 0; kcomp < m_numComponents; ++kcomp) { + if (m_phaseID[kcomp] == iph) { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- delete_multiphase One of the species is a component %d - %s with mole number %g\n", + kcomp, m_speciesName[kcomp].c_str(), m_molNumSpecies_old[kcomp]); + } +#endif + if (m_molNumSpecies_old[kcomp] != 0.0) { + deltaLarge = 0.0; + for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) { + irxn = kspec - m_numComponents; + if (m_phaseID[kspec] != iph) { + if (m_stoichCoeffRxnMatrix[irxn][kcomp] != 0.0) { + dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix[irxn][kcomp]; + if (dxWant + m_molNumSpecies_old[kspec] < 0.0) { + dxPerm = -m_molNumSpecies_old[kspec]; + } + for (jcomp = 0; kcomp < m_numComponents; ++kcomp) { + if (jcomp != kcomp) { + if (m_phaseID[jcomp] == iph) { + dxPerm = 0.0; + } else { + dj = dxWant * m_stoichCoeffRxnMatrix[irxn][jcomp]; + if (dj + m_molNumSpecies_old[kcomp] < 0.0) { + dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix[irxn][jcomp]; + } + if (fabs(dxPerm2) < fabs(dxPerm)) { + dxPerm = dxPerm2; + } + } + } + } + } + if (dxPerm != 0.0) { + delta_species(kspec, &dxPerm); + } + } + } + + } + if (m_molNumSpecies_old[kcomp] != 0.0) { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- delete_multiphase One of the species is a component %d - %s still with mole number %g\n", + kcomp, m_speciesName[kcomp].c_str(), m_molNumSpecies_old[kcomp]); + plogf(" --- zeroing it \n"); + } +#endif + m_molNumSpecies_old[kcomp] = 0.0; + } + } + } + + /* * Loop over all of the inactive species in the phase: * Right now we reinstate all species in a deleted multiphase. @@ -2561,12 +2593,18 @@ namespace VCSnonideal { } } } + + /* + * Zero out the total moles counters for the phase + */ + m_tPhaseMoles_old[iph] = 0.0; + m_tPhaseMoles_new[iph] = 0.0; + m_deltaPhaseMoles[iph] = 0.0; /* * Upload the state to the VP object */ - Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, - VCS_DATA_PTR(m_molNumSpecies_old), - VCS_DATA_PTR(m_tPhaseMoles_old)); + Vphase->setTotalMoles(0.0); + return successful; } /**********************************************************************************/ @@ -3275,8 +3313,8 @@ namespace VCSnonideal { ncTrial = m_numComponents; int numPreDeleted = m_numRxnTot - m_numRxnRdc; if (numPreDeleted != (m_numSpeciesTot - m_numSpeciesRdc)) { - plogf("we shouldn't be here\n"); - exit(-1); + plogf("vcs_basopt:: we shouldn't be here\n"); + std::exit(-1); } m_numRxnTot = m_numSpeciesTot - ncTrial; m_numRxnRdc = m_numRxnTot - numPreDeleted; @@ -4434,6 +4472,33 @@ namespace VCSnonideal { m_totalMolNum = sum; return m_totalMolNum; } + +#ifdef DEBUG_MODE + void VCS_SOLVE::check_tmoles() const { + int i; + double sum, m_tPhaseMoles_old_a; + //vcs_VolPhase *Vphase; + for (i = 0; i < m_numPhases; i++) { + m_tPhaseMoles_old_a = TPhInertMoles[i]; + + for (int k = 0; k < m_numSpeciesTot; k++) { + if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM) { + if (m_phaseID[k] == i) { + m_tPhaseMoles_old_a += m_molNumSpecies_old[k]; + } + } + } + sum += m_tPhaseMoles_old_a; + + double denom = m_tPhaseMoles_old[i]+ m_tPhaseMoles_old_a + 1.0E-19; + if (!vcs_doubleEqual(m_tPhaseMoles_old[i]/denom, m_tPhaseMoles_old_a/denom)) { + plogf("check_tmoles: we have found a problem with phase %d: %20.15g, %20.15g\n", + i, m_tPhaseMoles_old[i], m_tPhaseMoles_old_a); + //std::exit(-1); + } + } + } +#endif /*****************************************************************************/ // This routine uploads the state of the system into all of the @@ -5052,8 +5117,8 @@ namespace VCSnonideal { // If it isn't, we don't know what's happening if (m_molNumSpecies_old[kspec] != 0.0) { w_kspec = 0.0; - plogf("we shouldn't be here\n"); - exit(-1); + plogf("vcs_birthGuess:: we shouldn't be here\n"); + std::exit(-1); } #endif int ss = m_SSPhase[kspec];