ODE solvers: further rationalisation and the addition of the Runge-Kutta-Fehlberg method

applications/test/ODE/Test-ODE.C

@@ -109,8 +109,11 @@ int main(int argc, char *argv[])
     // Create the ODE system
     testODE ode;
 
+    dictionary dict;
+    dict.add("solver", args[1]);
+
     // Create the selected ODE system solver
-    autoPtr<ODESolver> odeSolver = ODESolver::New(args[1], ode);
+    autoPtr<ODESolver> odeSolver = ODESolver::New(ode, dict);
 
     // Initialise the ODE system fields
     scalar xStart = 1.0;
@@ -124,26 +127,26 @@ int main(int argc, char *argv[])
     scalarField dyStart(ode.nEqns());
     ode.derivatives(xStart, yStart, dyStart);
 
-    Info<< setw(10) << "eps" << setw(12) << "hEst";
+    Info<< setw(10) << "relTol" << setw(12) << "dxEst";
-    Info<< setw(13) << "hDid" << setw(14) << "hNext" << endl;
+    Info<< setw(13) << "dxDid" << setw(14) << "dxNext" << endl;
     Info<< setprecision(6);
 
     for (label i=0; i<15; i++)
     {
-        scalar eps = ::Foam::exp(-scalar(i + 1));
+        scalar relTol = ::Foam::exp(-scalar(i + 1));
 
         scalar x = xStart;
         scalarField y(yStart);
-        scalarField dydx(dyStart);
 
-        scalar hEst = 0.6;
+        scalar dxEst = 0.6;
-        scalar hDid, hNext;
+        scalar dxNext = dxEst;
 
-        odeSolver->solve(ode, x, y, dydx, eps, hEst, hDid, hNext);
+        odeSolver->relTol() = relTol;
+        odeSolver->solve(ode, x, y, dxNext);
 
-        Info<< scientific << setw(13) << eps;
+        Info<< scientific << setw(13) << relTol;
-        Info<< fixed << setw(11) << hEst;
+        Info<< fixed << setw(11) << dxEst;
-        Info<< setw(13) << hDid << setw(13) << hNext
+        Info<< setw(13) << x - xStart << setw(13) << dxNext
             << setw(13) << y[0] << setw(13) << y[1]
             << setw(13) << y[2] << setw(13) << y[3]
             << endl;
@@ -159,12 +162,13 @@ int main(int argc, char *argv[])
     yEnd[2] = ::Foam::jn(2, xEnd);
     yEnd[3] = ::Foam::jn(3, xEnd);
 
-    scalar hEst = 0.5;
+    scalar dxEst = 0.5;
 
-    odeSolver->solve(ode, x, xEnd, y, 1e-4, hEst);
+    odeSolver->relTol() = 1e-4;
+    odeSolver->solve(ode, x, xEnd, y, dxEst);
 
     Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
-    Info      << "Numerical:  y(2.0) = " << y << ", hEst = " << hEst << endl;
+    Info      << "Numerical:  y(2.0) = " << y << ", dxEst = " << dxEst << endl;
 
     Info<< "\nEnd\n" << endl;
 

src/ODE/Make/files

@@ -2,7 +2,8 @@ ODESolvers/ODESolver/ODESolver.C
 ODESolvers/ODESolver/ODESolverNew.C
 
 ODESolvers/adaptiveSolver/adaptiveSolver.C
-ODESolvers/RKCK5/RKCK5.C
+ODESolvers/RKF45/RKF45.C
+ODESolvers/RKCK45/RKCK45.C
 ODESolvers/KRR4/KRR4.C
 ODESolvers/SIBS/SIBS.C
 ODESolvers/SIBS/SIMPR.C

src/ODE/ODESolvers/KRR4/KRR4.C

@@ -31,7 +31,7 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(KRR4, 0);
-    addToRunTimeSelectionTable(ODESolver, KRR4, ODE);
+    addToRunTimeSelectionTable(ODESolver, KRR4, dictionary);
 
 const scalar
     KRR4::gamma = 1.0/2.0,
@@ -50,15 +50,16 @@ const scalar
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::KRR4::KRR4(const ODESystem& ode)
+Foam::KRR4::KRR4(const ODESystem& ode, const dictionary& dict)
 :
-    ODESolver(ode),
+    ODESolver(ode, dict),
-    adaptiveSolver(ode),
+    adaptiveSolver(ode, dict),
     g1_(n_),
     g2_(n_),
     g3_(n_),
     g4_(n_),
     err_(n_),
+    dydx_(n_),
     dfdx_(n_),
     dfdy_(n_, n_),
     a_(n_, n_),
@@ -75,8 +76,7 @@ Foam::scalar Foam::KRR4::solve
     const scalarField& y0,
     const scalarField& dydx0,
     const scalar dx,
-    scalarField& y,
+    scalarField& y
-    const scalar eps
 ) const
 {
     ode.jacobian(x0, y0, dfdx_, dfdy_);
@@ -142,7 +142,7 @@ Foam::scalar Foam::KRR4::solve
         err_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];
     }
 
