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

Changed the name of a variable.

Browse source
This commit is contained in:
Harry Moffat 2008-04-18 22:21:47 +00:00
parent a10cf07683
commit 69e42f89cc
8 changed files with 228 additions and 229 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

336
Cantera/src/equil/vcs_TP.cpp
View file

@ -18,179 +18,177 @@
namespace VCSnonideal {
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg)
int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg)
/**************************************************************************
*
* vcs_TP:
*
* Solve an equilibrium problem at a particular fixed temperature
* and pressure
*
* ipr = 1 -> Print results to standard output
* 0 -> don't report on anything
* ip1 = 1 -> Print intermediate results.
* maxit -> Maximum number of iterations for the algorithm
* T = Temperature (Kelvin)
* pres = Pressure (units given by if__ variable)
*
* Return Codes
* ------------------
* 0 = Equilibrium Achieved
* 1 = Range space error encountered. The element abundance criteria are
* only partially satisfied. Specifically, the first NC= (number of
* components) conditions are satisfied. However, the full NE
* (number of elements) conditions are not satisfied. The equilibrirum
* condition is returned.
* -1 = Maximum number of iterations is exceeded. Convergence was not
* found.
***************************************************************************/
{
int retn, iconv;
/*
* Store the temperature and pressure in the private global variables
*/
T = T_arg;
Pres = pres_arg;
/*
* Evaluate the standard state free energies
* at the current temperatures and pressures.
*/
iconv = vcs_evalSS_TP(ipr, ip1, T, pres_arg);
/*
* Prepare the problem data:
* ->nondimensionalize the free energies using
* the divisor, R * T
*/
vcs_nondim_TP();
/*
* Prep the fe field
*/
vcs_fePrep_TP();
/*
* Decide whether we need an initial estimate of the solution
* If so, go get one. If not, then
*/
if (iest == -1) {
retn = vcs_inest_TP();
if (retn != VCS_SUCCESS) {
plogf("vcs_inest_TP returned a failure flag\n");
}
}
/*
* Solve the problem at a fixed Temperature and Pressure
* (all information concerning Temperature and Pressure has already
* been derived. The free energies are now in dimensionless form.)
*/
iconv = vcs_solve_TP(ipr, ip1, maxit);
/*
* Redimensionalize the free energies using
* the reverse of vcs_nondim to add back units.
*/
vcs_redim_TP();
/*
* Return the convergence success flag.
*/
return iconv;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*ARGSUSED*/
int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
/**************************************************************************
*
* vcs_evalSS_TP:
*
* IPR = 1 -> Print results to standard output
* 0 -> don't report on anything
* IP1 = 1 -> Print intermediate results.
* T = Temperature (Kelvin)
* Pres = Pressure (units of if__ variable)
*
* Evaluate the standard state free energies at the current temperature
* and pressure. Ideal gas pressure contribution is added in here.
*
***************************************************************************/
{
// int i;
//double R;
/*
* At this level of the program, we are still using values
* for the free energies that have units.
*/
// R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
/*
* We need to special case VCS_UNITS_UNITLESS, here.
* cpc_ts_GStar_calc() returns units of Kelvin. Also, the temperature
* comes into play in calculating the ideal equation of state
* contributions, and other equations of state also. Therefore,
* we will emulate the VCS_UNITS_KELVIN case, here by chaning
* the initial gibbs free energy units to Kelvin before feeding
* them to the cpc_ts_GStar_calc() routine. Then, we will revert
* them back to unitless at the end of this routine.
*/
/*
* Loop over the species calculating the standard state Gibbs free
* energies. -> These are energies that only depend upon the Temperature
* and possibly on the pressure (i.e., ideal gas, etc).
*/
// HKM -> We can change this to looks over phases, calling the vcs_VolPhase
// object. Working to get rid of VCS_SPECIES_THERMO object
//for (i = 0; i < m_numSpeciesTot; ++i) {
// VCS_SPECIES_THERMO *spt = SpeciesThermo[i];
// ff[i] = R * spt->GStar_R_calc(i, Temp, pres);
//}
for (int iph = 0; iph < NPhase; iph++) {
vcs_VolPhase* vph = VPhaseList[iph];
vph->setState_TP(T, Pres);
vph->sendToVCSGStar(VCS_DATA_PTR(m_SSfeSpecies));
}
if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
for (int i = 0; i < m_numSpeciesTot; ++i) {
m_SSfeSpecies[i] /= Temp;
}
}
return VCS_SUCCESS;
} /***************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
void VCS_SOLVE::vcs_fePrep_TP(void)
/**************************************************************************
*
*
***************************************************************************/
{
int i;
for (i = 0; i < m_numSpeciesTot; ++i) {
/**************************************************************************
*
* vcs_TP:
*
* Solve an equilibrium problem at a particular fixed temperature
* and pressure
*
* ipr = 1 -> Print results to standard output
* 0 -> don't report on anything
* ip1 = 1 -> Print intermediate results.
* maxit -> Maximum number of iterations for the algorithm
* T = Temperature (Kelvin)
* pres = Pressure (units given by if__ variable)
*
* Return Codes
* ------------------
* 0 = Equilibrium Achieved
* 1 = Range space error encountered. The element abundance criteria are
* only partially satisfied. Specifically, the first NC= (number of
* components) conditions are satisfied. However, the full NE
* (number of elements) conditions are not satisfied. The equilibrirum
* condition is returned.
* -1 = Maximum number of iterations is exceeded. Convergence was not
* found.
***************************************************************************/
{
int retn, iconv;
/*
* For single species phases, initialize the chemical
* potential with the value of the standard state chemical
* potential. This value doesn't change during the calculation
* Store the temperature and pressure in the private global variables
*/
if (SSPhase[i]) {
m_gibbsSpecies[i] = m_SSfeSpecies[i];
T = T_arg;
Pres = pres_arg;
/*
* Evaluate the standard state free energies
* at the current temperatures and pressures.
*/
iconv = vcs_evalSS_TP(ipr, ip1, T, pres_arg);
/*
* Prepare the problem data:
* ->nondimensionalize the free energies using
* the divisor, R * T
*/
vcs_nondim_TP();
/*
* Prep the fe field
*/
vcs_fePrep_TP();
/*
* Decide whether we need an initial estimate of the solution
* If so, go get one. If not, then
*/
if (iest == -1) {
retn = vcs_inest_TP();
if (retn != VCS_SUCCESS) {
plogf("vcs_inest_TP returned a failure flag\n");
}
}
}
} /* vcs_fePrep_TP() ********************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*
* Solve the problem at a fixed Temperature and Pressure
* (all information concerning Temperature and Pressure has already
* been derived. The free energies are now in dimensionless form.)
*/
iconv = vcs_solve_TP(ipr, ip1, maxit);
/*
* Redimensionalize the free energies using
* the reverse of vcs_nondim to add back units.
*/
vcs_redim_TP();
/*
* Return the convergence success flag.
*/
return iconv;
}
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*ARGSUSED*/
int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
/**************************************************************************
*
* vcs_evalSS_TP:
*
* IPR = 1 -> Print results to standard output
* 0 -> don't report on anything
* IP1 = 1 -> Print intermediate results.
* T = Temperature (Kelvin)
* Pres = Pressure (units of if__ variable)
*
* Evaluate the standard state free energies at the current temperature
* and pressure. Ideal gas pressure contribution is added in here.
*
***************************************************************************/
{
// int i;
//double R;
/*
* At this level of the program, we are still using values
* for the free energies that have units.
*/
// R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
/*
* We need to special case VCS_UNITS_UNITLESS, here.
* cpc_ts_GStar_calc() returns units of Kelvin. Also, the temperature
* comes into play in calculating the ideal equation of state
* contributions, and other equations of state also. Therefore,
* we will emulate the VCS_UNITS_KELVIN case, here by chaning
* the initial gibbs free energy units to Kelvin before feeding
* them to the cpc_ts_GStar_calc() routine. Then, we will revert
* them back to unitless at the end of this routine.
*/
/*
* Loop over the species calculating the standard state Gibbs free
* energies. -> These are energies that only depend upon the Temperature
* and possibly on the pressure (i.e., ideal gas, etc).
*/
// HKM -> We can change this to looks over phases, calling the vcs_VolPhase
// object. Working to get rid of VCS_SPECIES_THERMO object
//for (i = 0; i < m_numSpeciesTot; ++i) {
// VCS_SPECIES_THERMO *spt = SpeciesThermo[i];
// ff[i] = R * spt->GStar_R_calc(i, Temp, pres);
//}
for (int iph = 0; iph < NPhase; iph++) {
vcs_VolPhase* vph = VPhaseList[iph];
vph->setState_TP(T, Pres);
vph->sendToVCSGStar(VCS_DATA_PTR(m_SSfeSpecies));
}
if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) {
for (int i = 0; i < m_numSpeciesTot; ++i) {
m_SSfeSpecies[i] /= Temp;
}
}
return VCS_SUCCESS;
} /***************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/
void VCS_SOLVE::vcs_fePrep_TP(void)
/**************************************************************************
*
*
***************************************************************************/
{
int i;
for (i = 0; i < m_numSpeciesTot; ++i) {
/*
* For single species phases, initialize the chemical
* potential with the value of the standard state chemical
* potential. This value doesn't change during the calculation
*/
if (SSPhase[i]) {
m_feSpecies_curr[i] = m_SSfeSpecies[i];
}
}
} /* vcs_fePrep_TP() ********************************************************/
}

