Cantera  2.0
solveProb.cpp
1 /*
2  * @file: solveSP.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"
16 
17 /* Standard include files */
18 #include <cstdio>
19 #include <cstdlib>
20 #include <cmath>
21 #include <vector>
22 
23 using namespace std;
24 namespace Cantera
25 {
26 
27 /***************************************************************************
28  * STATIC ROUTINES DEFINED IN THIS FILE
29  ***************************************************************************/
30 
31 static doublereal calcWeightedNorm(const doublereal [], const doublereal dx[], size_t);
32 
33 //================================================================================================
34 // Main constructor
35 solveProb::solveProb(ResidEval* resid) :
36  m_residFunc(resid),
37  m_neq(0),
38  m_atol(0),
39  m_rtol(1.0E-4),
40  m_maxstep(1000),
41  m_ioflag(0)
42 {
44 
45  // Dimension solution vector
46  size_t dim1 = std::max<size_t>(1, m_neq);
47 
48  m_atol.resize(dim1, 1.0E-9);
49  m_netProductionRatesSave.resize(dim1, 0.0);
50  m_numEqn1.resize(dim1, 0.0);
51  m_numEqn2.resize(dim1, 0.0);
52  m_CSolnSave.resize(dim1, 0.0);
53  m_CSolnSP.resize(dim1, 0.0);
54  m_CSolnSPInit.resize(dim1, 0.0);
55  m_CSolnSPOld.resize(dim1, 0.0);
56  m_wtResid.resize(dim1, 0.0);
57  m_wtSpecies.resize(dim1, 0.0);
58  m_resid.resize(dim1, 0.0);
59  m_ipiv.resize(dim1, 0);
60  m_topBounds.resize(dim1, 1.0);
61  m_botBounds.resize(dim1, 0.0);
62 
63  m_Jac.resize(dim1, dim1, 0.0);
64  m_JacCol.resize(dim1, 0);
65  for (size_t k = 0; k < dim1; k++) {
66  m_JacCol[k] = m_Jac.ptrColumn(k);
67  }
68 
69 }
70 //================================================================================================
71 // Empty destructor
73 {
74 }
75 //================================================================================================
76 /*
77  * The following calculation is a Newton's method to
78  * get the surface fractions of the surface and bulk species by
79  * requiring that the
80  * surface species production rate = 0 and that the bulk fractions are
81  * proportional to their production rates.
82  */
83 int solveProb::solve(int ifunc, doublereal time_scale,
84  doublereal reltol)
85 {
86  doublereal EXTRA_ACCURACY = 0.001;
87  if (ifunc == SOLVEPROB_JACOBIAN) {
88  EXTRA_ACCURACY *= 0.001;
89  }
90  int info = 0;
91  size_t label_t = npos; /* Species IDs for time control */
92  size_t label_d; /* Species IDs for damping control */
93  size_t label_t_old = npos;
94  doublereal label_factor = 1.0;
95  int iter=0; // iteration number on numlinear solver
96  int iter_max=1000; // maximum number of nonlinear iterations
97  int nrhs=1;
98  doublereal deltaT = 1.0E-10; // Delta time step
99  doublereal damp=1.0, tmp;
100  // Weighted L2 norm of the residual. Currently, this is only
101  // used for IO purposes. It doesn't control convergence.
102  // Therefore, it is turned off when DEBUG_SOLVEPROB isn't defined.
103  doublereal resid_norm;
104  doublereal inv_t = 0.0;
105  doublereal t_real = 0.0, update_norm = 1.0E6;
106 
107  bool do_time = false, not_converged = true;
108 
109 #ifdef DEBUG_SOLVEPROB
110 #ifdef DEBUG_SOLVEPROB_TIME
111  doublereal t1;
112 #endif
113 #else
114  if (m_ioflag > 1) {
115  m_ioflag = 1;
116  }
117 #endif
118 
119 #ifdef DEBUG_SOLVEPROB
120 #ifdef DEBUG_SOLVEPROB_TIME
121  Cantera::clockWC wc;
122  if (m_ioflag) {
123  t1 = wc.secondsWC();
124  }
125 #endif
126 #endif
127 
128  /*
129  * Set the initial value of the do_time parameter
130  */
131  if (ifunc == SOLVEPROB_INITIALIZE || ifunc == SOLVEPROB_TRANSIENT) {
132  do_time = true;
133  }
134 
135  /*
136  * upload the initial conditions
137  */
139 
140  /*
141  * Store the initial guess in the soln vector,
142  * CSolnSP, and in an separate vector CSolnSPInit.
143  */
144  std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPInit.begin());
145 
146 
147 
148  if (m_ioflag) {
149  print_header(m_ioflag, ifunc, time_scale, reltol,
151  }
152 
153  /*
154  * Quick return when there isn't a surface problem to solve
155  */
156  if (m_neq == 0) {
157  not_converged = false;
158  update_norm = 0.0;
159  }
160 
161  /* ------------------------------------------------------------------
162  * Start of Newton's method
163  * ------------------------------------------------------------------
164  */
165  while (not_converged && iter < iter_max) {
166  iter++;
167  /*
168  * Store previous iteration's solution in the old solution vector
169  */
170  std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPOld.begin());
171 
172  /*
173  * Evaluate the largest surface species for each surface phase every
174  * 5 iterations.
175  */
176  // if (iter%5 == 4) {
177  // evalSurfLarge(DATA_PTR(m_CSolnSP));
178  // }
179 
180  /*
181  * Calculate the value of the time step
182  * - heuristics to stop large oscillations in deltaT
183  */
184  if (do_time) {
185  /* don't hurry increase in time step at the same time as damping */
186  if (damp < 1.0) {
187  label_factor = 1.0;
188  }
190  &label_t, &label_t_old, &label_factor, m_ioflag);
191  if (iter < 10) {
192  inv_t = tmp;
193  } else if (tmp > 2.0*inv_t) {
194  inv_t = 2.0*inv_t;
195  } else {
196  inv_t = tmp;
197  }
198 
199  /*
200  * Check end condition
201  */
202 
203  if (ifunc == SOLVEPROB_TRANSIENT) {
204  tmp = t_real + 1.0/inv_t;
205  if (tmp > time_scale) {
206  inv_t = 1.0/(time_scale - t_real);
207  }
208  }
209  } else {
210  /* make steady state calc a step of 1 million seconds to
211  prevent singular jacobians for some pathological cases */
212  inv_t = 1.0e-6;
213  }
214  deltaT = 1.0/inv_t;
215 
216  /*
217  * Call the routine to numerically evaluation the jacobian
218  * and residual for the current iteration.
219  */
221  DATA_PTR(m_CSolnSPOld), do_time, deltaT);
222 
223  /*
224  * Calculate the weights. Make sure the calculation is carried
225  * out on the first iteration.
226  */
227  if (iter%4 == 1) {
230  }
231 
232  /*
233  * Find the weighted norm of the residual
234  */
235  resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid), DATA_PTR(m_resid), m_neq);
236 
237 #ifdef DEBUG_SOLVEPROB
238  if (m_ioflag > 1) {
239  printIterationHeader(m_ioflag, damp, inv_t, t_real, iter, do_time);
240  /*
241  * Print out the residual and jacobian
242  */
243  printResJac(m_ioflag, m_neq, m_Jac, DATA_PTR(m_resid),
244  DATA_PTR(m_wtResid), resid_norm);
245  }
246 #endif
247 
248  /*
249  * Solve Linear system (with LAPACK). The solution is in resid[]
250  */
251 
252  ct_dgetrf(m_neq, m_neq, m_JacCol[0], m_neq, DATA_PTR(m_ipiv), info);
253  if (info==0) {
254  ct_dgetrs(ctlapack::NoTranspose, m_neq, nrhs, m_JacCol[0],
256  info);
257  }
258  /*
259  * Force convergence if residual is small to avoid
260  * "nan" results from the linear solve.
261  */
262  else {
263  if (m_ioflag) {
264  printf("solveSurfSS: Zero pivot, assuming converged: %g (%d)\n",
265  resid_norm, info);
266  }
267  for (size_t jcol = 0; jcol < m_neq; jcol++) {
268  m_resid[jcol] = 0.0;
269  }
270 
271  /* print out some helpful info */
272  if (m_ioflag > 1) {
273  printf("-----\n");
274  printf("solveSurfProb: iter %d t_real %g delta_t %g\n\n",
275  iter,t_real, 1.0/inv_t);
276  printf("solveSurfProb: init guess, current concentration,"
277  "and prod rate:\n");
278 
279  printf("-----\n");
280  }
281  if (do_time) {
282  t_real += time_scale;
283  }
284 #ifdef DEBUG_SOLVEPROB
285  if (m_ioflag) {
286  printf("\nResidual is small, forcing convergence!\n");
287  }
288 #endif
289  }
290 
291  /*
292  * Calculate the Damping factor needed to keep all unknowns
293  * between 0 and 1, and not allow too large a change (factor of 2)
294  * in any unknown.
295  */
296 
297 
298  damp = calc_damping(DATA_PTR(m_CSolnSP), DATA_PTR(m_resid), m_neq, &label_d);
299 
300 
301  /*
302  * Calculate the weighted norm of the update vector
303  * Here, resid is the delta of the solution, in concentration
304  * units.
305  */
306  update_norm = calcWeightedNorm(DATA_PTR(m_wtSpecies),
308  /*
309  * Update the solution vector and real time
310  * Crop the concentrations to zero.
311  */
312  for (size_t irow = 0; irow < m_neq; irow++) {
313  m_CSolnSP[irow] -= damp * m_resid[irow];
314  }
315 
316 
317  if (do_time) {
318  t_real += damp/inv_t;
319  }
320 
321  if (m_ioflag) {
322  printIteration(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
323  update_norm, resid_norm,
326  DATA_PTR(m_wtSpecies), m_neq, do_time);
327  }
328 
329  if (ifunc == SOLVEPROB_TRANSIENT) {
330  not_converged = (t_real < time_scale);
331  } else {
332  if (do_time) {
333  if (t_real > time_scale ||
334  (resid_norm < 1.0e-7 &&
335  update_norm*time_scale/t_real < EXTRA_ACCURACY)) {
336  do_time = false;
337 #ifdef DEBUG_SOLVEPROB
338  if (m_ioflag > 1) {
339  printf("\t\tSwitching to steady solve.\n");
340  }
341 #endif
342  }
343  } else {
344  not_converged = ((update_norm > EXTRA_ACCURACY) ||
345  (resid_norm > EXTRA_ACCURACY));
346  }
347  }
348  } /* End of Newton's Method while statement */
349 
350  /*
351  * End Newton's method. If not converged, print error message and
352  * recalculate sdot's at equal site fractions.
353  */
354  if (not_converged) {
355  if (m_ioflag) {
356  printf("#$#$#$# Error in solveProb $#$#$#$ \n");
357  printf("Newton iter on surface species did not converge, "
358  "update_norm = %e \n", update_norm);
359  printf("Continuing anyway\n");
360  }
361  }
362 #ifdef DEBUG_SOLVEPROB
363 #ifdef DEBUG_SOLVEPROB_TIME
364  if (m_ioflag) {
365  printf("\nEnd of solve, time used: %e\n", wc.secondsWC()-t1);
366  }
367 #endif
368 #endif
369 
370  /*
371  * Decide on what to return in the solution vector
372  * - right now, will always return the last solution
373  * no matter how bad
374  */
375  if (m_ioflag) {
377  false, deltaT);
378  resid_norm = calcWeightedNorm(DATA_PTR(m_wtResid),
380  printFinal(m_ioflag, damp, label_d, label_t, inv_t, t_real, iter,
381  update_norm, resid_norm, DATA_PTR(m_netProductionRatesSave),
384  DATA_PTR(m_wtResid), m_neq, do_time);
385  }
386 
387  /*
388  * Return with the appropriate flag
389  */
390  if (update_norm > 1.0) {
391  return -1;
392  }
393  return 0;
394 }
395 //================================================================================================
396 /*
397  * Update the surface states of the surface phases.
398  */
399 void solveProb::reportState(doublereal* const CSolnSP) const
400 {
401  std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), CSolnSP);
402 }
403 //================================================================================================
404 /*
405  * This calculates the net production rates of all species
406  *
407  * This calculates the function eval.
408  * (should switch to special_species formulation for sum condition)
409  *
410  * @internal
411  * This routine uses the m_numEqn1 and m_netProductionRatesSave vectors
412  * as temporary internal storage.
413  */
414 void solveProb::fun_eval(doublereal* const resid, const doublereal* const CSoln,
415  const doublereal* const CSolnOld, const bool do_time,
416  const doublereal deltaT)
417 {
418  if (do_time) {
419  m_residFunc->evalSimpleTD(0.0, CSoln, CSolnOld, deltaT, resid);
420  } else {
421  m_residFunc->evalSS(0.0, CSoln, resid);
422  }
423 }
424 //================================================================================================
425 /*
426  * Calculate the Jacobian and residual
427  *
428  * @internal
429  * This routine uses the m_numEqn2 vector
430  * as temporary internal storage.
431  */
432 void solveProb::resjac_eval(std::vector<doublereal*> &JacCol,
433  doublereal resid[], doublereal CSoln[],
434  const doublereal CSolnOld[], const bool do_time,
435  const doublereal deltaT)
436 {
437  doublereal dc, cSave, sd;
438  doublereal* col_j;
439  /*
440  * Calculate the residual
441  */
442  fun_eval(resid, CSoln, CSolnOld, do_time, deltaT);
443  /*
444  * Now we will look over the columns perturbing each unknown.
445  */
446 
447  for (size_t kCol = 0; kCol < m_neq; kCol++) {
448  cSave = CSoln[kCol];
449  sd = fabs(cSave) + fabs(CSoln[kCol]) + m_atol[kCol] * 1.0E6;
450  if (sd < 1.0E-200) {
451  sd = 1.0E-4;
452  }
453  dc = std::max(1.0E-11 * sd, fabs(cSave) * 1.0E-6);
454  CSoln[kCol] += dc;
455  fun_eval(DATA_PTR(m_numEqn2), CSoln, CSolnOld, do_time, deltaT);
456  col_j = JacCol[kCol];
457  for (size_t i = 0; i < m_neq; i++) {
458  col_j[i] = (m_numEqn2[i] - resid[i])/dc;
459  }
460  CSoln[kCol] = cSave;
461  }
462 
463 }
464 //================================================================================================
465 #define APPROACH 0.50
466 // This function calculates a damping factor for the Newton iteration update
467 // vector, dxneg, to insure that all solution components stay within prescribed bounds
468 /*
469  * The default for this class is that all solution components are bounded between zero and one.
470  * this is because the original unknowns were mole fractions and surface site fractions.
471  *
472  * dxneg[] = negative of the update vector.
473  *
474  * The constant "APPROACH" sets the fraction of the distance to the boundary
475  * that the step can take. If the full step would not force any fraction
476  * outside of the bounds, then Newton's method is mostly allowed to operate normally.
477  * There is also some solution damping employed.
478  *
479  * @param x Vector of the current solution components
480  * @param dxneg Vector of the negative of the full solution update vector.
481  * @param dim Size of the solution vector
482  * @param label return int, stating which solution component caused the most damping.
483  */
484 doublereal solveProb::calc_damping(doublereal x[], doublereal dxneg[], size_t dim, size_t* label)
485 {
486  doublereal damp = 1.0, xnew, xtop, xbot;
487  static doublereal damp_old = 1.0;
488  *label = npos;
489 
490  for (size_t i = 0; i < dim; i++) {
491  doublereal topBounds = m_topBounds[i];
492  doublereal botBounds = m_botBounds[i];
493  /*
494  * Calculate the new suggested new value of x[i]
495  */
496  double delta_x = - dxneg[i];
497  xnew = x[i] - damp * dxneg[i];
498 
499  /*
500  * Calculate the allowed maximum and minimum values of x[i]
501  * - Only going to allow x[i] to converge to the top and bottom bounds by a
502  * single order of magnitude at one time
503  */
504  bool canCrossOrigin = false;
505  if (topBounds > 0.0 && botBounds < 0.0) {
506  canCrossOrigin = true;
507  }
508 
509  xtop = topBounds - 0.1 * fabs(topBounds - x[i]);
510 
511  xbot = botBounds + 0.1 * fabs(x[i] - botBounds);
512 
513  if (xnew > xtop) {
514  damp = - APPROACH * (xtop - x[i]) / dxneg[i];
515  *label = i;
516  } else if (xnew < xbot) {
517  damp = APPROACH * (x[i] - xbot) / dxneg[i];
518  *label = i;
519  }
520  // else if (fabs(xnew) > 2.0*MAX(fabs(x[i]), 1.0E-10)) {
521  // damp = 0.5 * MAX(fabs(x[i]), 1.0E-9)/ fabs(xnew);
522  // *label = i;
523  // }
524  double denom = fabs(x[i]) + 1.0E5 * m_atol[i];
525  if ((fabs(delta_x) / denom) > 0.3) {
526  double newdamp = 0.3 * denom / fabs(delta_x);
527  if (canCrossOrigin) {
528  if (xnew * x[i] < 0.0) {
529  if (fabs(x[i]) < 1.0E8 * m_atol[i]) {
530  newdamp = 2.0 * fabs(x[i]) / fabs(delta_x);
531  }
532  }
533  }
534  damp = std::min(damp, newdamp);
535  }
536 
537  }
538 
539  /*
540  * Only allow the damping parameter to increase by a factor of three each
541  * iteration. Heuristic to avoid oscillations in the value of damp
542  */
543  if (damp > damp_old*3) {
544  damp = damp_old*3;
545  *label = npos;
546  }
547 
548  /*
549  * Save old value of the damping parameter for use
550  * in subsequent calls.
551  */
552  damp_old = damp;
553  return damp;
554 
555 }
556 #undef APPROACH
557 //================================================================================================
558 /*
559  * This function calculates the norm of an update, dx[],
560  * based on the weighted values of x.
561  */
562 static doublereal calcWeightedNorm(const doublereal wtX[], const doublereal dx[], size_t dim)
563 {
564  doublereal norm = 0.0;
565  doublereal tmp;
566  if (dim == 0) {
567  return 0.0;
568  }
569  for (size_t i = 0; i < dim; i++) {
570  tmp = dx[i] / wtX[i];
571  norm += tmp * tmp;
572  }
573  return (sqrt(norm/dim));
574 }
575 //================================================================================================
576 /*
577  * Calculate the weighting factors for norms wrt both the species
578  * concentration unknowns and the residual unknowns.
579  *
580  */
581 void solveProb::calcWeights(doublereal wtSpecies[], doublereal wtResid[],
582  const doublereal CSoln[])
583 {
584  /*
585  * First calculate the weighting factor
586  */
587 
588  for (size_t k = 0; k < m_neq; k++) {
589  wtSpecies[k] = m_atol[k] + m_rtol * fabs(CSoln[k]);
590  }
591  /*
592  * Now do the residual Weights. Since we have the Jacobian, we
593  * will use it to generate a number based on the what a significant
594  * change in a solution variable does to each residual.
595  * This is a row sum scale operation.
596  */
597  for (size_t k = 0; k < m_neq; k++) {
598  wtResid[k] = 0.0;
599  for (size_t jcol = 0; jcol < m_neq; jcol++) {
600  wtResid[k] += fabs(m_Jac(k,jcol) * wtSpecies[jcol]);
601  }
602  }
603 }
604 //================================================================================================
605 /*
606  * This routine calculates a pretty conservative 1/del_t based
607  * on MAX_i(sdot_i/(X_i*SDen0)). This probably guarantees
608  * diagonal dominance.
609  *
610  * Small surface fractions are allowed to intervene in the del_t
611  * determination, no matter how small. This may be changed.
612  * Now minimum changed to 1.0e-12,
613  *
614  * Maximum time step set to time_scale.
615  */
616 doublereal solveProb::
617 calc_t(doublereal netProdRateSolnSP[], doublereal Csoln[],
618  size_t* label, size_t* label_old, doublereal* label_factor, int ioflag)
619 {
620  doublereal tmp, inv_timeScale=0.0;
621  for (size_t k = 0; k < m_neq; k++) {
622  if (Csoln[k] <= 1.0E-10) {
623  tmp = 1.0E-10;
624  } else {
625  tmp = Csoln[k];
626  }
627  tmp = fabs(netProdRateSolnSP[k]/ tmp);
628 
629 
630  if (netProdRateSolnSP[k]> 0.0) {
631  tmp /= 100.;
632  }
633  if (tmp > inv_timeScale) {
634  inv_timeScale = tmp;
635  *label = k;
636  }
637  }
638 
639  /*
640  * Increase time step exponentially as same species repeatedly
641  * controls time step
642  */
643  if (*label == *label_old) {
644  *label_factor *= 1.5;
645  } else {
646  *label_old = *label;
647  *label_factor = 1.0;
648  }
649  inv_timeScale = inv_timeScale / *label_factor;
650 #ifdef DEBUG_SOLVEPROB
651  if (ioflag > 1) {
652  if (*label_factor > 1.0) {
653  printf("Delta_t increase due to repeated controlling species = %e\n",
654  *label_factor);
655  }
656  int kkin = m_kinSpecIndex[*label];
657 
658  string sn = " "
659  printf("calc_t: spec=%d(%s) sf=%e pr=%e dt=%e\n",
660  *label, sn.c_str(), XMolSolnSP[*label],
661  netProdRateSolnSP[*label], 1.0/inv_timeScale);
662  }
663 #endif
664 
665  return (inv_timeScale);
666 
667 }
668 //====================================================================================================================
669 // Set the bottom and top bounds on the solution vector
670 /*
671  * The default is for the bottom is 0.0, while the default for the top is 1.0
672  *
673  * @param botBounds Vector of bottom bounds
674  * @param topBounds vector of top bounds
675  */
676 void solveProb::setBounds(const doublereal botBounds[], const doublereal topBounds[])
677 {
678  for (size_t k = 0; k < m_neq; k++) {
679  m_botBounds[k] = botBounds[k];
680  m_topBounds[k] = topBounds[k];
681  }
682 }
683 //====================================================================================================================
684 /*
685  * printResJac(): prints out the residual and Jacobian.
686  *
687  */
688 #ifdef DEBUG_SOLVEPROB
689 void solveProb::printResJac(int ioflag, int neq, const Array2D& Jac,
690  doublereal resid[], doublereal wtRes[],
691  doublereal norm)
692 {
693 
694 }
695 #endif
696 //================================================================================================
697 /*
698  * Optional printing at the start of the solveProb problem
699  */
700 void solveProb::print_header(int ioflag, int ifunc, doublereal time_scale,
701  doublereal reltol,
702  doublereal netProdRate[])
703 {
704  int damping = 1;
705  if (ioflag) {
706  printf("\n================================ SOLVEPROB CALL SETUP "
707  "========================================\n");
708  if (ifunc == SOLVEPROB_INITIALIZE) {
709  printf("\n SOLVEPROB Called with Initialization turned on\n");
710  printf(" Time scale input = %9.3e\n", time_scale);
711  } else if (ifunc == SOLVEPROB_RESIDUAL) {
712  printf("\n SOLVEPROB Called to calculate steady state residual\n");
713  printf(" from a good initial guess\n");
714  } else if (ifunc == SOLVEPROB_JACOBIAN) {
715  printf("\n SOLVEPROB Called to calculate steady state jacobian\n");
716  printf(" from a good initial guess\n");
717  } else if (ifunc == SOLVEPROB_TRANSIENT) {
718  printf("\n SOLVEPROB Called to integrate surface in time\n");
719  printf(" for a total of %9.3e sec\n", time_scale);
720  } else {
721  fprintf(stderr,"Unknown ifunc flag = %d\n", ifunc);
722  exit(EXIT_FAILURE);
723  }
724 
725 
726 
727  if (damping) {
728  printf(" Damping is ON \n");
729  } else {
730  printf(" Damping is OFF \n");
731  }
732 
733  printf(" Reltol = %9.3e, Abstol = %9.3e\n", reltol, m_atol[0]);
734  }
735 
736  /*
737  * Print out the initial guess
738  */
739 #ifdef DEBUG_SOLVEPROB
740  if (ioflag > 1) {
741  printf("\n================================ INITIAL GUESS "
742  "========================================\n");
743  int kindexSP = 0;
744  for (int isp = 0; isp < m_numSurfPhases; isp++) {
745  InterfaceKinetics* m_kin = m_objects[isp];
746  int surfIndex = m_kin->surfacePhaseIndex();
747  int nPhases = m_kin->nPhases();
748  m_kin->getNetProductionRates(netProdRate);
749  updateMFKinSpecies(XMolKinSpecies, isp);
750 
751  printf("\n IntefaceKinetics Object # %d\n\n", isp);
752 
753  printf("\t Number of Phases = %d\n", nPhases);
754  printf("\t Phase:SpecName Prod_Rate MoleFraction kindexSP\n");
755  printf("\t -------------------------------------------------------"
756  "----------\n");
757 
758  int kspindex = 0;
759  bool inSurfacePhase = false;
760  for (int ip = 0; ip < nPhases; ip++) {
761  if (ip == surfIndex) {
762  inSurfacePhase = true;
763  } else {
764  inSurfacePhase = false;
765  }
766  ThermoPhase& THref = m_kin->thermo(ip);
767  int nsp = THref.nSpecies();
768  string pname = THref.id();
769  for (int k = 0; k < nsp; k++) {
770  string sname = THref.speciesName(k);
771  string cname = pname + ":" + sname;
772  if (inSurfacePhase) {
773  printf("\t %-24s %10.3e %10.3e %d\n", cname.c_str(),
774  netProdRate[kspindex], XMolKinSpecies[kspindex],
775  kindexSP);
776  kindexSP++;
777  } else {
778  printf("\t %-24s %10.3e %10.3e\n", cname.c_str(),
779  netProdRate[kspindex], XMolKinSpecies[kspindex]);
780  }
781  kspindex++;
782  }
783  }
784  printf("=========================================================="
785  "=================================\n");
786  }
787  }
788 #endif
789  if (ioflag == 1) {
790  printf("\n\n\t Iter Time Del_t Damp DelX "
791  " Resid Name-Time Name-Damp\n");
792  printf("\t -----------------------------------------------"
793  "------------------------------------\n");
794  }
795 }
796 //================================================================================================
797 void solveProb::printIteration(int ioflag, doublereal damp, size_t label_d,
798  size_t label_t,
799  doublereal inv_t, doublereal t_real, int iter,
800  doublereal update_norm, doublereal resid_norm,
801  doublereal netProdRate[], doublereal CSolnSP[],
802  doublereal resid[],
803  doublereal wtSpecies[], size_t dim, bool do_time)
804 {
805  size_t i, k;
806  string nm;
807  if (ioflag == 1) {
808 
809  printf("\t%6d ", iter);
810  if (do_time) {
811  printf("%9.4e %9.4e ", t_real, 1.0/inv_t);
812  } else
813  for (i = 0; i < 22; i++) {
814  printf(" ");
815  }
816  if (damp < 1.0) {
817  printf("%9.4e ", damp);
818  } else
819  for (i = 0; i < 11; i++) {
820  printf(" ");
821  }
822  printf("%9.4e %9.4e", update_norm, resid_norm);
823  if (do_time) {
824  k = label_t;
825  printf(" %s", int2str(k).c_str());
826  } else {
827  for (i = 0; i < 16; i++) {
828  printf(" ");
829  }
830  }
831  if (label_d != npos) {
832  k = label_d;
833  printf(" %s", int2str(k).c_str());
834  }
835  printf("\n");
836  }
837 #ifdef DEBUG_SOLVEPROB
838  else if (ioflag > 1) {
839 
840  updateMFSolnSP(XMolSolnSP);
841  printf("\n\t Weighted norm of update = %10.4e\n", update_norm);
842 
843  printf("\t Name Prod_Rate XMol Conc "
844  " Conc_Old wtConc");
845  if (damp < 1.0) {
846  printf(" UnDamped_Conc");
847  }
848  printf("\n");
849  printf("\t---------------------------------------------------------"
850  "-----------------------------\n");
851  int kindexSP = 0;
852  for (int isp = 0; isp < m_numSurfPhases; isp++) {
853  int nsp = m_nSpeciesSurfPhase[isp];
854  InterfaceKinetics* m_kin = m_objects[isp];
855  //int surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp];
857  for (int k = 0; k < nsp; k++, kindexSP++) {
858  int kspIndex = m_kinSpecIndex[kindexSP];
859  nm = m_kin->kineticsSpeciesName(kspIndex);
860  printf("\t%-16s %10.3e %10.3e %10.3e %10.3e %10.3e ",
861  nm.c_str(),
862  m_numEqn1[kspIndex],
863  XMolSolnSP[kindexSP],
864  CSolnSP[kindexSP], CSolnSP[kindexSP]+damp*resid[kindexSP],
865  wtSpecies[kindexSP]);
866  if (damp < 1.0) {
867  printf("%10.4e ", CSolnSP[kindexSP]+(damp-1.0)*resid[kindexSP]);
868  if (label_d == kindexSP) {
869  printf(" Damp ");
870  }
871  }
872  if (label_t == kindexSP) {
873  printf(" Tctrl");
874  }
875  printf("\n");
876  }
877 
878  }
879 
880  printf("\t--------------------------------------------------------"
881  "------------------------------\n");
882  }
883 #endif
884 } /* printIteration */
885 
886 //================================================================================================
887 void solveProb::printFinal(int ioflag, doublereal damp, size_t label_d, size_t label_t,
888  doublereal inv_t, doublereal t_real, int iter,
889  doublereal update_norm, doublereal resid_norm,
890  doublereal netProdRateKinSpecies[], const doublereal CSolnSP[],
891  const doublereal resid[],
892  const doublereal wtSpecies[], const doublereal wtRes[],
893  size_t dim, bool do_time)
894 {
895  size_t i, k;
896  string nm;
897  if (ioflag == 1) {
898 
899  printf("\tFIN%3d ", iter);
900  if (do_time) {
901  printf("%9.4e %9.4e ", t_real, 1.0/inv_t);
902  } else
903  for (i = 0; i < 22; i++) {
904  printf(" ");
905  }
906  if (damp < 1.0) {
907  printf("%9.4e ", damp);
908  } else
909  for (i = 0; i < 11; i++) {
910  printf(" ");
911  }
912  printf("%9.4e %9.4e", update_norm, resid_norm);
913  if (do_time) {
914  k = label_t;
915  printf(" %s", int2str(k).c_str());
916  } else {
917  for (i = 0; i < 16; i++) {
918  printf(" ");
919  }
920  }
921  if (label_d != npos) {
922  k = label_d;
923 
924  printf(" %s", int2str(k).c_str());
925  }
926  printf(" -- success\n");
927  }
928 #ifdef DEBUG_SOLVEPROB
929  else if (ioflag > 1) {
930 
931 
932  printf("\n================================== FINAL RESULT ========="
933  "==================================================\n");
934 
935  printf("\n Weighted norm of solution update = %10.4e\n", update_norm);
936  printf(" Weighted norm of residual update = %10.4e\n\n", resid_norm);
937 
938  printf(" Name Prod_Rate XMol Conc "
939  " wtConc Resid Resid/wtResid wtResid");
940  if (damp < 1.0) {
941  printf(" UnDamped_Conc");
942  }
943  printf("\n");
944  printf("---------------------------------------------------------------"
945  "---------------------------------------------\n");
946 
947  for (int k = 0; k < m_neq; k++, k++) {
948  printf("%-16s %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e",
949  nm.c_str(),
950  m_numEqn1[k],
951  XMolSolnSP[k],
952  CSolnSP[k],
953  wtSpecies[k],
954  resid[k],
955  resid[k]/wtRes[k], wtRes[k]);
956  if (damp < 1.0) {
957  printf("%10.4e ", CSolnSP[k]+(damp-1.0)*resid[k]);
958  if (label_d == k) {
959  printf(" Damp ");
960  }
961  }
962  if (label_t == k) {
963  printf(" Tctrl");
964  }
965  printf("\n");
966  }
967 
968  printf("\n");
969  printf("==============================================================="
970  "============================================\n\n");
971  }
972 #endif
973 }
974 //================================================================================================
975 #ifdef DEBUG_SOLVEPROB
976 void solveProb::
977 printIterationHeader(int ioflag, doublereal damp,doublereal inv_t, doublereal t_real,
978  int iter, bool do_time)
979 {
980  if (ioflag > 1) {
981  printf("\n===============================Iteration %5d "
982  "=================================\n", iter);
983  if (do_time) {
984  printf(" Transient step with: Real Time_n-1 = %10.4e sec,", t_real);
985  printf(" Time_n = %10.4e sec\n", t_real + 1.0/inv_t);
986  printf(" Delta t = %10.4e sec", 1.0/inv_t);
987  } else {
988  printf(" Steady Solve ");
989  }
990  if (damp < 1.0) {
991  printf(", Damping value = %10.4e\n", damp);
992  } else {
993  printf("\n");
994  }
995  }
996 }
997 #endif
998 //================================================================================================
999 void solveProb::setAtol(const doublereal atol[])
1000 {
1001  for (size_t k = 0; k < m_neq; k++, k++) {
1002  m_atol[k] = atol[k];
1003  }
1004 }
1005 //================================================================================================
1006 void solveProb::setAtolConst(const doublereal atolconst)
1007 {
1008  for (size_t k = 0; k < m_neq; k++, k++) {
1009  m_atol[k] = atolconst;
1010  }
1011 }
1012 //================================================================================================
1013 
1014 
1015 }