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