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