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

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,8 +2,15 @@ ODESolvers/ODESolver/ODESolver.C
 ODESolvers/ODESolver/ODESolverNew.C
 
 ODESolvers/adaptiveSolver/adaptiveSolver.C
-ODESolvers/RKCK5/RKCK5.C
+ODESolvers/Euler/Euler.C
-ODESolvers/KRR4/KRR4.C
+ODESolvers/EulerSI/EulerSI.C
+ODESolvers/Trapezoid/Trapezoid.C
+ODESolvers/RKF45/RKF45.C
+ODESolvers/RKCK45/RKCK45.C
+ODESolvers/RKDP45/RKDP45.C
+ODESolvers/Rosenbrock21/Rosenbrock21.C
+ODESolvers/Rosenbrock43/Rosenbrock43.C
+ODESolvers/rodas43/rodas43.C
 ODESolvers/SIBS/SIBS.C
 ODESolvers/SIBS/SIMPR.C
 ODESolvers/SIBS/polyExtrapolate.C

src/ODE/ODESolvers/Euler/Euler.C

@@ -0,0 +1,88 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "Euler.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(Euler, 0);
+    addToRunTimeSelectionTable(ODESolver, Euler, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::Euler::Euler(const ODESystem& ode, const dictionary& dict)
+:
+    ODESolver(ode, dict),
+    adaptiveSolver(ode, dict),
+    err_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::Euler::solve
+(
+    const ODESystem& ode,
+    const scalar x0,
+    const scalarField& y0,
+    const scalarField& dydx0,
+    const scalar dx,
+    scalarField& y
+) const
+{
+    // Calculate error estimate from the change in state:
+    forAll(err_, i)
+    {
+        err_[i] = dx*dydx0[i];
+    }
+
+    // Update the state
+    forAll(y, i)
+    {
+        y[i] = y0[i] + err_[i];
+    }
+
+    return normalizeError(y0, y, err_);
+}
+
+
+void Foam::Euler::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/Euler/Euler.H

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::Euler
+
+Description
+    Euler ODE solver of order (0)1.
+
+    The method calculates the new state from:
+    \f[
+        y_{n+1} = y_n + \delta_x f
+    \f]
+    The error is estimated directly from the change in the solution,
+    i.e. the difference between the 0th and 1st order solutions:
+    \f[
+        err_{n+1} = y_{n+1} - y_n
+    \f]
+
+SourceFiles
+    Euler.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Euler_H
+#define Euler_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Euler Declaration
+\*---------------------------------------------------------------------------*/
+
+class Euler
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        mutable scalarField err_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Euler");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        Euler(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/EulerSI/EulerSI.C

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "EulerSI.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(EulerSI, 0);
+    addToRunTimeSelectionTable(ODESolver, EulerSI, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::EulerSI::EulerSI(const ODESystem& ode, const dictionary& dict)
+:
+    ODESolver(ode, dict),
+    adaptiveSolver(ode, dict),
+    err_(n_),
+    dydx_(n_),
+    dfdx_(n_),
+    dfdy_(n_, n_),
+    a_(n_, n_),
+    pivotIndices_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::EulerSI::solve
+(
+    const ODESystem& ode,
+    const scalar x0,
+    const scalarField& y0,
+    const scalarField& dydx0,
+    const scalar dx,
+    scalarField& y
+) const
+{
+    ode.jacobian(x0, y0, dfdx_, dfdy_);
+
+    for (register label i=0; i<n_; i++)
+    {
+        for (register label j=0; j<n_; j++)
+        {
+            a_[i][j] = -dfdy_[i][j];
+        }
+
+        a_[i][i] += 1.0/dx;
+    }
+
+    LUDecompose(a_, pivotIndices_);
+
+    // Calculate error estimate from the change in state:
+    forAll(err_, i)
+    {
+        err_[i] = dydx0[i] + dx*dfdx_[i];
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, err_);
+
+    forAll(y, i)
+    {
+        y[i] = y0[i] + err_[i];
+    }
+
+    return normalizeError(y0, y, err_);
+}
+
+
+void Foam::EulerSI::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/EulerSI/EulerSI.H

@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::EulerSI
+
+Description
+    Semi-implicit Euler ODE solver of order (0)1.
+
+    The method calculates the new state from:
+    \f[
+        y_{n+1} = y_n
+          + \delta_x\left[I - \delta_x\frac{\partial f}{\partial y}\right]^{-1}
+            \cdot \left[f(y_n) + \delta_x\frac{\partial f}{\partial x}\right]
+    \f]
+    The error is estimated directly from the change in the solution,
+    i.e. the difference between the 0th and 1st order solutions:
+    \f[
+        err_{n+1} = y_{n+1} - y_n
+    \f]
+
+SourceFiles
+    EulerSI.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef EulerSI_H
+#define EulerSI_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class EulerSI Declaration
+\*---------------------------------------------------------------------------*/
+
+class EulerSI
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        mutable scalarField err_;
+        mutable scalarField dydx_;
+        mutable scalarField dfdx_;
+        mutable scalarSquareMatrix dfdy_;
+        mutable scalarSquareMatrix a_;
+        mutable labelList pivotIndices_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("EulerSI");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        EulerSI(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/ODESolver/ODESolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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-4)),
+    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
@@ -36,13 +36,23 @@ Description
         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.,
+        Nørsett, S.,
+        Wanner, G.,
+        Springer-Verlag, Berlin. 1993, ISBN 3-540-56670-8.
+    \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 +63,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class RKCK5 Declaration
+                           Class RKCK45 Declaration
 \*---------------------------------------------------------------------------*/
 
-class RKCK5
+class RKCK45
 :
     public ODESolver,
     public adaptiveSolver
@@ -65,11 +75,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 +96,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 +115,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 +124,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/RKDP45/RKDP45.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 "RKDP45.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(RKDP45, 0);
+    addToRunTimeSelectionTable(ODESolver, RKDP45, dictionary);
+
+const scalar
+    RKDP45::c2 = 1.0/5.0,
+    RKDP45::c3 = 3.0/10.0,
+    RKDP45::c4 = 4.0/5.0,
+    RKDP45::c5 = 8.0/9.0,
+
+    RKDP45::a21 = 1.0/5.0,
+
+    RKDP45::a31 = 3.0/40.0,
+    RKDP45::a32 = 9.0/40.0,
+
+    RKDP45::a41 = 44.0/45.0,
+    RKDP45::a42 = -56.0/15.0,
+    RKDP45::a43 = 32.0/9.0,
+
+    RKDP45::a51 = 19372.0/6561.0,
+    RKDP45::a52 = -25360.0/2187.0,
+    RKDP45::a53 = 64448.0/6561.0,
+    RKDP45::a54 = -212.0/729.0,
+
+    RKDP45::a61 = 9017.0/3168.0,
+    RKDP45::a62 = -355.0/33.0,
+    RKDP45::a63 = 46732.0/5247.0,
+    RKDP45::a64 = 49.0/176.0,
+    RKDP45::a65 = -5103.0/18656.0,
+
+    RKDP45::b1 = 35.0/384.0,
+    RKDP45::b3 = 500.0/1113.0,
+    RKDP45::b4 = 125.0/192.0,
+    RKDP45::b5 = -2187.0/6784.0,
+    RKDP45::b6 = 11.0/84.0,
+
+    RKDP45::e1 = 5179.0/57600.0 - RKDP45::b1,
+    RKDP45::e3 = 7571.0/16695.0 - RKDP45::b3,
+    RKDP45::e4 = 393.0/640.0 - RKDP45::b4,
+    RKDP45::e5 = -92097.0/339200.0 - RKDP45::b5,
+    RKDP45::e6 = 187.0/2100.0 - RKDP45::b6,
+    RKDP45::e7 = 1.0/40.0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::RKDP45::RKDP45(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::RKDP45::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 + dx, yTemp_, k6_);
+
+    forAll(y, i)
+    {
+        y[i] = y0[i]
+          + dx*(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b5*k5_[i] + b6*k6_[i]);
+    }
+
+    // Reuse k2_ for the derivative of the new state
+    odes.derivatives(x0 + dx, y, k2_);
+
+    forAll(err_, i)
+    {
+        err_[i] =
+            dx
+           *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[i]
+             + e7*k2_[i]);
+    }
+
+    return normalizeError(y0, y, err_);
+}
+
+
+void Foam::RKDP45::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/RKDP45/RKDP45.H

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::RKDP45
+
+Description
+    4/5th Order Dormand–Prince Runge-Kutta ODE solver.
+
+    References:
+    \verbatim
+        "A family of embedded Runge-Kutta formulae"
+        Dormand, J. R.,
+        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.,
+        Nørsett, S.,
+        Wanner, G.,
+        Springer-Verlag, Berlin. 1993, ISBN 3-540-56670-8.
+    \endverbatim
+
+SeeAlso
+    Foam::RKF45
+    Foam::RKCK45
+
+SourceFiles
+    RKDP45.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RKDP45_H
+#define RKDP45_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class RKDP45 Declaration
+\*---------------------------------------------------------------------------*/
+
+class RKDP45
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        //- RKDP Constants
+        static const scalar
+            c2, c3, c4, c5,
+            a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
+            a61, a62, a63, a64, a65,
+            b1, b3, b4, b5, b6,
+            e1, e3, e4, e5, e6, e7;
+
+        // 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("RKDP45");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        RKDP45(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/RKF45/RKF45.C

@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::b1  = 16.0/135.0,
+    RKF45::b3  = 6656.0/12825.0,
+    RKF45::b4  = 28561.0/56430.0,
+    RKF45::b5  = -9.0/50.0,
+    RKF45::b6  = 2.0/55.0,
+
+    RKF45::e1  = 25.0/216.0 - RKF45::b1,
+    RKF45::e3  = 1408.0/2565.0 - RKF45::b3,
+    RKF45::e4  = 2197.0/4104.0 - RKF45::b4,
+    RKF45::e5  = -1.0/5.0 - RKF45::b5,
+    RKF45::e6  = -RKF45::b6;
+}
+
+
+// * * * * * * * * * * * * * * * * 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
+           *(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b5*k5_[i] + b6*k6_[i]);
+    }
+
+    // Calculate the error estimate from the difference between the
+    // 4th-order and 5th-order solutions
+    forAll(err_, i)
+    {
+        err_[i] =
+            dx
+           *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[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,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    Based on code from:
+    \verbatim
+        "Solving Ordinary Differential Equations I: Nonstiff Problems,
+         second edition",
+        Hairer, E.,
+        Nørsett, S.,
+        Wanner, G.,
+        Springer-Verlag, Berlin. 1993, ISBN 3-540-56670-8.
+    \endverbatim
+
+    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.
+
+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,
+            b1, b3, b4, b5, b6,
+            e1, e3, e4, e5, e6;
+
+
+        // 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/Rosenbrock21/Rosenbrock21.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,42 +23,40 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "KRR4.H"
+#include "Rosenbrock21.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(KRR4, 0);
+    defineTypeNameAndDebug(Rosenbrock21, 0);
-    addToRunTimeSelectionTable(ODESolver, KRR4, ODE);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock21, dictionary);
 
 const scalar
-    KRR4::gamma = 1.0/2.0,
+    Rosenbrock21::gamma = 1 + 1.0/std::sqrt(2.0),
-    KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
+    Rosenbrock21::a21 = 1.0/gamma,
-    KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0,
+    Rosenbrock21::c2 = 1.0,
-    KRR4::c41 = -112.0/125.0, KRR4::c42 = -54.0/125.0, KRR4::c43 = -2.0/5.0,
+    Rosenbrock21::c21 = -2.0/gamma,
-    KRR4::b1 = 19.0/9.0, KRR4::b2 = 1.0/2.0, KRR4::b3 = 25.0/108.0,
+    Rosenbrock21::b1 = (3.0/2.0)/gamma,
-    KRR4::b4 = 125.0/108.0,
+    Rosenbrock21::b2 = (1.0/2.0)/gamma,
-    KRR4::e1 = 17.0/54.0, KRR4::e2 = 7.0/36.0, KRR4::e3 = 0.0,
+    Rosenbrock21::e1 = b1 - 1.0/gamma,
-    KRR4::e4 = 125.0/108.0,
+    Rosenbrock21::e2 = b2,
-    KRR4::c1X = 1.0/2.0, KRR4::c2X = -3.0/2.0, KRR4::c3X = 121.0/50.0,
+    Rosenbrock21::d1 = 1,
-    KRR4::c4X = 29.0/250.0,
+    Rosenbrock21::d2 = -1;
-    KRR4::a2X = 1.0, KRR4::a3X = 3.0/5.0;
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::KRR4::KRR4(const ODESystem& ode)
+Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict)
 :
-    ODESolver(ode),
+    ODESolver(ode, dict),
-    adaptiveSolver(ode),
+    adaptiveSolver(ode, dict),
-    g1_(n_),
+    k1_(n_),
-    g2_(n_),
+    k2_(n_),
-    g3_(n_),
-    g4_(n_),
     err_(n_),
+    dydx_(n_),
     dfdx_(n_),
     dfdy_(n_, n_),
     a_(n_, n_),
@@ -68,15 +66,14 @@ Foam::KRR4::KRR4(const ODESystem& ode)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::KRR4::solve
+Foam::scalar Foam::Rosenbrock21::solve
 (
     const ODESystem& ode,
     const scalar x0,
     const scalarField& y0,
     const scalarField& dydx0,
     const scalar dx,
-    scalarField& y,
+    scalarField& y
-    const scalar eps
 ) const
 {
     ode.jacobian(x0, y0, dfdx_, dfdy_);
@@ -93,72 +90,49 @@ Foam::scalar Foam::KRR4::solve
 
     LUDecompose(a_, pivotIndices_);
 
-    forAll(g1_, i)
+    // Calculate k1:
+    forAll(k1_, i)
     {
-        g1_[i] = dydx0[i] + dx*c1X*dfdx_[i];
+        k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
     }
 
-    LUBacksubstitute(a_, pivotIndices_, g1_);
+    LUBacksubstitute(a_, pivotIndices_, k1_);
 
+    // Calculate k2:
     forAll(y, i)
     {
-        y[i] = y0[i] + a21*g1_[i];
+        y[i] = y0[i] + a21*k1_[i];
     }
 
-    ode.derivatives(x0 + a2X*dx, y, dydx_);
+    ode.derivatives(x0 + c2*dx, y, dydx_);
 
-    forAll(g2_, i)
+    forAll(k2_, i)
     {
-        g2_[i] = dydx_[i] + dx*c2X*dfdx_[i] + c21*g1_[i]/dx;
+        k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
     }
 
-    LUBacksubstitute(a_, pivotIndices_, g2_);
+    LUBacksubstitute(a_, pivotIndices_, k2_);
 
+    // Calculate error and update state:
     forAll(y, i)
     {
-        y[i] = y0[i] + a31*g1_[i] + a32*g2_[i];
+        y[i] = y0[i] + b1*k1_[i] + b2*k2_[i];
+        err_[i] = e1*k1_[i] + e2*k2_[i];
     }
 
-    ode.derivatives(x0 + a3X*dx, y, dydx_);
+    return normalizeError(y0, y, err_);
-
-    forAll(g3_, i)
-    {
-        g3_[i] = dydx0[i] + dx*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/dx;
-    }
-
-    LUBacksubstitute(a_, pivotIndices_, g3_);
-
-    forAll(g4_, i)
-    {
-        g4_[i] = dydx_[i] + dx*c4X*dfdx_[i]
-           + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/dx;
-    }
-
-    LUBacksubstitute(a_, pivotIndices_, g4_);
-
-    forAll(y, i)
-    {
-        y[i] = y0[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];
-        err_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];
-    }
-
-    return normalizeError(eps, y0, y, err_);
 }
 
 
-void Foam::KRR4::solve
+void Foam::Rosenbrock21::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/Rosenbrock21/Rosenbrock21.H

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::Rosenbrock21
+
+Description
+    Embedded Rosenbrock ODE solver of order (1)2.
+
+SourceFiles
+    Rosenbrock21.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Rosenbrock21_H
+#define Rosenbrock21_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Rosenbrock21 Declaration
+\*---------------------------------------------------------------------------*/
+
+class Rosenbrock21
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        mutable scalarField k1_;
+        mutable scalarField k2_;
+        mutable scalarField err_;
+        mutable scalarField dydx_;
+        mutable scalarField dfdx_;
+        mutable scalarSquareMatrix dfdy_;
+        mutable scalarSquareMatrix a_;
+        mutable labelList pivotIndices_;
+
+        static const scalar
+            a21,
+            c21,
+            b1, b2,
+            gamma,
+            c2,
+            e1, e2,
+            d1, d2;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Rosenbrock21");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        Rosenbrock21(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/Rosenbrock43/Rosenbrock43.C

@@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Rosenbrock43.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(Rosenbrock43, 0);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock43, dictionary);
+
+const scalar
+
+    Rosenbrock43::a21 = 2.0,
+    Rosenbrock43::a31 = 48.0/25.0,
+    Rosenbrock43::a32 = 6.0/25.0,
+
+    Rosenbrock43::c21 = -8.0,
+    Rosenbrock43::c31 = 372.0/25.0,
+    Rosenbrock43::c32 = 12.0/5.0,
+
+    Rosenbrock43::c41 = -112.0/125.0,
+    Rosenbrock43::c42 = -54.0/125.0,
+    Rosenbrock43::c43 = -2.0/5.0,
+
+    Rosenbrock43::b1 = 19.0/9.0,
+    Rosenbrock43::b2 = 1.0/2.0,
+    Rosenbrock43::b3 = 25.0/108.0,
+    Rosenbrock43::b4 = 125.0/108.0,
+
+    Rosenbrock43::e1 = 34.0/108.0,
+    Rosenbrock43::e2 = 7.0/36.0,
+    Rosenbrock43::e3 = 0.0,
+    Rosenbrock43::e4 = 125.0/108.0,
+
+    Rosenbrock43::gamma = 0.5,
+    Rosenbrock43::c2 = 1.0,
+    Rosenbrock43::c3  = 3.0/5.0,
+
+    Rosenbrock43::d1 = 1.0/2.0,
+    Rosenbrock43::d2 = -3.0/2.0,
+    Rosenbrock43::d3 = 605.0/250.0,
+    Rosenbrock43::d4 = 29.0/250.0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict)
+:
+    ODESolver(ode, dict),
+    adaptiveSolver(ode, dict),
+    k1_(n_),
+    k2_(n_),
+    k3_(n_),
+    k4_(n_),
+    err_(n_),
+    dydx_(n_),
+    dfdx_(n_),
+    dfdy_(n_, n_),
+    a_(n_, n_),
+    pivotIndices_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::Rosenbrock43::solve
+(
+    const ODESystem& ode,
+    const scalar x0,
+    const scalarField& y0,
+    const scalarField& dydx0,
+    const scalar dx,
+    scalarField& y
+) const
+{
+    ode.jacobian(x0, y0, dfdx_, dfdy_);
+
+    for (register label i=0; i<n_; i++)
+    {
+        for (register label j=0; j<n_; j++)
+        {
+            a_[i][j] = -dfdy_[i][j];
+        }
+
+        a_[i][i] += 1.0/(gamma*dx);
+    }
+
+    LUDecompose(a_, pivotIndices_);
+
+    // Calculate k1:
+    forAll(k1_, i)
+    {
+        k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k1_);
+
+    // Calculate k2:
+    forAll(y, i)
+    {
+        y[i] = y0[i] + a21*k1_[i];
+    }
+
+    ode.derivatives(x0 + c2*dx, y, dydx_);
+
+    forAll(k2_, i)
+    {
+        k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k2_);
+
+    // Calculate k3:
+    forAll(y, i)
+    {
+        y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
+    }
+
+    ode.derivatives(x0 + c3*dx, y, dydx_);
+
+    forAll(k3_, i)
+    {
+        k3_[i] = dydx_[i] + dx*d3*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k3_);
+
+    // Calculate k4:
+    forAll(k4_, i)
+    {
+        k4_[i] = dydx_[i] + dx*d4*dfdx_[i]
+          + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k4_);
+
+    // Calculate error and update state:
+    forAll(y, i)
+    {
+        y[i] = y0[i] + b1*k1_[i] + b2*k2_[i] + b3*k3_[i] + b4*k4_[i];
+        err_[i] = e1*k1_[i] + e2*k2_[i] + e4*k4_[i];
+    }
+
+    return normalizeError(y0, y, err_);
+}
+
+
+void Foam::Rosenbrock43::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.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,18 +22,19 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::KRR4
+    Foam::Rosenbrock43
 
 Description
-    4th Order Kaps Rentrop Rosenbrock stiff ODE solver
+    Embedded Rosenbrock ODE solver of order (3)4.
 
-    References:
+    Based on code from:
     \verbatim
-        "Generalized Runge-Kutta methods of order four with stepsize control
+        "Solving Ordinary Differential Equations II: Stiff
-         for stiff ordinary differential equations"
+         and Differential-Algebraic Problems, second edition",
-        Kaps, P.,
+        Hairer, E.,
-        Rentrop, P.,
+        Nørsett, S.,
-        Numerische Thematic, vol. 33, 1979, pp. 55–68.
+        Wanner, G.,
+        Springer-Verlag, Berlin. 1996.
     \endverbatim
 
     Constants from:
@@ -44,12 +45,12 @@ Description
     \endverbatim
 
 SourceFiles
-    KRR4.C
+    Rosenbrock43.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef KRR4_H
+#ifndef Rosenbrock43_H
-#define KRR4_H
+#define Rosenbrock43_H
 
 #include "ODESolver.H"
 #include "adaptiveSolver.H"
@@ -60,46 +61,48 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class KRR4 Declaration
+                           Class Rosenbrock43 Declaration
 \*---------------------------------------------------------------------------*/
 
