diff --git a/Cantera/src/equil/equil.h b/Cantera/src/equil/equil.h index ed1ff7801..c4211aff8 100644 --- a/Cantera/src/equil/equil.h +++ b/Cantera/src/equil/equil.h @@ -72,7 +72,7 @@ namespace Cantera { * @ingroup equil */ int equilibrate(thermo_t& s, const char* XY, - int solver = -1, doublereal rtol = 1.0e-9, int maxsteps = 1000, + int solver = -1, doublereal rtol = 1.0e-9, int maxsteps = 5000, int maxiter = 100, int loglevel = -99); //! Equilibrate a MultiPhase object @@ -101,7 +101,7 @@ namespace Cantera { * @ingroup equil */ doublereal equilibrate(MultiPhase& s, const char* XY, - doublereal rtol = 1.0e-9, int maxsteps = 1000, int maxiter = 100, + doublereal rtol = 1.0e-9, int maxsteps = 5000, int maxiter = 100, int loglevel = -99); } diff --git a/Cantera/src/equil/vcs_Gibbs.cpp b/Cantera/src/equil/vcs_Gibbs.cpp index b30597037..18c8dd29b 100644 --- a/Cantera/src/equil/vcs_Gibbs.cpp +++ b/Cantera/src/equil/vcs_Gibbs.cpp @@ -23,7 +23,8 @@ namespace VCSnonideal { /*****************************************************************************/ /*****************************************************************************/ -double VCS_SOLVE::vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles) +double VCS_SOLVE::vcs_Total_Gibbs(double *molesSp, double *chemPot, + double *tPhMoles) /************************************************************************* * @@ -37,9 +38,8 @@ double VCS_SOLVE::vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles) *************************************************************************/ { double g = 0.0; - int kspec; - int iph; - for (iph = 0; iph < NPhase; iph++) { + + for (int iph = 0; iph < NPhase; iph++) { vcs_VolPhase *Vphase = VPhaseList[iph]; if ((TPhInertMoles[iph] > 0.0) && (tPhMoles[iph] > 0.0)) { g += TPhInertMoles[iph] * @@ -49,8 +49,9 @@ double VCS_SOLVE::vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles) } } } - for (kspec = 0; kspec < m_numSpeciesRdc; ++kspec) { - g += w[kspec] * fe[kspec]; + + for (int kspec = 0; kspec < m_numSpeciesRdc; ++kspec) { + g += molesSp[kspec] * chemPot[kspec]; } return g; diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 725f3087f..c9748b6da 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -123,6 +123,7 @@ namespace VCSnonideal { dg.resize(nspecies0, 0.0); dgl.resize(nspecies0, 0.0); + m_deltaGRxn_tmp.resize(nspecies0, 0.0); ds.resize(nspecies0, 0.0); fel.resize(nspecies0, 0.0); diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 0c85cc6c8..de5d887e7 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -294,6 +294,7 @@ private: ); int force(int iti); + int globStepDamp(int iti); void vcs_switch2D(double * const * const Jac, int k1, int k2); double l2normdg(double dg[]); #ifdef DEBUG_MODE @@ -398,7 +399,15 @@ public: * handled by the alt_min treatment or * should be handled as a major species. */ - std::vector scSize; + std::vector scSize; + + //! Standard state chemical potentials for species K at the current + //! temperature and pressure. + /*! + * The first NC entries are for components. The following NR entries are + * for the current non-component species in the mechanism. + */ + std::vector ff; //! Dimensionless/Dimensional free energy for all the species in the mechanism at the //! current T, P, and mole numbers. @@ -409,13 +418,11 @@ public: */ std::vector m_gibbsSpecies; - //! Standard state chemical potentials for species K at the current - //! temperature and pressure. + //! Old free energy vector from the previous iteration /*! - * The first NC entries are for components. The following NR entries are - * for the current non-component species in the mechanism. + * fe[] is copied into fel[] */ - std::vector ff; + std::vector fel; //! Dimensionless trial free energy for all the species in the mechanism //! at the current T, P, and mole numbers. @@ -430,7 +437,7 @@ public: * Initial estimate: 0 user estimate * -1 machine estimate */ - int iest; + int iest; //! Total moles of the species /*! @@ -469,7 +476,7 @@ public: //! electric potential of the iph phase std::vector phasePhi; - //! Tentative value of the mole number vector. It's also used to store the + //! Tentative value of the mole number vector. It's also used to store the //! mole fraction vector. std::vector wt; @@ -487,7 +494,9 @@ public: std::vector dg; //! Last deltag[irxn] from the previous step - std::vector dgl; + std::vector dgl; + + std::vector m_deltaGRxn_tmp; //! Reaction Adjustments for each species /*! @@ -495,8 +504,7 @@ public: */ std::vector ds; - std::vector fel; /* fel[k] = Old Free Energy vector from the previous - * iteration. fe[] is copied into fel[] */ + std::vector ga; /* ga[j] = Element abundances for jth element from * estimate * -> this is calculated from the current mole diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 634427e84..fc57e78a2 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -1003,10 +1003,9 @@ namespace VCSnonideal { if (vcs_debug_print_lvl >= 2) { for (k = 0; k < m_numComponents; k++) { plogf(" --- "); plogf("%-12.12s", SpName[k].c_str()); - plogf(" c%11.4E%11.4E%11.4E |", + plogf(" c%11.4E%11.4E%11.4E |\n", soln[k], soln[k]+ds[k], ds[k]); } - plogendl(); plogf(" "); vcs_print_line("-", 80); plogf(" --- Finished Main Loop"); plogendl(); @@ -1105,7 +1104,7 @@ namespace VCSnonideal { /* * Print Intermediate results */ - // HKM Actually always need to calculate this + // Actually always need to calculate this // or else nonprintouts get different results and sometimes // fail in the line search algorithm -> Why is this? vcs_dfe(VCS_DATA_PTR(wt), 1, 1, 0, m_numSpeciesRdc); @@ -1168,7 +1167,10 @@ namespace VCSnonideal { plogf(" --- %18s = %15.7E\n", Vphase->PhaseName.c_str(), TPhMoles1[iph]); } plogf(" "); vcs_print_line("-", 103); - plogf(" --- Total Dimensionless Gibbs Free Energy = %15.7E", + plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n", + vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(fel), + VCS_DATA_PTR(TPhMoles))); + plogf(" --- Total New Dimensionless Gibbs Free Energy = %20.13E", vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_gibbsSpecies), VCS_DATA_PTR(TPhMoles1))); plogendl(); @@ -1187,50 +1189,28 @@ namespace VCSnonideal { /* *************************************************************** */ /* **** CONVERGENCE FORCER SECTION ******************************* */ /* *************************************************************** */ - /* - * Save the previous delta G in the old vector for - * printout purposes - */ - if (printDetails) { - vcs_dcopy(VCS_DATA_PTR(dgl), VCS_DATA_PTR(dg), m_numRxnRdc); - } - forced = FALSE; - // if (! im && ! MajorSpeciesHaveConverged) { - forced = force(iti); - //} + + forced = globStepDamp(iti); + /* * Print out the changes to the solution that FORCER produced */ if (printDetails && forced) { - if (iti != 0) { -#ifdef DEBUG_MODE - if (vcs_debug_print_lvl >= 3) { - plogf(" *** vcs_dfe for printout only:"); - } -#endif - vcs_updateVP(0); - vcs_dfe(VCS_DATA_PTR(soln), 0, 1, 0, m_numSpeciesRdc); -#ifdef DEBUG_MODE - if (vcs_debug_print_lvl >= 3) { - plogf(" *** vcs_deltag call for printouts only;"); - } -#endif - vcs_deltag(1, false); - } plogf(" -----------------------------------------------------\n"); plogf(" --- FORCER SUBROUTINE changed the solution:\n"); - plogf(" --- SPECIES Status TENT MOLES"); - plogf(" FINAL MOLES TENT_DEL_G/RT FINAL_DELTA_G/RT\n"); + plogf(" --- SPECIES Status INIT MOLES TENT_MOLES"); + plogf(" FINAL MOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n"); for (i = 0; i < m_numComponents; ++i) { plogf(" --- %-12.12s", SpName[i].c_str()); - plogf(" %14.6E%14.6E\n", wt[i], soln[i]); + plogf(" %14.6E %14.6E %14.6E\n", soln[i], soln[i] + ds[i], wt[i]); } for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) { irxn = kspec - m_numComponents; plogf(" --- %-12.12s", SpName[kspec].c_str()); - plogf(" %2d %14.6E%14.6E%14.6E%14.6E\n", spStatus[irxn], - wt[kspec], soln[kspec], dgl[irxn], dg[irxn]); + plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", spStatus[irxn], + soln[kspec], soln[kspec]+ds[kspec], wt[kspec], dgl[irxn], + m_deltaGRxn_tmp[irxn], dg[irxn]); } print_space(26); plogf("Norms of Delta G():%14.6E%14.6E\n", @@ -1242,9 +1222,9 @@ namespace VCSnonideal { } else { plogf(" Total moles of liquid = %15.7E\n", 0.0); } - plogf(" Total Dimensionless Gibbs Free Energy = %15.7E\n", - vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_gibbsSpecies), - VCS_DATA_PTR(TPhMoles))); + plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n", + vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_gibbsSpecies), + VCS_DATA_PTR(TPhMoles1))); plogf(" -----------------------------------------------------"); plogendl(); } @@ -1263,10 +1243,11 @@ namespace VCSnonideal { * we have already done this inside the FORCED * loop. */ - if (! forced) { - vcs_dcopy(VCS_DATA_PTR(TPhMoles), VCS_DATA_PTR(TPhMoles1), NPhase); - vcs_dcopy(VCS_DATA_PTR(soln), VCS_DATA_PTR(wt), m_numSpeciesRdc); - } + vcs_dcopy(VCS_DATA_PTR(TPhMoles), VCS_DATA_PTR(TPhMoles1), NPhase); + vcs_dcopy(VCS_DATA_PTR(soln), VCS_DATA_PTR(wt), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(dgl), VCS_DATA_PTR(dg), m_numRxnRdc); + vcs_dcopy(VCS_DATA_PTR(fel), VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesRdc); + vcs_updateVP(0); /* * Increment the iteration counters @@ -1275,8 +1256,9 @@ namespace VCSnonideal { ++it1; #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { - plogf(" --- Increment counter increased, step is accepted: %4d\n", + plogf(" --- Increment counter increased, step is accepted: %4d", m_VCount->Its); + plogendl(); } #endif /*************************************************************************/ @@ -2096,7 +2078,7 @@ namespace VCSnonideal { /*****************************************************************************/ /*****************************************************************************/ - int VCS_SOLVE::delta_species(int kspec, double *delta_ptr) + int VCS_SOLVE::delta_species(int kspec, double *delta_ptr) /************************************************************************ * @@ -2620,92 +2602,87 @@ namespace VCSnonideal { vcs_deltag(0, true); } - /*****************************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ - - int VCS_SOLVE::force(int iti) - - /************************************************************************** - * - * force: - * - * Convergence Forcer: - * - * This routine optimizes the minimization of the total gibbs free - * energy: - * Gibbs = sum_k( fe_k * w_k ) - * along the current direction ds[], by choosing a value, al: (0 0), - * does this code section kick in. It finds the point on the parabola - * where the slope is equal to zero. - * - * NOTE: The algorithm used to find the slope is not quite accurate. - * The term, sum_k( (fe_k_n - fe_k_n-1) * w_k_n-1 ) - * is dropped from s1, and, the term, - * sum_k( (fe_k_n - fe_k_n-1) * w_k_n ), is dropped from s2 - *************************************************************************/ - { + /* globalStepDamp + * + * Convergence Forcer: + * + * This routine optimizes the minimization of the total gibbs free + * energy by making sure the slope of the following functional stays + * negative: + * + * d_Gibbs/ds = sum_k( m_deltaGRxn * ds[k] ) + * + * along the current direction ds[], by choosing a value, al: (0 0), + * does this code section kick in. It finds the point on the parabola + * where the slope is equal to zero. + * + */ + int VCS_SOLVE::globStepDamp(int iti) { double s1, s2, al; - int i, iph; - double *dptr = VCS_DATA_PTR(m_gibbsSpecies); - //int numSpeciesRdc = m_numSpeciesRdc; + int irxn, kspec, iph; + double *dptr = VCS_DATA_PTR(dg); /* *************************************************** */ /* **** CALCULATE SLOPE AT END OF THE STEP ********** */ /* *************************************************** */ s2 = 0.0; - for (i = 0; i < m_numSpeciesRdc; ++i) { - s2 += dptr[i] * ds[i]; + for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { + kspec = irxn + m_numComponents; + s2 += dptr[irxn] * ds[kspec]; } -#ifdef DEBUG_NOT - if (s2 <= 0.0) { -#ifdef DEBUG_NOT - if (vcs_debug_print_lvl >= 2) { - plogf(" --- subroutine FORCE produced no adjustments,"); - plogf(" failed s2 test\n"); - } -#endif - return FALSE; - } -#endif -#ifdef DEBUG_MODE - if (vcs_debug_print_lvl >= 2) { - plogf(" --- subroutine FORCE: End Slope = %g\n", s2); - } -#endif + + /* *************************************************** */ /* **** CALCULATE ORIGINAL SLOPE ********************* */ /* ************************************************** */ s1 = 0.0; - dptr = VCS_DATA_PTR(fel); - for (i = 0; i < m_numSpeciesRdc; ++i) { - s1 += dptr[i] * ds[i]; + dptr = VCS_DATA_PTR(dgl); + for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { + kspec = irxn + m_numComponents; + s1 += dptr[irxn] * ds[kspec]; } -#ifdef DEBUG_NOT - if (s1 >= 0.0) { + +#ifdef DEBUG_MODE + if (vcs_debug_print_lvl >= 2) { + plogf(" --- subroutine FORCE: Beginning Slope = %g\n", s1); + plogf(" --- subroutine FORCE: End Slope = %g\n", s2); + } +#endif + + if (s1 > 0.0) { #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { plogf(" --- subroutine FORCE produced no adjustments,"); - plogf(" failed s1 test -PROBLEM!!\n"); + if (s1 < 1.0E-40) { + plogf(" s1 positive but really small"); + } else { + plogf(" failed s1 test"); + } + plogendl(); } #endif return FALSE; } -#endif + + if (s2 <= 0.0) { #ifdef DEBUG_MODE - if (vcs_debug_print_lvl >= 2) { - plogf(" --- subroutine FORCE: Beginning Slope = %g\n", s1); + if (vcs_debug_print_lvl >= 2) { + plogf(" --- subroutine FORCE produced no adjustments, s2 < 0"); + plogendl(); + } +#endif + return FALSE; } -#endif + /* *************************************************** */ - /* **** FIT PARABOLA ********************************* */ + /* **** FIT PCJ2822ARABOLA ********************************* */ /* *************************************************** */ al = 1.0; if (fabs(s1 -s2) > 1.0E-200) { @@ -2724,17 +2701,25 @@ namespace VCSnonideal { plogf(" --- subroutine FORCE produced a damping factor = %g\n", al); } #endif + /* *************************************************** */ /* **** ADJUST MOLE NUMBERS, CHEM. POT *************** */ /* *************************************************** */ - dptr = VCS_DATA_PTR(soln); - for (i = 0; i < m_numSpeciesRdc; ++i) { - dptr[i] += al * ds[i]; +#ifdef DEBUG_MODE + if (vcs_debug_print_lvl >= 2) { + vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_tmp), VCS_DATA_PTR(dg), + m_numRxnRdc); + } +#endif + + dptr = VCS_DATA_PTR(wt); + for (kspec = 0; kspec < m_numSpeciesRdc; ++kspec) { + wt[kspec] = soln[kspec] + al * ds[kspec]; } for (iph = 0; iph < NPhase; iph++) { - TPhMoles[iph] += al * DelTPhMoles[iph]; + TPhMoles1[iph] = TPhMoles[iph] + al * DelTPhMoles[iph]; } - vcs_updateVP(0); + vcs_updateVP(1); #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { @@ -2748,17 +2733,34 @@ namespace VCSnonideal { * only step is being carried out, then we don't need to * update the minor noncomponents. */ - vcs_dfe(dptr, 0, iti, 0, m_numSpeciesRdc); + // vcs_dfe(dptr, 1, iti, 0, m_numSpeciesRdc); + vcs_dfe(dptr, 1, 0, 0, m_numSpeciesRdc); /* * Evaluate DeltaG for all components if ITI=0, and for * major components only if ITI NE 0 */ - vcs_deltag(iti, false); + // vcs_deltag(iti, false); + vcs_deltag(0, false); + + dptr = VCS_DATA_PTR(dg); + s2 = 0.0; + for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { + kspec = irxn + m_numComponents; + s2 += dptr[irxn] * ds[kspec]; + } + + +#ifdef DEBUG_MODE + if (vcs_debug_print_lvl >= 2) { + plogf(" --- subroutine FORCE: Adj End Slope = %g", s2); + plogendl(); + } +#endif return TRUE; - } /* force() *****************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ + } + + + /* * vcs_RxnStepSizes(): *