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

[Equil] Replace PropertyCalculator with simple lambdas

Browse source
This commit is contained in:
Ray Speth 2016-04-18 20:18:11 -04:00
parent f4ad150905
commit 3566542c0b
3 changed files with 23 additions and 129 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

6
include/cantera/equil/ChemEquil.h
View file

@ -12,6 +12,7 @@
#include "cantera/base/ct_defs.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/thermo/ThermoPhase.h"
#include <functional>
namespace Cantera
{
@ -55,9 +56,6 @@ public:
bool contin;
};
template<class M>
class PropertyCalculator;
/**
* @defgroup equil Chemical Equilibrium
*/
@ -266,7 +264,7 @@ protected:
//! it is computed. It's initialized to #m_mm.
size_t m_nComponents;
std::unique_ptr<PropertyCalculator<thermo_t> > m_p1, m_p2;
std::function<double(ThermoPhase&)> m_p1, m_p2;
//! Current value of the mole fractions in the single phase. length = #m_kk.
vector_fp m_molefractions;

44
src/equil/ChemEquil.cpp
View file

@ -7,8 +7,6 @@
// Copyright 2001 California Institute of Technology
#include "cantera/equil/ChemEquil.h"
#include "PropertyCalculator.h"
#include "cantera/base/stringUtils.h"
#include "cantera/equil/MultiPhaseEquil.h"
@ -339,37 +337,37 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
switch (XY) {
case TP:
case PT:
m_p1.reset(new TemperatureCalculator<thermo_t>);
m_p2.reset(new PressureCalculator<thermo_t>);
m_p1 = [](ThermoPhase& s) { return s.temperature(); };
m_p2 = [](ThermoPhase& s) { return s.pressure(); };
break;
case HP:
case PH:
tempFixed = false;
m_p1.reset(new EnthalpyCalculator<thermo_t>);
m_p2.reset(new PressureCalculator<thermo_t>);
m_p1 = [](ThermoPhase& s) { return s.enthalpy_mass(); };
m_p2 = [](ThermoPhase& s) { return s.pressure(); };
break;
case SP:
case PS:
tempFixed = false;
m_p1.reset(new EntropyCalculator<thermo_t>);
m_p2.reset(new PressureCalculator<thermo_t>);
m_p1 = [](ThermoPhase& s) { return s.entropy_mass(); };
m_p2 = [](ThermoPhase& s) { return s.pressure(); };
break;
case SV:
case VS:
tempFixed = false;
m_p1.reset(new EntropyCalculator<thermo_t>);
m_p2.reset(new DensityCalculator<thermo_t>);
m_p1 = [](ThermoPhase& s) { return s.entropy_mass(); };
m_p2 = [](ThermoPhase& s) { return s.density(); };
break;
case TV:
case VT:
m_p1.reset(new TemperatureCalculator<thermo_t>);
m_p2.reset(new DensityCalculator<thermo_t>);
m_p1 = [](ThermoPhase& s) { return s.temperature(); };
m_p2 = [](ThermoPhase& s) { return s.density(); };
break;
case UV:
case VU:
tempFixed = false;
m_p1.reset(new IntEnergyCalculator<thermo_t>);
m_p2.reset(new DensityCalculator<thermo_t>);
m_p1 = [](ThermoPhase& s) { return s.intEnergy_mass(); };
m_p2 = [](ThermoPhase& s) { return s.density(); };
break;
default:
throw CanteraError("equilibrate","illegal property pair.");
@ -387,8 +385,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
// Before we do anything to change the ThermoPhase object, we calculate and
// store the two specified thermodynamic properties that we are after.
double xval = m_p1->value(s);
double yval = m_p2->value(s);
double xval = m_p1(s);
double yval = m_p2(s);
size_t mm = m_mm;
size_t nvar = mm + 1;
@ -447,11 +445,11 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
s.setTemperature(tmax);
setInitialMoles(s, elMolesGoal, loglevel - 1);
phigh = m_p1->value(s);
phigh = m_p1(s);
s.setTemperature(tmin);
setInitialMoles(s, elMolesGoal, loglevel - 1);
plow = m_p1->value(s);
plow = m_p1(s);
// start with T at the midpoint of the range
doublereal t0 = 0.5*(tmin + tmax);
@ -461,7 +459,7 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
for (int it = 0; it < 10; it++) {
// set the composition and get p1
setInitialMoles(s, elMolesGoal, loglevel - 1);
pval = m_p1->value(s);
pval = m_p1(s);
// If this value of p1 is greater than the specified property value,
// then the current temperature is too high. Use it as the new upper
@ -566,8 +564,8 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr,
// check for convergence.
equilResidual(s, x, elMolesGoal, res_trial, xval, yval);
double f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin());
double xx = m_p1->value(s);
double yy = m_p2->value(s);
double xx = m_p1(s);
double yy = m_p2(s);
double deltax = (xx - xval)/xval;
double deltay = (yy - yval)/yval;
bool passThis = true;
@ -785,8 +783,8 @@ void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x,
}
}
double xx = m_p1->value(s);
double yy = m_p2->value(s);
double xx = m_p1(s);
double yy = m_p2(s);
resid[m_mm] = xx/xval - 1.0;
resid[m_skip] = yy/yval - 1.0;

102
src/equil/PropertyCalculator.h
View file

@ -1,102 +0,0 @@
/**
* @file PropertyCalculator.h
*/
// Copyright 2001 California Institute of Technology
#ifndef CT_PROP_CALC_H
#define CT_PROP_CALC_H
#include "cantera/base/ct_defs.h"
namespace Cantera
{
/// Classes used by ChemEquil. These classes are used only by the
/// ChemEquil equilibrium solver. Each one returns a particular
/// property of the object supplied as the argument.
///
template<class M>
class PropertyCalculator
{
public:
virtual ~PropertyCalculator() {}
virtual doublereal value(const M& s) =0;
virtual std::string symbol() =0;
};
template<class M>
class EnthalpyCalculator : public PropertyCalculator<M>
{
public:
virtual doublereal value(const M& s) {
return s.enthalpy_mass();
}
virtual std::string symbol() {
return "H";
}
};
template<class M>
class EntropyCalculator : public PropertyCalculator<M>
{
public:
virtual doublereal value(const M& s) {
return s.entropy_mass();
}
virtual std::string symbol() {
return "S";
}
};
template<class M>
class TemperatureCalculator : public PropertyCalculator<M>
{
public:
virtual doublereal value(const M& s) {
return s.temperature();
}
virtual std::string symbol() {
return "T";
}
};
template<class M>
class PressureCalculator : public PropertyCalculator<M>
{
public:
virtual doublereal value(const M& s) {
return s.pressure();
}
virtual std::string symbol() {
return "P";
}
};
template<class M>
class DensityCalculator : public PropertyCalculator<M>
{
public:
virtual doublereal value(const M& s) {
return s.density();
}
virtual std::string symbol() {
return "V";
}
};
template<class M>
class IntEnergyCalculator : public PropertyCalculator<M>
{
public:
virtual doublereal value(const M& s) {
return s.intEnergy_mass();
}
virtual std::string symbol() {
return "U";
}
};
}
#endif