FreeFOAM The Cross-Platform CFD Toolkit
KRR4.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 "KRR4.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
32 
33 namespace Foam
34 {
35  addToRunTimeSelectionTable(ODESolver, KRR4, ODE);
36 
37 const scalar
38  KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25,
39  KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296;
40 
41 const scalar
42  KRR4::gamma = 1.0/2.0,
43  KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
44  KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0,
45  KRR4::c41 = -112.0/125.0, KRR4::c42 = -54.0/125.0, KRR4::c43 = -2.0/5.0,
46  KRR4::b1 = 19.0/9.0, KRR4::b2 = 1.0/2.0, KRR4::b3 = 25.0/108.0,
47  KRR4::b4 = 125.0/108.0,
48  KRR4::e1 = 17.0/54.0, KRR4::e2 = 7.0/36.0, KRR4::e3 = 0.0,
49  KRR4::e4 = 125.0/108.0,
50  KRR4::c1X = 1.0/2.0, KRR4::c2X = -3.0/2.0, KRR4::c3X = 121.0/50.0,
51  KRR4::c4X = 29.0/250.0,
52  KRR4::a2X = 1.0, KRR4::a3X = 3.0/5.0;
53 };
54 
55 
56 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
57 
59 :
60  ODESolver(ode),
61  yTemp_(n_, 0.0),
62  dydxTemp_(n_, 0.0),
63  g1_(n_, 0.0),
64  g2_(n_, 0.0),
65  g3_(n_, 0.0),
66  g4_(n_, 0.0),
67  yErr_(n_, 0.0),
68  dfdx_(n_, 0.0),
69  dfdy_(n_, n_, 0.0),
70  a_(n_, n_, 0.0),
71  pivotIndices_(n_, 0.0)
72 {}
73 
74 
75 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
76 
78 (
79  const ODE& ode,
80  scalar& x,
81  scalarField& y,
82  scalarField& dydx,
83  const scalar eps,
84  const scalarField& yScale,
85  const scalar hTry,
86  scalar& hDid,
87  scalar& hNext
88 ) const
89 {
90  scalar xTemp = x;
91  yTemp_ = y;
92  dydxTemp_ = dydx;
93 
94  ode.jacobian(xTemp, yTemp_, dfdx_, dfdy_);
95 
96  scalar h = hTry;
97 
98  for (register label jtry=0; jtry<maxtry; jtry++)
99  {
100  for (register label i=0; i<n_; i++)
101  {
102  for (register label j=0; j<n_; j++)
103  {
104  a_[i][j] = -dfdy_[i][j];
105  }
106 
107  a_[i][i] += 1.0/(gamma*h);
108  }
109 
110  LUDecompose(a_, pivotIndices_);
111 
112  for (register label i=0; i<n_; i++)
113  {
114  g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];
115  }
116 
117  LUBacksubstitute(a_, pivotIndices_, g1_);
118 
119  for (register label i=0; i<n_; i++)
120  {
121  y[i] = yTemp_[i] + a21*g1_[i];
122  }
123 
124  x = xTemp + a2X*h;
125  ode.derivatives(x, y, dydx_);
126 
127  for (register label i=0; i<n_; i++)
128  {
129  g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/h;
130  }
131 
132  LUBacksubstitute(a_, pivotIndices_, g2_);
133 
134  for (register label i=0; i<n_; i++)
135  {
136  y[i] = yTemp_[i] + a31*g1_[i] + a32*g2_[i];
137  }
138 
139  x = xTemp + a3X*h;
140  ode.derivatives(x, y, dydx_);
141 
142  for (register label i=0; i<n_; i++)
143  {
144  g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;
145  }
146 
147  LUBacksubstitute(a_, pivotIndices_, g3_);
148 
149  for (register label i=0; i<n_; i++)
150  {
151  g4_[i] = dydx_[i] + h*c4X*dfdx_[i]
152  + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;
153  }
154 
155  LUBacksubstitute(a_, pivotIndices_, g4_);
156 
157  for (register label i=0; i<n_; i++)
158  {
159  y[i] = yTemp_[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];
160  yErr_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];
161  }
162 
163  x = xTemp + h;
164 
165  if (x == xTemp)
166  {
167  FatalErrorIn("ODES::KRR4")
168  << "stepsize not significant"
169  << exit(FatalError);
170  }
171 
172  scalar maxErr = 0.0;
173  for (register label i=0; i<n_; i++)
174  {
175  maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
176  }
177  maxErr /= eps;
178 
179  if (maxErr <= 1.0)
180  {
181  hDid = h;
182  hNext = (maxErr > errcon ? safety*h*pow(maxErr, pgrow) : grow*h);
183  return;
184  }
185  else
186  {
187  hNext = safety*h*pow(maxErr, pshrink);
188  h = (h >= 0.0 ? max(hNext, shrink*h) : min(hNext, shrink*h));
189  }
190  }
191 
192  FatalErrorIn("ODES::KRR4")
193  << "exceeded maxtry"
194  << exit(FatalError);
195 }
196 
197 
198 // ************************ vim: set sw=4 sts=4 et: ************************ //