ODESolvers: Completed Rosenbrock methods and removed legacy KRR4

This commit is contained in:
Henry
2013-11-04 12:21:40 +00:00
parent 29b3e9adfe
commit d9cdb08934
7 changed files with 271 additions and 223 deletions

View File

@ -5,8 +5,8 @@ ODESolvers/adaptiveSolver/adaptiveSolver.C
ODESolvers/RKF45/RKF45.C ODESolvers/RKF45/RKF45.C
ODESolvers/RKCK45/RKCK45.C ODESolvers/RKCK45/RKCK45.C
ODESolvers/RKDP45/RKDP45.C ODESolvers/RKDP45/RKDP45.C
ODESolvers/KRR43/KRR43.C
ODESolvers/Rosenbrock43/Rosenbrock43.C ODESolvers/Rosenbrock43/Rosenbrock43.C
ODESolvers/rodas43/rodas43.C
ODESolvers/SIBS/SIBS.C ODESolvers/SIBS/SIBS.C
ODESolvers/SIBS/SIMPR.C ODESolvers/SIBS/SIMPR.C
ODESolvers/SIBS/polyExtrapolate.C ODESolvers/SIBS/polyExtrapolate.C

View File

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<n_; i++)
{
for (register label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
}
a_[i][i] += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
forAll(k1_, i)
{
k1_[i] = dydx0[i] + dx*c1X*dfdx_[i];
}
LUBacksubstitute(a_, pivotIndices_, k1_);
forAll(y, i)
{
y[i] = y0[i] + a21*k1_[i];
}
ode.derivatives(x0 + a2X*dx, y, dydx_);
forAll(k2_, i)
{
k2_[i] = dydx_[i] + dx*c2X*dfdx_[i] + c21*k1_[i]/dx;
}
LUBacksubstitute(a_, pivotIndices_, k2_);
forAll(y, i)
{
y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
}
ode.derivatives(x0 + a3X*dx, y, dydx_);
forAll(k3_, i)
{
k3_[i] = dydx_[i] + dx*c3X*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, k3_);
forAll(k4_, i)
{
k4_[i] = dydx_[i] + dx*c4X*dfdx_[i]
+ (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, k4_);
forAll(y, i)
{
y[i] = y0[i] + b1*k1_[i] + b2*k2_[i] + b3*k3_[i] + b4*k4_[i];
err_[i] = e1*k1_[i] + e2*k2_[i] + e3*k3_[i] + e4*k4_[i];
}
return normalizeError(y0, y, err_);
}
void Foam::KRR43::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
}
// ************************************************************************* //

View File

@ -25,15 +25,16 @@ Class
Foam::Rosenbrock43 Foam::Rosenbrock43
Description Description
4th Order Kaps Rentrop Rosenbrock stiff ODE solver Embedded Rosenbrock ODE solver of order (3)4.
References: Based on code from:
\verbatim \verbatim
"Generalized Runge-Kutta methods of order four with stepsize control "Solving Ordinary Differential Equations II: Stiff
for stiff ordinary differential equations" and Differential-Algebraic Problems, second edition",
Kaps, P., Hairer, E.,
Rentrop, P., Nørsett, S.,
Numerische Thematic, vol. 33, 1979, pp. 5568. Wanner, G.,
Springer-Verlag, Berlin. 1996.
\endverbatim \endverbatim
Constants from: Constants from:
@ -43,16 +44,6 @@ Description
ACM Transactions on Mathematical Software, vol. 8, 1982, pp. 93113. ACM Transactions on Mathematical Software, vol. 8, 1982, pp. 93113.
\endverbatim \endverbatim
Based on code from:
\verbatim
"Solving Ordinary Differential Equations I: Nonstiff Problems,
second edition",
Hairer, E.,
Nørsett, S.,
Wanner, G.,
Springer-Verlag, Berlin. 1993, ISBN 3-540-56670-8.
\endverbatim
SourceFiles SourceFiles
Rosenbrock43.C Rosenbrock43.C

View File

