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