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