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

Added stateCalc to some vcs_VolPhase functions. This is work in progress.

Browse source 
deltaGrxn() is now always based on irxn index.
This commit is contained in:
Harry Moffat 2008-06-16 16:25:01 +00:00
parent 045b4110e8
commit 0bd6ae596d
6 changed files with 96 additions and 49 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
Cantera/src/equil/vcs_VolPhase.cpp
View file

@ -146,7 +146,11 @@ namespace VCSnonideal {
if (&b != this) {
int old_num = NVolSpecies;
m_owningSolverObject = b.m_owningSolverObject;
// Note: we comment this out for the assignment operator
// specifically, because it isn't true for the assignment
// operator but is true for a copy constructor
// m_owningSolverObject = b.m_owningSolverObject;
VP_ID = b.VP_ID;
Domain_ID = b.Domain_ID;
SingleSpecies = b.SingleSpecies;
@ -688,7 +692,9 @@ namespace VCSnonideal {
* in all of the phases in a VCS problem. Only the
* entries for the current phase are filled in.
*/
void vcs_VolPhase::sendToVCS_ActCoeff(double * const AC) const {
void vcs_VolPhase::sendToVCS_ActCoeff(const int stateCalc,
double * const AC) {
updateFromVCS_MoleNumbers(stateCalc);
if (!m_UpToDate_AC) {
evaluateActCoeff();
}

2
Cantera/src/equil/vcs_VolPhase.h
View file

@ -226,7 +226,7 @@ namespace VCSnonideal {
* in all of the phases in a VCS problem. Only the
* entries for the current phase are filled in.
*/
void sendToVCS_ActCoeff(double * const AC) const;
void sendToVCS_ActCoeff(const int stateCalc, double * const AC);
//! set the electric potential of the phase
/*!

37
Cantera/src/equil/vcs_report.cpp
View file

@ -203,22 +203,21 @@ namespace VCSnonideal {
}
plogf(" | DG/RT Rxn |\n");
print_line("-", m_numComponents*10 + 45);
for (i = 0; i < m_numRxnTot; i++) {
int kspec = m_indexRxnToSpecies[i];
for (int irxn = 0; irxn < m_numRxnTot; irxn++) {
int kspec = m_indexRxnToSpecies[irxn];
plogf(" %3d ", kspec);
plogf("%-10.10s", m_speciesName[kspec].c_str());
plogf("|%10.3g |", m_molNumSpecies_old[kspec]);
for (j = 0; j < m_numComponents; j++) {
plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]);
plogf(" %6.2f", m_stoichCoeffRxnMatrix[irxn][j]);
}
// Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
plogf(" |%10.3g |", m_deltaGRxn_new[kspec]);
plogf(" |%10.3g |", m_deltaGRxn_new[irxn]);
plogf("\n");
}
print_line("-", m_numComponents*10 + 45);
plogf("\n");
/*
/*
* ------------------ TABLE OF PHASE INFORMATION ---------------------
*/
std::vector<double> gaPhase(m_numElemConstraints, 0.0);
@ -262,7 +261,7 @@ namespace VCSnonideal {
plogf(" %10.3g", gaPhase[j]);
gaTPhase[j] += gaPhase[j];
}
gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(m_molNumSpecies_old),
gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(m_feSpecies_old));
gibbsTotal += gibbsPhase;
plogf(" | %18.11E |\n", gibbsPhase);
@ -313,8 +312,12 @@ namespace VCSnonideal {
vcs_printChemPotUnits(m_VCS_UnitsFormat);
plogf(")\n");
plogf(" Name TKMoles StandStateChemPot "
" ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)\n");
print_line("-", 115);
" ln(AC) ln(X_i) | F z_i phi | ChemPot | (-lnMnaught)");
#ifdef DEBUG_MODE
plogf("| (MolNum ChemPot)|");
#endif
plogf("\n");
print_line("-", 147);
for (i = 0; i < nspecies; ++i) {
l = sortindex[i];
int pid = m_phaseID[l];
@ -346,12 +349,22 @@ namespace VCSnonideal {
exit(-1);
}
plogf(" %12.4E |", m_feSpecies_old[l]);
if (m_lnMnaughtSpecies[l] != 0.0) {
plogf(" (%14.7E)", - m_lnMnaughtSpecies[l]);
if( m_lnMnaughtSpecies[l] != 0.0) {
plogf("(%11.5E)", - m_lnMnaughtSpecies[l]);
} else {
plogf(" ");
}
#ifdef DEBUG_MODE
plogf("| %20.13E |", m_feSpecies_old[l] * m_molNumSpecies_old[l]);
#endif
plogf("\n");
}
print_line("-", 115);
#ifdef DEBUG_MODE
for (i = 0; i < 125; i++) plogf(" ");
plogf("%20.13E\n", g);
#endif
print_line("-", 147);
/*
* ------------- TABLE OF SOLUTION COUNTERS --------------------------

18
Cantera/src/equil/vcs_rxnadj.cpp
View file

@ -372,13 +372,14 @@ namespace VCSnonideal {
* Note, this is a dangerous routine that leaves the underlying objects in
* an unknown state.
*/
double VCS_SOLVE::deltaG_Recalc_Rxn(const int irxn, const double *const molNum,
double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc,
const int irxn, const double *const molNum,
double * const ac, double * const mu_i) {
int kspec = irxn + m_numComponents;
int *pp_ptr = m_phaseParticipation[irxn];
for (int iphase = 0; iphase < m_numPhases; iphase++) {
if (pp_ptr[iphase]) {
vcs_chemPotPhase(iphase, molNum, ac, mu_i);
vcs_chemPotPhase(stateCalc, iphase, molNum, ac, mu_i);
}
}
double deltaG = mu_i[kspec];
@ -424,7 +425,9 @@ namespace VCSnonideal {
/*
* Calculate the deltaG value at the dx = 0.0 point
*/
double deltaGOrig = deltaG_Recalc_Rxn(irxn, molNumBase, acBase,
vcs_setFlagsVolPhases(false, VCS_STATECALC_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) {
@ -461,8 +464,10 @@ namespace VCSnonideal {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx_orig;
molSum += molNumBase[k];
}
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
double deltaG1 = deltaG_Recalc_Rxn(irxn, VCS_DATA_PTR(m_molNumSpecies_new),
double deltaG1 = deltaG_Recalc_Rxn(VCS_STATECALC_NEW,
irxn, VCS_DATA_PTR(m_molNumSpecies_new),
ac, VCS_DATA_PTR(m_feSpecies_new));
/*
@ -500,7 +505,9 @@ namespace VCSnonideal {
for (k = 0; k < m_numComponents; k++) {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx;
}
double deltaG = deltaG_Recalc_Rxn(irxn, VCS_DATA_PTR(m_molNumSpecies_new),
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));
/*
* If deltaG hasn't switched signs when going the full distance
@ -524,6 +531,7 @@ namespace VCSnonideal {
}
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)",

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

@ -341,7 +341,7 @@ public:
* (VCS species order)
*
*/
void vcs_chemPotPhase(const int iph, const double *const molNum,
void vcs_chemPotPhase(const int stateCalc, const int iph, const double *const molNum,
double * const ac, double * const mu_i,
const bool do_deleted = false);
@ -1301,7 +1301,8 @@ private:
*
* @return Returns the dimensionless deltaG of the reaction
*/
double deltaG_Recalc_Rxn(const int irxn, const double *const molNum,
double deltaG_Recalc_Rxn(const int stateCalc,
const int irxn, const double *const molNum,
double * const ac, double * const mu_i);
//! Delete memory that isn't just resizeable STL containers
@ -1494,16 +1495,13 @@ public:
//std::vector<double> wt;
std::vector<double> m_molNumSpecies_new;
//! Delta G(I) for the noncomponent species in the mechanism.
//! Delta G(irxn) for the noncomponent species in the mechanism.
/*!
* Computed by the subroutine DELTAG. DG is the free
* energy change for the reaction which
* forms species K from the
* component species. This vector has length
* equal to the number of noncomponent
* species in the mechanism. It starts with
* the first current noncomponent species
* in the mechanism.
* Computed by the subroutine deltaG. m_deltaGRxn is the free
* energy change for the reaction which forms species K from the
* component species. This vector has length equal to the number
* of noncomponent species in the mechanism. It starts with
* the first current noncomponent species in the mechanism.
*/
std::vector<double> m_deltaGRxn_new;

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

@ -926,11 +926,28 @@ namespace VCSnonideal {
(doPhaseDeleteIph == -1) &&
(m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
double dx_old = dx;
#ifdef DEBUG_HKM_NOT
std::vector<double> feSpecies_tmp(m_numSpeciesTot);
feSpecies_tmp = m_feSpecies_old;
std::vector<double> molNumSpecies_tmp(m_numSpeciesTot);
molNumSpecies_tmp = m_molNumSpecies_old;
#endif
#ifdef DEBUG_MODE
dx = vcs_line_search(irxn, dx_old, ANOTE);
#else
dx = vcs_line_search(irxn, dx_old);
#endif
#ifdef DEBUG_HKM_NOT
for (int kk = 0; kk < m_numSpeciesTot; kk++) {
if (feSpecies_tmp[kk] != m_feSpecies_old[kk]) {
printf("we are here\n");
}
if (molNumSpecies_tmp[kk] != m_molNumSpecies_old[kk]) {
printf("we are here\n");
}
}
#endif
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
}
m_deltaMolNumSpecies[kspec] = dx;
@ -1942,14 +1959,14 @@ namespace VCSnonideal {
* Store the final Delta G values for each non-component species
* in the species slot rather than the reaction slot
*/
kspec = m_numSpeciesTot;
i = m_numRxnTot;
for (irxn = 0; irxn < m_numRxnTot; ++irxn) {
--kspec;
--i;
m_deltaGRxn_new[kspec] = m_deltaGRxn_new[i];
}
vcs_dzero(VCS_DATA_PTR(m_deltaGRxn_new), m_numComponents);
// kspec = m_numSpeciesTot;
// i = m_numRxnTot;
//for (irxn = 0; irxn < m_numRxnTot; ++irxn) {
// --kspec;
// --i;
// m_deltaGRxn_new[kspec] = m_deltaGRxn_new[i];
//}
// vcs_dzero(VCS_DATA_PTR(m_deltaGRxn_new), m_numComponents);
/*
* Evaluate the final mole fractions
* storring them in wt[]
@ -2182,13 +2199,18 @@ namespace VCSnonideal {
m_molNumSpecies_old[kspec] += dx;
int iph = m_phaseID[kspec];
m_tPhaseMoles_old[iph] += dx;
vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD);
for (j = 0; j < m_numComponents; ++j) {
iph = m_phaseID[j];
tmp = sc_irxn[j] * dx;
m_molNumSpecies_old[j] += tmp;
m_tPhaseMoles_old[iph] += tmp;
if (m_molNumSpecies_old[j] < 0.0) {
m_molNumSpecies_old[j] = 0.0;
if (tmp != 0.0) {
iph = m_phaseID[j];
m_molNumSpecies_old[j] += tmp;
m_tPhaseMoles_old[iph] += tmp;
vcs_setFlagsVolPhase(iph, false, VCS_STATECALC_OLD);
if (m_molNumSpecies_old[j] < 0.0) {
m_molNumSpecies_old[j] = 0.0;
}
}
}
}
@ -2936,7 +2958,7 @@ namespace VCSnonideal {
* only step is being carried out, then we don't need to
* update the minor noncomponents.
*/
vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
vcs_dfe(VCS_STATECALC_NEW, 0, 0, m_numSpeciesRdc);
/*
@ -4388,14 +4410,14 @@ namespace VCSnonideal {
* (VCS species order)
*
*/
void VCS_SOLVE::vcs_chemPotPhase(const int iph, const double *const molNum,
void VCS_SOLVE::vcs_chemPotPhase(const int stateCalc,
const int iph, const double *const molNum,
double * const ac, double * const mu_i,
const bool do_deleted) {
vcs_VolPhase *Vphase = m_VolPhaseList[iph];
int nkk = Vphase->NVolSpecies;
int k, kspec;
int stateCalc = VCS_STATECALC_OLD;
#ifdef DEBUG_MODE
//if (m_debug_print_lvl >= 2) {
@ -4414,7 +4436,7 @@ namespace VCSnonideal {
}
Vphase->setMolesFromVCS(stateCalc, molNum);
Vphase->sendToVCS_ActCoeff(ac);
Vphase->sendToVCS_ActCoeff(stateCalc, ac);
double phi = Vphase->electricPotential();
double Faraday_phi = m_Faraday_dim * phi;
@ -4695,7 +4717,7 @@ namespace VCSnonideal {
Vphase = m_VolPhaseList[iphase];
if (!Vphase->SingleSpecies) {
Vphase->setMolesFromVCS(stateCalc, molNum);
Vphase->sendToVCS_ActCoeff(VCS_DATA_PTR(actCoeff_ptr));
Vphase->sendToVCS_ActCoeff(stateCalc, VCS_DATA_PTR(actCoeff_ptr));
}
m_phasePhi[iphase] = Vphase->electricPotential();
m_phaseACAreCurrent[iphase] = 1;