diff --git a/Cantera/src/equil/vcs_elem.cpp b/Cantera/src/equil/vcs_elem.cpp index 0af328933..f92c646e4 100644 --- a/Cantera/src/equil/vcs_elem.cpp +++ b/Cantera/src/equil/vcs_elem.cpp @@ -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; diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index 139eec9a2..400788350 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -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(); } diff --git a/Cantera/src/equil/vcs_nondim.cpp b/Cantera/src/equil/vcs_nondim.cpp index 4d3723c52..f47ce211d 100644 --- a/Cantera/src/equil/vcs_nondim.cpp +++ b/Cantera/src/equil/vcs_nondim.cpp @@ -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) { diff --git a/Cantera/src/equil/vcs_prep.cpp b/Cantera/src/equil/vcs_prep.cpp index c39e9c3fe..012b2381e 100644 --- a/Cantera/src/equil/vcs_prep.cpp +++ b/Cantera/src/equil/vcs_prep.cpp @@ -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; diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 519e9cc9a..7fd33fe6b 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -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]; } diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index b9671da84..8703b281f 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -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; diff --git a/Cantera/src/equil/vcs_setMolesLinProg.cpp b/Cantera/src/equil/vcs_setMolesLinProg.cpp index 298d0ae70..1194d275f 100644 --- a/Cantera/src/equil/vcs_setMolesLinProg.cpp +++ b/Cantera/src/equil/vcs_setMolesLinProg.cpp @@ -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 diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 059af1121..42cf9b703 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -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); diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 89728449c..869d8a5aa 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -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 soln; + std::vector m_molNumSpecies_old; //! Specifies the species unknown type /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 0c93f14ef..3107fb71a 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -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; } diff --git a/Cantera/src/equil/vcs_species_thermo.cpp b/Cantera/src/equil/vcs_species_thermo.cpp index 7a30c159e..e142fb1a5 100644 --- a/Cantera/src/equil/vcs_species_thermo.cpp +++ b/Cantera/src/equil/vcs_species_thermo.cpp @@ -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);