Cantera  2.5.1
solveSP.cpp
1 /*
2  * @file: solveSP.cpp Implicit surface site concentration solver
3  */
4 
5 // This file is part of Cantera. See License.txt in the top-level directory or
6 // at https://cantera.org/license.txt for license and copyright information.
7 
11 
12 using namespace std;
13 namespace Cantera
14 {
15 
16 // STATIC ROUTINES DEFINED IN THIS FILE
17 
18 static doublereal calc_damping(doublereal* x, doublereal* dx, size_t dim, int*);
19 static doublereal calcWeightedNorm(const doublereal [], const doublereal dx[], size_t);
20 
21 // solveSP Class Definitions
22 
23 solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) :
24  m_SurfChemPtr(surfChemPtr),
25  m_objects(surfChemPtr->getObjects()),
26  m_neq(0),
27  m_bulkFunc(bulkFunc),
28  m_numSurfPhases(0),
29  m_numTotSurfSpecies(0),
30  m_numBulkPhasesSS(0),
31  m_numTotBulkSpeciesSS(0),
32  m_atol(1.0E-15),
33  m_rtol(1.0E-4),
34  m_maxstep(1000),
35  m_maxTotSpecies(0),
36  m_ioflag(0)
37 {
38  m_numSurfPhases = 0;
39  for (size_t n = 0; n < m_objects.size(); n++) {
40  InterfaceKinetics* kin = m_objects[n];
41  size_t surfPhaseIndex = kin->surfacePhaseIndex();
42  if (surfPhaseIndex != npos) {
44  m_indexKinObjSurfPhase.push_back(n);
45  m_kinObjPhaseIDSurfPhase.push_back(surfPhaseIndex);
46  } else {
47  throw CanteraError("solveSP::solveSP",
48  "InterfaceKinetics object has no surface phase");
49  }
50  SurfPhase* sp = dynamic_cast<SurfPhase*>(&kin->thermo(surfPhaseIndex));
51  if (!sp) {
52  throw CanteraError("solveSP::solveSP",
53  "Inconsistent ThermoPhase object within "
54  "InterfaceKinetics object");
55  }
56 
57  m_ptrsSurfPhase.push_back(sp);
58  size_t nsp = sp->nSpecies();
59  m_nSpeciesSurfPhase.push_back(nsp);
60  m_numTotSurfSpecies += nsp;
61  }
62  // We rely on ordering to figure things out
64 
65  if (bulkFunc == BULK_DEPOSITION) {
67  } else {
69  }
70 
71  m_maxTotSpecies = 0;
72  for (size_t n = 0; n < m_numSurfPhases; n++) {
73  size_t tsp = m_objects[n]->nTotalSpecies();
74  m_maxTotSpecies = std::max(m_maxTotSpecies, tsp);
75  }
77 
79  m_numEqn1.resize(m_maxTotSpecies, 0.0);
80  m_numEqn2.resize(m_maxTotSpecies, 0.0);
81  m_XMolKinSpecies.resize(m_maxTotSpecies, 0.0);
82  m_CSolnSave.resize(m_neq, 0.0);
83  m_spSurfLarge.resize(m_numSurfPhases, 0);
87 
88  size_t kindexSP = 0;
89  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
90  size_t iKinObject = m_indexKinObjSurfPhase[isp];
91  InterfaceKinetics* kin = m_objects[iKinObject];
92  size_t surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp];
93  size_t kstart = kin->kineticsSpeciesIndex(0, surfPhaseIndex);
94  size_t nsp = m_nSpeciesSurfPhase[isp];
95  m_eqnIndexStartSolnPhase[isp] = kindexSP;
96  for (size_t k = 0; k < nsp; k++, kindexSP++) {
97  m_kinSpecIndex[kindexSP] = kstart + k;
98  m_kinObjIndex[kindexSP] = isp;
99  }
100  }
101 
102  // Dimension solution vector
103  size_t dim1 = std::max<size_t>(1, m_neq);
104  m_CSolnSP.resize(dim1, 0.0);
105  m_CSolnSPInit.resize(dim1, 0.0);
106  m_CSolnSPOld.resize(dim1, 0.0);
107  m_wtResid.resize(dim1, 0.0);
108  m_wtSpecies.resize(dim1, 0.0);
109  m_resid.resize(dim1, 0.0);
110  m_Jac.resize(dim1, dim1, 0.0);
111 }
112 
113 int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
114  doublereal PGas, doublereal reltol, doublereal abstol)
115 {
116  doublereal EXTRA_ACCURACY = 0.001;
117  if (ifunc == SFLUX_JACOBIAN) {
118  EXTRA_ACCURACY *= 0.001;
119  }
120  int label_t=-1; // Species IDs for time control
121  int label_d = -1; // Species IDs for damping control
122  int label_t_old=-1;
123  doublereal label_factor = 1.0;
124  int iter=0; // iteration number on nonlinear solver
125  int iter_max=1000; // maximum number of nonlinear iterations
126  doublereal deltaT = 1.0E-10; // Delta time step
127  doublereal damp=1.0;
128  doublereal inv_t = 0.0;
129  doublereal t_real = 0.0, update_norm = 1.0E6;
130  bool do_time = false, not_converged = true;
131  m_ioflag = std::min(m_ioflag, 1);
132 
133  // Set the initial value of the do_time parameter
134  if (ifunc == SFLUX_INITIALIZE || ifunc == SFLUX_TRANSIENT) {
135  do_time = true;
136  }
137 
138  // Store the initial guess for the surface problem in the soln vector,
139  // CSoln, and in an separate vector CSolnInit.
140  size_t loc = 0;
141  for (size_t n = 0; n < m_numSurfPhases; n++) {
142  m_ptrsSurfPhase[n]->getConcentrations(m_numEqn1.data());
143  for (size_t k = 0; k < m_nSpeciesSurfPhase[n]; k++) {
144  m_CSolnSP[loc] = m_numEqn1[k];
145  loc++;
146  }
147  }
148 
150 
151  // Calculate the largest species in each phase
152  evalSurfLarge(m_CSolnSP.data());
153 
154  if (m_ioflag) {
155  print_header(m_ioflag, ifunc, time_scale, true, reltol, abstol);
156  }
157 
158  // Quick return when there isn't a surface problem to solve
159  if (m_neq == 0) {
160  not_converged = false;
161  update_norm = 0.0;
162  }
163 
164  // Start of Newton's method
165  while (not_converged && iter < iter_max) {
166  iter++;
167  // Store previous iteration's solution in the old solution vector
169 
170  // Evaluate the largest surface species for each surface phase every
171  // 5 iterations.
172  if (iter%5 == 4) {
173  evalSurfLarge(m_CSolnSP.data());
174  }
175 
176  // Calculate the value of the time step
177  // - heuristics to stop large oscillations in deltaT
178  if (do_time) {
179  // don't hurry increase in time step at the same time as damping
180  if (damp < 1.0) {
181  label_factor = 1.0;
182  }
183  double tmp = calc_t(m_netProductionRatesSave.data(),
184  m_XMolKinSpecies.data(),
185  &label_t, &label_t_old, &label_factor, m_ioflag);
186  if (iter < 10) {
187  inv_t = tmp;
188  } else if (tmp > 2.0*inv_t) {
189  inv_t = 2.0*inv_t;
190  } else {
191  inv_t = tmp;
192  }
193 
194  // Check end condition
195  if (ifunc == SFLUX_TRANSIENT) {
196  tmp = t_real + 1.0/inv_t;
197  if (tmp > time_scale) {
198  inv_t = 1.0/(time_scale - t_real);
199  }
200  }
201  } else {
202  // make steady state calc a step of 1 million seconds to prevent
203  // singular Jacobians for some pathological cases
204  inv_t = 1.0e-6;
205  }
206  deltaT = 1.0/inv_t;
207 
208  // Call the routine to numerically evaluation the Jacobian and residual
209  // for the current iteration.
210  resjac_eval(m_Jac, m_resid.data(), m_CSolnSP.data(),
211  m_CSolnSPOld.data(), do_time, deltaT);
212 
213  // Calculate the weights. Make sure the calculation is carried out on
214  // the first iteration.
215  if (iter%4 == 1) {
216  calcWeights(m_wtSpecies.data(), m_wtResid.data(),
217  m_Jac, m_CSolnSP.data(), abstol, reltol);
218  }
219 
220  // Find the weighted norm of the residual
221  double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
222 
223  // Solve Linear system. The solution is in m_resid
224  solve(m_Jac, m_resid.data());
225 
226  // Calculate the Damping factor needed to keep all unknowns between 0
227  // and 1, and not allow too large a change (factor of 2) in any unknown.
228  damp = calc_damping(m_CSolnSP.data(), m_resid.data(), m_neq, &label_d);
229 
230  // Calculate the weighted norm of the update vector Here, resid is the
231  // delta of the solution, in concentration units.
232  update_norm = calcWeightedNorm(m_wtSpecies.data(),
233  m_resid.data(), m_neq);
234 
235  // Update the solution vector and real time Crop the concentrations to
236  // zero.
237  for (size_t irow = 0; irow < m_neq; irow++) {
238  m_CSolnSP[irow] -= damp * m_resid[irow];
239  }
240  for (size_t irow = 0; irow < m_neq; irow++) {
241  m_CSolnSP[irow] = std::max(0.0, m_CSolnSP[irow]);
242  }
243  updateState(m_CSolnSP.data());
244 
245  if (do_time) {
246  t_real += damp/inv_t;
247  }
248 
249  if (m_ioflag) {
250  printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
251  update_norm, resid_norm, do_time);
252  }
253 
254  if (ifunc == SFLUX_TRANSIENT) {
255  not_converged = (t_real < time_scale);
256  } else {
257  if (do_time) {
258  if (t_real > time_scale ||
259  (resid_norm < 1.0e-7 &&
260  update_norm*time_scale/t_real < EXTRA_ACCURACY)) {
261  do_time = false;
262  }
263  } else {
264  not_converged = ((update_norm > EXTRA_ACCURACY) ||
265  (resid_norm > EXTRA_ACCURACY));
266  }
267  }
268  } // End of Newton's Method while statement
269 
270  // End Newton's method. If not converged, print error message and
271  // recalculate sdot's at equal site fractions.
272  if (not_converged && m_ioflag) {
273  writelog("#$#$#$# Error in solveSP $#$#$#$ \n");
274  writelogf("Newton iter on surface species did not converge, "
275  "update_norm = %e \n", update_norm);
276  writelog("Continuing anyway\n");
277  }
278 
279  // Decide on what to return in the solution vector. Right now, will always
280  // return the last solution no matter how bad
281  if (m_ioflag) {
282  fun_eval(m_resid.data(), m_CSolnSP.data(), m_CSolnSPOld.data(),
283  false, deltaT);
284  double resid_norm = calcWeightedNorm(m_wtResid.data(), m_resid.data(), m_neq);
285  printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
286  update_norm, resid_norm, do_time, true);
287  }
288 
289  // Return with the appropriate flag
290  if (update_norm > 1.0) {
291  return -1;
292  }
293  return 1;
294 }
295 
296 void solveSP::updateState(const doublereal* CSolnSP)
297 {
298  size_t loc = 0;
299  for (size_t n = 0; n < m_numSurfPhases; n++) {
300  m_ptrsSurfPhase[n]->setConcentrations(CSolnSP + loc);
301  loc += m_nSpeciesSurfPhase[n];
302  }
303 }
304 
305 void solveSP::updateMFSolnSP(doublereal* XMolSolnSP)
306 {
307  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
308  size_t keqnStart = m_eqnIndexStartSolnPhase[isp];
309  m_ptrsSurfPhase[isp]->getMoleFractions(XMolSolnSP + keqnStart);
310  }
311 }
312 
313 void solveSP::updateMFKinSpecies(doublereal* XMolKinSpecies, int isp)
314 {
315  InterfaceKinetics* kin = m_objects[isp];
316  for (size_t iph = 0; iph < kin->nPhases(); iph++) {
317  size_t ksi = kin->kineticsSpeciesIndex(0, iph);
318  kin->thermo(iph).getMoleFractions(XMolKinSpecies + ksi);
319  }
320 }
321 
322 void solveSP::evalSurfLarge(const doublereal* CSolnSP)
323 {
324  size_t kindexSP = 0;
325  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
326  doublereal Clarge = CSolnSP[kindexSP];
327  m_spSurfLarge[isp] = 0;
328  kindexSP++;
329  for (size_t k = 1; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) {
330  if (CSolnSP[kindexSP] > Clarge) {
331  Clarge = CSolnSP[kindexSP];
332  m_spSurfLarge[isp] = k;
333  }
334  }
335  }
336 }
337 
338 void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln,
339  const doublereal* CSolnOld, const bool do_time,
340  const doublereal deltaT)
341 {
342  size_t k;
343  doublereal lenScale = 1.0E-9;
344  if (m_numSurfPhases > 0) {
345  // update the surface concentrations with the input surface
346  // concentration vector
347  updateState(CSoln);
348 
349  // Get the net production rates of all of the species in the
350  // surface kinetics mechanism
351  //
352  // HKM Should do it here for all kinetics objects so that
353  // bulk will eventually work.
354  if (do_time) {
355  size_t kindexSP = 0;
356  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
357  size_t nsp = m_nSpeciesSurfPhase[isp];
358  InterfaceKinetics* kinPtr = m_objects[isp];
359  size_t surfIndex = kinPtr->surfacePhaseIndex();
360  size_t kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex);
361  size_t kins = kindexSP;
363  for (k = 0; k < nsp; k++, kindexSP++) {
364  resid[kindexSP] =
365  (CSoln[kindexSP] - CSolnOld[kindexSP]) / deltaT
366  - m_netProductionRatesSave[kstart + k];
367  }
368 
369  size_t kspecial = kins + m_spSurfLarge[isp];
370  double sd = m_ptrsSurfPhase[isp]->siteDensity();
371  resid[kspecial] = sd;
372  for (k = 0; k < nsp; k++) {
373  resid[kspecial] -= CSoln[kins + k];
374  }
375  }
376  } else {
377  size_t kindexSP = 0;
378  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
379  size_t nsp = m_nSpeciesSurfPhase[isp];
380  InterfaceKinetics* kinPtr = m_objects[isp];
381  size_t surfIndex = kinPtr->surfacePhaseIndex();
382  size_t kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex);
383  size_t kins = kindexSP;
385  for (k = 0; k < nsp; k++, kindexSP++) {
386  resid[kindexSP] = - m_netProductionRatesSave[kstart + k];
387  }
388  size_t kspecial = kins + m_spSurfLarge[isp];
389  double sd = m_ptrsSurfPhase[isp]->siteDensity();
390  resid[kspecial] = sd;
391  for (k = 0; k < nsp; k++) {
392  resid[kspecial] -= CSoln[kins + k];
393  }
394  }
395  }
396 
397  if (m_bulkFunc == BULK_DEPOSITION) {
398  size_t kindexSP = m_numTotSurfSpecies;
399  for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) {
400  doublereal* XBlk = m_numEqn1.data();
401  size_t nsp = m_nSpeciesSurfPhase[isp];
402  size_t surfPhaseIndex = m_indexKinObjSurfPhase[isp];
403  InterfaceKinetics* kin = m_objects[isp];
404  double grRate = 0.0;
405  size_t kstart = kin->kineticsSpeciesIndex(0, surfPhaseIndex);
406  for (k = 0; k < nsp; k++) {
407  if (m_netProductionRatesSave[kstart + k] > 0.0) {
408  grRate += m_netProductionRatesSave[kstart + k];
409  }
410  }
411  resid[kindexSP] = m_bulkPhasePtrs[isp]->molarDensity();
412  for (k = 0; k < nsp; k++) {
413  resid[kindexSP] -= CSoln[kindexSP + k];
414  }
415  if (grRate > 0.0) {
416  for (k = 1; k < nsp; k++) {
417  if (m_netProductionRatesSave[kstart + k] > 0.0) {
418  resid[kindexSP + k] = XBlk[k] * grRate
419  - m_netProductionRatesSave[kstart + k];
420  } else {
421  resid[kindexSP + k] = XBlk[k] * grRate;
422  }
423  }
424  } else {
425  grRate = 1.0E-6;
426  //! @todo the appearance of k in this formula is suspicious
427  grRate += fabs(m_netProductionRatesSave[kstart + k]);
428  for (k = 1; k < nsp; k++) {
429  resid[kindexSP + k] = grRate * (XBlk[k] - 1.0/nsp);
430  }
431  }
432  if (do_time) {
433  for (k = 1; k < nsp; k++) {
434  resid[kindexSP + k] +=
435  lenScale / deltaT *
436  (CSoln[kindexSP + k]- CSolnOld[kindexSP + k]);
437  }
438  }
439  kindexSP += nsp;
440  }
441  }
442  }
443 }
444 
446  doublereal resid[], doublereal CSoln[],
447  const doublereal CSolnOld[], const bool do_time,
448  const doublereal deltaT)
449 {
450  size_t kColIndex = 0;
451  // Calculate the residual
452  fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
453  // Now we will look over the columns perturbing each unknown.
454  for (size_t jsp = 0; jsp < m_numSurfPhases; jsp++) {
455  size_t nsp = m_nSpeciesSurfPhase[jsp];
456  double sd = m_ptrsSurfPhase[jsp]->siteDensity();
457  for (size_t kCol = 0; kCol < nsp; kCol++) {
458  double cSave = CSoln[kColIndex];
459  double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
460  CSoln[kColIndex] += dc;
461  fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT);
462  for (size_t i = 0; i < m_neq; i++) {
463  jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc;
464  }
465  CSoln[kColIndex] = cSave;
466  kColIndex++;
467  }
468  }
469 
470  if (m_bulkFunc == BULK_DEPOSITION) {
471  for (size_t jsp = 0; jsp < m_numBulkPhasesSS; jsp++) {
472  size_t nsp = m_numBulkSpecies[jsp];
473  double sd = m_bulkPhasePtrs[jsp]->molarDensity();
474  for (size_t kCol = 0; kCol < nsp; kCol++) {
475  double cSave = CSoln[kColIndex];
476  double dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
477  CSoln[kColIndex] += dc;
478  fun_eval(m_numEqn2.data(), CSoln, CSolnOld, do_time, deltaT);
479  for (size_t i = 0; i < m_neq; i++) {
480  jac(i, kColIndex) = (m_numEqn2[i] - resid[i])/dc;
481  }
482  CSoln[kColIndex] = cSave;
483  kColIndex++;
484  }
485  }
486  }
487 }
488 
489 /*!
490  * This function calculates a damping factor for the Newton iteration update
491  * vector, dxneg, to insure that all site and bulk fractions, x, remain
492  * bounded between zero and one.
493  *
494  * dxneg[] = negative of the update vector.
495  *
496  * The constant "APPROACH" sets the fraction of the distance to the boundary
497  * that the step can take. If the full step would not force any fraction
498  * outside of 0-1, then Newton's method is allowed to operate normally.
499  */
500 static doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, int* label)
501 {
502  const doublereal APPROACH = 0.80;
503  doublereal damp = 1.0;
504  static doublereal damp_old = 1.0; //! @todo this variable breaks thread safety
505  *label = -1;
506 
507  for (size_t i = 0; i < dim; i++) {
508  // Calculate the new suggested new value of x[i]
509  double xnew = x[i] - damp * dxneg[i];
510 
511  // Calculate the allowed maximum and minimum values of x[i]
512  // - Only going to allow x[i] to converge to zero by a
513  // single order of magnitude at a time
514  double xtop = 1.0 - 0.1*fabs(1.0-x[i]);
515  double xbot = fabs(x[i]*0.1) - 1.0e-16;
516  if (xnew > xtop) {
517  damp = - APPROACH * (1.0 - x[i]) / dxneg[i];
518  *label = int(i);
519  } else if (xnew < xbot) {
520  damp = APPROACH * x[i] / dxneg[i];
521  *label = int(i);
522  } else if (xnew > 3.0*std::max(x[i], 1.0E-10)) {
523  damp = - 2.0 * std::max(x[i], 1.0E-10) / dxneg[i];
524  *label = int(i);
525  }
526  }
527  damp = std::max(damp, 1e-2);
528 
529  // Only allow the damping parameter to increase by a factor of three each
530  // iteration. Heuristic to avoid oscillations in the value of damp
531  if (damp > damp_old*3) {
532  damp = damp_old*3;
533  *label = -1;
534  }
535 
536  // Save old value of the damping parameter for use in subsequent calls.
537  damp_old = damp;
538  return damp;
539 
540 } /* calc_damping */
541 
542 /*
543  * This function calculates the norm of an update, dx[], based on the
544  * weighted values of x.
545  */
546 static doublereal calcWeightedNorm(const doublereal wtX[], const doublereal dx[], size_t dim)
547 {
548  doublereal norm = 0.0;
549  if (dim == 0) {
550  return 0.0;
551  }
552  for (size_t i = 0; i < dim; i++) {
553  norm += pow(dx[i] / wtX[i], 2);
554  }
555  return sqrt(norm/dim);
556 }
557 
558 void solveSP::calcWeights(doublereal wtSpecies[], doublereal wtResid[],
559  const Array2D& Jac, const doublereal CSoln[],
560  const doublereal abstol, const doublereal reltol)
561 {
562  // First calculate the weighting factor for the concentrations of the
563  // surface species and bulk species.
564  size_t kindex = 0;
565  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
566  double sd = m_ptrsSurfPhase[isp]->siteDensity();
567  for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindex++) {
568  wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
569  }
570  }
571  if (m_bulkFunc == BULK_DEPOSITION) {
572  for (size_t isp = 0; isp < m_numBulkPhasesSS; isp++) {
573  double sd = m_bulkPhasePtrs[isp]->molarDensity();
574  for (size_t k = 0; k < m_numBulkSpecies[isp]; k++, kindex++) {
575  wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
576  }
577  }
578  }
579 
580  // Now do the residual Weights. Since we have the Jacobian, we will use it
581  // to generate a number based on the what a significant change in a solution
582  // variable does to each residual. This is a row sum scale operation.
583  for (size_t k = 0; k < m_neq; k++) {
584  wtResid[k] = 0.0;
585  for (size_t jcol = 0; jcol < m_neq; jcol++) {
586  wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]);
587  }
588  }
589 }
590 
591 doublereal solveSP::calc_t(doublereal netProdRateSolnSP[],
592  doublereal XMolSolnSP[],
593  int* label, int* label_old,
594  doublereal* label_factor, int ioflag)
595 {
596  doublereal inv_timeScale = 1.0E-10;
597  size_t kindexSP = 0;
598  *label = 0;
599  updateMFSolnSP(XMolSolnSP);
600  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
601  // Get the interface kinetics associated with this surface
602  InterfaceKinetics* kin = m_objects[isp];
603 
604  // Calculate the start of the species index for surfaces within
605  // the InterfaceKinetics object
606  size_t surfIndex = kin->surfacePhaseIndex();
607  size_t kstart = kin->kineticsSpeciesIndex(0, surfIndex);
608  kin->getNetProductionRates(m_numEqn1.data());
609  double sden = kin->thermo(surfIndex).molarDensity();
610  for (size_t k = 0; k < m_nSpeciesSurfPhase[isp]; k++, kindexSP++) {
611  size_t kspindex = kstart + k;
612  netProdRateSolnSP[kindexSP] = m_numEqn1[kspindex];
613  double tmp = std::max(XMolSolnSP[kindexSP], 1.0e-10);
614  tmp *= sden;
615  tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp);
616  if (netProdRateSolnSP[kindexSP]> 0.0) {
617  tmp /= 100.;
618  }
619  if (tmp > inv_timeScale) {
620  inv_timeScale = tmp;
621  *label = int(kindexSP);
622  }
623  }
624  }
625 
626  // Increase time step exponentially as same species repeatedly controls time
627  // step
628  if (*label == *label_old) {
629  *label_factor *= 1.5;
630  } else {
631  *label_old = *label;
632  *label_factor = 1.0;
633  }
634  return inv_timeScale / *label_factor;
635 } // calc_t
636 
637 void solveSP::print_header(int ioflag, int ifunc, doublereal time_scale,
638  int damping, doublereal reltol, doublereal abstol)
639 {
640  if (ioflag) {
641  writelog("\n================================ SOLVESP CALL SETUP "
642  "========================================\n");
643  if (ifunc == SFLUX_INITIALIZE) {
644  writelog("\n SOLVESP Called with Initialization turned on\n");
645  writelogf(" Time scale input = %9.3e\n", time_scale);
646  } else if (ifunc == SFLUX_RESIDUAL) {
647  writelog("\n SOLVESP Called to calculate steady state residual\n");
648  writelog(" from a good initial guess\n");
649  } else if (ifunc == SFLUX_JACOBIAN) {
650  writelog("\n SOLVESP Called to calculate steady state Jacobian\n");
651  writelog(" from a good initial guess\n");
652  } else if (ifunc == SFLUX_TRANSIENT) {
653  writelog("\n SOLVESP Called to integrate surface in time\n");
654  writelogf(" for a total of %9.3e sec\n", time_scale);
655  } else {
656  throw CanteraError("solveSP::print_header",
657  "Unknown ifunc flag = {}", ifunc);
658  }
659 
660  if (m_bulkFunc == BULK_DEPOSITION) {
661  writelog(" The composition of the Bulk Phases will be calculated\n");
662  } else if (m_bulkFunc == BULK_ETCH) {
663  writelog(" Bulk Phases have fixed compositions\n");
664  } else {
665  throw CanteraError("solveSP::print_header",
666  "Unknown bulkFunc flag = {}", m_bulkFunc);
667  }
668 
669  if (damping) {
670  writelog(" Damping is ON \n");
671  } else {
672  writelog(" Damping is OFF \n");
673  }
674 
675  writelogf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol);
676  }
677 
678  if (ioflag == 1) {
679  writelog("\n\n\t Iter Time Del_t Damp DelX "
680  " Resid Name-Time Name-Damp\n");
681  writelog("\t -----------------------------------------------"
682  "------------------------------------\n");
683  }
684 }
685 
686 void solveSP::printIteration(int ioflag, doublereal damp, int label_d,
687  int label_t, doublereal inv_t, doublereal t_real,
688  size_t iter, doublereal update_norm,
689  doublereal resid_norm, bool do_time, bool final)
690 {
691  if (ioflag == 1) {
692  if (final) {
693  writelogf("\tFIN%3d ", iter);
694  } else {
695  writelogf("\t%6d ", iter);
696  }
697  if (do_time) {
698  writelogf("%9.4e %9.4e ", t_real, 1.0/inv_t);
699  } else {
700  writeline(' ', 22, false);
701  }
702  if (damp < 1.0) {
703  writelogf("%9.4e ", damp);
704  } else {
705  writeline(' ', 11, false);
706  }
707  writelogf("%9.4e %9.4e", update_norm, resid_norm);
708  if (do_time) {
709  size_t k = m_kinSpecIndex[label_t];
710  size_t isp = m_kinObjIndex[label_t];
711  writelog(" %-16s", m_objects[isp]->kineticsSpeciesName(k));
712  } else {
713  writeline(' ', 16, false);
714  }
715  if (label_d >= 0) {
716  size_t k = m_kinSpecIndex[label_d];
717  size_t isp = m_kinObjIndex[label_d];
718  writelogf(" %-16s", m_objects[isp]->kineticsSpeciesName(k));
719  }
720  if (final) {
721  writelog(" -- success");
722  }
723  writelog("\n");
724  }
725 } // printIteration
726 
727 }
Declarations for the implicit integration of surface site density equations (see Kinetics Managers an...
Header for a simple thermodynamics model of a surface phase derived from ThermoPhase,...
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:32
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition: DenseMatrix.h:51
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the matrix.
Definition: DenseMatrix.cpp:69
Advances the surface coverages of the associated set of SurfacePhase objects in time.
A kinetics manager for heterogeneous reaction mechanisms.
thermo_t & thermo(size_t n=0)
This method returns a reference to the nth ThermoPhase object defined in this kinetics mechanism.
Definition: Kinetics.h:227
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:167
virtual void getNetProductionRates(doublereal *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition: Kinetics.cpp:433
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:261
size_t surfacePhaseIndex() const
This returns the integer index of the phase which has ThermoPhase type cSurf.
Definition: Kinetics.h:202
double molarDensity() const
Molar density (kmol/m^3).
Definition: Phase.cpp:700
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:285
void getMoleFractions(double *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:572
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition: SurfPhase.h:143
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:188
void fun_eval(doublereal *resid, const doublereal *CSolnSP, const doublereal *CSolnOldSP, const bool do_time, const doublereal deltaT)
Main Function evaluation.
Definition: solveSP.cpp:338
void print_header(int ioflag, int ifunc, doublereal time_scale, int damping, doublereal reltol, doublereal abstol)
Printing routine that optionally gets called at the start of every invocation.
Definition: solveSP.cpp:637
std::vector< size_t > m_indexKinObjSurfPhase
Mapping between the surface phases and the InterfaceKinetics objects.
Definition: solveSP.h:352
size_t m_numTotSurfSpecies
Total number of surface species in all surface phases.
Definition: solveSP.h:344
void updateMFSolnSP(doublereal *XMolSolnSP)
Update mole fraction vector consisting of unknowns in surface problem.
Definition: solveSP.cpp:305
std::vector< size_t > m_kinObjIndex
Index between the equation index and the index of the InterfaceKinetics object.
Definition: solveSP.h:433
vector_fp m_netProductionRatesSave
Temporary vector with length equal to max m_maxTotSpecies.
Definition: solveSP.h:458
void evalSurfLarge(const doublereal *CSolnSP)
Update the vector that keeps track of the largest species in each surface phase.
Definition: solveSP.cpp:322
void resjac_eval(DenseMatrix &jac, doublereal *resid, doublereal *CSolnSP, const doublereal *CSolnSPOld, const bool do_time, const doublereal deltaT)
Main routine that calculates the current residual and Jacobian.
Definition: solveSP.cpp:445
void calcWeights(doublereal wtSpecies[], doublereal wtResid[], const Array2D &Jac, const doublereal CSolnSP[], const doublereal abstol, const doublereal reltol)
Calculate the solution and residual weights.
Definition: solveSP.cpp:558
void printIteration(int ioflag, doublereal damp, int label_d, int label_t, doublereal inv_t, doublereal t_real, size_t iter, doublereal update_norm, doublereal resid_norm, bool do_time, bool final=false)
Printing routine that gets called after every iteration.
Definition: solveSP.cpp:686
int m_bulkFunc
This variable determines how the bulk phases are to be handled.
Definition: solveSP.h:329
std::vector< SurfPhase * > m_ptrsSurfPhase
Vector of surface phase pointers.
Definition: solveSP.h:366
vector_fp m_CSolnSP
Solution vector. length MAX(1, m_neq)
Definition: solveSP.h:470
vector_fp m_wtResid
Weights for the residual norm calculation. length MAX(1, m_neq)
Definition: solveSP.h:479
vector_fp m_resid
Residual for the surface problem.
Definition: solveSP.h:497
vector_fp m_CSolnSave
Temporary vector with length equal to max m_maxTotSpecies.
Definition: solveSP.h:467
vector_fp m_XMolKinSpecies
Vector of mole fractions. length m_maxTotSpecies.
Definition: solveSP.h:500
std::vector< ThermoPhase * > m_bulkPhasePtrs
Vector of bulk phase pointers, length is equal to m_numBulkPhases.
Definition: solveSP.h:416
std::vector< size_t > m_numBulkSpecies
Vector of number of species in the m_numBulkPhases phases.
Definition: solveSP.h:406
std::vector< size_t > m_kinSpecIndex
Index between the equation index and the position in the kinetic species array for the appropriate ki...
Definition: solveSP.h:426
doublereal calc_t(doublereal netProdRateSolnSP[], doublereal XMolSolnSP[], int *label, int *label_old, doublereal *label_factor, int ioflag)
Calculate a conservative delta T to use in a pseudo-steady state algorithm.
Definition: solveSP.cpp:591
std::vector< size_t > m_nSpeciesSurfPhase
Vector of length number of surface phases containing the number of surface species in each phase.
Definition: solveSP.h:359
vector_fp m_wtSpecies
Weights for the species concentrations norm calculation.
Definition: solveSP.h:485
const int BULK_ETCH
Etching of a bulk phase is to be expected.
Definition: solveSP.h:56
void updateMFKinSpecies(doublereal *XMolKinSp, int isp)
Update the mole fraction vector for a specific kinetic species vector corresponding to one InterfaceK...
Definition: solveSP.cpp:313
size_t m_numTotBulkSpeciesSS
Total number of species in all bulk phases.
Definition: solveSP.h:413
const int BULK_DEPOSITION
Deposition of a bulk phase is to be expected.
Definition: solveSP.h:51
std::vector< size_t > m_spSurfLarge
Vector containing the indices of the largest species in each surface phase.
Definition: solveSP.h:442
int solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, doublereal PGas, doublereal reltol, doublereal abstol)
Main routine that actually calculates the pseudo steady state of the surface problem.
Definition: solveSP.cpp:113
vector_fp m_numEqn2
Temporary vector with length equal to max m_maxTotSpecies.
Definition: solveSP.h:464
size_t m_numSurfPhases
Number of surface phases in the surface problem.
Definition: solveSP.h:336
std::vector< size_t > m_kinObjPhaseIDSurfPhase
Phase ID in the InterfaceKinetics object of the surface phase.
Definition: solveSP.h:387
vector_fp m_numEqn1
Temporary vector with length equal to max m_maxTotSpecies.
Definition: solveSP.h:461
DenseMatrix m_Jac
Jacobian.
Definition: solveSP.h:504
size_t m_numBulkPhasesSS
Total number of volumetric condensed phases included in the steady state problem handled by this rout...
Definition: solveSP.h:400
size_t m_neq
Total number of equations to solve in the implicit problem.
Definition: solveSP.h:323
std::vector< size_t > m_eqnIndexStartSolnPhase
Index of the start of the unknowns for each solution phase.
Definition: solveSP.h:378
std::vector< InterfaceKinetics * > & m_objects
Vector of interface kinetics objects.
Definition: solveSP.h:317
void updateState(const doublereal *cSurfSpec)
Update the surface states of the surface phases.
Definition: solveSP.cpp:296
size_t m_maxTotSpecies
Maximum number of species in any single kinetics operator -> also maxed wrt the total # of solution s...
Definition: solveSP.h:455
vector_fp m_CSolnSPOld
Saved solution vector at the old time step. length MAX(1, m_neq)
Definition: solveSP.h:476
vector_fp m_CSolnSPInit
Saved initial solution vector. length MAX(1, m_neq)
Definition: solveSP.h:473
const int SFLUX_TRANSIENT
The transient calculation is performed here for an amount of time specified by "time_scale".
Definition: solveSP.h:41
const int SFLUX_RESIDUAL
Need to solve the surface problem in order to calculate the surface fluxes of gas-phase species.
Definition: solveSP.h:28
const int SFLUX_JACOBIAN
Calculation of the surface problem is due to the need for a numerical Jacobian for the gas-problem.
Definition: solveSP.h:34
const int SFLUX_INITIALIZE
This assumes that the initial guess supplied to the routine is far from the correct one.
Definition: solveSP.h:21
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:158
void writelogf(const char *fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:179
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).