12
Cantera/src/equil/vcs_inest.cpp
View file

@ -195,14 +195,14 @@ namespace VCSnonideal {
wt[kspec] = 0.0;
}
}
vcs_dcopy(VCS_DATA_PTR(m_gibbsSpecies), VCS_DATA_PTR(m_SSfeSpecies),
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(m_SSfeSpecies),
nspecies);
for (kspec = 0; kspec < m_numComponents; ++kspec) {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
if (! SSPhase[kspec]) {
iph = PhaseID[kspec];
m_gibbsSpecies[kspec] += log(wt[kspec] / TPhMoles[iph]);
m_feSpecies_curr[kspec] += log(wt[kspec] / TPhMoles[iph]);
}
} else {
wt[kspec] = 0.0;
@ -214,11 +214,11 @@ namespace VCSnonideal {
for (kspec = 0; kspec < nspecies; ++kspec) {
plogf("%s", pprefix); plogf("%-12.12s", SpName[kspec].c_str());
if (kspec < m_numComponents)
plogf("fe* = %15.5g ff = %15.5g\n", m_gibbsSpecies[kspec],
plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_curr[kspec],
m_SSfeSpecies[kspec]);
else
plogf("fe* = %15.5g ff = %15.5g dg* = %15.5g\n",
m_gibbsSpecies[kspec], m_SSfeSpecies[kspec], dg[kspec-m_numComponents]);
m_feSpecies_curr[kspec], m_SSfeSpecies[kspec], dg[kspec-m_numComponents]);
}
}
#endif
@ -324,7 +324,7 @@ namespace VCSnonideal {
/* ******************************************* */
vcs_dfe(molNum, 0, 0, 0, nspecies);
for (kspec = 0, s = 0.0; kspec < nspecies; ++kspec) {
s += ds[kspec] * m_gibbsSpecies[kspec];
s += ds[kspec] * m_feSpecies_curr[kspec];
}
if (s == 0.0) {
finished = TRUE; continue;
@ -502,7 +502,7 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix,
vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_gibbsSpecies),
vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles)));
plogendl();
}

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

