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