RKDP45.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 "RKDP45.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
35 
36 const scalar
37  RKDP45::c2 = 1.0/5.0,
38  RKDP45::c3 = 3.0/10.0,
39  RKDP45::c4 = 4.0/5.0,
40  RKDP45::c5 = 8.0/9.0,
41 
42  RKDP45::a21 = 1.0/5.0,
43 
44  RKDP45::a31 = 3.0/40.0,
45  RKDP45::a32 = 9.0/40.0,
46 
47  RKDP45::a41 = 44.0/45.0,
48  RKDP45::a42 = -56.0/15.0,
49  RKDP45::a43 = 32.0/9.0,
50 
51  RKDP45::a51 = 19372.0/6561.0,
52  RKDP45::a52 = -25360.0/2187.0,
53  RKDP45::a53 = 64448.0/6561.0,
54  RKDP45::a54 = -212.0/729.0,
55 
56  RKDP45::a61 = 9017.0/3168.0,
57  RKDP45::a62 = -355.0/33.0,
58  RKDP45::a63 = 46732.0/5247.0,
59  RKDP45::a64 = 49.0/176.0,
60  RKDP45::a65 = -5103.0/18656.0,
61 
62  RKDP45::b1 = 35.0/384.0,
63  RKDP45::b3 = 500.0/1113.0,
64  RKDP45::b4 = 125.0/192.0,
65  RKDP45::b5 = -2187.0/6784.0,
66  RKDP45::b6 = 11.0/84.0,
67 
68  RKDP45::e1 = 5179.0/57600.0 - RKDP45::b1,
69  RKDP45::e3 = 7571.0/16695.0 - RKDP45::b3,
70  RKDP45::e4 = 393.0/640.0 - RKDP45::b4,
71  RKDP45::e5 = -92097.0/339200.0 - RKDP45::b5,
72  RKDP45::e6 = 187.0/2100.0 - RKDP45::b6,
73  RKDP45::e7 = 1.0/40.0;
74 }
75 
76 
77 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
78 
80 :
81  ODESolver(ode, dict),
83  yTemp_(n_),
84  k2_(n_),
85  k3_(n_),
86  k4_(n_),
87  k5_(n_),
88  k6_(n_),
89  err_(n_)
90 {}
91 
92 
93 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
94 
96 {
97  if (ODESolver::resize())
98  {
100 
101  resizeField(yTemp_);
102  resizeField(k2_);
103  resizeField(k3_);
104  resizeField(k4_);
105  resizeField(k5_);
106  resizeField(k6_);
107  resizeField(err_);
108 
109  return true;
110  }
111  else
112  {
113  return false;
114  }
115 }
116 
117 
119 (
120  const scalar x0,
121  const scalarField& y0,
122  const label li,
123  const scalarField& dydx0,
124  const scalar dx,
125  scalarField& y
126 ) const
127 {
128  forAll(yTemp_, i)
129  {
130  yTemp_[i] = y0[i] + a21*dx*dydx0[i];
131  }
132 
133  odes_.derivatives(x0 + c2*dx, yTemp_, li, k2_);
134 
135  forAll(yTemp_, i)
136  {
137  yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
138  }
139 
140  odes_.derivatives(x0 + c3*dx, yTemp_, li, k3_);
141 
142  forAll(yTemp_, i)
143  {
144  yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
145  }
146 
147  odes_.derivatives(x0 + c4*dx, yTemp_, li, k4_);
148 
149  forAll(yTemp_, i)
150  {
151  yTemp_[i] = y0[i]
152  + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
153  }
154 
155  odes_.derivatives(x0 + c5*dx, yTemp_, li, k5_);
156 
157  forAll(yTemp_, i)
158  {
159  yTemp_[i] = y0[i]
160  + dx
161  *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
162  }
163 
164  odes_.derivatives(x0 + dx, yTemp_, li, k6_);
165 
166  forAll(y, i)
167  {
168  y[i] = y0[i]
169  + dx*(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b5*k5_[i] + b6*k6_[i]);
170  }
171 
172  // Reuse k2_ for the derivative of the new state
173  odes_.derivatives(x0 + dx, y, li, k2_);
174 
175  forAll(err_, i)
176  {
177  err_[i] =
178  dx
179  *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[i]
180  + e7*k2_[i]);
181  }
182 
183  return normaliseError(y0, y, err_);
184 }
185 
186 
188 (
189  scalar& x,
190  scalarField& y,
191  const label li,
192  scalar& dxTry
193 ) const
194 {
195  adaptiveSolver::solve(odes_, x, y, li, dxTry);
196 }
197 
198 
199 // ************************************************************************* //
scalar y
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:51
virtual bool resize()=0
Resize the ODE solver.
Definition: ODESolver.C:89
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:47
4/5th Order Dormand–Prince Runge-Kutta ODE solver.
Definition: RKDP45.H:69
virtual bool resize()
Resize the ODE solver.
Definition: RKDP45.C:95
virtual void solve(scalar &x, scalarField &y, const label li, scalar &dxTry) const
Inherit solve from ODESolver.
Definition: ODESolver.C:116
RKDP45(const ODESystem &ode, const dictionary &dict)
Construct from ODESystem.
Definition: RKDP45.C:79
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.
bool resize(const label n)
Resize the ODE solver.
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
An ODE solver for chemistry.
Definition: ode.H:53
const dimensionedScalar c2
Second radiation constant: default SI units: [m K].
Namespace for OpenFOAM.
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
addToRunTimeSelectionTable(polyPatch, mergedCyclicPolyPatch, word)
dimensionedScalar y0(const dimensionedScalar &ds)
defineTypeNameAndDebug(combustionModel, 0)
dictionary dict