-class KRR4
+class Rosenbrock43
 :
     public ODESolver,
     public adaptiveSolver
 {
     // Private data
 
-        mutable scalarField g1_;
+        mutable scalarField k1_;
-        mutable scalarField g2_;
+        mutable scalarField k2_;
-        mutable scalarField g3_;
+        mutable scalarField k3_;
-        mutable scalarField g4_;
+        mutable scalarField k4_;
         mutable scalarField err_;
+        mutable scalarField dydx_;
         mutable scalarField dfdx_;
         mutable scalarSquareMatrix dfdy_;
         mutable scalarSquareMatrix a_;
         mutable labelList pivotIndices_;
 
         static const scalar
-            gamma,
             a21, a31, a32,
-            c21, c31, c32, c41, c42, c43,
+            c21, c31, c32,
+            c41, c42, c43,
             b1, b2, b3, b4,
             e1, e2, e3, e4,
-            c1X, c2X, c3X, c4X,
+            gamma,
-            a2X, a3X;
+            c2, c3,
+            d1, d2, d3, d4;
 
 
 public:
 
     //- Runtime type information
-    TypeName("KRR4");
+    TypeName("Rosenbrock43");
 
 
     // Constructors
 
         //- Construct from ODE
-        KRR4(const ODESystem& ode);
+        Rosenbrock43(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -112,8 +115,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 +124,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/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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,11 @@ Class
 Description
     Foam::SIBS
 
+    Bader, G. and Deuflhard, P.
+    "A Semi-Implicit Mid-Point Rule for
+     Stiff Systems of Ordinary Differential Equations."
+     Numer. Math. 41, 373-398, 1983.
+
 SourceFiles
     SIMPR.C
     polyExtrapolate.C
@@ -67,6 +72,7 @@ class SIBS
         mutable scalarField yTemp_;
         mutable scalarField ySeq_;
         mutable scalarField yErr_;
+        mutable scalarField dydx0_;
         mutable scalarField dfdx_;
         mutable scalarSquareMatrix dfdy_;
 
@@ -110,7 +116,7 @@ public:
     // Constructors
 
         //- Construct from ODE
-        SIBS(const ODESystem& ode);
+        SIBS(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
@@ -120,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/Trapezoid/Trapezoid.C

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "Trapezoid.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(Trapezoid, 0);
+    addToRunTimeSelectionTable(ODESolver, Trapezoid, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::Trapezoid::Trapezoid(const ODESystem& ode, const dictionary& dict)
+:
+    ODESolver(ode, dict),
+    adaptiveSolver(ode, dict),
+    err_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::Trapezoid::solve
+(
+    const ODESystem& ode,
+    const scalar x0,
+    const scalarField& y0,
+    const scalarField& dydx0,
+    const scalar dx,
+    scalarField& y
+) const
+{
+    // Predict the state using 1st-order Trapezoid method
+    forAll(y, i)
+    {
+        y[i] = y0[i] + dx*dydx0[i];
+    }
+
+    // Evaluate the system for the predicted state
+    ode.derivatives(x0 + dx, y, err_);
+
+    // Update the state as the average between the prediction and the correction
+    // and estimate the error from the difference
+    forAll(y, i)
+    {
+        y[i] = y0[i] + 0.5*dx*(dydx0[i] + err_[i]);
+        err_[i] = 0.5*dx*(err_[i] - dydx0[i]);
+    }
+
+    return normalizeError(y0, y, err_);
+}
+
+
+void Foam::Trapezoid::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/Trapezoid/Trapezoid.H

@@ -0,0 +1,104 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::Trapezoid
+
+Description
+    Trapezoidal ODE solver of order (1)2.
+
+SourceFiles
+    Trapezoid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Trapezoid_H
+#define Trapezoid_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Trapezoid Declaration
+\*---------------------------------------------------------------------------*/
+
+class Trapezoid
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        mutable scalarField err_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Trapezoid");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        Trapezoid(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/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,19 @@ 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;
+    if (err > pow(maxScale_/safeScale_, -1.0/alphaInc_))
+    {
+        dxTry =
+            min(max(safeScale_*pow(err, -alphaInc_), minScale_), maxScale_)*dx;
+    }
+    else
+    {
+        dxTry = safeScale_*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/ODE/ODESolvers/rodas43/rodas43.C

@@ -0,0 +1,232 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "rodas43.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(rodas43, 0);
+    addToRunTimeSelectionTable(ODESolver, rodas43, dictionary);
+
+const scalar
+    rodas43::c2 = 0.386,
+    rodas43::c3 = 0.21,
+    rodas43::c4 = 0.63,
+    rodas43::bet2p = 0.0317,
+    rodas43::bet3p = 0.0635,
+    rodas43::bet4p = 0.3438,
+    rodas43::d1 =  0.25,
+    rodas43::d2 = -0.1043,
+    rodas43::d3 =  0.1035,
+    rodas43::d4 = -0.3620000000000023e-01,
+    rodas43::a21 =  0.1544e1,
+    rodas43::a31 =  0.9466785280815826,
+    rodas43::a32 =  0.2557011698983284,
+    rodas43::a41 =  0.3314825187068521e1,
+    rodas43::a42 =  0.2896124015972201e1,
+    rodas43::a43 =  0.9986419139977817,
+    rodas43::a51 =  0.1221224509226641e1,
+    rodas43::a52 =  0.6019134481288629e1,
+    rodas43::a53 =  0.1253708332932087e2,
+    rodas43::a54 = -0.6878860361058950,
+    rodas43::c21 = -0.56688e1,
+    rodas43::c31 = -0.2430093356833875e1,
+    rodas43::c32 = -0.2063599157091915,
+    rodas43::c41 = -0.1073529058151375,
+    rodas43::c42 = -0.9594562251023355e1,
+    rodas43::c43 = -0.2047028614809616e2,
+    rodas43::c51 =  0.7496443313967647e1,
+    rodas43::c52 = -0.1024680431464352e2,
+    rodas43::c53 = -0.3399990352819905e2,
+    rodas43::c54 =  0.1170890893206160e2,
+    rodas43::c61 =  0.8083246795921522e1,
+    rodas43::c62 = -0.7981132988064893e1,
+    rodas43::c63 = -0.3152159432874371e2,
+    rodas43::c64 =  0.1631930543123136e2,
+    rodas43::c65 = -0.6058818238834054e1,
+    rodas43::gamma = 0.25;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict)
+:
+    ODESolver(ode, dict),
+    adaptiveSolver(ode, dict),
+    k1_(n_),
+    k2_(n_),
+    k3_(n_),
+    k4_(n_),
+    k5_(n_),
+    dy_(n_),
+    err_(n_),
+    dydx_(n_),
+    dfdx_(n_),
+    dfdy_(n_, n_),
+    a_(n_, n_),
+    pivotIndices_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::rodas43::solve
+(
+    const ODESystem& ode,
+    const scalar x0,
+    const scalarField& y0,
+    const scalarField& dydx0,
+    const scalar dx,
+    scalarField& y
+) const
+{
+    ode.jacobian(x0, y0, dfdx_, dfdy_);
+
+    for (register label i=0; i<n_; i++)
+    {
+        for (register label j=0; j<n_; j++)
+        {
+            a_[i][j] = -dfdy_[i][j];
+        }
+
+        a_[i][i] += 1.0/(gamma*dx);
+    }
+
+    LUDecompose(a_, pivotIndices_);
+
+    // Calculate k1:
+    forAll(k1_, i)
+    {
+        k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k1_);
+
+    // Calculate k2:
+    forAll(y, i)
+    {
+        y[i] = y0[i] + a21*k1_[i];
+    }
+
+    ode.derivatives(x0 + c2*dx, y, dydx_);
+
+    forAll(k2_, i)
+    {
+        k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k2_);
+
+    // Calculate k3:
+    forAll(y, i)
+    {
+        y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
+    }
+
+    ode.derivatives(x0 + c3*dx, y, dydx_);
+
+    forAll(k3_, i)
+    {
+        k3_[i] = dydx_[i] + dx*d3*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k3_);
+
+    // Calculate k4:
+    forAll(y, i)
+    {
+        y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
+    }
+
+    ode.derivatives(x0 + c4*dx, y, dydx_);
+
+    forAll(k4_, i)
+    {
+        k4_[i] = dydx_[i] + dx*d4*dfdx_[i]
+          + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k4_);
+
+    // Calculate k5:
+    forAll(y, i)
+    {
+        dy_[i] = a51*k1_[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i];
+        y[i] = y0[i] + dy_[i];
+    }
+
+    ode.derivatives(x0 + dx, y, dydx_);
+
+    forAll(k5_, i)
+    {
+        k5_[i] = dydx_[i]
+          + (c51*k1_[i] + c52*k2_[i] + c53*k3_[i] + c54*k4_[i])/dx;
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, k5_);
+
+    // Calculate new state and error
+    forAll(y, i)
+    {
+        dy_[i] += k5_[i];
+        y[i] = y0[i] + dy_[i];
+    }
+
+    ode.derivatives(x0 + dx, y, dydx_);
+
+    forAll(err_, i)
+    {
+        err_[i] = dydx_[i]
+          + (c61*k1_[i] + c62*k2_[i] + c63*k3_[i] + c64*k4_[i] + c65*k5_[i])/dx;
+    }
+
+    LUBacksubstitute(a_, pivotIndices_, err_);
+
+    forAll(y, i)
+    {
+        y[i] = y0[i] + dy_[i] + err_[i];
+    }
+
+    return normalizeError(y0, y, err_);
+}
+
+
+void Foam::rodas43::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/rodas43/rodas43.H

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::rodas43
+
+Description
+    Stiffly-stable embedded Rosenbrock ODE solver of order (3)4.
+
+    References:
+    \verbatim
+        "Solving Ordinary Differential Equations II: Stiff
+         and Differential-Algebraic Problems, second edition",
+        Hairer, E.,
+        Nørsett, S.,
+        Wanner, G.,
+        Springer-Verlag, Berlin. 1996.
+    \endverbatim
+
+SourceFiles
+    rodas43.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rodas43_H
+#define rodas43_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class rodas43 Declaration
+\*---------------------------------------------------------------------------*/
+
+class rodas43
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        mutable scalarField k1_;
+        mutable scalarField k2_;
+        mutable scalarField k3_;
+        mutable scalarField k4_;
+        mutable scalarField k5_;
+        mutable scalarField dy_;
+        mutable scalarField err_;
+        mutable scalarField dydx_;
+        mutable scalarField dfdx_;
+        mutable scalarSquareMatrix dfdy_;
+        mutable scalarSquareMatrix a_;
+        mutable labelList pivotIndices_;
+
+        static const scalar
+            c2, c3, c4,
+            bet2p, bet3p, bet4p,
+            d1, d2, d3, d4,
+            a21, a31, a32,
+            a41, a42, a43,
+            a51, a52, a53, a54,
+            c21, c31, c32,
+            c41, c42, c43,
+            c51, c52, c53, c54,
+            c61, c62, c63, c64, c65,
+            gamma;
+
+public:
+
+    //- Runtime type information
+    TypeName("rodas43");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        rodas43(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/OSspecific/POSIX/signals/sigFpe.C

@@ -45,6 +45,7 @@ License
 #endif
 
 #include <stdint.h>
+#include <limits>
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -53,26 +54,7 @@ struct sigaction Foam::sigFpe::oldAction_;
 
 void Foam::sigFpe::fillSignallingNan(UList<scalar>& lst)
 {
-#ifdef LINUX
+    lst = std::numeric_limits<scalar>::signaling_NaN();
-
-    // initialize to signalling NaN
-#   ifdef WM_SP
-
-    const uint32_t sNAN = 0x7ff7fffflu;
-    uint32_t* dPtr = reinterpret_cast<uint32_t*>(lst.begin());
-
-#   else
-
-    const uint64_t sNAN = 0x7ff7ffffffffffffllu;
-    uint64_t* dPtr = reinterpret_cast<uint64_t*>(lst.begin());
-
-#   endif
-
-    forAll(lst, i)
-    {
-        *dPtr++ = sNAN;
-    }
-#endif
 }
 
 

src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -198,7 +198,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const FixedList<T, Size>& L)
             // Write end delimiter
             os << token::END_BLOCK;
         }
-        else if (Size < 11 && contiguous<T>())
+        else if (Size <= 1 ||(Size < 11 && contiguous<T>()))
         {
             // Write start delimiter
             os << token::BEGIN_LIST;

src/OpenFOAM/containers/Lists/UIndirectList/UIndirectListIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,7 +67,7 @@ Foam::Ostream& Foam::operator<<
             // Write end delimiter
             os << token::END_BLOCK;
         }
-        else if (L.size() < 11 && contiguous<T>())
+        else if(L.size() <= 1 || (L.size() < 11 && contiguous<T>()))
         {
             // Write size and start delimiter
             os << L.size() << token::BEGIN_LIST;

src/OpenFOAM/containers/Lists/UList/UListIO.C

@@ -92,7 +92,7 @@ Foam::Ostream& Foam::operator<<(Foam::Ostream& os, const Foam::UList<T>& L)
             // Write end delimiter
             os << token::END_BLOCK;
         }
-        else if (L.size() == 1 || (L.size() < 11 && contiguous<T>()))
+        else if (L.size() <= 1 || (L.size() < 11 && contiguous<T>()))
         {
             // Write size and start delimiter
             os << L.size() << token::BEGIN_LIST;

src/OpenFOAM/meshes/polyMesh/polyMesh.C

@@ -1328,24 +1328,9 @@ void Foam::polyMesh::findTetFacePt
 {
     const polyMesh& mesh = *this;
 
-    tetFaceI = -1;
+    tetIndices tet(polyMeshTetDecomposition::findTet(mesh, cellI, pt));
-    tetPtI = -1;
+    tetFaceI = tet.face();
-
+    tetPtI = tet.tetPt();
-    List<tetIndices> cellTets =
-        polyMeshTetDecomposition::cellTetIndices(mesh, cellI);
-
-    forAll(cellTets, tetI)
-    {
-        const tetIndices& cellTetIs = cellTets[tetI];
-
-        if (cellTetIs.tet(mesh).inside(pt))
-        {
-            tetFaceI = cellTetIs.face();
-            tetPtI = cellTetIs.tetPt();
-
-            return;
-        }
-    }
 }
 
 

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C

@@ -311,15 +311,14 @@ Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts
         {
             FatalErrorIn
             (
-                "Foam::labelList"
+                "labelList"
-                "Foam::polyMeshTetDecomposition::findFaceBasePts"
+                "polyMeshTetDecomposition::findFaceBasePts"
                 "("
                     "const polyMesh& mesh, "
                     "scalar tol, "
                     "bool report"
                 ")"
-            )
+            )   << "Coupled face base point exchange failure for face "
-                << "Coupled face base point exchange failure for face "
                 << fI
                 << abort(FatalError);
         }
@@ -561,8 +560,8 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
         {
             WarningIn
             (
-                "Foam::List<Foam::tetIndices> "
+                "List<tetIndices> "
-                "Foam::polyMeshTetDecomposition::faceTetIndices"
+                "polyMeshTetDecomposition::faceTetIndices"
                 "("
                     "const polyMesh&, "
                     "label, "
@@ -613,6 +612,76 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::faceTetIndices
 }
 
 
+Foam::tetIndices Foam::polyMeshTetDecomposition::triangleTetIndices
+(
+    const polyMesh& mesh,
+    const label fI,
+    const label cI,
+    const label tetPtI
+)
+{
+    static label nWarnings = 0;
+    static const label maxWarnings = 100;
+
+    const face& f = mesh.faces()[fI];
+    bool own = (mesh.faceOwner()[fI] == cI);
+    label tetBasePtI = mesh.tetBasePtIs()[fI];
+    if (tetBasePtI == -1)
+    {
+        if (nWarnings < maxWarnings)
+        {
+            WarningIn
+            (
+                "tetIndices "
+                "polyMeshTetDecomposition::triangleTetIndices"
+                "("
+                    "const polyMesh&, "
+                    "label, "
+                    "label, "
+                    "label"
+                ")"
+            )   << "No base point for face " << fI << ", " << f
+                << ", produces a valid tet decomposition."
+                << endl;
+            nWarnings++;
+        }
+        if (nWarnings == maxWarnings)
+        {
+            Warning<< "Suppressing any further warnings." << endl;
+            nWarnings++;
+        }
+
+        tetBasePtI = 0;
+    }
+
+    tetIndices faceTetIs;
+
+    label facePtI = (tetPtI + tetBasePtI) % f.size();
+    label otherFacePtI = f.fcIndex(facePtI);
+
+    faceTetIs.cell() = cI;
+
+    faceTetIs.face() = fI;
+
+    faceTetIs.faceBasePt() = tetBasePtI;
+
+    if (own)
+    {
+        faceTetIs.facePtA() = facePtI;
+        faceTetIs.facePtB() = otherFacePtI;
+    }
+    else
+    {
+        faceTetIs.facePtA() = otherFacePtI;
+        faceTetIs.facePtB() = facePtI;
+    }
+
+    faceTetIs.tetPt() = tetPtI;
+
+    return faceTetIs;
+}
+
+
 Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
 (
     const polyMesh& mesh,
@@ -644,4 +713,56 @@ Foam::List<Foam::tetIndices> Foam::polyMeshTetDecomposition::cellTetIndices
 }
 
 
+Foam::tetIndices Foam::polyMeshTetDecomposition::findTet
+(
+    const polyMesh& mesh,
+    label cI,
+    const point& pt
+)
+{
+    const faceList& pFaces = mesh.faces();
+    const cellList& pCells = mesh.cells();
+
+    const cell& thisCell = pCells[cI];
+
+    tetIndices tetContainingPt;
+
+
+    forAll(thisCell, cFI)
+    {
+        label fI = thisCell[cFI];
+        const face& f = pFaces[fI];
+
+        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
+        {
+            // Get tetIndices of face triangle
+            tetIndices faceTetIs
+            (
+                triangleTetIndices
+                (
+                    mesh,
+                    fI,
+                    cI,
+                    tetPtI
+                )
+            );
+
+            // Check if inside
+            if (faceTetIs.tet(mesh).inside(pt))
+            {
+                tetContainingPt = faceTetIs;
+                break;
+            }
+        }
+
+        if (tetContainingPt.cell() != -1)
+        {
+            break;
+        }
+    }
+
+    return tetContainingPt;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H

@@ -130,6 +130,15 @@ public:
             label cI
         );
 
+        //- Return the tet decomposition of the given triangle of the given face
+        static tetIndices triangleTetIndices
+        (
+            const polyMesh& mesh,
+            label fI,
+            label cI,
+            const label tetPtI      // offset in face
+        );
+
         //- Return the tet decomposition of the given cell, see
         //  findFacePt for the meaning of the indices
         static List<tetIndices> cellTetIndices
@@ -137,6 +146,14 @@ public:
             const polyMesh& mesh,
             label cI
         );
+
+        //- Find the tet decomposition of the cell containing the given point
+        static tetIndices findTet
+        (
+            const polyMesh& mesh,
+            label cI,
+            const point& pt
+        );
 };
 
 

src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C

@@ -1784,7 +1784,7 @@ void Foam::faceCoupleInfo::subDivisionMatch
                     << " points:" << masterPoints[masterEdge[0]]
                     << ' ' << masterPoints[masterEdge[1]]
                     << " corresponding cut edge: (" << cutPoint0
-                    << ' ' << cutPoint1
+                    << ") " << cutPoint1
                     << abort(FatalError);
             }
 

src/dynamicMesh/polyMeshAdder/polyMeshAdder.C

@@ -38,51 +38,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Append all mapped elements of a list to a DynamicList
-void Foam::polyMeshAdder::append
-(
-    const labelList& map,
-    const labelList& lst,
-    DynamicList<label>& dynLst
-)
-{
-    dynLst.setCapacity(dynLst.size() + lst.size());
-
-    forAll(lst, i)
-    {
-        const label newElem = map[lst[i]];
-
-        if (newElem != -1)
-        {
-            dynLst.append(newElem);
-        }
-    }
-}
-
-
-// Append all mapped elements of a list to a DynamicList
-void Foam::polyMeshAdder::append
-(
-    const labelList& map,
-    const labelList& lst,
-    const SortableList<label>& sortedLst,
-    DynamicList<label>& dynLst
-)
-{
-    dynLst.setCapacity(dynLst.size() + lst.size());
-
-    forAll(lst, i)
-    {
-        const label newElem = map[lst[i]];
-
-        if (newElem != -1 && findSortedIndex(sortedLst, newElem) == -1)
-        {
-            dynLst.append(newElem);
-        }
-    }
-}
-
-
 // Get index of patch in new set of patchnames/types
 Foam::label Foam::polyMeshAdder::patchIndex
 (
@@ -913,6 +868,7 @@ void Foam::polyMeshAdder::mergePrimitives
 
 void Foam::polyMeshAdder::mergePointZones
 (
+    const label nAllPoints,
     const pointZoneMesh& pz0,
     const pointZoneMesh& pz1,
     const labelList& from0ToAllPoints,
@@ -934,60 +890,141 @@ void Foam::polyMeshAdder::mergePointZones
     }
     zoneNames.shrink();
 
-    // Point labels per merged zone
-    pzPoints.setSize(zoneNames.size());
 
+
+    // Zone(s) per point. Two levels: if only one zone
+    // stored in pointToZone. Any extra stored in additionalPointToZones.
+    // This is so we only allocate labelLists per point if absolutely
+    // necesary.
+    labelList pointToZone(nAllPoints, -1);
+    labelListList addPointToZones(nAllPoints);
+
+    // mesh0 zones kept
     forAll(pz0, zoneI)
     {
-        append(from0ToAllPoints, pz0[zoneI], pzPoints[zoneI]);
+        const pointZone& pz = pz0[zoneI];
-    }
 
-    // Get sorted zone contents for duplicate element recognition
+        forAll(pz, i)
-    PtrList<SortableList<label> > pzPointsSorted(pzPoints.size());
-    forAll(pzPoints, zoneI)
         {
-        pzPointsSorted.set
+            label point0 = pz[i];
-        (
+            label allPointI = from0ToAllPoints[point0];
-            zoneI,
+
-            new SortableList<label>(pzPoints[zoneI])
+            if (pointToZone[allPointI] == -1)
-        );
+            {
+                pointToZone[allPointI] = zoneI;
+            }
+            else if (pointToZone[allPointI] != zoneI)
+            {
+                labelList& pZones = addPointToZones[allPointI];
+                if (findIndex(pZones, zoneI) == -1)
+                {
+                    pZones.append(zoneI);
+                }
+            }
+        }
     }
 
-    // Now we have full addressing for points so do the pointZones of mesh1.
+    // mesh1 zones renumbered
     forAll(pz1, zoneI)
     {
-        // Relabel all points of zone and add to correct pzPoints.
+        const pointZone& pz = pz1[zoneI];
         const label allZoneI = from1ToAll[zoneI];
 
-        append
+        forAll(pz, i)
-        (
+        {
-            from1ToAllPoints,
+            label point1 = pz[i];
-            pz1[zoneI],
+            label allPointI = from1ToAllPoints[point1];
-            pzPointsSorted[allZoneI],
+
-            pzPoints[allZoneI]
+            if (pointToZone[allPointI] == -1)
-        );
+            {
+                pointToZone[allPointI] = allZoneI;
+            }
+            else if (pointToZone[allPointI] != allZoneI)
+            {
+                labelList& pZones = addPointToZones[allPointI];
+                if (findIndex(pZones, allZoneI) == -1)
+                {
+                    pZones.append(allZoneI);
+                }
+            }
+        }
     }
 
+    // Extract back into zones
+
+    // 1. Count
+    labelList nPoints(zoneNames.size(), 0);
+    forAll(pointToZone, allPointI)
+    {
+        label zoneI = pointToZone[allPointI];
+        if (zoneI != -1)
+        {
+            nPoints[zoneI]++;
+        }
+    }
+    forAll(addPointToZones, allPointI)
+    {
+        const labelList& pZones = addPointToZones[allPointI];
+        forAll(pZones, i)
+        {
+            nPoints[pZones[i]]++;
+        }
+    }
+
+    // 2. Fill
+    pzPoints.setSize(zoneNames.size());
+    forAll(pzPoints, zoneI)
+    {
+        pzPoints[zoneI].setCapacity(nPoints[zoneI]);
+    }
+    forAll(pointToZone, allPointI)
+    {
+        label zoneI = pointToZone[allPointI];
+        if (zoneI != -1)
+        {
+            pzPoints[zoneI].append(allPointI);
+        }
+    }
+    forAll(addPointToZones, allPointI)
+    {
+        const labelList& pZones = addPointToZones[allPointI];
+        forAll(pZones, i)
+        {
+            pzPoints[pZones[i]].append(allPointI);
+        }
+    }
     forAll(pzPoints, i)
     {
         pzPoints[i].shrink();
+        stableSort(pzPoints[i]);
     }
 }
 
 
 void Foam::polyMeshAdder::mergeFaceZones
 (
-    const faceZoneMesh& fz0,
+    const labelList& allOwner,
-    const faceZoneMesh& fz1,
+
+    const polyMesh& mesh0,
+    const polyMesh& mesh1,
+
     const labelList& from0ToAllFaces,
     const labelList& from1ToAllFaces,
 
+    const labelList& from1ToAllCells,
+
     DynamicList<word>& zoneNames,
     labelList& from1ToAll,
     List<DynamicList<label> >& fzFaces,
     List<DynamicList<bool> >& fzFlips
 )
 {
+    const faceZoneMesh& fz0 = mesh0.faceZones();
+    const labelList& owner0 = mesh0.faceOwner();
+    const faceZoneMesh& fz1 = mesh1.faceZones();
+    const labelList& owner1 = mesh1.faceOwner();
+
+
     zoneNames.setCapacity(fz0.size() + fz1.size());
     zoneNames.append(fz0.names());
 
@@ -1000,72 +1037,146 @@ void Foam::polyMeshAdder::mergeFaceZones
     zoneNames.shrink();
 
 
-    // Create temporary lists for faceZones.
+    // Zone(s) per face
-    fzFaces.setSize(zoneNames.size());
+    labelList faceToZone(allOwner.size(), -1);
-    fzFlips.setSize(zoneNames.size());
+    labelListList addFaceToZones(allOwner.size());
+    boolList faceToFlip(allOwner.size(), false);
+    boolListList addFaceToFlips(allOwner.size());
+
+    // mesh0 zones kept
     forAll(fz0, zoneI)
     {
-        DynamicList<label>& newZone = fzFaces[zoneI];
-        DynamicList<bool>& newFlip = fzFlips[zoneI];
-
-        newZone.setCapacity(fz0[zoneI].size());
-        newFlip.setCapacity(newZone.size());
-
         const labelList& addressing = fz0[zoneI];
         const boolList& flipMap = fz0[zoneI].flipMap();
 
         forAll(addressing, i)
         {
-            label faceI = addressing[i];
+            label face0 = addressing[i];
+            bool flip0 = flipMap[i];
 
-            if (from0ToAllFaces[faceI] != -1)
+            label allFaceI = from0ToAllFaces[face0];
+            if (allFaceI != -1)
             {
-                newZone.append(from0ToAllFaces[faceI]);
+                // Check if orientation same
-                newFlip.append(flipMap[i]);
+                label allCell0 = owner0[face0];
-            }
+                if (allOwner[allFaceI] != allCell0)
-        }
-    }
-
-    // Get sorted zone contents for duplicate element recognition
-    PtrList<SortableList<label> > fzFacesSorted(fzFaces.size());
-    forAll(fzFaces, zoneI)
                 {
-        fzFacesSorted.set
+                    flip0 = !flip0;
-        (
-            zoneI,
-            new SortableList<label>(fzFaces[zoneI])
-        );
                 }
 
-    // Now we have full addressing for faces so do the faceZones of mesh1.
+                if (faceToZone[allFaceI] == -1)
+                {
+                    faceToZone[allFaceI] = zoneI;
+                    faceToFlip[allFaceI] = flip0;
+                }
+                else if (faceToZone[allFaceI] != zoneI)
+                {
+                    labelList& fZones = addFaceToZones[allFaceI];
+                    boolList& flipZones = addFaceToFlips[allFaceI];
+
+                    if (findIndex(fZones, zoneI) == -1)
+                    {
+                        fZones.append(zoneI);
+                        flipZones.append(flip0);
+                    }
+                }
+            }
+        }
+    }
+
+    // mesh1 zones renumbered
     forAll(fz1, zoneI)
     {
-        label allZoneI = from1ToAll[zoneI];
-
-        DynamicList<label>& newZone = fzFaces[allZoneI];
-        const SortableList<label>& newZoneSorted = fzFacesSorted[allZoneI];
-        DynamicList<bool>& newFlip = fzFlips[allZoneI];
-
-        newZone.setCapacity(newZone.size() + fz1[zoneI].size());
-        newFlip.setCapacity(newZone.size());
-
         const labelList& addressing = fz1[zoneI];
         const boolList& flipMap = fz1[zoneI].flipMap();
 
+        const label allZoneI = from1ToAll[zoneI];
+
         forAll(addressing, i)
         {
-            label faceI = addressing[i];
+            label face1 = addressing[i];
-            label allFaceI = from1ToAllFaces[faceI];
+            bool flip1 = flipMap[i];
 
-            if
+            label allFaceI = from1ToAllFaces[face1];
-            (
+            if (allFaceI != -1)
-                allFaceI != -1
-             && findSortedIndex(newZoneSorted, allFaceI) == -1
-            )
             {
-                newZone.append(allFaceI);
+                // Check if orientation same
-                newFlip.append(flipMap[i]);
+                label allCell1 = from1ToAllCells[owner1[face1]];
+                if (allOwner[allFaceI] != allCell1)
+                {
+                    flip1 = !flip1;
                 }
+
+                if (faceToZone[allFaceI] == -1)
+                {
+                    faceToZone[allFaceI] = allZoneI;
+                    faceToFlip[allFaceI] = flip1;
+                }
+                else if (faceToZone[allFaceI] != allZoneI)
+                {
+                    labelList& fZones = addFaceToZones[allFaceI];
+                    boolList& flipZones = addFaceToFlips[allFaceI];
+
+                    if (findIndex(fZones, allZoneI) == -1)
+                    {
+                        fZones.append(allZoneI);
+                        flipZones.append(flip1);
+                    }
+                }
+            }
+        }
+    }
+
+
+    // Extract back into zones
+
+    // 1. Count
+    labelList nFaces(zoneNames.size(), 0);
+    forAll(faceToZone, allFaceI)
+    {
+        label zoneI = faceToZone[allFaceI];
+        if (zoneI != -1)
+        {
+            nFaces[zoneI]++;
+        }
+    }
+    forAll(addFaceToZones, allFaceI)
+    {
+        const labelList& fZones = addFaceToZones[allFaceI];
+        forAll(fZones, i)
+        {
+            nFaces[fZones[i]]++;
+        }
+    }
+
+    // 2. Fill
+    fzFaces.setSize(zoneNames.size());
+    fzFlips.setSize(zoneNames.size());
+    forAll(fzFaces, zoneI)
+    {
+        fzFaces[zoneI].setCapacity(nFaces[zoneI]);
+        fzFlips[zoneI].setCapacity(nFaces[zoneI]);
+    }
+    forAll(faceToZone, allFaceI)
+    {
+        label zoneI = faceToZone[allFaceI];
+        bool flip = faceToFlip[allFaceI];
+        if (zoneI != -1)
+        {
+            fzFaces[zoneI].append(allFaceI);
+            fzFlips[zoneI].append(flip);
+        }
+    }
+    forAll(addFaceToZones, allFaceI)
+    {
+        const labelList& fZones = addFaceToZones[allFaceI];
+        const boolList& flipZones = addFaceToFlips[allFaceI];
+
+        forAll(fZones, i)
+        {
+            label zoneI = fZones[i];
+            fzFaces[zoneI].append(allFaceI);
+            fzFlips[zoneI].append(flipZones[i]);
         }
     }
 