@ -0,0 +1,232 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<n_; i++)
{
for (register label j=0; j<n_; j++)
{
a_[i][j] = -dfdy_[i][j];
}
a_[i][i] += 1.0/(gamma*dx);
}
LUDecompose(a_, pivotIndices_);
// Calculate k1:
forAll(k1_, i)
{
k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
}
LUBacksubstitute(a_, pivotIndices_, k1_);
// Calculate k2:
forAll(y, i)
{
y[i] = y0[i] + a21*k1_[i];
}
ode.derivatives(x0 + c2*dx, y, dydx_);
forAll(k2_, i)
{
k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
}
LUBacksubstitute(a_, pivotIndices_, k2_);
// Calculate k3:
forAll(y, i)
{
y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
}
ode.derivatives(x0 + c3*dx, y, dydx_);
forAll(k3_, i)
{
k3_[i] = dydx_[i] + dx*d3*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, k3_);
// Calculate k4:
forAll(y, i)
{
y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
}
ode.derivatives(x0 + c4*dx, y, dydx_);
forAll(k4_, i)
{
k4_[i] = dydx_[i] + dx*d4*dfdx_[i]
+ (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, k4_);
// Calculate k5:
forAll(y, i)
{
dy_[i] = a51*k1_[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i];
y[i] = y0[i] + dy_[i];
}
ode.derivatives(x0 + dx, y, dydx_);
forAll(k5_, i)
{
k5_[i] = dydx_[i]
+ (c51*k1_[i] + c52*k2_[i] + c53*k3_[i] + c54*k4_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, k5_);
// Calculate new state and error
forAll(y, i)
{
dy_[i] += k5_[i];
y[i] = y0[i] + dy_[i];
}
ode.derivatives(x0 + dx, y, dydx_);
forAll(err_, i)
{
err_[i] = dydx_[i]
+ (c61*k1_[i] + c62*k2_[i] + c63*k3_[i] + c64*k4_[i] + c65*k5_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, err_);
forAll(y, i)
{
y[i] = y0[i] + dy_[i] + err_[i];
}
return normalizeError(y0, y, err_);
}
void Foam::rodas43::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -22,34 +22,28 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class Class
Foam::KRR43 Foam::rodas43
Description Description
4th Order Kaps Rentrop Rosenbrock stiff ODE solver Stiffly-stable embedded Rosenbrock ODE solver of order (3)4.
References: References:
\verbatim \verbatim
"Generalized Runge-Kutta methods of order four with stepsize control "Solving Ordinary Differential Equations II: Stiff
for stiff ordinary differential equations" and Differential-Algebraic Problems, second edition",
Kaps, P., Hairer, E.,
Rentrop, P., Nørsett, S.,
Numerische Thematic, vol. 33, 1979, pp. 5568. Wanner, G.,
\endverbatim Springer-Verlag, Berlin. 1996.
Constants from:
\verbatim
"Implementation of Rosenbrock Methods"
Shampine, L. F.,
ACM Transactions on Mathematical Software, vol. 8, 1982, pp. 93113.
\endverbatim \endverbatim
SourceFiles SourceFiles
KRR43.C rodas43.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#ifndef KRR43_H #ifndef rodas43_H
#define KRR43_H #define rodas43_H
#include "ODESolver.H" #include "ODESolver.H"
#include "adaptiveSolver.H" #include "adaptiveSolver.H"
@ -60,10 +54,10 @@ namespace Foam
{ {
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class KRR43 Declaration Class rodas43 Declaration
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
class KRR43 class rodas43
: :
public ODESolver, public ODESolver,
public adaptiveSolver public adaptiveSolver
@ -74,6 +68,8 @@ class KRR43
mutable scalarField k2_; mutable scalarField k2_;
mutable scalarField k3_; mutable scalarField k3_;
mutable scalarField k4_; mutable scalarField k4_;
mutable scalarField k5_;
mutable scalarField dy_;
mutable scalarField err_; mutable scalarField err_;
mutable scalarField dydx_; mutable scalarField dydx_;
mutable scalarField dfdx_; mutable scalarField dfdx_;
@ -82,25 +78,28 @@ class KRR43
mutable labelList pivotIndices_; mutable labelList pivotIndices_;
static const scalar static const scalar
gamma, c2, c3, c4,
bet2p, bet3p, bet4p,
d1, d2, d3, d4,
a21, a31, a32, a21, a31, a32,
c21, c31, c32, c41, c42, c43, a41, a42, a43,
b1, b2, b3, b4, a51, a52, a53, a54,
e1, e2, e3, e4, c21, c31, c32,
c1X, c2X, c3X, c4X, c41, c42, c43,
a2X, a3X; c51, c52, c53, c54,
c61, c62, c63, c64, c65,
gamma;
public: public:
//- Runtime type information //- Runtime type information
TypeName("KRR43"); TypeName("rodas43");
// Constructors // Constructors
//- Construct from ODE //- Construct from ODE
KRR43(const ODESystem& ode, const dictionary& dict); rodas43(const ODESystem& ode, const dictionary& dict);
// Member Functions // Member Functions

View File

@ -33,7 +33,7 @@ EulerImplicitCoeffs
odeCoeffs odeCoeffs
{ {
solver KRR4; solver Rosenbrock43;
absTol 1e-12; absTol 1e-12;
relTol 0.01; relTol 0.01;
} }

View File

@ -33,7 +33,7 @@ EulerImplicitCoeffs
odeCoeffs odeCoeffs
{ {
solver KRR4; solver Rosenbrock43;
absTol 1e-12; absTol 1e-12;
relTol 0.01; relTol 0.01;
} }