@ -103,7 +103,7 @@ void VCS_SOLVE::vcs_nondim_TP(void) {
* Thus, we may divide it by the temperature.
*/
m_SSfeSpecies[i] *= tf;
m_gibbsSpecies[i] *= tf;
m_feSpecies_curr[i] *= tf;
dg[i] *= tf;
dgl[i] *= tf;
m_feSpecies_old[i] *= tf;
@ -143,7 +143,7 @@ void VCS_SOLVE::vcs_redim_TP(void)
* FF(I), to make it have units, i.e. mu = RT * mu_star
*/
m_SSfeSpecies[i] *= tf;
m_gibbsSpecies[i] *= tf;
m_feSpecies_curr[i] *= tf;
dg[i] *= tf;
dgl[i] *= tf;
m_feSpecies_old[i] *= tf;

2
Cantera/src/equil/vcs_prep.cpp
View file

@ -284,7 +284,7 @@ int VCS_SOLVE::vcs_prep(void) {
/*
* Initialize various arrays in the data to zero
*/
vcs_dzero(VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesTot);
vcs_dzero(VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesTot);
vcs_vdzero(m_feSpecies_old, m_numSpeciesTot);
vcs_vdzero(wt, m_numSpeciesTot);
vcs_dzero(&(DnPhase[0][0]), m_numSpeciesTot*NPhase);

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

@ -122,7 +122,7 @@ int VCS_SOLVE::vcs_report(int iconv)
for (i = 0; i < m_numComponents; ++i) {
plogf(" %-12.12s", SpName[i].c_str());
print_space(13);
plogf("%14.7E %14.7E %12.4E", soln[i], wt[i], m_gibbsSpecies[i]);
plogf("%14.7E %14.7E %12.4E", soln[i], wt[i], m_feSpecies_curr[i]);
plogf(" %3d", SpeciesUnknownType[i]);
plogf("\n");
}
@ -132,10 +132,10 @@ int VCS_SOLVE::vcs_report(int iconv)
print_space(13);
if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_MOLNUM) {
plogf("%14.7E %14.7E %12.4E", soln[l], wt[l], m_gibbsSpecies[l]);
plogf("%14.7E %14.7E %12.4E", soln[l], wt[l], m_feSpecies_curr[l]);
plogf(" MolNum ");
} else if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" NA %14.7E %12.4E", 1.0, m_gibbsSpecies[l]);
plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_curr[l]);
plogf(" Voltage = %14.7E", soln[l]);
} else {
plogf("we have a problem\n");
@ -256,7 +256,7 @@ int VCS_SOLVE::vcs_report(int iconv)
gaTPhase[j] += gaPhase[j];
}
gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(soln),
VCS_DATA_PTR(m_gibbsSpecies));
VCS_DATA_PTR(m_feSpecies_curr));
gibbsTotal += gibbsPhase;
plogf(" | %18.11E |\n", gibbsPhase);
}
@ -280,7 +280,7 @@ int VCS_SOLVE::vcs_report(int iconv)
* energy of zero
*/
g = vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_gibbsSpecies),
g = vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles));
plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
if (inertYes)
@ -325,17 +325,17 @@ int VCS_SOLVE::vcs_report(int iconv)
if (tpmoles > 0.0 && soln[l] > 0.0) {
lx = log(soln[l]) - log(tpmoles);
} else {
lx = m_gibbsSpecies[l] - m_SSfeSpecies[l] - log(ActCoeff[l]) + SpecLnMnaught[l];
lx = m_feSpecies_curr[l] - m_SSfeSpecies[l] - log(ActCoeff[l]) + SpecLnMnaught[l];
}
}
plogf("%14.7E |", lx);
plogf("%14.7E | ", eContrib);
double tmp = m_SSfeSpecies[l] + log(ActCoeff[l]) + lx - SpecLnMnaught[l] + eContrib;
if (fabs(m_gibbsSpecies[l] - tmp) > 1.0E-8) {
if (fabs(m_feSpecies_curr[l] - tmp) > 1.0E-8) {
plogf("\n\t\twe have a problem - doesn't add up\n");
exit(-1);
}
plogf(" %12.4E |", m_gibbsSpecies[l]);
plogf(" %12.4E |", m_feSpecies_curr[l]);
if (SpecLnMnaught[l] != 0.0) {
plogf(" (%14.7E)", - SpecLnMnaught[l]);
}

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

@ -109,7 +109,7 @@ namespace VCSnonideal {
scSize.resize(nspecies0, 0.0);
m_spSize.resize(nspecies0, 1.0);
m_gibbsSpecies.resize(nspecies0, 0.0);
m_feSpecies_curr.resize(nspecies0, 0.0);
m_SSfeSpecies.resize(nspecies0, 0.0);
m_feSpecies_new.resize(nspecies0, 0.0);
soln.resize(nspecies0, 0.0);
@ -789,7 +789,7 @@ namespace VCSnonideal {
k = ind[kspec];
soln[kspec] = pub->w[k];
wt[kspec] = pub->mf[k];
m_gibbsSpecies[kspec] = pub->m_gibbsSpecies[k];
m_feSpecies_curr[kspec] = pub->m_gibbsSpecies[k];
}
/*
@ -925,7 +925,7 @@ namespace VCSnonideal {
plogf("voltage species = %g\n", soln[k1]);
}
pub->mf[i] = wt[k1];
pub->m_gibbsSpecies[i] = m_gibbsSpecies[k1];
pub->m_gibbsSpecies[i] = m_feSpecies_curr[k1];
pub->VolPM[i] = VolPM[k1];
}

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

@ -430,8 +430,9 @@ public:
* The first NC entries are for components. The following
* NR entries are for the current non-component species in the mechanism.
* The dimension of this vector is specified by the m_VCS_UnitsFormat variable.
* Length = number of species
*/
std::vector<double> m_gibbsSpecies;
std::vector<double> m_feSpecies_curr;
//! Free energy vector from the start of the current iteration
/*!

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

@ -442,8 +442,8 @@ namespace VCSnonideal {
}
}
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesRdc);
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(dgl), VCS_DATA_PTR(dg), m_numRxnRdc);
@ -903,7 +903,7 @@ namespace VCSnonideal {
* which have yet to be processed in the main loop
*/
for (ll = kspec+1; ll < m_numSpeciesRdc; ++ll) {
m_feSpecies_old[ll] = m_gibbsSpecies[ll];
m_feSpecies_old[ll] = m_feSpecies_curr[ll];
}
for (ll = irxn+1; ll < m_numRxnRdc; ++ll) {
dgl[ll] = dg[ll];
@ -1120,7 +1120,7 @@ namespace VCSnonideal {
vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_old),
VCS_DATA_PTR(TPhMoles)));
plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E",
vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_gibbsSpecies),
vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles1)));
plogendl();
}
@ -1158,7 +1158,7 @@ namespace VCSnonideal {
plogf(" Total moles of liquid = %15.7E\n", 0.0);
}
plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n",
vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_gibbsSpecies),
vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles1)));
plogf(" -----------------------------------------------------");
plogendl();
@ -1181,14 +1181,14 @@ namespace VCSnonideal {
for (i = 0; i < m_numComponents; ++i) {
plogf(" --- %-12.12s", SpName[i].c_str()); plogf(" ");
plogf("%14.6E%14.6E%14.6E%14.6E\n", soln[i],
wt[i], m_feSpecies_old[i], m_gibbsSpecies[i]);
wt[i], m_feSpecies_old[i], m_feSpecies_curr[i]);
}
for (i = m_numComponents; i < m_numSpeciesRdc; ++i) {
l1 = i - m_numComponents;
plogf(" --- %-12.12s", SpName[i].c_str());
plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
spStatus[l1], soln[i],
wt[i], m_feSpecies_old[i], m_gibbsSpecies[i],
wt[i], m_feSpecies_old[i], m_feSpecies_curr[i],
dgl[l1], dg[l1]);
}
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
@ -1196,7 +1196,7 @@ namespace VCSnonideal {
plogf(" --- %-12.12s", SpName[kspec].c_str());
plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
spStatus[l1], soln[kspec],
wt[kspec], m_feSpecies_old[kspec], m_gibbsSpecies[kspec],
wt[kspec], m_feSpecies_old[kspec], m_feSpecies_curr[kspec],
dgl[l1], dg[l1]);
}
plogf(" ---"); print_space(56);
@ -1216,7 +1216,7 @@ namespace VCSnonideal {
vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_old),
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_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles1)));
plogendl();
if (m_VCount->Its > 550) {
@ -1244,7 +1244,7 @@ namespace VCSnonideal {
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(m_feSpecies_old), VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc);
vcs_updateVP(0);
/*
@ -2260,7 +2260,7 @@ namespace VCSnonideal {
spStatus[irxn] = VCS_SPECIES_DELETED;
dg[irxn] = 0.0;
dgl[irxn] = 0.0;
m_gibbsSpecies[kspec] = 0.0;
m_feSpecies_curr[kspec] = 0.0;
m_feSpecies_old[kspec] = 0.0;
wt[kspec] = 0.0;
/*
@ -2525,7 +2525,7 @@ namespace VCSnonideal {
* for formation reactions
*/
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
/*
* Recalculate the DeltaG's of the formation reactions for the
@ -2613,7 +2613,7 @@ namespace VCSnonideal {
* ~ infinite dilution.
*/
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
/*
* Recalculate the DeltaG's of the formation reactions for the
@ -3155,10 +3155,10 @@ namespace VCSnonideal {
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
if (spStatus[irxn] != VCS_SPECIES_MINOR) {
icase = 0;
dg[irxn] = m_gibbsSpecies[ir[irxn]];
dg[irxn] = m_feSpecies_curr[ir[irxn]];
dtmp_ptr = sc[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
dg[irxn] += dtmp_ptr[kspec] * m_gibbsSpecies[kspec];
dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
icase = 1;
}
@ -3174,10 +3174,10 @@ namespace VCSnonideal {
/* ************************************************* */
for (irxn = 0; irxn < irxnl; ++irxn) {
icase = 0;
dg[irxn] = m_gibbsSpecies[ir[irxn]];
dg[irxn] = m_feSpecies_curr[ir[irxn]];
dtmp_ptr = sc[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
dg[irxn] += dtmp_ptr[kspec] * m_gibbsSpecies[kspec];
dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
icase = 1;
}
@ -3193,10 +3193,10 @@ namespace VCSnonideal {
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
if (spStatus[irxn] <= VCS_SPECIES_MINOR) {
icase = 0;
dg[irxn] = m_gibbsSpecies[ir[irxn]];
dg[irxn] = m_feSpecies_curr[ir[irxn]];
dtmp_ptr = sc[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
dg[irxn] += dtmp_ptr[kspec] * m_gibbsSpecies[kspec];
dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
icase = 1;
}
@ -4328,24 +4328,24 @@ namespace VCSnonideal {
exit(-1);
}
#endif
m_gibbsSpecies[kspec] =
m_feSpecies_curr[kspec] =
m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase];
} else {
if (SSPhase[kspec]) {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
} else {
if (z[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
iph = PhaseID[kspec];
if (tPhMoles_ptr[iph] > 0.0) {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]
+ log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ Charge[kspec] * Faraday_dim * phasePhi[iphase];
} else {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ Charge[kspec] * Faraday_dim * phasePhi[iphase];
}
@ -4371,24 +4371,24 @@ namespace VCSnonideal {
exit(-1);
}
#endif
m_gibbsSpecies[kspec] =
m_feSpecies_curr[kspec] =
m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase];
} else {
if (SSPhase[kspec]) {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
} else {
if (z[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
iph = PhaseID[kspec];
if (tPhMoles_ptr[iph] > 0.0) {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]
+ log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ Charge[kspec] * Faraday_dim * phasePhi[iphase]; ;
} else {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]
+ Charge[kspec] * Faraday_dim * phasePhi[iphase];
}
@ -4415,24 +4415,24 @@ namespace VCSnonideal {
exit(-1);
}
#endif
m_gibbsSpecies[kspec] =
m_feSpecies_curr[kspec] =
m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase]; ;
} else {
if (SSPhase[kspec]) {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
} else {
if (z[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) {
iph = PhaseID[kspec];
if (tPhMoles_ptr[iph] > 0.0) {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]
+ log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF)
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec];
} else {
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec];
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec];
}
} else {
st_ptr = SpeciesThermo[kspec];
m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec])
- tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec];
}
}
@ -4674,7 +4674,7 @@ namespace VCSnonideal {
SWAP(wt[k1], wt[k2], t1);
SWAP(m_SSfeSpecies[k1], m_SSfeSpecies[k2], t1);
SWAP(m_spSize[k1], m_spSize[k2], t1);
SWAP(m_gibbsSpecies[k1], m_gibbsSpecies[k2], t1);
SWAP(m_feSpecies_curr[k1], m_feSpecies_curr[k2], t1);
SWAP(ds[k1], ds[k2], t1);
SWAP(m_feSpecies_old[k1], m_feSpecies_old[k2], t1);
SWAP(m_feSpecies_new[k1], m_feSpecies_new[k2], t1);
@ -4783,10 +4783,10 @@ namespace VCSnonideal {
#endif
if (kspec >= m_numComponents) {
irxn = kspec - m_numComponents;
dg[irxn] = m_gibbsSpecies[kspec];
dg[irxn] = m_feSpecies_curr[kspec];
dtmp_ptr = sc[irxn];
for (kcomp = 0; kcomp < m_numComponents; ++kcomp) {
dg[irxn] += dtmp_ptr[kcomp] * m_gibbsSpecies[kcomp];
dg[irxn] += dtmp_ptr[kcomp] * m_feSpecies_curr[kcomp];
}
}
}
@ -4802,10 +4802,10 @@ namespace VCSnonideal {
iph = PhaseID[kspec];
if (iph == iphase ) {
if (soln[kspec] > 0.0) zeroedPhase = FALSE;
dg[irxn] = m_gibbsSpecies[kspec];
dg[irxn] = m_feSpecies_curr[kspec];
dtmp_ptr = sc[irxn];
for (kcomp = 0; kcomp < m_numComponents; ++kcomp) {
dg[irxn] += dtmp_ptr[kcomp] * m_gibbsSpecies[kcomp];
dg[irxn] += dtmp_ptr[kcomp] * m_feSpecies_curr[kcomp];
}
}
}