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

[Equil] Eliminate "total moles" scaling in VCS solver

Browse source
This commit is contained in:
Ray Speth 2017-08-23 14:34:45 -04:00
parent 1b97c49d8d
commit 3820ca3147
5 changed files with 16 additions and 113 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -1426,15 +1426,6 @@ public:
*/
char m_unitsState;
//! Multiplier for the mole numbers within the nondimensional formulation
/*!
* All numbers within the main routine are on an absolute basis. This
* presents some problems wrt very large and very small mole numbers. We get
* around this by using a multiplier coming into and coming out of the
* equilibrium routines
*/
double m_totalMoleScale;
//! specifies the activity convention of the phase containing the species
/*!
* * 0 = molar based

76
src/equil/vcs_nondim.cpp
View file

@ -30,60 +30,6 @@ void VCS_SOLVE::vcs_nondim_TP()
}
m_Faraday_dim = ElectronCharge * Avogadro / (m_temperature * GasConstant);
// Scale the total moles if necessary: First find out the total moles
double tmole_orig = vcs_tmoles();
// Then add in the total moles of elements that are goals. Either one or
// the other is specified here.
double esum = 0.0;
for (size_t i = 0; i < m_nelem; ++i) {
if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) {
esum += fabs(m_elemAbundancesGoal[i]);
}
}
tmole_orig += esum;
// Ok now test out the bounds on the total moles that this program can
// handle. These are a bit arbitrary. However, it would seem that any
// reasonable input would be between these two numbers below.
if (tmole_orig < 1.0E-200 || tmole_orig > 1.0E200) {
throw CanteraError("VCS_SOLVE::vcs_nondim_TP",
"Total input moles, {} is outside the range handled by vcs.\n",
tmole_orig);
}
// Determine the scale of the problem
if (tmole_orig > 1.0E4) {
m_totalMoleScale = tmole_orig / 1.0E4;
} else if (tmole_orig < 1.0E-4) {
m_totalMoleScale = tmole_orig / 1.0E-4;
} else {
m_totalMoleScale = 1.0;
}
if (m_totalMoleScale != 1.0) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice\n", m_totalMoleScale);
}
for (size_t i = 0; i < m_nsp; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale);
}
}
for (size_t i = 0; i < m_nelem; ++i) {
m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale);
}
for (size_t iph = 0; iph < m_numPhases; iph++) {
TPhInertMoles[iph] *= (1.0 / m_totalMoleScale);
if (TPhInertMoles[iph] != 0.0) {
vcs_VolPhase* vphase = m_VolPhaseList[iph].get();
vphase->setTotalMolesInert(TPhInertMoles[iph]);
}
}
vcs_tmoles();
}
}
}
@ -103,28 +49,6 @@ void VCS_SOLVE::vcs_redim_TP()
}
m_Faraday_dim *= tf;
}
if (m_totalMoleScale != 1.0) {
if (m_debug_print_lvl >= 2) {
plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g\n", m_totalMoleScale);
}
for (size_t i = 0; i < m_nsp; ++i) {
if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
m_molNumSpecies_old[i] *= m_totalMoleScale;
}
}
for (size_t i = 0; i < m_nelem; ++i) {
m_elemAbundancesGoal[i] *= m_totalMoleScale;
}
for (size_t iph = 0; iph < m_numPhases; iph++) {
TPhInertMoles[iph] *= m_totalMoleScale;
if (TPhInertMoles[iph] != 0.0) {
vcs_VolPhase* vphase = m_VolPhaseList[iph].get();
vphase->setTotalMolesInert(TPhInertMoles[iph]);
}
}
vcs_tmoles();
}
}
}

40
src/equil/vcs_report.cpp
View file

@ -11,7 +11,7 @@ namespace Cantera
{
int VCS_SOLVE::vcs_report(int iconv)
{
bool printActualMoles = true, inertYes = false;
bool inertYes = false;
char originalUnitsState = m_unitsState;
std::vector<size_t> sortindex(m_nsp, 0);
vector_fp xy(m_nsp, 0.0);
@ -38,10 +38,6 @@ int VCS_SOLVE::vcs_report(int iconv)
if (m_unitsState == VCS_DIMENSIONAL_G) {
vcs_nondim_TP();
}
double molScale = 1.0;
if (printActualMoles) {
molScale = m_totalMoleScale;
}
vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp);
@ -66,11 +62,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature);
plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA);
plogf("\t\ttotal Volume = %15.5g m**3\n", m_totalVol * molScale);
if (!printActualMoles) {
plogf("\t\tMole Scale = %15.5g kmol (all mole numbers and volumes are scaled by this value)\n",
molScale);
}
plogf("\t\ttotal Volume = %15.5g m**3\n", m_totalVol);
// TABLE OF SPECIES IN DECREASING MOLE NUMBERS
plogf("\n\n");
@ -81,8 +73,8 @@ int VCS_SOLVE::vcs_report(int iconv)
for (size_t i = 0; i < m_numComponents; ++i) {
plogf(" %-12.12s", m_speciesName[i]);
writeline(' ', 13, false);
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i] * molScale,
m_molNumSpecies_new[i] * molScale, m_feSpecies_old[i]);
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i],
m_molNumSpecies_new[i], m_feSpecies_old[i]);
plogf(" %3d", m_speciesUnknownType[i]);
plogf("\n");
}
@ -92,12 +84,12 @@ int VCS_SOLVE::vcs_report(int iconv)
writeline(' ', 13, false);
if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_MOLNUM) {
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[j] * molScale,
m_molNumSpecies_new[j] * molScale, m_feSpecies_old[j]);
plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[j],
m_molNumSpecies_new[j], m_feSpecies_old[j]);
plogf(" KMolNum ");
} else if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_old[j]);
plogf(" Voltage = %14.7E", m_molNumSpecies_old[j] * molScale);
plogf(" Voltage = %14.7E", m_molNumSpecies_old[j]);
} else {
throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem");
}
@ -112,7 +104,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" Inert Species in phase %16s ",
m_VolPhaseList[i]->PhaseName);
}
plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i] * molScale,
plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i],
TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0);
}
}
@ -122,8 +114,8 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %-12.12s", m_speciesName[kspec]);
// Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve
plogf(" %14.7E %14.7E %12.4E",
m_molNumSpecies_old[kspec]*molScale,
m_molNumSpecies_new[kspec]*molScale, m_deltaGRxn_new[kspec]);
m_molNumSpecies_old[kspec],
m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]);
if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) {
plogf(" KMol_Num");
} else if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
@ -151,7 +143,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" | |\n");
plogf(" NonComponent | Moles |");
for (size_t j = 0; j < m_numComponents; j++) {
plogf(" %10.3g", m_molNumSpecies_old[j] * molScale);
plogf(" %10.3g", m_molNumSpecies_old[j]);
}
plogf(" | DG/RT Rxn |\n");
writeline('-', m_numComponents*10 + 45);
@ -159,7 +151,7 @@ int VCS_SOLVE::vcs_report(int iconv)
size_t kspec = m_indexRxnToSpecies[irxn];
plogf(" %3d ", kspec);
plogf("%-10.10s", m_speciesName[kspec]);
plogf("|%10.3g |", m_molNumSpecies_old[kspec]*molScale);
plogf("|%10.3g |", m_molNumSpecies_old[kspec]);
for (size_t j = 0; j < m_numComponents; j++) {
plogf(" %6.2f", m_stoichCoeffRxnMatrix(j,irxn));
}
@ -198,7 +190,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" %3d ", iphase);
vcs_VolPhase* VPhase = m_VolPhaseList[iphase].get();
plogf("%-12.12s |",VPhase->PhaseName);
plogf("%10.3e |", m_tPhaseMoles_old[iphase]*molScale);
plogf("%10.3e |", m_tPhaseMoles_old[iphase]);
totalMoles += m_tPhaseMoles_old[iphase];
if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles() &&
!vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) {
@ -240,7 +232,7 @@ int VCS_SOLVE::vcs_report(int iconv)
for (size_t i = 0; i < m_nelem; ++i) {
writeline(' ', 26, false);
plogf("%-2.2s", m_elementName[i]);
plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale);
plogf("%20.12E %20.12E", m_elemAbundances[i], m_elemAbundancesGoal[i]);
plogf(" %3d %3d\n", m_elType[i], m_elementActive[i]);
}
plogf("\n");
@ -258,7 +250,7 @@ int VCS_SOLVE::vcs_report(int iconv)
size_t j = sortindex[i];
size_t pid = m_phaseID[j];
plogf(" %-12.12s", m_speciesName[j]);
plogf(" %14.7E ", m_molNumSpecies_old[j]*molScale);
plogf(" %14.7E ", m_molNumSpecies_old[j]);
plogf("%14.7E ", m_SSfeSpecies[j]);
plogf("%14.7E ", log(m_actCoeffSpecies_old[j]));
double tpmoles = m_tPhaseMoles_old[pid];
@ -291,7 +283,7 @@ int VCS_SOLVE::vcs_report(int iconv)
plogf(" ");
}
plogf("| %20.9E |", m_feSpecies_old[j] * m_molNumSpecies_old[j] * molScale);
plogf("| %20.9E |", m_feSpecies_old[j] * m_molNumSpecies_old[j]);
plogf("\n");
}
for (size_t i = 0; i < 125; i++) {

1
src/equil/vcs_solve.cpp
View file

@ -41,7 +41,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) :
m_tolmaj2(1.0E-10),
m_tolmin2(1.0E-8),
m_unitsState(VCS_DIMENSIONAL_G),
m_totalMoleScale(1.0),
m_useActCoeffJac(0),
m_totalVol(mphase->volume()),
m_Faraday_dim(ElectronCharge * Avogadro),

3
src/equil/vcs_solve_TP.cpp
View file

@ -942,9 +942,6 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1,
} else {
plogf(" (only major species):");
}
if (m_totalMoleScale != 1.0) {
plogf(" (Total Mole Scale = %g)", m_totalMoleScale);
}
plogf("\n");
plogf(" --- Species Status Initial_KMoles Final_KMoles Initial_Mu/RT");
plogf(" Mu/RT Init_Del_G/RT Delta_G/RT\n");