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