ODESolver.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) 2011-2015 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 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
38 
40 :
41  odes_(ode),
42  n_(ode.nEqns()),
43  absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
44  relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-4)),
45  maxSteps_(10000)
46 {}
47 
48 
50 (
51  const ODESystem& ode,
52  const scalarField& absTol,
53  const scalarField& relTol
54 )
55 :
56  odes_(ode),
57  n_(ode.nEqns()),
58  absTol_(absTol),
59  relTol_(relTol),
60  maxSteps_(10000)
61 {}
62 
63 
64 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
65 
67 (
68  const scalarField& y0,
69  const scalarField& y,
70  const scalarField& err
71 ) const
72 {
73  // Calculate the maximum error
74  scalar maxErr = 0.0;
75  forAll(err, i)
76  {
77  scalar tol = absTol_[i] + relTol_[i]*max(mag(y0[i]), mag(y[i]));
78  maxErr = max(maxErr, mag(err[i])/tol);
79  }
80 
81  return maxErr;
82 }
83 
84 
86 (
87  scalar& x,
88  scalarField& y,
89  stepState& step
90 ) const
91 {
92  scalar x0 = x;
93  solve(x, y, step.dxTry);
94  step.dxDid = x - x0;
95 }
96 
97 
99 (
100  const scalar xStart,
101  const scalar xEnd,
102  scalarField& y,
103  scalar& dxTry
104 ) const
105 {
106  stepState step(dxTry);
107  scalar x = xStart;
108 
109  for (label nStep=0; nStep<maxSteps_; nStep++)
110  {
111  // Store previous iteration dxTry
112  scalar dxTry0 = step.dxTry;
113 
114  step.reject = false;
115 
116  // Check if this is a truncated step and set dxTry to integrate to xEnd
117  if ((x + step.dxTry - xEnd)*(x + step.dxTry - xStart) > 0)
118  {
119  step.last = true;
120  step.dxTry = xEnd - x;
121  }
122 
123  // Integrate as far as possible up to step.dxTry
124  solve(x, y, step);
125 
126  // Check if reached xEnd
127  if ((x - xEnd)*(xEnd - xStart) >= 0)
128  {
129  if (nStep > 0 && step.last)
130  {
131  step.dxTry = dxTry0;
132  }
133 
134  dxTry = step.dxTry;
135 
136  return;
137  }
138 
139  step.first = false;
140 
141  // If the step.dxTry was reject set step.prevReject
142  if (step.reject)
143  {
144  step.prevReject = true;
145  }
146  }
147 
149  << "Integration steps greater than maximum " << maxSteps_
150  << "xStart = " << xStart << ", xEnd = " << xEnd
151  << ", x = " << x << ", dxDid = " << step.dxDid
152  << exit(FatalError);
153 }
154 
155 
156 // ************************************************************************* //
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const double e
Elementary charge.
Definition: doubleFloat.H:78
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
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:70
ODESolver(const ODESolver &)
Disallow default bitwise copy construct.
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
defineTypeNameAndDebug(combustionModel, 0)
virtual void solve(scalar &x, scalarField &y, scalar &dxTry) const
Solve the ODE system as far as possible upto dxTry.
Definition: ODESolver.H:178
scalarField relTol_
Relative convergence tolerance per step.
Definition: ODESolver.H:67
scalar normalizeError(const scalarField &y0, const scalarField &y, const scalarField &err) const
Return the nomalized scalar error.
Definition: ODESolver.C:67
scalarField absTol_
Absolute convergence tolerance per step.
Definition: ODESolver.H:64
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual label nEqns() const =0
Return the number of equations in the system.
label n_
Size of the ODESystem.
Definition: ODESolver.H:61
Namespace for OpenFOAM.