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-2021 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  const label li,
120  scalar& dxTry
121 ) const
122 {
123  stepState step(dxTry);
124  solve(x, y, li, step);
125  dxTry = step.dxTry;
126 }
127 
128 
130 (
131  scalar& x,
132  scalarField& y,
133  const label li,
134  stepState& step
135 ) const
136 {
137  scalar x0 = x;
138  solve(x, y, li, step.dxTry);
139  step.dxDid = x - x0;
140 }
141 
142 
144 (
145  const scalar xStart,
146  const scalar xEnd,
147  scalarField& y,
148  const label li,
149  scalar& dxTry
150 ) const
151 {
152  stepState step(dxTry);
153  scalar x = xStart;
154 
155  for (label nStep=0; nStep<maxSteps_; nStep++)
156  {
157  // Store previous iteration dxTry
158  scalar dxTry0 = step.dxTry;
159 
160  step.reject = false;
161 
162  // Check if this is a truncated step and set dxTry to integrate to xEnd
163  if ((x + step.dxTry - xEnd)*(x + step.dxTry - xStart) > 0)
164  {
165  step.last = true;
166  step.dxTry = xEnd - x;
167  }
168 
169  // Integrate as far as possible up to step.dxTry
170  solve(x, y, li, step);
171 
172  // Check if reached xEnd
173  if ((x - xEnd)*(xEnd - xStart) >= 0)
174  {
175  if (nStep > 0 && step.last)
176  {
177  step.dxTry = dxTry0;
178  }
179 
180  dxTry = step.dxTry;
181 
182  return;
183  }
184 
185  step.first = false;
186 
187  // If the step.dxTry was reject set step.prevReject
188  if (step.reject)
189  {
190  step.prevReject = true;
191  }
192  }
193 
195  << "Integration steps greater than maximum " << maxSteps_ << nl
196  << " xStart = " << xStart << ", xEnd = " << xEnd
197  << ", x = " << x << ", dxDid = " << step.dxDid << nl
198  << " y = " << y
199  << exit(FatalError);
200 }
201 
202 
203 // ************************************************************************* //
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:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
scalar normaliseError(const scalarField &y0, const scalarField &y, const scalarField &err) const
Return the normalised scalar error.
Definition: ODESolver.C:40
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:260
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
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: ODESolver.C:116
scalarField relTol_
Relative convergence tolerance per step.
Definition: ODESolver.H:70
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:105
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.