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

Rationalisation of the ODE solvers library

Browse Source
This commit is contained in:
Henry
2013-10-31 23:51:16 +00:00
parent 30db0da817
commit dedfb01ada
14 changed files with 598 additions and 438 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
applications/test/ODE/Test-ODE.C
View File

@ -29,7 +29,6 @@ Description
#include "IOmanip.H"
#include "ODESystem.H"
#include "ODESolver.H"
#include "RK.H"
using namespace Foam;
@ -137,11 +136,10 @@ int main(int argc, char *argv[])
scalarField y(yStart);
scalarField dydx(dyStart);
scalarField yScale(ode.nEqns(), 1.0);
scalar hEst = 0.6;
scalar hDid, hNext;
odeSolver->solve(ode, x, y, dydx, eps, yScale, hEst, hDid, hNext);
odeSolver->solve(ode, x, y, dydx, eps, hEst, hDid, hNext);
Info<< scientific << setw(13) << eps;
Info<< fixed << setw(11) << hEst;

24
src/ODE/Make/files
View File

@ -1,19 +1,11 @@
ODE = ODE
ODESolvers = ODESolvers
ODESolversODESolver = ODESolvers/ODESolver
ODESolversKRR4 = ODESolvers/KRR4
ODESolversRK = ODESolvers/RK
ODESolversSIBS = ODESolvers/SIBS
ODESolvers/ODESolver/ODESolver.C
ODESolvers/ODESolver/ODESolverNew.C
$(ODESolversODESolver)/ODESolver.C
$(ODESolversODESolver)/ODESolverNew.C
$(ODESolversRK)/RK.C
$(ODESolversKRR4)/KRR4.C
$(ODESolversSIBS)/SIBS.C
$(ODESolversSIBS)/SIMPR.C
$(ODESolversSIBS)/polyExtrapolate.C
ODESolvers/adaptiveSolver/adaptiveSolver.C
ODESolvers/RKCK5/RKCK5.C
ODESolvers/KRR4/KRR4.C
ODESolvers/SIBS/SIBS.C
ODESolvers/SIBS/SIMPR.C
ODESolvers/SIBS/polyExtrapolate.C
LIB = $(FOAM_LIBBIN)/libODE

247
src/ODE/ODESolvers/KRR4/KRR4.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -31,13 +31,8 @@ License
namespace Foam
{
defineTypeNameAndDebug(KRR4, 0);
addToRunTimeSelectionTable(ODESolver, KRR4, ODE);
const scalar
KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25,
KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296;
const scalar
KRR4::gamma = 1.0/2.0,
KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
@ -58,166 +53,112 @@ const scalar
Foam::KRR4::KRR4(const ODESystem& ode)
:
ODESolver(ode),
yTemp_(n_, 0.0),
dydxTemp_(n_, 0.0),
g1_(n_, 0.0),
g2_(n_, 0.0),
g3_(n_, 0.0),
g4_(n_, 0.0),
yErr_(n_, 0.0),
dfdx_(n_, 0.0),
dfdy_(n_, n_, 0.0),
a_(n_, n_, 0.0),
pivotIndices_(n_, 0.0)
adaptiveSolver(ode),
g1_(n_),
g2_(n_),
g3_(n_),
g4_(n_),
err_(n_),
dfdx_(n_),
dfdy_(n_, n_),
a_(n_, n_),
pivotIndices_(n_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::KRR4::solve
Foam::scalar Foam::KRR4::solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y,
const scalar eps
) 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(g1_, i)
{
g1_[i] = dydx0[i] + dx*c1X*dfdx_[i];
}
LUBacksubstitute(a_, pivotIndices_, g1_);
forAll(y, i)
{
y[i] = y0[i] + a21*g1_[i];
}
ode.derivatives(x0 + a2X*dx, y, dydx_);
forAll(g2_, i)
{
g2_[i] = dydx_[i] + dx*c2X*dfdx_[i] + c21*g1_[i]/dx;
}
LUBacksubstitute(a_, pivotIndices_, g2_);
forAll(y, i)
{
y[i] = y0[i] + a31*g1_[i] + a32*g2_[i];
}
ode.derivatives(x0 + a3X*dx, y, dydx_);
forAll(g3_, i)
{
g3_[i] = dydx0[i] + dx*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, g3_);
forAll(g4_, i)
{
g4_[i] = dydx_[i] + dx*c4X*dfdx_[i]
+ (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/dx;
}
LUBacksubstitute(a_, pivotIndices_, g4_);
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];
}
return normalizeError(eps, y0, y, err_);
}
void Foam::KRR4::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
scalar& hNext
const scalar dxTry,
scalar& dxDid,
scalar& dxNext
) const
{
scalar xTemp = x;
yTemp_ = y;
dydxTemp_ = dydx;
ode.jacobian(xTemp, yTemp_, dfdx_, dfdy_);
scalar h = hTry;
for (register label jtry=0; jtry<maxtry; jtry++)
{
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*h);
}
LUDecompose(a_, pivotIndices_);
for (register label i=0; i<n_; i++)
{
g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];
}
LUBacksubstitute(a_, pivotIndices_, g1_);
for (register label i=0; i<n_; i++)
{
y[i] = yTemp_[i] + a21*g1_[i];
}
x = xTemp + a2X*h;
ode.derivatives(x, y, dydx_);
for (register label i=0; i<n_; i++)
{
g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/h;
}
LUBacksubstitute(a_, pivotIndices_, g2_);
for (register label i=0; i<n_; i++)
{
y[i] = yTemp_[i] + a31*g1_[i] + a32*g2_[i];
}
x = xTemp + a3X*h;
ode.derivatives(x, y, dydx_);
for (register label i=0; i<n_; i++)
{
g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;
}
LUBacksubstitute(a_, pivotIndices_, g3_);
for (register label i=0; i<n_; i++)
{
g4_[i] = dydx_[i] + h*c4X*dfdx_[i]
+ (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;
}
LUBacksubstitute(a_, pivotIndices_, g4_);
for (register label i=0; i<n_; i++)
{
y[i] = yTemp_[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];
yErr_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];
}
x = xTemp + h;
if (x == xTemp)
{
FatalErrorIn
(
"void Foam::KRR4::solve"
"("
"const ODESystem&, "
"scalar&, "
"scalarField&, "
"scalarField&, "
"const scalar, "
"const scalarField&, "
"const scalar, "
"scalar&, "
"scalar&"
") const"
) << "solver stalled: step size = 0"
<< exit(FatalError);
}
scalar maxErr = 0.0;
for (register label i=0; i<n_; i++)
{
maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
}
maxErr /= eps;
if (maxErr <= 1.0)
{
hDid = h;
hNext = (maxErr > errcon ? safety*h*pow(maxErr, pgrow) : grow*h);
return;
}
else
{
hNext = safety*h*pow(maxErr, pshrink);
h = (h >= 0.0 ? max(hNext, shrink*h) : min(hNext, shrink*h));
}
}
FatalErrorIn
(
"void Foam::KRR4::solve"
"("
"const ODESystem&, "
"scalar&, "
"scalarField&, "
"scalarField&, "
"const scalar, "
"const scalarField&, "
"const scalar, "
"scalar&, "
"scalar&"
") const"
) << "Maximum number of solver iterations exceeded"
<< exit(FatalError);
adaptiveSolver::solve(odes, x, y, dydx, eps, dxTry, dxDid, dxNext);
}

55
src/ODE/ODESolvers/KRR4/KRR4.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -25,11 +25,26 @@ Class
Foam::KRR4
Description
Foam::KRR4
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
SourceFiles
KRR4CK.C
KRR4QS.C
KRR4.C
\*---------------------------------------------------------------------------*/
@ -37,6 +52,7 @@ SourceFiles
#define KRR4_H
#include "ODESolver.H"
#include "adaptiveSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -49,26 +65,21 @@ namespace Foam
class KRR4
:
public ODESolver
public ODESolver,
public adaptiveSolver
{
// Private data
mutable scalarField yTemp_;
mutable scalarField dydxTemp_;
mutable scalarField g1_;
mutable scalarField g2_;
mutable scalarField g3_;
mutable scalarField g4_;
mutable scalarField yErr_;
mutable scalarField err_;
mutable scalarField dfdx_;
mutable scalarSquareMatrix dfdy_;
mutable scalarSquareMatrix a_;
mutable labelList pivotIndices_;
static const int maxtry = 40;
static const scalar safety, grow, pgrow, shrink, pshrink, errcon;
static const scalar
gamma,
a21, a31, a32,
@ -93,6 +104,19 @@ public:
// 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 scalar eps
) const;
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
@ -100,10 +124,9 @@ public:
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
scalar& hNext
const scalar dxTry,
scalar& dxDid,
scalar& dxNext
) const;
};

