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

Variable name changes involving activity coefficients

Browse source 
Some changes from m_actCoefSpecies_new to m_actCoeffSpecies_old
This commit is contained in:
Harry Moffat 2008-05-14 20:03:57 +00:00
parent f741f96b3c
commit 66a341d868
5 changed files with 67 additions and 49 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -44,6 +44,7 @@ int VCS_SOLVE::vcs_report(int iconv)
* Print out a report on the state of the equilibrium problem to
* standard output.
* This prints out the current contents of the VCS_SOLVE class, V.
* The "old" solution vector is printed out.
***************************************************************************/
{
int i, j, l, k, inertYes = FALSE, kspec;
@ -315,7 +316,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %-12.12s", m_speciesName[l].c_str());
plogf(" %14.7E ", m_molNumSpecies_old[l]);
plogf("%14.7E ", m_SSfeSpecies[l]);
plogf("%14.7E ", log(ActCoeff[l]));
plogf("%14.7E ", log(m_actCoeffSpecies_old[l]));
double tpmoles = m_tPhaseMoles_old[pid];
double phi = m_phasePhi[pid];
double eContrib = phi * m_chargeSpecies[l] * Faraday_dim;
@ -326,12 +327,12 @@ int VCS_SOLVE::vcs_report(int iconv)
if (tpmoles > 0.0 && m_molNumSpecies_old[l] > 0.0) {
lx = log(m_molNumSpecies_old[l]) - log(tpmoles);
} else {
lx = m_feSpecies_curr[l] - m_SSfeSpecies[l] - log(ActCoeff[l]) + SpecLnMnaught[l];
lx = m_feSpecies_curr[l] - m_SSfeSpecies[l] - log(m_actCoeffSpecies_old[l]) + SpecLnMnaught[l];
}
}
plogf("%14.7E |", lx);
plogf("%14.7E | ", eContrib);
double tmp = m_SSfeSpecies[l] + log(ActCoeff[l]) + lx - SpecLnMnaught[l] + eContrib;
double tmp = m_SSfeSpecies[l] + log(m_actCoeffSpecies_old[l]) + lx - SpecLnMnaught[l] + eContrib;
if (fabs(m_feSpecies_curr[l] - tmp) > 1.0E-8) {
plogf("\n\t\twe have a problem - doesn't add up\n");
exit(-1);

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

@ -416,8 +416,8 @@ namespace VCSnonideal {
double dx = dx_orig;
double *sc_irxn = m_stoichCoeffRxnMatrix[irxn];
double *molNumBase = VCS_DATA_PTR(m_molNumSpecies_old);
double *acBase = VCS_DATA_PTR(ActCoeff0);
double *ac = VCS_DATA_PTR(ActCoeff);
double *acBase = VCS_DATA_PTR(m_actCoeffSpecies_old);
double *ac = VCS_DATA_PTR(m_actCoeffSpecies_new);
double molSum = 0.0;
double slope;
/*

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

@ -181,8 +181,8 @@ namespace VCSnonideal {
SpecActConvention.resize(nspecies0, 0);
PhaseActConvention.resize(nphase0, 0);
SpecLnMnaught.resize(nspecies0, 0.0);
ActCoeff.resize(nspecies0, 1.0);
ActCoeff0.resize(nspecies0, 1.0);
m_actCoeffSpecies_new.resize(nspecies0, 1.0);
m_actCoeffSpecies_old.resize(nspecies0, 1.0);
CurrPhAC.resize(nphase0, 0);
m_wtSpecies.resize(nspecies0, 0.0);
m_chargeSpecies.resize(nspecies0, 0.0);

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

@ -291,7 +291,24 @@ public:
void vcs_tmoles();
void vcs_deltag(int l, bool doDeleted);
void vcs_switch_pos(int ifunc, int k1, int k2);
//! Swaps the indecises for all of the global data for two species, k1
//! and k2.
/*!
*
* @param ifunc: If true, switch the species data and the noncomponent reaction
* data. This must be called for a non-component species only.
* If false, switch the species data only. Typically, we use this
* option when determining the component species and at the
* end of the calculation, when we want to return unscrambled
* results. All rxn data will be out-of-date.
*
* @param k1 First species index
*
* @param k2 Second species index
*/
void vcs_switch_pos(const int ifunc, const int k1, const int k2);
void vcs_deltag_Phase(int iphase, bool doDeleted);
//! birthGuess returns the number of moles of a species
@ -1256,14 +1273,15 @@ public:
*
* Length = number of species
*/
std::vector<double> ActCoeff;
std::vector<double> m_actCoeffSpecies_new;
//! Molar-based Activity Coefficients for Species
//! Molar-based Activity Coefficients for Species based on old mole numbers
/*!
* These activity coefficients are based on the m_molNumSpecies_old values
* Molar based activity coeffients.
* Length = number of species
*/
std::vector<double> ActCoeff0;
std::vector<double> m_actCoeffSpecies_old;
//! Change in activity coefficient with mole number
/*!

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

@ -400,7 +400,7 @@ namespace VCSnonideal {
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(ActCoeff0), VCS_DATA_PTR(ActCoeff), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_old), VCS_DATA_PTR(m_actCoeffSpecies_new), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc);
/* Go find a new reaction adjustment ->
@ -2027,7 +2027,7 @@ namespace VCSnonideal {
goto L_ZERO_SPECIES;
}
} else {
double ac0 = ActCoeff[kspec];
double ac0 = m_actCoeffSpecies_new[kspec];
double ac = ac0;
double w0 = w_kspec;
double dd = exp(-dg_irxn);
@ -2035,8 +2035,8 @@ namespace VCSnonideal {
wTrial = w0 * ac0 / ac * dd;
*wt_kspec = wTrial;
Vphase->setMolesFromVCS(VCS_DATA_PTR(m_molNumSpecies_new));
Vphase->sendToVCSActCoeff(VCS_DATA_PTR(ActCoeff));
double ac1 = ActCoeff[kspec];
Vphase->sendToVCSActCoeff(VCS_DATA_PTR(m_actCoeffSpecies_new));
double ac1 = m_actCoeffSpecies_new[kspec];
double acprime = 0.0;
if (fabs(wTrial - w0) > 1.0E-8 * w0) {
acprime = (ac1 - ac0) / (wTrial - w0);
@ -3249,7 +3249,7 @@ namespace VCSnonideal {
irxn = kspec - m_numComponents;
if (m_deltaGRxn_new[irxn] > 50.0) m_deltaGRxn_new[irxn] = 50.0;
if (m_deltaGRxn_new[irxn] < -50.0) m_deltaGRxn_new[irxn] = -50.0;
poly += exp(-m_deltaGRxn_new[irxn])/ActCoeff[kspec];
poly += exp(-m_deltaGRxn_new[irxn])/m_actCoeffSpecies_new[kspec];
}
/*
* Calculate m_deltaGRxn_new[] for each species in a zeroed multispecies phase.
@ -4159,7 +4159,7 @@ namespace VCSnonideal {
*
* m_feSpecies(I) = m_SSfeSpecies(I)
* + ln(ActCoeff_M[I] * m(I))
* + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
* + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]
* where m[I] is the molality of the ith solute
*
* m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units)
@ -4190,7 +4190,8 @@ namespace VCSnonideal {
* T = temperature
* V = potential of the interface = phi_electrode - phi_solution
*
* For these species, the solution vector unknown, z[I], is V, the phase voltage, in volts.
* For these species, the solution vector unknown, z[I], is V, the phase voltage,
* in volts.
*
* Input
* --------
@ -4224,6 +4225,7 @@ namespace VCSnonideal {
int l1, l2, iph, kspec, irxn;
int iphase;
double *tPhMoles_ptr;
double *actCoeff_ptr;
double *tlogMoles;
vcs_VolPhase *Vphase;
VCS_SPECIES_THERMO *st_ptr;
@ -4256,8 +4258,10 @@ namespace VCSnonideal {
#endif
if (kk <= VCS_STATECALC_OLD) {
tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_old);
actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_old);
} else {
tPhMoles_ptr = VCS_DATA_PTR(m_tPhaseMoles_new);
actCoeff_ptr = VCS_DATA_PTR(m_actCoeffSpecies_new);
}
tlogMoles = VCS_DATA_PTR(TmpPhase);
/*
@ -4313,7 +4317,7 @@ namespace VCSnonideal {
Vphase = VPhaseList[iphase];
if (!Vphase->SingleSpecies) {
Vphase->setMolesFromVCS(z);
Vphase->sendToVCSActCoeff(VCS_DATA_PTR(ActCoeff));
Vphase->sendToVCSActCoeff(VCS_DATA_PTR(actCoeff_ptr));
}
m_phasePhi[iphase] = Vphase->electricPotential();
CurrPhAC[iphase] = 1;
@ -4351,14 +4355,14 @@ namespace VCSnonideal {
iph = PhaseID[kspec];
if (tPhMoles_ptr[iph] > 0.0) {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]
+ log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
+ log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase];
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(actCoeff_ptr[kspec] * z[kspec])
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase];
}
@ -4394,14 +4398,14 @@ namespace VCSnonideal {
iph = PhaseID[kspec];
if (tPhMoles_ptr[iph] > 0.0) {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]
+ log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
+ log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase]; ;
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(actCoeff_ptr[kspec] * z[kspec])
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iphase];
}
@ -4438,14 +4442,14 @@ namespace VCSnonideal {
iph = PhaseID[kspec];
if (tPhMoles_ptr[iph] > 0.0) {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]
+ log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
+ log(actCoeff_ptr[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec];
} else {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
st_ptr = SpeciesThermo[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(actCoeff_ptr[kspec] * z[kspec])
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec];
}
}
@ -4580,7 +4584,7 @@ namespace VCSnonideal {
else {
plogf("we shouldn't be here");
plogendl();
exit(-1);
std::exit(-1);
}
#endif
}
@ -4608,27 +4612,22 @@ namespace VCSnonideal {
}
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
void VCS_SOLVE::vcs_switch_pos(int ifunc, int k1, int k2)
/**************************************************************************
*
* vcs_switch_pos:
*
* Swaps the indecises for all of the global data for two species, k1
* and k2.
*
* ifunc: If true, switch the species data and the noncomponent reaction
* data. This must be called for a non-component species only.
*
* If false, switch the species data only. Typically, we use this
* option when determining the component species and at the
* end of the calculation, when we want to return unscrambled
* results.
*************************************************************************/
{
// Swaps the indecises for all of the global data for two species, k1
// and k2.
/*
* @param ifunc: If true, switch the species data and the noncomponent reaction
* data. This must be called for a non-component species only.
* If false, switch the species data only. Typically, we use this
* option when determining the component species and at the
* end of the calculation, when we want to return unscrambled
* results. All rxn data will be out-of-date.
*
* @param k1 First species index
*
* @param k2 Second species index
*/
void VCS_SOLVE::vcs_switch_pos(const int ifunc, const int k1, const int k2) {
register int j;
register double t1 = 0.0;
int i1, i2, iph, kp1, kp2;
@ -4679,8 +4678,8 @@ namespace VCSnonideal {
SWAP(indPhSp[k1], indPhSp[k2], j);
SWAP(SpecActConvention[k1], SpecActConvention[k2], j);
SWAP(SpecLnMnaught[k1], SpecLnMnaught[k2], t1);
SWAP(ActCoeff[k1], ActCoeff[k2], t1);
SWAP(ActCoeff0[k1], ActCoeff0[k2], t1);
SWAP(m_actCoeffSpecies_new[k1], m_actCoeffSpecies_new[k2], t1);
SWAP(m_actCoeffSpecies_old[k1], m_actCoeffSpecies_old[k2], t1);
SWAP(m_wtSpecies[k1], m_wtSpecies[k2], t1);
SWAP(m_chargeSpecies[k1], m_chargeSpecies[k2], t1);
SWAP(SpeciesThermo[k1], SpeciesThermo[k2], st_tmp);
@ -4854,7 +4853,7 @@ namespace VCSnonideal {
if (iph == iphase) {
if (m_deltaGRxn_new[irxn] > 50.0) m_deltaGRxn_new[irxn] = 50.0;
if (m_deltaGRxn_new[irxn] < -50.0) m_deltaGRxn_new[irxn] = -50.0;
phaseDG -= exp(-m_deltaGRxn_new[irxn])/ActCoeff[kspec];
phaseDG -= exp(-m_deltaGRxn_new[irxn])/m_actCoeffSpecies_new[kspec];
}
}
/*