From 739d4e483011c940840134a0e84ae2796a5cf473 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 21 Aug 2017 17:59:29 -0400 Subject: [PATCH] [Equil] Eliminate redundant variables from VCS_SOLVE Most of these were originally members of VCS_PROB --- include/cantera/equil/vcs_solve.h | 23 +--- src/equil/vcs_MultiPhaseEquil.cpp | 57 +++------ src/equil/vcs_prob.cpp | 2 +- src/equil/vcs_solve.cpp | 118 +++--------------- .../LatticeSolid_LiSi/output_blessed.txt | 6 +- .../LatticeSolid_LiSi/verbose_blessed.txt | 6 +- 6 files changed, 39 insertions(+), 173 deletions(-) diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index c838f72f5..d18e9d81b 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -763,7 +763,7 @@ public: const double* const fe); //! Transfer the results of the equilibrium calculation back from VCS_SOLVE - int vcs_prob_update(); + void vcs_prob_update(); //! Fully specify the problem to be solved void vcs_prob_specifyFully(); @@ -1006,26 +1006,6 @@ private: vector_fp m_wx; public: - //! @{ Variables moved from VCS_PROB - - //! Vector of chemical potentials of the species. This is a calculated - //! output quantity. length = number of species. - vector_fp m_gibbsSpecies; - - //! Total number of moles of the kth species. - /*! - * This is both an input and an output variable. On input, this is an - * estimate of the mole numbers. The actual element abundance vector - * contains the problem specification. - * - * On output, this contains the solution for the total number of moles of - * the kth species. This vector contains the species in their original order. - */ - vector_fp w; - //! Mole fraction vector. This is a calculated vector, calculated from w[]. - //! length number of species. - vector_fp mf; - //! Print level for print routines int m_printLvl; @@ -1070,7 +1050,6 @@ public: * */ size_t addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT); - //! @} //! This routine resizes the number of elements in the VCS_SOLVE object by //! adding a new element to the end of the element list diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index 23ad912ea..2334f4dde 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -443,24 +443,10 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, } int iSuccess = m_vsolve.vcs(ipr, ip1, maxit); - // Transfer the information back to the MultiPhase object. Note we don't - // just call setMoles, because some multispecies solution phases may be - // zeroed out, and that would cause a problem for that routine. Also, the - // mole fractions of such zeroed out phases actually contain information - // about likely reemergent states. - m_mix->uploadMoleFractionsFromPhases(); - size_t kGlob = 0; - for (size_t ip = 0; ip < m_vsolve.m_numPhases; ip++) { - double phaseMole = 0.0; - ThermoPhase& tref = m_mix->phase(ip); - for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) { - phaseMole += m_vsolve.w[kGlob]; - } - m_mix->setPhaseMoles(ip, phaseMole); - } - double te = tickTock.secondsWC(); if (printLvl > 0) { + vector_fp mu(m_mix->nSpecies()); + m_mix->getChemPotentials(mu.data()); plogf("\n Results from vcs:\n"); if (iSuccess != 0) { plogf("\nVCS FAILED TO CONVERGE!\n"); @@ -478,20 +464,20 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, for (size_t i = 0; i < m_mix->nSpecies(); i++) { plogf("%-12s", m_mix->speciesName(i)); if (m_vsolve.m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" %15.3e %15.3e ", 0.0, m_vsolve.mf[i]); - plogf("%15.3e\n", m_vsolve.m_gibbsSpecies[i]); + plogf(" %15.3e %15.3e ", 0.0, m_mix->moleFraction(i)); + plogf("%15.3e\n", mu[i]); } else { - plogf(" %15.3e %15.3e ", m_vsolve.w[i], m_vsolve.mf[i]); - if (m_vsolve.w[i] <= 0.0) { + plogf(" %15.3e %15.3e ", m_mix->speciesMoles(i), m_mix->moleFraction(i)); + if (m_mix->speciesMoles(i) <= 0.0) { size_t iph = m_vsolve.m_phaseID[i]; vcs_VolPhase* VPhase = m_vsolve.m_VolPhaseList[iph].get(); if (VPhase->nSpecies() > 1) { plogf(" -1.000e+300\n"); } else { - plogf("%15.3e\n", m_vsolve.m_gibbsSpecies[i]); + plogf("%15.3e\n", mu[i]); } } else { - plogf("%15.3e\n", m_vsolve.m_gibbsSpecies[i]); + plogf("%15.3e\n", mu[i]); } } } @@ -513,8 +499,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) throw CanteraError("vcs_MultiPhaseEquil::reportCSV", "Failure to open file"); } - vector_fp& mf = m_vsolve.mf; - double* fe = &m_vsolve.m_gibbsSpecies[0]; vector_fp VolPM; vector_fp activity; vector_fp ac; @@ -524,7 +508,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) double vol = 0.0; for (size_t iphase = 0; iphase < nphase; iphase++) { - size_t istart = m_mix->speciesIndex(0, iphase); ThermoPhase& tref = m_mix->phase(iphase); size_t nSpecies = tref.nSpecies(); VolPM.resize(nSpecies, 0.0); @@ -534,7 +517,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) double TMolesPhase = volP->totalMoles(); double VolPhaseVolumes = 0.0; for (size_t k = 0; k < nSpecies; k++) { - VolPhaseVolumes += VolPM[k] * mf[istart + k]; + VolPhaseVolumes += VolPM[k] * tref.moleFraction(k); } VolPhaseVolumes *= TMolesPhase; vol += VolPhaseVolumes; @@ -549,7 +532,6 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) fprintf(FP,"Number VCS iterations = %d\n", m_vsolve.m_VCount->Its); for (size_t iphase = 0; iphase < nphase; iphase++) { - size_t istart = m_mix->speciesIndex(0, iphase); ThermoPhase& tref = m_mix->phase(iphase); string phaseName = tref.name(); vcs_VolPhase* volP = m_vsolve.m_VolPhaseList[iphase].get(); @@ -569,7 +551,7 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) tref.getChemPotentials(&mu[0]); double VolPhaseVolumes = 0.0; for (size_t k = 0; k < nSpecies; k++) { - VolPhaseVolumes += VolPM[k] * mf[istart + k]; + VolPhaseVolumes += VolPM[k] * tref.moleFraction(k); } VolPhaseVolumes *= TMolesPhase; vol += VolPhaseVolumes; @@ -594,9 +576,9 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) "%11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n", sName.c_str(), phaseName.c_str(), TMolesPhase, - mf[istart + k], molalities[k], ac[k], activity[k], + tref.moleFraction(k), molalities[k], ac[k], activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6, - mf[istart + k] * TMolesPhase, + tref.moleFraction(k) * TMolesPhase, VolPM[k], VolPhaseVolumes); } } else { @@ -618,23 +600,12 @@ void vcs_MultiPhaseEquil::reportCSV(const std::string& reportFile) "%11.3e, %11.3e,% 11.3e, %11.3e, %11.3e\n", sName.c_str(), phaseName.c_str(), TMolesPhase, - mf[istart + k], molalities[k], ac[k], + tref.moleFraction(k), molalities[k], ac[k], activity[k], mu0[k]*1.0E-6, mu[k]*1.0E-6, - mf[istart + k] * TMolesPhase, + tref.moleFraction(k) * TMolesPhase, VolPM[k], VolPhaseVolumes); } } - - // Check consistency: These should be equal - tref.getChemPotentials(fe+istart); - for (size_t k = 0; k < nSpecies; k++) { - if (!vcs_doubleEqual(fe[istart+k], mu[k])) { - fprintf(FP,"ERROR: incompatibility!\n"); - fclose(FP); - throw CanteraError("vcs_MultiPhaseEquil::reportCSV", "incompatibility!"); - } - } - } fclose(FP); } diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index f958c7da9..eba21578c 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -48,7 +48,7 @@ void VCS_SOLVE::prob_report(int print_lvl) plogf("%16s %5d %16s", m_mix->speciesName(i), m_phaseID[i], Vphase->PhaseName); if (m_doEstimateEquil >= 0) { - plogf(" %-10.5g", w[i]); + plogf(" %-10.5g", m_molNumSpecies_old[i]); } else { plogf(" N/A"); } diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index ab328cf45..8171656d7 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -49,9 +49,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : m_debug_print_lvl(0), m_timing_print_lvl(1) { - m_gibbsSpecies.resize(m_nsp, 0.0); - w.resize(m_nsp, 0.0); - mf.resize(m_nsp, 0.0); m_speciesThermoList.resize(m_nsp); for (size_t kspec = 0; kspec < m_nsp; kspec++) { m_speciesThermoList[kspec].reset(new VCS_SPECIES_THERMO()); @@ -215,9 +212,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : // object into the vprob object. addPhaseElements(VolPhase); VolPhase->setState_TP(m_temperature, m_pressurePA); - vector_fp muPhase(tPhase->nSpecies(),0.0); - tPhase->getChemPotentials(&muPhase[0]); - double tMoles = 0.0; // Loop through each species in the current phase for (size_t k = 0; k < nSpPhase; k++) { @@ -235,21 +229,14 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : m_speciesUnknownType[kT] = VolPhase->speciesUnknownType(k); if (m_speciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) { // Set the initial number of kmoles of the species - // and the mole fraction vector - w[kT] = mphase->speciesMoles(kT); - tMoles += w[kT]; - mf[kT] = mphase->moleFraction(kT); + m_molNumSpecies_old[kT] = mphase->speciesMoles(kT); } else if (m_speciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - w[kT] = tPhase->electricPotential(); - mf[kT] = mphase->moleFraction(kT); + m_molNumSpecies_old[kT] = tPhase->electricPotential(); } else { throw CanteraError(" vcs_Cantera_to_vsolve() ERROR", "Unknown species type: {}", m_speciesUnknownType[kT]); } - // transfer chemical potential vector - m_gibbsSpecies[kT] = muPhase[k]; - // Transfer the species information from the // volPhase structure to the VPROB structure // This includes: @@ -325,24 +312,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : kT++; } - // Now go back through the species in the phase and assign a valid mole - // fraction to all phases, even if the initial estimate of the total - // number of moles is zero. - if (tMoles > 0.0) { - for (size_t k = 0; k < nSpPhase; k++) { - size_t kTa = VolPhase->spGlobalIndexVCS(k); - mf[kTa] = w[kTa] / tMoles; - } - } else { - // Perhaps, we could do a more sophisticated treatment below. - // But, will start with this. - for (size_t k = 0; k < nSpPhase; k++) { - size_t kTa = VolPhase->spGlobalIndexVCS(k); - mf[kTa]= 1.0 / (double) nSpPhase; - } - } - - VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, &w[0]); + VolPhase->setMolesFromVCS(VCS_STATECALC_OLD, &m_molNumSpecies_old[0]); // Now, calculate a sample naught Gibbs free energy calculation // at the specified temperature. @@ -358,7 +328,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : for (size_t j = 0; j < m_nelem; j++) { for (size_t kspec = 0; kspec < m_nsp; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - m_elemAbundancesGoal[j] += m_formulaMatrix(kspec,j) * w[kspec]; + m_elemAbundancesGoal[j] += m_formulaMatrix(kspec,j) * m_molNumSpecies_old[kspec]; } } if (m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[j] < 1.0E-10) { @@ -383,9 +353,9 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase, VolPhase->PhaseName.c_str()); if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" Volts = %-10.5g\n", w[i]); + plogf(" Volts = %-10.5g\n", m_molNumSpecies_old[i]); } else { - plogf(" %-10.5g\n", w[i]); + plogf(" %-10.5g\n", m_molNumSpecies_old[i]); } } @@ -412,9 +382,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : plogf("\n"); } - // Copy the equilibrium mole number estimate - m_molNumSpecies_old = w; - // TPhInertMoles[] -> must be copied over here for (size_t iph = 0; iph < m_numPhases; iph++) { vcs_VolPhase* Vphase = m_VolPhaseList[iph].get(); @@ -573,7 +540,6 @@ int VCS_SOLVE::vcs(int ipr, int ip1, int maxit) void VCS_SOLVE::vcs_prob_specifyFully() { - size_t kT = 0; // Whether we have an estimate or not gets overwritten on // the call to the equilibrium solver. m_temperature = m_mix->temperature(); @@ -590,17 +556,9 @@ void VCS_SOLVE::vcs_prob_specifyFully() vcs_VolPhase* volPhase = m_VolPhaseList[iphase].get(); volPhase->setState_TP(m_temperature, m_pressurePA); - vector_fp muPhase(tPhase->nSpecies(),0.0); - tPhase->getChemPotentials(&muPhase[0]); // Loop through each species in the current phase size_t nSpPhase = tPhase->nSpecies(); - for (size_t k = 0; k < nSpPhase; k++) { - // transfer chemical potential vector - m_gibbsSpecies[invSpecies[kT]] = muPhase[k]; - - kT++; - } if ((nSpPhase == 1) && (volPhase->phiVarIndex() == 0)) { volPhase->setExistence(VCS_PHASE_EXIST_ALWAYS); } else if (volPhase->totalMoles() > 0.0) { @@ -627,9 +585,9 @@ void VCS_SOLVE::vcs_prob_specifyFully() plogf("%16s %5d %16s", m_speciesName[i].c_str(), iphase, VolPhase->PhaseName.c_str()); if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - plogf(" Volts = %-10.5g\n", w[i]); + plogf(" Volts = %-10.5g\n", m_molNumSpecies_old[i]); } else { - plogf(" %-10.5g\n", w[i]); + plogf(" %-10.5g\n", m_molNumSpecies_old[i]); } } @@ -668,59 +626,17 @@ void VCS_SOLVE::vcs_prob_specifyFully() m_numRxnRdc = m_numRxnTot; } -int VCS_SOLVE::vcs_prob_update() +void VCS_SOLVE::vcs_prob_update() { - vcs_tmoles(); - m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA, - &m_molNumSpecies_old[0], &m_PMVolumeSpecies[0]); - - for (size_t i = 0; i < m_nsp; ++i) { - // Find the index of I in the index vector, m_speciesIndexVector[]. Call - // it K1 and continue. - size_t k1 = 0; - for (size_t j = 0; j < m_nsp; ++j) { - k1 = j; - if (m_speciesMapIndex[j] == i) { - break; - } - } - - // Switch the species data back from K1 into I - if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - w[i] = m_molNumSpecies_old[k1]; - } else { - w[i] = 0.0; - } - m_gibbsSpecies[i] = m_feSpecies_old[k1]; + // Transfer the information back to the MultiPhase object. Note we don't + // just call setMoles, because some multispecies solution phases may be + // zeroed out, and that would cause a problem for that routine. Also, the + // mole fractions of such zeroed out phases actually contain information + // about likely reemergent states. + m_mix->uploadMoleFractionsFromPhases(); + for (size_t ip = 0; ip < m_numPhases; ip++) { + m_mix->setPhaseMoles(ip, m_VolPhaseList[ip]->totalMoles()); } - - for (size_t iph = 0; iph < m_numPhases; iph++) { - vcs_VolPhase* vPhase = m_VolPhaseList[iph].get(); - double sumMoles = vPhase->totalMolesInert(); - const vector_fp & mfVector = vPhase->moleFractions(); - for (size_t k = 0; k < vPhase->nSpecies(); k++) { - size_t kT = vPhase->spGlobalIndexVCS(k); - size_t kOrig = m_speciesMapIndex[kT]; - mf[kOrig] = mfVector[k]; - if (vPhase->phiVarIndex() == k) { - double tmp = m_molNumSpecies_old[vPhase->spGlobalIndexVCS(k)]; - if (!vcs_doubleEqual(vPhase->electricPotential(), tmp)) { - throw CanteraError("VCS_SOLVE::vcs_prob_update", - "We have an inconsistency in voltage, {} {}", - vPhase->electricPotential(), tmp); - } - } - if (vPhase->speciesUnknownType(k) != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - sumMoles += w[kOrig]; - } - } - if (! vcs_doubleEqual(sumMoles, vPhase->totalMoles())) { - throw CanteraError("VCS_SOLVE::vcs_prob_update", - "We have an inconsistency in total moles, {} {}", - sumMoles, vPhase->totalMoles()); - } - } - return VCS_SUCCESS; } void VCS_SOLVE::vcs_counters_init(int ifunc) diff --git a/test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt b/test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt index be8da8781..66b9cd125 100644 --- a/test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt +++ b/test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt @@ -231,9 +231,9 @@ Pressure = 101325 Pa ------------------------------------------------------------- LiCl(L) 7.000e+00 7.000e-01 -4.501e+08 KCl(L) 3.000e+00 3.000e-01 -4.992e+08 -Li7Si3(S) 5.000e-01 5.000e-01 -4.237e+08 -Li(i) 1.398e-13 1.398e-13 -1.691e+08 -V(i) 5.000e-01 5.000e-01 -3.513e+06 +Li7Si3(S) 5.000e-01 5.000e-01 -4.201e+08 +Li(i) 1.398e-13 1.398e-13 -1.655e+08 +V(i) 5.000e-01 5.000e-01 8.980e+04 LiFixed 1.000e+02 1.000e+00 -1.656e+08 ------------------------------------------------------------- VCS solver succeeded diff --git a/test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt b/test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt index 3a80177e1..a640bb11b 100644 --- a/test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt +++ b/test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt @@ -470,9 +470,9 @@ Pressure = 101325 Pa ------------------------------------------------------------- LiCl(L) 7.000e+00 7.000e-01 -4.501e+08 KCl(L) 3.000e+00 3.000e-01 -4.992e+08 -Li7Si3(S) 5.000e-01 5.000e-01 -4.237e+08 -Li(i) 1.398e-13 1.398e-13 -1.691e+08 -V(i) 5.000e-01 5.000e-01 -3.513e+06 +Li7Si3(S) 5.000e-01 5.000e-01 -4.201e+08 +Li(i) 1.398e-13 1.398e-13 -1.655e+08 +V(i) 5.000e-01 5.000e-01 8.980e+04 LiFixed 1.000e+02 1.000e+00 -1.656e+08 ------------------------------------------------------------- VCS solver succeeded