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

[Equil] Eliminate some redundant variables in VCS_SOLVE

Browse source
This commit is contained in:
Ray Speth 2017-08-18 17:54:25 -04:00
parent bba0d8edf0
commit da3ba35945
4 changed files with 38 additions and 155 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -1024,50 +1024,6 @@ public:
//! Mole fraction vector. This is a calculated vector, calculated from w[].
//! length number of species.
vector_fp mf;
//! Element abundances for jth element
/*!
* This is input from the input file and is considered a constant from
* thereon within the vcs_solve_TP().
*/
vector_fp gai;
//! Formula Matrix for the problem
/*!
* FormulaMatrix(kspec,j) = Number of elements, j, in the kspec species
*/
Array2D FormulaMatrix;
//! Specifies the species unknown type
/*!
* There are two types. One is the straightforward species, with the mole
* number w[k], as the unknown. The second is the an interfacial voltage
* where w[k] refers to the interfacial voltage in volts.
*
* These species types correspond to metallic electrons corresponding to
* electrodes. The voltage and other interfacial conditions sets up an
* interfacial current, which is set to zero in this initial treatment.
* Later we may have non-zero interfacial currents.
*/
vector_int SpeciesUnknownType;
//! Mapping between the species and the phases
std::vector<size_t> PhaseID;
//! Specifies whether an element constraint is active
/*!
* The default is true
* Length = nelements
*/
vector_int ElActive;
//! Molecular weight of species
/*!
* WtSpecies[k] = molecular weight of species in gm/mol
*/
vector_fp WtSpecies;
//! Charge of each species
vector_fp Charge;
//! Array of phase structures
std::vector<vcs_VolPhase*> VPhaseList;
@ -1126,9 +1082,6 @@ public:
size_t addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT);
//! @}
//! Calculate the element abundance vector from the mole numbers
void set_gai();
//! This routine resizes the number of elements in the VCS_SOLVE object by
//! adding a new element to the end of the element list
/*!

4
src/equil/vcs_MultiPhaseEquil.cpp
View file

@ -477,13 +477,13 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil,
"-----------\n");
for (size_t i = 0; i < m_mix->nSpecies(); i++) {
plogf("%-12s", m_mix->speciesName(i));
if (m_vsolve.SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
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]);
} else {
plogf(" %15.3e %15.3e ", m_vsolve.w[i], m_vsolve.mf[i]);
if (m_vsolve.w[i] <= 0.0) {
size_t iph = m_vsolve.PhaseID[i];
size_t iph = m_vsolve.m_phaseID[i];
vcs_VolPhase* VPhase = m_vsolve.VPhaseList[iph];
if (VPhase->nSpecies() > 1) {
plogf(" -1.000e+300\n");

31
src/equil/vcs_prob.cpp
View file

@ -22,18 +22,6 @@ using namespace std;
namespace Cantera
{
void VCS_SOLVE::set_gai()
{
gai.assign(gai.size(), 0.0);
for (size_t j = 0; j < m_nelem; j++) {
for (size_t kspec = 0; kspec < m_nsp; kspec++) {
if (SpeciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
gai[j] += FormulaMatrix(kspec,j) * w[kspec];
}
}
}
}
void VCS_SOLVE::prob_report(int print_lvl)
{
m_printLvl = print_lvl;
@ -56,17 +44,17 @@ void VCS_SOLVE::prob_report(int print_lvl)
plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_Moles Species_Type\n");
for (size_t i = 0; i < m_nsp; i++) {
vcs_VolPhase* Vphase = VPhaseList[PhaseID[i]];
plogf("%16s %5d %16s", m_mix->speciesName(i), PhaseID[i],
vcs_VolPhase* Vphase = VPhaseList[m_phaseID[i]];
plogf("%16s %5d %16s", m_mix->speciesName(i), m_phaseID[i],
Vphase->PhaseName);
if (m_doEstimateEquil >= 0) {
plogf(" %-10.5g", w[i]);
} else {
plogf(" N/A");
}
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" Mol_Num");
} else if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
} else if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Voltage");
} else {
plogf(" ");
@ -99,8 +87,8 @@ void VCS_SOLVE::prob_report(int print_lvl)
for (size_t i = 0; i < m_nelem; ++i) {
writeline(' ', 26, false);
plogf("%-2.2s", m_elementName[i]);
plogf("%20.12E ", gai[i]);
plogf("%3d %3d\n", m_elType[i], ElActive[i]);
plogf("%20.12E ", m_elemAbundancesGoal[i]);
plogf("%3d %3d\n", m_elType[i], m_elementActive[i]);
}
plogf("\nChemical Potentials: (J/kmol)\n");
@ -166,12 +154,9 @@ size_t VCS_SOLVE::addElement(const char* elNameNew, int elType, int elactive)
m_nelem++;
m_numComponents++;
gai.push_back(0.0);
FormulaMatrix.resize(m_nsp, m_nelem, 0.0);
m_formulaMatrix.resize(m_nsp, m_nelem);
m_formulaMatrix.resize(m_nsp, m_nelem, 0.0);
m_stoichCoeffRxnMatrix.resize(m_nelem, m_nsp, 0.0);
m_elType.push_back(elType);
ElActive.push_back(elactive);
m_elementActive.push_back(elactive);
m_elemAbundances.push_back(0.0);
m_elemAbundancesGoal.push_back(0.0);
@ -191,7 +176,7 @@ size_t VCS_SOLVE::addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT
size_t e = volPhase->elemGlobalIndex(eVP);
AssertThrowMsg(e != npos, "VCS_PROB::addOnePhaseSpecies",
"element not found");
FormulaMatrix(kT,e) = fm(k,eVP);
m_formulaMatrix(kT,e) = fm(k,eVP);
}
// Tell the phase object about the current position of the species within

111
src/equil/vcs_solve.cpp
View file

@ -53,10 +53,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_gibbsSpecies.resize(m_nsp, 0.0);
w.resize(m_nsp, 0.0);
mf.resize(m_nsp, 0.0);
SpeciesUnknownType.resize(m_nsp, VCS_SPECIES_TYPE_MOLNUM);
PhaseID.resize(m_nsp, npos);
WtSpecies.resize(m_nsp, 0.0);
Charge.resize(m_nsp, 0.0);
SpeciesThermo.resize(m_nsp,0);
for (size_t kspec = 0; kspec < m_nsp; kspec++) {
SpeciesThermo[kspec] = new VCS_SPECIES_THERMO(0, 0);
@ -115,7 +111,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_indexRxnToSpecies.resize(m_nsp, 0);
// Initialize all species to be major species
m_speciesStatus.resize(m_nsp, 1);
m_speciesStatus.resize(m_nsp, VCS_SPECIES_MAJOR);
m_SSPhase.resize(2*m_nsp, 0);
m_phaseID.resize(m_nsp, 0);
@ -233,28 +229,28 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
for (size_t k = 0; k < nSpPhase; k++) {
// Obtain the molecular weight of the species from the
// ThermoPhase object
WtSpecies[kT] = tPhase->molecularWeight(k);
m_wtSpecies[kT] = tPhase->molecularWeight(k);
// Obtain the charges of the species from the ThermoPhase object
Charge[kT] = tPhase->charge(k);
m_chargeSpecies[kT] = tPhase->charge(k);
// Set the phaseid of the species
PhaseID[kT] = iphase;
m_phaseID[kT] = iphase;
// Transfer the type of unknown
SpeciesUnknownType[kT] = VolPhase->speciesUnknownType(k);
if (SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_MOLNUM) {
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);
} else if (SpeciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
} else if (m_speciesUnknownType[kT] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
w[kT] = tPhase->electricPotential();
mf[kT] = mphase->moleFraction(kT);
} else {
throw CanteraError(" vcs_Cantera_to_vsolve() ERROR",
"Unknown species type: {}", SpeciesUnknownType[kT]);
"Unknown species type: {}", m_speciesUnknownType[kT]);
}
// transfer chemical potential vector
@ -278,7 +274,7 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
sProp->WtSpecies = tPhase->molecularWeight(k);
sProp->FormulaMatrixCol.resize(m_nelem, 0.0);
for (size_t e = 0; e < m_nelem; e++) {
sProp->FormulaMatrixCol[e] = FormulaMatrix(kT,e);
sProp->FormulaMatrixCol[e] = m_formulaMatrix(kT,e);
}
sProp->Charge = tPhase->charge(k);
sProp->SurfaceSpecies = false;
@ -364,10 +360,17 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
}
}
// Transfer initial element abundances to the vprob object.
// We have to find the mapping index from one to the other
gai.resize(m_nelem, 0.0);
set_gai();
// Transfer initial element abundances based on the species mole numbers
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];
}
}
if (m_elType[j] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[j] < 1.0E-10) {
m_elemAbundancesGoal[j] = 0.0;
}
}
// Printout the species information: PhaseID's and mole nums
if (m_printLvl > 1) {
@ -380,12 +383,12 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
plogf(" species phaseID phaseName ");
plogf(" Initial_Estimated_kMols\n");
for (size_t i = 0; i < m_nsp; i++) {
size_t iphase = PhaseID[i];
size_t iphase = m_phaseID[i];
vcs_VolPhase* VolPhase = VPhaseList[iphase];
plogf("%16s %5d %16s", mphase->speciesName(i).c_str(), iphase,
VolPhase->PhaseName.c_str());
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Volts = %-10.5g\n", w[i]);
} else {
plogf(" %-10.5g\n", w[i]);
@ -415,12 +418,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
plogf("\n");
}
// Copy over the species molecular weights
m_wtSpecies = WtSpecies;
// Copy over the charges
m_chargeSpecies = Charge;
// Copy the VCS_SPECIES_THERMO structures
for (size_t kspec = 0; kspec < m_nsp; kspec++) {
delete m_speciesThermoList[kspec];
@ -432,9 +429,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
}
}
// Copy the species unknown type
m_speciesUnknownType = SpeciesUnknownType;
// w[] -> Copy the equilibrium mole number estimate if it exists.
if (w.size() != 0) {
m_molNumSpecies_old = w;
@ -443,35 +437,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_molNumSpecies_old.assign(m_molNumSpecies_old.size(), 0.0);
}
// Formulate the Goal Element Abundance Vector
if (gai.size() != 0) {
for (size_t i = 0; i < m_nelem; i++) {
m_elemAbundancesGoal[i] = gai[i];
if (m_elType[i] == VCS_ELEM_TYPE_LATTICERATIO && m_elemAbundancesGoal[i] < 1.0E-10) {
m_elemAbundancesGoal[i] = 0.0;
}
}
} else {
if (m_doEstimateEquil == 0) {
double sum = 0;
for (size_t j = 0; j < m_nelem; j++) {
m_elemAbundancesGoal[j] = 0.0;
for (size_t kspec = 0; kspec < m_nsp; kspec++) {
if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
sum += m_molNumSpecies_old[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 * sum) {
m_elemAbundancesGoal[j] = 0.0;
}
}
} else {
throw CanteraError("VCS_SOLVE::VCS_SOLVE",
"Element Abundances, m_elemAbundancesGoal[], not specified");
}
}
// zero out values that will be filled in later
//
// TPhMoles[] -> Untouched here. These will be filled in vcs_prep.c
@ -495,24 +460,19 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_elementMapIndex[i] = i;
}
// Define all species to be major species, initially.
for (size_t i = 0; i < m_nsp; i++) {
m_speciesStatus[i] = VCS_SPECIES_MAJOR;
}
// PhaseID: Fill in the species to phase mapping. Check for bad values at
// the same time.
if (PhaseID.size() != 0) {
if (m_phaseID.size() != 0) {
std::vector<size_t> numPhSp(m_numPhases, 0);
for (size_t kspec = 0; kspec < m_nsp; kspec++) {
size_t iph = PhaseID[kspec];
size_t iph = m_phaseID[kspec];
if (iph >= m_numPhases) {
throw CanteraError("VCS_SOLVE::VCS_SOLVE",
"Species to Phase Mapping, PhaseID, has a bad value\n"
"\tPhaseID[{}] = {}\n"
"\tm_phaseID[{}] = {}\n"
"Allowed values: 0 to {}", kspec, iph, m_numPhases - 1);
}
m_phaseID[kspec] = PhaseID[kspec];
m_phaseID[kspec] = m_phaseID[kspec];
m_speciesLocalPhaseIndex[kspec] = numPhSp[iph];
numPhSp[iph]++;
}
@ -561,21 +521,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_speciesName[i] = m_mix->speciesName(i);
}
// FormulaMatrix[] -> Copy the formula matrix over
for (size_t i = 0; i < m_nsp; i++) {
bool nonzero = false;
for (size_t j = 0; j < m_nelem; j++) {
if (FormulaMatrix(i,j) != 0.0) {
nonzero = true;
}
m_formulaMatrix(i,j) = FormulaMatrix(i,j);
}
if (!nonzero) {
throw CanteraError("VCS_SOLVE::VCS_SOLVE",
"species {} {} has a zero formula matrix!", i, m_speciesName[i]);
}
}
// Copy over all of the phase information. Use the object's assignment
// operator
for (size_t iph = 0; iph < m_numPhases; iph++) {
@ -765,7 +710,7 @@ int VCS_SOLVE::vcs_prob_specifyFully()
vcs_VolPhase* VolPhase = m_VolPhaseList[iphase];
plogf("%16s %5d %16s", m_speciesName[i].c_str(), iphase,
VolPhase->PhaseName.c_str());
if (SpeciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" Volts = %-10.5g\n", w[i]);
} else {
plogf(" %-10.5g\n", w[i]);
@ -836,7 +781,7 @@ int VCS_SOLVE::vcs_prob_update()
}
// Switch the species data back from K1 into I
if (SpeciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
w[i] = m_molNumSpecies_old[k1];
} else {
w[i] = 0.0;