Cantera  2.0
newton_utils.cpp
Go to the documentation of this file.
1 /**
2  * @file newton_utils.cpp
3  */
4
5 #include "cantera/base/ct_defs.h"
7
8 #include <cstdio>
9
10 using namespace std;
11
12 namespace Cantera
13 {
14
15 class Indx
16 {
17 public:
18  Indx(size_t nv, size_t np) : m_nv(nv), m_np(np) {}
19  size_t m_nv, m_np;
20  size_t operator()(size_t m, size_t j) {
21  return j*m_nv + m;
22  }
23 };
24
25
26 /**
27  * Return a damping coefficient that keeps the solution after taking one
28  * Newton step between specified lower and upper bounds. This function only
29  * considers one domain.
30  */
31 doublereal bound_step(const doublereal* x, const doublereal* step,
32  Domain1D& r, int loglevel)
33 {
34
35  char buf[100];
36  size_t np = r.nPoints();
37  size_t nv = r.nComponents();
38  Indx index(nv, np);
39  doublereal above, below, val, newval;
40  size_t m, j;
41  doublereal fbound = 1.0;
42  bool wroteTitle = false;
43  for (m = 0; m < nv; m++) {
44  above = r.upperBound(m);
45  below = r.lowerBound(m);
46
47  for (j = 0; j < np; j++) {
48  val = x[index(m,j)];
49  if (loglevel > 0) {
50  if (val > above + 1.0e-12 || val < below - 1.0e-12) {
51  sprintf(buf, "domain %s: %20s(%s) = %10.3e (%10.3e, %10.3e)\n",
52  int2str(r.domainIndex()).c_str(),
53  r.componentName(m).c_str(), int2str(j).c_str(),
54  val, below, above);
55  writelog(string("\nERROR: solution out of bounds.\n")+buf);
56  }
57  }
58
59  newval = val + step[index(m,j)];
60
61  if (newval > above) {
62  fbound = std::max(0.0, std::min(fbound,
63  (above - val)/(newval - val)));
64  } else if (newval < below) {
65  fbound = std::min(fbound, (val - below)/(val - newval));
66  }
67
68  if (loglevel > 1 && (newval > above || newval < below)) {
69  if (!wroteTitle) {
70  writelog("\nNewton step takes solution out of bounds.\n\n");
71  sprintf(buf," %12s %12s %4s %10s %10s %10s %10s\n",
72  "domain","component","pt","value","step","min","max");
73  wroteTitle = true;
74  writelog(buf);
75  }
76  sprintf(buf, " %4s %12s %4s %10.3e %10.3e %10.3e %10.3e\n",
77  int2str(r.domainIndex()).c_str(),
78  r.componentName(m).c_str(), int2str(j).c_str(),
79  val, step[index(m,j)], below, above);
80  writelog(buf);
81  }
82  }
83  }
84  return fbound;
85 }
86
87
88
89 /**
90  * This function computes the square of a weighted norm of a step
91  * vector for one domain.
92  *
93  * @param x Solution vector for this domain.
94  * @param step Newton step vector for this domain.
95  * @param r Object representing the domain. Used to get tolerances,
96  * number of components, and number of points.
97  *
98  * The return value is
99  * \f[
100  * \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2
101  * \f]
102  * where the error weight for solution component \f$n\f$ is given by
103  * \f[
104  * w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
105  * \f]
106  * Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for
107  * component n, and multiplies the average magnitude of
108  * solution component n in the domain. The second term,
109  * \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component
110  * n.
111  *
112  */
113 doublereal norm_square(const doublereal* x,
114  const doublereal* step, Domain1D& r)
115 {
116  doublereal f, ewt, esum, sum = 0.0;
117  size_t n, j;
118  doublereal f2max = 0.0;
119  size_t nv = r.nComponents();
120  size_t np = r.nPoints();
121
122  for (n = 0; n < nv; n++) {
123  esum = 0.0;
124  for (j = 0; j < np; j++) {
125  esum += fabs(x[nv*j + n]);
126  }
127  ewt = r.rtol(n)*esum/np + r.atol(n);
128  for (j = 0; j < np; j++) {
129  f = step[nv*j + n]/ewt;
130  sum += f*f;
131  if (f*f > f2max) {
132  f2max = f*f;
133  }
134  }
135  }
136  return sum;
137 }
138 }