FreeFOAM The Cross-Platform CFD Toolkit
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) 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 <ODE/ODESolver.H>
27 
28 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
29 
31 namespace Foam
32 {
33  defineRunTimeSelectionTable(ODESolver, ODE);
34 };
35 
36 
37 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
38 
40 :
41  n_(ode.nEqns()),
42  yScale_(n_),
43  dydx_(n_)
44 {}
45 
46 
47 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
48 
50 (
51  const ODE& ode,
52  const scalar xStart,
53  const scalar xEnd,
54  scalarField& y,
55  const scalar eps,
56  scalar& hEst
57 ) const
58 {
59  const label MAXSTP = 10000;
60 
61  scalar x = xStart;
62  scalar h = hEst;
63  scalar hNext = 0;
64  scalar hPrev = 0;
65 
66  for (label nStep=0; nStep<MAXSTP; nStep++)
67  {
68  ode.derivatives(x, y, dydx_);
69 
70  for (label i=0; i<n_; i++)
71  {
72  yScale_[i] = mag(y[i]) + mag(dydx_[i]*h) + SMALL;
73  }
74 
75  if ((x + h - xEnd)*(x + h - xStart) > 0.0)
76  {
77  h = xEnd - x;
78  hPrev = hNext;
79  }
80 
81  hNext = 0;
82  scalar hDid;
83  solve(ode, x, y, dydx_, eps, yScale_, h, hDid, hNext);
84 
85  if ((x - xEnd)*(xEnd - xStart) >= 0.0)
86  {
87  if (hPrev != 0)
88  {
89  hEst = hPrev;
90  }
91  else
92  {
93  hEst = hNext;
94  }
95 
96  return;
97  }
98 
99  h = hNext;
100  }
101 
103  (
104  "ODESolver::solve"
105  "(const ODE& ode, const scalar xStart, const scalar xEnd,"
106  "scalarField& yStart, const scalar eps, scalar& hEst) const"
107  ) << "Too many integration steps"
108  << exit(FatalError);
109 }
110 
111 // ************************ vim: set sw=4 sts=4 et: ************************ //