rodas34.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2021 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 label li,
129  const scalarField& dydx0,
130  const scalar dx,
131  scalarField& y
132 ) const
133 {
134  odes_.jacobian(x0, y0, li, dfdx_, dfdy_);
135 
136  for (label i=0; i<n_; i++)
137  {
138  for (label j=0; j<n_; j++)
139  {
140  a_(i, j) = -dfdy_(i, j);
141  }
142 
143  a_(i, i) += 1.0/(gamma*dx);
144  }
145 
146  LUDecompose(a_, pivotIndices_);
147 
148  // Calculate k1:
149  forAll(k1_, i)
150  {
151  k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
152  }
153 
154  LUBacksubstitute(a_, pivotIndices_, k1_);
155 
156  // Calculate k2:
157  forAll(y, i)
158  {
159  y[i] = y0[i] + a21*k1_[i];
160  }
161 
162  odes_.derivatives(x0 + c2*dx, y, li, dydx_);
163 
164  forAll(k2_, i)
165  {
166  k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
167  }
168 
169  LUBacksubstitute(a_, pivotIndices_, k2_);
170 
171  // Calculate k3:
172  forAll(y, i)
173  {
174  y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
175  }
176 
177  odes_.derivatives(x0 + c3*dx, y, li, dydx_);
178 
179  forAll(k3_, i)
180  {
181  k3_[i] = dydx_[i] + dx*d3*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
182  }
183 
184  LUBacksubstitute(a_, pivotIndices_, k3_);
185 
186  // Calculate k4:
187  forAll(y, i)
188  {
189  y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
190  }
191 
192  odes_.derivatives(x0 + c4*dx, y, li, dydx_);
193 
194  forAll(k4_, i)
195  {
196  k4_[i] = dydx_[i] + dx*d4*dfdx_[i]
197  + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
198  }
199 
200  LUBacksubstitute(a_, pivotIndices_, k4_);
201 
202  // Calculate k5:
203  forAll(y, i)
204  {
205  dy_[i] = a51*k1_[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i];
206  y[i] = y0[i] + dy_[i];
207  }
208 
209  odes_.derivatives(x0 + dx, y, li, dydx_);
210 
211  forAll(k5_, i)
212  {
213  k5_[i] = dydx_[i]
214  + (c51*k1_[i] + c52*k2_[i] + c53*k3_[i] + c54*k4_[i])/dx;
215  }
216 
217  LUBacksubstitute(a_, pivotIndices_, k5_);
218 
219  // Calculate new state and error
220  forAll(y, i)
221  {
222  dy_[i] += k5_[i];
223  y[i] = y0[i] + dy_[i];
224  }
225 
226  odes_.derivatives(x0 + dx, y, li, dydx_);
227 
228  forAll(err_, i)
229  {
230  err_[i] = dydx_[i]
231  + (c61*k1_[i] + c62*k2_[i] + c63*k3_[i] + c64*k4_[i] + c65*k5_[i])/dx;
232  }
233 
234  LUBacksubstitute(a_, pivotIndices_, err_);
235 
236  forAll(y, i)
237  {
238  y[i] = y0[i] + dy_[i] + err_[i];
239  }
240 
241  return normaliseError(y0, y, err_);
242 }
243 
244 
246 (
247  scalar& x,
248  scalarField& y,
249  const label li,
250  scalar& dxTry
251 ) const
252 {
253  adaptiveSolver::solve(odes_, x, y, li, dxTry);
254 }
255 
256 
257 // ************************************************************************* //
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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:156
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
virtual void derivatives(const scalar x, const scalarField &y, const label li, scalarField &dydx) const =0
Calculate the derivatives in dydx.
An ODE solver for chemistry.
Definition: ode.H:50
scalar normaliseError(const scalarField &y0, const scalarField &y, const scalarField &err) const
Return the normalised scalar error.
Definition: ODESolver.C:40
Macros for easy insertion into run-time selection tables.
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)
virtual scalar solve(const scalar x0, const scalarField &y0, const label li, const scalarField &dydx0, const scalar dx, scalarField &y) const
Solve a single step dx and return the error.
Definition: rodas34.C:125
defineTypeNameAndDebug(combustionModel, 0)
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:50
virtual void jacobian(const scalar x, const scalarField &y, const label li, scalarField &dfdx, scalarSquareMatrix &dfdy) const =0
Calculate the Jacobian of the system.
static void resizeField(UList< Type > &f, const label n)
Definition: ODESolverI.H:48
virtual scalar solve(const scalar x0, const scalarField &y0, const label li, const scalarField &dydx0, const scalar dx, scalarField &y) const =0
Solve a single step dx and return the error.
label n_
Size of the ODESystem (adjustable)
Definition: ODESolver.H:64
virtual bool resize()
Resize the ODE solver.
Definition: rodas34.C:96
Namespace for OpenFOAM.
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution.