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