ODESolver.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-2018 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 "ODESolver.H"
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
30 namespace Foam
31 {
32  defineTypeNameAndDebug(ODESolver, 0);
33  defineRunTimeSelectionTable(ODESolver, dictionary);
34 }
35 
36 
37 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
38 
40 (
41  const scalarField& y0,
42  const scalarField& y,
43  const scalarField& err
44 ) const
45 {
46  // Calculate the maximum error
47  scalar maxErr = 0.0;
48  forAll(err, i)
49  {
50  scalar tol = absTol_[i] + relTol_[i]*max(mag(y0[i]), mag(y[i]));
51  maxErr = max(maxErr, mag(err[i])/tol);
52  }
53 
54  return maxErr;
55 }
56 
57 
58 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
59 
61 :
62  odes_(ode),
63  maxN_(ode.nEqns()),
64  n_(ode.nEqns()),
65  absTol_(n_, dict.lookupOrDefault<scalar>("absTol", small)),
66  relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-4)),
67  maxSteps_(dict.lookupOrDefault<scalar>("maxSteps", 10000))
68 {}
69 
70 
72 (
73  const ODESystem& ode,
74  const scalarField& absTol,
75  const scalarField& relTol
76 )
77 :
78  odes_(ode),
79  maxN_(ode.nEqns()),
80  n_(ode.nEqns()),
81  absTol_(absTol),
82  relTol_(relTol),
83  maxSteps_(10000)
84 {}
85 
86 
87 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
88 
90 {
91  if (odes_.nEqns() != n_)
92  {
93  if (odes_.nEqns() > maxN_)
94  {
96  << "Specified number of equations " << odes_.nEqns()
97  << " greater than maximum " << maxN_
98  << abort(FatalError);
99  }
100 
101  n_ = odes_.nEqns();
102 
105 
106  return true;
107  }
108  else
109  {
110  return false;
111  }
112 }
113 
114 
116 (
117  scalar& x,
118  scalarField& y,
119  scalar& dxTry
120 ) const
121 {
122  stepState step(dxTry);
123  solve(x, y, step);
124  dxTry = step.dxTry;
125 }
126 
127 
129 (
130  scalar& x,
131  scalarField& y,
132  stepState& step
133 ) const
134 {
135  scalar x0 = x;
136  solve(x, y, step.dxTry);
137  step.dxDid = x - x0;
138 }
139 
140 
142 (
143  const scalar xStart,
144  const scalar xEnd,
145  scalarField& y,
146  scalar& dxTry
147 ) const
148 {
149  stepState step(dxTry);
150  scalar x = xStart;
151 
152  for (label nStep=0; nStep<maxSteps_; nStep++)
153  {
154  // Store previous iteration dxTry
155  scalar dxTry0 = step.dxTry;
156 
157  step.reject = false;
158 
159  // Check if this is a truncated step and set dxTry to integrate to xEnd
160  if ((x + step.dxTry - xEnd)*(x + step.dxTry - xStart) > 0)
161  {
162  step.last = true;
163  step.dxTry = xEnd - x;
164  }
165 
166  // Integrate as far as possible up to step.dxTry
167  solve(x, y, step);
168 
169  // Check if reached xEnd
170  if ((x - xEnd)*(xEnd - xStart) >= 0)
171  {
172  if (nStep > 0 && step.last)
173  {
174  step.dxTry = dxTry0;
175  }
176 
177  dxTry = step.dxTry;
178 
179  return;
180  }
181 
182  step.first = false;
183 
184  // If the step.dxTry was reject set step.prevReject
185  if (step.reject)
186  {
187  step.prevReject = true;
188  }
189  }
190 
192  << "Integration steps greater than maximum " << maxSteps_ << nl
193  << " xStart = " << xStart << ", xEnd = " << xEnd
194  << ", x = " << x << ", dxDid = " << step.dxDid << nl
195  << " y = " << y
196  << exit(FatalError);
197 }
198 
199 
200 // ************************************************************************* //
virtual void solve(scalar &x, scalarField &y, scalar &dxTry) const
Solve the ODE system as far as possible up to dxTry.
Definition: ODESolver.C:116
ODESolver(const ODESystem &ode, const dictionary &dict)
Construct for given ODESystem.
Definition: ODESolver.C:60
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
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:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const ODESystem & odes_
Reference to ODESystem.
Definition: ODESolver.H:58
label maxSteps_
The maximum number of sub-steps allowed for the integration step.
Definition: ODESolver.H:73
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const label maxN_
Maximum size of the ODESystem.
Definition: ODESolver.H:61
static const char nl
Definition: Ostream.H:265
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
defineTypeNameAndDebug(combustionModel, 0)
scalarField relTol_
Relative convergence tolerance per step.
Definition: ODESolver.H:70
scalar normalizeError(const scalarField &y0, const scalarField &y, const scalarField &err) const
Return the nomalized scalar error.
Definition: ODESolver.C:40
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
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
virtual label nEqns() const =0
Return the number of equations in the system.
label n_
Size of the ODESystem (adjustable)
Definition: ODESolver.H:64
Namespace for OpenFOAM.