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

[VCS] Use local variable declarations where possible

Browse source
This commit is contained in:
Ray Speth 2014-05-23 22:47:03 +00:00
parent 217b396ec5
commit 4dcdc01fb8
18 changed files with 586 additions and 824 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

131
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -60,19 +60,15 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
{
// doublereal dt = 1.0e3;
doublereal Vtarget = m_mix->volume();
doublereal dVdP;
if ((XY != TV) && (XY != HV) && (XY != UV) && (XY != SV)) {
throw CanteraError("vcs_MultiPhaseEquil::equilibrate_TV",
"Wrong XY flag:" + int2str(XY));
}
int maxiter = 100;
int iSuccess = 0;
int innerXY;
double Pnow;
if (XY == TV) {
m_mix->setTemperature(xtarget);
}
double Pnew;
int strt = estimateEquil;
double P1 = 0.0;
double V1 = 0.0;
@ -83,24 +79,21 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
doublereal Vnow, Verr;
int printLvlSub = std::max(0, printLvl - 1);
for (int n = 0; n < maxiter; n++) {
Pnow = m_mix->pressure();
double Pnow = m_mix->pressure();
switch (XY) {
case TV:
iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
break;
case HV:
innerXY = HP;
iSuccess = equilibrate_HP(xtarget, innerXY, Tlow, Thigh, strt,
iSuccess = equilibrate_HP(xtarget, HP, Tlow, Thigh, strt,
printLvlSub, err, maxsteps, loglevel);
break;
case UV:
innerXY = UP;
iSuccess = equilibrate_HP(xtarget, innerXY, Tlow, Thigh, strt,
iSuccess = equilibrate_HP(xtarget, UP, Tlow, Thigh, strt,
printLvlSub, err, maxsteps, loglevel);
break;
case SV:
innerXY = SP;
iSuccess = equilibrate_SP(xtarget, Tlow, Thigh, strt,
printLvlSub, err, maxsteps, loglevel);
break;
@ -127,9 +120,10 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
if (Verr < err) {
goto done;
}
double Pnew;
// find dV/dP
if (n > 1) {
dVdP = (V2 - V1) / (P2 - P1);
double dVdP = (V2 - V1) / (P2 - P1);
if (dVdP == 0.0) {
throw CanteraError("vcs_MultiPhase::equilibrate_TV",
"dVdP == 0.0");
@ -145,7 +139,7 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget,
} else {
m_mix->setPressure(Pnow*1.01);
dVdP = (m_mix->volume() - Vnow)/(0.01*Pnow);
double dVdP = (m_mix->volume() - Vnow)/(0.01*Pnow);
Pnew = Pnow + 0.5*(Vtarget - Vnow)/dVdP;
if (Pnew < 0.5* Pnow) {
Pnew = 0.5 * Pnow;
@ -188,11 +182,9 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
Thigh = 2.0 * m_mix->maxTemp();
}
doublereal cpb = 1.0, dT, dTa, dTmax, Tnew;
doublereal Hnow;
doublereal cpb = 1.0, Tnew;
doublereal Hlow = Undef;
doublereal Hhigh = Undef;
doublereal Herr, HConvErr;
doublereal Tnow = m_mix->temperature();
int printLvlSub = std::max(printLvl - 1, 0);
@ -203,11 +195,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
Tnow = m_mix->temperature();
iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
strt = 0;
if (XY == UP) {
Hnow = m_mix->IntEnergy();
} else {
Hnow = m_mix->enthalpy();
}
double Hnow = (XY == UP) ? m_mix->IntEnergy() : m_mix->enthalpy();
double pmoles[10];
pmoles[0] = m_mix->phaseMoles(0);
double Tmoles = pmoles[0];
@ -237,6 +225,7 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
Hhigh = Hnow;
}
}
double dT, dTa, dTmax, Tnew;
if (Hlow != Undef && Hhigh != Undef) {
cpb = (Hhigh - Hlow)/(Thigh - Tlow);
dT = (Htarget - Hnow)/cpb;
@ -257,8 +246,8 @@ int vcs_MultiPhaseEquil::equilibrate_HP(doublereal Htarget,
}
double acpb = std::max(fabs(cpb), 1.0E-6);
double denom = std::max(fabs(Htarget), acpb);
Herr = Htarget - Hnow;
HConvErr = fabs((Herr)/denom);
double Herr = Htarget - Hnow;
double HConvErr = fabs((Herr)/denom);
if (printLvl > 0) {
plogf(" equilibrate_HP: It = %d, Tcurr = %g Hcurr = %g, Htarget = %g\n",
n, Tnow, Hnow, Htarget);
@ -320,10 +309,8 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
}
doublereal cpb = 1.0, dT, dTa, dTmax, Tnew;
doublereal Snow;
doublereal Slow = Undef;
doublereal Shigh = Undef;
doublereal Serr, SConvErr;
doublereal Tnow = m_mix->temperature();
if (Tnow < Tlow) {
Tlow = Tnow;
@ -341,7 +328,7 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
Tnow = m_mix->temperature();
int iSuccess = equilibrate_TP(strt, printLvlSub, err, maxsteps, loglevel);
strt = 0;
Snow = m_mix->entropy();
double Snow = m_mix->entropy();
double pmoles[10];
pmoles[0] = m_mix->phaseMoles(0);
double Tmoles = pmoles[0];
@ -396,8 +383,8 @@ int vcs_MultiPhaseEquil::equilibrate_SP(doublereal Starget,
double acpb = std::max(fabs(cpb), 1.0E-6);
double denom = std::max(fabs(Starget), acpb);
Serr = Starget - Snow;
SConvErr = fabs((Serr)/denom);
double Serr = Starget - Snow;
double SConvErr = fabs((Serr)/denom);
if (printLvl > 0) {
plogf(" equilibrate_SP: It = %d, Tcurr = %g Scurr = %g, Starget = %g\n",
n, Tnow, Snow, Starget);
@ -441,10 +428,9 @@ int vcs_MultiPhaseEquil::equilibrate(int XY, int estimateEquil,
int printLvl, doublereal err,
int maxsteps, int loglevel)
{
int iSuccess;
doublereal xtarget;
if (XY == TP) {
iSuccess = equilibrate_TP(estimateEquil, printLvl, err, maxsteps, loglevel);
return equilibrate_TP(estimateEquil, printLvl, err, maxsteps, loglevel);
} else if (XY == HP || XY == UP) {
if (XY == HP) {
xtarget = m_mix->enthalpy();
@ -453,36 +439,34 @@ int vcs_MultiPhaseEquil::equilibrate(int XY, int estimateEquil,
}
double Tlow = 0.5 * m_mix->minTemp();
double Thigh = 2.0 * m_mix->maxTemp();
iSuccess = equilibrate_HP(xtarget, XY, Tlow, Thigh,
estimateEquil, printLvl, err, maxsteps, loglevel);
return equilibrate_HP(xtarget, XY, Tlow, Thigh,
estimateEquil, printLvl, err, maxsteps, loglevel);
} else if (XY == SP) {
xtarget = m_mix->entropy();
double Tlow = 0.5 * m_mix->minTemp();
double Thigh = 2.0 * m_mix->maxTemp();
iSuccess = equilibrate_SP(xtarget, Tlow, Thigh,
estimateEquil, printLvl, err, maxsteps, loglevel);
return equilibrate_SP(xtarget, Tlow, Thigh,
estimateEquil, printLvl, err, maxsteps, loglevel);
} else if (XY == TV) {
xtarget = m_mix->temperature();
iSuccess = equilibrate_TV(XY, xtarget,
estimateEquil, printLvl, err, maxsteps, loglevel);
return equilibrate_TV(XY, xtarget,
estimateEquil, printLvl, err, maxsteps, loglevel);
} else if (XY == HV) {
xtarget = m_mix->enthalpy();
iSuccess = equilibrate_TV(XY, xtarget,
estimateEquil, printLvl, err, maxsteps, loglevel);
return equilibrate_TV(XY, xtarget,
estimateEquil, printLvl, err, maxsteps, loglevel);
} else if (XY == UV) {
xtarget = m_mix->IntEnergy();
iSuccess = equilibrate_TV(XY, xtarget,
estimateEquil, printLvl, err, maxsteps, loglevel);
return equilibrate_TV(XY, xtarget,
estimateEquil, printLvl, err, maxsteps, loglevel);
} else if (XY == SV) {
xtarget = m_mix->entropy();
iSuccess = equilibrate_TV(XY, xtarget, estimateEquil,
return equilibrate_TV(XY, xtarget, estimateEquil,
printLvl, err, maxsteps, loglevel);
} else {
throw CanteraError(" vcs_MultiPhaseEquil::equilibrate",
"Unsupported Option");
}
return iSuccess;
}
int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
@ -631,12 +615,7 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
{
size_t k;
size_t istart;
size_t nSpecies;
double vol = 0.0;
string sName;
size_t nphase = m_vprob.NPhase;
FILE* FP = fopen(reportFile.c_str(), "w");
@ -657,19 +636,18 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
std::vector<double> mu0;
std::vector<double> molalities;
vol = 0.0;
for (size_t iphase = 0; iphase < nphase; iphase++) {
istart = m_mix->speciesIndex(0, iphase);
size_t istart = m_mix->speciesIndex(0, iphase);
Cantera::ThermoPhase& tref = m_mix->phase(iphase);
nSpecies = tref.nSpecies();
size_t nSpecies = tref.nSpecies();
VolPM.resize(nSpecies, 0.0);
tref.getPartialMolarVolumes(VCS_DATA_PTR(VolPM));
vcs_VolPhase* volP = m_vprob.VPhaseList[iphase];
double TMolesPhase = volP->totalMoles();
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;
@ -685,14 +663,13 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
fprintf(FP,"Number VCS iterations = %d\n", m_vprob.m_Iterations);
for (size_t iphase = 0; iphase < nphase; iphase++) {
istart = m_mix->speciesIndex(0, iphase);
size_t istart = m_mix->speciesIndex(0, iphase);
Cantera::ThermoPhase& tref = m_mix->phase(iphase);
Cantera::ThermoPhase* tp = &tref;
string phaseName = tref.name();
vcs_VolPhase* volP = m_vprob.VPhaseList[iphase];
double TMolesPhase = volP->totalMoles();
//AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase));
nSpecies = tref.nSpecies();
size_t nSpecies = tref.nSpecies();
activity.resize(nSpecies, 0.0);
ac.resize(nSpecies, 0.0);
@ -701,15 +678,15 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
VolPM.resize(nSpecies, 0.0);
molalities.resize(nSpecies, 0.0);
int actConvention = tp->activityConvention();
tp->getActivities(VCS_DATA_PTR(activity));
tp->getActivityCoefficients(VCS_DATA_PTR(ac));
tp->getStandardChemPotentials(VCS_DATA_PTR(mu0));
int actConvention = tref.activityConvention();
tref.getActivities(VCS_DATA_PTR(activity));
tref.getActivityCoefficients(VCS_DATA_PTR(ac));
tref.getStandardChemPotentials(VCS_DATA_PTR(mu0));
tp->getPartialMolarVolumes(VCS_DATA_PTR(VolPM));
tp->getChemPotentials(VCS_DATA_PTR(mu));
tref.getPartialMolarVolumes(VCS_DATA_PTR(VolPM));
tref.getChemPotentials(VCS_DATA_PTR(mu));
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;
@ -717,9 +694,9 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
if (actConvention == 1) {
MolalityVPSSTP* mTP = static_cast<MolalityVPSSTP*>(tp);
MolalityVPSSTP* mTP = static_cast<MolalityVPSSTP*>(&tref);
mTP->getMolalities(VCS_DATA_PTR(molalities));
tp->getChemPotentials(VCS_DATA_PTR(mu));
tref.getChemPotentials(VCS_DATA_PTR(mu));
if (iphase == 0) {
fprintf(FP," Name, Phase, PhaseMoles, Mole_Fract, "
@ -730,8 +707,8 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
" , , ,"
" (J/kmol), (J/kmol), (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++) {
std::string sName = tref.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(),
@ -752,11 +729,11 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
" , , ,"
" (J/kmol), (J/kmol), (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++) {
std::string sName = tref.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(),
@ -772,8 +749,8 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile)
/*
* Check consistency: These should be equal
*/
tp->getChemPotentials(fe+istart);
for (k = 0; k < nSpecies; k++) {
tref.getChemPotentials(fe+istart);
for (size_t k = 0; k < nSpecies; k++) {
if (!vcs_doubleEqual(fe[istart+k], mu[k])) {
fprintf(FP,"ERROR: incompatibility!\n");
fclose(FP);
@ -828,9 +805,6 @@ int vcs_Cantera_to_vprob(Cantera::MultiPhase* mphase,
vprob->Vol = mphase->volume();
vprob->Title = "MultiPhase Object";
Cantera::ThermoPhase* tPhase = 0;
bool gasPhase;
int printLvl = vprob->m_printLvl;
/*
@ -842,7 +816,7 @@ int vcs_Cantera_to_vprob(Cantera::MultiPhase* mphase,
/*
* Get the thermophase object - assume volume phase
*/
tPhase = &(mphase->phase(iphase));
Cantera::ThermoPhase* tPhase = &(mphase->phase(iphase));
size_t nelem = tPhase->nElements();
/*
@ -850,11 +824,7 @@ int vcs_Cantera_to_vprob(Cantera::MultiPhase* mphase,
* current phase.
*/
int eos = tPhase->eosType();
if (eos == cIdealGas) {
gasPhase = true;
} else {
gasPhase = false;
}
bool gasPhase = (eos == cIdealGas);
/*
* Find out the number of species in the phase
@ -1226,10 +1196,9 @@ int vcs_Cantera_update_vprob(Cantera::MultiPhase* mphase,
vprob->T = mphase->temperature();
vprob->PresPA = mphase->pressure();
vprob->Vol = mphase->volume();
Cantera::ThermoPhase* tPhase = 0;
for (size_t iphase = 0; iphase < totNumPhases; iphase++) {
tPhase = &(mphase->phase(iphase));
Cantera::ThermoPhase* tPhase = &(mphase->phase(iphase));
vcs_VolPhase* volPhase = vprob->VPhaseList[iphase];
/*
* Set the electric potential of the volume phase from the

5
src/equil/vcs_TP.cpp
View file

@ -8,7 +8,6 @@ namespace VCSnonideal
{
int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg)
{
int retn, iconv;
/*
* Store the temperature and pressure in the private global variables
*/
@ -18,7 +17,7 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg
* Evaluate the standard state free energies
* at the current temperatures and pressures.
*/
iconv = vcs_evalSS_TP(ipr, ip1, m_temperature, pres_arg);
int iconv = vcs_evalSS_TP(ipr, ip1, m_temperature, pres_arg);
/*
* Prepare the problem data:
@ -35,7 +34,7 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg
* If so, go get one. If not, then
*/
if (m_doEstimateEquil) {
retn = vcs_inest_TP();
int retn = vcs_inest_TP();
if (retn != VCS_SUCCESS) {
plogf("vcs_inest_TP returned a failure flag\n");
}

40
src/equil/vcs_VolPhase.cpp
View file

@ -481,8 +481,6 @@ void vcs_VolPhase::setMoleFractionsState(const double totalMoles,
void vcs_VolPhase::setMolesFromVCS(const int stateCalc,
const double* molesSpeciesVCS)
{
size_t kglob;
double tmp;
v_totalMoles = m_totalMolesInert;
if (molesSpeciesVCS == 0) {
@ -524,15 +522,15 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc,
for (size_t k = 0; k < m_numSpecies; k++) {
if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
kglob = IndSpecies[k];
size_t kglob = IndSpecies[k];
v_totalMoles += std::max(0.0, molesSpeciesVCS[kglob]);
}
}
if (v_totalMoles > 0.0) {
for (size_t k = 0; k < m_numSpecies; k++) {
if (m_speciesUnknownType[k] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
kglob = IndSpecies[k];
tmp = std::max(0.0, molesSpeciesVCS[kglob]);
size_t kglob = IndSpecies[k];
double tmp = std::max(0.0, molesSpeciesVCS[kglob]);
Xmol_[k] = tmp / v_totalMoles;
}
}
@ -551,7 +549,7 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc,
* in the equation system
*/
if (m_phiVarIndex != npos) {
kglob = IndSpecies[m_phiVarIndex];
size_t kglob = IndSpecies[m_phiVarIndex];
if (m_numSpecies == 1) {
Xmol_[m_phiVarIndex] = 1.0;
} else {
@ -1200,9 +1198,6 @@ static bool chargeNeutralityElement(const Cantera::ThermoPhase* const tPhase)
size_t vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase* const tPhase)
{
size_t e, k, eT;
std::string ename;
size_t eFound = npos;
size_t nebase = tPhase->nElements();
size_t ne = nebase;
size_t ns = tPhase->nSpecies();
@ -1227,6 +1222,7 @@ size_t vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase* const tPhase
m_elementType[ChargeNeutralityElement] = VCS_ELEM_TYPE_CHARGENEUTRALITY;
}
size_t eFound = npos;
if (hasChargedSpecies(tPhase)) {
if (cne) {
/*
@ -1238,18 +1234,16 @@ size_t vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase* const tPhase
* by toggling the ElActive variable. If we find we need it
* later, we will retoggle ElActive to true.
*/
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
if (ename == "E") {
for (size_t eT = 0; eT < nebase; eT++) {
if (tPhase->elementName(eT) == "E") {
eFound = eT;
m_elementActive[eT] = 0;
m_elementType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
}
}
} else {
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
if (ename == "E") {
for (size_t eT = 0; eT < nebase; eT++) {
if (tPhase->elementName(eT) == "E") {
eFound = eT;
m_elementType[eT] = VCS_ELEM_TYPE_ELECTRONCHARGE;
}
@ -1273,10 +1267,9 @@ size_t vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase* const tPhase
elemResize(ne);
e = 0;
for (eT = 0; eT < nebase; eT++) {
ename = tPhase->elementName(eT);
m_elementNames[e] = ename;
size_t e = 0;
for (size_t eT = 0; eT < nebase; eT++) {
m_elementNames[e] = tPhase->elementName(eT);
m_elementType[e] = tPhase->elementType(eT);
e++;
}
@ -1288,15 +1281,14 @@ size_t vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase* const tPhase
sss << "phase" << VP_ID_;
pname = sss.str();
}
ename = "cn_" + pname;
e = ChargeNeutralityElement;
m_elementNames[e] = ename;
m_elementNames[e] = "cn_" + pname;
}
double* const* const fm = m_formulaMatrix.baseDataAddr();
for (k = 0; k < ns; k++) {
for (size_t k = 0; k < ns; k++) {
e = 0;
for (eT = 0; eT < nebase; eT++) {
for (size_t eT = 0; eT < nebase; eT++) {
fm[e][k] = tPhase->nAtoms(k, eT);
e++;
}
@ -1306,7 +1298,7 @@ size_t vcs_VolPhase::transferElementsFM(const Cantera::ThermoPhase* const tPhase
}
if (cne) {
for (k = 0; k < ns; k++) {
for (size_t k = 0; k < ns; k++) {
fm[ChargeNeutralityElement][k] = tPhase->charge(k);
}
}

44
src/equil/vcs_elem.cpp
View file

@ -30,9 +30,6 @@ void VCS_SOLVE::vcs_elab()
bool VCS_SOLVE::vcs_elabcheck(int ibound)
{
size_t top = m_numComponents;
double eval, scale;
int numNonZero;
bool multisign = false;
if (ibound) {
top = m_numElemConstraints;
}
@ -51,23 +48,21 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound)
"Problem with charge neutrality condition");
}
if (m_elemAbundancesGoal[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) {
scale = VCS_DELETE_MINORSPECIES_CUTOFF;
double scale = VCS_DELETE_MINORSPECIES_CUTOFF;
/*
* Find out if the constraint is a multisign constraint.
* If it is, then we have to worry about roundoff error
* in the addition of terms. We are limited to 13
* digits of finite arithmetic accuracy.
*/
numNonZero = 0;
multisign = false;
bool multisign = false;
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
eval = m_formulaMatrix[i][kspec];
double eval = m_formulaMatrix[i][kspec];
if (eval < 0.0) {
multisign = true;
}
if (eval != 0.0) {
scale = std::max(scale, fabs(eval * m_molNumSpecies_old[kspec]));
numNonZero++;
}
}
if (multisign) {
@ -112,12 +107,9 @@ void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase)
int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
{
int retn = 0, its;
bool goodSpec;
double xx, par, saveDir, dir;
int retn = 0;
#ifdef DEBUG_MODE
double l2before = 0.0, l2after = 0.0;
std::vector<double> ga_save(m_numElemConstraints, 0.0);
vcs_dcopy(VCS_DATA_PTR(ga_save), VCS_DATA_PTR(m_elemAbundances), m_numElemConstraints);
if (m_debug_print_lvl >= 2) {
@ -131,7 +123,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
for (size_t i = 0; i < m_numElemConstraints; ++i) {
x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
}
l2before = 0.0;
double l2before = 0.0;
for (size_t i = 0; i < m_numElemConstraints; ++i) {
l2before += x[i] * x[i];
}
@ -144,12 +136,10 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
* formula matrix, and no other species have zero values either.
*
*/
int numNonZero = 0;
bool changed = false;
bool multisign = false;
for (size_t i = 0; i < m_numElemConstraints; ++i) {
numNonZero = 0;
multisign = false;
int numNonZero = 0;
bool multisign = false;
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
double eval = m_formulaMatrix[i][kspec];
@ -286,10 +276,10 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
/*
* Now apply the new direction without creating negative species.
*/
par = 0.5;
double par = 0.5;
for (size_t i = 0; i < m_numComponents; ++i) {
if (m_molNumSpecies_old[i] > 0.0) {
xx = -x[i] / m_molNumSpecies_old[i];
double xx = -x[i] / m_molNumSpecies_old[i];
if (par < xx) {
par = xx;
}
@ -350,10 +340,10 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
continue;
}
saveDir = 0.0;
goodSpec = true;
double saveDir = 0.0;
bool goodSpec = true;
for (size_t i = 0; i < m_numComponents; ++i) {
dir = m_formulaMatrix[i][kspec] * (m_elemAbundancesGoal[i] - m_elemAbundances[i]);
double dir = m_formulaMatrix[i][kspec] * (m_elemAbundancesGoal[i] - m_elemAbundances[i]);
if (fabs(dir) > 1.0E-10) {
if (dir > 0.0) {
if (saveDir < 0.0) {
@ -375,8 +365,8 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
}
}
if (goodSpec) {
its = 0;
xx = 0.0;
int its = 0;
double xx = 0.0;
for (size_t i = 0; i < m_numComponents; ++i) {
if (m_formulaMatrix[i][kspec] != 0.0) {
xx += (m_elemAbundancesGoal[i] - m_elemAbundances[i]) / m_formulaMatrix[i][kspec];
@ -500,15 +490,15 @@ L_CLEANUP:
;
vcs_tmoles();
#ifdef DEBUG_MODE
l2after = 0.0;
for (int i = 0; i < m_numElemConstraints; ++i) {
double l2after = 0.0;
for (size_t i = 0; i < m_numElemConstraints; ++i) {
l2after += SQUARE(m_elemAbundances[i] - m_elemAbundancesGoal[i]);
}
l2after = sqrt(l2after/m_numElemConstraints);
if (m_debug_print_lvl >= 2) {
plogf(" --- Elem_Abund: Correct Initial "
" Final\n");
for (int i = 0; i < m_numElemConstraints; ++i) {
for (size_t i = 0; i < m_numElemConstraints; ++i) {
plogf(" --- ");
plogf("%-2.2s", m_elementName[i].c_str());
plogf(" %20.12E %20.12E %20.12E\n", m_elemAbundancesGoal[i], ga_save[i], m_elemAbundances[i]);

35
src/equil/vcs_elem_rearrange.cpp
View file

@ -19,14 +19,11 @@ namespace VCSnonideal
int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
double* const sm, double* const ss)
{
size_t j, k, l, i, jl, ml, jr, ielem;
bool lindep;
size_t ncomponents = m_numComponents;
double test = -1.0E10;
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
plogf(" ");
for (i=0; i<77; i++) {
for (size_t i=0; i<77; i++) {
plogf("-");
}
plogf("\n");
@ -41,10 +38,11 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
* Use a temporary work array for the element numbers
* Also make sure the value of test is unique.
*/
lindep = false;
bool lindep = false;
double test = -1.0E10;
do {
lindep = false;
for (i = 0; i < m_numElemConstraints; ++i) {
for (size_t i = 0; i < m_numElemConstraints; ++i) {
test -= 1.0;
aw[i] = m_elemAbundancesGoal[i];
if (test == aw[i]) {
@ -57,9 +55,10 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
* Top of a loop of some sort based on the index JR. JR is the
* current number independent elements found.
*/
jr = npos;
size_t jr = npos;
do {
++jr;
size_t k;
/*
* Top of another loop point based on finding a linearly
* independent species
@ -70,7 +69,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
* for the largest remaining species. Return its identity in K.
*/
k = m_numElemConstraints;
for (ielem = jr; ielem < m_numElemConstraints; ielem++) {
for (size_t ielem = jr; ielem < m_numElemConstraints; ielem++) {
if (m_elementActive[ielem]) {
if (aw[ielem] != test) {
k = ielem;
@ -98,13 +97,13 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
* 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
* from the current component.
*/
for (j = 0; j < ncomponents; ++j) {
for (size_t j = 0; j < ncomponents; ++j) {
sm[j + jr*ncomponents] = m_formulaMatrix[k][j];
}
if (jl > 0) {
@ -114,9 +113,9 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
* (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];
@ -125,8 +124,8 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
* Now make the new column, (*,JR), orthogonal to the
* previous columns
*/
for (j = 0; j < jl; ++j) {
for (l = 0; l < ncomponents; ++l) {
for (size_t j = 0; j < jl; ++j) {
for (size_t l = 0; l < ncomponents; ++l) {
sm[l + jr*ncomponents] -= ss[j] * sm[l + j*ncomponents];
}
}
@ -137,7 +136,7 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa,
* 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) {
sa[jr] += SQUARE(sm[ml + jr*ncomponents]);
}
/* **************************************************** */
@ -181,8 +180,6 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
if (ipos == jpos) {
return;
}
size_t j;
vcs_VolPhase* volPhase;
#ifdef DEBUG_MODE
if (ipos > (m_numElemConstraints - 1) ||
jpos > (m_numElemConstraints - 1)) {
@ -197,7 +194,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
* to reflect the switch in the element positions.
*/
for (size_t iph = 0; iph < m_numPhases; iph++) {
volPhase = m_VolPhaseList[iph];
vcs_VolPhase* volPhase = m_VolPhaseList[iph];
for (size_t e = 0; e < volPhase->nElemConstraints(); e++) {
if (volPhase->elemGlobalIndex(e) == ipos) {
volPhase->setElemGlobalIndex(e, jpos);
@ -212,7 +209,7 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos)
std::swap(m_elementMapIndex[ipos], m_elementMapIndex[jpos]);
std::swap(m_elType[ipos], m_elType[jpos]);
std::swap(m_elementActive[ipos], m_elementActive[jpos]);
for (j = 0; j < m_numSpeciesTot; ++j) {
for (size_t j = 0; j < m_numSpeciesTot; ++j) {
std::swap(m_formulaMatrix[ipos][j], m_formulaMatrix[jpos][j]);
}
std::swap(m_elementName[ipos], m_elementName[jpos]);

3
src/equil/vcs_equilibrate.cpp
View file

@ -36,7 +36,6 @@ int vcs_equilibrate(thermo_t& s, const char* XY,
{
MultiPhase* m = 0;
int retn = 1;
int retnSub = 0;
if (solver == 2) {
m = new MultiPhase;
@ -78,7 +77,7 @@ int vcs_equilibrate(thermo_t& s, const char* XY,
if (estimateEquil == 0) {
useThermoPhaseElementPotentials = true;
}
retnSub = e->equilibrate(s, XY,
int retnSub = e->equilibrate(s, XY,
useThermoPhaseElementPotentials, loglevel-1);
if (retnSub < 0) {
delete e;

77
src/equil/vcs_inest.cpp
View file

@ -22,22 +22,10 @@ static char pprefix[20] = " --- vcs_inest: ";
void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
double* const ss, double test)
{
size_t lt, ikl, kspec, iph, irxn;
double s;
double s1 = 0.0;
double xl, par;
bool finished, conv;
size_t nspecies = m_numSpeciesTot;
size_t nrxn = m_numRxnTot;
// double *molNum = VCS_DATA_PTR(m_molNumSpecies_old);
double TMolesMultiphase;
double* xtphMax = VCS_DATA_PTR(m_TmpPhase);
double* xtphMin = VCS_DATA_PTR(m_TmpPhase2);
ikl = 0;
lt = 0;
/*
* CALL ROUTINE TO SOLVE MAX(CC*molNum) SUCH THAT AX*molNum = BB
@ -52,7 +40,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
plogf("%s Mole Numbers returned from linear programming (vcs_inest initial guess):\n",
pprefix);
plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix);
for (kspec = 0; kspec < nspecies; ++kspec) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
plogf("%s ", pprefix);
plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]);
@ -62,31 +50,27 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
pprefix);
plogendl();
plogf("%s Element Goal Actual\n", pprefix);
int jj = 0;
for (size_t j = 0; j < m_numElemConstraints; j++) {
if (m_elementActive[j]) {
double tmp = 0.0;
for (kspec = 0; kspec < nspecies; ++kspec) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
tmp += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec];
}
plogf("%s ", pprefix);
plogf(" %-9.9s", (m_elementName[j]).c_str());
plogf(" %12.3g %12.3g\n", m_elemAbundancesGoal[j], tmp);
jj++;
}
}
plogendl();
}
#endif
/*
* Make sure all species have positive definite mole numbers
* Set voltages to zero for now, until we figure out what to do
*/
vcs_dzero(VCS_DATA_PTR(m_deltaMolNumSpecies), nspecies);
for (kspec = 0; kspec < nspecies; ++kspec) {
iph = m_phaseID[kspec];
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (m_molNumSpecies_old[kspec] <= 0.0) {
/*
@ -105,6 +89,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
* Now find the optimized basis that spans the stoichiometric
* coefficient matrix
*/
bool conv;
(void) vcs_basopt(false, aw, sa, sm, ss, test, &conv);
/* ***************************************************************** */
@ -118,22 +103,22 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
/*
* m_tPhaseMoles_new[] will consist of just the component moles
*/
for (iph = 0; iph < m_numPhases; iph++) {
for (size_t iph = 0; iph < m_numPhases; iph++) {
m_tPhaseMoles_new[iph] = TPhInertMoles[iph] + 1.0E-20;
}
for (kspec = 0; kspec < m_numComponents; ++kspec) {
for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
m_tPhaseMoles_new[m_phaseID[kspec]] += m_molNumSpecies_old[kspec];
}
}
TMolesMultiphase = 0.0;
for (iph = 0; iph < m_numPhases; iph++) {
double TMolesMultiphase = 0.0;
for (size_t iph = 0; iph < m_numPhases; iph++) {
if (! m_VolPhaseList[iph]->m_singleSpecies) {
TMolesMultiphase += m_tPhaseMoles_new[iph];
}
}
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_molNumSpecies_old), nspecies);
for (kspec = 0; kspec < m_numComponents; ++kspec) {
for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_MOLNUM) {
m_molNumSpecies_new[kspec] = 0.0;
}
@ -141,10 +126,10 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_SSfeSpecies),
nspecies);
for (kspec = 0; kspec < m_numComponents; ++kspec) {
for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
if (! m_SSPhase[kspec]) {
iph = m_phaseID[kspec];
size_t iph = m_phaseID[kspec];
m_feSpecies_new[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]);
}
} else {
@ -154,7 +139,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
vcs_deltag(0, true, VCS_STATECALC_NEW);
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
for (kspec = 0; kspec < nspecies; ++kspec) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
plogf("%s", pprefix);
plogf("%-12.12s", m_speciesName[kspec].c_str());
if (kspec < m_numComponents)
@ -169,13 +154,15 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
/* ********************************************************** */
/* **** ESTIMATE REACTION ADJUSTMENTS *********************** */
/* ********************************************************** */
double* xtphMax = VCS_DATA_PTR(m_TmpPhase);
double* xtphMin = VCS_DATA_PTR(m_TmpPhase2);
vcs_dzero(VCS_DATA_PTR(m_deltaPhaseMoles), m_numPhases);
for (iph = 0; iph < m_numPhases; iph++) {
for (size_t iph = 0; iph < m_numPhases; iph++) {
xtphMax[iph] = log(m_tPhaseMoles_new[iph] * 1.0E32);
xtphMin[iph] = log(m_tPhaseMoles_new[iph] * 1.0E-32);
}
for (irxn = 0; irxn < nrxn; ++irxn) {
kspec = m_indexRxnToSpecies[irxn];
for (size_t irxn = 0; irxn < nrxn; ++irxn) {
size_t kspec = m_indexRxnToSpecies[irxn];
/*
* For single species phases, we will not estimate the
* mole numbers. If the phase exists, it stays. If it
@ -183,7 +170,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
* existence here.
*/
if (! m_SSPhase[kspec]) {
iph = m_phaseID[kspec];
size_t iph = m_phaseID[kspec];
if (m_deltaGRxn_new[irxn] > xtphMax[iph]) {
m_deltaGRxn_new[irxn] = 0.8 * xtphMax[iph];
}
@ -213,7 +200,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
}
#ifdef DEBUG_MODE
if (m_debug_print_lvl >= 2) {
for (kspec = 0; kspec < nspecies; ++kspec) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf("%sdirection (", pprefix);
plogf("%-12.12s", m_speciesName[kspec].c_str());
@ -233,8 +220,8 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
/* *********************************************************** */
/* **** KEEP COMPONENT SPECIES POSITIVE ********************** */
/* *********************************************************** */
par = 0.5;
for (kspec = 0; kspec < m_numComponents; ++kspec) {
double par = 0.5;
for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (par < -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec]) {
par = -m_deltaMolNumSpecies[kspec] / m_molNumSpecies_new[kspec];
@ -250,16 +237,19 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
/* ******************************************** */
/* **** CALCULATE NEW MOLE NUMBERS ************ */
/* ******************************************** */
finished = false;
bool finished = false;
size_t lt = 0;
size_t ikl = 0;
double s1 = 0.0;
do {
for (kspec = 0; kspec < m_numComponents; ++kspec) {
for (size_t kspec = 0; kspec < m_numComponents; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[kspec] = m_molNumSpecies_new[kspec] + par * m_deltaMolNumSpecies[kspec];
} else {
m_deltaMolNumSpecies[kspec] = 0.0;
}
}
for (kspec = m_numComponents; kspec < nspecies; ++kspec) {
for (size_t kspec = m_numComponents; kspec < nspecies; ++kspec) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (m_deltaMolNumSpecies[kspec] != 0.0) {
m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par;
@ -279,7 +269,8 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
/* ******************************************* */
vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
vcs_dfe(VCS_STATECALC_OLD, 0, 0, nspecies);
for (kspec = 0, s = 0.0; kspec < nspecies; ++kspec) {
double s = 0.0;
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec];
}
if (s == 0.0) {
@ -304,7 +295,7 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm,
/* **************************************************** */
/* **** FIT PARABOLA THROUGH HALF AND FULL STEPS ****** */
/* **************************************************** */
xl = (1.0 - s / (s1 - s)) * 0.5;
double xl = (1.0 - s / (s1 - s)) * 0.5;
if (xl < 0.0) {
/* *************************************************** */
/* *** POOR DIRECTION, REDUCE STEP SIZE TO 0.2 ******* */
@ -332,7 +323,7 @@ finished:
plogf("%s Final Mole Numbers produced by inest:\n",
pprefix);
plogf("%s SPECIES MOLE_NUMBER\n", pprefix);
for (kspec = 0; kspec < nspecies; ++kspec) {
for (size_t kspec = 0; kspec < nspecies; ++kspec) {
plogf("%s ", pprefix);
plogf("%-12.12s", m_speciesName[kspec].c_str());
plogf(" %g", m_molNumSpecies_old[kspec]);
@ -345,9 +336,7 @@ finished:
int VCS_SOLVE::vcs_inest_TP()
{
int retn = 0;
double test;
Cantera::clockWC tickTock;
test = -1.0E20;
if (m_doEstimateEquil > 0) {
/*
@ -398,6 +387,7 @@ int VCS_SOLVE::vcs_inest_TP()
plogendl();
}
#endif
double test = -1.0E20;
vcs_inest(VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), VCS_DATA_PTR(sm),
VCS_DATA_PTR(ss), test);
/*
@ -472,8 +462,7 @@ int VCS_SOLVE::vcs_inest_TP()
/*
* Record time
*/
double tsecond = tickTock.secondsWC();
m_VCount->T_Time_inest += tsecond;
m_VCount->T_Time_inest += tickTock.secondsWC();
(m_VCount->T_Calls_Inest)++;
return retn;
}

34
src/equil/vcs_nondim.cpp
View file

@ -18,7 +18,6 @@ namespace VCSnonideal
{
double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const
{
double Farad;
if (TKelvin <= 0.0) {
TKelvin = 293.15;
}
@ -26,59 +25,47 @@ double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const
case VCS_UNITS_MKS:
case VCS_UNITS_KJMOL:
case VCS_UNITS_KCALMOL:
Farad = Cantera::ElectronCharge * Cantera::Avogadro /
return Cantera::ElectronCharge * Cantera::Avogadro /
(TKelvin * Cantera::GasConstant);
break;
case VCS_UNITS_UNITLESS:
Farad = Cantera::ElectronCharge * Cantera::Avogadro;
break;
return Cantera::ElectronCharge * Cantera::Avogadro;
case VCS_UNITS_KELVIN:
Farad = Cantera::ElectronCharge * Cantera::Avogadro/ TKelvin;
break;
return Cantera::ElectronCharge * Cantera::Avogadro/ TKelvin;
default:
plogf("vcs_nondim_Farad error: unknown units: %d\n", mu_units);
plogendl();
exit(EXIT_FAILURE);
}
return Farad;
}
double VCS_SOLVE::vcs_nondimMult_TP(int mu_units, double TKelvin) const
{
double rt;
if (TKelvin <= 0.0) {
TKelvin = 293.15;
}
switch (mu_units) {
case VCS_UNITS_KCALMOL:
rt = TKelvin * Cantera::GasConst_cal_mol_K * 1e-3;
break;
return TKelvin * Cantera::GasConst_cal_mol_K * 1e-3;
case VCS_UNITS_UNITLESS:
rt = 1.0;
break;
return 1.0;
case VCS_UNITS_KJMOL:
rt = TKelvin * Cantera::GasConstant * 1e-6;
break;
return TKelvin * Cantera::GasConstant * 1e-6;
case VCS_UNITS_KELVIN:
rt = TKelvin;
break;
return TKelvin;
case VCS_UNITS_MKS:
rt = TKelvin * Cantera::GasConstant;
break;
return TKelvin * Cantera::GasConstant;
default:
plogf("vcs_nondimMult_TP error: unknown units: %d\n", mu_units);
plogendl();
exit(EXIT_FAILURE);
}
return rt;
}
void VCS_SOLVE::vcs_nondim_TP()
{
double tf;
if (m_unitsState == VCS_DIMENSIONAL_G) {
m_unitsState = VCS_NONDIMENSIONAL_G;
tf = 1.0 / vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
double tf = 1.0 / vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
/*
* Modify the standard state and total chemical potential data,
@ -166,10 +153,9 @@ void VCS_SOLVE::vcs_nondim_TP()
void VCS_SOLVE::vcs_redim_TP(void)
{
double tf;
if (m_unitsState != VCS_DIMENSIONAL_G) {
m_unitsState = VCS_DIMENSIONAL_G;
tf = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
double tf = vcs_nondimMult_TP(m_VCS_UnitsFormat, m_temperature);
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
/*
* Modify the standard state and total chemical potential data,

125
src/equil/vcs_phaseStability.cpp
View file

@ -139,10 +139,6 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const
int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
{
int nfound = 0;
vcs_VolPhase* Vphase = 0;
double stoicC;
double molComp;
std::vector<int> linkedPhases;
phasePopProblemLists_.clear();
/*
@ -161,19 +157,17 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
*/
for (size_t j = 0; j < m_numComponents; j++) {
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
molComp = m_molNumSpecies_old[j];
if (molComp <= 0.0) {
if (m_molNumSpecies_old[j] <= 0.0) {
std::vector<size_t> &jList = zeroedComponentLinkedPhasePops[j];
size_t iph = m_phaseID[j];
jList.push_back(iph);
for (size_t irxn = 0; irxn < m_numRxnTot; irxn++) {
size_t kspec = irxn + m_numComponents;
iph = m_phaseID[kspec];
Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
int existence = Vphase->exists();
if (existence < 0) {
stoicC = m_stoichCoeffRxnMatrix[irxn][j];
if (stoicC > 0.0) {
if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) {
if (std::find(jList.begin(), jList.end(), iph) != jList.end()) {
jList.push_back(iph);
}
@ -190,15 +184,15 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
* Cut out components which have a pos stoichiometric value with another species in the phase.
*/
std::vector< std::vector<size_t> > zeroedPhaseLinkedZeroComponents(m_numPhases);
std::vector<int> linkedPhases;
/*
* The logic below calculates zeroedPhaseLinkedZeroComponents
*/
for (size_t iph = 0; iph < m_numPhases; iph++) {
std::vector<size_t> &iphList = zeroedPhaseLinkedZeroComponents[iph];
iphList.clear();
Vphase = m_VolPhaseList[iph];
int existence = Vphase->exists();
if (existence < 0) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
if (Vphase->exists() < 0) {
linkedPhases.clear();
size_t nsp = Vphase->nSpecies();
@ -208,10 +202,8 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
for (size_t j = 0; j < m_numComponents; j++) {
if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) {
molComp = m_molNumSpecies_old[j];
if (molComp <= 0.0) {
stoicC = m_stoichCoeffRxnMatrix[irxn][j];
if (stoicC < 0.0) {
if (m_molNumSpecies_old[j] <= 0.0) {
if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) {
bool foundPos = false;
for (size_t kk = 0; kk < nsp; kk++) {
size_t kkspec = Vphase->spGlobalIndexVCS(kk);
@ -240,9 +232,8 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
*
*/
for (size_t iph = 0; iph < m_numPhases; iph++) {
Vphase = m_VolPhaseList[iph];
int existence = Vphase->exists();
if (existence < 0) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
if (Vphase->exists() < 0) {
std::vector<size_t> &iphList = zeroedPhaseLinkedZeroComponents[iph];
std::vector<size_t> popProblem(0);
popProblem.push_back(iph);
@ -266,11 +257,8 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList()
size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
{
size_t iphasePop = npos;
size_t irxn, kspec;
doublereal FephaseMax = -1.0E30;
doublereal Fephase = -1.0E30;
vcs_VolPhase* Vphase = 0;
#ifdef DEBUG_MODE
char anote[128];
@ -281,7 +269,7 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
}
#endif
for (size_t iph = 0; iph < m_numPhases; iph++) {
Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
int existence = Vphase->exists();
#ifdef DEBUG_MODE
strcpy(anote, "");
@ -303,8 +291,8 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector<size_t> & phasePopPhaseIDs)
* Single Phase Stability Resolution
*
***********************************************************************/
kspec = Vphase->spGlobalIndexVCS(0);
irxn = kspec - m_numComponents;
size_t kspec = Vphase->spGlobalIndexVCS(0);
size_t irxn = kspec - m_numComponents;
doublereal deltaGRxn = m_deltaGRxn_old[irxn];
Fephase = exp(-deltaGRxn) - 1.0;
if (Fephase > 0.0) {
@ -403,16 +391,13 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
size_t irxn = kspec - m_numComponents;
std::vector<size_t> creationGlobalRxnNumbers;
doublereal s;
// Calculate the initial moles of the phase being born.
// Here we set it to 10x of the value which would cause the phase to be
// zeroed out within the algorithm. We may later adjust the value.
doublereal tPhaseMoles = 10. * m_totalMolNum * VCS_DELETE_PHASE_CUTOFF;
#ifdef DEBUG_MODE
int existence = Vphase->exists();
if (existence > 0) {
if (Vphase->exists() > 0) {
printf("ERROR vcs_popPhaseRxnStepSizes called for a phase that exists!");
exit(-1);
}
@ -424,7 +409,7 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop)
#endif
// Section for a single-species phase
if (Vphase->m_singleSpecies) {
s = 0.0;
double s = 0.0;
double* dnPhase_irxn = m_deltaMolNumPhase[irxn];
for (size_t j = 0; j < m_numComponents; ++j) {
if (!m_SSPhase[j]) {
@ -595,10 +580,8 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
/*
* We will use the _new state calc here
*/
size_t kspec, irxn, k, i, kc, kc_spec;
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
const size_t nsp = Vphase->nSpecies();
doublereal deltaGRxn;
int minNumberIterations = 3;
if (nsp <= 1) {
minNumberIterations = 1;
@ -635,14 +618,14 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
std::vector<size_t> componentList;
for (k = 0; k < nsp; k++) {
kspec = Vphase->spGlobalIndexVCS(k);
for (size_t k = 0; k < nsp; k++) {
size_t kspec = Vphase->spGlobalIndexVCS(k);
if (kspec < m_numComponents) {
componentList.push_back(k);
}
}
for (k = 0; k < m_numComponents; k++) {
for (size_t k = 0; k < m_numComponents; k++) {
m_feSpecies_Deficient[k] = m_feSpecies_old[k];
}
normUpdate = 0.1 * vcs_l2norm(fracDelta_new);
@ -663,7 +646,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
}
#endif
for (k = 0; k < nsp; k++) {
for (size_t k = 0; k < nsp; k++) {
if (fracDelta_new[k] < 1.0E-13) {
fracDelta_new[k] = 1.0E-13;
}
@ -678,14 +661,14 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
// Given a set of fracDelta's, we calculate the fracDelta's
// for the component species, if any
for (i = 0; i < componentList.size(); i++) {
kc = componentList[i];
kc_spec = Vphase->spGlobalIndexVCS(kc);
for (size_t i = 0; i < componentList.size(); i++) {
size_t kc = componentList[i];
size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
fracDelta_old[kc] = 0.0;
for (k = 0; k < nsp; k++) {
kspec = Vphase->spGlobalIndexVCS(k);
for (size_t k = 0; k < nsp; k++) {
size_t kspec = Vphase->spGlobalIndexVCS(k);
if (kspec >= m_numComponents) {
irxn = kspec - m_numComponents;
size_t irxn = kspec - m_numComponents;
fracDelta_old[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_old[k];
}
}
@ -693,7 +676,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
// Now, calculate the predicted mole fractions, X_est[k]
double sumFrac = 0.0;
for (k = 0; k < nsp; k++) {
for (size_t k = 0; k < nsp; k++) {
sumFrac += fracDelta_old[k];
}
// Necessary because this can be identically zero. -> we need to fix this algorithm!
@ -701,10 +684,9 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
sumFrac = 1.0;
}
double sum_Xcomp = 0.0;
for (k = 0; k < nsp; k++) {
for (size_t k = 0; k < nsp; k++) {
X_est[k] = fracDelta_old[k] / sumFrac;
kc_spec = Vphase->spGlobalIndexVCS(k);
if (kc_spec < m_numComponents) {
if (Vphase->spGlobalIndexVCS(k) < m_numComponents) {
sum_Xcomp += X_est[k];
}
}
@ -725,9 +707,9 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
* First calculate altered chemical potentials for component species
* belonging to this phase.
*/
for (i = 0; i < componentList.size(); i++) {
kc = componentList[i];
kc_spec = Vphase->spGlobalIndexVCS(kc);
for (size_t i = 0; i < componentList.size(); i++) {
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]
+ log(m_actCoeffSpecies_new[kc_spec] * X_est[kc]);
@ -737,14 +719,12 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
}
}
for (i = 0; i < componentList.size(); i++) {
kc = componentList[i];
kc_spec = Vphase->spGlobalIndexVCS(kc);
for (k = 0; k < Vphase->nSpecies(); k++) {
kspec = Vphase->spGlobalIndexVCS(k);
for (size_t i = 0; i < componentList.size(); i++) {
size_t kc_spec = Vphase->spGlobalIndexVCS(componentList[i]);
for (size_t k = 0; k < Vphase->nSpecies(); k++) {
size_t kspec = Vphase->spGlobalIndexVCS(k);
if (kspec >= m_numComponents) {
irxn = kspec - m_numComponents;
size_t irxn = kspec - m_numComponents;
if (i == 0) {
m_deltaGRxn_Deficient[irxn] = m_deltaGRxn_old[irxn];
}
@ -763,11 +743,11 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
*/
sum = 0.0;
funcPhaseStability = sum_Xcomp - 1.0;
for (k = 0; k < nsp; k++) {
kspec = Vphase->spGlobalIndexVCS(k);
for (size_t k = 0; k < nsp; k++) {
size_t kspec = Vphase->spGlobalIndexVCS(k);
if (kspec >= m_numComponents) {
irxn = kspec - m_numComponents;
deltaGRxn = m_deltaGRxn_Deficient[irxn];
size_t irxn = kspec - m_numComponents;
double deltaGRxn = m_deltaGRxn_Deficient[irxn];
if (deltaGRxn > 50.0) {
deltaGRxn = 50.0;
}
@ -785,11 +765,10 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
/*
* Calculate the raw estimate of the new fracs
*/
for (k = 0; k < nsp; k++) {
kspec = Vphase->spGlobalIndexVCS(k);
for (size_t k = 0; k < nsp; k++) {
size_t kspec = Vphase->spGlobalIndexVCS(k);
double b = E_phi[k] / sum * (1.0 - sum_Xcomp);
if (kspec >= m_numComponents) {
irxn = kspec - m_numComponents;
fracDelta_raw[k] = b;
}
}
@ -797,14 +776,14 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
// Given a set of fracDelta's, we calculate the fracDelta's
// for the component species, if any
for (i = 0; i < componentList.size(); i++) {
kc = componentList[i];
kc_spec = Vphase->spGlobalIndexVCS(kc);
for (size_t i = 0; i < componentList.size(); i++) {
size_t kc = componentList[i];
size_t kc_spec = Vphase->spGlobalIndexVCS(kc);
fracDelta_raw[kc] = 0.0;
for (k = 0; k < nsp; k++) {
kspec = Vphase->spGlobalIndexVCS(k);
for (size_t k = 0; k < nsp; k++) {
size_t kspec = Vphase->spGlobalIndexVCS(k);
if (kspec >= m_numComponents) {
irxn = kspec - m_numComponents;
size_t irxn = kspec - m_numComponents;
fracDelta_raw[kc] += m_stoichCoeffRxnMatrix[irxn][kc_spec] * fracDelta_raw[k];
}
}
@ -816,14 +795,14 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
* Now possibly dampen the estimate.
*/
doublereal sumADel = 0.0;
for (k = 0; k < nsp; k++) {
for (size_t k = 0; k < nsp; k++) {
delFrac[k] = fracDelta_raw[k] - fracDelta_old[k];
sumADel += fabs(delFrac[k]);
}
normUpdate = vcs_l2norm(delFrac);
dirProd = 0.0;
for (k = 0; k < nsp; k++) {
for (size_t k = 0; k < nsp; k++) {
dirProd += fracDelta_old[k] * delFrac[k];
}
bool crossedSign = false;
@ -850,7 +829,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
}
}
for (k = 0; k < nsp; k++) {
for (size_t k = 0; k < nsp; k++) {
if (fabs(damp * delFrac[k]) > 0.3*fabs(fracDelta_old[k])) {
damp = std::max(0.3*fabs(fracDelta_old[k]) / fabs(delFrac[k]),
1.0E-8/fabs(delFrac[k]));
@ -870,7 +849,7 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph)
damp = 0.000001;
}
for (k = 0; k < nsp; k++) {
for (size_t k = 0; k < nsp; k++) {
fracDelta_new[k] = fracDelta_old[k] + damp * (delFrac[k]);
}

47
src/equil/vcs_prep.cpp
View file

@ -19,9 +19,6 @@ namespace VCSnonideal
{
void VCS_SOLVE::vcs_SSPhase()
{
size_t iph;
vcs_VolPhase* Vphase;
std::vector<int> numPhSpecies(m_numPhases, 0);
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
@ -32,8 +29,8 @@ void VCS_SOLVE::vcs_SSPhase()
* has been earmarked as a multispecies phase.
* Treat that species as a single-species phase
*/
for (iph = 0; iph < m_numPhases; iph++) {
Vphase = m_VolPhaseList[iph];
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
Vphase->m_singleSpecies = false;
if (TPhInertMoles[iph] > 0.0) {
Vphase->setExistence(2);
@ -52,8 +49,8 @@ void VCS_SOLVE::vcs_SSPhase()
* single species phase or not.
*/
for (size_t kspec = 0; kspec < m_numSpeciesTot; kspec++) {
iph = m_phaseID[kspec];
Vphase = m_VolPhaseList[iph];
size_t iph = m_phaseID[kspec];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
if (Vphase->m_singleSpecies) {
m_SSPhase[kspec] = true;
} else {
@ -64,13 +61,7 @@ void VCS_SOLVE::vcs_SSPhase()
int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
{
size_t kspec, i;
int retn = VCS_SUCCESS;
double pres, test;
double* aw, *sa, *sm, *ss;
bool modifiedSoln = false;
bool conv;
m_debug_print_lvl = printLvl;
/*
@ -90,11 +81,11 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
}
m_numRxnRdc = m_numRxnTot;
m_numSpeciesRdc = m_numSpeciesTot;
for (i = 0; i < m_numRxnRdc; ++i) {
for (size_t i = 0; i < m_numRxnRdc; ++i) {
m_indexRxnToSpecies[i] = m_numElemConstraints + i;
}
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
size_t pID = m_phaseID[kspec];
size_t spPhIndex = m_speciesLocalPhaseIndex[kspec];
vcs_VolPhase* vPhase = m_VolPhaseList[pID];
@ -130,23 +121,20 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
*
* For voltage unknowns, set these to zero for the moment.
*/
test = -1.0e-10;
double test = -1.0e-10;
bool modifiedSoln = false;
if (m_doEstimateEquil < 0) {
double sum = 0.0;
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
sum += fabs(m_molNumSpecies_old[kspec]);
}
}
if (fabs(sum) < 1.0E-6) {
modifiedSoln = true;
if (m_pressurePA <= 0.0) {
pres = 1.01325E5;
} else {
pres = m_pressurePA;
}
double pres = (m_pressurePA <= 0.0) ? 1.01325E5 : m_pressurePA;
retn = vcs_evalSS_TP(0, 0, m_temperature, pres);
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec];
} else {
@ -163,14 +151,15 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
* reaction matrix.
*/
std::vector<double> awSpace(m_numSpeciesTot + (m_numElemConstraints + 2)*(m_numElemConstraints), 0.0);
aw = VCS_DATA_PTR(awSpace);
double* aw = VCS_DATA_PTR(awSpace);
if (aw == NULL) {
plogf("vcs_prep_oneTime: failed to get memory: global bailout\n");
return VCS_NOMEMORY;
}
sa = aw + m_numSpeciesTot;
sm = sa + m_numElemConstraints;
ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
double* sa = aw + m_numSpeciesTot;
double* sm = sa + m_numElemConstraints;
double* ss = sm + (m_numElemConstraints)*(m_numElemConstraints);
bool conv;
retn = vcs_basopt(true, aw, sa, sm, ss, test, &conv);
if (retn != VCS_SUCCESS) {
plogf("vcs_prep_oneTime:");
@ -182,7 +171,7 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
if (m_numSpeciesTot >= m_numComponents) {
m_numRxnTot = m_numRxnRdc = m_numSpeciesTot - m_numComponents;
for (i = 0; i < m_numRxnRdc; ++i) {
for (size_t i = 0; i < m_numRxnRdc; ++i) {
m_indexRxnToSpecies[i] = m_numComponents + i;
}
} else {
@ -209,7 +198,7 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
// If we mucked up the solution unknowns because they were all
// zero to start with, set them back to zero here
if (modifiedSoln) {
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
m_molNumSpecies_old[kspec] = 0.0;
}
}

60
src/equil/vcs_prob.cpp
View file

@ -184,7 +184,6 @@ static void print_char(const char letter, const int num)
void VCS_PROB::prob_report(int print_lvl)
{
m_printLvl = print_lvl;
vcs_VolPhase* Vphase = 0;
/*
* Printout the species information: PhaseID's and mole nums
*/
@ -215,7 +214,7 @@ void VCS_PROB::prob_report(int print_lvl)
plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_Moles Species_Type\n");
for (size_t i = 0; i < nspecies; i++) {
Vphase = VPhaseList[PhaseID[i]];
vcs_VolPhase* Vphase = VPhaseList[PhaseID[i]];
plogf("%16s %5d %16s", SpName[i].c_str(), PhaseID[i],
Vphase->PhaseName.c_str());
if (iest >= 0) {
@ -245,7 +244,7 @@ void VCS_PROB::prob_report(int print_lvl)
plogf(" TMolesInert TKmoles\n");
for (size_t iphase = 0; iphase < NPhase; iphase++) {
Vphase = VPhaseList[iphase];
vcs_VolPhase* Vphase = VPhaseList[iphase];
std::string EOS_cstr = string16_EOSType(Vphase->m_eqnState);
plogf("%16s %5d %5d %8d ", Vphase->PhaseName.c_str(),
Vphase->VP_ID_, Vphase->m_singleSpecies, Vphase->m_gasPhase);
@ -288,7 +287,7 @@ void VCS_PROB::prob_report(int print_lvl)
plogf(" Species (phase) "
" SS0ChemPot StarChemPot\n");
for (size_t iphase = 0; iphase < NPhase; iphase++) {
Vphase = VPhaseList[iphase];
vcs_VolPhase* Vphase = VPhaseList[iphase];
Vphase->setState_TP(T, PresPA);
for (size_t kindex = 0; kindex < Vphase->nSpecies(); kindex++) {
size_t kglob = Vphase->spGlobalIndexVCS(kindex);
@ -317,24 +316,20 @@ void VCS_PROB::prob_report(int print_lvl)
void VCS_PROB::addPhaseElements(vcs_VolPhase* volPhase)
{
size_t e, eVP;
size_t foundPos = npos;
size_t neVP = volPhase->nElemConstraints();
std::string en;
std::string enVP;
/*
* Loop through the elements in the vol phase object
*/
for (eVP = 0; eVP < neVP; eVP++) {
foundPos = npos;
enVP = volPhase->elementName(eVP);
for (size_t eVP = 0; eVP < neVP; eVP++) {
size_t foundPos = npos;
std::string enVP = volPhase->elementName(eVP);
/*
* Search for matches with the existing elements.
* If found, then fill in the entry in the global
* mapping array.
*/
for (e = 0; e < ne; e++) {
en = ElName[e];
for (size_t e = 0; e < ne; e++) {
std::string en = ElName[e];
if (!strcmp(enVP.c_str(), en.c_str())) {
volPhase->setElemGlobalIndex(eVP, e);
foundPos = e;
@ -343,7 +338,7 @@ void VCS_PROB::addPhaseElements(vcs_VolPhase* volPhase)
if (foundPos == npos) {
int elType = volPhase->elementType(eVP);
int elactive = volPhase->elementActive(eVP);
e = addElement(enVP.c_str(), elType, elactive);
size_t e = addElement(enVP.c_str(), elType, elactive);
volPhase->setElemGlobalIndex(eVP, e);
}
}
@ -366,7 +361,6 @@ size_t VCS_PROB::addElement(const char* elNameNew, int elType, int elactive)
size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT)
{
size_t e, eVP;
if (kT > nspecies) {
/*
* Need to expand the number of species here
@ -375,8 +369,8 @@ size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT)
exit(EXIT_FAILURE);
}
double const* const* const fm = volPhase->getFormulaMatrix();
for (eVP = 0; eVP < volPhase->nElemConstraints(); eVP++) {
e = volPhase->elemGlobalIndex(eVP);
for (size_t eVP = 0; eVP < volPhase->nElemConstraints(); eVP++) {
size_t e = volPhase->elemGlobalIndex(eVP);
#ifdef DEBUG_MODE
if (e == npos) {
exit(EXIT_FAILURE);
@ -394,18 +388,11 @@ size_t VCS_PROB::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT)
void VCS_PROB::reportCSV(const std::string& reportFile)
{
size_t k;
size_t istart;
double vol = 0.0;
string sName;
FILE* FP = fopen(reportFile.c_str(), "w");
if (!FP) {
plogf("Failure to open file\n");
exit(EXIT_FAILURE);
}
double Temp = T;
std::vector<double> volPM(nspecies, 0.0);
std::vector<double> activity(nspecies, 0.0);
@ -414,10 +401,10 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
std::vector<double> mu0(nspecies, 0.0);
std::vector<double> molalities(nspecies, 0.0);
vol = 0.0;
double vol = 0.0;
size_t iK = 0;
for (size_t iphase = 0; iphase < NPhase; iphase++) {
istart = iK;
size_t istart = iK;
vcs_VolPhase* volP = VPhaseList[iphase];
//const Cantera::ThermoPhase *tptr = volP->ptrThermoPhase();
size_t nSpeciesPhase = volP->nSpecies();
@ -426,7 +413,7 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
double TMolesPhase = volP->totalMoles();
double VolPhaseVolumes = 0.0;
for (k = 0; k < nSpeciesPhase; k++) {
for (size_t k = 0; k < nSpeciesPhase; k++) {
iK++;
VolPhaseVolumes += volPM[istart + k] * mf[istart + k];
}
@ -436,7 +423,7 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
fprintf(FP,"--------------------- VCS_MULTIPHASE_EQUIL FINAL REPORT"
" -----------------------------\n");
fprintf(FP,"Temperature = %11.5g kelvin\n", Temp);
fprintf(FP,"Temperature = %11.5g kelvin\n", T);
fprintf(FP,"Pressure = %11.5g Pascal\n", PresPA);
fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
fprintf(FP,"Number Basis optimizations = %d\n", m_NumBasisOptimizations);
@ -444,7 +431,7 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
iK = 0;
for (size_t iphase = 0; iphase < NPhase; iphase++) {
istart = iK;
size_t istart = iK;
vcs_VolPhase* volP = VPhaseList[iphase];
const Cantera::ThermoPhase* tp = volP->ptrThermoPhase();
@ -469,13 +456,12 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
tp->getPartialMolarVolumes(VCS_DATA_PTR(volPM));
tp->getChemPotentials(VCS_DATA_PTR(mu));
double VolPhaseVolumes = 0.0;
for (k = 0; k < nSpeciesPhase; k++) {
for (size_t k = 0; k < nSpeciesPhase; k++) {
VolPhaseVolumes += volPM[k] * mf[istart + k];
}
VolPhaseVolumes *= TMolesPhase;
vol += VolPhaseVolumes;
if (actConvention == 1) {
const Cantera::MolalityVPSSTP* mTP = static_cast<const Cantera::MolalityVPSSTP*>(tp);
tp->getChemPotentials(VCS_DATA_PTR(mu));
@ -491,8 +477,8 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
" , , ,"
" (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
}
for (k = 0; k < nSpeciesPhase; k++) {
sName = tp->speciesName(k);
for (size_t k = 0; k < nSpeciesPhase; k++) {
std::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(),
@ -513,11 +499,11 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
" , , ,"
" (J/kmol), (J/kmol), (kmol), (m**3/kmol), (m**3)\n");
}
for (k = 0; k < nSpeciesPhase; k++) {
for (size_t k = 0; k < nSpeciesPhase; k++) {
molalities[k] = 0.0;
}
for (k = 0; k < nSpeciesPhase; k++) {
sName = tp->speciesName(k);
for (size_t k = 0; k < nSpeciesPhase; k++) {
std::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(),
@ -534,7 +520,7 @@ void VCS_PROB::reportCSV(const std::string& reportFile)
* Check consistency: These should be equal
*/
tp->getChemPotentials(VCS_DATA_PTR(m_gibbsSpecies)+istart);
for (k = 0; k < nSpeciesPhase; k++) {
for (size_t k = 0; k < nSpeciesPhase; k++) {
if (!vcs_doubleEqual(m_gibbsSpecies[istart+k], mu[k])) {
fprintf(FP,"ERROR: incompatibility!\n");
fclose(FP);

7
src/equil/vcs_rearrange.cpp
View file

@ -15,17 +15,16 @@ namespace VCSnonideal
{
int VCS_SOLVE::vcs_rearrange()
{
size_t i, l, j;
size_t k1 = 0;
/* - Loop over all of the species */
for (i = 0; i < m_numSpeciesTot; ++i) {
for (size_t i = 0; i < m_numSpeciesTot; ++i) {
/*
* Find the index of I in the index vector m_speciesIndexVector[].
* Call it k1 and continue.
*/
for (j = 0; j < m_numSpeciesTot; ++j) {
l = m_speciesMapIndex[j];
for (size_t j = 0; j < m_numSpeciesTot; ++j) {
size_t l = m_speciesMapIndex[j];
k1 = j;
if (l == i) {
break;

53
src/equil/vcs_report.cpp
View file

@ -29,13 +29,10 @@ static void print_line(const std::string& schar, size_t num)
int VCS_SOLVE::vcs_report(int iconv)
{
bool printActualMoles = true, inertYes = false;
size_t i, j, l, k, kspec;
size_t nspecies = m_numSpeciesTot;
double g;
char originalUnitsState = m_unitsState;
std::vector<size_t> sortindex(nspecies,0);
std::vector<double> xy(nspecies,0.0);
@ -43,7 +40,7 @@ int VCS_SOLVE::vcs_report(int iconv)
/* **** SORT DEPENDENT SPECIES IN DECREASING ORDER OF MOLES ***** */
/* ************************************************************** */
for (i = 0; i < nspecies; ++i) {
for (size_t i = 0; i < nspecies; ++i) {
sortindex[i] = i;
xy[i] = m_molNumSpecies_old[i];
}
@ -52,8 +49,8 @@ int VCS_SOLVE::vcs_report(int iconv)
* and the sort index vector, sortindex, according to
* the magnitude of the mole fraction vector.
*/
for (l = m_numComponents; l < m_numSpeciesRdc; ++l) {
k = vcs_optMax(VCS_DATA_PTR(xy), 0, l, m_numSpeciesRdc);
for (size_t l = m_numComponents; l < m_numSpeciesRdc; ++l) {
size_t k = vcs_optMax(VCS_DATA_PTR(xy), 0, l, m_numSpeciesRdc);
if (k != l) {
std::swap(xy[k], xy[l]);
std::swap(sortindex[k], sortindex[l]);
@ -113,7 +110,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" Species Equilibrium kmoles ");
plogf("Mole Fraction ChemPot/RT SpecUnkType\n");
print_line("-", 80);
for (i = 0; i < m_numComponents; ++i) {
for (size_t i = 0; i < m_numComponents; ++i) {
plogf(" %-12.12s", m_speciesName[i].c_str());
print_space(13);
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i] * molScale,
@ -121,8 +118,8 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %3d", m_speciesUnknownType[i]);
plogf("\n");
}
for (i = m_numComponents; i < m_numSpeciesRdc; ++i) {
l = sortindex[i];
for (size_t i = m_numComponents; i < m_numSpeciesRdc; ++i) {
size_t l = sortindex[i];
plogf(" %-12.12s", m_speciesName[l].c_str());
print_space(13);
@ -139,7 +136,7 @@ int VCS_SOLVE::vcs_report(int iconv)
}
plogf("\n");
}
for (i = 0; i < m_numPhases; i++) {
for (size_t i = 0; i < m_numPhases; i++) {
if (TPhInertMoles[i] > 0.0) {
inertYes = true;
if (i == 0) {
@ -154,15 +151,15 @@ int VCS_SOLVE::vcs_report(int iconv)
}
if (m_numSpeciesRdc != nspecies) {
plogf("\n SPECIES WITH LESS THAN 1.0E-32 KMOLES:\n\n");
for (kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) {
for (size_t kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) {
plogf(" %-12.12s", m_speciesName[kspec].c_str());
// Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
plogf(" %14.7E %14.7E %12.4E",
m_molNumSpecies_old[kspec]*molScale,
m_molNumSpecies_new[kspec]*molScale, m_deltaGRxn_new[kspec]);
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" KMol_Num");
} else if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
} else if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Voltage");
} else {
plogf(" Unknown");
@ -180,17 +177,17 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("\n");
print_line("-", m_numComponents*10 + 45);
plogf(" |ComponentID|");
for (j = 0; j < m_numComponents; j++) {
for (size_t j = 0; j < m_numComponents; j++) {
plogf(" %3d", j);
}
plogf(" | |\n");
plogf(" | Components|");
for (j = 0; j < m_numComponents; j++) {
for (size_t j = 0; j < m_numComponents; j++) {
plogf(" %10.10s", m_speciesName[j].c_str());
}
plogf(" | |\n");
plogf(" NonComponent | Moles |");
for (j = 0; j < m_numComponents; j++) {
for (size_t j = 0; j < m_numComponents; j++) {
plogf(" %10.3g", m_molNumSpecies_old[j] * molScale);
}
plogf(" | DG/RT Rxn |\n");
@ -200,7 +197,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %3d ", kspec);
plogf("%-10.10s", m_speciesName[kspec].c_str());
plogf("|%10.3g |", m_molNumSpecies_old[kspec]*molScale);
for (j = 0; j < m_numComponents; j++) {
for (size_t j = 0; j < m_numComponents; j++) {
plogf(" %6.2f", m_stoichCoeffRxnMatrix[irxn][j]);
}
plogf(" |%10.3g |", m_deltaGRxn_new[irxn]);
@ -221,17 +218,17 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("\n");
print_line("-", m_numElemConstraints*10 + 58);
plogf(" | ElementID |");
for (j = 0; j < m_numElemConstraints; j++) {
for (size_t j = 0; j < m_numElemConstraints; j++) {
plogf(" %3d", j);
}
plogf(" | |\n");
plogf(" | Element |");
for (j = 0; j < m_numElemConstraints; j++) {
for (size_t j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.10s", (m_elementName[j]).c_str());
}
plogf(" | |\n");
plogf(" PhaseName |KMolTarget |");
for (j = 0; j < m_numElemConstraints; j++) {
for (size_t j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", m_elemAbundancesGoal[j]);
}
plogf(" | Gibbs Total |\n");
@ -249,7 +246,7 @@ int VCS_SOLVE::vcs_report(int iconv)
}
}
vcs_elabPhase(iphase, VCS_DATA_PTR(gaPhase));
for (j = 0; j < m_numElemConstraints; j++) {
for (size_t j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", gaPhase[j]);
gaTPhase[j] += gaPhase[j];
}
@ -260,7 +257,7 @@ int VCS_SOLVE::vcs_report(int iconv)
}
print_line("-", m_numElemConstraints*10 + 58);
plogf(" TOTAL |%10.3e |", totalMoles);
for (j = 0; j < m_numElemConstraints; j++) {
for (size_t j = 0; j < m_numElemConstraints; j++) {
plogf(" %10.3g", gaTPhase[j]);
}
plogf(" | %18.11E |\n", gibbsTotal);
@ -278,8 +275,8 @@ int VCS_SOLVE::vcs_report(int iconv)
* energy of zero
*/
g = vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
VCS_DATA_PTR(m_tPhaseMoles_old));
double g = vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
VCS_DATA_PTR(m_tPhaseMoles_old));
plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
if (inertYes) {
plogf("\t\t(Inert species have standard free energy of zero)\n");
@ -287,7 +284,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("\nElemental Abundances (kmol): ");
plogf(" Actual Target Type ElActive\n");
for (i = 0; i < m_numElemConstraints; ++i) {
for (size_t i = 0; i < m_numElemConstraints; ++i) {
print_space(26);
plogf("%-2.2s", (m_elementName[i]).c_str());
plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale);
@ -311,8 +308,8 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("| (MolNum ChemPot)|");
plogf("\n");
print_line("-", 147);
for (i = 0; i < nspecies; ++i) {
l = sortindex[i];
for (size_t i = 0; i < nspecies; ++i) {
size_t l = sortindex[i];
size_t pid = m_phaseID[l];
plogf(" %-12.12s", m_speciesName[l].c_str());
plogf(" %14.7E ", m_molNumSpecies_old[l]*molScale);
@ -351,7 +348,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("| %20.9E |", m_feSpecies_old[l] * m_molNumSpecies_old[l] * molScale);
plogf("\n");
}
for (i = 0; i < 125; i++) {
for (size_t i = 0; i < 125; i++) {
plogf(" ");
}
plogf(" %20.9E\n", g);

80
src/equil/vcs_rxnadj.cpp
View file

@ -19,12 +19,8 @@ namespace VCSnonideal
size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
{
size_t kspec, iph;
size_t iphDel = npos;
double s, xx, dss;
size_t k = 0;
vcs_VolPhase* Vphase = 0;
double* dnPhase_irxn;
#ifdef DEBUG_MODE
char ANOTE[128];
if (m_debug_print_lvl >= 2) {
@ -59,7 +55,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
sprintf(ANOTE, "Normal Calc");
#endif
kspec = m_indexRxnToSpecies[irxn];
size_t kspec = m_indexRxnToSpecies[irxn];
if (m_speciesStatus[kspec] == VCS_SPECIES_ZEROEDPHASE) {
m_deltaMolNumSpecies[kspec] = 0.0;
@ -68,7 +64,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
#endif
} else if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
dnPhase_irxn = m_deltaMolNumPhase[irxn];
double* dnPhase_irxn = m_deltaMolNumPhase[irxn];
if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) {
/********************************************************************/
@ -84,10 +80,10 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
* First decide if this species is part of a multiphase that
* is nontrivial in size.
*/
iph = m_phaseID[kspec];
size_t iph = m_phaseID[kspec];
double tphmoles = m_tPhaseMoles_old[iph];
double trphmoles = tphmoles / m_totalMolNum;
Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
if (Vphase->exists() && (trphmoles > VCS_DELETE_PHASE_CUTOFF)) {
m_deltaMolNumSpecies[kspec] = m_totalMolNum * VCS_SMALL_MULTIPHASE_SPECIES;
if (m_speciesStatus[kspec] == VCS_SPECIES_STOICHZERO) {
@ -169,6 +165,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
/*
* Start of the regular processing
*/
double s;
if (m_SSPhase[kspec]) {
s = 0.0;
} else {
@ -182,7 +179,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
}
}
for (size_t j = 0; j < m_numPhases; j++) {
Vphase = m_VolPhaseList[j];
vcs_VolPhase* Vphase = m_VolPhaseList[j];
if (!Vphase->m_singleSpecies) {
if (m_tPhaseMoles_old[j] > 0.0) {
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
@ -254,12 +251,13 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
* will zero out first.
* -> The species to be zeroed out will be "k".
*/
double dss;
if (m_deltaGRxn_new[irxn] > 0.0) {
dss = m_molNumSpecies_old[kspec];
k = kspec;
for (size_t j = 0; j < m_numComponents; ++j) {
if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) {
xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
double xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
@ -271,7 +269,7 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
dss = 1.0e10;
for (size_t j = 0; j < m_numComponents; ++j) {
if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) {
xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
double xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
@ -382,16 +380,11 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial)
int VCS_SOLVE::vcs_rxn_adj_cg()
{
size_t irxn, j;
size_t k = 0;
size_t kspec;
int soldel = 0;
double s, xx, dss;
double* dnPhase_irxn;
#ifdef DEBUG_MODE
char ANOTE[128];
plogf(" ");
for (j = 0; j < 77; j++) {
for (size_t j = 0; j < 77; j++) {
plogf("-");
}
plogf("\n --- Subroutine rxn_adj_cg() called\n");
@ -404,13 +397,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
* We also evaluate whether the matrix is appropriate for
* this algorithm. If not, we bail out.
*/
for (irxn = 0; irxn < m_numRxnRdc; ++irxn) {
for (size_t irxn = 0; irxn < m_numRxnRdc; ++irxn) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Normal Calc");
#endif
kspec = m_indexRxnToSpecies[irxn];
dnPhase_irxn = m_deltaMolNumPhase[irxn];
size_t kspec = m_indexRxnToSpecies[irxn];
double* dnPhase_irxn = m_deltaMolNumPhase[irxn];
if (m_molNumSpecies_old[kspec] == 0.0 && (!m_SSPhase[kspec])) {
/* *******************************************************************/
@ -471,17 +464,18 @@ 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];
}
for (j = 0; j < m_numComponents; ++j) {
for (size_t j = 0; j < m_numComponents; ++j) {
if (!m_SSPhase[j]) {
s += SQUARE(m_stoichCoeffRxnMatrix[irxn][j]) / m_molNumSpecies_old[j];
}
}
for (j = 0; j < m_numPhases; j++) {
for (size_t j = 0; j < m_numPhases; j++) {
if (!(m_VolPhaseList[j])->m_singleSpecies) {
if (m_tPhaseMoles_old[j] > 0.0) {
s -= SQUARE(dnPhase_irxn[j]) / m_tPhaseMoles_old[j];
@ -502,12 +496,14 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
* Then, we need to follow the reaction to see which species
* will zero out first.
*/
size_t k;
double dss;
if (m_deltaGRxn_new[irxn] > 0.0) {
dss = m_molNumSpecies_old[kspec];
k = kspec;
for (j = 0; j < m_numComponents; ++j) {
for (size_t j = 0; j < m_numComponents; ++j) {
if (m_stoichCoeffRxnMatrix[irxn][j] > 0.0) {
xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
double xx = m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
@ -517,9 +513,9 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
dss = -dss;
} else {
dss = 1.0e10;
for (j = 0; j < m_numComponents; ++j) {
for (size_t j = 0; j < m_numComponents; ++j) {
if (m_stoichCoeffRxnMatrix[irxn][j] < 0.0) {
xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
double xx = -m_molNumSpecies_old[j] / m_stoichCoeffRxnMatrix[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
@ -538,7 +534,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
if (dss != 0.0) {
m_molNumSpecies_old[kspec] += dss;
m_tPhaseMoles_old[m_phaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
for (size_t j = 0; j < m_numComponents; ++j) {
m_molNumSpecies_old[j] += dss * m_stoichCoeffRxnMatrix[irxn][j];
m_tPhaseMoles_old[m_phaseID[j]] += dss * m_stoichCoeffRxnMatrix[irxn][j];
}
@ -582,7 +578,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg()
#ifdef DEBUG_MODE
plogf(" ");
for (j = 0; j < 77; j++) {
for (size_t j = 0; j < 77; j++) {
plogf("-");
}
plogf("\n");
@ -610,27 +606,24 @@ double VCS_SOLVE::vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal)
double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn)
{
size_t kspec, k, l, kph;
double s;
double* sc_irxn;
kspec = m_indexRxnToSpecies[irxn];
kph = m_phaseID[kspec];
size_t kspec = m_indexRxnToSpecies[irxn];
size_t kph = m_phaseID[kspec];
double np_kspec = m_tPhaseMoles_old[kph];
if (np_kspec < 1.0E-13) {
np_kspec = 1.0E-13;
}
sc_irxn = m_stoichCoeffRxnMatrix[irxn];
double* sc_irxn = m_stoichCoeffRxnMatrix[irxn];
/*
* First the diagonal term of the Jacobian
*/
s = m_np_dLnActCoeffdMolNum[kspec][kspec] / np_kspec;
double s = m_np_dLnActCoeffdMolNum[kspec][kspec] / np_kspec;
/*
* Next, the other terms. Note this only a loop over the components
* So, it's not too expensive to calculate.
*/
for (l = 0; l < m_numComponents; l++) {
for (size_t l = 0; l < m_numComponents; l++) {
if (!m_SSPhase[l]) {
for (k = 0; k < m_numComponents; ++k) {
for (size_t k = 0; k < m_numComponents; ++k) {
if (m_phaseID[k] == m_phaseID[l]) {
double np = m_tPhaseMoles_old[m_phaseID[k]];
if (np > 0.0) {
@ -697,7 +690,6 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
#endif
{
int its = 0;
size_t k;
size_t kspec = m_indexRxnToSpecies[irxn];
const int MAXITS = 10;
double dx = dx_orig;
@ -705,8 +697,6 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
double* molNumBase = VCS_DATA_PTR(m_molNumSpecies_old);
double* acBase = VCS_DATA_PTR(m_actCoeffSpecies_old);
double* ac = VCS_DATA_PTR(m_actCoeffSpecies_new);
double molSum = 0.0;
double slope;
/*
* Calculate the deltaG value at the dx = 0.0 point
*/
@ -743,9 +733,9 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
}
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNumBase, m_numSpeciesRdc);
molSum = molNumBase[kspec];
double molSum = molNumBase[kspec];
m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx_orig;
for (k = 0; k < m_numComponents; k++) {
for (size_t k = 0; k < m_numComponents; k++) {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx_orig;
molSum += molNumBase[k];
}
@ -769,7 +759,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
*/
if (fabs(deltaG1) < 0.8 * forig) {
if (deltaG1 * deltaGOrig < 0.0) {
slope = (deltaG1 - deltaGOrig) / dx_orig;
double slope = (deltaG1 - deltaGOrig) / dx_orig;
dx = -deltaGOrig / slope;
} else {
dx = dx_orig;
@ -786,7 +776,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
*/
dx *= 0.5;
m_molNumSpecies_new[kspec] = molNumBase[kspec] + dx;
for (k = 0; k < m_numComponents; k++) {
for (size_t k = 0; k < m_numComponents; k++) {
m_molNumSpecies_new[k] = molNumBase[k] + sc_irxn[k] * dx;
}
vcs_setFlagsVolPhases(false, VCS_STATECALC_NEW);
@ -806,7 +796,7 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig)
*/
if (fabs(deltaG) / forig < (1.0 - 0.1 * dx / dx_orig)) {
if (deltaG * deltaGOrig < 0.0) {
slope = (deltaG - deltaGOrig) / dx;
double slope = (deltaG - deltaGOrig) / dx;
dx = -deltaGOrig / slope;
}
goto finalize;

34
src/equil/vcs_solve.cpp
View file

@ -222,15 +222,14 @@ VCS_SOLVE::~VCS_SOLVE()
void VCS_SOLVE::vcs_delete_memory()
{
size_t j;
size_t nspecies = m_numSpeciesTot;
for (j = 0; j < m_numPhases; j++) {
for (size_t j = 0; j < m_numPhases; j++) {
delete m_VolPhaseList[j];
m_VolPhaseList[j] = 0;
}
for (j = 0; j < nspecies; j++) {
for (size_t j = 0; j < nspecies; j++) {
delete m_speciesThermoList[j];
m_speciesThermoList[j] = 0;
}
@ -247,7 +246,6 @@ void VCS_SOLVE::vcs_delete_memory()
int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit)
{
int retn = 0, iconv = 0;
size_t nspecies0, nelements0, nphase0;
Cantera::clockWC tickTock;
int iprintTime = std::max(ipr, ip1);
@ -266,9 +264,9 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit)
* This function is called to create the private data
* using the public data.
*/
nspecies0 = vprob->nspecies + 10;
nelements0 = vprob->ne;
nphase0 = vprob->NPhase;
size_t nspecies0 = vprob->nspecies + 10;
size_t nelements0 = vprob->ne;
size_t nphase0 = vprob->NPhase;
vcs_initSizes(nspecies0, nelements0, nphase0);
@ -380,7 +378,6 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit)
int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
{
vcs_VolPhase* Vphase = 0;
const char* ser =
"vcs_pub_to_priv ERROR :ill defined interface -> bailout:\n\t";
@ -566,7 +563,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
* TPhInertMoles[] -> must be copied over here
*/
for (size_t iph = 0; iph < nph; iph++) {
Vphase = pub->VPhaseList[iph];
vcs_VolPhase* Vphase = pub->VPhaseList[iph];
TPhInertMoles[iph] = Vphase->totalMolesInert();
}
@ -626,7 +623,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
numPhSp[iph]++;
}
for (size_t iph = 0; iph < nph; iph++) {
Vphase = pub->VPhaseList[iph];
vcs_VolPhase* Vphase = pub->VPhaseList[iph];
if (numPhSp[iph] != Vphase->nSpecies()) {
plogf("%sNumber of species in phase %d, %s, doesn't match\n",
ser, iph, Vphase->PhaseName.c_str());
@ -703,7 +700,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
* It should point to the species thermo pointer in the private
* data space.
*/
Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
for (size_t k = 0; k < Vphase->nSpecies(); k++) {
vcs_SpeciesProperties* sProp = Vphase->speciesProperty(k);
size_t kT = Vphase->spGlobalIndexVCS(k);
@ -715,7 +712,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
* Specify the Activity Convention information
*/
for (size_t iph = 0; iph < nph; iph++) {
Vphase = m_VolPhaseList[iph];
vcs_VolPhase* Vphase = m_VolPhaseList[iph];
m_phaseActConvention[iph] = Vphase->p_activityConvention;
if (Vphase->p_activityConvention != 0) {
/*
@ -762,10 +759,8 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub)
int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub)
{
size_t kspec, k, i, j, iph;
string yo("vcs_prob_specify ERROR: ");
int retn = VCS_SUCCESS;
bool status_change = false;
m_temperature = pub->T;
m_pressurePA = pub->PresPA;
@ -779,8 +774,8 @@ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub)
m_tolmaj2 = 0.01 * m_tolmaj;
m_tolmin2 = 0.01 * m_tolmin;
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
k = m_speciesMapIndex[kspec];
for (size_t kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
size_t k = m_speciesMapIndex[kspec];
m_molNumSpecies_old[kspec] = pub->w[k];
m_molNumSpecies_new[kspec] = pub->mf[k];
m_feSpecies_old[kspec] = pub->m_gibbsSpecies[k];
@ -789,8 +784,8 @@ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub)
/*
* Transfer the element abundance goals to the solve object
*/
for (i = 0; i < m_numElemConstraints; i++) {
j = m_elementMapIndex[i];
for (size_t i = 0; i < m_numElemConstraints; i++) {
size_t j = m_elementMapIndex[i];
m_elemAbundancesGoal[i] = pub->gai[j];
}
@ -813,7 +808,8 @@ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub)
* condition.
*/
for (iph = 0; iph < m_numPhases; iph++) {
bool status_change = false;
for (size_t iph = 0; iph < m_numPhases; iph++) {
vcs_VolPhase* vPhase = m_VolPhaseList[iph];
vcs_VolPhase* pub_phase_ptr = pub->VPhaseList[iph];

531
src/equil/vcs_solve_TP.cpp
View file

File diff suppressed because it is too large Load diff

9
src/equil/vcs_species_thermo.cpp
View file

@ -124,9 +124,8 @@ double VCS_SPECIES_THERMO::GStar_R_calc(size_t kglob, double TKelvin,
double pres)
{
char yo[] = "VCS_SPECIES_THERMO::GStar_R_calc ";
double fe, T;
fe = G0_R_calc(kglob, TKelvin);
T = TKelvin;
double fe = G0_R_calc(kglob, TKelvin);
double T = TKelvin;
if (UseCanteraCalls) {
if (m_VCS_UnitsFormat != VCS_UNITS_MKS) {
throw Cantera::CanteraError("VCS_SPECIES_THERMO::GStar_R_calc",
@ -157,9 +156,9 @@ double VCS_SPECIES_THERMO::VolStar_calc(size_t kglob, double TKelvin,
double presPA)
{
char yo[] = "VCS_SPECIES_THERMO::VStar_calc ";
double vol, T;
double vol;
T = TKelvin;
double T = TKelvin;
if (UseCanteraCalls) {
if (m_VCS_UnitsFormat != VCS_UNITS_MKS) {
throw Cantera::CanteraError("VCS_SPECIES_THERMO::VolStar_calc",

95
src/equil/vcs_util.cpp
View file

@ -18,8 +18,7 @@ namespace VCSnonideal
#ifndef USE_MEMSET
void vcs_dzero(double* vector, int length)
{
int i;
for (i = 0; i < length; i++) {
for (int i = 0; i < length; i++) {
vector[i] = 0.0;
}
}
@ -28,8 +27,7 @@ void vcs_dzero(double* vector, int length)
#ifndef USE_MEMSET
void vcs_izero(int* vector, int length)
{
int i;
for (i = 0; i < length; i++) {
for (int i = 0; i < length; i++) {
vector[i] = 0;
}
}
@ -38,8 +36,7 @@ void vcs_izero(int* vector, int length)
#ifndef USE_MEMSET
void vcs_dcopy(double* const vec_to, const double* const vec_from, int length)
{
int i;
for (i = 0; i < length; i++) {
for (int i = 0; i < length; i++) {
vec_to[i] = vec_from[i];
}
}
@ -48,8 +45,7 @@ void vcs_dcopy(double* const vec_to, const double* const vec_from, int length)
#ifndef USE_MEMSET
void vcs_icopy(int* vec_to, int* vec_from, int length)
{
int i;
for (i = 0; i < length; i++) {
for (int i = 0; i < length; i++) {
vec_to[i] = vec_from[i];
}
}
@ -109,13 +105,12 @@ void vcs_vicopy(std::vector<int> &vec_to,
size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n)
{
size_t i;
size_t largest = j;
double big = x[j];
if (xSize) {
assert(xSize[j] > 0.0);
big *= xSize[j];
for (i = j + 1; i < n; ++i) {
for (size_t i = j + 1; i < n; ++i) {
assert(xSize[i] > 0.0);
if ((x[i] * xSize[i]) > big) {
largest = i;
@ -123,7 +118,7 @@ size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n)
}
}
} else {
for (i = j + 1; i < n; ++i) {
for (size_t i = j + 1; i < n; ++i) {
if (x[i] > big) {
largest = i;
big = x[i];
@ -135,12 +130,12 @@ size_t vcs_optMax(const double* x, const double* xSize, size_t j, size_t n)
int vcs_max_int(const int* vector, int length)
{
int i, retn;
int retn;
if (vector == NULL || length <= 0) {
return 0;
}
retn = vector[0];
for (i = 1; i < length; i++) {
for (int i = 1; i < length; i++) {
retn = std::max(retn, vector[i]);
}
return retn;
@ -148,30 +143,23 @@ int vcs_max_int(const int* vector, int length)
double vcsUtil_gasConstant(int mu_units)
{
double r;
switch (mu_units) {
case VCS_UNITS_KCALMOL:
r = Cantera::GasConst_cal_mol_K * 1e-3;
break;
return Cantera::GasConst_cal_mol_K * 1e-3;
case VCS_UNITS_UNITLESS:
r = 1.0;
break;
return 1.0;
case VCS_UNITS_KJMOL:
r = Cantera::GasConstant * 1e-6;
break;
return Cantera::GasConstant * 1e-6;
case VCS_UNITS_KELVIN:
r = 1.0;
break;
return 1.0;
case VCS_UNITS_MKS:
/* joules / kg-mol K = kg m2 / s2 kg-mol K */
r = Cantera::GasConstant;
break;
return Cantera::GasConstant;
default:
plogf("vcs_gasConstant error: uknown units: %d\n",
mu_units);
exit(EXIT_FAILURE);
}
return r;
}
void vcs_print_line(const char* string, int num)
@ -189,77 +177,66 @@ const char* vcs_speciesType_string(int speciesStatus, int length)
const char* sss;
switch (speciesStatus) {
case VCS_SPECIES_COMPONENT:
sss = "Component Species";
break;
return "Component Species";
case VCS_SPECIES_MAJOR:
sss = "Major Species";
break;
return "Major Species";
case VCS_SPECIES_MINOR:
sss = "Minor Species";
break;
return "Minor Species";
case VCS_SPECIES_ZEROEDPHASE:
if (length < 48) {
sss = "Set Zeroed-Phase";
return "Set Zeroed-Phase";
} else {
sss = "Purposely Zeroed-Phase Species (not in problem)";
return "Purposely Zeroed-Phase Species (not in problem)";
}
break;
case VCS_SPECIES_ZEROEDMS:
if (length < 23) {
sss = "Zeroed-MS Phase";
return "Zeroed-MS Phase";
} else {
sss = "Zeroed-MS Phase Species";
return "Zeroed-MS Phase Species";
}
break;
case VCS_SPECIES_ZEROEDSS:
if (length < 23) {
sss = "Zeroed-SS Phase";
return "Zeroed-SS Phase";
} else {
sss = "Zeroed-SS Phase Species";
return "Zeroed-SS Phase Species";
}
break;
case VCS_SPECIES_DELETED:
if (length < 22) {
sss = "Deleted Species";
return "Deleted Species";
} else if (length < 40) {
sss = "Deleted-Small Species";
return "Deleted-Small Species";
} else {
sss = "Deleted-Small Species in a MS phase";
return "Deleted-Small Species in a MS phase";
}
break;
case VCS_SPECIES_ACTIVEBUTZERO:
if (length < 47) {
sss = "Tmp Zeroed in MS";
return "Tmp Zeroed in MS";
} else {
sss = "Zeroed Species in an active MS phase (tmp)";
return "Zeroed Species in an active MS phase (tmp)";
}
break;
case VCS_SPECIES_STOICHZERO:
if (length < 56) {
sss = "Stoich Zeroed in MS";
return "Stoich Zeroed in MS";
} else {
sss = "Zeroed Species in an active MS phase (Stoich Constraint)";
return "Zeroed Species in an active MS phase (Stoich Constraint)";
}
break;
case VCS_SPECIES_INTERFACIALVOLTAGE:
if (length < 29) {
sss = "InterfaceVoltage";
return "InterfaceVoltage";
} else {
sss = "InterfaceVoltage Species";
return "InterfaceVoltage Species";
}
break;
default:
sss = "unknown species type";
return "unknown species type";
}
return sss;
}
void vcs_print_stringTrunc(const char* str, size_t space, int alignment)
{
size_t i, ls = 0, rs = 0;
size_t ls = 0, rs = 0;
size_t len = strlen(str);
if ((len) >= space) {
for (i = 0; i < space; i++) {
for (size_t i = 0; i < space; i++) {
plogf("%c", str[i]);
}
} else {
@ -272,13 +249,13 @@ void vcs_print_stringTrunc(const char* str, size_t space, int alignment)
rs = space - len - ls;
}
if (ls != 0) {
for (i = 0; i < ls; i++) {
for (size_t i = 0; i < ls; i++) {
plogf(" ");
}
}
plogf("%s", str);
if (rs != 0) {
for (i = 0; i < rs; i++) {
for (size_t i = 0; i < rs; i++) {
plogf(" ");
}
}