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

[Equil] Make better use of local variables

Browse source
This commit is contained in:
Ray Speth 2016-04-03 17:40:38 -04:00
parent 60eed786fa
commit a60217cfc6
15 changed files with 412 additions and 600 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -110,7 +110,7 @@ typedef double(*VCS_FUNC_PTR)(double xval, double Vtarget,
* @param vec vector of doubles
* @return Returns the l2 norm of the vector
*/
double vcs_l2norm(const vector_fp vec);
double vcs_l2norm(const vector_fp& vec);
//! Finds the location of the maximum component in a double vector
/*!

148
src/equil/BasisOptimize.cpp
View file

@ -10,6 +10,7 @@ using namespace std;
namespace Cantera
{
int BasisOptimize_print_lvl = 0;
static const double USEDBEFORE = -1;
//! Print a string within a given space limit.
/*!
@ -28,40 +29,27 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
std::vector<size_t>& orderVectorElements,
vector_fp& formRxnMatrix)
{
size_t j, jj, k=0, kk, i, jl, ml;
std::string ename;
std::string sname;
// Get the total number of elements defined in the multiphase object
size_t ne = mphase->nElements();
// Get the total number of species in the multiphase object
size_t nspecies = mphase->nSpecies();
doublereal tmp;
doublereal const USEDBEFORE = -1;
// Perhaps, initialize the element ordering
if (orderVectorElements.size() < ne) {
orderVectorElements.resize(ne);
for (j = 0; j < ne; j++) {
orderVectorElements[j] = j;
}
iota(orderVectorElements.begin(), orderVectorElements.end(), 0);
}
// Perhaps, initialize the species ordering
if (orderVectorSpecies.size() != nspecies) {
orderVectorSpecies.resize(nspecies);
for (k = 0; k < nspecies; k++) {
orderVectorSpecies[k] = k;
}
iota(orderVectorSpecies.begin(), orderVectorSpecies.end(), 0);
}
if (BasisOptimize_print_lvl >= 1) {
writelog(" ");
for (i=0; i<77; i++) {
writelog("-");
}
writelog("\n");
writeline('-', 77);
writelog(" --- Subroutine BASOPT called to ");
writelog("calculate the number of components and ");
writelog("evaluate the formation matrix\n");
@ -70,22 +58,22 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
writelog(" --- Formula Matrix used in BASOPT calculation\n");
writelog(" --- Species | Order | ");
for (j = 0; j < ne; j++) {
jj = orderVectorElements[j];
for (size_t j = 0; j < ne; j++) {
size_t jj = orderVectorElements[j];
writelog(" ");
ename = mphase->elementName(jj);
std::string ename = mphase->elementName(jj);
print_stringTrunc(ename.c_str(), 4, 1);
writelogf("(%1d)", j);
}
writelog("\n");
for (k = 0; k < nspecies; k++) {
kk = orderVectorSpecies[k];
for (size_t k = 0; k < nspecies; k++) {
size_t kk = orderVectorSpecies[k];
writelog(" --- ");
sname = mphase->speciesName(kk);
std::string sname = mphase->speciesName(kk);
print_stringTrunc(sname.c_str(), 11, 1);
writelogf(" | %4d |", k);
for (j = 0; j < ne; j++) {
jj = orderVectorElements[j];
for (size_t j = 0; j < ne; j++) {
size_t jj = orderVectorElements[j];
double num = mphase->nAtoms(kk,jj);
writelogf("%6.1g ", num);
}
@ -117,10 +105,8 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
}
// For debugging purposes keep an unmodified copy of the array.
vector_fp molNumBase;
molNumBase = molNum;
vector_fp molNumBase = molNum;
double molSave = 0.0;
size_t jr = 0;
// Top of a loop of some sort based on the index JR. JR is the current
@ -128,11 +114,13 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
while (jr < nComponents) {
// Top of another loop point based on finding a linearly independent
// species
size_t k = npos;
while (true) {
// Search the remaining part of the mole number vector, molNum for
// the largest remaining species. Return its identity. kk is the raw
// number. k is the orderVectorSpecies index.
kk = max_element(molNum.begin(), molNum.end()) - molNum.begin();
size_t kk = max_element(molNum.begin(), molNum.end()) - molNum.begin();
size_t j;
for (j = 0; j < nspecies; j++) {
if (orderVectorSpecies[j] == kk) {
k = j;
@ -162,9 +150,9 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// Modified Gram-Schmidt Method, p. 202 Dalquist
// QR factorization of a matrix without row pivoting.
jl = jr;
size_t jl = jr;
for (j = 0; j < ne; ++j) {
jj = orderVectorElements[j];
size_t jj = orderVectorElements[j];
sm[j + jr*ne] = mphase->nAtoms(kk,jj);
}
if (jl > 0) {
@ -173,7 +161,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// different than Dalquist) R_JA_JA = 1
for (j = 0; j < jl; ++j) {
ss[j] = 0.0;
for (i = 0; i < ne; ++i) {
for (size_t i = 0; i < ne; ++i) {
ss[j] += sm[i + jr*ne] * sm[i + j*ne];
}
ss[j] /= sa[j];
@ -191,8 +179,8 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// Find the new length of the new column in Q.
// It will be used in the denominator in future row calcs.
sa[jr] = 0.0;
for (ml = 0; ml < ne; ++ml) {
tmp = sm[ml + jr*ne];
for (size_t ml = 0; ml < ne; ++ml) {
double tmp = sm[ml + jr*ne];
sa[jr] += tmp * tmp;
}
@ -205,9 +193,9 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// REARRANGE THE DATA
if (jr != k) {
if (BasisOptimize_print_lvl >= 1) {
kk = orderVectorSpecies[k];
size_t kk = orderVectorSpecies[k];
writelogf(" --- %-12.12s", mphase->speciesName(kk));
jj = orderVectorSpecies[jr];
size_t jj = orderVectorSpecies[jr];
writelogf("(%9.2g) replaces %-12.12s",
molSave, mphase->speciesName(jj));
writelogf("(%9.2g) as component %3d\n", molNum[jj], jr);
@ -253,19 +241,19 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// Note the rearrangement of elements need only be done once in the problem.
// It's actually very similar to the top of this program with ne being the
// species and nc being the elements!!
for (k = 0; k < nComponents; ++k) {
kk = orderVectorSpecies[k];
for (j = 0; j < nComponents; ++j) {
jj = orderVectorElements[j];
for (size_t k = 0; k < nComponents; ++k) {
size_t kk = orderVectorSpecies[k];
for (size_t j = 0; j < nComponents; ++j) {
size_t jj = orderVectorElements[j];
sm[j + k*ne] = mphase->nAtoms(kk, jj);
}
}
for (i = 0; i < nNonComponents; ++i) {
k = nComponents + i;
kk = orderVectorSpecies[k];
for (j = 0; j < nComponents; ++j) {
jj = orderVectorElements[j];
for (size_t i = 0; i < nNonComponents; ++i) {
size_t k = nComponents + i;
size_t kk = orderVectorSpecies[k];
for (size_t j = 0; j < nComponents; ++j) {
size_t jj = orderVectorElements[j];
formRxnMatrix[j + i * ne] = - mphase->nAtoms(kk, jj);
}
}
@ -284,39 +272,36 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
writelogf(" --- Number of Components = %d\n", nComponents);
writelog(" --- Formula Matrix:\n");
writelog(" --- Components: ");
for (k = 0; k < nComponents; k++) {
kk = orderVectorSpecies[k];
for (size_t k = 0; k < nComponents; k++) {
size_t kk = orderVectorSpecies[k];
writelogf(" %3d (%3d) ", k, kk);
}
writelog("\n --- Components Moles: ");
for (k = 0; k < nComponents; k++) {
kk = orderVectorSpecies[k];
for (size_t k = 0; k < nComponents; k++) {
size_t kk = orderVectorSpecies[k];
writelogf("%-11.3g", molNumBase[kk]);
}
writelog("\n --- NonComponent | Moles | ");
for (i = 0; i < nComponents; i++) {
kk = orderVectorSpecies[i];
for (size_t i = 0; i < nComponents; i++) {
size_t kk = orderVectorSpecies[i];
writelogf("%-11.10s", mphase->speciesName(kk));
}
writelog("\n");
for (i = 0; i < nNonComponents; i++) {
k = i + nComponents;
kk = orderVectorSpecies[k];
for (size_t i = 0; i < nNonComponents; i++) {
size_t k = i + nComponents;
size_t kk = orderVectorSpecies[k];
writelogf(" --- %3d (%3d) ", k, kk);
writelogf("%-10.10s", mphase->speciesName(kk));
writelogf("|%10.3g|", molNumBase[kk]);
// Print the negative of formRxnMatrix[]; it's easier to interpret.
for (j = 0; j < nComponents; j++) {
for (size_t j = 0; j < nComponents; j++) {
writelogf(" %6.2f", - formRxnMatrix[j + i * ne]);
}
writelog("\n");
}
writelog(" ");
for (i=0; i<77; i++) {
writelog("-");
}
writelog("\n");
writeline('-', 77);
}
return nComponents;
@ -368,19 +353,13 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
std::vector<size_t>& orderVectorSpecies,
std::vector<size_t>& orderVectorElements)
{
size_t j, k, i, jl, ml, jr, ielem, jj, kk=0;
size_t nelements = mphase->nElements();
std::string ename;
// Get the total number of species in the multiphase object
size_t nspecies = mphase->nSpecies();
double test = -1.0E10;
if (BasisOptimize_print_lvl > 0) {
writelog(" ");
for (i=0; i<77; i++) {
writelog("-");
}
writelog("\n");
writeline('-', 77);
writelog(" --- Subroutine ElemRearrange() called to ");
writelog("check stoich. coefficient matrix\n");
writelog(" --- and to rearrange the element ordering once\n");
@ -389,7 +368,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// Perhaps, initialize the element ordering
if (orderVectorElements.size() < nelements) {
orderVectorElements.resize(nelements);
for (j = 0; j < nelements; j++) {
for (size_t j = 0; j < nelements; j++) {
orderVectorElements[j] = j;
}
}
@ -398,7 +377,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// this ordering is assumed to yield the component species for the problem
if (orderVectorSpecies.size() != nspecies) {
orderVectorSpecies.resize(nspecies);
for (k = 0; k < nspecies; k++) {
for (size_t k = 0; k < nspecies; k++) {
orderVectorSpecies[k] = k;
}
}
@ -408,9 +387,9 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// zero species to the end of the element ordering.
vector_fp eAbund(nelements,0.0);
if (elementAbundances.size() != nelements) {
for (j = 0; j < nelements; j++) {
for (size_t j = 0; j < nelements; j++) {
eAbund[j] = 0.0;
for (k = 0; k < nspecies; k++) {
for (size_t k = 0; k < nspecies; k++) {
eAbund[j] += fabs(mphase->nAtoms(k, j));
}
}
@ -425,25 +404,26 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// Top of a loop of some sort based on the index JR. JR is the current
// number independent elements found.
jr = 0;
size_t jr = 0;
while (jr < nComponents) {
// Top of another loop point based on finding a linearly independent
// element
size_t k = nelements;
while (true) {
// Search the element vector. We first locate elements that are
// present in any amount. Then, we locate elements that are not
// present in any amount. Return its identity in K.
k = nelements;
for (ielem = jr; ielem < nelements; ielem++) {
size_t kk;
for (size_t ielem = jr; ielem < nelements; ielem++) {
kk = orderVectorElements[ielem];
if (eAbund[kk] != test && eAbund[kk] > 0.0) {
if (eAbund[kk] != USEDBEFORE && eAbund[kk] > 0.0) {
k = ielem;
break;
}
}
for (ielem = jr; ielem < nelements; ielem++) {
for (size_t ielem = jr; ielem < nelements; ielem++) {
kk = orderVectorElements[ielem];
if (eAbund[kk] != test) {
if (eAbund[kk] != USEDBEFORE) {
k = ielem;
break;
}
@ -460,21 +440,21 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// Assign a large negative number to the element that we have
// just found, in order to take it out of further consideration.
eAbund[kk] = test;
eAbund[kk] = USEDBEFORE;
// CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX
// LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX
// Modified Gram-Schmidt Method, p. 202 Dalquist
// QR factorization of a matrix without row pivoting.
jl = jr;
size_t jl = jr;
// Fill in the row for the current element, k, under consideration
// The row will contain the Formula matrix value for that element
// with respect to the vector of component species. (note j and k
// indices are flipped compared to the previous routine)
for (j = 0; j < nComponents; ++j) {
jj = orderVectorSpecies[j];
for (size_t j = 0; j < nComponents; ++j) {
size_t jj = orderVectorSpecies[j];
kk = orderVectorElements[k];
sm[j + jr*nComponents] = mphase->nAtoms(jj,kk);
}
@ -482,9 +462,9 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// Compute the coefficients of JA column of the the upper
// triangular R matrix, SS(J) = R_J_JR (this is slightly
// different than Dalquist) R_JA_JA = 1
for (j = 0; j < jl; ++j) {
for (size_t j = 0; j < jl; ++j) {
ss[j] = 0.0;
for (i = 0; i < nComponents; ++i) {
for (size_t i = 0; i < nComponents; ++i) {
ss[j] += sm[i + jr*nComponents] * sm[i + j*nComponents];
}
ss[j] /= sa[j];
@ -492,7 +472,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// Now make the new column, (*,JR), orthogonal to the
// previous columns
for (j = 0; j < jl; ++j) {
for (size_t j = 0; j < jl; ++j) {
for (size_t i = 0; i < nComponents; ++i) {
sm[i + jr*nComponents] -= ss[j] * sm[i + j*nComponents];
}
@ -502,7 +482,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// Find the new length of the new column in Q.
// It will be used in the denominator in future row calcs.
sa[jr] = 0.0;
for (ml = 0; ml < nComponents; ++ml) {
for (size_t ml = 0; ml < nComponents; ++ml) {
double tmp = sm[ml + jr*nComponents];
sa[jr] += tmp * tmp;
}
@ -514,7 +494,7 @@ void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
// REARRANGE THE DATA
if (jr != k) {
if (BasisOptimize_print_lvl > 0) {
kk = orderVectorElements[k];
size_t kk = orderVectorElements[k];
writelog(" --- ");
writelogf("%-2.2s", mphase->elementName(kk));
writelog("replaces ");

213
src/equil/ChemEquil.cpp
View file

@ -96,15 +96,12 @@ void ChemEquil::initialize(thermo_t& s)
// set up elemental composition matrix
size_t mneg = npos;
doublereal na, ewt;
for (size_t m = 0; m < m_mm; m++) {
for (size_t k = 0; k < m_kk; k++) {
na = s.nAtoms(k,m);
// handle the case of negative atom numbers (used to
// represent positive ions, where the 'element' is an
// electron
if (na < 0.0) {
if (s.nAtoms(k,m) < 0.0) {
// if negative atom numbers have already been specified
// for some element other than this one, throw
// an exception
@ -113,11 +110,10 @@ void ChemEquil::initialize(thermo_t& s)
"negative atom numbers allowed for only one element");
}
mneg = m;
ewt = s.atomicWeight(m);
// the element should be an electron... if it isn't
// print a warning.
if (ewt > 1.0e-3) {
if (s.atomicWeight(m) > 1.0e-3) {
writelog("WARNING: species {} has {} atoms of element {},"
" but this element is not an electron.\n",
s.speciesName(k), s.nAtoms(k,m), s.elementName(m));
@ -256,18 +252,16 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT,
m_orderVectorSpecies, m_orderVectorElements);
s.getChemPotentials(mu_RT.data());
doublereal rrt = 1.0/(GasConstant* s.temperature());
scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt);
scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(),
1.0/(GasConstant* s.temperature()));
if (ChemEquil_print_lvl > 0) {
for (size_t m = 0; m < m_nComponents; m++) {
size_t isp = m_component[m];
writelogf("isp = %d, %s\n", isp, s.speciesName(isp));
}
double pres = s.pressure();
double temp = s.temperature();
writelogf("Pressure = %g\n", pres);
writelogf("Temperature = %g\n", temp);
writelogf("Pressure = %g\n", s.pressure());
writelogf("Temperature = %g\n", s.temperature());
writelog(" id Name MF mu/RT \n");
for (size_t n = 0; n < s.nSpecies(); n++) {
writelogf("%10d %15s %10.5g %10.5g\n",
@ -316,11 +310,9 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT,
int ChemEquil::equilibrate(thermo_t& s, const char* XY,
bool useThermoPhaseElementPotentials, int loglevel)
{
vector_fp elMolesGoal(s.nElements());
initialize(s);
update(s);
copy(m_elementmolefracs.begin(), m_elementmolefracs.end(),
elMolesGoal.begin());
vector_fp elMolesGoal = m_elementmolefracs;
return equilibrate(s, XY, elMolesGoal, useThermoPhaseElementPotentials,
loglevel-1);
}
@ -330,7 +322,6 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
bool useThermoPhaseElementPotentials,
int loglevel)
{
doublereal xval, yval, tmp;
int fail = 0;
bool tempFixed = true;
int XY = _equilflag(XYstr);
@ -396,8 +387,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
// Before we do anything to change the ThermoPhase object, we calculate and
// store the two specified thermodynamic properties that we are after.
xval = m_p1->value(s);
yval = m_p2->value(s);
double xval = m_p1->value(s);
double yval = m_p2->value(s);
size_t mm = m_mm;
size_t nvar = mm + 1;
@ -409,10 +400,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
// specified property calculation.
//
// We choose the equation of the element with the highest element abundance.
size_t m;
tmp = -1.0;
double tmp = -1.0;
for (size_t im = 0; im < m_nComponents; im++) {
m = m_orderVectorElements[im];
size_t m = m_orderVectorElements[im];
if (elMolesGoal[m] > tmp) {
m_skip = m;
tmp = elMolesGoal[m];
@ -520,7 +510,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
if (s.temperature() < 100.) {
writelog("we are here {:g}\n", s.temperature());
}
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
x[m] *= 1.0 / s.RT();
}
} else {
@ -555,7 +545,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
// containing that element.
vector_fp above(nvar);
vector_fp below(nvar);
for (m = 0; m < mm; m++) {
for (size_t m = 0; m < mm; m++) {
above[m] = 200.0;
below[m] = -2000.0;
if (elMolesGoal[m] < m_elemFracCutoff && m != m_eloc) {
@ -571,20 +561,17 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
vector_fp grad(nvar, 0.0); // gradient of f = F*F/2
vector_fp oldx(nvar, 0.0); // old solution
vector_fp oldresid(nvar, 0.0);
doublereal f, oldf;
doublereal fctr = 1.0, newval;
for (int iter = 0; iter < options.maxIterations; iter++) {
// check for convergence.
equilResidual(s, x, elMolesGoal, res_trial, xval, yval);
f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin());
doublereal xx, yy, deltax, deltay;
xx = m_p1->value(s);
yy = m_p2->value(s);
deltax = (xx - xval)/xval;
deltay = (yy - yval)/yval;
double f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin());
double xx = m_p1->value(s);
double yy = m_p2->value(s);
double deltax = (xx - xval)/xval;
double deltay = (yy - yval)/yval;
bool passThis = true;
for (m = 0; m < nvar; m++) {
for (size_t m = 0; m < nvar; m++) {
double tval = options.relTolerance;
if (m < mm) {
// Special case convergence requirements for electron element.
@ -612,7 +599,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
if (iter > 0 && passThis && fabs(deltax) < options.relTolerance
&& fabs(deltay) < options.relTolerance) {
options.iterations = iter;
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
m_lambda[m] = x[m]* s.RT();
}
@ -641,7 +628,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
if (ChemEquil_print_lvl > 0) {
writelogf("Jacobian matrix %d:\n", iter);
for (m = 0; m <= m_mm; m++) {
for (size_t m = 0; m <= m_mm; m++) {
writelog(" [ ");
for (size_t n = 0; n <= m_mm; n++) {
writelog("{:10.5g} ", jac(m,n));
@ -661,7 +648,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
}
oldx = x;
oldf = f;
double oldf = f;
scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0);
// Solve the system
@ -677,9 +664,9 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
// find the factor by which the Newton step can be multiplied
// to keep the solution within bounds.
fctr = 1.0;
for (m = 0; m < nvar; m++) {
newval = x[m] + res_trial[m];
double fctr = 1.0;
for (size_t m = 0; m < nvar; m++) {
double newval = x[m] + res_trial[m];
if (newval > above[m]) {
fctr = std::max(0.0,
std::min(fctr,0.8*(above[m] - x[m])/(newval - x[m])));
@ -729,11 +716,9 @@ int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx,
double oldf, vector_fp& grad, vector_fp& step, vector_fp& x,
double& f, vector_fp& elmols, double xval, double yval)
{
double damp;
// Carry out a delta damping approach on the dimensionless element
// potentials.
damp = 1.0;
double damp = 1.0;
for (size_t m = 0; m < m_mm; m++) {
if (m == m_eloc) {
if (step[m] > 1.25) {
@ -770,9 +755,7 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x,
const vector_fp& elmFracGoal, vector_fp& resid,
doublereal xval, doublereal yval, int loglevel)
{
doublereal xx, yy;
doublereal temp = exp(x[m_mm]);
setToEquilState(s, x, temp);
setToEquilState(s, x, exp(x[m_mm]));
// residuals are the total element moles
vector_fp& elmFrac = m_elementmolefracs;
@ -802,8 +785,8 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x,
}
}
xx = m_p1->value(s);
yy = m_p2->value(s);
double xx = m_p1->value(s);
double yy = m_p2->value(s);
resid[m_mm] = xx/xval - 1.0;
resid[m_skip] = yy/yval - 1.0;
@ -823,25 +806,23 @@ void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x,
size_t len = x.size();
r0.resize(len);
r1.resize(len);
size_t n, m;
doublereal rdx, dx, xsave;
doublereal atol = 1.e-10;
equilResidual(s, x, elmols, r0, xval, yval, loglevel-1);
m_doResPerturb = false;
for (n = 0; n < len; n++) {
xsave = x[n];
dx = std::max(atol, fabs(xsave) * 1.0E-7);
for (size_t n = 0; n < len; n++) {
double xsave = x[n];
double dx = std::max(atol, fabs(xsave) * 1.0E-7);
x[n] = xsave + dx;
dx = x[n] - xsave;
rdx = 1.0/dx;
double rdx = 1.0/dx;
// calculate perturbed residual
equilResidual(s, x, elmols, r1, xval, yval, loglevel-1);
// compute nth column of Jacobian
for (m = 0; m < x.size(); m++) {
for (size_t m = 0; m < x.size(); m++) {
jac(m, n) = (r1[m] - r0[m])*rdx;
}
x[n] = xsave;
@ -855,7 +836,6 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t,
double pressureConst)
{
double n_t_calc = 0.0;
double tmp;
// Calculate the activity coefficients of the solution, at the previous
// solution state.
@ -865,7 +845,7 @@ double ChemEquil::calcEmoles(thermo_t& s, vector_fp& x, const double& n_t,
s.getActivityCoefficients(actCoeff.data());
for (size_t k = 0; k < m_kk; k++) {
tmp = - (m_muSS_RT[k] + log(actCoeff[k]));
double tmp = - (m_muSS_RT[k] + log(actCoeff[k]));
for (size_t m = 0; m < m_mm; m++) {
tmp += nAtoms(k,m) * x[m];
}
@ -893,11 +873,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
// things go drastically wrong, we will restore the saved state.
vector_fp state;
s.saveState(state);
double tmp, sum;
bool modifiedMatrix = false;
size_t neq = m_mm+1;
int retn = 1;
size_t m, n, k, im;
DenseMatrix a1(neq, neq, 0.0);
vector_fp b(neq, 0.0);
vector_fp n_i(m_kk,0.0);
@ -920,14 +898,14 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
vector_fp eMolesCalc(m_mm, 0.0);
vector_fp eMolesFix(m_mm, 0.0);
double elMolesTotal = 0.0;
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
elMolesTotal += elMoles[m];
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
eMolesFix[m] += nAtoms(k,m) * n_i[k];
}
}
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
if (elMoles[m] > 1.0E-70) {
x[m] = clip(x[m], -100.0, 50.0);
} else {
@ -936,16 +914,15 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
double n_t = 0.0;
double sum2 = 0.0;
double nAtomsMax = 1.0;
s.setMoleFractions(Xmol_i_calc.data());
s.setPressure(pressureConst);
s.getActivityCoefficients(actCoeff.data());
for (k = 0; k < m_kk; k++) {
tmp = - (m_muSS_RT[k] + log(actCoeff[k]));
sum2 = 0.0;
for (m = 0; m < m_mm; m++) {
sum = nAtoms(k,m);
for (size_t k = 0; k < m_kk; k++) {
double tmp = - (m_muSS_RT[k] + log(actCoeff[k]));
double sum2 = 0.0;
for (size_t m = 0; m < m_mm; m++) {
double sum = nAtoms(k,m);
tmp += sum * x[m];
sum2 += sum;
nAtomsMax = std::max(nAtomsMax, sum2);
@ -959,25 +936,23 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
if (ChemEquil_print_lvl > 0) {
writelog("estimateEP_Brinkley::\n\n");
double temp = s.temperature();
double pres = s.pressure();
writelogf("temp = %g\n", temp);
writelogf("pres = %g\n", pres);
writelogf("temp = %g\n", s.temperature());
writelogf("pres = %g\n", s.pressure());
writelog("Initial mole numbers and mu_SS:\n");
writelog(" Name MoleNum mu_SS actCoeff\n");
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
writelogf("%15s %13.5g %13.5g %13.5g\n",
s.speciesName(k), n_i[k], m_muSS_RT[k], actCoeff[k]);
}
writelogf("Initial n_t = %10.5g\n", n_t);
writelog("Comparison of Goal Element Abundance with Initial Guess:\n");
writelog(" eName eCurrent eGoal\n");
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
writelogf("%5s %13.5g %13.5g\n",
s.elementName(m), eMolesFix[m], elMoles[m]);
}
}
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
if (m != m_eloc && elMoles[m] <= options.absElemTol) {
x[m] = -200.;
}
@ -986,7 +961,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
// Main Loop.
for (int iter = 0; iter < 20* options.maxIterations; iter++) {
// Save the old solution
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
x_old[m] = x[m];
}
x_old[m_mm] = n_t;
@ -998,29 +973,28 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
double n_t_calc = calcEmoles(s, x, n_t, Xmol_i_calc, eMolesCalc, n_i_calc,
pressureConst);
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
Xmol_i_calc[k] = n_i_calc[k]/n_t_calc;
}
if (ChemEquil_print_lvl > 0) {
writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n");
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
writelogf("%15s: %10.5g %10.5g\n",
s.speciesName(k), n_i_calc[k], Xmol_i_calc[k]);
}
writelogf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc);
writelogf("(iter %d) element moles bal: Goal Calculated\n", iter);
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
writelogf(" %8s: %10.5g %10.5g \n",
s.elementName(m), elMoles[m], eMolesCalc[m]);
}
}
double nCutoff;
bool normalStep = true;
// Decide if we are to do a normal step or a modified step
size_t iM = npos;
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
if (elMoles[m] > 0.001 * elMolesTotal) {
if (eMolesCalc[m] > 1000. * elMoles[m]) {
normalStep = false;
@ -1038,8 +1012,8 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
if (!normalStep) {
beta = 1.0;
resid[m_mm] = 0.0;
for (im = 0; im < m_mm; im++) {
m = m_orderVectorElements[im];
for (size_t im = 0; im < m_mm; im++) {
size_t m = m_orderVectorElements[im];
resid[m] = 0.0;
if (im < m_nComponents && elMoles[m] > 0.001 * elMolesTotal) {
if (eMolesCalc[m] > 1000. * elMoles[m]) {
@ -1082,25 +1056,25 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
// Jordan factorization scheme in the future. For Example the scheme
// below would fail for the set: HCl NH4Cl, NH3. Hopefully, it's
// caught by the equal rows logic below.
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
lumpSum[m] = 1;
}
nCutoff = 1.0E-9 * n_t_calc;
double nCutoff = 1.0E-9 * n_t_calc;
if (ChemEquil_print_lvl > 0) {
writelog(" Lump Sum Elements Calculation: \n");
}
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
size_t kMSp = npos;
size_t kMSp2 = npos;
int nSpeciesWithElem = 0;
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (n_i_calc[k] > nCutoff && fabs(nAtoms(k,m)) > 0.001) {
nSpeciesWithElem++;
if (kMSp != npos) {
kMSp2 = k;
double factor = fabs(nAtoms(kMSp,m) / nAtoms(kMSp2,m));
for (n = 0; n < m_mm; n++) {
for (size_t n = 0; n < m_mm; n++) {
if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) {
lumpSum[m] = 0;
break;
@ -1118,19 +1092,19 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
// Formulate the matrix.
for (im = 0; im < m_mm; im++) {
m = m_orderVectorElements[im];
for (size_t im = 0; im < m_mm; im++) {
size_t m = m_orderVectorElements[im];
if (im < m_nComponents) {
for (n = 0; n < m_mm; n++) {
for (size_t n = 0; n < m_mm; n++) {
a1(m,n) = 0.0;
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
a1(m,n) += nAtoms(k,m) * nAtoms(k,n) * n_i_calc[k];
}
}
a1(m,m_mm) = eMolesCalc[m];
a1(m_mm, m) = eMolesCalc[m];
} else {
for (n = 0; n <= m_mm; n++) {
for (size_t n = 0; n <= m_mm; n++) {
a1(m,n) = 0.0;
}
a1(m,m) = 1.0;
@ -1140,9 +1114,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
// Formulate the residual, resid, and the estimate for the
// convergence criteria, sum
sum = 0.0;
for (im = 0; im < m_mm; im++) {
m = m_orderVectorElements[im];
double sum = 0.0;
for (size_t im = 0; im < m_mm; im++) {
size_t m = m_orderVectorElements[im];
if (im < m_nComponents) {
resid[m] = elMoles[m] - eMolesCalc[m];
} else {
@ -1154,6 +1128,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
// criteria by a condition limited by finite precision of
// inverting a matrix. Other equations with just positive
// coefficients aren't limited by this.
double tmp;
if (m == m_eloc) {
tmp = resid[m] / (elMoles[m] + elMolesTotal*1.0E-6 + options.absElemTol);
} else {
@ -1162,12 +1137,12 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
sum += tmp * tmp;
}
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
if (a1(m,m) < 1.0E-50) {
if (ChemEquil_print_lvl > 0) {
writelogf(" NOTE: Diagonalizing the analytical Jac row %d\n", m);
}
for (n = 0; n < m_mm; n++) {
for (size_t n = 0; n < m_mm; n++) {
a1(m,n) = 0.0;
}
a1(m,m) = 1.0;
@ -1185,17 +1160,16 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
if (ChemEquil_print_lvl > 0) {
writelog("Matrix:\n");
for (m = 0; m <= m_mm; m++) {
for (size_t m = 0; m <= m_mm; m++) {
writelog(" [");
for (n = 0; n <= m_mm; n++) {
for (size_t n = 0; n <= m_mm; n++) {
writelogf(" %10.5g", a1(m,n));
}
writelogf("] = %10.5g\n", resid[m]);
}
}
tmp = resid[m_mm] /(n_t + 1.0E-15);
sum += tmp * tmp;
sum += pow(resid[m_mm] /(n_t + 1.0E-15), 2);
if (ChemEquil_print_lvl > 0) {
writelogf("(it %d) Convergence = %g\n", iter, sum);
}
@ -1210,16 +1184,16 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
// Row Sum scaling
for (m = 0; m <= m_mm; m++) {
tmp = 0.0;
for (n = 0; n <= m_mm; n++) {
for (size_t m = 0; m <= m_mm; m++) {
double tmp = 0.0;
for (size_t n = 0; n <= m_mm; n++) {
tmp += fabs(a1(m,n));
}
if (m < m_mm && tmp < 1.0E-30) {
if (ChemEquil_print_lvl > 0) {
writelogf(" NOTE: Diagonalizing row %d\n", m);
}
for (n = 0; n <= m_mm; n++) {
for (size_t n = 0; n <= m_mm; n++) {
if (n != m) {
a1(m,n) = 0.0;
a1(n,m) = 0.0;
@ -1227,7 +1201,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
tmp = 1.0/tmp;
for (n = 0; n <= m_mm; n++) {
for (size_t n = 0; n <= m_mm; n++) {
a1(m,n) *= tmp;
}
resid[m] *= tmp;
@ -1235,9 +1209,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
if (ChemEquil_print_lvl > 0) {
writelog("Row Summed Matrix:\n");
for (m = 0; m <= m_mm; m++) {
for (size_t m = 0; m <= m_mm; m++) {
writelog(" [");
for (n = 0; n <= m_mm; n++) {
for (size_t n = 0; n <= m_mm; n++) {
writelogf(" %10.5g", a1(m,n));
}
writelogf("] = %10.5g\n", resid[m]);
@ -1264,11 +1238,11 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
// 1.0E-3. If two rows are anywhere close to being equivalent, the
// algorithm can get stuck in an oscillatory mode.
modifiedMatrix = false;
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
size_t sameAsRow = npos;
for (size_t im = 0; im < m; im++) {
bool theSame = true;
for (n = 0; n < m_mm; n++) {
for (size_t n = 0; n < m_mm; n++) {
if (fabs(a1(m,n) - a1(im,n)) > 1.0E-7) {
theSame = false;
break;
@ -1287,7 +1261,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
modifiedMatrix = true;
for (n = 0; n < m_mm; n++) {
for (size_t n = 0; n < m_mm; n++) {
if (n != m) {
a1(m,m) += fabs(a1(m,n));
a1(m,n) = 0.0;
@ -1298,9 +1272,9 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
if (ChemEquil_print_lvl > 0 && modifiedMatrix) {
writelog("Row Summed, MODIFIED Matrix:\n");
for (m = 0; m <= m_mm; m++) {
for (size_t m = 0; m <= m_mm; m++) {
writelog(" [");
for (n = 0; n <= m_mm; n++) {
for (size_t n = 0; n <= m_mm; n++) {
writelogf(" %10.5g", a1(m,n));
}
writelogf("] = %10.5g\n", resid[m]);
@ -1320,7 +1294,7 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
// Figure out the damping coefficient: Use a delta damping
// coefficient formulation: magnitude of change is capped to exp(1).
beta = 1.0;
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
if (resid[m] > 1.0) {
beta = std::min(beta, 1.0 / resid[m]);
}
@ -1333,14 +1307,14 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x,
}
}
// Update the solution vector
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
x[m] += beta * resid[m];
}
n_t *= exp(beta * resid[m_mm]);
if (ChemEquil_print_lvl > 0) {
writelogf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter);
for (m = 0; m < m_mm; m++) {
for (size_t m = 0; m < m_mm; m++) {
writelogf(" %5s %10.5g %10.5g %10.5g\n",
s.elementName(m), x_old[m], x[m], resid[m]);
}
@ -1372,13 +1346,12 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal)
return;
}
s.getMoleFractions(m_molefractions.data());
size_t k;
size_t maxPosEloc = npos;
size_t maxNegEloc = npos;
double maxPosVal = -1.0;
double maxNegVal = -1.0;
if (ChemEquil_print_lvl > 0) {
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (nAtoms(k,m_eloc) > 0.0 && m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) {
maxPosVal = m_molefractions[k];
maxPosEloc = k;
@ -1392,7 +1365,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal)
double sumPos = 0.0;
double sumNeg = 0.0;
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (nAtoms(k,m_eloc) > 0.0) {
sumPos += nAtoms(k,m_eloc) * m_molefractions[k];
}
@ -1412,7 +1385,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal)
s.speciesName(maxPosEloc),
m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor);
}
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (nAtoms(k,m_eloc) > 0.0) {
m_molefractions[k] *= factor;
}
@ -1424,7 +1397,7 @@ void ChemEquil::adjustEloc(thermo_t& s, vector_fp& elMolesGoal)
s.speciesName(maxNegEloc),
m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor);
}
for (k = 0; k < m_kk; k++) {
for (size_t k = 0; k < m_kk; k++) {
if (nAtoms(k,m_eloc) < 0.0) {
m_molefractions[k] *= factor;
}

200
src/equil/MultiPhase.cpp
View file

@ -68,8 +68,7 @@ MultiPhase& MultiPhase::operator=(const MultiPhase& right)
void MultiPhase::addPhases(MultiPhase& mix)
{
size_t n;
for (n = 0; n < mix.nPhases(); n++) {
for (size_t n = 0; n < mix.nPhases(); n++) {
addPhase(mix.m_phase[n], mix.m_moles[n]);
}
}
@ -77,9 +76,7 @@ void MultiPhase::addPhases(MultiPhase& mix)
void MultiPhase::addPhases(std::vector<ThermoPhase*>& phases,
const vector_fp& phaseMoles)
{
size_t np = phases.size();
size_t n;
for (n = 0; n < np; n++) {
for (size_t n = 0; n < phases.size(); n++) {
addPhase(phases[n], phaseMoles[n]);
}
init();
@ -105,11 +102,9 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
// determine if this phase has new elements for each new element, add an
// entry in the map from names to index number + 1:
string ename;
// iterate over the elements in this phase
size_t m, nel = p->nElements();
for (m = 0; m < nel; m++) {
ename = p->elementName(m);
for (size_t m = 0; m < p->nElements(); m++) {
string ename = p->elementName(m);
// if no entry is found for this element name, then it is a new element.
// In this case, add the name to the list of names, increment the
@ -154,9 +149,6 @@ void MultiPhase::init()
if (m_init) {
return;
}
size_t ip, kp, k = 0, nsp, m;
size_t mlocal;
string sym;
// allocate space for the atomic composition matrix
m_atoms.resize(m_nel, m_nsp, 0.0);
@ -165,15 +157,14 @@ void MultiPhase::init()
// iterate over the elements
// -> fill in m_atoms(m,k), m_snames(k), m_spphase(k), m_spstart(ip)
for (m = 0; m < m_nel; m++) {
sym = m_enames[m];
k = 0;
for (size_t m = 0; m < m_nel; m++) {
size_t k = 0;
// iterate over the phases
for (ip = 0; ip < nPhases(); ip++) {
for (size_t ip = 0; ip < nPhases(); ip++) {
ThermoPhase* p = m_phase[ip];
nsp = p->nSpecies();
mlocal = p->elementIndex(sym);
for (kp = 0; kp < nsp; kp++) {
size_t nsp = p->nSpecies();
size_t mlocal = p->elementIndex(m_enames[m]);
for (size_t kp = 0; kp < nsp; kp++) {
if (mlocal != npos) {
m_atoms(m, k) = p->nAtoms(kp, mlocal);
}
@ -189,18 +180,6 @@ void MultiPhase::init()
}
}
if (m_eloc != npos) {
doublereal esum;
for (k = 0; k < m_nsp; k++) {
esum = 0.0;
for (m = 0; m < m_nel; m++) {
if (m != m_eloc) {
esum += m_atoms(m,k) * m_atomicNumber[m];
}
}
}
}
// set the initial composition within each phase to the
// mole fractions stored in the phase objects
m_init = true;
@ -241,13 +220,12 @@ doublereal MultiPhase::speciesMoles(size_t k) const
doublereal MultiPhase::elementMoles(size_t m) const
{
doublereal sum = 0.0, phasesum;
size_t i, k = 0, ik, nsp;
for (i = 0; i < nPhases(); i++) {
phasesum = 0.0;
nsp = m_phase[i]->nSpecies();
for (ik = 0; ik < nsp; ik++) {
k = speciesIndex(ik, i);
doublereal sum = 0.0;
for (size_t i = 0; i < nPhases(); i++) {
double phasesum = 0.0;
size_t nsp = m_phase[i]->nSpecies();
for (size_t ik = 0; ik < nsp; ik++) {
size_t k = speciesIndex(ik, i);
phasesum += m_atoms(m,k)*m_moleFractions[k];
}
sum += phasesum * m_moles[i];
@ -258,8 +236,7 @@ doublereal MultiPhase::elementMoles(size_t m) const
doublereal MultiPhase::charge() const
{
doublereal sum = 0.0;
size_t i;
for (i = 0; i < nPhases(); i++) {
for (size_t i = 0; i < nPhases(); i++) {
sum += phaseCharge(i);
}
return sum;
@ -284,9 +261,9 @@ size_t MultiPhase::speciesIndex(const std::string& speciesName, const std::strin
doublereal MultiPhase::phaseCharge(size_t p) const
{
doublereal phasesum = 0.0;
size_t ik, k, nsp = m_phase[p]->nSpecies();
for (ik = 0; ik < nsp; ik++) {
k = speciesIndex(ik, p);
size_t nsp = m_phase[p]->nSpecies();
for (size_t ik = 0; ik < nsp; ik++) {
size_t k = speciesIndex(ik, p);
phasesum += m_phase[p]->charge(ik)*m_moleFractions[k];
}
return Faraday*phasesum*m_moles[p];
@ -294,9 +271,9 @@ doublereal MultiPhase::phaseCharge(size_t p) const
void MultiPhase::getChemPotentials(doublereal* mu) const
{
size_t i, loc = 0;
updatePhases();
for (i = 0; i < nPhases(); i++) {
size_t loc = 0;
for (size_t i = 0; i < nPhases(); i++) {
m_phase[i]->getChemPotentials(mu + loc);
loc += m_phase[i]->nSpecies();
}
@ -305,10 +282,10 @@ void MultiPhase::getChemPotentials(doublereal* mu) const
void MultiPhase::getValidChemPotentials(doublereal not_mu,
doublereal* mu, bool standard) const
{
size_t i, loc = 0;
updatePhases();
// iterate over the phases
for (i = 0; i < nPhases(); i++) {
size_t loc = 0;
for (size_t i = 0; i < nPhases(); i++) {
if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
if (!standard) {
m_phase[i]->getChemPotentials(mu + loc);
@ -333,10 +310,9 @@ bool MultiPhase::solutionSpecies(size_t k) const
doublereal MultiPhase::gibbs() const
{
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < nPhases(); i++) {
for (size_t i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->gibbs_mole() * m_moles[i];
}
@ -346,10 +322,9 @@ doublereal MultiPhase::gibbs() const
doublereal MultiPhase::enthalpy() const
{
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < nPhases(); i++) {
for (size_t i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->enthalpy_mole() * m_moles[i];
}
@ -359,10 +334,9 @@ doublereal MultiPhase::enthalpy() const
doublereal MultiPhase::IntEnergy() const
{
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < nPhases(); i++) {
for (size_t i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->intEnergy_mole() * m_moles[i];
}
@ -372,10 +346,9 @@ doublereal MultiPhase::IntEnergy() const
doublereal MultiPhase::entropy() const
{
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < nPhases(); i++) {
for (size_t i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->entropy_mole() * m_moles[i];
}
@ -385,10 +358,9 @@ doublereal MultiPhase::entropy() const
doublereal MultiPhase::cp() const
{
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < nPhases(); i++) {
for (size_t i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->cp_mole() * m_moles[i];
}
@ -430,13 +402,12 @@ void MultiPhase::getMoles(doublereal* molNum) const
{
// First copy in the mole fractions
copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
size_t ik;
doublereal* dtmp = molNum;
for (size_t ip = 0; ip < nPhases(); ip++) {
doublereal phasemoles = m_moles[ip];
ThermoPhase* p = m_phase[ip];
size_t nsp = p->nSpecies();
for (ik = 0; ik < nsp; ik++) {
for (size_t ik = 0; ik < nsp; ik++) {
*(dtmp++) *= phasemoles;
}
}
@ -447,14 +418,13 @@ void MultiPhase::setMoles(const doublereal* n)
if (!m_init) {
init();
}
size_t ip, loc = 0;
size_t ik, k = 0, nsp;
doublereal phasemoles;
for (ip = 0; ip < nPhases(); ip++) {
size_t loc = 0;
size_t k = 0;
for (size_t ip = 0; ip < nPhases(); ip++) {
ThermoPhase* p = m_phase[ip];
nsp = p->nSpecies();
phasemoles = 0.0;
for (ik = 0; ik < nsp; ik++) {
size_t nsp = p->nSpecies();
double phasemoles = 0.0;
for (size_t ik = 0; ik < nsp; ik++) {
phasemoles += n[k];
k++;
}
@ -502,9 +472,8 @@ void MultiPhase::setState_TPMoles(const doublereal T, const doublereal Pres,
void MultiPhase::getElemAbundances(doublereal* elemAbundances) const
{
size_t eGlobal;
calcElemAbundances();
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
elemAbundances[eGlobal] = m_elemAbundances[eGlobal];
}
}
@ -512,20 +481,18 @@ void MultiPhase::getElemAbundances(doublereal* elemAbundances) const
void MultiPhase::calcElemAbundances() const
{
size_t loc = 0;
size_t eGlobal;
size_t ik, kGlobal;
doublereal spMoles;
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
m_elemAbundances[eGlobal] = 0.0;
}
for (size_t ip = 0; ip < nPhases(); ip++) {
ThermoPhase* p = m_phase[ip];
size_t nspPhase = p->nSpecies();
doublereal phasemoles = m_moles[ip];
for (ik = 0; ik < nspPhase; ik++) {
kGlobal = loc + ik;
for (size_t ik = 0; ik < nspPhase; ik++) {
size_t kGlobal = loc + ik;
spMoles = m_moleFractions[kGlobal] * phasemoles;
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
for (size_t eGlobal = 0; eGlobal < m_nel; eGlobal++) {
m_elemAbundances[eGlobal] += m_atoms(eGlobal, kGlobal) * spMoles;
}
}
@ -535,9 +502,8 @@ void MultiPhase::calcElemAbundances() const
doublereal MultiPhase::volume() const
{
int i;
doublereal sum = 0;
for (i = 0; i < int(nPhases()); i++) {
for (size_t i = 0; i < nPhases(); i++) {
double vol = 1.0/m_phase[i]->molarDensity();
sum += m_moles[i] * vol;
}
@ -549,14 +515,7 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
int loglevel)
{
bool strt = false;
doublereal dt;
doublereal h0;
int n;
doublereal hnow, herr = 1.0;
doublereal snow, s0;
doublereal Tlow = -1.0, Thigh = -1.0;
doublereal Hlow = Undef, Hhigh = Undef, tnew;
doublereal dta=0.0, dtmax, cpb;
doublereal dta = 0.0;
if (!m_init) {
init();
}
@ -566,10 +525,11 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
MultiPhaseEquil e(this);
return e.equilibrate(XY, err, maxsteps, loglevel);
} else if (XY == HP) {
h0 = enthalpy();
Tlow = 0.5*m_Tmin; // lower bound on T
Thigh = 2.0*m_Tmax; // upper bound on T
for (n = 0; n < maxiter; n++) {
double h0 = enthalpy();
double Tlow = 0.5*m_Tmin; // lower bound on T
double Thigh = 2.0*m_Tmax; // upper bound on T
doublereal Hlow = Undef, Hhigh = Undef;
for (int n = 0; n < maxiter; n++) {
// if 'strt' is false, the current composition will be used as
// the starting estimate; otherwise it will be estimated
MultiPhaseEquil e(this, strt);
@ -578,7 +538,7 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
try {
e.equilibrate(TP, err, maxsteps, loglevel);
hnow = enthalpy();
double hnow = enthalpy();
// the equilibrium enthalpy monotonically increases with T;
// if the current value is below the target, the we know the
// current temperature is too low. Set
@ -595,25 +555,26 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
Hhigh = hnow;
}
}
double dt;
if (Hlow != Undef && Hhigh != Undef) {
cpb = (Hhigh - Hlow)/(Thigh - Tlow);
double cpb = (Hhigh - Hlow)/(Thigh - Tlow);
dt = (h0 - hnow)/cpb;
dta = fabs(dt);
dtmax = 0.5*fabs(Thigh - Tlow);
double dtmax = 0.5*fabs(Thigh - Tlow);
if (dta > dtmax) {
dt *= dtmax/dta;
}
} else {
tnew = sqrt(Tlow*Thigh);
double tnew = sqrt(Tlow*Thigh);
dt = tnew - m_temp;
}
herr = fabs((h0 - hnow)/h0);
double herr = fabs((h0 - hnow)/h0);
if (herr < err) {
return err;
}
tnew = m_temp + dt;
double tnew = m_temp + dt;
if (tnew < 0.0) {
tnew = 0.5*m_temp;
}
@ -629,7 +590,7 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
if (!strt) {
strt = true;
} else {
tnew = 0.5*(m_temp + Thigh);
double tnew = 0.5*(m_temp + Thigh);
if (fabs(tnew - m_temp) < 1.0) {
tnew = m_temp + 1.0;
}
@ -640,31 +601,31 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
throw CanteraError("MultiPhase::equilibrate_MultiPhaseEquil",
"No convergence for T");
} else if (XY == SP) {
s0 = entropy();
Tlow = 1.0; // lower bound on T
Thigh = 1.0e6; // upper bound on T
for (n = 0; n < maxiter; n++) {
double s0 = entropy();
double Tlow = 1.0; // lower bound on T
double Thigh = 1.0e6; // upper bound on T
for (int n = 0; n < maxiter; n++) {
MultiPhaseEquil e(this, strt);
try {
e.equilibrate(TP, err, maxsteps, loglevel);
snow = entropy();
double snow = entropy();
if (snow < s0) {
Tlow = std::max(Tlow, m_temp);
} else {
Thigh = std::min(Thigh, m_temp);
}
dt = (s0 - snow)*m_temp/cp();
dtmax = 0.5*fabs(Thigh - Tlow);
double dt = (s0 - snow)*m_temp/cp();
double dtmax = 0.5*fabs(Thigh - Tlow);
dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
dta = fabs(dt);
if (dta > dtmax) {
dt *= dtmax/dta;
}
if (herr < err || dta < 1.0e-4) {
if (dta < 1.0e-4) {
return err;
}
tnew = m_temp + dt;
double tnew = m_temp + dt;
setTemperature(tnew);
// if the size of Delta T is not too large, use
@ -676,7 +637,7 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
if (!strt) {
strt = true;
} else {
tnew = 0.5*(m_temp + Thigh);
double tnew = 0.5*(m_temp + Thigh);
setTemperature(tnew);
}
}
@ -685,25 +646,22 @@ double MultiPhase::equilibrate_MultiPhaseEquil(int XY, doublereal err,
"No convergence for T");
} else if (XY == TV) {
doublereal v0 = volume();
doublereal dVdP;
int n;
bool start = true;
doublereal vnow, pnow, verr;
for (n = 0; n < maxiter; n++) {
pnow = pressure();
for (int n = 0; n < maxiter; n++) {
double pnow = pressure();
MultiPhaseEquil e(this, start);
start = false;
e.equilibrate(TP, err, maxsteps, loglevel);
vnow = volume();
verr = fabs((v0 - vnow)/v0);
double vnow = volume();
double verr = fabs((v0 - vnow)/v0);
if (verr < err) {
return err;
}
// find dV/dP
setPressure(pnow*1.01);
dVdP = (volume() - vnow)/(0.01*pnow);
double dVdP = (volume() - vnow)/(0.01*pnow);
setPressure(pnow + 0.5*(v0 - vnow)/dVdP);
}
} else {
@ -852,11 +810,8 @@ std::string MultiPhase::phaseName(const size_t iph) const
int MultiPhase::phaseIndex(const std::string& pName) const
{
std::string tmp;
for (int iph = 0; iph < (int) nPhases(); iph++) {
const ThermoPhase* tptr = m_phase[iph];
tmp = tptr->id();
if (tmp == pName) {
if (m_phase[iph]->id() == pName) {
return iph;
}
}
@ -890,8 +845,8 @@ bool MultiPhase::tempOK(const size_t p) const
void MultiPhase::uploadMoleFractionsFromPhases()
{
size_t ip, loc = 0;
for (ip = 0; ip < nPhases(); ip++) {
size_t loc = 0;
for (size_t ip = 0; ip < nPhases(); ip++) {
ThermoPhase* p = m_phase[ip];
p->getMoleFractions(&m_moleFractions[loc]);
loc += p->nSpecies();
@ -901,11 +856,10 @@ void MultiPhase::uploadMoleFractionsFromPhases()
void MultiPhase::updatePhases() const
{
size_t p, nsp, loc = 0;
for (p = 0; p < nPhases(); p++) {
nsp = m_phase[p]->nSpecies();
size_t loc = 0;
for (size_t p = 0; p < nPhases(); p++) {
m_phase[p]->setState_TPX(m_temp, m_press, &m_moleFractions[loc]);
loc += nsp;
loc += m_phase[p]->nSpecies();
m_temp_OK[p] = true;
if (m_temp < m_phase[p]->minTemp() || m_temp > m_phase[p]->maxTemp()) {
m_temp_OK[p] = false;

199
src/equil/MultiPhaseEquil.cpp
View file

@ -20,14 +20,13 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
m_press = mix->pressure();
m_temp = mix->temperature();
size_t m, k;
m_force = true;
m_nel = 0;
m_nsp = 0;
m_eloc = 1000;
m_incl_species.resize(m_nsp_mix,1);
m_incl_element.resize(m_nel_mix,1);
for (m = 0; m < m_nel_mix; m++) {
for (size_t m = 0; m < m_nel_mix; m++) {
string enm = mix->elementName(m);
// element 'E' or 'e' represents an electron; this requires special
// handling, so save its index for later use
@ -40,7 +39,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
// number of 'atoms' of electrons (positive ions).
if (m_mix->elementMoles(m) <= 0.0 && m != m_eloc) {
m_incl_element[m] = 0;
for (k = 0; k < m_nsp_mix; k++) {
for (size_t k = 0; k < m_nsp_mix; k++) {
if (m_mix->nAtoms(k,m) != 0.0) {
m_incl_species[k] = 0;
}
@ -55,7 +54,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
m_nel++;
}
// add the included elements other than electrons
for (m = 0; m < m_nel_mix; m++) {
for (size_t m = 0; m < m_nel_mix; m++) {
if (m_incl_element[m] == 1 && m != m_eloc) {
m_nel++;
m_element.push_back(m);
@ -71,9 +70,8 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
// only extend to 273.15 K, and give unphysical results above this
// temperature, leading (incorrectly) to Gibbs free energies at high
// temperature lower than for liquid water.
size_t ip;
for (k = 0; k < m_nsp_mix; k++) {
ip = m_mix->speciesPhaseIndex(k);
for (size_t k = 0; k < m_nsp_mix; k++) {
size_t ip = m_mix->speciesPhaseIndex(k);
if (!m_mix->solutionSpecies(k) &&
!m_mix->tempOK(ip)) {
m_incl_species[k] = 0;
@ -88,7 +86,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
}
// Now build the list of all species to be included in the calculation.
for (k = 0; k < m_nsp_mix; k++) {
for (size_t k = 0; k < m_nsp_mix; k++) {
if (m_incl_species[k] ==1) {
m_nsp++;
m_species.push_back(k);
@ -107,8 +105,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
m_dxi.resize(nFree());
// initialize the mole numbers to the mixture composition
size_t ik;
for (ik = 0; ik < m_nsp; ik++) {
for (size_t ik = 0; ik < m_nsp; ik++) {
m_moles[ik] = m_mix->speciesMoles(m_species[ik]);
}
@ -121,9 +118,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
m_A.resize(m_nel, m_nsp, 0.0);
m_N.resize(m_nsp, nFree());
m_order.resize(std::max(m_nsp, m_nel), 0);
for (k = 0; k < m_nsp; k++) {
m_order[k] = k;
}
iota(m_order.begin(), m_order.begin() + m_nsp, 0);
// if the 'start' flag is set, estimate the initial mole numbers by doing a
// linear Gibbs minimization. In this case, only the elemental composition
@ -141,7 +136,7 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
unsort(m_work);
}
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
m_moles[k] += m_work[k];
m_lastmoles[k] = m_moles[k];
if (m_mix->solutionSpecies(m_species[k])) {
@ -181,8 +176,7 @@ doublereal MultiPhaseEquil::equilibrate(int XY, doublereal err,
void MultiPhaseEquil::updateMixMoles()
{
fill(m_work3.begin(), m_work3.end(), 0.0);
size_t k;
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
m_work3[m_species[k]] = m_moles[k];
}
m_mix->setMoles(m_work3.data());
@ -191,8 +185,7 @@ void MultiPhaseEquil::updateMixMoles()
void MultiPhaseEquil::finish()
{
fill(m_work3.begin(), m_work3.end(), 0.0);
size_t k;
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
m_work3[m_species[k]] = (m_moles[k] > 0.0 ? m_moles[k] : 0.0);
}
m_mix->setMoles(m_work3.data());
@ -200,13 +193,9 @@ void MultiPhaseEquil::finish()
int MultiPhaseEquil::setInitialMoles(int loglevel)
{
size_t ik, j;
double not_mu = 1.0e12;
m_mix->getValidChemPotentials(not_mu, m_mu.data(), true);
doublereal dg_rt;
int idir;
double nu;
double delta_xi, dxi_min = 1.0e10;
double dxi_min = 1.0e10;
bool redo = true;
int iter = 0;
@ -220,18 +209,18 @@ int MultiPhaseEquil::setInitialMoles(int loglevel)
}
// loop over all reactions
for (j = 0; j < nFree(); j++) {
dg_rt = 0.0;
for (size_t j = 0; j < nFree(); j++) {
double dg_rt = 0.0;
dxi_min = 1.0e10;
for (ik = 0; ik < m_nsp; ik++) {
for (size_t ik = 0; ik < m_nsp; ik++) {
dg_rt += mu(ik) * m_N(ik,j);
}
// fwd or rev direction
idir = (dg_rt < 0.0 ? 1 : -1);
int idir = (dg_rt < 0.0 ? 1 : -1);
for (ik = 0; ik < m_nsp; ik++) {
nu = m_N(ik, j);
for (size_t ik = 0; ik < m_nsp; ik++) {
double nu = m_N(ik, j);
// set max change in progress variable by
// non-negativity requirement
@ -239,7 +228,7 @@ int MultiPhaseEquil::setInitialMoles(int loglevel)
// isn't zero. This causes differences between
// optimized and debug versions of the code
if (nu*idir < 0) {
delta_xi = fabs(0.99*moles(ik)/nu);
double delta_xi = fabs(0.99*moles(ik)/nu);
// if a component has nearly zero moles, redo
// with a new set of components
if (!redo && delta_xi < 1.0e-10 && ik < m_nel) {
@ -249,7 +238,7 @@ int MultiPhaseEquil::setInitialMoles(int loglevel)
}
}
// step the composition by dxi_min
for (ik = 0; ik < m_nsp; ik++) {
for (size_t ik = 0; ik < m_nsp; ik++) {
moles(ik) += m_N(ik, j) * idir*dxi_min;
}
}
@ -261,36 +250,33 @@ int MultiPhaseEquil::setInitialMoles(int loglevel)
void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
{
size_t m, k, j;
// if the input species array has the wrong size, ignore it
// and consider the species for components in declaration order.
if (order.size() != m_nsp) {
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
m_order[k] = k;
}
} else {
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
m_order[k] = order[k];
}
}
size_t nRows = m_nel;
size_t nColumns = m_nsp;
doublereal fctr;
// set up the atomic composition matrix
for (m = 0; m < nRows; m++) {
for (k = 0; k < nColumns; k++) {
for (size_t m = 0; m < nRows; m++) {
for (size_t k = 0; k < nColumns; k++) {
m_A(m, k) = m_mix->nAtoms(m_species[m_order[k]], m_element[m]);
}
}
// Do Gaussian elimination
for (m = 0; m < nRows; m++) {
for (size_t m = 0; m < nRows; m++) {
// Check for rows that are zero
bool isZeroRow = true;
for (k = m; k < nColumns; k++) {
for (size_t k = m; k < nColumns; k++) {
if (fabs(m_A(m,k)) > sqrt(Tiny)) {
isZeroRow = false;
break;
@ -301,7 +287,7 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
size_t n = nRows - 1;
bool foundSwapCandidate = false;
for (; n > m; n--) {
for (k = m; k < nColumns; k++) {
for (size_t k = m; k < nColumns; k++) {
if (fabs(m_A(n,k)) > sqrt(Tiny)) {
foundSwapCandidate = true;
break;
@ -313,7 +299,7 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
}
if (m != n) {
// Swap this row with the last non-zero row
for (k = 0; k < nColumns; k++) {
for (size_t k = 0; k < nColumns; k++) {
std::swap(m_A(n,k), m_A(m,k));
}
} else {
@ -334,7 +320,7 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
// satisfies these criteria.
doublereal maxmoles = -999.0;
size_t kmax = 0;
for (k = m+1; k < nColumns; k++) {
for (size_t k = m+1; k < nColumns; k++) {
if (m_A(m,k) != 0.0 && fabs(m_moles[m_order[k]]) > maxmoles) {
kmax = k;
maxmoles = fabs(m_moles[m_order[k]]);
@ -352,8 +338,8 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
}
// scale row m so that the diagonal element is unity
fctr = 1.0/m_A(m,m);
for (k = 0; k < nColumns; k++) {
double fctr = 1.0/m_A(m,m);
for (size_t k = 0; k < nColumns; k++) {
m_A(m,k) *= fctr;
}
@ -361,7 +347,7 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
// * (row m) from row n, so that A(n,m) = 0.
for (size_t n = m+1; n < m_nel; n++) {
fctr = m_A(n,m)/m_A(m,m);
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
m_A(n,k) -= m_A(m,k)*fctr;
}
}
@ -369,11 +355,11 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
// The left m_nel columns of A are now upper-diagonal. Now
// reduce the m_nel columns to diagonal form by back-solving
for (m = std::min(nRows,nColumns)-1; m > 0; m--) {
for (size_t m = std::min(nRows,nColumns)-1; m > 0; m--) {
for (size_t n = m-1; n != npos; n--) {
if (m_A(n,m) != 0.0) {
fctr = m_A(n,m);
for (k = m; k < m_nsp; k++) {
double fctr = m_A(n,m);
for (size_t k = m; k < m_nsp; k++) {
m_A(n,k) -= fctr*m_A(m,k);
}
}
@ -383,11 +369,11 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
// create stoichiometric coefficient matrix.
for (size_t n = 0; n < m_nsp; n++) {
if (n < m_nel) {
for (k = 0; k < nFree(); k++) {
for (size_t k = 0; k < nFree(); k++) {
m_N(n, k) = -m_A(n, k + m_nel);
}
} else {
for (k = 0; k < nFree(); k++) {
for (size_t k = 0; k < nFree(); k++) {
m_N(n, k) = 0.0;
}
m_N(n, n - m_nel) = 1.0;
@ -395,9 +381,9 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
}
// find reactions involving solution phase species
for (j = 0; j < nFree(); j++) {
for (size_t j = 0; j < nFree(); j++) {
m_solnrxn[j] = false;
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
if (m_N(k, j) != 0 && m_mix->solutionSpecies(m_species[m_order[k]])) {
m_solnrxn[j] = true;
}
@ -408,8 +394,7 @@ void MultiPhaseEquil::getComponents(const std::vector<size_t>& order)
void MultiPhaseEquil::unsort(vector_fp& x)
{
m_work2 = x;
size_t k;
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
x[m_order[k]] = m_work2[k];
}
}
@ -417,19 +402,18 @@ void MultiPhaseEquil::unsort(vector_fp& x)
void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN,
int loglevel)
{
size_t k, ik;
if (omega < 0.0) {
throw CanteraError("MultiPhaseEquil::step","negative omega");
}
for (ik = 0; ik < m_nel; ik++) {
k = m_order[ik];
for (size_t ik = 0; ik < m_nel; ik++) {
size_t k = m_order[ik];
m_lastmoles[k] = m_moles[k];
m_moles[k] += omega * deltaN[k];
}
for (ik = m_nel; ik < m_nsp; ik++) {
k = m_order[ik];
for (size_t ik = m_nel; ik < m_nsp; ik++) {
size_t k = m_order[ik];
m_lastmoles[k] = m_moles[k];
if (m_majorsp[k]) {
m_moles[k] += omega * deltaN[k];
@ -444,7 +428,6 @@ void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN,
doublereal MultiPhaseEquil::stepComposition(int loglevel)
{
m_iter++;
size_t ik, k = 0;
doublereal grad0 = computeReactionSteps(m_dxi);
// compute the mole fraction changes.
@ -458,9 +441,9 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
// scale omega to keep the major species non-negative
doublereal FCTR = 0.99;
const doublereal MAJOR_THRESHOLD = 1.0e-12;
doublereal omega = 1.0, omax, omegamax = 1.0;
for (ik = 0; ik < m_nsp; ik++) {
k = m_order[ik];
double omegamax = 1.0;
for (size_t ik = 0; ik < m_nsp; ik++) {
size_t k = m_order[ik];
if (ik < m_nel) {
FCTR = 0.99;
if (m_moles[k] < MAJOR_THRESHOLD) {
@ -477,7 +460,7 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
if (m_moles[k] < MAJOR_THRESHOLD) {
m_force = true;
}
omax = m_moles[k]*FCTR/(fabs(m_work[k]) + Tiny);
double omax = m_moles[k]*FCTR/(fabs(m_work[k]) + Tiny);
if (m_work[k] < 0.0 && omax < omegamax) {
omegamax = omax;
if (omegamax < 1.0e-5) {
@ -490,7 +473,7 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
}
} else {
if (m_work[k] < 0.0 && m_moles[k] > 0.0) {
omax = -m_moles[k]/m_work[k];
double omax = -m_moles[k]/m_work[k];
if (omax < omegamax) {
omegamax = omax;
if (omegamax < 1.0e-5) {
@ -510,14 +493,14 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
doublereal not_mu = 1.0e12;
m_mix->getValidChemPotentials(not_mu, m_mu.data());
doublereal grad1 = 0.0;
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
grad1 += m_work[k] * m_mu[m_species[k]];
}
omega = omegamax;
double omega = omegamax;
if (grad1 > 0.0) {
omega *= fabs(grad0) / (grad1 + fabs(grad0));
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
m_moles[k] = m_lastmoles[k];
}
step(omega, m_work);
@ -527,8 +510,6 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel)
doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi)
{
size_t j, k, ik, kc, ip;
doublereal stoich, nmoles, csum, term1, fctr, rfctr;
vector_fp nu;
doublereal grad = 0.0;
dxi.resize(nFree());
@ -536,24 +517,23 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi)
doublereal not_mu = 1.0e12;
m_mix->getValidChemPotentials(not_mu, m_mu.data());
for (j = 0; j < nFree(); j++) {
for (size_t j = 0; j < nFree(); j++) {
// get stoichiometric vector
getStoichVector(j, nu);
// compute Delta G
doublereal dg_rt = 0.0;
for (k = 0; k < m_nsp; k++) {
for (size_t k = 0; k < m_nsp; k++) {
dg_rt += m_mu[m_species[k]] * nu[k];
}
dg_rt /= (m_temp * GasConstant);
m_deltaG_RT[j] = dg_rt;
fctr = 1.0;
double fctr = 1.0;
// if this is a formation reaction for a single-component phase,
// check whether reaction should be included
ik = j + m_nel;
k = m_order[ik];
size_t k = m_order[j + m_nel];
if (!m_dsoln[k]) {
if (m_moles[k] <= 0.0 && dg_rt > 0.0) {
fctr = 0.0;
@ -564,36 +544,34 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi)
fctr = 1.0;
} else {
// component sum
csum = 0.0;
double csum = 0.0;
for (k = 0; k < m_nel; k++) {
kc = m_order[k];
stoich = nu[kc];
nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny;
csum += stoich*stoich*m_dsoln[kc]/nmoles;
size_t kc = m_order[k];
double nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny;
csum += pow(nu[kc], 2)*m_dsoln[kc]/nmoles;
}
// noncomponent term
kc = m_order[j + m_nel];
nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny;
term1 = m_dsoln[kc]/nmoles;
size_t kc = m_order[j + m_nel];
double nmoles = fabs(m_mix->speciesMoles(m_species[kc])) + Tiny;
double term1 = m_dsoln[kc]/nmoles;
// sum over solution phases
doublereal sum = 0.0, psum;
for (ip = 0; ip < m_mix->nPhases(); ip++) {
doublereal sum = 0.0;
for (size_t ip = 0; ip < m_mix->nPhases(); ip++) {
ThermoPhase& p = m_mix->phase(ip);
if (p.nSpecies() > 1) {
psum = 0.0;
double psum = 0.0;
for (k = 0; k < m_nsp; k++) {
kc = m_species[k];
if (m_mix->speciesPhaseIndex(kc) == ip) {
stoich = nu[k];
psum += stoich * stoich;
psum += pow(nu[k], 2);
}
}
sum -= psum / (fabs(m_mix->phaseMoles(ip)) + Tiny);
}
}
rfctr = term1 + csum + sum;
double rfctr = term1 + csum + sum;
if (fabs(rfctr) < Tiny) {
fctr = 1.0;
} else {
@ -602,8 +580,7 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi)
}
dxi[j] = -fctr*dg_rt;
size_t m;
for (m = 0; m < m_nel; m++) {
for (size_t m = 0; m < m_nel; m++) {
if (m_moles[m_order[m]] <= 0.0 && (m_N(m, j)*dxi[j] < 0.0)) {
dxi[j] = 0.0;
}
@ -628,7 +605,6 @@ void MultiPhaseEquil::computeN()
m_sortindex[k] = moleFractions[k].second;
}
bool ok;
for (size_t m = 0; m < m_nel; m++) {
size_t k = 0;
for (size_t ik = 0; ik < m_nsp; ik++) {
@ -637,7 +613,7 @@ void MultiPhaseEquil::computeN()
break;
}
}
ok = false;
bool ok = false;
for (size_t ij = 0; ij < m_nel; ij++) {
if (k == m_order[ij]) {
ok = true;
@ -685,17 +661,10 @@ double MultiPhaseEquil::phaseMoles(size_t iph) const
void MultiPhaseEquil::reportCSV(const std::string& reportFile)
{
size_t k;
size_t istart;
size_t nSpecies;
double vol = 0.0;
string sName;
FILE* FP = fopen(reportFile.c_str(), "w");
if (!FP) {
throw CanteraError("MultiPhaseEquil::reportCSV", "Failure to open file");
}
double Temp = m_mix->temperature();
double pres = m_mix->pressure();
vector_fp mf(m_nsp_mix, 1.0);
vector_fp fe(m_nsp_mix, 0.0);
vector_fp VolPM;
@ -705,18 +674,18 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
vector_fp mu0;
vector_fp molalities;
vol = 0.0;
double vol = 0.0;
for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) {
istart = m_mix->speciesIndex(0, iphase);
size_t istart = m_mix->speciesIndex(0, iphase);
ThermoPhase& tref = m_mix->phase(iphase);
nSpecies = tref.nSpecies();
size_t nSpecies = tref.nSpecies();
VolPM.resize(nSpecies, 0.0);
tref.getMoleFractions(&mf[istart]);
tref.getPartialMolarVolumes(VolPM.data());
double TMolesPhase = phaseMoles(iphase);
double VolPhaseVolumes = 0.0;
for (k = 0; k < nSpecies; k++) {
for (size_t k = 0; k < nSpecies; k++) {
VolPhaseVolumes += VolPM[k] * mf[istart + k];
}
VolPhaseVolumes *= TMolesPhase;
@ -724,18 +693,18 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
}
fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
" -----------------------------\n");
fprintf(FP,"Temperature = %11.5g kelvin\n", Temp);
fprintf(FP,"Pressure = %11.5g Pascal\n", pres);
fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature());
fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure());
fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) {
istart = m_mix->speciesIndex(0, iphase);
size_t istart = m_mix->speciesIndex(0, iphase);
ThermoPhase& tref = m_mix->phase(iphase);
ThermoPhase* tp = &tref;
tp->getMoleFractions(&mf[istart]);
string phaseName = tref.name();
double TMolesPhase = phaseMoles(iphase);
nSpecies = tref.nSpecies();
size_t nSpecies = tref.nSpecies();
activity.resize(nSpecies, 0.0);
ac.resize(nSpecies, 0.0);
mu0.resize(nSpecies, 0.0);
@ -749,7 +718,7 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
tp->getPartialMolarVolumes(VolPM.data());
tp->getChemPotentials(mu.data());
double VolPhaseVolumes = 0.0;
for (k = 0; k < nSpecies; k++) {
for (size_t k = 0; k < nSpecies; k++) {
VolPhaseVolumes += VolPM[k] * mf[istart + k];
}
VolPhaseVolumes *= TMolesPhase;
@ -768,8 +737,8 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
", , ,"
" (kJ/gmol), (kJ/gmol), (kmol), (m**3/kmol), (m**3)\n");
}
for (k = 0; k < nSpecies; k++) {
sName = tp->speciesName(k);
for (size_t k = 0; k < nSpecies; k++) {
string sName = tp->speciesName(k);
fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e,"
"%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
sName.c_str(),
@ -789,11 +758,11 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
", , ,"
" (kJ/gmol), (kJ/gmol), (kmol), (m**3/kmol), (m**3)\n");
}
for (k = 0; k < nSpecies; k++) {
for (size_t k = 0; k < nSpecies; k++) {
molalities[k] = 0.0;
}
for (k = 0; k < nSpecies; k++) {
sName = tp->speciesName(k);
for (size_t k = 0; k < nSpecies; k++) {
string sName = tp->speciesName(k);
fprintf(FP,"%12s, %11s, %11.3e, %11.3e, %11.3e, %11.3e, %11.3e, "
"%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n",
sName.c_str(),

48
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -68,7 +68,6 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
double P2 = 0.0;
doublereal Tlow = 0.5 * m_mix->minTemp();
doublereal Thigh = 2.0 * m_mix->maxTemp();
doublereal Vnow, Verr;
int printLvlSub = std::max(0, printLvl - 1);
for (int n = 0; n < maxiter; n++) {
double Pnow = m_mix->pressure();
@ -93,7 +92,7 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
break;
}
strt = false;
Vnow = m_mix->volume();
double Vnow = m_mix->volume();
if (n == 0) {
V2 = Vnow;
P2 = Pnow;
@ -107,7 +106,7 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
V1 = Vnow;
}
Verr = fabs((Vtarget - Vnow)/Vtarget);
double Verr = fabs((Vtarget - Vnow)/Vtarget);
if (Verr < err) {
goto done;
}
@ -171,7 +170,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
Thigh = 2.0 * m_mix->maxTemp();
}
doublereal cpb = 1.0, Tnew;
doublereal cpb = 1.0;
doublereal Hlow = Undef;
doublereal Hhigh = Undef;
doublereal Tnow = m_mix->temperature();
@ -211,17 +210,17 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
Hhigh = Hnow;
}
}
double dT, dTa, dTmax, Tnew;
double dT;
if (Hlow != Undef && Hhigh != Undef) {
cpb = (Hhigh - Hlow)/(Thigh - Tlow);
dT = (Htarget - Hnow)/cpb;
dTa = fabs(dT);
dTmax = 0.5*fabs(Thigh - Tlow);
double dTa = fabs(dT);
double dTmax = 0.5*fabs(Thigh - Tlow);
if (dTa > dTmax) {
dT *= dTmax/dTa;
}
} else {
Tnew = sqrt(Tlow*Thigh);
double Tnew = sqrt(Tlow*Thigh);
dT = clip(Tnew - Tnow, -200.0, 200.0);
}
double acpb = std::max(fabs(cpb), 1.0E-6);
@ -244,7 +243,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
}
goto done;
}
Tnew = Tnow + dT;
double Tnew = Tnow + dT;
if (Tnew < 0.0) {
Tnew = 0.5*Tnow;
}
@ -253,7 +252,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
if (!estimateEquil) {
strt = -1;
} else {
Tnew = 0.5*(Tnow + Thigh);
double Tnew = 0.5*(Tnow + Thigh);
if (fabs(Tnew - Tnow) < 1.0) {
Tnew = Tnow + 1.0;
}
@ -286,7 +285,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
Thigh = 2.0 * m_mix->maxTemp();
}
doublereal cpb = 1.0, dT, dTa, dTmax, Tnew;
doublereal cpb = 1.0, dT;
doublereal Slow = Undef;
doublereal Shigh = Undef;
doublereal Tnow = m_mix->temperature();
@ -338,9 +337,9 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
if (Slow != Undef && Shigh != Undef) {
cpb = (Shigh - Slow)/(Thigh - Tlow);
dT = (Starget - Snow)/cpb;
Tnew = Tnow + dT;
dTa = fabs(dT);
dTmax = 0.5*fabs(Thigh - Tlow);
double Tnew = Tnow + dT;
double dTa = fabs(dT);
double dTmax = 0.5*fabs(Thigh - Tlow);
if (Tnew > Thigh || Tnew < Tlow) {
dTmax = 1.5*fabs(Thigh - Tlow);
}
@ -349,7 +348,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
dT *= dTmax/dTa;
}
} else {
Tnew = sqrt(Tlow*Thigh);
double Tnew = sqrt(Tlow*Thigh);
dT = Tnew - Tnow;
}
@ -373,7 +372,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
}
return iSuccess;
}
Tnew = Tnow + dT;
double Tnew = Tnow + dT;
if (Tnew < 0.0) {
Tnew = 0.5*Tnow;
}
@ -382,7 +381,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
if (!estimateEquil) {
strt = -1;
} else {
Tnew = 0.5*(Tnow + Thigh);
double Tnew = 0.5*(Tnow + Thigh);
if (fabs(Tnew - Tnow) < 1.0) {
Tnew = Tnow + 1.0;
}
@ -464,13 +463,11 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
// Check obvious bounds on the temperature and pressure NOTE, we may want to
// do more here with the real bounds given by the ThermoPhase objects.
double T = m_mix->temperature();
if (T <= 0.0) {
if (m_mix->temperature() <= 0.0) {
throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
"Temperature less than zero on input");
}
double pres = m_mix->pressure();
if (pres <= 0.0) {
if (m_mix->pressure() <= 0.0) {
throw CanteraError("vcs_MultiPhaseEquil::equilibrate",
"Pressure less than zero on input");
}
@ -568,7 +565,6 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
{
double vol = 0.0;
size_t nphase = m_vprob.NPhase;
FILE* FP = fopen(reportFile.c_str(), "w");
@ -576,8 +572,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
throw CanteraError("vcs_MultiPhaseEquil::reportCSV",
"Failure to open file");
}
double Temp = m_mix->temperature();
double pres = m_mix->pressure();
vector_fp& mf = m_vprob.mf;
double* fe = &m_vprob.m_gibbsSpecies[0];
vector_fp VolPM;
@ -587,7 +581,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
vector_fp mu0;
vector_fp molalities;
vol = 0.0;
double vol = 0.0;
for (size_t iphase = 0; iphase < nphase; iphase++) {
size_t istart = m_mix->speciesIndex(0, iphase);
ThermoPhase& tref = m_mix->phase(iphase);
@ -607,8 +601,8 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
" -----------------------------\n");
fprintf(FP,"Temperature = %11.5g kelvin\n", Temp);
fprintf(FP,"Pressure = %11.5g Pascal\n", pres);
fprintf(FP,"Temperature = %11.5g kelvin\n", m_mix->temperature());
fprintf(FP,"Pressure = %11.5g Pascal\n", m_mix->pressure());
fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
fprintf(FP,"Number Basis optimizations = %d\n", m_vprob.m_NumBasisOptimizations);
fprintf(FP,"Number VCS iterations = %d\n", m_vprob.m_Iterations);

4
src/equil/vcs_VolPhase.cpp
View file

@ -57,9 +57,7 @@ vcs_VolPhase::vcs_VolPhase(VCS_SOLVE* owningSolverObject) :
vcs_VolPhase::~vcs_VolPhase()
{
for (size_t k = 0; k < m_numSpecies; k++) {
vcs_SpeciesProperties* sp = ListSpeciesPtr[k];
delete sp;
sp = 0;
delete ListSpeciesPtr[k];
}
}

36
src/equil/vcs_phaseStability.cpp
View file

@ -142,7 +142,6 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
// out components which have a pos stoichiometric value with another species
// in the phase.
std::vector< std::vector<size_t> > zeroedPhaseLinkedZeroComponents(m_numPhases);
vector_int linkedPhases;
// The logic below calculates zeroedPhaseLinkedZeroComponents
for (size_t iph = 0; iph < m_numPhases; iph++) {
@ -150,7 +149,6 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
iphList.clear();
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
if (Vphase->exists() < 0) {
linkedPhases.clear();
size_t nsp = Vphase->nSpecies();
for (size_t k = 0; k < nsp; k++) {
size_t kspec = Vphase->spGlobalIndexVCS(k);
@ -447,26 +445,18 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
vector_fp X_est(nsp, 0.0);
vector_fp delFrac(nsp, 0.0);
vector_fp E_phi(nsp, 0.0);
vector_fp fracDelta_new(nsp, 0.0);
vector_fp fracDelta_old(nsp, 0.0);
vector_fp fracDelta_raw(nsp, 0.0);
vector<size_t> creationGlobalRxnNumbers(nsp, npos);
m_deltaGRxn_Deficient = m_deltaGRxn_old;
vector_fp m_feSpecies_Deficient(m_numComponents, 0.0);
doublereal damp = 1.0;
doublereal dampOld = 1.0;
doublereal normUpdate = 1.0;
doublereal normUpdateOld = 1.0;
doublereal sum = 0.0;
doublereal dirProd = 0.0;
doublereal dirProdOld = 0.0;
vector_fp feSpecies_Deficient = m_feSpecies_old;
// get the activity coefficients
Vphase->sendToVCS_ActCoeff(VCS_STATECALC_OLD, &m_actCoeffSpecies_new[0]);
// Get the stored estimate for the composition of the phase if
// it gets created
fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
vector_fp fracDelta_new = Vphase->creationMoleNumbers(creationGlobalRxnNumbers);
std::vector<size_t> componentList;
for (size_t k = 0; k < nsp; k++) {
@ -476,11 +466,8 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
}
}
for (size_t k = 0; k < m_numComponents; k++) {
m_feSpecies_Deficient[k] = m_feSpecies_old[k];
}
normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
damp = 1.0E-2;
double normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
double damp = 1.0E-2;
if (doSuccessiveSubstitution) {
int KP = 0;
@ -500,11 +487,12 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
}
}
bool converged = false;
double dirProd = 0.0;
for (int its = 0; its < 200 && (!converged); its++) {
dampOld = damp;
normUpdateOld = normUpdate;
double dampOld = damp;
double normUpdateOld = normUpdate;
fracDelta_old = fracDelta_new;
dirProdOld = dirProd;
double dirProdOld = dirProd;
// Given a set of fracDelta's, we calculate the fracDelta's
// for the component species, if any
@ -552,10 +540,10 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
size_t kc = componentList[i];
size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
if (X_est[kc] > VCS_DELETE_MINORSPECIES_CUTOFF) {
m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
+ log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
} else {
m_feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
feSpecies_Deficient[kc_spec] = m_feSpecies_old[kc_spec]
+ log(m_actCoeffSpecies_new[kc_spec] * VCS_DELETE_MINORSPECIES_CUTOFF);
}
}
@ -571,14 +559,14 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
}
if (m_stoichCoeffRxnMatrix(kc_spec,irxn) != 0.0) {
m_deltaGRxn_Deficient[irxn] +=
m_stoichCoeffRxnMatrix(kc_spec,irxn) * (m_feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
m_stoichCoeffRxnMatrix(kc_spec,irxn) * (feSpecies_Deficient[kc_spec]- m_feSpecies_old[kc_spec]);
}
}
}
}
// Calculate the E_phi's
sum = 0.0;
double sum = 0.0;
funcPhaseStability = sum_Xcomp - 1.0;
for (size_t k = 0; k < nsp; k++) {
size_t kspec = Vphase->spGlobalIndexVCS(k);

2
src/equil/vcs_prob.cpp
View file

@ -99,11 +99,9 @@ VCS_PROB::~VCS_PROB()
{
for (size_t i = 0; i < nspecies; i++) {
delete SpeciesThermo[i];
SpeciesThermo[i] = 0;
}
for (size_t iph = 0; iph < NPhase; iph++) {
delete VPhaseList[iph];
VPhaseList[iph] = 0;
}
}

7
src/equil/vcs_rxnadj.cpp
View file

@ -358,12 +358,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
}
// Start of the regular processing
double s;
if (m_SSPhase[kspec]) {
s = 0.0;
} else {
s = 1.0 / m_molNumSpecies_old[kspec];
}
double s = (m_SSPhase[kspec]) ? 0.0 : 1.0 / m_molNumSpecies_old[kspec];
for (size_t j = 0; j < m_numComponents; ++j) {
if (!m_SSPhase[j]) {
s += pow(m_stoichCoeffRxnMatrix(j,irxn), 2) / m_molNumSpecies_old[j];

22
src/equil/vcs_setMolesLinProg.cpp
View file

@ -31,18 +31,13 @@ static void printProgress(const vector<string> &spName,
int VCS_SOLVE::vcs_setMolesLinProg()
{
size_t ik, irxn;
double test = -1.0E-10;
if (m_debug_print_lvl >= 2) {
plogf(" --- call setInitialMoles\n");
}
double dg_rt;
int idir;
double nu;
double delta_xi, dxi_min = 1.0e10;
bool redo = true;
double dxi_min = 1.0e10;
int retn;
int iter = 0;
bool abundancesOK = true;
@ -53,7 +48,7 @@ int VCS_SOLVE::vcs_setMolesLinProg()
vector_fp wx(m_numElemConstraints, 0.0);
vector_fp aw(m_numSpeciesTot, 0.0);
for (ik = 0; ik < m_numSpeciesTot; ik++) {
for (size_t ik = 0; ik < m_numSpeciesTot; ik++) {
if (m_speciesUnknownType[ik] != VCS_SPECIES_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[ik] = max(0.0, m_molNumSpecies_old[ik]);
}
@ -63,6 +58,7 @@ int VCS_SOLVE::vcs_setMolesLinProg()
printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies);
}
bool redo = true;
while (redo) {
if (!vcs_elabcheck(0)) {
if (m_debug_print_lvl >= 2) {
@ -98,10 +94,10 @@ int VCS_SOLVE::vcs_setMolesLinProg()
}
// loop over all reactions
for (irxn = 0; irxn < m_numRxnTot; irxn++) {
for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
// dg_rt is the Delta_G / RT value for the reaction
ik = m_numComponents + irxn;
dg_rt = m_SSfeSpecies[ik];
size_t ik = m_numComponents + irxn;
double dg_rt = m_SSfeSpecies[ik];
dxi_min = 1.0e10;
const double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
for (size_t jcomp = 0; jcomp < m_numElemConstraints; jcomp++) {
@ -110,17 +106,17 @@ int VCS_SOLVE::vcs_setMolesLinProg()
// fwd or rev direction.
// idir > 0 implies increasing the current species
// idir < 0 implies decreasing the current species
idir = (dg_rt < 0.0 ? 1 : -1);
int idir = (dg_rt < 0.0 ? 1 : -1);
if (idir < 0) {
dxi_min = m_molNumSpecies_old[ik];
}
for (size_t jcomp = 0; jcomp < m_numComponents; jcomp++) {
nu = sc_irxn[jcomp];
double nu = sc_irxn[jcomp];
// set max change in progress variable by
// non-negativity requirement
if (nu*idir < 0) {
delta_xi = fabs(m_molNumSpecies_old[jcomp]/nu);
double delta_xi = fabs(m_molNumSpecies_old[jcomp]/nu);
// if a component has nearly zero moles, redo
// with a new set of components
if (!redo && delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) {

96
src/equil/vcs_solve_TP.cpp
View file

@ -48,10 +48,9 @@ void VCS_SOLVE::checkDelta1(double* const dsLocal,
int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
{
int stage = MAIN;
int solveFail;
bool allMinorZeroedSpecies = false;
size_t it1 = 0;
size_t npb, iti;
size_t iti;
int rangeErrorFound = 0;
bool giveUpOnElemAbund = false;
int finalElemAbundAttempts = 0;
@ -75,7 +74,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
m_aw.assign(m_numSpeciesTot, 0.0);
m_wx.assign(m_numElemConstraints, 0.0);
solveFail = false;
int solveFail = false;
// Evaluate the elemental composition
vcs_elab();
@ -216,7 +215,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
// reason why we are here is because all of the non-component
// species in the problem have been eliminated for one reason or
// another.
npb = vcs_recheck_deleted();
size_t npb = vcs_recheck_deleted();
// If we haven't found any species that needed adding we are done.
if (npb <= 0) {
@ -232,7 +231,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
}
} else if (stage == RETURN_A) {
// CLEANUP AND RETURN BLOCK
npb = vcs_recheck_deleted();
size_t npb = vcs_recheck_deleted();
// If we haven't found any species that needed adding we are done.
if (npb > 0) {
@ -249,7 +248,7 @@ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit)
} else if (stage == RETURN_B) {
// Add back deleted species in non-zeroed phases. Estimate their
// mole numbers.
npb = vcs_add_all_deleted();
size_t npb = vcs_add_all_deleted();
if (npb > 0) {
iti = 0;
if (m_debug_print_lvl >= 1) {
@ -1427,12 +1426,8 @@ void VCS_SOLVE::solve_tp_elem_abund_check(size_t& iti, int& stage, bool& lec,
double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete,
char* ANOTE) const
{
double dx = 0.0, a;
double w_kspec = m_molNumSpecies_old[kspec];
double molNum_kspec_new;
double wTrial, tmp;
double dg_irxn = m_deltaGRxn_old[irxn];
doublereal s;
size_t iph = m_phaseID[kspec];
*do_delete = false;
@ -1445,20 +1440,22 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete,
sprintf(ANOTE,"minor species alternative calc");
}
if (dg_irxn >= 23.0) {
molNum_kspec_new = w_kspec * 1.0e-10;
double molNum_kspec_new = w_kspec * 1.0e-10;
if (w_kspec < VCS_DELETE_MINORSPECIES_CUTOFF) {
goto L_ZERO_SPECIES;
// delete the species from the current list of active species,
// because its concentration has gotten too small.
*do_delete = true;
return - w_kspec;
}
return molNum_kspec_new - w_kspec;
} else {
if (fabs(dg_irxn) <= m_tolmin2) {
molNum_kspec_new = w_kspec;
return 0.0;
}
}
// get the diagonal of the activity coefficient Jacobian
s = m_np_dLnActCoeffdMolNum(kspec,kspec) / m_tPhaseMoles_old[iph];
double s = m_np_dLnActCoeffdMolNum(kspec,kspec) / m_tPhaseMoles_old[iph];
// We fit it to a power law approximation of the activity coefficient
//
@ -1469,10 +1466,10 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete,
// We then solve the resulting calculation:
//
// gamma * x = gamma_0 * x0 exp (-deltaG/RT);
a = clip(w_kspec * s, -1.0+1e-8, 100.0);
tmp = clip(-dg_irxn / (1.0 + a), -200.0, 200.0);
wTrial = w_kspec * exp(tmp);
molNum_kspec_new = wTrial;
double a = clip(w_kspec * s, -1.0+1e-8, 100.0);
double tmp = clip(-dg_irxn / (1.0 + a), -200.0, 200.0);
double wTrial = w_kspec * exp(tmp);
double molNum_kspec_new = wTrial;
if (wTrial > 100. * w_kspec) {
double molNumMax = 0.0001 * m_tPhaseMoles_old[iph];
@ -1491,39 +1488,34 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete,
}
if ((molNum_kspec_new) < VCS_DELETE_MINORSPECIES_CUTOFF) {
goto L_ZERO_SPECIES;
// delete the species from the current list of active species,
// because its concentration has gotten too small.
*do_delete = true;
return - w_kspec;
}
return molNum_kspec_new - w_kspec;
// Alternate return based for cases where we need to delete the species
// from the current list of active species, because its concentration
// has gotten too small.
L_ZERO_SPECIES:
;
*do_delete = true;
return - w_kspec;
} else {
// Voltage calculation
// Need to check the sign -> This is good for electrons
dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim;
double dx = m_deltaGRxn_old[irxn]/ m_Faraday_dim;
if (ANOTE) {
sprintf(ANOTE,"voltage species alternative calc");
}
return dx;
}
return dx;
}
int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr)
{
size_t irxn = kspec - m_numComponents;
int retn = 1;
double delta = *delta_ptr;
AssertThrowMsg(kspec >= m_numComponents, "VCS_SOLVE::delta_species",
"delete_species() ERROR: called for a component {}", kspec);
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
// Attempt the given dx. If it doesn't work, try to see if a smaller one
// would work,
double dx = delta;
double dx = *delta_ptr;
double* sc_irxn = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
for (size_t j = 0; j < m_numComponents; ++j) {
if (m_molNumSpecies_old[j] > 0.0) {
@ -1742,7 +1734,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
}
}
double dj, dxWant, dxPerm = 0.0, dxPerm2 = 0.0;
double dxPerm = 0.0, dxPerm2 = 0.0;
for (size_t kcomp = 0; kcomp < m_numComponents; ++kcomp) {
if (m_phaseID[kcomp] == iph) {
if (m_debug_print_lvl >= 2) {
@ -1754,7 +1746,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
size_t irxn = kspec - m_numComponents;
if (m_phaseID[kspec] != iph) {
if (m_stoichCoeffRxnMatrix(kcomp,irxn) != 0.0) {
dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(kcomp,irxn);
double dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(kcomp,irxn);
if (dxWant + m_molNumSpecies_old[kspec] < 0.0) {
dxPerm = -m_molNumSpecies_old[kspec];
}
@ -1763,7 +1755,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph)
if (m_phaseID[jcomp] == iph) {
dxPerm = 0.0;
} else {
dj = dxWant * m_stoichCoeffRxnMatrix(jcomp,irxn);
double dj = dxWant * m_stoichCoeffRxnMatrix(jcomp,irxn);
if (dj + m_molNumSpecies_old[kcomp] < 0.0) {
dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix(jcomp,irxn);
}
@ -1938,7 +1930,6 @@ bool VCS_SOLVE::recheck_deleted_phase(const int iphase)
size_t VCS_SOLVE::vcs_add_all_deleted()
{
size_t retn;
if (m_numSpeciesRdc == m_numSpeciesTot) {
return 0;
}
@ -1981,7 +1972,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted()
size_t iph = m_phaseID[kspec];
if (m_tPhaseMoles_old[iph] > 0.0) {
double dx = m_molNumSpecies_new[kspec];
retn = delta_species(kspec, &dx);
size_t retn = delta_species(kspec, &dx);
if (retn == 0) {
if (m_debug_print_lvl) {
plogf(" --- add_deleted(): delta_species() failed for species %s (%d) with mol number %g\n",
@ -2012,7 +2003,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted()
vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot);
vcs_deltag(0, true, VCS_STATECALC_OLD);
retn = 0;
size_t retn = 0;
for (size_t irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) {
size_t kspec = m_indexRxnToSpecies[irxn];
size_t iph = m_phaseID[kspec];
@ -2149,7 +2140,6 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
size_t k;
size_t juse = npos;
size_t jlose = npos;
double* scrxn_ptr;
clockWC tickTock;
if (m_debug_print_lvl >= 2) {
plogf(" ");
@ -2197,7 +2187,6 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
m_numComponents = ncTrial;
*usedZeroedSpecies = false;
vector_int ipiv(ncTrial);
int info;
// Use a temporary work array for the mole numbers, aw[]
std::copy(m_molNumSpecies_old.begin(),
@ -2457,6 +2446,7 @@ L_END_LOOP:
}
// Solve the linear system to calculate the reaction matrix,
// m_stoichCoeffRxnMatrix.
int info;
ct_dgetrf(ncTrial, ncTrial, sm, m_numElemConstraints, &ipiv[0], info);
if (info) {
plogf("vcs_solve_TP ERROR: Error factorizing stoichiometric coefficient matrix\n");
@ -2618,7 +2608,7 @@ L_END_LOOP:
// m_deltaMolNumPhase(iphase,irxn), and the phase participation array,
// PhaseParticipation[irxn][iphase]
for (size_t irxn = 0; irxn < m_numRxnTot; ++irxn) {
scrxn_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
double* scrxn_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
size_t kspec = m_indexRxnToSpecies[irxn];
size_t iph = m_phaseID[kspec];
m_deltaMolNumPhase(iph,irxn) = 1.0;
@ -2702,8 +2692,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const
size_t iph = m_phaseID[kspec];
int irxn = int(kspec) - int(m_numComponents);
vcs_VolPhase* VPhase = m_VolPhaseList[iph];
int phaseExist = VPhase->exists();
int phaseExist = m_VolPhaseList[iph]->exists();
// Treat zeroed out species first
if (m_molNumSpecies_old[kspec] <= 0.0) {
@ -2714,8 +2703,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const
// see if the species has an element which is so low that species will
// always be zero
for (size_t j = 0; j < m_numElemConstraints; ++j) {
int elType = m_elType[j];
if (elType == VCS_ELEM_TYPE_ABSPOS) {
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
double atomComp = m_formulaMatrix(kspec,j);
if (atomComp > 0.0) {
double maxPermissible = m_elemAbundancesGoal[j] / atomComp;
@ -2758,9 +2746,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const
}
}
} else if (negChangeComp < 0.0) {
size_t jph = m_phaseID[j];
vcs_VolPhase* jVPhase = m_VolPhaseList[jph];
if (jVPhase->exists() <= 0) {
if (m_VolPhaseList[m_phaseID[j]]->exists() <= 0) {
if (m_debug_print_lvl >= 2) {
plogf(" --- %s is prevented from popping into existence because"
" a needed component %s is in a zeroed-phase that would be "
@ -3395,7 +3381,6 @@ bool VCS_SOLVE::vcs_evaluate_speciesType()
void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted,
const int vcsState, const bool alterZeroedPhases)
{
int icase = 0;
size_t irxnl = m_numRxnRdc;
if (doDeleted) {
irxnl = m_numRxnTot;
@ -3435,7 +3420,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted,
for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
size_t kspec = irxn + m_numComponents;
if (m_speciesStatus[kspec] != VCS_SPECIES_MINOR) {
icase = 0;
int icase = 0;
deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
for (kspec = 0; kspec < m_numComponents; ++kspec) {
@ -3452,7 +3437,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted,
} else if (L == 0) {
// ALL REACTIONS
for (size_t irxn = 0; irxn < irxnl; ++irxn) {
icase = 0;
int icase = 0;
deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
@ -3471,7 +3456,7 @@ void VCS_SOLVE::vcs_deltag(const int L, const bool doDeleted,
for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
size_t kspec = irxn + m_numComponents;
if (m_speciesStatus[kspec] <= VCS_SPECIES_MINOR) {
icase = 0;
int icase = 0;
deltaGRxn[irxn] = feSpecies[m_indexRxnToSpecies[irxn]];
double* dtmp_ptr = m_stoichCoeffRxnMatrix.ptrColumn(irxn);
for (kspec = 0; kspec < m_numComponents; ++kspec) {
@ -3719,10 +3704,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted,
throw CanteraError("VCS_SOLVE::vcs_deltag_Phase", "bad stateCalc");
}
size_t irxnl = m_numRxnRdc;
if (doDeleted) {
irxnl = m_numRxnTot;
}
size_t irxnl = (doDeleted) ? m_numRxnTot : m_numRxnRdc;
vcs_VolPhase* vPhase = m_VolPhaseList[iphase];
if (m_debug_print_lvl >= 2) {
@ -3900,19 +3882,17 @@ double VCS_SOLVE::vcs_birthGuess(const int kspec)
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
return dx;
}
double w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF;
// Check to make sure that species is zero in the solution vector
// If it isn't, we don't know what's happening
AssertThrowMsg(m_molNumSpecies_old[kspec] == 0.0,
"VCS_SOLVE::vcs_birthGuess", "we shouldn't be here");
int ss = m_SSPhase[kspec];
if (!ss) {
if (!m_SSPhase[kspec]) {
// Logic to handle species in multiple species phases. We cap the moles
// here at 1.0E-15 kmol.
bool soldel_ret;
double dxm = vcs_minor_alt_calc(kspec, irxn, &soldel_ret);
dx = std::min(w_kspec + dxm, 1e-15);
dx = std::min(VCS_DELETE_MINORSPECIES_CUTOFF + dxm, 1e-15);
} else {
// Logic to handle single species phases. There is no real way to
// estimate the moles. So we set it to a small number.

10
src/equil/vcs_solve_phaseStability.cpp
View file

@ -18,7 +18,6 @@ int VCS_SOLVE::vcs_PS(VCS_PROB* vprob, int iphase, int printLvl, double& feStabl
{
// ifunc determines the problem type
int ifunc = 0;
int iStab = 0;
// This function is called to create the private data using the public data.
size_t nspecies0 = vprob->nspecies + 10;
@ -91,7 +90,7 @@ int VCS_SOLVE::vcs_PS(VCS_PROB* vprob, int iphase, int printLvl, double& feStabl
// 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.)
iStab = vcs_solve_phaseStability(iphase, ifunc, feStable, printLvl);
int iStab = vcs_solve_phaseStability(iphase, ifunc, feStable, printLvl);
// Redimensionalize the free energies using the reverse of vcs_nondim to add
// back units.
@ -109,7 +108,6 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc,
{
double test = -1.0E-10;
bool usedZeroedSpecies;
int iStab = 0;
vector_fp sm(m_numElemConstraints*m_numElemConstraints, 0.0);
vector_fp ss(m_numElemConstraints, 0.0);
@ -133,12 +131,10 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc,
m_deltaGRxn_Deficient = m_deltaGRxn_old;
funcVal = vcs_phaseStabilityTest(iph);
if (funcVal > 0.0) {
iStab = 1;
return 1;
} else {
iStab = 0;
return 0;
}
return iStab;
}
}

18
src/equil/vcs_species_thermo.cpp
View file

@ -94,26 +94,21 @@ double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin,
throw CanteraError("VCS_SPECIES_THERMO::GStar_R_calc",
"Possible inconsistency");
}
size_t kspec = IndexSpeciesPhase;
OwningPhase->setState_TP(TKelvin, pres);
fe = OwningPhase->GStar_calc_one(kspec);
fe = OwningPhase->GStar_calc_one(IndexSpeciesPhase);
double R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
fe /= R;
return fe;
return fe / R;
}
double VCS_SPECIES_THERMO::VolStar_calc(size_t kglob, double TKelvin,
double presPA)
{
double vol;
if (m_VCS_UnitsFormat != VCS_UNITS_MKS) {
throw CanteraError("VCS_SPECIES_THERMO::VolStar_calc",
"Possible inconsistency");
}
size_t kspec = IndexSpeciesPhase;
OwningPhase->setState_TP(TKelvin, presPA);
vol = OwningPhase->VolStar_calc_one(kspec);
return vol;
return OwningPhase->VolStar_calc_one(IndexSpeciesPhase);
}
double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin)
@ -128,9 +123,8 @@ double VCS_SPECIES_THERMO::G0_R_calc(size_t kglob, double TKelvin)
throw CanteraError("VCS_SPECIES_THERMO::G0_R_calc",
"Possible inconsistency");
}
size_t kspec = IndexSpeciesPhase;
OwningPhase->setState_T(TKelvin);
double fe = OwningPhase->G0_calc_one(kspec);
double fe = OwningPhase->G0_calc_one(IndexSpeciesPhase);
double R = vcsUtil_gasConstant(m_VCS_UnitsFormat);
fe /= R;
SS0_feSave = fe;
@ -144,9 +138,7 @@ double VCS_SPECIES_THERMO::eval_ac(size_t kglob)
// they are, then the currPhAC[] boolean may be used to reduce repeated
// work. Just set currPhAC[iph], when the activity coefficients for all
// species in the phase are reevaluated.
size_t kspec = IndexSpeciesPhase;
double ac = OwningPhase->AC_calc_one(kspec);
return ac;
return OwningPhase->AC_calc_one(IndexSpeciesPhase);
}
}

7
src/equil/vcs_util.cpp
View file

@ -20,17 +20,16 @@ using namespace std;
namespace Cantera
{
double vcs_l2norm(const vector_fp vec)
double vcs_l2norm(const vector_fp& vec)
{
size_t len = vec.size();
if (len == 0) {
if (vec.empty()) {
return 0.0;
}
double sum = 0.0;
for (const auto& val : vec) {
sum += val * val;
}
return std::sqrt(sum / len);
return std::sqrt(sum / vec.size());
}
size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n)