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

[Thermo] 'report' only prints species with non-trivial mole fractions

Browse source 
The threshold for printing species is a new, optional parameter to the
'report' function. If any species are excluded, the total number of
minor species and their aggregate mass and mole fraction are printed.
This commit is contained in:
Ray Speth 2014-06-03 16:54:15 +00:00
parent dbfe428437
commit ea3f74eb27
21 changed files with 120 additions and 61 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
include/cantera/thermo/MolalityVPSSTP.h
View file

@ -698,8 +698,11 @@ public:
/*!
* @param show_thermo If true, extra information is printed out
* about the thermodynamic state of the system.
* @param threshold Show information about species with mole fractions
* greater than *threshold*.
*/
virtual std::string report(bool show_thermo = true) const;
virtual std::string report(bool show_thermo=true,
doublereal threshold=1e-14) const;
protected:

5
include/cantera/thermo/MolarityIonicVPSSTP.h
View file

@ -276,8 +276,11 @@ public:
/*!
* @param show_thermo If true, extra information is printed out
* about the thermodynamic state of the system.
* @param threshold Show information about species with mole fractions
* greater than *threshold*.
*/
virtual std::string report(bool show_thermo = true) const;
virtual std::string report(bool show_thermo=true,
doublereal threshold=1e-14) const;
private:
//! Initialize lengths of local variables after all species have been identified.

3
include/cantera/thermo/PseudoBinaryVPSSTP.h
View file

@ -187,7 +187,8 @@ public:
* @param show_thermo If true, extra information is printed out
* about the thermodynamic state of the system.
*/
virtual std::string report(bool show_thermo = true) const;
virtual std::string report(bool show_thermo=true,
doublereal threshold=1e-14) const;
private:
//! Initialize lengths of local variables after all species have

4
include/cantera/thermo/PureFluidPhase.h
View file

@ -499,8 +499,10 @@ public:
/*!
* @param show_thermo If true, extra information is printed out
* about the thermodynamic state of the system.
* @param threshold Unused in this subclass
*/
virtual std::string report(bool show_thermo = true) const;
virtual std::string report(bool show_thermo=true,
doublereal threshold=1e-14) const;
protected:

5
include/cantera/thermo/ThermoPhase.h
View file

@ -1535,8 +1535,11 @@ public:
/*!
* @param show_thermo If true, extra information is printed out
* about the thermodynamic state of the system.
* @param threshold Show information about species with mole fractions
* greater than *threshold*.
*/
virtual std::string report(bool show_thermo = true) const;
virtual std::string report(bool show_thermo=true,
doublereal threshold=1e-14) const;
//! returns a summary of the state of the phase to a comma separated file.
//! To customize the data included in the report, derived classes should

2
interfaces/cython/cantera/_cantera.pxd
View file

@ -46,7 +46,7 @@ cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera":
# miscellaneous
int eosType()
string report(cbool) except +
string report(cbool, double) except +
string name()
void setName(string)
string id()

4
interfaces/cython/cantera/mixture.pyx
View file

@ -48,7 +48,7 @@ cdef class Mixture:
def __dealloc__(self):
del self.mix
def report(self):
def report(self, threshold=1e-14):
"""
Generate a report describing the thermodynamic state of this mixture. To
print the report to the screen, simply call the mixture object. The
@ -61,7 +61,7 @@ cdef class Mixture:
for i,phase in enumerate(self._phases):
s.append('************ Phase {0} ************'.format(phase.name))
s.append('Moles: {0}'.format(self.phase_moles(i)))
s.append(phase.report())
s.append(phase.report(threshold=threshold))
return '\n'.join(s)

20
interfaces/cython/cantera/test/test_thermo.py
View file

@ -114,12 +114,22 @@ class TestThermoPhase(utilities.CanteraTest):
with self.assertRaises(ValueError):
self.phase['H2','O2'].Y = [0.1, 0.2, 0.3]
def test_report(self):
report = self.phase.report()
self.assertTrue(self.phase.name in report)
self.assertTrue('temperature' in report)
def test_full_report(self):
report = self.phase.report(threshold=0.0)
self.assertIn(self.phase.name, report)
self.assertIn('temperature', report)
self.assertNotIn('minor', report)
for name in self.phase.species_names:
self.assertTrue(name in report)
self.assertIn(name, report)
def test_default_report(self):
self.phase.X = 'H2:0.1, O2:0.9, HO2:1e-10, H2O2:1e-20'
report = self.phase.report()
self.assertIn('minor', report)
for name in (' H2 ', ' O2 ', ' HO2 '):
self.assertIn(name, report)
for name in (' H2O2 ', ' OH ', ' AR '):
self.assertNotIn(name, report)
def test_name(self):
self.assertEqual(self.phase.name, 'ohmech')

8
interfaces/cython/cantera/thermo.pyx
View file

@ -27,7 +27,7 @@ cdef class ThermoPhase(_SolutionBase):
if 'source' not in kwargs:
self.thermo_basis = mass_basis
def report(self, show_thermo=True):
def report(self, show_thermo=True, float threshold=1e-14):
"""
Generate a report describing the thermodynamic state of this phase. To
print the report to the terminal, simply call the phase object. The
@ -36,10 +36,10 @@ cdef class ThermoPhase(_SolutionBase):
>>> phase()
>>> print(phase.report())
"""
return pystr(self.thermo.report(bool(show_thermo)))
return pystr(self.thermo.report(bool(show_thermo), threshold))
def __call__(self):
print(self.report())
def __call__(self, *args, **kwargs):
print(self.report(*args, **kwargs))
property name:
"""

