Cantera  2.0
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  *
52  */
53 class ResidJacEval : public ResidEval
54 {
55 
56 public:
57 
58  //!Default constructor
59  /*!
60  * @param atol Initial value of the global tolerance (defaults to 1.0E-13)
61  */
62  ResidJacEval(doublereal atol = 1.0e-13);
63 
64  //!Copy Constructor for the %ResidJacEval object
65  /*!
66  * @param right Item to be copied
67  */
68  ResidJacEval(const ResidJacEval& right);
69 
70  /// Destructor. Does nothing.
71  virtual ~ResidJacEval();
72 
73  //! Assignment operator
74  /*!
75  * This is NOT a virtual function.
76  *
77  * @param right Reference to %ResidJacEval object to be copied into the
78  * current one.
79  */
80  ResidJacEval& operator=(const ResidJacEval& right);
81 
82  //! Duplication routine for objects which inherit from
83  //! residJacEval
84  /*!
85  * This virtual routine can be used to duplicate %ResidJacEval objects
86  * inherited from %ResidJacEval even if the application only has
87  * a pointer to %ResidJacEval to work with.
88  *
89  * These routines are basically wrappers around the derived copy
90  * constructor.
91  */
92  virtual ResidJacEval* duplMyselfAsResidJacEval() const;
93 
94  //! Return the number of equations in the equation system
95  virtual int nEquations() const;
96 
97  //! Evaluate the residual function
98  /*!
99  * @param t Time (input)
100  * @param delta_t The current value of the time step (input)
101  * @param y Solution vector (input, do not modify)
102  * @param ydot Rate of change of solution vector. (input, do not modify)
103  * @param resid Value of the residual that is computed (output)
104  * @param evalType Type of the residual being computed (defaults to Base_ResidEval)
105  * @param id_x Index of the variable that is being numerically differenced to find
106  * the jacobian (defaults to -1, which indicates that no variable is being
107  * differenced or that the residual doesn't take this issue into account)
108  * @param delta_x Value of the delta used in the numerical differencing
109  *
110  * @return Returns a flag to indicate that operation is successful.
111  * 1 Means a successful operation
112  * -0 or neg value Means an unsuccessful operation
113  */
114  virtual int evalResidNJ(const doublereal t, const doublereal delta_t,
115  const doublereal* const y,
116  const doublereal* const ydot,
117  doublereal* const resid,
118  const ResidEval_Type_Enum evalType = Base_ResidEval,
119  const int id_x = -1,
120  const doublereal delta_x = 0.0);
121 
122 
123  /**
124  * Evaluate the residual function. Called by the
125  * integrator.
126  * @param t time. (input)
127  * @param y solution vector. (input)
128  * @param ydot rate of change of solution vector. (input)
129  * @param r residual vector (output)
130  */
131  virtual int eval(const doublereal t, const doublereal* const y,
132  const doublereal* const ydot,
133  doublereal* const r);
134 
135 
136 
137  //! Fill in the initial conditions
138  /*!
139  * Values for both the solution and the value of ydot may be provided.
140  *
141  * @param t0 Time (input)
142  * @param y Solution vector (output)
143  * @param ydot Rate of change of solution vector. (output)
144  *
145  * @return Returns a flag to indicate that operation is successful.
146  * 1 Means a successful operation
147  * -0 or neg value Means an unsuccessful operation
148  */
149  virtual int getInitialConditions(const doublereal t0, doublereal* const y, doublereal* const ydot);
150 
151  //! Filter the solution predictions
152  /*!
153  * Codes might provide a predicted step change. This routine filters the predicted
154  * solution vector eliminating illegal directions.
155  *
156  * @param t Time (input)
157  * @param ybase Solution vector (input, output)
158  * @param step Proposed step in the solution that will be cropped
159  *
160  * @return Return the norm of the amount of filtering
161  */
162  virtual doublereal filterNewStep(const doublereal t, const doublereal* const ybase,
163  doublereal* const step);
164 
165  //! Filter the solution predictions
166  /*!
167  * Codes might provide a predicted solution vector. This routine filters the predicted
168  * solution vector.
169  *
170  * @param t Time (input)
171  * @param y Solution vector (input, output)
172  *
173  * @return Return the norm of the amount of filtering
174  */
175  virtual doublereal filterSolnPrediction(const doublereal t, doublereal* const y);
176 
177  //! Set a global value of the absolute tolerance
178  /*!
179  * @param atol Value of atol
180  */
181  void setAtol(doublereal atol);
182 
183  //! Evaluate the time tracking equations, if any
184  /*!
185  * Evaluate time integrated quantities that are calculated at the
186  * end of every successful time step. This call is made once at the end of every successful
187  * time step that advances the time. It's also made once at the start of the time stepping.
188  *
189  * @param t Time (input)
190  * @param delta_t The current value of the time step (input)
191  * @param y Solution vector (input, do not modify)
192  * @param ydot Rate of change of solution vector. (input, do not modify)
193  *
194  * @return Returns a flag to indicate that operation is successful.
195  * 1 Means a successful operation
196  * -0 or neg value Means an unsuccessful operation
197  */
198  virtual int evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal* const y,
199  const doublereal* const ydot);
200 
201  //! Evaluate any stopping criteria other than a final time limit
202  /*!
203  * If we are to stop the time integration for any reason other than reaching a final time limit, tout,
204  * provide a test here. This call is made at the end of every successful time step iteration
205  *
206  * @return If true, the the time stepping is stopped. If false, then time stepping is stopped if t >= tout
207  * Defaults to false.
208  *
209  * @param t Time (input)
210  * @param delta_t The current value of the time step (input)
211  * @param y Solution vector (input, do not modify)
212  * @param ydot Rate of change of solution vector. (input, do not modify)
213  */
214  virtual bool evalStoppingCritera(const doublereal t,
215  const doublereal delta_t,
216  const doublereal* const y,
217  const doublereal* const ydot);
218 
219  //! Return a vector of delta y's for calculation of the numerical Jacobian
220  /*!
221  * There is a default algorithm provided.
222  *
223  * delta_y[i] = atol[i] + 1.0E-6 ysoln[i]
224  * delta_y[i] = atol[i] + MAX(1.0E-6 ysoln[i] * 0.01 * solnWeights[i])
225  *
226  * @param t Time (input)
227  * @param y Solution vector (input, do not modify)
228  * @param ydot Rate of change of solution vector. (input, do not modify)
229  * @param delta_y Value of the delta to be used in calculating the numerical jacobian
230  * @param solnWeights Value of the solution weights that are used in determining convergence (default = 0)
231  *
232  * @return Returns a flag to indicate that operation is successful.
233  * 1 Means a successful operation
234  * -0 or neg value Means an unsuccessful operation
235  */
236  virtual int
237  calcDeltaSolnVariables(const doublereal t,
238  const doublereal* const y,
239  const doublereal* const ydot,
240  doublereal* const delta_y,
241  const doublereal* const solnWeights = 0);
242 
243  //! Returns a vector of column scale factors that can be used to column scale Jacobians.
244  /*!
245  * Default to yScales[] = 1.0
246  *
247  * @param t Time (input)
248  * @param y Solution vector (input, do not modify)
249  * @param y_old Old Solution vector (input, do not modify)
250  * @param yScales Value of the column scales
251  */
252  virtual void calcSolnScales(const doublereal t, const doublereal* const y,
253  const doublereal* const y_old, doublereal* const yScales);
254 
255  //! This function may be used to create output at various points in the execution of an application.
256  /*!
257  *
258  * @param ifunc identity of the call
259  * 0 Initial call
260  * 1 Called at the end of every successful time step
261  * -1 Called at the end of every unsuccessful time step
262  * 2 Called at the end of every call to integrateRJE()
263  *
264  * @param t Time (input)
265  * @param delta_t The current value of the time step (input)
266  * @param y Solution vector (input, do not modify)
267  * @param ydot Rate of change of solution vector. (input)
268  */
269  virtual void user_out2(const int ifunc, const doublereal t,
270  const doublereal delta_t,
271  const doublereal* const y,
272  const doublereal* const ydot);
273 
274  //! This function may be used to create output at various points in the execution of an application.
275  /*!
276  * This routine calls user_out2().
277  *
278  * @param ifunc identity of the call
279  * @param t Time (input)
280  * @param y Solution vector (input, do not modify)
281  * @param ydot Rate of change of solution vector. (input)
282  */
283  virtual void user_out(const int ifunc, const doublereal t,
284  const doublereal* y,
285  const doublereal* ydot);
286 
287  //! Multiply the matrix by another matrix that leads to better conditioning
288  /*!
289  * Provide a left sided matrix that will multiply the current jacobian, after scaling
290  * and lead to a better conditioned system.
291  * This routine is called just before the matrix is factored.
292  *
293  * Original Problem:
294  * J delta_x = - Resid
295  *
296  * New problem:
297  * M (J delta_x) = - M Resid
298  *
299  * @param matrix Pointer to the current jacobian (if zero, it's already been factored)
300  * @param nrows offsets for the matrix
301  * @param rhs residual vector. This also needs to be lhs multiplied by M
302  *
303  * @return Returns a flag to indicate that operation is successful.
304  * 1 Means a successful operation
305  * -0 or neg value Means an unsuccessful operation
306  */
307  virtual int matrixConditioning(doublereal* const matrix, const int nrows,
308  doublereal* const rhs);
309 
310  //! Calculate an analytical jacobian and the residual at the current time and values.
311  /*!
312  * Only called if the jacFormation method is set to analytical
313  *
314  * @param t Time (input)
315  * @param delta_t The current value of the time step (input)
316  * @param cj Coefficient of yprime used in the evaluation of the jacobian
317  * @param y Solution vector (input, do not modify)
318  * @param ydot Rate of change of solution vector. (input, do not modify)
319  * @param J Reference to the SquareMatrix object to be calculated (output)
320  * @param resid Value of the residual that is computed (output)
321  *
322  * @return Returns a flag to indicate that operation is successful.
323  * 1 Means a successful operation
324  * -0 or neg value Means an unsuccessful operation
325  */
326  virtual int evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
327  const doublereal* const y, const doublereal* const ydot,
328  GeneralMatrix& J, doublereal* const resid);
329 
330 
331  //! Calculate an analytical jacobian and the residual at the current time and values.
332  /*!
333  * Only called if the jacFormation method is set to analytical
334  *
335  * @param t Time (input)
336  * @param delta_t The current value of the time step (input)
337  * @param cj Coefficient of yprime used in the evaluation of the jacobian
338  * @param y Solution vector (input, do not modify)
339  * @param ydot Rate of change of solution vector. (input, do not modify)
340  * @param jacobianColPts Pointer to the vector of pts to columns of the SquareMatrix
341  * object to be calculated (output)
342  * @param resid Value of the residual that is computed (output)
343  *
344  * @return Returns a flag to indicate that operation is successful.
345  * 1 Means a successful operation
346  * -0 or neg value Means an unsuccessful operation
347  */
348  virtual int evalJacobianDP(const doublereal t, const doublereal delta_t, doublereal cj,
349  const doublereal* const y,
350  const doublereal* const ydot,
351  doublereal* const* jacobianColPts,
352  doublereal* const resid);
353 
354 protected:
355 
356  //! constant value of atol
357  doublereal m_atol;
358 
359  //! Number of equations
360  int neq_;
361 };
362 }
363 
364 #endif
365