From dedfb01ada052c831b9c67c1e96e4855b9b8cc99 Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 31 Oct 2013 23:51:16 +0000 Subject: [PATCH] Rationalisation of the ODE solvers library --- applications/test/ODE/Test-ODE.C | 4 +- src/ODE/Make/files | 24 +- src/ODE/ODESolvers/KRR4/KRR4.C | 247 +++++++----------- src/ODE/ODESolvers/KRR4/KRR4.H | 55 ++-- src/ODE/ODESolvers/ODESolver/ODESolver.C | 8 +- src/ODE/ODESolvers/ODESolver/ODESolver.H | 2 - src/ODE/ODESolvers/RK/RK.C | 204 --------------- src/ODE/ODESolvers/RKCK5/RKCK5.C | 155 +++++++++++ src/ODE/ODESolvers/{RK/RK.H => RKCK5/RKCK5.H} | 77 +++--- src/ODE/ODESolvers/SIBS/SIBS.C | 7 +- src/ODE/ODESolvers/SIBS/SIBS.H | 1 - .../adaptiveSolver/adaptiveSolver.C | 123 +++++++++ .../adaptiveSolver/adaptiveSolver.H | 127 +++++++++ src/ODE/ODESystem/ODESystem.H | 2 +- 14 files changed, 598 insertions(+), 438 deletions(-) delete mode 100644 src/ODE/ODESolvers/RK/RK.C create mode 100644 src/ODE/ODESolvers/RKCK5/RKCK5.C rename src/ODE/ODESolvers/{RK/RK.H => RKCK5/RKCK5.H} (63%) create mode 100644 src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C create mode 100644 src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C index e5d3026828..7848ca0d8b 100644 --- a/applications/test/ODE/Test-ODE.C +++ b/applications/test/ODE/Test-ODE.C @@ -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; diff --git a/src/ODE/Make/files b/src/ODE/Make/files index 0e6c49a915..cb569e78f5 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -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 diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C index e9218641a7..6f0b938782 100644 --- a/src/ODE/ODESolvers/KRR4/KRR4.C +++ b/src/ODE/ODESolvers/KRR4/KRR4.C @@ -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 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); } diff --git a/src/ODE/ODESolvers/KRR4/KRR4.H b/src/ODE/ODESolvers/KRR4/KRR4.H index 1eb434a0ba..710651fde9 100644 --- a/src/ODE/ODESolvers/KRR4/KRR4.H +++ b/src/ODE/ODESolvers/KRR4/KRR4.H @@ -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; }; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index 5f7188a079..60a599e8d2 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.C +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C @@ -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 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) { diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H index 04315ea7dc..a293088685 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H @@ -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 diff --git a/src/ODE/ODESolvers/RK/RK.C b/src/ODE/ODESolvers/RK/RK.C deleted file mode 100644 index 06c5e6ae81..0000000000 --- a/src/ODE/ODESolvers/RK/RK.C +++ /dev/null @@ -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 . - -\*---------------------------------------------------------------------------*/ - -#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= 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; - } -} - - -// ************************************************************************* // diff --git a/src/ODE/ODESolvers/RKCK5/RKCK5.C b/src/ODE/ODESolvers/RKCK5/RKCK5.C new file mode 100644 index 0000000000..5050c0e7ca --- /dev/null +++ b/src/ODE/ODESolvers/RKCK5/RKCK5.C @@ -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 . + +\*---------------------------------------------------------------------------*/ + +#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); +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/RK/RK.H b/src/ODE/ODESolvers/RKCK5/RKCK5.H similarity index 63% rename from src/ODE/ODESolvers/RK/RK.H rename to src/ODE/ODESolvers/RKCK5/RKCK5.H index 20fa3aa789..73a33e36bd 100644 --- a/src/ODE/ODESolvers/RK/RK.H +++ b/src/ODE/ODESolvers/RKCK5/RKCK5.H @@ -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 . 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; }; diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C index 95ab44776c..964e288949 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.C +++ b/src/ODE/ODESolvers/SIBS/SIBS.C @@ -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. + +\*---------------------------------------------------------------------------*/ + +#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; +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H new file mode 100644 index 0000000000..40f07c2d93 --- /dev/null +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H @@ -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 . + +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 + +// ************************************************************************* // diff --git a/src/ODE/ODESystem/ODESystem.H b/src/ODE/ODESystem/ODESystem.H index 90da37c427..50b3fbffed 100644 --- a/src/ODE/ODESystem/ODESystem.H +++ b/src/ODE/ODESystem/ODESystem.H @@ -51,7 +51,7 @@ public: // Constructors - //- Construct null8 + //- Construct null ODESystem() {}