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