Cantera  2.1.2
solveSP.cpp
1 /*
2  * @file: solveSP.cpp Implicit surface site concentration solver
3  */
4 /*
5  * Copyright 2004 Sandia Corporation. Under the terms of Contract
6  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
7  * retains certain rights in this software.
8  * See file License.txt for licensing information.
9  */
10
12 #include "cantera/base/clockWC.h"
14
15 #include <cstdio>
16
17 using namespace std;
18 namespace Cantera
19 {
20
21 /***************************************************************************
22  * STATIC ROUTINES DEFINED IN THIS FILE
23  ***************************************************************************/
24
25 static doublereal calc_damping(doublereal* x, doublereal* dx, size_t dim, int*);
26 static doublereal calcWeightedNorm(const doublereal [], const doublereal dx[], size_t);
27
28 /***************************************************************************
29  * solveSP Class Definitions
30  ***************************************************************************/
31
32 solveSP::solveSP(ImplicitSurfChem* surfChemPtr, int bulkFunc) :
33  m_SurfChemPtr(surfChemPtr),
34  m_objects(surfChemPtr->getObjects()),
35  m_neq(0),
36  m_bulkFunc(bulkFunc),
37  m_numSurfPhases(0),
38  m_numTotSurfSpecies(0),
39  m_numBulkPhasesSS(0),
40  m_numTotBulkSpeciesSS(0),
41  m_atol(1.0E-15),
42  m_rtol(1.0E-4),
43  m_maxstep(1000),
44  m_maxTotSpecies(0),
45  m_ioflag(0)
46 {
47  m_numSurfPhases = 0;
48  size_t numPossibleSurfPhases = m_objects.size();
49  for (size_t n = 0; n < numPossibleSurfPhases; n++) {
50  InterfaceKinetics* m_kin = m_objects[n];
51  size_t surfPhaseIndex = m_kin->surfacePhaseIndex();
52  if (surfPhaseIndex != npos) {
54  m_indexKinObjSurfPhase.push_back(n);
55  m_kinObjPhaseIDSurfPhase.push_back(surfPhaseIndex);
56  } else {
57  throw CanteraError("solveSP",
58  "InterfaceKinetics object has no surface phase");
59  }
60  ThermoPhase* tp = &(m_kin->thermo(surfPhaseIndex));
61  SurfPhase* sp = dynamic_cast<SurfPhase*>(tp);
62  if (!sp) {
63  throw CanteraError("solveSP",
64  "Inconsistent ThermoPhase object within "
65  "InterfaceKinetics object");
66  }
67
68  m_ptrsSurfPhase.push_back(sp);
69  size_t nsp = sp->nSpecies();
70  m_nSpeciesSurfPhase.push_back(nsp);
71  m_numTotSurfSpecies += nsp;
72
73  }
74  /*
75  * We rely on ordering to figure things out
76  */
78
79  if (bulkFunc == BULK_DEPOSITION) {
81  } else {
83  }
84
85  m_maxTotSpecies = 0;
86  for (size_t n = 0; n < m_numSurfPhases; n++) {
87  size_t tsp = m_objects[n]->nTotalSpecies();
88  m_maxTotSpecies = std::max(m_maxTotSpecies, tsp);
89  }
91
92
94  m_numEqn1.resize(m_maxTotSpecies, 0.0);
95  m_numEqn2.resize(m_maxTotSpecies, 0.0);
96  m_XMolKinSpecies.resize(m_maxTotSpecies, 0.0);
97  m_CSolnSave.resize(m_neq, 0.0);
98
99  m_spSurfLarge.resize(m_numSurfPhases, 0);
100
103  m_eqnIndexStartSolnPhase.resize(m_numSurfPhases + m_numBulkPhasesSS, 0);
104
105  size_t kindexSP = 0;
106  size_t isp, k, nsp, kstart;
107  for (isp = 0; isp < m_numSurfPhases; isp++) {
108  size_t iKinObject = m_indexKinObjSurfPhase[isp];
109  InterfaceKinetics* m_kin = m_objects[iKinObject];
110  size_t surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp];
111  kstart = m_kin->kineticsSpeciesIndex(0, surfPhaseIndex);
112  nsp = m_nSpeciesSurfPhase[isp];
113  m_eqnIndexStartSolnPhase[isp] = kindexSP;
114  for (k = 0; k < nsp; k++, kindexSP++) {
115  m_kinSpecIndex[kindexSP] = kstart + k;
116  m_kinObjIndex[kindexSP] = isp;
117  }
118  }
119
120  // Dimension solution vector
121  size_t dim1 = std::max<size_t>(1, m_neq);
122  m_CSolnSP.resize(dim1, 0.0);
123  m_CSolnSPInit.resize(dim1, 0.0);
124  m_CSolnSPOld.resize(dim1, 0.0);
125  m_wtResid.resize(dim1, 0.0);
126  m_wtSpecies.resize(dim1, 0.0);
127  m_resid.resize(dim1, 0.0);
128  m_ipiv.resize(dim1, 0);
129
130  m_Jac.resize(dim1, dim1, 0.0);
131  m_JacCol.resize(dim1, 0);
132  for (size_t k = 0; k < dim1; k++) {
133  m_JacCol[k] = m_Jac.ptrColumn(k);
134  }
135 }
136
138 {
139 }
140
141 int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin,
142  doublereal PGas, doublereal reltol, doublereal abstol)
143 {
144  doublereal EXTRA_ACCURACY = 0.001;
145  if (ifunc == SFLUX_JACOBIAN) {
146  EXTRA_ACCURACY *= 0.001;
147  }
148  int info = 0;
149  int label_t=-1; /* Species IDs for time control */
150  int label_d = -1; /* Species IDs for damping control */
151  int label_t_old=-1;
152  doublereal label_factor = 1.0;
153  int iter=0; // iteration number on numlinear solver
154  int iter_max=1000; // maximum number of nonlinear iterations
155  int nrhs=1;
156  doublereal deltaT = 1.0E-10; // Delta time step
157  doublereal damp=1.0, tmp;
158  // Weighted L2 norm of the residual. Currently, this is only
159  // used for IO purposes. It doesn't control convergence.
160  doublereal resid_norm;
161  doublereal inv_t = 0.0;
162  doublereal t_real = 0.0, update_norm = 1.0E6;
163
164  bool do_time = false, not_converged = true;
165
166  if (m_ioflag > 1) {
167  m_ioflag = 1;
168  }
169
170  /*
171  * Set the initial value of the do_time parameter
172  */
173  if (ifunc == SFLUX_INITIALIZE || ifunc == SFLUX_TRANSIENT) {
174  do_time = true;
175  }
176
177  /*
178  * Store the initial guess for the surface problem in the soln vector,
179  * CSoln, and in an separate vector CSolnInit.
180  */
181  size_t loc = 0;
182  for (size_t n = 0; n < m_numSurfPhases; n++) {
183  SurfPhase* sf_ptr = m_ptrsSurfPhase[n];
185  size_t nsp = m_nSpeciesSurfPhase[n];
186  for (size_t k = 0; k <nsp; k++) {
187  m_CSolnSP[loc] = m_numEqn1[k];
188  loc++;
189  }
190  }
191
192  std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPInit.begin());
193
194  // Calculate the largest species in each phase
196  /*
197  * Get the net production rate of all species in the kinetics manager.
198  */
199  // m_kin->getNetProductionRates(DATA_PTR(m_netProductionRatesSave));
200
201  if (m_ioflag) {
202  print_header(m_ioflag, ifunc, time_scale, true, reltol, abstol,
203  TKelvin, PGas, DATA_PTR(m_netProductionRatesSave),
205  }
206
207  /*
208  * Quick return when there isn't a surface problem to solve
209  */
210  if (m_neq == 0) {
211  not_converged = false;
212  update_norm = 0.0;
213  }
214
215  /* ------------------------------------------------------------------
216  * Start of Newton's method
217  * ------------------------------------------------------------------
218  */
219  while (not_converged && iter < iter_max) {
220  iter++;
221  /*
222  * Store previous iteration's solution in the old solution vector
223  */
224  std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPOld.begin());
225
226  /*
227  * Evaluate the largest surface species for each surface phase every
228  * 5 iterations.
229  */
230  if (iter%5 == 4) {
232  }
233
234  /*
235  * Calculate the value of the time step
236  * - heuristics to stop large oscillations in deltaT
237  */
238  if (do_time) {
239  /* don't hurry increase in time step at the same time as damping */
240  if (damp < 1.0) {
241  label_factor = 1.0;
242  }
245  &label_t, &label_t_old, &label_factor, m_ioflag);
246  if (iter < 10) {
247  inv_t = tmp;
248  } else if (tmp > 2.0*inv_t) {
249  inv_t = 2.0*inv_t;
250  } else {
251  inv_t = tmp;
252  }
253
254  /*
255  * Check end condition
256  */
257
258  if (ifunc == SFLUX_TRANSIENT) {
259  tmp = t_real + 1.0/inv_t;
260  if (tmp > time_scale) {
261  inv_t = 1.0/(time_scale - t_real);
262  }
263  }
264  } else {
265  /* make steady state calc a step of 1 million seconds to
266  prevent singular jacobians for some pathological cases */
267  inv_t = 1.0e-6;
268  }
269  deltaT = 1.0/inv_t;
270
271  /*
272  * Call the routine to numerically evaluation the jacobian
273  * and residual for the current iteration.
274  */
276  DATA_PTR(m_CSolnSPOld), do_time, deltaT);
277
278  /*
279  * Calculate the weights. Make sure the calculation is carried
280  * out on the first iteration.
281  */
282  if (iter%4 == 1) {
284  m_Jac, DATA_PTR(m_CSolnSP), abstol, reltol);
285  }
286
287  /*
288  * Find the weighted norm of the residual
289  */
290  resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid),
292
293  /*
294  * Solve Linear system (with LAPACK). The solution is in resid[]
295  */
296  ct_dgetrf(m_neq, m_neq, m_JacCol[0], m_neq, DATA_PTR(m_ipiv), info);
297  if (info==0) {
298  ct_dgetrs(ctlapack::NoTranspose, m_neq, nrhs, m_JacCol[0],
300  info);
301  }
302  /*
303  * Force convergence if residual is small to avoid
304  * "nan" results from the linear solve.
305  */
306  else {
307  if (m_ioflag) {
308  printf("solveSurfSS: Zero pivot, assuming converged: %g (%d)\n",
309  resid_norm, info);
310  }
311  for (size_t jcol = 0; jcol < m_neq; jcol++) {
312  m_resid[jcol] = 0.0;
313  }
314
315  /* print out some helpful info */
316  if (m_ioflag > 1) {
317  printf("-----\n");
318  printf("solveSurfProb: iter %d t_real %g delta_t %g\n\n",
319  iter,t_real, 1.0/inv_t);
320  printf("solveSurfProb: init guess, current concentration,"
321  "and prod rate:\n");
322  for (size_t jcol = 0; jcol < m_neq; jcol++) {
323  printf("\t%s %g %g %g\n", int2str(jcol).c_str(),
324  m_CSolnSPInit[jcol], m_CSolnSP[jcol],
326  }
327  printf("-----\n");
328  }
329  if (do_time) {
330  t_real += time_scale;
331  }
332  }
333
334  /*
335  * Calculate the Damping factor needed to keep all unknowns
336  * between 0 and 1, and not allow too large a change (factor of 2)
337  * in any unknown.
338  */
339
340  damp = calc_damping(DATA_PTR(m_CSolnSP), DATA_PTR(m_resid), m_neq, &label_d);
341
342  /*
343  * Calculate the weighted norm of the update vector
344  * Here, resid is the delta of the solution, in concentration
345  * units.
346  */
347  update_norm = calcWeightedNorm(DATA_PTR(m_wtSpecies),
349  /*
350  * Update the solution vector and real time
351  * Crop the concentrations to zero.
352  */
353  for (size_t irow = 0; irow < m_neq; irow++) {
354  m_CSolnSP[irow] -= damp * m_resid[irow];
355  }
356  for (size_t irow = 0; irow < m_neq; irow++) {
357  m_CSolnSP[irow] = std::max(0.0, m_CSolnSP[irow]);
358  }
360
361  if (do_time) {
362  t_real += damp/inv_t;
363  }
364
365  if (m_ioflag) {
366  printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
367  update_norm, resid_norm,
371  m_neq, do_time);
372  }
373
374  if (ifunc == SFLUX_TRANSIENT) {
375  not_converged = (t_real < time_scale);
376  } else {
377  if (do_time) {
378  if (t_real > time_scale ||
379  (resid_norm < 1.0e-7 &&
380  update_norm*time_scale/t_real < EXTRA_ACCURACY)) {
381  do_time = false;
382  }
383  } else {
384  not_converged = ((update_norm > EXTRA_ACCURACY) ||
385  (resid_norm > EXTRA_ACCURACY));
386  }
387  }
388  } /* End of Newton's Method while statement */
389
390  /*
391  * End Newton's method. If not converged, print error message and
392  * recalculate sdot's at equal site fractions.
393  */
394  if (not_converged) {
395  if (m_ioflag) {
396  printf("#$#$#$# Error in solveSP$#$#$#\$ \n");
397  printf("Newton iter on surface species did not converge, "
398  "update_norm = %e \n", update_norm);
399  printf("Continuing anyway\n");
400  }
401  }
402
403  /*
404  * Decide on what to return in the solution vector
405  * - right now, will always return the last solution
406  * no matter how bad
407  */
408  if (m_ioflag) {
410  false, deltaT);
411  resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid),
413  printFinal(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
414  update_norm, resid_norm, DATA_PTR(m_netProductionRatesSave),
417  DATA_PTR(m_wtResid), m_neq, do_time,
418  TKelvin, PGas);
419  }
420
421  /*
422  * Return with the appropriate flag
423  */
424  if (update_norm > 1.0) {
425  return -1;
426  }
427  return 1;
428 }
429
430 void solveSP::updateState(const doublereal* CSolnSP)
431 {
432  size_t loc = 0;
433  for (size_t n = 0; n < m_numSurfPhases; n++) {
434  m_ptrsSurfPhase[n]->setConcentrations(CSolnSP + loc);
435  loc += m_nSpeciesSurfPhase[n];
436  }
437 }
438
439 void solveSP::updateMFSolnSP(doublereal* XMolSolnSP)
440 {
441  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
442  size_t keqnStart = m_eqnIndexStartSolnPhase[isp];
443  m_ptrsSurfPhase[isp]->getMoleFractions(XMolSolnSP + keqnStart);
444  }
445 }
446
447 void solveSP::updateMFKinSpecies(doublereal* XMolKinSpecies, int isp)
448 {
449  InterfaceKinetics* m_kin = m_objects[isp];
450  size_t nph = m_kin->nPhases();
451  for (size_t iph = 0; iph < nph; iph++) {
452  size_t ksi = m_kin->kineticsSpeciesIndex(0, iph);
453  ThermoPhase& thref = m_kin->thermo(iph);
454  thref.getMoleFractions(XMolKinSpecies + ksi);
455  }
456 }
457
458 void solveSP::evalSurfLarge(const doublereal* CSolnSP)
459 {
460  size_t kindexSP = 0;
461  for (size_t isp = 0; isp < m_numSurfPhases; isp++) {
462  size_t nsp = m_nSpeciesSurfPhase[isp];
463  doublereal Clarge = CSolnSP[kindexSP];
464  m_spSurfLarge[isp] = 0;
465  kindexSP++;
466  for (size_t k = 1; k < nsp; k++, kindexSP++) {
467  if (CSolnSP[kindexSP] > Clarge) {
468  Clarge = CSolnSP[kindexSP];
469  m_spSurfLarge[isp] = k;
470  }
471  }
472  }
473 }
474
475 void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln,
476  const doublereal* CSolnOld, const bool do_time,
477  const doublereal deltaT)
478 {
479  size_t isp, nsp, kstart, k, kindexSP, kins, kspecial;
480  doublereal lenScale = 1.0E-9;
481  doublereal sd = 0.0;
482  doublereal grRate;
483  if (m_numSurfPhases > 0) {
484  /*
485  * update the surface concentrations with the input surface
486  * concentration vector
487  */
488  updateState(CSoln);
489  /*
490  * Get the net production rates of all of the species in the
491  * surface kinetics mechanism
492  *
493  * HKM Should do it here for all kinetics objects so that
494  * bulk will eventually work.
495  */
496
497  if (do_time) {
498  kindexSP = 0;
499  for (isp = 0; isp < m_numSurfPhases; isp++) {
500  nsp = m_nSpeciesSurfPhase[isp];
501  InterfaceKinetics* kinPtr = m_objects[isp];
502  size_t surfIndex = kinPtr->surfacePhaseIndex();
503  kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex);
504  kins = kindexSP;
506  for (k = 0; k < nsp; k++, kindexSP++) {
507  resid[kindexSP] =
508  (CSoln[kindexSP] - CSolnOld[kindexSP]) / deltaT
509  - m_netProductionRatesSave[kstart + k];
510  }
511
512  kspecial = kins + m_spSurfLarge[isp];
513  sd = m_ptrsSurfPhase[isp]->siteDensity();
514  resid[kspecial] = sd;
515  for (k = 0; k < nsp; k++) {
516  resid[kspecial] -= CSoln[kins + k];
517  }
518  }
519  } else {
520  kindexSP = 0;
521  for (isp = 0; isp < m_numSurfPhases; isp++) {
522  nsp = m_nSpeciesSurfPhase[isp];
523  InterfaceKinetics* kinPtr = m_objects[isp];
524  size_t surfIndex = kinPtr->surfacePhaseIndex();
525  kstart = kinPtr->kineticsSpeciesIndex(0, surfIndex);
526  kins = kindexSP;
528  for (k = 0; k < nsp; k++, kindexSP++) {
529  resid[kindexSP] = - m_netProductionRatesSave[kstart + k];
530  }
531  kspecial = kins + m_spSurfLarge[isp];
532  sd = m_ptrsSurfPhase[isp]->siteDensity();
533  resid[kspecial] = sd;
534  for (k = 0; k < nsp; k++) {
535  resid[kspecial] -= CSoln[kins + k];
536  }
537  }
538  }
539
540  if (m_bulkFunc == BULK_DEPOSITION) {
541  kindexSP = m_numTotSurfSpecies;
542  for (isp = 0; isp < m_numBulkPhasesSS; isp++) {
543  doublereal* XBlk = DATA_PTR(m_numEqn1);
544  //ThermoPhase *THptr = m_bulkPhasePtrs[isp];
545  //THptr->getMoleFractions(XBlk);
546  nsp = m_nSpeciesSurfPhase[isp];
547  size_t surfPhaseIndex = m_indexKinObjSurfPhase[isp];
548  InterfaceKinetics* m_kin = m_objects[isp];
549  grRate = 0.0;
550  kstart = m_kin->kineticsSpeciesIndex(0, surfPhaseIndex);
551  for (k = 0; k < nsp; k++) {
552  if (m_netProductionRatesSave[kstart + k] > 0.0) {
553  grRate += m_netProductionRatesSave[kstart + k];
554  }
555  }
556  resid[kindexSP] = m_bulkPhasePtrs[isp]->molarDensity();
557  for (k = 0; k < nsp; k++) {
558  resid[kindexSP] -= CSoln[kindexSP + k];
559  }
560  if (grRate > 0.0) {
561  for (k = 1; k < nsp; k++) {
562  if (m_netProductionRatesSave[kstart + k] > 0.0) {
563  resid[kindexSP + k] = XBlk[k] * grRate
564  - m_netProductionRatesSave[kstart + k];
565  } else {
566  resid[kindexSP + k] = XBlk[k] * grRate;
567  }
568  }
569  } else {
570  grRate = 1.0E-6;
571  grRate += fabs(m_netProductionRatesSave[kstart + k]);
572  for (k = 1; k < nsp; k++) {
573  resid[kindexSP + k] = grRate * (XBlk[k] - 1.0/nsp);
574  }
575  }
576  if (do_time) {
577  for (k = 1; k < nsp; k++) {
578  resid[kindexSP + k] +=
579  lenScale / deltaT *
580  (CSoln[kindexSP + k]- CSolnOld[kindexSP + k]);
581  }
582  }
583  kindexSP += nsp;
584  }
585  }
586  }
587 }
588
589 void solveSP::resjac_eval(std::vector<doublereal*> &JacCol,
590  doublereal resid[], doublereal CSoln[],
591  const doublereal CSolnOld[], const bool do_time,
592  const doublereal deltaT)
593 {
594  size_t kColIndex = 0, nsp, jsp, i, kCol;
595  doublereal dc, cSave, sd;
596  doublereal* col_j;
597  /*
598  * Calculate the residual
599  */
600  fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
601  /*
602  * Now we will look over the columns perturbing each unknown.
603  */
604  for (jsp = 0; jsp < m_numSurfPhases; jsp++) {
605  nsp = m_nSpeciesSurfPhase[jsp];
606  sd = m_ptrsSurfPhase[jsp]->siteDensity();
607  for (kCol = 0; kCol < nsp; kCol++) {
608  cSave = CSoln[kColIndex];
609  dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
610  CSoln[kColIndex] += dc;
611  fun_eval(DATA_PTR(m_numEqn2), CSoln, CSolnOld, do_time, deltaT);
612  col_j = JacCol[kColIndex];
613  for (i = 0; i < m_neq; i++) {
614  col_j[i] = (m_numEqn2[i] - resid[i])/dc;
615  }
616  CSoln[kColIndex] = cSave;
617  kColIndex++;
618  }
619  }
620
621  if (m_bulkFunc == BULK_DEPOSITION) {
622  for (jsp = 0; jsp < m_numBulkPhasesSS; jsp++) {
623  nsp = m_numBulkSpecies[jsp];
624  sd = m_bulkPhasePtrs[jsp]->molarDensity();
625  for (kCol = 0; kCol < nsp; kCol++) {
626  cSave = CSoln[kColIndex];
627  dc = std::max(1.0E-10 * sd, fabs(cSave) * 1.0E-7);
628  CSoln[kColIndex] += dc;
629  fun_eval(DATA_PTR(m_numEqn2), CSoln, CSolnOld, do_time, deltaT);
630  col_j = JacCol[kColIndex];
631  for (i = 0; i < m_neq; i++) {
632  col_j[i] = (m_numEqn2[i] - resid[i])/dc;
633  }
634  CSoln[kColIndex] = cSave;
635  kColIndex++;
636  }
637  }
638  }
639 }
640
641 /*!
642  * This function calculates a damping factor for the Newton iteration update
643  * vector, dxneg, to insure that all site and bulk fractions, x, remain
644  * bounded between zero and one.
645  *
646  * dxneg[] = negative of the update vector.
647  *
648  * The constant "APPROACH" sets the fraction of the distance to the boundary
649  * that the step can take. If the full step would not force any fraction
650  * outside of 0-1, then Newton's method is allowed to operate normally.
651  */
652 static doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, int* label)
653 {
654  const doublereal APPROACH = 0.80;
655  doublereal damp = 1.0, xnew, xtop, xbot;
656  static doublereal damp_old = 1.0;
657
658  *label = -1;
659
660  for (size_t i = 0; i < dim; i++) {
661
662  /*
663  * Calculate the new suggested new value of x[i]
664  */
665
666  xnew = x[i] - damp * dxneg[i];
667
668  /*
669  * Calculate the allowed maximum and minimum values of x[i]
670  * - Only going to allow x[i] to converge to zero by a
671  * single order of magnitude at a time
672  */
673
674  xtop = 1.0 - 0.1*fabs(1.0-x[i]);
675  xbot = fabs(x[i]*0.1) - 1.0e-16;
676  if (xnew > xtop) {
677  damp = - APPROACH * (1.0 - x[i]) / dxneg[i];
678  *label = int(i);
679  } else if (xnew < xbot) {
680  damp = APPROACH * x[i] / dxneg[i];
681  *label = int(i);
682  } else if (xnew > 3.0*std::max(x[i], 1.0E-10)) {
683  damp = - 2.0 * std::max(x[i], 1.0E-10) / dxneg[i];
684  *label = int(i);
685  }
686  }
687
688  if (damp < 1.0e-2) {
689  damp = 1.0e-2;
690  }
691  /*
692  * Only allow the damping parameter to increase by a factor of three each
693  * iteration. Heuristic to avoid oscillations in the value of damp
694  */
695
696  if (damp > damp_old*3) {
697  damp = damp_old*3;
698  *label = -1;
699  }
700
701  /*
702  * Save old value of the damping parameter for use
703  * in subsequent calls.
704  */
705
706  damp_old = damp;
707  return damp;
708
709 } /* calc_damping */
710
711 /*
712  * This function calculates the norm of an update, dx[],
713  * based on the weighted values of x.
714  */
715 static doublereal calcWeightedNorm(const doublereal wtX[], const doublereal dx[], size_t dim)
716 {
717  doublereal norm = 0.0;
718  doublereal tmp;
719  if (dim == 0) {
720  return 0.0;
721  }
722  for (size_t i = 0; i < dim; i++) {
723  tmp = dx[i] / wtX[i];
724  norm += tmp * tmp;
725  }
726  return sqrt(norm/dim);
727 }
728
729 void solveSP::calcWeights(doublereal wtSpecies[], doublereal wtResid[],
730  const Array2D& Jac, const doublereal CSoln[],
731  const doublereal abstol, const doublereal reltol)
732 {
733  size_t k, jcol, kindex, isp, nsp;
734  doublereal sd;
735  /*
736  * First calculate the weighting factor for the concentrations of
737  * the surface species and bulk species.
738  */
739  kindex = 0;
740  for (isp = 0; isp < m_numSurfPhases; isp++) {
741  nsp = m_nSpeciesSurfPhase[isp];
742  sd = m_ptrsSurfPhase[isp]->siteDensity();
743  for (k = 0; k < nsp; k++, kindex++) {
744  wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
745  }
746  }
747  if (m_bulkFunc == BULK_DEPOSITION) {
748  for (isp = 0; isp < m_numBulkPhasesSS; isp++) {
749  nsp = m_numBulkSpecies[isp];
750  sd = m_bulkPhasePtrs[isp]->molarDensity();
751  for (k = 0; k < nsp; k++, kindex++) {
752  wtSpecies[kindex] = abstol * sd + reltol * fabs(CSoln[kindex]);
753  }
754  }
755  }
756  /*
757  * Now do the residual Weights. Since we have the Jacobian, we
758  * will use it to generate a number based on the what a significant
759  * change in a solution variable does to each residual.
760  * This is a row sum scale operation.
761  */
762  for (k = 0; k < m_neq; k++) {
763  wtResid[k] = 0.0;
764  for (jcol = 0; jcol < m_neq; jcol++) {
765  wtResid[k] += fabs(Jac(k,jcol) * wtSpecies[jcol]);
766  }
767  }
768 }
769
770 doublereal solveSP::
771 calc_t(doublereal netProdRateSolnSP[], doublereal XMolSolnSP[],
772  int* label, int* label_old, doublereal* label_factor, int ioflag)
773 {
774  size_t k, isp, nsp, kstart;
775  doublereal inv_timeScale = 1.0E-10;
776  doublereal sden, tmp;
777  size_t kindexSP = 0;
778  *label = 0;
779  updateMFSolnSP(XMolSolnSP);
780  for (isp = 0; isp < m_numSurfPhases; isp++) {
781  nsp = m_nSpeciesSurfPhase[isp];
782
783  // Get the interface kinetics associated with this surface
784  InterfaceKinetics* m_kin = m_objects[isp];
785
786  // Calculate the start of the species index for surfaces within
787  // the InterfaceKinetics object
788  size_t surfIndex = m_kin->surfacePhaseIndex();
789  kstart = m_kin->kineticsSpeciesIndex(0, surfIndex);
790  ThermoPhase& THref = m_kin->thermo(surfIndex);
791
793
794  sden = THref.molarDensity();
795  for (k = 0; k < nsp; k++, kindexSP++) {
796  size_t kspindex = kstart + k;
797  netProdRateSolnSP[kindexSP] = m_numEqn1[kspindex];
798  if (XMolSolnSP[kindexSP] <= 1.0E-10) {
799  tmp = 1.0E-10;
800  } else {
801  tmp = XMolSolnSP[kindexSP];
802  }
803  tmp *= sden;
804  tmp = fabs(netProdRateSolnSP[kindexSP]/ tmp);
805  if (netProdRateSolnSP[kindexSP]> 0.0) {
806  tmp /= 100.;
807  }
808  if (tmp > inv_timeScale) {
809  inv_timeScale = tmp;
810  *label = int(kindexSP);
811  }
812  }
813  }
814
815  /*
816  * Increase time step exponentially as same species repeatedly
817  * controls time step
818  */
819  if (*label == *label_old) {
820  *label_factor *= 1.5;
821  } else {
822  *label_old = *label;
823  *label_factor = 1.0;
824  }
825  return inv_timeScale / *label_factor;
826 } /* calc_t */
827
828 void solveSP::print_header(int ioflag, int ifunc, doublereal time_scale,
829  int damping, doublereal reltol, doublereal abstol,
830  doublereal TKelvin,
831  doublereal PGas, doublereal netProdRate[],
832  doublereal XMolKinSpecies[])
833 {
834  if (ioflag) {
835  printf("\n================================ SOLVESP CALL SETUP "
836  "========================================\n");
837  if (ifunc == SFLUX_INITIALIZE) {
838  printf("\n SOLVESP Called with Initialization turned on\n");
839  printf(" Time scale input = %9.3e\n", time_scale);
840  } else if (ifunc == SFLUX_RESIDUAL) {
841  printf("\n SOLVESP Called to calculate steady state residual\n");
842  printf(" from a good initial guess\n");
843  } else if (ifunc == SFLUX_JACOBIAN) {
844  printf("\n SOLVESP Called to calculate steady state jacobian\n");
845  printf(" from a good initial guess\n");
846  } else if (ifunc == SFLUX_TRANSIENT) {
847  printf("\n SOLVESP Called to integrate surface in time\n");
848  printf(" for a total of %9.3e sec\n", time_scale);
849  } else {
850  fprintf(stderr,"Unknown ifunc flag = %d\n", ifunc);
851  exit(EXIT_FAILURE);
852  }
853
854  if (m_bulkFunc == BULK_DEPOSITION) {
855  printf(" The composition of the Bulk Phases will be calculated\n");
856  } else if (m_bulkFunc == BULK_ETCH) {
857  printf(" Bulk Phases have fixed compositions\n");
858  } else {
859  fprintf(stderr,"Unknown bulkFunc flag = %d\n", m_bulkFunc);
860  exit(EXIT_FAILURE);
861  }
862
863  if (damping) {
864  printf(" Damping is ON \n");
865  } else {
866  printf(" Damping is OFF \n");
867  }
868
869  printf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, abstol);
870  }
871
872  if (ioflag == 1) {
873  printf("\n\n\t Iter Time Del_t Damp DelX "
874  " Resid Name-Time Name-Damp\n");
875  printf("\t -----------------------------------------------"
876  "------------------------------------\n");
877  }
878 }
879
880 void solveSP::printIteration(int ioflag, doublereal damp, int label_d,
881  int label_t,
882  doublereal inv_t, doublereal t_real, size_t iter,
883  doublereal update_norm, doublereal resid_norm,
884  doublereal netProdRate[], doublereal CSolnSP[],
885  doublereal resid[], doublereal XMolSolnSP[],
886  doublereal wtSpecies[], size_t dim, bool do_time)
887 {
888  size_t i, k;
889  string nm;
890  if (ioflag == 1) {
891
892  printf("\t%6s ", int2str(iter).c_str());
893  if (do_time) {
894  printf("%9.4e %9.4e ", t_real, 1.0/inv_t);
895  } else
896  for (i = 0; i < 22; i++) {
897  printf(" ");
898  }
899  if (damp < 1.0) {
900  printf("%9.4e ", damp);
901  } else
902  for (i = 0; i < 11; i++) {
903  printf(" ");
904  }
905  printf("%9.4e %9.4e", update_norm, resid_norm);
906  if (do_time) {
907  k = m_kinSpecIndex[label_t];
908  size_t isp = m_kinObjIndex[label_t];
909  InterfaceKinetics* m_kin = m_objects[isp];
910  nm = m_kin->kineticsSpeciesName(k);
911  printf(" %-16s", nm.c_str());
912  } else {
913  for (i = 0; i < 16; i++) {
914  printf(" ");
915  }
916  }
917  if (label_d >= 0) {
918  k = m_kinSpecIndex[label_d];
919  size_t isp = m_kinObjIndex[label_d];
920  InterfaceKinetics* m_kin = m_objects[isp];
921  nm = m_kin->kineticsSpeciesName(k);
922  printf(" %-16s", nm.c_str());
923  }
924  printf("\n");
925  }
926 } /* printIteration */
927
928
929 void solveSP::printFinal(int ioflag, doublereal damp, int label_d, int label_t,
930  doublereal inv_t, doublereal t_real, size_t iter,
931  doublereal update_norm, doublereal resid_norm,
932  doublereal netProdRateKinSpecies[], const doublereal CSolnSP[],
933  const doublereal resid[], doublereal XMolSolnSP[],
934  const doublereal wtSpecies[], const doublereal wtRes[],
935  size_t dim, bool do_time,
936  doublereal TKelvin, doublereal PGas)
937 {
938  size_t i, k;
939  string nm;
940  if (ioflag == 1) {
941
942  printf("\tFIN%3s ", int2str(iter).c_str());
943  if (do_time) {
944  printf("%9.4e %9.4e ", t_real, 1.0/inv_t);
945  } else
946  for (i = 0; i < 22; i++) {
947  printf(" ");
948  }
949  if (damp < 1.0) {
950  printf("%9.4e ", damp);
951  } else
952  for (i = 0; i < 11; i++) {
953  printf(" ");
954  }
955  printf("%9.4e %9.4e", update_norm, resid_norm);
956  if (do_time) {
957  k = m_kinSpecIndex[label_t];
958  size_t isp = m_kinObjIndex[label_t];
959  InterfaceKinetics* m_kin = m_objects[isp];
960  nm = m_kin->kineticsSpeciesName(k);
961  printf(" %-16s", nm.c_str());
962  } else {
963  for (i = 0; i < 16; i++) {
964  printf(" ");
965  }
966  }
967  if (label_d >= 0) {
968  k = m_kinSpecIndex[label_d];
969  size_t isp = m_kinObjIndex[label_d];
970  InterfaceKinetics* m_kin = m_objects[isp];
971  nm = m_kin->kineticsSpeciesName(k);
972  printf(" %-16s", nm.c_str());
973  }
974  printf(" -- success\n");
975  }
976 }
977
978 }
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
int m_bulkFunc
This variable determines how the bulk phases are to be handled.
Definition: solveSP.h:358
virtual void getNetProductionRates(doublereal *net)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
std::vector< InterfaceKinetics * > & m_objects
Vector of interface kinetics objects.
Definition: solveSP.h:346
void resize(size_t n, size_t m, doublereal v=0.0)
Resize the array, and fill the new entries with 'v'.
Definition: Array.h:128
size_t m_neq
Total number of equations to solve in the implicit problem.
Definition: solveSP.h:352
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
std::vector< doublereal * > m_JacCol
Vector of pointers to the top of the columns of the Jacobian.
Definition: solveSP.h:544
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:288
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
std::string kineticsSpeciesName(size_t k) const
Return the name of the kth species in the kinetics manager.
Definition: Kinetics.cpp:186
vector_fp m_CSolnSPInit
Saved initial solution vector. length MAX(1, m_neq)
Definition: solveSP.h:506
size_t m_numBulkPhasesSS
Total number of volumetric condensed phases included in the steady state problem handled by this rout...
Definition: solveSP.h:429
void updateMFSolnSP(doublereal *XMolSolnSP)
Update mole fraction vector consisting of unknowns in surface problem.
Definition: solveSP.cpp:439
vector_fp m_CSolnSP
Solution vector. length MAX(1, m_neq)
Definition: solveSP.h:503
doublereal molarDensity() const
Molar density (kmol/m^3).
Definition: Phase.cpp:597
doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are contiguous in memory.
Definition: Array.h:363
const int BULK_ETCH
Etching of a bulk phase is to be expected.
Definition: solveSP.h:60
void getMoleFractions(doublereal *const x) const
Get the species mole fraction vector.
Definition: Phase.cpp:519
void getConcentrations(doublereal *const c) const
Get the species concentrations (kmol/m^3).
Definition: Phase.cpp:572
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:33
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:101
void evalSurfLarge(const doublereal *CSolnSP)
Update the vector that keeps track of the largest species in each surface phase.
Definition: solveSP.cpp:458
A simple thermodynamic model for a surface phase, assuming an ideal solution model.
Definition: SurfPhase.h:143
size_t surfacePhaseIndex()
This returns the integer index of the phase which has ThermoPhase type cSurf.
Definition: Kinetics.h:263
size_t m_numTotSurfSpecies
Total number of surface species in all surface phases.
Definition: solveSP.h:373
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:141
const int SFLUX_INITIALIZE
This assumes that the initial guess supplied to the routine is far from the correct one...
Definition: solveSP.h:25
void print_header(int ioflag, int ifunc, doublereal time_scale, int damping, doublereal reltol, doublereal abstol, doublereal TKelvin, doublereal PGas, doublereal netProdRate[], doublereal XMolKinSpecies[])
Printing routine that optionally gets called at the start of every invocation.
Definition: solveSP.cpp:828
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:228
std::vector< size_t > m_eqnIndexStartSolnPhase
Index of the start of the unknowns for each solution phase.
Definition: solveSP.h:407
Array2D m_Jac
Jacobian.
Definition: solveSP.h:548
A kinetics manager for heterogeneous reaction mechanisms.
vector_fp m_wtResid
Weights for the residual norm calculation. length MAX(1, m_neq)
Definition: solveSP.h:512
const int SFLUX_RESIDUAL
Need to solve the surface problem in order to calculate the surface fluxes of gas-phase species...
Definition: solveSP.h:32
vector_fp m_XMolKinSpecies
Vector of mole fractions. length m_maxTotSpecies.
Definition: solveSP.h:534
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:771
vector_fp m_numEqn1
Temporary vector with length equal to max m_maxTotSpecies.
Definition: solveSP.h:494
const int BULK_DEPOSITION
Deposition of a bulk phase is to be expected.
Definition: solveSP.h:55
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:38
std::vector< size_t > m_spSurfLarge
Vector containing the indices of the largest species in each surface phase.
Definition: solveSP.h:475
vector_fp m_CSolnSave
Temporary vector with length equal to max m_maxTotSpecies.
Definition: solveSP.h:500
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:388
vector_int m_ipiv
pivots. length MAX(1, m_neq)
Definition: solveSP.h:537
void updateMFKinSpecies(doublereal *XMolKinSp, int isp)
Update the mole fraction vector for a specific kinetic species vector corresponding to one InterfaceK...
Definition: solveSP.cpp:447
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
vector_fp m_numEqn2
Temporary vector with length equal to max m_maxTotSpecies.
Definition: solveSP.h:497
const int SFLUX_TRANSIENT
The transient calculation is performed here for an amount of time specified by "time_scale".
Definition: solveSP.h:45
size_t kineticsSpeciesIndex(size_t k, size_t n) const
The location of species k of phase n in species arrays.
Definition: Kinetics.h:331
std::vector< SurfPhase * > m_ptrsSurfPhase
Vector of surface phase pointers.
Definition: solveSP.h:395
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:252
std::vector< size_t > m_kinObjPhaseIDSurfPhase
Phase ID in the InterfaceKinetics object of the surface phase.
Definition: solveSP.h:416
vector_fp m_netProductionRatesSave
Temporary vector with length equal to max m_maxTotSpecies.
Definition: solveSP.h:491
std::vector< size_t > m_kinObjIndex
Index between the equation index and the index of the InterfaceKinetics object.
Definition: solveSP.h:466
Advances the surface coverages of the associated set of SurfacePhase objects in time.
size_t m_numSurfPhases
Number of surface phases in the surface problem.
Definition: solveSP.h:365
std::vector< ThermoPhase * > m_bulkPhasePtrs
Vector of bulk phase pointers, length is equal to m_numBulkPhases.
Definition: solveSP.h:449
~solveSP()
Destructor. Deletes the integrator.
Definition: solveSP.cpp:137
std::vector< size_t > m_numBulkSpecies
Vector of number of species in the m_numBulkPhases phases.
Definition: solveSP.h:435
vector_fp m_CSolnSPOld
Saved solution vector at the old time step. length MAX(1, m_neq)
Definition: solveSP.h:509
vector_fp m_resid
Residual for the surface problem.
Definition: solveSP.h:531
Header file for implicit surface problem solver (see Chemical Kinetics and class solveSP).
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
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:459
size_t m_maxTotSpecies
Maximum number of species in any single kinetics operator -> also maxed wrt the total # of solution s...
Definition: solveSP.h:488
std::vector< size_t > m_indexKinObjSurfPhase
Mapping between the surface phases and the InterfaceKinetics objects.
Definition: solveSP.h:381
vector_fp m_wtSpecies
Weights for the species concentrations norm calculation.
Definition: solveSP.h:518
void printFinal(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, doublereal netProdRateKinSpecies[], const doublereal CSolnSP[], const doublereal resid[], doublereal XMolSolnSP[], const doublereal wtSpecies[], const doublereal wtRes[], size_t dim, bool do_time, doublereal TKelvin, doublereal PGas)
Print a summary of the solution.
Definition: solveSP.cpp:929
void resjac_eval(std::vector< doublereal * > &JacCol, 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:589
void fun_eval(doublereal *resid, const doublereal *CSolnSP, const doublereal *CSolnOldSP, const bool do_time, const doublereal deltaT)
Main Function evaluation.
Definition: solveSP.cpp:475
size_t m_numTotBulkSpeciesSS
Total number of species in all bulk phases.
Definition: solveSP.h:446
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, doublereal netProdRate[], doublereal CSolnSP[], doublereal resid[], doublereal XMolSolnSP[], doublereal wtSpecies[], size_t dim, bool do_time)
Printing routine that gets called after every iteration.
Definition: solveSP.cpp:880
void updateState(const doublereal *cSurfSpec)
Update the surface states of the surface phases.
Definition: solveSP.cpp:430
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:729