Cantera  2.2.1
solveProb.cpp
Go to the documentation of this file.
1 /**
2  * @file: solveProb.cpp Implicit solver for nonlinear problems
3  */
4
5 /*
6  * Copyright 2004 Sandia Corporation. Under the terms of Contract
7  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
8  * retains certain rights in this software.
9  * See file License.txt for licensing information.
10  */
11
13 #include "cantera/base/clockWC.h"
15 #include "cantera/base/global.h"
16
17 using namespace std;
18 namespace Cantera
19 {
20
21 /***************************************************************************
22  * STATIC ROUTINES DEFINED IN THIS FILE
23  ***************************************************************************/
24
25 static doublereal calcWeightedNorm(const doublereal [], const doublereal dx[], size_t);
26
27 solveProb::solveProb(ResidEval* resid) :
28  m_residFunc(resid),
29  m_neq(0),
30  m_atol(0),
31  m_rtol(1.0E-4),
32  m_maxstep(1000),
33  m_ioflag(0)
34 {
35  warn_deprecated("class solveProb", "To be removed after Cantera 2.2.");
37
38  // Dimension solution vector
39  size_t dim1 = std::max<size_t>(1, m_neq);
40
41  m_atol.resize(dim1, 1.0E-9);
42  m_netProductionRatesSave.resize(dim1, 0.0);
43  m_numEqn1.resize(dim1, 0.0);
44  m_numEqn2.resize(dim1, 0.0);
45  m_CSolnSave.resize(dim1, 0.0);
46  m_CSolnSP.resize(dim1, 0.0);
47  m_CSolnSPInit.resize(dim1, 0.0);
48  m_CSolnSPOld.resize(dim1, 0.0);
49  m_wtResid.resize(dim1, 0.0);
50  m_wtSpecies.resize(dim1, 0.0);
51  m_resid.resize(dim1, 0.0);
52  m_topBounds.resize(dim1, 1.0);
53  m_botBounds.resize(dim1, 0.0);
54
55  m_Jac.resize(dim1, dim1, 0.0);
56  m_JacCol.resize(dim1, 0);
57  for (size_t k = 0; k < dim1; k++) {
58  m_JacCol[k] = m_Jac.ptrColumn(k);
59  }
60
61 }
62
63 int solveProb::solve(int ifunc, doublereal time_scale,
64  doublereal reltol)
65 {
66  /*
67  * The following calculation is a Newton's method to get the surface fractions
68  * of the surface and bulk species by requiring that the surface species
69  * production rate = 0 and that the bulk fractions are proportional to their
70  * production rates.
71  */
72  doublereal EXTRA_ACCURACY = 0.001;
73  if (ifunc == SOLVEPROB_JACOBIAN) {
74  EXTRA_ACCURACY *= 0.001;
75  }
76  int info = 0;
77  size_t label_t = npos; /* Species IDs for time control */
78  size_t label_d; /* Species IDs for damping control */
79  size_t label_t_old = npos;
80  doublereal label_factor = 1.0;
81  int iter=0; // iteration number on numlinear solver
82  int iter_max=1000; // maximum number of nonlinear iterations
83  doublereal deltaT = 1.0E-10; // Delta time step
84  doublereal damp=1.0, tmp;
85  // Weighted L2 norm of the residual. Currently, this is only
86  // used for IO purposes. It doesn't control convergence.
87  // Therefore, it is turned off when DEBUG_SOLVEPROB isn't defined.
88  doublereal resid_norm;
89  doublereal inv_t = 0.0;
90  doublereal t_real = 0.0, update_norm = 1.0E6;
91
92  bool do_time = false, not_converged = true;
93
94 #ifdef DEBUG_SOLVEPROB
95 #ifdef DEBUG_SOLVEPROB_TIME
96  doublereal t1;
97 #endif
98 #else
99  if (m_ioflag > 1) {
100  m_ioflag = 1;
101  }
102 #endif
103
104 #ifdef DEBUG_SOLVEPROB
105 #ifdef DEBUG_SOLVEPROB_TIME
106  Cantera::clockWC wc;
107  if (m_ioflag) {
108  t1 = wc.secondsWC();
109  }
110 #endif
111 #endif
112
113  /*
114  * Set the initial value of the do_time parameter
115  */
116  if (ifunc == SOLVEPROB_INITIALIZE || ifunc == SOLVEPROB_TRANSIENT) {
117  do_time = true;
118  }
119
120  /*
121  * upload the initial conditions
122  */
124
125  /*
126  * Store the initial guess in the soln vector,
127  * CSolnSP, and in an separate vector CSolnSPInit.
128  */
129  std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPInit.begin());
130
131
132
133  if (m_ioflag) {
136  }
137
138  /*
139  * Quick return when there isn't a surface problem to solve
140  */
141  if (m_neq == 0) {
142  not_converged = false;
143  update_norm = 0.0;
144  }
145
146  /* ------------------------------------------------------------------
147  * Start of Newton's method
148  * ------------------------------------------------------------------
149  */
150  while (not_converged && iter < iter_max) {
151  iter++;
152  /*
153  * Store previous iteration's solution in the old solution vector
154  */
155  std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPOld.begin());
156
157  /*
158  * Calculate the value of the time step
159  * - heuristics to stop large oscillations in deltaT
160  */
161  if (do_time) {
162  /* don't hurry increase in time step at the same time as damping */
163  if (damp < 1.0) {
164  label_factor = 1.0;
165  }
167  &label_t, &label_t_old, &label_factor, m_ioflag);
168  if (iter < 10) {
169  inv_t = tmp;
170  } else if (tmp > 2.0*inv_t) {
171  inv_t = 2.0*inv_t;
172  } else {
173  inv_t = tmp;
174  }
175
176  /*
177  * Check end condition
178  */
179
180  if (ifunc == SOLVEPROB_TRANSIENT) {
181  tmp = t_real + 1.0/inv_t;
182  if (tmp > time_scale) {
183  inv_t = 1.0/(time_scale - t_real);
184  }
185  }
186  } else {
187  /* make steady state calc a step of 1 million seconds to
188  prevent singular Jacobians for some pathological cases */
189  inv_t = 1.0e-6;
190  }
191  deltaT = 1.0/inv_t;
192
193  /*
194  * Call the routine to numerically evaluation the Jacobian
195  * and residual for the current iteration.
196  */
198  DATA_PTR(m_CSolnSPOld), do_time, deltaT);
199
200  /*
201  * Calculate the weights. Make sure the calculation is carried
202  * out on the first iteration.
203  */
204  if (iter%4 == 1) {
207  }
208
209  /*
210  * Find the weighted norm of the residual
211  */
212  resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid), DATA_PTR(m_resid), m_neq);
213
214 #ifdef DEBUG_SOLVEPROB
215  if (m_ioflag > 1) {
216  printIterationHeader(m_ioflag, damp, inv_t, t_real, iter, do_time);
217  /*
218  * Print out the residual and Jacobian
219  */
220  printResJac(m_ioflag, m_neq, m_Jac, DATA_PTR(m_resid),
221  DATA_PTR(m_wtResid), resid_norm);
222  }
223 #endif
224
225  /*
226  * Solve Linear system (with LAPACK). The solution is in resid[]
227  */
228  info = m_Jac.factor();
229  if (info==0) {
230  m_Jac.solve(&m_resid[0]);
231  }
232  /*
233  * Force convergence if residual is small to avoid
234  * "nan" results from the linear solve.
235  */
236  else {
237  if (m_ioflag) {
238  printf("solveSurfSS: Zero pivot, assuming converged: %g (%d)\n",
239  resid_norm, info);
240  }
241  for (size_t jcol = 0; jcol < m_neq; jcol++) {
242  m_resid[jcol] = 0.0;
243  }
244
245  /* print out some helpful info */
246  if (m_ioflag > 1) {
247  printf("-----\n");
248  printf("solveSurfProb: iter %d t_real %g delta_t %g\n\n",
249  iter,t_real, 1.0/inv_t);
250  printf("solveSurfProb: init guess, current concentration,"
251  "and prod rate:\n");
252
253  printf("-----\n");
254  }
255  if (do_time) {
256  t_real += time_scale;
257  }
258 #ifdef DEBUG_SOLVEPROB
259  if (m_ioflag) {
260  printf("\nResidual is small, forcing convergence!\n");
261  }
262 #endif
263  }
264
265  /*
266  * Calculate the Damping factor needed to keep all unknowns
267  * between 0 and 1, and not allow too large a change (factor of 2)
268  * in any unknown.
269  */
270
271
272  damp = calc_damping(DATA_PTR(m_CSolnSP), DATA_PTR(m_resid), m_neq, &label_d);
273
274
275  /*
276  * Calculate the weighted norm of the update vector
277  * Here, resid is the delta of the solution, in concentration
278  * units.
279  */
280  update_norm = calcWeightedNorm(DATA_PTR(m_wtSpecies),
282  /*
283  * Update the solution vector and real time
284  * Crop the concentrations to zero.
285  */
286  for (size_t irow = 0; irow < m_neq; irow++) {
287  m_CSolnSP[irow] -= damp * m_resid[irow];
288  }
289
290
291  if (do_time) {
292  t_real += damp/inv_t;
293  }
294
295  if (m_ioflag) {
296  printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
297  update_norm, resid_norm,
300  DATA_PTR(m_wtSpecies), m_neq, do_time);
301  }
302
303  if (ifunc == SOLVEPROB_TRANSIENT) {
304  not_converged = (t_real < time_scale);
305  } else {
306  if (do_time) {
307  if (t_real > time_scale ||
308  (resid_norm < 1.0e-7 &&
309  update_norm*time_scale/t_real < EXTRA_ACCURACY)) {
310  do_time = false;
311 #ifdef DEBUG_SOLVEPROB
312  if (m_ioflag > 1) {
314  }
315 #endif
316  }
317  } else {
318  not_converged = ((update_norm > EXTRA_ACCURACY) ||
319  (resid_norm > EXTRA_ACCURACY));
320  }
321  }
322  } /* End of Newton's Method while statement */
323
324  /*
325  * End Newton's method. If not converged, print error message and
326  * recalculate sdot's at equal site fractions.
327  */
328  if (not_converged) {
329  if (m_ioflag) {
330  printf("#$#$#$# Error in solveProb$#$#$#\$ \n");
331  printf("Newton iter on surface species did not converge, "
332  "update_norm = %e \n", update_norm);
333  printf("Continuing anyway\n");
334  }
335  }
336 #ifdef DEBUG_SOLVEPROB
337 #ifdef DEBUG_SOLVEPROB_TIME
338  if (m_ioflag) {
339  printf("\nEnd of solve, time used: %e\n", wc.secondsWC()-t1);
340  }
341 #endif
342 #endif
343
344  /*
345  * Decide on what to return in the solution vector
346  * - right now, will always return the last solution
347  * no matter how bad
348  */
349  if (m_ioflag) {
351  false, deltaT);
352  resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid),
354  printFinal(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
355  update_norm, resid_norm, DATA_PTR(m_netProductionRatesSave),
358  DATA_PTR(m_wtResid), m_neq, do_time);
359  }
360
361  /*
362  * Return with the appropriate flag
363  */
364  if (update_norm > 1.0) {
365  return -1;
366  }
367  return 0;
368 }
369
370 void solveProb::reportState(doublereal* const CSolnSP) const
371 {
372  std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), CSolnSP);
373 }
374
375 void solveProb::fun_eval(doublereal* const resid, const doublereal* const CSoln,
376  const doublereal* const CSolnOld, const bool do_time,
377  const doublereal deltaT)
378 {
379  /*
380  * This routine uses the m_numEqn1 and m_netProductionRatesSave vectors
381  * as temporary internal storage.
382  */
383  if (do_time) {
384  m_residFunc->evalSimpleTD(0.0, CSoln, CSolnOld, deltaT, resid);
385  } else {
386  m_residFunc->evalSS(0.0, CSoln, resid);
387  }
388 }
389
390 void solveProb::resjac_eval(std::vector<doublereal*> &JacCol,
391  doublereal resid[], doublereal CSoln[],
392  const doublereal CSolnOld[], const bool do_time,
393  const doublereal deltaT)
394 {
395  doublereal dc, cSave, sd;
396  doublereal* col_j;
397  /*
398  * Calculate the residual
399  */
400  fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
401  /*
402  * Now we will look over the columns perturbing each unknown.
403  */
404
405  for (size_t kCol = 0; kCol < m_neq; kCol++) {
406  cSave = CSoln[kCol];
407  sd = fabs(cSave) + fabs(CSoln[kCol]) + m_atol[kCol] * 1.0E6;
408  if (sd < 1.0E-200) {
409  sd = 1.0E-4;
410  }
411  dc = std::max(1.0E-11 * sd, fabs(cSave) * 1.0E-6);
412  CSoln[kCol] += dc;
413  // Use the m_numEqn2 vector as temporary internal storage.
414  fun_eval(DATA_PTR(m_numEqn2), CSoln, CSolnOld, do_time, deltaT);
415  col_j = JacCol[kCol];
416  for (size_t i = 0; i < m_neq; i++) {
417  col_j[i] = (m_numEqn2[i] - resid[i])/dc;
418  }
419  CSoln[kCol] = cSave;
420  }
421
422 }
423
424 doublereal solveProb::calc_damping(doublereal x[], doublereal dxneg[], size_t dim, size_t* label)
425 {
426  const doublereal APPROACH = 0.50;
427  doublereal damp = 1.0, xnew, xtop, xbot;
428  static doublereal damp_old = 1.0;
429  *label = npos;
430
431  for (size_t i = 0; i < dim; i++) {
432  doublereal topBounds = m_topBounds[i];
433  doublereal botBounds = m_botBounds[i];
434  /*
435  * Calculate the new suggested new value of x[i]
436  */
437  double delta_x = - dxneg[i];
438  xnew = x[i] - damp * dxneg[i];
439
440  /*
441  * Calculate the allowed maximum and minimum values of x[i]
442  * - Only going to allow x[i] to converge to the top and bottom bounds by a
443  * single order of magnitude at one time
444  */
445  bool canCrossOrigin = false;
446  if (topBounds > 0.0 && botBounds < 0.0) {
447  canCrossOrigin = true;
448  }
449
450  xtop = topBounds - 0.1 * fabs(topBounds - x[i]);
451
452  xbot = botBounds + 0.1 * fabs(x[i] - botBounds);
453
454  if (xnew > xtop) {
455  damp = - APPROACH * (xtop - x[i]) / dxneg[i];
456  *label = i;
457  } else if (xnew < xbot) {
458  damp = APPROACH * (x[i] - xbot) / dxneg[i];
459  *label = i;
460  }
461  double denom = fabs(x[i]) + 1.0E5 * m_atol[i];
462  if ((fabs(delta_x) / denom) > 0.3) {
463  double newdamp = 0.3 * denom / fabs(delta_x);
464  if (canCrossOrigin) {
465  if (xnew * x[i] < 0.0) {
466  if (fabs(x[i]) < 1.0E8 * m_atol[i]) {
467  newdamp = 2.0 * fabs(x[i]) / fabs(delta_x);
468  }
469  }
470  }
471  damp = std::min(damp, newdamp);
472  }
473
474  }
475
476  /*
477  * Only allow the damping parameter to increase by a factor of three each
478  * iteration. Heuristic to avoid oscillations in the value of damp
479  */
480  if (damp > damp_old*3) {
481  damp = damp_old*3;
482  *label = npos;
483  }
484
485  /*
486  * Save old value of the damping parameter for use
487  * in subsequent calls.
488  */
489  damp_old = damp;
490  return damp;
491
492 }
493
494 /*
495  * This function calculates the norm of an update, dx[],
496  * based on the weighted values of x.
497  */
498 static doublereal calcWeightedNorm(const doublereal wtX[], const doublereal dx[], size_t dim)
499 {
500  doublereal norm = 0.0;
501  doublereal tmp;
502  if (dim == 0) {
503  return 0.0;
504  }
505  for (size_t i = 0; i < dim; i++) {
506  tmp = dx[i] / wtX[i];
507  norm += tmp * tmp;
508  }
509  return sqrt(norm/dim);
510 }
511
512 void solveProb::calcWeights(doublereal wtSpecies[], doublereal wtResid[],
513  const doublereal CSoln[])
514 {
515  /*
516  * First calculate the weighting factor
517  */
518
519  for (size_t k = 0; k < m_neq; k++) {
520  wtSpecies[k] = m_atol[k] + m_rtol * fabs(CSoln[k]);
521  }
522  /*
523  * Now do the residual Weights. Since we have the Jacobian, we
524  * will use it to generate a number based on the what a significant
525  * change in a solution variable does to each residual.
526  * This is a row sum scale operation.
527  */
528  for (size_t k = 0; k < m_neq; k++) {
529  wtResid[k] = 0.0;
530  for (size_t jcol = 0; jcol < m_neq; jcol++) {
531  wtResid[k] += fabs(m_Jac(k,jcol) * wtSpecies[jcol]);
532  }
533  }
534 }
535
536 doublereal solveProb::calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[],
537  size_t* label, size_t* label_old,
538  doublereal* label_factor, int ioflag)
539 {
540  doublereal tmp, inv_timeScale=0.0;
541  for (size_t k = 0; k < m_neq; k++) {
542  if (Csoln[k] <= 1.0E-10) {
543  tmp = 1.0E-10;
544  } else {
545  tmp = Csoln[k];
546  }
547  tmp = fabs(netProdRateSolnSP[k]/ tmp);
548
549
550  if (netProdRateSolnSP[k]> 0.0) {
551  tmp /= 100.;
552  }
553  if (tmp > inv_timeScale) {
554  inv_timeScale = tmp;
555  *label = k;
556  }
557  }
558
559  /*
560  * Increase time step exponentially as same species repeatedly
561  * controls time step
562  */
563  if (*label == *label_old) {
564  *label_factor *= 1.5;
565  } else {
566  *label_old = *label;
567  *label_factor = 1.0;
568  }
569  inv_timeScale = inv_timeScale / *label_factor;
570 #ifdef DEBUG_SOLVEPROB
571  if (ioflag > 1) {
572  if (*label_factor > 1.0) {
573  printf("Delta_t increase due to repeated controlling species = %e\n",
574  *label_factor);
575  }
576  int kkin = m_kinSpecIndex[*label];
577
578  string sn = " "
579  printf("calc_t: spec=%d(%s) sf=%e pr=%e dt=%e\n",
580  *label, sn.c_str(), XMolSolnSP[*label],
581  netProdRateSolnSP[*label], 1.0/inv_timeScale);
582  }
583 #endif
584
585  return inv_timeScale;
586
587 }
588
589 void solveProb::setBounds(const doublereal botBounds[], const doublereal topBounds[])
590 {
591  for (size_t k = 0; k < m_neq; k++) {
592  m_botBounds[k] = botBounds[k];
593  m_topBounds[k] = topBounds[k];
594  }
595 }
596
597 #ifdef DEBUG_SOLVEPROB
598 void solveProb::printResJac(int ioflag, int neq, const Array2D& Jac,
599  doublereal resid[], doublereal wtRes[],
600  doublereal norm)
601 {
602
603 }
604 #endif
605
606 void solveProb::print_header(int ioflag, int ifunc, doublereal time_scale,
607  doublereal reltol,
608  doublereal netProdRate[])
609 {
610  int damping = 1;
611  if (ioflag) {
612  printf("\n================================ SOLVEPROB CALL SETUP "
613  "========================================\n");
614  if (ifunc == SOLVEPROB_INITIALIZE) {
615  printf("\n SOLVEPROB Called with Initialization turned on\n");
616  printf(" Time scale input = %9.3e\n", time_scale);
617  } else if (ifunc == SOLVEPROB_RESIDUAL) {
618  printf("\n SOLVEPROB Called to calculate steady state residual\n");
619  printf(" from a good initial guess\n");
620  } else if (ifunc == SOLVEPROB_JACOBIAN) {
621  printf("\n SOLVEPROB Called to calculate steady state Jacobian\n");
622  printf(" from a good initial guess\n");
623  } else if (ifunc == SOLVEPROB_TRANSIENT) {
624  printf("\n SOLVEPROB Called to integrate surface in time\n");
625  printf(" for a total of %9.3e sec\n", time_scale);
626  } else {
628  "Unknown ifunc flag = " + int2str(ifunc));
629  }
630
631
632
633  if (damping) {
634  printf(" Damping is ON \n");
635  } else {
636  printf(" Damping is OFF \n");
637  }
638
639  printf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, m_atol[0]);
640  }
641
642  /*
643  * Print out the initial guess
644  */
645 #ifdef DEBUG_SOLVEPROB
646  if (ioflag > 1) {
647  printf("\n================================ INITIAL GUESS "
648  "========================================\n");
649  int kindexSP = 0;
650  for (int isp = 0; isp < m_numSurfPhases; isp++) {
651  InterfaceKinetics* m_kin = m_objects[isp];
652  int surfIndex = m_kin->surfacePhaseIndex();
653  int nPhases = m_kin->nPhases();
654  m_kin->getNetProductionRates(netProdRate);
655  updateMFKinSpecies(XMolKinSpecies, isp);
656
657  printf("\n IntefaceKinetics Object # %d\n\n", isp);
658
659  printf("\t Number of Phases = %d\n", nPhases);
660  printf("\t Phase:SpecName Prod_Rate MoleFraction kindexSP\n");
661  printf("\t -------------------------------------------------------"
662  "----------\n");
663
664  int kspindex = 0;
665  bool inSurfacePhase = false;
666  for (int ip = 0; ip < nPhases; ip++) {
667  if (ip == surfIndex) {
668  inSurfacePhase = true;
669  } else {
670  inSurfacePhase = false;
671  }
672  ThermoPhase& THref = m_kin->thermo(ip);
673  int nsp = THref.nSpecies();
674  string pname = THref.id();
675  for (int k = 0; k < nsp; k++) {
676  string sname = THref.speciesName(k);
677  string cname = pname + ":" + sname;
678  if (inSurfacePhase) {
679  printf("\t %-24s %10.3e %10.3e %d\n", cname.c_str(),
680  netProdRate[kspindex], XMolKinSpecies[kspindex],
681  kindexSP);
682  kindexSP++;
683  } else {
684  printf("\t %-24s %10.3e %10.3e\n", cname.c_str(),
685  netProdRate[kspindex], XMolKinSpecies[kspindex]);
686  }
687  kspindex++;
688  }
689  }
690  printf("=========================================================="
691  "=================================\n");
692  }
693  }
694 #endif
695  if (ioflag == 1) {
696  printf("\n\n\t Iter Time Del_t Damp DelX "
697  " Resid Name-Time Name-Damp\n");
698  printf("\t -----------------------------------------------"
699  "------------------------------------\n");
700  }
701 }
702
703 void solveProb::printIteration(int ioflag, doublereal damp, size_t label_d,
704  size_t label_t,
705  doublereal inv_t, doublereal t_real, int iter,
706  doublereal update_norm, doublereal resid_norm,
707  doublereal netProdRate[], doublereal CSolnSP[],
708  doublereal resid[],
709  doublereal wtSpecies[], size_t dim, bool do_time)
710 {
711  size_t i, k;
712  string nm;
713  if (ioflag == 1) {
714
715  printf("\t%6d ", iter);
716  if (do_time) {
717  printf("%9.4e %9.4e ", t_real, 1.0/inv_t);
718  } else
719  for (i = 0; i < 22; i++) {
720  printf(" ");
721  }
722  if (damp < 1.0) {
723  printf("%9.4e ", damp);
724  } else
725  for (i = 0; i < 11; i++) {
726  printf(" ");
727  }
728  printf("%9.4e %9.4e", update_norm, resid_norm);
729  if (do_time) {
730  k = label_t;
731  printf(" %s", int2str(k).c_str());
732  } else {
733  for (i = 0; i < 16; i++) {
734  printf(" ");
735  }
736  }
737  if (label_d != npos) {
738  k = label_d;
739  printf(" %s", int2str(k).c_str());
740  }
741  printf("\n");
742  }
743 #ifdef DEBUG_SOLVEPROB
744  else if (ioflag > 1) {
745
746  updateMFSolnSP(XMolSolnSP);
747  printf("\n\t Weighted norm of update = %10.4e\n", update_norm);
748
749  printf("\t Name Prod_Rate XMol Conc "
750  " Conc_Old wtConc");
751  if (damp < 1.0) {
752  printf(" UnDamped_Conc");
753  }
754  printf("\n");
755  printf("\t---------------------------------------------------------"
756  "-----------------------------\n");
757  int kindexSP = 0;
758  for (int isp = 0; isp < m_numSurfPhases; isp++) {
759  int nsp = m_nSpeciesSurfPhase[isp];
760  InterfaceKinetics* m_kin = m_objects[isp];
762  for (int k = 0; k < nsp; k++, kindexSP++) {
763  int kspIndex = m_kinSpecIndex[kindexSP];
764  nm = m_kin->kineticsSpeciesName(kspIndex);
765  printf("\t%-16s %10.3e %10.3e %10.3e %10.3e %10.3e ",
766  nm.c_str(),
767  m_numEqn1[kspIndex],
768  XMolSolnSP[kindexSP],
769  CSolnSP[kindexSP], CSolnSP[kindexSP]+damp*resid[kindexSP],
770  wtSpecies[kindexSP]);
771  if (damp < 1.0) {
772  printf("%10.4e ", CSolnSP[kindexSP]+(damp-1.0)*resid[kindexSP]);
773  if (label_d == kindexSP) {
774  printf(" Damp ");
775  }
776  }
777  if (label_t == kindexSP) {
778  printf(" Tctrl");
779  }
780  printf("\n");
781  }
782
783  }
784
785  printf("\t--------------------------------------------------------"
786  "------------------------------\n");
787  }
788 #endif
789 }
790
791 void solveProb::printFinal(int ioflag, doublereal damp, size_t label_d, size_t label_t,
792  doublereal inv_t, doublereal t_real, int iter,
793  doublereal update_norm, doublereal resid_norm,
794  doublereal netProdRateKinSpecies[], const doublereal CSolnSP[],
795  const doublereal resid[],
796  const doublereal wtSpecies[], const doublereal wtRes[],
797  size_t dim, bool do_time)
798 {
799  size_t i, k;
800  string nm;
801  if (ioflag == 1) {
802
803  printf("\tFIN%3d ", iter);
804  if (do_time) {
805  printf("%9.4e %9.4e ", t_real, 1.0/inv_t);
806  } else
807  for (i = 0; i < 22; i++) {
808  printf(" ");
809  }
810  if (damp < 1.0) {
811  printf("%9.4e ", damp);
812  } else
813  for (i = 0; i < 11; i++) {
814  printf(" ");
815  }
816  printf("%9.4e %9.4e", update_norm, resid_norm);
817  if (do_time) {
818  k = label_t;
819  printf(" %s", int2str(k).c_str());
820  } else {
821  for (i = 0; i < 16; i++) {
822  printf(" ");
823  }
824  }
825  if (label_d != npos) {
826  k = label_d;
827
828  printf(" %s", int2str(k).c_str());
829  }
830  printf(" -- success\n");
831  }
832 #ifdef DEBUG_SOLVEPROB
833  else if (ioflag > 1) {
834  printf("\n================================== FINAL RESULT ========="
835  "==================================================\n");
836
837  printf("\n Weighted norm of solution update = %10.4e\n", update_norm);
838  printf(" Weighted norm of residual update = %10.4e\n\n", resid_norm);
839
840  printf(" Name Prod_Rate XMol Conc "
841  " wtConc Resid Resid/wtResid wtResid");
842  if (damp < 1.0) {
843  printf(" UnDamped_Conc");
844  }
845  printf("\n");
846  printf("---------------------------------------------------------------"
847  "---------------------------------------------\n");
848
849  for (int k = 0; k < m_neq; k++, k++) {
850  printf("%-16s %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e",
851  nm.c_str(),
852  m_numEqn1[k],
853  XMolSolnSP[k],
854  CSolnSP[k],
855  wtSpecies[k],
856  resid[k],
857  resid[k]/wtRes[k], wtRes[k]);
858  if (damp < 1.0) {
859  printf("%10.4e ", CSolnSP[k]+(damp-1.0)*resid[k]);
860  if (label_d == k) {
861  printf(" Damp ");
862  }
863  }
864  if (label_t == k) {
865  printf(" Tctrl");
866  }
867  printf("\n");
868  }
869
870  printf("\n");
871  printf("==============================================================="
872  "============================================\n\n");
873  }
874 #endif
875 }
876
877 #ifdef DEBUG_SOLVEPROB
878 void solveProb::printIterationHeader(int ioflag, doublereal damp,
879  doublereal inv_t, doublereal t_real,
880  int iter, bool do_time)
881 {
882  if (ioflag > 1) {
883  printf("\n===============================Iteration %5d "
884  "=================================\n", iter);
885  if (do_time) {
886  printf(" Transient step with: Real Time_n-1 = %10.4e sec,", t_real);
887  printf(" Time_n = %10.4e sec\n", t_real + 1.0/inv_t);
888  printf(" Delta t = %10.4e sec", 1.0/inv_t);
889  } else {
891  }
892  if (damp < 1.0) {
893  printf(", Damping value = %10.4e\n", damp);
894  } else {
895  printf("\n");
896  }
897  }
898 }
899 #endif
900
901 void solveProb::setAtol(const doublereal atol[])
902 {
903  for (size_t k = 0; k < m_neq; k++, k++) {
904  m_atol[k] = atol[k];
905  }
906 }
907
908 void solveProb::setAtolConst(const doublereal atolconst)
909 {
910  for (size_t k = 0; k < m_neq; k++, k++) {
911  m_atol[k] = atolconst;
912  }
913 }
914
915 }
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
virtual void setBounds(const doublereal botBounds[], const doublereal topBounds[])
Set the bottom and top bounds on the solution vector.
Definition: solveProb.cpp:589
const int SOLVEPROB_INITIALIZE
Solution Methods.
Definition: solveProb.h:52
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
virtual doublereal calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[], size_t *label, size_t *label_old, doublereal *label_factor, int ioflag)
Calculate a conservative delta T to use in a pseudo-steady state algorithm.
Definition: solveProb.cpp:536
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:285
virtual doublereal calc_damping(doublereal x[], doublereal dxneg[], size_t dim, size_t *label)
This function calculates a damping factor for the Newton iteration update vector, dxneg...
Definition: solveProb.cpp:424
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
std::string kineticsSpeciesName(size_t k) const
Return the name of the kth species in the kinetics manager.
Definition: Kinetics.cpp:359
vector_fp m_wtSpecies
Weights for the species concentrations norm calculation.
Definition: solveProb.h:404
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:78
The class provides the wall clock timer in seconds.
Definition: clockWC.h:51
int solve(int ifunc, doublereal time_scale, doublereal reltol)
Main routine that actually calculates the pseudo steady state of the surface problem.
Definition: solveProb.cpp:63
virtual void printIteration(int ioflag, doublereal damp, size_t label_d, size_t label_t, doublereal inv_t, doublereal t_real, int iter, doublereal update_norm, doublereal resid_norm, doublereal netProdRate[], doublereal CSolnSP[], doublereal resid[], doublereal wtSpecies[], size_t dim, bool do_time)
Printing routine that gets called after every iteration.
Definition: solveProb.cpp:703
This file contains definitions for utility functions and text for modules, inputfiles, logs, textlogs, (see Input File Handling, Diagnostic Output, and Writing messages to the screen).
size_t m_neq
Total number of equations to solve in the implicit problem.
Definition: solveProb.h:350
vector_fp m_numEqn1
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:368
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:29
virtual void getNetProductionRates(doublereal *wdot)
Species net production rates [kmol/m^3/s or kmol/m^2/s].
Definition: Kinetics.cpp:501
ResidEval * m_residFunc
residual function pointer to be solved.
Definition: solveProb.h:344
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:97
vector_fp m_botBounds
Bottom bounds for the solution vector.
Definition: solveProb.h:442
vector_fp m_CSolnSPOld
Saved solution vector at the old time step.
Definition: solveProb.h:392
size_t surfacePhaseIndex()
This returns the integer index of the phase which has ThermoPhase type cSurf.
Definition: Kinetics.h:260
virtual void reportState(doublereal *const CSoln) const
Report the current state of the solution.
Definition: solveProb.cpp:370
virtual int getInitialConditions(const doublereal t0, doublereal *const y, doublereal *const ydot)
Fill in the initial conditions.
Definition: ResidEval.h:119
Virtual base class for DAE residual function evaluators.
Definition: ResidEval.h:34
size_t nPhases() const
The number of phases participating in the reaction mechanism.
Definition: Kinetics.h:225
A kinetics manager for heterogeneous reaction mechanisms.
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Definition: clockWC.cpp:35
int solve(doublereal *b, size_t nrhs=1, size_t ldb=0)
Solves the Ax = b system returning x in the b spot.
doublereal m_rtol
m_rtol is the relative error tolerance.
Definition: solveProb.h:356
vector_fp m_atol
m_atol is the absolute tolerance in real units.
Definition: solveProb.h:353
vector_fp m_resid
Residual for the surface problem.
Definition: solveProb.h:417
int factor()
Factors the A matrix, overwriting A.
std::string id() const
Return the string id for the phase.
Definition: Phase.cpp:147
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
vector_fp m_CSolnSP
Solution vector.
Definition: solveProb.h:380
virtual int nEquations() const =0
Return the number of equations in the equation system.
size_t nSpecies() const
Returns the number of species in the phase.
Definition: Phase.h:265
SquareMatrix m_Jac
Jacobian.
Definition: solveProb.h:430
virtual void fun_eval(doublereal *const resid, const doublereal *const CSolnSP, const doublereal *const CSolnSPOld, const bool do_time, const doublereal deltaT)
Main Function evaluation.
Definition: solveProb.cpp:375
virtual 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: solveProb.cpp:390
std::vector< doublereal * > m_JacCol
Vector of pointers to the top of the columns of the Jacobians.
Definition: solveProb.h:424
vector_fp m_CSolnSPInit
Saved initial solution vector.
Definition: solveProb.h:386
Header file for implicit nonlinear solver with the option of a pseudotransient (see Numerical Utiliti...
Contains declarations for string manipulation functions within Cantera.
virtual void calcWeights(doublereal wtSpecies[], doublereal wtResid[], const doublereal CSolnSP[])
Calculate the solution and residual weights.
Definition: solveProb.cpp:512
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
vector_fp m_netProductionRatesSave
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:365
virtual doublereal * ptrColumn(size_t j)
Return a pointer to the top of column j, columns are assumed to be contiguous in memory.
virtual void print_header(int ioflag, int ifunc, doublereal time_scale, doublereal reltol, doublereal netProdRate[])
Printing routine that gets called at the start of every invocation.
Definition: solveProb.cpp:606
vector_fp m_CSolnSave
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:374
vector_fp m_numEqn2
Temporary vector with length MAX(1, m_neq)
Definition: solveProb.h:371
virtual void printFinal(int ioflag, doublereal damp, size_t label_d, size_t label_t, doublereal inv_t, doublereal t_real, int iter, doublereal update_norm, doublereal resid_norm, doublereal netProdRateKinSpecies[], const doublereal CSolnSP[], const doublereal resid[], const doublereal wtSpecies[], const doublereal wtRes[], size_t dim, bool do_time)
Print a summary of the solution.
Definition: solveProb.cpp:791
vector_fp m_wtResid
Weights for the residual norm calculation.
Definition: solveProb.h:398
vector_fp m_topBounds
Top bounds for the solution vector.
Definition: solveProb.h:436
std::string speciesName(size_t k) const
Name of the species with index k.
Definition: Phase.cpp:272