-    return normalizeError(eps, y0, y, err_);
+    return normalizeError(y0, y, err_);
 }
 
 
@@ -151,14 +151,10 @@ void Foam::KRR4::solve
     const ODESystem& odes,
     scalar& x,
     scalarField& y,
-    scalarField& dydx,
+    scalar& dxTry
-    const scalar eps,
-    const scalar dxTry,
-    scalar& dxDid,
-    scalar& dxNext
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dydx, eps, dxTry, dxDid, dxNext);
+    adaptiveSolver::solve(odes, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/KRR4/KRR4.H

@@ -75,6 +75,7 @@ class KRR4
         mutable scalarField g3_;
         mutable scalarField g4_;
         mutable scalarField err_;
+        mutable scalarField dydx_;
         mutable scalarField dfdx_;
         mutable scalarSquareMatrix dfdy_;
         mutable scalarSquareMatrix a_;
@@ -99,7 +100,7 @@ public:
     // Constructors
 
         //- Construct from ODE
-        KRR4(const ODESystem& ode);
+        KRR4(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -112,8 +113,7 @@ public:
             const scalarField& y0,
             const scalarField& dydx0,
             const scalar dx,
-            scalarField& y,
+            scalarField& y
-            const scalar eps
         ) const;
 
         //- Solve the ODE system and the update the state
@@ -122,11 +122,7 @@ public:
             const ODESystem& ode,
             scalar& x,
             scalarField& y,
-            scalarField& dydx,
+            scalar& dxTry
-            const scalar eps,
-            const scalar dxTry,
-            scalar& dxDid,
-            scalar& dxNext
         ) const;
 };
 

src/ODE/ODESolvers/ODESolver/ODESolver.C

@@ -30,76 +30,103 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(ODESolver, 0);
-    defineRunTimeSelectionTable(ODESolver, ODE);
+    defineRunTimeSelectionTable(ODESolver, dictionary);
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::ODESolver::ODESolver(const ODESystem& ode)
+Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict)
 :
     n_(ode.nEqns()),
-    dydx_(n_)
+    absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
+    relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-6)),
+    maxSteps_(10000)
+{}
+
+
+Foam::ODESolver::ODESolver
+(
+    const ODESystem& ode,
+    const scalarField& absTol,
+    const scalarField& relTol
+)
+:
+    n_(ode.nEqns()),
+    absTol_(absTol),
+    relTol_(relTol),
+    maxSteps_(10000)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::scalar Foam::ODESolver::normalizeError
+(
+    const scalarField& y0,
+    const scalarField& y,
+    const scalarField& err
+) const
+{
+    // Calculate the maximum error
+    scalar maxErr = 0.0;
+    forAll(err, i)
+    {
+        scalar tol = absTol_[i] + relTol_[i]*max(mag(y0[i]), mag(y[i]));
+        maxErr = max(maxErr, mag(err[i])/tol);
+    }
+
+    return maxErr;
+}
+
+
 void Foam::ODESolver::solve
 (
     const ODESystem& ode,
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
-    const scalar eps,
+    scalar& dxEst
-    scalar& hEst
 ) const
 {
-    const label MAXSTP = 10000;
-
     scalar x = xStart;
-    scalar h = hEst;
+    bool truncated = false;
-    scalar hNext = 0;
-    scalar hPrev = 0;
 
-    for (label nStep=0; nStep<MAXSTP; nStep++)
+    for (label nStep=0; nStep<maxSteps_; nStep++)
     {
-        ode.derivatives(x, y, dydx_);
+        // Store previous iteration dxEst
+        scalar dxEst0 = dxEst;
 
-        if ((x + h - xEnd)*(x + h - xStart) > 0.0)
+        // Check if this is a truncated step and set dxEst to integrate to xEnd
+        if ((x + dxEst - xEnd)*(x + dxEst - xStart) > 0)
         {
-            h = xEnd - x;
+            truncated = true;
-            hPrev = hNext;
+            dxEst = xEnd - x;
         }
 
-        hNext = 0;
+        // Integrate as far as possible up to dxEst
-        scalar hDid;
+        solve(ode, x, y, dxEst);
-        solve(ode, x, y, dydx_, eps, h, hDid, hNext);
 
-        if ((x - xEnd)*(xEnd - xStart) >= 0.0)
+        // Check if reached xEnd
+        if ((x - xEnd)*(xEnd - xStart) >= 0)
         {
-            if (hPrev != 0)
+            if (nStep > 0 && truncated)
             {
-                hEst = hPrev;
+                dxEst = dxEst0;
-            }
-            else
-            {
-                hEst = hNext;
             }
 
             return;
         }
-
-        h = hNext;
     }
 
     FatalErrorIn
     (
         "ODESolver::solve"
         "(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
-        "scalarField& yStart, const scalar eps, scalar& hEst) const"
+        "scalarField& y, scalar& dxEst) const"
-    )   << "Too many integration steps"
+    )   << "Integration steps greater than maximum " << maxSteps_
         << exit(FatalError);
 }
 
+
 // ************************************************************************* //

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -53,13 +53,30 @@ class ODESolver
 
 protected:
 
-    // Private data
+    // Protected data
 
+        //- Size of the ODESystem
         label n_;
-        mutable scalarField dydx_;
+
+        //- Absolute convergence tolerance per step
+        scalarField absTol_;
+
+        //- Relative convergence tolerance per step
+        scalarField relTol_;
+
+        //- The maximum number of sub-steps allowed for the integration step
+        label maxSteps_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
+
+        //- Return the nomalized scalar error
+        scalar normalizeError
+        (
+            const scalarField& y0,
+            const scalarField& y,
+            const scalarField& err
+        ) const;
 
         //- Disallow default bitwise copy construct
         ODESolver(const ODESolver&);
@@ -80,16 +97,24 @@ public:
         (
             autoPtr,
             ODESolver,
-            ODE,
+            dictionary,
-            (const ODESystem& ode),
+            (const ODESystem& ode, const dictionary& dict),
-            (ode)
+            (ode, dict)
         );
 
 
     // Constructors
 
-        //- Construct for given ODE
+        //- Construct for given ODESystem
-        ODESolver(const ODESystem& ode);
+        ODESolver(const ODESystem& ode, const dictionary& dict);
+
+        //- Construct for given ODESystem specifying tolerances
+        ODESolver
+        (
+            const ODESystem& ode,
+            const scalarField& absTol,
+            const scalarField& relTol
+        );
 
 
     // Selectors
@@ -97,8 +122,8 @@ public:
         //- Select null constructed
         static autoPtr<ODESolver> New
         (
-            const word& ODESolverTypeName,
+            const ODESystem& ode,
-            const ODESystem& ode
+            const dictionary& dict
         );
 
 
@@ -109,27 +134,37 @@ public:
 
     // Member Functions
 
+        scalarField& absTol()
+        {
+            return absTol_;
+        }
+
+        scalarField& relTol()
+        {
+            return relTol_;
+        }
+
+        //- Solve the ODE system as far as possible upto dxTry
+        //  adjusting the step as necessary to provide a solution within
+        //  the specified tolerance.
+        //  Update the state and return an estimate for the next step in dxTry
         virtual void solve
         (
             const ODESystem& ode,
             scalar& x,
             scalarField& y,
-            scalarField& dydx,
+            scalar& dxTry
-            const scalar eps,
-            const scalar hTry,
-            scalar& hDid,
-            scalar& hNext
         ) const = 0;
 
-
+        //- Solve the ODE system from xStart to xEnd, update the state
+        //  and return an estimate for the next step in dxTry
         virtual void solve
         (
             const ODESystem& ode,
             const scalar xStart,
             const scalar xEnd,
             scalarField& y,
-            const scalar eps,
+            scalar& dxEst
-            scalar& hEst
         ) const;
 };
 

src/ODE/ODESolvers/ODESolver/ODESolverNew.C

@@ -29,29 +29,30 @@ License
 
 Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
 (
-    const Foam::word& ODESolverTypeName,
+    const ODESystem& odes,
-    const Foam::ODESystem& ode
+    const dictionary& dict
 )
 {
+    word ODESolverTypeName(dict.lookup("solver"));
     Info<< "Selecting ODE solver " << ODESolverTypeName << endl;
 
-    ODEConstructorTable::iterator cstrIter =
+    dictionaryConstructorTable::iterator cstrIter =
-        ODEConstructorTablePtr_->find(ODESolverTypeName);
+        dictionaryConstructorTablePtr_->find(ODESolverTypeName);
 
-    if (cstrIter == ODEConstructorTablePtr_->end())
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
         FatalErrorIn
         (
             "ODESolver::New"
-            "(const word& ODESolverTypeName, const ODESystem& ode)"
+            "(const dictionary& dict, const ODESystem& odes)"
         )   << "Unknown ODESolver type "
             << ODESolverTypeName << nl << nl
             << "Valid ODESolvers are : " << endl
-            << ODEConstructorTablePtr_->sortedToc()
+            << dictionaryConstructorTablePtr_->sortedToc()
             << exit(FatalError);
     }
 
-    return autoPtr<ODESolver>(cstrIter()(ode));
+    return autoPtr<ODESolver>(cstrIter()(odes, dict));
 }
 
 

src/ODE/ODESolvers/RKCK45/RKCK45.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,41 +23,58 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "RKCK5.H"
+#include "RKCK45.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(RKCK5, 0);
+    defineTypeNameAndDebug(RKCK45, 0);
-    addToRunTimeSelectionTable(ODESolver, RKCK5, ODE);
+    addToRunTimeSelectionTable(ODESolver, RKCK45, dictionary);
 
 const scalar
