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

Eliminate unnecessary counter member variables

Browse source
This commit is contained in:
Ray Speth 2015-12-01 15:54:55 -05:00
parent 2dfa43ffa7
commit 4887775109
21 changed files with 141 additions and 222 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
include/cantera/equil/MultiPhase.h
View file

@ -426,7 +426,7 @@ public:
//! Number of phases.
size_t nPhases() const {
return m_np;
return m_phase.size();
}
//! Return true is species \a kGlob is a species in a multicomponent
@ -617,9 +617,6 @@ private:
*/
std::map<std::string, size_t> m_enamemap;
//! Number of phases in the MultiPhase object
size_t m_np;
//! Current value of the temperature (kelvin)
doublereal m_temp;

2
include/cantera/equil/MultiPhaseEquil.h
View file

@ -161,7 +161,7 @@ protected:
return (m_nsp > m_nel) ? m_nsp - m_nel : 0;
}
size_t m_nel_mix, m_nsp_mix, m_np;
size_t m_nel_mix, m_nsp_mix;
size_t m_nel, m_nsp;
size_t m_eloc;
int m_iter;

2
include/cantera/kinetics/GasKinetics.h
View file

@ -70,8 +70,6 @@ public:
virtual void update_rates_C();
protected:
size_t m_nfall;
//! Reaction index of each falloff reaction
std::vector<size_t> m_fallindx;

10
include/cantera/kinetics/ImplicitSurfChem.h
View file

@ -227,22 +227,12 @@ protected:
std::vector<size_t> m_specStartIndex;
//! Total number of surface phases.
/*!
* This is also equal to the number of InterfaceKinetics objects
* as there is a 1-1 correspondence between InterfaceKinetics objects
* and surface phases.
*/
size_t m_nsurf;
//! Total number of surface species in all surface phases
/*!
* This is the total number of unknowns in m_mode 0 problem
*/
size_t m_nv;
size_t m_numBulkPhases;
std::vector<size_t> m_nspBulkPhases;
size_t m_numTotalBulkSpecies;
size_t m_numTotalSpecies;

6
include/cantera/kinetics/InterfaceKinetics.h
View file

