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

Name changes on variables, only

Browse source
This commit is contained in:
Harry Moffat 2008-05-17 23:52:13 +00:00
parent 431eb16260
commit 7a99ca115f
6 changed files with 165 additions and 86 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -106,12 +106,12 @@ namespace VCSnonideal {
* Calculate some quantities that may need updating
*/
vcs_tmoles();
Vol = vcs_VolTotal(m_temperature, m_pressurePA,
VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM));
m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA,
VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_PMVolumeSpecies));
plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
plogf("\t\tVolume = %15.5g m**3\n", Vol);
plogf("\t\tVolume = %15.5g m**3\n", m_totalVol);
/*
* -------- TABLE OF SPECIES IN DECREASING MOLE NUMBERS --------------

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

@ -285,7 +285,7 @@ namespace VCSnonideal {
/*
* First the diagonal term of the Jacobian
*/
s = dLnActCoeffdMolNum[kspec][kspec];
s = m_dLnActCoeffdMolNum[kspec][kspec];
/*
* Next, the other terms. Note this only a loop over the components
* So, it's not too expensive to calculate.
@ -294,11 +294,11 @@ namespace VCSnonideal {
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] * dLnActCoeffdMolNum[k][l];
s += sc_irxn[k] * sc_irxn[l] * m_dLnActCoeffdMolNum[k][l];
}
}
if (kph == m_phaseID[l]) {
s += sc_irxn[l] * (dLnActCoeffdMolNum[kspec][l] + dLnActCoeffdMolNum[l][kspec]);
s += sc_irxn[l] * (m_dLnActCoeffdMolNum[kspec][l] + m_dLnActCoeffdMolNum[l][kspec]);
}
}
}
@ -334,7 +334,7 @@ namespace VCSnonideal {
* -> This scatter calculation is carried out in the
* vcs_VolPhase object.
*/
Vphase->sendToVCSLnActCoeffJac(dLnActCoeffdMolNum.baseDataAddr());
Vphase->sendToVCSLnActCoeffJac(m_dLnActCoeffdMolNum.baseDataAddr());
}
}
}

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

