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