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

Fixed a serious error in InterfaceKinetics that had crept into this branch.

Browse source 
Started upgrade/update of equilibrium solver
       - Fixed numerical jacobian calculation of the ln activity coefficients.
This commit is contained in:
Harry Moffat 2012-11-29 02:18:05 +00:00
parent 64d252fdbf
commit c6f7abbf2f
13 changed files with 730 additions and 652 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
ext/SConscript
View file

@ -26,8 +26,7 @@ def prep_f2c(env):
if env['VERBOSE'] :
print "INFO 2: prep_f2c: adding MSDOS to CPPDEFINES"
if env['VERBOSE'] :
localenv.Append(CPPDEFINES=['added'])
if env['VERBOSE'] :
print "INFO 2: CPPDEFINES ", localenv['CPPDEFINES']
#print "localenv", localenv.Dump()
#exit(0)

6
include/cantera/equil/vcs_VolPhase.h
View file

@ -916,13 +916,13 @@ private:
mutable std::vector<double> ActCoeff;
//! Vector of the derivatives of the ln activity coefficient wrt to the
//! current mole number
//! current mole number multiplied by the current phase moles
/*!
* dLnActCoeffdMolNumber[j][k];
* np_dLnActCoeffdMolNumber[j][k];
* j = id of the species mole number
* k = id of the species activity coefficient
*/
mutable DoubleStarStar dLnActCoeffdMolNumber;
mutable DoubleStarStar np_dLnActCoeffdMolNumber;
//! Status
/*!

17
include/cantera/equil/vcs_internal.h
View file

@ -521,6 +521,23 @@ void vcs_print_stringTrunc(const char* str, size_t space, int alignment);
*/
bool vcs_doubleEqual(double d1, double d2);
//! Sorts a vector of ints in place from lowest to the highest values
/*!
* The vector is returned sorted from lowest to highest.
*
* @param x Reference to a vector of ints.
*/
void vcs_heapsort(std::vector<int> &x);
//! Sorts a vector of ints and eliminates duplicates from the resulting list
/*!
* @param xOrderedUnique Ordered vector of unique ints that were part of the original list
* @param x Reference to a constant vector of ints.
*/
void vcs_orderedUnique(std::vector<int> & xOrderedUnique, const std::vector<int> & x);
}
#endif

11
include/cantera/equil/vcs_solve.h
View file

