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

[Kinetics] Make better use of local variables

Browse source
This commit is contained in:
Ray Speth 2016-05-11 16:22:41 -04:00
parent 248874de83
commit 920ff1d897
8 changed files with 217 additions and 319 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

11
src/kinetics/Falloff.cpp
View file

@ -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);
}

11
src/kinetics/Group.cpp
View file

@ -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<std::string>& 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 << "-";

18
src/kinetics/ImplicitSurfChem.cpp
View file

@ -28,24 +28,23 @@ ImplicitSurfChem::ImplicitSurfChem(vector<InterfaceKinetics*> 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];

26
src/kinetics/InterfaceKinetics.cpp
View file

@ -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<InterfaceKinetics*> k;
k.push_back(this);
vector<InterfaceKinetics*> 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<InterfaceKinetics*> k;
k.push_back(this);
vector<InterfaceKinetics*> k{this};
m_integrator = new ImplicitSurfChem(k);
m_integrator->initialize();
}

20
src/kinetics/Kinetics.cpp
View file

@ -121,12 +121,11 @@ void Kinetics::checkSpeciesArraySize(size_t kk) const
void Kinetics::assignShallowPointers(const std::vector<thermo_t*> & 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<Reaction> 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);
}

225
src/kinetics/ReactionPath.cpp
View file

@ -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<size_t> 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<string>& 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<size_t, int> 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<size_t, int> 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<size_t, int> 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<string>& in_nodes = r.included();
vector<string>& 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<size_t, map<size_t, Group> >& 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]);

19
src/kinetics/importKinetics.cpp
View file

@ -154,20 +154,15 @@ bool importKinetics(const XML_Node& phase, std::vector<ThermoPhase*> th,
}
phase_ids.push_back(owning_phase);
int np = static_cast<int>(phase_ids.size());
int nt = static_cast<int>(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<ThermoPhase*> 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<ThermoPhase*> phases(1);
phases[0] = th;
std::vector<ThermoPhase*> 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

206
src/kinetics/solveSP.cpp
View file

@ -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<SurfPhase*>(tp);
SurfPhase* sp = dynamic_cast<SurfPhase*>(&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 <nsp; k++) {
m_ptrsSurfPhase[n]->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");