Cantera  2.2.1
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"
14
15 using namespace std;
16
17 namespace Cantera
18 {
19 ResidJacEval::ResidJacEval(doublereal atol) :
20  m_atol(atol)
21 {
22 }
23
25 {
26  *this = right;
27 }
28
30 {
31  if (this == &right) {
32  return *this;
33  }
34
35  ResidEval::operator=(right);
36
37  m_atol = right.m_atol;
38  neq_ = right.neq_;
39
40  return *this;
41 }
42
44 {
45  return new ResidJacEval(*this);
46 }
47
49 {
50  return neq_;
51 }
52
53 void ResidJacEval::setAtol(doublereal atol)
54 {
55  m_atol = atol;
56  if (m_atol <= 0.0) {
57  throw CanteraError("ResidJacEval::setAtol",
58  "atol must be greater than zero");
59  }
60 }
61
62 int ResidJacEval::getInitialConditions(doublereal t0, doublereal* const y,
63  doublereal* const ydot)
64 {
65  for (int i = 0; i < neq_; i++) {
66  y[i] = 0.0;
67  }
68  if (ydot) {
69  for (int i = 0; i < neq_; i++) {
70  ydot[i] = 0.0;
71  }
72  }
73  return 1;
74 }
75
76 void ResidJacEval::user_out2(const int ifunc, const doublereal t,
77  const doublereal deltaT, const doublereal* y,
78  const doublereal* ydot)
79 {
80 }
81
82 void ResidJacEval::user_out(const int ifunc, const doublereal t,
83  const doublereal* y, const doublereal* ydot)
84 {
85  user_out2(ifunc, t, 0.0, y, ydot);
86 }
87
88 int ResidJacEval::evalTimeTrackingEqns(const doublereal t,
89  const doublereal delta_t,
90  const doublereal* y,
91  const doublereal* ydot)
92 {
93  return 1;
94 }
95
96 int ResidJacEval::calcDeltaSolnVariables(const doublereal t,
97  const doublereal* const ySoln,
98  const doublereal* const ySolnDot,
99  doublereal* const deltaYSoln,
100  const doublereal* const solnWeights)
101 {
102  if (!solnWeights) {
103  for (int i = 0; i < neq_; i++) {
104  deltaYSoln[i] = m_atol + fabs(1.0E-6 * ySoln[i]);
105  }
106  } else {
107  for (int i = 0; i < neq_; i++) {
108  deltaYSoln[i] = std::max(1.0E-2 * solnWeights[i], 1.0E-6 * fabs(ySoln[i]));
109  }
110  }
111  return 1;
112 }
113
114 void ResidJacEval::calcSolnScales(const doublereal t,
115  const doublereal* const ysoln,
116  const doublereal* const ysolnOld,
117  doublereal* const ysolnScales)
118 {
119  if (ysolnScales) {
120  if (ysolnScales[0] == 0.0) {
121  for (int i = 0; i < neq_; i++) {
122  ysolnScales[i] = 1.0;
123  }
124  }
125  }
126 }
127
128 doublereal ResidJacEval::filterNewStep(doublereal t, const doublereal* const ybase, doublereal* const step)
129 {
130  return 0.0;
131 }
132
133 doublereal ResidJacEval::filterSolnPrediction(doublereal t, doublereal* const y)
134 {
135  return 0.0;
136 }
137
138 bool ResidJacEval::evalStoppingCritera(const doublereal t,
139  const doublereal delta_t,
140  const doublereal* const y,
141  const doublereal* const ydot)
142 {
143  return false;
144 }
145
146 int ResidJacEval::matrixConditioning(doublereal* const matrix, const int nrows,
147  doublereal* const rhs)
148 {
149  return 1;
150 }
151
152 int ResidJacEval::evalResidNJ(const doublereal t, const doublereal deltaT,
153  const doublereal* y, const doublereal* ydot,
154  doublereal* const resid,
155  const ResidEval_Type_Enum evalType,
156  const int id_x, const doublereal delta_x)
157 {
158  throw CanteraError("ResidJacEval::evalResidNJ()", "Not implemented\n");
159  return 1;
160 }
161
162 int ResidJacEval::eval(const doublereal t, const doublereal* const y, const doublereal* const ydot,
163  doublereal* const r)
164 {
165  double deltaT = -1.0;
166  return evalResidNJ(t, deltaT, y, ydot, r);
167 }
168
169 int ResidJacEval::evalJacobian(const doublereal t, const doublereal delta_t,
170  doublereal cj, const doublereal* const y,
171  const doublereal* const ydot, GeneralMatrix& J,
172  doublereal* const resid)
173 {
174  doublereal* const* jac_colPts = J.colPts();
175  return evalJacobianDP(t, delta_t, cj, y, ydot, jac_colPts, resid);
176 }
177
178 int ResidJacEval::evalJacobianDP(const doublereal t, const doublereal delta_t,
179  const doublereal c_j,
180  const doublereal* const y,
181  const doublereal* const ydot,
182  doublereal* const* jac_colPts,
183  doublereal* const resid)
184 {
185  throw CanteraError("ResidJacEval::evalJacobianDP()", "Not implemented\n");
186  return 1;
187 }
188
189 }
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:24
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.
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:99
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.