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

[Equil] Eliminate redundant variables from VCS_SOLVE

Browse source 
Most of these were originally members of VCS_PROB
This commit is contained in:
Ray Speth 2017-08-21 17:59:29 -04:00
parent 7eb939dc5f
commit 739d4e4830
6 changed files with 39 additions and 173 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

23
include/cantera/equil/vcs_solve.h
View file

@ -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

57
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -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);
}

2
src/equil/vcs_prob.cpp
View file

@ -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");
}

118
src/equil/vcs_solve.cpp
View file

@ -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)

6
test_problems/VCSnonideal/LatticeSolid_LiSi/output_blessed.txt
View file

@ -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

6
test_problems/VCSnonideal/LatticeSolid_LiSi/verbose_blessed.txt
View file

@ -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