8
src/ODE/ODESolvers/ODESolver/ODESolver.C
View File

@ -39,7 +39,6 @@ namespace Foam
Foam::ODESolver::ODESolver(const ODESystem& ode)
:
n_(ode.nEqns()),
yScale_(n_),
dydx_(n_)
{}
@ -67,11 +66,6 @@ void Foam::ODESolver::solve
{
ode.derivatives(x, y, dydx_);
for (label i=0; i<n_; i++)
{
yScale_[i] = mag(y[i]) + mag(dydx_[i]*h) + SMALL;
}
if ((x + h - xEnd)*(x + h - xStart) > 0.0)
{
h = xEnd - x;
@ -80,7 +74,7 @@ void Foam::ODESolver::solve
hNext = 0;
scalar hDid;
solve(ode, x, y, dydx_, eps, yScale_, h, hDid, hNext);
solve(ode, x, y, dydx_, eps, h, hDid, hNext);
if ((x - xEnd)*(xEnd - xStart) >= 0.0)
{

2
src/ODE/ODESolvers/ODESolver/ODESolver.H
View File

@ -56,7 +56,6 @@ protected:
// Private data
label n_;
mutable scalarField yScale_;
mutable scalarField dydx_;
@ -117,7 +116,6 @@ public:
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
scalar& hNext

204
src/ODE/ODESolvers/RK/RK.C
View File

@ -1,204 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 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 "RK.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(RK, 0);
addToRunTimeSelectionTable(ODESolver, RK, ODE);
const scalar
RK::safety=0.9, RK::pGrow=-0.2, RK::pShrink=-0.25, RK::errCon=1.89e-4;
const scalar
RK::a2 = 0.2, RK::a3 = 0.3, RK::a4 = 0.6, RK::a5 = 1.0, RK::a6 = 0.875,
RK::b21 = 0.2, RK::b31 = 3.0/40.0, RK::b32 = 9.0/40.0,
RK::b41 = 0.3, RK::b42 = -0.9, RK::b43 = 1.2,
RK::b51 = -11.0/54.0, RK::b52 = 2.5, RK::b53 = -70.0/27.0,
RK::b54 = 35.0/27.0,
RK::b61 = 1631.0/55296.0, RK::b62 = 175.0/512.0, RK::b63 = 575.0/13824.0,
RK::b64 = 44275.0/110592.0, RK::b65 = 253.0/4096.0,
RK::c1 = 37.0/378.0, RK::c3 = 250.0/621.0,
RK::c4 = 125.0/594.0, RK::c6 = 512.0/1771.0,
RK::dc1 = RK::c1 - 2825.0/27648.0, RK::dc3 = RK::c3 - 18575.0/48384.0,
RK::dc4 = RK::c4 - 13525.0/55296.0, RK::dc5 = -277.00/14336.0,
RK::dc6 = RK::c6 - 0.25;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::RK::RK(const ODESystem& ode)
:
ODESolver(ode),
yTemp_(n_, 0.0),
ak2_(n_, 0.0),
ak3_(n_, 0.0),
ak4_(n_, 0.0),
ak5_(n_, 0.0),
ak6_(n_, 0.0),
yErr_(n_, 0.0),
yTemp2_(n_, 0.0)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::RK::solve
(
const ODESystem& ode,
const scalar x,
const scalarField& y,
const scalarField& dydx,
const scalar h,
scalarField& yout,
scalarField& yerr
) const
{
forAll(yTemp_, i)
{
yTemp_[i] = y[i] + b21*h*dydx[i];
}
ode.derivatives(x + a2*h, yTemp_, ak2_);
forAll(yTemp_, i)
{
yTemp_[i] = y[i] + h*(b31*dydx[i] + b32*ak2_[i]);
}
ode.derivatives(x + a3*h, yTemp_, ak3_);
forAll(yTemp_, i)
{
yTemp_[i] = y[i] + h*(b41*dydx[i] + b42*ak2_[i] + b43*ak3_[i]);
}
ode.derivatives(x + a4*h, yTemp_, ak4_);
forAll(yTemp_, i)
{
yTemp_[i] = y[i]
+ h*(b51*dydx[i] + b52*ak2_[i] + b53*ak3_[i] + b54*ak4_[i]);
}
ode.derivatives(x + a5*h, yTemp_, ak5_);
forAll(yTemp_, i)
{
yTemp_[i] = y[i]
+ h*
(
b61*dydx[i] + b62*ak2_[i] + b63*ak3_[i]
+ b64*ak4_[i] + b65*ak5_[i]
);
}
ode.derivatives(x + a6*h, yTemp_, ak6_);
forAll(yout, i)
{
yout[i] = y[i]
+ h*(c1*dydx[i] + c3*ak3_[i] + c4*ak4_[i] + c6*ak6_[i]);
}
forAll(yerr, i)
{
yerr[i] =
h*
(
dc1*dydx[i] + dc3*ak3_[i] + dc4*ak4_[i]
+ dc5*ak5_[i] + dc6*ak6_[i]
);
}
}
void Foam::RK::solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
scalar& hNext
) const
{
scalar h = hTry;
scalar maxErr = 0.0;
for (;;)
{
solve(ode, x, y, dydx, h, yTemp2_, yErr_);
maxErr = 0.0;
for (register label i=0; i<n_; i++)
{
maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
}
maxErr /= eps;
if (maxErr <= 1.0)
{
break;
}
{
scalar hTemp = safety*h*pow(maxErr, pShrink);
h = (h >= 0.0 ? max(hTemp, 0.1*h) : min(hTemp, 0.1*h));
}
if (h < VSMALL)
{
FatalErrorIn("RK::solve")
<< "stepsize underflow"
<< exit(FatalError);
}
}
hDid = h;
x += h;
y = yTemp2_;
if (maxErr > errCon)
{
hNext = safety*h*pow(maxErr, pGrow);
}
else
{
hNext = 5.0*h;
}
}
// ************************************************************************* //

155
src/ODE/ODESolvers/RKCK5/RKCK5.C Normal file
View File

@ -0,0 +1,155 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "RKCK5.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(RKCK5, 0);
addToRunTimeSelectionTable(ODESolver, RKCK5, ODE);
const scalar
RKCK5::a2 = 0.2, RKCK5::a3 = 0.3, RKCK5::a4 = 0.6, RKCK5::a5 = 1.0,
RKCK5::a6 = 0.875,
RKCK5::b21 = 0.2, RKCK5::b31 = 3.0/40.0, RKCK5::b32 = 9.0/40.0,
RKCK5::b41 = 0.3, RKCK5::b42 = -0.9, RKCK5::b43 = 1.2,
RKCK5::b51 = -11.0/54.0, RKCK5::b52 = 2.5, RKCK5::b53 = -70.0/27.0,
RKCK5::b54 = 35.0/27.0,
RKCK5::b61 = 1631.0/55296.0, RKCK5::b62 = 175.0/512.0,
RKCK5::b63 = 575.0/13824.0, RKCK5::b64 = 44275.0/110592.0,
RKCK5::b65 = 253.0/4096.0,
RKCK5::c1 = 37.0/378.0, RKCK5::c3 = 250.0/621.0,
RKCK5::c4 = 125.0/594.0, RKCK5::c6 = 512.0/1771.0,
RKCK5::dc1 = RKCK5::c1 - 2825.0/27648.0,
RKCK5::dc3 = RKCK5::c3 - 18575.0/48384.0,
RKCK5::dc4 = RKCK5::c4 - 13525.0/55296.0, RKCK5::dc5 = -277.00/14336.0,
RKCK5::dc6 = RKCK5::c6 - 0.25;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::RKCK5::RKCK5(const ODESystem& ode)
:
ODESolver(ode),
adaptiveSolver(ode),
yTemp_(n_),
k2_(n_),
k3_(n_),
k4_(n_),
k5_(n_),
k6_(n_),
err_(n_)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::RKCK5::solve
(
const ODESystem& odes,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y,
const scalar eps
) const
{
forAll(yTemp_, i)
{
yTemp_[i] = y0[i] + b21*dx*dydx0[i];
}
odes.derivatives(x0 + a2*dx, yTemp_, k2_);
forAll(yTemp_, i)
{
yTemp_[i] = y0[i] + dx*(b31*dydx0[i] + b32*k2_[i]);
}
odes.derivatives(x0 + a3*dx, yTemp_, k3_);
forAll(yTemp_, i)
{
yTemp_[i] = y0[i] + dx*(b41*dydx0[i] + b42*k2_[i] + b43*k3_[i]);
}
odes.derivatives(x0 + a4*dx, yTemp_, k4_);
forAll(yTemp_, i)
{
yTemp_[i] = y0[i]
+ dx*(b51*dydx0[i] + b52*k2_[i] + b53*k3_[i] + b54*k4_[i]);
}
odes.derivatives(x0 + a5*dx, yTemp_, k5_);
forAll(yTemp_, i)
{
yTemp_[i] = y0[i]
+ dx
*(b61*dydx0[i] + b62*k2_[i] + b63*k3_[i] + b64*k4_[i] + b65*k5_[i]);
}
odes.derivatives(x0 + a6*dx, yTemp_, k6_);
forAll(y, i)
{
y[i] = y0[i]
+ dx*(c1*dydx0[i] + c3*k3_[i] + c4*k4_[i] + c6*k6_[i]);
}
forAll(err_, i)
{
err_[i] =
dx
*(dc1*dydx0[i] + dc3*k3_[i] + dc4*k4_[i] + dc5*k5_[i] + dc6*k6_[i]);
}
return normalizeError(eps, y0, y, err_);
}
void Foam::RKCK5::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalar dxTry,
scalar& dxDid,
scalar& dxNext
) const
{
adaptiveSolver::solve(odes, x, y, dydx, eps, dxTry, dxDid, dxNext);
}
// ************************************************************************* //

77
src/ODE/ODESolvers/RK/RK.H → src/ODE/ODESolvers/RKCK5/RKCK5.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,21 +22,30 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::RK
Foam::RKCK5
Description
Runge-Kutta ODE solver
5th Order Cash-Karp Runge-Kutta ODE solver
References:
\verbatim
"A variable order Runge-Kutta method for initial value problems
with rapidly varying right-hand sides"
Cash, J.R.,
Karp, A.H.
ACM Transactions on Mathematical Software, vol. 16, 1990, pp. 201–222.
\endverbatim
SourceFiles
RKCK.C
RKQS.C
RKCK5.C
\*---------------------------------------------------------------------------*/
#ifndef RK_H
#define RK_H
#ifndef RKCK5_H
#define RKCK5_H
#include "ODESolver.H"
#include "adaptiveSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -44,16 +53,17 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class RK Declaration
Class RKCK5 Declaration
\*---------------------------------------------------------------------------*/
class RK
class RKCK5
:
public ODESolver
public ODESolver,
public adaptiveSolver
{
// Private data
static const scalar safety, pGrow, pShrink, errCon;
//- RKCK Constants
static const scalar
a2, a3, a4, a5, a6,
b21, b31, b32, b41, b42, b43,
@ -61,43 +71,45 @@ class RK
c1, c3, c4, c6,
dc1, dc3, dc4, dc5, dc6;
// Temporary fields
mutable scalarField yTemp_;
mutable scalarField ak2_;
mutable scalarField ak3_;
mutable scalarField ak4_;
mutable scalarField ak5_;
mutable scalarField ak6_;
mutable scalarField k2_;
mutable scalarField k3_;
mutable scalarField k4_;
mutable scalarField k5_;
mutable scalarField k6_;
//- Error-estimate field
mutable scalarField err_;
mutable scalarField yErr_;
mutable scalarField yTemp2_;
public:
//- Runtime type information
TypeName("RK");
TypeName("RKCK5");
// Constructors
//- Construct from ODE
RK(const ODESystem& ode);
RKCK5(const ODESystem& ode);
// Member Functions
void solve
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x,
const scalarField& y,
const scalarField& dydx,
const scalar h,
scalarField& yout,
scalarField& yerr
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y,
const scalar eps
) const;
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
@ -105,10 +117,9 @@ public:
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
scalar& hNext
const scalar dxTry,
scalar& dxDid,
scalar& dxNext
) const;
};