@@ -1073,12 +1184,19 @@ void Foam::polyMeshAdder::mergeFaceZones
     {
         fzFaces[i].shrink();
         fzFlips[i].shrink();
+
+        labelList order;
+        sortedOrder(fzFaces[i], order);
+        fzFaces[i] = UIndirectList<label>(fzFaces[i], order)();
+        fzFlips[i] = UIndirectList<bool>(fzFlips[i], order)();
     }
 }
 
 
 void Foam::polyMeshAdder::mergeCellZones
 (
+    const label nAllCells,
+
     const cellZoneMesh& cz0,
     const cellZoneMesh& cz1,
     const labelList& from1ToAllCells,
@@ -1099,32 +1217,118 @@ void Foam::polyMeshAdder::mergeCellZones
     zoneNames.shrink();
 
 
-    // Create temporary lists for cellZones.
+    // Zone(s) per cell. Two levels: if only one zone
-    czCells.setSize(zoneNames.size());
+    // stored in cellToZone. Any extra stored in additionalCellToZones.
+    // This is so we only allocate labelLists per cell if absolutely
+    // necesary.
+    labelList cellToZone(nAllCells, -1);
+    labelListList addCellToZones(nAllCells);
+
+    // mesh0 zones kept
     forAll(cz0, zoneI)
     {
-        // Insert mesh0 cells
+        const cellZone& cz = cz0[zoneI];
-        czCells[zoneI].append(cz0[zoneI]);
+        forAll(cz, i)
+        {
+            label cell0 = cz[i];
+
+            if (cellToZone[cell0] == -1)
+            {
+                cellToZone[cell0] = zoneI;
+            }
+            else if (cellToZone[cell0] != zoneI)
+            {
+                labelList& cZones = addCellToZones[cell0];
+                if (findIndex(cZones, zoneI) == -1)
+                {
+                    cZones.append(zoneI);
+                }
+            }
+        }
     }
 
-
+    // mesh1 zones renumbered
-    // Cell mapping is trivial.
     forAll(cz1, zoneI)
     {
+        const cellZone& cz = cz1[zoneI];
         const label allZoneI = from1ToAll[zoneI];
+        forAll(cz, i)
+        {
+            label cell1 = cz[i];
+            label allCellI = from1ToAllCells[cell1];
 
-        append(from1ToAllCells, cz1[zoneI], czCells[allZoneI]);
+            if (cellToZone[allCellI] == -1)
+            {
+                cellToZone[allCellI] = allZoneI;
+            }
+            else if (cellToZone[allCellI] != allZoneI)
+            {
+                labelList& cZones = addCellToZones[allCellI];
+                if (findIndex(cZones, allZoneI) == -1)
+                {
+                    cZones.append(allZoneI);
+                }
+            }
+        }
     }
 
+    // Extract back into zones
+
+    // 1. Count
+    labelList nCells(zoneNames.size(), 0);
+    forAll(cellToZone, allCellI)
+    {
+        label zoneI = cellToZone[allCellI];
+        if (zoneI != -1)
+        {
+            nCells[zoneI]++;
+        }
+    }
+    forAll(addCellToZones, allCellI)
+    {
+        const labelList& cZones = addCellToZones[allCellI];
+        forAll(cZones, i)
+        {
+            nCells[cZones[i]]++;
+        }
+    }
+
+    // 2. Fill
+    czCells.setSize(zoneNames.size());
+    forAll(czCells, zoneI)
+    {
+        czCells[zoneI].setCapacity(nCells[zoneI]);
+    }
+    forAll(cellToZone, allCellI)
+    {
+        label zoneI = cellToZone[allCellI];
+        if (zoneI != -1)
+        {
+            czCells[zoneI].append(allCellI);
+        }
+    }
+    forAll(addCellToZones, allCellI)
+    {
+        const labelList& cZones = addCellToZones[allCellI];
+        forAll(cZones, i)
+        {
+            czCells[cZones[i]].append(allCellI);
+        }
+    }
     forAll(czCells, i)
     {
         czCells[i].shrink();
+        stableSort(czCells[i]);
     }
 }
 
 
 void Foam::polyMeshAdder::mergeZones
 (
+    const label nAllPoints,
+    const labelList& allOwner,
+    const label nAllCells,
+
     const polyMesh& mesh0,
     const polyMesh& mesh1,
     const labelList& from0ToAllPoints,
@@ -1148,6 +1352,7 @@ void Foam::polyMeshAdder::mergeZones
     labelList from1ToAllPZones;
     mergePointZones
     (
+        nAllPoints,
         mesh0.pointZones(),
         mesh1.pointZones(),
         from0ToAllPoints,
@@ -1161,10 +1366,12 @@ void Foam::polyMeshAdder::mergeZones
     labelList from1ToAllFZones;
     mergeFaceZones
     (
-        mesh0.faceZones(),
+        allOwner,
-        mesh1.faceZones(),
+        mesh0,
+        mesh1,
         from0ToAllFaces,
         from1ToAllFaces,
+        from1ToAllCells,
 
         faceZoneNames,
         from1ToAllFZones,
@@ -1175,6 +1382,7 @@ void Foam::polyMeshAdder::mergeZones
     labelList from1ToAllCZones;
     mergeCellZones
     (
+        nAllCells,
         mesh0.cellZones(),
         mesh1.cellZones(),
         from1ToAllCells,
@@ -1353,6 +1561,10 @@ Foam::autoPtr<Foam::polyMesh> Foam::polyMeshAdder::add
 
     mergeZones
     (
+        allPoints.size(),
+        allOwner,
+        nCells,
+
         mesh0,
         mesh1,
 
@@ -1566,6 +1778,10 @@ Foam::autoPtr<Foam::mapAddedPolyMesh> Foam::polyMeshAdder::add
 
     mergeZones
     (
+        allPoints.size(),
+        allOwner,
+        nCells,
+
         mesh0,
         mesh1,
 

src/dynamicMesh/polyMeshAdder/polyMeshAdder.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,29 +61,9 @@ class polyMeshAdder
 {
 
 private:
-    // Private data
-
 
     // Private Member Functions
 
-        //- Append all mapped elements of a list to a DynamicList
-        static void append
-        (
-            const labelList& map,
-            const labelList& lst,
-            DynamicList<label>&
-        );
-
-        //- Append all mapped elements of a list to a DynamicList that are
-        //  not already present in the sorted list.
-        static void append
-        (
-            const labelList& map,
-            const labelList& lst,
-            const SortableList<label>& sortedLst,
-            DynamicList<label>&
-        );
-
         //- Index of patch in allPatches. Add if nonexisting.
         static label patchIndex
         (
@@ -180,6 +160,8 @@ private:
         //- Merge point zones
         static void mergePointZones
         (
+            const label nAllPoints,
+
             const pointZoneMesh& pz0,
             const pointZoneMesh& pz1,
             const labelList& from0ToAllPoints,
@@ -193,10 +175,13 @@ private:
         //- Merge face zones
         static void mergeFaceZones
         (
-            const faceZoneMesh& fz0,
+            const labelList& allOwner,
-            const faceZoneMesh& fz1,
+
+            const polyMesh& mesh0,
+            const polyMesh& mesh1,
             const labelList& from0ToAllFaces,
             const labelList& from1ToAllFaces,
+            const labelList& from1ToAllCells,
 
             DynamicList<word>& zoneNames,
             labelList& from1ToAll,
@@ -207,6 +192,8 @@ private:
         //- Merge cell zones
         static void mergeCellZones
         (
+            const label nAllCells,
+
             const cellZoneMesh& cz0,
             const cellZoneMesh& cz1,
             const labelList& from1ToAllCells,
@@ -219,6 +206,10 @@ private:
         //- Merge point/face/cell zone information
         static void mergeZones
         (
+            const label nAllPoints,
+            const labelList& allOwner,
+            const label nAllCells,
+
             const polyMesh& mesh0,
             const polyMesh& mesh1,
             const labelList& from0ToAllPoints,

src/finiteVolume/finiteVolume/ddtSchemes/CoEulerDdtScheme/CoEulerDdtScheme.C

@@ -184,7 +184,8 @@ CoEulerDdtScheme<Type>::fvcDdt
         );
 
         tdtdt().internalField() =
-            rDeltaT.internalField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
+            rDeltaT.internalField()*dt.value()
+           *(1.0 - mesh().Vsc0()/mesh().Vsc());
 
         return tdtdt;
     }
@@ -237,7 +238,7 @@ CoEulerDdtScheme<Type>::fvcDdt
                 rDeltaT.internalField()*
                 (
                     vf.internalField()
-                  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                  - vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.boundaryField()*
                 (
@@ -289,7 +290,7 @@ CoEulerDdtScheme<Type>::fvcDdt
                 rDeltaT.internalField()*rho.value()*
                 (
                     vf.internalField()
-                  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                  - vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.boundaryField()*rho.value()*
                 (
@@ -342,7 +343,7 @@ CoEulerDdtScheme<Type>::fvcDdt
                 (
                     rho.internalField()*vf.internalField()
                   - rho.oldTime().internalField()
-                   *vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                   *vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.boundaryField()*
                 (
@@ -387,15 +388,15 @@ CoEulerDdtScheme<Type>::fvmDdt
 
     scalarField rDeltaT(CorDeltaT()().internalField());
 
-    fvm.diag() = rDeltaT*mesh().V();
+    fvm.diag() = rDeltaT*mesh().Vsc();
 
     if (mesh().moving())
     {
-        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
+        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
-        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
+        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;
@@ -422,17 +423,17 @@ CoEulerDdtScheme<Type>::fvmDdt
 
     scalarField rDeltaT(CorDeltaT()().internalField());
 
-    fvm.diag() = rDeltaT*rho.value()*mesh().V();
+    fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
 
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
-            *rho.value()*vf.oldTime().internalField()*mesh().V0();
+            *rho.value()*vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
         fvm.source() = rDeltaT
-            *rho.value()*vf.oldTime().internalField()*mesh().V();
+            *rho.value()*vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;
@@ -459,19 +460,19 @@ CoEulerDdtScheme<Type>::fvmDdt
 
     scalarField rDeltaT(CorDeltaT()().internalField());
 
-    fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
+    fvm.diag() = rDeltaT*rho.internalField()*mesh().Vsc();
 
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
             *rho.oldTime().internalField()
-            *vf.oldTime().internalField()*mesh().V0();
+            *vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
         fvm.source() = rDeltaT
             *rho.oldTime().internalField()
-            *vf.oldTime().internalField()*mesh().V();
+            *vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;

src/finiteVolume/finiteVolume/ddtSchemes/EulerDdtScheme/EulerDdtScheme.C

@@ -74,7 +74,7 @@ EulerDdtScheme<Type>::fvcDdt
         );
 
         tdtdt().internalField() =
-            rDeltaT.value()*dt.value()*(1.0 - mesh().V0()/mesh().V());
+            rDeltaT.value()*dt.value()*(1.0 - mesh().Vsc0()/mesh().Vsc());
 
         return tdtdt;
     }
@@ -127,7 +127,7 @@ EulerDdtScheme<Type>::fvcDdt
                 rDeltaT.value()*
                 (
                     vf.internalField()
-                  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                  - vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.value()*
                 (
@@ -179,7 +179,7 @@ EulerDdtScheme<Type>::fvcDdt
                 rDeltaT.value()*rho.value()*
                 (
                     vf.internalField()
-                  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                  - vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.value()*rho.value()*
                 (
@@ -232,7 +232,7 @@ EulerDdtScheme<Type>::fvcDdt
                 (
                     rho.internalField()*vf.internalField()
                   - rho.oldTime().internalField()
-                   *vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                   *vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.value()*
                 (
@@ -277,15 +277,15 @@ EulerDdtScheme<Type>::fvmDdt
 
     scalar rDeltaT = 1.0/mesh().time().deltaTValue();
 
-    fvm.diag() = rDeltaT*mesh().V();
+    fvm.diag() = rDeltaT*mesh().Vsc();
 
     if (mesh().moving())
     {
-        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
+        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
-        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
+        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;
@@ -312,17 +312,17 @@ EulerDdtScheme<Type>::fvmDdt
 
     scalar rDeltaT = 1.0/mesh().time().deltaTValue();
 
-    fvm.diag() = rDeltaT*rho.value()*mesh().V();
+    fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
 
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
-            *rho.value()*vf.oldTime().internalField()*mesh().V0();
+            *rho.value()*vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
         fvm.source() = rDeltaT
-            *rho.value()*vf.oldTime().internalField()*mesh().V();
+            *rho.value()*vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;
@@ -349,19 +349,19 @@ EulerDdtScheme<Type>::fvmDdt
 
     scalar rDeltaT = 1.0/mesh().time().deltaTValue();
 
-    fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
+    fvm.diag() = rDeltaT*rho.internalField()*mesh().Vsc();
 
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
             *rho.oldTime().internalField()
-            *vf.oldTime().internalField()*mesh().V0();
+            *vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
         fvm.source() = rDeltaT
             *rho.oldTime().internalField()
-            *vf.oldTime().internalField()*mesh().V();
+            *vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;

src/finiteVolume/finiteVolume/ddtSchemes/localEulerDdtScheme/localEulerDdtScheme.C

@@ -81,7 +81,8 @@ localEulerDdtScheme<Type>::fvcDdt
         );
 
         tdtdt().internalField() =
-            rDeltaT.internalField()*dt.value()*(1.0 - mesh().V0()/mesh().V());
+            rDeltaT.internalField()*dt.value()
+           *(1.0 - mesh().Vsc0()/mesh().Vsc());
 
         return tdtdt;
     }
@@ -134,7 +135,7 @@ localEulerDdtScheme<Type>::fvcDdt
                 rDeltaT.internalField()*
                 (
                     vf.internalField()
-                  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                  - vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.boundaryField()*
                 (
@@ -186,7 +187,7 @@ localEulerDdtScheme<Type>::fvcDdt
                 rDeltaT.internalField()*rho.value()*
                 (
                     vf.internalField()
-                  - vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                  - vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.boundaryField()*rho.value()*
                 (
@@ -239,7 +240,7 @@ localEulerDdtScheme<Type>::fvcDdt
                 (
                     rho.internalField()*vf.internalField()
                   - rho.oldTime().internalField()
-                   *vf.oldTime().internalField()*mesh().V0()/mesh().V()
+                   *vf.oldTime().internalField()*mesh().Vsc0()/mesh().Vsc()
                 ),
                 rDeltaT.boundaryField()*
                 (
@@ -284,15 +285,15 @@ localEulerDdtScheme<Type>::fvmDdt
 
     const scalarField& rDeltaT = localRDeltaT().internalField();
 
-    fvm.diag() = rDeltaT*mesh().V();
+    fvm.diag() = rDeltaT*mesh().Vsc();
 
     if (mesh().moving())
     {
-        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V0();
+        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
-        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().V();
+        fvm.source() = rDeltaT*vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;
@@ -319,17 +320,17 @@ localEulerDdtScheme<Type>::fvmDdt
 
     const scalarField& rDeltaT = localRDeltaT().internalField();
 
-    fvm.diag() = rDeltaT*rho.value()*mesh().V();
+    fvm.diag() = rDeltaT*rho.value()*mesh().Vsc();
 
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
-            *rho.value()*vf.oldTime().internalField()*mesh().V0();
+            *rho.value()*vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
         fvm.source() = rDeltaT
-            *rho.value()*vf.oldTime().internalField()*mesh().V();
+            *rho.value()*vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;
@@ -356,19 +357,19 @@ localEulerDdtScheme<Type>::fvmDdt
 
     const scalarField& rDeltaT = localRDeltaT().internalField();
 
-    fvm.diag() = rDeltaT*rho.internalField()*mesh().V();
+    fvm.diag() = rDeltaT*rho.internalField()*mesh().Vsc();
 
     if (mesh().moving())
     {
         fvm.source() = rDeltaT
             *rho.oldTime().internalField()
-            *vf.oldTime().internalField()*mesh().V0();
+            *vf.oldTime().internalField()*mesh().Vsc0();
     }
     else
     {
         fvm.source() = rDeltaT
             *rho.oldTime().internalField()
-            *vf.oldTime().internalField()*mesh().V();
+            *vf.oldTime().internalField()*mesh().Vsc();
     }
 
     return tfvm;

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C

@@ -461,7 +461,8 @@ void Foam::meshRefinement::markFeatureCellLevel
     // Track all particles to their end position (= starting feature point)
     // Note that the particle might have started on a different processor
     // so this will transfer across nicely until we can start tracking proper.
-    startPointCloud.move(td, GREAT);
+    scalar maxTrackLen = 2.0*mesh_.bounds().mag();
+    startPointCloud.move(td, maxTrackLen);
 
 
     // Reset level
@@ -531,7 +532,7 @@ void Foam::meshRefinement::markFeatureCellLevel
     while (true)
     {
         // Track all particles to their end position.
-        cloud.move(td, GREAT);
+        cloud.move(td, maxTrackLen);
 
 
         label nParticles = 0;

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

@@ -0,0 +1,226 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "findCellParticle.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::findCellParticle::findCellParticle
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const point& end,
+    const label data
+)
+:
+    particle(mesh, position, cellI, tetFaceI, tetPtI),
+    end_(end),
+    data_(data)
+{}
+
+
+Foam::findCellParticle::findCellParticle
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    particle(mesh, is, readFields)
+{
+    if (readFields)
+    {
+        if (is.format() == IOstream::ASCII)
+        {
+            is >> end_;
+            data_ = readLabel(is);
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&end_),
+                sizeof(end_) + sizeof(data_)
+            );
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "findCellParticle::findCellParticle"
+        "(const Cloud<findCellParticle>&, Istream&, bool)"
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::findCellParticle::move
+(
+    trackingData& td,
+    const scalar maxTrackLen
+)
+{
+    td.switchProcessor = false;
+    td.keepParticle = true;
+
+    scalar tEnd = (1.0 - stepFraction())*maxTrackLen;
+    scalar dtMax = tEnd;
+
+    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
+    {
+        // set the lagrangian time-step
+        scalar dt = min(dtMax, tEnd);
+
+        dt *= trackToFace(end_, td);
+
+        tEnd -= dt;
+        stepFraction() = 1.0 - tEnd/maxTrackLen;
+    }
+
+    if (tEnd < SMALL || !td.keepParticle)
+    {
+        // Hit endpoint or patch. If patch hit could do fancy stuff but just
+        // to use the patch point is good enough for now.
+        td.cellToData()[cell()].append(data());
+        td.cellToEnd()[cell()].append(position());
+    }
+
+    return td.keepParticle;
+}
+
+
+bool Foam::findCellParticle::hitPatch
+(
+    const polyPatch&,
+    trackingData& td,
+    const label patchI,
+    const scalar trackFraction,
+    const tetIndices& tetIs
+)
+{
+    return false;
+}
+
+
+void Foam::findCellParticle::hitWedgePatch
+(
+    const wedgePolyPatch&,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::findCellParticle::hitSymmetryPatch
+(
+    const symmetryPolyPatch&,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::findCellParticle::hitCyclicPatch
+(
+    const cyclicPolyPatch&,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::findCellParticle::hitProcessorPatch
+(
+    const processorPolyPatch&,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.switchProcessor = true;
+}
+
+
+void Foam::findCellParticle::hitWallPatch
+(
+    const wallPolyPatch& wpp,
+    trackingData& td,
+    const tetIndices&
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+void Foam::findCellParticle::hitPatch
+(
+    const polyPatch& wpp,
+    trackingData& td
+)
+{
+    // Remove particle
+    td.keepParticle = false;
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const findCellParticle& p)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << static_cast<const particle&>(p)
+            << token::SPACE << p.end_
+            << token::SPACE << p.data_;
+    }
+    else
+    {
+        os  << static_cast<const particle&>(p);
+        os.write
+        (
+            reinterpret_cast<const char*>(&p.end_),
+            sizeof(p.end_) + sizeof(p.data_)
+        );
+    }
+
+    // Check state of Ostream
+    os.check("Ostream& operator<<(Ostream&, const findCellParticle&)");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/postProcessing/functionObjects/field/nearWallFields/findCellParticle.H

@@ -0,0 +1,259 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::findCellParticle
+
+Description
+    Particle class that finds cells by tracking
+
+SourceFiles
+    findCellParticle.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef findCellParticle_H
+#define findCellParticle_H
+
+#include "particle.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class findCellParticleCloud;
+
+/*---------------------------------------------------------------------------*\
+                     Class findCellParticle Declaration
+\*---------------------------------------------------------------------------*/
+
+class findCellParticle
+:
+    public particle
+{
+    // Private data
+
+        //- end point to track to
+        point end_;
+
+        //- passive data
+        label data_;
+
+public:
+
+    friend class Cloud<findCellParticle>;
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackingData
+    :
+        public particle::TrackingData<Cloud<findCellParticle> >
+    {
+        labelListList& cellToData_;
+        List<List<point> >& cellToEnd_;
+
+    public:
+
+
+        // Constructors
+
+            trackingData
+            (
+                Cloud<findCellParticle>& cloud,
+                labelListList& cellToData,
+                List<List<point> >& cellToEnd
+            )
+            :
+                particle::TrackingData<Cloud<findCellParticle> >(cloud),
+                cellToData_(cellToData),
+                cellToEnd_(cellToEnd)
+            {}
+
+
+        // Member functions
+
+            labelListList& cellToData()
+            {
+                return cellToData_;
+            }
+            List<List<point> >& cellToEnd()
+            {
+                return cellToEnd_;
+            }
+    };
+
+
+
+    // Constructors
+
+        //- Construct from components
+        findCellParticle
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const point& end,
+            const label data
+        );
+
+        //- Construct from Istream
+        findCellParticle
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct and return a clone
+        autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new findCellParticle(*this));
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<findCellParticle> operator()(Istream& is) const
+            {
+                return autoPtr<findCellParticle>
+                (
+                    new findCellParticle(mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        //- point to track to
+        const point& end() const
+        {
+            return end_;
+        }
+
+        //- transported label
+        label data() const
+        {
+            return data_;
+        }
+
+
+
+        // Tracking
+
+            //- Track all particles to their end point
+            bool move(trackingData&, const scalar);
+
+
+            //- Overridable function to handle the particle hitting a patch
+            //  Executed before other patch-hitting functions
+            bool hitPatch
+            (
+                const polyPatch&,
+                trackingData& td,
+                const label patchI,
+                const scalar trackFraction,
+                const tetIndices& tetIs
+            );
+
+            //- Overridable function to handle the particle hitting a wedge
+            void hitWedgePatch
+            (
+                const wedgePolyPatch&,
+                trackingData& td
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //  symmetryPlane
+            void hitSymmetryPatch
+            (
+                const symmetryPolyPatch&,
+                trackingData& td
+            );
+
+            //- Overridable function to handle the particle hitting a cyclic
+            void hitCyclicPatch
+            (
+                const cyclicPolyPatch&,
+                trackingData& td
+            );
+
+            //- Overridable function to handle the particle hitting a
+            //- processorPatch
+            void hitProcessorPatch
+            (
+                const processorPolyPatch&,
+                trackingData& td
+            );
+
+            //- Overridable function to handle the particle hitting a wallPatch
+            void hitWallPatch
+            (
+                const wallPolyPatch&,
+                trackingData& td,
+                const tetIndices&
+            );
+
+            //- Overridable function to handle the particle hitting a polyPatch
+            void hitPatch
+            (
+                const polyPatch&,
+                trackingData& td
+            );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<(Ostream&, const findCellParticle&);
+};
+
+
+template<>
+inline bool contiguous<findCellParticle>()
+{
+    return true;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "findCellParticle.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTemplateTypeNameAndDebug(Cloud<findCellParticle>, 0);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -25,6 +25,9 @@ License
 
 #include "nearWallFields.H"
 #include "wordReList.H"
+#include "findCellParticle.H"
+#include "mappedPatchBase.H"
+#include "OBJstream.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -34,6 +37,189 @@ defineTypeNameAndDebug(nearWallFields, 0);
 }
 
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::nearWallFields::calcAddressing()
+{
+    const fvMesh& mesh = refCast<const fvMesh>(obr_);
+
+    // Count number of faces
+    label nPatchFaces = 0;
+    forAllConstIter(labelHashSet, patchSet_, iter)
+    {
+        label patchI = iter.key();
+        nPatchFaces += mesh.boundary()[patchI].size();
+    }
+
+    // Global indexing
+    globalIndex globalCells(mesh.nCells());
+    globalIndex globalWalls(nPatchFaces);
+
+    if (debug)
+    {
+        Info<< "nearWallFields::calcAddressing() :"
+            << " nPatchFaces:" << globalWalls.size() << endl;
+    }
+
+    // Construct cloud
+    Cloud<findCellParticle> cloud(mesh, IDLList<findCellParticle>());
+
+    // Add particles to track to sample locations
+    nPatchFaces = 0;
+
+    forAllConstIter(labelHashSet, patchSet_, iter)
+    {
+        label patchI = iter.key();
+        const fvPatch& patch = mesh.boundary()[patchI];
+
+        vectorField nf = patch.nf();
+        vectorField faceCellCentres = patch.patch().faceCellCentres();
+
+        forAll(patch, patchFaceI)
+        {
+            label meshFaceI = patch.start()+patchFaceI;
+
+            // Find starting point on face (since faceCentre might not
+            // be on face-diagonal decomposition)
+            pointIndexHit startInfo
+            (
+                mappedPatchBase::facePoint
+                (
+                    mesh,
+                    meshFaceI,
+                    polyMesh::FACEDIAGTETS
+                )
+            );
+
+
+            point start;
+            if (startInfo.hit())
+            {
+                start = startInfo.hitPoint();
+            }
+            else
+            {
+                // Fallback: start tracking from neighbouring cell centre
+                start = faceCellCentres[patchFaceI];
+            }
+
+            const point end = start-distance_*nf[patchFaceI];
+
+            // Find tet for starting location
+            label cellI = -1;
+            label tetFaceI = -1;
+            label tetPtI = -1;
+            mesh.findCellFacePt(start, cellI, tetFaceI, tetPtI);
+
+            // Add to cloud. Add originating face as passive data
+            cloud.addParticle
+            (
+                new findCellParticle
+                (
+                    mesh,
+                    start,
+                    cellI,
+                    tetFaceI,
+                    tetPtI,
+                    end,
+                    globalWalls.toGlobal(nPatchFaces)    // passive data
+                )
+            );
+
+            nPatchFaces++;
+        }
+    }
+
+
+
+    if (debug)
+    {
+        // Dump particles
+        OBJstream str
+        (
+            mesh.time().path()
+           /"wantedTracks_" + mesh.time().timeName() + ".obj"
+        );
+        Info<< "nearWallFields::calcAddressing() :"
+            << "Dumping tracks to " << str.name() << endl;
+
+        forAllConstIter(Cloud<findCellParticle>, cloud, iter)
+        {
+            const findCellParticle& tp = iter();
+            str.write(linePointRef(tp.position(), tp.end()));
+        }
+    }
+
+
+
+    // Per cell: empty or global wall index and end location
+    cellToWalls_.setSize(mesh.nCells());
+    cellToSamples_.setSize(mesh.nCells());
+
+    // Database to pass into findCellParticle::move
+    findCellParticle::trackingData td(cloud, cellToWalls_, cellToSamples_);
+
+    // Track all particles to their end position.
+    scalar maxTrackLen = 2.0*mesh.bounds().mag();
+
+
+    //Debug: collect start points
+    pointField start;
+    if (debug)
+    {
+        start.setSize(nPatchFaces);
+        nPatchFaces = 0;
+        forAllConstIter(Cloud<findCellParticle>, cloud, iter)
+        {
+            const findCellParticle& tp = iter();
+            start[nPatchFaces++] = tp.position();
+        }
+    }
+
+
+    cloud.move(td, maxTrackLen);
+
+
+    // Rework cell-to-globalpatchface into a map
+    List<Map<label> > compactMap;
+    getPatchDataMapPtr_.reset
+    (
+        new mapDistribute
+        (
+            globalWalls,
+            cellToWalls_,
+            compactMap
+        )
+    );
+
+
+    // Debug: dump resulting tracks
+    if (debug)
+    {
+        getPatchDataMapPtr_().distribute(start);
+        {
+            OBJstream str
+            (
+                mesh.time().path()
+               /"obtainedTracks_" + mesh.time().timeName() + ".obj"
+            );
+            Info<< "nearWallFields::calcAddressing() :"
+                << "Dumping obtained to " << str.name() << endl;
+
+            forAll(cellToWalls_, cellI)
+            {
+                const List<point>& ends = cellToSamples_[cellI];
+                const labelList& cData = cellToWalls_[cellI];
+                forAll(cData, i)
+                {
+                    str.write(linePointRef(ends[i], start[cData[i]]));
+                }
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::nearWallFields::nearWallFields
@@ -127,16 +313,11 @@ void Foam::nearWallFields::read(const dictionary& dict)
             reverseFieldMap_.insert(sampleFldName, fldName);
         }
 
-        Info<< type() << " " << name_ << ": Creating " << fieldMap_.size()
+        Info<< type() << " " << name_ << ": Sampling " << fieldMap_.size()
             << " fields" << endl;
 
-        createFields(vsf_);
+        // Do analysis
-        createFields(vvf_);
+        calcAddressing();
-        createFields(vSpheretf_);
-        createFields(vSymmtf_);
-        createFields(vtf_);
-
-        Info<< endl;
     }
 }
 
@@ -147,15 +328,6 @@ void Foam::nearWallFields::execute()
     {
         Info<< "nearWallFields:execute()" << endl;
     }
-
-    //if (active_)
-    //{
-    //    sampleFields(vsf_);
-    //    sampleFields(vvf_);
-    //    sampleFields(vSpheretf_);
-    //    sampleFields(vSymmtf_);
-    //    sampleFields(vtf_);
-    //}
 }
 
 
@@ -165,8 +337,6 @@ void Foam::nearWallFields::end()
     {
         Info<< "nearWallFields:end()" << endl;
     }
-    // Update fields
-    execute();
 }
 
 
@@ -186,6 +356,28 @@ void Foam::nearWallFields::write()
     // Do nothing
     if (active_)
     {
+        if
+        (
+            fieldMap_.size()
+         && vsf_.empty()
+         && vvf_.empty()
+         && vSpheretf_.empty()
+         && vSymmtf_.empty()
+         && vtf_.empty()
+        )
+        {
+            Info<< type() << " " << name_ << ": Creating " << fieldMap_.size()
+                << " fields" << endl;
+
+            createFields(vsf_);
+            createFields(vvf_);
+            createFields(vSpheretf_);
+            createFields(vSymmtf_);
+            createFields(vtf_);
+
+            Info<< endl;
+        }
+
         Info<< type() << " " << name_ << " output:" << nl;
 
         Info<< "    Writing sampled fields to " << obr_.time().timeName()

src/postProcessing/functionObjects/field/nearWallFields/nearWallFields.H

@@ -76,6 +76,7 @@ SourceFiles
 #include "OFstream.H"
 #include "volFields.H"
 #include "Tuple2.H"
+#include "interpolationCellPoint.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -122,7 +123,21 @@ protected:
             //- From resulting back to original field
             HashTable<word> reverseFieldMap_;
 
-        //- Locally constructed fields
+
+        // Calculated addressing
+
+            //- From cell to seed patch faces
+            labelListList cellToWalls_;
+
+            //- From cell to tracked end point
+            List<List<point> > cellToSamples_;
+
+            //- Map from cell based data back to patch based data
+            autoPtr<mapDistribute> getPatchDataMapPtr_;
+
+
+        // Locally constructed fields
+
             PtrList<volScalarField> vsf_;
             PtrList<volVectorField> vvf_;
             PtrList<volSphericalTensorField> vSpheretf_;
@@ -132,11 +147,8 @@ protected:
 
     // Protected Member Functions
 
-        //- Disallow default bitwise copy construct
+        //- Calculate addressing from cells back to patch faces
-        nearWallFields(const nearWallFields&);
+        void calcAddressing();
-
-        //- Disallow default bitwise assignment
-        void operator=(const nearWallFields&);
 
         template<class Type>
         void createFields
@@ -144,6 +156,14 @@ protected:
             PtrList<GeometricField<Type, fvPatchField, volMesh> >&
         ) const;
 
+        //- Override boundary fields with sampled values
+        template<class Type>
+        void sampleBoundaryField
+        (
+            const interpolationCellPoint<Type>& interpolator,
+            GeometricField<Type, fvPatchField, volMesh>& fld
+        ) const;
+
         template<class Type>
         void sampleFields
         (
@@ -151,6 +171,14 @@ protected:
         ) const;
 
 
+private:
+
+        //- Disallow default bitwise copy construct
+        nearWallFields(const nearWallFields&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const nearWallFields&);
+
 public:
 
     //- Runtime type information

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "nearWallFields.H"
-#include "mappedFieldFvPatchFields.H"
-#include "interpolationCellPoint.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,38 +62,58 @@ void Foam::nearWallFields::createFields
                 io.rename(sampleFldName);
 
                 sflds.set(sz, new vfType(io, fld));
-                vfType& sampleFld = sflds[sz];
 
-                // Reset the bcs to be mapped
+                Info<< "    created " << sflds[sz].name() << " to sample "
+                    << fld.name() << endl;
+            }
+        }
+    }
+}
+
+
+template<class Type>
+void Foam::nearWallFields::sampleBoundaryField
+(
+    const interpolationCellPoint<Type>& interpolator,
+    GeometricField<Type, fvPatchField, volMesh>& fld
+) const
+{
+    // Construct flat fields for all patch faces to be sampled
+    Field<Type> sampledValues(getPatchDataMapPtr_().constructSize());
+
+    forAll(cellToWalls_, cellI)
+    {
+        const labelList& cData = cellToWalls_[cellI];
+
+        forAll(cData, i)
+        {
+            const point& samplePt = cellToSamples_[cellI][i];
+            sampledValues[cData[i]] = interpolator.interpolate(samplePt, cellI);
+        }
+    }
+
+    // Send back sampled values to patch faces
+    getPatchDataMapPtr_().reverseDistribute
+    (
+        getPatchDataMapPtr_().constructSize(),
+        sampledValues
+    );
+
+    // Pick up data
+    label nPatchFaces = 0;
     forAllConstIter(labelHashSet, patchSet_, iter)
     {
         label patchI = iter.key();
 
-                    sampleFld.boundaryField().set
+        fvPatchField<Type>& pfld = fld.boundaryField()[patchI];
-                    (
-                        patchI,
-                        new mappedFieldFvPatchField<Type>
-                        (
-                            sampleFld.mesh().boundary()[patchI],
-                            sampleFld.dimensionedInternalField(),
 
-                            sampleFld.mesh().name(),
+        Field<Type> newFld(pfld.size());
-                            mappedPatchBase::NEARESTCELL,
+        forAll(pfld, i)
-                            word::null,     // samplePatch
+        {
-                            -distance_,
+            newFld[i] = sampledValues[nPatchFaces++];
-
-                            sampleFld.name(),       // fieldName
-                            false,                  // setAverage
-                            pTraits<Type>::zero,    // average
-                            interpolationCellPoint<Type>::typeName
-                        )
-                    );
         }
 
-                Info<< "    created " << sampleFld.name() << " to sample "
+        pfld == newFld;
-                    << fld.name() << endl;
-            }
-        }
     }
 }
 
@@ -115,8 +133,12 @@ void Foam::nearWallFields::sampleFields
 
         // Take over internal and boundary values
         sflds[i] == fld;
-        // Evaluate to update the mapped
+
-        sflds[i].correctBoundaryConditions();
+        // Construct interpolation method
+        interpolationCellPoint<Type> interpolator(fld);
+
+        // Override sampled values
+        sampleBoundaryField(interpolator, sflds[i]);
     }
 }
 

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_;
 

tutorials/combustion/LTSReactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties

@@ -17,7 +17,7 @@ FoamFile
 
 chemistryType
 {
-    chemistrySolver   ode;
+    chemistrySolver   EulerImplicit;
     chemistryThermo   psi;
 }
 
@@ -33,8 +33,9 @@ EulerImplicitCoeffs
 
 odeCoeffs
 {
-    solver          KRR4;
+    solver          Rosenbrock43;
-    eps             0.05;
+    absTol          1e-12;
+    relTol          0.01;
 }
 
 // ************************************************************************* //

tutorials/combustion/chemFoam/gri/constant/chemistryProperties

@@ -33,8 +33,9 @@ EulerImplicitCoeffs
 
 odeCoeffs
 {
-    solver          SIBS;
+    solver          Rosenbrock43;
-    eps             0.01;
+    absTol          1e-12;
+    relTol          1e-2;
 }
 
 

tutorials/combustion/chemFoam/h2/constant/chemistryProperties

@@ -28,7 +28,8 @@ initialChemicalTimeStep 1e-10;
 odeCoeffs
 {
     solver          SIBS;
-    eps             0.001;
+    absTol          1e-12;
+    relTol          1e-2;
 }
 
 

tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties

@@ -28,7 +28,8 @@ initialChemicalTimeStep 1e-10;
 odeCoeffs
 {
     solver          SIBS;
-    eps             1e-03;
+    absTol          1e-12;
+    relTol          1e-3;
 }
 
 

tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties

@@ -28,7 +28,8 @@ initialChemicalTimeStep 1e-10;
 odeCoeffs
 {
     solver          SIBS;
-    eps             1e-04;
+    absTol          1e-14;
+    relTol          1e-4;
 }
 
 

tutorials/combustion/fireFoam/les/oppositeBurningPanels/constant/panelRegion/chemistryProperties

@@ -28,7 +28,8 @@ initialChemicalTimeStep 1e-07;
 odeCoeffs
 {
     solver          SIBS;
-    eps             0.05;
+    absTol          1e-12;
+    relTol          0.01;
 }
 
 

tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/constant/chemistryProperties

@@ -33,8 +33,9 @@ EulerImplicitCoeffs
 
 odeCoeffs
 {
-    solver          KRR4;
+    solver          Rosenbrock43;
-    eps             0.05;
+    absTol          1e-12;
+    relTol          0.01;
 }
 
 // ************************************************************************* //