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

Changed names of a variable.

Browse source
This commit is contained in:
Harry Moffat 2008-04-21 19:27:07 +00:00
parent 75c11d40b3
commit 5ee0f16d88
11 changed files with 241 additions and 241 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

82
Cantera/src/equil/vcs_elem.cpp
View file

@ -29,7 +29,7 @@ void VCS_SOLVE::vcs_elab(void)
ga[j] = 0.0;
for (int i = 0; i < m_numSpeciesTot; ++i) {
if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
ga[j] += FormulaMatrix[j][i] * soln[i];
ga[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i];
}
}
}
@ -102,7 +102,7 @@ int VCS_SOLVE::vcs_elabcheck(int ibound) {
multisign = true;
}
if (eval != 0.0) {
scale = MAX(scale, fabs(eval * soln[kspec]));
scale = MAX(scale, fabs(eval * m_molNumSpecies_old[kspec]));
numNonZero++;
}
}
@ -144,7 +144,7 @@ void VCS_SOLVE::vcs_elabPhase(int iphase, double * const elemAbundPhase)
* Computes the elemental abundances vector for a single phase,
* elemAbundPhase[], and returns it through the argument list.
* The mole numbers of species are taken from the current value
* in soln[].
* in m_molNumSpecies_old[].
*************************************************************************/
{
int i, j;
@ -153,7 +153,7 @@ void VCS_SOLVE::vcs_elabPhase(int iphase, double * const elemAbundPhase)
for (i = 0; i < m_numSpeciesTot; ++i) {
if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (PhaseID[i] == iphase) {
elemAbundPhase[j] += FormulaMatrix[j][i] * soln[i];
elemAbundPhase[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i];
}
}
}
@ -194,7 +194,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
*
* ga Current element abundances
* gai Required elemental abundances
* soln Current mole number of species.
* m_molNumSpecies_old Current mole number of species.
* FormulaMatrix[][] Formular matrix of the species
* ne Number of elements
* nc Number of components.
@ -261,7 +261,7 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
double eval = FormulaMatrix[i][kspec];
if (eval > 0.0) {
soln[kspec] = gai[i] / eval;
m_molNumSpecies_old[kspec] = gai[i] / eval;
changed = true;
}
}
@ -283,9 +283,9 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
for (kspec = m_numComponents; kspec < m_numSpeciesTot; kspec++) {
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
double eval = FormulaMatrix[i][kspec];
diff -= eval * soln[kspec];
diff -= eval * m_molNumSpecies_old[kspec];
}
soln[compID] = MAX(0.0,diff/FormulaMatrix[i][compID]);
m_molNumSpecies_old[compID] = MAX(0.0,diff/FormulaMatrix[i][compID]);
changed = true;
}
}
@ -314,18 +314,18 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
double atomComp = FormulaMatrix[i][kspec];
if (atomComp > 0.0) {
double maxPermissible = gai[i] / atomComp;
if (soln[kspec] > maxPermissible) {
if (m_molNumSpecies_old[kspec] > maxPermissible) {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 3) {
plogf(" --- vcs_elcorr: Reduced species %s from %g to %g due to %s max bounds constraint\n",
SpName[kspec].c_str(), soln[kspec], maxPermissible, ElName[i].c_str());
SpName[kspec].c_str(), m_molNumSpecies_old[kspec], maxPermissible, ElName[i].c_str());
}
#endif
soln[kspec] = maxPermissible;
m_molNumSpecies_old[kspec] = maxPermissible;
changed = true;
if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF) {
soln[kspec] = 0.0;
if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF) {
m_molNumSpecies_old[kspec] = 0.0;
if (SSPhase[kspec]) {
spStatus[kspec] = VCS_SPECIES_ZEROEDSS;
} else {
@ -372,8 +372,8 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
*/
par = 0.5;
for (i = 0; i < m_numComponents; ++i) {
if (soln[i] > 0.0) {
xx = -x[i] / soln[i];
if (m_molNumSpecies_old[i] > 0.0) {
xx = -x[i] / m_molNumSpecies_old[i];
if (par < xx) par = xx;
}
}
@ -385,27 +385,27 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
retn = 2;
par *= 0.9999;
for (i = 0; i < m_numComponents; ++i) {
double tmp = soln[i] + par * x[i];
double tmp = m_molNumSpecies_old[i] + par * x[i];
if (tmp > 0.0) {
soln[i] = tmp;
m_molNumSpecies_old[i] = tmp;
} else {
if (SSPhase[i]) {
soln[i] = 0.0;
m_molNumSpecies_old[i] = 0.0;
} else {
soln[i] = soln[i] * 0.0001;
m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001;
}
}
}
} else {
for (i = 0; i < m_numComponents; ++i) {
double tmp = soln[i] + x[i];
double tmp = m_molNumSpecies_old[i] + x[i];
if (tmp > 0.0) {
soln[i] = tmp;
m_molNumSpecies_old[i] = tmp;
} else {
if (SSPhase[i]) {
soln[i] = 0.0;
m_molNumSpecies_old[i] = 0.0;
} else {
soln[i] = soln[i] * 0.0001;
m_molNumSpecies_old[i] = m_molNumSpecies_old[i] * 0.0001;
}
}
}
@ -466,15 +466,15 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
}
}
if (its > 0) xx /= its;
soln[kspec] += xx;
soln[kspec] = MAX(soln[kspec], 1.0E-10);
m_molNumSpecies_old[kspec] += xx;
m_molNumSpecies_old[kspec] = MAX(m_molNumSpecies_old[kspec], 1.0E-10);
/*
* If we are dealing with a deleted species, then
* we need to reinsert it into the active list.
*/
if (kspec >= m_numSpeciesRdc) {
vcs_reinsert_deleted(kspec);
soln[m_numSpeciesRdc - 1] = xx;
m_molNumSpecies_old[m_numSpeciesRdc - 1] = xx;
vcs_elab();
goto L_CLEANUP;
}
@ -493,9 +493,9 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
if (ga[i] > 0.0) {
if (FormulaMatrix[i][kspec] < 0.0) {
soln[kspec] -= ga[i] / FormulaMatrix[i][kspec] ;
if (soln[kspec] < 0.0) {
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] -= ga[i] / FormulaMatrix[i][kspec] ;
if (m_molNumSpecies_old[kspec] < 0.0) {
m_molNumSpecies_old[kspec] = 0.0;
}
vcs_elab();
break;
@ -503,9 +503,9 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
}
if (ga[i] < 0.0) {
if (FormulaMatrix[i][kspec] > 0.0) {
soln[kspec] -= ga[i] / FormulaMatrix[i][kspec];
if (soln[kspec] < 0.0) {
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] -= ga[i] / FormulaMatrix[i][kspec];
if (m_molNumSpecies_old[kspec] < 0.0) {
m_molNumSpecies_old[kspec] = 0.0;
}
vcs_elab();
break;
@ -531,26 +531,26 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
if (dev < 0.0) {
if (FormulaMatrix[i][kspec] < 0.0) {
if (soln[kspec] > 0.0) {
if (m_molNumSpecies_old[kspec] > 0.0) {
useZeroed = false;
}
}
} else {
if (FormulaMatrix[i][kspec] > 0.0) {
if (soln[kspec] > 0.0) {
if (m_molNumSpecies_old[kspec] > 0.0) {
useZeroed = false;
}
}
}
}
for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
if (soln[kspec] > 0.0 || useZeroed) {
if (m_molNumSpecies_old[kspec] > 0.0 || useZeroed) {
if (dev < 0.0) {
if (FormulaMatrix[i][kspec] < 0.0) {
double delta = dev / FormulaMatrix[i][kspec] ;
soln[kspec] += delta;
if (soln[kspec] < 0.0) {
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] += delta;
if (m_molNumSpecies_old[kspec] < 0.0) {
m_molNumSpecies_old[kspec] = 0.0;
}
vcs_elab();
break;
@ -559,9 +559,9 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
if (dev > 0.0) {
if (FormulaMatrix[i][kspec] > 0.0) {
double delta = dev / FormulaMatrix[i][kspec] ;
soln[kspec] += delta;
if (soln[kspec] < 0.0) {
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] += delta;
if (m_molNumSpecies_old[kspec] < 0.0) {
m_molNumSpecies_old[kspec] = 0.0;
}
vcs_elab();
break;

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

@ -52,7 +52,7 @@ namespace VCSnonideal {
int nrxn = m_numRxnTot;
vcs_VolPhase *Vphase = 0;
double *molNum = VCS_DATA_PTR(soln);
double *molNum = VCS_DATA_PTR(m_molNumSpecies_old);
double TMolesMultiphase;
double *xtphMax = VCS_DATA_PTR(TmpPhase);
double *xtphMin = VCS_DATA_PTR(TmpPhase2);
@ -504,7 +504,7 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix,
vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_curr),
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles)));
plogendl();
}

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

@ -113,7 +113,7 @@ void VCS_SOLVE::vcs_nondim_TP(void) {
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
for (i = 0; i < m_numSpeciesTot; ++i) {
if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
soln[i] *= 1.0E3;
m_molNumSpecies_old[i] *= 1.0E3;
}
}
for (i = 0; i < m_numElemConstraints; ++i) {
@ -153,7 +153,7 @@ void VCS_SOLVE::vcs_redim_TP(void)
if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
for (i = 0; i < m_numSpeciesTot; ++i) {
if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
soln[i] /= 1.0E3;
m_molNumSpecies_old[i] /= 1.0E3;
}
}
for (i = 0; i < m_numElemConstraints; ++i) {

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

@ -186,7 +186,7 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
double sum = 0.0;
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
sum += fabs(soln[kspec]);
sum += fabs(m_molNumSpecies_old[kspec]);
}
}
if (fabs(sum) < 1.0E-6) {
@ -196,9 +196,9 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
retn = vcs_evalSS_TP(0, 0, T, pres);
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
soln[kspec] = - m_SSfeSpecies[kspec];
m_molNumSpecies_old[kspec] = - m_SSfeSpecies[kspec];
} else {
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] = 0.0;
}
}
}
@ -256,7 +256,7 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl)
// zero to start with, set them back to zero here
if (modifiedSoln) {
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] = 0.0;
}
}
return VCS_SUCCESS;

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

