src/ODE/ODESolvers/SIBS/SIBS.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Description
  26     Semi-Implicit Burlisch-Stoer scheme for stiff equation sets
  27     Numerical Recipes in C, Secn 16.6 page 742-747
  28     Alternative reference in Numerical Recipes in C++
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "SIBS.H"
  33 #include "addToRunTimeSelectionTable.H"
  34
  35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  36
  37 defineTypeNameAndDebug(Foam::SIBS, 0);
  38
  39 namespace Foam
  40 {
  41     addToRunTimeSelectionTable(ODESolver, SIBS, ODE);
  42
  43     const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70};
  44
  45     const scalar
  46         SIBS::safe1 = 0.25, SIBS::safe2 = 0.7,
  47         SIBS::redMax = 1.0e-5, SIBS::redMin = 0.0, SIBS::scaleMX = 0.1;
  48 };
  49
  50
  51 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  52
  53 Foam::SIBS::SIBS(ODE& ode)
  54 :
  55     ODESolver(ode),
  56     a_(iMaxX_),
  57     alpha_(kMaxX_, kMaxX_),
  58     d_p_(ode_.nEqns(), kMaxX_),
  59     x_p_(kMaxX_),
  60     err_(kMaxX_),
  61
  62     yTemp_(ode_.nEqns()),
  63     ySeq_(ode_.nEqns()),
  64     yErr_(ode_.nEqns()),
  65     dfdx_(ode_.nEqns()),
  66     dfdy_(ode_.nEqns(), ode_.nEqns()),
  67     first_(1),
  68     epsOld_(-1.0)
  69 {}
  70
  71
  72 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  73
  74 void Foam::SIBS::solve
  75 (
  76     scalar& x,
  77     scalarField& y,
  78     scalarField& dydx,
  79     const scalar eps,
  80     const scalarField& yScale,
  81     const scalar hTry,
  82     scalar& hDid,
  83     scalar& hNext
  84 ) const
  85 {
  86     bool exitflag = false;
  87     const label nEqns = ode_.nEqns();
  88
  89     if (eps != epsOld_)
  90     {
  91         hNext = xNew_ = -GREAT;
  92         scalar eps1 = safe1*eps;
  93         a_[0] = nSeq_[0] + 1;
  94
  95         for (register label k=0; k<kMaxX_; k++)
  96         {
  97             a_[k + 1] = a_[k] + nSeq_[k + 1];
  98         }
  99
 100         for (register label iq = 1; iq<kMaxX_; iq++)
 101         {
 102             for (register label k=0; k<iq; k++)
 103             {
 104                 alpha_[k][iq] =
 105                     pow(eps1, (a_[k + 1] - a_[iq + 1])
 106                    /((a_[iq + 1] - a_[0] + 1.0)*(2*k + 3)));
 107             }
 108         }
 109
 110         epsOld_ = eps;
 111         a_[0] += nEqns;
 112
 113         for (register label k=0; k<kMaxX_; k++)
 114         {
 115             a_[k + 1] = a_[k] + nSeq_[k + 1];
 116         }
 117
 118         for (kOpt_ = 1; kOpt_<kMaxX_ - 1; kOpt_++)
 119         {
 120             if (a_[kOpt_ + 1] > a_[kOpt_]*alpha_[kOpt_ - 1][kOpt_])
 121             {
 122                 break;
 123             }
 124         }
 125
 126         kMax_ = kOpt_;
 127     }
 128
 129     label k=0;
 130     scalar h = hTry;
 131     yTemp_ = y;
 132
 133     ode_.jacobian(x, y, dfdx_, dfdy_);
 134
 135     if (x != xNew_ || h != hNext)
 136     {
 137         first_ = 1;
 138         kOpt_ = kMax_;
 139     }
 140
 141     label km=0;
 142     label reduct=0;
 143     scalar maxErr = SMALL;
 144
 145     for (;;)
 146     {
 147         scalar red = 1.0;
 148
 149         for (k=0; k <= kMax_; k++)
 150         {
 151             xNew_ = x + h;
 152
 153             if (xNew_ == x)
 154             {
 155                 FatalErrorIn("ODES::SIBS")
 156                     << "step size underflow"
 157                     << exit(FatalError);
 158             }
 159
 160             SIMPR(x, yTemp_, dydx, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
 161             scalar xest = sqr(h/nSeq_[k]);
 162
 163             polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);
 164
 165             if (k != 0)
 166             {
 167                 maxErr = SMALL;
 168                 for (register label i=0; i<nEqns; i++)
 169                 {
 170                     maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
 171                 }
 172                 maxErr /= eps;
 173                 km = k - 1;
 174                 err_[km] = pow(maxErr/safe1, 1.0/(2*km + 3));
 175             }
 176
 177             if (k != 0 && (k >= kOpt_ - 1 || first_))
 178             {
 179                 if (maxErr < 1.0)
 180                 {
 181                     exitflag = true;
 182                     break;
 183                 }
 184
 185                 if (k == kMax_ || k == kOpt_ + 1)
 186                 {
 187                     red = safe2/err_[km];
 188                     break;
 189                 }
 190                 else if (k == kOpt_ && alpha_[kOpt_ - 1][kOpt_] < err_[km])
 191                 {
 192                     red = 1.0/err_[km];
 193                     break;
 194                 }
 195                 else if (kOpt_ == kMax_ && alpha_[km][kMax_ - 1] < err_[km])
 196                 {
 197                     red = alpha_[km][kMax_ - 1]*safe2/err_[km];
 198                     break;
 199                 }
 200                 else if (alpha_[km][kOpt_] < err_[km])
 201                 {
 202                     red = alpha_[km][kOpt_ - 1]/err_[km];
 203                     break;
 204                 }
 205             }
 206         }
 207
 208         if (exitflag)
 209         {
 210             break;
 211         }
 212
 213         red = min(red, redMin);
 214         red = max(red, redMax);
 215         h *= red;
 216         reduct = 1;
 217     }
 218
 219     x = xNew_;
 220     hDid = h;
 221     first_=0;
 222     scalar wrkmin = GREAT;
 223     scalar scale = 1.0;
 224
 225     for (register label kk=0; kk<=km; kk++)
 226     {
 227         scalar fact = max(err_[kk], scaleMX);
 228         scalar work = fact*a_[kk + 1];
 229         if (work < wrkmin)
 230         {
 231             scale = fact;
 232             wrkmin = work;
 233             kOpt_ = kk + 1;
 234         }
 235     }
 236
 237     hNext = h/scale;
 238
 239     if (kOpt_ >= k && kOpt_ != kMax_ && !reduct)
 240     {
 241         scalar fact = max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX);
 242         if (a_[kOpt_ + 1]*fact <= wrkmin)
 243         {
 244             hNext = h/fact;
 245             kOpt_++;
 246         }
 247     }
 248 }
 249
 250
 251 // ************************************************************************* //