Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/test/ODE/Test-ODE.C

@@ -142,7 +142,7 @@ int main(int argc, char *argv[])
         scalar dxNext = dxEst;
 
         odeSolver->relTol() = relTol;
-        odeSolver->solve(ode, x, y, dxNext);
+        odeSolver->solve(x, y, dxNext);
 
         Info<< scientific << setw(13) << relTol;
         Info<< fixed << setw(11) << dxEst;
@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
     scalar dxEst = 0.5;
 
     odeSolver->relTol() = 1e-4;
-    odeSolver->solve(ode, x, xEnd, y, dxEst);
+    odeSolver->solve(x, xEnd, y, dxEst);
 
     Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
     Info      << "Numerical:  y(2.0) = " << y << ", dxEst = " << dxEst << endl;

src/ODE/Make/files

@@ -8,11 +8,11 @@ ODESolvers/Trapezoid/Trapezoid.C
 ODESolvers/RKF45/RKF45.C
 ODESolvers/RKCK45/RKCK45.C
 ODESolvers/RKDP45/RKDP45.C
-ODESolvers/Rosenbrock21/Rosenbrock21.C
+ODESolvers/Rosenbrock12/Rosenbrock12.C
-ODESolvers/Rosenbrock32/Rosenbrock32.C
+ODESolvers/Rosenbrock23/Rosenbrock23.C
-ODESolvers/Rosenbrock43/Rosenbrock43.C
+ODESolvers/Rosenbrock34/Rosenbrock34.C
-ODESolvers/rodas32/rodas32.C
+ODESolvers/rodas23/rodas23.C
-ODESolvers/rodas43/rodas43.C
+ODESolvers/rodas34/rodas34.C
 ODESolvers/SIBS/SIBS.C
 ODESolvers/SIBS/SIMPR.C
 ODESolvers/SIBS/polyExtrapolate.C

src/ODE/ODESolvers/Euler/Euler.C

@@ -49,7 +49,6 @@ Foam::Euler::Euler(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Euler::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -75,13 +74,12 @@ Foam::scalar Foam::Euler::solve
 
 void Foam::Euler::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Euler/Euler.H

@@ -84,7 +84,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -95,7 +94,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/EulerSI/EulerSI.C

@@ -54,7 +54,6 @@ Foam::EulerSI::EulerSI(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::EulerSI::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -62,7 +61,7 @@ Foam::scalar Foam::EulerSI::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -95,13 +94,12 @@ Foam::scalar Foam::EulerSI::solve
 
 void Foam::EulerSI::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/EulerSI/EulerSI.H

@@ -91,7 +91,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -102,7 +101,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/ODESolver/ODESolver.C

@@ -38,6 +38,7 @@ namespace Foam
 
 Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict)
 :
+    odes_(ode),
     n_(ode.nEqns()),
     absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
     relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-4)),
@@ -52,6 +53,7 @@ Foam::ODESolver::ODESolver
     const scalarField& relTol
 )
 :
+    odes_(ode),
     n_(ode.nEqns()),
     absTol_(absTol),
     relTol_(relTol),
@@ -82,48 +84,70 @@ Foam::scalar Foam::ODESolver::normalizeError
 
 void Foam::ODESolver::solve
 (
-    const ODESystem& ode,
+    scalar& x,
+    scalarField& y,
+    stepState& step
+) const
+{
+    solve(x, y, step.dxTry);
+}
+
+
+void Foam::ODESolver::solve
+(
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
-    scalar& dxEst
+    scalar& dxTry
 ) const
 {
+    stepState step(dxTry);
     scalar x = xStart;
-    bool truncated = false;
 
     for (label nStep=0; nStep<maxSteps_; nStep++)
     {
-        // Store previous iteration dxEst
+        // Store previous iteration dxTry
-        scalar dxEst0 = dxEst;
+        scalar dxTry0 = step.dxTry;
 
-        // Check if this is a truncated step and set dxEst to integrate to xEnd
+        step.reject = false;
-        if ((x + dxEst - xEnd)*(x + dxEst - xStart) > 0)
+
+        // Check if this is a truncated step and set dxTry to integrate to xEnd
+        if ((x + step.dxTry - xEnd)*(x + step.dxTry - xStart) > 0)
         {
-            truncated = true;
+            step.last = true;
-            dxEst = xEnd - x;
+            step.dxTry = xEnd - x;
         }
 
-        // Integrate as far as possible up to dxEst
+        // Integrate as far as possible up to step.dxTry
-        solve(ode, x, y, dxEst);
+        solve(x, y, step);
 
         // Check if reached xEnd
         if ((x - xEnd)*(xEnd - xStart) >= 0)
         {
-            if (nStep > 0 && truncated)
+            if (nStep > 0 && step.last)
             {
-                dxEst = dxEst0;
+                step.dxTry = dxTry0;
             }
 
+            dxTry = step.dxTry;
+
             return;
         }
+
+        step.first = false;
+
+        // If the step.dxTry was reject set step.prevReject
+        if (step.reject)
+        {
+            step.prevReject = true;
+        }
     }
 
     FatalErrorIn
     (
         "ODESolver::solve"
-        "(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
+        "(const scalar xStart, const scalar xEnd,"
-        "scalarField& y, scalar& dxEst) const"
+        "scalarField& y, scalar& dxTry) const"
     )   << "Integration steps greater than maximum " << maxSteps_
         << exit(FatalError);
 }

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -55,6 +55,9 @@ protected:
 
     // Protected data
 
+        //- Reference to ODESystem
+        const ODESystem& odes_;
+
         //- Size of the ODESystem
         label n_;
 
@@ -90,6 +93,30 @@ public:
     //- Runtime type information
     TypeName("ODESolver");
 
+    class stepState
+    {
+        public:
+
+        const bool forward;
+        scalar dxTry;
+        scalar dxDid;
+        bool first;
+        bool last;
+        bool reject;
+        bool prevReject;
+
+        stepState(const scalar dx)
+        :
+            forward(dx > 0 ? true : false),
+            dxTry(dx),
+            dxDid(0),
+            first(true),
+            last(false),
+            reject(false),
+            prevReject(false)
+        {}
+    };
+
 
     // Declare run-time constructor selection table
 
@@ -150,17 +177,31 @@ public:
         //  Update the state and return an estimate for the next step in dxTry
         virtual void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry
-        ) const = 0;
+        ) const //= 0;
+        {
+            stepState step(dxTry);
+            solve(x, y, step);
+            dxTry = step.dxTry;
+        }
+
+        //- 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
+        (
+            scalar& x,
+            scalarField& y,
+            stepState& step
+        ) const;
 
         //- 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,

src/ODE/ODESolvers/RKCK45/RKCK45.C

@@ -89,7 +89,6 @@ Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKCK45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -102,21 +101,21 @@ Foam::scalar Foam::RKCK45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    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_);
+    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_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -124,7 +123,7 @@ Foam::scalar Foam::RKCK45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -133,7 +132,7 @@ Foam::scalar Foam::RKCK45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + c6*dx, yTemp_, k6_);
+    odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
 
     forAll(y, i)
     {
@@ -154,13 +153,12 @@ Foam::scalar Foam::RKCK45::solve
 
 void Foam::RKCK45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/RKCK45/RKCK45.H

@@ -34,10 +34,7 @@ Description
         Cash, J.R.,
         Karp, A.H.
         ACM Transactions on Mathematical Software, vol. 16, 1990, pp. 201–222.
-    \endverbatim
 
-    Based on code from:
-    \verbatim
         "Solving Ordinary Differential Equations I: Nonstiff Problems,
          second edition",
         Hairer, E.,
@@ -110,7 +107,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -121,7 +117,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/RKDP45/RKDP45.C

@@ -94,7 +94,6 @@ Foam::RKDP45::RKDP45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKDP45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -107,21 +106,21 @@ Foam::scalar Foam::RKDP45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    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_);
+    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_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -129,7 +128,7 @@ Foam::scalar Foam::RKDP45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -138,7 +137,7 @@ Foam::scalar Foam::RKDP45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + dx, yTemp_, k6_);
+    odes_.derivatives(x0 + dx, yTemp_, k6_);
 
     forAll(y, i)
     {
@@ -147,7 +146,7 @@ Foam::scalar Foam::RKDP45::solve
     }
 
     // Reuse k2_ for the derivative of the new state
-    odes.derivatives(x0 + dx, y, k2_);
+    odes_.derivatives(x0 + dx, y, k2_);
 
     forAll(err_, i)
     {
@@ -163,13 +162,12 @@ Foam::scalar Foam::RKDP45::solve
 
 void Foam::RKDP45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/RKDP45/RKDP45.H

@@ -34,10 +34,7 @@ Description
         Prince, P. J.,
         Journal of Computational and Applied Mathematics,
         6 (1), 1980: pp. 19-26.
-    \endverbatim
 
-    Based on code from:
-    \verbatim
         "Solving Ordinary Differential Equations I: Nonstiff Problems,
          second edition",
         Hairer, E.,
@@ -114,7 +111,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -125,7 +121,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/RKF45/RKF45.C

@@ -90,7 +90,6 @@ Foam::RKF45::RKF45(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::RKF45::solve
 (
-    const ODESystem& odes,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -103,21 +102,21 @@ Foam::scalar Foam::RKF45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes.derivatives(x0 + c2*dx, yTemp_, k2_);
+    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_);
+    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_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
 
     forAll(yTemp_, i)
     {
@@ -125,7 +124,7 @@ Foam::scalar Foam::RKF45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
 
     forAll(yTemp_, i)
     {
@@ -134,7 +133,7 @@ Foam::scalar Foam::RKF45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes.derivatives(x0 + c6*dx, yTemp_, k6_);
+    odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
 
     // Calculate the 5th-order solution
     forAll(y, i)
@@ -159,13 +158,12 @@ Foam::scalar Foam::RKF45::solve
 
 void Foam::RKF45::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/RKF45/RKF45.H

@@ -34,8 +34,6 @@ Description
         Fehlberg, E.,
         NASA Technical Report 315, 1969.
 
-    Based on code from:
-    \verbatim
         "Solving Ordinary Differential Equations I: Nonstiff Problems,
          second edition",
         Hairer, E.,
@@ -48,31 +46,6 @@ Description
     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.
-
 SourceFiles
     RKF45.C
 
@@ -138,7 +111,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -149,7 +121,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C

@@ -23,33 +23,33 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Rosenbrock21.H"
+#include "Rosenbrock12.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(Rosenbrock21, 0);
+    defineTypeNameAndDebug(Rosenbrock12, 0);
-    addToRunTimeSelectionTable(ODESolver, Rosenbrock21, dictionary);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock12, dictionary);
 
 const scalar
-    Rosenbrock21::gamma = 1 + 1.0/std::sqrt(2.0),
+    Rosenbrock12::gamma = 1 + 1.0/std::sqrt(2.0),
-    Rosenbrock21::a21 = 1.0/gamma,
+    Rosenbrock12::a21 = 1.0/gamma,
-    Rosenbrock21::c2 = 1.0,
+    Rosenbrock12::c2 = 1.0,
-    Rosenbrock21::c21 = -2.0/gamma,
+    Rosenbrock12::c21 = -2.0/gamma,
-    Rosenbrock21::b1 = (3.0/2.0)/gamma,
+    Rosenbrock12::b1 = (3.0/2.0)/gamma,
-    Rosenbrock21::b2 = (1.0/2.0)/gamma,
+    Rosenbrock12::b2 = (1.0/2.0)/gamma,
-    Rosenbrock21::e1 = b1 - 1.0/gamma,
+    Rosenbrock12::e1 = b1 - 1.0/gamma,
-    Rosenbrock21::e2 = b2,
+    Rosenbrock12::e2 = b2,
-    Rosenbrock21::d1 = gamma,
+    Rosenbrock12::d1 = gamma,
-    Rosenbrock21::d2 = -gamma;
+    Rosenbrock12::d2 = -gamma;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict)
+Foam::Rosenbrock12::Rosenbrock12(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -66,9 +66,8 @@ Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::Rosenbrock21::solve
+Foam::scalar Foam::Rosenbrock12::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -76,7 +75,7 @@ Foam::scalar Foam::Rosenbrock21::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -104,7 +103,7 @@ Foam::scalar Foam::Rosenbrock21::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -124,15 +123,14 @@ Foam::scalar Foam::Rosenbrock21::solve
 }
 
 
-void Foam::Rosenbrock21::solve
+void Foam::Rosenbrock12::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::Rosenbrock21
+    Foam::Rosenbrock12
 
 Description
     L-stable embedded Rosenbrock ODE solver of order (1)2.
@@ -39,12 +39,12 @@ Description
     \endverbatim
 
 SourceFiles
-    Rosenbrock21.C
+    Rosenbrock12.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Rosenbrock21_H
+#ifndef Rosenbrock12_H
-#define Rosenbrock21_H
+#define Rosenbrock12_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -55,10 +55,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class Rosenbrock21 Declaration
+                           Class Rosenbrock12 Declaration
 \*---------------------------------------------------------------------------*/
 
-class Rosenbrock21
+class Rosenbrock12
 :
     public ODESolver,
     public adaptiveSolver
@@ -87,13 +87,13 @@ class Rosenbrock21
 public:
 
     //- Runtime type information
-    TypeName("Rosenbrock21");
+    TypeName("Rosenbrock12");
 
 
     // Constructors
 
         //- Construct from ODE
-        Rosenbrock21(const ODESystem& ode, const dictionary& dict);
+        Rosenbrock12(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -101,7 +101,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -112,7 +111,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C

@@ -23,45 +23,45 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Rosenbrock32.H"
+#include "Rosenbrock23.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(Rosenbrock32, 0);
+    defineTypeNameAndDebug(Rosenbrock23, 0);
-    addToRunTimeSelectionTable(ODESolver, Rosenbrock32, dictionary);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock23, dictionary);
 
 const scalar
-    Rosenbrock32::a21 = 1,
+    Rosenbrock23::a21 = 1,
-    Rosenbrock32::a31 = 1,
+    Rosenbrock23::a31 = 1,
-    Rosenbrock32::a32 = 0,
+    Rosenbrock23::a32 = 0,
 
-    Rosenbrock32::c21 = -1.0156171083877702091975600115545,
+    Rosenbrock23::c21 = -1.0156171083877702091975600115545,
-    Rosenbrock32::c31 = 4.0759956452537699824805835358067,
+    Rosenbrock23::c31 = 4.0759956452537699824805835358067,
-    Rosenbrock32::c32 = 9.2076794298330791242156818474003,
+    Rosenbrock23::c32 = 9.2076794298330791242156818474003,
 
-    Rosenbrock32::b1 = 1,
+    Rosenbrock23::b1 = 1,
-    Rosenbrock32::b2 = 6.1697947043828245592553615689730,
+    Rosenbrock23::b2 = 6.1697947043828245592553615689730,
-    Rosenbrock32::b3 = -0.4277225654321857332623837380651,
+    Rosenbrock23::b3 = -0.4277225654321857332623837380651,
 
-    Rosenbrock32::e1 = 0.5,
+    Rosenbrock23::e1 = 0.5,
-    Rosenbrock32::e2 = -2.9079558716805469821718236208017,
+    Rosenbrock23::e2 = -2.9079558716805469821718236208017,
-    Rosenbrock32::e3 = 0.2235406989781156962736090927619,
+    Rosenbrock23::e3 = 0.2235406989781156962736090927619,
 
-    Rosenbrock32::gamma = 0.43586652150845899941601945119356,
+    Rosenbrock23::gamma = 0.43586652150845899941601945119356,
-    Rosenbrock32::c2 = 0.43586652150845899941601945119356,
+    Rosenbrock23::c2 = 0.43586652150845899941601945119356,
 
-    Rosenbrock32::d1 = 0.43586652150845899941601945119356,
+    Rosenbrock23::d1 = 0.43586652150845899941601945119356,
-    Rosenbrock32::d2 = 0.24291996454816804366592249683314,
+    Rosenbrock23::d2 = 0.24291996454816804366592249683314,
-    Rosenbrock32::d3 = 2.1851380027664058511513169485832;
+    Rosenbrock23::d3 = 2.1851380027664058511513169485832;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict)
+Foam::Rosenbrock23::Rosenbrock23(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -79,9 +79,8 @@ Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::Rosenbrock32::solve
+Foam::scalar Foam::Rosenbrock23::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -89,7 +88,7 @@ Foam::scalar Foam::Rosenbrock32::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -117,7 +116,7 @@ Foam::scalar Foam::Rosenbrock32::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -146,15 +145,14 @@ Foam::scalar Foam::Rosenbrock32::solve
 }
 
 
-void Foam::Rosenbrock32::solve
+void Foam::Rosenbrock23::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::Rosenbrock32
+    Foam::Rosenbrock23
 
 Description
     L-stable embedded Rosenbrock ODE solver of order (3)4.
@@ -42,12 +42,12 @@ Description
     \endverbatim
 
 SourceFiles
-    Rosenbrock32.C
+    Rosenbrock23.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Rosenbrock32_H
+#ifndef Rosenbrock23_H
-#define Rosenbrock32_H
+#define Rosenbrock23_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -58,10 +58,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class Rosenbrock32 Declaration
+                           Class Rosenbrock23 Declaration
 \*---------------------------------------------------------------------------*/
 
-class Rosenbrock32
+class Rosenbrock23
 :
     public ODESolver,
     public adaptiveSolver
@@ -91,13 +91,13 @@ class Rosenbrock32
 public:
 
     //- Runtime type information
-    TypeName("Rosenbrock32");
+    TypeName("Rosenbrock23");
 
 
     // Constructors
 
         //- Construct from ODE
-        Rosenbrock32(const ODESystem& ode, const dictionary& dict);
+        Rosenbrock23(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -105,7 +105,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -116,7 +115,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C

@@ -23,89 +23,89 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Rosenbrock43.H"
+#include "Rosenbrock34.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(Rosenbrock43, 0);
+    defineTypeNameAndDebug(Rosenbrock34, 0);
-    addToRunTimeSelectionTable(ODESolver, Rosenbrock43, dictionary);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock34, dictionary);
 
 const scalar
-    // L-Stable constants from Hairer et. al.
-    Rosenbrock43::a21 = 2,
-    Rosenbrock43::a31 = 1.867943637803922,
-    Rosenbrock43::a32 = 0.2344449711399156,
-
-    Rosenbrock43::c21 = -7.137615036412310,
-    Rosenbrock43::c31 = 2.580708087951457,
-    Rosenbrock43::c32 = 0.6515950076447975,
-    Rosenbrock43::c41 = -2.137148994382534,
-    Rosenbrock43::c42 = -0.3214669691237626,
-    Rosenbrock43::c43 = -0.6949742501781779,
-
-    Rosenbrock43::b1 = 2.255570073418735,
-    Rosenbrock43::b2 = 0.2870493262186792,
-    Rosenbrock43::b3 = 0.435317943184018,
-    Rosenbrock43::b4 = 1.093502252409163,
-
-    Rosenbrock43::e1 = -0.2815431932141155,
-    Rosenbrock43::e2 = -0.0727619912493892,
-    Rosenbrock43::e3 = -0.1082196201495311,
-    Rosenbrock43::e4 = -1.093502252409163,
-
-    Rosenbrock43::gamma = 0.57282,
-    Rosenbrock43::c2 = 1.14564,
-    Rosenbrock43::c3 = 0.65521686381559,
-
-    Rosenbrock43::d1 = 0.57282,
-    Rosenbrock43::d2 = -1.769193891319233,
-    Rosenbrock43::d3 = 0.7592633437920482,
-    Rosenbrock43::d4 = -0.1049021087100450;
-
     // Constants by Shampine
     // More accurate than the L-Stable coefficients for small step-size
     // but less stable for large step-size
+    Rosenbrock34::a21 = 2,
+    Rosenbrock34::a31 = 48.0/25.0,
+    Rosenbrock34::a32 = 6.0/25.0,
+
+    Rosenbrock34::c21 = -8,
+    Rosenbrock34::c31 = 372.0/25.0,
+    Rosenbrock34::c32 = 12.0/5.0,
+
+    Rosenbrock34::c41 = -112.0/125.0,
+    Rosenbrock34::c42 = -54.0/125.0,
+    Rosenbrock34::c43 = -2.0/5.0,
+
+    Rosenbrock34::b1 = 19.0/9.0,
+    Rosenbrock34::b2 = 1.0/2.0,
+    Rosenbrock34::b3 = 25.0/108.0,
+    Rosenbrock34::b4 = 125.0/108.0,
+
+    Rosenbrock34::e1 = 34.0/108.0,
+    Rosenbrock34::e2 = 7.0/36.0,
+    Rosenbrock34::e3 = 0,
+    Rosenbrock34::e4 = 125.0/108.0,
+
+    Rosenbrock34::gamma = 1.0/2.0,
+    Rosenbrock34::c2 = 1,
+    Rosenbrock34::c3  = 3.0/5.0,
+
+    Rosenbrock34::d1 = 1.0/2.0,
+    Rosenbrock34::d2 = -3.0/2.0,
+    Rosenbrock34::d3 = 605.0/250.0,
+    Rosenbrock34::d4 = 29.0/250.0;
+
     /*
-    Rosenbrock43::a21 = 2,
+    // L-Stable constants from Hairer et. al.
-    Rosenbrock43::a31 = 48.0/25.0,
+    Rosenbrock34::a21 = 2,
-    Rosenbrock43::a32 = 6.0/25.0,
+    Rosenbrock34::a31 = 1.867943637803922,
+    Rosenbrock34::a32 = 0.2344449711399156,
 
-    Rosenbrock43::c21 = -8,
+    Rosenbrock34::c21 = -7.137615036412310,
-    Rosenbrock43::c31 = 372.0/25.0,
+    Rosenbrock34::c31 = 2.580708087951457,
-    Rosenbrock43::c32 = 12.0/5.0,
+    Rosenbrock34::c32 = 0.6515950076447975,
+    Rosenbrock34::c41 = -2.137148994382534,
+    Rosenbrock34::c42 = -0.3214669691237626,
+    Rosenbrock34::c43 = -0.6949742501781779,
 
-    Rosenbrock43::c41 = -112.0/125.0,
+    Rosenbrock34::b1 = 2.255570073418735,
-    Rosenbrock43::c42 = -54.0/125.0,
+    Rosenbrock34::b2 = 0.2870493262186792,
-    Rosenbrock43::c43 = -2.0/5.0,
+    Rosenbrock34::b3 = 0.435317943184018,
+    Rosenbrock34::b4 = 1.093502252409163,
 
-    Rosenbrock43::b1 = 19.0/9.0,
+    Rosenbrock34::e1 = -0.2815431932141155,
-    Rosenbrock43::b2 = 1.0/2.0,
+    Rosenbrock34::e2 = -0.0727619912493892,
-    Rosenbrock43::b3 = 25.0/108.0,
+    Rosenbrock34::e3 = -0.1082196201495311,
-    Rosenbrock43::b4 = 125.0/108.0,
+    Rosenbrock34::e4 = -1.093502252409163,
 
-    Rosenbrock43::e1 = 34.0/108.0,
+    Rosenbrock34::gamma = 0.57282,
-    Rosenbrock43::e2 = 7.0/36.0,
+    Rosenbrock34::c2 = 1.14564,
-    Rosenbrock43::e3 = 0,
+    Rosenbrock34::c3 = 0.65521686381559,
-    Rosenbrock43::e4 = 125.0/108.0,
 
-    Rosenbrock43::gamma = 1.0/2.0,
+    Rosenbrock34::d1 = 0.57282,
-    Rosenbrock43::c2 = 1,
+    Rosenbrock34::d2 = -1.769193891319233,
-    Rosenbrock43::c3  = 3.0/5.0,
+    Rosenbrock34::d3 = 0.7592633437920482,
-
+    Rosenbrock34::d4 = -0.1049021087100450;
-    Rosenbrock43::d1 = 1.0/2.0,
-    Rosenbrock43::d2 = -3.0/2.0,
-    Rosenbrock43::d3 = 605.0/250.0,
-    Rosenbrock43::d4 = 29.0/250.0;
     */
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict)
+Foam::Rosenbrock34::Rosenbrock34(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -124,9 +124,8 @@ Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::Rosenbrock43::solve
+Foam::scalar Foam::Rosenbrock34::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -134,7 +133,7 @@ Foam::scalar Foam::Rosenbrock43::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -162,7 +161,7 @@ Foam::scalar Foam::Rosenbrock43::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -177,7 +176,7 @@ Foam::scalar Foam::Rosenbrock43::solve
         y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
     }
 
-    ode.derivatives(x0 + c3*dx, y, dydx_);
+    odes_.derivatives(x0 + c3*dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -206,15 +205,14 @@ Foam::scalar Foam::Rosenbrock43::solve
 }
 
 
-void Foam::Rosenbrock43::solve
+void Foam::Rosenbrock34::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H

@@ -22,12 +22,11 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::Rosenbrock43
+    Foam::Rosenbrock34
 
 Description
     L-stable embedded Rosenbrock ODE solver of order (3)4.
 
-    Based on code from:
     \verbatim
         "Solving Ordinary Differential Equations II: Stiff
          and Differential-Algebraic Problems, second edition",
@@ -37,7 +36,7 @@ Description
         Springer-Verlag, Berlin. 1996.
     \endverbatim
 
-    Also provided are the constants from:
+    The default constants are from:
     \verbatim
         "Implementation of Rosenbrock Methods"
         Shampine, L. F.,
@@ -46,13 +45,15 @@ Description
     with which the scheme is more accurate than with the L-Stable coefficients
     for small step-size but less stable for large step-size.
 
+    The L-Stable scheme constants are provided commented-out in Rosenbrock34.C
+
 SourceFiles
-    Rosenbrock43.C
+    Rosenbrock34.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Rosenbrock43_H
+#ifndef Rosenbrock34_H
-#define Rosenbrock43_H
+#define Rosenbrock34_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -63,10 +64,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class Rosenbrock43 Declaration
+                           Class Rosenbrock34 Declaration
 \*---------------------------------------------------------------------------*/
 
-class Rosenbrock43
+class Rosenbrock34
 :
     public ODESolver,
     public adaptiveSolver
@@ -98,13 +99,13 @@ class Rosenbrock43
 public:
 
     //- Runtime type information
-    TypeName("Rosenbrock43");
+    TypeName("Rosenbrock34");
 
 
     // Constructors
 
         //- Construct from ODE
-        Rosenbrock43(const ODESystem& ode, const dictionary& dict);
+        Rosenbrock34(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -112,7 +113,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -123,7 +123,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/SIBS/SIBS.C

@@ -70,13 +70,12 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
 
 void Foam::SIBS::solve
 (
-    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    ode.derivatives(x, y, dydx0_);
+    odes_.derivatives(x, y, dydx0_);
 
     scalar h = dxTry;
     bool exitflag = false;
@@ -124,7 +123,7 @@ void Foam::SIBS::solve
     label k = 0;
     yTemp_ = y;
 
-    ode.jacobian(x, y, dfdx_, dfdy_);
+    odes_.jacobian(x, y, dfdx_, dfdy_);
 
     if (x != xNew_ || h != dxTry)
     {
@@ -151,7 +150,7 @@ void Foam::SIBS::solve
                     << exit(FatalError);
             }
 
-            SIMPR(ode, x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
+            SIMPR(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_);

src/ODE/ODESolvers/SIBS/SIBS.H

@@ -84,7 +84,6 @@ class SIBS
 
         void SIMPR
         (
-            const ODESystem& ode,
             const scalar xStart,
             const scalarField& y,
             const scalarField& dydx,
@@ -123,7 +122,6 @@ public:
 
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/SIBS/SIMPR.C

@@ -29,7 +29,6 @@ License
 
 void Foam::SIBS::SIMPR
 (
-    const ODESystem& ode,
     const scalar xStart,
     const scalarField& y,
     const scalarField& dydx,
@@ -72,7 +71,7 @@ void Foam::SIBS::SIMPR
 
     scalar x = xStart + h;
 
-    ode.derivatives(x, ytemp, yEnd);
+    odes_.derivatives(x, ytemp, yEnd);
 
     for (register label nn=2; nn<=nSteps; nn++)
     {
@@ -90,7 +89,7 @@ void Foam::SIBS::SIMPR
 
         x += h;
 
-        ode.derivatives(x, ytemp, yEnd);
+        odes_.derivatives(x, ytemp, yEnd);
     }
     for (register label i=0; i<n_; i++)
     {

src/ODE/ODESolvers/Trapezoid/Trapezoid.C

@@ -49,7 +49,6 @@ Foam::Trapezoid::Trapezoid(const ODESystem& ode, const dictionary& dict)
 
 Foam::scalar Foam::Trapezoid::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -64,7 +63,7 @@ Foam::scalar Foam::Trapezoid::solve
     }
 
     // Evaluate the system for the predicted state
-    ode.derivatives(x0 + dx, y, err_);
+    odes_.derivatives(x0 + dx, y, err_);
 
     // Update the state as the average between the prediction and the correction
     // and estimate the error from the difference
@@ -80,13 +79,12 @@ Foam::scalar Foam::Trapezoid::solve
 
 void Foam::Trapezoid::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/Trapezoid/Trapezoid.H

@@ -74,7 +74,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -85,7 +84,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C

@@ -64,7 +64,7 @@ void Foam::adaptiveSolver::solve
     do
     {
         // Solve step and provide error estimate
-        err = solve(odes, x, y, dydx0_, dx, yTemp_);
+        err = solve(x, y, dydx0_, dx, yTemp_);
 
         // If error is large reduce dx
         if (err > 1)

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H

@@ -77,7 +77,6 @@ public:
         //- Solve a single step dx and return the error
         virtual scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,

src/ODE/ODESolvers/rodas23/rodas23.C

@@ -23,35 +23,35 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "rodas32.H"
+#include "rodas23.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(rodas32, 0);
+    defineTypeNameAndDebug(rodas23, 0);
-    addToRunTimeSelectionTable(ODESolver, rodas32, dictionary);
+    addToRunTimeSelectionTable(ODESolver, rodas23, dictionary);
 
 const scalar
-    rodas32::c3 = 1,
+    rodas23::c3 = 1,
-    rodas32::d1 = 1.0/2.0,
+    rodas23::d1 = 1.0/2.0,
-    rodas32::d2 = 3.0/2.0,
+    rodas23::d2 = 3.0/2.0,
-    rodas32::a31 = 2,
+    rodas23::a31 = 2,
-    rodas32::a41 = 2,
+    rodas23::a41 = 2,
-    rodas32::c21 = 4,
+    rodas23::c21 = 4,
-    rodas32::c31 = 1,
+    rodas23::c31 = 1,
-    rodas32::c32 = -1,
+    rodas23::c32 = -1,
-    rodas32::c41 = 1,
+    rodas23::c41 = 1,
-    rodas32::c42 = -1,
+    rodas23::c42 = -1,
-    rodas32::c43 = -8.0/3.0,
+    rodas23::c43 = -8.0/3.0,
-    rodas32::gamma = 1.0/2.0;
+    rodas23::gamma = 1.0/2.0;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict)
+Foam::rodas23::rodas23(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -70,9 +70,8 @@ Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::rodas32::solve
+Foam::scalar Foam::rodas23::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -80,7 +79,7 @@ Foam::scalar Foam::rodas32::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -117,7 +116,7 @@ Foam::scalar Foam::rodas32::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -133,7 +132,7 @@ Foam::scalar Foam::rodas32::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(err_, i)
     {
@@ -151,15 +150,14 @@ Foam::scalar Foam::rodas32::solve
 }
 
 
-void Foam::rodas32::solve
+void Foam::rodas23::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/rodas23/rodas23.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::rodas32
+    Foam::rodas23
 
 Description
     L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (2)3.
@@ -42,12 +42,12 @@ Description
     \endverbatim
 
 SourceFiles
-    rodas32.C
+    rodas23.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef rodas32_H
+#ifndef rodas23_H
-#define rodas32_H
+#define rodas23_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -58,10 +58,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class rodas32 Declaration
+                           Class rodas23 Declaration
 \*---------------------------------------------------------------------------*/
 
-class rodas32
+class rodas23
 :
     public ODESolver,
     public adaptiveSolver
@@ -91,13 +91,13 @@ class rodas32
 public:
 
     //- Runtime type information
-    TypeName("rodas32");
+    TypeName("rodas23");
 
 
     // Constructors
 
         //- Construct from ODE
-        rodas32(const ODESystem& ode, const dictionary& dict);
+        rodas23(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -105,7 +105,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -116,7 +115,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/rodas34/rodas34.C

@@ -23,56 +23,56 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "rodas43.H"
+#include "rodas34.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(rodas43, 0);
+    defineTypeNameAndDebug(rodas34, 0);
-    addToRunTimeSelectionTable(ODESolver, rodas43, dictionary);
+    addToRunTimeSelectionTable(ODESolver, rodas34, dictionary);
 
 const scalar
-    rodas43::c2 = 0.386,
+    rodas34::c2 = 0.386,
-    rodas43::c3 = 0.21,
+    rodas34::c3 = 0.21,
-    rodas43::c4 = 0.63,
+    rodas34::c4 = 0.63,
-    rodas43::d1 =  0.25,
+    rodas34::d1 =  0.25,
-    rodas43::d2 = -0.1043,
+    rodas34::d2 = -0.1043,
-    rodas43::d3 =  0.1035,
+    rodas34::d3 =  0.1035,
-    rodas43::d4 = -0.3620000000000023e-01,
+    rodas34::d4 = -0.3620000000000023e-01,
-    rodas43::a21 =  0.1544e1,
+    rodas34::a21 =  0.1544e1,
-    rodas43::a31 =  0.9466785280815826,
+    rodas34::a31 =  0.9466785280815826,
-    rodas43::a32 =  0.2557011698983284,
+    rodas34::a32 =  0.2557011698983284,
-    rodas43::a41 =  0.3314825187068521e1,
+    rodas34::a41 =  0.3314825187068521e1,
-    rodas43::a42 =  0.2896124015972201e1,
+    rodas34::a42 =  0.2896124015972201e1,
-    rodas43::a43 =  0.9986419139977817,
+    rodas34::a43 =  0.9986419139977817,
-    rodas43::a51 =  0.1221224509226641e1,
+    rodas34::a51 =  0.1221224509226641e1,
-    rodas43::a52 =  0.6019134481288629e1,
+    rodas34::a52 =  0.6019134481288629e1,
-    rodas43::a53 =  0.1253708332932087e2,
+    rodas34::a53 =  0.1253708332932087e2,
-    rodas43::a54 = -0.6878860361058950,
+    rodas34::a54 = -0.6878860361058950,
-    rodas43::c21 = -0.56688e1,
+    rodas34::c21 = -0.56688e1,
-    rodas43::c31 = -0.2430093356833875e1,
+    rodas34::c31 = -0.2430093356833875e1,
-    rodas43::c32 = -0.2063599157091915,
+    rodas34::c32 = -0.2063599157091915,
-    rodas43::c41 = -0.1073529058151375,
+    rodas34::c41 = -0.1073529058151375,
-    rodas43::c42 = -0.9594562251023355e1,
+    rodas34::c42 = -0.9594562251023355e1,
-    rodas43::c43 = -0.2047028614809616e2,
+    rodas34::c43 = -0.2047028614809616e2,
-    rodas43::c51 =  0.7496443313967647e1,
+    rodas34::c51 =  0.7496443313967647e1,
-    rodas43::c52 = -0.1024680431464352e2,
+    rodas34::c52 = -0.1024680431464352e2,
-    rodas43::c53 = -0.3399990352819905e2,
+    rodas34::c53 = -0.3399990352819905e2,
-    rodas43::c54 =  0.1170890893206160e2,
+    rodas34::c54 =  0.1170890893206160e2,
-    rodas43::c61 =  0.8083246795921522e1,
+    rodas34::c61 =  0.8083246795921522e1,
-    rodas43::c62 = -0.7981132988064893e1,
+    rodas34::c62 = -0.7981132988064893e1,
-    rodas43::c63 = -0.3152159432874371e2,
+    rodas34::c63 = -0.3152159432874371e2,
-    rodas43::c64 =  0.1631930543123136e2,
+    rodas34::c64 =  0.1631930543123136e2,
-    rodas43::c65 = -0.6058818238834054e1,
+    rodas34::c65 = -0.6058818238834054e1,
-    rodas43::gamma = 0.25;
+    rodas34::gamma = 0.25;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict)
+Foam::rodas34::rodas34(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
     adaptiveSolver(ode, dict),
@@ -93,9 +93,8 @@ Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::rodas43::solve
+Foam::scalar Foam::rodas34::solve
 (
-    const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
@@ -103,7 +102,7 @@ Foam::scalar Foam::rodas43::solve
     scalarField& y
 ) const
 {
-    ode.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, dfdx_, dfdy_);
 
     for (register label i=0; i<n_; i++)
     {
@@ -131,7 +130,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, dydx_);
 
     forAll(k2_, i)
     {
@@ -146,7 +145,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
     }
 
-    ode.derivatives(x0 + c3*dx, y, dydx_);
+    odes_.derivatives(x0 + c3*dx, y, dydx_);
 
     forAll(k3_, i)
     {
@@ -161,7 +160,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
     }
 
-    ode.derivatives(x0 + c4*dx, y, dydx_);
+    odes_.derivatives(x0 + c4*dx, y, dydx_);
 
     forAll(k4_, i)
     {
@@ -178,7 +177,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(k5_, i)
     {
@@ -195,7 +194,7 @@ Foam::scalar Foam::rodas43::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    ode.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, dydx_);
 
     forAll(err_, i)
     {
@@ -214,15 +213,14 @@ Foam::scalar Foam::rodas43::solve
 }
 
 
-void Foam::rodas43::solve
+void Foam::rodas34::solve
 (
-    const ODESystem& odes,
     scalar& x,
     scalarField& y,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, dxTry);
 }
 
 

src/ODE/ODESolvers/rodas34/rodas34.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::rodas43
+    Foam::rodas34
 
 Description
     L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (3)4.
@@ -38,12 +38,12 @@ Description
     \endverbatim
 
 SourceFiles
-    rodas43.C
+    rodas34.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef rodas43_H
+#ifndef rodas34_H
-#define rodas43_H
+#define rodas34_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -54,10 +54,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class rodas43 Declaration
+                           Class rodas34 Declaration
 \*---------------------------------------------------------------------------*/
 
-class rodas43
+class rodas34
 :
     public ODESolver,
     public adaptiveSolver
@@ -92,13 +92,13 @@ class rodas43
 public:
 
     //- Runtime type information
-    TypeName("rodas43");
+    TypeName("rodas34");
 
 
     // Constructors
 
         //- Construct from ODE
-        rodas43(const ODESystem& ode, const dictionary& dict);
+        rodas34(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -106,7 +106,6 @@ public:
         //- Solve a single step dx and return the error
         scalar solve
         (
-            const ODESystem& ode,
             const scalar x0,
             const scalarField& y0,
             const scalarField& dydx0,
@@ -117,7 +116,6 @@ public:
         //- Solve the ODE system and the update the state
         void solve
         (
-            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalar& dxTry

src/ODE/ODESolvers/seulex/seulex.C

@@ -35,14 +35,13 @@ namespace Foam
     addToRunTimeSelectionTable(ODESolver, seulex, dictionary);
 
     const scalar
-        seulex::EPS = VSMALL,
+        seulex::stepFactor1_ = 0.6,
-        seulex::STEPFAC1 = 0.6,
+        seulex::stepFactor2_ = 0.93,
-        seulex::STEPFAC2 = 0.93,
+        seulex::stepFactor3_ = 0.1,
-        seulex::STEPFAC3 = 0.1,
+        seulex::stepFactor4_ = 4.0,
-        seulex::STEPFAC4 = 4.0,
+        seulex::stepFactor5_ = 0.5,
-        seulex::STEPFAC5 = 0.5,
+        seulex::kFactor1_ = 0.7,
-        seulex::KFAC1 = 0.7,
+        seulex::kFactor2_ = 0.9;
-        seulex::KFAC2 = 0.9;
 }
 
 
@@ -51,150 +50,222 @@ namespace Foam
 Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
 :
     ODESolver(ode, dict),
-    ode_(ode),
+    jacRedo_(min(1.0e-4, min(relTol_))),
-    nSeq_(IMAXX),
+    nSeq_(iMaxx_),
-    cost_(IMAXX),
+    cpu_(iMaxx_),
+    coeff_(iMaxx_, iMaxx_),
+    theta_(2.0*jacRedo_),
+    table_(kMaxx_, n_),
     dfdx_(n_),
-    factrl_(IMAXX),
-    table_(KMAXX,n_),
-    fSave_((IMAXX-1)*(IMAXX+1)/2 + 2,n_),
     dfdy_(n_),
-    calcJac_(false),
-    dense_(false),
     a_(n_),
-    coeff_(IMAXX,IMAXX),
+    pivotIndices_(n_),
-    pivotIndices_(n_, 0.0),
+    dxOpt_(iMaxx_),
-    hOpt_(IMAXX),
+    temp_(iMaxx_),
-    work_(IMAXX),
+    y0_(n_),
-    ySaved_(n_),
     ySequence_(n_),
     scale_(n_),
-    del_(n_),
+    dy_(n_),
     yTemp_(n_),
-    dyTemp_(n_),
+    dydx_(n_)
-    delTemp_(n_)
 {
-    const scalar costfunc = 1.0, costjac = 5.0, costlu = 1.0, costsolve = 1.0;
+    // The CPU time factors for the major parts of the algorithm
+    const scalar cpuFunc = 1.0, cpuJac = 5.0, cpuLU = 1.0, cpuSolve = 1.0;
 
-    jacRedo_ = min(1.0e-4, min(relTol_));
-    theta_ = 2.0*jacRedo_;
     nSeq_[0] = 2;
     nSeq_[1] = 3;
 
-    for (label i = 2; i < IMAXX; i++)
+    for (int i=2; i<iMaxx_; i++)
     {
         nSeq_[i] = 2*nSeq_[i-2];
     }
-    cost_[0] = costjac + costlu + nSeq_[0]*(costfunc + costsolve);
+    cpu_[0] = cpuJac + cpuLU + nSeq_[0]*(cpuFunc + cpuSolve);
 
-    for (label k=0;k<KMAXX;k++)
+    for (int k=0; k<kMaxx_; k++)
     {
-        cost_[k+1] = cost_[k] + (nSeq_[k+1]-1)*(costfunc + costsolve) + costlu;
+        cpu_[k+1] = cpu_[k] + (nSeq_[k+1]-1)*(cpuFunc + cpuSolve) + cpuLU;
     }
 
-    hnext_ =- VGREAT;
+    // Set the extrapolation coefficients array
-
+    for (int k=0; k<iMaxx_; k++)
-    //NOTE: the first element of relTol_ and absTol_ are used here.
-    scalar logfact =- log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
-
-    kTarg_ = min(1, min(KMAXX - 1, label(logfact)));
-
-    for (label k = 0; k < IMAXX; k++)
     {
-        for (label l = 0; l < k; l++)
+        for (int l=0; l<k; l++)
         {
             scalar ratio = scalar(nSeq_[k])/nSeq_[l];
             coeff_[k][l] = 1.0/(ratio - 1.0);
         }
     }
-
-    factrl_[0] = 1.0;
-    for (label k = 0; k < IMAXX - 1; k++)
-    {
-        factrl_[k+1] = (k+1)*factrl_[k];
-    }
 }
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-
+bool Foam::seulex::seul
-void Foam::seulex::solve
 (
-    const ODESystem& ode,
+    const scalar x0,
-    scalar& x,
+    const scalarField& y0,
+    const scalar dxTot,
+    const label k,
     scalarField& y,
-    scalar& dxTry
+    const scalarField& scale
 ) const
 {
-    step(dxTry, x, y);
+    label nSteps = nSeq_[k];
-    dxTry = hnext_;
+    scalar dx = dxTot/nSteps;
+
+    for (label i=0; i<n_; i++)
+    {
+        for (label j=0; j<n_; j++)
+        {
+            a_[i][j] = -dfdy_[i][j];
+        }
+
+        a_[i][i] += 1.0/dx;
+    }
+
+    LUDecompose(a_, pivotIndices_);
+
+    scalar xnew = x0 + dx;
+    odes_.derivatives(xnew, y0, dy_);
+    LUBacksubstitute(a_, pivotIndices_, dy_);
+
+    yTemp_ = y0;
+
+    for (label nn=1; nn<nSteps; nn++)
+    {
+        yTemp_ += dy_;
+        xnew += dx;
+
+        if (nn == 1 && k<=1)
+        {
+            scalar dy1 = 0.0;
+            for (label i=0; i<n_; i++)
+            {
+                dy1 += sqr(dy_[i]/scale[i]);
+            }
+            dy1 = sqrt(dy1);
+
+            odes_.derivatives(x0 + dx, yTemp_, dydx_);
+            for (label i=0; i<n_; i++)
+            {
+                dy_[i] = dydx_[i] - dy_[i]/dx;
+            }
+
+            LUBacksubstitute(a_, pivotIndices_, dy_);
+
+            scalar dy2 = 0.0;
+            for (label i=0; i<n_; i++)
+            {
+                dy2 += sqr(dy_[i]/scale[i]);
+            }
+            dy2 = sqrt(dy2);
+            theta_ = dy2/min(1.0, dy1 + SMALL);
+
+            if (theta_ > 1.0)
+            {
+                return false;
+            }
+        }
+
+        odes_.derivatives(xnew, yTemp_, dy_);
+        LUBacksubstitute(a_, pivotIndices_, dy_);
+    }
+
+    for (label i=0; i<n_; i++)
+    {
+        y[i] = yTemp_[i] + dy_[i];
+    }
+
+    return true;
 }
 
 
-void Foam::seulex::step(const scalar& htry, scalar& x,  scalarField& y) const
+void Foam::seulex::extrapolate
+(
+    const label k,
+    scalarRectangularMatrix& table,
+    scalarField& y
+) const
 {
-    static bool first_step = true, last_step = false;
+    for (int j=k-1; j>0; j--)
-    static bool forward, reject = false, prev_reject = false;
-    static scalar errold;
-    label i, k;
-    scalar fac, h, hnew, err;
-    bool firstk;
-
-    work_[0] = 1.e30;
-    h = htry;
-    forward = h > 0 ? true : false;
-
-    ySaved_ = y;
-
-    if (h != hnext_ && !first_step)
     {
-        last_step = true;
+        for (label i=0; i<n_; i++)
+        {
+            table[j-1][i] =
+                table[j][i] + coeff_[k][j]*(table[j][i] - table[j-1][i]);
+        }
     }
 
-    if (reject)
+    for (int i=0; i<n_; i++)
     {
-        prev_reject = true;
+        y[i] = table[0][i] + coeff_[k][0]*(table[0][i] - y[i]);
-        last_step = false;
+    }
-        theta_=2.0*jacRedo_;
+}
+
+
+void Foam::seulex::solve
+(
+    scalar& x,
+    scalarField& y,
+    stepState& step
+) const
+{
+    temp_[0] = GREAT;
+    scalar dx = step.dxTry;
+    y0_ = y;
+    dxOpt_[0] = mag(0.1*dx);
+
+    if (step.first || step.prevReject)
+    {
+        theta_ = 2.0*jacRedo_;
     }
 
-    for (i = 0; i < n_; i++)
+    if (step.first)
+    {
+        // NOTE: the first element of relTol_ and absTol_ are used here.
+        scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
+        kTarg_ = max(1, min(kMaxx_ - 1, int(logTol)));
+    }
+
+    forAll(scale_, i)
     {
         scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
     }
 
-    reject = false;
+    bool jacUpdated = false;
-    firstk = true;
-    hnew = mag(h);
 
-    if (theta_ > jacRedo_ && !calcJac_)
+    if (theta_ > jacRedo_)
     {
-        ode_.jacobian(x, y, dfdx_, dfdy_);
+        odes_.jacobian(x, y, dfdx_, dfdy_);
-        calcJac_ = true;
+        jacUpdated = true;
     }
 
-    while (firstk || reject)
+    int k;
+    scalar dxNew = mag(dx);
+    bool firstk = true;
+
+    while (firstk || step.reject)
     {
-        h = forward ? hnew : -hnew;
+        dx = step.forward ? dxNew : -dxNew;
         firstk = false;
-        reject = false;
+        step.reject = false;
 
-        if (mag(h) <= mag(x)*EPS)
+        if (mag(dx) <= mag(x)*SMALL)
         {
-             WarningIn("seulex::step(const scalar")
+             WarningIn("seulex::solve(scalar& x, scalarField& y, stepState&")
-                    << "step size underflow in step :"  << h << endl;
+                    << "step size underflow :"  << dx << endl;
         }
-        label ipt=-1;
 
-        for (k = 0; k <= kTarg_ + 1; k++)
+        scalar errOld;
+
+        for (k=0; k<=kTarg_+1; k++)
         {
-            bool success = dy(x, ySaved_, h, k, ySequence_, ipt, scale_);
+            bool success = seul(x, y0_, dx, k, ySequence_, scale_);
 
             if (!success)
             {
-                reject = true;
+                step.reject = true;
-                hnew = mag(h)*STEPFAC5;
+                dxNew = mag(dx)*stepFactor5_;
                 break;
             }
 
@@ -204,45 +275,55 @@ void Foam::seulex::step(const scalar& htry, scalar& x,  scalarField& y) const
             }
             else
             {
-                for (i = 0; i < n_; i++)
+                forAll(ySequence_, i)
+                {
                     table_[k-1][i] = ySequence_[i];
+                }
             }
+
             if (k != 0)
             {
-                polyextr(k, table_, y);
+                extrapolate(k, table_, y);
-                err = 0.0;
+                scalar err = 0.0;
-                for (i = 0; i < n_; i++)
+                forAll(scale_, i)
                 {
-                    scale_[i] = absTol_[i] + relTol_[i]*mag(ySaved_[i]);
+                    scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]);
                     err += sqr((y[i] - table_[0][i])/scale_[i]);
                 }
                 err = sqrt(err/n_);
-                if (err > 1.0/EPS || (k > 1 && err >= errold))
+                if (err > 1.0/SMALL || (k > 1 && err >= errOld))
                 {
-                    reject = true;
+                    step.reject = true;
-                    hnew = mag(h)*STEPFAC5;
+                    dxNew = mag(dx)*stepFactor5_;
                     break;
                 }
-                errold = min(4.0*err, 1.0);
+                errOld = min(4*err, 1);
                 scalar expo = 1.0/(k + 1);
-                scalar facmin = pow(STEPFAC3, expo);
+                scalar facmin = pow(stepFactor3_, expo);
+                scalar fac;
                 if (err == 0.0)
                 {
                     fac = 1.0/facmin;
                 }
                 else
                 {
-                    fac = STEPFAC2/pow(err/STEPFAC1, expo);
+                    fac = stepFactor2_/pow(err/stepFactor1_, expo);
-                    fac = max(facmin/STEPFAC4, min(1.0/facmin, fac));
+                    fac = max(facmin/stepFactor4_, min(1.0/facmin, fac));
                 }
-                hOpt_[k] = mag(h*fac);
+                dxOpt_[k] = mag(dx*fac);
-                work_[k] = cost_[k]/hOpt_[k];
+                temp_[k] = cpu_[k]/dxOpt_[k];
 
-                if ((first_step || last_step) && err <= 1.0)
+                if ((step.first || step.last) && err <= 1.0)
                 {
                     break;
                 }
-                if (k == kTarg_-1 && !prev_reject && !first_step && !last_step)
+
+                if
+                (
+                    k == kTarg_ - 1
+                 && !step.prevReject
+                 && !step.first && !step.last
+                )
                 {
                     if (err <= 1.0)
                     {
@@ -250,16 +331,17 @@ void Foam::seulex::step(const scalar& htry, scalar& x,  scalarField& y) const
                     }
                     else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0)
                     {
-                        reject = true;
+                        step.reject = true;
                         kTarg_ = k;
-                        if (kTarg_>1 && work_[k-1] < KFAC1*work_[k])
+                        if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
                         {
                             kTarg_--;
                         }
-                        hnew = hOpt_[kTarg_];
+                        dxNew = dxOpt_[kTarg_];
                         break;
                     }
                 }
+
                 if (k == kTarg_)
                 {
                     if (err <= 1.0)
@@ -268,51 +350,55 @@ void Foam::seulex::step(const scalar& htry, scalar& x,  scalarField& y) const
                     }
                     else if (err > nSeq_[k + 1]*2.0)
                     {
-                        reject = true;
+                        step.reject = true;
-                        if (kTarg_>1 && work_[k-1] < KFAC1*work_[k])
+                        if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
                         {
                             kTarg_--;
                         }
-                        hnew = hOpt_[kTarg_];
+                        dxNew = dxOpt_[kTarg_];
                         break;
                     }
                 }
+
                 if (k == kTarg_+1)
                 {
                     if (err > 1.0)
                     {
-                        reject = true;
+                        step.reject = true;
-                        if (kTarg_ > 1 && work_[kTarg_-1] < KFAC1*work_[kTarg_])
+                        if
+                        (
+                            kTarg_ > 1
+                         && temp_[kTarg_-1] < kFactor1_*temp_[kTarg_]
+                        )
                         {
                             kTarg_--;
                         }
-                        hnew = hOpt_[kTarg_];
+                        dxNew = dxOpt_[kTarg_];
                     }
                     break;
                 }
             }
         }
-        if (reject)
+        if (step.reject)
         {
-            prev_reject = true;
+            step.prevReject = true;
-            if (!calcJac_)
+            if (!jacUpdated)
             {
                 theta_ = 2.0*jacRedo_;
 
-                if (theta_ > jacRedo_ && !calcJac_)
+                if (theta_ > jacRedo_ && !jacUpdated)
                 {
-                    ode_.jacobian(x, y, dfdx_, dfdy_);
+                    odes_.jacobian(x, y, dfdx_, dfdy_);
-                    calcJac_ = true;
+                    jacUpdated = true;
                 }
             }
         }
     }
 
-    calcJac_ = false;
+    jacUpdated = false;
 
-    x += h;
+    step.dxDid = dx;
-    hdid_ = h;
+    x += dx;
-    first_step = false;
 
     label kopt;
     if (k == 1)
@@ -322,180 +408,56 @@ void Foam::seulex::step(const scalar& htry, scalar& x,  scalarField& y) const
     else if (k <= kTarg_)
     {
         kopt=k;
-        if (work_[k-1] < KFAC1*work_[k])
+        if (temp_[k-1] < kFactor1_*temp_[k])
         {
             kopt = k - 1;
         }
-        else if (work_[k] < KFAC2*work_[k - 1])
+        else if (temp_[k] < kFactor2_*temp_[k - 1])
         {
-            kopt = min(k + 1, KMAXX - 1);
+            kopt = min(k + 1, kMaxx_ - 1);
         }
     }
     else
     {
         kopt = k - 1;
-        if (k > 2 && work_[k-2] < KFAC1*work_[k - 1])
+        if (k > 2 && temp_[k-2] < kFactor1_*temp_[k - 1])
         {
             kopt = k - 2;
         }
-        if (work_[k] < KFAC2*work_[kopt])
+        if (temp_[k] < kFactor2_*temp_[kopt])
         {
-            kopt = min(k, KMAXX - 1);
+            kopt = min(k, kMaxx_ - 1);
         }
     }
 
-    if (prev_reject)
+    if (step.prevReject)
     {
         kTarg_ = min(kopt, k);
-        hnew = min(mag(h), hOpt_[kTarg_]);
+        dxNew = min(mag(dx), dxOpt_[kTarg_]);
-        prev_reject = false;
+        step.prevReject = false;
     }
     else
     {
         if (kopt <= k)
         {
-            hnew = hOpt_[kopt];
+            dxNew = dxOpt_[kopt];
         }
         else
         {
-            if (k < kTarg_ && work_[k] < KFAC2*work_[k - 1])
+            if (k < kTarg_ && temp_[k] < kFactor2_*temp_[k - 1])
             {
-                hnew = hOpt_[k]*cost_[kopt + 1]/cost_[k];
+                dxNew = dxOpt_[k]*cpu_[kopt + 1]/cpu_[k];
             }
             else
             {
-                hnew = hOpt_[k]*cost_[kopt]/cost_[k];
+                dxNew = dxOpt_[k]*cpu_[kopt]/cpu_[k];
             }
         }
         kTarg_ = kopt;
     }
 
-    if (forward)
+    step.dxTry = step.forward ? dxNew : -dxNew;
-    {
-        hnext_ = hnew;
-    }
-    else
-    {
-        hnext_ =- hnew;
-    }
 }
 
 
-bool Foam::seulex::dy
-(
-    const scalar& x,
-    scalarField& y,
-    const scalar htot,
-    const label k,
-    scalarField& yend,
-    label& ipt,
-    scalarField& scale
-) const
-{
-    label nstep = nSeq_[k];
-    scalar h = htot/nstep;
-    for (label i = 0; i < n_; i++)
-    {
-        for (label j = 0; j < n_; j++) a_[i][j] = -dfdy_[i][j];
-        a_[i][i] += 1.0/h;
-    }
-    LUDecompose(a_, pivotIndices_);
-    scalar xnew = x + h;
-
-    ode_.derivatives(xnew, y, del_);
-
-    yTemp_ = y;
-
-    LUBacksubstitute(a_, pivotIndices_, del_);
-
-    if (dense_ && nstep == k + 1)
-    {
-        ipt++;
-        for (label i = 0; i < n_; i++)
-        {
-            fSave_[ipt][i] = del_[i];
-        }
-    }
-    for (label nn = 1; nn < nstep; nn++)
-    {
-        for (label i=0;i<n_;i++)
-        {
-            yTemp_[i] += del_[i];
-        }
-        xnew += h;
-
-        ode_.derivatives(xnew, yTemp_, yend);
-
-        if (nn == 1 && k<=1)
-        {
-            scalar del1=0.0;
-            for (label i = 0; i < n_; i++)
-            {
-                del1 += sqr(del_[i]/scale[i]);
-            }
-            del1 = sqrt(del1);
-
-            ode_.derivatives(x+h, yTemp_, dyTemp_);
-            for (label i=0;i<n_;i++)
-            {
-                del_[i] = dyTemp_[i] - del_[i]/h;
-            }
-
-            LUBacksubstitute(a_, pivotIndices_, del_);
-
-            scalar del2 = 0.0;
-            for (label i = 0; i <n_ ; i++)
-            {
-                del2 += sqr(del_[i]/scale[i]);
-            }
-            del2 = sqrt(del2);
-            theta_ = del2 / min(1.0, del1 + SMALL);
-
-            if (theta_ > 1.0)
-            {
-                return false;
-            }
-        }
-
-        delTemp_ = yend;
-        LUBacksubstitute(a_, pivotIndices_, delTemp_);
-        del_ = delTemp_;
-
-        if (dense_ && nn >= nstep-k-1)
-        {
-            ipt++;
-            for (label i=0;i<n_;i++)
-            {
-                fSave_[ipt][i]=del_[i];
-            }
-        }
-    }
-
-    yend = yTemp_ + del_;
-    return true;
-}
-
-
-void Foam::seulex::polyextr
-(
-    const label k,
-    scalarRectangularMatrix& table,
-    scalarField& last
-) const
-{
-    label l=last.size();
-    for (label j = k - 1; j > 0; j--)
-    {
-        for (label i=0; i<l; i++)
-        {
-            table[j-1][i] =
-                table[j][i] + coeff_[k][j]*(table[j][i] - table[j-1][i]);
-        }
-    }
-    for (label i = 0; i < l; i++)
-    {
-        last[i] = table[0][i] + coeff_[k][0]*(table[0][i] - last[i]);
-    }
-}
-
 // ************************************************************************* //

src/ODE/ODESolvers/seulex/seulex.H

@@ -28,7 +28,7 @@ Description
     An extrapolation-algorithm, based on the linearly implicit Euler method
     with step size control and order selection.
 
-    The implementation is based on the SEULEX code in
+    The implementation is based on the SEULEX algorithm in
     \verbatim
         "Solving Ordinary Differential Equations II: Stiff
          and Differential-Algebraic Problems, second edition",
@@ -47,8 +47,6 @@ SourceFiles
 #define seulex_H
 
 #include "ODESolver.H"
-
-#include "scalarFieldField.H"
 #include "scalarMatrices.H"
 #include "labelField.H"
 
@@ -65,64 +63,65 @@ class seulex
 :
     public ODESolver
 {
-
     // Private data
 
-        //- seulex constants
+        // Static constants
-        static const scalar EPS;
 
-        static const label KMAXX=12, IMAXX = KMAXX + 1;
+            static const label kMaxx_ = 12;
+            static const label iMaxx_ = kMaxx_ + 1;
 
-        static const scalar
+            static const scalar
-            STEPFAC1, STEPFAC2, STEPFAC3, STEPFAC4, STEPFAC5, KFAC1, KFAC2;
+                stepFactor1_, stepFactor2_, stepFactor3_,
+                stepFactor4_, stepFactor5_,
+                kFactor1_, kFactor2_;
 
+        // Evaluated constants
 
-        //- Reference to ODESystem
+            scalar jacRedo_;
-        const ODESystem& ode_;
+            labelField nSeq_;
+            scalarField cpu_;
+            scalarSquareMatrix coeff_;
 
-        // Temporary fields
+        // Temporary storage
-        mutable label kTarg_;
+        // held to avoid dynamic memory allocation between calls
-        mutable labelField nSeq_;
+        // and to transfer internal values between functions
-        mutable scalarField cost_, dfdx_, factrl_;
-        mutable scalarRectangularMatrix table_, fSave_;
-        mutable scalarSquareMatrix dfdy_;
-        mutable scalar jacRedo_, theta_;
-        mutable bool calcJac_, dense_;
-        mutable scalarSquareMatrix a_, coeff_;
 
+            mutable scalar theta_;
+            mutable label kTarg_;
+            mutable scalarRectangularMatrix table_;
 
-        mutable scalar hnext_, hdid_;
+            mutable scalarField dfdx_;
+            mutable scalarSquareMatrix dfdy_;
+            mutable scalarSquareMatrix a_;
+            mutable labelList pivotIndices_;
 
-        mutable labelList pivotIndices_;
+            // Fields space for "solve" function
+            mutable scalarField dxOpt_, temp_;
+            mutable scalarField y0_, ySequence_, scale_;
 
-        // Fields space for "step" function
+            // Fields used in "seul" function
-        mutable scalarField hOpt_ ,work_, ySaved_, ySequence_, scale_;
+            mutable scalarField dy_, yTemp_, dydx_;
-
-        // Fields used in "dy" function
-        mutable scalarField del_, yTemp_, dyTemp_, delTemp_;
 
 
     // Private Member Functions
 
-        void step(const scalar& htry, scalar& x,  scalarField& y) const;
+        //- Computes the j-th line of the extrapolation table
-
+        bool seul
-        bool dy
         (
-            const scalar& x,
+            const scalar x0,
-            scalarField& y,
+            const scalarField& y0,
-            const scalar htot,
+            const scalar dxTot,
             const label k,
-            scalarField& yend,
+            scalarField& y,
-            label& ipt,
+            const scalarField& scale
-            scalarField& scale
         ) const;
 
-
+        //- Polynomial extrpolation
-        void polyextr
+        void extrapolate
         (
             const label k,
             scalarRectangularMatrix& table,
-            scalarField& last
+            scalarField& y
         ) const;
 
 
@@ -139,13 +138,14 @@ public:
 
 
     // Member Functions
-    void solve
+
-    (
+        //- Solve the ODE system and the update the state
-        const ODESystem& ode,
+        void solve
-        scalar& x,
+        (
-        scalarField& y,
+            scalar& x,
-        scalar& dxTry
+            scalarField& y,
-    ) const;
+            stepState& step
+        ) const;
 };
 
 

src/postProcessing/functionObjects/field/Make/options

@@ -6,15 +6,16 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
+    -lsurfMesh \
     -llagrangian \
     -lfileFormats \
     -lsampling \
     -lincompressibleTransportModels \
     -lcompressibleTurbulenceModel \
     -lincompressibleTurbulenceModel
-

src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.C

@@ -33,7 +33,7 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(nearWallFields, 0);
+    defineTypeNameAndDebug(nearWallFields, 0);
 }
 
 
@@ -72,8 +72,8 @@ void Foam::nearWallFields::calcAddressing()
         label patchI = iter.key();
         const fvPatch& patch = mesh.boundary()[patchI];
 
-        vectorField nf = patch.nf();
+        vectorField nf(patch.nf());
-        vectorField faceCellCentres = patch.patch().faceCellCentres();
+        vectorField faceCellCentres(patch.patch().faceCellCentres());
 
         forAll(patch, patchFaceI)
         {

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

@@ -70,14 +70,7 @@ void Foam::ode<ChemistryModel>::solve
     cTp_[nSpecie] = T;
     cTp_[nSpecie+1] = p;
 
-    odeSolver_->solve
+    odeSolver_->solve(0, deltaT, cTp_, subDeltaT);
-    (
-        *this,
-        0,
-        deltaT,
-        cTp_,
-        subDeltaT
-    );
 
     for (register int i=0; i<nSpecie; i++)
     {

tutorials/combustion/chemFoam/README

@@ -8,4 +8,4 @@ ic8h18 : iso-Octane combustion, 874 species, 3796 reactions
 Results interpreted in 'validation' sub-folder, where OpenFOAM results
 are compared against those predicted by CHEMKIN II.
 
-Overall the best performing ODE solver is seulex followed closely by rodas32.
+Overall the best performing ODE solver is seulex followed closely by rodas23.

tutorials/combustion/chemFoam/h2/constant/chemistryProperties

@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 
 odeCoeffs
 {
-    solver          seulex; //rodas32;
+    solver          seulex;
     absTol          1e-12;
     relTol          1e-1;
 }

tutorials/multiphase/interDyMFoam/ras/floatingObject/system/decomposeParDict

@@ -15,9 +15,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains 4;
+numberOfSubdomains 8;
 
-method          hierarchical;
+method          scotch;
 
 simpleCoeffs
 {