Cantera  2.5.1
MultiNewton.cpp
Go to the documentation of this file.
1 //! @file MultiNewton.cpp Damped Newton solver for 1D multi-domain problems
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
8 
9 #include <ctime>
10 
11 using namespace std;
12 
13 namespace Cantera
14 {
15 
16 // unnamed-namespace for local helpers
17 namespace
18 {
19 
20 class Indx
21 {
22 public:
23  Indx(size_t nv, size_t np) : m_nv(nv), m_np(np) {}
24  size_t m_nv, m_np;
25  size_t operator()(size_t m, size_t j) {
26  return j*m_nv + m;
27  }
28 };
29 
30 /**
31  * Return a damping coefficient that keeps the solution after taking one
32  * Newton step between specified lower and upper bounds. This function only
33  * considers one domain.
34  */
35 doublereal bound_step(const doublereal* x, const doublereal* step,
36  Domain1D& r, int loglevel)
37 {
38  size_t np = r.nPoints();
39  size_t nv = r.nComponents();
40  Indx index(nv, np);
41  doublereal fbound = 1.0;
42  bool wroteTitle = false;
43  for (size_t m = 0; m < nv; m++) {
44  double above = r.upperBound(m);
45  double below = r.lowerBound(m);
46 
47  for (size_t j = 0; j < np; j++) {
48  double val = x[index(m,j)];
49  if (loglevel > 0 && (val > above + 1.0e-12 || val < below - 1.0e-12)) {
50  writelog("\nERROR: solution out of bounds.\n");
51  writelog("domain {:d}: {:>20s}({:d}) = {:10.3e} ({:10.3e}, {:10.3e})\n",
52  r.domainIndex(), r.componentName(m), j, val, below, above);
53  }
54 
55  double newval = val + step[index(m,j)];
56 
57  if (newval > above) {
58  fbound = std::max(0.0, std::min(fbound,
59  (above - val)/(newval - val)));
60  } else if (newval < below) {
61  fbound = std::min(fbound, (val - below)/(val - newval));
62  }
63 
64  if (loglevel > 1 && (newval > above || newval < below)) {
65  if (!wroteTitle) {
66  writelog("\nNewton step takes solution out of bounds.\n\n");
67  writelog(" {:>12s} {:>12s} {:>4s} {:>10s} {:>10s} {:>10s} {:>10s}\n",
68  "domain","component","pt","value","step","min","max");
69  wroteTitle = true;
70  }
71  writelog(" {:4d} {:>12s} {:4d} {:10.3e} {:10.3e} {:10.3e} {:10.3e}\n",
72  r.domainIndex(), r.componentName(m), j,
73  val, step[index(m,j)], below, above);
74  }
75  }
76  }
77  return fbound;
78 }
79 
80 /**
81  * This function computes the square of a weighted norm of a step vector for one
82  * domain.
83  *
84  * @param x Solution vector for this domain.
85  * @param step Newton step vector for this domain.
86  * @param r Object representing the domain. Used to get tolerances,
87  * number of components, and number of points.
88  *
89  * The return value is
90  * \f[
91  * \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2
92  * \f]
93  * where the error weight for solution component \f$n\f$ is given by
94  * \f[
95  * w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
96  * \f]
97  * Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for component n,
98  * and multiplies the average magnitude of solution component n in the domain.
99  * The second term, \f$\epsilon_{a,n}\f$, is the absolute error tolerance for
100  * component n.
101  */
102 doublereal norm_square(const doublereal* x,
103  const doublereal* step, Domain1D& r)
104 {
105  double sum = 0.0;
106  doublereal f2max = 0.0;
107  size_t nv = r.nComponents();
108  size_t np = r.nPoints();
109 
110  for (size_t n = 0; n < nv; n++) {
111  double esum = 0.0;
112  for (size_t j = 0; j < np; j++) {
113  esum += fabs(x[nv*j + n]);
114  }
115  double ewt = r.rtol(n)*esum/np + r.atol(n);
116  for (size_t j = 0; j < np; j++) {
117  double f = step[nv*j + n]/ewt;
118  sum += f*f;
119  f2max = std::max(f*f, f2max);
120  }
121  }
122  return sum;
123 }
124 
125 } // end unnamed-namespace
126 
127 
128 // constants
129 const doublereal DampFactor = sqrt(2.0);
130 const size_t NDAMP = 7;
131 
132 // ---------------- MultiNewton methods ----------------
133 
134 MultiNewton::MultiNewton(int sz)
135  : m_maxAge(5)
136 {
137  m_n = sz;
138  m_elapsed = 0.0;
139 }
140 
141 void MultiNewton::resize(size_t sz)
142 {
143  m_n = sz;
144  m_x.resize(m_n);
145  m_stp.resize(m_n);
146  m_stp1.resize(m_n);
147 }
148 
149 doublereal MultiNewton::norm2(const doublereal* x,
150  const doublereal* step, OneDim& r) const
151 {
152  double sum = 0.0;
153  size_t nd = r.nDomains();
154  for (size_t n = 0; n < nd; n++) {
155  double f = norm_square(x + r.start(n), step + r.start(n), r.domain(n));
156  sum += f;
157  }
158  sum /= r.size();
159  return sqrt(sum);
160 }
161 
162 void MultiNewton::step(doublereal* x, doublereal* step,
163  OneDim& r, MultiJac& jac, int loglevel)
164 {
165  r.eval(npos, x, step);
166  for (size_t n = 0; n < r.size(); n++) {
167  step[n] = -step[n];
168  }
169 
170  try {
171  jac.solve(step, step);
172  } catch (CanteraError&) {
173  int iok = jac.info() - 1;
174  if (iok >= 0) {
175  size_t nd = r.nDomains();
176  size_t n;
177  for (n = nd-1; n != npos; n--) {
178  if (iok >= static_cast<int>(r.start(n))) {
179  break;
180  }
181  }
182  Domain1D& dom = r.domain(n);
183  size_t offset = iok - r.start(n);
184  size_t pt = offset/dom.nComponents();
185  size_t comp = offset - pt*dom.nComponents();
186  throw CanteraError("MultiNewton::step",
187  "Jacobian is singular for domain {}, component {} at point {}\n"
188  "(Matrix row {})",
189  dom.id(), dom.componentName(comp), pt, iok);
190  } else {
191  throw;
192  }
193  }
194 }
195 
196 doublereal MultiNewton::boundStep(const doublereal* x0,
197  const doublereal* step0, const OneDim& r, int loglevel)
198 {
199  doublereal fbound = 1.0;
200  for (size_t i = 0; i < r.nDomains(); i++) {
201  fbound = std::min(fbound,
202  bound_step(x0 + r.start(i), step0 + r.start(i),
203  r.domain(i), loglevel));
204  }
205  return fbound;
206 }
207 
208 int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
209  doublereal* x1, doublereal* step1, doublereal& s1,
210  OneDim& r, MultiJac& jac, int loglevel, bool writetitle)
211 {
212  // write header
213  if (loglevel > 0 && writetitle) {
214  writelog("\n\nDamped Newton iteration:\n");
215  writeline('-', 65, false);
216 
217  writelog("\n{} {:>9s} {:>9s} {:>9s} {:>9s} {:>9s} {:>5s} {:>5s}\n",
218  "m","F_damp","F_bound","log10(ss)",
219  "log10(s0)","log10(s1)","N_jac","Age");
220  writeline('-', 65);
221  }
222 
223  // compute the weighted norm of the undamped step size step0
224  doublereal s0 = norm2(x0, step0, r);
225 
226  // compute the multiplier to keep all components in bounds
227  doublereal fbound = boundStep(x0, step0, r, loglevel-1);
228 
229  // if fbound is very small, then x0 is already close to the boundary and
230  // step0 points out of the allowed domain. In this case, the Newton
231  // algorithm fails, so return an error condition.
232  if (fbound < 1.e-10) {
233  debuglog("\nAt limits.\n", loglevel);
234  return -3;
235  }
236 
237  // ---------- Attempt damped step ----------
238 
239  // damping coefficient starts at 1.0
240  doublereal damp = 1.0;
241  size_t m;
242  for (m = 0; m < NDAMP; m++) {
243  double ff = fbound*damp;
244 
245  // step the solution by the damped step size
246  for (size_t j = 0; j < m_n; j++) {
247  x1[j] = ff*step0[j] + x0[j];
248  }
249 
250  // compute the next undamped step that would result if x1 is accepted
251  step(x1, step1, r, jac, loglevel-1);
252 
253  // compute the weighted norm of step1
254  s1 = norm2(x1, step1, r);
255 
256  // write log information
257  if (loglevel > 0) {
258  doublereal ss = r.ssnorm(x1,step1);
259  writelog("\n{:d} {:9.5f} {:9.5f} {:9.5f} {:9.5f} {:9.5f} {:4d} {:d}/{:d}",
260  m, damp, fbound, log10(ss+SmallNumber),
261  log10(s0+SmallNumber), log10(s1+SmallNumber),
262  jac.nEvals(), jac.age(), m_maxAge);
263  }
264 
265  // if the norm of s1 is less than the norm of s0, then accept this
266  // damping coefficient. Also accept it if this step would result in a
267  // converged solution. Otherwise, decrease the damping coefficient and
268  // try again.
269  if (s1 < 1.0 || s1 < s0) {
270  break;
271  }
272  damp /= DampFactor;
273  }
274 
275  // If a damping coefficient was found, return 1 if the solution after
276  // stepping by the damped step would represent a converged solution, and
277  // return 0 otherwise. If no damping coefficient could be found, return -2.
278  if (m < NDAMP) {
279  if (s1 > 1.0) {
280  return 0;
281  } else {
282  return 1;
283  }
284  } else {
285  return -2;
286  }
287 }
288 
289 int MultiNewton::solve(doublereal* x0, doublereal* x1,
290  OneDim& r, MultiJac& jac, int loglevel)
291 {
292  clock_t t0 = clock();
293  int m = 0;
294  bool forceNewJac = false;
295  doublereal s1=1.e30;
296 
297  copy(x0, x0 + m_n, &m_x[0]);
298 
299  bool frst = true;
300  doublereal rdt = r.rdt();
301  int j0 = jac.nEvals();
302  int nJacReeval = 0;
303 
304  while (true) {
305  // Check whether the Jacobian should be re-evaluated.
306  if (jac.age() > m_maxAge) {
307  if (loglevel > 0) {
308  writelog("\nMaximum Jacobian age reached ({})\n", m_maxAge);
309  }
310  forceNewJac = true;
311  }
312 
313  if (forceNewJac) {
314  r.eval(npos, &m_x[0], &m_stp[0], 0.0, 0);
315  jac.eval(&m_x[0], &m_stp[0], 0.0);
316  jac.updateTransient(rdt, r.transientMask().data());
317  forceNewJac = false;
318  }
319 
320  // compute the undamped Newton step
321  step(&m_x[0], &m_stp[0], r, jac, loglevel-1);
322 
323  // increment the Jacobian age
324  jac.incrementAge();
325 
326  // damp the Newton step
327  m = dampStep(&m_x[0], &m_stp[0], x1, &m_stp1[0], s1, r, jac, loglevel-1, frst);
328  if (loglevel == 1 && m >= 0) {
329  if (frst) {
330  writelog("\n\n {:>10s} {:>10s} {:>5s}",
331  "log10(ss)","log10(s1)","N_jac");
332  writelog("\n ------------------------------------");
333  }
334  doublereal ss = r.ssnorm(&m_x[0], &m_stp[0]);
335  writelog("\n {:10.4f} {:10.4f} {:d}",
336  log10(ss),log10(s1),jac.nEvals());
337  }
338  frst = false;
339 
340  // Successful step, but not converged yet. Take the damped step, and try
341  // again.
342  if (m == 0) {
343  copy(x1, x1 + m_n, m_x.begin());
344  } else if (m == 1) {
345  // convergence
346  if (rdt == 0) {
347  jac.setAge(0); // for efficient sensitivity analysis
348  }
349  break;
350  } else if (m < 0) {
351  // If dampStep fails, first try a new Jacobian if an old one was
352  // being used. If it was a new Jacobian, then return -1 to signify
353  // failure.
354  if (jac.age() > 1) {
355  forceNewJac = true;
356  if (nJacReeval > 3) {
357  break;
358  }
359  nJacReeval++;
360  debuglog("\nRe-evaluating Jacobian, since no damping "
361  "coefficient\ncould be found with this Jacobian.\n",
362  loglevel);
363  } else {
364  break;
365  }
366  }
367  }
368 
369  if (m < 0) {
370  copy(m_x.begin(), m_x.end(), x1);
371  }
372  if (m > 0 && jac.nEvals() == j0) {
373  m = 100;
374  }
375  m_elapsed += (clock() - t0)/(1.0*CLOCKS_PER_SEC);
376  return m;
377 }
378 
379 } // end namespace Cantera
int solve(const doublereal *const b, doublereal *const x)
Solve the matrix problem Ax = b.
Definition: BandMatrix.cpp:275
int info() const
LAPACK "info" flag after last factor/solve operation.
Definition: BandMatrix.h:261
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:61
Base class for one-dimensional domains.
Definition: Domain1D.h:39
size_t nComponents() const
Number of components at each grid point.
Definition: Domain1D.h:134
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition: Domain1D.cpp:56
Class MultiJac evaluates the Jacobian of a system of equations defined by a residual function supplie...
Definition: MultiJac.h:23
int nEvals() const
Number of Jacobian evaluations.
Definition: MultiJac.h:41
void setAge(int age)
Set the Jacobian age.
Definition: MultiJac.h:58
void eval(doublereal *x0, doublereal *resid0, double rdt)
Evaluate the Jacobian at x0.
Definition: MultiJac.cpp:43
void incrementAge()
Increment the Jacobian age.
Definition: MultiJac.h:51
int age() const
Number of times 'incrementAge' has been called since the last evaluation.
Definition: MultiJac.h:46
Container class for multiple-domain 1D problems.
Definition: OneDim.h:26
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition: OneDim.h:85
size_t size() const
Total solution vector length;.
Definition: OneDim.h:96
doublereal ssnorm(doublereal *x, doublereal *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x.
Definition: OneDim.cpp:290
size_t nDomains() const
Number of domains.
Definition: OneDim.h:54
void eval(size_t j, double *x, double *r, doublereal rdt=-1.0, int count=1)
Evaluate the multi-domain residual function.
Definition: OneDim.cpp:258
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition: OneDim.h:59
doublereal rdt() const
Reciprocal of the time step.
Definition: OneDim.h:150
void debuglog(const std::string &msg, int loglevel)
Write a message to the log only if loglevel > 0.
Definition: global.h:140
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:188
const double SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:149
void writelog(const std::string &fmt, const Args &... args)
Write a formatted message to the screen.
Definition: global.h:158
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
int solve(DenseMatrix &A, double *b, size_t nrhs, size_t ldb)
Solve Ax = b. Array b is overwritten on exit with x.
Various templated functions that carry out common vector operations (see Templated Utility Functions)...