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-2015 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 
96 Foam::scalar Foam::rodas34::solve
97 (
98  const scalar x0,
99  const scalarField& y0,
100  const scalarField& dydx0,
101  const scalar dx,
102  scalarField& y
103 ) const
104 {
105  odes_.jacobian(x0, y0, dfdx_, dfdy_);
106 
107  for (label i=0; i<n_; i++)
108  {
109  for (label j=0; j<n_; j++)
110  {
111  a_[i][j] = -dfdy_[i][j];
112  }
113 
114  a_[i][i] += 1.0/(gamma*dx);
115  }
116 
117  LUDecompose(a_, pivotIndices_);
118 
119  // Calculate k1:
120  forAll(k1_, i)
121  {
122  k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
123  }
124 
125  LUBacksubstitute(a_, pivotIndices_, k1_);
126 
127  // Calculate k2:
128  forAll(y, i)
129  {
130  y[i] = y0[i] + a21*k1_[i];
131  }
132 
133  odes_.derivatives(x0 + c2*dx, y, dydx_);
134 
135  forAll(k2_, i)
136  {
137  k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
138  }
139 
140  LUBacksubstitute(a_, pivotIndices_, k2_);
141 
142  // Calculate k3:
143  forAll(y, i)
144  {
145  y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
146  }
147 
148  odes_.derivatives(x0 + c3*dx, y, dydx_);
149 
150  forAll(k3_, i)
151  {
152  k3_[i] = dydx_[i] + dx*d3*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
153  }
154 
155  LUBacksubstitute(a_, pivotIndices_, k3_);
156 
157  // Calculate k4:
158  forAll(y, i)
159  {
160  y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
161  }
162 
163  odes_.derivatives(x0 + c4*dx, y, dydx_);
164 
165  forAll(k4_, i)
166  {
167  k4_[i] = dydx_[i] + dx*d4*dfdx_[i]
168  + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
169  }
170 
171  LUBacksubstitute(a_, pivotIndices_, k4_);
172 
173  // Calculate k5:
174  forAll(y, i)
175  {
176  dy_[i] = a51*k1_[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i];
177  y[i] = y0[i] + dy_[i];
178  }
179 
180  odes_.derivatives(x0 + dx, y, dydx_);
181 
182  forAll(k5_, i)
183  {
184  k5_[i] = dydx_[i]
185  + (c51*k1_[i] + c52*k2_[i] + c53*k3_[i] + c54*k4_[i])/dx;
186  }
187 
188  LUBacksubstitute(a_, pivotIndices_, k5_);
189 
190  // Calculate new state and error
191  forAll(y, i)
192  {
193  dy_[i] += k5_[i];
194  y[i] = y0[i] + dy_[i];
195  }
196 
197  odes_.derivatives(x0 + dx, y, dydx_);
198 
199  forAll(err_, i)
200  {
201  err_[i] = dydx_[i]
202  + (c61*k1_[i] + c62*k2_[i] + c63*k3_[i] + c64*k4_[i] + c65*k5_[i])/dx;
203  }
204 
205  LUBacksubstitute(a_, pivotIndices_, err_);
206 
207  forAll(y, i)
208  {
209  y[i] = y0[i] + dy_[i] + err_[i];
210  }
211 
212  return normalizeError(y0, y, err_);
213 }
214 
215 
217 (
218  scalar& x,
219  scalarField& y,
220  scalar& dxTry
221 ) const
222 {
223  adaptiveSolver::solve(odes_, x, y, dxTry);
224 }
225 
226 
227 // ************************************************************************* //
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution.
label n_
Size of the ODESystem.
Definition: ODESolver.H:61
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.
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
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
virtual void derivatives(const scalar x, const scalarField &y, scalarField &dydx) const =0
Calculate the derivatives in dydx.
Namespace for OpenFOAM.
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:46
dictionary dict
An ODE solver for chemistry.
Definition: ode.H:50
#define forAll(list, i)
Definition: UList.H:421
const ODESystem & odes_
Reference to ODESystem.
Definition: ODESolver.H:58
Macros for easy insertion into run-time selection tables.
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:97
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:50
virtual void jacobian(const scalar x, const scalarField &y, scalarField &dfdx, scalarSquareMatrix &dfdy) const =0
Calculate the Jacobian of the system.
rodas34(const ODESystem &ode, const dictionary &dict)
Construct from ODE.
Definition: rodas34.C:75
scalar normalizeError(const scalarField &y0, const scalarField &y, const scalarField &err) const
Return the nomalized scalar error.
Definition: ODESolver.C:67
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)