diff --git a/Cantera/src/equil/vcs_TP.cpp b/Cantera/src/equil/vcs_TP.cpp index aaecb9989..b7ca363ba 100644 --- a/Cantera/src/equil/vcs_TP.cpp +++ b/Cantera/src/equil/vcs_TP.cpp @@ -18,179 +18,177 @@ namespace VCSnonideal { -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ -int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg) + int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg) - /************************************************************************** - * - * vcs_TP: - * - * Solve an equilibrium problem at a particular fixed temperature - * and pressure - * - * ipr = 1 -> Print results to standard output - * 0 -> don't report on anything - * ip1 = 1 -> Print intermediate results. - * maxit -> Maximum number of iterations for the algorithm - * T = Temperature (Kelvin) - * pres = Pressure (units given by if__ variable) - * - * Return Codes - * ------------------ - * 0 = Equilibrium Achieved - * 1 = Range space error encountered. The element abundance criteria are - * only partially satisfied. Specifically, the first NC= (number of - * components) conditions are satisfied. However, the full NE - * (number of elements) conditions are not satisfied. The equilibrirum - * condition is returned. - * -1 = Maximum number of iterations is exceeded. Convergence was not - * found. - ***************************************************************************/ -{ - int retn, iconv; - /* - * Store the temperature and pressure in the private global variables - */ - T = T_arg; - Pres = pres_arg; - /* - * Evaluate the standard state free energies - * at the current temperatures and pressures. - */ - iconv = vcs_evalSS_TP(ipr, ip1, T, pres_arg); - - /* - * Prepare the problem data: - * ->nondimensionalize the free energies using - * the divisor, R * T - */ - vcs_nondim_TP(); - /* - * Prep the fe field - */ - vcs_fePrep_TP(); - /* - * Decide whether we need an initial estimate of the solution - * If so, go get one. If not, then - */ - if (iest == -1) { - retn = vcs_inest_TP(); - if (retn != VCS_SUCCESS) { - plogf("vcs_inest_TP returned a failure flag\n"); - } - } - /* - * Solve the problem at a fixed Temperature and Pressure - * (all information concerning Temperature and Pressure has already - * been derived. The free energies are now in dimensionless form.) - */ - iconv = vcs_solve_TP(ipr, ip1, maxit); - - /* - * Redimensionalize the free energies using - * the reverse of vcs_nondim to add back units. - */ - vcs_redim_TP(); - /* - * Return the convergence success flag. - */ - return iconv; -} -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ -/*ARGSUSED*/ -int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) - - /************************************************************************** - * - * vcs_evalSS_TP: - * - * IPR = 1 -> Print results to standard output - * 0 -> don't report on anything - * IP1 = 1 -> Print intermediate results. - * T = Temperature (Kelvin) - * Pres = Pressure (units of if__ variable) - * - * Evaluate the standard state free energies at the current temperature - * and pressure. Ideal gas pressure contribution is added in here. - * - ***************************************************************************/ -{ - // int i; - //double R; - /* - * At this level of the program, we are still using values - * for the free energies that have units. - */ - // R = vcsUtil_gasConstant(m_VCS_UnitsFormat); - - /* - * We need to special case VCS_UNITS_UNITLESS, here. - * cpc_ts_GStar_calc() returns units of Kelvin. Also, the temperature - * comes into play in calculating the ideal equation of state - * contributions, and other equations of state also. Therefore, - * we will emulate the VCS_UNITS_KELVIN case, here by chaning - * the initial gibbs free energy units to Kelvin before feeding - * them to the cpc_ts_GStar_calc() routine. Then, we will revert - * them back to unitless at the end of this routine. - */ - - - /* - * Loop over the species calculating the standard state Gibbs free - * energies. -> These are energies that only depend upon the Temperature - * and possibly on the pressure (i.e., ideal gas, etc). - */ - // HKM -> We can change this to looks over phases, calling the vcs_VolPhase - // object. Working to get rid of VCS_SPECIES_THERMO object - //for (i = 0; i < m_numSpeciesTot; ++i) { - // VCS_SPECIES_THERMO *spt = SpeciesThermo[i]; - // ff[i] = R * spt->GStar_R_calc(i, Temp, pres); - //} - - for (int iph = 0; iph < NPhase; iph++) { - vcs_VolPhase* vph = VPhaseList[iph]; - vph->setState_TP(T, Pres); - vph->sendToVCSGStar(VCS_DATA_PTR(m_SSfeSpecies)); - } - - if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { - for (int i = 0; i < m_numSpeciesTot; ++i) { - m_SSfeSpecies[i] /= Temp; - } - } - return VCS_SUCCESS; -} /***************************************************************************/ - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - -void VCS_SOLVE::vcs_fePrep_TP(void) - - /************************************************************************** - * - * - ***************************************************************************/ -{ - int i; - for (i = 0; i < m_numSpeciesTot; ++i) { + /************************************************************************** + * + * vcs_TP: + * + * Solve an equilibrium problem at a particular fixed temperature + * and pressure + * + * ipr = 1 -> Print results to standard output + * 0 -> don't report on anything + * ip1 = 1 -> Print intermediate results. + * maxit -> Maximum number of iterations for the algorithm + * T = Temperature (Kelvin) + * pres = Pressure (units given by if__ variable) + * + * Return Codes + * ------------------ + * 0 = Equilibrium Achieved + * 1 = Range space error encountered. The element abundance criteria are + * only partially satisfied. Specifically, the first NC= (number of + * components) conditions are satisfied. However, the full NE + * (number of elements) conditions are not satisfied. The equilibrirum + * condition is returned. + * -1 = Maximum number of iterations is exceeded. Convergence was not + * found. + ***************************************************************************/ + { + int retn, iconv; /* - * For single species phases, initialize the chemical - * potential with the value of the standard state chemical - * potential. This value doesn't change during the calculation + * Store the temperature and pressure in the private global variables */ - if (SSPhase[i]) { - m_gibbsSpecies[i] = m_SSfeSpecies[i]; + T = T_arg; + Pres = pres_arg; + /* + * Evaluate the standard state free energies + * at the current temperatures and pressures. + */ + iconv = vcs_evalSS_TP(ipr, ip1, T, pres_arg); + + /* + * Prepare the problem data: + * ->nondimensionalize the free energies using + * the divisor, R * T + */ + vcs_nondim_TP(); + /* + * Prep the fe field + */ + vcs_fePrep_TP(); + /* + * Decide whether we need an initial estimate of the solution + * If so, go get one. If not, then + */ + if (iest == -1) { + retn = vcs_inest_TP(); + if (retn != VCS_SUCCESS) { + plogf("vcs_inest_TP returned a failure flag\n"); + } } - } -} /* vcs_fePrep_TP() ********************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ + /* + * Solve the problem at a fixed Temperature and Pressure + * (all information concerning Temperature and Pressure has already + * been derived. The free energies are now in dimensionless form.) + */ + iconv = vcs_solve_TP(ipr, ip1, maxit); + + /* + * Redimensionalize the free energies using + * the reverse of vcs_nondim to add back units. + */ + vcs_redim_TP(); + /* + * Return the convergence success flag. + */ + return iconv; + } + /*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ + /*ARGSUSED*/ + int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) + + /************************************************************************** + * + * vcs_evalSS_TP: + * + * IPR = 1 -> Print results to standard output + * 0 -> don't report on anything + * IP1 = 1 -> Print intermediate results. + * T = Temperature (Kelvin) + * Pres = Pressure (units of if__ variable) + * + * Evaluate the standard state free energies at the current temperature + * and pressure. Ideal gas pressure contribution is added in here. + * + ***************************************************************************/ + { + // int i; + //double R; + /* + * At this level of the program, we are still using values + * for the free energies that have units. + */ + // R = vcsUtil_gasConstant(m_VCS_UnitsFormat); + + /* + * We need to special case VCS_UNITS_UNITLESS, here. + * cpc_ts_GStar_calc() returns units of Kelvin. Also, the temperature + * comes into play in calculating the ideal equation of state + * contributions, and other equations of state also. Therefore, + * we will emulate the VCS_UNITS_KELVIN case, here by chaning + * the initial gibbs free energy units to Kelvin before feeding + * them to the cpc_ts_GStar_calc() routine. Then, we will revert + * them back to unitless at the end of this routine. + */ + + + /* + * Loop over the species calculating the standard state Gibbs free + * energies. -> These are energies that only depend upon the Temperature + * and possibly on the pressure (i.e., ideal gas, etc). + */ + // HKM -> We can change this to looks over phases, calling the vcs_VolPhase + // object. Working to get rid of VCS_SPECIES_THERMO object + //for (i = 0; i < m_numSpeciesTot; ++i) { + // VCS_SPECIES_THERMO *spt = SpeciesThermo[i]; + // ff[i] = R * spt->GStar_R_calc(i, Temp, pres); + //} + + for (int iph = 0; iph < NPhase; iph++) { + vcs_VolPhase* vph = VPhaseList[iph]; + vph->setState_TP(T, Pres); + vph->sendToVCSGStar(VCS_DATA_PTR(m_SSfeSpecies)); + } + + if (m_VCS_UnitsFormat == VCS_UNITS_UNITLESS) { + for (int i = 0; i < m_numSpeciesTot; ++i) { + m_SSfeSpecies[i] /= Temp; + } + } + return VCS_SUCCESS; + } /***************************************************************************/ + + /*****************************************************************************/ + /*****************************************************************************/ + /*****************************************************************************/ + + void VCS_SOLVE::vcs_fePrep_TP(void) + + /************************************************************************** + * + * + ***************************************************************************/ + { + int i; + for (i = 0; i < m_numSpeciesTot; ++i) { + /* + * For single species phases, initialize the chemical + * potential with the value of the standard state chemical + * potential. This value doesn't change during the calculation + */ + if (SSPhase[i]) { + m_feSpecies_curr[i] = m_SSfeSpecies[i]; + } + } + } /* vcs_fePrep_TP() ********************************************************/ + } diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index 5df9477c2..3d6ee1984 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -195,14 +195,14 @@ namespace VCSnonideal { wt[kspec] = 0.0; } } - vcs_dcopy(VCS_DATA_PTR(m_gibbsSpecies), VCS_DATA_PTR(m_SSfeSpecies), + vcs_dcopy(VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(m_SSfeSpecies), nspecies); for (kspec = 0; kspec < m_numComponents; ++kspec) { if (SpeciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { if (! SSPhase[kspec]) { iph = PhaseID[kspec]; - m_gibbsSpecies[kspec] += log(wt[kspec] / TPhMoles[iph]); + m_feSpecies_curr[kspec] += log(wt[kspec] / TPhMoles[iph]); } } else { wt[kspec] = 0.0; @@ -214,11 +214,11 @@ namespace VCSnonideal { for (kspec = 0; kspec < nspecies; ++kspec) { plogf("%s", pprefix); plogf("%-12.12s", SpName[kspec].c_str()); if (kspec < m_numComponents) - plogf("fe* = %15.5g ff = %15.5g\n", m_gibbsSpecies[kspec], + plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_curr[kspec], m_SSfeSpecies[kspec]); else plogf("fe* = %15.5g ff = %15.5g dg* = %15.5g\n", - m_gibbsSpecies[kspec], m_SSfeSpecies[kspec], dg[kspec-m_numComponents]); + m_feSpecies_curr[kspec], m_SSfeSpecies[kspec], dg[kspec-m_numComponents]); } } #endif @@ -324,7 +324,7 @@ namespace VCSnonideal { /* ******************************************* */ vcs_dfe(molNum, 0, 0, 0, nspecies); for (kspec = 0, s = 0.0; kspec < nspecies; ++kspec) { - s += ds[kspec] * m_gibbsSpecies[kspec]; + s += ds[kspec] * m_feSpecies_curr[kspec]; } if (s == 0.0) { finished = TRUE; continue; @@ -502,7 +502,7 @@ namespace VCSnonideal { #ifdef DEBUG_MODE if (vcs_debug_print_lvl >= 2) { plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix, - vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_gibbsSpecies), + vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(TPhMoles))); plogendl(); } diff --git a/Cantera/src/equil/vcs_nondim.cpp b/Cantera/src/equil/vcs_nondim.cpp index 83863942f..4d3723c52 100644 --- a/Cantera/src/equil/vcs_nondim.cpp +++ b/Cantera/src/equil/vcs_nondim.cpp @@ -103,7 +103,7 @@ void VCS_SOLVE::vcs_nondim_TP(void) { * Thus, we may divide it by the temperature. */ m_SSfeSpecies[i] *= tf; - m_gibbsSpecies[i] *= tf; + m_feSpecies_curr[i] *= tf; dg[i] *= tf; dgl[i] *= tf; m_feSpecies_old[i] *= tf; @@ -143,7 +143,7 @@ void VCS_SOLVE::vcs_redim_TP(void) * FF(I), to make it have units, i.e. mu = RT * mu_star */ m_SSfeSpecies[i] *= tf; - m_gibbsSpecies[i] *= tf; + m_feSpecies_curr[i] *= tf; dg[i] *= tf; dgl[i] *= tf; m_feSpecies_old[i] *= tf; diff --git a/Cantera/src/equil/vcs_prep.cpp b/Cantera/src/equil/vcs_prep.cpp index b3f8e47c4..e4b2656a2 100644 --- a/Cantera/src/equil/vcs_prep.cpp +++ b/Cantera/src/equil/vcs_prep.cpp @@ -284,7 +284,7 @@ int VCS_SOLVE::vcs_prep(void) { /* * Initialize various arrays in the data to zero */ - vcs_dzero(VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesTot); + vcs_dzero(VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesTot); vcs_vdzero(m_feSpecies_old, m_numSpeciesTot); vcs_vdzero(wt, m_numSpeciesTot); vcs_dzero(&(DnPhase[0][0]), m_numSpeciesTot*NPhase); diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index e5ced8093..e5324b37b 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -122,7 +122,7 @@ int VCS_SOLVE::vcs_report(int iconv) for (i = 0; i < m_numComponents; ++i) { plogf(" %-12.12s", SpName[i].c_str()); print_space(13); - plogf("%14.7E %14.7E %12.4E", soln[i], wt[i], m_gibbsSpecies[i]); + plogf("%14.7E %14.7E %12.4E", soln[i], wt[i], m_feSpecies_curr[i]); plogf(" %3d", SpeciesUnknownType[i]); plogf("\n"); } @@ -132,10 +132,10 @@ int VCS_SOLVE::vcs_report(int iconv) print_space(13); if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_MOLNUM) { - plogf("%14.7E %14.7E %12.4E", soln[l], wt[l], m_gibbsSpecies[l]); + plogf("%14.7E %14.7E %12.4E", soln[l], wt[l], m_feSpecies_curr[l]); plogf(" MolNum "); } else if (SpeciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" NA %14.7E %12.4E", 1.0, m_gibbsSpecies[l]); + plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_curr[l]); plogf(" Voltage = %14.7E", soln[l]); } else { plogf("we have a problem\n"); @@ -256,7 +256,7 @@ int VCS_SOLVE::vcs_report(int iconv) gaTPhase[j] += gaPhase[j]; } gibbsPhase = vcs_GibbsPhase(iphase, VCS_DATA_PTR(soln), - VCS_DATA_PTR(m_gibbsSpecies)); + VCS_DATA_PTR(m_feSpecies_curr)); gibbsTotal += gibbsPhase; plogf(" | %18.11E |\n", gibbsPhase); } @@ -280,7 +280,7 @@ int VCS_SOLVE::vcs_report(int iconv) * energy of zero */ - g = vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_gibbsSpecies), + g = vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(TPhMoles)); plogf("\n\tTotal Dimensionless Gibbs Free Energy = G/RT = %15.7E\n", g); if (inertYes) @@ -325,17 +325,17 @@ int VCS_SOLVE::vcs_report(int iconv) if (tpmoles > 0.0 && soln[l] > 0.0) { lx = log(soln[l]) - log(tpmoles); } else { - lx = m_gibbsSpecies[l] - m_SSfeSpecies[l] - log(ActCoeff[l]) + SpecLnMnaught[l]; + lx = m_feSpecies_curr[l] - m_SSfeSpecies[l] - log(ActCoeff[l]) + SpecLnMnaught[l]; } } plogf("%14.7E |", lx); plogf("%14.7E | ", eContrib); double tmp = m_SSfeSpecies[l] + log(ActCoeff[l]) + lx - SpecLnMnaught[l] + eContrib; - if (fabs(m_gibbsSpecies[l] - tmp) > 1.0E-8) { + if (fabs(m_feSpecies_curr[l] - tmp) > 1.0E-8) { plogf("\n\t\twe have a problem - doesn't add up\n"); exit(-1); } - plogf(" %12.4E |", m_gibbsSpecies[l]); + plogf(" %12.4E |", m_feSpecies_curr[l]); if (SpecLnMnaught[l] != 0.0) { plogf(" (%14.7E)", - SpecLnMnaught[l]); } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 508bcc923..5fe5907aa 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -109,7 +109,7 @@ namespace VCSnonideal { scSize.resize(nspecies0, 0.0); m_spSize.resize(nspecies0, 1.0); - m_gibbsSpecies.resize(nspecies0, 0.0); + m_feSpecies_curr.resize(nspecies0, 0.0); m_SSfeSpecies.resize(nspecies0, 0.0); m_feSpecies_new.resize(nspecies0, 0.0); soln.resize(nspecies0, 0.0); @@ -789,7 +789,7 @@ namespace VCSnonideal { k = ind[kspec]; soln[kspec] = pub->w[k]; wt[kspec] = pub->mf[k]; - m_gibbsSpecies[kspec] = pub->m_gibbsSpecies[k]; + m_feSpecies_curr[kspec] = pub->m_gibbsSpecies[k]; } /* @@ -925,7 +925,7 @@ namespace VCSnonideal { plogf("voltage species = %g\n", soln[k1]); } pub->mf[i] = wt[k1]; - pub->m_gibbsSpecies[i] = m_gibbsSpecies[k1]; + pub->m_gibbsSpecies[i] = m_feSpecies_curr[k1]; pub->VolPM[i] = VolPM[k1]; } diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 8553639d7..c94b1053b 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -430,8 +430,9 @@ public: * The first NC entries are for components. The following * NR entries are for the current non-component species in the mechanism. * The dimension of this vector is specified by the m_VCS_UnitsFormat variable. + * Length = number of species */ - std::vector m_gibbsSpecies; + std::vector m_feSpecies_curr; //! Free energy vector from the start of the current iteration /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 379dbad5c..a814d96cd 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -442,8 +442,8 @@ namespace VCSnonideal { } } - vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesRdc); - vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc); vcs_dcopy(VCS_DATA_PTR(ActCoeff0), VCS_DATA_PTR(ActCoeff), m_numSpeciesRdc); vcs_dcopy(VCS_DATA_PTR(dgl), VCS_DATA_PTR(dg), m_numRxnRdc); @@ -903,7 +903,7 @@ namespace VCSnonideal { * which have yet to be processed in the main loop */ for (ll = kspec+1; ll < m_numSpeciesRdc; ++ll) { - m_feSpecies_old[ll] = m_gibbsSpecies[ll]; + m_feSpecies_old[ll] = m_feSpecies_curr[ll]; } for (ll = irxn+1; ll < m_numRxnRdc; ++ll) { dgl[ll] = dg[ll]; @@ -1120,7 +1120,7 @@ namespace VCSnonideal { vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(TPhMoles))); plogf(" --- Total tentative Dimensionless Gibbs Free Energy = %20.13E", - vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_gibbsSpecies), + vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(TPhMoles1))); plogendl(); } @@ -1158,7 +1158,7 @@ namespace VCSnonideal { plogf(" Total moles of liquid = %15.7E\n", 0.0); } plogf(" Total New Dimensionless Gibbs Free Energy = %20.13E\n", - vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_gibbsSpecies), + vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(TPhMoles1))); plogf(" -----------------------------------------------------"); plogendl(); @@ -1181,14 +1181,14 @@ namespace VCSnonideal { for (i = 0; i < m_numComponents; ++i) { plogf(" --- %-12.12s", SpName[i].c_str()); plogf(" "); plogf("%14.6E%14.6E%14.6E%14.6E\n", soln[i], - wt[i], m_feSpecies_old[i], m_gibbsSpecies[i]); + wt[i], m_feSpecies_old[i], m_feSpecies_curr[i]); } for (i = m_numComponents; i < m_numSpeciesRdc; ++i) { l1 = i - m_numComponents; plogf(" --- %-12.12s", SpName[i].c_str()); plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", spStatus[l1], soln[i], - wt[i], m_feSpecies_old[i], m_gibbsSpecies[i], + wt[i], m_feSpecies_old[i], m_feSpecies_curr[i], dgl[l1], dg[l1]); } for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { @@ -1196,7 +1196,7 @@ namespace VCSnonideal { plogf(" --- %-12.12s", SpName[kspec].c_str()); plogf(" %2d %14.6E%14.6E%14.6E%14.6E%14.6E%14.6E\n", spStatus[l1], soln[kspec], - wt[kspec], m_feSpecies_old[kspec], m_gibbsSpecies[kspec], + wt[kspec], m_feSpecies_old[kspec], m_feSpecies_curr[kspec], dgl[l1], dg[l1]); } plogf(" ---"); print_space(56); @@ -1216,7 +1216,7 @@ namespace VCSnonideal { vcs_Total_Gibbs(VCS_DATA_PTR(soln), VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(TPhMoles))); plogf(" --- Total New Dimensionless Gibbs Free Energy = %20.13E", - vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_gibbsSpecies), + vcs_Total_Gibbs(VCS_DATA_PTR(wt), VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(TPhMoles1))); plogendl(); if (m_VCount->Its > 550) { @@ -1244,7 +1244,7 @@ namespace VCSnonideal { vcs_dcopy(VCS_DATA_PTR(TPhMoles), VCS_DATA_PTR(TPhMoles1), NPhase); vcs_dcopy(VCS_DATA_PTR(soln), VCS_DATA_PTR(wt), m_numSpeciesRdc); vcs_dcopy(VCS_DATA_PTR(dgl), VCS_DATA_PTR(dg), m_numRxnRdc); - vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_gibbsSpecies), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc); vcs_updateVP(0); /* @@ -2260,7 +2260,7 @@ namespace VCSnonideal { spStatus[irxn] = VCS_SPECIES_DELETED; dg[irxn] = 0.0; dgl[irxn] = 0.0; - m_gibbsSpecies[kspec] = 0.0; + m_feSpecies_curr[kspec] = 0.0; m_feSpecies_old[kspec] = 0.0; wt[kspec] = 0.0; /* @@ -2525,7 +2525,7 @@ namespace VCSnonideal { * for formation reactions */ for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]; + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } /* * Recalculate the DeltaG's of the formation reactions for the @@ -2613,7 +2613,7 @@ namespace VCSnonideal { * ~ infinite dilution. */ for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]; + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } /* * Recalculate the DeltaG's of the formation reactions for the @@ -3155,10 +3155,10 @@ namespace VCSnonideal { for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { if (spStatus[irxn] != VCS_SPECIES_MINOR) { icase = 0; - dg[irxn] = m_gibbsSpecies[ir[irxn]]; + dg[irxn] = m_feSpecies_curr[ir[irxn]]; dtmp_ptr = sc[irxn]; for (kspec = 0; kspec < m_numComponents; ++kspec) { - dg[irxn] += dtmp_ptr[kspec] * m_gibbsSpecies[kspec]; + dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec]; if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) { icase = 1; } @@ -3174,10 +3174,10 @@ namespace VCSnonideal { /* ************************************************* */ for (irxn = 0; irxn < irxnl; ++irxn) { icase = 0; - dg[irxn] = m_gibbsSpecies[ir[irxn]]; + dg[irxn] = m_feSpecies_curr[ir[irxn]]; dtmp_ptr = sc[irxn]; for (kspec = 0; kspec < m_numComponents; ++kspec) { - dg[irxn] += dtmp_ptr[kspec] * m_gibbsSpecies[kspec]; + dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec]; if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) { icase = 1; } @@ -3193,10 +3193,10 @@ namespace VCSnonideal { for (irxn = 0; irxn < m_numRxnRdc; ++irxn) { if (spStatus[irxn] <= VCS_SPECIES_MINOR) { icase = 0; - dg[irxn] = m_gibbsSpecies[ir[irxn]]; + dg[irxn] = m_feSpecies_curr[ir[irxn]]; dtmp_ptr = sc[irxn]; for (kspec = 0; kspec < m_numComponents; ++kspec) { - dg[irxn] += dtmp_ptr[kspec] * m_gibbsSpecies[kspec]; + dg[irxn] += dtmp_ptr[kspec] * m_feSpecies_curr[kspec]; if (soln[kspec] < VCS_DELETE_MINORSPECIES_CUTOFF && dtmp_ptr[kspec] < 0.0) { icase = 1; } @@ -4328,24 +4328,24 @@ namespace VCSnonideal { exit(-1); } #endif - m_gibbsSpecies[kspec] = + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase]; } else { if (SSPhase[kspec]) { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]; + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } else { if (z[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) { iph = PhaseID[kspec]; if (tPhMoles_ptr[iph] > 0.0) { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF) - tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase]; } else { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]; + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } } else { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec]) + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec]) - tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase]; } @@ -4371,24 +4371,24 @@ namespace VCSnonideal { exit(-1); } #endif - m_gibbsSpecies[kspec] = + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase]; } else { if (SSPhase[kspec]) { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]; + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } else { if (z[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) { iph = PhaseID[kspec]; if (tPhMoles_ptr[iph] > 0.0) { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF) - tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase]; ; } else { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]; + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } } else { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec]) + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec]) - tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase]; } @@ -4415,24 +4415,24 @@ namespace VCSnonideal { exit(-1); } #endif - m_gibbsSpecies[kspec] = + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + Charge[kspec] * Faraday_dim * phasePhi[iphase]; ; } else { if (SSPhase[kspec]) { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]; + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } else { if (z[kspec] <= VCS_DELETE_MINORSPECIES_CUTOFF) { iph = PhaseID[kspec]; if (tPhMoles_ptr[iph] > 0.0) { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * VCS_DELETE_MINORSPECIES_CUTOFF) - tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]; } else { - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec]; + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec]; } } else { st_ptr = SpeciesThermo[kspec]; - m_gibbsSpecies[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec]) + m_feSpecies_curr[kspec] = m_SSfeSpecies[kspec] + log(ActCoeff[kspec] * z[kspec]) - tlogMoles[PhaseID[kspec]] - SpecLnMnaught[kspec]; } } @@ -4674,7 +4674,7 @@ namespace VCSnonideal { SWAP(wt[k1], wt[k2], t1); SWAP(m_SSfeSpecies[k1], m_SSfeSpecies[k2], t1); SWAP(m_spSize[k1], m_spSize[k2], t1); - SWAP(m_gibbsSpecies[k1], m_gibbsSpecies[k2], t1); + SWAP(m_feSpecies_curr[k1], m_feSpecies_curr[k2], t1); SWAP(ds[k1], ds[k2], t1); SWAP(m_feSpecies_old[k1], m_feSpecies_old[k2], t1); SWAP(m_feSpecies_new[k1], m_feSpecies_new[k2], t1); @@ -4783,10 +4783,10 @@ namespace VCSnonideal { #endif if (kspec >= m_numComponents) { irxn = kspec - m_numComponents; - dg[irxn] = m_gibbsSpecies[kspec]; + dg[irxn] = m_feSpecies_curr[kspec]; dtmp_ptr = sc[irxn]; for (kcomp = 0; kcomp < m_numComponents; ++kcomp) { - dg[irxn] += dtmp_ptr[kcomp] * m_gibbsSpecies[kcomp]; + dg[irxn] += dtmp_ptr[kcomp] * m_feSpecies_curr[kcomp]; } } } @@ -4802,10 +4802,10 @@ namespace VCSnonideal { iph = PhaseID[kspec]; if (iph == iphase ) { if (soln[kspec] > 0.0) zeroedPhase = FALSE; - dg[irxn] = m_gibbsSpecies[kspec]; + dg[irxn] = m_feSpecies_curr[kspec]; dtmp_ptr = sc[irxn]; for (kcomp = 0; kcomp < m_numComponents; ++kcomp) { - dg[irxn] += dtmp_ptr[kcomp] * m_gibbsSpecies[kcomp]; + dg[irxn] += dtmp_ptr[kcomp] * m_feSpecies_curr[kcomp]; } } }