@ -1958,16 +1958,15 @@ public:
*/
std::vector<double> m_actCoeffSpecies_old;
//! Change in activity coefficient with mole number
//! Change in the log of the activity coefficient with respect to the mole number
//! multiplied by the phase mole number
/*!
* length = [nspecies][nspecies]
*
* (This is a temporary array that
* gets regenerated every time it's
* needed. It is not swapped wrt species
* (unused atm)
* This is a temporary array that gets regenerated every time it's
* needed. It is not swapped wrt species.
*/
DoubleStarStar m_dLnActCoeffdMolNum;
DoubleStarStar m_np_dLnActCoeffdMolNum;
//! Molecular weight of each species
/*!

2
include/cantera/kinetics/InterfaceKinetics.h
View file

@ -556,7 +556,7 @@ public:
* @param iphase Index of the phase. This is the order within the internal thermo vector object
* @param exists Boolean indicating whether the phase exists or not
*/
void setPhaseExistence(const size_t iphase, const bool exists);
void setPhaseExistence(const size_t iphase, const int exists);
//! Set the stability of a phase in the reaction object
/*!

5
include/cantera/thermo/Phase.h
View file

@ -543,9 +543,8 @@ public:
//! Note the density of a phase is an independent variable.
//! @param[in] density density (kg/m^3).
virtual void setDensity(const doublereal density_) {
if (density_ <= 0) {
throw CanteraError("Phase::setDensity",
"density must be positive");
if (density_ <= 0.0) {
throw CanteraError("Phase::setDensity()", "density must be positive");
}
m_dens = density_;
}

43
src/equil/vcs_VolPhase.cpp
View file

@ -213,7 +213,7 @@ vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b)
StarMolarVol = b.StarMolarVol;
PartialMolarVol = b.PartialMolarVol;
ActCoeff = b.ActCoeff;
dLnActCoeffdMolNumber = b.dLnActCoeffdMolNumber;
np_dLnActCoeffdMolNumber = b.np_dLnActCoeffdMolNumber;
m_vcsStateStatus = b.m_vcsStateStatus;
m_phi = b.m_phi;
m_UpToDate = false;
@ -311,7 +311,7 @@ void vcs_VolPhase::resize(const size_t phaseNum, const size_t nspecies,
StarMolarVol.resize(nspecies, -1.0);
PartialMolarVol.resize(nspecies, -1.0);
ActCoeff.resize(nspecies, 1.0);
dLnActCoeffdMolNumber.resize(nspecies, nspecies, 0.0);
np_dLnActCoeffdMolNumber.resize(nspecies, nspecies, 0.0);
m_speciesUnknownType.resize(nspecies, VCS_SPECIES_TYPE_MOLNUM);
@ -914,6 +914,11 @@ double vcs_VolPhase::_updateVolPM() const
void vcs_VolPhase::_updateLnActCoeffJac()
{
double phaseTotalMoles = v_totalMoles;
if (phaseTotalMoles < 1.0E-14) {
phaseTotalMoles = 1.0;
}
/*
* Evaluate the current base activity coefficients if necessary
*/
@ -923,15 +928,15 @@ void vcs_VolPhase::_updateLnActCoeffJac()
if (!TP_ptr) {
return;
}
TP_ptr->getdlnActCoeffdlnN(m_numSpecies, &dLnActCoeffdMolNumber[0][0]);
TP_ptr->getdlnActCoeffdlnN(m_numSpecies, &np_dLnActCoeffdMolNumber[0][0]);
for (size_t j = 0; j < m_numSpecies; j++) {
double moles_j_base = v_totalMoles * Xmol_[j];
double* const lnActCoeffCol = dLnActCoeffdMolNumber[j];
double moles_j_base = phaseTotalMoles * Xmol_[j];
double* const np_lnActCoeffCol = np_dLnActCoeffdMolNumber[j];
if (moles_j_base < 1.0E-200) {
moles_j_base = 1.0E-7 * moles_j_base + 1.0E-20 * v_totalMoles + 1.0E-150;
moles_j_base = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150;
}
for (size_t k = 0; k < m_numSpecies; k++) {
lnActCoeffCol[k] /= moles_j_base;
np_lnActCoeffCol[k] = np_lnActCoeffCol[k] * phaseTotalMoles / moles_j_base;
}
}
@ -939,7 +944,7 @@ void vcs_VolPhase::_updateLnActCoeffJac()
// Make copies of ActCoeff and Xmol_ for use in taking differences
std::vector<double> ActCoeff_Base(ActCoeff);
std::vector<double> Xmol_Base(Xmol_);
double TMoles_base = v_totalMoles;
double TMoles_base = phaseTotalMoles;
/*
* Loop over the columns species to be deltad
@ -950,17 +955,17 @@ void vcs_VolPhase::_updateLnActCoeffJac()
* -> Note Xmol_[] and Tmoles are always positive or zero
* quantities.
*/
double moles_j_base = v_totalMoles * Xmol_Base[j];
deltaMoles_j = 1.0E-7 * moles_j_base + 1.0E-13 * v_totalMoles + 1.0E-150;
double moles_j_base = phaseTotalMoles * Xmol_Base[j];
deltaMoles_j = 1.0E-7 * moles_j_base + 1.0E-13 * phaseTotalMoles + 1.0E-150;
/*
* Now, update the total moles in the phase and all of the
* mole fractions based on this.
*/
v_totalMoles = TMoles_base + deltaMoles_j;
phaseTotalMoles = TMoles_base + deltaMoles_j;
for (size_t k = 0; k < m_numSpecies; k++) {
Xmol_[k] = Xmol_Base[k] * TMoles_base / v_totalMoles;
Xmol_[k] = Xmol_Base[k] * TMoles_base / phaseTotalMoles;
}
Xmol_[j] = (moles_j_base + deltaMoles_j) / v_totalMoles;
Xmol_[j] = (moles_j_base + deltaMoles_j) / phaseTotalMoles;
/*
* Go get new values for the activity coefficients.
@ -971,12 +976,12 @@ void vcs_VolPhase::_updateLnActCoeffJac()
/*
* Calculate the column of the matrix
*/
double* const lnActCoeffCol = dLnActCoeffdMolNumber[j];
double* const np_lnActCoeffCol = np_dLnActCoeffdMolNumber[j];
for (size_t k = 0; k < m_numSpecies; k++) {
double tmp;
tmp = (ActCoeff[k] - ActCoeff_Base[k]) /
((ActCoeff[k] + ActCoeff_Base[k]) * 0.5 * deltaMoles_j);
if (fabs(tmp - lnActCoeffCol[k]) > 1.0E-4 * fabs(tmp) + fabs(lnActCoeffCol[k])) {
if (fabs(tmp - np_lnActCoeffCol[k]) > 1.0E-4 * fabs(tmp) + fabs(np_lnActCoeffCol[k])) {
// printf(" we have an error\n");
}
@ -1014,7 +1019,7 @@ void vcs_VolPhase::_updateLnActCoeffJac()
* k = id of the species activity coefficient
*/
void
vcs_VolPhase::sendToVCS_LnActCoeffJac(double* const* const LnACJac_VCS)
vcs_VolPhase::sendToVCS_LnActCoeffJac(double* const* const np_LnACJac_VCS)
{
/*
* update the Ln Act Coeff jacobian entries with respect to the
@ -1028,11 +1033,11 @@ vcs_VolPhase::sendToVCS_LnActCoeffJac(double* const* const LnACJac_VCS)
*/
for (size_t j = 0; j < m_numSpecies; j++) {
size_t jglob = IndSpecies[j];
double* const lnACJacVCS_col = LnACJac_VCS[jglob];
const double* const lnACJac_col = dLnActCoeffdMolNumber[j];
double* const np_lnACJacVCS_col = np_LnACJac_VCS[jglob];
const double* const np_lnACJac_col = np_dLnActCoeffdMolNumber[j];
for (size_t k = 0; k < m_numSpecies; k++) {
size_t kglob = IndSpecies[k];
lnACJacVCS_col[kglob] = lnACJac_col[k];
np_lnACJacVCS_col[kglob] = np_lnACJac_col[k];
}
}
}

2
src/equil/vcs_root1d.cpp
View file

@ -487,7 +487,7 @@ QUAD_BAIL:
}
#ifdef DEBUG_MODE
if (printLvl >= 3) {
fprintf(fp, "\nvcs_root1d failure in %d its\n", its);
fprintf(fp, "\nvcs_root1d failure in %lu its\n", its);
}
#endif
}

