RKF45.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2013-2016 OpenFOAM Foundation
9 Copyright (C) 2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "RKF45.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
34namespace Foam
35{
38
39const scalar
40 RKF45::c2 = 1.0/4.0,
41 RKF45::c3 = 3.0/8.0,
42 RKF45::c4 = 12.0/13.0,
43 RKF45::c5 = 1.0,
44 RKF45::c6 = 1.0/2.0,
45
46 RKF45::a21 = 1.0/4.0,
47 RKF45::a31 = 3.0/32.0,
48 RKF45::a32 = 9.0/32.0,
49 RKF45::a41 = 1932.0/2197.0,
50 RKF45::a42 = -7200.0/2197.0,
51 RKF45::a43 = 7296.0/2197.0,
52 RKF45::a51 = 439.0/216.0,
53 RKF45::a52 = -8.0,
54 RKF45::a53 = 3680.0/513.0,
55 RKF45::a54 = -845.0/4104.0,
56 RKF45::a61 = -8.0/27.0,
57 RKF45::a62 = 2.0,
58 RKF45::a63 = -3544.0/2565.0,
59 RKF45::a64 = 1859.0/4104.0,
60 RKF45::a65 = -11.0/40.0,
61
62 RKF45::b1 = 16.0/135.0,
63 RKF45::b3 = 6656.0/12825.0,
64 RKF45::b4 = 28561.0/56430.0,
65 RKF45::b5 = -9.0/50.0,
66 RKF45::b6 = 2.0/55.0,
67
68 RKF45::e1 = 25.0/216.0 - RKF45::b1,
69 RKF45::e3 = 1408.0/2565.0 - RKF45::b3,
70 RKF45::e4 = 2197.0/4104.0 - RKF45::b4,
71 RKF45::e5 = -1.0/5.0 - RKF45::b5,
72 RKF45::e6 = -RKF45::b6;
73}
74
75
76// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
77
79:
82 yTemp_(n_),
83 k2_(n_),
84 k3_(n_),
85 k4_(n_),
86 k5_(n_),
87 k6_(n_),
88 err_(n_)
89{}
90
91
92// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
93
95{
97 {
99
100 resizeField(yTemp_);
101 resizeField(k2_);
102 resizeField(k3_);
103 resizeField(k4_);
104 resizeField(k5_);
105 resizeField(k6_);
106 resizeField(err_);
107
108 return true;
109 }
110
111 return false;
112}
113
114
116(
117 const scalar x0,
118 const scalarField& y0,
119 const scalarField& dydx0,
120 const scalar dx,
122) const
123{
124 forAll(yTemp_, i)
125 {
126 yTemp_[i] = y0[i] + a21*dx*dydx0[i];
127 }
128
129 odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
130
131 forAll(yTemp_, i)
132 {
133 yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
134 }
135
136 odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
137
138 forAll(yTemp_, i)
139 {
140 yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
141 }
142
143 odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
144
145 forAll(yTemp_, i)
146 {
147 yTemp_[i] = y0[i]
148 + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
149 }
150
151 odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
152
153 forAll(yTemp_, i)
154 {
155 yTemp_[i] = y0[i]
156 + dx
157 *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
158 }
159
160 odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
161
162 // Calculate the 5th-order solution
163 forAll(y, i)
164 {
165 y[i] = y0[i]
166 + dx
167 *(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b5*k5_[i] + b6*k6_[i]);
168 }
169
170 // Calculate the error estimate from the difference between the
171 // 4th-order and 5th-order solutions
172 forAll(err_, i)
173 {
174 err_[i] =
175 dx
176 *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[i]);
177 }
178
179 return normalizeError(y0, y, err_);
180}
181
182
184(
185 scalar& x,
186 scalarField& y,
187 scalar& dxTry
188) const
189{
190 adaptiveSolver::solve(odes_, x, y, dxTry);
191}
192
193
194// ************************************************************************* //
scalar y
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
virtual bool resize()
Resize the ODE solver.
Definition: Euler.C:53
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:57
virtual bool resize()=0
Resize the ODE solver.
Definition: ODESolver.C:92
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:50
4/5th Order Runge-Kutta-Fehlberg ODE solver
Definition: RKF45.H:76
virtual bool resize()
Resize the ODE solver.
Definition: RKF45.C:94
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
An ODE solver for chemistry.
Definition: ode.H:55
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
Namespace for OpenFOAM.
dimensionedScalar y0(const dimensionedScalar &ds)
dictionary dict
CEqn solve()
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333