Cantera  2.2.1
NonlinearSolver.cpp
Go to the documentation of this file.
1 /**
2  * @file NonlinearSolver.cpp Damped Newton solver for 0D and 1D problems
3  */
4
5 /*
6  * Copyright 2004 Sandia Corporation. Under the terms of Contract
7  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
8  * retains certain rights in this software.
9  * See file License.txt for licensing information.
10  */
11
14
15 #include "cantera/base/clockWC.h"
17
18 #include <limits>
19
20 using namespace std;
21
22 namespace Cantera
23 {
24
25 //-----------------------------------------------------------
26 // Constants
27 //-----------------------------------------------------------
28 //! Dampfactor is the factor by which the damping factor is reduced by when a
29 //! reduction in step length is warranted
30 const doublereal DampFactor = 4.0;
31 //! Number of damping steps that are carried out before the solution is deemed
32 //! a failure
33 const int NDAMP = 7;
34
35 bool NonlinearSolver::s_TurnOffTiming(false);
36
37 #ifdef DEBUG_NUMJAC
38 bool NonlinearSolver::s_print_NumJac(true);
39 #else
40 bool NonlinearSolver::s_print_NumJac(false);
41 #endif
42
43 // Turn off printing of dogleg information
44 bool NonlinearSolver::s_print_DogLeg(false);
45
46 // Turn off solving the system twice and comparing the answer.
47 /*
48  * Turn this on if you want to compare the Hessian and Newton solve results.
49  */
50 bool NonlinearSolver::s_doBothSolvesAndCompare(false);
51
52 // This toggle turns off the use of the Hessian when it is warranted by the condition number.
53 /*
54  * This is a debugging option.
55  */
56 bool NonlinearSolver::s_alwaysAssumeNewtonGood(false);
57
58 NonlinearSolver::NonlinearSolver(ResidJacEval* func) :
59  m_func(func),
61  neq_(0),
62  m_manualDeltaStepSet(0),
63  m_normResid_0(0.0),
64  m_normResid_Bound(0.0),
65  m_normResid_1(0.0),
66  m_normDeltaSoln_Newton(0.0),
67  m_normDeltaSoln_CP(0.0),
68  m_normResidTrial(0.0),
69  m_resid_scaled(false),
70  m_dampBound(1.0),
71  m_dampRes(1.0),
72  delta_t_n(-1.0),
73  m_nfe(0),
74  m_colScaling(0),
75  m_rowScaling(0),
76  m_numTotalLinearSolves(0),
77  m_numTotalNewtIts(0),
78  m_min_newt_its(0),
79  maxNewtIts_(100),
80  m_jacFormMethod(NSOLN_JAC_NUM),
81  m_nJacEval(0),
82  time_n(0.0),
83  m_matrixConditioning(0),
84  m_order(1),
85  rtol_(1.0E-3),
86  atolBase_(1.0E-10),
87  userResidRtol_(1.0E-3),
88  checkUserResidualTols_(0),
89  m_print_flag(0),
90  m_ScaleSolnNormToResNorm(0.001),
91  jacCopyPtr_(0),
92  HessianPtr_(0),
93  residNorm2Cauchy_(0.0),
94  dogLegID_(0),
95  dogLegAlpha_(1.0),
96  RJd_norm_(0.0),
97  lambdaStar_(0.0),
98  Jd_(0),
99  deltaX_trust_(0),
100  norm_deltaX_trust_(0.0),
101  trustDelta_(1.0),
102  trustRegionInitializationMethod_(2),
103  trustRegionInitializationFactor_(1.0),
104  Nuu_(0.0),
105  dist_R0_(0.0),
106  dist_R1_(0.0),
107  dist_R2_(0.0),
108  dist_Total_(0.0),
109  JdJd_norm_(0.0),
110  normTrust_Newton_(0.0),
111  normTrust_CP_(0.0),
112  doDogLeg_(0),
113  doAffineSolve_(0) ,
114  CurrentTrustFactor_(1.0),
115  NextTrustFactor_(1.0),
116  ResidWtsReevaluated_(false),
117  ResidDecreaseSDExp_(0.0),
118  ResidDecreaseSD_(0.0),
119  ResidDecreaseNewtExp_(0.0),
120  ResidDecreaseNewt_(0.0)
121 {
122  warn_deprecated("class NonlinearSolver", "To be removed after Cantera 2.2.");
123  neq_ = m_func->nEquations();
124
125  m_ewt.resize(neq_, rtol_);
126  m_deltaStepMinimum.resize(neq_, 0.001);
127  m_deltaStepMaximum.resize(neq_, 1.0E10);
128  m_y_n_curr.resize(neq_, 0.0);
129  m_ydot_n_curr.resize(neq_, 0.0);
130  m_y_nm1.resize(neq_, 0.0);
131  m_ydot_nm1.resize(neq_, 0.0);
132  m_y_n_trial.resize(neq_, 0.0);
133  m_ydot_trial.resize(neq_, 0.0);
134  m_colScales.resize(neq_, 1.0);
135  m_rowScales.resize(neq_, 1.0);
136  m_rowWtScales.resize(neq_, 1.0);
137  m_resid.resize(neq_, 0.0);
138  m_wksp.resize(neq_, 0.0);
139  m_wksp_2.resize(neq_, 0.0);
140  m_residWts.resize(neq_, 0.0);
141  atolk_.resize(neq_, atolBase_);
142  deltaX_Newton_.resize(neq_, 0.0);
143  m_step_1.resize(neq_, 0.0);
144  m_y_n_trial.resize(neq_, 0.0);
145  doublereal hb = std::numeric_limits<double>::max();
146  m_y_high_bounds.resize(neq_, hb);
147  m_y_low_bounds.resize(neq_, -hb);
148
149  for (size_t i = 0; i < neq_; i++) {
150  atolk_[i] = atolBase_;
151  m_ewt[i] = atolk_[i];
152  }
153
154  deltaX_CP_.resize(neq_, 0.0);
155  Jd_.resize(neq_, 0.0);
156  deltaX_trust_.resize(neq_, 1.0);
157 }
158
160  m_func(right.m_func),
162  neq_(0),
163  m_manualDeltaStepSet(0),
164  m_normResid_0(0.0),
165  m_normResid_Bound(0.0),
166  m_normResid_1(0.0),
167  m_normDeltaSoln_Newton(0.0),
168  m_normDeltaSoln_CP(0.0),
169  m_normResidTrial(0.0),
170  m_resid_scaled(false),
171  m_dampBound(1.0),
172  m_dampRes(1.0),
173  delta_t_n(-1.0),
174  m_nfe(0),
175  m_colScaling(0),
176  m_rowScaling(0),
177  m_numTotalLinearSolves(0),
178  m_numTotalNewtIts(0),
179  m_min_newt_its(0),
180  maxNewtIts_(100),
181  m_jacFormMethod(NSOLN_JAC_NUM),
182  m_nJacEval(0),
183  time_n(0.0),
184  m_matrixConditioning(0),
185  m_order(1),
186  rtol_(1.0E-3),
187  atolBase_(1.0E-10),
188  userResidRtol_(1.0E-3),
189  checkUserResidualTols_(0),
190  m_print_flag(0),
191  m_ScaleSolnNormToResNorm(0.001),
192  jacCopyPtr_(0),
193  HessianPtr_(0),
194  residNorm2Cauchy_(0.0),
195  dogLegID_(0),
196  dogLegAlpha_(1.0),
197  RJd_norm_(0.0),
198  lambdaStar_(0.0),
199  norm_deltaX_trust_(0.0),
200  trustDelta_(1.0),
201  trustRegionInitializationMethod_(2),
202  trustRegionInitializationFactor_(1.0),
203  Nuu_(0.0),
204  dist_R0_(0.0),
205  dist_R1_(0.0),
206  dist_R2_(0.0),
207  dist_Total_(0.0),
208  JdJd_norm_(0.0),
209  normTrust_Newton_(0.0),
210  normTrust_CP_(0.0),
211  doDogLeg_(0),
212  doAffineSolve_(0),
213  CurrentTrustFactor_(1.0),
214  NextTrustFactor_(1.0),
215  ResidWtsReevaluated_(false),
216  ResidDecreaseSDExp_(0.0),
217  ResidDecreaseSD_(0.0),
218  ResidDecreaseNewtExp_(0.0),
219  ResidDecreaseNewt_(0.0)
220 {
221  *this =operator=(right);
222 }
223
225 {
226  delete jacCopyPtr_;
227  delete HessianPtr_;
228 }
229
231 {
232  if (this == &right) {
233  return *this;
234  }
235  // rely on the ResidJacEval duplMyselfAsresidJacEval() function to
236  // create a deep copy
238
239  solnType_ = right.solnType_;
240  neq_ = right.neq_;
241  m_ewt = right.m_ewt;
244  m_y_n_curr = right.m_y_n_curr;
246  m_y_nm1 = right.m_y_nm1;
247  m_ydot_nm1 = right.m_ydot_nm1;
248  m_y_n_trial = right.m_y_n_trial;
249  m_ydot_trial = right.m_ydot_trial;
250  m_step_1 = right.m_step_1;
251  m_colScales = right.m_colScales;
252  m_rowScales = right.m_rowScales;
254  m_resid = right.m_resid;
255  m_wksp = right.m_wksp;
256  m_wksp_2 = right.m_wksp_2;
257  m_residWts = right.m_residWts;
267  m_dampBound = right.m_dampBound;
268  m_dampRes = right.m_dampRes;
269  delta_t_n = right.delta_t_n;
270  m_nfe = right.m_nfe;
271  m_colScaling = right.m_colScaling;
272  m_rowScaling = right.m_rowScaling;
276  maxNewtIts_ = right.maxNewtIts_;
278  m_nJacEval = right.m_nJacEval;
279  time_n = right.time_n;
281  m_order = right.m_order;
282  rtol_ = right.rtol_;
283  atolBase_ = right.atolBase_;
284  atolk_ = right.atolk_;
288  m_print_flag = right.m_print_flag;
290
291  delete jacCopyPtr_;
292  jacCopyPtr_ = (right.jacCopyPtr_)->duplMyselfAsGeneralMatrix();
293  delete HessianPtr_;
294  HessianPtr_ = (right.HessianPtr_)->duplMyselfAsGeneralMatrix();
295
296  deltaX_CP_ = right.deltaX_CP_;
299  dogLegID_ = right.dogLegID_;
300  dogLegAlpha_ = right.dogLegAlpha_;
301  RJd_norm_ = right.RJd_norm_;
302  lambdaStar_ = right.lambdaStar_;
303  Jd_ = right.Jd_;
306  trustDelta_ = right.trustDelta_;
309  Nuu_ = right.Nuu_;
310  dist_R0_ = right.dist_R0_;
311  dist_R1_ = right.dist_R1_;
312  dist_R2_ = right.dist_R2_;
313  dist_Total_ = right.dist_Total_;
314  JdJd_norm_ = right.JdJd_norm_;
317  doDogLeg_ = right.doDogLeg_;
321
327
328  return *this;
329 }
330
331 void NonlinearSolver::createSolnWeights(const doublereal* const y)
332 {
333  for (size_t i = 0; i < neq_; i++) {
334  m_ewt[i] = rtol_ * fabs(y[i]) + atolk_[i];
335  if (DEBUG_MODE_ENABLED && m_ewt[i] <= 0.0) {
336  throw CanteraError(" NonlinearSolver::createSolnWeights()", "ewts <= 0.0");
337  }
338  }
339 }
340
341 void NonlinearSolver::setBoundsConstraints(const doublereal* const y_low_bounds,
342  const doublereal* const y_high_bounds)
343 {
344  for (size_t i = 0; i < neq_; i++) {
345  m_y_low_bounds[i] = y_low_bounds[i];
346  m_y_high_bounds[i] = y_high_bounds[i];
347  }
348 }
349
350 void NonlinearSolver::setSolverScheme(int doDogLeg, int doAffineSolve)
351 {
352  doDogLeg_ = doDogLeg;
353  doAffineSolve_ = doAffineSolve;
354 }
355
357 {
358  return m_y_low_bounds;
359 }
360
362 {
363  return m_y_high_bounds;
364 }
365
366 doublereal NonlinearSolver::solnErrorNorm(const doublereal* const delta_y, const char* title, int printLargest,
367  const doublereal dampFactor) const
368 {
369  doublereal sum_norm = 0.0, error;
370  for (size_t i = 0; i < neq_; i++) {
371  error = delta_y[i] / m_ewt[i];
372  sum_norm += (error * error);
373  }
374  sum_norm = sqrt(sum_norm / neq_);
375  if (printLargest) {
376  if ((printLargest == 1) || (m_print_flag >= 4 && m_print_flag <= 5)) {
377
378  printf("\t\t solnErrorNorm(): ");
379  if (title) {
380  printf("%s", title);
381  } else {
382  printf(" Delta soln norm ");
383  }
384  printf(" = %-11.4E\n", sum_norm);
385  } else if (m_print_flag >= 6) {
386
387
388
389  const int num_entries = printLargest;
390  printf("\t\t ");
391  writeline('-', 90);
392  printf("\t\t solnErrorNorm(): ");
393  if (title) {
394  printf("%s", title);
395  } else {
396  printf(" Delta soln norm ");
397  }
398  printf(" = %-11.4E\n", sum_norm);
399
400  doublereal dmax1, normContrib;
401  int j;
402  std::vector<size_t> imax(num_entries, npos);
403  printf("\t\t Printout of Largest Contributors: (damp = %g)\n", dampFactor);
404  printf("\t\t I weightdeltaY/sqtN| deltaY "
405  "ysolnOld ysolnNew Soln_Weights\n");
406  printf("\t\t ");
407  writeline('-', 88);
408
409  for (int jnum = 0; jnum < num_entries; jnum++) {
410  dmax1 = -1.0;
411  for (size_t i = 0; i < neq_; i++) {
412  bool used = false;
413  for (j = 0; j < jnum; j++) {
414  if (imax[j] == i) {
415  used = true;
416  }
417  }
418  if (!used) {
419  error = delta_y[i] / m_ewt[i];
420  normContrib = sqrt(error * error);
421  if (normContrib > dmax1) {
422  imax[jnum] = i;
423  dmax1 = normContrib;
424  }
425  }
426  }
427  size_t i = imax[jnum];
428  if (i != npos) {
429  error = delta_y[i] / m_ewt[i];
430  normContrib = sqrt(error * error);
431  printf("\t\t %4s %12.4e | %12.4e %12.4e %12.4e %12.4e\n", int2str(i).c_str(), normContrib/sqrt((double)neq_),
432  delta_y[i], m_y_n_curr[i], m_y_n_curr[i] + dampFactor * delta_y[i], m_ewt[i]);
433
434  }
435  }
436  printf("\t\t ");
437  writeline('-', 90);
438  }
439  }
440  return sum_norm;
441 }
442
443 doublereal NonlinearSolver::residErrorNorm(const doublereal* const resid, const char* title, const int printLargest,
444  const doublereal* const y) const
445 {
446  doublereal sum_norm = 0.0, error;
447
448  for (size_t i = 0; i < neq_; i++) {
449  if (DEBUG_MODE_ENABLED) {
450  checkFinite(resid[i]);
451  }
452  error = resid[i] / m_residWts[i];
453  if (DEBUG_MODE_ENABLED) {
455  }
456  sum_norm += (error * error);
457  }
458  sum_norm = sqrt(sum_norm / neq_);
459  if (DEBUG_MODE_ENABLED) {
460  checkFinite(sum_norm);
461  }
462  if (printLargest) {
463  const int num_entries = printLargest;
464  doublereal dmax1, normContrib;
465  int j;
466  std::vector<size_t> imax(num_entries, npos);
467
468  if (m_print_flag >= 4 && m_print_flag <= 5) {
469  printf("\t\t residErrorNorm():");
470  if (title) {
471  printf(" %s ", title);
472  } else {
473  printf(" residual L2 norm ");
474  }
475  printf("= %12.4E\n", sum_norm);
476  }
477  if (m_print_flag >= 6) {
478  printf("\t\t ");
479  writeline('-', 90);
480  printf("\t\t residErrorNorm(): ");
481  if (title) {
482  printf(" %s ", title);
483  } else {
484  printf(" residual L2 norm ");
485  }
486  printf("= %12.4E\n", sum_norm);
487  printf("\t\t Printout of Largest Contributors to norm:\n");
488  printf("\t\t I |Resid/ResWt| UnsclRes ResWt | y_curr\n");
489  printf("\t\t ");
490  writeline('-', 88);
491  for (int jnum = 0; jnum < num_entries; jnum++) {
492  dmax1 = -1.0;
493  for (size_t i = 0; i < neq_; i++) {
494  bool used = false;
495  for (j = 0; j < jnum; j++) {
496  if (imax[j] == i) {
497  used = true;
498  }
499  }
500  if (!used) {
501  error = resid[i] / m_residWts[i];
502  normContrib = sqrt(error * error);
503  if (normContrib > dmax1) {
504  imax[jnum] = i;
505  dmax1 = normContrib;
506  }
507  }
508  }
509  size_t i = imax[jnum];
510  if (i != npos) {
511  error = resid[i] / m_residWts[i];
512  normContrib = sqrt(error * error);
513  printf("\t\t %4s %12.4e %12.4e %12.4e | %12.4e\n", int2str(i).c_str(), normContrib, resid[i], m_residWts[i], y[i]);
514  }
515  }
516
517  printf("\t\t ");
518  writeline('-', 90);
519  }
520  }
521  return sum_norm;
522 }
523
524 void NonlinearSolver::setColumnScaling(bool useColScaling, const double* const scaleFactors)
525 {
526  if (useColScaling) {
527  if (scaleFactors) {
528  m_colScaling = 2;
529  for (size_t i = 0; i < neq_; i++) {
530  m_colScales[i] = scaleFactors[i];
531  if (m_colScales[i] <= 1.0E-200) {
532  throw CanteraError("NonlinearSolver::setColumnScaling() ERROR", "Bad column scale factor");
533  }
534  }
535  } else {
536  m_colScaling = 1;
537  }
538  } else {
539  m_colScaling = 0;
540  }
541 }
542
543 void NonlinearSolver::setRowScaling(bool useRowScaling)
544 {
545  m_rowScaling = useRowScaling;
546 }
547
549 {
550  if (m_colScaling == 1) {
551  for (size_t i = 0; i < neq_; i++) {
552  m_colScales[i] = m_ewt[i];
553  }
554  } else {
555  for (size_t i = 0; i < neq_; i++) {
556  m_colScales[i] = 1.0;
557  }
558  }
559  if (m_colScaling) {
561  }
562 }
563
564 int NonlinearSolver::doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal* const y_curr,
565  const doublereal* const ydot_curr, const ResidEval_Type_Enum evalType) const
566 {
567  int retn = m_func->evalResidNJ(time_curr, delta_t_n, y_curr, ydot_curr, DATA_PTR(m_resid), evalType);
568  m_nfe++;
569  m_resid_scaled = false;
570  return retn;
571 }
572
573 void NonlinearSolver::scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, doublereal* const ydot_comm,
574  doublereal time_curr, int num_newt_its)
575 {
576  int irow, jcol;
577  size_t ivec[2];
578  jac.nRowsAndStruct(ivec);
579  double* colP_j;
580
581  /*
582  * Column scaling -> We scale the columns of the Jacobian
583  * by the nominal important change in the solution vector
584  */
585  if (m_colScaling) {
586  if (!jac.factored()) {
587  if (jac.matrixType_ == 0) {
588  /*
589  * Go get new scales -> Took this out of this inner loop.
590  * Needs to be done at a larger scale.
591  */
592  // setColumnScales();
593
594  /*
595  * Scale the new Jacobian
596  */
597  doublereal* jptr = &(*(jac.begin()));
598  for (jcol = 0; jcol < (int) neq_; jcol++) {
599  for (irow = 0; irow < (int) neq_; irow++) {
600  *jptr *= m_colScales[jcol];
601  jptr++;
602  }
603  }
604  } else if (jac.matrixType_ == 1) {
605  int kl = static_cast<int>(ivec[0]);
606  int ku = static_cast<int>(ivec[1]);
607  for (jcol = 0; jcol < (int) neq_; jcol++) {
608  colP_j = (doublereal*) jac.ptrColumn(jcol);
609  for (irow = jcol - ku; irow <= jcol + kl; irow++) {
610  if (irow >= 0 && irow < (int) neq_) {
611  colP_j[kl + ku + irow - jcol] *= m_colScales[jcol];
612  }
613  }
614  }
615  }
616  }
617  }
618  /*
619  * row sum scaling -> Note, this is an unequivocal success
620  * at keeping the small numbers well balanced and nonnegative.
621  */
622  if (! jac.factored()) {
623  /*
624  * Ok, this is ugly. jac.begin() returns an vector<double> iterator
625  * to the first data location.
626  * Then &(*()) reverts it to a doublereal *.
627  */
628  doublereal* jptr = &(*(jac.begin()));
629  for (irow = 0; irow < (int) neq_; irow++) {
630  m_rowScales[irow] = 0.0;
631  m_rowWtScales[irow] = 0.0;
632  }
633  if (jac.matrixType_ == 0) {
634  for (jcol = 0; jcol < (int) neq_; jcol++) {
635  for (irow = 0; irow < (int) neq_; irow++) {
636  if (m_rowScaling) {
637  m_rowScales[irow] += fabs(*jptr);
638  }
639  if (m_colScaling) {
640  // This is needed in order to mitgate the change in J_ij carried out just above this loop.
641  // Alternatively, we could move this loop up to the top
642  m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol] / m_colScales[jcol];
643  } else {
644  m_rowWtScales[irow] += fabs(*jptr) * m_ewt[jcol];
645  }
646  if (DEBUG_MODE_ENABLED) {
647  checkFinite(m_rowWtScales[irow]);
648  }
649  jptr++;
650  }
651  }
652  } else if (jac.matrixType_ == 1) {
653  int kl = static_cast<int>(ivec[0]);
654  int ku = static_cast<int>(ivec[1]);
655  for (jcol = 0; jcol < (int) neq_; jcol++) {
656  colP_j = (doublereal*) jac.ptrColumn(jcol);
657  for (irow = jcol - ku; irow <= jcol + kl; irow++) {
658  if (irow >= 0 && irow < (int) neq_) {
659  double vv = fabs(colP_j[kl + ku + irow - jcol]);
660  if (m_rowScaling) {
661  m_rowScales[irow] += vv;
662  }
663  if (m_colScaling) {
664  // This is needed in order to mitgate the change in J_ij carried out just above this loop.
665  // Alternatively, we could move this loop up to the top
666  m_rowWtScales[irow] += vv * m_ewt[jcol] / m_colScales[jcol];
667  } else {
668  m_rowWtScales[irow] += vv * m_ewt[jcol];
669  }
670  if (DEBUG_MODE_ENABLED) {
671  checkFinite(m_rowWtScales[irow]);
672  }
673  }
674  }
675  }
676  }
677  if (m_rowScaling) {
678  for (irow = 0; irow < (int) neq_; irow++) {
679  m_rowScales[irow] = 1.0/m_rowScales[irow];
680  }
681  } else {
682  for (irow = 0; irow < (int) neq_; irow++) {
683  m_rowScales[irow] = 1.0;
684  }
685  }
686  // What we have defined is a maximum value that the residual can be and still pass.
687  // This isn't sufficient.
688
689  if (m_rowScaling) {
690  if (jac.matrixType_ == 0) {
691  jptr = &(*(jac.begin()));
692  for (jcol = 0; jcol < (int) neq_; jcol++) {
693  for (irow = 0; irow < (int) neq_; irow++) {
694  *jptr *= m_rowScales[irow];
695  jptr++;
696  }
697  }
698  } else if (jac.matrixType_ == 1) {
699  int kl = static_cast<int>(ivec[0]);
700  int ku = static_cast<int>(ivec[1]);
701  for (jcol = 0; jcol < (int) neq_; jcol++) {
702  colP_j = (doublereal*) jac.ptrColumn(jcol);
703  for (irow = jcol - ku; irow <= jcol + kl; irow++) {
704  if (irow >= 0 && irow < (int) neq_) {
705  colP_j[kl + ku + irow - jcol] *= m_rowScales[irow];
706  }
707  }
708  }
709  }
710  }
711
712  if (num_newt_its % 5 == 1) {
713  computeResidWts();
714  }
715
716  }
717 }
718
720 {
721  if (! jacCopyPtr_->factored()) {
722
723  if (checkUserResidualTols_ != 1) {
724  doublereal sum = 0.0;
725  for (size_t irow = 0; irow < neq_; irow++) {
726  m_residWts[irow] = m_rowWtScales[irow] / neq_;
727  sum += m_residWts[irow];
728  }
729  sum /= neq_;
730  for (size_t irow = 0; irow < neq_; irow++) {
731  m_residWts[irow] = (m_residWts[irow] + atolBase_ * atolBase_ * sum);
732  }
733  if (checkUserResidualTols_ == 2) {
734  for (size_t irow = 0; irow < neq_; irow++) {
735  m_residWts[irow] = std::min(m_residWts[irow], userResidAtol_[irow] + userResidRtol_ * m_rowWtScales[irow] / neq_);
736  }
737  }
738  } else {
739  for (size_t irow = 0; irow < neq_; irow++) {
740  m_residWts[irow] = userResidAtol_[irow] + userResidRtol_ * m_rowWtScales[irow] / neq_;
741  }
742  }
743
744
745  for (size_t irow = 0; irow < neq_; irow++) {
746  m_wksp[irow] = 0.0;
747  }
748  doublereal* jptr = &((*jacCopyPtr_)(0,0));
749  for (size_t jcol = 0; jcol < neq_; jcol++) {
750  for (size_t irow = 0; irow < neq_; irow++) {
751  m_wksp[irow] += (*jptr) * m_ewt[jcol];
752  jptr++;
753  }
754  }
755  doublereal resNormOld = 0.0;
756  doublereal error;
757
758  for (size_t irow = 0; irow < neq_; irow++) {
759  error = m_wksp[irow] / m_residWts[irow];
760  resNormOld += error * error;
761  }
762  resNormOld = sqrt(resNormOld / neq_);
763
764  if (resNormOld > 0.0) {
765  m_ScaleSolnNormToResNorm = resNormOld;
766  }
768
769  // Recalculate the residual weights now that we know the value of m_ScaleSolnNormToResNorm
770  computeResidWts();
771  } else {
772  throw CanteraError("NonlinearSolver::calcSolnToResNormVector()" , "Logic error");
773  }
774 }
775
776 int NonlinearSolver::doNewtonSolve(const doublereal time_curr, const doublereal* const y_curr,
777  const doublereal* const ydot_curr, doublereal* const delta_y,
778  GeneralMatrix& jac)
779 {
780  // multiply the residual by -1
781  if (m_rowScaling && !m_resid_scaled) {
782  for (size_t n = 0; n < neq_; n++) {
783  delta_y[n] = -m_rowScales[n] * m_resid[n];
784  }
785  m_resid_scaled = true;
786  } else {
787  for (size_t n = 0; n < neq_; n++) {
788  delta_y[n] = -m_resid[n];
789  }
790  }
791
792
793
794  /*
795  * Solve the system -> This also involves inverting the
796  * matrix
797  */
798  int info = jac.solve(DATA_PTR(delta_y));
799
800
801  /*
802  * reverse the column scaling if there was any.
803  */
804  if (m_colScaling) {
805  for (size_t irow = 0; irow < neq_; irow++) {
806  delta_y[irow] = delta_y[irow] * m_colScales[irow];
807  }
808  }
809
810 #ifdef DEBUG_JAC
811  if (printJacContributions) {
812  for (size_t iNum = 0; iNum < numRows; iNum++) {
813  if (iNum > 0) {
814  focusRow++;
815  }
816  doublereal dsum = 0.0;
817  vector_fp& Jdata = jacBack.data();
818  doublereal dRow = Jdata[neq_ * focusRow + focusRow];
819  printf("\n Details on delta_Y for row %d \n", focusRow);
820  printf(" Value before = %15.5e, delta = %15.5e,"
821  "value after = %15.5e\n", y_curr[focusRow],
822  delta_y[focusRow], y_curr[focusRow] + delta_y[focusRow]);
823  if (!freshJac) {
824  printf(" Old Jacobian\n");
825  }
826  printf(" col delta_y aij "
827  "contrib \n");
828  printf("--------------------------------------------------"
829  "---------------------------------------------\n");
830  printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n",
831  focusRow, delta_y[focusRow],
832  dRow, RRow[iNum] / dRow, RRow[iNum]);
833  dsum += RRow[iNum] / dRow;
834  for (size_t ii = 0; ii < neq_; ii++) {
835  if (ii != focusRow) {
836  doublereal aij = Jdata[neq_ * ii + focusRow];
837  doublereal contrib = aij * delta_y[ii] * (-1.0) / dRow;
838  dsum += contrib;
839  if (fabs(contrib) > Pcutoff) {
840  printf("%6d %15.5e %15.5e %15.5e\n", ii,
841  delta_y[ii] , aij, contrib);
842  }
843  }
844  }
845  printf("--------------------------------------------------"
846  "---------------------------------------------\n");
847  printf(" %15.5e %15.5e\n",
848  delta_y[focusRow], dsum);
849  }
850  }
851
852 #endif
853
856  return info;
857 }
858
859 int NonlinearSolver::doAffineNewtonSolve(const doublereal* const y_curr, const doublereal* const ydot_curr,
860  doublereal* const delta_y, GeneralMatrix& jac)
861 {
862  /*
863  * Notes on banded Hessian solve: The matrix for jT j has a larger band
864  * width. Both the top and bottom band widths are doubled, going from KU
865  * to KU+KL and KL to KU+KL in size. This is not an impossible increase
866  * in cost, but has to be considered.
867  */
868  bool newtonGood = true;
869  // We can default to QR here ( or not )
870  jac.useFactorAlgorithm(1);
871  int useQR = jac.factorAlgorithm();
872  // multiply the residual by -1
873  // Scale the residual if there is row scaling. Note, the matrix has already been scaled
874  if (m_rowScaling && !m_resid_scaled) {
875  for (size_t n = 0; n < neq_; n++) {
876  delta_y[n] = -m_rowScales[n] * m_resid[n];
877  }
878  m_resid_scaled = true;
879  } else {
880  for (size_t n = 0; n < neq_; n++) {
881  delta_y[n] = -m_resid[n];
882  }
883  }
884
885  // Factor the matrix using a standard Newton solve
886  m_conditionNumber = 1.0E300;
887  int info = 0;
888  if (!jac.factored()) {
889  if (useQR) {
890  info = jac.factorQR();
891  } else {
892  info = jac.factor();
893  }
894  }
895  /*
896  * Find the condition number of the matrix
897  * If we have failed to factor, we will fall back to calculating and factoring a modified Hessian
898  */
899  if (info == 0) {
900  doublereal rcond = 0.0;
901  if (useQR) {
902  rcond = jac.rcondQR();
903  } else {
904  doublereal a1norm = jac.oneNorm();
905  rcond = jac.rcond(a1norm);
906  }
907  if (rcond > 0.0) {
908  m_conditionNumber = 1.0 / rcond;
909  }
910  } else {
911  m_conditionNumber = 1.0E300;
912  newtonGood = false;
913  if (m_print_flag >= 1) {
914  printf("\t\t doAffineNewtonSolve: ");
915  if (useQR) {
916  printf("factorQR()");
917  } else {
918  printf("factor()");
919  }
920  printf(" returned with info = %d, indicating a zero row or column\n", info);
921  }
922  }
923  bool doHessian = false;
925  doHessian = true;
926  }
927  if (m_conditionNumber < 1.0E7) {
928  if (m_print_flag >= 4) {
929  printf("\t\t doAffineNewtonSolve: Condition number = %g during regular solve\n", m_conditionNumber);
930  }
931
932  /*
933  * Solve the system -> This also involves inverting the matrix
934  */
935  int info = jac.solve(DATA_PTR(delta_y));
936  if (info) {
937  if (m_print_flag >= 2) {
938  printf("\t\t doAffineNewtonSolve() ERROR: QRSolve returned INFO = %d. Switching to Hessian solve\n", info);
939  }
940  doHessian = true;
941  newtonGood = false;
942  }
943  /*
944  * reverse the column scaling if there was any on a successful solve
945  */
946  if (m_colScaling) {
947  for (size_t irow = 0; irow < neq_; irow++) {
948  delta_y[irow] = delta_y[irow] * m_colScales[irow];
949  }
950  }
951
952  } else {
953  if (jac.matrixType_ == 1) {
954  newtonGood = true;
955  if (m_print_flag >= 3) {
956  printf("\t\t doAffineNewtonSolve() WARNING: Condition number too large, %g, But Banded Hessian solve "
957  "not implemented yet \n", m_conditionNumber);
958  }
959  } else {
960  doHessian = true;
961  newtonGood = false;
962  if (m_print_flag >= 3) {
963  printf("\t\t doAffineNewtonSolve() WARNING: Condition number too large, %g. Doing a Hessian solve \n", m_conditionNumber);
964  }
965  }
966  }
967
968  if (doHessian) {
969  // Store the old value for later comparison
970  vector_fp delyNewton(neq_);
971  for (size_t irow = 0; irow < neq_; irow++) {
972  delyNewton[irow] = delta_y[irow];
973  }
974
975  // Get memory if not done before
976  if (HessianPtr_ == 0) {
978  }
979
980  /*
981  * Calculate the symmetric Hessian
982  */
983  GeneralMatrix& hessian = *HessianPtr_;
984  GeneralMatrix& jacCopy = *jacCopyPtr_;
985  hessian.zero();
986  if (m_rowScaling) {
987  for (size_t i = 0; i < neq_; i++) {
988  for (size_t j = i; j < neq_; j++) {
989  for (size_t k = 0; k < neq_; k++) {
990  hessian(i,j) += jacCopy(k,i) * jacCopy(k,j) * m_rowScales[k] * m_rowScales[k];
991  }
992  hessian(j,i) = hessian(i,j);
993  }
994  }
995  } else {
996  for (size_t i = 0; i < neq_; i++) {
997  for (size_t j = i; j < neq_; j++) {
998  for (size_t k = 0; k < neq_; k++) {
999  hessian(i,j) += jacCopy(k,i) * jacCopy(k,j);
1000  }
1001  hessian(j,i) = hessian(i,j);
1002  }
1003  }
1004  }
1005
1006  /*
1007  * Calculate the matrix norm of the Hessian
1008  */
1009  doublereal hnorm = 0.0;
1010  doublereal hcol = 0.0;
1011  if (m_colScaling) {
1012  for (size_t i = 0; i < neq_; i++) {
1013  for (size_t j = i; j < neq_; j++) {
1014  hcol += fabs(hessian(j,i)) * m_colScales[j];
1015  }
1016  for (size_t j = i+1; j < neq_; j++) {
1017  hcol += fabs(hessian(i,j)) * m_colScales[j];
1018  }
1019  hcol *= m_colScales[i];
1020  if (hcol > hnorm) {
1021  hnorm = hcol;
1022  }
1023  }
1024  } else {
1025  for (size_t i = 0; i < neq_; i++) {
1026  for (size_t j = i; j < neq_; j++) {
1027  hcol += fabs(hessian(j,i));
1028  }
1029  for (size_t j = i+1; j < neq_; j++) {
1030  hcol += fabs(hessian(i,j));
1031  }
1032  if (hcol > hnorm) {
1033  hnorm = hcol;
1034  }
1035  }
1036  }
1037  /*
1038  * Add junk to the Hessian diagonal
1039  * -> Note, testing indicates that this will get too big for ill-conditioned systems.
1040  */
1041  hcol = sqrt(static_cast<double>(neq_)) * 1.0E-7 * hnorm;
1042 #ifdef DEBUG_HKM_NOT
1043  if (hcol > 1.0) {
1044  hcol = 1.0E1;
1045  }
1046 #endif
1047  if (m_colScaling) {
1048  for (size_t i = 0; i < neq_; i++) {
1049  hessian(i,i) += hcol / (m_colScales[i] * m_colScales[i]);
1050  }
1051  } else {
1052  for (size_t i = 0; i < neq_; i++) {
1053  hessian(i,i) += hcol;
1054  }
1055  }
1056
1057  /*
1058  * Factor the Hessian
1059  */
1060  int info = 0;
1061  ct_dpotrf(ctlapack::UpperTriangular, neq_, &(*(HessianPtr_->begin())), neq_, info);
1062  if (info) {
1063  if (m_print_flag >= 2) {
1064  printf("\t\t doAffineNewtonSolve() ERROR: Hessian isn't positive definite DPOTRF returned INFO = %d\n", info);
1065  }
1066  return info;
1067  }
1068
1069  vector_fp delyH(neq_);
1070  // First recalculate the scaled residual. It got wiped out doing the Newton solve
1071  if (m_rowScaling) {
1072  for (size_t n = 0; n < neq_; n++) {
1073  delyH[n] = -m_rowScales[n] * m_resid[n];
1074  }
1075  } else {
1076  for (size_t n = 0; n < neq_; n++) {
1077  delyH[n] = -m_resid[n];
1078  }
1079  }
1080
1081  if (m_rowScaling) {
1082  for (size_t j = 0; j < neq_; j++) {
1083  delta_y[j] = 0.0;
1084  for (size_t i = 0; i < neq_; i++) {
1085  delta_y[j] += delyH[i] * jacCopy(i,j) * m_rowScales[i];
1086  }
1087  }
1088  } else {
1089  for (size_t j = 0; j < neq_; j++) {
1090  delta_y[j] = 0.0;
1091  for (size_t i = 0; i < neq_; i++) {
1092  delta_y[j] += delyH[i] * jacCopy(i,j);
1093  }
1094  }
1095  }
1096
1097
1098  /*
1099  * Solve the factored Hessian System
1100  */
1101  ct_dpotrs(ctlapack::UpperTriangular, neq_, 1,&(*(hessian.begin())), neq_, delta_y, neq_, info);
1102  if (info) {
1103  if (m_print_flag >= 2) {
1104  printf("\t\t NonlinearSolver::doAffineNewtonSolve() ERROR: DPOTRS returned INFO = %d\n", info);
1105  }
1106  return info;
1107  }
1108  /*
1109  * reverse the column scaling if there was any.
1110  */
1111  if (m_colScaling) {
1112  for (size_t irow = 0; irow < neq_; irow++) {
1113  delta_y[irow] = delta_y[irow] * m_colScales[irow];
1114  }
1115  }
1116
1117
1118  if (doDogLeg_ && m_print_flag > 7) {
1119  double normNewt = solnErrorNorm(&delyNewton[0]);
1120  double normHess = solnErrorNorm(delta_y);
1121  printf("\t\t doAffineNewtonSolve(): Printout Comparison between Hessian deltaX and Newton deltaX\n");
1122
1123  printf("\t\t I Hessian+Junk Newton");
1124  if (newtonGood || s_alwaysAssumeNewtonGood) {
1125  printf(" (USING NEWTON DIRECTION)\n");
1126  } else {
1127  printf(" (USING HESSIAN DIRECTION)\n");
1128  }
1129  printf("\t\t Norm: %12.4E %12.4E\n", normHess, normNewt);
1130
1131  printf("\t\t --------------------------------------------------------\n");
1132  for (size_t i = 0; i < neq_; i++) {
1133  printf("\t\t %3s %13.5E %13.5E\n", int2str(i).c_str(), delta_y[i], delyNewton[i]);
1134  }
1135  printf("\t\t --------------------------------------------------------\n");
1136  } else if (doDogLeg_ && m_print_flag >= 4) {
1137  double normNewt = solnErrorNorm(&delyNewton[0]);
1138  double normHess = solnErrorNorm(delta_y);
1139  printf("\t\t doAffineNewtonSolve(): Hessian update norm = %12.4E \n"
1140  "\t\t Newton update norm = %12.4E \n", normHess, normNewt);
1141  if (newtonGood || s_alwaysAssumeNewtonGood) {
1142  printf("\t\t (USING NEWTON DIRECTION)\n");
1143  } else {
1144  printf("\t\t (USING HESSIAN DIRECTION)\n");
1145  }
1146  }
1147
1148  /*
1149  * Choose the delta_y to use
1150  */
1151  if (newtonGood || s_alwaysAssumeNewtonGood) {
1152  copy(delyNewton.begin(), delyNewton.end(), delta_y);
1153  }
1154  }
1155
1156 #ifdef DEBUG_JAC
1157  if (printJacContributions) {
1158  for (int iNum = 0; iNum < numRows; iNum++) {
1159  if (iNum > 0) {
1160  focusRow++;
1161  }
1162  doublereal dsum = 0.0;
1163  vector_fp& Jdata = jacBack.data();
1164  doublereal dRow = Jdata[neq_ * focusRow + focusRow];
1165  printf("\n Details on delta_Y for row %d \n", focusRow);
1166  printf(" Value before = %15.5e, delta = %15.5e,"
1167  "value after = %15.5e\n", y_curr[focusRow],
1168  delta_y[focusRow], y_curr[focusRow] + delta_y[focusRow]);
1169  if (!freshJac) {
1170  printf(" Old Jacobian\n");
1171  }
1172  printf(" col delta_y aij "
1173  "contrib \n");
1174  printf("-----------------------------------------------------------------------------------------------\n");
1175  printf(" Res(%d) %15.5e %15.5e %15.5e (Res = %g)\n",
1176  focusRow, delta_y[focusRow],
1177  dRow, RRow[iNum] / dRow, RRow[iNum]);
1178  dsum += RRow[iNum] / dRow;
1179  for (int ii = 0; ii < neq_; ii++) {
1180  if (ii != focusRow) {
1181  doublereal aij = Jdata[neq_ * ii + focusRow];
1182  doublereal contrib = aij * delta_y[ii] * (-1.0) / dRow;
1183  dsum += contrib;
1184  if (fabs(contrib) > Pcutoff) {
1185  printf("%6d %15.5e %15.5e %15.5e\n", ii,
1186  delta_y[ii] , aij, contrib);
1187  }
1188  }
1189  }
1190  printf("-----------------------------------------------------------------------------------------------\n");
1191  printf(" %15.5e %15.5e\n",
1192  delta_y[focusRow], dsum);
1193  }
1194  }
1195
1196 #endif
1197
1200  return info;
1201
1202 }
1203
1205 {
1206  doublereal rowFac = 1.0;
1207  doublereal colFac = 1.0;
1208  doublereal normSoln;
1209  // Calculate the descent direction
1210  /*
1211  * For confirmation of the scaling factors, see Dennis and Schnabel p, 152, p, 156 and my notes
1212  *
1213  * The colFac and rowFac values are used to eliminate the scaling of the matrix from the
1214  * actual equation
1215  *
1216  * Here we calculate the steepest descent direction. This is equation (11) in the notes. It is
1217  * stored in deltaX_CP_[].The value corresponds to d_descent[].
1218  */
1219  for (size_t j = 0; j < neq_; j++) {
1220  deltaX_CP_[j] = 0.0;
1221  if (m_colScaling) {
1222  colFac = 1.0 / m_colScales[j];
1223  }
1224  for (size_t i = 0; i < neq_; i++) {
1225  if (m_rowScaling) {
1226  rowFac = 1.0 / m_rowScales[i];
1227  }
1228  deltaX_CP_[j] -= m_resid[i] * jac(i,j) * colFac * rowFac * m_ewt[j] * m_ewt[j]
1229  / (m_residWts[i] * m_residWts[i]);
1230  if (DEBUG_MODE_ENABLED) {
1231  checkFinite(deltaX_CP_[j]);
1232  }
1233  }
1234  }
1235
1236  /*
1237  * Calculate J_hat d_y_descent. This is formula 18 in the notes.
1238  */
1239  for (size_t i = 0; i < neq_; i++) {
1240  Jd_[i] = 0.0;
1241  if (m_rowScaling) {
1242  rowFac = 1.0 / m_rowScales[i];
1243  } else {
1244  rowFac = 1.0;
1245  }
1246  for (size_t j = 0; j < neq_; j++) {
1247  if (m_colScaling) {
1248  colFac = 1.0 / m_colScales[j];
1249  }
1250  Jd_[i] += deltaX_CP_[j] * jac(i,j) * rowFac * colFac / m_residWts[i];
1251  }
1252  }
1253
1254  /*
1255  * Calculate the distance along the steepest descent until the Cauchy point
1256  * This is Eqn. 17 in the notes.
1257  */
1258  RJd_norm_ = 0.0;
1259  JdJd_norm_ = 0.0;
1260  for (size_t i = 0; i < neq_; i++) {
1261  RJd_norm_ += m_resid[i] * Jd_[i] / m_residWts[i];
1262  JdJd_norm_ += Jd_[i] * Jd_[i];
1263  }
1264  if (fabs(JdJd_norm_) < 1.0E-290) {
1265  if (fabs(RJd_norm_) < 1.0E-300) {
1266  lambdaStar_ = 0.0;
1267  } else {
1268  throw CanteraError("NonlinearSolver::doCauchyPointSolve()", "Unexpected condition: norms are zero");
1269  }
1270  } else {
1272  }
1273
1274  /*
1275  * Now we modify the steepest descent vector such that its length is equal to the
1276  * Cauchy distance. From now on, if we want to recreate the descent vector, we have
1277  * to unnormalize it by dividing by lambdaStar_.
1278  */
1279  for (size_t i = 0; i < neq_; i++) {
1280  deltaX_CP_[i] *= lambdaStar_;
1281  }
1282
1283
1284  doublereal normResid02 = m_normResid_0 * m_normResid_0 * neq_;
1285
1286  /*
1287  * Calculate the expected square of the risdual at the Cauchy point if the linear model is correct
1288  */
1289  if (fabs(JdJd_norm_) < 1.0E-290) {
1290  residNorm2Cauchy_ = normResid02;
1291  } else {
1292  residNorm2Cauchy_ = normResid02 - RJd_norm_ * RJd_norm_ / (JdJd_norm_);
1293  }
1294
1295
1296  // Extra printout section
1297  if (m_print_flag > 2) {
1298  // Calculate the expected residual at the Cauchy point if the linear model is correct
1299  doublereal residCauchy = 0.0;
1300  if (residNorm2Cauchy_ > 0.0) {
1301  residCauchy = sqrt(residNorm2Cauchy_ / neq_);
1302  } else {
1303  if (fabs(JdJd_norm_) < 1.0E-290) {
1304  residCauchy = m_normResid_0;
1305  } else {
1306  residCauchy = m_normResid_0 - sqrt(RJd_norm_ * RJd_norm_ / (JdJd_norm_));
1307  }
1308  }
1309  // Compute the weighted norm of the undamped step size descentDir_[]
1310  if ((s_print_DogLeg || doDogLeg_) && m_print_flag >= 6) {
1311  normSoln = solnErrorNorm(DATA_PTR(deltaX_CP_), "SteepestDescentDir", 10);
1312  } else {
1313  normSoln = solnErrorNorm(DATA_PTR(deltaX_CP_), "SteepestDescentDir", 0);
1314  }
1315  if ((s_print_DogLeg || doDogLeg_) && m_print_flag >= 5) {
1316  printf("\t\t doCauchyPointSolve: Steepest descent to Cauchy point: \n");
1317  printf("\t\t\t R0 = %g \n", m_normResid_0);
1318  printf("\t\t\t Rpred = %g\n", residCauchy);
1319  printf("\t\t\t Rjd = %g\n", RJd_norm_);
1320  printf("\t\t\t JdJd = %g\n", JdJd_norm_);
1321  printf("\t\t\t deltaX = %g\n", normSoln);
1322  printf("\t\t\t lambda = %g\n", lambdaStar_);
1323  }
1324  } else {
1325  // Calculate the norm of the Cauchy solution update in any case
1326  normSoln = solnErrorNorm(DATA_PTR(deltaX_CP_), "SteepestDescentDir", 0);
1327  }
1328  return normSoln;
1329 }
1330
1331 void NonlinearSolver::descentComparison(doublereal time_curr, doublereal* ydot0, doublereal* ydot1, int& numTrials)
1332 {
1333  doublereal ff = 1.0E-5;
1334  doublereal ffNewt = 1.0E-5;
1335  doublereal* y_n_1 = DATA_PTR(m_wksp);
1336  doublereal cauchyDistanceNorm = solnErrorNorm(DATA_PTR(deltaX_CP_));
1337  if (cauchyDistanceNorm < 1.0E-2) {
1338  ff = 1.0E-9 / cauchyDistanceNorm;
1339  if (ff > 1.0E-2) {
1340  ff = 1.0E-2;
1341  }
1342  }
1343  for (size_t i = 0; i < neq_; i++) {
1344  y_n_1[i] = m_y_n_curr[i] + ff * deltaX_CP_[i];
1345  }
1346  /*
1347  * Calculate the residual that would result if y1[] were the new solution vector
1348  * -> m_resid[] contains the result of the residual calculation
1349  */
1351  doResidualCalc(time_curr, solnType_, y_n_1, ydot1, Base_LaggedSolutionComponents);
1352  } else {
1353  doResidualCalc(time_curr, solnType_, y_n_1, ydot0, Base_LaggedSolutionComponents);
1354  }
1355
1356  doublereal normResid02 = m_normResid_0 * m_normResid_0 * neq_;
1357  doublereal residSteep = residErrorNorm(DATA_PTR(m_resid));
1358  doublereal residSteep2 = residSteep * residSteep * neq_;
1359  doublereal funcDecreaseSD = 0.5 * (residSteep2 - normResid02) / (ff * cauchyDistanceNorm);
1360
1361  doublereal sNewt = solnErrorNorm(DATA_PTR(deltaX_Newton_));
1362  if (sNewt > 1.0) {
1363  ffNewt = ffNewt / sNewt;
1364  }
1365  for (size_t i = 0; i < neq_; i++) {
1366  y_n_1[i] = m_y_n_curr[i] + ffNewt * deltaX_Newton_[i];
1367  }
1368  /*
1369  * Calculate the residual that would result if y1[] were the new solution vector.
1370  * Here we use the lagged solution components in the residual calculation as well. We are
1371  * interested in the linear model and its agreement with the nonlinear model.
1372  *
1373  * -> m_resid[] contains the result of the residual calculation
1374  */
1376  doResidualCalc(time_curr, solnType_, y_n_1, ydot1, Base_LaggedSolutionComponents);
1377  } else {
1378  doResidualCalc(time_curr, solnType_, y_n_1, ydot0, Base_LaggedSolutionComponents);
1379  }
1380  doublereal residNewt = residErrorNorm(DATA_PTR(m_resid));
1381  doublereal residNewt2 = residNewt * residNewt * neq_;
1382
1383  doublereal funcDecreaseNewt2 = 0.5 * (residNewt2 - normResid02) / (ffNewt * sNewt);
1384
1385  // This is the expected initial rate of decrease in the Cauchy direction.
1386  // -> This is Eqn. 29 = Rhat dot Jhat dy / || d ||
1387  doublereal funcDecreaseSDExp = RJd_norm_ / cauchyDistanceNorm * lambdaStar_;
1388
1389  doublereal funcDecreaseNewtExp2 = - normResid02 / sNewt;
1390
1391  if (m_normResid_0 > 1.0E-100) {
1392  ResidDecreaseSDExp_ = funcDecreaseSDExp / neq_ / m_normResid_0;
1393  ResidDecreaseSD_ = funcDecreaseSD / neq_ / m_normResid_0;
1394  ResidDecreaseNewtExp_ = funcDecreaseNewtExp2 / neq_ / m_normResid_0;
1395  ResidDecreaseNewt_ = funcDecreaseNewt2 / neq_ / m_normResid_0;
1396  } else {
1397  ResidDecreaseSDExp_ = 0.0;
1398  ResidDecreaseSD_ = funcDecreaseSD / neq_;
1399  ResidDecreaseNewtExp_ = 0.0;
1400  ResidDecreaseNewt_ = funcDecreaseNewt2 / neq_;
1401  }
1402  numTrials += 2;
1403
1404  /*
1405  * HKM These have been shown to exactly match up.
1406  * The steepest direction is always largest even when there are variable solution weights
1407  *
1408  * HKM When a Hessian is used with junk on the diagonal, funcDecreaseNewtExp2 is no longer accurate as the
1409  * direction gets significantly shorter with increasing condition number. This suggests an algorithm where the
1410  * Newton step from the Hessian should be increased so as to match funcDecreaseNewtExp2 = funcDecreaseNewt2.
1411  * This roughly equals the ratio of the norms of the Hessian and Newton steps. This increased Newton step can
1412  * then be used with the trust region double dogleg algorithm.
1413  */
1414  if ((s_print_DogLeg && m_print_flag >= 3) || (doDogLeg_ && m_print_flag >= 5)) {
1415  printf("\t\t descentComparison: initial rate of decrease of func in Cauchy dir (expected) = %g\n", funcDecreaseSDExp);
1416  printf("\t\t descentComparison: initial rate of decrease of func in Cauchy dir = %g\n", funcDecreaseSD);
1417  printf("\t\t descentComparison: initial rate of decrease of func in Newton dir (expected) = %g\n", funcDecreaseNewtExp2);
1418  printf("\t\t descentComparison: initial rate of decrease of func in Newton dir = %g\n", funcDecreaseNewt2);
1419  }
1420  if ((s_print_DogLeg && m_print_flag >= 3) || (doDogLeg_ && m_print_flag >= 4)) {
1421  printf("\t\t descentComparison: initial rate of decrease of Resid in Cauchy dir (expected) = %g\n", ResidDecreaseSDExp_);
1422  printf("\t\t descentComparison: initial rate of decrease of Resid in Cauchy dir = %g\n", ResidDecreaseSD_);
1423  printf("\t\t descentComparison: initial rate of decrease of Resid in Newton dir (expected) = %g\n", ResidDecreaseNewtExp_);
1424  printf("\t\t descentComparison: initial rate of decrease of Resid in Newton dir = %g\n", ResidDecreaseNewt_);
1425  }
1426
1427  if ((s_print_DogLeg && m_print_flag >= 5) || (doDogLeg_ && m_print_flag >= 5)) {
1428  if (funcDecreaseNewt2 >= 0.0) {
1429  printf("\t\t %13.5E %22.16E\n", funcDecreaseNewtExp2, m_normResid_0);
1430  double ff = ffNewt * 1.0E-5;
1431  for (int ii = 0; ii < 13; ii++) {
1432  ff *= 10.;
1433  if (ii == 12) {
1434  ff = ffNewt;
1435  }
1436  for (size_t i = 0; i < neq_; i++) {
1437  y_n_1[i] = m_y_n_curr[i] + ff * deltaX_Newton_[i];
1438  }
1439  numTrials += 1;
1441  doResidualCalc(time_curr, solnType_, y_n_1, ydot1, Base_LaggedSolutionComponents);
1442  } else {
1443  doResidualCalc(time_curr, solnType_, y_n_1, ydot0, Base_LaggedSolutionComponents);
1444  }
1445  residNewt = residErrorNorm(DATA_PTR(m_resid));
1446  residNewt2 = residNewt * residNewt * neq_;
1447  funcDecreaseNewt2 = 0.5 * (residNewt2 - normResid02) / (ff * sNewt);
1448  printf("\t\t %10.3E %13.5E %22.16E\n", ff, funcDecreaseNewt2, residNewt);
1449  }
1450  }
1451  }
1452 }
1453
1455 {
1456  /*
1457  * Gamma = ||grad f ||**4
1458  * ---------------------------------------------
1460  */
1461  /*
1462  * This hasn't worked. so will do it heuristically. One issue is that the Newton
1463  * direction is not the inverse of the Hessian times the gradient. The Hessian
1464  * is the matrix squared. Until I have the inverse of the Hessian from QR factorization
1465  * I may not be able to do it this way.
1466  */
1467
1468  /*
1469  * Heuristic algorithm - Find out where on the Newton line the residual is the same
1470  * as the residual at the Cauchy point. Then, go halfway to
1471  * the Newton point and call that Nuu.
1472  * Maybe we need to check that the linearized residual is
1473  * monotonic along that line. However, we haven't needed to yet.
1474  */
1475  doublereal residSteepLin = expectedResidLeg(0, 1.0);
1476  doublereal Nres2CP = residSteepLin * residSteepLin * neq_;
1477  doublereal Nres2_o = m_normResid_0 * m_normResid_0 * neq_;
1478  doublereal a = Nres2CP / Nres2_o;
1479  doublereal betaEqual = (2.0 - sqrt(4.0 - 4 * (1.0 - a))) / 2.0;
1480  doublereal beta = (1.0 + betaEqual) / 2.0;
1481
1482
1483  Nuu_ = beta;
1484
1486  for (size_t i = 0; i < neq_; i++) {
1487  m_wksp[i] = Nuu_ * deltaX_Newton_[i] - deltaX_CP_[i];
1488  }
1490  dist_R2_ = (1.0 - Nuu_) * m_normDeltaSoln_Newton;
1492
1493  /*
1494  * Calculate the trust distances
1495  */
1498 }
1499
1500 int NonlinearSolver::lambdaToLeg(const doublereal lambda, doublereal& alpha) const
1501 {
1502
1503  if (lambda < dist_R0_ / dist_Total_) {
1504  alpha = lambda * dist_Total_ / dist_R0_;
1505  return 0;
1506  } else if (lambda < ((dist_R0_ + dist_R1_)/ dist_Total_)) {
1507  alpha = (lambda * dist_Total_ - dist_R0_) / dist_R1_;
1508  return 1;
1509  }
1510  alpha = (lambda * dist_Total_ - dist_R0_ - dist_R1_) / dist_R2_;
1511  return 2;
1512 }
1513
1514 doublereal NonlinearSolver::expectedResidLeg(int leg, doublereal alpha) const
1515 {
1516  doublereal resD2, res2, resNorm;
1517  doublereal normResid02 = m_normResid_0 * m_normResid_0 * neq_;
1518
1519  if (leg == 0) {
1520  /*
1521  * We are on the steepest descent line
1522  * along that line
1523  * R2 = R2 + 2 lambda R dot Jd + lambda**2 Jd dot Jd
1524  */
1525
1526  doublereal tmp = - 2.0 * alpha + alpha * alpha;
1527  doublereal tmp2 = - RJd_norm_ * lambdaStar_;
1528  resD2 = tmp2 * tmp;
1529
1530  } else if (leg == 1) {
1531
1532  /*
1533  * Same formula as above for lambda=1.
1534  */
1535  doublereal tmp2 = - RJd_norm_ * lambdaStar_;
1536  doublereal RdotJS = - tmp2;
1537  doublereal JsJs = tmp2;
1538
1539
1540  doublereal res0_2 = m_normResid_0 * m_normResid_0 * neq_;
1541
1542  res2 = res0_2 + (1.0 - alpha) * 2 * RdotJS - 2 * alpha * Nuu_ * res0_2
1543  + (1.0 - alpha) * (1.0 - alpha) * JsJs
1544  + alpha * alpha * Nuu_ * Nuu_ * res0_2
1545  - 2 * alpha * Nuu_ * (1.0 - alpha) * RdotJS;
1546
1547  return sqrt(res2 / neq_);
1548
1549  } else {
1550  doublereal beta = Nuu_ + alpha * (1.0 - Nuu_);
1551  doublereal tmp2 = normResid02;
1552  doublereal tmp = 1.0 - 2.0 * beta + 1.0 * beta * beta - 1.0;
1553  resD2 = tmp * tmp2;
1554  }
1555
1556  res2 = m_normResid_0 * m_normResid_0 * neq_ + resD2;
1557  if (res2 < 0.0) {
1558  resNorm = m_normResid_0 - sqrt(resD2/neq_);
1559  } else {
1560  resNorm = sqrt(res2 / neq_);
1561  }
1562
1563  return resNorm;
1564
1565 }
1566
1567 void NonlinearSolver::residualComparisonLeg(const doublereal time_curr, const doublereal* const ydot0, int& legBest,
1568  doublereal& alphaBest) const
1569 {
1570  doublereal* y1 = DATA_PTR(m_wksp);
1571  doublereal* ydot1 = DATA_PTR(m_wksp_2);
1572  doublereal sLen;
1573  doublereal alpha = 0;
1574
1575  doublereal residSteepBest = 1.0E300;
1576  doublereal residSteepLinBest = 0.0;
1577  if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
1578  printf("\t\t residualComparisonLeg() \n");
1579  printf("\t\t Point StepLen Residual_Actual Residual_Linear RelativeMatch\n");
1580  }
1581  // First compare at 1/4 along SD curve
1582  std::vector<doublereal> alphaT;
1583  alphaT.push_back(0.00);
1584  alphaT.push_back(0.01);
1585  alphaT.push_back(0.1);
1586  alphaT.push_back(0.25);
1587  alphaT.push_back(0.50);
1588  alphaT.push_back(0.75);
1589  alphaT.push_back(1.0);
1590  for (size_t iteration = 0; iteration < alphaT.size(); iteration++) {
1591  alpha = alphaT[iteration];
1592  for (size_t i = 0; i < neq_; i++) {
1593  y1[i] = m_y_n_curr[i] + alpha * deltaX_CP_[i];
1594  }
1596  calc_ydot(m_order, y1, ydot1);
1597  }
1598  sLen = alpha * solnErrorNorm(DATA_PTR(deltaX_CP_));
1599  /*
1600  * Calculate the residual that would result if y1[] were the new solution vector
1601  * -> m_resid[] contains the result of the residual calculation
1602  */
1604  doResidualCalc(time_curr, solnType_, y1, ydot1, Base_LaggedSolutionComponents);
1605  } else {
1606  doResidualCalc(time_curr, solnType_, y1, ydot0, Base_LaggedSolutionComponents);
1607  }
1608
1609
1610  doublereal residSteep = residErrorNorm(DATA_PTR(m_resid));
1611  doublereal residSteepLin = expectedResidLeg(0, alpha);
1612  if (residSteep < residSteepBest) {
1613  legBest = 0;
1614  alphaBest = alpha;
1615  residSteepBest = residSteep;
1616  residSteepLinBest = residSteepLin;
1617  }
1618
1619  doublereal relFit = (residSteep - residSteepLin) / (fabs(residSteepLin) + 1.0E-10);
1620  if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
1621  printf("\t\t (%2d - % 10.3g) % 15.8E % 15.8E % 15.8E % 15.8E\n", 0, alpha, sLen, residSteep, residSteepLin , relFit);
1622  }
1623  }
1624
1625  for (size_t iteration = 0; iteration < alphaT.size(); iteration++) {
1626  doublereal alpha = alphaT[iteration];
1627  for (size_t i = 0; i < neq_; i++) {
1628  y1[i] = m_y_n_curr[i] + (1.0 - alpha) * deltaX_CP_[i];
1629  y1[i] += alpha * Nuu_ * deltaX_Newton_[i];
1630  }
1632  calc_ydot(m_order, y1, ydot1);
1633  }
1634  /*
1635  * Calculate the residual that would result if y1[] were the new solution vector
1636  * -> m_resid[] contains the result of the residual calculation
1637  */
1639  doResidualCalc(time_curr, solnType_, y1, ydot1, Base_LaggedSolutionComponents);
1640  } else {
1641  doResidualCalc(time_curr, solnType_, y1, ydot0, Base_LaggedSolutionComponents);
1642  }
1643
1644  for (size_t i = 0; i < neq_; i++) {
1645  y1[i] -= m_y_n_curr[i];
1646  }
1647  sLen = solnErrorNorm(DATA_PTR(y1));
1648
1649  doublereal residSteep = residErrorNorm(DATA_PTR(m_resid));
1650  doublereal residSteepLin = expectedResidLeg(1, alpha);
1651  if (residSteep < residSteepBest) {
1652  legBest = 1;
1653  alphaBest = alpha;
1654  residSteepBest = residSteep;
1655  residSteepLinBest = residSteepLin;
1656  }
1657
1658  doublereal relFit = (residSteep - residSteepLin) / (fabs(residSteepLin) + 1.0E-10);
1659  if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
1660  printf("\t\t (%2d - % 10.3g) % 15.8E % 15.8E % 15.8E % 15.8E\n", 1, alpha, sLen, residSteep, residSteepLin , relFit);
1661  }
1662  }
1663
1664  for (size_t iteration = 0; iteration < alphaT.size(); iteration++) {
1665  doublereal alpha = alphaT[iteration];
1666  for (size_t i = 0; i < neq_; i++) {
1667  y1[i] = m_y_n_curr[i] + (Nuu_ + alpha * (1.0 - Nuu_))* deltaX_Newton_[i];
1668  }
1670  calc_ydot(m_order, y1, ydot1);
1671  }
1672  sLen = (Nuu_ + alpha * (1.0 - Nuu_)) * solnErrorNorm(DATA_PTR(deltaX_Newton_));
1673  /*
1674  * Calculate the residual that would result if y1[] were the new solution vector
1675  * -> m_resid[] contains the result of the residual calculation
1676  */
1678  doResidualCalc(time_curr, solnType_, y1, ydot1, Base_LaggedSolutionComponents);
1679  } else {
1680  doResidualCalc(time_curr, solnType_, y1, ydot0, Base_LaggedSolutionComponents);
1681  }
1682
1683
1684
1685  doublereal residSteep = residErrorNorm(DATA_PTR(m_resid));
1686  doublereal residSteepLin = expectedResidLeg(2, alpha);
1687  if (residSteep < residSteepBest) {
1688  legBest = 2;
1689  alphaBest = alpha;
1690  residSteepBest = residSteep;
1691  residSteepLinBest = residSteepLin;
1692  }
1693  doublereal relFit = (residSteep - residSteepLin) / (fabs(residSteepLin) + 1.0E-10);
1694  if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
1695  printf("\t\t (%2d - % 10.3g) % 15.8E % 15.8E % 15.8E % 15.8E\n", 2, alpha, sLen, residSteep, residSteepLin , relFit);
1696  }
1697  }
1698  if (s_print_DogLeg || (doDogLeg_ && m_print_flag > 6)) {
1699  printf("\t\t Best Result: \n");
1700  doublereal relFit = (residSteepBest - residSteepLinBest) / (fabs(residSteepLinBest) + 1.0E-10);
1701  if (m_print_flag <= 6) {
1702  printf("\t\t Leg %2d alpha %5g: NonlinResid = %g LinResid = %g, relfit = %g\n",
1703  legBest, alphaBest, residSteepBest, residSteepLinBest, relFit);
1704  } else {
1705  if (legBest == 0) {
1706  sLen = alpha * solnErrorNorm(DATA_PTR(deltaX_CP_));
1707  } else if (legBest == 1) {
1708  for (size_t i = 0; i < neq_; i++) {
1709  y1[i] = (1.0 - alphaBest) * deltaX_CP_[i];
1710  y1[i] += alphaBest * Nuu_ * deltaX_Newton_[i];
1711  }
1712  sLen = solnErrorNorm(DATA_PTR(y1));
1713  } else {
1714  sLen = (Nuu_ + alpha * (1.0 - Nuu_)) * solnErrorNorm(DATA_PTR(deltaX_Newton_));
1715  }
1716  printf("\t\t (%2d - % 10.3g) % 15.8E % 15.8E % 15.8E % 15.8E\n", legBest, alphaBest, sLen,
1717  residSteepBest, residSteepLinBest , relFit);
1718  }
1719  }
1720
1721 }
1722
1724 {
1727 }
1728
1730 {
1731  for (size_t i = 0; i < neq_; i++) {
1732  m_deltaStepMinimum[i] = 1000. * atolk_[i];
1733  m_deltaStepMinimum[i] = std::max(m_deltaStepMinimum[i], 0.1 * fabs(m_y_n_curr[i]));
1734  }
1735 }
1736
1738 {
1739  for (size_t i = 0; i < neq_; i++) {
1740  doublereal goodVal = deltaX_trust_[i] * trustDelta_;
1741  if (deltaX_trust_[i] * trustDelta_ > m_deltaStepMinimum[i]) {
1742  m_deltaStepMinimum[i] = 1.1 * goodVal;
1743  }
1744
1745  }
1746 }
1747
1748 void NonlinearSolver::setDeltaBoundsMagnitudes(const doublereal* const deltaStepMinimum)
1749 {
1750  for (size_t i = 0; i < neq_; i++) {
1751  m_deltaStepMinimum[i] = deltaStepMinimum[i];
1752  }
1754 }
1755
1756 double
1757 NonlinearSolver::deltaBoundStep(const doublereal* const y_n_curr, const doublereal* const step_1)
1758 {
1759  size_t i_fbounds = 0;
1760  int ifbd = 0;
1761  int i_fbd = 0;
1762  doublereal UPFAC = 2.0;
1763
1764  doublereal sameSign = 0.0;
1765  doublereal ff;
1766  doublereal f_delta_bounds = 1.0;
1767  doublereal ff_alt;
1768  for (size_t i = 0; i < neq_; i++) {
1769  doublereal y_new = y_n_curr[i] + step_1[i];
1770  sameSign = y_new * y_n_curr[i];
1771
1772  /*
1773  * Now do a delta bounds
1774  * Increase variables by a factor of UPFAC only
1775  * decrease variables by a factor of 2 only
1776  */
1777  ff = 1.0;
1778
1779
1780  if (sameSign >= 0.0) {
1781  if ((fabs(y_new) > UPFAC * fabs(y_n_curr[i])) &&
1782  (fabs(y_new - y_n_curr[i]) > m_deltaStepMinimum[i])) {
1783  ff = (UPFAC - 1.0) * fabs(y_n_curr[i]/(y_new - y_n_curr[i]));
1784  ff_alt = fabs(m_deltaStepMinimum[i] / (y_new - y_n_curr[i]));
1785  ff = std::max(ff, ff_alt);
1786  ifbd = 1;
1787  }
1788  if ((fabs(2.0 * y_new) < fabs(y_n_curr[i])) &&
1789  (fabs(y_new - y_n_curr[i]) > m_deltaStepMinimum[i])) {
1790  ff = y_n_curr[i]/(y_new - y_n_curr[i]) * (1.0 - 2.0)/2.0;
1791  ff_alt = fabs(m_deltaStepMinimum[i] / (y_new - y_n_curr[i]));
1792  ff = std::max(ff, ff_alt);
1793  ifbd = 0;
1794  }
1795  } else {
1796  /*
1797  * This handles the case where the value crosses the origin.
1798  * - First we don't let it cross the origin until its shrunk to the size of m_deltaStepMinimum[i]
1799  */
1800  if (fabs(y_n_curr[i]) > m_deltaStepMinimum[i]) {
1801  ff = y_n_curr[i]/(y_new - y_n_curr[i]) * (1.0 - 2.0)/2.0;
1802  ff_alt = fabs(m_deltaStepMinimum[i] / (y_new - y_n_curr[i]));
1803  ff = std::max(ff, ff_alt);
1804  if (y_n_curr[i] >= 0.0) {
1805  ifbd = 0;
1806  } else {
1807  ifbd = 1;
1808  }
1809  }
1810  /*
1811  * Second when it does cross the origin, we make sure that its magnitude is only 50% of the previous value.
1812  */
1813  else if (fabs(y_new) > 0.5 * fabs(y_n_curr[i])) {
1814  ff = y_n_curr[i]/(y_new - y_n_curr[i]) * (-1.5);
1815  ff_alt = fabs(m_deltaStepMinimum[i] / (y_new - y_n_curr[i]));
1816  ff = std::max(ff, ff_alt);
1817  ifbd = 0;
1818  }
1819  }
1820
1821  if (ff < f_delta_bounds) {
1822  f_delta_bounds = ff;
1823  i_fbounds = i;
1824  i_fbd = ifbd;
1825  }
1826
1827
1828  }
1829
1830
1831  /*
1832  * Report on any corrections
1833  */
1834  if (m_print_flag >= 3) {
1835  if (f_delta_bounds < 1.0) {
1836  if (i_fbd) {
1837  printf("\t\tdeltaBoundStep: Increase of Variable %s causing "
1838  "delta damping of %g: origVal = %10.3g, undampedNew = %10.3g, dampedNew = %10.3g\n",
1839  int2str(i_fbounds).c_str(), f_delta_bounds, y_n_curr[i_fbounds], y_n_curr[i_fbounds] + step_1[i_fbounds],
1840  y_n_curr[i_fbounds] + f_delta_bounds * step_1[i_fbounds]);
1841  } else {
1842  printf("\t\tdeltaBoundStep: Decrease of variable %s causing"
1843  "delta damping of %g: origVal = %10.3g, undampedNew = %10.3g, dampedNew = %10.3g\n",
1844  int2str(i_fbounds).c_str(), f_delta_bounds, y_n_curr[i_fbounds], y_n_curr[i_fbounds] + step_1[i_fbounds],
1845  y_n_curr[i_fbounds] + f_delta_bounds * step_1[i_fbounds]);
1846  }
1847  }
1848  }
1849
1850
1851  return f_delta_bounds;
1852 }
1853
1855 {
1856  doublereal trustDeltaOld = trustDelta_;
1857  doublereal wtSum = 0.0;
1858  for (size_t i = 0; i < neq_; i++) {
1859  wtSum += m_ewt[i];
1860  }
1861  wtSum /= neq_;
1862  doublereal trustNorm = solnErrorNorm(DATA_PTR(deltaX_trust_));
1863  doublereal deltaXSizeOld = trustNorm;
1864  doublereal trustNormGoal = trustNorm * trustDelta_;
1865
1866  doublereal oldVal;
1867  doublereal fabsy;
1868  // we use the old value of the trust region as an indicator
1869  for (size_t i = 0; i < neq_; i++) {
1870  oldVal = deltaX_trust_[i];
1871  fabsy = fabs(m_y_n_curr[i]);
1872  // First off make sure that each trust region vector is 1/2 the size of each variable or smaller
1873  // unless overridden by the deltaStepMininum value.
1874  doublereal newValue = trustNormGoal * m_ewt[i];
1875  if (newValue > 0.5 * fabsy) {
1876  if (fabsy * 0.5 > m_deltaStepMinimum[i]) {
1877  deltaX_trust_[i] = 0.5 * fabsy;
1878  } else {
1880  }
1881  } else {
1882  if (newValue > 4.0 * oldVal) {
1883  newValue = 4.0 * oldVal;
1884  } else if (newValue < 0.25 * oldVal) {
1885  newValue = 0.25 * oldVal;
1886  if (deltaX_trust_[i] < m_deltaStepMinimum[i]) {
1887  newValue = m_deltaStepMinimum[i];
1888  }
1889  }
1890  deltaX_trust_[i] = newValue;
1891  if (deltaX_trust_[i] > 0.75 * m_deltaStepMaximum[i]) {
1892  deltaX_trust_[i] = 0.75 * m_deltaStepMaximum[i];
1893  }
1894  }
1895  }
1896
1897
1898  // Final renormalization.
1900  doublereal sum = trustNormGoal / trustNorm;
1901  for (size_t i = 0; i < neq_; i++) {
1902  deltaX_trust_[i] = deltaX_trust_[i] * sum;
1903  }
1905  trustDelta_ = trustNormGoal / norm_deltaX_trust_;
1906
1907  if (doDogLeg_ && m_print_flag >= 4) {
1908  printf("\t\t reajustTrustVector(): Trust size = %11.3E: Old deltaX size = %11.3E trustDelta_ = %11.3E\n"
1909  "\t\t new deltaX size = %11.3E trustdelta_ = %11.3E\n",
1910  trustNormGoal, deltaXSizeOld, trustDeltaOld, norm_deltaX_trust_, trustDelta_);
1911  }
1912 }
1913
1915 {
1917  return;
1918  }
1920  for (size_t i = 0; i < neq_; i++) {
1922  }
1923  trustDelta_ = 1.0;
1924  }
1926  for (size_t i = 0; i < neq_; i++) {
1928  }
1929  doublereal cpd = calcTrustDistance(deltaX_CP_);
1930  if ((doDogLeg_ && m_print_flag >= 4)) {
1931  printf("\t\t initializeTrustRegion(): Relative Distance of Cauchy Vector wrt Trust Vector = %g\n", cpd);
1932  }
1936  if ((doDogLeg_ && m_print_flag >= 4)) {
1937  printf("\t\t initializeTrustRegion(): Relative Distance of Cauchy Vector wrt Trust Vector = %g\n", cpd);
1938  }
1939  }
1941  for (size_t i = 0; i < neq_; i++) {
1943  }
1944  doublereal cpd = calcTrustDistance(deltaX_Newton_);
1945  if ((doDogLeg_ && m_print_flag >= 4)) {
1946  printf("\t\t initializeTrustRegion(): Relative Distance of Newton Vector wrt Trust Vector = %g\n", cpd);
1947  }
1948  trustDelta_ = trustDelta_ * cpd;
1951  if ((doDogLeg_ && m_print_flag >= 4)) {
1952  printf("\t\t initializeTrustRegion(): Relative Distance of Newton Vector wrt Trust Vector = %g\n", cpd);
1953  }
1954  }
1955 }
1956
1957 void NonlinearSolver::fillDogLegStep(int leg, doublereal alpha, std::vector<doublereal> & deltaX) const
1958 {
1959  if (leg == 0) {
1960  for (size_t i = 0; i < neq_; i++) {
1961  deltaX[i] = alpha * deltaX_CP_[i];
1962  }
1963  } else if (leg == 2) {
1964  for (size_t i = 0; i < neq_; i++) {
1965  deltaX[i] = (alpha + (1.0 - alpha) * Nuu_) * deltaX_Newton_[i];
1966  }
1967  } else {
1968  for (size_t i = 0; i < neq_; i++) {
1969  deltaX[i] = deltaX_CP_[i] * (1.0 - alpha) + alpha * Nuu_ * deltaX_Newton_[i];
1970  }
1971  }
1972 }
1973
1974 doublereal NonlinearSolver::calcTrustDistance(std::vector<doublereal> const& deltaX) const
1975 {
1976  doublereal sum = 0.0;
1977  doublereal tmp = 0.0;
1978  for (size_t i = 0; i < neq_; i++) {
1979  tmp = deltaX[i] / deltaX_trust_[i];
1980  sum += tmp * tmp;
1981  }
1982  return sqrt(sum / neq_) / trustDelta_;
1983 }
1984
1985 int NonlinearSolver::calcTrustIntersection(doublereal trustDelta, doublereal& lambda, doublereal& alpha) const
1986 {
1987  doublereal dist;
1988  if (normTrust_Newton_ < trustDelta) {
1989  lambda = 1.0;
1990  alpha = 1.0;
1991  return 2;
1992  }
1993
1994  if (normTrust_Newton_ * Nuu_ < trustDelta) {
1995  alpha = (trustDelta - normTrust_Newton_ * Nuu_) / (normTrust_Newton_ - normTrust_Newton_ * Nuu_);
1996  dist = dist_R0_ + dist_R1_ + alpha * dist_R2_;
1997  lambda = dist / dist_Total_;
1998  return 2;
1999  }
2000  if (normTrust_CP_ > trustDelta) {
2001  lambda = 1.0;
2002  dist = dist_R0_ * trustDelta / normTrust_CP_;
2003  lambda = dist / dist_Total_;
2004  alpha = trustDelta / normTrust_CP_;
2005  return 0;
2006  }
2007  doublereal sumv = 0.0;
2008  for (size_t i = 0; i < neq_; i++) {
2009  sumv += (deltaX_Newton_[i] / deltaX_trust_[i]) * (deltaX_CP_[i] / deltaX_trust_[i]);
2010  }
2011
2012  doublereal a = normTrust_Newton_ * normTrust_Newton_ * Nuu_ * Nuu_;
2013  doublereal b = 2.0 * Nuu_ * sumv;
2014  doublereal c = normTrust_CP_ * normTrust_CP_ - trustDelta * trustDelta;
2015
2016  alpha =(-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
2017
2018
2019  dist = dist_R0_ + alpha * dist_R1_;
2020  lambda = dist / dist_Total_;
2021  return 1;
2022 }
2023
2024 doublereal NonlinearSolver::boundStep(const doublereal* const y, const doublereal* const step0)
2025 {
2026  size_t i_lower = npos;
2027  doublereal f_bounds = 1.0;
2028  doublereal ff, y_new;
2029
2030  for (size_t i = 0; i < neq_; i++) {
2031  y_new = y[i] + step0[i];
2032  /*
2033  * Force the step to only take 80% a step towards the lower bounds
2034  */
2035  if (step0[i] < 0.0) {
2036  if (y_new < (y[i] + 0.8 * (m_y_low_bounds[i] - y[i]))) {
2037  doublereal legalDelta = 0.8*(m_y_low_bounds[i] - y[i]);
2038  ff = legalDelta / step0[i];
2039  if (ff < f_bounds) {
2040  f_bounds = ff;
2041  i_lower = i;
2042  }
2043  }
2044  }
2045  /*
2046  * Force the step to only take 80% a step towards the high bounds
2047  */
2048  if (step0[i] > 0.0) {
2049  if (y_new > (y[i] + 0.8 * (m_y_high_bounds[i] - y[i]))) {
2050  doublereal legalDelta = 0.8*(m_y_high_bounds[i] - y[i]);
2051  ff = legalDelta / step0[i];
2052  if (ff < f_bounds) {
2053  f_bounds = ff;
2054  i_lower = i;
2055  }
2056  }
2057  }
2058
2059  }
2060
2061  /*
2062  * Report on any corrections
2063  */
2064  if (m_print_flag >= 3) {
2065  if (f_bounds != 1.0) {
2066  printf("\t\tboundStep: Variable %s causing bounds damping of %g\n", int2str(i_lower).c_str(), f_bounds);
2067  }
2068  }
2069
2070  doublereal f_delta_bounds = deltaBoundStep(y, step0);
2071  return std::min(f_bounds, f_delta_bounds);
2072 }
2073
2074 int NonlinearSolver::dampStep(const doublereal time_curr, const doublereal* const y_n_curr,
2075  const doublereal* const ydot_n_curr, doublereal* const step_1,
2076  doublereal* const y_n_1, doublereal* const ydot_n_1, doublereal* const step_2,
2077  doublereal& stepNorm_2, GeneralMatrix& jac, bool writetitle, int& num_backtracks)
2078 {
2079  int m;
2080  int info = 0;
2081  int retnTrial = NSOLN_RETN_FAIL_DAMPSTEP;
2082  // Compute the weighted norm of the undamped step size step_1
2083  doublereal stepNorm_1 = solnErrorNorm(step_1);
2084
2085  doublereal* step_1_orig = DATA_PTR(m_wksp);
2086  for (size_t j = 0; j < neq_; j++) {
2087  step_1_orig[j] = step_1[j];
2088  }
2089
2090
2091  // Compute the multiplier to keep all components in bounds.A value of one indicates that there is no limitation
2092  // on the current step size in the nonlinear method due to bounds constraints (either negative values of delta
2093  // bounds constraints.
2094  m_dampBound = boundStep(y_n_curr, step_1);
2095
2096  // If fbound is very small, then y0 is already close to the boundary and step0 points out of the allowed domain. In
2097  // this case, the Newton algorithm fails, so return an error condition.
2098  if (m_dampBound < 1.e-30) {
2099  if (m_print_flag > 1) {
2100  printf("\t\t\tdampStep(): At limits.\n");
2101  }
2102  return -3;
2103  }
2104
2105  //--------------------------------------------
2106  // Attempt damped step
2107  //--------------------------------------------
2108
2109  // damping coefficient starts at 1.0
2110  m_dampRes = 1.0;
2111
2112  doublereal ff = m_dampBound;
2113  num_backtracks = 0;
2114  for (m = 0; m < NDAMP; m++) {
2115
2116  ff = m_dampBound * m_dampRes;
2117
2118  // step the solution by the damped step size
2119  /*
2120  * Whenever we update the solution, we must also always
2121  * update the time derivative.
2122  */
2123  for (size_t j = 0; j < neq_; j++) {
2124  step_1[j] = ff * step_1_orig[j];
2125  y_n_1[j] = y_n_curr[j] + step_1[j];
2126  }
2127
2129  calc_ydot(m_order, y_n_1, ydot_n_1);
2130  }
2131  /*
2132  * Calculate the residual that would result if y1[] were the new solution vector
2133  * -> m_resid[] contains the result of the residual calculation
2134  */
2136  info = doResidualCalc(time_curr, solnType_, y_n_1, ydot_n_1, Base_LaggedSolutionComponents);
2137  } else {
2138  info = doResidualCalc(time_curr, solnType_, y_n_1, ydot_n_curr, Base_LaggedSolutionComponents);
2139  }
2140  if (info != 1) {
2141  if (m_print_flag > 0) {
2142  printf("\t\t\tdampStep(): current trial step and damping led to Residual Calc ERROR %d. Bailing\n", info);
2143  }
2144  return -1;
2145  }
2148  if (m == 0) {
2150  }
2151
2152  bool steepEnough = (m_normResidTrial < m_normResid_0 * (0.9 * (1.0 - ff) * (1.0 - ff)* (1.0 - ff) + 0.1));
2153
2154  if (m_normResidTrial < 1.0 || steepEnough) {
2155  if (m_print_flag >= 5) {
2156  if (m_normResidTrial < 1.0) {
2157  printf("\t dampStep(): Current trial step and damping"
2158  " coefficient accepted because residTrial test step < 1:\n");
2159  printf("\t resid0 = %g, residTrial = %g\n", m_normResid_0, m_normResidTrial);
2160  } else if (steepEnough) {
2161  printf("\t dampStep(): Current trial step and damping"
2162  " coefficient accepted because resid0 > residTrial and steep enough:\n");
2163  printf("\t resid0 = %g, residTrial = %g\n", m_normResid_0, m_normResidTrial);
2164  } else {
2165  printf("\t dampStep(): Current trial step and damping"
2166  " coefficient accepted because residual solution damping is turned off:\n");
2167  printf("\t resid0 = %g, residTrial = %g\n", m_normResid_0, m_normResidTrial);
2168  }
2169  }
2170  /*
2171  * We aren't going to solve the system if we don't need to. Therefore, return an estimate
2172  * of the next solution update based on the ratio of the residual reduction.
2173  */
2174  if (m_normResid_0 > 0.0) {
2175  stepNorm_2 = stepNorm_1 * m_normResidTrial / m_normResid_0;
2176  } else {
2177  stepNorm_2 = 0;
2178  }
2179  if (m_normResidTrial < 1.0) {
2180  retnTrial = 3;
2181  } else {
2182  retnTrial = 4;
2183  }
2184  break;
2185  }
2186
2187  // Compute the next undamped step, step1[], that would result if y1[] were accepted.
2188  // We now have two steps that we have calculated step0[] and step1[]
2190  info = doNewtonSolve(time_curr, y_n_1, ydot_n_1, step_2, jac);
2191  } else {
2192  info = doNewtonSolve(time_curr, y_n_1, ydot_n_curr, step_2, jac);
2193  }
2194  if (info) {
2195  if (m_print_flag > 0) {
2196  printf("\t\t\tdampStep: current trial step and damping led to LAPACK ERROR %d. Bailing\n", info);
2197  }
2198  return -1;
2199  }
2200
2201  // compute the weighted norm of step1
2202  stepNorm_2 = solnErrorNorm(step_2);
2203
2205  if (m_print_flag >= 5) {
2206  print_solnDelta_norm_contrib((const doublereal*) step_1_orig, "DeltaSoln",
2207  (const doublereal*) step_2, "DeltaSolnTrial",
2208  "dampNewt: Important Entries for Weighted Soln Updates:",
2209  y_n_curr, y_n_1, ff, 5);
2210  }
2211  if (m_print_flag >= 4) {
2212  printf("\t\t\tdampStep(): s1 = %g, s2 = %g, dampBound = %g,"
2213  "dampRes = %g\n", stepNorm_1, stepNorm_2, m_dampBound, m_dampRes);
2214  }
2215
2216
2217  // if the norm of s1 is less than the norm of s0, then
2218  // accept this damping coefficient. Also accept it if this
2219  // step would result in a converged solution. Otherwise,
2220  // decrease the damping coefficient and try again.
2221
2222  if (stepNorm_2 < 0.8 || stepNorm_2 < stepNorm_1) {
2223  if (stepNorm_2 < 1.0) {
2224  if (m_print_flag >= 3) {
2225  if (stepNorm_2 < 1.0) {
2226  printf("\t\t\tdampStep: current trial step and damping coefficient accepted because test step < 1\n");
2227  printf("\t\t\t s2 = %g, s1 = %g\n", stepNorm_2, stepNorm_1);
2228  }
2229  }
2230  retnTrial = 2;
2231  } else {
2232  retnTrial = 1;
2233  }
2234  break;
2235  } else {
2236  if (m_print_flag > 1) {
2237  printf("\t\t\tdampStep: current step rejected: (s1 = %g > "
2238  "s0 = %g)", stepNorm_2, stepNorm_1);
2239  if (m < (NDAMP-1)) {
2240  printf(" Decreasing damping factor and retrying");
2241  } else {
2242  printf(" Giving up!!!");
2243  }
2244  printf("\n");
2245  }
2246  }
2247  num_backtracks++;
2248  m_dampRes /= DampFactor;
2249  }
2250
2251  // If a damping coefficient was found, return 1 if the
2252  // solution after stepping by the damped step would represent
2253  // a converged solution, and return 0 otherwise. If no damping
2254  // coefficient could be found, return NSOLN_RETN_FAIL_DAMPSTEP.
2255  if (m < NDAMP) {
2256  if (m_print_flag >= 4) {
2257  printf("\t dampStep(): current trial step accepted retnTrial = %d, its = %d, damp = %g\n", retnTrial, m+1, ff);
2258  }
2259  return retnTrial;
2260  } else {
2261  if (stepNorm_2 < 0.5 && (stepNorm_1 < 0.5)) {
2262  if (m_print_flag >= 4) {
2263  printf("\t dampStep(): current trial step accepted kindof retnTrial = %d, its = %d, damp = %g\n", 2, m+1, ff);
2264  }
2265  return 2;
2266  }
2267  if (stepNorm_2 < 1.0) {
2268  if (m_print_flag >= 4) {
2269  printf("\t dampStep(): current trial step accepted and soln converged retnTrial ="
2270  "%d, its = %d, damp = %g\n", 0, m+1, ff);
2271  }
2272  return 0;
2273  }
2274  }
2275  if (m_print_flag >= 4) {
2276  printf("\t dampStep(): current direction is rejected! retnTrial = %d, its = %d, damp = %g\n",
2277  NSOLN_RETN_FAIL_DAMPSTEP, m+1, ff);
2278  }
2279  return NSOLN_RETN_FAIL_DAMPSTEP;
2280 }
2281
2282 int NonlinearSolver::dampDogLeg(const doublereal time_curr, const doublereal* y_n_curr,
2283  const doublereal* ydot_n_curr, std::vector<doublereal> & step_1,
2284  doublereal* const y_n_1, doublereal* const ydot_n_1,
2285  doublereal& stepNorm_1, doublereal& stepNorm_2, GeneralMatrix& jac, int& numTrials)
2286 {
2287  doublereal lambda;
2288  int info;
2289
2290  bool success = false;
2291  bool haveASuccess = false;
2292  doublereal trustDeltaOld = trustDelta_;
2293  doublereal* stepLastGood = DATA_PTR(m_wksp);
2294  //--------------------------------------------
2295  // Attempt damped step
2296  //--------------------------------------------
2297
2298  // damping coefficient starts at 1.0
2299  m_dampRes = 1.0;
2300  int m;
2301  doublereal tlen;
2302
2303
2304  for (m = 0; m < NDAMP; m++) {
2305  numTrials++;
2306  /*
2307  * Find the initial value of lambda that satisfies the trust distance, trustDelta_
2308  */
2310  if (m_print_flag >= 4) {
2311  tlen = trustRegionLength();
2312  printf("\t\t dampDogLeg: trust region with length %13.5E has intersection at leg = %d, alpha = %g, lambda = %g\n",
2313  tlen, dogLegID_, dogLegAlpha_, lambda);
2314  }
2315  /*
2316  * Figure out the new step vector, step_1, based on (leg, alpha). Here we are using the
2317  * intersection of the trust oval with the dog-leg curve.
2318  */
2320
2321  /*
2322  * OK, now that we have step0, Bound the step
2323  */
2324  m_dampBound = boundStep(y_n_curr, DATA_PTR(step_1));
2325  /*
2326  * Decrease the step length if we are bound
2327  */
2328  if (m_dampBound < 1.0) {
2329  for (size_t j = 0; j < neq_; j++) {
2330  step_1[j] = step_1[j] * m_dampBound;
2331  }
2332  }
2333  /*
2334  * Calculate the new solution value y1[] given the step size
2335  */
2336  for (size_t j = 0; j < neq_; j++) {
2337  y_n_1[j] = y_n_curr[j] + step_1[j];
2338  }
2339  /*
2340  * Calculate the new solution time derivative given the step size
2341  */
2343  calc_ydot(m_order, y_n_1, ydot_n_1);
2344  }
2345  /*
2346  * OK, we have the step0. Now, ask the question whether it satisfies the acceptance criteria
2347  * as a good step. The overall outcome is returned in the variable info.
2348  */
2349  info = decideStep(time_curr, dogLegID_, dogLegAlpha_, y_n_curr, ydot_n_curr, step_1,
2350  y_n_1, ydot_n_1, trustDeltaOld);
2352
2353  /*
2354  * The algorithm failed to find a solution vector sufficiently different than the current point
2355  */
2356  if (info == -1) {
2357
2358  if (m_print_flag >= 1) {
2359  doublereal stepNorm = solnErrorNorm(DATA_PTR(step_1));
2360  printf("\t\t dampDogLeg: Current direction rejected, update became too small %g\n", stepNorm);
2361  success = false;
2362  break;
2363  }
2364  }
2365  if (info == -2) {
2366  if (m_print_flag >= 1) {
2367  printf("\t\t dampDogLeg: current trial step and damping led to LAPACK ERROR %d. Bailing\n", info);
2368  success = false;
2369  break;
2370  }
2371  }
2372  if (info == 0) {
2373  success = true;
2374  break;
2375  }
2376  if (info == 3) {
2377
2378  haveASuccess = true;
2379  // Store the good results in stepLastGood
2380  copy(step_1.begin(), step_1.end(), stepLastGood);
2381  // Within the program decideStep(), we have already increased the value of trustDelta_. We store the
2382  // value of step0 in step1, recalculate a larger step0 in the next fillDogLegStep(),
2383  // and then attempt to see if the larger step works in the next iteration
2384  }
2385  if (info == 2) {
2386  // Step was a failure. If we had a previous success with a smaller stepsize, haveASuccess is true
2387  // and we execute the next block and break. If we didn't have a previous success, trustDelta_ has
2388  // already been decreased in the decideStep() routine. We go back and try another iteration with
2389  // a smaller trust region.
2390  if (haveASuccess) {
2391  copy(stepLastGood, stepLastGood+neq_, step_1.begin());
2392  for (size_t j = 0; j < neq_; j++) {
2393  y_n_1[j] = y_n_curr[j] + step_1[j];
2394  }
2396  calc_ydot(m_order, y_n_1, ydot_n_1);
2397  }
2398  success = true;
2399  break;
2400  } else {
2401
2402  }
2403  }
2404  }
2405
2406  /*
2407  * Estimate s1, the norm after the next step
2408  */
2409  stepNorm_1 = solnErrorNorm(DATA_PTR(step_1));
2410  stepNorm_2 = stepNorm_1;
2411  if (m_dampBound < 1.0) {
2412  stepNorm_2 /= m_dampBound;
2413  }
2414  stepNorm_2 /= lambda;
2415  stepNorm_2 *= m_normResidTrial / m_normResid_0;
2416
2417
2418  if (success) {
2419  if (m_normResidTrial < 1.0) {
2420  if (normTrust_Newton_ < trustDelta_ && m_dampBound == 1.0) {
2421  return 1;
2422  } else {
2423  return 0;
2424  }
2425  }
2426  return 0;
2427  }
2428  return NSOLN_RETN_FAIL_DAMPSTEP;
2429 }
2430
2431 int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal alpha,
2432  const doublereal* const y_n_curr,
2433  const doublereal* const ydot_n_curr, const std::vector<doublereal> & step_1,
2434  const doublereal* const y_n_1, const doublereal* const ydot_n_1,
2435  doublereal trustDeltaOld)
2436 {
2437  int retn = 2;
2438  int info;
2439  doublereal ll;
2440  // Calculate the solution step length
2441  doublereal stepNorm = solnErrorNorm(DATA_PTR(step_1));
2442
2443  // Calculate the initial (R**2 * neq) value for the old function
2444  doublereal normResid0_2 = m_normResid_0 * m_normResid_0 * neq_;
2445
2446  // Calculate the distance to the Cauchy point
2447  doublereal cauchyDistanceNorm = solnErrorNorm(DATA_PTR(deltaX_CP_));
2448
2449  // This is the expected initial rate of decrease in the Cauchy direction.
2450  // -> This is Eqn. 29 = Rhat dot Jhat dy / || d ||
2451  doublereal funcDecreaseSDExp = RJd_norm_ / cauchyDistanceNorm * lambdaStar_;
2452  if (funcDecreaseSDExp > 0.0) {
2453  if (m_print_flag >= 5) {
2454  printf("\t\tdecideStep(): Unexpected condition -> Cauchy slope is positive\n");
2455  }
2456  }
2457
2458  /*
2459  * Calculate the residual that would result if y1[] were the new solution vector.
2460  * The Lagged solution components are kept lagged here. Unfortunately, it just doesn't work in some cases to use a
2461  * Jacobian from a lagged state and then use a residual from an unlagged condition. The linear model doesn't
2462  * agree with the nonlinear model.
2463  * -> m_resid[] contains the result of the residual calculation
2464  */
2466  info = doResidualCalc(time_curr, solnType_, y_n_1, ydot_n_1, Base_LaggedSolutionComponents);
2467  } else {
2468  info = doResidualCalc(time_curr, solnType_, y_n_1, ydot_n_curr, Base_LaggedSolutionComponents);
2469  }
2470
2471  if (info != 1) {
2472  if (m_print_flag >= 2) {
2473  printf("\t\tdecideStep: current trial step and damping led to Residual Calc ERROR %d. Bailing\n", info);
2474  }
2475  return -2;
2476  }
2477  /*
2478  * Ok we have a successful new residual. Calculate the normalized residual value and store it in
2479  * m_normResidTrial
2480  */
2482  doublereal normResidTrial_2 = neq_ * m_normResidTrial * m_normResidTrial;
2483
2484  /*
2485  * We have a minimal acceptance test for passage. deltaf < 1.0E-4 (CauchySlope) (deltS)
2486  * This is the condition that D&S use in 6.4.5
2487  */
2488  doublereal funcDecrease = 0.5 * (normResidTrial_2 - normResid0_2);
2489  doublereal acceptableDelF = funcDecreaseSDExp * stepNorm * 1.0E-4;
2490  if (funcDecrease < acceptableDelF) {
2493  retn = 0;
2494  if (m_print_flag >= 4) {
2495  printf("\t\t decideStep: Norm Residual(leg=%1d, alpha=%10.2E) = %11.4E passes\n",
2496  dogLegID_, dogLegAlpha_, m_normResidTrial);
2497  }
2498  } else {
2499  if (m_print_flag >= 4) {
2500  printf("\t\t decideStep: Norm Residual(leg=%1d, alpha=%10.2E) = %11.4E failes\n",
2501  dogLegID_, dogLegAlpha_, m_normResidTrial);
2502  }
2503  trustDelta_ *= 0.33;
2504  CurrentTrustFactor_ *= 0.33;
2505  retn = 2;
2506  // error condition if step is getting too small
2507  if (rtol_ * stepNorm < 1.0E-6) {
2508  retn = -1;
2509  }
2510  return retn;
2511  }
2512  /*
2513  * Figure out the next trust region. We are here iff retn = 0
2514  *
2515  * If we had to bounds delta the update, decrease the trust region
2516  */
2517  if (m_dampBound >= 1.0) {
2518  retn = 0;
2519  /*
2520  * Calculate the expected residual from the quadratic model
2521  */
2522  doublereal expectedNormRes = expectedResidLeg(leg, alpha);
2523  doublereal expectedFuncDecrease = 0.5 * (neq_ * expectedNormRes * expectedNormRes - normResid0_2);
2524  if (funcDecrease > 0.1 * expectedFuncDecrease) {
2525  if ((m_normResidTrial > 0.5 * m_normResid_0) && (m_normResidTrial > 0.1)) {
2526  trustDelta_ *= 0.5;
2527  NextTrustFactor_ *= 0.5;
2528  ll = trustRegionLength();
2529  if (m_print_flag >= 4) {
2530  printf("\t\t decideStep: Trust region decreased from %g to %g due to bad quad approximation\n",
2531  ll*2, ll);
2532  }
2533  }
2534  } else {
2535  /*
2536  * If we are doing well, consider increasing the trust region and recalculating
2537  */
2538  if (funcDecrease < 0.8 * expectedFuncDecrease || (m_normResidTrial < 0.33 * m_normResid_0)) {
2539  if (trustDelta_ <= trustDeltaOld && (leg != 2 || alpha < 0.75)) {
2540  trustDelta_ *= 2.0;
2541  CurrentTrustFactor_ *= 2;
2543  ll = trustRegionLength();
2544  if (m_print_flag >= 4) {
2545  if (m_normResidTrial < 0.33 * m_normResid_0) {
2546  printf("\t\t decideStep: Redo line search with trust region increased from %g to %g due to good nonlinear behavior\n",
2547  ll*0.5, ll);
2548  } else {
2549  printf("\t\t decideStep: Redi line search with trust region increased from %g to %g due to good linear model approximation\n",
2550  ll*0.5, ll);
2551  }
2552  }
2553  retn = 3;
2554  } else {
2555  /*
2556  * Increase the size of the trust region for the next calculation
2557  */
2558  if (m_normResidTrial < 0.99 * expectedNormRes || (m_normResidTrial < 0.20 * m_normResid_0) ||
2559  (funcDecrease < -1.0E-50 && (funcDecrease < 0.9 *expectedFuncDecrease))) {
2560  if (leg == 2 && alpha == 1.0) {
2561  ll = trustRegionLength();
2562  if (ll < 2.0 * m_normDeltaSoln_Newton) {
2563  trustDelta_ *= 2.0;
2564  NextTrustFactor_ *= 2.0;
2566  ll = trustRegionLength();
2567  if (m_print_flag >= 4) {
2568  printf("\t\t decideStep: Trust region further increased from %g to %g next step due to good linear model behavior\n",
2569  ll*0.5, ll);
2570  }
2571  }
2572  } else {
2573  ll = trustRegionLength();
2574  trustDelta_ *= 2.0;
2575  NextTrustFactor_ *= 2.0;
2577  ll = trustRegionLength();
2578  if (m_print_flag >= 4) {
2579  printf("\t\t decideStep: Trust region further increased from %g to %g next step due to good linear model behavior\n",
2580  ll*0.5, ll);
2581  }
2582  }
2583  }
2584  }
2585  }
2586  }
2587  }
2588  return retn;
2589 }
2590
2591 int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_comm, doublereal* const ydot_comm,
2592  doublereal CJ, doublereal time_curr, GeneralMatrix& jac,
2593  int& num_newt_its, int& num_linear_solves,
2594  int& num_backtracks, int loglevelInput)
2595 {
2596  clockWC wc;
2597  int convRes = 0;
2598  solnType_ = SolnType;
2599  int info = 0;
2600  if (neq_ <= 0) {
2601  return 1;
2602  }
2603
2604  num_linear_solves -= m_numTotalLinearSolves;
2605  int retnDamp = 0;
2606  int retnCode = 0;
2607  bool forceNewJac = false;
2608
2609  delete jacCopyPtr_;
2611
2612  doublereal stepNorm_1;
2613  doublereal stepNorm_2;
2614  int legBest;
2615  doublereal alphaBest;
2616  bool trInit = false;
2617
2618  copy(y_comm, y_comm + neq_, m_y_n_curr.begin());
2619
2620  if (SolnType != NSOLN_TYPE_STEADY_STATE || ydot_comm) {
2621  copy(ydot_comm, ydot_comm + neq_, m_ydot_n_curr.begin());
2622  copy(ydot_comm, ydot_comm + neq_, m_ydot_trial.begin());
2623  }
2624  // Redo the solution weights every time we enter the function
2626  m_normDeltaSoln_Newton = 1.0E1;
2627  bool frst = true;
2628  num_newt_its = 0;
2629  num_backtracks = 0;
2630  int i_numTrials;
2631  m_print_flag = loglevelInput;
2632
2634  trInit = true;
2635  } else if (trustRegionInitializationMethod_ == 1) {
2636  trInit = true;
2638  } else {
2639  deltaX_trust_.assign(neq_, 1.0);
2640  trustDelta_ = 1.0;
2641  }
2642
2643  if (m_print_flag == 2 || m_print_flag == 3) {
2644  printf("\tsolve_nonlinear_problem():\n\n");
2645  if (doDogLeg_) {
2646  printf("\tWt Iter Resid NewJac log(CN)| dRdS_CDexp dRdS_CD dRdS_Newtexp dRdS_Newt |"
2647  "DS_Cauchy DS_Newton DS_Trust | legID legAlpha Fbound | CTF NTF | nTr|"
2648  "DS_Final ResidLag ResidFull\n");
2649  printf("\t---------------------------------------------------------------------------------------------------"
2650  "--------------------------------------------------------------------------------\n");
2651  } else {
2652  printf("\t Wt Iter Resid NewJac | Fbound ResidBound | DampIts Fdamp DS_Step1 DS_Step2"
2653  "ResidLag | DS_Damp DS_Newton ResidFull\n");
2654  printf("\t--------------------------------------------------------------------------------------------------"
2655  "----------------------------------\n");
2656  }
2657  }
2658
2659  while (1 > 0) {
2660
2661  CurrentTrustFactor_ = 1.0;
2662  NextTrustFactor_ = 1.0;
2663  ResidWtsReevaluated_ = false;
2664  i_numTrials = 0;
2665  /*
2666  * Increment Newton Solve counter
2667  */
2669  num_newt_its++;
2671
2672  if (m_print_flag > 3) {
2673  printf("\t");
2674  writeline('=', 119);
2675  printf("\tsolve_nonlinear_problem(): iteration %d:\n",
2676  num_newt_its);
2677  }
2678  /*
2679  * If we are far enough away from the solution, redo the solution weights and the trust vectors.
2680  */
2681  if (m_normDeltaSoln_Newton > 1.0E2 || (num_newt_its % 5) == 1) {
2683  if (trInit && (DEBUG_MODE_ENABLED || doDogLeg_)) {
2685  }
2686  }
2687
2688  /*
2689  * Set default values of Delta bounds constraints
2690  */
2691  if (!m_manualDeltaStepSet) {
2693  }
2694
2695  // Check whether the Jacobian should be re-evaluated.
2696
2697  forceNewJac = true;
2698
2699  if (forceNewJac) {
2700  if (m_print_flag > 3) {
2701  printf("\t solve_nonlinear_problem(): Getting a new Jacobian\n");
2702  }
2703  info = beuler_jac(jac, DATA_PTR(m_resid), time_curr, CJ, DATA_PTR(m_y_n_curr),
2704  DATA_PTR(m_ydot_n_curr), num_newt_its);
2705  if (info != 1) {
2706  if (m_print_flag > 0) {
2707  printf("\t solve_nonlinear_problem(): Jacobian Formation Error: %d Bailing\n", info);
2708  }
2710  goto done;
2711  }
2712  } else {
2713  if (m_print_flag > 1) {
2714  printf("\t solve_nonlinear_problem(): Solving system with old Jacobian\n");
2715  }
2716  }
2717  /*
2718  * Go get new scales
2719  */
2720  calcColumnScales();
2721
2722
2723  /*
2724  * Calculate the base residual
2725  */
2726  if (m_print_flag >= 6) {
2727  printf("\t solve_nonlinear_problem(): Calculate the base residual\n");
2728  }
2730  if (info != 1) {
2731  if (m_print_flag > 0) {
2732  printf("\t solve_nonlinear_problem(): Residual Calc ERROR %d. Bailing\n", info);
2733  }
2735  goto done;
2736  }
2737
2738  /*
2739  * Scale the matrix and the RHS, if they aren't already scaled
2740  * Figure out and store the residual scaling factors.
2741  */
2742  scaleMatrix(jac, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr), time_curr, num_newt_its);
2743
2744
2745  /*
2746  * Optional print out the initial residual
2747  */
2748  if (m_print_flag >= 6) {
2749  m_normResid_0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 10, DATA_PTR(m_y_n_curr));
2750  } else {
2751  m_normResid_0 = residErrorNorm(DATA_PTR(m_resid), "Initial norm of the residual", 0, DATA_PTR(m_y_n_curr));
2752  if (m_print_flag == 4 || m_print_flag == 5) {
2753  printf("\t solve_nonlinear_problem(): Initial Residual Norm = %13.4E\n", m_normResid_0);
2754  }
2755  }
2756
2757  if (DEBUG_MODE_ENABLED || doDogLeg_) {
2758  if (m_print_flag > 3) {
2759  printf("\t solve_nonlinear_problem(): Calculate the steepest descent direction and Cauchy Point\n");
2760  }
2762  }
2763
2764  // compute the undamped Newton step
2765  if (doAffineSolve_) {
2766  if (m_print_flag >= 4) {
2767  printf("\t solve_nonlinear_problem(): Calculate the Newton direction via an Affine solve\n");
2768  }
2770  } else {
2771  if (m_print_flag >= 4) {
2772  printf("\t solve_nonlinear_problem(): Calculate the Newton direction via a Newton solve\n");
2773  }
2775  }
2776
2777  if (info) {
2779  if (m_print_flag > 0) {
2780  printf("\t solve_nonlinear_problem(): Matrix Inversion Error: %d Bailing\n", info);
2781  }
2782  goto done;
2783  }
2785
2786  if (m_print_flag >= 6) {
2787  m_normDeltaSoln_Newton = solnErrorNorm(DATA_PTR(deltaX_Newton_), "Initial Undamped Newton Step of the iteration", 10);
2788  } else {
2789  m_normDeltaSoln_Newton = solnErrorNorm(DATA_PTR(deltaX_Newton_), "Initial Undamped Newton Step of the iteration", 0);
2790  }
2791
2792  if (m_numTotalNewtIts == 1) {
2794  if (m_print_flag > 3) {
2796  printf("\t solve_nonlinear_problem(): Initialize the trust region size as the length of the Cauchy Vector times %f\n",
2798  } else {
2799  printf("\t solve_nonlinear_problem(): Initialize the trust region size as the length of the Newton Vector times %f\n",
2801  }
2802  }
2804  trInit = true;
2805  }
2806  }
2807
2808  if (DEBUG_MODE_ENABLED && doDogLeg_ && m_print_flag >= 4) {
2809  doublereal trustD = calcTrustDistance(m_step_1);
2810  if (trustD > trustDelta_) {
2811  printf("\t\t Newton's method step size, %g trustVectorUnits, larger than trust region, %g trustVectorUnits\n",
2812  trustD, trustDelta_);
2813  printf("\t\t Newton's method step size, %g trustVectorUnits, larger than trust region, %g trustVectorUnits\n",
2814  trustD, trustDelta_);
2815  } else {
2816  printf("\t\t Newton's method step size, %g trustVectorUnits, smaller than trust region, %g trustVectorUnits\n",
2817  trustD, trustDelta_);
2818  }
2819  }
2820
2821  /*
2822  * Filter out bad directions
2823  */
2825
2826
2827
2828  if (s_print_DogLeg && m_print_flag >= 4) {
2829  printf("\t solve_nonlinear_problem(): Compare descent rates for Cauchy and Newton directions\n");
2830  descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(m_ydot_trial), i_numTrials);
2831  } else {
2832  if (doDogLeg_) {
2833  descentComparison(time_curr, DATA_PTR(m_ydot_n_curr), DATA_PTR(m_ydot_trial), i_numTrials);
2834  }
2835  }
2836
2837
2838
2839  if (doDogLeg_) {
2841  if (DEBUG_MODE_ENABLED && s_print_DogLeg && m_print_flag >= 5) {
2842  printf("\t solve_nonlinear_problem(): Compare Linear and nonlinear residuals along double dog-leg path\n");
2843  residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr), legBest, alphaBest);
2844  }
2845  if (m_print_flag >= 4) {
2846  printf("\t solve_nonlinear_problem(): Calculate damping along dog-leg path to ensure residual decrease\n");
2847  }
2848  retnDamp = dampDogLeg(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr),
2849  m_step_1, DATA_PTR(m_y_n_trial), DATA_PTR(m_ydot_trial), stepNorm_1, stepNorm_2, jac, i_numTrials);
2850  } else if (DEBUG_MODE_ENABLED && s_print_DogLeg && m_print_flag >= 5) {
2851  printf("\t solve_nonlinear_problem(): Compare Linear and nonlinear residuals along double dog-leg path\n");
2852  residualComparisonLeg(time_curr, DATA_PTR(m_ydot_n_curr), legBest, alphaBest);
2853  }
2854
2855  // Damp the Newton step
2856  /*
2857  * On return the recommended new solution and derivatives is located in:
2858  * y_new
2859  * y_dot_new
2860  * The update delta vector is located in
2861  * stp1
2862  * The estimate of the solution update norm for the next step is located in
2863  * s1
2864  */
2865  if (!doDogLeg_) {
2866  retnDamp = dampStep(time_curr, DATA_PTR(m_y_n_curr), DATA_PTR(m_ydot_n_curr),
2868  DATA_PTR(m_wksp_2), stepNorm_2, jac, frst, i_numTrials);
2869  frst = false;
2870  num_backtracks += i_numTrials;
2871  stepNorm_1 = solnErrorNorm(DATA_PTR(m_step_1));
2872  }
2873
2874
2875  /*
2876  * Impose the minimum number of Newton iterations criteria
2877  */
2878  if (num_newt_its < m_min_newt_its) {
2879  if (retnDamp > NSOLN_RETN_CONTINUE) {
2880  if (m_print_flag > 2) {
2881  printf("\t solve_nonlinear_problem(): Damped Newton successful (m=%d) but minimum Newton"
2882  "iterations not attained. Resolving ...\n", retnDamp);
2883  }
2884  retnDamp = NSOLN_RETN_CONTINUE;
2885  }
2886  }
2887
2888  /*
2889  * Impose max Newton iteration
2890  */
2891  if (num_newt_its > maxNewtIts_) {
2893  if (m_print_flag > 1) {
2894  printf("\t solve_nonlinear_problem(): Damped Newton unsuccessful (max newts exceeded) sfinal = %g\n",
2895  stepNorm_1);
2896  }
2897  }
2898
2899  /*
2900  * Do a full residual calculation with the unlagged solution components.
2901  * Then get the norm of the residual
2902  */
2904  if (info != 1) {
2905  if (m_print_flag > 0) {
2906  printf("\t solve_nonlinear_problem(): current trial step and damping led to Residual Calc "
2907  "ERROR %d. Bailing\n", info);
2908  }
2910  goto done;
2911  }
2912  if (m_print_flag >= 4) {
2913  m_normResid_full = residErrorNorm(DATA_PTR(m_resid), " Resulting full residual norm", 10, DATA_PTR(m_y_n_trial));
2914  if (fabs(m_normResid_full - m_normResid_1) > 1.0E-3 * (m_normResid_1 + m_normResid_full + 1.0E-4)) {
2915  if (m_print_flag >= 4) {
2916  printf("\t solve_nonlinear_problem(): Full residual norm changed from %g to %g due to "
2917  "lagging of components\n", m_normResid_1, m_normResid_full);
2918  }
2919  }
2920  } else {
2922  }
2923
2924  /*
2925  * Check the convergence criteria
2926  */
2927  convRes = 0;
2928  if (retnDamp > NSOLN_RETN_CONTINUE) {
2929  convRes = convergenceCheck(retnDamp, stepNorm_1);
2930  }
2931
2932
2933
2934
2935  bool m_filterIntermediate = false;
2936  if (m_filterIntermediate) {
2937  if (retnDamp == NSOLN_RETN_CONTINUE) {
2939  }
2940  }
2941
2942  // Exchange new for curr solutions
2943  if (retnDamp >= NSOLN_RETN_CONTINUE) {
2945
2948  }
2949  }
2950
2951  if (m_print_flag == 2 || m_print_flag == 3) {
2952  if (ResidWtsReevaluated_) {
2953  printf("\t*");
2954  } else {
2955  printf("\t ");
2956  }
2957  printf(" %3d %11.3E", num_newt_its, m_normResid_0);
2958  bool m_jacAge = false;
2959  if (!m_jacAge) {
2960  printf(" Y ");
2961  } else {
2962  printf(" N ");
2963  }
2964  if (doDogLeg_) {
2965  printf("%5.1F |", log10(m_conditionNumber));
2966  // printf("\t Iter Resid NewJac | DS_Cauchy DS_Newton DS_Trust | legID legAlpha Fbound | | DS_F ResidFinal \n");
2967  printf("%10.3E %10.3E %10.3E %10.3E|", ResidDecreaseSDExp_, ResidDecreaseSD_,
2969  printf("%10.3E %10.3E %10.3E|", m_normDeltaSoln_CP , m_normDeltaSoln_Newton, norm_deltaX_trust_ * trustDelta_);
2970  printf("%2d %10.2E %10.2E", dogLegID_ , dogLegAlpha_, m_dampBound);
2971  printf("| %3.2f %3.2f |", CurrentTrustFactor_, NextTrustFactor_);
2972  printf(" %2d ", i_numTrials);
2973  printf("| %10.3E %10.3E %10.3E", stepNorm_1, m_normResid_1, m_normResid_full);
2974  } else {
2975  printf(" |");
2976  printf("%10.2E %10.3E |", m_dampBound, m_normResid_Bound);
2977  printf("%2d %10.2E %10.3E %10.3E %10.3E", i_numTrials + 1, m_dampRes,
2978  stepNorm_1 / (m_dampRes * m_dampBound), stepNorm_2, m_normResid_1);
2979  printf("| %10.3E %10.3E %10.3E", stepNorm_1, m_normDeltaSoln_Newton, m_normResid_full);
2980  }
2981  printf("\n");
2982
2983  }
2984  if (m_print_flag >= 4) {
2985  if (doDogLeg_) {
2986  if (convRes > 0) {
2987  printf("\t solve_nonlinear_problem(): Problem Converged, stepNorm = %11.3E, reduction of res from %11.3E to %11.3E\n",
2988  stepNorm_1, m_normResid_0, m_normResid_full);
2989  printf("\t");
2990  writeline('=', 119);
2991  } else {
2992  printf("\t solve_nonlinear_problem(): Successfull step taken with stepNorm = %11.3E, reduction of res from %11.3E to %11.3E\n",
2993  stepNorm_1, m_normResid_0, m_normResid_full);
2994  }
2995  } else {
2996  if (convRes > 0) {
2997  printf("\t solve_nonlinear_problem(): Damped Newton iteration successful, nonlin "
2998  "converged, final estimate of the next solution update norm = %-12.4E\n", stepNorm_2);
2999  printf("\t");
3000  writeline('=', 119);
3001  } else if (retnDamp >= NSOLN_RETN_CONTINUE) {
3002  printf("\t solve_nonlinear_problem(): Damped Newton iteration successful, "
3003  "estimate of the next solution update norm = %-12.4E\n", stepNorm_2);
3004  } else {
3005  printf("\t solve_nonlinear_problem(): Damped Newton unsuccessful, final estimate "
3006  "of the next solution update norm = %-12.4E\n", stepNorm_2);
3007  }
3008  }
3009  }
3010
3011  /*
3012  * Do a full residual calculation with the unlagged solution components calling ShowSolution to perhaps print out the solution.
3013  */
3014  if (m_print_flag >= 4) {
3015  if (convRes > 0 || m_print_flag >= 6) {
3017  if (info != 1) {
3018  if (m_print_flag > 0) {
3019  printf("\t solve_nonlinear_problem(): Final ShowSolution residual eval returned an error! "
3020  "ERROR %d. Bailing\n", info);
3021  }
3023  goto done;
3024  }
3025  }
3026  }
3027
3028  // convergence
3029  if (convRes) {
3030  goto done;
3031  }
3032
3033  // If dampStep fails, first try a new Jacobian if an old
3034  // one was being used. If it was a new Jacobian, then
3035  // return -1 to signify failure.
3036  else if (retnDamp < NSOLN_RETN_CONTINUE) {
3037  goto done;
3038  }
3039  }
3040
3041 done:
3042
3043
3044  if (m_print_flag == 2 || m_print_flag == 3) {
3045  if (convRes > 0) {
3046  if (doDogLeg_) {
3047  if (convRes == 3) {
3048  printf("\t | | "
3049  " | | converged = 3 |(%11.3E) \n", stepNorm_2);
3050  } else {
3051  printf("\t | | "
3052  " | | converged = %1d | %10.3E %10.3E\n", convRes,
3053  stepNorm_2, m_normResidTrial);
3054  }
3055  printf("\t-----------------------------------------------------------------------------------------------------"
3056  "------------------------------------------------------------------------------\n");
3057  } else {
3058  if (convRes == 3) {
3059  printf("\t | "
3060  " | converged = 3 | (%11.3E) \n", stepNorm_2);
3061  } else {
3062  printf("\t | "
3063  " | converged = %1d | %10.3E %10.3E\n", convRes,
3064  stepNorm_2, m_normResidTrial);
3065  }
3066  printf("\t------------------------------------------------------------------------------------"
3067  "-----------------------------------------------\n");
3068  }
3069  }
3070
3071
3072
3073  }
3074
3075  copy(m_y_n_curr.begin(), m_y_n_curr.end(), y_comm);
3077  copy(m_ydot_n_curr.begin(), m_ydot_n_curr.end(), ydot_comm);
3078  }
3079
3080  num_linear_solves += m_numTotalLinearSolves;
3081
3082  doublereal time_elapsed = wc.secondsWC();
3083  if (m_print_flag > 1) {
3084  if (retnDamp > 0) {
3086  printf("\tNonlinear problem solved successfully in %d its\n",
3087  num_newt_its);
3088  } else {
3089  printf("\tNonlinear problem solved successfully in %d its, time elapsed = %g sec\n",
3090  num_newt_its, time_elapsed);
3091  }
3092  } else {
3093  printf("\tNonlinear problem failed to solve after %d its\n", num_newt_its);
3094  }
3095  }
3096  retnCode = retnDamp;
3097  if (retnDamp > 0) {
3098  retnCode = NSOLN_RETN_SUCCESS;
3099  }
3100
3101
3102  return retnCode;
3103 }
3104
3105 void NonlinearSolver::setPreviousTimeStep(const std::vector<doublereal>& y_nm1,
3106  const std::vector<doublereal>& ydot_nm1)
3107 {
3108  m_y_nm1 = y_nm1;
3109  m_ydot_nm1 = ydot_nm1;
3110 }
3111
3112 void NonlinearSolver::print_solnDelta_norm_contrib(const doublereal* const step_1,
3113  const char* const stepNorm_1,
3114  const doublereal* const step_2,
3115  const char* const stepNorm_2,
3116  const char* const title,
3117  const doublereal* const y_n_curr,
3118  const doublereal* const y_n_1,
3119  doublereal damp,
3120  size_t num_entries)
3121 {
3122  bool used;
3123  doublereal dmax0, dmax1, error, rel_norm;
3124  printf("\t\t%s currentDamp = %g\n", title, damp);
3125  printf("\t\t I ysolnOld %13s ysolnNewRaw | ysolnNewTrial "
3126  "%10s ysolnNewTrialRaw | solnWeight wtDelSoln wtDelSolnTrial\n", stepNorm_1, stepNorm_2);
3127  std::vector<size_t> imax(num_entries, npos);
3128  printf("\t\t ");
3129  writeline('-', 125);
3130  for (size_t jnum = 0; jnum < num_entries; jnum++) {
3131  dmax1 = -1.0;
3132  for (size_t i = 0; i < neq_; i++) {
3133  used = false;
3134  for (size_t j = 0; j < jnum; j++) {
3135  if (imax[j] == i) {
3136  used = true;
3137  }
3138  }
3139  if (!used) {
3140  error = step_1[i] / m_ewt[i];
3141  rel_norm = sqrt(error * error);
3142  error = step_2[i] / m_ewt[i];
3143  rel_norm += sqrt(error * error);
3144  if (rel_norm > dmax1) {
3145  imax[jnum] = i;
3146  dmax1 = rel_norm;
3147  }
3148  }
3149  }
3150  if (imax[jnum] != npos) {
3151  size_t i = imax[jnum];
3152  error = step_1[i] / m_ewt[i];
3153  dmax0 = sqrt(error * error);
3154  error = step_2[i] / m_ewt[i];
3155  dmax1 = sqrt(error * error);
3156  printf("\t\t %4s %12.4e %12.4e %12.4e | %12.4e %12.4e %12.4e |%12.4e %12.4e %12.4e\n",
3157  int2str(i).c_str(), y_n_curr[i], step_1[i], y_n_curr[i] + step_1[i], y_n_1[i],
3158  step_2[i], y_n_1[i]+ step_2[i], m_ewt[i], dmax0, dmax1);
3159  }
3160  }
3161  printf("\t\t ");
3162  writeline('-', 125);
3163 }
3164
3165 //! This routine subtracts two numbers for one another
3166 /*!
3167  * This routine subtracts 2 numbers. If the difference is less
3168  * than 1.0E-14 times the magnitude of the smallest number, then diff returns an exact zero.
3169  * It also returns an exact zero if the difference is less than
3170  * 1.0E-300.
3171  *
3172  * returns: a - b
3173  *
3174  * This routine is used in numerical differencing schemes in order
3175  * to avoid roundoff errors resulting in creating Jacobian terms.
3176  * Note: This is a slow routine. However, Jacobian errors may cause
3177  * loss of convergence. Therefore, in practice this routine has proved cost-effective.
3178  *
3179  * @param a Value of a
3180  * @param b value of b
3181  *
3182  * @return returns the difference between a and b
3183  */
3184 static inline doublereal subtractRD(doublereal a, doublereal b)
3185 {
3186  doublereal diff = a - b;
3187  doublereal d = std::min(fabs(a), fabs(b));
3188  d *= 1.0E-14;
3190  if (ad < 1.0E-300) {
3191  diff = 0.0;
3192  }
3193  if (ad < d) {
3194  diff = 0.0;
3195  }
3196  return diff;
3197 }
3198
3199 int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f,
3200  doublereal time_curr, doublereal CJ,
3201  doublereal* const y, doublereal* const ydot,
3202  int num_newt_its)
3203 {
3204  double* col_j;
3205  int info;
3206  doublereal ysave, dy;
3207  doublereal ydotsave = 0;
3208  int retn = 1;
3209
3210  /*
3211  * Clear the factor flag
3212  */
3213  J.clearFactorFlag();
3215  /********************************************************************
3216  * Call the function to get a Jacobian.
3217  */
3218  info = m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f);
3219  m_nJacEval++;
3220  m_nfe++;
3221  if (info != 1) {
3222  return info;
3223  }
3224  } else {
3225  if (J.matrixType_ == 0) {
3226  /*******************************************************************
3227  * Generic algorithm to calculate a numerical Jacobian
3228  */
3229  /*
3230  * Calculate the current value of the RHS given the
3231  * current conditions.
3232  */
3233
3234  info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
3235  m_nfe++;
3236  if (info != 1) {
3237  return info;
3238  }
3239  m_nJacEval++;
3240  if (DEBUG_MODE_ENABLED) {
3241  for (size_t ii = 0; ii < neq_; ii++) {
3242  checkFinite(f[ii]);
3243  }
3244  }
3245
3246  /*
3247  * Malloc a vector and call the function object to return a set of
3248  * deltaY's that are appropriate for calculating the numerical
3249  * derivative.
3250  */
3251  vector_fp dyVector(neq_);
3252  retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, &dyVector[0], DATA_PTR(m_ewt));
3253  if (s_print_NumJac) {
3254  if (m_print_flag >= 7) {
3255  if (retn != 1) {
3256  printf("\t\t beuler_jac ERROR! calcDeltaSolnVariables() returned an error flag\n");
3257  printf("\t\t We will bail from the nonlinear solver after calculating the Jacobian");
3258  }
3259  if (neq_ < 20) {
3260  printf("\t\tUnk m_ewt y dyVector ResN\n");
3261  for (size_t iii = 0; iii < neq_; iii++) {
3262  printf("\t\t %4s %16.8e %16.8e %16.8e %16.8e \n",
3263  int2str(iii).c_str(), m_ewt[iii], y[iii], dyVector[iii], f[iii]);
3264  }
3265  }
3266  }
3267  }
3268
3269  /*
3270  * Loop over the variables, formulating a numerical derivative
3271  * of the dense matrix.
3272  * For the delta in the variable, we will use a variety of approaches
3273  * The original approach was to use the error tolerance amount.
3274  * This may not be the best approach, as it could be overly large in
3275  * some instances and overly small in others.
3276  * We will first protect from being overly small, by using the usual
3277  * sqrt of machine precision approach, i.e., 1.0E-7,
3278  * to bound the lower limit of the delta.
3279  */
3280  for (size_t j = 0; j < neq_; j++) {
3281
3282
3283  /*
3284  * Get a pointer into the column of the matrix
3285  */
3286
3287
3288  col_j = (doublereal*) J.ptrColumn(j);
3289  ysave = y[j];
3290  dy = dyVector[j];
3291
3292  y[j] = ysave + dy;
3293  dy = y[j] - ysave;
3295  ydotsave = ydot[j];
3296  ydot[j] += dy * CJ;
3297  }
3298  /*
3299  * Call the function
3300  */
3301
3302
3303  info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp),
3304  JacDelta_ResidEval, static_cast<int>(j), dy);
3305  m_nfe++;
3306
3307  if (DEBUG_MODE_ENABLED) {
3308  if (fabs(dy) < 1.0E-300) {
3309  throw CanteraError("NonlinearSolver::beuler_jac", "dy is equal to zero");
3310  }
3311  for (size_t ii = 0; ii < neq_; ii++) {
3312  checkFinite(m_wksp[ii]);
3313  }
3314  }
3315
3316  if (info != 1) {
3317  return info;
3318  }
3319
3320  doublereal diff;
3321  for (size_t i = 0; i < neq_; i++) {
3322  diff = subtractRD(m_wksp[i], f[i]);
3323  col_j[i] = diff / dy;
3324  }
3325  y[j] = ysave;
3327  ydot[j] = ydotsave;
3328  }
3329
3330  }
3331  } else if (J.matrixType_ == 1) {
3332  int ku, kl;
3333  size_t ivec[2];
3334  size_t n = J.nRowsAndStruct(ivec);
3335  kl = static_cast<int>(ivec[0]);
3336  ku = static_cast<int>(ivec[1]);
3337  if (n != neq_) {
3338  printf("we have probs\n");
3339  exit(-1);
3340  }
3341
3342  // --------------------------------- BANDED MATRIX BRAIN DEAD ---------------------------------------------------
3343  info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
3344  m_nfe++;
3345  if (info != 1) {
3346  return info;
3347  }
3348  m_nJacEval++;
3349
3350
3351  vector_fp dyVector(neq_);
3352  retn = m_func->calcDeltaSolnVariables(time_curr, y, ydot, &dyVector[0], DATA_PTR(m_ewt));
3353  if (s_print_NumJac) {
3354  if (m_print_flag >= 7) {
3355  if (retn != 1) {
3356  printf("\t\t beuler_jac ERROR! calcDeltaSolnVariables() returned an error flag\n");
3357  printf("\t\t We will bail from the nonlinear solver after calculating the Jacobian");
3358  }
3359  if (neq_ < 20) {
3360  printf("\t\tUnk m_ewt y dyVector ResN\n");
3361  for (size_t iii = 0; iii < neq_; iii++) {
3362  printf("\t\t %4s %16.8e %16.8e %16.8e %16.8e \n",
3363  int2str(iii).c_str(), m_ewt[iii], y[iii], dyVector[iii], f[iii]);
3364  }
3365  }
3366  }
3367  }
3368
3369
3370  for (size_t j = 0; j < neq_; j++) {
3371
3372
3373  col_j = (doublereal*) J.ptrColumn(j);
3374  ysave = y[j];
3375  dy = dyVector[j];
3376
3377
3378  y[j] = ysave + dy;
3379  dy = y[j] - ysave;
3381  ydotsave = ydot[j];
3382  ydot[j] += dy * CJ;
3383  }
3384
3385  info = m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_wksp), JacDelta_ResidEval, static_cast<int>(j), dy);
3386  m_nfe++;
3387  if (DEBUG_MODE_ENABLED) {
3388  if (fabs(dy) < 1.0E-300) {
3389  throw CanteraError("NonlinearSolver::beuler_jac", "dy is equal to zero");
3390  }
3391  for (size_t ii = 0; ii < neq_; ii++) {
3392  checkFinite(m_wksp[ii]);
3393  }
3394  }
3395  if (info != 1) {
3396  return info;
3397  }
3398
3399  doublereal diff;
3400
3401
3402
3403  for (int i = (int) j - ku; i <= (int) j + kl; i++) {
3404  if (i >= 0 && i < (int) neq_) {
3405  diff = subtractRD(m_wksp[i], f[i]);
3406  col_j[kl + ku + i - j] = diff / dy;
3407  }
3408  }
3409  y[j] = ysave;
3411  ydot[j] = ydotsave;
3412  }
3413
3414  }
3415
3416  double vSmall;
3417  size_t ismall = J.checkRows(vSmall);
3418  if (vSmall < 1.0E-100) {
3419  printf("WE have a zero row, %s\n", int2str(ismall).c_str());
3420  exit(-1);
3421  }
3422  ismall = J.checkColumns(vSmall);
3423  if (vSmall < 1.0E-100) {
3424  printf("WE have a zero column, %s\n", int2str(ismall).c_str());
3425  exit(-1);
3426  }
3427
3428  // ---------------------BANDED MATRIX BRAIN DEAD -----------------------
3429  }
3430  }
3431
3432  if (m_print_flag >= 7 && s_print_NumJac) {
3433  if (neq_ < 30) {
3434  printf("\t\tCurrent Matrix and Residual:\n");
3435  printf("\t\t I,J | ");
3436  for (size_t j = 0; j < neq_; j++) {
3437  printf(" %5s ", int2str(j).c_str());
3438  }
3439  printf("| Residual \n");
3440  printf("\t\t --");
3441  for (size_t j = 0; j < neq_; j++) {
3442  printf("------------");
3443  }
3444  printf("| -----------\n");
3445
3446
3447  for (size_t i = 0; i < neq_; i++) {
3448  printf("\t\t %4s |", int2str(i).c_str());
3449  for (size_t j = 0; j < neq_; j++) {
3450  printf(" % 11.4E", J(i,j));
3451  }
3452  printf(" | % 11.4E\n", f[i]);
3453  }
3454
3455  printf("\t\t --");
3456  for (size_t j = 0; j < neq_; j++) {
3457  printf("------------");
3458  }
3459  printf("--------------\n");
3460  }
3461  }
3462  /*
3463  * Make a copy of the data. Note, this Jacobian copy occurs before any matrix scaling operations.
3464  * It's the raw matrix producted by this routine.
3465  */
3466  *jacCopyPtr_ = J;
3467
3468  return retn;
3469 }
3470
3471 void NonlinearSolver::calc_ydot(const int order,
3472  const doublereal* const y_curr,
3473  doublereal* const ydot_curr) const
3474 {
3475  if (!ydot_curr) {
3476  return;
3477  }
3478  doublereal c1;
3479  switch (order) {
3480  case 0:
3481  case 1: /* First order forward Euler/backward Euler */
3482  c1 = 1.0 / delta_t_n;
3483  for (size_t i = 0; i < neq_; i++) {
3484  ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]);
3485  }
3486  return;
3487  case 2: /* Second order Adams-Bashforth / Trapezoidal Rule */
3488  c1 = 2.0 / delta_t_n;
3489  for (size_t i = 0; i < neq_; i++) {
3490  ydot_curr[i] = c1 * (y_curr[i] - m_y_nm1[i]) - m_ydot_nm1[i];
3491  }
3492
3493  return;
3494  default:
3495  throw CanteraError("calc_ydot()", "Case not covered");
3496  }
3497 }
3498
3499 doublereal NonlinearSolver::filterNewStep(const doublereal timeCurrent,
3500  const doublereal* const ybase, doublereal* const step0)
3501 {
3502  return m_func->filterNewStep(timeCurrent, ybase, step0);
3503 }
3504
3505 doublereal NonlinearSolver::filterNewSolution(const doublereal timeCurrent,
3506  doublereal* const y_current, doublereal* const ydot_current)
3507 {
3508  return m_func->filterSolnPrediction(timeCurrent, y_current);
3509 }
3510
3512 {
3513  ResidWtsReevaluated_ = true;
3514  if (checkUserResidualTols_ == 1) {
3515  for (size_t i = 0; i < neq_; i++) {
3517  if (DEBUG_MODE_ENABLED) {
3518  checkFinite(m_residWts[i]);
3519  }
3520  }
3521  } else {
3522  doublereal sum = 0.0;
3523  for (size_t i = 0; i < neq_; i++) {
3524  m_residWts[i] = m_rowWtScales[i];
3525  if (DEBUG_MODE_ENABLED) {
3526  checkFinite(m_residWts[i]);
3527  }
3528  sum += m_residWts[i];
3529  }
3530  sum /= neq_;
3531  for (size_t i = 0; i < neq_; i++) {
3533  }
3534  if (checkUserResidualTols_ == 2) {
3535  for (size_t i = 0; i < neq_; i++) {
3536  double uR = userResidAtol_[i] + userResidRtol_ * m_rowWtScales[i] / rtol_;
3537  m_residWts[i] = std::min(m_residWts[i], uR);
3538  }
3539  }
3540  }
3541 }
3542
3543 void NonlinearSolver::getResidWts(doublereal* const residWts) const
3544 {
3545  for (size_t i = 0; i < neq_; i++) {
3546  residWts[i] = (m_residWts)[i];
3547  }
3548 }
3549
3550 int NonlinearSolver::convergenceCheck(int dampCode, doublereal s1)
3551 {
3552  int retn = 0;
3553  if (m_dampBound < 0.9999) {
3554  return retn;
3555  }
3556  if (m_dampRes < 0.9999) {
3557  return retn;
3558  }
3559  if (dampCode <= 0) {
3560  return retn;
3561  }
3562  if (dampCode == 3) {
3563  if (s1 < 1.0E-2) {
3564  if (m_normResid_full < 1.0E-1) {
3565  return 3;
3566  }
3567  }
3568  if (s1 < 0.8) {
3569  if (m_normDeltaSoln_Newton < 1.0) {
3570  if (m_normResid_full < 1.0E-1) {
3571  return 2;
3572  }
3573  }
3574  }
3575  }
3576  if (dampCode == 4) {
3577  if (s1 < 1.0E-2) {
3578  if (m_normResid_full < 1.0E-1) {
3579  return 3;
3580  }
3581  }
3582  }
3583
3584  if (s1 < 0.8) {
3585  if (m_normDeltaSoln_Newton < 1.0) {
3586  if (m_normResid_full < 1.0E-1) {
3587  return 2;
3588  }
3589  }
3590  }
3591  if (dampCode == 1 || dampCode == 2) {
3592  if (s1 < 1.0) {
3593  if (m_normResid_full < 1.0E-1) {
3594  return 1;
3595  }
3596  }
3597  }
3598  return retn;
3599 }
3600
3601 void NonlinearSolver::setAtol(const doublereal* const atol)
3602 {
3603  for (size_t i = 0; i < neq_; i++) {
3604  if (atol[i] <= 0.0) {
3605  throw CanteraError("NonlinearSolver::setAtol()",
3606  "Atol is less than or equal to zero");
3607  }
3608  atolk_[i]= atol[i];
3609  }
3610 }
3611
3612 void NonlinearSolver::setRtol(const doublereal rtol)
3613 {
3614  if (rtol <= 0.0) {
3615  throw CanteraError("NonlinearSolver::setRtol()",
3616  "Rtol is <= zero");
3617  }
3618  rtol_ = rtol;
3619 }
3620
3621 void NonlinearSolver::setResidualTols(double residRtol, double* residATol, int residNormHandling)
3622 {
3623  if (residNormHandling < 0 || residNormHandling > 2) {
3624  throw CanteraError("NonlinearSolver::setResidualTols()",
3625  "Unknown int for residNormHandling");
3626  }
3627  checkUserResidualTols_ = residNormHandling;
3628  userResidRtol_ = residRtol;
3629  if (residATol) {
3630  userResidAtol_.resize(neq_);
3631  for (size_t i = 0; i < neq_; i++) {
3632  userResidAtol_[i] = residATol[i];
3633  }
3634  } else {
3635  if (residNormHandling ==1 || residNormHandling == 2) {
3636  throw CanteraError("NonlinearSolver::setResidualTols()",
3637  "Must set residATol vector");
3638  }
3639  }
3640 }
3641
3643 {
3644  m_print_flag = printLvl;
3645 }
3646
3647 }
static bool s_alwaysAssumeNewtonGood
This toggle turns off the use of the Hessian when it is warranted by the condition number...
virtual doublereal * ptrColumn(size_t j)=0
Return a pointer to the top of column j, columns are assumed to be contiguous in memory.
int matrixType_
Matrix type.
int trustRegionInitializationMethod_
Method for handling the trust region initialization.
void computeResidWts()
Compute the Residual Weights.
std::vector< doublereal > m_y_low_bounds
Lower bounds vector for each species.
doublereal ResidDecreaseSDExp_
Expected DResid_dS for the steepest descent path - output variable.
std::string int2str(const int n, const std::string &fmt)
Convert an int to a string using a format converter.
Definition: stringUtils.cpp:39
#define NSOLN_RETN_SUCCESS
The nonlinear solve is successful.
int beuler_jac(GeneralMatrix &J, doublereal *const f, doublereal time_curr, doublereal CJ, doublereal *const y, doublereal *const ydot, int num_newt_its)
Function called to evaluate the Jacobian matrix and the current residual vector at the current time s...
int m_matrixConditioning
Boolean indicating matrix conditioning.
Base residual calculation for the showSolution routine.
Definition: ResidJacEval.h:35
void checkFinite(const double tmp)
Check to see that a number is finite (not NaN, +Inf or -Inf)
Definition: checkFinite.cpp:43
doublereal Nuu_
Relative distance down the Newton step that the second dogleg starts.
int m_rowScaling
int indicating whether row scaling is turned on (1) or not (0)
virtual int nEquations() const
Return the number of equations in the equation system.
int convergenceCheck(int dampCode, doublereal s1)
Check to see if the nonlinear problem has converged.
doublereal solnErrorNorm(const doublereal *const delta_y, const char *title=0, int printLargest=0, const doublereal dampFactor=1.0) const
L2 norm of the delta of the solution vector.
std::vector< doublereal > atolk_
absolute tolerance in the solution unknown
doublereal boundStep(const doublereal *const y, const doublereal *const step0)
Bound the step.
doublereal m_normResidTrial
Norm of the residual for a trial calculation which may or may not be used.
doublereal lambdaStar_
Value of lambdaStar_ which is used to calculate the Cauchy point.
doublereal ResidDecreaseNewtExp_
Expected DResid_dS for the Newton path - output variable.
std::vector< doublereal > m_colScales
Vector of column scaling factors.
doublereal normTrust_Newton_
Norm of the Newton Step wrt trust region.
std::vector< doublereal > userResidAtol_
absolute tolerance in the unscaled solution unknowns
Declarations for a simple class that implements an Ansi C wall clock timer (see Cantera::clockWC).
virtual GeneralMatrix * duplMyselfAsGeneralMatrix() const =0
Duplicator member function.
doublereal userResidRtol_
absolute tolerance in the unscaled solution unknowns
doublereal m_normResid_Bound
Norm of the residual after it has been bounded.
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:165
virtual doublereal oneNorm() const =0
Calculate the one norm of the matrix.
virtual void clearFactorFlag()
clear the factored flag
virtual void setPreviousTimeStep(const std::vector< doublereal > &y_nm1, const std::vector< doublereal > &ydot_nm1)
Set the values for the previous time step.
#define NSOLN_RETN_MATRIXINVERSIONERROR
The nonlinear problem's Jacobian is singular.
int doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal *const y_curr, const doublereal *const ydot_curr, const ResidEval_Type_Enum evalType=Base_ResidEval) const
Compute the current residual.
virtual int evalResidNJ(const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot, doublereal *const resid, const ResidEval_Type_Enum evalType=Base_ResidEval, const int id_x=-1, const doublereal delta_x=0.0)
Evaluate the residual function.
int dogLegID_
Current leg.
virtual doublereal rcond(doublereal a1norm)=0
Returns an estimate of the inverse of the condition number for the matrix.
virtual int evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj, const doublereal *const y, const doublereal *const ydot, GeneralMatrix &J, doublereal *const resid)
Calculate an analytical Jacobian and the residual at the current time and values. ...
std::vector< doublereal > m_ydot_nm1
Vector containing the solution derivative at the previous time step.
void setupDoubleDogleg()
Setup the parameters for the double dog leg.
#define NSOLN_RETN_CONTINUE
Problem isn't solved yet.
virtual vector_fp::iterator begin()=0
Return an iterator pointing to the first element.
void calcSolnToResNormVector()
Calculate the scaling factor for translating residual norms into solution norms.
The nonlinear problem is part of a steady state calculation.
void warn_deprecated(const std::string &method, const std::string &extra)
Print a warning indicating that method is deprecated.
Definition: global.cpp:78
std::vector< doublereal > deltaX_trust_
Vector of trust region values.
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
int dampStep(const doublereal time_curr, const doublereal *const y_n_curr, const doublereal *const ydot_n_curr, doublereal *const step_1, doublereal *const y_n_1, doublereal *const ydot_n_1, doublereal *step_2, doublereal &stepNorm_2, GeneralMatrix &jac, bool writetitle, int &num_backtracks)
Find a damping coefficient through a look-ahead mechanism.
Generic matrix.
Definition: GeneralMatrix.h:24
Base residual calculation for the Jacobian calculation.
Definition: ResidJacEval.h:28
#define NSOLN_RETN_FAIL_DAMPSTEP
The nonlinear problem didn't solve the problem.
virtual int factorQR()
Factors the A matrix using the QR algorithm, overwriting A.
Definition: GeneralMatrix.h:81
The class provides the wall clock timer in seconds.
Definition: clockWC.h:51
int m_print_flag
Determines the level of printing for each time step.
int m_nfe
Counter for the total number of function evaluations.
virtual int factorAlgorithm() const =0
Return the factor algorithm used.
NonlinearSolver & operator=(const NonlinearSolver &right)
Assignment operator.
doublereal m_normDeltaSoln_CP
Norm of the distance to the Cauchy point using the solution norm.
std::vector< doublereal > deltaX_Newton_
Newton Step - This is the Newton step determined from the straight Jacobian.
doublereal residErrorNorm(const doublereal *const resid, const char *title=0, const int printLargest=0, const doublereal *const y=0) const
L2 norm of the residual of the equation system.
doublereal NextTrustFactor_
Factor indicating how much trust region has been changed next iteration - output variable.
doublereal filterNewSolution(const doublereal timeCurrent, doublereal *const y_current, doublereal *const ydot_current)
Apply a filter to the solution.
#define NSOLN_JAC_ANAL
The Jacobian is calculated from an analytical function.
virtual int factor()=0
Factors the A matrix, overwriting A.
#define NSOLN_RETN_MAXIMUMITERATIONSEXCEEDED
The nonlinear problem's max number of iterations has been exceeded.
double subtractRD(double a, double b)
This routine subtracts two numbers for one another.
Definition: BEulerInt.cpp:468
std::vector< doublereal > m_rowScales
Weights for normalizing the values of the residuals.
virtual ResidJacEval * duplMyselfAsResidJacEval() const
Duplication routine for objects derived from residJacEval.
std::vector< doublereal > m_ewt
Soln error weights.
ResidEval_Type_Enum
Differentiates the type of residual evaluations according to functionality.
Definition: ResidJacEval.h:24
doublereal dogLegAlpha_
Current Alpha param along the leg.
int m_min_newt_its
Minimum number of Newton iterations to use.
void setDefaultDeltaBoundsMagnitudes()
Set default deulta bounds amounts.
std::vector< doublereal > Jd_
Jacobian times the steepest descent direction in the normalized coordinates.
doublereal m_normDeltaSoln_Newton
Norm of the solution update created by the iteration in its raw, undamped form, using the solution no...
static bool s_print_NumJac
Turn on or off printing of the Jacobian.
int checkUserResidualTols_
Check the residual tolerances explicitly against user input.
NonlinearSolver(ResidJacEval *func)
Default constructor.
Wrappers for the function evaluators for Nonlinear solvers and Time steppers.
Definition: ResidJacEval.h:51
std::vector< doublereal > deltaX_CP_
Steepest descent direction. This is also the distance to the Cauchy Point.
doublereal m_dampBound
Damping factor imposed by hard bounds and by delta bounds.
doublereal trustDelta_
int m_colScaling
The type of column scaling used in the matrix inversion of the problem.
virtual int calcDeltaSolnVariables(const doublereal t, const doublereal *const y, const doublereal *const ydot, doublereal *const delta_y, const doublereal *const solnWeights=0)
Return a vector of delta y's for calculation of the numerical Jacobian.
void createSolnWeights(const doublereal *const y)
Create solution weights for convergence criteria.
std::vector< doublereal > m_deltaStepMaximum
Value of the delta step magnitudes.
int m_numLocalLinearSolves
Number of local linear solves done during the current iteration.
virtual void zero()=0
Zero the matrix elements.
void print_solnDelta_norm_contrib(const doublereal *const step_1, const char *const stepNorm_1, const doublereal *const step_2, const char *const stepNorm_2, const char *const title, const doublereal *const y_n_curr, const doublereal *const y_n_1, doublereal damp, size_t num_entries)
Print solution norm contribution.
int doNewtonSolve(const doublereal time_curr, const doublereal *const y_curr, const doublereal *const ydot_curr, doublereal *const delta_y, GeneralMatrix &jac)
Compute the undamped Newton step.
doublereal calcTrustDistance(std::vector< doublereal > const &deltaX) const
Calculate the trust distance of a step in the solution variables.
Class that calculates the solution to a nonlinear, dense, set of equations (see Numerical Utilities w...
std::vector< doublereal > m_residWts
Vector of residual weights.
doublereal trustRegionLength() const
Calculate the length of the current trust region in terms of the solution error norm.
void calcColumnScales()
Set the column scaling vector at the current time.
bool m_resid_scaled
Boolean indicating whether we should scale the residual.
doublereal ResidDecreaseNewt_
Actual DResid_dS for the Newton path - output variable.
void setAtol(const doublereal *const atol)
Set the absolute tolerances for the solution variables.
int dampDogLeg(const doublereal time_curr, const doublereal *y_n_curr, const doublereal *ydot_n_curr, std::vector< doublereal > &step_1, doublereal *const y_n_1, doublereal *const ydot_n_1, doublereal &stepNorm_1, doublereal &stepNorm_2, GeneralMatrix &jac, int &num_backtracks)
Damp using the dog leg approach.
double secondsWC()
Returns the wall clock time in seconds since the last reset.
Definition: clockWC.cpp:35
void error(const std::string &msg)
Write an error message and quit.
Definition: global.cpp:72
void initializeTrustRegion()
Initialize the size of the trust vector.
int m_nJacEval
Number of Jacobian evaluations.
std::vector< doublereal > m_y_n_curr
Vector containing the current solution vector within the nonlinear solver.
doublereal ResidDecreaseSD_
Actual DResid_dS for the steepest descent path - output variable.
void setDeltaBoundsMagnitudes(const doublereal *const deltaBoundsMagnitudes)
Set the delta Bounds magnitudes by hand.
size_t neq_
Local copy of the number of equations.
int doAffineNewtonSolve(const doublereal *const y_curr, const doublereal *const ydot_curr, doublereal *const delta_y, GeneralMatrix &jac)
Compute the Newton step, either by direct Newton's or by solving a close problem that is represented ...
std::vector< doublereal > m_ydot_trial
Value of the solution time derivative at the new point that is to be considered.
void fillDogLegStep(int leg, doublereal alpha, std::vector< doublereal > &deltaX) const
Fill a dogleg solution step vector.
virtual size_t checkRows(doublereal &valueSmall) const =0
Check to see if we have any zero rows in the Jacobian.
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:99
Class that calculates the solution to a nonlinear system.
doublereal m_dampRes
Additional damping factor due to bounds on the residual and solution norms.
#define NSOLN_RETN_JACOBIANFORMATIONERROR
The nonlinear problem's Jacobian formation produced an error.
void setResidualTols(double residRtol, double *residAtol, int residNormHandling=2)
Set the relative and absolute tolerances for the Residual norm comparisons, if used.
std::vector< doublereal > m_y_high_bounds
Bounds vector for each species.
doublereal rtol_
value of the relative tolerance to use in solving the equation set
const int NDAMP
Number of damping steps that are carried out before the solution is deemed a failure.
Definition: BEulerInt.cpp:1289
static bool s_print_DogLeg
Turn on extra printing of dogleg information.
doublereal dist_R2_
Distance of the second leg of the dogleg in terms of the solution norm.
doublereal dist_R1_
Distance of the first leg of the dogleg in terms of the solution norm.
void scaleMatrix(GeneralMatrix &jac, doublereal *const y_comm, doublereal *const ydot_comm, doublereal time_curr, int num_newt_its)
Scale the matrix.
bool ResidWtsReevaluated_
Boolean indicating that the residual weights have been reevaluated this iteration - output variable...
doublereal RJd_norm_
Residual dot Jd norm.
void setRowScaling(bool useRowScaling)
Set the rowscaling that are used for the inversion of the matrix.
std::vector< doublereal > m_ydot_n_curr
Vector containing the time derivative of the current solution vector within the nonlinear solver (whe...
void descentComparison(doublereal time_curr, doublereal *ydot0, doublereal *ydot1, int &numTrials)
This is a utility routine that can be used to print out the rates of the initial residual decline...
doublereal filterNewStep(const doublereal timeCurrent, const doublereal *const ybase, doublereal *const step0)
Apply a filtering process to the step.
static bool s_TurnOffTiming
Turn off printing of time.
int m_numTotalLinearSolves
Total number of linear solves taken by the solver object.
virtual void calcSolnScales(const doublereal t, const doublereal *const y, const doublereal *const y_old, doublereal *const yScales)
Returns a vector of column scale factors that can be used to column scale Jacobians.
int m_manualDeltaStepSet
Boolean indicating whether a manual delta bounds has been input.
Cantera::GeneralMatrix * jacCopyPtr_
Copy of the Jacobian that doesn't get overwritten when the inverse is determined. ...
#define NSOLN_RETN_RESIDUALFORMATIONERROR
The nonlinear problem's base residual produced an error.
virtual ~NonlinearSolver()
Destructor.
virtual bool factored() const
true if the current factorization is up to date with the matrix
std::vector< doublereal > m_resid
Value of the residual for the nonlinear problem.
int m_jacFormMethod
Jacobian formation method.
int solve_nonlinear_problem(int SolnType, doublereal *const y_comm, doublereal *const ydot_comm, doublereal CJ, doublereal time_curr, GeneralMatrix &jac, int &num_newt_its, int &num_linear_solves, int &num_backtracks, int loglevelInput)
Find the solution to F(X) = 0 by damped Newton iteration.
int lambdaToLeg(const doublereal lambda, doublereal &alpha) const
Change the global lambda coordinate into the (leg,alpha) coordinate for the double dogleg...
std::vector< doublereal > m_deltaStepMinimum
Soln Delta bounds magnitudes.
std::vector< doublereal > m_y_n_trial
Vector containing the solution at the new point that is to be considered.
std::vector< doublereal > m_y_nm1
Vector containing the solution at the previous time step.
doublereal m_normResid_full
Norm of the residual at the end of the first leg of the current iteration.
doublereal m_normResid_0
Norm of the residual at the start of each nonlinear iteration.
doublereal time_n
Current system time.
virtual int solve(doublereal *b, size_t nrhs=1, size_t ldb=0)=0
Solves the Ax = b system returning x in the b spot.
int decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal *const y_n_curr, const doublereal *const ydot_n_curr, const std::vector< doublereal > &step_1, const doublereal *const y_n_1, const doublereal *const ydot_n_1, doublereal trustDeltaOld)
Decide whether the current step is acceptable and adjust the trust region size.
int m_numTotalNewtIts
Total number of Newton iterations.
virtual size_t nRowsAndStruct(size_t *const iStruct=0) const =0
Return the size and structure of the matrix.
int m_order
Order of the time step method = 1.
Cantera::GeneralMatrix * HessianPtr_
Hessian.
virtual size_t checkColumns(doublereal &valueSmall) const =0
Check to see if we have any zero columns in the Jacobian.
doublereal delta_t_n
Delta t for the current step.
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
int solnType_
Solution type.
void residualComparisonLeg(const doublereal time_curr, const doublereal *const ydot0, int &legBest, doublereal &alphaBest) const
Here we print out the residual at various points along the double dogleg, comparing against the quadr...
doublereal expectedResidLeg(int leg, doublereal alpha) const
Calculated the expected residual along the double dogleg curve.
void setRtol(const doublereal rtol)
Set the relative tolerances for the solution variables.
void setBoundsConstraints(const doublereal *const y_low_bounds, const doublereal *const y_high_bounds)
Set bounds constraints for all variables in the problem.
Delta residual calculation for the Jacbobian calculation.
Definition: ResidJacEval.h:30
void calc_ydot(const int order, const doublereal *const y_curr, doublereal *const ydot_curr) const
Internal function to calculate the time derivative of the solution at the new step.
doublereal normTrust_CP_
Norm of the Cauchy Step direction wrt trust region.
doublereal CurrentTrustFactor_
Factor indicating how much trust region has been changed this iteration - output variable.
doublereal JdJd_norm_
Dot product of the Jd_ variable defined above with itself.
Contains declarations for string manipulation functions within Cantera.
void getResidWts(doublereal *const residWts) const
Return the residual weights.
#define DATA_PTR(vec)
Creates a pointer to the start of the raw data for a vector.
Definition: ct_defs.h:36
int calcTrustIntersection(doublereal trustVal, doublereal &lambda, doublereal &alpha) const
Given a trust distance, this routine calculates the intersection of the this distance with the double...
doublereal m_conditionNumber
Condition number of the matrix.
std::vector< doublereal > m_wksp_2
Workspace of length neq_.
virtual doublereal filterSolnPrediction(const doublereal t, doublereal *const y)
Filter the solution predictions.
void setSolverScheme(int doDogLeg, int doAffineSolve)
Parameter to turn on solution solver schemes.
doublereal norm_deltaX_trust_
Current norm of the vector deltaX_trust_ in terms of the solution norm.
int doAffineSolve_
General toggle for turning on Affine solve with Hessian.
ResidJacEval * m_func
Pointer to the residual and Jacobian evaluator for the function.
doublereal deltaBoundStep(const doublereal *const y, const doublereal *const step0)
Return the factor by which the undamped Newton step 'step0' must be multiplied in order to keep the u...
doublereal trustRegionInitializationFactor_
Factor used to set the initial trust region.
doublereal m_normResid_1
Norm of the residual at the end of the first leg of the current iteration.
std::vector< doublereal > m_rowWtScales
Weights for normalizing the values of the residuals.
void setPrintLvl(int printLvl)
Set the print level from the nonlinear solver.
static bool s_doBothSolvesAndCompare
Turn on solving both the Newton and Hessian systems and comparing the results.
doublereal residNorm2Cauchy_
Expected value of the residual norm at the Cauchy point if the quadratic model were valid...
doublereal dist_R0_
Distance of the zeroeth leg of the dogleg in terms of the solution norm.
doublereal atolBase_
Base value of the absolute tolerance.
std::vector< doublereal > m_wksp
Workspace of length neq_.
int maxNewtIts_
Maximum number of Newton iterations.
doublereal m_ScaleSolnNormToResNorm
Scale factor for turning residual norms into solution norms.
virtual doublereal filterNewStep(const doublereal t, const doublereal *const ybase, doublereal *const step)
Filter the solution predictions.
Base residual calculation containing any lagged components.
Definition: ResidJacEval.h:42
virtual doublereal rcondQR()
Returns an estimate of the inverse of the condition number for the matrix.
Definition: GeneralMatrix.h:91
int doDogLeg_
General toggle for turning on dog leg damping.
doublereal dist_Total_
Distance of the sum of all legs of the doglegs in terms of the solution norm.
std::vector< doublereal > & lowBoundsConstraintVector()
Return an editable vector of the low bounds constraints.
std::vector< doublereal > m_step_1
Value of the step to be taken in the solution.
std::vector< doublereal > & highBoundsConstraintVector()
Return an editable vector of the high bounds constraints.
virtual void useFactorAlgorithm(int fAlgorithm)
Change the way the matrix is factored.
doublereal doCauchyPointSolve(GeneralMatrix &jac)
Calculate the steepest descent direction and the Cauchy Point where the quadratic formulation of the ...
void setColumnScaling(bool useColScaling, const double *const scaleFactors=0)
Set the column scaling that are used for the inversion of the matrix.