diff --git a/src/ODE/Make/files b/src/ODE/Make/files index 0c869a7b84..1ef8f8ab81 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -5,8 +5,8 @@ ODESolvers/adaptiveSolver/adaptiveSolver.C ODESolvers/RKF45/RKF45.C ODESolvers/RKCK45/RKCK45.C ODESolvers/RKDP45/RKDP45.C -ODESolvers/KRR43/KRR43.C ODESolvers/Rosenbrock43/Rosenbrock43.C +ODESolvers/rodas43/rodas43.C ODESolvers/SIBS/SIBS.C ODESolvers/SIBS/SIMPR.C ODESolvers/SIBS/polyExtrapolate.C diff --git a/src/ODE/ODESolvers/KRR43/KRR43.C b/src/ODE/ODESolvers/KRR43/KRR43.C deleted file mode 100644 index 60b6795f93..0000000000 --- a/src/ODE/ODESolvers/KRR43/KRR43.C +++ /dev/null @@ -1,174 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see . - -\*---------------------------------------------------------------------------*/ - -#include "KRR43.H" -#include "addToRunTimeSelectionTable.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(KRR43, 0); - addToRunTimeSelectionTable(ODESolver, KRR43, dictionary); - -const scalar - KRR43::gamma = 1.0/2.0, - KRR43::a21 = 2.0, - KRR43::a31 = 48.0/25.0, - KRR43::a32 = 6.0/25.0, - KRR43::c21 = -8.0, - KRR43::c31 = 372.0/25.0, - KRR43::c32 = 12.0/5.0, - KRR43::c41 = -112.0/125.0, - KRR43::c42 = -54.0/125.0, - KRR43::c43 = -2.0/5.0, - KRR43::b1 = 19.0/9.0, - KRR43::b2 = 1.0/2.0, - KRR43::b3 = 25.0/108.0, - KRR43::b4 = 125.0/108.0, - KRR43::e1 = 17.0/54.0, - KRR43::e2 = 7.0/36.0, - KRR43::e3 = 0.0, - KRR43::e4 = 125.0/108.0, - KRR43::c1X = 1.0/2.0, - KRR43::c2X = -3.0/2.0, - KRR43::c3X = 121.0/50.0, - KRR43::c4X = 29.0/250.0, - KRR43::a2X = 1.0, - KRR43::a3X = 3.0/5.0; -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::KRR43::KRR43(const ODESystem& ode, const dictionary& dict) -: - ODESolver(ode, dict), - adaptiveSolver(ode, dict), - k1_(n_), - k2_(n_), - k3_(n_), - k4_(n_), - err_(n_), - dydx_(n_), - dfdx_(n_), - dfdy_(n_, n_), - a_(n_, n_), - pivotIndices_(n_) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -Foam::scalar Foam::KRR43::solve -( - const ODESystem& ode, - const scalar x0, - const scalarField& y0, - const scalarField& dydx0, - const scalar dx, - scalarField& y -) const -{ - ode.jacobian(x0, y0, dfdx_, dfdy_); - - for (register label i=0; i. + +\*---------------------------------------------------------------------------*/ + +#include "rodas43.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(rodas43, 0); + addToRunTimeSelectionTable(ODESolver, rodas43, dictionary); + +const scalar + rodas43::c2 = 0.386, + rodas43::c3 = 0.21, + rodas43::c4 = 0.63, + rodas43::bet2p = 0.0317, + rodas43::bet3p = 0.0635, + rodas43::bet4p = 0.3438, + rodas43::d1 = 0.25, + rodas43::d2 = -0.1043, + rodas43::d3 = 0.1035, + rodas43::d4 = -0.3620000000000023e-01, + rodas43::a21 = 0.1544e1, + rodas43::a31 = 0.9466785280815826, + rodas43::a32 = 0.2557011698983284, + rodas43::a41 = 0.3314825187068521e1, + rodas43::a42 = 0.2896124015972201e1, + rodas43::a43 = 0.9986419139977817, + rodas43::a51 = 0.1221224509226641e1, + rodas43::a52 = 0.6019134481288629e1, + rodas43::a53 = 0.1253708332932087e2, + rodas43::a54 = -0.6878860361058950, + rodas43::c21 = -0.56688e1, + rodas43::c31 = -0.2430093356833875e1, + rodas43::c32 = -0.2063599157091915, + rodas43::c41 = -0.1073529058151375, + rodas43::c42 = -0.9594562251023355e1, + rodas43::c43 = -0.2047028614809616e2, + rodas43::c51 = 0.7496443313967647e1, + rodas43::c52 = -0.1024680431464352e2, + rodas43::c53 = -0.3399990352819905e2, + rodas43::c54 = 0.1170890893206160e2, + rodas43::c61 = 0.8083246795921522e1, + rodas43::c62 = -0.7981132988064893e1, + rodas43::c63 = -0.3152159432874371e2, + rodas43::c64 = 0.1631930543123136e2, + rodas43::c65 = -0.6058818238834054e1, + rodas43::gamma = 0.25; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict) +: + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + k1_(n_), + k2_(n_), + k3_(n_), + k4_(n_), + k5_(n_), + dy_(n_), + err_(n_), + dydx_(n_), + dfdx_(n_), + dfdy_(n_, n_), + a_(n_, n_), + pivotIndices_(n_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::scalar Foam::rodas43::solve +( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y +) const +{ + ode.jacobian(x0, y0, dfdx_, dfdy_); + + for (register label i=0; i. Class - Foam::KRR43 + Foam::rodas43 Description - 4th Order Kaps Rentrop Rosenbrock stiff ODE solver + Stiffly-stable embedded Rosenbrock ODE solver of order (3)4. References: \verbatim - "Generalized Runge-Kutta methods of order four with stepsize control - for stiff ordinary differential equations" - Kaps, P., - Rentrop, P., - Numerische Thematic, vol. 33, 1979, pp. 55–68. - \endverbatim - - Constants from: - \verbatim - "Implementation of Rosenbrock Methods" - Shampine, L. F., - ACM Transactions on Mathematical Software, vol. 8, 1982, pp. 93–113. + "Solving Ordinary Differential Equations II: Stiff + and Differential-Algebraic Problems, second edition", + Hairer, E., + Nørsett, S., + Wanner, G., + Springer-Verlag, Berlin. 1996. \endverbatim SourceFiles - KRR43.C + rodas43.C \*---------------------------------------------------------------------------*/ -#ifndef KRR43_H -#define KRR43_H +#ifndef rodas43_H +#define rodas43_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -60,10 +54,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class KRR43 Declaration + Class rodas43 Declaration \*---------------------------------------------------------------------------*/ -class KRR43 +class rodas43 : public ODESolver, public adaptiveSolver @@ -74,6 +68,8 @@ class KRR43 mutable scalarField k2_; mutable scalarField k3_; mutable scalarField k4_; + mutable scalarField k5_; + mutable scalarField dy_; mutable scalarField err_; mutable scalarField dydx_; mutable scalarField dfdx_; @@ -82,25 +78,28 @@ class KRR43 mutable labelList pivotIndices_; static const scalar - gamma, + c2, c3, c4, + bet2p, bet3p, bet4p, + d1, d2, d3, d4, a21, a31, a32, - c21, c31, c32, c41, c42, c43, - b1, b2, b3, b4, - e1, e2, e3, e4, - c1X, c2X, c3X, c4X, - a2X, a3X; - + a41, a42, a43, + a51, a52, a53, a54, + c21, c31, c32, + c41, c42, c43, + c51, c52, c53, c54, + c61, c62, c63, c64, c65, + gamma; public: //- Runtime type information - TypeName("KRR43"); + TypeName("rodas43"); // Constructors //- Construct from ODE - KRR43(const ODESystem& ode, const dictionary& dict); + rodas43(const ODESystem& ode, const dictionary& dict); // Member Functions diff --git a/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties index 2c59e83d27..0fc57b4bba 100644 --- a/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties +++ b/tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties @@ -33,7 +33,7 @@ EulerImplicitCoeffs odeCoeffs { - solver KRR4; + solver Rosenbrock43; absTol 1e-12; relTol 0.01; } diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties index 41e24d715b..93d33d5232 100644 --- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties +++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties @@ -33,7 +33,7 @@ EulerImplicitCoeffs odeCoeffs { - solver KRR4; + solver Rosenbrock43; absTol 1e-12; relTol 0.01; }