research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Cleaned up the global damping algorithm. made sure the total gibbs free

Browse source 
energy printouts are monotonically decreasing.
This commit is contained in:
Harry Moffat 2008-04-11 22:44:49 +00:00
parent e463b77fa2
commit 6c7e375020
5 changed files with 150 additions and 138 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
Cantera/src/equil/equil.h
View file

@ -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);
}

13
Cantera/src/equil/vcs_Gibbs.cpp
View file

@ -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;

1
Cantera/src/equil/vcs_solve.cpp
View file

@ -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);

30
Cantera/src/equil/vcs_solve.h
View file

@ -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<double> scSize;
std::vector<double> 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<double> 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<double> 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<double> ff;
std::vector<double> 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<double> 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<double> wt;
@ -487,7 +494,9 @@ public:
std::vector<double> dg;
//! Last deltag[irxn] from the previous step
std::vector<double> dgl;
std::vector<double> dgl;
std::vector<double> m_deltaGRxn_tmp;
//! Reaction Adjustments for each species
/*!
@ -495,8 +504,7 @@ public:
*/
std::vector<double> ds;
std::vector<double> fel; /* fel[k] = Old Free Energy vector from the previous
* iteration. fe[] is copied into fel[] */
std::vector<double> ga; /* ga[j] = Element abundances for jth element from
* estimate
* -> this is calculated from the current mole

240
Cantera/src/equil/vcs_solve_TP.cpp
View file

@ -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<al<1)
* such that the a parabola approximation to Gibbs(al) fit to the
* end points al = 0 and al = 1 is minimizied.
* s1 = slope of Gibbs function at al = 0, which is the previous
* solution = d(Gibbs)/d(al).
* s2 = slope of Gibbs function at al = 1, which is the current
* solution = d(Gibbs)/d(al).
* Only if there has been an inflection point (i.e., s1 < 0 and s2 > 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<al<1)
* such that the a parabola approximation to Gibbs(al) fit to the
* end points al = 0 and al = 1 is minimizied.
* s1 = slope of Gibbs function at al = 0, which is the previous
* solution = d(Gibbs)/d(al).
* s2 = slope of Gibbs function at al = 1, which is the current
* solution = d(Gibbs)/d(al).
* Only if there has been an inflection point (i.e., s1 < 0 and s2 > 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():
*