zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ODE solvers: Added Rosenbrock43 (equivalent to KRR43 but more consistent)

Browse Source
This commit is contained in:
Henry
2013-11-03 23:21:32 +00:00
parent 16d07b77da
commit 734d88a875
8 changed files with 407 additions and 61 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
src/ODE/Make/files
View File

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

89
src/ODE/ODESolvers/KRR4/KRR4.C → src/ODE/ODESolvers/KRR43/KRR43.C
View File

@ -23,41 +23,54 @@ License
\*---------------------------------------------------------------------------*/
#include "KRR4.H"
#include "KRR43.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(KRR4, 0);
addToRunTimeSelectionTable(ODESolver, KRR4, dictionary);
defineTypeNameAndDebug(KRR43, 0);
addToRunTimeSelectionTable(ODESolver, KRR43, dictionary);
const scalar
KRR4::gamma = 1.0/2.0,
KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0,
KRR4::c41 = -112.0/125.0, KRR4::c42 = -54.0/125.0, KRR4::c43 = -2.0/5.0,
KRR4::b1 = 19.0/9.0, KRR4::b2 = 1.0/2.0, KRR4::b3 = 25.0/108.0,
KRR4::b4 = 125.0/108.0,
KRR4::e1 = 17.0/54.0, KRR4::e2 = 7.0/36.0, KRR4::e3 = 0.0,
KRR4::e4 = 125.0/108.0,
KRR4::c1X = 1.0/2.0, KRR4::c2X = -3.0/2.0, KRR4::c3X = 121.0/50.0,
KRR4::c4X = 29.0/250.0,
KRR4::a2X = 1.0, KRR4::a3X = 3.0/5.0;
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::KRR4::KRR4(const ODESystem& ode, const dictionary& dict)
Foam::KRR43::KRR43(const ODESystem& ode, const dictionary& dict)
:
ODESolver(ode, dict),
adaptiveSolver(ode, dict),
g1_(n_),
g2_(n_),
g3_(n_),
g4_(n_),
k1_(n_),
k2_(n_),
k3_(n_),
k4_(n_),
err_(n_),
dydx_(n_),
dfdx_(n_),
@ -69,7 +82,7 @@ Foam::KRR4::KRR4(const ODESystem& ode, const dictionary& dict)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::KRR4::solve
Foam::scalar Foam::KRR43::solve
(
const ODESystem& ode,
const scalar x0,
@ -93,60 +106,60 @@ Foam::scalar Foam::KRR4::solve
LUDecompose(a_, pivotIndices_);
forAll(g1_, i)
forAll(k1_, i)
{
g1_[i] = dydx0[i] + dx*c1X*dfdx_[i];
k1_[i] = dydx0[i] + dx*c1X*dfdx_[i];
}
LUBacksubstitute(a_, pivotIndices_, g1_);
LUBacksubstitute(a_, pivotIndices_, k1_);
forAll(y, i)
{
y[i] = y0[i] + a21*g1_[i];
y[i] = y0[i] + a21*k1_[i];
}
ode.derivatives(x0 + a2X*dx, y, dydx_);
forAll(g2_, i)
forAll(k2_, i)
{
g2_[i] = dydx_[i] + dx*c2X*dfdx_[i] + c21*g1_[i]/dx;
k2_[i] = dydx_[i] + dx*c2X*dfdx_[i] + c21*k1_[i]/dx;
}
LUBacksubstitute(a_, pivotIndices_, g2_);
LUBacksubstitute(a_, pivotIndices_, k2_);
forAll(y, i)
{
y[i] = y0[i] + a31*g1_[i] + a32*g2_[i];
y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
}
ode.derivatives(x0 + a3X*dx, y, dydx_);
forAll(g3_, i)
forAll(k3_, i)
{
g3_[i] = dydx0[i] + dx*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/dx;
k3_[i] = dydx_[i] + dx*c3X*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, g3_);
LUBacksubstitute(a_, pivotIndices_, k3_);
forAll(g4_, i)
forAll(k4_, i)
{
g4_[i] = dydx_[i] + dx*c4X*dfdx_[i]
+ (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/dx;
k4_[i] = dydx_[i] + dx*c4X*dfdx_[i]
+ (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, g4_);
LUBacksubstitute(a_, pivotIndices_, k4_);
forAll(y, i)
{
y[i] = y0[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];
err_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[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::KRR4::solve
void Foam::KRR43::solve
(
const ODESystem& odes,
scalar& x,

24
src/ODE/ODESolvers/KRR4/KRR4.H → src/ODE/ODESolvers/KRR43/KRR43.H
View File

@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::KRR4
Foam::KRR43
Description
4th Order Kaps Rentrop Rosenbrock stiff ODE solver
@ -44,12 +44,12 @@ Description
\endverbatim
SourceFiles
KRR4.C
KRR43.C
\*---------------------------------------------------------------------------*/
#ifndef KRR4_H
#define KRR4_H
#ifndef KRR43_H
#define KRR43_H
#include "ODESolver.H"
#include "adaptiveSolver.H"
@ -60,20 +60,20 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class KRR4 Declaration
Class KRR43 Declaration
\*---------------------------------------------------------------------------*/
class KRR4
class KRR43
:
public ODESolver,
public adaptiveSolver
{
// Private data
mutable scalarField g1_;
mutable scalarField g2_;
mutable scalarField g3_;
mutable scalarField g4_;
mutable scalarField k1_;
mutable scalarField k2_;
mutable scalarField k3_;
mutable scalarField k4_;
mutable scalarField err_;
mutable scalarField dydx_;
mutable scalarField dfdx_;
@ -94,13 +94,13 @@ class KRR4
public:
//- Runtime type information
TypeName("KRR4");
TypeName("KRR43");
// Constructors
//- Construct from ODE
KRR4(const ODESystem& ode, const dictionary& dict);
KRR43(const ODESystem& ode, const dictionary& dict);
// Member Functions

3
src/ODE/ODESolvers/RKCK45/RKCK45.H
View File

@ -34,7 +34,10 @@ Description
Cash, J.R.,
Karp, A.H.
ACM Transactions on Mathematical Software, vol. 16, 1990, pp. 201–222.
\endverbatim
Based on code from:
\verbatim
"Solving Ordinary Differential Equations I: Nonstiff Problems,
second edition",
Hairer, E.,

3
src/ODE/ODESolvers/RKDP45/RKDP45.H
View File

@ -34,7 +34,10 @@ Description
Prince, P. J.,
Journal of Computational and Applied Mathematics,
6 (1), 1980: pp. 19-26.
\endverbatim
Based on code from:
\verbatim
"Solving Ordinary Differential Equations I: Nonstiff Problems,
second edition",
Hairer, E.,

13
src/ODE/ODESolvers/RKF45/RKF45.H
View File

@ -34,6 +34,8 @@ Description
Fehlberg, E.,
NASA Technical Report 315, 1969.
Based on code from:
\verbatim
"Solving Ordinary Differential Equations I: Nonstiff Problems,
second edition",
Hairer, E.,
@ -69,16 +71,7 @@ Description
The difference between the 5th and 4th order (\f$\epsilon = \left|\Delta_5
- \Delta_4\right|\f$) can be used to determine if the stepsize \f$h\f$ needs
to be adjusted. The step-size does not need to be reduced if the following
inequation holds true:
\f[
\left|\Delta_5 - \Delta_4\right| < \epsilon_{abs} +
\epsilon_{rel}\left|y\right|
\f]
where \f$\left|y\right|\f$ is the element wise maximum of
\f$\left|y_k\right|\f$ and \f$\left|y_{k+1}\right|\f$.
to be adjusted.
SourceFiles
RKF45.C

184
src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C Normal file
View File

@ -0,0 +1,184 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "Rosenbrock43.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(Rosenbrock43, 0);
addToRunTimeSelectionTable(ODESolver, Rosenbrock43, dictionary);
const scalar
Rosenbrock43::a21 = 2.0,
Rosenbrock43::a31 = 48.0/25.0,
Rosenbrock43::a32 = 6.0/25.0,
Rosenbrock43::c21 = -8.0,
Rosenbrock43::c31 = 372.0/25.0,
Rosenbrock43::c32 = 12.0/5.0,
Rosenbrock43::c41 = -112.0/125.0,
Rosenbrock43::c42 = -54.0/125.0,
Rosenbrock43::c43 = -2.0/5.0,
Rosenbrock43::b1 = 19.0/9.0,
Rosenbrock43::b2 = 1.0/2.0,
Rosenbrock43::b3 = 25.0/108.0,
Rosenbrock43::b4 = 125.0/108.0,
Rosenbrock43::e1 = 34.0/108.0,
Rosenbrock43::e2 = 7.0/36.0,
Rosenbrock43::e3 = 0.0,
Rosenbrock43::e4 = 125.0/108.0,
Rosenbrock43::gamma = 0.5,
Rosenbrock43::c2 = 1.0,
Rosenbrock43::c3 = 3.0/5.0,
Rosenbrock43::d1 = 1.0/2.0,
Rosenbrock43::d2 = -3.0/2.0,
Rosenbrock43::d3 = 605.0/250.0,
Rosenbrock43::d4 = 29.0/250.0;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::Rosenbrock43::Rosenbrock43(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::Rosenbrock43::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_);
// compute k1:
forAll(k1_, i)
{
k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
}
LUBacksubstitute(a_, pivotIndices_, k1_);
// compute 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_);
// compute 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_);
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_);
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] + e4*k4_[i];
}
return normalizeError(y0, y, err_);
}
void Foam::Rosenbrock43::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
}
// ************************************************************************* //

149
src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.H Normal file
View File

@ -0,0 +1,149 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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/>.
Class
Foam::Rosenbrock43
Description
4th Order Kaps Rentrop Rosenbrock stiff ODE solver
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.
\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
Rosenbrock43.C
\*---------------------------------------------------------------------------*/
#ifndef Rosenbrock43_H
#define Rosenbrock43_H
#include "ODESolver.H"
#include "adaptiveSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class Rosenbrock43 Declaration
\*---------------------------------------------------------------------------*/
class Rosenbrock43
:
public ODESolver,
public adaptiveSolver
{
// Private data
mutable scalarField k1_;
mutable scalarField k2_;
mutable scalarField k3_;
mutable scalarField k4_;
mutable scalarField err_;
mutable scalarField dydx_;
mutable scalarField dfdx_;
mutable scalarSquareMatrix dfdy_;
mutable scalarSquareMatrix a_;
mutable labelList pivotIndices_;
static const scalar
a21, a31, a32,
c21, c31, c32,
c41, c42, c43,
b1, b2, b3, b4,
e1, e2, e3, e4,
gamma,
c2, c3,
d1, d2, d3, d4;
public:
//- Runtime type information
TypeName("Rosenbrock43");
// Constructors
//- Construct from ODE
Rosenbrock43(const ODESystem& ode, const dictionary& dict);
// Member Functions
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y
) const;
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //