Cantera  2.1.2
MultiNewton.cpp
Go to the documentation of this file.
1 /**
2  * @file MultiNewton.cpp Damped Newton solver for 1D multi-domain problems
3  */
4
5 /*
6  * Copyright 2001 California Institute of Technology
7  */
8
10
14
15 #include <cstdio>
16
17 using namespace std;
18
19 namespace Cantera
20 {
21
22 // unnamed-namespace for local helpers
23 namespace
24 {
25
26 class Indx
27 {
28 public:
29  Indx(size_t nv, size_t np) : m_nv(nv), m_np(np) {}
30  size_t m_nv, m_np;
31  size_t operator()(size_t m, size_t j) {
32  return j*m_nv + m;
33  }
34 };
35
36 /**
37  * Return a damping coefficient that keeps the solution after taking one
38  * Newton step between specified lower and upper bounds. This function only
39  * considers one domain.
40  */
41 doublereal bound_step(const doublereal* x, const doublereal* step,
42  Domain1D& r, int loglevel)
43 {
44  char buf;
45  size_t np = r.nPoints();
46  size_t nv = r.nComponents();
47  Indx index(nv, np);
48  doublereal above, below, val, newval;
49  size_t m, j;
50  doublereal fbound = 1.0;
51  bool wroteTitle = false;
52  for (m = 0; m < nv; m++) {
53  above = r.upperBound(m);
54  below = r.lowerBound(m);
55
56  for (j = 0; j < np; j++) {
57  val = x[index(m,j)];
58  if (loglevel > 0) {
59  if (val > above + 1.0e-12 || val < below - 1.0e-12) {
60  sprintf(buf, "domain %s: %20s(%s) = %10.3e (%10.3e, %10.3e)\n",
61  int2str(r.domainIndex()).c_str(),
62  r.componentName(m).c_str(), int2str(j).c_str(),
63  val, below, above);
64  writelog(string("\nERROR: solution out of bounds.\n")+buf);
65  }
66  }
67
68  newval = val + step[index(m,j)];
69
70  if (newval > above) {
71  fbound = std::max(0.0, std::min(fbound,
72  (above - val)/(newval - val)));
73  } else if (newval < below) {
74  fbound = std::min(fbound, (val - below)/(val - newval));
75  }
76
77  if (loglevel > 1 && (newval > above || newval < below)) {
78  if (!wroteTitle) {
79  writelog("\nNewton step takes solution out of bounds.\n\n");
80  sprintf(buf," %12s %12s %4s %10s %10s %10s %10s\n",
81  "domain","component","pt","value","step","min","max");
82  wroteTitle = true;
83  writelog(buf);
84  }
85  sprintf(buf, " %4s %12s %4s %10.3e %10.3e %10.3e %10.3e\n",
86  int2str(r.domainIndex()).c_str(),
87  r.componentName(m).c_str(), int2str(j).c_str(),
88  val, step[index(m,j)], below, above);
89  writelog(buf);
90  }
91  }
92  }
93  return fbound;
94 }
95
96 /**
97  * This function computes the square of a weighted norm of a step
98  * vector for one domain.
99  *
100  * @param x Solution vector for this domain.
101  * @param step Newton step vector for this domain.
102  * @param r Object representing the domain. Used to get tolerances,
103  * number of components, and number of points.
104  *
105  * The return value is
106  * \f[
107  * \sum_{n,j} \left(\frac{s_{n,j}}{w_n}\right)^2
108  * \f]
109  * where the error weight for solution component \f$n\f$ is given by
110  * \f[
111  * w_n = \epsilon_{r,n} \frac{\sum_j |x_{n,j}|}{J} + \epsilon_{a,n}.
112  * \f]
113  * Here \f$\epsilon_{r,n} \f$ is the relative error tolerance for
114  * component n, and multiplies the average magnitude of
115  * solution component n in the domain. The second term,
116  * \f$\epsilon_{a,n}\f$, is the absolute error tolerance for component
117  * n.
118  */
119 doublereal norm_square(const doublereal* x,
120  const doublereal* step, Domain1D& r)
121 {
122  doublereal f, ewt, esum, sum = 0.0;
123  size_t n, j;
124  doublereal f2max = 0.0;
125  size_t nv = r.nComponents();
126  size_t np = r.nPoints();
127
128  for (n = 0; n < nv; n++) {
129  esum = 0.0;
130  for (j = 0; j < np; j++) {
131  esum += fabs(x[nv*j + n]);
132  }
133  ewt = r.rtol(n)*esum/np + r.atol(n);
134  for (j = 0; j < np; j++) {
135  f = step[nv*j + n]/ewt;
136  sum += f*f;
137  if (f*f > f2max) {
138  f2max = f*f;
139  }
140  }
141  }
142  return sum;
143 }
144
145 } // end unnamed-namespace
146
147 //-----------------------------------------------------------
148 // constants
149 //-----------------------------------------------------------
150
151 const string dashedline =
152  "-----------------------------------------------------------------";
153
154 const doublereal DampFactor = sqrt(2.0);
155 const size_t NDAMP = 7;
156
157 //-----------------------------------------------------------
158 // MultiNewton methods
159 //-----------------------------------------------------------
160
161 MultiNewton::MultiNewton(int sz)
162  : m_maxAge(5)
163 {
164  m_n = sz;
165  m_elapsed = 0.0;
166 }
167
168 MultiNewton::~MultiNewton()
169 {
170  for (size_t i = 0; i < m_workarrays.size(); i++) {
171  delete[] m_workarrays[i];
172  }
173 }
174
175 void MultiNewton::resize(size_t sz)
176 {
177  m_n = sz;
178  for (size_t i = 0; i < m_workarrays.size(); i++) {
179  delete[] m_workarrays[i];
180  }
181  m_workarrays.clear();
182 }
183
184 doublereal MultiNewton::norm2(const doublereal* x,
185  const doublereal* step, OneDim& r) const
186 {
187  doublereal f, sum = 0.0;//, fmx = 0.0;
188  size_t nd = r.nDomains();
189  for (size_t n = 0; n < nd; n++) {
190  f = norm_square(x + r.start(n), step + r.start(n),
191  r.domain(n));
192  sum += f;
193  }
194  sum /= r.size();
195  return sqrt(sum);
196 }
197
198 void MultiNewton::step(doublereal* x, doublereal* step,
199  OneDim& r, MultiJac& jac, int loglevel)
200 {
201  size_t iok;
202  size_t sz = r.size();
203  r.eval(npos, x, step);
204 #undef DEBUG_STEP
205 #ifdef DEBUG_STEP
206  vector_fp ssave(sz, 0.0);
207  for (size_t n = 0; n < sz; n++) {
208  step[n] = -step[n];
209  ssave[n] = step[n];
210  }
211 #else
212  for (size_t n = 0; n < sz; n++) {
213  step[n] = -step[n];
214  }
215 #endif
216
217  iok = jac.solve(step, step);
218
219  // if iok is non-zero, then solve failed
220  if (iok != 0) {
221  iok--;
222  size_t nd = r.nDomains();
223  size_t n;
224  for (n = nd-1; n != npos; n--)
225  if (iok >= r.start(n)) {
226  break;
227  }
228  Domain1D& dom = r.domain(n);
229  size_t offset = iok - r.start(n);
230  size_t pt = offset/dom.nComponents();
231  size_t comp = offset - pt*dom.nComponents();
232  throw CanteraError("MultiNewton::step",
233  "Jacobian is singular for domain "+
234  dom.id() + ", component "
235  +dom.componentName(comp)+" at point "
236  +int2str(pt)+"\n(Matrix row "
237  +int2str(iok)+") \nsee file bandmatrix.csv\n");
238  } else if (int(iok) < 0)
239  throw CanteraError("MultiNewton::step",
240  "iok = "+int2str(iok));
241
242 #ifdef DEBUG_STEP
243  bool ok = false;
244  Domain1D* d;
245  if (!ok) {
246  for (size_t n = 0; n < sz; n++) {
247  d = r.pointDomain(n);
248  int nvd = d->nComponents();
249  int pt = (n - d->loc())/nvd;
250  cout << "step: " << pt << " " <<
251  r.pointDomain(n)->componentName(n - d->loc() - nvd*pt)
252  << " " << x[n] << " " << ssave[n] << " " << step[n] << endl;
253  }
254  }
255 #endif
256 }
257
258 doublereal MultiNewton::boundStep(const doublereal* x0,
259  const doublereal* step0, const OneDim& r, int loglevel)
260 {
261  doublereal fbound = 1.0;
262  for (size_t i = 0; i < r.nDomains(); i++) {
263  fbound = std::min(fbound,
264  bound_step(x0 + r.start(i), step0 + r.start(i),
265  r.domain(i), loglevel));
266  }
267  return fbound;
268 }
269
270 int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
271  doublereal* x1, doublereal* step1, doublereal& s1,
272  OneDim& r, MultiJac& jac, int loglevel, bool writetitle)
273 {
275  if (loglevel > 0 && writetitle) {
276  writelog("\n\nDamped Newton iteration:\n");
277  writelog(dashedline);
278
279  sprintf(m_buf,"\n%s %9s %9s %9s %9s %9s %5s %5s\n",
280  "m","F_damp","F_bound","log10(ss)",
281  "log10(s0)","log10(s1)","N_jac","Age");
282  writelog(m_buf);
283  writelog(dashedline+"\n");
284  }
285
286  // compute the weighted norm of the undamped step size step0
287  doublereal s0 = norm2(x0, step0, r);
288
289  // compute the multiplier to keep all components in bounds
290  doublereal fbound = boundStep(x0, step0, r, loglevel-1);
291
292  // if fbound is very small, then x0 is already close to the
293  // boundary and step0 points out of the allowed domain. In
294  // this case, the Newton algorithm fails, so return an error
295  // condition.
296  if (fbound < 1.e-10) {
297  writelog("\nAt limits.\n", loglevel);
298  return -3;
299  }
300
301
302  //--------------------------------------------
303  // Attempt damped step
304  //--------------------------------------------
305
306  // damping coefficient starts at 1.0
307  doublereal damp = 1.0;
308
309  doublereal ff;
310
311  size_t m;
312  for (m = 0; m < NDAMP; m++) {
313
314  ff = fbound*damp;
315
316  // step the solution by the damped step size
317  for (size_t j = 0; j < m_n; j++) {
318  x1[j] = ff*step0[j] + x0[j];
319  }
320
321  // compute the next undamped step that would result if x1
322  // is accepted
323  step(x1, step1, r, jac, loglevel-1);
324
325  // compute the weighted norm of step1
326  s1 = norm2(x1, step1, r);
327
329  if (loglevel > 0) {
330  doublereal ss = r.ssnorm(x1,step1);
331  sprintf(m_buf,"\n%s %9.5f %9.5f %9.5f %9.5f %9.5f %4d %d/%d",
332  int2str(m).c_str(), damp, fbound, log10(ss+SmallNumber),
333  log10(s0+SmallNumber),
334  log10(s1+SmallNumber),
335  jac.nEvals(), jac.age(), m_maxAge);
336  writelog(m_buf);
337  }
338
339  // if the norm of s1 is less than the norm of s0, then
340  // accept this damping coefficient. Also accept it if this
341  // step would result in a converged solution. Otherwise,
342  // decrease the damping coefficient and try again.
343
344  if (s1 < 1.0 || s1 < s0) {
345  break;
346  }
347  damp /= DampFactor;
348  }
349
350  // If a damping coefficient was found, return 1 if the
351  // solution after stepping by the damped step would represent
352  // a converged solution, and return 0 otherwise. If no damping
353  // coefficient could be found, return -2.
354  if (m < NDAMP) {
355  if (s1 > 1.0) {
356  return 0;
357  } else {
358  return 1;
359  }
360  } else {
361  return -2;
362  }
363 }
364
365 int MultiNewton::solve(doublereal* x0, doublereal* x1,
366  OneDim& r, MultiJac& jac, int loglevel)
367 {
368  clock_t t0 = clock();
369  int m = 0;
370  bool forceNewJac = false;
371  doublereal s1=1.e30;
372
373  doublereal* x = getWorkArray();
374  doublereal* stp = getWorkArray();
375  doublereal* stp1 = getWorkArray();
376
377  copy(x0, x0 + m_n, x);
378
379  bool frst = true;
380  doublereal rdt = r.rdt();
381  int j0 = jac.nEvals();
382  int nJacReeval = 0;
383
384  while (1 > 0) {
385
386  // Check whether the Jacobian should be re-evaluated.
387  if (jac.age() > m_maxAge) {
388  writelog("\nMaximum Jacobian age reached ("+int2str(m_maxAge)+")\n", loglevel);
389  forceNewJac = true;
390  }
391
392  if (forceNewJac) {
393  r.eval(npos, x, stp, 0.0, 0);
394  jac.eval(x, stp, 0.0);
396  forceNewJac = false;
397  }
398
399  // compute the undamped Newton step
400  step(x, stp, r, jac, loglevel-1);
401
402  // increment the Jacobian age
403  jac.incrementAge();
404
405  // damp the Newton step
406  m = dampStep(x, stp, x1, stp1, s1, r, jac, loglevel-1, frst);
407  if (loglevel == 1 && m >= 0) {
408  if (frst) {
409  sprintf(m_buf,"\n\n %10s %10s %5s ",
410  "log10(ss)","log10(s1)","N_jac");
411  writelog(m_buf);
412  sprintf(m_buf,"\n ------------------------------------");
413  writelog(m_buf);
414  }
415  doublereal ss = r.ssnorm(x, stp);
416  sprintf(m_buf,"\n %10.4f %10.4f %d ",
417  log10(ss),log10(s1),jac.nEvals());
418  writelog(m_buf);
419  }
420  frst = false;
421
422  // Successful step, but not converged yet. Take the damped
423  // step, and try again.
424  if (m == 0) {
425  copy(x1, x1 + m_n, x);
426  }
427
428  // convergence
429  else if (m == 1) {
430  break;
431  }
432
433  // If dampStep fails, first try a new Jacobian if an old
434  // one was being used. If it was a new Jacobian, then
435  // return -1 to signify failure.
436  else if (m < 0) {
437  if (jac.age() > 1) {
438  forceNewJac = true;
439  if (nJacReeval > 3) {
440  break;
441  }
442  nJacReeval++;
443  writelog("\nRe-evaluating Jacobian, since no damping "
444  "coefficient\ncould be found with this Jacobian.\n",
445  loglevel);
446  } else {
447  break;
448  }
449  }
450  }
451
452  if (m < 0) {
453  copy(x, x + m_n, x1);
454  }
455  if (m > 0 && jac.nEvals() == j0) {
456  m = 100;
457  }
458  releaseWorkArray(x);
459  releaseWorkArray(stp);
460  releaseWorkArray(stp1);
461  m_elapsed += (clock() - t0)/(1.0*CLOCKS_PER_SEC);
462  return m;
463 }
464
466 {
467  doublereal* w = 0;
468
469  if (!m_workarrays.empty()) {
470  w = m_workarrays.back();
471  m_workarrays.pop_back();
472  } else {
473  w = new doublereal[m_n];
474  }
475  return w;
476 }
477
478 void MultiNewton::releaseWorkArray(doublereal* work)
479 {
480  m_workarrays.push_back(work);
481 }
482
483 } // end namespace Cantera
Container class for multiple-domain 1D problems.
Definition: OneDim.h:22
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:40
int nEvals() const
Number of Jacobian evaluations.
Definition: MultiJac.h:45
int age() const
Number of times 'incrementAge' has been called since the last evaluation.
Definition: MultiJac.h:51
doublereal ssnorm(doublereal *x, doublereal *r)
Steady-state max norm (infinity norm) of the residual evaluated using solution x. ...
Definition: OneDim.cpp:268
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:236
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:173
void step(doublereal *x, doublereal *step, OneDim &r, MultiJac &jac, int loglevel)
Compute the undamped Newton step.
virtual std::string componentName(size_t n) const
Name of the nth component. May be overloaded.
Definition: Domain1D.h:208
virtual size_t loc(size_t j=0) const
Location of the start of the local solution vector in the global solution vector,.
Definition: Domain1D.h:461
int solve(const doublereal *const b, doublereal *const x)
Solve the matrix problem Ax = b.
Definition: BandMatrix.cpp:267
void incrementAge()
Increment the Jacobian age.
Definition: MultiJac.h:56
const double DampFactor
Dampfactor is the factor by which the damping factor is reduced by when a reduction in step length is...
Definition: BEulerInt.cpp:1327
size_t m_n
Number of rows and columns of the matrix.
Definition: BandMatrix.h:350
doublereal rdt() const
Reciprocal of the time step.
Definition: OneDim.h:135
size_t nComponents() const
Number of components at each grid point.
Definition: Domain1D.h:164
Base class for one-dimensional domains.
Definition: Domain1D.h:37
int solve(doublereal *x0, doublereal *x1, OneDim &r, MultiJac &jac, int loglevel)
Find the solution to F(X) = 0 by damped Newton iteration.
size_t start(size_t i) const
The index of the start of domain i in the solution vector.
Definition: OneDim.h:78
Domain1D & domain(size_t i) const
Return a reference to domain i.
Definition: OneDim.h:54
doublereal norm2(const doublereal *x, const doublereal *step, OneDim &r) const
Compute the weighted 2-norm of step.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
const int NDAMP
Number of damping steps that are carried out before the solution is deemed a failure.
Definition: BEulerInt.cpp:1328
void resize(size_t points)
Change the problem size.
size_t nDomains() const
Number of domains.
Definition: OneDim.h:49
void eval(doublereal *x0, doublereal *resid0, double rdt)
Evaluate the Jacobian at x0.
Definition: MultiJac.cpp:48
Class MultiJac evaluates the Jacobian of a system of equations defined by a residual function supplie...
Definition: MultiJac.h:26
doublereal boundStep(const doublereal *x0, const doublereal *step0, const OneDim &r, int loglevel)
Return the factor by which the undamped Newton step 'step0' must be multiplied in order to keep all s...
const doublereal SmallNumber
smallest number to compare to zero.
Definition: ct_defs.h:139
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:165
Templates for operations on vector-like objects.
Contains declarations for string manipulation functions within Cantera.
int dampStep(const doublereal *x0, const doublereal *step0, doublereal *x1, doublereal *step1, doublereal &s1, OneDim &r, MultiJac &jac, int loglevel, bool writetitle)
On entry, step0 must contain an undamped Newton step for the solution x0.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
void releaseWorkArray(doublereal *work)
Release a work array by pushing its pointer onto the stack of available arrays.
void writelog(const std::string &msg)
Write a message to the screen.
Definition: global.cpp:43
size_t size() const
Total solution vector length;.
Definition: OneDim.h:83
Domain1D * pointDomain(size_t i)
Return a pointer to the domain global point i belongs to.
Definition: OneDim.cpp:224
doublereal * getWorkArray()
Get a pointer to an array of length m_n for temporary work space.
Definitions for the classes that are thrown when Cantera experiences an error condition (also contain...