@ -62,7 +62,7 @@ int VCS_SOLVE::vcs_report(int iconv)
for (i = 0; i < nspecies; ++i) {
sortindex[i] = i;
xy[i] = soln[i];
xy[i] = m_molNumSpecies_old[i];
}
/*
* Sort the XY vector, the mole fraction vector,
@ -105,7 +105,7 @@ int VCS_SOLVE::vcs_report(int iconv)
* Calculate some quantities that may need updating
*/
vcs_tmoles();
Vol = vcs_VolTotal(T, Pres, VCS_DATA_PTR(soln), VCS_DATA_PTR(VolPM));
Vol = vcs_VolTotal(T, Pres, VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM));
plogf("\t\tTemperature = %15.2g Kelvin\n", T);
plogf("\t\tPressure = %15.5g Atmos\n", Pres);
@ -122,7 +122,7 @@ int VCS_SOLVE::vcs_report(int iconv)
for (i = 0; i < m_numComponents; ++i) {
plogf(" %-12.12s", SpName[i].c_str());
print_space(13);
plogf("%14.7E %14.7E %12.4E", soln[i], m_molNumSpecies_new[i], m_feSpecies_curr[i]);
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i], m_molNumSpecies_new[i], m_feSpecies_curr[i]);
plogf(" %3d", SpeciesUnknownType[i]);
plogf("\n");
}
@ -132,11 +132,11 @@ int VCS_SOLVE::vcs_report(int iconv)
print_space(13);
if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_MOLNUM) {
plogf("%14.7E %14.7E %12.4E", soln[l], m_molNumSpecies_new[l], m_feSpecies_curr[l]);
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[l], m_molNumSpecies_new[l], m_feSpecies_curr[l]);
plogf(" MolNum ");
} else if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_curr[l]);
plogf(" Voltage = %14.7E", soln[l]);
plogf(" Voltage = %14.7E", m_molNumSpecies_old[l]);
} else {
plogf("we have a problem\n");
exit(-1);
@ -161,7 +161,7 @@ int VCS_SOLVE::vcs_report(int iconv)
for (kspec = m_numSpeciesRdc; kspec < nspecies; ++kspec) {
plogf(" %-12.12s", SpName[kspec].c_str());
plogf(" %14.7E %14.7E %12.4E",
soln[kspec], m_molNumSpecies_new[kspec], dg[kspec]);
m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], dg[kspec]);
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" Mol_Num");
} else if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
@ -193,7 +193,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" | |\n");
plogf(" NonComponent | Moles |");
for (j = 0; j < m_numComponents; j++) {
plogf(" %10.3g", soln[j]);
plogf(" %10.3g", m_molNumSpecies_old[j]);
}
plogf(" | DG/RT Rxn |\n");
print_line("-", m_numComponents*10 + 45);
@ -201,7 +201,7 @@ int VCS_SOLVE::vcs_report(int iconv)
int irxn = ir[i];
plogf(" %3d ", irxn);
plogf("%-10.10s", SpName[irxn].c_str());
plogf("|%10.3g |", soln[irxn]);
plogf("|%10.3g |", m_molNumSpecies_old[irxn]);
for (j = 0; j < m_numComponents; j++) {
plogf(" %6.2f", sc[i][j]);
}
@ -255,7 +255,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %10.3g", gaPhase[j]);
gaTPhase[j] += gaPhase[j];
}
gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(soln),
gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(m_feSpecies_curr));
gibbsTotal += gibbsPhase;
plogf(" | %18.11E |\n", gibbsPhase);
@ -280,7 +280,7 @@ int VCS_SOLVE::vcs_report(int iconv)
* energy of zero
*/
g = vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_curr),
g = vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr),
VCS_DATA_PTR(TPhMoles));
plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g);
if (inertYes)
@ -312,7 +312,7 @@ int VCS_SOLVE::vcs_report(int iconv)
l = sortindex[i];
int pid = PhaseID[l];
plogf(" %-12.12s", SpName[l].c_str());
plogf(" %14.7E ", soln[l]);
plogf(" %14.7E ", m_molNumSpecies_old[l]);
plogf("%14.7E ", m_SSfeSpecies[l]);
plogf("%14.7E ", log(ActCoeff[l]));
double tpmoles = TPhMoles[pid];
@ -322,8 +322,8 @@ int VCS_SOLVE::vcs_report(int iconv)
if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
lx = 0.0;
} else {
if (tpmoles > 0.0 && soln[l] > 0.0) {
lx = log(soln[l]) - log(tpmoles);
if (tpmoles > 0.0 && m_molNumSpecies_old[l] > 0.0) {
lx = log(m_molNumSpecies_old[l]) - log(tpmoles);
} else {
lx = m_feSpecies_curr[l] - m_SSfeSpecies[l] - log(ActCoeff[l]) + SpecLnMnaught[l];
}

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

@ -78,7 +78,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
kspec = ir[irxn];
dnPhase_irxn = DnPhase[irxn];
if (soln[kspec] == 0.0 && (! SSPhase[kspec])) {
if (m_molNumSpecies_old[kspec] == 0.0 && (! SSPhase[kspec])) {
/* *******************************************************************/
/* **** MULTISPECIES PHASE WITH total moles equal to zero ************/
/* *******************************************************************/
@ -114,7 +114,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
#ifdef DEBUG_MODE
sprintf(ANOTE,"Skipped: converged DG = %11.3E\n", dg[irxn]);
plogf(" --- "); plogf("%-12.12s", SpName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", soln[kspec], ds[kspec], ANOTE);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], ds[kspec], ANOTE);
#endif
continue;
}
@ -127,7 +127,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
sprintf(ANOTE,"Skipped: IC = %3d and DG >0: %11.3E\n",
spStatus[irxn], dg[irxn]);
plogf(" --- "); plogf("%-12.12s", SpName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", soln[kspec], ds[kspec], ANOTE);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], ds[kspec], ANOTE);
#endif
continue;
}
@ -135,9 +135,9 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
* Start of the regular processing
*/
if (SSPhase[kspec]) s = 0.0;
else s = 1.0 / soln[kspec];
else s = 1.0 / m_molNumSpecies_old[kspec];
for (j = 0; j < m_numComponents; ++j) {
if (! SSPhase[j]) s += SQUARE(sc[irxn][j]) / soln[j];
if (! SSPhase[j]) s += SQUARE(sc[irxn][j]) / m_molNumSpecies_old[j];
}
for (j = 0; j < NPhase; j++) {
if (! (VPhaseList[j])->SingleSpecies) {
@ -160,11 +160,11 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
* will zero out first.
*/
if (dg[irxn] > 0.0) {
dss = soln[kspec];
dss = m_molNumSpecies_old[kspec];
k = kspec;
for (j = 0; j < m_numComponents; ++j) {
if (sc[irxn][j] > 0.0) {
xx = soln[j] / sc[irxn][j];
xx = m_molNumSpecies_old[j] / sc[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
@ -176,7 +176,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
dss = 1.0e10;
for (j = 0; j < m_numComponents; ++j) {
if (sc[irxn][j] < 0.0) {
xx = -soln[j] / sc[irxn][j];
xx = -m_molNumSpecies_old[j] / sc[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
@ -193,13 +193,13 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
* added back into the component species.
*/
if (dss != 0.0) {
soln[kspec] += dss;
m_molNumSpecies_old[kspec] += dss;
TPhMoles[PhaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
soln[j] += dss * sc[irxn][j];
m_molNumSpecies_old[j] += dss * sc[irxn][j];
TPhMoles[PhaseID[j]] += dss * sc[irxn][j];
}
soln[k] = 0.0;
m_molNumSpecies_old[k] = 0.0;
TPhMoles[PhaseID[k]] = 0.0;
#ifdef DEBUG_MODE
plogf(" --- vcs_st2 Special section to delete ");
@ -219,7 +219,7 @@ int VCS_SOLVE::vcs_rxn_adj_cg(void)
} /* End of regular processing */
#ifdef DEBUG_MODE
plogf(" --- "); plogf("%-12.12s", SpName[kspec].c_str());
plogf(" %12.4E %12.4E | %s\n", soln[kspec], ds[kspec], ANOTE);
plogf(" %12.4E %12.4E | %s\n", m_molNumSpecies_old[kspec], ds[kspec], ANOTE);
#endif
} /* End of loop over non-component stoichiometric formation reactions */
@ -409,7 +409,7 @@ double VCS_SOLVE::vcs_line_search(int irxn, double dx_orig)
const int MAXITS = 10;
double dx = dx_orig;
double *sc_irxn = sc[irxn];
double *molNumBase = VCS_DATA_PTR(soln);
double *molNumBase = VCS_DATA_PTR(m_molNumSpecies_old);
double *acBase = VCS_DATA_PTR(ActCoeff0);
double *ac = VCS_DATA_PTR(ActCoeff);
double molSum = 0.0;

30
Cantera/src/equil/vcs_setMolesLinProg.cpp
View file

@ -92,13 +92,13 @@ int VCS_SOLVE::vcs_setMolesLinProg() {
for (ik = 0; ik < m_numSpeciesTot; ik++) {
if (SpeciesUnknownType[ik] != VCS_SPECIES_INTERFACIALVOLTAGE) {
soln[ik] = MAX(0.0, soln[ik]);
m_molNumSpecies_old[ik] = MAX(0.0, m_molNumSpecies_old[ik]);
}
}
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
printProgress(SpName, soln, m_SSfeSpecies);
printProgress(SpName, m_molNumSpecies_old, m_SSfeSpecies);
}
#endif
@ -122,7 +122,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() {
}
/*
* Now find the optimized basis that spans the stoichiometric
* coefficient matrix, based on the current composition, soln[]
* coefficient matrix, based on the current composition, m_molNumSpecies_old[]
* We also calculate sc[][], the reaction matrix.
*/
retn = vcs_basopt(FALSE, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa),
@ -155,7 +155,7 @@ int VCS_SOLVE::vcs_setMolesLinProg() {
// idir < 0 implies decreasing the current species
idir = (dg_rt < 0.0 ? 1 : -1);
if (idir < 0) {
dxi_min = soln[ik];
dxi_min = m_molNumSpecies_old[ik];
}
for (jcomp = 0; jcomp < m_numComponents; jcomp++) {
@ -164,11 +164,11 @@ int VCS_SOLVE::vcs_setMolesLinProg() {
// set max change in progress variable by
// non-negativity requirement
if (nu*idir < 0) {
delta_xi = fabs(soln[jcomp]/nu);
delta_xi = fabs(m_molNumSpecies_old[jcomp]/nu);
// if a component has nearly zero moles, redo
// with a new set of components
if (!redo) {
if (delta_xi < 1.0e-10 && (soln[ik] >= 1.0E-10)) {
if (delta_xi < 1.0e-10 && (m_molNumSpecies_old[ik] >= 1.0E-10)) {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf(" --- Component too small: %s\n", SpName[jcomp].c_str());
@ -185,17 +185,17 @@ int VCS_SOLVE::vcs_setMolesLinProg() {
// we are zeroing components and species on every step.
// Redo the iteration, if a component went from positive to zero on this step.
double dsLocal = idir*dxi_min;
soln[ik] += dsLocal;
soln[ik] = MAX(0.0, soln[ik]);
m_molNumSpecies_old[ik] += dsLocal;
m_molNumSpecies_old[ik] = MAX(0.0, m_molNumSpecies_old[ik]);
for (jcomp = 0; jcomp < m_numComponents; jcomp++) {
bool full = false;
if (soln[jcomp] > 1.0E-15) {
if (m_molNumSpecies_old[jcomp] > 1.0E-15) {
full = true;
}
soln[jcomp] += sc_irxn[jcomp] * dsLocal;
soln[jcomp] = MAX(0.0, soln[jcomp]);
m_molNumSpecies_old[jcomp] += sc_irxn[jcomp] * dsLocal;
m_molNumSpecies_old[jcomp] = MAX(0.0, m_molNumSpecies_old[jcomp]);
if (full) {
if (soln[jcomp] < 1.0E-60) {
if (m_molNumSpecies_old[jcomp] < 1.0E-60) {
redo = true;
}
}
@ -204,18 +204,18 @@ int VCS_SOLVE::vcs_setMolesLinProg() {
// set the moles of the phase objects to match
// updateMixMoles();
// Update the phase objects with the contents of the soln vector
// Update the phase objects with the contents of the m_molNumSpecies_old vector
// vcs_updateVP(0);
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
printProgress(SpName, soln, m_SSfeSpecies);
printProgress(SpName, m_molNumSpecies_old, m_SSfeSpecies);
}
#endif
}
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl == 1) {
printProgress(SpName, soln, m_SSfeSpecies);
printProgress(SpName, m_molNumSpecies_old, m_SSfeSpecies);
plogf(" --- setInitialMoles end\n");
}
#endif

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

@ -112,7 +112,7 @@ namespace VCSnonideal {
m_feSpecies_curr.resize(nspecies0, 0.0);
m_SSfeSpecies.resize(nspecies0, 0.0);
m_feSpecies_new.resize(nspecies0, 0.0);
soln.resize(nspecies0, 0.0);
m_molNumSpecies_old.resize(nspecies0, 0.0);
SpeciesUnknownType.resize(nspecies0, VCS_SPECIES_TYPE_MOLNUM);
@ -543,10 +543,10 @@ namespace VCSnonideal {
* w[] -> Copy the equilibrium mole number estimate if it exists.
*/
if (pub->w.size() != 0) {
vcs_vdcopy(soln, pub->w, nspecies);
vcs_vdcopy(m_molNumSpecies_old, pub->w, nspecies);
} else {
iest = -1;
vcs_dzero(VCS_DATA_PTR(soln), nspecies);
vcs_dzero(VCS_DATA_PTR(m_molNumSpecies_old), nspecies);
}
/*
@ -560,7 +560,7 @@ namespace VCSnonideal {
gai[j] = 0.0;
for (kspec = 0; kspec < nspecies; kspec++) {
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
gai[j] += FormulaMatrix[j][kspec] * soln[kspec];
gai[j] += FormulaMatrix[j][kspec] * m_molNumSpecies_old[kspec];
}
}
}
@ -787,7 +787,7 @@ namespace VCSnonideal {
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
k = ind[kspec];
soln[kspec] = pub->w[k];
m_molNumSpecies_old[kspec] = pub->w[k];
m_molNumSpecies_new[kspec] = pub->mf[k];
m_feSpecies_curr[kspec] = pub->m_gibbsSpecies[k];
}
@ -903,7 +903,7 @@ namespace VCSnonideal {
int k1 = 0;
vcs_tmoles();
Vol = vcs_VolTotal(T, Pres, VCS_DATA_PTR(soln), VCS_DATA_PTR(VolPM));
Vol = vcs_VolTotal(T, Pres, VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(VolPM));
for (i = 0; i < m_numSpeciesTot; ++i) {
/*
@ -919,10 +919,10 @@ namespace VCSnonideal {
* - Switch the species data back from K1 into I
*/
if (pub->SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
pub->w[i] = soln[k1];
pub->w[i] = m_molNumSpecies_old[k1];
} else {
pub->w[i] = 0.0;
plogf("voltage species = %g\n", soln[k1]);
plogf("voltage species = %g\n", m_molNumSpecies_old[k1]);
}
pub->mf[i] = m_molNumSpecies_new[k1];
pub->m_gibbsSpecies[i] = m_feSpecies_curr[k1];
@ -956,7 +956,7 @@ namespace VCSnonideal {
if (pubPhase->m_phiVarIndex == k) {
k1 = vPhase->IndSpecies[k];
double tmp = soln[k1];
double tmp = m_molNumSpecies_old[k1];
if (! vcs_doubleEqual( pubPhase->electricPotential() , tmp)) {
plogf("We have an inconsistency in voltage, %g, %g\n",
pubPhase->electricPotential(), tmp);

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

@ -459,10 +459,10 @@ public:
//! Total moles of the species
/*!
* soln[k] = Total number of moles of the kth species.
* Total number of moles of the kth species.
* Length = Total number of species = m
*/
std::vector<double> soln;
std::vector<double> m_molNumSpecies_old;
//! Specifies the species unknown type
/*!

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

@ -294,7 +294,7 @@ namespace VCSnonideal {
plogf(" 0");
}
print_space(47-m_numElemConstraints*3);
plogf("%12.5E %12.5E", RT * m_SSfeSpecies[i], soln[i]);
plogf("%12.5E %12.5E", RT * m_SSfeSpecies[i], m_molNumSpecies_old[i]);
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" Mol_Num");
} else if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
@ -307,12 +307,12 @@ namespace VCSnonideal {
}
for (i = 0; i < m_numSpeciesTot; ++i) {
if (soln[i] < 0.0) {
if (m_molNumSpecies_old[i] < 0.0) {
plogf("On Input species %-12s has a "
"negative MF, setting it small",
SpName[i].c_str());
plogendl();
soln[i] = VCS_DELETE_SPECIES_CUTOFF;
m_molNumSpecies_old[i] = VCS_DELETE_SPECIES_CUTOFF;
}
}
@ -325,7 +325,7 @@ namespace VCSnonideal {
/* ******************************************* */
/* **** EVALUATE ALL CHEMICAL POTENTIALS ***** */
/* ******************************************* */
vcs_dfe(VCS_DATA_PTR(soln), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
/*
* HKM -> If there was a machine estimate, we used to branch
* to the code segment which determined whether we needed a
@ -386,7 +386,7 @@ namespace VCSnonideal {
}
#endif
vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx));
vcs_dfe(VCS_DATA_PTR(soln), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
}
#ifdef DEBUG_MODE
else {
@ -422,7 +422,7 @@ namespace VCSnonideal {
* We have already evaluated the major non-components
*/
if (uptodate_minors == FALSE) {
vcs_dfe(VCS_DATA_PTR(soln), 0, 1, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 1, 0, m_numSpeciesRdc);
vcs_deltag(1, false);
}
uptodate_minors = TRUE;
@ -561,7 +561,7 @@ namespace VCSnonideal {
plogf(" --- %s currently zeroed (SpStatus=%-2d):",
SpName[kspec].c_str(), spStatus[irxn]);
plogf("%3d DG = %11.4E WT = %11.4E W = %11.4E DS = %11.4E\n",
irxn, dg[irxn], m_molNumSpecies_new[kspec], soln[kspec], ds[kspec]);
irxn, dg[irxn], m_molNumSpecies_new[kspec], m_molNumSpecies_old[kspec], ds[kspec]);
}
#endif
// HKM Alternative is to not allow ds[] = 0.0 phases
@ -572,7 +572,7 @@ namespace VCSnonideal {
// This could change in the future.
//if (dg[irxn] >= 0.0 || ds[kspec] <= 0.0) {
if (dg[irxn] >= 0.0 ) {
m_molNumSpecies_new[kspec] = soln[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
ds[kspec] = 0.0;
resurrect = false;
#ifdef DEBUG_MODE
@ -626,17 +626,17 @@ namespace VCSnonideal {
if (ds[kspec] > 0.0) {
dx = ds[kspec] * 0.01;
m_molNumSpecies_new[kspec] = soln[kspec] + dx;
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
} else {
m_molNumSpecies_new[kspec] = TMoles * VCS_DELETE_PHASE_CUTOFF * 10.;
dx = m_molNumSpecies_new[kspec] - soln[kspec];
dx = m_molNumSpecies_new[kspec] - m_molNumSpecies_old[kspec];
}
ds[kspec] = dx;
#ifdef DEBUG_MODE
sprintf(ANOTE, "Born:IC=-1 to IC=1:DG=%11.4E", dg[irxn]);
#endif
} else {
m_molNumSpecies_new[kspec] = soln[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
ds[kspec] = 0.0;
dx = 0.0;
}
@ -649,7 +649,7 @@ namespace VCSnonideal {
* to minor species.
*/
if (iti != 0) {
m_molNumSpecies_new[kspec] = soln[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
ds[kspec] = 0.0;
dx = 0.0;
#ifdef DEBUG_MODE
@ -657,7 +657,7 @@ namespace VCSnonideal {
if (vcs_debug_print_lvl >= 2) {
plogf(" --- "); plogf("%-12s", SpName[kspec].c_str());
plogf("%3d%11.4E%11.4E%11.4E | %s",
spStatus[irxn], soln[kspec], m_molNumSpecies_new[kspec],
spStatus[irxn], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
ds[kspec], ANOTE);
plogendl();
}
@ -733,7 +733,7 @@ namespace VCSnonideal {
* irxn loop if it is superconverged.
*/
if (fabs(dg[irxn]) <= tolmaj2) {
m_molNumSpecies_new[kspec] = soln[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
ds[kspec] = 0.0;
dx = 0.0;
#ifdef DEBUG_MODE
@ -741,7 +741,7 @@ namespace VCSnonideal {
if (vcs_debug_print_lvl >= 2) {
plogf(" --- "); plogf("%-12s", SpName[kspec].c_str());
plogf("%3d%11.4E%11.4E%11.4E | %s",
spStatus[irxn], soln[kspec], m_molNumSpecies_new[kspec],
spStatus[irxn], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
ds[kspec], ANOTE);
plogendl();
}
@ -770,7 +770,7 @@ namespace VCSnonideal {
/*
* Form a tentative value of the new species moles
*/
m_molNumSpecies_new[kspec] = soln[kspec] + dx;
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
/*
* Check for non-positive mole fraction of major species.
* If we find one, we branch to a section below. Then,
@ -799,9 +799,9 @@ namespace VCSnonideal {
* major species in multispecies phases.
* Decrease its concentration by a factor of 10.
*/
dx = -0.9 * soln[kspec];
dx = -0.9 * m_molNumSpecies_old[kspec];
ds[kspec] = dx;
m_molNumSpecies_new[kspec] = soln[kspec] + dx;
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
/*
* Change major to minor if the current species
* has a mole number that is less than 1/100 of the
@ -832,7 +832,7 @@ namespace VCSnonideal {
* Calculate a dx that will wipe out the
* moles in the phase.
*/
dx = -soln[kspec];
dx = -m_molNumSpecies_old[kspec];
/*
* Calculate an update that doesn't create a negative mole
* number for a component species. Actually, restrict this
@ -841,15 +841,15 @@ namespace VCSnonideal {
*/
for (j = 0; j < m_numComponents; ++j) {
if (sc_irxn[j] != 0.0) {
wx[j] = soln[j] + sc_irxn[j] * dx;
if (wx[j] <= soln[j] * 0.01 - 1.0E-150) {
dx = MAX(dx, soln[j] * -0.99 / sc_irxn[j]);
wx[j] = m_molNumSpecies_old[j] + sc_irxn[j] * dx;
if (wx[j] <= m_molNumSpecies_old[j] * 0.01 - 1.0E-150) {
dx = MAX(dx, m_molNumSpecies_old[j] * -0.99 / sc_irxn[j]);
}
} else {
wx[j] = soln[j];
wx[j] = m_molNumSpecies_old[j];
}
}
m_molNumSpecies_new[kspec] = soln[kspec] + dx;
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx;
if (m_molNumSpecies_new[kspec] > 0.0) {
ds[kspec] = dx;
#ifdef DEBUG_MODE
@ -874,9 +874,9 @@ namespace VCSnonideal {
* to eliminate the phase in this special section
* outside the main loop.
*/
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] = 0.0;
for (j = 0; j < m_numComponents; ++j) {
soln[j] = wx[j];
m_molNumSpecies_old[j] = wx[j];
}
/*
* Change the total number of moles in all phases due to
@ -896,7 +896,7 @@ namespace VCSnonideal {
* set of reactions being considered. The set of reactions
* is determined by the value of iti.
*/
vcs_dfe(VCS_DATA_PTR(soln), 0, iti, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, iti, 0, m_numSpeciesRdc);
vcs_deltag(iti, false);
/*
* Redefine the starting conditions for noncomponents
@ -923,7 +923,7 @@ namespace VCSnonideal {
if (im && iti != 0) {
goto L_EQUILIB_CHECK;
}
m_molNumSpecies_new[kspec] = soln[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec];
ds[kspec] = 0.0;
dx = 0.0;
}
@ -935,7 +935,7 @@ namespace VCSnonideal {
/*
* Skip the line search if we are birthing a species
*/
if (dx != 0.0 && (soln[kspec] > 0.0) &&
if (dx != 0.0 && (m_molNumSpecies_old[kspec] > 0.0) &&
(SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
double dx_old = dx;
#ifdef DEBUG_MODE
@ -988,10 +988,10 @@ namespace VCSnonideal {
#endif
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
m_molNumSpecies_new[kspec] = soln[kspec] + ds[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + ds[kspec];
plogf(" --- "); plogf("%-12.12s", SpName[kspec].c_str());
plogf("%3d%11.4E%11.4E%11.4E | %s",
spStatus[irxn], soln[kspec], m_molNumSpecies_new[kspec],
spStatus[irxn], m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec],
ds[kspec], ANOTE);
plogendl();
}
@ -1004,7 +1004,7 @@ namespace VCSnonideal {
for (k = 0; k < m_numComponents; k++) {
plogf(" --- "); plogf("%-12.12s", SpName[k].c_str());
plogf(" c%11.4E%11.4E%11.4E |\n",
soln[k], soln[k]+ds[k], ds[k]);
m_molNumSpecies_old[k], m_molNumSpecies_old[k]+ds[k], ds[k]);
}
plogf(" "); vcs_print_line("-", 80);
plogf(" --- Finished Main Loop");
@ -1020,8 +1020,8 @@ namespace VCSnonideal {
*/
par = 0.5;
for (k = 0; k < m_numComponents; ++k) {
if (soln[k] > 0.0) {
xx = -ds[k] / soln[k];
if (m_molNumSpecies_old[k] > 0.0) {
xx = -ds[k] / m_molNumSpecies_old[k];
if (par < xx) {
par = xx;
#ifdef DEBUG_MODE
@ -1072,7 +1072,7 @@ namespace VCSnonideal {
* consistent with a new estimate of the state of the system.
*/
for (kspec = 0; kspec < m_numSpeciesTot; ++kspec) {
m_molNumSpecies_new[kspec] = soln[kspec] + ds[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + ds[kspec];
if (m_molNumSpecies_new[kspec] < 0.0 && (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE)) {
plogf("vcs_solve_TP: ERROR on step change wt[%d:%s]: %g < 0.0",
kspec, SpName[kspec].c_str(), m_molNumSpecies_new[kspec]);
@ -1117,7 +1117,7 @@ namespace VCSnonideal {
/* *************************************************************** */
if (printDetails) {
plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_old),
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
VCS_DATA_PTR(TPhMoles)));
plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E",
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr),
@ -1138,13 +1138,13 @@ namespace VCSnonideal {
plogf(" FINAL MOLES INIT_DEL_G/RT TENT_DEL_G/RT FINAL_DELTA_G/RT\n");
for (i = 0; i < m_numComponents; ++i) {
plogf(" --- %-12.12s", SpName[i].c_str());
plogf(" %14.6E %14.6E %14.6E\n", soln[i], soln[i] + ds[i], m_molNumSpecies_new[i]);
plogf(" %14.6E %14.6E %14.6E\n", m_molNumSpecies_old[i], m_molNumSpecies_old[i] + ds[i], m_molNumSpecies_new[i]);
}
for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
irxn = kspec - m_numComponents;
plogf(" --- %-12.12s", SpName[kspec].c_str());
plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", spStatus[irxn],
soln[kspec], soln[kspec]+ds[kspec], m_molNumSpecies_new[kspec], dgl[irxn],
m_molNumSpecies_old[kspec], m_molNumSpecies_old[kspec]+ds[kspec], m_molNumSpecies_new[kspec], dgl[irxn],
m_deltaGRxn_tmp[irxn], dg[irxn]);
}
print_space(26);
@ -1180,14 +1180,14 @@ namespace VCSnonideal {
plogf(" Mu/RT Init_Del_G/RT Delta_G/RT\n");
for (i = 0; i < m_numComponents; ++i) {
plogf(" --- %-12.12s", SpName[i].c_str()); plogf(" ");
plogf("%14.6E%14.6E%14.6E%14.6E\n", soln[i],
plogf("%14.6E%14.6E%14.6E%14.6E\n", m_molNumSpecies_old[i],
m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_curr[i]);
}
for (i = m_numComponents; i < m_numSpeciesRdc; ++i) {
l1 = i - m_numComponents;
plogf(" --- %-12.12s", SpName[i].c_str());
plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
spStatus[l1], soln[i],
spStatus[l1], m_molNumSpecies_old[i],
m_molNumSpecies_new[i], m_feSpecies_old[i], m_feSpecies_curr[i],
dgl[l1], dg[l1]);
}
@ -1195,7 +1195,7 @@ namespace VCSnonideal {
l1 = kspec - m_numComponents;
plogf(" --- %-12.12s", SpName[kspec].c_str());
plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n",
spStatus[l1], soln[kspec],
spStatus[l1], m_molNumSpecies_old[kspec],
m_molNumSpecies_new[kspec], m_feSpecies_old[kspec], m_feSpecies_curr[kspec],
dgl[l1], dg[l1]);
}
@ -1213,7 +1213,7 @@ namespace VCSnonideal {
}
plogf(" "); vcs_print_line("-", 103);
plogf(" --- Total Old Dimensionless Gibbs Free Energy = %20.13E\n",
vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_old),
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_old),
VCS_DATA_PTR(TPhMoles)));
plogf(" --- Total New Dimensionless Gibbs Free Energy = %20.13E",
vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_feSpecies_curr),
@ -1242,7 +1242,7 @@ namespace VCSnonideal {
* loop.
*/
vcs_dcopy(VCS_DATA_PTR(TPhMoles), VCS_DATA_PTR(TPhMoles1), NPhase);
vcs_dcopy(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_molNumSpecies_new), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_molNumSpecies_new), m_numSpeciesRdc);
vcs_dcopy(VCS_DATA_PTR(dgl), VCS_DATA_PTR(dg), m_numRxnRdc);
vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc);
@ -1277,17 +1277,17 @@ namespace VCSnonideal {
TPhMoles[iph]/TMoles <= VCS_DELETE_PHASE_CUTOFF) {
soldel = 1;
for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
if (PhaseID[kspec] == iph && soln[kspec] > 0.0) {
if (PhaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) {
irxn = kspec - m_numComponents;
if (kspec < m_numComponents) {
if (soln[kspec] > VCS_DELETE_SPECIES_CUTOFF) {
if (m_molNumSpecies_old[kspec] > VCS_DELETE_SPECIES_CUTOFF) {
soldel = 0;
break;
}
} else {
for (k = 0; k < m_numComponents; k++) {
if (sc[irxn][k] != 0.0) {
if (soln[kspec]/soln[k] > VCS_DELETE_PHASE_CUTOFF) {
if (m_molNumSpecies_old[kspec]/m_molNumSpecies_old[k] > VCS_DELETE_PHASE_CUTOFF) {
soldel = 0;
break;
}
@ -1321,7 +1321,7 @@ namespace VCSnonideal {
VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test,
&usedZeroedSpecies);
if (retn != VCS_SUCCESS) return retn;
vcs_dfe(VCS_DATA_PTR(soln), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, true);
uptodate_minors = TRUE;
if (conv) {
@ -1354,7 +1354,7 @@ namespace VCSnonideal {
}
#endif
vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx));
vcs_dfe(VCS_DATA_PTR(soln), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, true);
uptodate_minors = TRUE;
}
@ -1382,7 +1382,7 @@ namespace VCSnonideal {
*/
dofast = (m_numComponents != 1);
for (i = 1; i < m_numComponents; ++i) {
if ((soln[i - 1] * m_spSize[i-1]) < (soln[i] * m_spSize[i])) {
if ((m_molNumSpecies_old[i - 1] * m_spSize[i-1]) < (m_molNumSpecies_old[i] * m_spSize[i])) {
dofast = FALSE;
break;
}
@ -1394,20 +1394,20 @@ namespace VCSnonideal {
for (j = m_numComponents - 1; j >= 0; j--) {
bool doSwap = false;
if (SSPhase[j]) {
doSwap = (soln[l] * m_spSize[l]) > (soln[j] * m_spSize[j] * 1.01);
doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) > (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
if (!SSPhase[i]) {
if (doSwap) {
doSwap = (soln[l]) > (soln[j] * 1.01);
doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
}
}
} else {
if (SSPhase[i]) {
doSwap = (soln[l] * m_spSize[l]) > (soln[j] * m_spSize[j] * 1.01);
doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) > (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
if (!doSwap) {
doSwap = (soln[l]) > (soln[j] * 1.01);
doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
}
} else {
doSwap = (soln[l] * m_spSize[l]) > (soln[j] * m_spSize[j] * 1.01);
doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) > (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
}
}
if (doSwap) {
@ -1428,7 +1428,7 @@ namespace VCSnonideal {
}
#ifdef DEBUG_NOT
if (spStatus[i] == VCS_SPECIES_ZEROEDMS) {
if (soln[j] == 0.0) {
if (m_molNumSpecies_old[j] == 0.0) {
if (sc[i][j] != 0.0) {
if (dg[i] < 0.0) {
#ifdef DEBUG_MODE
@ -1454,20 +1454,20 @@ namespace VCSnonideal {
for (j = 0; j < m_numComponents; ++j) {
bool doSwap = false;
if (SSPhase[j]) {
doSwap = (soln[l] * m_spSize[l]) > (soln[j] * m_spSize[j] * 1.01);
doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) > (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
if (!SSPhase[l]) {
if (doSwap) {
doSwap = (soln[l]) > (soln[j] * 1.01);
doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
}
}
} else {
if (SSPhase[l]) {
doSwap = (soln[l] * m_spSize[l]) > (soln[j] * m_spSize[j] * 1.01);
doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) > (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
if (!doSwap) {
doSwap = (soln[l]) > (soln[j] * 1.01);
doSwap = (m_molNumSpecies_old[l]) > (m_molNumSpecies_old[j] * 1.01);
}
} else {
doSwap = (soln[l] * m_spSize[l]) > (soln[j] * m_spSize[j] * 1.01);
doSwap = (m_molNumSpecies_old[l] * m_spSize[l]) > (m_molNumSpecies_old[j] * m_spSize[j] * 1.01);
}
}
if (doSwap) {
@ -1488,7 +1488,7 @@ namespace VCSnonideal {
}
#ifdef DEBUG_NOT
if (spStatus[i] == VCS_SPECIES_ZEROEDMS) {
if (soln[j] == 0.0) {
if (m_molNumSpecies_old[j] == 0.0) {
if (sc[i][j] != 0.0) {
if (dg[i] < 0.0) {
#ifdef DEBUG_MODE
@ -1579,7 +1579,7 @@ namespace VCSnonideal {
* For this special case, we must reevaluate thermo functions
*/
if (iti != 0) {
vcs_dfe(VCS_DATA_PTR(soln), 0, 0, kspec, kspec+1);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, kspec, kspec+1);
vcs_deltag(0, false);
}
}
@ -1660,7 +1660,7 @@ namespace VCSnonideal {
* for minor species, if needed.
*/
if (iti != 0) {
vcs_dfe(VCS_DATA_PTR(soln), 0, 1, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 1, 0, m_numSpeciesRdc);
vcs_deltag(1, false);
uptodate_minors = TRUE;
}
@ -1770,7 +1770,7 @@ namespace VCSnonideal {
/*
* Go back to evaluate the total moles of gas and liquid.
*/
vcs_dfe(VCS_DATA_PTR(soln), 0, 0, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
/*
*
@ -1860,7 +1860,7 @@ namespace VCSnonideal {
* for minor species and go back to do a full iteration
*/
MajorSpeciesHaveConverged = true;
vcs_dfe(VCS_DATA_PTR(soln), 0, 1, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 1, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
iti = 0;
goto L_MAINLOOP_ALL_SPECIES;
@ -1879,7 +1879,7 @@ namespace VCSnonideal {
* for minor species and go back to do a full iteration
*/
MajorSpeciesHaveConverged = true;
vcs_dfe(VCS_DATA_PTR(soln), 0, 1, 0, m_numSpeciesRdc);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 1, 0, m_numSpeciesRdc);
vcs_deltag(0, false);
iti = 0;
goto L_MAINLOOP_ALL_SPECIES;
@ -1921,7 +1921,7 @@ namespace VCSnonideal {
} else {
iph = PhaseID[kspec];
if (TPhMoles[iph] != 0.0) {
m_molNumSpecies_new[kspec] = soln[kspec] / TPhMoles[iph];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] / TPhMoles[iph];
} else {
/*
* For MultiSpecies phases that are zeroed out,
@ -2010,7 +2010,7 @@ namespace VCSnonideal {
*************************************************************************/
{
double dx;
double w_kspec = soln[kspec];
double w_kspec = m_molNumSpecies_old[kspec];
double *wt_kspec = VCS_DATA_PTR(m_molNumSpecies_new) + kspec;
double wTrial;
double *ds_kspec = VCS_DATA_PTR(ds) + kspec;
@ -2140,18 +2140,18 @@ namespace VCSnonideal {
double dx = delta;
double *sc_irxn = sc[irxn];
for (j = 0; j < m_numComponents; ++j) {
if (soln[j] > 0.0) {
if (m_molNumSpecies_old[j] > 0.0) {
tmp = sc_irxn[j] * dx;
if (-tmp > soln[j]) {
if (-tmp > m_molNumSpecies_old[j]) {
retn = 0;
dx = MIN(dx, - soln[j] / sc_irxn[j]);
dx = MIN(dx, - m_molNumSpecies_old[j] / sc_irxn[j]);
}
}
/*
* If the component has a zero concentration and is a reactant
* in the formation reaction, then dx == 0.0, and we just return.
*/
if (soln[j] <= 0.0) {
if (m_molNumSpecies_old[j] <= 0.0) {
if (sc_irxn[j] < 0.0) {
*delta_ptr = 0.0;
return 0;
@ -2162,16 +2162,16 @@ namespace VCSnonideal {
* ok, we found a positive dx. implement it.
*/
*delta_ptr = dx;
soln[kspec] += dx;
m_molNumSpecies_old[kspec] += dx;
int iph = PhaseID[kspec];
TPhMoles[iph] += dx;
for (j = 0; j < m_numComponents; ++j) {
iph = PhaseID[j];
tmp = sc_irxn[j] * dx;
soln[j] += tmp;
m_molNumSpecies_old[j] += tmp;
TPhMoles[iph] += tmp;
if (soln[j] < 0.0) {
soln[j] = 0.0;
if (m_molNumSpecies_old[j] < 0.0) {
m_molNumSpecies_old[j] = 0.0;
}
}
}
@ -2203,13 +2203,13 @@ namespace VCSnonideal {
* Calculate a delta that will eliminate the species.
*/
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
double dx = -(soln[kspec]);
double dx = -(m_molNumSpecies_old[kspec]);
if (dx != 0.0) {
retn = delta_species(kspec, &dx);
if (!retn) {
plogf("zero_species: Couldn't zero the species %d, "
"did delta of %g. orig conc of %g\n",
kspec, dx, soln[kspec] + dx);
kspec, dx, m_molNumSpecies_old[kspec] + dx);
}
}
}
@ -2273,7 +2273,7 @@ namespace VCSnonideal {
/*
* Adjust the total moles in a phase downwards.
*/
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(soln), VCS_DATA_PTR(TPhMoles));
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles));
/*
* Adjust the current number of active species and reactions counters
@ -2291,7 +2291,7 @@ namespace VCSnonideal {
for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (PhaseID[kspec] == iph) {
if (soln[kspec] > 0.0) {
if (m_molNumSpecies_old[kspec] > 0.0) {
Vphase->Existence = 1;
break;
}
@ -2338,7 +2338,7 @@ namespace VCSnonideal {
#endif
/*
* Set the species back to minor species status
* this adjusts soln[] and TPhMoles[]
* this adjusts m_molNumSpecies_old[] and TPhMoles[]
* HKM -> make this a relative mole number!
*/
dx = VCS_DELETE_SPECIES_CUTOFF * 10.;
@ -2351,7 +2351,7 @@ namespace VCSnonideal {
}
int iph = PhaseID[kspec];
vcs_VolPhase *Vphase = VPhaseList[iph];
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(soln), VCS_DATA_PTR(TPhMoles));
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles));
/*
* We may have popped a multispecies phase back
* into existence. If we did, we have to check
@ -2433,11 +2433,11 @@ namespace VCSnonideal {
/*
* calculate an extent of rxn, dx, that zeroes out the species.
*/
dx = - (soln[kspec]);
dx = - (m_molNumSpecies_old[kspec]);
/*
* Set the mole numbers of that species to zero.
*/
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] = 0.0;
m_molNumSpecies_new[kspec] = 0.0;
ds[kspec] = 0.0;
/*
@ -2452,9 +2452,9 @@ namespace VCSnonideal {
* HKM -> note, this will cause a loss of moles!
*/
for (j = 0; j < m_numComponents; ++j) {
soln[j] += sc[irxn][j] * dx;
if (soln[j] < 0.0) {
soln[j] = 0.0;
m_molNumSpecies_old[j] += sc[irxn][j] * dx;
if (m_molNumSpecies_old[j] < 0.0) {
m_molNumSpecies_old[j] = 0.0;
}
}
}
@ -2472,7 +2472,7 @@ namespace VCSnonideal {
for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) {
if (PhaseID[kspec] == iph) {
irxn = kspec - m_numComponents;
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] = 0.0;
m_molNumSpecies_new[kspec] = 0.0;
ds[kspec] = 0.0;
spStatus[irxn] = VCS_SPECIES_ZEROEDPHASE;
@ -2497,7 +2497,7 @@ namespace VCSnonideal {
/*
* Upload the state to the VP object
*/
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(soln), VCS_DATA_PTR(TPhMoles), iph);
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(TPhMoles), iph);
} /* delete_multiphase() *****************************************************/
@ -2574,7 +2574,7 @@ namespace VCSnonideal {
vcs_reinsert_deleted(kspec);
npb++;
} else {
soln[kspec] = 0.0;
m_molNumSpecies_old[kspec] = 0.0;
}
} else if (TPhMoles[iph] > 0.0) {
if (dg[irxn] < xtcutoff[iph]) {
@ -2632,7 +2632,7 @@ namespace VCSnonideal {
}
}
vcs_dfe(VCS_DATA_PTR(soln), 0, 0, 0, m_numSpeciesTot);
vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), 0, 0, 0, m_numSpeciesTot);
vcs_deltag(0, true);
}
@ -2748,7 +2748,7 @@ namespace VCSnonideal {
dptr = VCS_DATA_PTR(m_molNumSpecies_new);
for (kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
m_molNumSpecies_new[kspec] = soln[kspec] + al * ds[kspec];
m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + al * ds[kspec];
}
for (iph = 0; iph < NPhase; iph++) {
TPhMoles1[iph] = TPhMoles[iph] + al * DelTPhMoles[iph];
@ -2834,7 +2834,7 @@ namespace VCSnonideal {
* top of the loop, when necessary
*/
if (UseActCoeffJac) {
vcs_CalcLnActCoeffJac(VCS_DATA_PTR(soln));
vcs_CalcLnActCoeffJac(VCS_DATA_PTR(m_molNumSpecies_old));
}
/************************************************************************
******** LOOP OVER THE FORMATION REACTIONS *****************************
@ -2851,7 +2851,7 @@ namespace VCSnonideal {
dnPhase_irxn = DnPhase[irxn];
if (soln[kspec] == 0.0 && (! SSPhase[kspec])) {
if (m_molNumSpecies_old[kspec] == 0.0 && (! SSPhase[kspec])) {
/********************************************************************/
/******* MULTISPECIES PHASE WITH total moles equal to zero *********/
/*******************************************************************/
@ -2907,7 +2907,7 @@ namespace VCSnonideal {
if (vcs_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", SpName[kspec].c_str());
plogf(" %12.4E %12.4E %12.4E | %s\n",
soln[kspec], ds[kspec], dg[irxn], ANOTE);
m_molNumSpecies_old[kspec], ds[kspec], dg[irxn], ANOTE);
}
#endif
continue;
@ -2923,7 +2923,7 @@ namespace VCSnonideal {
if (vcs_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", SpName[kspec].c_str());
plogf(" %12.4E %12.4E %12.4E | %s\n",
soln[kspec], ds[kspec], dg[irxn], ANOTE);
m_molNumSpecies_old[kspec], ds[kspec], dg[irxn], ANOTE);
}
#endif
continue;
@ -2934,12 +2934,12 @@ namespace VCSnonideal {
if (SSPhase[kspec]) {
s = 0.0;
} else {
s = 1.0 / soln[kspec] ;
s = 1.0 / m_molNumSpecies_old[kspec] ;
}
for (j = 0; j < m_numComponents; ++j) {
if (!SSPhase[j]) {
if (soln[j] > 0.0) {
s += SQUARE(sc[irxn][j]) / soln[j];
if (m_molNumSpecies_old[j] > 0.0) {
s += SQUARE(sc[irxn][j]) / m_molNumSpecies_old[j];
}
}
}
@ -2976,19 +2976,19 @@ namespace VCSnonideal {
double stoicC = sc[irxn][j];
if (stoicC != 0.0) {
double negChangeComp = - stoicC * ds[kspec];
if (negChangeComp > soln[j]) {
if (soln[j] > 0.0) {
if (negChangeComp > m_molNumSpecies_old[j]) {
if (m_molNumSpecies_old[j] > 0.0) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %d (%10s) going neg", ds[kspec],
-soln[j]/stoicC, j, SpName[j].c_str());
-m_molNumSpecies_old[j]/stoicC, j, SpName[j].c_str());
#endif
ds[kspec] = - soln[j] / stoicC;
ds[kspec] = - m_molNumSpecies_old[j] / stoicC;
} else {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to component %d (%10s) zero", ds[kspec],
-soln[j]/stoicC, j, SpName[j].c_str());
-m_molNumSpecies_old[j]/stoicC, j, SpName[j].c_str());
#endif
ds[kspec] = 0.0;
}
@ -2996,13 +2996,13 @@ namespace VCSnonideal {
}
}
// Implement a damping term that limits ds to the size of the mole number
if (-ds[kspec] > soln[kspec]) {
if (-ds[kspec] > m_molNumSpecies_old[kspec]) {
#ifdef DEBUG_MODE
sprintf(ANOTE, "Delta damped from %g "
"to %g due to %s going negative", ds[kspec],
-soln[kspec], SpName[kspec].c_str());
-m_molNumSpecies_old[kspec], SpName[kspec].c_str());
#endif
ds[kspec] = -soln[kspec];
ds[kspec] = -m_molNumSpecies_old[kspec];
}
} else {
@ -3019,11 +3019,11 @@ namespace VCSnonideal {
* -> The species to be zeroed out will be "k".
*/
if (dg[irxn] > 0.0) {
dss = soln[kspec];
dss = m_molNumSpecies_old[kspec];
k = kspec;
for (j = 0; j < m_numComponents; ++j) {
if (sc[irxn][j] > 0.0) {
xx = soln[j] / sc[irxn][j];
xx = m_molNumSpecies_old[j] / sc[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
@ -3035,7 +3035,7 @@ namespace VCSnonideal {
dss = 1.0e10;
for (j = 0; j < m_numComponents; ++j) {
if (sc[irxn][j] < 0.0) {
xx = -soln[j] / sc[irxn][j];
xx = -m_molNumSpecies_old[j] / sc[irxn][j];
if (xx < dss) {
dss = xx;
k = j;
@ -3052,13 +3052,13 @@ namespace VCSnonideal {
* added back into the component species.
*/
if (dss != 0.0) {
soln[kspec] += dss;
m_molNumSpecies_old[kspec] += dss;
TPhMoles[PhaseID[kspec]] += dss;
for (j = 0; j < m_numComponents; ++j) {
soln[j] += dss * sc[irxn][j];
m_molNumSpecies_old[j] += dss * sc[irxn][j];
TPhMoles[PhaseID[j]] += dss * sc[irxn][j];
}
soln[k] = 0.0;
m_molNumSpecies_old[k] = 0.0;
iph = PhaseID[k];
Vphase = VPhaseList[iph];
Vphase->Existence = 0;
@ -3085,7 +3085,7 @@ namespace VCSnonideal {
if (vcs_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", SpName[kspec].c_str());
plogf(" %12.4E %12.4E %12.4E | %s\n",
soln[kspec], ds[kspec], dg[irxn], ANOTE);
m_molNumSpecies_old[kspec], ds[kspec], dg[irxn], ANOTE);
}
#endif
} /* End of loop over SpeciesUnknownType */
@ -3159,7 +3159,7 @@ namespace VCSnonideal {
dtmp_ptr = sc[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
icase = 1;
}
}
@ -3178,7 +3178,7 @@ namespace VCSnonideal {
dtmp_ptr = sc[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
icase = 1;
}
}
@ -3197,7 +3197,7 @@ namespace VCSnonideal {
dtmp_ptr = sc[irxn];
for (kspec = 0; kspec < m_numComponents; ++kspec) {
dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec];
if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
if (m_molNumSpecies_old[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) {
icase = 1;
}
}
@ -3258,7 +3258,7 @@ namespace VCSnonideal {
for (k = 0; k < Vphase->NVolSpecies; k++) {
kspec = Vphase->IndSpecies[k];
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
sum += soln[kspec];
sum += m_molNumSpecies_old[kspec];
}
if (sum > 0.0) break;
}
@ -3396,7 +3396,7 @@ namespace VCSnonideal {
/*
* Use a temporary work array for the mole numbers, aw[]
*/
vcs_dcopy(aw, VCS_DATA_PTR(soln), m_numSpeciesTot);
vcs_dcopy(aw, VCS_DATA_PTR(m_molNumSpecies_old), m_numSpeciesTot);
/*
* Take out the Voltage unknowns from consideration
*/
@ -3609,8 +3609,8 @@ namespace VCSnonideal {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", (SpName[k]).c_str());
plogf("(%9.2g) replaces %-12.12s", soln[k], SpName[jr].c_str());
plogf("(%9.2g) as component %3d\n", soln[jr], jr);
plogf("(%9.2g) replaces %-12.12s", m_molNumSpecies_old[k], SpName[jr].c_str());
plogf("(%9.2g) as component %3d\n", m_molNumSpecies_old[jr], jr);
}
#endif
vcs_switch_pos(FALSE, jr, k);
@ -3620,7 +3620,7 @@ namespace VCSnonideal {
else {
if (vcs_debug_print_lvl >= 2) {
plogf(" --- %-12.12s", SpName[k].c_str());
plogf("(%9.2g) remains ", soln[k]);
plogf("(%9.2g) remains ", m_molNumSpecies_old[k]);
plogf(" as component %3d\n", jr);
}
}
@ -3765,7 +3765,7 @@ namespace VCSnonideal {
}
plogf("\n --- Components Moles:");
for (j = 0; j < ncTrial; j++) {
plogf("%10.3g", soln[j]);
plogf("%10.3g", m_molNumSpecies_old[j]);
}
plogf("\n --- NonComponent| Moles | ");
for (j = 0; j < ncTrial; j++) {
@ -3776,7 +3776,7 @@ namespace VCSnonideal {
for (i = 0; i < m_numRxnTot; i++) {
plogf(" --- %3d ", ir[i]);
plogf("%-10.10s", SpName[ir[i]].c_str());
plogf("|%10.3g|", soln[ir[i]]);
plogf("|%10.3g|", m_molNumSpecies_old[ir[i]]);
for (j = 0; j < ncTrial; j++) {
plogf(" %6.2f", sc[i][j]);
}
@ -3892,7 +3892,7 @@ namespace VCSnonideal {
return VCS_SPECIES_INTERFACIALVOLTAGE;
}
iph = PhaseID[kspec];
if (soln[kspec] <= 0.0) {
if (m_molNumSpecies_old[kspec] <= 0.0) {
if (dg[irxn] >= 0.0) {
/*
* We are here when the species is or should be zeroed out
@ -3920,7 +3920,7 @@ namespace VCSnonideal {
if (stoicC != 0.0) {
double negChangeComp = - stoicC;
if (negChangeComp > 0.0) {
if (soln[j] < 1.0E-60) {
if (m_molNumSpecies_old[j] < 1.0E-60) {
#ifdef DEBUG_MODE
if (vcs_debug_print_lvl >= 2) {
plogf(" --- %s would have popped back into existance but"
@ -3972,7 +3972,7 @@ namespace VCSnonideal {
* Check to see whether the current species is a major component
* of its phase. If it is, it is a major component
*/
if (soln[kspec] > (TPhMoles[iph] * 0.1)) return VCS_SPECIES_MAJOR;
if (m_molNumSpecies_old[kspec] > (TPhMoles[iph] * 0.1)) return VCS_SPECIES_MAJOR;
/*
* Main check in the loop:
* Check to see if there is a component with a mole number that is
@ -3985,7 +3985,7 @@ namespace VCSnonideal {
for (k = 0; k < m_numComponents; ++k) {
if (!(SSPhase[k])) {
if (sc[irxn][k] != 0.0) {
if (soln[kspec] * szAdj >= soln[k] * 0.01) {
if (m_molNumSpecies_old[kspec] * szAdj >= m_molNumSpecies_old[k] * 0.01) {
return VCS_SPECIES_MAJOR;
}
}
@ -4462,7 +4462,7 @@ namespace VCSnonideal {
for (j = 0; j < m_numElemConstraints; ++j) {
for (i = 0; i < m_numSpeciesTot; ++i) {
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
eav[j] += FormulaMatrix[j][i] * soln[i];
eav[j] += FormulaMatrix[j][i] * m_molNumSpecies_old[i];
}
}
}
@ -4538,7 +4538,7 @@ namespace VCSnonideal {
}
for (i = 0; i < m_numSpeciesTot; i++) {
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
TPhMoles[PhaseID[i]] += soln[i];
TPhMoles[PhaseID[i]] += m_molNumSpecies_old[i];
}
}
sum = 0.0;
@ -4571,7 +4571,7 @@ namespace VCSnonideal {
* This routine uploads the state of the system into all of the
* VolumePhase objects in the current problem.
* place
* 0 -> from soln
* 0 -> from m_molNumSpecies_old
* 1 -> from wt
*************************************************************************/
{
@ -4579,7 +4579,7 @@ namespace VCSnonideal {
for (int i = 0; i < NPhase; i++) {
Vphase = VPhaseList[i];
if (place == 0) {
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(soln),
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old),
VCS_DATA_PTR(TPhMoles), i);
} else if (place == 1) {
Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_new),
@ -4669,7 +4669,7 @@ namespace VCSnonideal {
pv2->IndSpecies[kp2] = k1;
vcsUtil_stsw(SpName, k1, k2);
SWAP(soln[k1], soln[k2], t1);
SWAP(m_molNumSpecies_old[k1], m_molNumSpecies_old[k2], t1);
SWAP(SpeciesUnknownType[k1], SpeciesUnknownType[k2], j);
SWAP(m_molNumSpecies_new[k1], m_molNumSpecies_new[k2], t1);
SWAP(m_SSfeSpecies[k1], m_SSfeSpecies[k2], t1);
@ -4801,7 +4801,7 @@ namespace VCSnonideal {
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
iph = PhaseID[kspec];
if (iph == iphase ) {
if (soln[kspec] > 0.0) zeroedPhase = FALSE;
if (m_molNumSpecies_old[kspec] > 0.0) zeroedPhase = FALSE;
dg[irxn] = m_feSpecies_curr[kspec];
dtmp_ptr = sc[irxn];
for (kcomp = 0; kcomp < m_numComponents; ++kcomp) {
@ -4896,7 +4896,7 @@ namespace VCSnonideal {
double w_kspec = VCS_DELETE_SPECIES_CUTOFF;
// Check to make sure that species is zero in the solution vector
// If it isn't, we don't know what's happening
if (soln[kspec] != 0.0) {
if (m_molNumSpecies_old[kspec] != 0.0) {
w_kspec = 0.0;
plogf("we shouldn't be here\n");
exit(-1);
@ -4936,13 +4936,13 @@ namespace VCSnonideal {
for (int j = 0; j < m_numComponents; ++j) {
// Only loop over element contraints that involve positive def. constraints
if (SpeciesUnknownType[j] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (soln[j] > 0.0) {
if (m_molNumSpecies_old[j] > 0.0) {
double tmp = sc_irxn[j] * dx;
if (3.0*(-tmp) > soln[j]) {
dx = MIN(dx, - 0.3333* soln[j] / sc_irxn[j]);
if (3.0*(-tmp) > m_molNumSpecies_old[j]) {
dx = MIN(dx, - 0.3333* m_molNumSpecies_old[j] / sc_irxn[j]);
}
}
if (soln[j] <= 0.0) {
if (m_molNumSpecies_old[j] <= 0.0) {
if (sc_irxn[j] < 0.0) {
dx = 0.0;
}

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

@ -446,7 +446,7 @@ double VCS_SOLVE::vcs_Gxs_calc(int iphase)
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
ts_ptr = SpeciesThermo[kspec];
ac = ts_ptr->eval_ac(kspec);
Gxs += soln[kspec]/totmol * log(ac);
Gxs += m_molNumSpecies_old[kspec]/totmol * log(ac);
} else {
plogf("FILL IN\n");
exit(-1);