Cantera  2.0
NonlinearSolver.h
Go to the documentation of this file.
1 /**
2  * @file NonlinearSolver.h
3  * Class that calculates the solution to a nonlinear, dense, set
4  * of equations (see \ref numerics
5  * and class \link Cantera::NonlinearSolver NonlinearSolver\endlink).
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 #ifndef CT_NONLINEARSOLVER_H
16 #define CT_NONLINEARSOLVER_H
17 
18 #include "ResidJacEval.h"
19 #include "SquareMatrix.h"
20 
21 namespace Cantera
22 {
23 
24 //@{
25 /// @name Constant which determines the type of the nonlinear solve
26 /*!
27  * I think steady state is the only option I'm gunning for
28  */
29 //! The nonlinear problem is part of a pseudo time dependent calculation (NOT TESTED)
30 #define NSOLN_TYPE_PSEUDO_TIME_DEPENDENT 2
31 //! The nonlinear problem is part of a time dependent calculation
32 #define NSOLN_TYPE_TIME_DEPENDENT 1
33 //! The nonlinear problem is part of a steady state calculation
34 #define NSOLN_TYPE_STEADY_STATE 0
35 //@}
36 
37 
38 //@{
39 /// @name Constant which determines the Return int from the nonlinear solver
40 /*!
41  * This int is returned from the nonlinear solver
42  */
43 //! The nonlinear solve is successful.
44 #define NSOLN_RETN_SUCCESS 1
45 //! Problem isn't solved yet
46 #define NSOLN_RETN_CONTINUE 0
47 //! The nonlinear problem started to take too small an update step. This indicates that either the
48 //! Jacobian is bad, or a constraint is being bumped up against.
49 #define NSOLN_RETN_FAIL_STEPTOOSMALL -1
50 //! The nonlinear problem didn't solve the problem
51 #define NSOLN_RETN_FAIL_DAMPSTEP -2
52 //! The nonlinear problem's jacobian is singular
53 #define NSOLN_RETN_MATRIXINVERSIONERROR -3
54 //! The nonlinear problem's jacobian formation produced an error
55 #define NSOLN_RETN_JACOBIANFORMATIONERROR -4
56 //! The nonlinear problem's base residual produced an error
57 #define NSOLN_RETN_RESIDUALFORMATIONERROR -5
58 //! The nonlinear problem's max number of iterations has been exceeded
59 #define NSOLN_RETN_MAXIMUMITERATIONSEXCEEDED -7
60 //@}
61 //@}
62 
63 //@{
64 /// @name Constant which determines the type of the Jacobian
65 //! The jacobian will be calculated from a numerical method
66 #define NSOLN_JAC_NUM 1
67 //! The jacobian is calculated from an analytical function
68 #define NSOLN_JAC_ANAL 2
69 //@}
70 
71 
72 //! Class that calculates the solution to a nonlinear system
73 /*!
74  * This is a small nonlinear solver that can solve highly nonlinear problems that
75  * must use a dense matrix to relax the system.
76  *
77  * Newton's method is used.
78  *
79  * Damping is used extensively when relaxing the system.
80  *
81  *
82  * The basic idea is that we predict a direction that is parameterized by an overall coordinate
83  * value, beta, from zero to one, This may or may not be the same as the value, damp,
84  * depending upon whether the direction is straight.
85  *
86  *
87  * TIME STEP TYPE
88  *
89  * The code solves a nonlinear problem. Frequently the nonlinear problem is created from time-dependent
90  * residual. Whenever you change the solution vector, you are also changing the derivative of the
91  * solution vector. Therefore, the code has the option of altering ydot, a vector of time derivatives
92  * of the solution in tandem with the solution vector and then feeding a residual and Jacobian routine
93  * with the time derivatives as well as the solution. The code has support for a backwards euler method
94  * and a second order Adams-Bashforth or Trapezoidal Rule.
95  *
96  * In order to use these methods, the solver must be initialized with delta_t and m_y_nm1[i] to specify
97  * the conditions at the previous time step. For second order methods, the time derivative at t_nm1 must
98  * also be supplied, m_ydot_nm1[i]. Then the solution type NSOLN_TYPE_TIME_DEPENDENT may be used to
99  * solve the problem.
100  *
101  * For steady state problem whose residual doesn't have a solution time derivative in it, you should
102  * use the NSOLN_TYPE_STEADY_STATE problem type.
103  *
104  * We have a NSOLN_TYPE_PSEUDO_TIME_DEPENDENT defined. However, this is not implemented yet. This would
105  * be a pseudo time dependent calculation, where an optional time derivative could be added in order to
106  * help equilibrate a nonlinear steady state system. The time transient is not important in and of
107  * itself. Many physical systems have a time dependence to them that provides a natural way to relax
108  * the nonlinear system.
109  *
110  * MATRIX SCALING
111  *
112  *
113  *
114  *
115  * @code
116  *
117  *
118  * NonlinearSolver *nls = new NonlinearSolver(&r1);
119  *
120  * int solnType = NSOLN_TYPE_STEADY_STATE ;
121  *
122  * nls->setDeltaBoundsMagnitudes(deltaBounds);
123  *
124  * nls->solve_nonlinear_problem(solnType, y_comm, ydot_comm, CJ, time_curr, jac,
125  * num_newt_its, num_linear_solves, numBacktracks,
126  * loglevelInput);
127  *
128  * @endcode
129  *
130  *
131  * @ingroup numerics
132  */
134 {
135 
136 public:
137  //! Default constructor
138  /*!
139  * @param func Residual and jacobian evaluator function object
140  */
142 
143  //!Copy Constructor for the %ThermoPhase object.
144  /*!
145  * @param right Item to be copied
146  */
147  NonlinearSolver(const NonlinearSolver& right);
148 
149  //! Destructor
151 
152  //! Assignment operator
153  /*!
154  * This is NOT a virtual function.
155  *
156  * @param right Reference to %NonlinearSolver object to be
157  * copied into the
158  * current one.
159  */
161 
162  //! Create solution weights for convergence criteria
163  /*!
164  * We create soln weights from the following formula
165  *
166  * wt[i] = rtol * abs(y[i]) + atol[i]
167  *
168  * The program always assumes that atol is specific
169  * to the solution component
170  *
171  * @param y vector of the current solution values
172  */
173  void createSolnWeights(const doublereal* const y);
174 
175  //! L2 norm of the delta of the solution vector
176  /*!
177  * calculate the norm of the solution vector. This will
178  * involve the column scaling of the matrix
179  *
180  * The third argument has a default of false. However,
181  * if true, then a table of the largest values is printed
182  * out to standard output.
183  *
184  * @param delta_y Vector to take the norm of
185  * @param title Optional title to be printed out
186  * @param printLargest int indicating how many specific lines should be printed out
187  * @param dampFactor Current value of the damping factor. Defaults to 1.
188  * only used for printout out a table.
189  *
190  * @return Returns the L2 norm of the delta
191  */
192  doublereal solnErrorNorm(const doublereal* const delta_y, const char* title = 0, int printLargest = 0,
193  const doublereal dampFactor = 1.0) const;
194 
195  //! L2 norm of the residual of the equation system
196  /*!
197  * Calculate the norm of the residual vector. This may
198  * involve using the row sum scaling from the matrix problem.
199  *
200  * The second argument has a default of false. However,
201  * if true, then a table of the largest values is printed
202  * out to standard output.
203  *
204  * @param resid Vector of the residuals
205  * @param title Optional title to be printed out
206  * @param printLargest Number of specific entries to be printed
207  * @param y Current value of y - only used for printouts
208  *
209  *
210  * @return Returns the L2 norm of the delta
211  */
212  doublereal residErrorNorm(const doublereal* const resid, const char* title = 0, const int printLargest = 0,
213  const doublereal* const y = 0) const;
214 
215  //! Compute the current residual
216  /*!
217  * The current value of the residual is stored in the internal work array m_resid, which is defined
218  * as mutable
219  *
220  * @param time_curr Value of the time
221  * @param typeCalc Type of the calculation
222  * @param y_curr Current value of the solution vector
223  * @param ydot_curr Current value of the time derivative of the solution vector
224  * @param evalType Base evaluation type
225  * Defaults to Base_ResidEval
226  *
227  * @return Returns a flag to indicate that operation is successful.
228  * 1 Means a successful operation
229  * -0 or neg value Means an unsuccessful operation
230  */
231  int doResidualCalc(const doublereal time_curr, const int typeCalc, const doublereal* const y_curr,
232  const doublereal* const ydot_curr,
233  const ResidEval_Type_Enum evalType = Base_ResidEval) const;
234 
235  //! Compute the undamped Newton step
236  /*!
237  *
238  * Compute the undamped Newton step. The residual function is
239  * evaluated at the current time, t_n, at the current values of the
240  * solution vector, m_y_n, and the solution time derivative, m_ydot_n.
241  * The Jacobian is not recomputed.
242  *
243  * A factored jacobian is reused, if available. If a factored jacobian
244  * is not available, then the jacobian is factored. Before factoring,
245  * the jacobian is row and column-scaled. Column scaling is not
246  * recomputed. The row scales are recomputed here, after column
247  * scaling has been implemented.
248  *
249  * @param time_curr Current value of the time
250  * @param y_curr Current value of the solution
251  * @param ydot_curr Current value of the solution derivative.
252  * @param delta_y return value of the raw change in y
253  * @param jac Jacobian
254  *
255  * @return Returns the result code from lapack. A zero means success. Anything
256  * else indicates a failure.
257  */
258  int doNewtonSolve(const doublereal time_curr, const doublereal* const y_curr,
259  const doublereal* const ydot_curr, doublereal* const delta_y,
260  GeneralMatrix& jac);
261 
262  //! Compute the newton step, either by direct newton's or by solving a close problem that is represented
263  //! by a Hessian (
264  /*!
265  * This is algorith A.6.5.1 in Dennis / Schnabel
266  *
267  * Compute the QR decomposition
268  *
269  * Compute the undamped Newton step. The residual function is
270  * evaluated at the current time, t_n, at the current values of the
271  * solution vector, m_y_n, and the solution time derivative, m_ydot_n.
272  * The Jacobian is not recomputed.
273  *
274  * A factored jacobian is reused, if available. If a factored jacobian
275  * is not available, then the jacobian is factored. Before factoring,
276  * the jacobian is row and column-scaled. Column scaling is not
277  * recomputed. The row scales are recomputed here, after column
278  * scaling has been implemented.
279  *
280  * @param y_curr Current value of the solution
281  * @param ydot_curr Current value of the solution derivative.
282  * @param delta_y return value of the raw change in y
283  * @param jac Jacobian
284  *
285  * Internal input
286  * ---------------
287  * internal m_resid Stored residual is used as input
288  *
289  *
290  * @return Returns the result code from lapack. A zero means success. Anything
291  * else indicates a failure.
292  */
293  int doAffineNewtonSolve(const doublereal* const y_curr, const doublereal* const ydot_curr,
294  doublereal* const delta_y, GeneralMatrix& jac);
295 
296  //! Calculate the length of the current trust region in terms of the solution error norm
297  /*!
298  * We carry out a norm of deltaX_trust_ first. Then, we multiply that value
299  * by trustDelta_
300  */
301  doublereal trustRegionLength() const;
302 
303  //! Set default deulta bounds amounts
304  /*!
305  * Delta bounds are set to 0.01 for all unknowns arbitrarily and capriciously
306  * Then, for each call to the nonlinear solver
307  * Then, they are increased to 1000 x atol
308  * then, they are increased to 0.1 fab(y[i])
309  */
311 
312  //! Adjust the step minimums
313  void adjustUpStepMinimums();
314 
315  //! Set the delta Bounds magnitudes by hand
316  /*!
317  * @param deltaBoundsMagnitudes set the deltaBoundsMagnitude vector
318  */
319  void setDeltaBoundsMagnitudes(const doublereal* const deltaBoundsMagnitudes);
320 
321 protected:
322 
323  //! Readjust the trust region vectors
324  /*!
325  * The trust region is made up of the trust region vector calculation and the trustDelta_ value
326  * We periodically recalculate the trustVector_ values so that they renormalize to the
327  * correct length. We change the trustDelta_ values regularly
328  *
329  * The trust region calculate is based on
330  *
331  * || delta_x dot 1/trustDeltaX_ || <= trustDelta_
332  *
333  */
334  void readjustTrustVector();
335 
336  //! Fill a dogleg solution step vector
337  /*!
338  * Previously, we have filled up deltaX_Newton_[], deltaX_CP_[], and Nuu_, so that
339  * this routine is straightforward.
340  *
341  * @param leg Leg of the dog leg you are on (0, 1, or 2)
342  * @param alpha Relative length along the dog length that you are on.
343  * @param deltaX Vector to be filled up
344  */
345  void fillDogLegStep(int leg, doublereal alpha, std::vector<doublereal> & deltaX) const;
346 
347  //! Calculate the trust distance of a step in the solution variables
348  /*!
349  * The trust distance is defined as the length of the step according to the norm wrt to the trust region.
350  * We calculate the trust distance by the following method.
351  *
352  * trustDist = || delta_x dot 1/trustDeltaX_ ||
353  *
354  * @param deltaX Current value of deltaX
355  */
356  doublereal calcTrustDistance(std::vector<doublereal> const& deltaX) const;
357 
358 
359 
360 public:
361  //! Bound the step
362  /*!
363  *
364  * Return the factor by which the undamped Newton step 'step0'
365  * must be multiplied in order to keep all solution components in
366  * all domains between their specified lower and upper bounds.
367  * Other bounds may be applied here as well.
368  *
369  * Currently the bounds are hard coded into this routine:
370  *
371  * Minimum value for all variables: - 0.01 * m_ewt[i]
372  * Maximum value = none.
373  *
374  * Thus, this means that all solution components are expected
375  * to be numerical greater than zero in the limit of time step
376  * truncation errors going to zero.
377  *
378  * Delta bounds: The idea behind these is that the Jacobian
379  * couldn't possibly be representative if the
380  * variable is changed by a lot. (true for
381  * nonlinear systems, false for linear systems)
382  * Maximum increase in variable in any one newton iteration:
383  * factor of 2
384  * Maximum decrease in variable in any one newton iteration:
385  * factor of 5
386  *
387  * @param y Current solution value of the old step
388  * @param step0 Proposed step change in the solution
389  *
390  * @return Returns the damping factor determined by the bounds calculation
391  */
392  doublereal boundStep(const doublereal* const y, const doublereal* const step0);
393 
394  //! Set bounds constraints for all variables in the problem
395  /*!
396  *
397  * @param y_low_bounds Vector of lower bounds
398  * @param y_high_bounds Vector of high bounds
399  */
400  void setBoundsConstraints(const doublereal* const y_low_bounds,
401  const doublereal* const y_high_bounds);
402 
403  //! Return an editable vector of the low bounds constraints
404  std::vector<doublereal> & lowBoundsConstraintVector();
405 
406  //! Return an editable vector of the high bounds constraints
407  std::vector<doublereal> & highBoundsConstraintVector();
408 
409  //! Internal function to calculate the time derivative of the solution at the new step
410  /*!
411  * Previously, the user must have supplied information about the previous time step for this routine to
412  * work as intended.
413  *
414  * @param order of the BDF method
415  * @param y_curr current value of the solution
416  * @param ydot_curr Calculated value of the solution derivative that is consistent with y_curr
417  */
418  void calc_ydot(const int order, const doublereal* const y_curr, doublereal* const ydot_curr) const;
419 
420  //! Function called to evaluate the jacobian matrix and the current
421  //! residual vector at the current time step
422  /*!
423  *
424  *
425  * @param J Jacobian matrix to be filled in
426  * @param f Right hand side. This routine returns the current
427  * value of the rhs (output), so that it does
428  * not have to be computed again.
429  * @param time_curr Current time
430  * @param CJ inverse of the value of deltaT
431  * @param y value of the solution vector
432  * @param ydot value of the time derivative of the solution vector
433  * @param num_newt_its Number of newton iterations
434  *
435  * @return Returns a flag to indicate that operation is successful.
436  * 1 Means a successful operation
437  * 0 Means an unsuccessful operation
438  */
439  int beuler_jac(GeneralMatrix& J, doublereal* const f,
440  doublereal time_curr, doublereal CJ, doublereal* const y,
441  doublereal* const ydot, int num_newt_its);
442 
443  //! Apply a filtering process to the step
444  /*!
445  * @param timeCurrent Current value of the time
446  * @param ybase current value of the solution
447  * @param step0 Proposed step change in the solution
448  *
449  * @return Returns the norm of the value of the amount filtered
450  */
451  doublereal filterNewStep(const doublereal timeCurrent, const doublereal* const ybase, doublereal* const step0);
452 
453  //! Apply a filter to the solution
454  /*!
455  * @param timeCurrent Current value of the time
456  * @param y_current current value of the solution
457  * @param ydot_current Current value of the solution derivative.
458  *
459  * @return Returns the norm of the value of the amount filtered
460  */
461  doublereal filterNewSolution(const doublereal timeCurrent, doublereal* const y_current,
462  doublereal* const ydot_current);
463 
464  //! Return the factor by which the undamped Newton step 'step0'
465  //! must be multiplied in order to keep the update within the bounds of an accurate jacobian.
466  /*!
467  *
468  * The idea behind these is that the Jacobian couldn't possibly be representative, if the
469  * variable is changed by a lot. (true for nonlinear systems, false for linear systems)
470  * Maximum increase in variable in any one newton iteration:
471  * factor of 1.5
472  * Maximum decrease in variable in any one newton iteration:
473  * factor of 2
474  *
475  * @param y Initial value of the solution vector
476  * @param step0 initial proposed step size
477  *
478  * @return returns the damping factor
479  */
480  doublereal deltaBoundStep(const doublereal* const y, const doublereal* const step0);
481 
482  //! Find a damping coefficient through a look-ahead mechanism
483  /*!
484  * On entry, step_1 must contain an undamped Newton step for the
485  * solution y_n_curr. This method attempts to find a damping coefficient
486  * such that all components stay in bounds, and the next
487  * undamped step would have a norm smaller than
488  * that of step_1. If successful, the new solution after taking the
489  * damped step is returned in y_n_1, and the undamped step at y_n_1 is
490  * returned in step_2.
491  *
492  * @param time_curr Current physical time
493  * @param y_n_curr Base value of the solution before any steps
494  * are taken
495  * @param ydot_n_curr Base value of the time derivative of the
496  * solution
497  * @param step_1 Initial step suggested.
498  * @param y_n_1 Value of y1, the suggested solution after damping
499  * @param ydot_n_1 Value of the time derivative of the solution at y_n_1
500  * @param step_2 Value of the step change from y_n_1 to y_n_2
501  * @param stepNorm_2 norm of the step change in going from y_n_1 to y_n_2
502  * @param jac Jacobian
503  * @param writetitle Write a title line
504  * @param num_backtracks Number of backtracks taken
505  *
506  * @return returns an integer indicating what happened.
507  */
508  int dampStep(const doublereal time_curr, const doublereal* const y_n_curr,
509  const doublereal* const ydot_n_curr, doublereal* const step_1,
510  doublereal* const y_n_1, doublereal* const ydot_n_1, doublereal* step_2,
511  doublereal& stepNorm_2, GeneralMatrix& jac, bool writetitle,
512  int& num_backtracks);
513 
514  //! Find the solution to F(X) = 0 by damped Newton iteration.
515  /*!
516  * On
517  * entry, x0 contains an initial estimate of the solution. On
518  * successful return, x1 contains the converged solution.
519  *
520  * SolnType = TRANSIENT -> we will assume we are relaxing a transient
521  * equation system for now. Will make it more general later,
522  * if an application comes up.
523  *
524  * @param SolnType Solution type
525  * @param y_comm Initial value of the solution. On return this is the converged
526  * value of the solution
527  * @param ydot_comm Initial value of the solution derivative. On return this is the
528  * converged value of the solution derivative.
529  * @param CJ Inverse of the value of deltaT
530  * @param time_curr Current value of the time
531  * @param jac Matrix that will be used to store the jacobian
532  * @param num_newt_its Number of newton iterations taken
533  * @param num_linear_solves Number of linear solves taken
534  * @param num_backtracks Number of backtracking steps taken
535  * @param loglevelInput Input log level determines the amount of printing.
536  *
537  *
538  * @return A positive value indicates a successful convergence
539  * -1 Failed convergence
540  */
541  int solve_nonlinear_problem(int SolnType, doublereal* const y_comm, doublereal* const ydot_comm, doublereal CJ,
542  doublereal time_curr, GeneralMatrix& jac, int& num_newt_its,
543  int& num_linear_solves, int& num_backtracks, int loglevelInput);
544 
545 private:
546  //! Set the column scales
547  void calcColumnScales();
548 
549 public:
550 
551  //! Set the column scaling that are used for the inversion of the matrix
552  /*!
553  * There are three ways to do this.
554  *
555  * The first method is to set the bool useColScaling to true, leaving the scaling factors unset.
556  * Then, the column scales will be set to the solution error weighting factors. This has the
557  * effect of ensuring that all delta variables will have the same order of magnitude at convergence
558  * end.
559  *
560  * The second way is the explicity set the column factors in the second argument of this function call.
561  *
562  * The final way to input the scales is to override the ResidJacEval member function call,
563  *
564  * calcSolnScales(double time_n, const double *m_y_n_curr, const double *m_y_nm1, double *m_colScales)
565  *
566  * Overriding this function call will trump all other ways to specify the column scaling factors.
567  *
568  * @param useColScaling Turn this on if you want to use column scaling in the calculations
569  * @param scaleFactors A vector of doubles that specifies the column factors.
570  */
571  void setColumnScaling(bool useColScaling, const double* const scaleFactors = 0);
572 
573 
574  //! Set the rowscaling that are used for the inversion of the matrix
575  /*!
576  * Row scaling is set here. Right now the row scaling is set internally in the code.
577  *
578  * @param useRowScaling Turn row scaling on or off.
579  */
580  void setRowScaling(bool useRowScaling);
581 
582  //! Scale the matrix
583  /*!
584  * @param jac Jacobian
585  * @param y_comm Current value of the solution vector
586  * @param ydot_comm Current value of the time derivative of the solution vector
587  * @param time_curr current value of the time
588  * @param num_newt_its Current value of the number of newt its
589  */
590  void scaleMatrix(GeneralMatrix& jac, doublereal* const y_comm, doublereal* const ydot_comm,
591  doublereal time_curr, int num_newt_its);
592 
593  //! Print solution norm contribution
594  /*!
595  * Prints out the most important entries to the update to the solution vector for the current step
596  *
597  * @param step_1 Raw update vector for the current nonlinear step
598  * @param stepNorm_1 Norm of the vector step_1
599  * @param step_2 Raw update vector for the next solution value based on the old matrix
600  * @param stepNorm_2 Norm of the vector step_2
601  * @param title title of the printout
602  * @param y_n_curr Old value of the solution
603  * @param y_n_1 New value of the solution after damping corrections
604  * @param damp Value of the damping factor
605  * @param num_entries Number of entries to print out
606  */
607  void
608  print_solnDelta_norm_contrib(const doublereal* const step_1, const char* const stepNorm_1,
609  const doublereal* const step_2, const char* const stepNorm_2,
610  const char* const title, const doublereal* const y_n_curr,
611  const doublereal* const y_n_1, doublereal damp, size_t num_entries);
612 
613  //! Compute the Residual Weights
614  /*!
615  * The residual weights are defined here to be equal to the inverse of the row scaling factors used to
616  * row scale the matrix, after column scaling is used. They are multiplied by 10-3 because the column
617  * weights are also multiplied by that same quantity.
618  *
619  * The basic idea is that a change in the solution vector on the order of the convergence tolerance
620  * multiplied by [RJC] which is of order one after row scaling should give you the relative weight
621  * of the row. Values of the residual for that row can then be normalized by the value of this weight.
622  * When the tolerance in delta x is achieved, the tolerance in the residual is also achieved.
623  */
624  void computeResidWts();
625 
626  //! Return the residual weights
627  /*!
628  * @param residWts Vector of length neq_
629  */
630  void getResidWts(doublereal* const residWts) const;
631 
632 
633 
634  //! Check to see if the nonlinear problem has converged
635  /*!
636  *
637  * @param dampCode Code from the damping routine
638  * @param s1 Value of the norm of the step change
639  *
640  * @return integer is returned. If positive, then the problem has converged
641  * 1 Successful step was taken: Next step's norm is less than 1.0.
642  * The final residual norm is less than 1.0.
643  * 2 Successful step: Next step's norm is less than 0.8.
644  * This step's norm is less than 1.0.
645  * The residual norm can be anything.
646  * 3 Success: The final residual is less than 1.0
647  * The predicted deltaSoln is below 1.0.
648  * 0 Not converged yet
649  */
650  int convergenceCheck(int dampCode, doublereal s1);
651 
652  //! Set the absolute tolerances for the solution variables
653  /*!
654  * Set the absolute tolerances used in the calculation
655  *
656  * @param atol Vector of length neq_ that contains the tolerances to be used for the solution variables
657  */
658  void setAtol(const doublereal* const atol);
659 
660  //! Set the relative tolerances for the solution variables
661  /*!
662  * Set the relative tolerances used in the calculation
663  *
664  * @param rtol single double
665  */
666  void setRtol(const doublereal rtol);
667 
668  //! Set the relative and absolute tolerances for the Residual norm comparisons, if used
669  /*!
670  * Residual norms are used to calculate convergence within the nonlinear solver, since
671  * these are the norms that are associated with convergence proofs, especially for ill-conditioned systems.
672  * Usually the residual weights for each row are calculated by the program such that they
673  * correlate with the convergence requirements on the solution variables input by the user using
674  * the routines setAtol() and setRtol().
675  * The residual weights are essentially calculated from the value
676  *
677  * residWeightNorm[i] = m_ScaleSolnNormToResNorm * sum_j ( fabs(A_i,j) ewt(j))
678  *
679  * The factor, m_ScaleSolnNormToResNorm, is computed periodically to ensure that the solution norms
680  * and the residual norms are converging at the same time and thus accounts for some-illconditioning issues
681  * but not all.
682  *
683  * The user specified tolerance for the residual is given by the following quantity
684  *
685  * residWeightNorm[i] = residAtol[i] + residRtol * m_rowWtScales[i] / neq
686  *
687  * @param residNormHandling Parameter that sets the default handling of the residual norms
688  * 0 The residual weighting vector is calculated to make sure that the solution
689  * norms are roughly 1 when the residual norm is roughly 1.
690  * This is the default if this routine is not called.
691  * 1 Use the user residual norm specified by the parameters in this routine
692  * 2 Use the minimum value of the residual weights calculcated by method 1 and 2.
693  * This is the default if this routine is called and this parameter isn't specified.
694  */
695  void setResidualTols(double residRtol, double* residATol, int residNormHandling = 2);
696 
697  //! Set the value of the maximum # of newton iterations
698  /*!
699  * @param maxNewtIts Maximum number of newton iterations
700  */
701  void setMaxNewtIts(const int maxNewtIts);
702 
703  //! Calculate the scaling factor for translating residual norms into solution norms.
704  /*!
705  * This routine calls computeResidWts() a couple of times in the calculation of m_ScaleSolnNormToResNorm.
706  * A more sophisticated routine may do more with signs to get a better value. Perhaps, a series of calculations
707  * with different signs attached may be in order. Then, m_ScaleSolnNormToResNorm would be calculated
708  * as the minimum of a series of calculations.
709  */
711 
712  //! Calculate the steepest descent direction and the Cauchy Point where the quadratic formulation
713  //! of the nonlinear problem expects a minimum along the descent direction.
714  /*!
715  * @param jac Jacobian matrix: must be unfactored.
716  *
717  * @return Returns the norm of the solution update
718  */
719  doublereal doCauchyPointSolve(GeneralMatrix& jac);
720 
721  //! This is a utility routine that can be used to print out the rates of the initial residual decline
722  /*!
723  * The residual**2 decline for various directions is printed out. The rate of decline of the
724  * square of the residuals multiplied by the number of equations along each direction is printed out
725  * This quantity can be directly related to the theory, and may be calculated from derivatives at the
726  * original point.
727  *
728  * ( (r)**2 * neq - (r0)**2 * neq ) / distance
729  *
730  * What's printed out:
731  *
732  * The theoretical linearized residual decline
733  * The actual residual decline in the steepest descent direction determined by numerical differencing
734  * The actual residual decline in the newton direction determined by numerical differencing
735  *
736  * This routine doesn't need to be called for the solution of the nonlinear problem.
737  *
738  * @param time_curr Current time
739  * @param ydot0 INPUT Current value of the derivative of the solution vector
740  * @param ydot1 INPUT Time derivatives of solution at the conditions which are evaluated for success
741  * @param numTrials OUTPUT Counter for the number of residual evaluations
742  */
743  void descentComparison(doublereal time_curr ,doublereal* ydot0, doublereal* ydot1, int& numTrials);
744 
745 
746  //! Setup the parameters for the double dog leg
747  /*!
748  * The calls to the doCauchySolve() and doNewtonSolve() routines are done at the main level. This routine comes
749  * after those calls. We calculate the point Nuu_ here, the distances of the dog-legs,
750  * and the norms of the CP and Newton points in terms of the trust vectors.
751  */
752  void setupDoubleDogleg();
753 
754  //! Change the global lambda coordinate into the (leg,alpha) coordinate for the double dogleg
755  /*!
756  * @param lambda Global value of the distance along the double dogleg
757  * @param alpha relative value along the particular leg
758  *
759  * @return Returns the leg number ( 0, 1, or 2).
760  */
761  int lambdaToLeg(const doublereal lambda, doublereal& alpha) const;
762 
763  //! Given a trust distance, this routine calculates the intersection of the this distance with the
764  //! double dogleg curve
765  /*!
766  * @param trustVal (INPUT) Value of the trust distance
767  * @param lambda (OUTPUT) Returns the internal coordinate of the double dogleg
768  * @param alpha (OUTPUT) Returns the relative distance along the appropriate leg
769  * @return leg (OUTPUT) Returns the leg ID (0, 1, or 2)
770  */
771  int calcTrustIntersection(doublereal trustVal, doublereal& lambda, doublereal& alpha) const;
772 
773  //! Initialize the size of the trust vector.
774  /*!
775  * The algorithm we use is to set it equal to the length of the Distance to the Cauchy point.
776  */
777  void initializeTrustRegion();
778 
779  //! Set Trust region initialization strategy
780  /*!
781  * The default is use method 2 with a factor of 1.
782  * Then, on subsequent invocations of solve_nonlinear_problem() the strategy flips to method 0.
783  *
784  * @param method Method to set the strategy
785  * 0 No strategy - Use the previous strategy
786  * 1 Factor of the solution error weights
787  * 2 Factor of the first Cauchy Point distance
788  * 3 Factor of the first Newton step distance
789  *
790  * @param factor Factor to use in combination with the method
791  *
792  */
793  void setTrustRegionInitializationMethod(int method, doublereal factor);
794 
795 
796  //! Damp using the dog leg approach
797  /*!
798  *
799  * @param time_curr INPUT Current value of the time
800  * @param y_n_curr INPUT Current value of the solution vector
801  * @param ydot_n_curr INPUT Current value of the derivative of the solution vector
802  * @param step_1 INPUT First trial step for the first iteration
803  * @param y_n_1 INPUT First trial value of the solution vector
804  * @param ydot_n_1 INPUT First trial value of the derivative of the solution vector
805  * @param stepNorm_1 OUTPUT Norm of the vector step_1
806  * @param stepNorm_2 OUTPUT Estimated norm of the vector step_2
807  * @param jac INPUT jacobian
808  * @param num_backtracks OUTPUT number of backtracks taken in the current damping step
809  *
810  * @return 1 Successful step was taken. The predicted residual norm is less than one
811  * 2 Successful step: Next step's norm is less than 0.8
812  * 3 Success: The final residual is less than 1.0
813  * A predicted deltaSoln1 is not produced however. s1 is estimated.
814  * 4 Success: The final residual is less than the residual
815  * from the previous step.
816  * A predicted deltaSoln1 is not produced however. s1 is estimated.
817  * 0 Uncertain Success: s1 is about the same as s0
818  * -2 Unsuccessful step.
819  */
820  int dampDogLeg(const doublereal time_curr, const doublereal* y_n_curr,
821  const doublereal* ydot_n_curr, std::vector<doublereal> & step_1,
822  doublereal* const y_n_1, doublereal* const ydot_n_1,
823  doublereal& stepNorm_1, doublereal& stepNorm_2, GeneralMatrix& jac, int& num_backtracks);
824 
825  //! Decide whether the current step is acceptable and adjust the trust region size
826  /*!
827  * This is an extension of algorithm 6.4.5 of Dennis and Schnabel.
828  *
829  * Here we decide whether to accept the current step
830  * At the end of the calculation a new estimate of the trust region is calculated
831  *
832  * @param time_curr INPUT Current value of the time
833  * @param leg INPUT Leg of the dogleg that we are on
834  * @param alpha INPUT Distance down that leg that we are on
835  * @param y_n_curr INPUT Current value of the solution vector
836  * @param ydot_n_curr INPUT Current value of the derivative of the solution vector
837  * @param step_1 INPUT Trial step
838  * @param y_n_1 OUTPUT Solution values at the conditions which are evaluated for success
839  * @param ydot_n_1 OUTPUT Time derivatives of solution at the conditions which are evaluated for success
840  * @param trustDeltaOld INPUT Value of the trust length at the old conditions
841  *
842  *
843  * @return This function returns a code which indicates whether the step will be accepted or not.
844  * 3 Step passed with flying colors. Try redoing the calculation with a bigger trust region.
845  * 2 Step didn't pass deltaF requirement. Decrease the size of the next trust region for a retry and return
846  * 0 The step passed.
847  * -1 The step size is now too small (||d || < 0.1). A really small step isn't decreasing the function.
848  * This is an error condition.
849  * -2 Current value of the solution vector caused a residual error in its evaluation.
850  * Step is a failure, and the step size must be reduced in order to proceed further.
851  */
852  int decideStep(const doublereal time_curr, int leg, doublereal alpha, const doublereal* const y_n_curr,
853  const doublereal* const ydot_n_curr,
854  const std::vector<doublereal> & step_1,
855  const doublereal* const y_n_1, const doublereal* const ydot_n_1, doublereal trustDeltaOld);
856 
857  //! Calculated the expected residual along the double dogleg curve.
858  /*!
859  * @param leg 0, 1, or 2 representing the curves of the dogleg
860  * @param alpha Relative distance along the particular curve.
861  *
862  * @return Returns the expected value of the residual at that point according to the quadratic model.
863  * The residual at the newton point will always be zero.
864  */
865  doublereal expectedResidLeg(int leg, doublereal alpha) const;
866 
867  //! Here we print out the residual at various points along the double dogleg, comparing against the quadratic model
868  //! in a table format
869  /*!
870  * @param time_curr INPUT current time
871  * @param ydot0 INPUT Current value of the derivative of the solution vector for non-time dependent
872  * determinations
873  * @param legBest OUTPUT leg of the dogleg that gives the lowest residual
874  * @param alphaBest OUTPUT distance along dogleg for best result.
875  */
876  void residualComparisonLeg(const doublereal time_curr, const doublereal* const ydot0, int& legBest,
877  doublereal& alphaBest) const;
878 
879  //! Set the print level from the nonlinear solver
880  /*!
881  *
882  * 0 -> absolutely nothing is printed for a single time step.
883  * 1 -> One line summary per solve_nonlinear call
884  * 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
885  * 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table
886  * 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed
887  * 5 -> Algorithm information on the nonlinear iterates are printed out
888  * 6 -> Additional info on the nonlinear iterates are printed out
889  * 7 -> Additional info on the linear solve is printed out.
890  * 8 -> Info on a per iterate of the linear solve is printed out.
891  *
892  * @param printLvl integer value
893  */
894  void setPrintLvl(int printLvl);
895 
896  //! Parameter to turn on solution solver schemes
897  /*!
898  * @param doDogLeg Parameter to turn on the double dog leg scheme
899  * Default is to always use a damping scheme in the Newton Direction.
900  * When this is nonzero, a model trust region approach is used using a double dog leg
901  * with the steepest descent direction used for small step sizes.
902  *
903  * @param doAffineSolve Parameter to turn on or off the solution of the system using a Hessian
904  * if the matrix has a bad condition number.
905  */
906  void setSolverScheme(int doDogLeg, int doAffineSolve);
907 
908 
909 
910  /*
911  * -----------------------------------------------------------------------------------------------------------------
912  * MEMBER DATA
913  * ------------------------------------------------------------------------------------------------
914  */
915 private:
916 
917  //! Pointer to the residual and jacobian evaluator for the
918  //! function
919  /*!
920  * See ResidJacEval.h for an evaluator.
921  */
923 
924  //! Solution type
926 
927  //! Local copy of the number of equations
928  size_t neq_;
929 
930  //! Soln error weights
931  std::vector<doublereal> m_ewt;
932 
933  //! Boolean indicating whether a manual delta bounds has been input.
935 
936  //! Soln Delta bounds magnitudes
937  std::vector<doublereal> m_deltaStepMinimum;
938 
939  //! Value of the delta step magnitudes
940  std::vector<doublereal> m_deltaStepMaximum;
941 
942  //! Vector containing the current solution vector within the nonlinear solver
943  std::vector<doublereal> m_y_n_curr;
944 
945  //! Vector containing the time derivative of the current solution vector within the nonlinear solver
946  //! (where applicable)
947  std::vector<doublereal> m_ydot_n_curr;
948 
949  //! Vector containing the solution at the previous time step
950  std::vector<doublereal> m_y_nm1;
951 
952  //! Vector containing the solution at the previous time step
953  std::vector<doublereal> m_y_n_1;
954 
955  //! Value of the solution time derivative at the new point that is to be considered
956  std::vector<doublereal> m_ydot_n_1;
957 
958  //! Value of the step to be taken in the solution
959  std::vector<doublereal> m_step_1;
960 
961  //! Vector of column scaling factors
962  std::vector<doublereal> m_colScales;
963 
964  //! Weights for normalizing the values of the residuals
965  /*!
966  * These are computed if row scaling, m_rowScaling, is turned on. They are calculated currently as the
967  * sum of the absolute values of the rows of the jacobian.
968  */
969  std::vector<doublereal> m_rowScales;
970 
971  //! Weights for normalizing the values of the residuals
972  /*!
973  * They are calculated as the sum of the absolute values of the jacobian
974  * multiplied by the solution weight function.
975  * This is carried out in scaleMatrix().
976  */
977  std::vector<doublereal> m_rowWtScales;
978 
979  //! Value of the residual for the nonlinear problem
980  mutable std::vector<doublereal> m_resid;
981 
982  //! Workspace of length neq_
983  mutable std::vector<doublereal> m_wksp;
984 
985  //! Workspace of length neq_
986  mutable std::vector<doublereal> m_wksp_2;
987 
988  /*****************************************************************************************
989  * INTERNAL WEIGHTS FOR TAKING SOLUTION NORMS
990  ******************************************************************************************/
991  //! Vector of residual weights
992  /*!
993  * These are used to establish useful and informative weighted norms of the residual vector.
994  */
995  std::vector<doublereal> m_residWts;
996 
997  //! Norm of the residual at the start of each nonlinear iteration
998  doublereal m_normResid_0;
999 
1000  //! Norm of the residual after it has been bounded
1001  doublereal m_normResid_Bound;
1002 
1003  //! Norm of the residual at the end of the first leg of the current iteration
1004  doublereal m_normResid_1;
1005 
1006  //! Norm of the residual at the end of the first leg of the current iteration
1007  doublereal m_normResid_full;
1008 
1009  //! Norm of the solution update created by the iteration in its raw, undamped form, using the solution norm
1011 
1012  //! Norm of the distance to the cauchy point using the solution norm
1014 
1015  //! Norm of the residual for a trial calculation which may or may not be used
1016  doublereal m_normResidTrial;
1017 
1018  //! Vector of the norm
1019  doublereal m_normResidPoints[15];
1020 
1021  //! Boolean indicating whether we should scale the residual
1022  mutable bool m_resid_scaled;
1023 
1024  /*****************************************************************************************
1025  * INTERNAL BOUNDARY INFO FOR SOLUTIONS
1026  *****************************************************************************************/
1027 
1028  //! Bounds vector for each species
1029  std::vector<doublereal> m_y_high_bounds;
1030 
1031  //! Lower bounds vector for each species
1032  std::vector<doublereal> m_y_low_bounds;
1033 
1034  //! Damping factor imposed by hard bounds and by delta bounds
1035  doublereal m_dampBound;
1036 
1037  //! Additional damping factor due to bounds on the residual and solution norms
1038  doublereal m_dampRes;
1039 
1040  //! Delta t for the current step
1041  doublereal delta_t_n;
1042 
1043  //! Counter for the total number of function evaluations
1044  mutable int m_nfe;
1045 
1046  /***********************************************************************************************
1047  * MATRIX INFORMATION
1048  **************************************************************************************/
1049 
1050  //! The type of column scaling used in the matrix inversion of the problem
1051  /*!
1052  * If 1 then colScaling = m_ewt[]
1053  * If 2 then colScaling = user set
1054  * if 0 then colScaling = 1.0
1055  */
1057 
1058  //! int indicating whether row scaling is turned on (1) or not (0)
1060 
1061  //! Total number of linear solves taken by the solver object
1063 
1064  //! Number of local linear solves done during the current iteration
1066 
1067  //! Total number of newton iterations
1069 
1070  //! Minimum number of newton iterations to use
1072 
1073  //! Maximum number of newton iterations
1075 
1076  //! Jacobian formation method
1077  /*!
1078  * 1 = numerical (default)
1079  * 2 = analytical
1080  */
1082 
1083  //! Number of Jacobian evaluations
1085 
1086  //! Current system time
1087  /*!
1088  * Note, we assume even for steady state problems that the residual
1089  * is a function of a system time.
1090  */
1091  doublereal time_n;
1092 
1093  //! Boolean indicating matrix conditioning
1095 
1096  //! Order of the time step method = 1
1097  int m_order;
1098 
1099  //! value of the relative tolerance to use in solving the equation set
1100  doublereal rtol_;
1101 
1102  //! Base value of the absolute tolerance
1103  doublereal atolBase_;
1104 
1105  //! Pointer containing the solution derivative at the previous time step
1106  doublereal* m_ydot_nm1;
1107 
1108  //! absolute tolerance in the solution unknown
1109  /*!
1110  * This is used to evaluating the weighting factor
1111  */
1112  std::vector<doublereal> atolk_;
1113 
1114  //! absolute tolerance in the unscaled solution unknowns
1115  std::vector<doublereal> userResidAtol_;
1116 
1117  //! absolute tolerance in the unscaled solution unknowns
1118  doublereal userResidRtol_;
1119 
1120  //! Check the residual tolerances explictly against user input
1121  /*!
1122  * 0 Don't calculate residual weights from residual tolerance inputs
1123  * 1 Calculate residual weights from residual tolerance inputs only
1124  * 2 Calculate residual weights from a minimum of the solution error weights process and the direct residual tolerance inputs
1125  */
1127 
1128  //! Determines the level of printing for each time step.
1129  /*!
1130  * 0 -> absolutely nothing is printed for a single time step.
1131  * 1 -> One line summary per solve_nonlinear call
1132  * 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
1133  * 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table
1134  * 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed
1135  * 5 -> Algorithm information on the nonlinear iterates are printed out
1136  * 6 -> Additional info on the nonlinear iterates are printed out
1137  * 7 -> Additional info on the linear solve is printed out.
1138  * 8 -> Info on a per iterate of the linear solve is printed out.
1139  */
1141 
1142  //! Scale factor for turning residual norms into solution norms
1144 
1145  //! Copy of the jacobian that doesn't get overwritten when the inverse is determined
1146  /*!
1147  * The jacobian stored here is the raw matrix, before any row or column scaling is carried out
1148  */
1150 
1151  //! Hessian
1153 
1154  /*********************************************************************************************
1155  * VARIABLES ASSOCIATED WITH STEPS AND ASSOCIATED DOUBLE DOGLEG PARAMETERS
1156  *********************************************************************************************/
1157 
1158  //! Steepest descent direction. This is also the distance to the Cauchy Point
1159  std::vector<doublereal> deltaX_CP_;
1160 
1161  //! Newton Step - This is the newton step determined from the straight Jacobian
1162  /*
1163  * Newton step for the current step only
1164  */
1165  std::vector<doublereal> deltaX_Newton_;
1166 
1167  //! Expected value of the residual norm at the Cauchy point if the quadratic model
1168  //! were valid
1169  doublereal residNorm2Cauchy_;
1170 
1171  //! Current leg
1173 
1174  //! Current Alpha param along the leg
1175  doublereal dogLegAlpha_;
1176 
1177  //! Residual dot Jd norm
1178  /*!
1179  * This is equal to R_hat dot J_hat d_y_descent
1180  */
1181  doublereal RJd_norm_;
1182 
1183  //! Value of lambdaStar_ which is used to calculate the Cauchy point
1184  doublereal lambdaStar_;
1185 
1186  //! Jacobian times the steepest descent direction in the normalized coordinates.
1187  /*!
1188  * This is equal to [ Jhat d^y_{descent} ] in the notes, Eqn. 18.
1189  */
1190  std::vector<doublereal> Jd_;
1191 
1192  //! Vector of trust region values.
1193  std::vector<doublereal> deltaX_trust_;
1194 
1195  //! Current norm of the vector deltaX_trust_ in terms of the solution norm
1196  mutable doublereal norm_deltaX_trust_;
1197 
1198  //! Current value of trust radius. This is used with deltaX_trust_ to
1199  //! calculate the max step size.
1200  doublereal trustDelta_;
1201 
1202  //! Method for handling the trust region initialization
1203  /*!
1204  * Then, on subsequent invocations of solve_nonlinear_problem() the strategy flips to method 0.
1205  *
1206  * method Method to set the strategy
1207  * 0 No strategy - Use the previous strategy
1208  * 1 Factor of the solution error weights
1209  * 2 Factor of the first Cauchy Point distance
1210  * 3 Factor of the first Newton step distance
1211  */
1213 
1214  //! Factor used to set the initial trust region
1216 
1217  //! Relative distance down the Newton step that the second dogleg starts
1218  doublereal Nuu_;
1219 
1220  //! Distance of the zeroeth leg of the dogleg in terms of the solution norm
1221  doublereal dist_R0_;
1222 
1223  //! Distance of the first leg of the dogleg in terms of the solution norm
1224  doublereal dist_R1_;
1225 
1226  //! Distance of the second leg of the dogleg in terms of the solution norm
1227  doublereal dist_R2_;
1228 
1229  //! Distance of the sum of all legs of the doglegs in terms of the solution norm
1230  doublereal dist_Total_;
1231 
1232  //! Dot product of the Jd_ variable defined above with itself.
1233  doublereal JdJd_norm_;
1234 
1235  //! Norm of the Newton Step wrt trust region
1236  doublereal normTrust_Newton_;
1237 
1238  //! Norm of the Cauchy Step direction wrt trust region
1239  doublereal normTrust_CP_;
1240 
1241  //! General toggle for turning on dog leg damping.
1243 
1244  //! General toggle for turning on Affine solve with Hessian
1246 
1247  //! Condition number of the matrix
1248  doublereal m_conditionNumber;
1249 
1250  //! Factor indicating how much trust region has been changed this iteration - output variable
1252 
1253  //! Factor indicating how much trust region has been changed next iteration - output variable
1254  doublereal NextTrustFactor_;
1255 
1256  //! Boolean indicating that the residual weights have been reevaluated this iteration - output variable
1258 
1259  //! Expected DResid_dS for the steepest descent path - output variable
1261 
1262  //! Actual DResid_dS for the steepest descent path - output variable
1263  doublereal ResidDecreaseSD_;
1264 
1265  //! Expected DResid_dS for the Newton path - output variable
1267 
1268  //! Actual DResid_dS for the newton path - output variable
1270 
1271  /*******************************************************************************************
1272  * STATIC VARIABLES
1273  *****************************************************************************************/
1274 
1275 public:
1276  //! Turn off printing of time
1277  /*!
1278  * Necessary to do for test suites
1279  */
1280  static bool s_TurnOffTiming;
1281 
1282  //! Turn on or off printing of the Jacobian
1283  static bool s_print_NumJac;
1284 
1285  //! Turn on extra printing of dogleg information
1286  static bool s_print_DogLeg;
1287 
1288  //! Turn on solving both the Newton and Hessian systems and comparing the results
1289  /*!
1290  * This is off by default
1291  */
1293 
1294  //! This toggle turns off the use of the Hessian when it is warranted by the condition number.
1295  /*!
1296  * This is a debugging option.
1297  */
1299 
1300 };
1301 
1302 }
1303 
1304 #endif