Cantera  2.0
MultiNewton.cpp
Go to the documentation of this file.
1 /**
2  * @file MultiNewton.cpp
3  *
4  * Damped Newton solver for 1D multi-domain problems
5  */
6
7 /*
8  * Copyright 2001 California Institute of Technology
9  */
10
11 #include <vector>
12 using namespace std;
13
15
19
20 #include <cstdio>
21 #include <cmath>
22 #include <ctime>
23
24 using namespace std;
25
26 namespace Cantera
27 {
28
29 //----------------------------------------------------------
30 // function declarations
31 //----------------------------------------------------------
32
33 // declarations for functions in newton_utils.h
34 doublereal bound_step(const doublereal* x,
35  const doublereal* step, Domain1D& r, int loglevel=0);
36 doublereal norm_square(const doublereal* x,
37  const doublereal* step, Domain1D& r);
38
39
40
41 //-----------------------------------------------------------
42 // constants
43 //-----------------------------------------------------------
44
45 const string dashedline =
46  "-----------------------------------------------------------------";
47
48 const doublereal DampFactor = sqrt(2.0);
49 const size_t NDAMP = 7;
50
51
52
53 //-----------------------------------------------------------
54 // MultiNewton methods
55 //-----------------------------------------------------------
56
57
58 MultiNewton::MultiNewton(int sz)
59  : m_maxAge(5)
60 {
61  m_n = sz;
62  m_elapsed = 0.0;
63 }
64
65 MultiNewton::~MultiNewton()
66 {
67  for (size_t i = 0; i < m_workarrays.size(); i++) {
68  delete[] m_workarrays[i];
69  }
70 }
71
72 /**
73  * Prepare for a new solution vector length.
74  */
75 void MultiNewton::resize(size_t sz)
76 {
77  m_n = sz;
78  for (size_t i = 0; i < m_workarrays.size(); i++) {
79  delete[] m_workarrays[i];
80  }
81  m_workarrays.clear();
82 }
83
84
85 /**
86  * Compute the weighted 2-norm of 'step'.
87  */
88 doublereal MultiNewton::norm2(const doublereal* x,
89  const doublereal* step, OneDim& r) const
90 {
91  doublereal f, sum = 0.0;//, fmx = 0.0;
92  size_t nd = r.nDomains();
93  for (size_t n = 0; n < nd; n++) {
94  f = norm_square(x + r.start(n), step + r.start(n),
95  r.domain(n));
96  sum += f;
97  }
98  sum /= r.size();
99  return sqrt(sum);
100 }
101
102
103 /**
104  * Compute the undamped Newton step. The residual function is
105  * evaluated at x, but the Jacobian is not recomputed.
106  */
107 void MultiNewton::step(doublereal* x, doublereal* step,
108  OneDim& r, MultiJac& jac, int loglevel)
109 {
110  size_t iok;
111  size_t sz = r.size();
112  r.eval(npos, x, step);
113 #undef DEBUG_STEP
114 #ifdef DEBUG_STEP
115  vector_fp ssave(sz, 0.0);
116  for (size_t n = 0; n < sz; n++) {
117  step[n] = -step[n];
118  ssave[n] = step[n];
119  }
120 #else
121  for (size_t n = 0; n < sz; n++) {
122  step[n] = -step[n];
123  }
124 #endif
125
126  iok = jac.solve(step, step);
127
128  // if iok is non-zero, then solve failed
129  if (iok != 0) {
130  iok--;
131  size_t nd = r.nDomains();
132  size_t n;
133  for (n = nd-1; n != npos; n--)
134  if (iok >= r.start(n)) {
135  break;
136  }
137  Domain1D& dom = r.domain(n);
138  size_t offset = iok - r.start(n);
139  size_t pt = offset/dom.nComponents();
140  size_t comp = offset - pt*dom.nComponents();
141  throw CanteraError("MultiNewton::step",
142  "Jacobian is singular for domain "+
143  dom.id() + ", component "
144  +dom.componentName(comp)+" at point "
145  +int2str(pt)+"\n(Matrix row "
146  +int2str(iok)+") \nsee file bandmatrix.csv\n");
147  } else if (int(iok) < 0)
148  throw CanteraError("MultiNewton::step",
149  "iok = "+int2str(iok));
150
151 #ifdef DEBUG_STEP
152  bool ok = false;
153  Domain1D* d;
154  if (!ok) {
155  for (size_t n = 0; n < sz; n++) {
156  d = r.pointDomain(n);
157  int nvd = d->nComponents();
158  int pt = (n - d->loc())/nvd;
159  cout << "step: " << pt << " " <<
160  r.pointDomain(n)->componentName(n - d->loc() - nvd*pt)
161  << " " << x[n] << " " << ssave[n] << " " << step[n] << endl;
162  }
163  }
164 #endif
165 }
166
167
168 /**
169  * Return the factor by which the undamped Newton step 'step0'
170  * must be multiplied in order to keep all solution components in
171  * all domains between their specified lower and upper bounds.
172  */
173 doublereal MultiNewton::boundStep(const doublereal* x0,
174  const doublereal* step0, const OneDim& r, int loglevel)
175 {
176  doublereal fbound = 1.0;
177  for (size_t i = 0; i < r.nDomains(); i++) {
178  fbound = std::min(fbound,
179  bound_step(x0 + r.start(i), step0 + r.start(i),
180  r.domain(i), loglevel));
181  }
182  return fbound;
183 }
184
185
186 /**
187  * On entry, step0 must contain an undamped Newton step for the
188  * solution x0. This method attempts to find a damping coefficient
189  * such that the next undamped step would have a norm smaller than
190  * that of step0. If successful, the new solution after taking the
191  * damped step is returned in x1, and the undamped step at x1 is
192  * returned in step1.
193  */
194 int MultiNewton::dampStep(const doublereal* x0, const doublereal* step0,
195  doublereal* x1, doublereal* step1, doublereal& s1,
196  OneDim& r, MultiJac& jac, int loglevel, bool writetitle)
197 {
198
200  if (loglevel > 0 && writetitle) {
201  writelog("\n\nDamped Newton iteration:\n");
202  writelog(dashedline);
203
204  sprintf(m_buf,"\n%s %9s %9s %9s %9s %9s %5s %5s\n",
205  "m","F_damp","F_bound","log10(ss)",
206  "log10(s0)","log10(s1)","N_jac","Age");
207  writelog(m_buf);
208  writelog(dashedline+"\n");
209  }
210
211  // compute the weighted norm of the undamped step size step0
212  doublereal s0 = norm2(x0, step0, r);
213
214  // compute the multiplier to keep all components in bounds
215  doublereal fbound = boundStep(x0, step0, r, loglevel-1);
216
217  // if fbound is very small, then x0 is already close to the
218  // boundary and step0 points out of the allowed domain. In
219  // this case, the Newton algorithm fails, so return an error
220  // condition.
221  if (fbound < 1.e-10) {
222  if (loglevel > 0) {
223  writelog("\nAt limits.\n");
224  }
225  return -3;
226  }
227
228
229  //--------------------------------------------
230  // Attempt damped step
231  //--------------------------------------------
232
233  // damping coefficient starts at 1.0
234  doublereal damp = 1.0;
235
236  doublereal ff;
237
238  size_t m;
239  for (m = 0; m < NDAMP; m++) {
240
241  ff = fbound*damp;
242
243  // step the solution by the damped step size
244  for (size_t j = 0; j < m_n; j++) {
245  x1[j] = ff*step0[j] + x0[j];
246  }
247
248  // compute the next undamped step that would result if x1
249  // is accepted
250  step(x1, step1, r, jac, loglevel-1);
251
252  // compute the weighted norm of step1
253  s1 = norm2(x1, step1, r);
254
256  if (loglevel > 0) {
257  doublereal ss = r.ssnorm(x1,step1);
258  sprintf(m_buf,"\n%s %9.5f %9.5f %9.5f %9.5f %9.5f %4d %d/%d",
259  int2str(m).c_str(), damp, fbound, log10(ss+SmallNumber),
260  log10(s0+SmallNumber),
261  log10(s1+SmallNumber),
262  jac.nEvals(), jac.age(), m_maxAge);
263  writelog(m_buf);
264  }
265
266  // if the norm of s1 is less than the norm of s0, then
267  // accept this damping coefficient. Also accept it if this
268  // step would result in a converged solution. Otherwise,
269  // decrease the damping coefficient and try again.
270
271  if (s1 < 1.0 || s1 < s0) {
272  break;
273  }
274  damp /= DampFactor;
275  }
276
277  // If a damping coefficient was found, return 1 if the
278  // solution after stepping by the damped step would represent
279  // a converged solution, and return 0 otherwise. If no damping
280  // coefficient could be found, return -2.
281  if (m < NDAMP) {
282  if (s1 > 1.0) {
283  return 0;
284  } else {
285  return 1;
286  }
287  } else {
288  return -2;
289  }
290 }
291
292
293 /**
294  * Find the solution to F(X) = 0 by damped Newton iteration. On
295  * entry, x0 contains an initial estimate of the solution. On
296  * successful return, x1 contains the converged solution.
297  */
298 int MultiNewton::solve(doublereal* x0, doublereal* x1,
299  OneDim& r, MultiJac& jac, int loglevel)
300 {
301  clock_t t0 = clock();
302  int m = 0;
303  bool forceNewJac = false;
304  doublereal s1=1.e30;
305
306  doublereal* x = getWorkArray();
307  doublereal* stp = getWorkArray();
308  doublereal* stp1 = getWorkArray();
309
310  copy(x0, x0 + m_n, x);
311
312  bool frst = true;
313  doublereal rdt = r.rdt();
314  int j0 = jac.nEvals();
315  int nJacReeval = 0;
316
317  while (1 > 0) {
318
319  // Check whether the Jacobian should be re-evaluated.
320  if (jac.age() > m_maxAge) {
321  if (loglevel > 0) {
322  writelog("\nMaximum Jacobian age reached ("+int2str(m_maxAge)+")\n");
323  }
324  forceNewJac = true;
325  }
326
327  if (forceNewJac) {
328  r.eval(npos, x, stp, 0.0, 0);
329  jac.eval(x, stp, 0.0);
331  forceNewJac = false;
332  }
333
334  // compute the undamped Newton step
335  step(x, stp, r, jac, loglevel-1);
336
337  // increment the Jacobian age
338  jac.incrementAge();
339
340  // damp the Newton step
341  m = dampStep(x, stp, x1, stp1, s1, r, jac, loglevel-1, frst);
342  if (loglevel == 1 && m >= 0) {
343  if (frst) {
344  sprintf(m_buf,"\n\n %10s %10s %5s ",
345  "log10(ss)","log10(s1)","N_jac");
346  writelog(m_buf);
347  sprintf(m_buf,"\n ------------------------------------");
348  writelog(m_buf);
349  }
350  doublereal ss = r.ssnorm(x, stp);
351  sprintf(m_buf,"\n %10.4f %10.4f %d ",
352  log10(ss),log10(s1),jac.nEvals());
353  writelog(m_buf);
354  }
355  frst = false;
356
357  // Successful step, but not converged yet. Take the damped
358  // step, and try again.
359  if (m == 0) {
360  copy(x1, x1 + m_n, x);
361  }
362
363  // convergence
364  else if (m == 1) {
365  goto done;
366  }
367
368  // If dampStep fails, first try a new Jacobian if an old
369  // one was being used. If it was a new Jacobian, then
370  // return -1 to signify failure.
371  else if (m < 0) {
372  if (jac.age() > 1) {
373  forceNewJac = true;
374  if (nJacReeval > 3) {
375  goto done;
376  }
377  nJacReeval++;
378  if (loglevel > 0)
379  writelog("\nRe-evaluating Jacobian, since no damping "
380  "coefficient\ncould be found with this Jacobian.\n");
381  } else {
382  goto done;
383  }
384  }
385  }
386
387 done:
388  if (m < 0) {
389  copy(x, x + m_n, x1);
390  }
391  if (m > 0 && jac.nEvals() == j0) {
392  m = 100;
393  }
394  releaseWorkArray(x);
395  releaseWorkArray(stp);
396  releaseWorkArray(stp1);
397  m_elapsed += (clock() - t0)/(1.0*CLOCKS_PER_SEC);
398  return m;
399 }
400
401
402 /**
403  * Get a pointer to an array of length m_n for temporary work
404  * space.
405  */
407 {
408  doublereal* w = 0;
409
410  if (!m_workarrays.empty()) {
411  w = m_workarrays.back();
412  m_workarrays.pop_back();
413  } else {
414  w = new doublereal[m_n];
415  }
416  return w;
417 }
418
419 /**
420  * Release a work array by pushing its pointer onto the stack of
421  * available arrays.
422  */
423 void MultiNewton::releaseWorkArray(doublereal* work)
424 {
425  m_workarrays.push_back(work);
426 }
427 }
428
429
430 // \$Log: Newton.cpp,v