applications/test/ODETest/ODETest.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 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
  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 Description
  25
  26 \*---------------------------------------------------------------------------*/
  27
  28 #include "argList.H"
  29 #include "IOmanip.H"
  30 #include "ODE.H"
  31 #include "ODESolver.H"
  32 #include "RK.H"
  33
  34 using namespace Foam;
  35
  36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  37
  38 class testODE
  39 :
  40     public ODE
  41 {
  42
  43 public:
  44
  45     testODE()
  46     {}
  47
  48     label nEqns() const
  49     {
  50         return 4;
  51     }
  52
  53     void derivatives
  54     (
  55         const scalar x,
  56         const scalarField& y,
  57         scalarField& dydx
  58     ) const
  59     {
  60         dydx[0] = -y[1];
  61         dydx[1] = y[0] - (1.0/x)*y[1];
  62         dydx[2] = y[1] - (2.0/x)*y[2];
  63         dydx[3] = y[2] - (3.0/x)*y[3];
  64     }
  65
  66     void jacobian
  67     (
  68         const scalar x,
  69         const scalarField& y,
  70         scalarField& dfdx,
  71         scalarSquareMatrix& dfdy
  72     ) const
  73     {
  74         dfdx[0] = 0.0;
  75         dfdx[1] = (1.0/sqr(x))*y[1];
  76         dfdx[2] = (2.0/sqr(x))*y[2];
  77         dfdx[3] = (3.0/sqr(x))*y[3];
  78
  79         dfdy[0][0] = 0.0;
  80         dfdy[0][1] = -1.0;
  81         dfdy[0][2] = 0.0;
  82         dfdy[0][3] = 0.0;
  83
  84         dfdy[1][0] = 1.0;
  85         dfdy[1][1] = -1.0/x;
  86         dfdy[1][2] = 0.0;
  87         dfdy[1][3] = 0.0;
  88
  89         dfdy[2][0] = 0.0;
  90         dfdy[2][1] = 1.0;
  91         dfdy[2][2] = -2.0/x;
  92         dfdy[2][3] = 0.0;
  93
  94         dfdy[3][0] = 0.0;
  95         dfdy[3][1] = 0.0;
  96         dfdy[3][2] = 1.0;
  97         dfdy[3][3] = -3.0/x;
  98     }
  99 };
 100
 101
 102 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 103 // Main program:
 104
 105 int main(int argc, char *argv[])
 106 {
 107     argList::validArgs.clear();
 108     argList::validArgs.append("ODESolver");
 109     argList args(argc, argv);
 110
 111     word ODESolverName(args.additionalArgs()[0]);
 112
 113     testODE ode;
 114     autoPtr<ODESolver> odeSolver = ODESolver::New(ODESolverName, ode);
 115
 116     scalar xStart = 1.0;
 117     scalarField yStart(ode.nEqns());
 118     yStart[0] = ::Foam::j0(xStart);
 119     yStart[1] = ::Foam::j1(xStart);
 120     yStart[2] = ::Foam::jn(2, xStart);
 121     yStart[3] = ::Foam::jn(3, xStart);
 122
 123     scalarField dyStart(ode.nEqns());
 124     ode.derivatives(xStart, yStart, dyStart);
 125
 126     Info<< setw(10) << "eps" << setw(12) << "hEst";
 127     Info<< setw(13) << "hDid" << setw(14) << "hNext" << endl;
 128     Info<< setprecision(6);
 129
 130     for (label i=0; i<15; i++)
 131     {
 132         scalar eps = ::Foam::exp(-scalar(i + 1));
 133
 134         scalar x = xStart;
 135         scalarField y = yStart;
 136         scalarField dydx = dyStart;
 137
 138         scalarField yScale(ode.nEqns(), 1.0);
 139         scalar hEst = 0.6;
 140         scalar hDid, hNext;
 141
 142         odeSolver->solve(ode, x, y, dydx, eps, yScale, hEst, hDid, hNext);
 143
 144         Info<< scientific << setw(13) << eps;
 145         Info<< fixed << setw(11) << hEst;
 146         Info<< setw(13) << hDid << setw(13) << hNext
 147             << setw(13) << y[0] << setw(13) << y[1]
 148             << setw(13) << y[2] << setw(13) << y[3]
 149             << endl;
 150     }
 151
 152     scalar x = xStart;
 153     scalar xEnd = x + 1.0;
 154     scalarField y = yStart;
 155
 156     scalarField yEnd(ode.nEqns());
 157     yEnd[0] = ::Foam::j0(xEnd);
 158     yEnd[1] = ::Foam::j1(xEnd);
 159     yEnd[2] = ::Foam::jn(2, xEnd);
 160     yEnd[3] = ::Foam::jn(3, xEnd);
 161
 162     scalar hEst = 0.5;
 163
 164     odeSolver->solve(ode, x, xEnd, y, 1e-4, hEst);
 165
 166     Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
 167     Info      << "Numerical:  y(2.0) = " << y << ", hEst = " << hEst << endl;
 168
 169     Info << "\nEnd\n" << endl;
 170
 171     return 0;
 172 }
 173
 174
 175 // ************************************************************************* //