FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

SIMPR.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "SIBS.H"
00027 
00028 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00029 
00030 void Foam::SIBS::SIMPR
00031 (
00032     const ODE& ode,
00033     const scalar xStart,
00034     const scalarField& y,
00035     const scalarField& dydx,
00036     const scalarField& dfdx,
00037     const scalarSquareMatrix& dfdy,
00038     const scalar deltaX,
00039     const label nSteps,
00040     scalarField& yEnd
00041 ) const
00042 {
00043     scalar h = deltaX/nSteps;
00044 
00045     scalarSquareMatrix a(n_);
00046     for (register label i=0; i<n_; i++)
00047     {
00048         for (register label j=0; j<n_; j++)
00049         {
00050             a[i][j] = -h*dfdy[i][j];
00051         }
00052         ++a[i][i];
00053     }
00054 
00055     labelList pivotIndices(n_);
00056     LUDecompose(a, pivotIndices);
00057 
00058     for (register label i=0; i<n_; i++)
00059     {
00060         yEnd[i] = h*(dydx[i] + h*dfdx[i]);
00061     }
00062 
00063     LUBacksubstitute(a, pivotIndices, yEnd);
00064 
00065     scalarField del(yEnd);
00066     scalarField ytemp(n_);
00067 
00068     for (register label i=0; i<n_; i++)
00069     {
00070         ytemp[i] = y[i] + del[i];
00071     }
00072 
00073     scalar x = xStart + h;
00074 
00075     ode.derivatives(x, ytemp, yEnd);
00076 
00077     for (register label nn=2; nn<=nSteps; nn++)
00078     {
00079         for (register label i=0; i<n_; i++)
00080         {
00081             yEnd[i] = h*yEnd[i] - del[i];
00082         }
00083 
00084         LUBacksubstitute(a, pivotIndices, yEnd);
00085 
00086         for (register label i=0; i<n_; i++)
00087         {
00088             ytemp[i] += (del[i] += 2.0*yEnd[i]);
00089         }
00090 
00091         x += h;
00092 
00093         ode.derivatives(x, ytemp, yEnd);
00094     }
00095     for (register label i=0; i<n_; i++)
00096     {
00097         yEnd[i] = h*yEnd[i] - del[i];
00098     }
00099 
00100     LUBacksubstitute(a, pivotIndices, yEnd);
00101 
00102     for (register label i=0; i<n_; i++)
00103     {
00104         yEnd[i] += ytemp[i];
00105     }
00106 }
00107 
00108 
00109 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines