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