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