7
src/ODE/ODESolvers/SIBS/SIBS.C
View File

@ -75,7 +75,6 @@ void Foam::SIBS::solve
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
scalar& hNext
@ -164,7 +163,11 @@ void Foam::SIBS::solve
maxErr = SMALL;
for (register label i=0; i<n_; i++)
{
maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
maxErr = max
(
maxErr,
mag(yErr_[i])/(mag(yTemp_[i]) + SMALL)
);
}
maxErr /= eps;
km = k - 1;

1
src/ODE/ODESolvers/SIBS/SIBS.H
View File

@ -122,7 +122,6 @@ public:
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
scalar& hNext

123
src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C Normal file
View File

@ -0,0 +1,123 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "adaptiveSolver.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
const scalar
adaptiveSolver::safeScale=0.9,
adaptiveSolver::alphaInc=0.2,
adaptiveSolver::alphaDec=0.25,
adaptiveSolver::minScale=0.2,
adaptiveSolver::maxScale=10;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::adaptiveSolver::adaptiveSolver(const ODESystem& ode)
:
yTemp_(ode.nEqns())
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::adaptiveSolver::normalizeError
(
const scalar eps,
const scalarField& y0,
const scalarField& y,
const scalarField& err
) const
{
scalar epsAbs = 1e-10;
scalar epsRel = eps;
// Calculate the maximum error
scalar maxErr = 0.0;
forAll(err, i)
{
scalar eps = epsAbs + epsRel*max(mag(y0[i]), mag(y[i]));
maxErr = max(maxErr, mag(err[i])/eps);
}
return maxErr;
}
void Foam::adaptiveSolver::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalar dxTry,
scalar& dxDid,
scalar& dxNext
) const
{
scalar dx = dxTry;
scalar err = 0.0;
// Loop over solver and adjust step-size as necessary
// to achieve desired error
do
{
// Solve step and provide error estimate
err = solve(odes, x, y, dydx, dx, yTemp_, eps);
// If error is large reduce dx
if (err > 1)
{
scalar scale = max(safeScale*pow(err, -alphaDec), minScale);
dx *= scale;
if (dx < VSMALL)
{
FatalErrorIn("adaptiveSolver::solve")
<< "stepsize underflow"
<< exit(FatalError);
}
}
} while (err > 1);
// Update the state
dxDid = dx;
x += dx;
y = yTemp_;
// If the error is small increase the step-size
dxNext = min(max(safeScale*pow(err, -alphaInc), minScale), maxScale)*dx;
}
// ************************************************************************* //

127
src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H Normal file
View File

@ -0,0 +1,127 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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::adaptiveSolver
Description
5th Order Cash-Karp Runge-Kutta ODE solver
References:
\verbatim
"A variable order Runge-Kutta method for initial value problems
with rapidly varying right-hand sides"
Cash, J.R.,
Karp, A.H.
ACM Transactions on Mathematical Software, vol. 16, 1990, pp. 201–222.
\endverbatim
SourceFiles
adaptiveSolver.C
\*---------------------------------------------------------------------------*/
#ifndef adaptiveSolver_H
#define adaptiveSolver_H
#include "ODESolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class adaptiveSolver Declaration
\*---------------------------------------------------------------------------*/
class adaptiveSolver
{
// Private data
//- Step-size adjustment controls
static const scalar safeScale, alphaInc, alphaDec, minScale, maxScale;
mutable scalarField yTemp_;
public:
// Constructors
//- Construct from ODESystem
adaptiveSolver(const ODESystem& ode);
//- Destructor
virtual ~adaptiveSolver()
{}
// Member Functions
//- Return the nomalized scalar error
scalar normalizeError
(
const scalar eps,
const scalarField& y0,
const scalarField& y,
const scalarField& err
) const;
//- Solve a single step dx and return the error
virtual scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
const scalar dx,
scalarField& y,
const scalar eps
) const = 0;
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalarField& dydx,
const scalar eps,
const scalar dxTry,
scalar& dxDid,
scalar& dxNext
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/ODE/ODESystem/ODESystem.H
View File

@ -51,7 +51,7 @@ public:
// Constructors
//- Construct null8
//- Construct null
ODESystem()
{}