198
src/equil/vcs_rxnadj.cpp
View file

@ -16,8 +16,7 @@
#include <cstdlib>
#include <cmath>
namespace VCSnonideal
{
namespace VCSnonideal {
// Calculates formation reaction step sizes.
/*
@ -64,7 +63,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
}
plogf("\n");
plogf(" --- Species KMoles Rxn_Adjustment DeltaG"
" | Comment\n");
" | Comment\n");
}
#endif
/*
@ -74,13 +73,13 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
if (m_useActCoeffJac) {
vcs_CalcLnActCoeffJac(VCS_DATA_PTR(m_molNumSpecies_old));
}
/************************************************************************
******** LOOP OVER THE FORMATION REACTIONS *****************************
************************************************************************/
/************************************************************************
******** LOOP OVER THE FORMATION REACTIONS *****************************
************************************************************************/
for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Normal Calc");
sprintf(ANOTE, "Normal Calc");
#endif
kspec = m_indexRxnToSpecies[irxn];
@ -94,7 +93,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
dnPhase_irxn = m_deltaMolNumPhase[irxn];
if (m_molNumSpecies_old[kspec] == 0.0 && (! m_SSPhase[kspec])) {
if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) {
/********************************************************************/
/******* MULTISPECIES PHASE WITH total moles equal to zero *********/
/*******************************************************************/
@ -111,44 +110,45 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
iph = m_phaseID[kspec];
double tphmoles = m_tPhaseMoles_old[iph];
double trphmoles = tphmoles / m_totalMolNum;
if (trphmoles > VCS_DELETE_PHASE_CUTOFF) {
Vphase = m_VolPhaseList[iph];
if (Vphase->exists() && (trphmoles > VCS_DELETE_PHASE_CUTOFF)) {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
m_deltaMolNumSpecies[kspec] = 0.0;
#ifdef DEBUG_MODE
sprintf(ANOTE,
"MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "MultSpec (%s): Species not born due to STOICH/PHASEPOP even though DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
#endif
} else {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.0;
#ifdef DEBUG_MODE
sprintf(ANOTE,
"MultSpec (%s): small species born again DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "MultSpec (%s): small species born again DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
#endif
}
} else {
#ifdef DEBUG_MODE
sprintf(ANOTE, "MultSpec (%s): phase come alive DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "MultSpec (%s):still dead, no phase pop, even though DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
#endif
Vphase = m_VolPhaseList[iph];
size_t numSpPhase = Vphase->nSpecies();
m_deltaMolNumSpecies[kspec] =
m_totalMolNum * 10.0 * VCS_DELETE_PHASE_CUTOFF / numSpPhase;
m_deltaMolNumSpecies[kspec] = 0.0;
if (Vphase->exists() > 0 && trphmoles > 0.0) {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES * 10.;
#ifdef DEBUG_MODE
sprintf(ANOTE,
"MultSpec (%s): birthed species because it was zero in a small existing phase with DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15), m_deltaGRxn_new[irxn]);
#endif
}
}
} else {
#ifdef DEBUG_MODE
sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E",
vcs_speciesType_string(m_speciesStatus[kspec], 15),
sprintf(ANOTE, "MultSpec (%s): still dead DG = %11.3E", vcs_speciesType_string(m_speciesStatus[kspec], 15),
m_deltaGRxn_new[irxn]);
#endif
m_deltaMolNumSpecies[kspec] = 0.0;
}
} else {
/********************************************************************/
@ -163,12 +163,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
*/
if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "Skipped: superconverged DG = %11.3E", m_deltaGRxn_new[irxn]);
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);
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
#endif
continue;
@ -179,13 +179,12 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
*/
if ((m_speciesStatus[kspec] != VCS_SPECIES_MAJOR) && (m_deltaGRxn_new[irxn] >= 0.0)) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E",
m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
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);
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
#endif
continue;
@ -196,7 +195,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
if (m_SSPhase[kspec]) {
s = 0.0;
} else {
s = 1.0 / m_molNumSpecies_old[kspec] ;
s = 1.0 / m_molNumSpecies_old[kspec];
}
for (size_t j = 0; j < m_numComponents; ++j) {
if (!m_SSPhase[j]) {
@ -207,7 +206,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
}
for (size_t j = 0; j < m_numPhases; j++) {
Vphase = m_VolPhaseList[j];
if (! Vphase->m_singleSpecies) {
if (!Vphase->m_singleSpecies) {
if (m_tPhaseMoles_old[j] > 0.0) {
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
}
@ -225,7 +224,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
#ifdef DEBUG_MODE
if (s_old != s) {
sprintf(ANOTE, "Normal calc: diag adjusted from %g "
"to %g due to act coeff", s_old, s);
"to %g due to act coeff", s_old, s);
}
#endif
}
@ -235,20 +234,20 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
for (size_t j = 0; j < m_numComponents; ++j) {
double stoicC = m_stoichCoeffRxnMatrix[irxn][j];
if (stoicC != 0.0) {
double negChangeComp = - stoicC * m_deltaMolNumSpecies[kspec];
double negChangeComp = -stoicC * m_deltaMolNumSpecies[kspec];
if (negChangeComp > m_molNumSpecies_old[j]) {
if (m_molNumSpecies_old[j] > 0.0) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %d (%10s) going neg", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
"to %g due to component %lu (%10s) going neg", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str());
#endif
m_deltaMolNumSpecies[kspec] = - m_molNumSpecies_old[j] / stoicC;
m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[j] / stoicC;
} else {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %d (%10s) zero", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j]/stoicC, j, m_speciesName[j].c_str());
"to %g due to component %lu (%10s) zero", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[j] / stoicC, j, m_speciesName[j].c_str());
#endif
m_deltaMolNumSpecies[kspec] = 0.0;
}
@ -259,8 +258,8 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
if (-m_deltaMolNumSpecies[kspec] > m_molNumSpecies_old[kspec]) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to %s going negative", m_deltaMolNumSpecies[kspec],
-m_molNumSpecies_old[kspec], m_speciesName[kspec].c_str());
"to %g due to %s going negative", m_deltaMolNumSpecies[kspec], -m_molNumSpecies_old[kspec],
m_speciesName[kspec].c_str());
#endif
m_deltaMolNumSpecies[kspec] = -m_molNumSpecies_old[kspec];
}
@ -321,17 +320,16 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
* 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];
#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);
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
#endif
continue;
@ -352,18 +350,18 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
#ifdef DEBUG_MODE
if (k != kspec) {
sprintf(ANOTE, "Delete component SS phase %d named %s - SS phases only",
iphDel, m_speciesName[k].c_str());
sprintf(ANOTE, "Delete component SS phase %lu named %s - SS phases only", iphDel,
m_speciesName[k].c_str());
} else {
sprintf(ANOTE, "Delete this SS phase %d - SS components only", iphDel);
sprintf(ANOTE, "Delete this SS phase %lu - SS components only", iphDel);
}
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);
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());
m_speciesName[k].c_str());
plogendl();
}
#endif
@ -390,8 +388,8 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
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);
m_molNumSpecies_old[kspec], m_deltaMolNumSpecies[kspec],
m_deltaGRxn_new[irxn], ANOTE);
}
#endif
} /* End of loop over m_speciesUnknownType */
@ -457,13 +455,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
*/
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Normal Calc");
sprintf(ANOTE, "Normal Calc");
#endif
kspec = m_indexRxnToSpecies[irxn];
dnPhase_irxn = m_deltaMolNumPhase[irxn];
if (m_molNumSpecies_old[kspec] == 0.0 && (! m_SSPhase[kspec])) {
if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) {
/* *******************************************************************/
/* **** MULTISPECIES PHASE WITH total moles equal to zero ************/
/* *******************************************************************/
@ -497,11 +495,11 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
*/
if (fabs(m_deltaGRxn_new[irxn]) <= m_tolmaj2) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "Skipped: converged DG = %11.3E\n", m_deltaGRxn_new[irxn]);
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
#endif
continue;
}
@ -511,12 +509,11 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
*/
if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR && m_deltaGRxn_new[irxn] >= 0.0) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E\n",
m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
sprintf(ANOTE, "Skipped: IC = %3d and DG >0: %11.3E\n", m_speciesStatus[kspec], m_deltaGRxn_new[irxn]);
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
m_deltaMolNumSpecies[kspec], ANOTE);
#endif
continue;
}
@ -529,7 +526,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
s = 1.0 / m_molNumSpecies_old[kspec];
}
for (j = 0; j < m_numComponents; ++j) {
if (! m_SSPhase[j]) {
if (!m_SSPhase[j]) {
s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
}
}
@ -589,10 +586,10 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
*/
if (dss != 0.0) {
m_molNumSpecies_old[kspec] += dss;
m_tPhaseMoles_old[m_phaseID[kspec]] += dss;
m_tPhaseMoles_old[m_phaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j];
m_tPhaseMoles_old[m_phaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
m_tPhaseMoles_old[m_phaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
}
m_molNumSpecies_old[k] = 0.0;
m_tPhaseMoles_old[m_phaseID[k]] = 0.0;
@ -619,7 +616,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
plogf(" --- ");
plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec],
m_deltaMolNumSpecies[kspec], ANOTE);
m_deltaMolNumSpecies[kspec], ANOTE);
#endif
} /* End of loop over non-component stoichiometric formation reactions */
@ -664,10 +661,10 @@ double VCS_SOLVE::vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal)
}
if (hessActCoef >= 0.0) {
diag += hessActCoef;
} else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) {
} else if (fabs(hessActCoef) < 0.6666 * hessianDiag_Ideal) {
diag += hessActCoef;
} else {
diag -= 0.6666 * hessianDiag_Ideal;
diag -= 0.6666 * hessianDiag_Ideal;
}
return diag;
}
@ -687,11 +684,15 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn)
double* sc_irxn;
kspec = m_indexRxnToSpecies[irxn];
kph = m_phaseID[kspec];
double np_kspec = m_tPhaseMoles_old[kph];
if (np_kspec < 1.0E-13) {
np_kspec = 1.0E-13;
}
sc_irxn = m_stoichCoeffRxnMatrix[irxn];
/*
* First the diagonal term of the Jacobian
*/
s = m_dLnActCoeffdMolNum[kspec][kspec];
s = m_np_dLnActCoeffdMolNum[kspec][kspec] / np_kspec;
/*
* Next, the other terms. Note this only a loop over the components
* So, it's not too expensive to calculate.
@ -700,12 +701,16 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn)
if (!m_SSPhase[l]) {
for (k = 0; k < m_numComponents; ++k) {
if (m_phaseID[k] == m_phaseID[l]) {
s += sc_irxn[k] * sc_irxn[l] * m_dLnActCoeffdMolNum[k][l];
double np = m_tPhaseMoles_old[m_phaseID[k]];
if (np > 0.0) {
s += sc_irxn[k] * sc_irxn[l] * m_np_dLnActCoeffdMolNum[k][l] / np;
}
}
}
if (kph == m_phaseID[l]) {
s += sc_irxn[l] * (m_dLnActCoeffdMolNum[kspec][l] + m_dLnActCoeffdMolNum[l][kspec]);
s += sc_irxn[l] * (m_np_dLnActCoeffdMolNum[kspec][l] + m_np_dLnActCoeffdMolNum[l][kspec]) / np_kspec;
}
}
}
return s;
@ -739,7 +744,7 @@ void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS)
* -> This scatter calculation is carried out in the
* vcs_VolPhase object.
*/
Vphase->sendToVCS_LnActCoeffJac(m_dLnActCoeffdMolNum.baseDataAddr());
Vphase->sendToVCS_LnActCoeffJac(m_np_dLnActCoeffdMolNum.baseDataAddr());
}
}
}
@ -778,9 +783,8 @@ void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS)
* Note, this is a dangerous routine that leaves the underlying objects in
* an unknown state.
*/
double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc,
const size_t irxn, const double* const molNum,
double* const ac, double* const mu_i)
double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc, const size_t irxn, const double* const molNum, double* const ac,
double* const mu_i)
{
size_t kspec = irxn + m_numComponents;
int* pp_ptr = m_phaseParticipation[irxn];
@ -812,8 +816,7 @@ double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc,
*
* @return Returns the optimized step length found by the search
*/
double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig,
char* const ANOTE)
double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char* const ANOTE)
#else
double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
#endif
@ -833,9 +836,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
* Calculate the deltaG value at the dx = 0.0 point
*/
vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
double deltaGOrig = deltaG_Recalc_Rxn(VCS_STATECALC_OLD,
irxn, molNumBase, acBase,
VCS_DATA_PTR(m_feSpecies_old));
double deltaGOrig = deltaG_Recalc_Rxn(VCS_STATECALC_OLD, irxn, molNumBase, acBase, VCS_DATA_PTR(m_feSpecies_old));
double forig = fabs(deltaGOrig) + 1.0E-15;
if (deltaGOrig > 0.0) {
if (dx_orig > 0.0) {
@ -844,7 +845,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
if (m_debug_print_lvl >= 2) {
//plogf(" --- %s :Warning possible error dx>0 dg > 0\n", SpName[kspec]);
}
sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx>0 dg > 0");
sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx>0 dg > 0");
#endif
return dx;
}
@ -855,7 +856,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
if (m_debug_print_lvl >= 2) {
//plogf(" --- %s :Warning possible error dx<0 dg < 0\n", SpName[kspec]);
}
sprintf(ANOTE,"Rxn reduced to zero step size in line search: dx<0 dg < 0");
sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx<0 dg < 0");
#endif
return dx;
}
@ -867,17 +868,16 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
}
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNumBase, m_numSpeciesRdc);
molSum = molNumBase[kspec];
molSum = molNumBase[kspec];
m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx_orig;
for (k = 0; k < m_numComponents; k++) {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx_orig;
molSum += molNumBase[k];
molSum += molNumBase[k];
}
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
double deltaG1 = deltaG_Recalc_Rxn(VCS_STATECALC_NEW,
irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
double deltaG1 = deltaG_Recalc_Rxn(VCS_STATECALC_NEW, irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
/*
* If deltaG hasn't switched signs when going the full distance
@ -892,7 +892,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
* If we have decreased somewhat, the deltaG return after finding
* a better estimate for the line search.
*/
if (fabs(deltaG1) < 0.8*forig) {
if (fabs(deltaG1) < 0.8 * forig) {
if (deltaG1 * deltaGOrig < 0.0) {
slope = (deltaG1 - deltaGOrig) / dx_orig;
dx = -deltaGOrig / slope;
@ -915,9 +915,8 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx;
}
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
double deltaG = deltaG_Recalc_Rxn(VCS_STATECALC_NEW,
irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
double deltaG = deltaG_Recalc_Rxn(VCS_STATECALC_NEW, irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
/*
* If deltaG hasn't switched signs when going the full distance
* then we are heading in the appropriate direction, and
@ -939,19 +938,16 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
}
}
finalize:
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
finalize: vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
if (its >= MAXITS) {
#ifdef DEBUG_MODE
sprintf(ANOTE,"Rxn reduced to zero step size from %g to %g (MAXITS)",
dx_orig, dx);
sprintf(ANOTE, "Rxn reduced to zero step size from %g to %g (MAXITS)", dx_orig, dx);
return dx;
#endif
}
#ifdef DEBUG_MODE
if (dx != dx_orig) {
sprintf(ANOTE,"Line Search reduced step size from %g to %g",
dx_orig, dx);
sprintf(ANOTE, "Line Search reduced step size from %g to %g", dx_orig, dx);
}
#endif

2
src/equil/vcs_solve.cpp
View file

@ -218,7 +218,7 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements,
*/
m_useActCoeffJac = true;
if (m_useActCoeffJac) {
m_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0);
m_np_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0);
}
m_PMVolumeSpecies.resize(nspecies0, 0.0);

4
src/equil/vcs_solve_TP.cpp
View file

@ -2075,7 +2075,7 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete
/*
* get the diagonal of the activity coefficient jacobian
*/
s = m_dLnActCoeffdMolNum[kspec][kspec];
s = m_np_dLnActCoeffdMolNum[kspec][kspec] / (m_tPhaseMoles_old[iph]);
// s *= (m_tPhaseMoles_old[iph]);
/*
* We fit it to a power law approximation of the activity coefficient
@ -5463,7 +5463,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k
std::swap(m_formulaMatrix[j][k1], m_formulaMatrix[j][k2]);
}
if (m_useActCoeffJac) {
vcs_switch2D(m_dLnActCoeffdMolNum.baseDataAddr(), k1, k2);
vcs_switch2D(m_np_dLnActCoeffdMolNum.baseDataAddr(), k1, k2);
}
std::swap(m_speciesStatus[k1], m_speciesStatus[k2]);
/*

1068
src/equil/vcs_util.cpp
View file

File diff suppressed because it is too large Load diff

21
src/kinetics/InterfaceKinetics.cpp
View file

@ -304,6 +304,8 @@ void InterfaceKinetics::getActivityConcentrations(doublereal* const conc)
* Update the equilibrium constants in molar units for all
* reversible reactions. Irreversible reactions have their
* equilibrium constant set to zero.
* For reactions involving charged species the equilibrium
* constant is adjusted according to the electrostatic potential.
*/
void InterfaceKinetics::updateKc()
{
@ -311,7 +313,10 @@ void InterfaceKinetics::updateKc()
//static vector_fp mu(nTotalSpecies());
if (m_nrev > 0) {
/*
* Get the vector of standard state electrochemical potentials for species in the Interfacial
* kinetics object and store it in m_mu0[]
*/
size_t nsp, ik = 0;
doublereal rt = GasConstant*thermo(0).temperature();
doublereal rrt = 1.0 / rt;
@ -327,8 +332,7 @@ void InterfaceKinetics::updateKc()
}
// compute Delta mu^0 for all reversible reactions
m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0),
DATA_PTR(m_rkcn));
m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), DATA_PTR(m_rkcn));
for (size_t i = 0; i < m_nrev; i++) {
size_t irxn = m_revindex[i];
@ -336,6 +340,7 @@ void InterfaceKinetics::updateKc()
throw CanteraError("InterfaceKinetics",
"illegal value: irxn = "+int2str(irxn));
}
// WARNING this may overflow HKM
m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
}
for (size_t i = 0; i != m_nirrev; ++i) {
@ -680,7 +685,7 @@ void InterfaceKinetics::updateROP()
if (m_rxnPhaseIsReactant[j][rp]) {
if (! m_phaseExists[rp]) {
m_ropnet[j] = 0.0;
m_ropr[j] = m_ropr[j] = 0.0;
m_ropr[j] = m_ropf[j] = 0.0;
}
}
}
@ -753,7 +758,7 @@ InterfaceKinetics::adjustRatesForIntermediatePhases()
}
#endif
//=================================================================================================
//=================================================================================================
/*
*
* getDeltaGibbs():
@ -1335,7 +1340,7 @@ solvePseudoSteadyStateProblem(int ifuncOverride, doublereal timeScaleOverride)
}
//================================================================================================
void InterfaceKinetics::setPhaseExistence(const size_t iphase, const bool exists)
void InterfaceKinetics::setPhaseExistence(const size_t iphase, const int exists)
{
if (iphase >= m_thermo.size()) {
throw CanteraError("InterfaceKinetics:setPhaseExistence", "out of bounds");
@ -1343,6 +1348,9 @@ void InterfaceKinetics::setPhaseExistence(const size_t iphase, const bool exists
if (exists) {
if (!m_phaseExists[iphase]) {
m_phaseExistsCheck--;
if (m_phaseExistsCheck < 0) {
m_phaseExistsCheck = 0;
}
m_phaseExists[iphase] = true;
}
m_phaseIsStable[iphase] = true;
@ -1353,6 +1361,7 @@ void InterfaceKinetics::setPhaseExistence(const size_t iphase, const bool exists
}
m_phaseIsStable[iphase] = false;
}
}
//================================================================================================
// Gets the phase existence int for the ith phase