src/ODE/ODESolvers/KRR4/KRR4.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     Fourth-order Kaps-Rentrop scheme with adjustive time-step size
  26     Numerical Recipes in C, Secn 16.6 page 739-742
  27     Alternative reference in Numerical Recipes in C++
  28
  29 \*---------------------------------------------------------------------------*/
  30
  31 #include "KRR4.H"
  32 #include "simpleMatrix.H"
  33 #include "addToRunTimeSelectionTable.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(Foam::KRR4, 0);
  38
  39 namespace Foam
  40 {
  41     addToRunTimeSelectionTable(ODESolver, KRR4, ODE);
  42
  43 const scalar
  44     KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25,
  45     KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296;
  46
  47 const scalar
  48     KRR4::gamma = 1.0/2.0,
  49     KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
  50     KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0,
  51     KRR4::c41 = -112.0/125.0, KRR4::c42 = -54.0/125.0, KRR4::c43 = -2.0/5.0,
  52     KRR4::b1 = 19.0/9.0, KRR4::b2 = 1.0/2.0, KRR4::b3 = 25.0/108.0,
  53     KRR4::b4 = 125.0/108.0,
  54     KRR4::e1 = 17.0/54.0, KRR4::e2 = 7.0/36.0, KRR4::e3 = 0.0,
  55     KRR4::e4 = 125.0/108.0,
  56     KRR4::c1X = 1.0/2.0, KRR4::c2X = -3.0/2.0, KRR4::c3X = 121.0/50.0,
  57     KRR4::c4X = 29.0/250.0,
  58     KRR4::a2X = 1.0, KRR4::a3X = 3.0/5.0;
  59 };
  60
  61
  62 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  63
  64 Foam::KRR4::KRR4(ODE& ode)
  65 :
  66     ODESolver(ode),
  67     yTemp_(ode_.nEqns()),
  68     dydxTemp_(ode_.nEqns()),
  69     g1_(ode_.nEqns()),
  70     g2_(ode_.nEqns()),
  71     g3_(ode_.nEqns()),
  72     g4_(ode_.nEqns()),
  73     yErr_(ode_.nEqns()),
  74     dfdx_(ode_.nEqns()),
  75     dfdy_(ode_.nEqns(), ode_.nEqns()),
  76     a_(ode_.nEqns(), ode_.nEqns()),
  77     pivotIndices_(ode_.nEqns())
  78 {}
  79
  80
  81 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  82
  83 void Foam::KRR4::solve
  84 (
  85     scalar& x,
  86     scalarField& y,
  87     scalarField& dydx,
  88     const scalar eps,
  89     const scalarField& yScale,
  90     const scalar hTry,
  91     scalar& hDid,
  92     scalar& hNext
  93 ) const
  94 {
  95     scalar xTemp = x;
  96     yTemp_ = y;
  97     dydxTemp_ = dydx;
  98     const label nEqns = ode_.nEqns();
  99
 100     ode_.jacobian(xTemp, yTemp_, dfdx_, dfdy_);
 101
 102     scalar h = hTry;
 103
 104     for (register label jtry=0; jtry<maxtry; jtry++)
 105     {
 106         for (register label i=0; i<nEqns; i++)
 107         {
 108             for (register label j=0; j<nEqns; j++)
 109             {
 110                 a_[i][j] = -dfdy_[i][j];
 111             }
 112
 113             a_[i][i] += 1.0/(gamma*h);
 114         }
 115
 116         simpleMatrix<scalar>::LUDecompose(a_, pivotIndices_);
 117
 118         for (register label i=0; i<nEqns; i++)
 119         {
 120             g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];
 121         }
 122
 123         simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g1_);
 124
 125         for (register label i=0; i<nEqns; i++)
 126         {
 127             y[i] = yTemp_[i] + a21*g1_[i];
 128         }
 129
 130         x = xTemp + a2X*h;
 131         ode_.derivatives(x, y, dydx_);
 132
 133         for (register label i=0; i<nEqns; i++)
 134         {
 135             g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/h;
 136         }
 137
 138         simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g2_);
 139
 140         for (register label i=0; i<nEqns; i++)
 141         {
 142             y[i] = yTemp_[i] + a31*g1_[i] + a32*g2_[i];
 143         }
 144
 145         x = xTemp + a3X*h;
 146         ode_.derivatives(x, y, dydx_);
 147
 148         for (register label i=0; i<nEqns; i++)
 149         {
 150             g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;
 151         }
 152
 153         simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g3_);
 154
 155         for (register label i=0; i<nEqns; i++)
 156         {
 157             g4_[i] = dydx_[i] + h*c4X*dfdx_[i]
 158                 + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;
 159         }
 160
 161         simpleMatrix<scalar>::LUBacksubstitute(a_, pivotIndices_, g4_);
 162
 163         for (register label i=0; i<nEqns; i++)
 164         {
 165             y[i] = yTemp_[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];
 166             yErr_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];
 167         }
 168
 169         x = xTemp + h;
 170
 171         if (x == xTemp)
 172         {
 173             FatalErrorIn("ODES::KRR4")
 174                 << "stepsize not significant"
 175                 << exit(FatalError);
 176         }
 177
 178         scalar maxErr = 0.0;
 179         for (register label i=0; i<nEqns; i++)
 180         {
 181             maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
 182         }
 183         maxErr /= eps;
 184
 185         if (maxErr <= 1.0)
 186         {
 187             hDid = h;
 188             hNext = (maxErr > errcon ? safety*h*pow(maxErr, pgrow) : grow*h);
 189             return;
 190         }
 191         else
 192         {
 193             hNext = safety*h*pow(maxErr, pshrink);
 194             h = (h >= 0.0 ? max(hNext, shrink*h) : min(hNext, shrink*h));
 195         }
 196     }
 197
 198     FatalErrorIn("ODES::KRR4")
 199         << "exceeded maxtry"
 200         << exit(FatalError);
 201 }
 202
 203
 204 // ************************************************************************* //