From 13f151a3586671f136a4f5cd0d8906cf275c7aa8 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Mon, 25 Oct 2010 22:12:55 +0000 Subject: [PATCH] Modification to the structure of solve_TP() Took the special case of SS species deletion in RxnStep and moved it into the main routine. It now counts as an interation step. I'm trying to simplify the structure of the program. --- Cantera/src/equil/vcs_VolPhase.cpp | 13 +- Cantera/src/equil/vcs_rxnadj.cpp | 90 +- Cantera/src/equil/vcs_solve.h | 12 +- Cantera/src/equil/vcs_solve_TP.cpp | 1302 +++++++++++++++------------- 4 files changed, 756 insertions(+), 661 deletions(-) diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index a6ca97434..51d4515a0 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -539,7 +539,16 @@ namespace VCSnonideal { m_existence = VCS_PHASE_EXIST_NO; } } + double fractotal = 1.0; v_totalMoles = totalMoles; + if (m_totalMolesInert > 0.0) { + if (m_totalMolesInert > v_totalMoles) { + printf("vcs_VolPhase::setMolesFractionsState: inerts greater than total: %g %g\n", + v_totalMoles, m_totalMolesInert); + exit(EXIT_FAILURE); + } + fractotal = 1.0 - m_totalMolesInert/v_totalMoles; + } double sum = 0.0; for (int k = 0; k < m_numSpecies; k++) { Xmol_[k] = moleFractions[k]; @@ -549,9 +558,9 @@ namespace VCSnonideal { printf("vcs_VolPhase::setMolesFractionsState: inappropriate usage\n"); exit(EXIT_FAILURE); } - if (sum != 1.0) { + if (sum != fractotal) { for (int k = 0; k < m_numSpecies; k++) { - Xmol_[k] /= sum; + Xmol_[k] *= (fractotal /sum); } } _updateMoleFractionDependencies(); diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index 7e68813a4..da1bd9bba 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -43,8 +43,9 @@ namespace VCSnonideal { * in this routine. The species is a noncomponent * - 2 : Same as one but, the zeroed species is a component. */ - int VCS_SOLVE::vcs_RxnStepSizes() { - int j, irxn, kspec, soldel = 0, iph; + int VCS_SOLVE::vcs_RxnStepSizes(int & forceComponentCalc, int &kSpecial) { + int j, irxn, kspec, iph; + int iphDel = -1; double s, xx, dss; int k = 0; vcs_VolPhase *Vphase = 0; @@ -315,24 +316,25 @@ namespace VCSnonideal { * and the code will recover. */ #ifdef DEBUG_MODE - sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", - m_deltaMolNumSpecies[kspec], - -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str()); + sprintf(ANOTE, "Delta damped from %g to %g due to delete %s", + m_deltaMolNumSpecies[kspec], + -m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str()); #endif - m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec]; + m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec]; #ifdef DEBUG_MODE - if (m_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); - } + if (m_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); + } #endif - continue; + continue; } /* * Delete the single species phase */ +#ifdef OLDSTUFF m_molNumSpecies_old[kspec] += dss; m_tPhaseMoles_old[m_phaseID[kspec]] += dss; for (j = 0; j < m_numComponents; ++j) { @@ -351,29 +353,53 @@ namespace VCSnonideal { exit(EXIT_FAILURE); } } +#else + + for (j = 0; j < m_numSpeciesTot; j++) { + m_deltaMolNumSpecies[j] = 0.0; + } + m_deltaMolNumSpecies[kspec] = dss; + for (j = 0; j < m_numComponents; ++j) { + m_deltaMolNumSpecies[j] = dss * m_stoichCoeffRxnMatrix[irxn][j]; + } + + iphDel = m_phaseID[k]; + kSpecial = k; + #ifdef DEBUG_MODE + if (k != kspec) { + sprintf(ANOTE, "Delete component SS phase %d named %s - SS phases only", + iphDel, m_speciesName[k].c_str()); + } else { + sprintf(ANOTE, "Delete this SS phase %d - SS components only", iphDel); + } if (m_debug_print_lvl >= 2) { - plogf(" --- vcs_RxnStepSizes Special section to delete %s", + 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); + plogf(" --- vcs_RxnStepSizes Special section to set up to delete %s", m_speciesName[k].c_str()); plogendl(); } + #endif - /* - * We need to immediately recompute the - * component basis, because we just zeroed - * it out. - */ - soldel = 1; if (k != kspec) { - soldel = 2; + forceComponentCalc = 1; #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { - plogf(" --- Immediate return to get new basis - Restart iteration\n"); + plogf(" --- Force a component recalculation \n"); plogendl(); } #endif - return soldel; } +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" "); vcs_print_line("-", 82); + } +#endif + return iphDel; +#endif } } } /* End of regular processing */ @@ -392,13 +418,12 @@ namespace VCSnonideal { plogf(" "); vcs_print_line("-", 82); } #endif - return soldel; + return iphDel; } - /*****************************************************************************/ - - - //! Calculates reaction adjustments using a full Hessian approximation - /*! + + //==================================================================================================================== + // Calculates reaction adjustments using a full Hessian approximation + /* * Calculates reaction adjustments. This does what equation 6.4-16, p. 143 * in Smith and Missen is suppose to do. However, a full matrix is * formed and then solved via a conjugate gradient algorithm. No @@ -613,8 +638,8 @@ namespace VCSnonideal { #endif return soldel; } - /*****************************************************************************/ + //==================================================================================================================== // Calculates the diagonal contribution to the Hessian due to // the dependence of the activity coefficients on the mole numbers. /* @@ -642,8 +667,8 @@ namespace VCSnonideal { } return diag; } - /*****************************************************************************/ + //==================================================================================================================== //! Calculates the diagonal contribution to the Hessian due to //! the dependence of the activity coefficients on the mole numbers. /*! @@ -680,8 +705,7 @@ namespace VCSnonideal { } return s; } - /*****************************************************************************/ - + //==================================================================================================================== // Recalculate all of the activity coefficients in all of the phases // based on input mole numbers /* diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 6100de9e5..2b43c55a6 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -549,13 +549,13 @@ public: * Special branching occurs sometimes. This causes the component basis * to be reevaluated * - * @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. + * @param forceComponentRecalc integer flagging whether a component recalculation needs + * to be carried out. + * @param kSpecial species number of phase being zeroed. + * + * @return Returns an int representing which phase may need to be zeroed */ - int vcs_RxnStepSizes(); + int vcs_RxnStepSizes(int & forceComponentCalc, int & kSpecial); //! Calculates the total number of moles of species in all phases. /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 8348a9ebb..556cd0b93 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -122,6 +122,8 @@ namespace VCSnonideal { double *dnPhase_irxn; double atomComp; int iphasePop; + int forceComponentCalc = 1; + int iphaseDelete; /* integer that determines which phase is being deleted */ std::vector phasePopPhaseIDs(0); #ifdef DEBUG_MODE char ANOTE[128]; @@ -285,6 +287,8 @@ namespace VCSnonideal { // Update the phase objects with the contents of the soln vector vcs_updateVP(VCS_STATECALC_OLD); vcs_deltag(0, false, VCS_STATECALC_OLD); + // Turn off the force componentCalc flag + forceComponentCalc = 0; if (conv) { goto L_RETURN_BLOCK; @@ -377,6 +381,9 @@ namespace VCSnonideal { check_tmoles(); #endif vcs_tmoles(); + /*************************************************************************/ + /************** COPY OLD into NEW and ZERO VECTORS ***********************/ + /*************************************************************************/ /* * Copy the old solution into the new solution as an initial guess */ @@ -387,15 +394,16 @@ namespace VCSnonideal { vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_new), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); - /* Go find a new reaction adjustment -> - * i.e., change in extent of reaction for each reaction. - * + /* * Zero out the entire vector of updates. We sometimes would * query these values below, and we want to be sure that no * information is left from previous iterations. */ vcs_dzero(VCS_DATA_PTR(m_deltaMolNumSpecies), m_numSpeciesTot); + /*************************************************************************/ + /************** DETERMINE IF DEAD PHASES POP INTO EXISTENCE **************/ + /*************************************************************************/ /* * First step is a major branch in the algorithm. * We first determine if a phase pops into existence. @@ -418,18 +426,21 @@ namespace VCSnonideal { #endif } } + + /*************************************************************************/ + /* DETERMINE THE REACTION STEP SIZES FOR MAIN STEP AND IF PHASES DIE *****/ + /*************************************************************************/ + /* + * Don't do this step if there is a phase pop + */ + iphaseDelete = -1; if (iphasePop < 0) { /* * Figure out the new reaction step sizes * for the major species (do minor species in the future too) */ - - soldel = vcs_RxnStepSizes(); - - if (soldel == 2) { - goto L_COMPONENT_CALC; - } - + kspec = -1; + iphaseDelete = vcs_RxnStepSizes(forceComponentCalc, kspec); } #ifdef DEBUG_MODE else { @@ -478,657 +489,702 @@ namespace VCSnonideal { * * * - */ + */ + if (iphaseDelete >= 0) { #ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf(" --- Main Loop Treatment of each non-component species "); - if (iti == 0) plogf("- Full Calculation:\n"); - else plogf("- Major Components Calculation:\n"); - plogf(" --- Species IC "); - plogf(" KMoles Tent_KMoles Rxn_Adj | Comment \n"); - } + if (m_debug_print_lvl >= 2) { + plogf(" --- Main Loop Treatment -> Circumvented due to Phase Deletion "); + plogendl(); + } #endif - - for (irxn = 0; irxn < m_numRxnRdc; irxn++) { - kspec = m_indexRxnToSpecies[irxn]; - sc_irxn = m_stoichCoeffRxnMatrix[irxn]; - iph = m_phaseID[kspec]; - Vphase = m_VolPhaseList[iph]; -#ifdef DEBUG_MODE - ANOTE[0] = '\0'; -#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 { - dx = 0.0; - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; + + for (k = 0; k < m_numSpeciesTot; k++) { + m_molNumSpecies_new[k] = m_molNumSpecies_old[k] + m_deltaMolNumSpecies[k]; + iph = m_phaseID[k]; + m_tPhaseMoles_new[iph] += m_deltaMolNumSpecies[k]; + } + if (kspec >= m_numComponents) { + if (m_molNumSpecies_new[k] != 0.0) { + printf("vcs_solve_tp:: we shouldn't be here!\n"); + exit(EXIT_FAILURE); } - } else { + if (m_SSPhase[kspec] == 1) { + m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; + } else { + printf("vcs_solve_tp:: we shouldn't be here!\n"); + exit(EXIT_FAILURE); + } + ++m_numRxnMinorZeroed; + allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc); + } + /* + * Set the flags indicating the mole numbers in the vcs_VolPhase + * objects are out of date. + */ + vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); + + /* + * Calculate the new chemical potentials using the tentative + * solution values. We only calculate a subset of these, because + * we have only updated a subset of the W(). + */ + vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); + + /* + * Evaluate DeltaG for all components if ITI=0, and for + * major components only if ITI NE 0 + */ + vcs_deltag(0, false, VCS_STATECALC_NEW); + } else { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- Main Loop Treatment of each non-component species "); + if (iti == 0) plogf("- Full Calculation:\n"); + else plogf("- Major Components Calculation:\n"); + plogf(" --- Species IC "); + plogf(" KMoles Tent_KMoles Rxn_Adj | Comment \n"); + } +#endif + for (irxn = 0; irxn < m_numRxnRdc; irxn++) { + kspec = m_indexRxnToSpecies[irxn]; + sc_irxn = m_stoichCoeffRxnMatrix[irxn]; + iph = m_phaseID[kspec]; + Vphase = m_VolPhaseList[iph]; +#ifdef DEBUG_MODE + ANOTE[0] = '\0'; +#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 { + dx = 0.0; + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; + } + } else { - if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) { - /********************************************************************/ - /************************ VOLTAGE SPECIES ***************************/ - /********************************************************************/ + if (m_speciesStatus[kspec] == VCS_SPECIES_INTERFACIALVOLTAGE) { + /********************************************************************/ + /************************ VOLTAGE SPECIES ***************************/ + /********************************************************************/ #ifdef DEBUG_MODE - dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE); + dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE); #else - dx = vcs_minor_alt_calc(kspec, irxn, &soldel); + 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 = (m_deltaMolNumSpecies[kspec] > 0.0); -#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]); + m_deltaMolNumSpecies[kspec] = dx; } + else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) { + /********************************************************************/ + /********************** ZEROED OUT SPECIES **************************/ + /********************************************************************/ + bool resurrect = (m_deltaMolNumSpecies[kspec] > 0.0); +#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 - if (m_deltaGRxn_new[irxn] >= 0.0 || !resurrect) { - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec]; - m_deltaMolNumSpecies[kspec] = 0.0; - resurrect = false; + if (m_deltaGRxn_new[irxn] >= 0.0 || !resurrect) { + 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) { - if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) { - sprintf(ANOTE, "Species stays zeroed even though dg neg due to " - "STOICH/PHASEPOP constraint: DG = %11.4E", - m_deltaGRxn_new[irxn]); - } else { - sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed", - m_deltaGRxn_new[irxn]); - } - } -#endif - } else { - for (int j = 0; j < m_numElemConstraints; ++j) { - int elType = m_elType[j]; - if (elType == VCS_ELEM_TYPE_ABSPOS) { - 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; - } - - 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"); -#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]); - } + sprintf(ANOTE, "Species stays zeroed: DG = %11.4E", m_deltaGRxn_new[irxn]); + if (m_deltaGRxn_new[irxn] < 0.0) { + if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) { + sprintf(ANOTE, "Species stays zeroed even though dg neg due to " + "STOICH/PHASEPOP constraint: DG = %11.4E", + m_deltaGRxn_new[irxn]); } else { - wx[j] = m_molNumSpecies_old[j]; + sprintf(ANOTE, "Species stays zeroed even though dg neg: DG = %11.4E, ds zeroed", + m_deltaGRxn_new[irxn]); } } - 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; + } else { + for (int j = 0; j < m_numElemConstraints; ++j) { + int elType = m_elType[j]; + if (elType == VCS_ELEM_TYPE_ABSPOS) { + 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; + } + 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_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"); +#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 - m_speciesStatus[kspec] = VCS_SPECIES_ZEROEDSS; - ++m_numRxnMinorZeroed; - allMinorZeroedSpecies = (m_numRxnMinorZeroed == m_numRxnRdc); + /* + * 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; - 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; +#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 (k = 0; k < m_numComponents; ++k) { - m_deltaMolNumSpecies[k] = 0.0; - } - for (iph = 0; iph < m_numPhases; iph++) { - m_deltaPhaseMoles[iph] = 0.0; - } + 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 */ - } - /***********************************************************************/ - /****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/ - /***********************************************************************/ - if (dx != 0.0 && (m_speciesUnknownType[kspec] != - VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { - /* - * Change the amount of the component compounds according - * to the reaction delta that we just computed. - * This should keep the amount of material constant. - */ + }/* End of Loop on ic[irxn] -> the type of species */ + } + /***********************************************************************/ + /****** CALCULATE KMOLE NUMBER CHANGE FOR THE COMPONENT BASIS **********/ + /***********************************************************************/ + if (dx != 0.0 && (m_speciesUnknownType[kspec] != + VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { + /* + * Change the amount of the component compounds according + * to the reaction delta that we just computed. + * This should keep the amount of material constant. + */ #ifdef DEBUG_MODE - if (fabs(m_deltaMolNumSpecies[kspec] -dx) > - 1.0E-14*(fabs(m_deltaMolNumSpecies[kspec]) + fabs(dx) + 1.0E-32)) { - plogf(" ds[kspec] = %20.16g dx = %20.16g , kspec = %d\n", - m_deltaMolNumSpecies[kspec], dx, kspec); - plogf("we have a problem!"); + if (fabs(m_deltaMolNumSpecies[kspec] -dx) > + 1.0E-14*(fabs(m_deltaMolNumSpecies[kspec]) + fabs(dx) + 1.0E-32)) { + plogf(" ds[kspec] = %20.16g dx = %20.16g , kspec = %d\n", + m_deltaMolNumSpecies[kspec], dx, kspec); + plogf("we have a problem!"); + plogendl(); + exit(EXIT_FAILURE); + } +#endif + for (k = 0; k < m_numComponents; ++k) { + m_deltaMolNumSpecies[k] += sc_irxn[k] * dx; + } + /* + * Calculate the tentative change in the total number of + * moles in all of the phases + */ + dnPhase_irxn = m_deltaMolNumPhase[irxn]; + for (iph = 0; iph < m_numPhases; iph++) { + m_deltaPhaseMoles[iph] += dx * dnPhase_irxn[iph]; + } + } + +#ifdef DEBUG_MODE + checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), + VCS_DATA_PTR(m_deltaPhaseMoles), kspec+1); +#endif + /* + * Branch point for returning - + */ +#ifndef DEBUG_MODE + L_MAIN_LOOP_END: ; +#endif +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; + plogf(" --- "); plogf("%-12.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(); + } + L_MAIN_LOOP_END_NO_PRINT: ; +#endif + if (doPhaseDeleteIph != -1) { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- "); + plogf("%-12.12s Main Loop Special Case deleting phase with species: ", + m_speciesName[doPhaseDeleteKspec].c_str()); + plogendl(); + } +#endif + break; + } + } /**************** END OF MAIN LOOP OVER FORMATION REACTIONS ************/ + +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + for (k = 0; k < m_numComponents; k++) { + plogf(" --- "); plogf("%-12.12s", m_speciesName[k].c_str()); + plogf(" c%11.4E%11.4E%11.4E |\n", + m_molNumSpecies_old[k], + m_molNumSpecies_old[k]+m_deltaMolNumSpecies[k], m_deltaMolNumSpecies[k]); + } + plogf(" "); vcs_print_line("-", 80); + plogf(" --- Finished Main Loop"); + plogendl(); + } +#endif + + /*************************************************************************/ + /*********** LIMIT REDUCTION OF BASIS SPECIES TO 99% *********************/ + /*************************************************************************/ + /* + * We have a tentative m_deltaMolNumSpecies[]. Now apply other criteria + * to limit it's magnitude. + * + * + */ + par = 0.5; + for (k = 0; k < m_numComponents; ++k) { + if (m_molNumSpecies_old[k] > 0.0) { + xx = -m_deltaMolNumSpecies[k] / m_molNumSpecies_old[k]; + if (par < xx) { + par = xx; +#ifdef DEBUG_MODE + ll = k; +#endif + } + } else { + if (m_deltaMolNumSpecies[k] < 0.0) { + /* + * If we are here, we then do a step which violates element + * conservation. + */ + iph = m_phaseID[k]; + m_deltaPhaseMoles[iph] -= m_deltaMolNumSpecies[k]; + m_deltaMolNumSpecies[k] = 0.0; + } + } + } + par = 1.0 / par; + if (par <= 1.01 && par > 0.0) { + /* Reduce the size of the step by the multiplicative factor, par */ + par *= 0.99; +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- Reduction in step size due to component "); + plogf("%s", m_speciesName[ll].c_str()); + plogf(" going negative = %11.3E", par); + plogendl(); + } +#endif + for (i = 0; i < m_numSpeciesTot; ++i) { + m_deltaMolNumSpecies[i] *= par; + } + for (iph = 0; iph < m_numPhases; iph++) { + m_deltaPhaseMoles[iph] *= par; + } + } else { + par = 1.0; + } +#ifdef DEBUG_MODE + checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), + VCS_DATA_PTR(m_deltaPhaseMoles), m_numSpeciesTot); +#endif + + /* + * Now adjust the wt[kspec]'s so that the reflect the decrease in + * the overall length of m_deltaMolNumSpecies[kspec] just calculated. At the end + * of this section wt[], m_deltaMolNumSpecies[], tPhMoles, and tPhMoles1 should all be + * consistent with a new estimate of the state of the system. + */ + for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) { + m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; + if (m_molNumSpecies_new[kspec] < 0.0 && (m_speciesUnknownType[kspec] + != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { + plogf("vcs_solve_TP: ERROR on step change wt[%d:%s]: %g < 0.0", + kspec, m_speciesName[kspec].c_str(), m_molNumSpecies_new[kspec]); plogendl(); exit(EXIT_FAILURE); } -#endif - for (k = 0; k < m_numComponents; ++k) { - m_deltaMolNumSpecies[k] += sc_irxn[k] * dx; - } - /* - * Calculate the tentative change in the total number of - * moles in all of the phases - */ - dnPhase_irxn = m_deltaMolNumPhase[irxn]; - for (iph = 0; iph < m_numPhases; iph++) { - m_deltaPhaseMoles[iph] += dx * dnPhase_irxn[iph]; - } } - -#ifdef DEBUG_MODE - checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), - VCS_DATA_PTR(m_deltaPhaseMoles), kspec+1); -#endif + /* - * Branch point for returning - + * Calculate the tentative total mole numbers for each phase */ -#ifndef DEBUG_MODE - L_MAIN_LOOP_END: ; -#endif -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; - plogf(" --- "); plogf("%-12.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(); - } - L_MAIN_LOOP_END_NO_PRINT: ; -#endif - if (doPhaseDeleteIph != -1) { -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf(" --- "); - plogf("%-12.12s Main Loop Special Case deleting phase with species: ", - m_speciesName[doPhaseDeleteKspec].c_str()); - plogendl(); - } -#endif - break; - } - } /**************** END OF MAIN LOOP OVER FORMATION REACTIONS ************/ - -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - for (k = 0; k < m_numComponents; k++) { - plogf(" --- "); plogf("%-12.12s", m_speciesName[k].c_str()); - plogf(" c%11.4E%11.4E%11.4E |\n", - m_molNumSpecies_old[k], - m_molNumSpecies_old[k]+m_deltaMolNumSpecies[k], m_deltaMolNumSpecies[k]); - } - plogf(" "); vcs_print_line("-", 80); - plogf(" --- Finished Main Loop"); - plogendl(); - } -#endif - - /*************************************************************************/ - /*********** LIMIT REDUCTION OF BASIS SPECIES TO 99% *********************/ - /*************************************************************************/ - /* - * We have a tentative m_deltaMolNumSpecies[]. Now apply other criteria - * to limit it's magnitude. - * - * - */ - par = 0.5; - for (k = 0; k < m_numComponents; ++k) { - if (m_molNumSpecies_old[k] > 0.0) { - xx = -m_deltaMolNumSpecies[k] / m_molNumSpecies_old[k]; - if (par < xx) { - par = xx; -#ifdef DEBUG_MODE - ll = k; -#endif - } - } else { - if (m_deltaMolNumSpecies[k] < 0.0) { - /* - * If we are here, we then do a step which violates element - * conservation. - */ - iph = m_phaseID[k]; - m_deltaPhaseMoles[iph] -= m_deltaMolNumSpecies[k]; - m_deltaMolNumSpecies[k] = 0.0; - } - } - } - par = 1.0 / par; - if (par <= 1.01 && par > 0.0) { - /* Reduce the size of the step by the multiplicative factor, par */ - par *= 0.99; -#ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf(" --- Reduction in step size due to component "); - plogf("%s", m_speciesName[ll].c_str()); - plogf(" going negative = %11.3E", par); - plogendl(); - } -#endif - for (i = 0; i < m_numSpeciesTot; ++i) { - m_deltaMolNumSpecies[i] *= par; - } for (iph = 0; iph < m_numPhases; iph++) { - m_deltaPhaseMoles[iph] *= par; + m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + m_deltaPhaseMoles[iph]; } - } else { - par = 1.0; - } -#ifdef DEBUG_MODE - checkDelta1(VCS_DATA_PTR(m_deltaMolNumSpecies), - VCS_DATA_PTR(m_deltaPhaseMoles), m_numSpeciesTot); -#endif - - /* - * Now adjust the wt[kspec]'s so that the reflect the decrease in - * the overall length of m_deltaMolNumSpecies[kspec] just calculated. At the end - * of this section wt[], m_deltaMolNumSpecies[], tPhMoles, and tPhMoles1 should all be - * consistent with a new estimate of the state of the system. - */ - for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) { - m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + m_deltaMolNumSpecies[kspec]; - if (m_molNumSpecies_new[kspec] < 0.0 && (m_speciesUnknownType[kspec] - != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) { - plogf("vcs_solve_TP: ERROR on step change wt[%d:%s]: %g < 0.0", - kspec, m_speciesName[kspec].c_str(), m_molNumSpecies_new[kspec]); + + /* + * Set the flags indicating the mole numbers in the vcs_VolPhase + * objects are out of date. + */ + vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); + + /* + * Calculate the new chemical potentials using the tentative + * solution values. We only calculate a subset of these, because + * we have only updated a subset of the W(). + */ + vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); + + /* + * Evaluate DeltaG for all components if ITI=0, and for + * major components only if ITI NE 0 + */ + vcs_deltag(0, false, VCS_STATECALC_NEW); + + /* *************************************************************** */ + /* **** CONVERGENCE FORCER SECTION ******************************* */ + /* *************************************************************** */ + if (printDetails) { + plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n", + vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old), + VCS_DATA_PTR(m_tPhaseMoles_old))); + plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E", + vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new), + VCS_DATA_PTR(m_tPhaseMoles_new))); plogendl(); - exit(EXIT_FAILURE); } - } - - /* - * Calculate the tentative total mole numbers for each phase - */ - for (iph = 0; iph < m_numPhases; iph++) { - m_tPhaseMoles_new[iph] = m_tPhaseMoles_old[iph] + m_deltaPhaseMoles[iph]; - } - /* - * Set the flags indicating the mole numbers in the vcs_VolPhase - * objects are out of date. - */ - vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW); + forced = vcs_globStepDamp(); - /* - * Calculate the new chemical potentials using the tentative - * solution values. We only calculate a subset of these, because - * we have only updated a subset of the W(). - */ - vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesTot); - - /* - * Evaluate DeltaG for all components if ITI=0, and for - * major components only if ITI NE 0 - */ - vcs_deltag(0, false, VCS_STATECALC_NEW); - - /* *************************************************************** */ - /* **** CONVERGENCE FORCER SECTION ******************************* */ - /* *************************************************************** */ - if (printDetails) { - plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n", - vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old), - VCS_DATA_PTR(m_tPhaseMoles_old))); - plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E", - vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new), - VCS_DATA_PTR(m_tPhaseMoles_new))); - plogendl(); - } - - forced = vcs_globStepDamp(); - - /* - * Print out the changes to the solution that FORCER produced - */ - if (printDetails && forced) { + /* + * Print out the changes to the solution that FORCER produced + */ + if (printDetails && forced) { - plogf(" -----------------------------------------------------\n"); - plogf(" --- FORCER SUBROUTINE changed the solution:\n"); - plogf(" --- SPECIES Status INIT MOLES TENT_MOLES"); - plogf(" FINAL KMOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n"); - for (i = 0; i < m_numComponents; ++i) { - plogf(" --- %-12.12s", m_speciesName[i].c_str()); - plogf(" %14.6E %14.6E %14.6E\n", m_molNumSpecies_old[i], - m_molNumSpecies_old[i] + m_deltaMolNumSpecies[i], m_molNumSpecies_new[i]); + plogf(" -----------------------------------------------------\n"); + plogf(" --- FORCER SUBROUTINE changed the solution:\n"); + plogf(" --- SPECIES Status INIT MOLES TENT_MOLES"); + plogf(" FINAL KMOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n"); + for (i = 0; i < m_numComponents; ++i) { + plogf(" --- %-12.12s", m_speciesName[i].c_str()); + plogf(" %14.6E %14.6E %14.6E\n", m_molNumSpecies_old[i], + m_molNumSpecies_old[i] + m_deltaMolNumSpecies[i], m_molNumSpecies_new[i]); + } + for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) { + irxn = kspec - m_numComponents; + plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); + plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", m_speciesStatus[kspec], + m_molNumSpecies_old[kspec], + m_molNumSpecies_old[kspec]+m_deltaMolNumSpecies[kspec], + m_molNumSpecies_new[kspec], m_deltaGRxn_old[irxn], + m_deltaGRxn_tmp[irxn], m_deltaGRxn_new[irxn]); + } + print_space(26); + plogf("Norms of Delta G():%14.6E%14.6E\n", + l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)), + l2normdg(VCS_DATA_PTR(m_deltaGRxn_new))); + plogf(" Total kmoles of gas = %15.7E\n", m_tPhaseMoles_old[0]); + if ((m_numPhases > 1) && (! (m_VolPhaseList[1])->m_singleSpecies)) { + plogf(" Total kmoles of liquid = %15.7E\n", m_tPhaseMoles_old[1]); + } else { + plogf(" Total kmoles of liquid = %15.7E\n", 0.0); + } + plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n", + vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new), + VCS_DATA_PTR(m_tPhaseMoles_new))); + plogf(" -----------------------------------------------------"); + plogendl(); } - for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) { - irxn = kspec - m_numComponents; - plogf(" --- %-12.12s", m_speciesName[kspec].c_str()); - plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", m_speciesStatus[kspec], - m_molNumSpecies_old[kspec], - m_molNumSpecies_old[kspec]+m_deltaMolNumSpecies[kspec], - m_molNumSpecies_new[kspec], m_deltaGRxn_old[irxn], - m_deltaGRxn_tmp[irxn], m_deltaGRxn_new[irxn]); - } - print_space(26); - plogf("Norms of Delta G():%14.6E%14.6E\n", - l2normdg(VCS_DATA_PTR(m_deltaGRxn_old)), - l2normdg(VCS_DATA_PTR(m_deltaGRxn_new))); - plogf(" Total kmoles of gas = %15.7E\n", m_tPhaseMoles_old[0]); - if ((m_numPhases > 1) && (! (m_VolPhaseList[1])->m_singleSpecies)) { - plogf(" Total kmoles of liquid = %15.7E\n", m_tPhaseMoles_old[1]); - } else { - plogf(" Total kmoles of liquid = %15.7E\n", 0.0); - } - plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n", - vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_new), - VCS_DATA_PTR(m_tPhaseMoles_new))); - plogf(" -----------------------------------------------------"); - plogendl(); } - /* *************************************************************** */ /* **** ITERATION SUMMARY PRINTOUT SECTION *********************** */ /* *************************************************************** */ @@ -1595,7 +1651,10 @@ namespace VCSnonideal { MajorSpeciesHaveConverged = false; /* * Go back and do another iteration with variable ITI - */ + */ + if (forceComponentCalc) { + goto L_COMPONENT_CALC; + } goto L_MAINLOOP_MM4_SPECIES; } } @@ -1661,6 +1720,9 @@ namespace VCSnonideal { * to the main loop to do another iteration. */ iti = 0; + if (forceComponentCalc) { + goto L_COMPONENT_CALC; + } goto L_MAINLOOP_ALL_SPECIES; } } @@ -2098,9 +2160,9 @@ namespace VCSnonideal { double delta = *delta_ptr; #ifdef DEBUG_MODE if (irxn < 0) { - plogf(" --- delete_species() ERROR: called for a component %d", kspec); - plogendl(); - exit(EXIT_FAILURE); + plogf(" --- delete_species() ERROR: called for a component %d", kspec); + plogendl(); + exit(EXIT_FAILURE); } #endif if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { @@ -2585,8 +2647,8 @@ namespace VCSnonideal { for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { iph = m_phaseID[kspec]; m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) - - m_lnMnaughtSpecies[kspec] - + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]); + - m_lnMnaughtSpecies[kspec] + + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]); } /* @@ -2747,8 +2809,8 @@ namespace VCSnonideal { for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { iph = m_phaseID[kspec]; m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) - - m_lnMnaughtSpecies[kspec] - + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]); + - m_lnMnaughtSpecies[kspec] + + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]); } /* * Recalculate the DeltaG's of the formation reactions for the @@ -2769,7 +2831,7 @@ namespace VCSnonideal { if (m_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); + m_speciesName[kspec].c_str(), kspec, dx); } #endif if (dx > 1.0E-50) {