SIBS.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) 2011-2019 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 "SIBS.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(SIBS, 0);
34  addToRunTimeSelectionTable(ODESolver, SIBS, dictionary);
35 
36  const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70};
37 
38  const scalar
39  SIBS::safe1 = 0.25,
40  SIBS::safe2 = 0.7,
41  SIBS::redMax = 1.0e-5,
42  SIBS::redMin = 0.7,
43  SIBS::scaleMX = 0.1;
44 }
45 
46 
47 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
48 
50 :
51  ODESolver(ode, dict),
52  a_(iMaxX_, 0.0),
53  alpha_(kMaxX_, 0.0),
54  d_p_(n_, kMaxX_, 0.0),
55  x_p_(kMaxX_, 0.0),
56  err_(kMaxX_, 0.0),
57 
58  yTemp_(n_, 0.0),
59  ySeq_(n_, 0.0),
60  yErr_(n_, 0.0),
61  dydx0_(n_),
62  dfdx_(n_, 0.0),
63  dfdy_(n_, 0.0),
64  first_(1),
65  epsOld_(-1.0)
66 {}
67 
68 
69 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
70 
72 {
73  if (ODESolver::resize())
74  {
75  resizeField(yTemp_);
76  resizeField(ySeq_);
77  resizeField(yErr_);
78  resizeField(dydx0_);
79  resizeField(dfdx_);
80  resizeMatrix(dfdy_);
81 
82  return true;
83  }
84  else
85  {
86  return false;
87  }
88 }
89 
90 
92 (
93  scalar& x,
94  scalarField& y,
95  const label li,
96  scalar& dxTry
97 ) const
98 {
99  odes_.derivatives(x, y, li, dydx0_);
100 
101  scalar h = dxTry;
102  bool exitflag = false;
103 
104  if (relTol_[0] != epsOld_)
105  {
106  dxTry = xNew_ = -great;
107  scalar eps1 = safe1*relTol_[0];
108  a_[0] = nSeq_[0] + 1;
109 
110  for (label k=0; k<kMaxX_; k++)
111  {
112  a_[k + 1] = a_[k] + nSeq_[k + 1];
113  }
114 
115  for (label iq = 1; iq<kMaxX_; iq++)
116  {
117  for (label k=0; k<iq; k++)
118  {
119  alpha_[k][iq] =
120  pow(eps1, (a_[k + 1] - a_[iq + 1])
121  /((a_[iq + 1] - a_[0] + 1.0)*(2*k + 3)));
122  }
123  }
124 
125  epsOld_ = relTol_[0];
126  a_[0] += n_;
127 
128  for (label k=0; k<kMaxX_; k++)
129  {
130  a_[k + 1] = a_[k] + nSeq_[k + 1];
131  }
132 
133  for (kOpt_ = 1; kOpt_<kMaxX_ - 1; kOpt_++)
134  {
135  if (a_[kOpt_ + 1] > a_[kOpt_]*alpha_[kOpt_ - 1][kOpt_])
136  {
137  break;
138  }
139  }
140 
141  kMax_ = kOpt_;
142  }
143 
144  label k = 0;
145  yTemp_ = y;
146 
147  odes_.jacobian(x, y, li, dfdx_, dfdy_);
148 
149  if (x != xNew_ || h != dxTry)
150  {
151  first_ = 1;
152  kOpt_ = kMax_;
153  }
154 
155  label km=0;
156  label reduct=0;
157  scalar maxErr = small;
158 
159  for (;;)
160  {
161  scalar red = 1.0;
162 
163  for (k=0; k <= kMax_; k++)
164  {
165  xNew_ = x + h;
166 
167  if (xNew_ == x)
168  {
170  << "step size underflow"
171  << exit(FatalError);
172  }
173 
174  SIMPR(x, yTemp_, li, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
175  scalar xest = sqr(h/nSeq_[k]);
176 
177  polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);
178 
179  if (k != 0)
180  {
181  maxErr = small;
182  for (label i=0; i<n_; i++)
183  {
184  maxErr = max
185  (
186  maxErr,
187  mag(yErr_[i])/(absTol_[i] + relTol_[i]*mag(yTemp_[i]))
188  );
189  }
190  km = k - 1;
191  err_[km] = pow(maxErr/safe1, 1.0/(2*km + 3));
192  }
193 
194  if (k != 0 && (k >= kOpt_ - 1 || first_))
195  {
196  if (maxErr < 1.0)
197  {
198  exitflag = true;
199  break;
200  }
201 
202  if (k == kMax_ || k == kOpt_ + 1)
203  {
204  red = safe2/err_[km];
205  break;
206  }
207  else if (k == kOpt_ && alpha_[kOpt_ - 1][kOpt_] < err_[km])
208  {
209  red = 1.0/err_[km];
210  break;
211  }
212  else if (kOpt_ == kMax_ && alpha_[km][kMax_ - 1] < err_[km])
213  {
214  red = alpha_[km][kMax_ - 1]*safe2/err_[km];
215  break;
216  }
217  else if (alpha_[km][kOpt_] < err_[km])
218  {
219  red = alpha_[km][kOpt_ - 1]/err_[km];
220  break;
221  }
222  }
223  }
224 
225  if (exitflag)
226  {
227  break;
228  }
229 
230  red = min(red, redMin);
231  red = max(red, redMax);
232  h *= red;
233  reduct = 1;
234  }
235 
236  x = xNew_;
237  first_ = 0;
238  scalar wrkmin = great;
239  scalar scale = 1.0;
240 
241  for (label kk=0; kk<=km; kk++)
242  {
243  scalar fact = max(err_[kk], scaleMX);
244  scalar work = fact*a_[kk + 1];
245  if (work < wrkmin)
246  {
247  scale = fact;
248  wrkmin = work;
249  kOpt_ = kk + 1;
250  }
251  }
252 
253  dxTry = h/scale;
254 
255  if (kOpt_ >= k && kOpt_ != kMax_ && !reduct)
256  {
257  scalar fact = max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX);
258  if (a_[kOpt_ + 1]*fact <= wrkmin)
259  {
260  dxTry = h/fact;
261  kOpt_++;
262  }
263  }
264 }
265 
266 
267 // ************************************************************************* //
dictionary dict
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:46
virtual bool resize()=0
Resize the ODE solver.
Definition: ODESolver.C:89
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
dimensionedSymmTensor sqr(const dimensionedVector &dv)
virtual void derivatives(const scalar x, const scalarField &y, const label li, scalarField &dydx) const =0
Calculate the derivatives in dydx.
label k
Boltzmann constant.
An ODE solver for chemistry.
Definition: ode.H:50
const dimensionedScalar h
Planck constant.
Macros for easy insertion into run-time selection tables.
scalar y
void resizeMatrix(scalarSquareMatrix &m) const
Definition: ODESolverI.H:61
const ODESystem & odes_
Reference to ODESystem.
Definition: ODESolver.H:58
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
virtual void solve(scalar &x, scalarField &y, const label li, scalar &dxTry) const
Solve the ODE system from the current state xStart, y.
Definition: SIBS.C:92
SIBS(const ODESystem &ode, const dictionary &dict)
Construct from ODE system.
Definition: SIBS.C:49
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
scalarField relTol_
Relative convergence tolerance per step.
Definition: ODESolver.H:70
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:50
virtual void jacobian(const scalar x, const scalarField &y, const label li, scalarField &dfdx, scalarSquareMatrix &dfdy) const =0
Calculate the Jacobian of the system.
scalarField absTol_
Absolute convergence tolerance per step.
Definition: ODESolver.H:67
dimensioned< scalar > mag(const dimensioned< Type > &)
static void resizeField(UList< Type > &f, const label n)
Definition: ODESolverI.H:48
virtual bool resize()
Resize the ODE solver.
Definition: SIBS.C:71
label n_
Size of the ODESystem (adjustable)
Definition: ODESolver.H:64
Namespace for OpenFOAM.