Cantera  2.1.2
ResidJacEval.h
Go to the documentation of this file.
1 /**
2  * @file ResidJacEval.h
3  *
4  * Dense, Square (not sparse) matrices.
5  */
6 
7 /*
8  * Copyright 2004 Sandia Corporation. Under the terms of Contract
9  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
10  * retains certain rights in this software.
11  * See file License.txt for licensing information.
12  */
13 
14 #ifndef CT_RESIDJACEVAL_H
15 #define CT_RESIDJACEVAL_H
16 
17 #include "ResidEval.h"
18 #include "GeneralMatrix.h"
19 
20 namespace Cantera
21 {
22 
23 //! Differentiates the type of residual evaluations according to functionality
25  //! Base residual calculation for the time-stepping function
27  //! Base residual calculation for the Jacobian calculation
29  //! Delta residual calculation for the Jacbobian calculation
31  //! Base residual calculation for the showSolution routine
32  /*!
33  * We calculate this when we want to display a solution
34  */
36  //! Base residual calculation containing any lagged components
37  /*!
38  * We use this to calculate residuals when doing line searches along
39  * directions determined by Jacobians that are missing contributions
40  * from lagged entries.
41  */
43 };
44 
45 //! Wrappers for the function evaluators for Nonlinear solvers and Time steppers
46 /*!
47  * A class for full (non-sparse dense matrices with Fortran-compatible data storage.
48  * The class adds support for identifying what types of calls are made to the residual
49  * evaluator by adding the ResidEval_Type_Enum class.
50  */
51 class ResidJacEval : public ResidEval
52 {
53 public:
54  //!Default constructor
55  /*!
56  * @param atol Initial value of the global tolerance (defaults to 1.0E-13)
57  */
58  ResidJacEval(doublereal atol = 1.0e-13);
59 
60  //!Copy Constructor
61  ResidJacEval(const ResidJacEval& right);
62 
63  //! Assignment operator
64  ResidJacEval& operator=(const ResidJacEval& right);
65 
66  //! Duplication routine for objects derived from residJacEval
67  /*!
68  * This virtual routine can be used to duplicate objects which inherit
69  * from ResidJacEval even if the application only has a pointer to
70  * ResidJacEval to work with.
71  *
72  * These routines are basically wrappers around the derived copy
73  * constructor.
74  */
75  virtual ResidJacEval* duplMyselfAsResidJacEval() const;
76 
77  //! Return the number of equations in the equation system
78  virtual int nEquations() const;
79 
80  //! Evaluate the residual function
81  /*!
82  * @param t Time (input)
83  * @param delta_t The current value of the time step (input)
84  * @param y Solution vector (input, do not modify)
85  * @param ydot Rate of change of solution vector. (input, do not modify)
86  * @param resid Value of the residual that is computed (output)
87  * @param evalType Type of the residual being computed (defaults to Base_ResidEval)
88  * @param id_x Index of the variable that is being numerically differenced to find
89  * the jacobian (defaults to -1, which indicates that no variable is being
90  * differenced or that the residual doesn't take this issue into account)
91  * @param delta_x Value of the delta used in the numerical differencing
92  *
93  * @return Returns a flag to indicate that operation is successful.
94  * 1 Means a successful operation
95  * -0 or neg value Means an unsuccessful operation
96  */
97  virtual int evalResidNJ(const doublereal t, const doublereal delta_t,
98  const doublereal* const y,
99  const doublereal* const ydot,
100  doublereal* const resid,
101  const ResidEval_Type_Enum evalType = Base_ResidEval,
102  const int id_x = -1,
103  const doublereal delta_x = 0.0);
104 
105  virtual int eval(const doublereal t, const doublereal* const y,
106  const doublereal* const ydot,
107  doublereal* const r);
108 
109  virtual int getInitialConditions(const doublereal t0, doublereal* const y, doublereal* const ydot);
110 
111  //! Filter the solution predictions
112  /*!
113  * Codes might provide a predicted step change. This routine filters the predicted
114  * solution vector eliminating illegal directions.
115  *
116  * @param t Time (input)
117  * @param ybase Solution vector (input, output)
118  * @param step Proposed step in the solution that will be cropped
119  *
120  * @return Return the norm of the amount of filtering
121  */
122  virtual doublereal filterNewStep(const doublereal t, const doublereal* const ybase,
123  doublereal* const step);
124 
125  //! Filter the solution predictions
126  /*!
127  * Codes might provide a predicted solution vector. This routine filters the predicted
128  * solution vector.
129  *
130  * @param t Time (input)
131  * @param y Solution vector (input, output)
132  *
133  * @return Return the norm of the amount of filtering
134  */
135  virtual doublereal filterSolnPrediction(const doublereal t, doublereal* const y);
136 
137  //! Set a global value of the absolute tolerance
138  /*!
139  * @param atol Value of atol
140  */
141  void setAtol(doublereal atol);
142 
143  //! Evaluate the time tracking equations, if any
144  /*!
145  * Evaluate time integrated quantities that are calculated at the
146  * end of every successful time step. This call is made once at the end of every successful
147  * time step that advances the time. It's also made once at the start of the time stepping.
148  *
149  * @param t Time (input)
150  * @param delta_t The current value of the time step (input)
151  * @param y Solution vector (input, do not modify)
152  * @param ydot Rate of change of solution vector. (input, do not modify)
153  *
154  * @return Returns a flag to indicate that operation is successful.
155  * 1 Means a successful operation
156  * -0 or neg value Means an unsuccessful operation
157  */
158  virtual int evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal* const y,
159  const doublereal* const ydot);
160 
161  //! Evaluate any stopping criteria other than a final time limit
162  /*!
163  * If we are to stop the time integration for any reason other than reaching a final time limit, tout,
164  * provide a test here. This call is made at the end of every successful time step iteration
165  *
166  * @return If true, the the time stepping is stopped. If false, then time stepping is stopped if t >= tout
167  * Defaults to false.
168  *
169  * @param t Time (input)
170  * @param delta_t The current value of the time step (input)
171  * @param y Solution vector (input, do not modify)
172  * @param ydot Rate of change of solution vector. (input, do not modify)
173  */
174  virtual bool evalStoppingCritera(const doublereal t,
175  const doublereal delta_t,
176  const doublereal* const y,
177  const doublereal* const ydot);
178 
179  //! Return a vector of delta y's for calculation of the numerical Jacobian
180  /*!
181  * There is a default algorithm provided.
182  *
183  * delta_y[i] = atol[i] + 1.0E-6 ysoln[i]
184  * delta_y[i] = atol[i] + MAX(1.0E-6 ysoln[i] * 0.01 * solnWeights[i])
185  *
186  * @param t Time (input)
187  * @param y Solution vector (input, do not modify)
188  * @param ydot Rate of change of solution vector. (input, do not modify)
189  * @param delta_y Value of the delta to be used in calculating the numerical jacobian
190  * @param solnWeights Value of the solution weights that are used in determining convergence (default = 0)
191  *
192  * @return Returns a flag to indicate that operation is successful.
193  * 1 Means a successful operation
194  * -0 or neg value Means an unsuccessful operation
195  */
196  virtual int
197  calcDeltaSolnVariables(const doublereal t,
198  const doublereal* const y,
199  const doublereal* const ydot,
200  doublereal* const delta_y,
201  const doublereal* const solnWeights = 0);
202 
203  //! Returns a vector of column scale factors that can be used to column scale Jacobians.
204  /*!
205  * Default to yScales[] = 1.0
206  *
207  * @param t Time (input)
208  * @param y Solution vector (input, do not modify)
209  * @param y_old Old Solution vector (input, do not modify)
210  * @param yScales Value of the column scales
211  */
212  virtual void calcSolnScales(const doublereal t, const doublereal* const y,
213  const doublereal* const y_old, doublereal* const yScales);
214 
215  //! This function may be used to create output at various points in the execution of an application.
216  /*!
217  * @param ifunc identity of the call
218  * 0 Initial call
219  * 1 Called at the end of every successful time step
220  * -1 Called at the end of every unsuccessful time step
221  * 2 Called at the end of every call to integrateRJE()
222  *
223  * @param t Time (input)
224  * @param delta_t The current value of the time step (input)
225  * @param y Solution vector (input, do not modify)
226  * @param ydot Rate of change of solution vector. (input)
227  */
228  virtual void user_out2(const int ifunc, const doublereal t,
229  const doublereal delta_t,
230  const doublereal* const y,
231  const doublereal* const ydot);
232 
233  //! This function may be used to create output at various points in the execution of an application.
234  /*!
235  * This routine calls user_out2().
236  *
237  * @param ifunc identity of the call
238  * @param t Time (input)
239  * @param y Solution vector (input, do not modify)
240  * @param ydot Rate of change of solution vector. (input)
241  */
242  virtual void user_out(const int ifunc, const doublereal t,
243  const doublereal* y,
244  const doublereal* ydot);
245 
246  //! Multiply the matrix by another matrix that leads to better conditioning
247  /*!
248  * Provide a left sided matrix that will multiply the current jacobian, after scaling
249  * and lead to a better conditioned system.
250  * This routine is called just before the matrix is factored.
251  *
252  * Original Problem:
253  * J delta_x = - Resid
254  *
255  * New problem:
256  * M (J delta_x) = - M Resid
257  *
258  * @param matrix Pointer to the current jacobian (if zero, it's already been factored)
259  * @param nrows offsets for the matrix
260  * @param rhs residual vector. This also needs to be lhs multiplied by M
261  *
262  * @return Returns a flag to indicate that operation is successful.
263  * 1 Means a successful operation
264  * -0 or neg value Means an unsuccessful operation
265  */
266  virtual int matrixConditioning(doublereal* const matrix, const int nrows,
267  doublereal* const rhs);
268 
269  //! Calculate an analytical jacobian and the residual at the current time and values.
270  /*!
271  * Only called if the jacFormation method is set to analytical
272  *
273  * @param t Time (input)
274  * @param delta_t The current value of the time step (input)
275  * @param cj Coefficient of yprime used in the evaluation of the jacobian
276  * @param y Solution vector (input, do not modify)
277  * @param ydot Rate of change of solution vector. (input, do not modify)
278  * @param J Reference to the SquareMatrix object to be calculated (output)
279  * @param resid Value of the residual that is computed (output)
280  *
281  * @return Returns a flag to indicate that operation is successful.
282  * 1 Means a successful operation
283  * -0 or neg value Means an unsuccessful operation
284  */
285  virtual int evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
286  const doublereal* const y, const doublereal* const ydot,
287  GeneralMatrix& J, doublereal* const resid);
288 
289  //! Calculate an analytical jacobian and the residual at the current time and values.
290  /*!
291  * Only called if the jacFormation method is set to analytical
292  *
293  * @param t Time (input)
294  * @param delta_t The current value of the time step (input)
295  * @param cj Coefficient of yprime used in the evaluation of the jacobian
296  * @param y Solution vector (input, do not modify)
297  * @param ydot Rate of change of solution vector. (input, do not modify)
298  * @param jacobianColPts Pointer to the vector of pts to columns of the SquareMatrix
299  * object to be calculated (output)
300  * @param resid Value of the residual that is computed (output)
301  *
302  * @return Returns a flag to indicate that operation is successful.
303  * 1 Means a successful operation
304  * -0 or neg value Means an unsuccessful operation
305  */
306  virtual int evalJacobianDP(const doublereal t, const doublereal delta_t, doublereal cj,
307  const doublereal* const y,
308  const doublereal* const ydot,
309  doublereal* const* jacobianColPts,
310  doublereal* const resid);
311 
312 protected:
313  //! constant value of atol
314  doublereal m_atol;
315 
316  //! Number of equations
317  int neq_;
318 };
319 }
320 
321 #endif
virtual int evalJacobianDP(const doublereal t, const doublereal delta_t, doublereal cj, const doublereal *const y, const doublereal *const ydot, doublereal *const *jacobianColPts, doublereal *const resid)
Calculate an analytical jacobian and the residual at the current time and values. ...
virtual void user_out2(const int ifunc, const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot)
This function may be used to create output at various points in the execution of an application...
Base residual calculation for the showSolution routine.
Definition: ResidJacEval.h:35
virtual int nEquations() const
Return the number of equations in the equation system.
virtual int evalResidNJ(const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot, doublereal *const resid, const ResidEval_Type_Enum evalType=Base_ResidEval, const int id_x=-1, const doublereal delta_x=0.0)
Evaluate the residual function.
virtual int evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj, const doublereal *const y, const doublereal *const ydot, GeneralMatrix &J, doublereal *const resid)
Calculate an analytical jacobian and the residual at the current time and values. ...
int neq_
Number of equations.
Definition: ResidJacEval.h:317
virtual int eval(const doublereal t, const doublereal *const y, const doublereal *const ydot, doublereal *const r)
Evaluate the residual function.
Generic matrix.
Definition: GeneralMatrix.h:23
Base residual calculation for the Jacobian calculation.
Definition: ResidJacEval.h:28
virtual int getInitialConditions(const doublereal t0, doublereal *const y, doublereal *const ydot)
Fill in the initial conditions.
ResidJacEval(doublereal atol=1.0e-13)
Default constructor.
virtual void user_out(const int ifunc, const doublereal t, const doublereal *y, const doublereal *ydot)
This function may be used to create output at various points in the execution of an application...
virtual ResidJacEval * duplMyselfAsResidJacEval() const
Duplication routine for objects derived from residJacEval.
ResidEval_Type_Enum
Differentiates the type of residual evaluations according to functionality.
Definition: ResidJacEval.h:24
Wrappers for the function evaluators for Nonlinear solvers and Time steppers.
Definition: ResidJacEval.h:51
virtual int calcDeltaSolnVariables(const doublereal t, const doublereal *const y, const doublereal *const ydot, doublereal *const delta_y, const doublereal *const solnWeights=0)
Return a vector of delta y's for calculation of the numerical Jacobian.
virtual int matrixConditioning(doublereal *const matrix, const int nrows, doublereal *const rhs)
Multiply the matrix by another matrix that leads to better conditioning.
Virtual base class for DAE residual function evaluators.
Definition: ResidEval.h:33
doublereal m_atol
constant value of atol
Definition: ResidJacEval.h:314
virtual int evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot)
Evaluate the time tracking equations, if any.
virtual void calcSolnScales(const doublereal t, const doublereal *const y, const doublereal *const y_old, doublereal *const yScales)
Returns a vector of column scale factors that can be used to column scale Jacobians.
ResidJacEval & operator=(const ResidJacEval &right)
Assignment operator.
Delta residual calculation for the Jacbobian calculation.
Definition: ResidJacEval.h:30
Declarations for the class GeneralMatrix which is a virtual base class for matrices handled by solver...
virtual doublereal filterSolnPrediction(const doublereal t, doublereal *const y)
Filter the solution predictions.
virtual bool evalStoppingCritera(const doublereal t, const doublereal delta_t, const doublereal *const y, const doublereal *const ydot)
Evaluate any stopping criteria other than a final time limit.
void setAtol(doublereal atol)
Set a global value of the absolute tolerance.
virtual doublereal filterNewStep(const doublereal t, const doublereal *const ybase, doublereal *const step)
Filter the solution predictions.
Base residual calculation for the time-stepping function.
Definition: ResidJacEval.h:26
Base residual calculation containing any lagged components.
Definition: ResidJacEval.h:42