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