rodas34.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "rodas34.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(rodas34, 0);
34  addToRunTimeSelectionTable(ODESolver, rodas34, dictionary);
35 
36 const scalar
37  rodas34::c2 = 0.386,
38  rodas34::c3 = 0.21,
39  rodas34::c4 = 0.63,
40  rodas34::d1 = 0.25,
41  rodas34::d2 = -0.1043,
42  rodas34::d3 = 0.1035,
43  rodas34::d4 = -0.3620000000000023e-01,
44  rodas34::a21 = 0.1544e1,
45  rodas34::a31 = 0.9466785280815826,
46  rodas34::a32 = 0.2557011698983284,
47  rodas34::a41 = 0.3314825187068521e1,
48  rodas34::a42 = 0.2896124015972201e1,
49  rodas34::a43 = 0.9986419139977817,
50  rodas34::a51 = 0.1221224509226641e1,
51  rodas34::a52 = 0.6019134481288629e1,
52  rodas34::a53 = 0.1253708332932087e2,
53  rodas34::a54 = -0.6878860361058950,
54  rodas34::c21 = -0.56688e1,
55  rodas34::c31 = -0.2430093356833875e1,
56  rodas34::c32 = -0.2063599157091915,
57  rodas34::c41 = -0.1073529058151375,
58  rodas34::c42 = -0.9594562251023355e1,
59  rodas34::c43 = -0.2047028614809616e2,
60  rodas34::c51 = 0.7496443313967647e1,
61  rodas34::c52 = -0.1024680431464352e2,
62  rodas34::c53 = -0.3399990352819905e2,
63  rodas34::c54 = 0.1170890893206160e2,
64  rodas34::c61 = 0.8083246795921522e1,
65  rodas34::c62 = -0.7981132988064893e1,
66  rodas34::c63 = -0.3152159432874371e2,
67  rodas34::c64 = 0.1631930543123136e2,
68  rodas34::c65 = -0.6058818238834054e1,
69  rodas34::gamma = 0.25;
70 }
71 
72 
73 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
74 
76 :
77  ODESolver(ode, dict),
78  adaptiveSolver(ode, dict),
79  k1_(n_),
80  k2_(n_),
81  k3_(n_),
82  k4_(n_),
83  k5_(n_),
84  dy_(n_),
85  err_(n_),
86  dydx_(n_),
87  dfdx_(n_),
88  dfdy_(n_, n_),
89  a_(n_, n_),
90  pivotIndices_(n_)
91 {}
92 
93 
94 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
95 
97 {
98  if (ODESolver::resize())
99  {
101 
102  resizeField(k1_);
103  resizeField(k2_);
104  resizeField(k3_);
105  resizeField(k4_);
106  resizeField(k5_);
107  resizeField(dy_);
108  resizeField(err_);
109  resizeField(dydx_);
110  resizeField(dfdx_);
111  resizeMatrix(dfdy_);
112  resizeMatrix(a_);
113  resizeField(pivotIndices_);
114 
115  return true;
116  }
117  else
118  {
119  return false;
120  }
121 }
122 
123 
124 Foam::scalar Foam::rodas34::solve
125 (
126  const scalar x0,
127  const scalarField& y0,
128  const scalarField& dydx0,
129  const scalar dx,
130  scalarField& y
131 ) const
132 {
133  odes_.jacobian(x0, y0, dfdx_, dfdy_);
134 
135  for (label i=0; i<n_; i++)
136  {
137  for (label j=0; j<n_; j++)
138  {
139  a_(i, j) = -dfdy_(i, j);
140  }
141 
142  a_(i, i) += 1.0/(gamma*dx);
143  }
144 
145  LUDecompose(a_, pivotIndices_);
146 
147  // Calculate k1:
148  forAll(k1_, i)
149  {
150  k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
151  }
152 
153  LUBacksubstitute(a_, pivotIndices_, k1_);
154 
155  // Calculate k2:
156  forAll(y, i)
157  {
158  y[i] = y0[i] + a21*k1_[i];
159  }
160 
161  odes_.derivatives(x0 + c2*dx, y, dydx_);
162 
163  forAll(k2_, i)
164  {
165  k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
166  }
167 
168  LUBacksubstitute(a_, pivotIndices_, k2_);
169 
170  // Calculate k3:
171  forAll(y, i)
172  {
173  y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
174  }
175 
176  odes_.derivatives(x0 + c3*dx, y, dydx_);
177 
178  forAll(k3_, i)
179  {
180  k3_[i] = dydx_[i] + dx*d3*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
181  }
182 
183  LUBacksubstitute(a_, pivotIndices_, k3_);
184 
185  // Calculate k4:
186  forAll(y, i)
187  {
188  y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
189  }
190 
191  odes_.derivatives(x0 + c4*dx, y, dydx_);
192 
193  forAll(k4_, i)
194  {
195  k4_[i] = dydx_[i] + dx*d4*dfdx_[i]
196  + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
197  }
198 
199  LUBacksubstitute(a_, pivotIndices_, k4_);
200 
201  // Calculate k5:
202  forAll(y, i)
203  {
204  dy_[i] = a51*k1_[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i];
205  y[i] = y0[i] + dy_[i];
206  }
207 
208  odes_.derivatives(x0 + dx, y, dydx_);
209 
210  forAll(k5_, i)
211  {
212  k5_[i] = dydx_[i]
213  + (c51*k1_[i] + c52*k2_[i] + c53*k3_[i] + c54*k4_[i])/dx;
214  }
215 
216  LUBacksubstitute(a_, pivotIndices_, k5_);
217 
218  // Calculate new state and error
219  forAll(y, i)
220  {
221  dy_[i] += k5_[i];
222  y[i] = y0[i] + dy_[i];
223  }
224 
225  odes_.derivatives(x0 + dx, y, dydx_);
226 
227  forAll(err_, i)
228  {
229  err_[i] = dydx_[i]
230  + (c61*k1_[i] + c62*k2_[i] + c63*k3_[i] + c64*k4_[i] + c65*k5_[i])/dx;
231  }
232 
233  LUBacksubstitute(a_, pivotIndices_, err_);
234 
235  forAll(y, i)
236  {
237  y[i] = y0[i] + dy_[i] + err_[i];
238  }
239 
240  return normalizeError(y0, y, err_);
241 }
242 
243 
245 (
246  scalar& x,
247  scalarField& y,
248  scalar& dxTry
249 ) const
250 {
251  adaptiveSolver::solve(odes_, x, y, dxTry);
252 }
253 
254 
255 // ************************************************************************* //
virtual scalar solve(const scalar x0, const scalarField &y0, const scalarField &dydx0, const scalar dx, scalarField &y) const =0
Solve a single step dx and return the error.
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:46
virtual bool resize()=0
Resize the ODE solver.
Definition: ODESolver.C:89
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
An ODE solver for chemistry.
Definition: ode.H:50
Macros for easy insertion into run-time selection tables.
virtual void jacobian(const scalar x, const scalarField &y, scalarField &dfdx, scalarSquareMatrix &dfdy) const =0
Calculate the Jacobian of the system.
bool resize(const label n)
Resize the ODE solver.
void resizeMatrix(scalarSquareMatrix &m) const
Definition: ODESolverI.H:61
const ODESystem & odes_
Reference to ODESystem.
Definition: ODESolver.H:58
rodas34(const ODESystem &ode, const dictionary &dict)
Construct from ODESystem.
Definition: rodas34.C:75
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
scalar normalizeError(const scalarField &y0, const scalarField &y, const scalarField &err) const
Return the nomalized scalar error.
Definition: ODESolver.C:40
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:50
static void resizeField(UList< Type > &f, const label n)
Definition: ODESolverI.H:48
virtual scalar solve(const scalar x0, const scalarField &y0, const scalarField &dydx0, const scalar dx, scalarField &y) const
Solve a single step dx and return the error.
Definition: rodas34.C:125
label n_
Size of the ODESystem (adjustable)
Definition: ODESolver.H:64
virtual bool resize()
Resize the ODE solver.
Definition: rodas34.C:96
virtual void derivatives(const scalar x, const scalarField &y, scalarField &dydx) const =0
Calculate the derivatives in dydx.
Namespace for OpenFOAM.
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution.