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