-    RKCK5::a2 = 0.2, RKCK5::a3 = 0.3, RKCK5::a4 = 0.6, RKCK5::a5 = 1.0,
+    RKCK45::c2 = 1.0/5.0,
-    RKCK5::a6 = 0.875,
+    RKCK45::c3 = 3.0/10.0,
-    RKCK5::b21 = 0.2, RKCK5::b31 = 3.0/40.0, RKCK5::b32 = 9.0/40.0,
+    RKCK45::c4 = 3.0/5.0,
-    RKCK5::b41 = 0.3, RKCK5::b42 = -0.9, RKCK5::b43 = 1.2,
+    RKCK45::c5 = 1.0,
-    RKCK5::b51 = -11.0/54.0, RKCK5::b52 = 2.5, RKCK5::b53 = -70.0/27.0,
+    RKCK45::c6 = 7.0/8.0,
-    RKCK5::b54 = 35.0/27.0,
+
-    RKCK5::b61 = 1631.0/55296.0, RKCK5::b62 = 175.0/512.0,
+    RKCK45::a21 = 1.0/5.0,
-    RKCK5::b63 = 575.0/13824.0, RKCK5::b64 = 44275.0/110592.0,
+    RKCK45::a31 = 3.0/40.0,
-    RKCK5::b65 = 253.0/4096.0,
+    RKCK45::a32 = 9.0/40.0,
-    RKCK5::c1 = 37.0/378.0, RKCK5::c3 = 250.0/621.0,
+    RKCK45::a41 = 3.0/10.0,
-    RKCK5::c4 = 125.0/594.0, RKCK5::c6 = 512.0/1771.0,
+    RKCK45::a42 = -9.0/10.0,
-    RKCK5::dc1 = RKCK5::c1 - 2825.0/27648.0,
+    RKCK45::a43 = 6.0/5.0,
-    RKCK5::dc3 = RKCK5::c3 - 18575.0/48384.0,
+    RKCK45::a51 = -11.0/54.0,
-    RKCK5::dc4 = RKCK5::c4 - 13525.0/55296.0, RKCK5::dc5 = -277.00/14336.0,
+    RKCK45::a52 = 5.0/2.0,
-    RKCK5::dc6 = RKCK5::c6 - 0.25;
+    RKCK45::a53 = -70.0/27.0,
+    RKCK45::a54 = 35.0/27.0,
+    RKCK45::a61 = 1631.0/55296.0,
+    RKCK45::a62 = 175.0/512.0,
+    RKCK45::a63 = 575.0/13824.0,
+    RKCK45::a64 = 44275.0/110592.0,
+    RKCK45::a65 = 253.0/4096.0,
+
+    RKCK45::b1 = 37.0/378.0,
+    RKCK45::b3 = 250.0/621.0,
+    RKCK45::b4 = 125.0/594.0,
+    RKCK45::b6 = 512.0/1771.0,
+
+    RKCK45::e1 = RKCK45::b1 - 2825.0/27648.0,
+    RKCK45::e3 = RKCK45::b3 - 18575.0/48384.0,
+    RKCK45::e4 = RKCK45::b4 - 13525.0/55296.0,
+    RKCK45::e5 = -277.00/14336.0,
+    RKCK45::e6 = RKCK45::b6 - 0.25;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::RKCK5::RKCK5(const ODESystem& ode)
+Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict)
 :
-    ODESolver(ode),
+    ODESolver(ode, dict),
-    adaptiveSolver(ode),
+    adaptiveSolver(ode, dict),
     yTemp_(n_),
     k2_(n_),
     k3_(n_),
@@ -70,85 +87,80 @@ Foam::RKCK5::RKCK5(const ODESystem& ode)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::RKCK5::solve
+Foam::scalar Foam::RKCK45::solve
 (
     const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
     const scalar dx,
-    scalarField& y,
+    scalarField& y
-    const scalar eps
 ) const
 {
     forAll(yTemp_, i)
     {
-        yTemp_[i] = y0[i] + b21*dx*dydx0[i];
+        yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + a2*dx, yTemp_, k2_);
+    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
 
     forAll(yTemp_, i)
     {
-        yTemp_[i] = y0[i] + dx*(b31*dydx0[i] + b32*k2_[i]);
+        yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes.derivatives(x0 + a3*dx, yTemp_, k3_);
+    odes.derivatives(x0 + c3*dx, yTemp_, k3_);
 
     forAll(yTemp_, i)
     {
-        yTemp_[i] = y0[i] + dx*(b41*dydx0[i] + b42*k2_[i] + b43*k3_[i]);
+        yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes.derivatives(x0 + a4*dx, yTemp_, k4_);
+    odes.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i]
-          + dx*(b51*dydx0[i] + b52*k2_[i] + b53*k3_[i] + b54*k4_[i]);
+          + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + a5*dx, yTemp_, k5_);
+    odes.derivatives(x0 + c5*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]);
+           *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + a6*dx, yTemp_, k6_);
+    odes.derivatives(x0 + c6*dx, yTemp_, k6_);
 
     forAll(y, i)
     {
         y[i] = y0[i]
-          + dx*(c1*dydx0[i] + c3*k3_[i] + c4*k4_[i] + c6*k6_[i]);
+          + dx*(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b6*k6_[i]);
     }
 
     forAll(err_, i)
     {
         err_[i] =
             dx
-           *(dc1*dydx0[i] + dc3*k3_[i] + dc4*k4_[i] + dc5*k5_[i] + dc6*k6_[i]);
+           *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[i]);
     }
 
-    return normalizeError(eps, y0, y, err_);
+    return normalizeError(y0, y, err_);
 }
 
 
-void Foam::RKCK5::solve
+void Foam::RKCK45::solve
 (
     const ODESystem& odes,
     scalar& x,
     scalarField& y,
-    scalarField& dydx,
+    scalar& dxTry
-    const scalar eps,
-    const scalar dxTry,
-    scalar& dxDid,
-    scalar& dxNext
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dydx, eps, dxTry, dxDid, dxNext);
+    adaptiveSolver::solve(odes, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/RKCK45/RKCK45.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,10 +22,10 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::RKCK5
+    Foam::RKCK45
 
 Description
-    5th Order Cash-Karp Runge-Kutta ODE solver
+    4/5th Order Cash-Karp Runge-Kutta ODE solver
 
     References:
     \verbatim
@@ -37,12 +37,12 @@ Description
     \endverbatim
 
 SourceFiles
-    RKCK5.C
+    RKCK45.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef RKCK5_H
+#ifndef RKCK45_H
-#define RKCK5_H
+#define RKCK45_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -53,10 +53,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class RKCK5 Declaration
+                           Class RKCK45 Declaration
 \*---------------------------------------------------------------------------*/
 
-class RKCK5
+class RKCK45
 :
     public ODESolver,
     public adaptiveSolver
@@ -65,11 +65,11 @@ class RKCK5
 
         //- RKCK Constants
         static const scalar
-            a2, a3, a4, a5, a6,
+            c2, c3, c4, c5, c6,
-            b21, b31, b32, b41, b42, b43,
+            a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
-            b51, b52, b53, b54, b61, b62, b63, b64, b65,
+            a61, a62, a63, a64, a65,
-            c1, c3, c4, c6,
+            b1, b3, b4, b6,
-            dc1, dc3, dc4, dc5, dc6;
+            e1, e3, e4, e5, e6;
 
         // Temporary fields
         mutable scalarField yTemp_;
@@ -86,13 +86,13 @@ class RKCK5
 public:
 
     //- Runtime type information
-    TypeName("RKCK5");
+    TypeName("RKCK45");
 
 
     // Constructors
 
         //- Construct from ODE
-        RKCK5(const ODESystem& ode);
+        RKCK45(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -105,8 +105,7 @@ public:
             const scalarField& y0,
             const scalarField& dydx0,
             const scalar dx,
-            scalarField& y,
+            scalarField& y
-            const scalar eps
         ) const;
 
         //- Solve the ODE system and the update the state
@@ -115,11 +114,7 @@ public:
             const ODESystem& ode,
             scalar& x,
             scalarField& y,
-            scalarField& dydx,
+            scalar& dxTry
-            const scalar eps,
-            const scalar dxTry,
-            scalar& dxDid,
-            scalar& dxNext
         ) const;
 };
 

src/ODE/ODESolvers/RKF45/RKF45.C

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "RKF45.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(RKF45, 0);
+    addToRunTimeSelectionTable(ODESolver, RKF45, dictionary);
+
+const scalar
+    RKF45::c2  = 1.0/4.0,
+    RKF45::c3  = 3.0/8.0,
+    RKF45::c4  = 12.0/13.0,
+    RKF45::c5  = 1.0,
+    RKF45::c6  = 1.0/2.0,
+
+    RKF45::a21 = 1.0/4.0,
+    RKF45::a31 = 3.0/32.0,
+    RKF45::a32 = 9.0/32.0,
+    RKF45::a41 = 1932.0/2197.0,
+    RKF45::a42 = -7200.0/2197.0,
+    RKF45::a43 = 7296.0/2197.0,
+    RKF45::a51 = 439.0/216.0,
+    RKF45::a52 = -8.0,
+    RKF45::a53 = 3680.0/513.0,
+    RKF45::a54 = -845.0/4104.0,
+    RKF45::a61 = -8.0/27.0,
+    RKF45::a62 = 2.0,
+    RKF45::a63 = -3544.0/2565.0,
+    RKF45::a64 = 1859.0/4104.0,
+    RKF45::a65 = -11.0/40.0,
+
+    RKF45::b41  = 25.0/216.0,
+    RKF45::b42  = 0.0,
+    RKF45::b43  = 1408.0/2565.0,
+    RKF45::b44  = 2197.0/4104.0,
+    RKF45::b45  = -1.0/5.0,
+    RKF45::b46  = 0.0,
+
+    RKF45::b51  = 16.0/135.0,
+    RKF45::b52  = 0.0,
+    RKF45::b53  = 6656.0/12825.0,
+    RKF45::b54  = 28561.0/56430.0,
+    RKF45::b55  = -9.0/50.0,
+    RKF45::b56  = 2.0/55.0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::RKF45::RKF45(const ODESystem& ode, const dictionary& dict)
+:
+    ODESolver(ode, dict),
+    adaptiveSolver(ode, dict),
+    yTemp_(n_),
+    k2_(n_),
+    k3_(n_),
+    k4_(n_),
+    k5_(n_),
+    k6_(n_),
+    err_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::RKF45::solve
+(
+    const ODESystem& odes,
+    const scalar x0,
+    const scalarField& y0,
+    const scalarField& dydx0,
+    const scalar dx,
+    scalarField& y
+) const
+{
+    forAll(yTemp_, i)
+    {
+        yTemp_[i] = y0[i] + a21*dx*dydx0[i];
+    }
+
+    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+
+    forAll(yTemp_, i)
+    {
+        yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
+    }
+
+    odes.derivatives(x0 + c3*dx, yTemp_, k3_);
+
+    forAll(yTemp_, i)
+    {
+        yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
+    }
+
+    odes.derivatives(x0 + c4*dx, yTemp_, k4_);
+
+    forAll(yTemp_, i)
+    {
+        yTemp_[i] = y0[i]
+          + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
+    }
+
+    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+
+    forAll(yTemp_, i)
+    {
+        yTemp_[i] = y0[i]
+          + dx
+           *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
+    }
+
+    odes.derivatives(x0 + c6*dx, yTemp_, k6_);
+
+    // Calculate the 5th-order solution
+    forAll(y, i)
+    {
+        y[i] = y0[i]
+          + dx
+           *(b51*dydx0[i] + b53*k3_[i] + b54*k4_[i] + b55*k5_[i] + b56*k6_[i]);
+    }
+
+    // Calculate the error estimate from the difference between the
+    // 4th-order and 5th-order solutions
+    forAll(err_, i)
+    {
+        err_[i] =
+            dx
+           *(b41*dydx0[i] + b43*k3_[i] + b44*k4_[i] + b45*k5_[i])
+          - (y[i] - y0[i])
+        ;
+    }
+
+    return normalizeError(y0, y, err_);
+}
+
+
+void Foam::RKF45::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/RKF45/RKF45.H

@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::RKF45
+
+Description
+    4/5th Order Runge-Kutta-Fehlberg ODE solver
+
+    References:
+    \verbatim
+        "Low-order classical Runge-Kutta formulas with step size control
+         and their application to some heat transfer problems."
+         Fehlberg, E.,
+         NASA Technical Report 315 1969
+
+         "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
+
+    This method embedds the 4-th order integration step into the 5-th order step
+    and allows to perform an adapdive step-size control using these two order
+    without the need of re-evaluation.
+
+    \f{align}{
+     k_1 &= h f(t_k, y_k) \\
+     k_2 &= h f(t_k + \frac{1}{4}h, y_k + \frac{1}{4}k_1) \\
+     k_3 &= h f(t_k + \frac{3}{8}h, y_k + \frac{3}{32}k_1 + \frac{9}{32}k_2) \\
+     k_4 &= h f(t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}k_1
+              \frac{7200}{2197}k_2 + \frac{7296}{2197}k_3) - \\
+     k_5 &= h f(t_k + h, y_k + \frac{439}{216}k_1 - 8k_2 + \frac{3680}{513}k_3 -
+           \frac{845}{4104}k_4) \\
+     k_6 &= h f(t_k + \frac{1}{2}h, y_k - \frac{8}{27}k_1 + 2k_2 -
+           \frac{3544}{2565}k_3 + \frac{1859}{4104}k_4 - \frac{11}{40}k_5) \\
+    \f}
+
+    Which yields the following update-steps for the 4th and 5th order:
+
+    \f{align}{
+     \Delta_4 &= \frac{25}{216}k_1 + \frac{1408}{2565}k_3 +
+                 \frac{2197}{4101}k_4 - \frac{1}{5}k_5 \\
+     \Delta_5 &= \frac{16}{135}k_1 + \frac{6656}{12825}k_3 +
+                 \frac{9}{50}k_5 + \frac{2}{55}k_6
+    \f}
+
+    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$.
+
+SourceFiles
+    RKF45.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RKF45_H
+#define RKF45_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class RKF45 Declaration
+\*---------------------------------------------------------------------------*/
+
+class RKF45
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        //- RKF45 Constants
+        static const scalar
+            c2, c3, c4, c5, c6,
+            a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
+            a61, a62, a63, a64, a65,
+            b41, b42, b43, b44, b45, b46,
+            b51, b52, b53, b54, b55, b56;
+
+
+        // Temporary fields
+        mutable scalarField yTemp_;
+        mutable scalarField k2_;
+        mutable scalarField k3_;
+        mutable scalarField k4_;
+        mutable scalarField k5_;
+        mutable scalarField k6_;
+
+        //- Error-estimate field
+        mutable scalarField err_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("RKF45");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        RKF45(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
+
+// ************************************************************************* //

src/ODE/ODESolvers/SIBS/SIBS.C

@@ -31,8 +31,7 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(SIBS, 0);
-
+    addToRunTimeSelectionTable(ODESolver, SIBS, dictionary);
-    addToRunTimeSelectionTable(ODESolver, SIBS, ODE);
 
     const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70};
 