7
interfaces/matlab/toolbox/@ThermoPhase/display.m
View file

@ -1,2 +1,5 @@
function display(self)
phase_get(thermo_hndl(self), 15, 1);
function display(self, threshold)
if nargin < 2
threshold = 1e-14
end
phase_get(thermo_hndl(self), 15, 1, threshold);

4
src/clib/ct.cpp
View file

@ -1336,11 +1336,11 @@ extern "C" {
}
}
int write_phase(int nth, int show_thermo)
int write_phase(int nth, int show_thermo, double threshold)
{
try {
bool stherm = (show_thermo != 0);
writelog(ThermoCabinet::item(nth).report(stherm)+"\n");
writelog(ThermoCabinet::item(nth).report(stherm, threshold)+"\n");
return 0;
} catch (...) {
return handleAllExceptions(-1, ERR);

2
src/clib/ct.h
View file

@ -38,7 +38,7 @@ extern "C" {
CANTERA_CAPI size_t phase_speciesIndex(int n, char* nm);
CANTERA_CAPI int phase_report(int nth,
int ibuf, char* buf, int show_thermo);
CANTERA_CAPI int write_phase(int nth, int show_thermo);
CANTERA_CAPI int write_phase(int nth, int show_thermo, double threshold);
CANTERA_CAPI double phase_nAtoms(int n, size_t k, size_t m);

4
src/matlab/phasemethods.cpp
View file

@ -110,6 +110,7 @@ void phasemethods(int nlhs, mxArray* plhs[],
}
else if (job < 20) {
double threshold;
switch (job) {
case 0:
@ -152,7 +153,8 @@ void phasemethods(int nlhs, mxArray* plhs[],
break;
case 15:
show_thermo = getInt(prhs[3]);
vv = write_phase(ph,show_thermo);
threshold = getDouble(prhs[4]);
vv = write_phase(ph,show_thermo,threshold);
break;
default:
mexErrMsgTxt("Unknown job number");

36
src/thermo/MolalityVPSSTP.cpp
View file

@ -502,7 +502,7 @@ void MolalityVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
/**
* Format a summary of the mixture state for output.
*/
std::string MolalityVPSSTP::report(bool show_thermo) const
std::string MolalityVPSSTP::report(bool show_thermo, doublereal threshold) const
{
char p[800];
string s = "";
@ -580,6 +580,8 @@ std::string MolalityVPSSTP::report(bool show_thermo) const
sprintf(p, " \n");
s += p;
int nMinor = 0;
doublereal xMinor = 0.0;
if (show_thermo) {
sprintf(p, " X "
" Molalities Chem.Pot. ChemPotSS ActCoeffMolal\n");
@ -591,14 +593,19 @@ std::string MolalityVPSSTP::report(bool show_thermo) const
" ------------ ------------ ------------ ------------\n");
s += p;
for (size_t k = 0; k < kk; k++) {
if (x[k] > SmallNumber) {
sprintf(p, "%18s %12.6g %12.6g %12.6g %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], molal[k], mu[k], muss[k], acMolal[k]);
if (x[k] > threshold) {
if (x[k] > SmallNumber) {
sprintf(p, "%18s %12.6g %12.6g %12.6g %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], molal[k], mu[k], muss[k], acMolal[k]);
} else {
sprintf(p, "%18s %12.6g %12.6g N/A %12.6g %12.6g \n",
speciesName(k).c_str(), x[k], molal[k], muss[k], acMolal[k]);
}
s += p;
} else {
sprintf(p, "%18s %12.6g %12.6g N/A %12.6g %12.6g \n",
speciesName(k).c_str(), x[k], molal[k], muss[k], acMolal[k]);
nMinor++;
xMinor += x[k];
}
s += p;
}
} else {
sprintf(p, " X"
@ -608,11 +615,20 @@ std::string MolalityVPSSTP::report(bool show_thermo) const
" ------------\n");
s += p;
for (size_t k = 0; k < kk; k++) {
sprintf(p, "%18s %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], molal[k]);
s += p;
if (x[k] > threshold) {
sprintf(p, "%18s %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], molal[k]);
s += p;
} else {
nMinor++;
xMinor += x[k];
}
}
}
if (nMinor) {
sprintf(p, " [%+5i minor] %12.6g\n", nMinor, xMinor);
s += p;
}
} catch (CanteraError& err) {
err.save();
}

2
src/thermo/MolarityIonicVPSSTP.cpp
View file

@ -457,7 +457,7 @@ void MolarityIonicVPSSTP::readXMLBinarySpecies(XML_Node& xmLBinarySpecies)
std::string xname = xmLBinarySpecies.name();
}
std::string MolarityIonicVPSSTP::report(bool show_thermo) const
std::string MolarityIonicVPSSTP::report(bool show_thermo, doublereal threshold) const
{
char p[800];
string s = "";

2
src/thermo/PseudoBinaryVPSSTP.cpp
View file

@ -175,7 +175,7 @@ void PseudoBinaryVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& i
GibbsExcessVPSSTP::initThermoXML(phaseNode, id_);
}
std::string PseudoBinaryVPSSTP::report(bool show_thermo) const
std::string PseudoBinaryVPSSTP::report(bool show_thermo, doublereal threshold) const
{
char p[800];
string s = "";

2
src/thermo/PureFluidPhase.cpp
View file

@ -362,7 +362,7 @@ void PureFluidPhase::setState_Psat(doublereal p, doublereal x)
setDensity(1.0/m_sub->v());
}
std::string PureFluidPhase::report(bool show_thermo) const
std::string PureFluidPhase::report(bool show_thermo, doublereal threshold) const
{
char p[800];
string s = "";

40
src/thermo/ThermoPhase.cpp
View file

@ -845,7 +845,7 @@ void ThermoPhase::getdlnActCoeffdlnN_numderiv(const size_t ld, doublereal* const
setState_PX(pres, DATA_PTR(Xmol_Base));
}
std::string ThermoPhase::report(bool show_thermo) const
std::string ThermoPhase::report(bool show_thermo, doublereal threshold) const
{
char p[800];
string s = "";
@ -910,6 +910,9 @@ std::string ThermoPhase::report(bool show_thermo) const
getChemPotentials(&mu[0]);
doublereal rt = GasConstant * temperature();
int nMinor = 0;
doublereal xMinor = 0.0;
doublereal yMinor = 0.0;
if (show_thermo) {
sprintf(p, " \n X "
" Y Chem. Pot. / RT \n");
@ -918,14 +921,20 @@ std::string ThermoPhase::report(bool show_thermo) const
"------------ ------------\n");
s += p;
for (size_t k = 0; k < kk; k++) {
if (x[k] > SmallNumber) {
sprintf(p, "%18s %12.6g %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], y[k], mu[k]/rt);
if (x[k] >= threshold) {
if (x[k] > SmallNumber) {
sprintf(p, "%18s %12.6g %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], y[k], mu[k]/rt);
} else {
sprintf(p, "%18s %12.6g %12.6g \n",
speciesName(k).c_str(), x[k], y[k]);
}
s += p;
} else {
sprintf(p, "%18s %12.6g %12.6g \n",
speciesName(k).c_str(), x[k], y[k]);
nMinor++;
xMinor += x[k];
yMinor += y[k];
}
s += p;
}
} else {
sprintf(p, " \n X "
@ -935,11 +944,22 @@ std::string ThermoPhase::report(bool show_thermo) const
" ------------\n");
s += p;
for (size_t k = 0; k < kk; k++) {
sprintf(p, "%18s %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], y[k]);
s += p;
if (x[k] >= threshold) {
sprintf(p, "%18s %12.6g %12.6g\n",
speciesName(k).c_str(), x[k], y[k]);
s += p;
} else {
nMinor++;
xMinor += x[k];
yMinor += y[k];
}
}
}
if (nMinor) {
sprintf(p, " [%+5i minor] %12.6g %12.6g\n",
nMinor, xMinor, yMinor);
s += p;
}
}
catch (CanteraError& err) {
err.save();

18
test_problems/ChemEquil_gri_matrix/gri_matrix.cpp
View file

@ -39,7 +39,7 @@ int main(int argc, char** argv)
g.setState_TPX(T, pres, DATA_PTR(Xmol));
try {
retnSub = equilibrate(g, "TP", solver);
cout << g;
cout << g.report(true, 0.0);
if (retnSub != 1) {
cerr << "ERROR: ChemEquil equilibration step failed at "
<< " T = " << T
@ -52,7 +52,7 @@ int main(int argc, char** argv)
exit(-1);
}
} catch (CanteraError& err) {
cout << g;
cout << g.report(true, 0.0);
std::cerr << err.what() << std::endl;
cerr << "ERROR: equilibration step failed at "
<< " T = " << T
@ -75,7 +75,7 @@ int main(int argc, char** argv)
g.setState_TPX(T, pres, DATA_PTR(Xmol));
try {
retnSub = equilibrate(g, "TP", solver);
cout << g;
cout << g.report(true, 0.0);
if (retnSub != 1) {
cerr << "ERROR: ChemEquil equilibration step failed at "
<< " T = " << T
@ -88,7 +88,7 @@ int main(int argc, char** argv)
exit(-1);
}
} catch (CanteraError& err) {
cout << g;
cout << g.report(true, 0.0);
std::cerr << err.what() << std::endl;
cerr << "ERROR: equilibration step failed at "
<< " T = " << T
@ -113,7 +113,7 @@ int main(int argc, char** argv)
g.setState_TPX(T, pres, DATA_PTR(Xmol));
try {
retnSub = equilibrate(g, "TP", solver);
cout << g;
cout << g.report(true, 0.0);
if (retnSub != 1) {
cerr << "ERROR: ChemEquil equilibration step failed at "
<< " T = " << T
@ -126,7 +126,7 @@ int main(int argc, char** argv)
exit(-1);
}
} catch (CanteraError& err) {
cout << g;
cout << g.report(true, 0.0);
std::cerr << err.what() << std::endl;
cerr << "ERROR: equilibration step failed at "
<< " T = " << T
@ -151,7 +151,7 @@ int main(int argc, char** argv)
g.setState_TPX(T, pres, DATA_PTR(Xmol));
try {
retnSub = equilibrate(g, "TP", solver);
cout << g;
cout << g.report(true, 0.0);
if (retnSub != 1) {
cerr << "ERROR: ChemEquil equilibration step failed at "
<< " T = " << T
@ -164,7 +164,7 @@ int main(int argc, char** argv)
exit(-1);
}
} catch (CanteraError& err) {
cout << g;
cout << g.report(true, 0.0);
std::cerr << err.what() << std::endl;
cerr << "ERROR: equilibration step failed at "
<< " T = " << T
@ -219,7 +219,7 @@ int main(int argc, char** argv)
<< endl;
}
if (showResults) {
cout << g;
cout << g.report(true, 0.0);
}
if (retnSub ==1) {
numSucc++;

4
test_problems/VCSnonideal/NaCl_equil/good_out.txt
View file

@ -359,13 +359,11 @@ Moles: 4.09718
X Y Chem. Pot. / RT
------------- ------------ ------------
O2 2.18341e-69 2.51533e-69 -182.771
H2 4.57268e-08 3.31864e-09 -32.6178
CO2 0 0
H2O 0.0237153 0.0153814 -124.003
NaCl 9.76108e-33 2.05377e-32 -174.517
N2 0.976285 0.984619 -23.0692
OH 9.14536e-08 5.59966e-08 -107.694
[ +3 minor] 9.76108e-33 2.05377e-32
*************** NaCl(S) *****************
Moles: 4.7877

4
test_problems/VPsilane_test/silane_equil.cpp
View file

@ -23,9 +23,7 @@ int main(int argc, char** argv)
g->setState_TPX(1500.0, 100.0, "SIH4:0.01, H2:0.99");
//g.setState_TPX(1500.0, 1.0132E5, "SIH4:0.01, H2:0.99");
equilibrate(*g, "TP");
std::string r = g->report(true);
cout << r;
cout << endl;
cout << g->report(true, 0.0) << endl;
return 0;
} catch (CanteraError& err) {
std::cerr << err.what() << std::endl;