src/ODE/ODESolvers/RK/RK.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     Fifth-order Cash-Karp embedded Runge-Kutta scheme with error control and
  26     adjustive time-step size
  27     Numerical Recipes in C, Secn 16.2 page 714-722
  28     Alternative reference in Numerical Recipes in C++
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "RK.H"
  33 #include "addToRunTimeSelectionTable.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(Foam::RK, 0);
  38
  39 namespace Foam
  40 {
  41     addToRunTimeSelectionTable(ODESolver, RK, ODE);
  42
  43 const scalar
  44     RK::safety=0.9, RK::pGrow=-0.2, RK::pShrink=-0.25, RK::errCon=1.89e-4;
  45
  46 const scalar
  47     RK::a2 = 0.2, RK::a3 = 0.3, RK::a4 = 0.6, RK::a5 = 1.0, RK::a6 = 0.875,
  48     RK::b21 = 0.2, RK::b31 = 3.0/40.0, RK::b32 = 9.0/40.0,
  49     RK::b41 = 0.3, RK::b42 = -0.9, RK::b43 = 1.2,
  50     RK::b51 = -11.0/54.0, RK::b52 = 2.5, RK::b53 = -70.0/27.0,
  51     RK::b54 = 35.0/27.0,
  52     RK::b61 = 1631.0/55296.0, RK::b62 = 175.0/512.0, RK::b63 = 575.0/13824.0,
  53     RK::b64 = 44275.0/110592.0, RK::b65 = 253.0/4096.0,
  54     RK::c1 = 37.0/378.0, RK::c3 = 250.0/621.0,
  55     RK::c4 = 125.0/594.0, RK::c6 = 512.0/1771.0,
  56     RK::dc1 = RK::c1 - 2825.0/27648.0, RK::dc3 = RK::c3 - 18575.0/48384.0,
  57     RK::dc4 = RK::c4 - 13525.0/55296.0, RK::dc5 = -277.00/14336.0,
  58     RK::dc6 = RK::c6 - 0.25;
  59 };
  60
  61
  62 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  63
  64 Foam::RK::RK(ODE& ode)
  65 :
  66     ODESolver(ode),
  67     yTemp_(ode.nEqns()),
  68     ak2_(ode.nEqns()),
  69     ak3_(ode.nEqns()),
  70     ak4_(ode.nEqns()),
  71     ak5_(ode.nEqns()),
  72     ak6_(ode.nEqns()),
  73     yErr_(ode.nEqns()),
  74     yTemp2_(ode.nEqns())
  75 {}
  76
  77
  78 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  79
  80 void Foam::RK::solve
  81 (
  82     const scalar x,
  83     const scalarField& y,
  84     const scalarField& dydx,
  85     const scalar h,
  86     scalarField& yout,
  87     scalarField& yerr
  88 ) const
  89 {
  90     forAll(yTemp_, i)
  91     {
  92         yTemp_[i] = y[i] + b21*h*dydx[i];
  93     }
  94
  95     ode_.derivatives(x + a2*h, yTemp_, ak2_);
  96
  97     forAll(yTemp_, i)
  98     {
  99         yTemp_[i] = y[i] + h*(b31*dydx[i] + b32*ak2_[i]);
 100     }
 101
 102     ode_.derivatives(x + a3*h, yTemp_, ak3_);
 103
 104     forAll(yTemp_, i)
 105     {
 106         yTemp_[i] = y[i] + h*(b41*dydx[i] + b42*ak2_[i] + b43*ak3_[i]);
 107     }
 108
 109     ode_.derivatives(x + a4*h, yTemp_, ak4_);
 110
 111     forAll(yTemp_, i)
 112     {
 113         yTemp_[i] = y[i]
 114           + h*(b51*dydx[i] + b52*ak2_[i] + b53*ak3_[i] + b54*ak4_[i]);
 115     }
 116
 117     ode_.derivatives(x + a5*h, yTemp_, ak5_);
 118
 119     forAll(yTemp_, i)
 120     {
 121         yTemp_[i] = y[i]
 122           + h*
 123             (
 124                 b61*dydx[i] + b62*ak2_[i] + b63*ak3_[i]
 125               + b64*ak4_[i] + b65*ak5_[i]
 126             );
 127     }
 128
 129     ode_.derivatives(x + a6*h, yTemp_, ak6_);
 130
 131     forAll(yout, i)
 132     {
 133         yout[i] = y[i]
 134           + h*(c1*dydx[i] + c3*ak3_[i] + c4*ak4_[i] + c6*ak6_[i]);
 135     }
 136
 137     forAll(yerr, i)
 138     {
 139         yerr[i] =
 140             h*
 141             (
 142                 dc1*dydx[i] + dc3*ak3_[i] + dc4*ak4_[i]
 143               + dc5*ak5_[i] + dc6*ak6_[i]
 144             );
 145     }
 146 }
 147
 148
 149 void Foam::RK::solve
 150 (
 151     scalar& x,
 152     scalarField& y,
 153     scalarField& dydx,
 154     const scalar eps,
 155     const scalarField& yScale,
 156     const scalar hTry,
 157     scalar& hDid,
 158     scalar& hNext
 159 ) const
 160 {
 161     scalar h = hTry;
 162     scalar maxErr = 0.0;
 163
 164     for (;;)
 165     {
 166         solve(x, y, dydx, h, yTemp2_, yErr_);
 167
 168         maxErr = 0.0;
 169         for (register label i=0; i<ode_.nEqns(); i++)
 170         {
 171             maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
 172         }
 173         maxErr /= eps;
 174
 175         if (maxErr <= 1.0)
 176         {
 177             break;
 178         }
 179
 180         {
 181             scalar hTemp = safety*h*pow(maxErr, pShrink);
 182             h = (h >= 0.0 ? max(hTemp, 0.1*h) : min(hTemp, 0.1*h));
 183         }
 184
 185         if (h < VSMALL)
 186         {
 187             FatalErrorIn("RK::solve")
 188                 << "stepsize underflow"
 189                 << exit(FatalError);
 190         }
 191     }
 192
 193     hDid = h;
 194
 195     x += h;
 196     y = yTemp2_;
 197
 198     if (maxErr > errCon)
 199     {
 200         hNext = safety*h*pow(maxErr, pGrow);
 201     }
 202     else
 203     {
 204         hNext = 5.0*h;
 205     }
 206 }
 207
 208
 209 // ************************************************************************* //