@@ -47,9 +46,9 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::SIBS::SIBS(const ODESystem& ode)
+Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
 :
-    ODESolver(ode),
+    ODESolver(ode, dict),
     a_(iMaxX_, 0.0),
     alpha_(kMaxX_, kMaxX_, 0.0),
     d_p_(n_, kMaxX_, 0.0),
@@ -59,6 +58,7 @@ Foam::SIBS::SIBS(const ODESystem& ode)
     yTemp_(n_, 0.0),
     ySeq_(n_, 0.0),
     yErr_(n_, 0.0),
+    dydx0_(n_),
     dfdx_(n_, 0.0),
     dfdy_(n_, n_, 0.0),
     first_(1),
@@ -73,19 +73,18 @@ void Foam::SIBS::solve
     const ODESystem& ode,
     scalar& x,
     scalarField& y,
-    scalarField& dydx,
+    scalar& dxTry
-    const scalar eps,
-    const scalar hTry,
-    scalar& hDid,
-    scalar& hNext
 ) const
 {
+    ode.derivatives(x, y, dydx0_);
+
+    scalar h = dxTry;
     bool exitflag = false;
 
-    if (eps != epsOld_)
+    if (relTol_[0] != epsOld_)
     {
-        hNext = xNew_ = -GREAT;
+        dxTry = xNew_ = -GREAT;
-        scalar eps1 = safe1*eps;
+        scalar eps1 = safe1*relTol_[0];
         a_[0] = nSeq_[0] + 1;
 
         for (register label k=0; k<kMaxX_; k++)
@@ -103,7 +102,7 @@ void Foam::SIBS::solve
             }
         }
 
-        epsOld_ = eps;
+        epsOld_ = relTol_[0];
         a_[0] += n_;
 
         for (register label k=0; k<kMaxX_; k++)
@@ -123,12 +122,11 @@ void Foam::SIBS::solve
     }
 
     label k = 0;
-    scalar h = hTry;
     yTemp_ = y;
 
     ode.jacobian(x, y, dfdx_, dfdy_);
 
-    if (x != xNew_ || h != hNext)
+    if (x != xNew_ || h != dxTry)
     {
         first_ = 1;
         kOpt_ = kMax_;
@@ -153,7 +151,7 @@ void Foam::SIBS::solve
                     << exit(FatalError);
             }
 
-            SIMPR(ode, x, yTemp_, dydx, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
+            SIMPR(ode, x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
             scalar xest = sqr(h/nSeq_[k]);
 
             polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);
@@ -166,10 +164,9 @@ void Foam::SIBS::solve
                     maxErr = max
                     (
                         maxErr,
-                        mag(yErr_[i])/(mag(yTemp_[i]) + SMALL)
+                        mag(yErr_[i])/(absTol_[i] + relTol_[i]*mag(yTemp_[i]))
                     );
                 }
-                maxErr /= eps;
                 km = k - 1;
                 err_[km] = pow(maxErr/safe1, 1.0/(2*km + 3));
             }
@@ -217,7 +214,6 @@ void Foam::SIBS::solve
     }
 
     x = xNew_;
-    hDid = h;
     first_ = 0;
     scalar wrkmin = GREAT;
     scalar scale = 1.0;
@@ -234,14 +230,14 @@ void Foam::SIBS::solve
         }
     }
 
-    hNext = h/scale;
+    dxTry = h/scale;
 
     if (kOpt_ >= k && kOpt_ != kMax_ && !reduct)
     {
         scalar fact = max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX);
         if (a_[kOpt_ + 1]*fact <= wrkmin)
         {
-            hNext = h/fact;
+            dxTry = h/fact;
             kOpt_++;
         }
     }

src/ODE/ODESolvers/SIBS/SIBS.H

@@ -72,6 +72,7 @@ class SIBS
         mutable scalarField yTemp_;
         mutable scalarField ySeq_;
         mutable scalarField yErr_;
+        mutable scalarField dydx0_;
         mutable scalarField dfdx_;
         mutable scalarSquareMatrix dfdy_;
 
@@ -115,7 +116,7 @@ public:
     // Constructors
 
         //- Construct from ODE
-        SIBS(const ODESystem& ode);
+        SIBS(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -125,11 +126,7 @@ public:
             const ODESystem& ode,
             scalar& x,
             scalarField& y,
-            scalarField& dydx,
+            scalar& dxTry
-            const scalar eps,
-            const scalar hTry,
-            scalar& hDid,
-            scalar& hNext
         ) const;
 };
 

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C

@@ -26,79 +26,50 @@ License
 #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)
+Foam::adaptiveSolver::adaptiveSolver
+(
+    const ODESystem& ode,
+    const dictionary& dict
+)
 :
+    safeScale_(dict.lookupOrDefault<scalar>("safeScale", 0.9)),
+    alphaInc_(dict.lookupOrDefault<scalar>("alphaIncrease", 0.2)),
+    alphaDec_(dict.lookupOrDefault<scalar>("alphaDecrease", 0.25)),
+    minScale_(dict.lookupOrDefault<scalar>("minScale", 0.2)),
+    maxScale_(dict.lookupOrDefault<scalar>("maxScale", 10)),
+    dydx0_(ode.nEqns()),
     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,
+    scalar& dxTry
-    const scalar eps,
-    const scalar dxTry,
-    scalar& dxDid,
-    scalar& dxNext
 ) const
 {
     scalar dx = dxTry;
     scalar err = 0.0;
 
+    odes.derivatives(x, y, dydx0_);
+
     // 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);
+        err = solve(odes, x, y, dydx0_, dx, yTemp_);
 
         // If error is large reduce dx
         if (err > 1)
         {
-            scalar scale = max(safeScale*pow(err, -alphaDec), minScale);
+            scalar scale = max(safeScale_*pow(err, -alphaDec_), minScale_);
             dx *= scale;
 
             if (dx < VSMALL)
@@ -111,12 +82,11 @@ void Foam::adaptiveSolver::solve
     } 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;
+    dxTry = min(max(safeScale_*pow(err, -alphaInc_), minScale_), maxScale_)*dx;
 }
 
 

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,16 +25,6 @@ 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
@@ -60,8 +50,12 @@ class adaptiveSolver
     // Private data
 
         //- Step-size adjustment controls
-        static const scalar safeScale, alphaInc, alphaDec, minScale, maxScale;
+        scalar safeScale_, alphaInc_, alphaDec_, minScale_, maxScale_;
 
+        //- Cache for dydx at the initial time
+        mutable scalarField dydx0_;
+
+        //- Temprorary for the test-step solution
         mutable scalarField yTemp_;
 
 
@@ -70,7 +64,7 @@ public:
     // Constructors
 
         //- Construct from ODESystem
-        adaptiveSolver(const ODESystem& ode);
+        adaptiveSolver(const ODESystem& ode, const dictionary& dict);
 
 
     //- Destructor
@@ -80,15 +74,6 @@ public:
 
     // 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
         (
@@ -97,8 +82,7 @@ public:
             const scalarField& y0,
             const scalarField& dydx0,
             const scalar dx,
-            scalarField& y,
+            scalarField& y
-            const scalar eps
         ) const = 0;
 
         //- Solve the ODE system and the update the state
@@ -107,11 +91,7 @@ public:
             const ODESystem& ode,
             scalar& x,
             scalarField& y,
-            scalarField& dydx,
+            scalar& dxTry
-            const scalar eps,
-            const scalar dxTry,
-            scalar& dxDid,
-            scalar& dxNext
         ) const;
 };
 

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C

@@ -36,9 +36,7 @@ Foam::ode<ChemistryModel>::ode
 :
     chemistrySolver<ChemistryModel>(mesh),
     coeffsDict_(this->subDict("odeCoeffs")),
-    solverName_(coeffsDict_.lookup("solver")),
+    odeSolver_(ODESolver::New(*this, coeffsDict_)),
-    odeSolver_(ODESolver::New(solverName_, *this)),
-    eps_(readScalar(coeffsDict_.lookup("eps"))),
     cTp_(this->nEqns())
 {}
 
@@ -78,7 +76,6 @@ void Foam::ode<ChemistryModel>::solve
         0,
         deltaT,
         cTp_,
-        eps_,
         subDeltaT
     );
 

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H

@@ -55,13 +55,8 @@ class ode
     // Private data
 
         dictionary coeffsDict_;
-        const word solverName_;
         autoPtr<ODESolver> odeSolver_;
 
-        // Model constants
-
-            scalar eps_;
-
         // Solver data
         mutable scalarField cTp_;