Cantera  2.1.2
ResidJacEval.cpp
Go to the documentation of this file.
1 /**
2  * @file ResidJacEval.cpp
3  */
4 
5 /*
6  * Copyright 2004 Sandia Corporation. Under the terms of Contract
7  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
8  * retains certain rights in this software.
9  * See file License.txt for licensing information.
10  */
11 
12 #include "cantera/base/ct_defs.h"
15 
16 using namespace std;
17 
18 namespace Cantera
19 {
20 ResidJacEval::ResidJacEval(doublereal atol) :
21  ResidEval(),
22  m_atol(atol)
23 {
24 }
25 
27  ResidEval()
28 {
29  *this = operator=(right);
30 }
31 
33 {
34  if (this == &right) {
35  return *this;
36  }
37 
38  ResidEval::operator=(right);
39 
40  m_atol = right.m_atol;
41  neq_ = right.neq_;
42 
43  return *this;
44 }
45 
47 {
48  return new ResidJacEval(*this);
49 }
50 
52 {
53  return neq_;
54 }
55 
56 void ResidJacEval::setAtol(doublereal atol)
57 {
58  m_atol = atol;
59  if (m_atol <= 0.0) {
60  throw CanteraError("ResidJacEval::setAtol",
61  "atol must be greater than zero");
62  }
63 }
64 
66 getInitialConditions(doublereal t0, doublereal* const y, doublereal* const ydot)
67 {
68  for (int i = 0; i < neq_; i++) {
69  y[i] = 0.0;
70  }
71  if (ydot) {
72  for (int i = 0; i < neq_; i++) {
73  ydot[i] = 0.0;
74  }
75  }
76  return 1;
77 }
78 
79 void ResidJacEval::
80 user_out2(const int ifunc, const doublereal t, const doublereal deltaT,
81  const doublereal* y, const doublereal* ydot)
82 {
83 }
84 
85 void ResidJacEval::
86 user_out(const int ifunc, const doublereal t,
87  const doublereal* y, const doublereal* ydot)
88 {
89  user_out2(ifunc, t, 0.0, y, ydot);
90 }
91 
93 evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal* y,
94  const doublereal* ydot)
95 {
96  return 1;
97 }
98 
100 calcDeltaSolnVariables(const doublereal t, const doublereal* const ySoln,
101  const doublereal* const ySolnDot, doublereal* const deltaYSoln,
102  const doublereal* const solnWeights)
103 {
104  if (!solnWeights) {
105  for (int i = 0; i < neq_; i++) {
106  deltaYSoln[i] = m_atol + fabs(1.0E-6 * ySoln[i]);
107  }
108  } else {
109  for (int i = 0; i < neq_; i++) {
110  deltaYSoln[i] = std::max(1.0E-2 * solnWeights[i], 1.0E-6 * fabs(ySoln[i]));
111  }
112  }
113  return 1;
114 }
115 
116 void ResidJacEval::
117 calcSolnScales(const doublereal t, const doublereal* const ysoln, const doublereal* const ysolnOld,
118  doublereal* const ysolnScales)
119 {
120  if (ysolnScales) {
121  if (ysolnScales[0] == 0.0) {
122  for (int i = 0; i < neq_; i++) {
123  ysolnScales[i] = 1.0;
124  }
125  }
126  }
127 }
128 
129 doublereal ResidJacEval::filterNewStep(doublereal t, const doublereal* const ybase, doublereal* const step)
130 {
131  return 0.0;
132 }
133 
134 doublereal ResidJacEval::filterSolnPrediction(doublereal t, doublereal* const y)
135 {
136  return 0.0;
137 }
138 
139 bool ResidJacEval::
140 evalStoppingCritera(const doublereal t,
141  const doublereal delta_t,
142  const doublereal* const y,
143  const doublereal* const ydot)
144 {
145  return false;
146 }
147 
148 int ResidJacEval::
149 matrixConditioning(doublereal* const matrix, const int nrows, doublereal* const rhs)
150 {
151  return 1;
152 }
153 
154 int ResidJacEval::
155 evalResidNJ(const doublereal t, const doublereal deltaT, const doublereal* y,
156  const doublereal* ydot, doublereal* const resid, const ResidEval_Type_Enum evalType,
157  const int id_x, const doublereal delta_x)
158 {
159  throw CanteraError("ResidJacEval::evalResidNJ()", "Not implemented\n");
160  return 1;
161 }
162 
163 int ResidJacEval::eval(const doublereal t, const doublereal* const y, const doublereal* const ydot,
164  doublereal* const r)
165 {
166  double deltaT = -1.0;
167  return evalResidNJ(t, deltaT, y, ydot, r);
168 }
169 
170 int ResidJacEval::
171 evalJacobian(const doublereal t, const doublereal delta_t, doublereal cj,
172  const doublereal* const y,
173  const doublereal* const ydot,
174  GeneralMatrix& J,
175  doublereal* const resid)
176 {
177  doublereal* const* jac_colPts = J.colPts();
178  return evalJacobianDP(t, delta_t, cj, y, ydot, jac_colPts, resid);
179 }
180 
181 int ResidJacEval::
182 evalJacobianDP(const doublereal t, const doublereal delta_t,
183  const doublereal c_j,
184  const doublereal* const y,
185  const doublereal* const ydot,
186  doublereal* const* jac_colPts,
187  doublereal* const resid)
188 {
189  throw CanteraError("ResidJacEval::evalJacobianDP()", "Not implemented\n");
190  return 1;
191 }
192 
193 }
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...
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.
This file contains definitions of terms that are used in internal routines and are unlikely to need m...
Generic matrix.
Definition: GeneralMatrix.h:23
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
Dense, Square (not sparse) matrices.
doublereal m_atol
constant value of atol
Definition: ResidJacEval.h:314
Base class for exceptions thrown by Cantera classes.
Definition: ctexceptions.h:68
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.
virtual doublereal *const * colPts()=0
Return a vector of const pointers to the columns.
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.