FreeFOAM The Cross-Platform CFD Toolkit
RK.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) 1991-2010 OpenCFD Ltd.
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 "RK.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
32 
33 namespace Foam
34 {
35  addToRunTimeSelectionTable(ODESolver, RK, ODE);
36 
37 const scalar
38  RK::safety=0.9, RK::pGrow=-0.2, RK::pShrink=-0.25, RK::errCon=1.89e-4;
39 
40 const scalar
41  RK::a2 = 0.2, RK::a3 = 0.3, RK::a4 = 0.6, RK::a5 = 1.0, RK::a6 = 0.875,
42  RK::b21 = 0.2, RK::b31 = 3.0/40.0, RK::b32 = 9.0/40.0,
43  RK::b41 = 0.3, RK::b42 = -0.9, RK::b43 = 1.2,
44  RK::b51 = -11.0/54.0, RK::b52 = 2.5, RK::b53 = -70.0/27.0,
45  RK::b54 = 35.0/27.0,
46  RK::b61 = 1631.0/55296.0, RK::b62 = 175.0/512.0, RK::b63 = 575.0/13824.0,
47  RK::b64 = 44275.0/110592.0, RK::b65 = 253.0/4096.0,
48  RK::c1 = 37.0/378.0, RK::c3 = 250.0/621.0,
49  RK::c4 = 125.0/594.0, RK::c6 = 512.0/1771.0,
50  RK::dc1 = RK::c1 - 2825.0/27648.0, RK::dc3 = RK::c3 - 18575.0/48384.0,
51  RK::dc4 = RK::c4 - 13525.0/55296.0, RK::dc5 = -277.00/14336.0,
52  RK::dc6 = RK::c6 - 0.25;
53 };
54 
55 
56 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
57 
59 :
60  ODESolver(ode),
61  yTemp_(n_, 0.0),
62  ak2_(n_, 0.0),
63  ak3_(n_, 0.0),
64  ak4_(n_, 0.0),
65  ak5_(n_, 0.0),
66  ak6_(n_, 0.0),
67  yErr_(n_, 0.0),
68  yTemp2_(n_, 0.0)
69 {}
70 
71 
72 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
73 
74 void Foam::RK::solve
75 (
76  const ODE& ode,
77  const scalar x,
78  const scalarField& y,
79  const scalarField& dydx,
80  const scalar h,
81  scalarField& yout,
82  scalarField& yerr
83 ) const
84 {
85  forAll(yTemp_, i)
86  {
87  yTemp_[i] = y[i] + b21*h*dydx[i];
88  }
89 
90  ode.derivatives(x + a2*h, yTemp_, ak2_);
91 
92  forAll(yTemp_, i)
93  {
94  yTemp_[i] = y[i] + h*(b31*dydx[i] + b32*ak2_[i]);
95  }
96 
97  ode.derivatives(x + a3*h, yTemp_, ak3_);
98 
99  forAll(yTemp_, i)
100  {
101  yTemp_[i] = y[i] + h*(b41*dydx[i] + b42*ak2_[i] + b43*ak3_[i]);
102  }
103 
104  ode.derivatives(x + a4*h, yTemp_, ak4_);
105 
106  forAll(yTemp_, i)
107  {
108  yTemp_[i] = y[i]
109  + h*(b51*dydx[i] + b52*ak2_[i] + b53*ak3_[i] + b54*ak4_[i]);
110  }
111 
112  ode.derivatives(x + a5*h, yTemp_, ak5_);
113 
114  forAll(yTemp_, i)
115  {
116  yTemp_[i] = y[i]
117  + h*
118  (
119  b61*dydx[i] + b62*ak2_[i] + b63*ak3_[i]
120  + b64*ak4_[i] + b65*ak5_[i]
121  );
122  }
123 
124  ode.derivatives(x + a6*h, yTemp_, ak6_);
125 
126  forAll(yout, i)
127  {
128  yout[i] = y[i]
129  + h*(c1*dydx[i] + c3*ak3_[i] + c4*ak4_[i] + c6*ak6_[i]);
130  }
131 
132  forAll(yerr, i)
133  {
134  yerr[i] =
135  h*
136  (
137  dc1*dydx[i] + dc3*ak3_[i] + dc4*ak4_[i]
138  + dc5*ak5_[i] + dc6*ak6_[i]
139  );
140  }
141 }
142 
143 
144 void Foam::RK::solve
145 (
146  const ODE& ode,
147  scalar& x,
148  scalarField& y,
149  scalarField& dydx,
150  const scalar eps,
151  const scalarField& yScale,
152  const scalar hTry,
153  scalar& hDid,
154  scalar& hNext
155 ) const
156 {
157  scalar h = hTry;
158  scalar maxErr = 0.0;
159 
160  for (;;)
161  {
162  solve(ode, x, y, dydx, h, yTemp2_, yErr_);
163 
164  maxErr = 0.0;
165  for (register label i=0; i<n_; i++)
166  {
167  maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
168  }
169  maxErr /= eps;
170 
171  if (maxErr <= 1.0)
172  {
173  break;
174  }
175 
176  {
177  scalar hTemp = safety*h*pow(maxErr, pShrink);
178  h = (h >= 0.0 ? max(hTemp, 0.1*h) : min(hTemp, 0.1*h));
179  }
180 
181  if (h < VSMALL)
182  {
183  FatalErrorIn("RK::solve")
184  << "stepsize underflow"
185  << exit(FatalError);
186  }
187  }
188 
189  hDid = h;
190 
191  x += h;
192  y = yTemp2_;
193 
194  if (maxErr > errCon)
195  {
196  hNext = safety*h*pow(maxErr, pGrow);
197  }
198  else
199  {
200  hNext = 5.0*h;
201  }
202 }
203 
204 
205 // ************************ vim: set sw=4 sts=4 et: ************************ //