@ -50,8 +50,8 @@ namespace VCSnonideal {
m_tolmaj2(0.0),
m_tolmin2(0.0),
UnitsState(VCS_DIMENSIONAL_G),
UseActCoeffJac(0),
Vol(0.0),
m_useActCoeffJac(0),
m_totalVol(0.0),
Faraday_dim(1.602e-19 * 6.022136736e26),
m_VCount(0),
vcs_debug_print_lvl(0),
@ -186,7 +186,7 @@ namespace VCSnonideal {
CurrPhAC.resize(nphase0, 0);
m_wtSpecies.resize(nspecies0, 0.0);
m_chargeSpecies.resize(nspecies0, 0.0);
SpeciesThermo.resize(nspecies0, (VCS_SPECIES_THERMO *)0);
m_speciesThermoList.resize(nspecies0, (VCS_SPECIES_THERMO *)0);
/*
* Malloc Phase Info
@ -199,12 +199,12 @@ namespace VCSnonideal {
/*
* For Future expansion
*/
UseActCoeffJac = true;
if (UseActCoeffJac ) {
dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0);
m_useActCoeffJac = true;
if (m_useActCoeffJac) {
m_dLnActCoeffdMolNum.resize(nspecies0, nspecies0, 0.0);
}
VolPM.resize(nspecies0, 0.0);
m_PMVolumeSpecies.resize(nspecies0, 0.0);
/*
* Malloc space for counters kept within vcs
@ -241,8 +241,8 @@ namespace VCSnonideal {
}
for (j = 0; j < nspecies; j++) {
delete SpeciesThermo[j];
SpeciesThermo[j] = 0;
delete m_speciesThermoList[j];
m_speciesThermoList[j] = 0;
}
delete m_VCount; m_VCount = 0;
@ -516,12 +516,12 @@ namespace VCSnonideal {
*
*/
for (kspec = 0; kspec < nspecies; kspec++) {
if (SpeciesThermo[kspec] != NULL) {
delete SpeciesThermo[kspec];
if (m_speciesThermoList[kspec] != NULL) {
delete m_speciesThermoList[kspec];
}
VCS_SPECIES_THERMO *spf = pub->SpeciesThermo[kspec];
SpeciesThermo[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO();
if (SpeciesThermo[kspec] == NULL) {
m_speciesThermoList[kspec] = spf->duplMyselfAsVCS_SPECIES_THERMO();
if (m_speciesThermoList[kspec] == NULL) {
plogf(" duplMyselfAsVCS_SPECIES_THERMO returned an error!\n");
return VCS_PUB_BAD;
}
@ -705,7 +705,7 @@ namespace VCSnonideal {
for (int k = 0; k < Vphase->NVolSpecies; k++) {
vcs_SpeciesProperties *sProp = Vphase->ListSpeciesPtr[k];
int kT = Vphase->IndSpecies[k];
sProp->SpeciesThermo = SpeciesThermo[kT];
sProp->SpeciesThermo = m_speciesThermoList[kT];
}
}
@ -747,9 +747,9 @@ namespace VCSnonideal {
/*
* Copy the volume info
*/
Vol = pub->Vol;
if (VolPM.size() != 0) {
vcs_dcopy(VCS_DATA_PTR(VolPM), VCS_DATA_PTR(pub->VolPM), nspecies);
m_totalVol = pub->Vol;
if (m_PMVolumeSpecies.size() != 0) {
vcs_dcopy(VCS_DATA_PTR(m_PMVolumeSpecies), VCS_DATA_PTR(pub->VolPM), nspecies);
}
/*
@ -767,7 +767,7 @@ namespace VCSnonideal {
*
* It's assumed we are solving the same problem.
*
* @param pub Pointer to VCS_PROB that will be used to
* @param pub Pointer to VCS_PROdB that will be used to
* initialize the current equilibrium problem
*/
int VCS_SOLVE::vcs_prob_specify(const VCS_PROB *pub) {
@ -781,7 +781,7 @@ namespace VCSnonideal {
m_VCS_UnitsFormat = pub->m_VCS_UnitsFormat;
m_doEstimateEquil = pub->iest;
Vol = pub->Vol;
m_totalVol = pub->Vol;
m_tolmaj = pub->tolmaj;
m_tolmin = pub->tolmin;
@ -906,8 +906,8 @@ namespace VCSnonideal {
int k1 = 0;
vcs_tmoles();
Vol = vcs_VolTotal(m_temperature, m_pressurePA,
VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM));
m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA,
VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_PMVolumeSpecies));
for (i = 0; i < m_numSpeciesTot; ++i) {
/*
@ -930,12 +930,12 @@ namespace VCSnonideal {
}
pub->mf[i] = m_molNumSpecies_new[k1];
pub->m_gibbsSpecies[i] = m_feSpecies_curr[k1];
pub->VolPM[i] = VolPM[k1];
pub->VolPM[i] = m_PMVolumeSpecies[k1];
}
pub->T = m_temperature;
pub->PresPA = m_pressurePA;
pub->Vol = Vol;
pub->Vol = m_totalVol;
int kT = 0;
for (int iph = 0; iph < pub->NPhase; iph++) {
vcs_VolPhase *pubPhase = pub->VPhaseList[iph];

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

@ -724,7 +724,19 @@ private:
int delete_species(int kspec);
void delete_multiphase(int iph);
int delta_species(int kspec, double *delta_ptr);
void add_deleted(void);
//! Provide an estimate for the deleted species in phases that
//! are not zeroed out
/*!
* Try to add back in all deleted species. An estimate of the kmol numbers
* are obtained and the species is added back into the equation system,
* into the old state vector.
*
* This routine is called at the end of the calculation, just before
* returning to the user.
*/
int vcs_add_all_deleted();
int recheck_deleted(void);
//! Alternative treatment for the update of a minor species
@ -1311,10 +1323,10 @@ public:
* needed. It is not swapped wrt species
* (unused atm)
*/
DoubleStarStar dLnActCoeffdMolNum;
DoubleStarStar m_dLnActCoeffdMolNum;
//! This boolean indicates whether the activity coefficients for a phase
//! are current.
//! This boolean indicates whether the activity coefficients for a phase
//! are current.
std::vector<int> CurrPhAC;
//! Molecular weight of each species
@ -1338,7 +1350,7 @@ public:
/*!
* SpeciesThermo[k] pointer to the thermo information for the kth species
*/
std::vector<VCS_SPECIES_THERMO *> SpeciesThermo;
std::vector<VCS_SPECIES_THERMO *> m_speciesThermoList;
//! Choice of Hessians
/*!
@ -1346,16 +1358,20 @@ public:
* Hessian based on Jacobian of the ln(ActCoeff) with respect to mole
* numbers
*/
int UseActCoeffJac;
int m_useActCoeffJac;
double Vol; /* Vol = Volume (m^3) */
//! Total volume of all phases
/*!
* units are m^3
*/
double m_totalVol;
//! Partialm molar volumes of the species
//! Partial molar volumes of the species
/*!
* units = mks (m^3/kmol) -determined by m_VCS_UnitsFormat
* Length = number of species
*/
std::vector<double> VolPM;
std::vector<double> m_PMVolumeSpecies;
//! dimensionless value of Faraday's constant
/*!

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

@ -1867,7 +1867,19 @@ namespace VCSnonideal {
* Add back deleted species in non-zeroed phases. Estimate their
* mole numbers.
*/
add_deleted();
npb = vcs_add_all_deleted();
if (npb > 0) {
MajorSpeciesHaveConverged = true;
iti = 0;
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 1) {
plogf(" --- add_all_deleted(): some rxns not converged. RETURNING TO LOOP!");
plogendl();
}
#endif
goto L_MAINLOOP_ALL_SPECIES;
}
/*
* Make sure the volume phase objects hold the same state and
* information as the vcs object. This also update the Cantera objects
@ -2083,9 +2095,6 @@ namespace VCSnonideal {
}
return dx;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
int VCS_SOLVE::delta_species(int kspec, double *delta_ptr)
@ -2284,9 +2293,10 @@ namespace VCSnonideal {
* have to signal the calling code
*/
return (m_numRxnRdc == 0);
} /* delete_species() ********************************************************/
}
/***************************************************************************/
/****************************************************************************
/*
*
* reinsert_deleted():
*
@ -2560,33 +2570,32 @@ namespace VCSnonideal {
}
}
return npb;
} /* recheck_deleted() *******************************************************/
}
/***********************************************************************************/
void VCS_SOLVE::add_deleted(void)
/*************************************************************************
*
* Provide an estimate for the deleted species in phases that
* are not zeroed out
*
*************************************************************************/
{
// Provide an estimate for the deleted species in phases that
// are not zeroed out
/*
* Try to add back in all deleted species. An estimate of the kmol numbers
* are obtained and the species is added back into the equation system,
* into the old state vector.
*/
int VCS_SOLVE::vcs_add_all_deleted() {
int iph, kspec, retn;
if (m_numSpeciesRdc == m_numSpeciesTot) return;
if (m_numSpeciesRdc == m_numSpeciesTot) return 0;
/*
* Use the standard chemical potentials for the chemical potentials
* of deleted species. Then, calculate Delta G for
* for formation reactions
*
* HKM Note: We need to update this step for nonunity activity
* coefficients.
* The formula will be fe = ff + RT * ln(actCoeff)
* where the activity coefficient is evaluated at
* ~ infinite dilution.
* for formation reactions.
* We are relying here on a old saved value of m_actCoeffSpecies_old[kspec]
* being sufficiently good. Note, we will recalculate everything at the
* end of the routine.
*/
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
iph = m_phaseID[kspec];
m_feSpecies_curr[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec])
- SpecLnMnaught[kspec]
+ m_chargeSpecies[kspec] * Faraday_dim * m_phasePhi[iph]);
}
/*
* Recalculate the DeltaG's of the formation reactions for the
@ -2594,20 +2603,73 @@ namespace VCSnonideal {
*/
vcs_deltag(0, true);
for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] > 0.0) {
double maxDG = MIN(m_deltaGRxn_new[irxn], 300);
double maxDG = MIN(m_deltaGRxn_new[irxn], 300.0);
double dx = m_tPhaseMoles_old[iph] * exp(- maxDG);
retn = delta_species(kspec, &dx);
if (retn == 0) {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl) {
plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
m_speciesName[kspec].c_str(), kspec, dx);
}
#endif
if (dx > 1.0E-50) {
dx = 1.0E-50;
retn = delta_species(kspec, &dx);
#ifdef DEBUG_MODE
if (retn == 0) {
if (vcs_debug_print_lvl) {
plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
m_speciesName[kspec].c_str(), kspec, dx);
}
}
#endif
}
}
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
if (retn != 0) {
plogf(" --- add_deleted(): species %s added back in with mol number %g",
m_speciesName[kspec].c_str(), dx);
plogendl();
} else {
plogf(" --- add_deleted(): species %s failed to be added back in");
plogendl();
}
}
#endif
}
}
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot);
vcs_deltag(0, true);
retn = 0;
for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
kspec = m_indexRxnToSpecies[irxn];
iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] > 0.0) {
if (fabs(m_deltaGRxn_old[irxn]) > m_tolmin) {
retn++;
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g",
m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec],
m_deltaGRxn_old[irxn]);
plogendl();
}
#endif
}
}
}
return retn;
}
/***********************************************************************************/
/* globalStepDamp
*
@ -2806,7 +2868,7 @@ namespace VCSnonideal {
* We update the matrix dlnActCoeffdmolNumber[][] at the
* top of the loop, when necessary
*/
if (UseActCoeffJac) {
if (m_useActCoeffJac) {
vcs_CalcLnActCoeffJac(VCS_DATA_PTR(m_molNumSpecies_old));
}
/************************************************************************
@ -2855,7 +2917,8 @@ namespace VCSnonideal {
#endif
Vphase = m_VolPhaseList[iph];
int numSpPhase = Vphase->NVolSpecies;
m_deltaMolNumSpecies[kspec] = m_totalMolNum * 10.0 * VCS_DELETE_PHASE_CUTOFF / numSpPhase;
m_deltaMolNumSpecies[kspec] =
m_totalMolNum * 10.0 * VCS_DELETE_PHASE_CUTOFF / numSpPhase;
}
--(m_numRxnMinorZeroed);
} else {
@ -2881,7 +2944,8 @@ namespace VCSnonideal {
if (vcs_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;
@ -2897,7 +2961,8 @@ namespace VCSnonideal {
if (vcs_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;
@ -2913,7 +2978,7 @@ namespace VCSnonideal {
for (j = 0; j < m_numComponents; ++j) {
if (!m_SSPhase[j]) {
if (m_molNumSpecies_old[j] > 0.0) {
s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
}
}
}
@ -2930,7 +2995,7 @@ namespace VCSnonideal {
* derivatives of the activity coefficients with respect to the
* mole numbers, even in our diagonal approximation.
*/
if (UseActCoeffJac) {
if (m_useActCoeffJac) {
double s_old = s;
s = vcs_Hessian_diag_adj(irxn, s_old);
#ifdef DEBUG_MODE
@ -3074,7 +3139,7 @@ namespace VCSnonideal {
}
/*****************************************************************************/
/**************************************************************************
/*
*
* vcs_deltag:
*
@ -3086,7 +3151,7 @@ namespace VCSnonideal {
* species I in this reaction.
*
* INPUT
* L = < 0 : Calculate reactions corresponding to
* L < 0 : Calculate reactions corresponding to
* major noncomponent and zeroed species only
* L = 0 : Do all noncomponent reactions, i, between
* 0 <= i < irxnl
@ -3185,7 +3250,7 @@ namespace VCSnonideal {
}
}
/* ************************************************* */
/* **** MULTISPECIES PHASES WITH ZERO MOLES************ */
/* **** MULTISPECIES PHASES WITH ZERO MOLES ******** */
/* ************************************************* */
/*
* Massage the free energies for species with zero mole fractions
@ -3273,9 +3338,7 @@ namespace VCSnonideal {
checkFinite(m_deltaGRxn_new[irxn]);
}
#endif
} /* vcs_deltag() ************************************************************/
/*****************************************************************************/
/*****************************************************************************/
}
/*****************************************************************************/
int VCS_SOLVE::vcs_basopt(int ifirst, double aw[], double sa[], double sm[],
@ -4450,7 +4513,7 @@ namespace VCSnonideal {
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
st_ptr = SpeciesThermo[kspec];
st_ptr = m_speciesThermoList[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(actCoeff_ptr[kspec] * z[kspec])
- tlogMoles[m_phaseID[kspec]] - SpecLnMnaught[kspec];
}
@ -4684,14 +4747,14 @@ namespace VCSnonideal {
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);
SWAP(VolPM[k1], VolPM[k2], t1);
SWAP(m_speciesThermoList[k1], m_speciesThermoList[k2], st_tmp);
SWAP(m_PMVolumeSpecies[k1], m_PMVolumeSpecies[k2], t1);
for (j = 0; j < m_numElemConstraints; ++j) {
SWAP(m_formulaMatrix[j][k1], m_formulaMatrix[j][k2], t1);
}
if (UseActCoeffJac) {
vcs_switch2D(dLnActCoeffdMolNum.baseDataAddr(), k1, k2);
if (m_useActCoeffJac) {
vcs_switch2D(m_dLnActCoeffdMolNum.baseDataAddr(), k1, k2);
}
/*

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

@ -389,7 +389,7 @@ double VCS_SOLVE::vcs_Gxs_phase_calc(vcs_VolPhase *Vphase, double *mf_PO)
if (Vphase->Activity_Coeff_Model != VCS_AC_CONSTANT) {
for (kspec = 0; kspec < Vphase->NVolSpecies; kspec++) {
kglob = Vphase->IndSpecies[kspec];
ts_ptr = SpeciesThermo[kglob];
ts_ptr = m_speciesThermoList[kglob];
ac = ts_ptr->eval_ac(kspec);
Gxs += mf_PO[kspec] * log(ac);
}
@ -434,7 +434,7 @@ double VCS_SOLVE::vcs_Gxs_calc(int iphase)
for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
if (m_phaseID[kspec] == iphase) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
ts_ptr = SpeciesThermo[kspec];
ts_ptr = m_speciesThermoList[kspec];
ac = ts_ptr->eval_ac(kspec);
Gxs += m_molNumSpecies_old[kspec]/totmol * log(ac);
} else {