src/ODE/ODESolvers/SIBS/SIMPR.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 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 "SIBS.H"
  28
  29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  30
  31 void Foam::SIBS::SIMPR
  32 (
  33     const ODE& ode,
  34     const scalar xStart,
  35     const scalarField& y,
  36     const scalarField& dydx,
  37     const scalarField& dfdx,
  38     const scalarSquareMatrix& dfdy,
  39     const scalar deltaX,
  40     const label nSteps,
  41     scalarField& yEnd
  42 ) const
  43 {
  44     scalar h = deltaX/nSteps;
  45
  46     scalarSquareMatrix a(n_);
  47     for (register label i=0; i<n_; i++)
  48     {
  49         for (register label j=0; j<n_; j++)
  50         {
  51             a[i][j] = -h*dfdy[i][j];
  52         }
  53         ++a[i][i];
  54     }
  55
  56     labelList pivotIndices(n_);
  57     LUDecompose(a, pivotIndices);
  58
  59     for (register label i=0; i<n_; i++)
  60     {
  61         yEnd[i] = h*(dydx[i] + h*dfdx[i]);
  62     }
  63
  64     LUBacksubstitute(a, pivotIndices, yEnd);
  65
  66     scalarField del(yEnd);
  67     scalarField ytemp(n_);
  68
  69     for (register label i=0; i<n_; i++)
  70     {
  71         ytemp[i] = y[i] + del[i];
  72     }
  73
  74     scalar x = xStart + h;
  75
  76     ode.derivatives(x, ytemp, yEnd);
  77
  78     for (register label nn=2; nn<=nSteps; nn++)
  79     {
  80         for (register label i=0; i<n_; i++)
  81         {
  82             yEnd[i] = h*yEnd[i] - del[i];
  83         }
  84
  85         LUBacksubstitute(a, pivotIndices, yEnd);
  86
  87         for (register label i=0; i<n_; i++)
  88         {
  89             ytemp[i] += (del[i] += 2.0*yEnd[i]);
  90         }
  91
  92         x += h;
  93
  94         ode.derivatives(x, ytemp, yEnd);
  95     }
  96     for (register label i=0; i<n_; i++)
  97     {
  98         yEnd[i] = h*yEnd[i] - del[i];
  99     }
 100
 101     LUBacksubstitute(a, pivotIndices, yEnd);
 102
 103     for (register label i=0; i<n_; i++)
 104     {
 105         yEnd[i] += ytemp[i];
 106     }
 107 }
 108
 109
 110 // ************************************************************************* //