FreeFOAM The Cross-Platform CFD Toolkit
SIBS.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 "SIBS.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
32 
33 namespace Foam
34 {
35  addToRunTimeSelectionTable(ODESolver, SIBS, ODE);
36 
37  const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70};
38 
39  const scalar
40  SIBS::safe1 = 0.25,
41  SIBS::safe2 = 0.7,
42  SIBS::redMax = 1.0e-5,
43  SIBS::redMin = 0.7,
44  SIBS::scaleMX = 0.1;
45 };
46 
47 
48 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
49 
51 :
52  ODESolver(ode),
53  a_(iMaxX_, 0.0),
54  alpha_(kMaxX_, kMaxX_, 0.0),
55  d_p_(n_, kMaxX_, 0.0),
56  x_p_(kMaxX_, 0.0),
57  err_(kMaxX_, 0.0),
58 
59  yTemp_(n_, 0.0),
60  ySeq_(n_, 0.0),
61  yErr_(n_, 0.0),
62  dfdx_(n_, 0.0),
63  dfdy_(n_, n_, 0.0),
64  first_(1),
65  epsOld_(-1.0)
66 {}
67 
68 
69 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
70 
72 (
73  const ODE& ode,
74  scalar& x,
75  scalarField& y,
76  scalarField& dydx,
77  const scalar eps,
78  const scalarField& yScale,
79  const scalar hTry,
80  scalar& hDid,
81  scalar& hNext
82 ) const
83 {
84  bool exitflag = false;
85 
86  if (eps != epsOld_)
87  {
88  hNext = xNew_ = -GREAT;
89  scalar eps1 = safe1*eps;
90  a_[0] = nSeq_[0] + 1;
91 
92  for (register label k=0; k<kMaxX_; k++)
93  {
94  a_[k + 1] = a_[k] + nSeq_[k + 1];
95  }
96 
97  for (register label iq = 1; iq<kMaxX_; iq++)
98  {
99  for (register label k=0; k<iq; k++)
100  {
101  alpha_[k][iq] =
102  pow(eps1, (a_[k + 1] - a_[iq + 1])
103  /((a_[iq + 1] - a_[0] + 1.0)*(2*k + 3)));
104  }
105  }
106 
107  epsOld_ = eps;
108  a_[0] += n_;
109 
110  for (register label k=0; k<kMaxX_; k++)
111  {
112  a_[k + 1] = a_[k] + nSeq_[k + 1];
113  }
114 
115  for (kOpt_ = 1; kOpt_<kMaxX_ - 1; kOpt_++)
116  {
117  if (a_[kOpt_ + 1] > a_[kOpt_]*alpha_[kOpt_ - 1][kOpt_])
118  {
119  break;
120  }
121  }
122 
123  kMax_ = kOpt_;
124  }
125 
126  label k=0;
127  scalar h = hTry;
128  yTemp_ = y;
129 
130  ode.jacobian(x, y, dfdx_, dfdy_);
131 
132  if (x != xNew_ || h != hNext)
133  {
134  first_ = 1;
135  kOpt_ = kMax_;
136  }
137 
138  label km=0;
139  label reduct=0;
140  scalar maxErr = SMALL;
141 
142  for (;;)
143  {
144  scalar red = 1.0;
145 
146  for (k=0; k <= kMax_; k++)
147  {
148  xNew_ = x + h;
149 
150  if (xNew_ == x)
151  {
152  FatalErrorIn("ODES::SIBS")
153  << "step size underflow"
154  << exit(FatalError);
155  }
156 
157  SIMPR(ode, x, yTemp_, dydx, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
158  scalar xest = sqr(h/nSeq_[k]);
159 
160  polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);
161 
162  if (k != 0)
163  {
164  maxErr = SMALL;
165  for (register label i=0; i<n_; i++)
166  {
167  maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
168  }
169  maxErr /= eps;
170  km = k - 1;
171  err_[km] = pow(maxErr/safe1, 1.0/(2*km + 3));
172  }
173 
174  if (k != 0 && (k >= kOpt_ - 1 || first_))
175  {
176  if (maxErr < 1.0)
177  {
178  exitflag = true;
179  break;
180  }
181 
182  if (k == kMax_ || k == kOpt_ + 1)
183  {
184  red = safe2/err_[km];
185  break;
186  }
187  else if (k == kOpt_ && alpha_[kOpt_ - 1][kOpt_] < err_[km])
188  {
189  red = 1.0/err_[km];
190  break;
191  }
192  else if (kOpt_ == kMax_ && alpha_[km][kMax_ - 1] < err_[km])
193  {
194  red = alpha_[km][kMax_ - 1]*safe2/err_[km];
195  break;
196  }
197  else if (alpha_[km][kOpt_] < err_[km])
198  {
199  red = alpha_[km][kOpt_ - 1]/err_[km];
200  break;
201  }
202  }
203  }
204 
205  if (exitflag)
206  {
207  break;
208  }
209 
210  red = min(red, redMin);
211  red = max(red, redMax);
212  h *= red;
213  reduct = 1;
214  }
215 
216  x = xNew_;
217  hDid = h;
218  first_=0;
219  scalar wrkmin = GREAT;
220  scalar scale = 1.0;
221 
222  for (register label kk=0; kk<=km; kk++)
223  {
224  scalar fact = max(err_[kk], scaleMX);
225  scalar work = fact*a_[kk + 1];
226  if (work < wrkmin)
227  {
228  scale = fact;
229  wrkmin = work;
230  kOpt_ = kk + 1;
231  }
232  }
233 
234  hNext = h/scale;
235 
236  if (kOpt_ >= k && kOpt_ != kMax_ && !reduct)
237  {
238  scalar fact = max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX);
239  if (a_[kOpt_ + 1]*fact <= wrkmin)
240  {
241  hNext = h/fact;
242  kOpt_++;
243  }
244  }
245 }
246 
247 
248 // ************************ vim: set sw=4 sts=4 et: ************************ //