@ -404,12 +404,6 @@ protected:
*/
std::vector<size_t> m_irrev;
//! Number of irreversible reactions in the mechanism
size_t m_nirrev;
//! Number of reversible reactions in the mechanism
size_t m_nrev;
//! Array of concentrations for each species in the kinetics mechanism
/*!
* An array of generalized concentrations \f$ C_k \f$ that are defined

5
include/cantera/kinetics/RxnRates.h
View file

@ -129,12 +129,12 @@ public:
m_mcov = 0.0;
size_t k;
doublereal th;
for (size_t n = 0; n < m_ncov; n++) {
for (size_t n = 0; n < m_ac.size(); n++) {
k = m_sp[n];
m_acov += m_ac[n] * theta[k];
m_ecov += m_ec[n] * theta[k];
}
for (size_t n = 0; n < m_nmcov; n++) {
for (size_t n = 0; n < m_mc.size(); n++) {
k = m_msp[n];
th = std::max(theta[k], Tiny);
m_mcov += m_mc[n]*std::log(th);
@ -177,7 +177,6 @@ protected:
doublereal m_acov, m_ecov, m_mcov;
std::vector<size_t> m_sp, m_msp;
vector_fp m_ac, m_ec, m_mc;
size_t m_ncov, m_nmcov;
};

21
include/cantera/numerics/Func1.h
View file

@ -820,9 +820,8 @@ class Poly1 : public Func1
public:
Poly1(size_t n, doublereal* c) :
Func1() {
m_n = n+1;
m_cpoly.resize(n+1);
std::copy(c, c+m_n, m_cpoly.begin());
std::copy(c, c+m_cpoly.size(), m_cpoly.begin());
}
Poly1(const Poly1& b) :
@ -836,7 +835,6 @@ public:
}
Func1::operator=(right);
m_cpoly = right.m_cpoly;
m_n = right.m_n;
m_parent = 0;
return *this;
}
@ -847,16 +845,15 @@ public:
}
virtual doublereal eval(doublereal t) const {
doublereal r = m_cpoly[m_n-1];
for (size_t n = 1; n < m_n; n++) {
doublereal r = m_cpoly[m_cpoly.size()-1];
for (size_t n = 1; n < m_cpoly.size(); n++) {
r *= t;
r += m_cpoly[m_n - n - 1];
r += m_cpoly[m_cpoly.size() - n - 1];
}
return r;
}
protected:
size_t m_n;
vector_fp m_cpoly;
};
@ -875,7 +872,6 @@ public:
Fourier1(size_t n, doublereal omega, doublereal a0,
doublereal* a, doublereal* b) :
Func1() {
m_n = n;
m_omega = omega;
m_a0_2 = 0.5*a0;
m_ccos.resize(n);
@ -898,7 +894,6 @@ public:
m_a0_2 = right.m_a0_2;
m_ccos = right.m_ccos;
m_csin = right.m_csin;
m_n = right.m_n;
m_parent = 0;
return *this;
}
@ -911,7 +906,7 @@ public:
virtual doublereal eval(doublereal t) const {
size_t n, nn;
doublereal sum = m_a0_2;
for (n = 0; n < m_n; n++) {
for (n = 0; n < m_ccos.size(); n++) {
nn = n + 1;
sum += m_ccos[n]*std::cos(m_omega*nn*t)
+ m_csin[n]*std::sin(m_omega*nn*t);
@ -920,7 +915,6 @@ public:
}
protected:
size_t m_n;
doublereal m_omega, m_a0_2;
vector_fp m_ccos, m_csin;
};
@ -937,7 +931,6 @@ class Arrhenius1 : public Func1
public:
Arrhenius1(size_t n, doublereal* c) :
Func1() {
m_n = n;
m_A.resize(n);
m_b.resize(n);
m_E.resize(n);
@ -959,7 +952,6 @@ public:
return *this;
}
Func1::operator=(right);
m_n = right.m_n;
m_A = right.m_A;
m_b = right.m_b;
m_E = right.m_E;
@ -974,14 +966,13 @@ public:
virtual doublereal eval(doublereal t) const {
doublereal sum = 0.0;
for (size_t n = 0; n < m_n; n++) {
for (size_t n = 0; n < m_A.size(); n++) {
sum += m_A[n]*std::pow(t,m_b[n])*std::exp(-m_E[n]/t);
}
return sum;
}
protected:
size_t m_n;
vector_fp m_A, m_b, m_E;
};

13
include/cantera/oneD/OneDim.h
View file

@ -47,7 +47,7 @@ public:
/// Number of domains.
size_t nDomains() const {
return m_nd;
return m_dom.size();
}
/// Return a reference to domain i.
@ -60,8 +60,8 @@ public:
//! Check that the specified domain index is in range.
//! Throws an exception if n is greater than nDomains()-1
void checkDomainIndex(size_t n) const {
if (n >= m_nd) {
throw IndexError("checkDomainIndex", "domains", n, m_nd-1);
if (n >= m_dom.size()) {
throw IndexError("checkDomainIndex", "domains", n, m_dom.size()-1);
}
}
@ -69,8 +69,8 @@ public:
//! Throws an exception if nn is less than nDomains(). Used before calls
//! which take an array pointer.
void checkDomainArraySize(size_t nn) const {
if (m_nd > nn) {
throw ArraySizeError("checkDomainArraySize", nn, m_nd);
if (m_dom.size() > nn) {
throw ArraySizeError("checkDomainArraySize", nn, m_dom.size());
}
}
@ -264,9 +264,6 @@ protected:
doublereal m_rdt; //!< reciprocal of time step
bool m_jac_ok; //!< if true, Jacobian is current
//! number of domains
size_t m_nd;
size_t m_bw; //!< Jacobian bandwidth
size_t m_size; //!< solution vector size

18
include/cantera/thermo/AdsorbateThermo.h
View file

@ -30,9 +30,7 @@ class Adsorbate : public SpeciesThermoInterpType
{
public:
//! Empty constructor
Adsorbate() :
m_nFreqs(0) {
}
Adsorbate() {}
//! Full Constructor
/*!
@ -43,10 +41,9 @@ public:
Adsorbate(double tlow, double thigh, double pref, const double* coeffs)
: SpeciesThermoInterpType(tlow, thigh, pref)
{
m_nFreqs = int(coeffs[0]);
m_freq.resize(int(coeffs[0]));
m_be = coeffs[1];
m_freq.resize(m_nFreqs);
std::copy(coeffs+2, coeffs + 2 + m_nFreqs, m_freq.begin());
std::copy(coeffs+2, coeffs + 2 + m_freq.size(), m_freq.begin());
}
virtual SpeciesThermoInterpType*
@ -76,15 +73,14 @@ public:
tlow = m_lowT;
thigh = m_highT;
pref = m_Pref;
coeffs[0] = static_cast<double>(m_nFreqs);
coeffs[0] = static_cast<double>(m_freq.size());
coeffs[1] = m_be;
for (size_t i = 2; i < m_nFreqs+2; i++) {
for (size_t i = 2; i < m_freq.size()+2; i++) {
coeffs[i] = m_freq[i-2];
}
}
protected:
size_t m_nFreqs;
//! array of vib frequencies
vector_fp m_freq;
doublereal m_be;
@ -92,7 +88,7 @@ protected:
doublereal _energy_RT(double T) const {
doublereal x, hnu_kt, hnu, sum = 0.0;
doublereal kt = T*Boltzmann;
for (size_t i = 0; i < m_nFreqs; i++) {
for (size_t i = 0; i < m_freq.size(); i++) {
hnu = Planck * m_freq[i];
hnu_kt = hnu/kt;
x = exp(-hnu_kt);
@ -104,7 +100,7 @@ protected:
doublereal _free_energy_RT(double T) const {
doublereal x, hnu_kt, sum = 0.0;
doublereal kt = T*Boltzmann;
for (size_t i = 0; i < m_nFreqs; i++) {
for (size_t i = 0; i < m_freq.size(); i++) {
hnu_kt = Planck * m_freq[i] / kt;
x = exp(-hnu_kt);
sum += log(1.0 - x);

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

@ -453,9 +453,6 @@ protected:
//! Current value of the molar density
doublereal m_molar_density;
//! Number of sublattice phases
size_t m_nlattice;
//! Vector of sublattic ThermoPhase objects
std::vector<LatticePhase*> m_lattice;

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

@ -98,9 +98,6 @@ public:
virtual void modifyParameters(doublereal* coeffs);
protected:
//! Number of temperature regions
size_t m_numTempRegions;
//! Lower boundaries of each temperature regions
vector_fp m_lowerTempBounds;

43
src/equil/MultiPhase.cpp
View file

@ -16,7 +16,6 @@ namespace Cantera
{
MultiPhase::MultiPhase() :
m_np(0),
m_temp(298.15),
m_press(OneBar),
m_nel(0),
@ -29,7 +28,6 @@ MultiPhase::MultiPhase() :
}
MultiPhase::MultiPhase(const MultiPhase& right) :
m_np(0),
m_temp(298.15),
m_press(OneBar),
m_nel(0),
@ -54,7 +52,6 @@ MultiPhase& MultiPhase::operator=(const MultiPhase& right)
m_spstart = right.m_spstart;
m_enames = right.m_enames;
m_enamemap = right.m_enamemap;
m_np = right.m_np;
m_temp = right.m_temp;
m_press = right.m_press;
m_nel = right.m_nel;
@ -72,7 +69,7 @@ MultiPhase& MultiPhase::operator=(const MultiPhase& right)
void MultiPhase::addPhases(MultiPhase& mix)
{
size_t n;
for (n = 0; n < mix.m_np; n++) {
for (n = 0; n < mix.nPhases(); n++) {
addPhase(mix.m_phase[n], mix.m_moles[n]);
}
}
@ -102,9 +99,7 @@ void MultiPhase::addPhase(ThermoPhase* p, doublereal moles)
m_moles.push_back(moles);
m_temp_OK.push_back(true);
// update the number of phases and the total number of
// species
m_np = m_phase.size();
// update the total number of species
m_nsp += p->nSpecies();
// determine if this phase has new elements for each new element, add an
@ -174,7 +169,7 @@ void MultiPhase::init()
sym = m_enames[m];
k = 0;
// iterate over the phases
for (ip = 0; ip < m_np; ip++) {
for (ip = 0; ip < nPhases(); ip++) {
ThermoPhase* p = m_phase[ip];
nsp = p->nSpecies();
mlocal = p->elementIndex(sym);
@ -248,7 +243,7 @@ doublereal MultiPhase::elementMoles(size_t m) const
{
doublereal sum = 0.0, phasesum;
size_t i, k = 0, ik, nsp;
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
phasesum = 0.0;
nsp = m_phase[i]->nSpecies();
for (ik = 0; ik < nsp; ik++) {
@ -264,7 +259,7 @@ doublereal MultiPhase::charge() const
{
doublereal sum = 0.0;
size_t i;
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
sum += phaseCharge(i);
}
return sum;
@ -301,7 +296,7 @@ void MultiPhase::getChemPotentials(doublereal* mu) const
{
size_t i, loc = 0;
updatePhases();
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
m_phase[i]->getChemPotentials(mu + loc);
loc += m_phase[i]->nSpecies();
}
@ -313,7 +308,7 @@ void MultiPhase::getValidChemPotentials(doublereal not_mu,
size_t i, loc = 0;
updatePhases();
// iterate over the phases
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
if (tempOK(i) || m_phase[i]->nSpecies() > 1) {
if (!standard) {
m_phase[i]->getChemPotentials(mu + loc);
@ -341,7 +336,7 @@ doublereal MultiPhase::gibbs() const
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->gibbs_mole() * m_moles[i];
}
@ -354,7 +349,7 @@ doublereal MultiPhase::enthalpy() const
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->enthalpy_mole() * m_moles[i];
}
@ -367,7 +362,7 @@ doublereal MultiPhase::IntEnergy() const
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->intEnergy_mole() * m_moles[i];
}
@ -380,7 +375,7 @@ doublereal MultiPhase::entropy() const
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->entropy_mole() * m_moles[i];
}
@ -393,7 +388,7 @@ doublereal MultiPhase::cp() const
size_t i;
doublereal sum = 0.0;
updatePhases();
for (i = 0; i < m_np; i++) {
for (i = 0; i < nPhases(); i++) {
if (m_moles[i] > 0.0) {
sum += m_phase[i]->cp_mole() * m_moles[i];
}
@ -437,7 +432,7 @@ void MultiPhase::getMoles(doublereal* molNum) const
copy(m_moleFractions.begin(), m_moleFractions.end(), molNum);
size_t ik;
doublereal* dtmp = molNum;
for (size_t ip = 0; ip < m_np; ip++) {
for (size_t ip = 0; ip < nPhases(); ip++) {
doublereal phasemoles = m_moles[ip];
ThermoPhase* p = m_phase[ip];
size_t nsp = p->nSpecies();
@ -455,7 +450,7 @@ void MultiPhase::setMoles(const doublereal* n)
size_t ip, loc = 0;
size_t ik, k = 0, nsp;
doublereal phasemoles;
for (ip = 0; ip < m_np; ip++) {
for (ip = 0; ip < nPhases(); ip++) {
ThermoPhase* p = m_phase[ip];
nsp = p->nSpecies();
phasemoles = 0.0;
@ -523,7 +518,7 @@ void MultiPhase::calcElemAbundances() const
for (eGlobal = 0; eGlobal < m_nel; eGlobal++) {
m_elemAbundances[eGlobal] = 0.0;
}
for (size_t ip = 0; ip < m_np; ip++) {
for (size_t ip = 0; ip < nPhases(); ip++) {
ThermoPhase* p = m_phase[ip];
size_t nspPhase = p->nSpecies();
doublereal phasemoles = m_moles[ip];
@ -542,7 +537,7 @@ doublereal MultiPhase::volume() const
{
int i;
doublereal sum = 0;
for (i = 0; i < int(m_np); i++) {
for (i = 0; i < int(nPhases()); i++) {
double vol = 1.0/m_phase[i]->molarDensity();
sum += m_moles[i] * vol;
}
@ -858,7 +853,7 @@ std::string MultiPhase::phaseName(const size_t iph) const
int MultiPhase::phaseIndex(const std::string& pName) const
{
std::string tmp;
for (int iph = 0; iph < (int) m_np; iph++) {
for (int iph = 0; iph < (int) nPhases(); iph++) {
const ThermoPhase* tptr = m_phase[iph];
tmp = tptr->id();
if (tmp == pName) {
@ -896,7 +891,7 @@ bool MultiPhase::tempOK(const size_t p) const
void MultiPhase::uploadMoleFractionsFromPhases()
{
size_t ip, loc = 0;
for (ip = 0; ip < m_np; ip++) {
for (ip = 0; ip < nPhases(); ip++) {
ThermoPhase* p = m_phase[ip];
p->getMoleFractions(&m_moleFractions[loc]);
loc += p->nSpecies();
@ -907,7 +902,7 @@ void MultiPhase::uploadMoleFractionsFromPhases()
void MultiPhase::updatePhases() const
{
size_t p, nsp, loc = 0;
for (p = 0; p < m_np; p++) {
for (p = 0; p < nPhases(); p++) {
nsp = m_phase[p]->nSpecies();
m_phase[p]->setState_TPX(m_temp, m_press, &m_moleFractions[loc]);
loc += nsp;

8
src/equil/MultiPhaseEquil.cpp
View file

@ -17,7 +17,6 @@ MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_
// store some mixture parameters locally
m_nel_mix = mix->nElements();
m_nsp_mix = mix->nSpecies();
m_np = mix->nPhases();
m_press = mix->pressure();
m_temp = mix->temperature();
@ -580,7 +579,7 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi)
// sum over solution phases
doublereal sum = 0.0, psum;
for (ip = 0; ip < m_np; ip++) {
for (ip = 0; ip < m_mix->nPhases(); ip++) {
ThermoPhase& p = m_mix->phase(ip);
if (p.nSpecies() > 1) {
psum = 0.0;
@ -691,7 +690,6 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
size_t nSpecies;
double vol = 0.0;
string sName;
size_t nphase = m_np;
FILE* FP = fopen(reportFile.c_str(), "w");
if (!FP) {
throw CanteraError("MultiPhaseEquil::reportCSV", "Failure to open file");
@ -708,7 +706,7 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
vector_fp molalities;
vol = 0.0;
for (size_t iphase = 0; iphase < nphase; iphase++) {
for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) {
istart = m_mix->speciesIndex(0, iphase);
ThermoPhase& tref = m_mix->phase(iphase);
nSpecies = tref.nSpecies();
@ -730,7 +728,7 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile)
fprintf(FP,"Pressure = %11.5g Pascal\n", pres);
fprintf(FP,"Total Volume = %11.5g m**3\n", vol);
for (size_t iphase = 0; iphase < nphase; iphase++) {
for (size_t iphase = 0; iphase < m_mix->nPhases(); iphase++) {
istart = m_mix->speciesIndex(0, iphase);
ThermoPhase& tref = m_mix->phase(iphase);
ThermoPhase* tp = &tref;

25
src/kinetics/GasKinetics.cpp
View file

@ -12,7 +12,6 @@ namespace Cantera
{
GasKinetics::GasKinetics(thermo_t* thermo) :
BulkKinetics(thermo),
m_nfall(0),
m_logp_ref(0.0),
m_logc_ref(0.0),
m_logStandConc(0.0),
@ -139,7 +138,7 @@ void GasKinetics::processFalloffReactions()
// use m_ropr for temporary storage of reduced pressure
vector_fp& pr = m_ropr;
for (size_t i = 0; i < m_nfall; i++) {
for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) {
pr[i] = concm_falloff_values[i] * m_rfn_low[i] / (m_rfn_high[i] + SmallNumber);
AssertFinite(pr[i], "GasKinetics::processFalloffReactions",
"pr[{}] is not finite.", i);
@ -147,7 +146,7 @@ void GasKinetics::processFalloffReactions()
m_falloffn.pr_to_falloff(pr.data(), falloff_work.data());
for (size_t i = 0; i < m_nfall; i++) {
for (size_t i = 0; i < m_falloff_low_rates.nReactions(); i++) {
if (reactionType(m_fallindx[i]) == FALLOFF_RXN) {
pr[i] *= m_rfn_high[i];
} else { // CHEMACT_RXN
@ -155,7 +154,7 @@ void GasKinetics::processFalloffReactions()
}
}
scatter_copy(pr.begin(), pr.begin() + m_nfall,
scatter_copy(pr.begin(), pr.begin() + m_falloff_low_rates.nReactions(),
m_ropf.begin(), m_fallindx.begin());
}
@ -175,7 +174,7 @@ void GasKinetics::updateROP()
m_3b_concm.multiply(m_ropf.data(), concm_3b_values.data());
}
if (m_nfall) {
if (m_falloff_high_rates.nReactions()) {
processFalloffReactions();
}
@ -223,7 +222,7 @@ void GasKinetics::getFwdRateConstants(doublereal* kfwd)
m_3b_concm.multiply(m_ropf.data(), concm_3b_values.data());
}
if (m_nfall) {
if (m_falloff_high_rates.nReactions()) {
processFalloffReactions();
}
@ -271,14 +270,15 @@ void GasKinetics::addFalloffReaction(FalloffReaction& r)
{
// install high and low rate coeff calculators and extend the high and low
// rate coeff value vectors
m_falloff_high_rates.install(m_nfall, r.high_rate);
size_t nfall = m_falloff_high_rates.nReactions();
m_falloff_high_rates.install(nfall, r.high_rate);
m_rfn_high.push_back(0.0);
m_falloff_low_rates.install(m_nfall, r.low_rate);
m_falloff_low_rates.install(nfall, r.low_rate);
m_rfn_low.push_back(0.0);
// add this reaction number to the list of falloff reactions
m_fallindx.push_back(nReactions()-1);
m_rfallindx[nReactions()-1] = m_nfall;
m_rfallindx[nReactions()-1] = nfall;
// install the enhanced third-body concentration calculator
map<size_t, double> efficiencies;
@ -292,14 +292,11 @@ void GasKinetics::addFalloffReaction(FalloffReaction& r)
"' while adding reaction '" + r.equation() + "'");
}
}
m_falloff_concm.install(m_nfall, efficiencies,
m_falloff_concm.install(nfall, efficiencies,
r.third_body.default_efficiency);
// install the falloff function calculator for this reaction
m_falloffn.install(m_nfall, r.reaction_type, r.falloff);
// increment the falloff reaction counter
++m_nfall;
m_falloffn.install(nfall, r.reaction_type, r.falloff);
}
void GasKinetics::addThreeBodyReaction(ThreeBodyReaction& r)

25
src/kinetics/ImplicitSurfChem.cpp
View file

@ -16,9 +16,7 @@ namespace Cantera
{
ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
m_nsurf(0),
m_nv(0),
m_numBulkPhases(0),
m_numTotalBulkSpecies(0),
m_numTotalSpecies(0),
m_atol(1.e-14),
@ -30,12 +28,11 @@ ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
m_commonTempPressForPhases(true),
m_ioFlag(0)
{
m_nsurf = k.size();
size_t ns, nsp;
size_t nt, ntmax = 0;
size_t kinSpIndex = 0;
// Loop over the number of surface kinetics objects
for (size_t n = 0; n < m_nsurf; n++) {
for (size_t n = 0; n < k.size(); n++) {
InterfaceKinetics* kinPtr = k[n];
m_vecKinPtrs.push_back(kinPtr);
ns = k[n]->surfacePhaseIndex();
@ -60,9 +57,7 @@ ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> k) :
ThermoPhase* thPtr = & kinPtr->thermo(ip);
if ((imatch = checkMatch(m_bulkPhases, thPtr)) == npos) {
m_bulkPhases.push_back(thPtr);
m_numBulkPhases++;
nsp = thPtr->nSpecies();
m_nspBulkPhases.push_back(nsp);
m_numTotalBulkSpecies += nsp;
imatch = m_bulkPhases.size() - 1;
}
@ -110,7 +105,7 @@ void ImplicitSurfChem::getInitialConditions(doublereal t0, size_t lenc,
void ImplicitSurfChem::getState(doublereal* c)
{
size_t loc = 0;
for (size_t n = 0; n < m_nsurf; n++) {
for (size_t n = 0; n < m_surf.size(); n++) {
m_surf[n]->getCoverages(c + loc);
loc += m_nsp[n];
}
@ -139,7 +134,7 @@ void ImplicitSurfChem::integrate0(doublereal t0, doublereal t1)
void ImplicitSurfChem::updateState(doublereal* c)
{
size_t loc = 0;
for (size_t n = 0; n < m_nsurf; n++) {
for (size_t n = 0; n < m_surf.size(); n++) {
m_surf[n]->setCoverages(c + loc);
loc += m_nsp[n];
}
@ -151,7 +146,7 @@ void ImplicitSurfChem::eval(doublereal time, doublereal* y,
updateState(y); // synchronize the surface state(s) with y
doublereal rs0, sum;
size_t loc = 0, kstart;
for (size_t n = 0; n < m_nsurf; n++) {
for (size_t n = 0; n < m_surf.size(); n++) {
rs0 = 1.0/m_surf[n]->siteDensity();
m_vecKinPtrs[n]->getNetProductionRates(m_work.data());
kstart = m_vecKinPtrs[n]->kineticsSpeciesIndex(0,m_surfindex[n]);
@ -244,13 +239,13 @@ void ImplicitSurfChem::solvePseudoSteadyStateProblem(int ifuncOverride,
void ImplicitSurfChem::getConcSpecies(doublereal* const vecConcSpecies) const
{
size_t kstart;
for (size_t ip = 0; ip < m_nsurf; ip++) {
for (size_t ip = 0; ip < m_surf.size(); ip++) {
ThermoPhase* TP_ptr = m_surf[ip];
kstart = m_specStartIndex[ip];
TP_ptr->getConcentrations(vecConcSpecies + kstart);
}
kstart = m_nv;
for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
ThermoPhase* TP_ptr = m_bulkPhases[ip];
TP_ptr->getConcentrations(vecConcSpecies + kstart);
kstart += TP_ptr->nSpecies();
@ -260,13 +255,13 @@ void ImplicitSurfChem::getConcSpecies(doublereal* const vecConcSpecies) const
void ImplicitSurfChem::setConcSpecies(const doublereal* const vecConcSpecies)
{
size_t kstart;
for (size_t ip = 0; ip < m_nsurf; ip++) {
for (size_t ip = 0; ip < m_surf.size(); ip++) {
ThermoPhase* TP_ptr = m_surf[ip];
kstart = m_specStartIndex[ip];
TP_ptr->setConcentrations(vecConcSpecies + kstart);
}
kstart = m_nv;
for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
ThermoPhase* TP_ptr = m_bulkPhases[ip];
TP_ptr->setConcentrations(vecConcSpecies + kstart);
kstart += TP_ptr->nSpecies();
@ -275,11 +270,11 @@ void ImplicitSurfChem::setConcSpecies(const doublereal* const vecConcSpecies)
void ImplicitSurfChem::setCommonState_TP(doublereal TKelvin, doublereal PresPa)
{
for (size_t ip = 0; ip < m_nsurf; ip++) {
for (size_t ip = 0; ip < m_surf.size(); ip++) {
ThermoPhase* TP_ptr = m_surf[ip];
TP_ptr->setState_TP(TKelvin, PresPa);
}
for (size_t ip = 0; ip < m_numBulkPhases; ip++) {
for (size_t ip = 0; ip < m_bulkPhases.size(); ip++) {
ThermoPhase* TP_ptr = m_bulkPhases[ip];
TP_ptr->setState_TP(TKelvin, PresPa);
}

16
src/kinetics/InterfaceKinetics.cpp
View file

@ -19,8 +19,6 @@ namespace Cantera
InterfaceKinetics::InterfaceKinetics(thermo_t* thermo) :
m_redo_rates(false),
m_nirrev(0),
m_nrev(0),
m_surf(0),
m_integrator(0),
m_logp0(0.0),
@ -65,8 +63,6 @@ InterfaceKinetics& InterfaceKinetics::operator=(const InterfaceKinetics& right)
m_rates = right.m_rates;
m_redo_rates = right.m_redo_rates;
m_irrev = right.m_irrev;
m_nirrev = right.m_nirrev;
m_nrev = right.m_nrev;
m_conc = right.m_conc;
m_actConc = right.m_actConc;
m_mu0 = right.m_mu0;
@ -198,7 +194,7 @@ void InterfaceKinetics::updateKc()
{
fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
if (m_nrev > 0) {
if (m_revindex.size() > 0) {
/*
* Get the vector of standard state electrochemical potentials for
* species in the Interfacial kinetics object and store it in m_mu0[]
@ -210,7 +206,7 @@ void InterfaceKinetics::updateKc()
// compute Delta mu^0 for all reversible reactions
getRevReactionDelta(m_mu0_Kc.data(), m_rkcn.data());
for (size_t i = 0; i < m_nrev; i++) {
for (size_t i = 0; i < m_revindex.size(); i++) {
size_t irxn = m_revindex[i];
if (irxn == npos || irxn >= nReactions()) {
throw CanteraError("InterfaceKinetics", "illegal value: irxn = {}", irxn);
@ -218,7 +214,7 @@ void InterfaceKinetics::updateKc()
// WARNING this may overflow HKM
m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt);
}
for (size_t i = 0; i != m_nirrev; ++i) {
for (size_t i = 0; i != m_irrev.size(); ++i) {
m_rkcn[ m_irrev[i] ] = 0.0;
}
}
@ -250,7 +246,7 @@ void InterfaceKinetics::checkPartialEquil()
vector_fp dmu(nTotalSpecies(), 0.0);
vector_fp rmu(std::max<size_t>(nReactions(), 1), 0.0);
if (m_nrev > 0) {
if (m_revindex.size() > 0) {
cout << "T = " << thermo(0).temperature() << " " << thermo(0).RT() << endl;
size_t nsp, ik=0;
doublereal delta;
@ -267,7 +263,7 @@ void InterfaceKinetics::checkPartialEquil()
// compute Delta mu^ for all reversible reactions
getRevReactionDelta(dmu.data(), rmu.data());
updateROP();
for (size_t i = 0; i < m_nrev; i++) {
for (size_t i = 0; i < m_revindex.size(); i++) {
size_t irxn = m_revindex[i];
writelog("Reaction {} {}\n",
reactionString(irxn), rmu[irxn]/thermo(0).RT());
@ -695,10 +691,8 @@ bool InterfaceKinetics::addReaction(shared_ptr<Reaction> r_base)
if (r.reversible) {
m_revindex.push_back(i);
m_nrev++;
} else {
m_irrev.push_back(i);
m_nirrev++;
}
m_rxnPhaseIsReactant.emplace_back(nPhases(), false);

6
src/kinetics/RxnRates.cpp
View file

@ -32,8 +32,6 @@ SurfaceArrhenius::SurfaceArrhenius()
, m_acov(0.0)
, m_ecov(0.0)
, m_mcov(0.0)
, m_ncov(0)
, m_nmcov(0)
{
}
@ -44,22 +42,18 @@ SurfaceArrhenius::SurfaceArrhenius(double A, double b, double Ta)
, m_acov(0.0)
, m_ecov(0.0)
, m_mcov(0.0)
, m_ncov(0)
, m_nmcov(0)
{
}
void SurfaceArrhenius::addCoverageDependence(size_t k, doublereal a,
doublereal m, doublereal e)
{
m_ncov++;
m_sp.push_back(k);
m_ac.push_back(a);
m_ec.push_back(e);
if (m != 0.0) {
m_msp.push_back(k);
m_mc.push_back(m);
m_nmcov++;
}
}

13
src/oneD/OneDim.cpp
View file

@ -15,7 +15,7 @@ namespace Cantera
OneDim::OneDim()
: m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5),
m_rdt(0.0), m_jac_ok(false),
m_nd(0), m_bw(0), m_size(0),
m_bw(0), m_size(0),
m_init(false), m_pts(0), m_solve_time(0.0),
m_ss_jac_age(10), m_ts_jac_age(20),
m_interrupt(0), m_nevals(0), m_evaltime(0.0)
@ -26,7 +26,7 @@ OneDim::OneDim()
OneDim::OneDim(vector<Domain1D*> domains) :
m_tmin(1.0e-16), m_tmax(10.0), m_tfactor(0.5),
m_rdt(0.0), m_jac_ok(false),
m_nd(0), m_bw(0), m_size(0),
m_bw(0), m_size(0),
m_init(false), m_solve_time(0.0),
m_ss_jac_age(10), m_ts_jac_age(20),
m_interrupt(0), m_nevals(0), m_evaltime(0.0)
@ -46,7 +46,7 @@ OneDim::~OneDim()
size_t OneDim::domainIndex(const std::string& name)
{
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < m_dom.size(); n++) {
if (domain(n).id() == name) {
return n;
}
@ -72,8 +72,7 @@ void OneDim::addDomain(Domain1D* d)
// add it also to the global domain list, and set its container and position
m_dom.push_back(d);
d->setContainer(this, m_nd);
m_nd++;
d->setContainer(this, m_dom.size()-1);
resize();
}
@ -140,7 +139,7 @@ void OneDim::resize()
// save the statistics for the last grid
saveStats();
m_pts = 0;
for (size_t i = 0; i < m_nd; i++) {
for (size_t i = 0; i < nDomains(); i++) {
Domain1D* d = m_dom[i];
size_t np = d->nPoints();
@ -181,7 +180,7 @@ void OneDim::resize()
m_jac.reset(new MultiJac(*this));
m_jac_ok = false;
for (size_t i = 0; i < m_nd; i++) {
for (size_t i = 0; i < nDomains(); i++) {
m_dom[i]->setJac(m_jac.get());
}
}

34
src/oneD/Sim1D.cpp
View file

@ -22,7 +22,7 @@ Sim1D::Sim1D(vector<Domain1D*>& domains) :
// domain-specific initialization of the solution vector.
m_x.resize(size(), 0.0);
m_xnew.resize(size(), 0.0);
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
domain(n)._getInitialSoln(&m_x[start(n)]);
}
@ -33,7 +33,7 @@ Sim1D::Sim1D(vector<Domain1D*>& domains) :
void Sim1D::setInitialGuess(const std::string& component, vector_fp& locs, vector_fp& vals)
{
for (size_t dom=0; dom<m_nd; dom++) {
for (size_t dom=0; dom<nDomains(); dom++) {
Domain1D& d = domain(dom);
size_t ncomp = d.nComponents();
for (size_t comp=0; comp<ncomp; comp++) {
@ -116,13 +116,13 @@ void Sim1D::restore(const std::string& fname, const std::string& id,
}
vector<XML_Node*> xd = f->getChildren("domain");
if (xd.size() != m_nd) {
if (xd.size() != nDomains()) {
throw CanteraError("Sim1D::restore", "Solution does not contain the "
" correct number of domains. Found {} expected {}.\n",
xd.size(), m_nd);
xd.size(), nDomains());
}
size_t sz = 0;
for (size_t m = 0; m < m_nd; m++) {
for (size_t m = 0; m < nDomains(); m++) {
if (loglevel > 0 && xd[m]->attrib("id") != domain(m).id()) {
writelog("Warning: domain names do not match: '" +
(*xd[m])["id"] + + "' and '" + domain(m).id() + "'\n");
@ -131,7 +131,7 @@ void Sim1D::restore(const std::string& fname, const std::string& id,
}
m_x.resize(sz);
m_xnew.resize(sz);
for (size_t m = 0; m < m_nd; m++) {
for (size_t m = 0; m < nDomains(); m++) {
domain(m).restore(*xd[m], &m_x[domain(m).loc()], loglevel);
}
resize();
@ -149,7 +149,7 @@ void Sim1D::setFlatProfile(size_t dom, size_t comp, doublereal v)
void Sim1D::showSolution(ostream& s)
{
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
if (domain(n).domainType() != cEmptyType) {
domain(n).showSolution_s(s, &m_x[start(n)]);
}
@ -158,7 +158,7 @@ void Sim1D::showSolution(ostream& s)
void Sim1D::showSolution()
{
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
if (domain(n).domainType() != cEmptyType) {
writelog("\n\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> "+domain(n).id()
+" <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n");
@ -169,14 +169,14 @@ void Sim1D::showSolution()
void Sim1D::getInitialSoln()
{
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
domain(n)._getInitialSoln(&m_x[start(n)]);
}
}
void Sim1D::finalize()
{
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
domain(n)._finalize(&m_x[start(n)]);
}
}
@ -319,7 +319,7 @@ int Sim1D::refine(int loglevel)
doublereal xmid, zmid;
std::vector<size_t> dsize;
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
Domain1D& d = domain(n);
Refiner& r = d.refiner();
@ -377,7 +377,7 @@ int Sim1D::refine(int loglevel)
// Now update each domain with the new grid.
size_t gridstart = 0, gridsize;
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
Domain1D& d = domain(n);
gridsize = dsize[n];
d.setupGrid(gridsize, &znew[gridstart]);
@ -406,7 +406,7 @@ int Sim1D::setFixedTemperature(doublereal t)
size_t m1 = 0;
std::vector<size_t> dsize;
for (n = 0; n < m_nd; n++) {
for (n = 0; n < nDomains(); n++) {
bool addnewpt=false;
Domain1D& d = domain(n);
size_t comp = d.nComponents();
@ -469,7 +469,7 @@ int Sim1D::setFixedTemperature(doublereal t)
// been constructed, but the domains themselves have not yet been modified.
// Now update each domain with the new grid.
size_t gridstart = 0, gridsize;
for (n = 0; n < m_nd; n++) {
for (n = 0; n < nDomains(); n++) {
Domain1D& d = domain(n);
gridsize = dsize[n];
d.setupGrid(gridsize, &znew[gridstart]);
@ -493,7 +493,7 @@ void Sim1D::setRefineCriteria(int dom, doublereal ratio,
Refiner& r = domain(dom).refiner();
r.setCriteria(ratio, slope, curve, prune);
} else {
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
Refiner& r = domain(n).refiner();
r.setCriteria(ratio, slope, curve, prune);
}
@ -506,7 +506,7 @@ void Sim1D::setGridMin(int dom, double gridmin)
Refiner& r = domain(dom).refiner();
r.setGridMin(gridmin);
} else {
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
Refiner& r = domain(n).refiner();
r.setGridMin(gridmin);
}
@ -519,7 +519,7 @@ void Sim1D::setMaxGridPoints(int dom, int npoints)
Refiner& r = domain(dom).refiner();
r.setMaxPoints(npoints);
} else {
for (size_t n = 0; n < m_nd; n++) {
for (size_t n = 0; n < nDomains(); n++) {
Refiner& r = domain(n).refiner();
r.setMaxPoints(npoints);
}

76
src/thermo/LatticeSolidPhase.cpp
View file

@ -20,15 +20,13 @@ namespace Cantera
{
LatticeSolidPhase::LatticeSolidPhase() :
m_press(-1.0),
m_molar_density(0.0),
m_nlattice(0)
m_molar_density(0.0)
{
}
LatticeSolidPhase::LatticeSolidPhase(const LatticeSolidPhase& right) :
m_press(-1.0),
m_molar_density(0.0),
m_nlattice(0)
m_molar_density(0.0)
{
*this = right;
}
@ -40,7 +38,6 @@ LatticeSolidPhase& LatticeSolidPhase::operator=(const LatticeSolidPhase& right)
m_tlast = right.m_tlast;
m_press = right.m_press;
m_molar_density = right.m_molar_density;
m_nlattice = right.m_nlattice;
deepStdVectorPointerCopy<LatticePhase>(right.m_lattice, m_lattice);
m_x = right.m_x;
theta_ = right.theta_;
@ -52,7 +49,7 @@ LatticeSolidPhase& LatticeSolidPhase::operator=(const LatticeSolidPhase& right)
LatticeSolidPhase::~LatticeSolidPhase()
{
// We own the sublattices. So we have to delete the sublattices
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
delete m_lattice[n];
m_lattice[n] = 0;
}
@ -66,14 +63,14 @@ ThermoPhase* LatticeSolidPhase::duplMyselfAsThermoPhase() const
doublereal LatticeSolidPhase::minTemp(size_t k) const
{
if (k != npos) {
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
if (lkstart_[n+1] < k) {
return m_lattice[n]->minTemp(k-lkstart_[n]);
}
}
}
doublereal mm = 1.0E300;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
double ml = m_lattice[n]->minTemp();
mm = std::min(mm, ml);
}
@ -83,14 +80,14 @@ doublereal LatticeSolidPhase::minTemp(size_t k) const
doublereal LatticeSolidPhase::maxTemp(size_t k) const
{
if (k != npos) {
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
if (lkstart_[n+1] < k) {
return (m_lattice[n])->maxTemp(k - lkstart_[n]);
}
}
}
doublereal mm = -1.0E300;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
double ml = m_lattice[n]->maxTemp();
mm = std::max(mm, ml);
}
@ -106,7 +103,7 @@ doublereal LatticeSolidPhase::enthalpy_mole() const
{
_updateThermo();
doublereal sum = 0.0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
sum += theta_[n] * m_lattice[n]->enthalpy_mole();
}
return sum;
@ -116,7 +113,7 @@ doublereal LatticeSolidPhase::intEnergy_mole() const
{
_updateThermo();
doublereal sum = 0.0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
sum += theta_[n] * m_lattice[n]->intEnergy_mole();
}
return sum;
@ -126,7 +123,7 @@ doublereal LatticeSolidPhase::entropy_mole() const
{
_updateThermo();
doublereal sum = 0.0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
sum += theta_[n] * m_lattice[n]->entropy_mole();
}
return sum;
@ -136,7 +133,7 @@ doublereal LatticeSolidPhase::gibbs_mole() const
{
_updateThermo();
doublereal sum = 0.0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
sum += theta_[n] * m_lattice[n]->gibbs_mole();
}
return sum;
@ -146,7 +143,7 @@ doublereal LatticeSolidPhase::cp_mole() const
{
_updateThermo();
doublereal sum = 0.0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
sum += theta_[n] * m_lattice[n]->cp_mole();
}
return sum;
@ -156,7 +153,7 @@ void LatticeSolidPhase::getActivityConcentrations(doublereal* c) const
{
_updateThermo();
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
m_lattice[n]->getMoleFractions(c+strt);
strt += m_lattice[n]->nSpecies();
}
@ -182,7 +179,7 @@ doublereal LatticeSolidPhase::logStandardConc(size_t k) const
void LatticeSolidPhase::setPressure(doublereal p)
{
m_press = p;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
m_lattice[n]->setPressure(m_press);
}
calcDensity();
@ -191,7 +188,7 @@ void LatticeSolidPhase::setPressure(doublereal p)
doublereal LatticeSolidPhase::calcDensity()
{
double sum = 0.0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
sum += theta_[n] * m_lattice[n]->density();
}
Phase::setDensity(sum);
@ -201,13 +198,13 @@ doublereal LatticeSolidPhase::calcDensity()
void LatticeSolidPhase::setMoleFractions(const doublereal* const x)
{
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nsp = m_lattice[n]->nSpecies();
m_lattice[n]->setMoleFractions(x + strt);
strt += nsp;
}
for (size_t k = 0; k < strt; k++) {
m_x[k] = x[k] / m_nlattice;
m_x[k] = x[k] / m_lattice.size();
}
Phase::setMoleFractions(m_x.data());
calcDensity();
@ -218,7 +215,7 @@ void LatticeSolidPhase::getMoleFractions(doublereal* const x) const
size_t strt = 0;
// the ifdef block should be the way we calculate this.!!!!!
Phase::getMoleFractions(x);
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nsp = m_lattice[n]->nSpecies();
double sum = 0.0;
for (size_t k = 0; k < nsp; k++) {
@ -245,7 +242,7 @@ void LatticeSolidPhase::getChemPotentials(doublereal* mu) const
{
_updateThermo();
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getChemPotentials(mu+strt);
strt += nlsp;
@ -256,7 +253,7 @@ void LatticeSolidPhase::getPartialMolarEnthalpies(doublereal* hbar) const
{
_updateThermo();
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getPartialMolarEnthalpies(hbar + strt);
strt += nlsp;
@ -267,7 +264,7 @@ void LatticeSolidPhase::getPartialMolarEntropies(doublereal* sbar) const
{
_updateThermo();
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getPartialMolarEntropies(sbar + strt);
strt += nlsp;
@ -278,7 +275,7 @@ void LatticeSolidPhase::getPartialMolarCp(doublereal* cpbar) const
{
_updateThermo();
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getPartialMolarCp(cpbar + strt);
strt += nlsp;
@ -289,7 +286,7 @@ void LatticeSolidPhase::getPartialMolarVolumes(doublereal* vbar) const
{
_updateThermo();
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nlsp = m_lattice[n]->nSpecies();
m_lattice[n]->getPartialMolarVolumes(vbar + strt);
strt += nlsp;
@ -300,7 +297,7 @@ void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const
{
_updateThermo();
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
m_lattice[n]->getStandardChemPotentials(mu0+strt);
strt += m_lattice[n]->nSpecies();
}
@ -309,7 +306,7 @@ void LatticeSolidPhase::getStandardChemPotentials(doublereal* mu0) const
void LatticeSolidPhase::getGibbs_RT_ref(doublereal* grt) const
{
_updateThermo();
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
m_lattice[n]->getGibbs_RT_ref(grt + lkstart_[n]);
}
}
@ -330,7 +327,7 @@ void LatticeSolidPhase::installSlavePhases(XML_Node* phaseNode)
XML_Node& la = phaseNode->child("thermo").child("LatticeArray");
std::vector<XML_Node*> lattices = la.getChildren("phase");
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
LatticePhase* lp = m_lattice[n];
vector_fp constArr(lp->nElements());
const vector_fp& aws = lp->atomicWeights();
@ -365,11 +362,11 @@ void LatticeSolidPhase::initThermo()
{
initLengths();
size_t loc = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nsp = m_lattice[n]->nSpecies();
lkstart_[n] = loc;
for (size_t k = 0; k < nsp; k++) {
m_x[loc] =m_lattice[n]->moleFraction(k) / (double) m_nlattice;
m_x[loc] =m_lattice[n]->moleFraction(k) / (double) m_lattice.size();
loc++;
}
lkstart_[n+1] = loc;
@ -380,8 +377,8 @@ void LatticeSolidPhase::initThermo()
void LatticeSolidPhase::initLengths()
{
theta_.resize(m_nlattice,0);
lkstart_.resize(m_nlattice+1);
theta_.resize(m_lattice.size(), 0);
lkstart_.resize(m_lattice.size() + 1);
m_x.resize(m_kk, 0.0);
tmpV_.resize(m_kk, 0.0);
}
@ -392,7 +389,7 @@ void LatticeSolidPhase::_updateThermo() const
if (m_tlast != tnow) {
getMoleFractions(m_x.data());
size_t strt = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
m_lattice[n]->setTemperature(tnow);
m_lattice[n]->setMoleFractions(&m_x[strt]);
m_lattice[n]->setPressure(m_press);
@ -406,7 +403,7 @@ void LatticeSolidPhase::setLatticeMoleFractionsByName(int nn, const std::string&
{
m_lattice[nn]->setMoleFractionsByName(x);
size_t loc = 0;
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
size_t nsp = m_lattice[n]->nSpecies();
double ndens = m_lattice[n]->molarDensity();
for (size_t k = 0; k < nsp; k++) {
@ -422,18 +419,17 @@ void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata)
eosdata._require("model","LatticeSolid");
XML_Node& la = eosdata.child("LatticeArray");
std::vector<XML_Node*> lattices = la.getChildren("phase");
m_nlattice = lattices.size();
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < lattices.size(); n++) {
m_lattice.push_back((LatticePhase*)newPhase(*lattices[n]));
}
std::vector<string> pnam;
std::vector<string> pval;
int np = getPairs(eosdata.child("LatticeStoichiometry"), pnam, pval);
theta_.resize(m_nlattice);
theta_.resize(m_lattice.size());
for (int i = 0; i < np; i++) {
double val = fpValueCheck(pval[i]);
bool found = false;
for (size_t j = 0; j < m_nlattice; j++) {
for (size_t j = 0; j < m_lattice.size(); j++) {
ThermoPhase& tp = *m_lattice[j];
string idj = tp.id();
if (idj == pnam[i]) {
@ -450,7 +446,7 @@ void LatticeSolidPhase::setParametersFromXML(const XML_Node& eosdata)
void LatticeSolidPhase::modifyOneHf298SS(const size_t k, const doublereal Hf298New)
{
for (size_t n = 0; n < m_nlattice; n++) {
for (size_t n = 0; n < m_lattice.size(); n++) {
if (lkstart_[n+1] < k) {
size_t kk = k-lkstart_[n];
SpeciesThermo& l_spthermo = m_lattice[n]->speciesThermo();

29
src/thermo/Nasa9PolyMultiTempRegion.cpp
View file

@ -20,25 +20,22 @@ using namespace std;
namespace Cantera
{
Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion() :
m_numTempRegions(0),
m_currRegion(0)
{
}
Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(vector<Nasa9Poly1*>& regionPts) :
m_numTempRegions(0),
m_currRegion(0)
{
m_numTempRegions = regionPts.size();
// From now on, we own these pointers
for (Nasa9Poly1* region : regionPts) {
m_regionPts.emplace_back(region);
}
m_lowerTempBounds.resize(m_numTempRegions);
m_lowerTempBounds.resize(regionPts.size());
m_lowT = m_regionPts[0]->minTemp();
m_highT = m_regionPts[m_numTempRegions-1]->maxTemp();
m_highT = m_regionPts[m_regionPts.size()-1]->maxTemp();
m_Pref = m_regionPts[0]->refPressure();
for (size_t i = 0; i < m_numTempRegions; i++) {
for (size_t i = 0; i < m_regionPts.size(); i++) {
m_lowerTempBounds[i] = m_regionPts[i]->minTemp();
if (fabs(m_regionPts[i]->refPressure() - m_Pref) > 0.0001) {
throw CanteraError("Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion",
@ -59,12 +56,11 @@ Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(vector<Nasa9Poly1*>& regionPt
Nasa9PolyMultiTempRegion::Nasa9PolyMultiTempRegion(const Nasa9PolyMultiTempRegion& b) :
SpeciesThermoInterpType(b),
m_numTempRegions(b.m_numTempRegions),
m_lowerTempBounds(b.m_lowerTempBounds),
m_currRegion(b.m_currRegion)
{
m_regionPts.resize(m_numTempRegions);
for (size_t i = 0; i < m_numTempRegions; i++) {
m_regionPts.resize(b.m_regionPts.size());
for (size_t i = 0; i < m_regionPts.size(); i++) {
m_regionPts[i].reset(new Nasa9Poly1(*b.m_regionPts[i]));
}
}
@ -74,11 +70,10 @@ Nasa9PolyMultiTempRegion::operator=(const Nasa9PolyMultiTempRegion& b)
{
if (&b != this) {
SpeciesThermoInterpType::operator=(b);
m_numTempRegions = b.m_numTempRegions;
m_lowerTempBounds = b.m_lowerTempBounds;
m_currRegion = b.m_currRegion;
m_regionPts.resize(m_numTempRegions);
for (size_t i = 0; i < m_numTempRegions; i++) {
m_regionPts.resize(b.m_regionPts.size());
for (size_t i = 0; i < m_regionPts.size(); i++) {
m_regionPts[i].reset(new Nasa9Poly1(*b.m_regionPts[i]));
}
}
@ -117,7 +112,7 @@ void Nasa9PolyMultiTempRegion::updateProperties(const doublereal* tt,
doublereal* s_R) const
{
m_currRegion = 0;
for (size_t i = 1; i < m_numTempRegions; i++) {
for (size_t i = 1; i < m_regionPts.size(); i++) {
if (tt[0] < m_lowerTempBounds[i]) {
break;
}
@ -133,7 +128,7 @@ void Nasa9PolyMultiTempRegion::updatePropertiesTemp(const doublereal temp,
{
// Now find the region
m_currRegion = 0;
for (size_t i = 1; i < m_numTempRegions; i++) {
for (size_t i = 1; i < m_regionPts.size(); i++) {
if (temp < m_lowerTempBounds[i]) {
break;
}
@ -154,12 +149,12 @@ void Nasa9PolyMultiTempRegion::reportParameters(size_t& n, int& type,
thigh = m_highT;
pref = m_Pref;
double ctmp[12];
coeffs[0] = double(m_numTempRegions);
coeffs[0] = double(m_regionPts.size());
int index = 1;
size_t n_tmp = 0;
int type_tmp = 0;
double pref_tmp = 0.0;
for (size_t iReg = 0; iReg < m_numTempRegions; iReg++) {
for (size_t iReg = 0; iReg < m_regionPts.size(); iReg++) {
m_regionPts[iReg]->reportParameters(n_tmp, type_tmp,
coeffs[index], coeffs[index+1],
pref_tmp, ctmp);
@ -173,7 +168,7 @@ void Nasa9PolyMultiTempRegion::reportParameters(size_t& n, int& type,
void Nasa9PolyMultiTempRegion::modifyParameters(doublereal* coeffs)
{
int index = 3;
for (size_t iReg = 0; iReg < m_numTempRegions; iReg++) {
for (size_t iReg = 0; iReg < m_regionPts.size(); iReg++) {
m_regionPts[iReg]->modifyParameters(coeffs + index);
index += 11;
}