applications/test/ODE/Test-ODE.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-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.append("ODESolver");
 108     argList args(argc, argv);
 109
 110     testODE ode;
 111     autoPtr<ODESolver> odeSolver = ODESolver::New(args[1], ode);
 112
 113     scalar xStart = 1.0;
 114     scalarField yStart(ode.nEqns());
 115     yStart[0] = ::Foam::j0(xStart);
 116     yStart[1] = ::Foam::j1(xStart);
 117     yStart[2] = ::Foam::jn(2, xStart);
 118     yStart[3] = ::Foam::jn(3, xStart);
 119
 120     scalarField dyStart(ode.nEqns());
 121     ode.derivatives(xStart, yStart, dyStart);
 122
 123     Info<< setw(10) << "eps" << setw(12) << "hEst";
 124     Info<< setw(13) << "hDid" << setw(14) << "hNext" << endl;
 125     Info<< setprecision(6);
 126
 127     for (label i=0; i<15; i++)
 128     {
 129         scalar eps = ::Foam::exp(-scalar(i + 1));
 130
 131         scalar x = xStart;
 132         scalarField y(yStart);
 133         scalarField dydx(dyStart);
 134
 135         scalarField yScale(ode.nEqns(), 1.0);
 136         scalar hEst = 0.6;
 137         scalar hDid, hNext;
 138
 139         odeSolver->solve(ode, x, y, dydx, eps, yScale, hEst, hDid, hNext);
 140
 141         Info<< scientific << setw(13) << eps;
 142         Info<< fixed << setw(11) << hEst;
 143         Info<< setw(13) << hDid << setw(13) << hNext
 144             << setw(13) << y[0] << setw(13) << y[1]
 145             << setw(13) << y[2] << setw(13) << y[3]
 146             << endl;
 147     }
 148
 149     scalar x = xStart;
 150     scalar xEnd = x + 1.0;
 151     scalarField y(yStart);
 152
 153     scalarField yEnd(ode.nEqns());
 154     yEnd[0] = ::Foam::j0(xEnd);
 155     yEnd[1] = ::Foam::j1(xEnd);
 156     yEnd[2] = ::Foam::jn(2, xEnd);
 157     yEnd[3] = ::Foam::jn(3, xEnd);
 158
 159     scalar hEst = 0.5;
 160
 161     odeSolver->solve(ode, x, xEnd, y, 1e-4, hEst);
 162
 163     Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
 164     Info      << "Numerical:  y(2.0) = " << y << ", hEst = " << hEst << endl;
 165
 166     Info<< "\nEnd\n" << endl;
 167
 168     return 0;
 169 }
 170
 171
 172 // ************************************************************************* //