src/ODE/ODESolvers/RK/RK.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "RK.H"
  28 #include "addToRunTimeSelectionTable.H"
  29
  30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  31
  32 defineTypeNameAndDebug(Foam::RK, 0);
  33
  34 namespace Foam
  35 {
  36     addToRunTimeSelectionTable(ODESolver, RK, ODE);
  37
  38 const scalar
  39     RK::safety=0.9, RK::pGrow=-0.2, RK::pShrink=-0.25, RK::errCon=1.89e-4;
  40
  41 const scalar
  42     RK::a2 = 0.2, RK::a3 = 0.3, RK::a4 = 0.6, RK::a5 = 1.0, RK::a6 = 0.875,
  43     RK::b21 = 0.2, RK::b31 = 3.0/40.0, RK::b32 = 9.0/40.0,
  44     RK::b41 = 0.3, RK::b42 = -0.9, RK::b43 = 1.2,
  45     RK::b51 = -11.0/54.0, RK::b52 = 2.5, RK::b53 = -70.0/27.0,
  46     RK::b54 = 35.0/27.0,
  47     RK::b61 = 1631.0/55296.0, RK::b62 = 175.0/512.0, RK::b63 = 575.0/13824.0,
  48     RK::b64 = 44275.0/110592.0, RK::b65 = 253.0/4096.0,
  49     RK::c1 = 37.0/378.0, RK::c3 = 250.0/621.0,
  50     RK::c4 = 125.0/594.0, RK::c6 = 512.0/1771.0,
  51     RK::dc1 = RK::c1 - 2825.0/27648.0, RK::dc3 = RK::c3 - 18575.0/48384.0,
  52     RK::dc4 = RK::c4 - 13525.0/55296.0, RK::dc5 = -277.00/14336.0,
  53     RK::dc6 = RK::c6 - 0.25;
  54 };
  55
  56
  57 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  58
  59 Foam::RK::RK(const ODE& ode)
  60 :
  61     ODESolver(ode),
  62     yTemp_(n_),
  63     ak2_(n_),
  64     ak3_(n_),
  65     ak4_(n_),
  66     ak5_(n_),
  67     ak6_(n_),
  68     yErr_(n_),
  69     yTemp2_(n_)
  70 {}
  71
  72
  73 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  74
  75 void Foam::RK::solve
  76 (
  77     const ODE& ode,
  78     const scalar x,
  79     const scalarField& y,
  80     const scalarField& dydx,
  81     const scalar h,
  82     scalarField& yout,
  83     scalarField& yerr
  84 ) const
  85 {
  86     forAll(yTemp_, i)
  87     {
  88         yTemp_[i] = y[i] + b21*h*dydx[i];
  89     }
  90
  91     ode.derivatives(x + a2*h, yTemp_, ak2_);
  92
  93     forAll(yTemp_, i)
  94     {
  95         yTemp_[i] = y[i] + h*(b31*dydx[i] + b32*ak2_[i]);
  96     }
  97
  98     ode.derivatives(x + a3*h, yTemp_, ak3_);
  99
 100     forAll(yTemp_, i)
 101     {
 102         yTemp_[i] = y[i] + h*(b41*dydx[i] + b42*ak2_[i] + b43*ak3_[i]);
 103     }
 104
 105     ode.derivatives(x + a4*h, yTemp_, ak4_);
 106
 107     forAll(yTemp_, i)
 108     {
 109         yTemp_[i] = y[i]
 110           + h*(b51*dydx[i] + b52*ak2_[i] + b53*ak3_[i] + b54*ak4_[i]);
 111     }
 112
 113     ode.derivatives(x + a5*h, yTemp_, ak5_);
 114
 115     forAll(yTemp_, i)
 116     {
 117         yTemp_[i] = y[i]
 118           + h*
 119             (
 120                 b61*dydx[i] + b62*ak2_[i] + b63*ak3_[i]
 121               + b64*ak4_[i] + b65*ak5_[i]
 122             );
 123     }
 124
 125     ode.derivatives(x + a6*h, yTemp_, ak6_);
 126
 127     forAll(yout, i)
 128     {
 129         yout[i] = y[i]
 130           + h*(c1*dydx[i] + c3*ak3_[i] + c4*ak4_[i] + c6*ak6_[i]);
 131     }
 132
 133     forAll(yerr, i)
 134     {
 135         yerr[i] =
 136             h*
 137             (
 138                 dc1*dydx[i] + dc3*ak3_[i] + dc4*ak4_[i]
 139               + dc5*ak5_[i] + dc6*ak6_[i]
 140             );
 141     }
 142 }
 143
 144
 145 void Foam::RK::solve
 146 (
 147     const ODE& ode,
 148     scalar& x,
 149     scalarField& y,
 150     scalarField& dydx,
 151     const scalar eps,
 152     const scalarField& yScale,
 153     const scalar hTry,
 154     scalar& hDid,
 155     scalar& hNext
 156 ) const
 157 {
 158     scalar h = hTry;
 159     scalar maxErr = 0.0;
 160
 161     for (;;)
 162     {
 163         solve(ode, x, y, dydx, h, yTemp2_, yErr_);
 164
 165         maxErr = 0.0;
 166         for (register label i=0; i<n_; i++)
 167         {
 168             maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
 169         }
 170         maxErr /= eps;
 171
 172         if (maxErr <= 1.0)
 173         {
 174             break;
 175         }
 176
 177         {
 178             scalar hTemp = safety*h*pow(maxErr, pShrink);
 179             h = (h >= 0.0 ? max(hTemp, 0.1*h) : min(hTemp, 0.1*h));
 180         }
 181
 182         if (h < VSMALL)
 183         {
 184             FatalErrorIn("RK::solve")
 185                 << "stepsize underflow"
 186                 << exit(FatalError);
 187         }
 188     }
 189
 190     hDid = h;
 191
 192     x += h;
 193     y = yTemp2_;
 194
 195     if (maxErr > errCon)
 196     {
 197         hNext = safety*h*pow(maxErr, pGrow);
 198     }
 199     else
 200     {
 201         hNext = 5.0*h;
 202     }
 203 }
 204
 205
 206 // ************************************************************************* //