From 920ff1d89773c6250c3a13795da2e25b64b620ed Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 11 May 2016 16:22:41 -0400 Subject: [PATCH] [Kinetics] Make better use of local variables --- src/kinetics/Falloff.cpp | 11 +- src/kinetics/Group.cpp | 11 +- src/kinetics/ImplicitSurfChem.cpp | 18 +-- src/kinetics/InterfaceKinetics.cpp | 26 ++-- src/kinetics/Kinetics.cpp | 20 +-- src/kinetics/ReactionPath.cpp | 225 ++++++++++++----------------- src/kinetics/importKinetics.cpp | 19 +-- src/kinetics/solveSP.cpp | 206 +++++++++++--------------- 8 files changed, 217 insertions(+), 319 deletions(-) diff --git a/src/kinetics/Falloff.cpp b/src/kinetics/Falloff.cpp index b24621ff2..e49537cf1 100644 --- a/src/kinetics/Falloff.cpp +++ b/src/kinetics/Falloff.cpp @@ -53,12 +53,11 @@ void Troe::updateTemp(double T, double* work) const double Troe::F(double pr, const double* work) const { - double lpr,f1,lgf, cc, nn; - lpr = log10(std::max(pr,SmallNumber)); - cc = -0.4 - 0.67 * (*work); - nn = 0.75 - 1.27 * (*work); - f1 = (lpr + cc)/ (nn - 0.14 * (lpr + cc)); - lgf = (*work) / (1.0 + f1 * f1); + double lpr = log10(std::max(pr,SmallNumber)); + double cc = -0.4 - 0.67 * (*work); + double nn = 0.75 - 1.27 * (*work); + double f1 = (lpr + cc)/ (nn - 0.14 * (lpr + cc)); + double lgf = (*work) / (1.0 + f1 * f1); return pow(10.0, lgf); } diff --git a/src/kinetics/Group.cpp b/src/kinetics/Group.cpp index a5f8bcf86..b66c81d61 100644 --- a/src/kinetics/Group.cpp +++ b/src/kinetics/Group.cpp @@ -13,8 +13,6 @@ namespace Cantera void Group::validate() { - size_t n = m_comp.size(); - // if already checked and not valid, return if (m_sign == -999) { return; @@ -22,7 +20,7 @@ void Group::validate() m_sign = 0; bool ok = true; - for (size_t m = 0; m < n; m++) { + for (size_t m = 0; m < m_comp.size(); m++) { if (m_comp[m] != 0) { if (m_sign == 0) { m_sign = m_comp[m]/abs(m_comp[m]); @@ -34,7 +32,6 @@ void Group::validate() } if (!ok) { m_sign = -999; - m_comp.resize(n,0); } } @@ -42,11 +39,9 @@ std::ostream& Group::fmt(std::ostream& s, const std::vector& esymbols) const { s << "("; - int nm; bool first = true; - size_t n = m_comp.size(); - for (size_t m = 0; m < n; m++) { - nm = m_comp[m]; + for (size_t m = 0; m < m_comp.size(); m++) { + int nm = m_comp[m]; if (nm != 0) { if (!first) { s << "-"; diff --git a/src/kinetics/ImplicitSurfChem.cpp b/src/kinetics/ImplicitSurfChem.cpp index 3559d0e42..67b6c8dfa 100644 --- a/src/kinetics/ImplicitSurfChem.cpp +++ b/src/kinetics/ImplicitSurfChem.cpp @@ -28,24 +28,23 @@ ImplicitSurfChem::ImplicitSurfChem(vector k) : m_commonTempPressForPhases(true), m_ioFlag(0) { - size_t ns, nsp; - size_t nt, ntmax = 0; + size_t ntmax = 0; size_t kinSpIndex = 0; // Loop over the number of surface kinetics objects for (size_t n = 0; n < k.size(); n++) { InterfaceKinetics* kinPtr = k[n]; m_vecKinPtrs.push_back(kinPtr); - ns = k[n]->surfacePhaseIndex(); + size_t ns = k[n]->surfacePhaseIndex(); if (ns == npos) { throw CanteraError("ImplicitSurfChem", "kinetics manager contains no surface phase"); } m_surfindex.push_back(ns); m_surf.push_back((SurfPhase*)&k[n]->thermo(ns)); - nsp = m_surf.back()->nSpecies(); + size_t nsp = m_surf.back()->nSpecies(); m_nsp.push_back(nsp); m_nv += m_nsp.back(); - nt = k[n]->nTotalSpecies(); + size_t nt = k[n]->nTotalSpecies(); ntmax = std::max(nt, ntmax); m_specStartIndex.push_back(kinSpIndex); kinSpIndex += nsp; @@ -144,13 +143,12 @@ void ImplicitSurfChem::eval(doublereal time, doublereal* y, doublereal* ydot, doublereal* p) { updateState(y); // synchronize the surface state(s) with y - doublereal rs0, sum; - size_t loc = 0, kstart; + size_t loc = 0; for (size_t n = 0; n < m_surf.size(); n++) { - rs0 = 1.0/m_surf[n]->siteDensity(); + double rs0 = 1.0/m_surf[n]->siteDensity(); m_vecKinPtrs[n]->getNetProductionRates(m_work.data()); - kstart = m_vecKinPtrs[n]->kineticsSpeciesIndex(0,m_surfindex[n]); - sum = 0.0; + size_t kstart = m_vecKinPtrs[n]->kineticsSpeciesIndex(0,m_surfindex[n]); + double sum = 0.0; for (size_t k = 1; k < m_nsp[n]; k++) { ydot[k + loc] = m_work[kstart + k] * rs0 * m_surf[n]->size(k); sum -= ydot[k]; diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 58cd21d47..661307ae9 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -226,12 +226,10 @@ void InterfaceKinetics::updateMu0() _update_rates_phi(); updateExchangeCurrentQuantities(); - size_t nsp, ik = 0; - size_t np = nPhases(); - for (size_t n = 0; n < np; n++) { + size_t ik = 0; + for (size_t n = 0; n < nPhases(); n++) { thermo(n).getStandardChemPotentials(m_mu0.data() + m_start[n]); - nsp = thermo(n).nSpecies(); - for (size_t k = 0; k < nsp; k++) { + for (size_t k = 0; k < thermo(n).nSpecies(); k++) { m_mu0_Kc[ik] = m_mu0[ik] + Faraday * m_phi[n] * thermo(n).charge(k); m_mu0_Kc[ik] -= thermo(0).RT() * thermo(n).logStandardConc(k); ik++; @@ -340,14 +338,13 @@ void InterfaceKinetics::applyVoltageKfwdCorrection(doublereal* const kf) // NOTE, there is some discussion about this point. Should we decrease the // activation energy below zero? I don't think this has been decided in any // definitive way. The treatment below is numerically more stable, however. - doublereal eamod; for (size_t i = 0; i < m_beta.size(); i++) { size_t irxn = m_ctrxn[i]; // If we calculate the BV form directly, we don't add the voltage // correction to the forward reaction rate constants. if (m_ctrxn_BVform[i] == 0) { - eamod = m_beta[i] * deltaElectricEnergy_[irxn]; + double eamod = m_beta[i] * deltaElectricEnergy_[irxn]; if (eamod != 0.0) { kf[irxn] *= exp(-eamod/thermo(0).RT()); } @@ -376,8 +373,7 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con if (m_ctrxn_BVform[i] == 0) { // Calculate the term and modify the forward reaction double tmp = exp(- m_beta[i] * m_deltaG0[irxn] / thermo(0).RT()); - double tmp2 = m_ProdStanConcReac[irxn]; - tmp *= 1.0 / tmp2 / Faraday; + tmp *= 1.0 / m_ProdStanConcReac[irxn] / Faraday; kfwd[irxn] *= tmp; } // If BVform is nonzero we don't need to do anything. @@ -392,8 +388,7 @@ void InterfaceKinetics::convertExchangeCurrentDensityFormulation(doublereal* con // constant so that it's in the exchange current density // formulation format double tmp = exp(m_beta[i] * m_deltaG0[irxn] * thermo(0).RT()); - double tmp2 = m_ProdStanConcReac[irxn]; - tmp *= Faraday * tmp2; + tmp *= Faraday * m_ProdStanConcReac[irxn]; kfwd[irxn] *= tmp; } } @@ -550,8 +545,7 @@ void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) void InterfaceKinetics::getDeltaElectrochemPotentials(doublereal* deltaM) { // Get the chemical potentials of the species - size_t np = nPhases(); - for (size_t n = 0; n < np; n++) { + for (size_t n = 0; n < nPhases(); n++) { thermo(n).getElectrochemPotentials(m_grt.data() + m_start[n]); } @@ -873,8 +867,7 @@ doublereal InterfaceKinetics::electrochem_beta(size_t irxn) const void InterfaceKinetics::advanceCoverages(doublereal tstep) { if (m_integrator == 0) { - vector k; - k.push_back(this); + vector k{this}; m_integrator = new ImplicitSurfChem(k); m_integrator->initialize(); } @@ -888,8 +881,7 @@ void InterfaceKinetics::solvePseudoSteadyStateProblem( { // create our own solver object if (m_integrator == 0) { - vector k; - k.push_back(this); + vector k{this}; m_integrator = new ImplicitSurfChem(k); m_integrator->initialize(); } diff --git a/src/kinetics/Kinetics.cpp b/src/kinetics/Kinetics.cpp index fae7da91b..55a8824ae 100644 --- a/src/kinetics/Kinetics.cpp +++ b/src/kinetics/Kinetics.cpp @@ -121,12 +121,11 @@ void Kinetics::checkSpeciesArraySize(size_t kk) const void Kinetics::assignShallowPointers(const std::vector & tpVector) { - size_t ns = tpVector.size(); - if (ns != m_thermo.size()) { + if (tpVector.size() != m_thermo.size()) { throw CanteraError(" Kinetics::assignShallowPointers", " Number of ThermoPhase objects arent't the same"); } - for (size_t i = 0; i < ns; i++) { + for (size_t i = 0; i < tpVector.size(); i++) { ThermoPhase* ntp = tpVector[i]; ThermoPhase* otp = m_thermo[i]; if (ntp->id() != otp->id()) { @@ -370,17 +369,13 @@ size_t Kinetics::kineticsSpeciesIndex(const std::string& nm, thermo_t& Kinetics::speciesPhase(const std::string& nm) { - size_t np = m_thermo.size(); - size_t k; - string id; - for (size_t n = 0; n < np; n++) { - k = thermo(n).speciesIndex(nm); + for (size_t n = 0; n < m_thermo.size(); n++) { + size_t k = thermo(n).speciesIndex(nm); if (k != npos) { return thermo(n); } } throw CanteraError("speciesPhase", "unknown species "+nm); - return thermo(0); } size_t Kinetics::speciesPhaseIndex(size_t k) @@ -391,7 +386,6 @@ size_t Kinetics::speciesPhaseIndex(size_t k) } } throw CanteraError("speciesPhaseIndex", "illegal species index: {}", k); - return npos; } double Kinetics::reactantStoichCoeff(size_t kSpec, size_t irxn) const @@ -576,14 +570,12 @@ bool Kinetics::addReaction(shared_ptr r) vector_fp rstoich, pstoich; for (const auto& sp : r->reactants) { - size_t k = kineticsSpeciesIndex(sp.first); - rk.push_back(k); + rk.push_back(kineticsSpeciesIndex(sp.first)); rstoich.push_back(sp.second); } for (const auto& sp : r->products) { - size_t k = kineticsSpeciesIndex(sp.first); - pk.push_back(k); + pk.push_back(kineticsSpeciesIndex(sp.first)); pstoich.push_back(sp.second); } diff --git a/src/kinetics/ReactionPath.cpp b/src/kinetics/ReactionPath.cpp index d9df10cde..eee0c3077 100644 --- a/src/kinetics/ReactionPath.cpp +++ b/src/kinetics/ReactionPath.cpp @@ -52,14 +52,13 @@ void Path::addReaction(size_t rxnNumber, doublereal value, void Path::writeLabel(ostream& s, doublereal threshold) { - size_t nn = m_label.size(); - if (nn == 0) { + if (m_label.size() == 0) { return; } doublereal v; for (const auto& label : m_label) { v = label.second/m_total; - if (nn == 1) { + if (m_label.size() == 1) { s << label.first << "\\l"; } else if (v > threshold) { s << label.first; @@ -113,16 +112,15 @@ ReactionPathDiagram::~ReactionPathDiagram() vector_int ReactionPathDiagram::reactions() { - size_t i, npaths = nPaths(); - doublereal flmax = 0.0, flxratio; - for (i = 0; i < npaths; i++) { + doublereal flmax = 0.0; + for (size_t i = 0; i < nPaths(); i++) { Path* p = path(i); flmax = std::max(p->flow(), flmax); } m_rxns.clear(); - for (i = 0; i < npaths; i++) { + for (size_t i = 0; i < nPaths(); i++) { for (const auto& rxn : path(i)->reactionMap()) { - flxratio = rxn.second/flmax; + double flxratio = rxn.second/flmax; if (flxratio > threshold) { m_rxns[rxn.first] = 1; } @@ -137,13 +135,10 @@ vector_int ReactionPathDiagram::reactions() void ReactionPathDiagram::add(ReactionPathDiagram& d) { - size_t np = nPaths(); - size_t n, k1, k2; - Path* p = 0; - for (n = 0; n < np; n++) { - p = path(n); - k1 = p->begin()->number; - k2 = p->end()->number; + for (size_t n = 0; n < nPaths(); n++) { + Path* p = path(n); + size_t k1 = p->begin()->number; + size_t k2 = p->end()->number; p->setFlow(p->flow() + d.flow(k1,k2)); } } @@ -151,22 +146,20 @@ void ReactionPathDiagram::add(ReactionPathDiagram& d) void ReactionPathDiagram::findMajorPaths(doublereal athreshold, size_t lda, doublereal* a) { - size_t nn = nNodes(); - size_t n, m, k1, k2; - doublereal fl, netmax = 0.0; - for (n = 0; n < nn; n++) { - for (m = n+1; m < nn; m++) { - k1 = m_speciesNumber[n]; - k2 = m_speciesNumber[m]; - fl = fabs(netFlow(k1,k2)); + double netmax = 0.0; + for (size_t n = 0; n < nNodes(); n++) { + for (size_t m = n+1; m < nNodes(); m++) { + size_t k1 = m_speciesNumber[n]; + size_t k2 = m_speciesNumber[m]; + double fl = fabs(netFlow(k1,k2)); netmax = std::max(fl, netmax); } } - for (n = 0; n < nn; n++) { - for (m = n+1; m < nn; m++) { - k1 = m_speciesNumber[n]; - k2 = m_speciesNumber[m]; - fl = fabs(netFlow(k1,k2)); + for (size_t n = 0; n < nNodes(); n++) { + for (size_t m = n+1; m < nNodes(); m++) { + size_t k1 = m_speciesNumber[n]; + size_t k2 = m_speciesNumber[m]; + double fl = fabs(netFlow(k1,k2)); if (fl > athreshold*netmax) { a[lda*k1 + k2] = 1; } @@ -176,21 +169,18 @@ void ReactionPathDiagram::findMajorPaths(doublereal athreshold, size_t lda, void ReactionPathDiagram::writeData(ostream& s) { - doublereal f1, f2; - size_t nnodes = nNodes(); - size_t i1, i2, k1, k2; s << title << endl; - for (i1 = 0; i1 < nnodes; i1++) { - k1 = m_speciesNumber[i1]; + for (size_t i1 = 0; i1 < nNodes(); i1++) { + size_t k1 = m_speciesNumber[i1]; s << m_nodes[k1]->name << " "; } s << endl; - for (i1 = 0; i1 < nnodes; i1++) { - k1 = m_speciesNumber[i1]; - for (i2 = i1+1; i2 < nnodes; i2++) { - k2 = m_speciesNumber[i2]; - f1 = flow(k1, k2); - f2 = flow(k2, k1); + for (size_t i1 = 0; i1 < nNodes(); i1++) { + size_t k1 = m_speciesNumber[i1]; + for (size_t i2 = i1+1; i2 < nNodes(); i2++) { + size_t k2 = m_speciesNumber[i2]; + double f1 = flow(k1, k2); + double f2 = flow(k2, k1); s << m_nodes[k1]->name << " " << m_nodes[k2]->name << " " << f1 << " " << -f2 << endl; } @@ -199,7 +189,7 @@ void ReactionPathDiagram::writeData(ostream& s) void ReactionPathDiagram::exportToDot(ostream& s) { - doublereal flxratio, flmax = 0.0, lwidth; + doublereal flmax = 0.0; s.precision(3); // a directed graph @@ -220,21 +210,17 @@ void ReactionPathDiagram::exportToDot(ostream& s) s << dot_options << endl; } - Path* p; - size_t kbegin, kend, i1, i2, k1, k2; - doublereal flx; - // draw paths representing net flows if (flow_type == NetFlow) { // if no scale was specified, normalize net flows by the maximum net // flow if (scale <= 0.0) { - for (i1 = 0; i1 < nNodes(); i1++) { - k1 = m_speciesNumber[i1]; + for (size_t i1 = 0; i1 < nNodes(); i1++) { + size_t k1 = m_speciesNumber[i1]; node(k1)->visible = false; - for (i2 = i1+1; i2 < nNodes(); i2++) { - k2 = m_speciesNumber[i2]; - flx = netFlow(k1, k2); + for (size_t i2 = i1+1; i2 < nNodes(); i2++) { + size_t k2 = m_speciesNumber[i2]; + double flx = netFlow(k1, k2); if (flx < 0.0) { flx = -flx; } @@ -247,15 +233,17 @@ void ReactionPathDiagram::exportToDot(ostream& s) flmax = std::max(flmax, 1e-10); // loop over all unique pairs of nodes - for (i1 = 0; i1 < nNodes(); i1++) { - k1 = m_speciesNumber[i1]; - for (i2 = i1+1; i2 < nNodes(); i2++) { - k2 = m_speciesNumber[i2]; - flx = netFlow(k1, k2); + for (size_t i1 = 0; i1 < nNodes(); i1++) { + size_t k1 = m_speciesNumber[i1]; + for (size_t i2 = i1+1; i2 < nNodes(); i2++) { + size_t k2 = m_speciesNumber[i2]; + double flx = netFlow(k1, k2); if (m_local != npos && k1 != m_local && k2 != m_local) { flx = 0.0; } if (flx != 0.0) { + double flxratio; + size_t kbegin, kend; // set beginning and end of the path based on the sign of // the net flow if (flx > 0.0) { @@ -279,7 +267,7 @@ void ReactionPathDiagram::exportToDot(ostream& s) s << "[fontname=\""+m_font+"\", style=\"setlinewidth("; if (arrow_width < 0) { - lwidth = 1.0 - 4.0 + double lwidth = 1.0 - 4.0 * log10(flxratio/threshold)/log10(threshold) + 1.0; s << lwidth << ")\""; s << ", arrowsize=" @@ -318,13 +306,12 @@ void ReactionPathDiagram::exportToDot(ostream& s) } } else { for (size_t i = 0; i < nPaths(); i++) { - p = path(i); - flmax = std::max(p->flow(), flmax); + flmax = std::max(path(i)->flow(), flmax); } for (size_t i = 0; i < nPaths(); i++) { - p = path(i); - flxratio = p->flow()/flmax; + Path* p = path(i); + double flxratio = p->flow()/flmax; if (m_local != npos) { if (p->begin()->number != m_local && p->end()->number != m_local) { @@ -338,7 +325,7 @@ void ReactionPathDiagram::exportToDot(ostream& s) << " -> s" << p->end()->number; if (arrow_width < 0) { - lwidth = 1.0 - 4.0 * log10(flxratio/threshold)/log10(threshold) + double lwidth = 1.0 - 4.0 * log10(flxratio/threshold)/log10(threshold) + 1.0; s << "[fontname=\""+m_font+"\", style=\"setlinewidth(" << lwidth @@ -394,11 +381,9 @@ void ReactionPathDiagram::addNode(size_t k, const string& nm, doublereal x) void ReactionPathDiagram::linkNodes(size_t k1, size_t k2, size_t rxn, doublereal value, string legend) { - SpeciesNode* begin = m_nodes[k1]; - SpeciesNode* end = m_nodes[k2]; Path* ff = m_paths[k1][k2]; if (!ff) { - ff= new Path(begin, end); + ff= new Path(m_nodes[k1], m_nodes[k2]); m_paths[k1][k2] = ff; m_pathlist.push_back(ff); } @@ -419,8 +404,6 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s) logfile << endl << "Reaction " << i+1 << ": " << s.reactionString(i); - size_t nrnet = m_reac[i].size(); - size_t npnet = m_prod[i].size(); std::vector r, p; for (size_t k = 0; k < s.nTotalSpecies(); k++) { if (s.reactantStoichCoeff(k,i)) { @@ -431,12 +414,9 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s) } } - Group b0, b1, bb; - vector& e = m_elementSymbols; - if (m_determinate[i]) { logfile << " ... OK." << endl; - } else if (nrnet == 2 && npnet == 2) { + } else if (m_reac[i].size() == 2 && m_prod[i].size() == 2) { // indices for the two reactants size_t kr0 = m_reac[i][0]; size_t kr1 = m_reac[i][1]; @@ -453,8 +433,8 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s) const Group* group_a0=0, *group_b0=0, *group_c0=0, *group_a1=0, *group_b1=0, *group_c1=0; - b0 = p0 - r0; - b1 = p1 - r0; + Group b0 = p0 - r0; + Group b1 = p1 - r0; if (b0.valid() && b1.valid()) { logfile << " ... ambiguous." << endl; } else if (!b0.valid() && !b1.valid()) { @@ -481,15 +461,15 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s) m_transfer[i][kr0][kp0] = p0; } logfile << " "; - group_a0->fmt(logfile,e); + group_a0->fmt(logfile, m_elementSymbols); logfile << " + "; - group_b0->fmt(logfile,e); - group_c0->fmt(logfile,e); + group_b0->fmt(logfile,m_elementSymbols); + group_c0->fmt(logfile, m_elementSymbols); logfile << " = "; - group_a0->fmt(logfile,e); - group_b0->fmt(logfile,e); + group_a0->fmt(logfile, m_elementSymbols); + group_b0->fmt(logfile, m_elementSymbols); logfile << " + "; - group_c0->fmt(logfile,e); + group_c0->fmt(logfile, m_elementSymbols); if (b1.valid()) { logfile << " [<= default] " << endl; } else { @@ -519,15 +499,15 @@ int ReactionPathBuilder::findGroups(ostream& logfile, Kinetics& s) } } logfile << " "; - group_a1->fmt(logfile,e); + group_a1->fmt(logfile, m_elementSymbols); logfile << " + "; - group_b1->fmt(logfile,e); - group_c1->fmt(logfile,e); + group_b1->fmt(logfile, m_elementSymbols); + group_c1->fmt(logfile, m_elementSymbols); logfile << " = "; - group_a1->fmt(logfile,e); - group_b1->fmt(logfile,e); + group_a1->fmt(logfile, m_elementSymbols); + group_b1->fmt(logfile, m_elementSymbols); logfile << " + "; - group_c1->fmt(logfile,e); + group_c1->fmt(logfile, m_elementSymbols); logfile << endl; } } else { @@ -544,17 +524,13 @@ void ReactionPathBuilder::writeGroup(ostream& out, const Group& g) void ReactionPathBuilder::findElements(Kinetics& kin) { - string ename; m_enamemap.clear(); m_nel = 0; - size_t np = kin.nPhases(); - ThermoPhase* p; - for (size_t i = 0; i < np; i++) { - p = &kin.thermo(i); + for (size_t i = 0; i < kin.nPhases(); i++) { + ThermoPhase* p = &kin.thermo(i); // iterate over the elements in this phase - size_t nel = p->nElements(); - for (size_t m = 0; m < nel; m++) { - ename = p->elementName(m); + for (size_t m = 0; m < p->nElements(); m++) { + string ename = p->elementName(m); // if no entry is found for this element name, then it is a new // element. In this case, add the name to the list of names, @@ -568,17 +544,14 @@ void ReactionPathBuilder::findElements(Kinetics& kin) } } m_atoms.resize(kin.nTotalSpecies(), m_nel, 0.0); - string sym; // iterate over the elements for (size_t m = 0; m < m_nel; m++) { - sym = m_elementSymbols[m]; size_t k = 0; // iterate over the phases - for (size_t ip = 0; ip < np; ip++) { + for (size_t ip = 0; ip < kin.nPhases(); ip++) { ThermoPhase* p = &kin.thermo(ip); - size_t nsp = p->nSpecies(); - size_t mlocal = p->elementIndex(sym); - for (size_t kp = 0; kp < nsp; kp++) { + size_t mlocal = p->elementIndex(m_elementSymbols[m]); + for (size_t kp = 0; kp < p->nSpecies(); kp++) { if (mlocal != npos) { m_atoms(k, m) = p->nAtoms(kp, mlocal); } @@ -621,18 +594,14 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) m_x.resize(m_ns); // not currently used ? m_elatoms.resize(m_nel, m_nr); - size_t nr, np, n, k; - size_t nmol; - map net; - for (size_t i = 0; i < m_nr; i++) { // construct the lists of reactant and product indices, not including // molecules that appear on both sides. m_reac[i].clear(); m_prod[i].clear(); - net.clear(); - nr = allReactants[i].size(); - np = allProducts[i].size(); + map net; + size_t nr = allReactants[i].size(); + size_t np = allProducts[i].size(); for (size_t ir = 0; ir < nr; ir++) { net[allReactants[i][ir]]--; } @@ -640,14 +609,14 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) net[allProducts[i][ip]]++; } - for (k = 0; k < m_ns; k++) { + for (size_t k = 0; k < m_ns; k++) { if (net[k] < 0) { - nmol = -net[k]; + size_t nmol = -net[k]; for (size_t jr = 0; jr < nmol; jr++) { m_reac[i].push_back(k); } } else if (net[k] > 0) { - nmol = net[k]; + size_t nmol = net[k]; for (size_t jp = 0; jp < nmol; jp++) { m_prod[i].push_back(k); } @@ -659,8 +628,8 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) // compute number of atoms of each element in each reaction, excluding // molecules that appear on both sides of the reaction. We only need to // compute this for the reactants, since the elements are conserved. - for (n = 0; n < nrnet; n++) { - k = m_reac[i][n]; + for (size_t n = 0; n < nrnet; n++) { + size_t k = m_reac[i][n]; for (size_t m = 0; m < m_nel; m++) { m_elatoms(m,i) += m_atoms(k,m); } @@ -683,14 +652,13 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) // element, *and* more than one product contains the element. In this case, // additional information is needed to determine the partitioning of the // reactant atoms of that element among the products. - int nar, nap; for (size_t i = 0; i < m_nr; i++) { - nr = m_reac[i].size(); - np = m_prod[i].size(); + size_t nr = m_reac[i].size(); + size_t np = m_prod[i].size(); m_determinate[i] = true; for (size_t m = 0; m < m_nel; m++) { - nar = 0; - nap = 0; + int nar = 0; + int nap = 0; for (size_t j = 0; j < nr; j++) { if (m_atoms(m_reac[i][j],m) > 0) { nar++; @@ -732,11 +700,8 @@ string reactionLabel(size_t i, size_t kr, size_t nr, int ReactionPathBuilder::build(Kinetics& s, const string& element, ostream& output, ReactionPathDiagram& r, bool quiet) { - doublereal f, ropf, ropr, fwd, rev; - string fwdlabel, revlabel; map warn; doublereal threshold = 0.0; - bool fwd_incl, rev_incl, force_incl; size_t m = m_enamemap[element]-1; r.element = element; if (m == npos) { @@ -750,8 +715,7 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, vector& in_nodes = r.included(); vector& out_nodes = r.excluded(); - vector_int status; - status.resize(s.nTotalSpecies(), 0); + vector_int status(s.nTotalSpecies(), 0); for (size_t ni = 0; ni < in_nodes.size(); ni++) { status[s.kineticsSpeciesIndex(in_nodes[ni])] = 1; } @@ -760,8 +724,8 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, } for (size_t i = 0; i < m_nr; i++) { - ropf = m_ropf[i]; - ropr = m_ropr[i]; + double ropf = m_ropf[i]; + double ropr = m_ropr[i]; // loop over reactions involving element m if (m_elatoms(m, i) > 0) { @@ -770,11 +734,11 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, for (size_t kr = 0; kr < nr; kr++) { size_t kkr = m_reac[i][kr]; - fwdlabel = reactionLabel(i, kr, nr, m_reac[i], s); + string fwdlabel = reactionLabel(i, kr, nr, m_reac[i], s); for (size_t kp = 0; kp < np; kp++) { size_t kkp = m_prod[i][kp]; - revlabel = ""; + string revlabel = ""; for (size_t j = 0; j < np; j++) { if (j != kp) { revlabel += " + "+ s.kineticsSpeciesName(m_prod[i][j]); @@ -797,6 +761,7 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, // the type of reaction to determine which reactant // species was the source of a given m-atom in the // product + double f; if ((m_atoms(kkp,m) < m_elatoms(m, i)) && (m_atoms(kkr,m) < m_elatoms(m, i))) { map >& g = m_transfer[i]; @@ -827,14 +792,14 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, f = m_atoms(kkp,m) * m_atoms(kkr,m) / m_elatoms(m, i); } - fwd = ropf*f; - rev = ropr*f; - force_incl = ((status[kkr] == 1) || (status[kkp] == 1)); + double fwd = ropf*f; + double rev = ropr*f; + bool force_incl = ((status[kkr] == 1) || (status[kkp] == 1)); - fwd_incl = ((fwd > threshold) || - (fwd > 0.0 && force_incl)); - rev_incl = ((rev > threshold) || - (rev > 0.0 && force_incl)); + bool fwd_incl = ((fwd > threshold) || + (fwd > 0.0 && force_incl)); + bool rev_incl = ((rev > threshold) || + (rev > 0.0 && force_incl)); if (fwd_incl || rev_incl) { if (!r.hasNode(kkr)) { r.addNode(kkr, s.kineticsSpeciesName(kkr), m_x[kkr]); diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index bd030e338..18acd4b47 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -154,20 +154,15 @@ bool importKinetics(const XML_Node& phase, std::vector th, } phase_ids.push_back(owning_phase); - int np = static_cast(phase_ids.size()); - int nt = static_cast(th.size()); - // for each referenced phase, attempt to find its id among those // phases specified. - bool phase_ok; - string phase_id; string msg = ""; - for (int n = 0; n < np; n++) { - phase_id = phase_ids[n]; - phase_ok = false; + for (size_t n = 0; n < phase_ids.size(); n++) { + string phase_id = phase_ids[n]; + bool phase_ok = false; // loop over the supplied 'ThermoPhase' objects representing // phases, to find an object with the same id. - for (int m = 0; m < nt; m++) { + for (size_t m = 0; m < th.size(); m++) { if (th[m]->id() == phase_id) { phase_ok = true; // if no phase with this id has been added to @@ -195,8 +190,7 @@ bool importKinetics(const XML_Node& phase, std::vector th, bool buildSolutionFromXML(XML_Node& root, const std::string& id, const std::string& nm, ThermoPhase* th, Kinetics* kin) { - XML_Node* x; - x = get_XML_NameID(nm, string("#")+id, &root); + XML_Node* x = get_XML_NameID(nm, string("#")+id, &root); if (!x) { return false; } @@ -207,8 +201,7 @@ bool buildSolutionFromXML(XML_Node& root, const std::string& id, // Create a vector of ThermoPhase pointers of length 1 having the current th // ThermoPhase as the entry. - std::vector phases(1); - phases[0] = th; + std::vector phases{th}; // Fill in the kinetics object k, by querying the const XML_Node tree // located by x. The source terms and eventually the source term vector will diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp index 8e4b3997b..965d01c5f 100644 --- a/src/kinetics/solveSP.cpp +++ b/src/kinetics/solveSP.cpp @@ -39,10 +39,9 @@ solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) : m_ioflag(0) { m_numSurfPhases = 0; - size_t numPossibleSurfPhases = m_objects.size(); - for (size_t n = 0; n < numPossibleSurfPhases; n++) { - InterfaceKinetics* m_kin = m_objects[n]; - size_t surfPhaseIndex = m_kin->surfacePhaseIndex(); + for (size_t n = 0; n < m_objects.size(); n++) { + InterfaceKinetics* kin = m_objects[n]; + size_t surfPhaseIndex = kin->surfacePhaseIndex(); if (surfPhaseIndex != npos) { m_numSurfPhases++; m_indexKinObjSurfPhase.push_back(n); @@ -51,8 +50,7 @@ solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) : throw CanteraError("solveSP", "InterfaceKinetics object has no surface phase"); } - ThermoPhase* tp = &m_kin->thermo(surfPhaseIndex); - SurfPhase* sp = dynamic_cast(tp); + SurfPhase* sp = dynamic_cast(&kin->thermo(surfPhaseIndex)); if (!sp) { throw CanteraError("solveSP", "Inconsistent ThermoPhase object within " @@ -91,15 +89,14 @@ solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) : m_eqnIndexStartSolnPhase.resize(m_numSurfPhases + m_numBulkPhasesSS, 0); size_t kindexSP = 0; - size_t isp, k, nsp, kstart; - for (isp = 0; isp < m_numSurfPhases; isp++) { + for (size_t isp = 0; isp < m_numSurfPhases; isp++) { size_t iKinObject = m_indexKinObjSurfPhase[isp]; - InterfaceKinetics* m_kin = m_objects[iKinObject]; + InterfaceKinetics* kin = m_objects[iKinObject]; size_t surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp]; - kstart = m_kin->kineticsSpeciesIndex(0, surfPhaseIndex); - nsp = m_nSpeciesSurfPhase[isp]; + size_t kstart = kin->kineticsSpeciesIndex(0, surfPhaseIndex); + size_t nsp = m_nSpeciesSurfPhase[isp]; m_eqnIndexStartSolnPhase[isp] = kindexSP; - for (k = 0; k < nsp; k++, kindexSP++) { + for (size_t k = 0; k < nsp; k++, kindexSP++) { m_kinSpecIndex[kindexSP] = kstart + k; m_kinObjIndex[kindexSP] = isp; } @@ -123,7 +120,6 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, if (ifunc == SFLUX_JACOBIAN) { EXTRA_ACCURACY *= 0.001; } - int info = 0; int label_t=-1; // Species IDs for time control int label_d = -1; // Species IDs for damping control int label_t_old=-1; @@ -131,10 +127,7 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, int iter=0; // iteration number on numlinear solver int iter_max=1000; // maximum number of nonlinear iterations doublereal deltaT = 1.0E-10; // Delta time step - doublereal damp=1.0, tmp; - // Weighted L2 norm of the residual. Currently, this is only used for IO - // purposes. It doesn't control convergence. - doublereal resid_norm; + doublereal damp=1.0; doublereal inv_t = 0.0; doublereal t_real = 0.0, update_norm = 1.0E6; bool do_time = false, not_converged = true; @@ -149,10 +142,8 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, // CSoln, and in an separate vector CSolnInit. size_t loc = 0; for (size_t n = 0; n < m_numSurfPhases; n++) { - SurfPhase* sf_ptr = m_ptrsSurfPhase[n]; - sf_ptr->getConcentrations(m_numEqn1.data()); - size_t nsp = m_nSpeciesSurfPhase[n]; - for (size_t k = 0; k getConcentrations(m_numEqn1.data()); + for (size_t k = 0; k < m_nSpeciesSurfPhase[n]; k++) { m_CSolnSP[loc] = m_numEqn1[k]; loc++; } @@ -192,7 +183,7 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, if (damp < 1.0) { label_factor = 1.0; } - tmp = calc_t(m_netProductionRatesSave.data(), + double tmp = calc_t(m_netProductionRatesSave.data(), m_XMolKinSpecies.data(), &label_t, &label_t_old, &label_factor, m_ioflag); if (iter < 10) { @@ -230,10 +221,10 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, } // Find the weighted norm of the residual - resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq); + double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq); // Solve Linear system. The solution is in resid[] - info = m_Jac.factor(); + int info = m_Jac.factor(); if (info==0) { m_Jac.solve(&m_resid[0]); } else { @@ -324,7 +315,7 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, if (m_ioflag) { fun_eval(m_resid.data(), m_CSolnSP.data(), m_CSolnSPOld.data(), false, deltaT); - resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq); + double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq); printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter, update_norm, resid_norm, do_time, true); } @@ -355,12 +346,10 @@ void solveSP::updateMFSolnSP(doublereal* XMolSolnSP) void solveSP::updateMFKinSpecies(doublereal* XMolKinSpecies, int isp) { - InterfaceKinetics* m_kin = m_objects[isp]; - size_t nph = m_kin->nPhases(); - for (size_t iph = 0; iph < nph; iph++) { - size_t ksi = m_kin->kineticsSpeciesIndex(0, iph); - ThermoPhase& thref = m_kin->thermo(iph); - thref.getMoleFractions(XMolKinSpecies + ksi); + InterfaceKinetics* kin = m_objects[isp]; + for (size_t iph = 0; iph < kin->nPhases(); iph++) { + size_t ksi = kin->kineticsSpeciesIndex(0, iph); + kin->thermo(iph).getMoleFractions(XMolKinSpecies + ksi); } } @@ -368,11 +357,10 @@ void solveSP::evalSurfLarge(const doublereal* CSolnSP) { size_t kindexSP = 0; for (size_t isp = 0; isp < m_numSurfPhases; isp++) { - size_t nsp = m_nSpeciesSurfPhase[isp]; doublereal Clarge = CSolnSP[kindexSP]; m_spSurfLarge[isp] = 0; kindexSP++; - for (size_t k = 1; k < nsp; k++, kindexSP++) { + for (size_t k = 1; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) { if (CSolnSP[kindexSP] > Clarge) { Clarge = CSolnSP[kindexSP]; m_spSurfLarge[isp] = k; @@ -385,10 +373,8 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln, const doublereal* CSolnOld, const bool do_time, const doublereal deltaT) { - size_t isp, nsp, kstart, k, kindexSP, kins, kspecial; + size_t k; doublereal lenScale = 1.0E-9; - doublereal sd = 0.0; - doublereal grRate; if (m_numSurfPhases > 0) { // update the surface concentrations with the input surface // concentration vector @@ -400,13 +386,13 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln, // HKM Should do it here for all kinetics objects so that // bulk will eventually work. if (do_time) { - kindexSP = 0; - for (isp = 0; isp < m_numSurfPhases; isp++) { - nsp = m_nSpeciesSurfPhase[isp]; + size_t kindexSP = 0; + for (size_t isp = 0; isp < m_numSurfPhases; isp++) { + size_t nsp = m_nSpeciesSurfPhase[isp]; InterfaceKinetics* kinPtr = m_objects[isp]; size_t surfIndex = kinPtr->surfacePhaseIndex(); - kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex); - kins = kindexSP; + size_t kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex); + size_t kins = kindexSP; kinPtr->getNetProductionRates(m_netProductionRatesSave.data()); for (k = 0; k < nsp; k++, kindexSP++) { resid[kindexSP] = @@ -414,27 +400,27 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln, - m_netProductionRatesSave[kstart + k]; } - kspecial = kins + m_spSurfLarge[isp]; - sd = m_ptrsSurfPhase[isp]->siteDensity(); + size_t kspecial = kins + m_spSurfLarge[isp]; + double sd = m_ptrsSurfPhase[isp]->siteDensity(); resid[kspecial] = sd; for (k = 0; k < nsp; k++) { resid[kspecial] -= CSoln[kins + k]; } } } else { - kindexSP = 0; - for (isp = 0; isp < m_numSurfPhases; isp++) { - nsp = m_nSpeciesSurfPhase[isp]; + size_t kindexSP = 0; + for (size_t isp = 0; isp < m_numSurfPhases; isp++) { + size_t nsp = m_nSpeciesSurfPhase[isp]; InterfaceKinetics* kinPtr = m_objects[isp]; size_t surfIndex = kinPtr->surfacePhaseIndex(); - kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex); - kins = kindexSP; + size_t kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex); + size_t kins = kindexSP; kinPtr->getNetProductionRates(m_netProductionRatesSave.data()); for (k = 0; k < nsp; k++, kindexSP++) { resid[kindexSP] = - m_netProductionRatesSave[kstart + k]; } - kspecial = kins + m_spSurfLarge[isp]; - sd = m_ptrsSurfPhase[isp]->siteDensity(); + size_t kspecial = kins + m_spSurfLarge[isp]; + double sd = m_ptrsSurfPhase[isp]->siteDensity(); resid[kspecial] = sd; for (k = 0; k < nsp; k++) { resid[kspecial] -= CSoln[kins + k]; @@ -443,14 +429,14 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln, } if (m_bulkFunc == BULK_DEPOSITION) { - kindexSP = m_numTotSurfSpecies; - for (isp = 0; isp < m_numBulkPhasesSS; isp++) { + size_t kindexSP = m_numTotSurfSpecies; + for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) { doublereal* XBlk = m_numEqn1.data(); - nsp = m_nSpeciesSurfPhase[isp]; + size_t nsp = m_nSpeciesSurfPhase[isp]; size_t surfPhaseIndex = m_indexKinObjSurfPhase[isp]; - InterfaceKinetics* m_kin = m_objects[isp]; - grRate = 0.0; - kstart = m_kin->kineticsSpeciesIndex(0, surfPhaseIndex); + InterfaceKinetics* kin = m_objects[isp]; + double grRate = 0.0; + size_t kstart = kin->kineticsSpeciesIndex(0, surfPhaseIndex); for (k = 0; k < nsp; k++) { if (m_netProductionRatesSave[kstart + k] > 0.0) { grRate += m_netProductionRatesSave[kstart + k]; @@ -471,6 +457,7 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln, } } else { grRate = 1.0E-6; + //! @todo the appearance of k in this formula is suspicious grRate += fabs(m_netProductionRatesSave[kstart + k]); for (k = 1; k < nsp; k++) { resid[kindexSP + k] = grRate * (XBlk[k] - 1.0/nsp); @@ -494,20 +481,19 @@ void solveSP::resjac_eval(SquareMatrix& jac, const doublereal CSolnOld[], const bool do_time, const doublereal deltaT) { - size_t kColIndex = 0, nsp, jsp, i, kCol; - doublereal dc, cSave, sd; + size_t kColIndex = 0; // Calculate the residual fun_eval(resid, CSoln, CSolnOld, do_time, deltaT); // Now we will look over the columns perturbing each unknown. - for (jsp = 0; jsp < m_numSurfPhases; jsp++) { - nsp = m_nSpeciesSurfPhase[jsp]; - sd = m_ptrsSurfPhase[jsp]->siteDensity(); - for (kCol = 0; kCol < nsp; kCol++) { - cSave = CSoln[kColIndex]; - dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7); + for (size_t jsp = 0; jsp < m_numSurfPhases; jsp++) { + size_t nsp = m_nSpeciesSurfPhase[jsp]; + double sd = m_ptrsSurfPhase[jsp]->siteDensity(); + for (size_t kCol = 0; kCol < nsp; kCol++) { + double cSave = CSoln[kColIndex]; + double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7); CSoln[kColIndex] += dc; fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT); - for (i = 0; i < m_neq; i++) { + for (size_t i = 0; i < m_neq; i++) { jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc; } CSoln[kColIndex] = cSave; @@ -516,15 +502,15 @@ void solveSP::resjac_eval(SquareMatrix& jac, } if (m_bulkFunc == BULK_DEPOSITION) { - for (jsp = 0; jsp < m_numBulkPhasesSS; jsp++) { - nsp = m_numBulkSpecies[jsp]; - sd = m_bulkPhasePtrs[jsp]->molarDensity(); - for (kCol = 0; kCol < nsp; kCol++) { - cSave = CSoln[kColIndex]; - dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7); + for (size_t jsp = 0; jsp < m_numBulkPhasesSS; jsp++) { + size_t nsp = m_numBulkSpecies[jsp]; + double sd = m_bulkPhasePtrs[jsp]->molarDensity(); + for (size_t kCol = 0; kCol < nsp; kCol++) { + double cSave = CSoln[kColIndex]; + double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7); CSoln[kColIndex] += dc; fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT); - for (i = 0; i < m_neq; i++) { + for (size_t i = 0; i < m_neq; i++) { jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc; } CSoln[kColIndex] = cSave; @@ -548,19 +534,19 @@ void solveSP::resjac_eval(SquareMatrix& jac, static doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, int* label) { const doublereal APPROACH = 0.80; - doublereal damp = 1.0, xnew, xtop, xbot; - static doublereal damp_old = 1.0; + doublereal damp = 1.0; + static doublereal damp_old = 1.0; //! @todo this variable breaks thread safety *label = -1; for (size_t i = 0; i < dim; i++) { // Calculate the new suggested new value of x[i] - xnew = x[i] - damp * dxneg[i]; + double xnew = x[i] - damp * dxneg[i]; // Calculate the allowed maximum and minimum values of x[i] // - Only going to allow x[i] to converge to zero by a // single order of magnitude at a time - xtop = 1.0 - 0.1*fabs(1.0-x[i]); - xbot = fabs(x[i]*0.1) - 1.0e-16; + double xtop = 1.0 - 0.1*fabs(1.0-x[i]); + double xbot = fabs(x[i]*0.1) - 1.0e-16; if (xnew > xtop) { damp = - APPROACH * (1.0 - x[i]) / dxneg[i]; *label = int(i); @@ -594,13 +580,11 @@ static doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, i static doublereal calcWeightedNorm(const doublereal wtX[], const doublereal dx[], size_t dim) { doublereal norm = 0.0; - doublereal tmp; if (dim == 0) { return 0.0; } for (size_t i = 0; i < dim; i++) { - tmp = dx[i] / wtX[i]; - norm += tmp * tmp; + norm += pow(dx[i] / wtX[i], 2); } return sqrt(norm/dim); } @@ -609,24 +593,19 @@ void solveSP::calcWeights(doublereal wtSpecies[], doublereal wtResid[], const Array2D& Jac, const doublereal CSoln[], const doublereal abstol, const doublereal reltol) { - size_t k, jcol, kindex, isp, nsp; - doublereal sd; - // First calculate the weighting factor for the concentrations of the // surface species and bulk species. - kindex = 0; - for (isp = 0; isp < m_numSurfPhases; isp++) { - nsp = m_nSpeciesSurfPhase[isp]; - sd = m_ptrsSurfPhase[isp]->siteDensity(); - for (k = 0; k < nsp; k++, kindex++) { + size_t kindex = 0; + for (size_t isp = 0; isp < m_numSurfPhases; isp++) { + double sd = m_ptrsSurfPhase[isp]->siteDensity(); + for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindex++) { wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]); } } if (m_bulkFunc == BULK_DEPOSITION) { - for (isp = 0; isp < m_numBulkPhasesSS; isp++) { - nsp = m_numBulkSpecies[isp]; - sd = m_bulkPhasePtrs[isp]->molarDensity(); - for (k = 0; k < nsp; k++, kindex++) { + for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) { + double sd = m_bulkPhasePtrs[isp]->molarDensity(); + for (size_t k = 0; k < m_numBulkSpecies[isp]; k++, kindex++) { wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]); } } @@ -635,9 +614,9 @@ void solveSP::calcWeights(doublereal wtSpecies[], doublereal wtResid[], // Now do the residual Weights. Since we have the Jacobian, we will use it // to generate a number based on the what a significant change in a solution // variable does to each residual. This is a row sum scale operation. - for (k = 0; k < m_neq; k++) { + for (size_t k = 0; k < m_neq; k++) { wtResid[k] = 0.0; - for (jcol = 0; jcol < m_neq; jcol++) { + for (size_t jcol = 0; jcol < m_neq; jcol++) { wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]); } } @@ -648,33 +627,24 @@ doublereal solveSP::calc_t(doublereal netProdRateSolnSP[], int* label, int* label_old, doublereal* label_factor, int ioflag) { - size_t k, isp, nsp, kstart; doublereal inv_timeScale = 1.0E-10; - doublereal sden, tmp; size_t kindexSP = 0; *label = 0; updateMFSolnSP(XMolSolnSP); - for (isp = 0; isp < m_numSurfPhases; isp++) { - nsp = m_nSpeciesSurfPhase[isp]; - + for (size_t isp = 0; isp < m_numSurfPhases; isp++) { // Get the interface kinetics associated with this surface - InterfaceKinetics* m_kin = m_objects[isp]; + InterfaceKinetics* kin = m_objects[isp]; // Calculate the start of the species index for surfaces within // the InterfaceKinetics object - size_t surfIndex = m_kin->surfacePhaseIndex(); - kstart = m_kin->kineticsSpeciesIndex(0, surfIndex); - ThermoPhase& THref = m_kin->thermo(surfIndex); - m_kin->getNetProductionRates(m_numEqn1.data()); - sden = THref.molarDensity(); - for (k = 0; k < nsp; k++, kindexSP++) { + size_t surfIndex = kin->surfacePhaseIndex(); + size_t kstart = kin->kineticsSpeciesIndex(0, surfIndex); + kin->getNetProductionRates(m_numEqn1.data()); + double sden = kin->thermo(surfIndex).molarDensity(); + for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) { size_t kspindex = kstart + k; netProdRateSolnSP[kindexSP] = m_numEqn1[kspindex]; - if (XMolSolnSP[kindexSP] <= 1.0E-10) { - tmp = 1.0E-10; - } else { - tmp = XMolSolnSP[kindexSP]; - } + double tmp = std::max(XMolSolnSP[kindexSP], 1.0e-10); tmp *= sden; tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp); if (netProdRateSolnSP[kindexSP]> 0.0) { @@ -752,8 +722,6 @@ void solveSP::printIteration(int ioflag, doublereal damp, int label_d, size_t iter, doublereal update_norm, doublereal resid_norm, bool do_time, bool final) { - size_t k; - string nm; if (ioflag == 1) { if (final) { writelogf("\tFIN%3d ", iter); @@ -772,20 +740,16 @@ void solveSP::printIteration(int ioflag, doublereal damp, int label_d, } writelogf("%9.4e %9.4e", update_norm, resid_norm); if (do_time) { - k = m_kinSpecIndex[label_t]; + size_t k = m_kinSpecIndex[label_t]; size_t isp = m_kinObjIndex[label_t]; - InterfaceKinetics* m_kin = m_objects[isp]; - nm = m_kin->kineticsSpeciesName(k); - writelog(" %-16s", nm); + writelog(" %-16s", m_objects[isp]->kineticsSpeciesName(k)); } else { writeline(' ', 16, false); } if (label_d >= 0) { - k = m_kinSpecIndex[label_d]; + size_t k = m_kinSpecIndex[label_d]; size_t isp = m_kinObjIndex[label_d]; - InterfaceKinetics* m_kin = m_objects[isp]; - nm = m_kin->kineticsSpeciesName(k); - writelogf(" %-16s", nm); + writelogf(" %-16s", m_objects[isp]->kineticsSpeciesName(k)); } if (final) { writelog(" -- success");