GIT: resolve conflict

applications/test/ODE/Test-ODE.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
@@ -29,7 +29,6 @@ Description
 #include "IOmanip.H"
 #include "ODESystem.H"
 #include "ODESolver.H"
-#include "RK.H"
 
 using namespace Foam;
 
@@ -110,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;
@@ -125,27 +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);
 
-        scalarField yScale(ode.nEqns(), 1.0);
+        scalar dxEst = 0.6;
-        scalar hEst = 0.6;
+        scalar dxNext = dxEst;
-        scalar hDid, hNext;
 
-        odeSolver->solve(ode, x, y, dydx, eps, yScale, 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;
@@ -161,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

@@ -1,19 +1,18 @@
-ODE  =  ODE
+ODESolvers/ODESolver/ODESolver.C
-ODESolvers  =  ODESolvers
+ODESolvers/ODESolver/ODESolverNew.C
-ODESolversODESolver  =  ODESolvers/ODESolver
-ODESolversKRR4  =  ODESolvers/KRR4
-ODESolversRK  =  ODESolvers/RK
-ODESolversSIBS  =  ODESolvers/SIBS
 
-$(ODESolversODESolver)/ODESolver.C
+ODESolvers/adaptiveSolver/adaptiveSolver.C
-$(ODESolversODESolver)/ODESolverNew.C
+ODESolvers/Euler/Euler.C
-
+ODESolvers/EulerSI/EulerSI.C
-$(ODESolversRK)/RK.C
+ODESolvers/Trapezoid/Trapezoid.C
-
+ODESolvers/RKF45/RKF45.C
-$(ODESolversKRR4)/KRR4.C
+ODESolvers/RKCK45/RKCK45.C
-
+ODESolvers/RKDP45/RKDP45.C
-$(ODESolversSIBS)/SIBS.C
+ODESolvers/Rosenbrock21/Rosenbrock21.C
-$(ODESolversSIBS)/SIMPR.C
+ODESolvers/Rosenbrock43/Rosenbrock43.C
-$(ODESolversSIBS)/polyExtrapolate.C
+ODESolvers/rodas43/rodas43.C
+ODESolvers/SIBS/SIBS.C
+ODESolvers/SIBS/SIMPR.C
+ODESolvers/SIBS/polyExtrapolate.C
 
 LIB = $(FOAM_LIBBIN)/libODE

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

@@ -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
@@ -22,21 +22,33 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::KRR4
+    Foam::EulerSI
 
 Description
-    Foam::KRR4
+    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
-    KRR4CK.C
+    EulerSI.C
-    KRR4QS.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef KRR4_H
+#ifndef EulerSI_H
-#define KRR4_H
+#define EulerSI_H
 
 #include "ODESolver.H"
+#include "adaptiveSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,66 +56,56 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class KRR4 Declaration
+                           Class EulerSI Declaration
 \*---------------------------------------------------------------------------*/
 
-class KRR4
+class EulerSI
 :
-    public ODESolver
+    public ODESolver,
+    public adaptiveSolver
 {
     // Private data
 
-        mutable scalarField yTemp_;
+        mutable scalarField err_;
-        mutable scalarField dydxTemp_;
+        mutable scalarField dydx_;
-        mutable scalarField g1_;
-        mutable scalarField g2_;
-        mutable scalarField g3_;
-        mutable scalarField g4_;
-        mutable scalarField yErr_;
         mutable scalarField dfdx_;
         mutable scalarSquareMatrix dfdy_;
         mutable scalarSquareMatrix a_;
         mutable labelList pivotIndices_;
 
-        static const int maxtry = 40;
-
-        static const scalar safety, grow, pgrow, shrink, pshrink, errcon;
-
-        static const scalar
-            gamma,
-            a21, a31, a32,
-            c21, c31, c32, c41, c42, c43,
-            b1, b2, b3, b4,
-            e1, e2, e3, e4,
-            c1X, c2X, c3X, c4X,
-            a2X, a3X;
-
 
 public:
 
     //- Runtime type information
-    TypeName("KRR4");
+    TypeName("EulerSI");
 
 
     // Constructors
 
         //- Construct from ODE
-        KRR4(const ODESystem& 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,
-            scalarField& dydx,
+            scalar& dxTry
-            const scalar eps,
-            const scalarField& yScale,
-            const scalar hTry,
-            scalar& hDid,
-            scalar& hNext
         ) const;
 };
 

src/ODE/ODESolvers/KRR4/KRR4.C

@@ -1,224 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "KRR4.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(KRR4, 0);
-
-    addToRunTimeSelectionTable(ODESolver, KRR4, ODE);
-
-const scalar
-    KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25,
-    KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296;
-
-const scalar
-    KRR4::gamma = 1.0/2.0,
-    KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0,
-    KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0,
-    KRR4::c41 = -112.0/125.0, KRR4::c42 = -54.0/125.0, KRR4::c43 = -2.0/5.0,
-    KRR4::b1 = 19.0/9.0, KRR4::b2 = 1.0/2.0, KRR4::b3 = 25.0/108.0,
-    KRR4::b4 = 125.0/108.0,
-    KRR4::e1 = 17.0/54.0, KRR4::e2 = 7.0/36.0, KRR4::e3 = 0.0,
-    KRR4::e4 = 125.0/108.0,
-    KRR4::c1X = 1.0/2.0, KRR4::c2X = -3.0/2.0, KRR4::c3X = 121.0/50.0,
-    KRR4::c4X = 29.0/250.0,
-    KRR4::a2X = 1.0, KRR4::a3X = 3.0/5.0;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::KRR4::KRR4(const ODESystem& ode)
-:
-    ODESolver(ode),
-    yTemp_(n_, 0.0),
-    dydxTemp_(n_, 0.0),
-    g1_(n_, 0.0),
-    g2_(n_, 0.0),
-    g3_(n_, 0.0),
-    g4_(n_, 0.0),
-    yErr_(n_, 0.0),
-    dfdx_(n_, 0.0),
-    dfdy_(n_, n_, 0.0),
-    a_(n_, n_, 0.0),
-    pivotIndices_(n_, 0.0)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::KRR4::solve
-(
-    const ODESystem& ode,
-    scalar& x,
-    scalarField& y,
-    scalarField& dydx,
-    const scalar eps,
-    const scalarField& yScale,
-    const scalar hTry,
-    scalar& hDid,
-    scalar& hNext
-) const
-{
-    scalar xTemp = x;
-    yTemp_ = y;
-    dydxTemp_ = dydx;
-
-    ode.jacobian(xTemp, yTemp_, dfdx_, dfdy_);
-
-    scalar h = hTry;
-
-    for (register label jtry=0; jtry<maxtry; jtry++)
-    {
-        for (register label i=0; i<n_; i++)
-        {
-            for (register label j=0; j<n_; j++)
-            {
-                a_[i][j] = -dfdy_[i][j];
-            }
-
-            a_[i][i] += 1.0/(gamma*h);
-        }
-
-        LUDecompose(a_, pivotIndices_);
-
-        for (register label i=0; i<n_; i++)
-        {
-            g1_[i] = dydxTemp_[i] + h*c1X*dfdx_[i];
-        }
-
-        LUBacksubstitute(a_, pivotIndices_, g1_);
-
-        for (register label i=0; i<n_; i++)
-        {
-            y[i] = yTemp_[i] + a21*g1_[i];
-        }
-
-        x = xTemp + a2X*h;
-        ode.derivatives(x, y, dydx_);
-
-        for (register label i=0; i<n_; i++)
-        {
-            g2_[i] = dydx_[i] + h*c2X*dfdx_[i] + c21*g1_[i]/h;
-        }
-
-        LUBacksubstitute(a_, pivotIndices_, g2_);
-
-        for (register label i=0; i<n_; i++)
-        {
-            y[i] = yTemp_[i] + a31*g1_[i] + a32*g2_[i];
-        }
-
-        x = xTemp + a3X*h;
-        ode.derivatives(x, y, dydx_);
-
-        for (register label i=0; i<n_; i++)
-        {
-            g3_[i] = dydx[i] + h*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/h;
-        }
-
-        LUBacksubstitute(a_, pivotIndices_, g3_);
-
-        for (register label i=0; i<n_; i++)
-        {
-            g4_[i] = dydx_[i] + h*c4X*dfdx_[i]
-                + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/h;
-        }
-
-        LUBacksubstitute(a_, pivotIndices_, g4_);
-
-        for (register label i=0; i<n_; i++)
-        {
-            y[i] = yTemp_[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i];
-            yErr_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i];
-        }
-
-        x = xTemp + h;
-
-        if (x == xTemp)
-        {
-            FatalErrorIn
-            (
-                "void Foam::KRR4::solve"
-                "("
-                    "const ODESystem&, "
-                    "scalar&, "
-                    "scalarField&, "
-                    "scalarField&, "
-                    "const scalar, "
-                    "const scalarField&, "
-                    "const scalar, "
-                    "scalar&, "
-                    "scalar&"
-                ") const"
-            )   << "solver stalled: step size = 0"
-                << exit(FatalError);
-        }
-
-        scalar maxErr = 0.0;
-        for (register label i=0; i<n_; i++)
-        {
-            maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
-        }
-        maxErr /= eps;
-
-        if (maxErr <= 1.0)
-        {
-            hDid = h;
-            hNext = (maxErr > errcon ? safety*h*pow(maxErr, pgrow) : grow*h);
-            return;
-        }
-        else
-        {
-            hNext = safety*h*pow(maxErr, pshrink);
-            h = (h >= 0.0 ? max(hNext, shrink*h) : min(hNext, shrink*h));
-        }
-    }
-
-    FatalErrorIn
-    (
-        "void Foam::KRR4::solve"
-        "("
-            "const ODESystem&, "
-            "scalar&, "
-            "scalarField&, "
-            "scalarField&, "
-            "const scalar, "
-            "const scalarField&, "
-            "const scalar, "
-            "scalar&, "
-            "scalar&"
-        ") const"
-    )   << "Maximum number of solver iterations exceeded"
-        << exit(FatalError);
-}
-
-
-// ************************************************************************* //

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,82 +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()),
-    yScale_(n_),
+    absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
-    dydx_(n_)
+    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;
 
-        for (label i=0; i<n_; i++)
+        // Check if this is a truncated step and set dxEst to integrate to xEnd
+        if ((x + dxEst - xEnd)*(x + dxEst - xStart) > 0)
         {
-            yScale_[i] = mag(y[i]) + mag(dydx_[i]*h) + SMALL;
+            truncated = true;
+            dxEst = xEnd - x;
         }
 
-        if ((x + h - xEnd)*(x + h - xStart) > 0.0)
+        // Integrate as far as possible up to dxEst
-        {
+        solve(ode, x, y, dxEst);
-            h = xEnd - x;
-            hPrev = hNext;
-        }
 
-        hNext = 0;
+        // Check if reached xEnd
-        scalar hDid;
+        if ((x - xEnd)*(xEnd - xStart) >= 0)
-        solve(ode, x, y, dydx_, eps, yScale_, h, hDid, hNext);
-
-        if ((x - xEnd)*(xEnd - xStart) >= 0.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,14 +53,30 @@ class ODESolver
 
 protected:
 
-    // Private data
+    // Protected data
 
+        //- Size of the ODESystem
         label n_;
-        mutable scalarField yScale_;
+
-        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&);
@@ -81,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
@@ -98,8 +122,8 @@ public:
         //- Select null constructed
         static autoPtr<ODESolver> New
         (
-            const word& ODESolverTypeName,
+            const ODESystem& ode,
-            const ODESystem& ode
+            const dictionary& dict
         );
 
 
@@ -110,28 +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 scalarField& yScale,
-            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/RK/RK.C

@@ -1,204 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "RK.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(RK, 0);
-    addToRunTimeSelectionTable(ODESolver, RK, ODE);
-
-const scalar
-    RK::safety=0.9, RK::pGrow=-0.2, RK::pShrink=-0.25, RK::errCon=1.89e-4;
-
-const scalar
-    RK::a2 = 0.2, RK::a3 = 0.3, RK::a4 = 0.6, RK::a5 = 1.0, RK::a6 = 0.875,
-    RK::b21 = 0.2, RK::b31 = 3.0/40.0, RK::b32 = 9.0/40.0,
-    RK::b41 = 0.3, RK::b42 = -0.9, RK::b43 = 1.2,
-    RK::b51 = -11.0/54.0, RK::b52 = 2.5, RK::b53 = -70.0/27.0,
-    RK::b54 = 35.0/27.0,
-    RK::b61 = 1631.0/55296.0, RK::b62 = 175.0/512.0, RK::b63 = 575.0/13824.0,
-    RK::b64 = 44275.0/110592.0, RK::b65 = 253.0/4096.0,
-    RK::c1 = 37.0/378.0, RK::c3 = 250.0/621.0,
-    RK::c4 = 125.0/594.0, RK::c6 = 512.0/1771.0,
-    RK::dc1 = RK::c1 - 2825.0/27648.0, RK::dc3 = RK::c3 - 18575.0/48384.0,
-    RK::dc4 = RK::c4 - 13525.0/55296.0, RK::dc5 = -277.00/14336.0,
-    RK::dc6 = RK::c6 - 0.25;
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::RK::RK(const ODESystem& ode)
-:
-    ODESolver(ode),
-    yTemp_(n_, 0.0),
-    ak2_(n_, 0.0),
-    ak3_(n_, 0.0),
-    ak4_(n_, 0.0),
-    ak5_(n_, 0.0),
-    ak6_(n_, 0.0),
-    yErr_(n_, 0.0),
-    yTemp2_(n_, 0.0)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::RK::solve
-(
-    const ODESystem& ode,
-    const scalar x,
-    const scalarField& y,
-    const scalarField& dydx,
-    const scalar h,
-    scalarField& yout,
-    scalarField& yerr
-) const
-{
-    forAll(yTemp_, i)
-    {
-        yTemp_[i] = y[i] + b21*h*dydx[i];
-    }
-
-    ode.derivatives(x + a2*h, yTemp_, ak2_);
-
-    forAll(yTemp_, i)
-    {
-        yTemp_[i] = y[i] + h*(b31*dydx[i] + b32*ak2_[i]);
-    }
-
-    ode.derivatives(x + a3*h, yTemp_, ak3_);
-
-    forAll(yTemp_, i)
-    {
-        yTemp_[i] = y[i] + h*(b41*dydx[i] + b42*ak2_[i] + b43*ak3_[i]);
-    }
-
-    ode.derivatives(x + a4*h, yTemp_, ak4_);
-
-    forAll(yTemp_, i)
-    {
-        yTemp_[i] = y[i]
-          + h*(b51*dydx[i] + b52*ak2_[i] + b53*ak3_[i] + b54*ak4_[i]);
-    }
-
-    ode.derivatives(x + a5*h, yTemp_, ak5_);
-
-    forAll(yTemp_, i)
-    {
-        yTemp_[i] = y[i]
-          + h*
-            (
-                b61*dydx[i] + b62*ak2_[i] + b63*ak3_[i]
-              + b64*ak4_[i] + b65*ak5_[i]
-            );
-    }
-
-    ode.derivatives(x + a6*h, yTemp_, ak6_);
-
-    forAll(yout, i)
-    {
-        yout[i] = y[i]
-          + h*(c1*dydx[i] + c3*ak3_[i] + c4*ak4_[i] + c6*ak6_[i]);
-    }
-
-    forAll(yerr, i)
-    {
-        yerr[i] =
-            h*
-            (
-                dc1*dydx[i] + dc3*ak3_[i] + dc4*ak4_[i]
-              + dc5*ak5_[i] + dc6*ak6_[i]
-            );
-    }
-}
-
-
-void Foam::RK::solve
-(
-    const ODESystem& ode,
-    scalar& x,
-    scalarField& y,
-    scalarField& dydx,
-    const scalar eps,
-    const scalarField& yScale,
-    const scalar hTry,
-    scalar& hDid,
-    scalar& hNext
-) const
-{
-    scalar h = hTry;
-    scalar maxErr = 0.0;
-
-    for (;;)
-    {
-        solve(ode, x, y, dydx, h, yTemp2_, yErr_);
-
-        maxErr = 0.0;
-        for (register label i=0; i<n_; i++)
-        {
-            maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
-        }
-        maxErr /= eps;
-
-        if (maxErr <= 1.0)
-        {
-            break;
-        }
-
-        {
-            scalar hTemp = safety*h*pow(maxErr, pShrink);
-            h = (h >= 0.0 ? max(hTemp, 0.1*h) : min(hTemp, 0.1*h));
-        }
-
-        if (h < VSMALL)
-        {
-            FatalErrorIn("RK::solve")
-                << "stepsize underflow"
-                << exit(FatalError);
-        }
-    }
-
-    hDid = h;
-
-    x += h;
-    y = yTemp2_;
-
-    if (maxErr > errCon)
-    {
-        hNext = safety*h*pow(maxErr, pGrow);
-    }
-    else
-    {
-        hNext = 5.0*h;
-    }
-}
-
-
-// ************************************************************************* //

src/ODE/ODESolvers/RKCK45/RKCK45.C

@@ -0,0 +1,167 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "RKCK45.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(RKCK45, 0);
+    addToRunTimeSelectionTable(ODESolver, RKCK45, dictionary);
+
+const scalar
+    RKCK45::c2 = 1.0/5.0,
+    RKCK45::c3 = 3.0/10.0,
+    RKCK45::c4 = 3.0/5.0,
+    RKCK45::c5 = 1.0,
+    RKCK45::c6 = 7.0/8.0,
+
+    RKCK45::a21 = 1.0/5.0,
+    RKCK45::a31 = 3.0/40.0,
+    RKCK45::a32 = 9.0/40.0,
+    RKCK45::a41 = 3.0/10.0,
+    RKCK45::a42 = -9.0/10.0,
+    RKCK45::a43 = 6.0/5.0,
+    RKCK45::a51 = -11.0/54.0,
+    RKCK45::a52 = 5.0/2.0,
+    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::RKCK45::RKCK45(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::RKCK45::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_);
+
+    forAll(y, i)
+    {
+        y[i] = y0[i]
+          + dx*(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b6*k6_[i]);
+    }
+
+    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::RKCK45::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/RKCK45/RKCK45.H

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::RKCK45
+
+Description
+    4/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
+
+    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
+    RKCK45.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef RKCK45_H
+#define RKCK45_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class RKCK45 Declaration
+\*---------------------------------------------------------------------------*/
+
+class RKCK45
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        //- RKCK 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, 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("RKCK45");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        RKCK45(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/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

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "Rosenbrock21.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(Rosenbrock21, 0);
+    addToRunTimeSelectionTable(ODESolver, Rosenbrock21, dictionary);
+
+const scalar
+    Rosenbrock21::gamma = 1 + 1.0/std::sqrt(2.0),
+    Rosenbrock21::a21 = 1.0/gamma,
+    Rosenbrock21::c2 = 1.0,
+    Rosenbrock21::c21 = -2.0/gamma,
+    Rosenbrock21::b1 = (3.0/2.0)/gamma,
+    Rosenbrock21::b2 = (1.0/2.0)/gamma,
+    Rosenbrock21::e1 = b1 - 1.0/gamma,
+    Rosenbrock21::e2 = b2,
+    Rosenbrock21::d1 = 1,
+    Rosenbrock21::d2 = -1;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict)
+:
+    ODESolver(ode, dict),
+    adaptiveSolver(ode, dict),
+    k1_(n_),
+    k2_(n_),
+    err_(n_),
+    dydx_(n_),
+    dfdx_(n_),
+    dfdy_(n_, n_),
+    a_(n_, n_),
+    pivotIndices_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::Rosenbrock21::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 error and update state:
+    forAll(y, i)
+    {
+        y[i] = y0[i] + b1*k1_[i] + b2*k2_[i];
+        err_[i] = e1*k1_[i] + e2*k2_[i];
+    }
+
+    return normalizeError(y0, y, err_);
+}
+
+
+void Foam::Rosenbrock21::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    L-stable embedded Rosenbrock ODE solver of order (1)2.
+
+    References:
+    \verbatim
+        "A second-order Rosenbrock method applied to
+         photochemical dispersion problems",
+        J. G. Verwer,
+        E. J. Spee,
+        J. G. Blom,
+        W. Hundsdorfer,
+        Siam Journal on Scientific Computing 01/1999; 20(4):1456-1480.
+    \endverbatim
+
+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,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    // L-Stable constants from Hairer et. al.
+    Rosenbrock43::a21 = 2,
+    Rosenbrock43::a31 = 1.867943637803922,
+    Rosenbrock43::a32 = 0.2344449711399156,
+
+    Rosenbrock43::c21 = -7.137615036412310,
+    Rosenbrock43::c31 = 2.580708087951457,
+    Rosenbrock43::c32 = 0.6515950076447975,
+    Rosenbrock43::c41 = -2.137148994382534,
+    Rosenbrock43::c42 = -0.3214669691237626,
+    Rosenbrock43::c43 = -0.6949742501781779,
+
+    Rosenbrock43::b1 = 2.255570073418735,
+    Rosenbrock43::b2 = 0.2870493262186792,
+    Rosenbrock43::b3 = 0.435317943184018,
+    Rosenbrock43::b4 = 1.093502252409163,
+
+    Rosenbrock43::e1 = -0.2815431932141155,
+    Rosenbrock43::e2 = -0.0727619912493892,
+    Rosenbrock43::e3 = -0.1082196201495311,
+    Rosenbrock43::e4 = -1.093502252409163,
+
+    Rosenbrock43::gamma = 0.57282,
+    Rosenbrock43::c2 = 1.14564,
+    Rosenbrock43::c3 = 0.65521686381559,
+
+    Rosenbrock43::d1 = 0.57282,
+    Rosenbrock43::d2 = -1.769193891319233,
+    Rosenbrock43::d3 = 0.7592633437920482,
+    Rosenbrock43::d4 = -0.1049021087100450;
+
+    // Constants by Shampine
+    // More accurate than the L-Stable coefficients for small step-size
+    // but less stable for large step-size
+    /*
+    Rosenbrock43::a21 = 2,
+    Rosenbrock43::a31 = 48.0/25.0,
+    Rosenbrock43::a32 = 6.0/25.0,
+
+    Rosenbrock43::c21 = -8,
+    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,
+    Rosenbrock43::e4 = 125.0/108.0,
+
+    Rosenbrock43::gamma = 0.5,
+    Rosenbrock43::c2 = 1,
+    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

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::Rosenbrock43
+
+Description
+    L-stable embedded Rosenbrock ODE solver of order (3)4.
+
+    Based on code from:
+    \verbatim
+        "Solving Ordinary Differential Equations II: Stiff
+         and Differential-Algebraic Problems, second edition",
+        Hairer, E.,
+        Nørsett, S.,
+        Wanner, G.,
+        Springer-Verlag, Berlin. 1996.
+    \endverbatim
+
+    Also provided are the constants from:
+    \verbatim
+        "Implementation of Rosenbrock Methods"
+        Shampine, L. F.,
+        ACM Transactions on Mathematical Software, vol. 8, 1982, pp. 93–113.
+    \endverbatim
+    with which the scheme is more accurate than with the L-Stable coefficients
+    for small step-size but less stable for large step-size.
+
+SourceFiles
+    Rosenbrock43.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Rosenbrock43_H
+#define Rosenbrock43_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Rosenbrock43 Declaration
+\*---------------------------------------------------------------------------*/
+
+class Rosenbrock43
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        mutable scalarField k1_;
+        mutable scalarField k2_;
+        mutable scalarField k3_;
+        mutable scalarField k4_;
+        mutable scalarField err_;
+        mutable scalarField dydx_;
+        mutable scalarField dfdx_;
+        mutable scalarSquareMatrix dfdy_;
+        mutable scalarSquareMatrix a_;
+        mutable labelList pivotIndices_;
+
+        static const scalar
+            a21, a31, a32,
+            c21, c31, c32,
+            c41, c42, c43,
+            b1, b2, b3, b4,
+            e1, e2, e3, e4,
+            gamma,
+            c2, c3,
+            d1, d2, d3, d4;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Rosenbrock43");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        Rosenbrock43(const ODESystem& ode, const dictionary& dict);
+
+
+    // Member Functions
+
+        //- Solve a single step dx and return the error
+        scalar solve
+        (
+            const ODESystem& ode,
+            const scalar x0,
+            const scalarField& y0,
+            const scalarField& dydx0,
+            const scalar dx,
+            scalarField& y
+        ) const;
+
+        //- Solve the ODE system and the update the state
+        void solve
+        (
+            const ODESystem& ode,
+            scalar& x,
+            scalarField& y,
+            scalar& dxTry
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/ODESolvers/SIBS/SIBS.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,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,20 +73,18 @@ void Foam::SIBS::solve
     const ODESystem& ode,
     scalar& x,
     scalarField& y,
-    scalarField& dydx,
+    scalar& dxTry
-    const scalar eps,
-    const scalarField& yScale,
-    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++)
@@ -104,7 +102,7 @@ void Foam::SIBS::solve
             }
         }
 
-        epsOld_ = eps;
+        epsOld_ = relTol_[0];
         a_[0] += n_;
 
         for (register label k=0; k<kMaxX_; k++)
@@ -124,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_;
@@ -154,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_);
@@ -164,9 +161,12 @@ void Foam::SIBS::solve
                 maxErr = SMALL;
                 for (register label i=0; i<n_; i++)
                 {
-                    maxErr = max(maxErr, mag(yErr_[i]/yScale[i]));
+                    maxErr = max
+                    (
+                        maxErr,
+                        mag(yErr_[i])/(absTol_[i] + relTol_[i]*mag(yTemp_[i]))
+                    );
                 }
-                maxErr /= eps;
                 km = k - 1;
                 err_[km] = pow(maxErr/safe1, 1.0/(2*km + 3));
             }
@@ -214,7 +214,6 @@ void Foam::SIBS::solve
     }
 
     x = xNew_;
-    hDid = h;
     first_ = 0;
     scalar wrkmin = GREAT;
     scalar scale = 1.0;
@@ -231,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,12 +126,7 @@ public:
             const ODESystem& ode,
             scalar& x,
             scalarField& y,
-            scalarField& dydx,
+            scalar& dxTry
-            const scalar eps,
-            const scalarField& yScale,
-            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

@@ -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
@@ -22,21 +22,21 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::RK
+    Foam::Trapezoid
 
 Description
-    Runge-Kutta ODE solver
+    Trapezoidal ODE solver of order (1)2.
 
 SourceFiles
-    RKCK.C
+    Trapezoid.C
-    RKQS.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef RK_H
+#ifndef Trapezoid_H
-#define RK_H
+#define Trapezoid_H
 
 #include "ODESolver.H"
+#include "adaptiveSolver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,71 +44,51 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class RK Declaration
+                           Class Trapezoid Declaration
 \*---------------------------------------------------------------------------*/
 
-class RK
+class Trapezoid
 :
-    public ODESolver
+    public ODESolver,
+    public adaptiveSolver
 {
     // Private data
-        static const scalar safety, pGrow, pShrink, errCon;
 
-        static const scalar
+        mutable scalarField err_;
-            a2, a3, a4, a5, a6,
-            b21, b31, b32, b41, b42, b43,
-            b51, b52, b53, b54, b61, b62, b63, b64, b65,
-            c1, c3, c4, c6,
-            dc1, dc3, dc4, dc5, dc6;
 
 
-        mutable scalarField yTemp_;
-        mutable scalarField ak2_;
-        mutable scalarField ak3_;
-        mutable scalarField ak4_;
-        mutable scalarField ak5_;
-        mutable scalarField ak6_;
-
-        mutable scalarField yErr_;
-        mutable scalarField yTemp2_;
-
 public:
 
     //- Runtime type information
-    TypeName("RK");
+    TypeName("Trapezoid");
 
 
     // Constructors
 
         //- Construct from ODE
-        RK(const ODESystem& ode);
+        Trapezoid(const ODESystem& ode, const dictionary& dict);
 
 
     // Member Functions
 
-        void solve
+        //- Solve a single step dx and return the error
+        scalar solve
         (
             const ODESystem& ode,
-            const scalar x,
+            const scalar x0,
-            const scalarField& y,
+            const scalarField& y0,
-            const scalarField& dydx,
+            const scalarField& dydx0,
-            const scalar h,
+            const scalar dx,
-            scalarField& yout,
+            scalarField& y
-            scalarField& yerr
         ) const;
 
-
+        //- Solve the ODE system and the update the state
         void solve
         (
             const ODESystem& ode,
             scalar& x,
             scalarField& y,
-            scalarField& dydx,
+            scalar& dxTry
-            const scalar eps,
-            const scalarField& yScale,
-            const scalar hTry,
-            scalar& hDid,
-            scalar& hNext
         ) const;
 };
 

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "adaptiveSolver.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+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  * * * * * * * * * * * * * //
+
+void Foam::adaptiveSolver::solve
+(
+    const ODESystem& odes,
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) 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, dydx0_, dx, yTemp_);
+
+        // If error is large reduce dx
+        if (err > 1)
+        {
+            scalar scale = max(safeScale_*pow(err, -alphaDec_), minScale_);
+            dx *= scale;
+
+            if (dx < VSMALL)
+            {
+                FatalErrorIn("adaptiveSolver::solve")
+                    << "stepsize underflow"
+                    << exit(FatalError);
+            }
+        }
+    } while (err > 1);
+
+    // Update the state
+    x += dx;
+    y = yTemp_;
+
+    // If the error is small increase the step-size
+    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

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::adaptiveSolver
+
+Description
+
+SourceFiles
+    adaptiveSolver.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef adaptiveSolver_H
+#define adaptiveSolver_H
+
+#include "ODESolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class adaptiveSolver Declaration
+\*---------------------------------------------------------------------------*/
+
+class adaptiveSolver
+{
+    // Private data
+
+        //- Step-size adjustment controls
+        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_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from ODESystem
+        adaptiveSolver(const ODESystem& ode, const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~adaptiveSolver()
+    {}
+
+
+    // Member Functions
+
+        //- Solve a single step dx and return the error
+        virtual scalar solve
+        (
+            const ODESystem& ode,
+            const scalar x0,
+            const scalarField& y0,
+            const scalarField& dydx0,
+            const scalar dx,
+            scalarField& y
+        ) const = 0;
+
+        //- 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/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/ODE/ODESystem/ODESystem.H

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

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/extrudePatchMesh/extrudePatchMesh.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
@@ -47,14 +47,15 @@ extrudePatchMesh::extrudePatchMesh
 (
     const fvMesh& mesh,
     const fvPatch& patch,
-    const dictionary& dict
+    const dictionary& dict,
+    const word regionName
 )
 :
     fvMesh
     (
         IOobject
         (
-            dict.lookup("region"),
+            regionName,
             mesh.facesInstance(),
             mesh,
             IOobject::READ_IF_PRESENT,

src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.H

@@ -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
@@ -25,12 +25,11 @@ Class
     extrudePatchMesh
 
 Description
-    Mesh at a patch created on the fly. The following entried should be used
+    Mesh at a patch created on the fly. The following entry should be used
     on the field boundary dictionary:
 
         // New Shell mesh data
 
-        region          "regionMesh";
         extrudeModel    linearNormal;
         linearNormalCoeffs
         {
@@ -120,7 +119,8 @@ public:
         (
             const fvMesh&,
             const fvPatch&,
-            const dictionary&
+            const dictionary&,
+            const word
         );
 
 

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];
+
+        forAll(pz, i)
+        {
+            label point0 = pz[i];
+            label allPointI = from0ToAllPoints[point0];
+
+            if (pointToZone[allPointI] == -1)
+            {
+                pointToZone[allPointI] = zoneI;
+            }
+            else if (pointToZone[allPointI] != zoneI)
+            {
+                labelList& pZones = addPointToZones[allPointI];
+                if (findIndex(pZones, zoneI) == -1)
+                {
+                    pZones.append(zoneI);
+                }
+            }
+        }
     }
 
-    // Get sorted zone contents for duplicate element recognition
+    // mesh1 zones renumbered
-    PtrList<SortableList<label> > pzPointsSorted(pzPoints.size());
-    forAll(pzPoints, zoneI)
-    {
-        pzPointsSorted.set
-        (
-            zoneI,
-            new SortableList<label>(pzPoints[zoneI])
-        );
-    }
-
-    // Now we have full addressing for points so do the pointZones of mesh1.
     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,85 +1037,166 @@ 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)
+                {
+                    flip0 = !flip0;
+                }
+
+                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);
+                    }
+                }
             }
         }
     }
 
-    // Get sorted zone contents for duplicate element recognition
+    // mesh1 zones renumbered
-    PtrList<SortableList<label> > fzFacesSorted(fzFaces.size());
-    forAll(fzFaces, zoneI)
-    {
-        fzFacesSorted.set
-        (
-            zoneI,
-            new SortableList<label>(fzFaces[zoneI])
-        );
-    }
-
-    // Now we have full addressing for faces so do the faceZones of mesh1.
     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]);
+        }
+    }
+
     forAll(fzFaces, i)
     {
         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,21 +123,32 @@ protected:
             //- From resulting back to original field
             HashTable<word> reverseFieldMap_;
 
-        //- Locally constructed fields
+
-        PtrList<volScalarField> vsf_;
+        // Calculated addressing
-        PtrList<volVectorField> vvf_;
+
-        PtrList<volSphericalTensorField> vSpheretf_;
+            //- From cell to seed patch faces
-        PtrList<volSymmTensorField> vSymmtf_;
+            labelListList cellToWalls_;
-        PtrList<volTensorField> vtf_;
+
+            //- 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_;
+            PtrList<volSymmTensorField> vSymmtf_;
+            PtrList<volTensorField> vtf_;
 
 
     // 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,35 +62,8 @@ 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 "
-                forAllConstIter(labelHashSet, patchSet_, iter)
-                {
-                    label patchI = iter.key();
-
-                    sampleFld.boundaryField().set
-                    (
-                        patchI,
-                        new mappedFieldFvPatchField<Type>
-                        (
-                            sampleFld.mesh().boundary()[patchI],
-                            sampleFld.dimensionedInternalField(),
-
-                            sampleFld.mesh().name(),
-                            mappedPatchBase::NEARESTCELL,
-                            word::null,     // samplePatch
-                            -distance_,
-
-                            sampleFld.name(),       // fieldName
-                            false,                  // setAverage
-                            pTraits<Type>::zero,    // average
-                            interpolationCellPoint<Type>::typeName
-                        )
-                    );
-                }
-
-                Info<< "    created " << sampleFld.name() << " to sample "
                     << fld.name() << endl;
             }
         }
@@ -100,6 +71,53 @@ void Foam::nearWallFields::createFields
 }
 
 
+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();
+
+        fvPatchField<Type>& pfld = fld.boundaryField()[patchI];
+
+        Field<Type> newFld(pfld.size());
+        forAll(pfld, i)
+        {
+            newFld[i] = sampledValues[nPatchFaces++];
+        }
+
+        pfld == newFld;
+    }
+}
+
+
 template<class Type>
 void Foam::nearWallFields::sampleFields
 (
@@ -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/regionModels/thermalBaffleModels/Make/options

@@ -7,7 +7,9 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/mesh/extrudeModel/lnInclude
 
 LIB_LIBS = \
     -lregionModels \
@@ -16,4 +18,5 @@ LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lOpenFOAM \
-    -lradiationModels
+    -lradiationModels \
+    -ldynamicMesh

src/regionModels/thermalBaffleModels/derivedFvPatchFields/thermalBaffle/thermalBaffleFvPatchScalarField.C

@@ -45,7 +45,8 @@ thermalBaffleFvPatchScalarField
     turbulentTemperatureRadCoupledMixedFvPatchScalarField(p, iF),
     owner_(false),
     baffle_(),
-    dict_(dictionary::null)
+    dict_(dictionary::null),
+    extrudeMeshPtr_()
 {}
 
 
@@ -66,8 +67,9 @@ thermalBaffleFvPatchScalarField
         mapper
     ),
     owner_(ptf.owner_),
-    baffle_(ptf.baffle_),
+    baffle_(),
-    dict_(ptf.dict_)
+    dict_(ptf.dict_),
+    extrudeMeshPtr_()
 {}
 
 
@@ -82,47 +84,47 @@ thermalBaffleFvPatchScalarField
     turbulentTemperatureRadCoupledMixedFvPatchScalarField(p, iF, dict),
     owner_(false),
     baffle_(),
-    dict_(dict)
+    dict_(dict),
+    extrudeMeshPtr_()
 {
-    if (!isA<mappedPatchBase>(patch().patch()))
-    {
-        FatalErrorIn
-        (
-            "thermalBaffleFvPatchScalarField::"
-            "thermalBaffleFvPatchScalarField\n"
-            "(\n"
-            "    const fvPatch& p,\n"
-            "    const DimensionedField<scalar, volMesh>& iF,\n"
-            "    const dictionary& dict\n"
-            ")\n"
-        )   << "\n    patch type '" << patch().type()
-            << "' not type '" << mappedPatchBase::typeName << "'"
-            << "\n    for patch " << patch().name()
-            << " of field " << dimensionedInternalField().name()
-            << " in file " << dimensionedInternalField().objectPath()
-            << exit(FatalError);
-    }
-
-    const mappedPatchBase& mpp =
-        refCast<const mappedPatchBase>(patch().patch());
-
-    const word nbrMesh = mpp.sampleRegion();
 
     const fvMesh& thisMesh = patch().boundaryMesh().mesh();
 
     typedef regionModels::thermalBaffleModels::thermalBaffleModel baffle;
 
-    if
+    if (thisMesh.name() == polyMesh::defaultRegion)
-    (
-        thisMesh.name() == polyMesh::defaultRegion
-     && !thisMesh.foundObject<baffle>(nbrMesh)
-     && !owner_
-    )
     {
-        Info << "Creating thermal baffle" <<  nbrMesh << endl;
+        const word regionName =
-        baffle_.reset(baffle::New(thisMesh, dict).ptr());
+            dict_.lookupOrDefault<word>("regionName", "none");
-        owner_ = true;
+
-        baffle_->rename(nbrMesh);
+        const word baffleName("3DBaffle" + regionName);
+
+        if
+        (
+            !thisMesh.time().foundObject<fvMesh>(regionName)
+         && regionName != "none"
+        )
+        {
+            if (extrudeMeshPtr_.empty())
+            {
+                createPatchMesh();
+            }
+
+            baffle_.reset(baffle::New(thisMesh, dict).ptr());
+            owner_ = true;
+            baffle_->rename(baffleName);
+        }
+        else if //Backwards compatibility (if region exists)
+        (
+            thisMesh.time().foundObject<fvMesh>(regionName)
+         && baffle_.empty()
+         && regionName != "none"
+        )
+        {
+            baffle_.reset(baffle::New(thisMesh, dict).ptr());
+            owner_ = true;
+            baffle_->rename(baffleName);
+        }
     }
 }
 
@@ -136,8 +138,9 @@ thermalBaffleFvPatchScalarField
 :
     turbulentTemperatureRadCoupledMixedFvPatchScalarField(ptf, iF),
     owner_(ptf.owner_),
-    baffle_(ptf.baffle_),
+    baffle_(),
-    dict_(ptf.dict_)
+    dict_(ptf.dict_),
+    extrudeMeshPtr_()
 {}
 
 
@@ -163,6 +166,36 @@ void thermalBaffleFvPatchScalarField::rmap
 }
 
 
+void thermalBaffleFvPatchScalarField::createPatchMesh()
+{
+    const fvMesh& defaultRegion =
+        db().time().lookupObject<fvMesh>(fvMesh::defaultRegion);
+
+    word regionName = dict_.lookup("regionName");
+
+    extrudeMeshPtr_.reset
+    (
+        new extrudePatchMesh
+        (
+            defaultRegion,
+            patch(),
+            dict_,
+            regionName
+        )
+    );
+
+    if (extrudeMeshPtr_.empty())
+    {
+        WarningIn
+        (
+            "thermalBaffleFvPatchScalarField::createPatchMesh()\n"
+        )   << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
+            << " patchMeshPtr not set."
+            << endl;
+    }
+}
+
+
 void thermalBaffleFvPatchScalarField::updateCoeffs()
 {
     if (this->updated())
@@ -189,27 +222,35 @@ void thermalBaffleFvPatchScalarField::write(Ostream& os) const
 
     if (thisMesh.name() == polyMesh::defaultRegion && owner_)
     {
-        word thermoModel = dict_.lookup("thermalBaffleModel");
 
-        os.writeKeyword("thermalBaffleModel")
+        os.writeKeyword("extrudeModel");
-            << thermoModel
+        os << word(dict_.lookup("extrudeModel"))
-            << token::END_STATEMENT << nl;
+           << token::END_STATEMENT << nl;
+
+        os.writeKeyword("nLayers");
+        os << readLabel(dict_.lookup("nLayers"))
+           << token::END_STATEMENT << nl;
+
+        os.writeKeyword("expansionRatio");
+        os << readScalar(dict_.lookup("expansionRatio"))
+           << token::END_STATEMENT << nl;
+
+        os.writeKeyword("columnCells");
+        os << readBool(dict_.lookup("columnCells"))
+           << token::END_STATEMENT << nl;
+
+        word extrudeModel(word(dict_.lookup("extrudeModel")) + "Coeffs");
+        os.writeKeyword(extrudeModel);
+        os << dict_.subDict(extrudeModel) << nl;
 
         word regionName = dict_.lookup("regionName");
         os.writeKeyword("regionName") << regionName
             << token::END_STATEMENT << nl;
 
-        bool infoOutput = readBool(dict_.lookup("infoOutput"));
-        os.writeKeyword("infoOutput") << infoOutput
-            << token::END_STATEMENT << nl;
-
         bool active = readBool(dict_.lookup("active"));
         os.writeKeyword("active") <<  active
             << token::END_STATEMENT << nl;
 
-        os.writeKeyword(word(thermoModel + "Coeffs"));
-        os << dict_.subDict(thermoModel + "Coeffs") << nl;
-
         os.writeKeyword("thermoType");
         os << dict_.subDict("thermoType") << nl;
 

src/regionModels/thermalBaffleModels/derivedFvPatchFields/thermalBaffle/thermalBaffleFvPatchScalarField.H

@@ -42,78 +42,111 @@ Description
     dictionary entries.
 
     \heading Patch usage
-
     Example of the boundary condition specification:
     \verbatim
     myPatch
     {
         type                compressible::temperatureThermoBaffle;
 
-        // Coupled boundary condition
+        // Underlaying coupled boundary condition
-        Tnbr                T;
+            Tnbr               T;
-        kappa               fluidThermo; // or solidThermo
+            kappa              fluidThermo; // or solidThermo
-        KappaName           none;
+            KappaName          none;
-        QrNbr               Qr; // or none.Name of Qr field on neighbour region
+            QrNbr              Qr;//or none.Name of Qr field on neighbourregion
-        Qr                  Qr; // or none.Name of Qr field on local region
+            Qr                 none;// or none.Name of Qr field on localregion
-
+            value              uniform 300;
 
         // Thermo baffle model
-        thermalBaffleModel   thermalBaffle;
+            regionName          baffleRegion; // solid region name
-        regionName          baffleRegion;
+            infoOutput          yes;
-        infoOutput          yes;
+            active              yes;
-        active              yes;
-        thermalBaffleCoeffs
-        {
-        }
 
-
+        // Solid thermo in solid region
-        // Solid thermo
+            thermoType
-        thermoType
-        {
-            type            heSolidThermo;
-            mixture         pureSolidMixture;
-            transport       constIso;
-            thermo          hConst;
-            equationOfState rhoConst;
-            specie          specie;
-            energy          sensibleEnthalpy;
-        }
-
-        mixture
-        {
-            specie
             {
-                nMoles          1;
+                type            heSolidThermo;
-                molWeight       20;
+                mixture         pureSolidMixture;
+                transport       constIso;
+                thermo          hConst;
+                equationOfState rhoConst;
+                specie          specie;
+                energy          sensibleEnthalpy;
             }
-            transport
-            {
-                kappa           0.01;
-            }
-            thermodynamics
-            {
-                Hf              0;
-                Cp              15;
-            }
-            density
-            {
-                rho             80;
-            }
-        }
 
-        radiation
+            mixture
-        {
+            {
-            radiationModel  opaqueSolid;
+                specie
-            absorptionEmissionModel none;
+                {
-            scatterModel    none;
+                    nMoles          1;
-        }
+                    molWeight       20;
+                }
+                transport
+                {
+                    kappa           0.01;
+                }
+                thermodynamics
+                {
+                    Hf              0;
+                    Cp              15;
+                }
+                density
+                {
+                    rho             80;
+                }
+            }
 
-        value               uniform 300;
+            radiation
+            {
+                radiationModel  opaqueSolid;
+                absorptionEmissionModel none;
+                scatterModel    none;
+            }
+
+
+        // Extrude model for new region
+
+            extrudeModel        linearNormal;
+            nLayers             50;
+            expansionRatio      1;
+            columnCells         false; //3D or 1D
+            linearNormalCoeffs
+            {
+                thickness           0.02;
+            }
+
+
+        // New mesh polyPatch information
+
+            bottomCoeffs
+            {
+                name        "bottom";
+                type        mappedWall;
+                sampleMode  nearestPatchFace;
+                samplePatch baffle3DWall_master;
+                offsetMode  uniform;
+                offset      (0 0 0);
+            }
+
+            topCoeffs
+            {
+                name        "top";
+                type        mappedWall;
+                sampleMode  nearestPatchFace;
+                samplePatch baffle3DWall_slave;
+                offsetMode  uniform;
+                offset      (0 0 0);
+            }
+
+            sideCoeffs
+            {
+                name        "side";
+                type        patch;
+            }
     }
     \endverbatim
 
 SeeAlso
-    Foam::turbulentTemperatureRadCoupledMixedFvPatchScalarField
+    Foam::turbulentTemperatureCoupledBaffleMixedFvPatchScalarField
     Foam::regionModels::thermalBaffleModels::thermalBaffleModel
 
 SourceFiles
@@ -128,6 +161,7 @@ SourceFiles
 #include "autoPtr.H"
 #include "regionModel.H"
 #include "thermalBaffleModel.H"
+#include "extrudePatchMesh.H"
 #include "turbulentTemperatureRadCoupledMixedFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -157,6 +191,14 @@ class thermalBaffleFvPatchScalarField
         //- Dictionary
         dictionary dict_;
 
+        //- Auto pointer to extrapolated mesh from patch
+        autoPtr<extrudePatchMesh> extrudeMeshPtr_;
+
+
+    // Private member functions
+
+        //- Extrude mesh
+        void createPatchMesh();
 
 public:
 

src/regionModels/thermalBaffleModels/thermalBaffleModel/thermalBaffleModelNew.C

@@ -53,11 +53,14 @@ autoPtr<thermalBaffleModel> thermalBaffleModel::New(const fvMesh& mesh)
             )
         );
 
-        thermalBafflePropertiesDict.lookup("thermalBaffleModel") >> modelType;
+        word modelType =
+            thermalBafflePropertiesDict.lookupOrDefault<word>
+            (
+                "thermalBaffleModel",
+                "thermalBaffle"
+            );
     }
 
-    Info<< "Selecting baffle model " << modelType << endl;
-
     meshConstructorTable::iterator cstrIter =
         meshConstructorTablePtr_->find(modelType);
 
@@ -82,9 +85,8 @@ autoPtr<thermalBaffleModel> thermalBaffleModel::New
     const dictionary& dict
 )
 {
-    word modelType = dict.lookup("thermalBaffleModel");
+    word modelType =
-
+        dict.lookupOrDefault<word>("thermalBaffleModel", "thermalBaffle");
-    Info<< "Selecting baffle model " << modelType << endl;
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(modelType);

src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C

@@ -310,10 +310,7 @@ void Foam::chemistryModel<CompType, ThermoType>::derivatives
     }
     dT /= rho*cp;
 
-    // limit the time-derivative, this is more stable for the ODE
+    dcdt[nSpecie_] = -dT;
-    // solver when calculating the allowed time step
-    const scalar dTLimited = min(500.0, mag(dT));
-    dcdt[nSpecie_] = -dT*dTLimited/(mag(dT) + 1.0e-10);
 
     // dp/dt = ...
     dcdt[nSpecie_ + 1] = 0.0;
@@ -743,7 +740,7 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
             RR_[i][celli] = (c[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
         }
     }
-
+    Info << "deltaTMin " << deltaTMin  << endl;
     return deltaTMin;
 }
 

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

src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.C

@@ -46,13 +46,17 @@ thermalBaffle1DFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
+    mappedPatchBase(p.patch()),
     mixedFvPatchScalarField(p, iF),
     TName_("T"),
     baffleActivated_(true),
     thickness_(p.size()),
     Qs_(p.size()),
     solidDict_(),
-    solidPtr_(NULL)
+    solidPtr_(NULL),
+    QrPrevious_(p.size()),
+    QrRelaxation_(0),
+    QrName_("undefined-Qr")
 {}
 
 
@@ -66,13 +70,17 @@ thermalBaffle1DFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
+    mappedPatchBase(p.patch(), ptf),
     mixedFvPatchScalarField(ptf, p, iF, mapper),
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
     Qs_(ptf.Qs_),
     solidDict_(ptf.solidDict_),
-    solidPtr_(ptf.solidPtr_)
+    solidPtr_(ptf.solidPtr_),
+    QrPrevious_(ptf.QrPrevious_),
+    QrRelaxation_(ptf.QrRelaxation_),
+    QrName_(ptf.QrName_)
 {}
 
 
@@ -85,33 +93,25 @@ thermalBaffle1DFvPatchScalarField
     const dictionary& dict
 )
 :
+    mappedPatchBase
+    (
+        p.patch(),
+        p.boundaryMesh().mesh().name(),
+        NEARESTPATCHFACE,
+        dict.lookup("samplePatch"),
+        0.0
+    ),
     mixedFvPatchScalarField(p, iF),
     TName_("T"),
-    baffleActivated_(readBool(dict.lookup("baffleActivated"))),
+    baffleActivated_(dict.lookupOrDefault<bool>("baffleActivated", true)),
-    thickness_(scalarField("thickness", dict, p.size())),
+    thickness_(),
-    Qs_(scalarField("Qs", dict, p.size())),
+    Qs_(),
     solidDict_(dict),
-    solidPtr_(new solidType(dict))
+    solidPtr_(),
+    QrPrevious_(p.size(), 0.0),
+    QrRelaxation_(dict.lookupOrDefault<scalar>("relaxation", 0)),
+    QrName_(dict.lookupOrDefault<word>("Qr", "none"))
 {
-    if (!isA<mappedPatchBase>(this->patch().patch()))
-    {
-        FatalErrorIn
-        (
-            "thermalBaffle1DFvPatchScalarField::"
-            "thermalBaffle1DFvPatchScalarField"
-            "("
-                "const fvPatch&,\n"
-                "const DimensionedField<scalar, volMesh>&, "
-                "const dictionary&"
-            ")"
-        )   << "\n    patch type '" << patch().type()
-            << "' not type '" << mappedPatchBase::typeName << "'"
-            << "\n    for patch " << patch().name()
-            << " of field " << dimensionedInternalField().name()
-            << " in file " << dimensionedInternalField().objectPath()
-            << exit(FatalError);
-    }
-
     fvPatchScalarField::operator=(scalarField("value", dict, p.size()));
 
     if (dict.found("refValue") && baffleActivated_)
@@ -139,13 +139,17 @@ thermalBaffle1DFvPatchScalarField
     const thermalBaffle1DFvPatchScalarField& ptf
 )
 :
+    mappedPatchBase(ptf.patch().patch(), ptf),
     mixedFvPatchScalarField(ptf),
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
     Qs_(ptf.Qs_),
     solidDict_(ptf.solidDict_),
-    solidPtr_(ptf.solidPtr_)
+    solidPtr_(ptf.solidPtr_),
+    QrPrevious_(ptf.QrPrevious_),
+    QrRelaxation_(ptf.QrRelaxation_),
+    QrName_(ptf.QrName_)
 {}
 
 
@@ -157,29 +161,118 @@ thermalBaffle1DFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
+    mappedPatchBase(ptf.patch().patch(), ptf),
     mixedFvPatchScalarField(ptf, iF),
     TName_(ptf.TName_),
     baffleActivated_(ptf.baffleActivated_),
     thickness_(ptf.thickness_),
     Qs_(ptf.Qs_),
     solidDict_(ptf.solidDict_),
-    solidPtr_(ptf.solidPtr_)
+    solidPtr_(ptf.solidPtr_),
+    QrPrevious_(ptf.QrPrevious_),
+    QrRelaxation_(ptf.QrRelaxation_),
+    QrName_(ptf.QrName_)
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class solidType>
-const solidType& thermalBaffle1DFvPatchScalarField<solidType>::solidPtr() const
+bool thermalBaffle1DFvPatchScalarField<solidType>::owner() const
 {
-    if (!solidPtr_.empty())
+    const label patchi = patch().index();
+
+    const label nbrPatchi = samplePolyPatch().index();
+
+    return (patchi < nbrPatchi);
+}
+
+
+template<class solidType>
+const solidType& thermalBaffle1DFvPatchScalarField<solidType>::solid() const
+{
+    if (this->owner())
     {
+        if (solidPtr_.empty())
+        {
+            solidPtr_.reset(new solidType(solidDict_));
+        }
         return solidPtr_();
     }
     else
     {
-        solidPtr_.reset(new solidType(solidDict_));
+        const fvPatch& nbrPatch =
-        return solidPtr_();
+            patch().boundaryMesh()[samplePolyPatch().index()];
+
+        const thermalBaffle1DFvPatchScalarField& nbrField =
+        refCast<const thermalBaffle1DFvPatchScalarField>
+        (
+            nbrPatch.template lookupPatchField<volScalarField, scalar>(TName_)
+        );
+
+        return nbrField.solid();
+    }
+}
+
+
+template<class solidType>
+const scalarField& thermalBaffle1DFvPatchScalarField<solidType>::
+baffleThickness() const
+{
+    if (this->owner())
+    {
+        if (thickness_.size() > 0)
+        {
+            return thickness_;
+        }
+        else
+        {
+            thickness_ = scalarField("thickness", solidDict_, patch().size());
+            return thickness_;
+        }
+    }
+    else
+    {
+        const fvPatch& nbrPatch =
+            patch().boundaryMesh()[samplePolyPatch().index()];
+        const thermalBaffle1DFvPatchScalarField& nbrField =
+        refCast<const thermalBaffle1DFvPatchScalarField>
+        (
+            nbrPatch.template lookupPatchField<volScalarField, scalar>(TName_)
+        );
+
+        return nbrField.thickness_;
+    }
+}
+
+
+template<class solidType>
+const scalarField& thermalBaffle1DFvPatchScalarField<solidType>::Qs() const
+{
+    if (this->owner())
+    {
+        if (Qs_.size() > 0)
+        {
+            return Qs_;
+        }
+        else
+        {
+            Qs_ = scalarField("Qs", solidDict_, patch().size());
+            return Qs_;
+        }
+    }
+    else
+    {
+        const fvPatch& nbrPatch =
+            patch().boundaryMesh()[samplePolyPatch().index()];
+
+        const thermalBaffle1DFvPatchScalarField& nbrField =
+        refCast<const thermalBaffle1DFvPatchScalarField>
+        (
+            nbrPatch.template lookupPatchField<volScalarField, scalar>(TName_)
+        );
+
+        return nbrField.Qs_;
     }
 }
 
@@ -219,18 +312,14 @@ void thermalBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
     {
         return;
     }
-
     // Since we're inside initEvaluate/evaluate there might be processor
     // comms underway. Change the tag we use.
     int oldTag = UPstream::msgType();
     UPstream::msgType() = oldTag+1;
 
-    const mappedPatchBase& mpp =
-        refCast<const mappedPatchBase>(patch().patch());
-
     const label patchi = patch().index();
 
-    const label nbrPatchi = mpp.samplePolyPatch().index();
+    const label nbrPatchi = samplePolyPatch().index();
 
     if (baffleActivated_)
     {
@@ -243,108 +332,49 @@ void thermalBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
             );
 
         // local properties
-
         const scalarField kappaw(turbModel.kappaEff(patchi));
 
         const fvPatchScalarField& Tp =
             patch().template lookupPatchField<volScalarField, scalar>(TName_);
 
-        const scalarField qDot(kappaw*Tp.snGrad());
+
+        scalarField Qr(Tp.size(), 0.0);
+
+        if (QrName_ != "none")
+        {
+            Qr = patch().template lookupPatchField<volScalarField, scalar>
+                (QrName_);
+
+            Qr = QrRelaxation_*Qr + (1.0 - QrRelaxation_)*QrPrevious_;
+            QrPrevious_ = Qr;
+        }
 
         tmp<scalarField> Ti = patchInternalField();
 
-        scalarField myh(patch().deltaCoeffs()*kappaw);
+        scalarField myKDelta(patch().deltaCoeffs()*kappaw);
 
-        // nbr properties
+        // nrb properties
-
+        const fvPatchScalarField& nbrTp =
-        const scalarField nbrKappaw(turbModel.kappaEff(nbrPatchi));
-
-        const fvPatchScalarField& nbrTw =
             turbModel.thermo().T().boundaryField()[nbrPatchi];
 
-        scalarField nbrQDot(nbrKappaw*nbrTw.snGrad());
+        // solid properties
-        mpp.map().distribute(nbrQDot);
+        scalarField kappas(patch().size(), 0.0);
-
+        forAll(kappas, i)
-        const thermalBaffle1DFvPatchScalarField& nbrField =
-        refCast<const thermalBaffle1DFvPatchScalarField>
-        (
-            nbrPatch.template lookupPatchField<volScalarField, scalar>(TName_)
-        );
-
-        scalarField nbrTi(nbrField.patchInternalField());
-        mpp.map().distribute(nbrTi);
-
-        scalarField nbrTp =
-           nbrPatch.template lookupPatchField<volScalarField, scalar>(TName_);
-        mpp.map().distribute(nbrTp);
-
-        scalarField nbrh(nbrPatch.deltaCoeffs()*nbrKappaw);
-        mpp.map().distribute(nbrh);
-
-
-        // heat source
-        const scalarField Q(Qs_/thickness_);
-
-        tmp<scalarField> tKDeltaw(new scalarField(patch().size()));
-        scalarField KDeltaw = tKDeltaw();
-
-        // Create fields for solid properties (p paramater not used)
-        forAll(KDeltaw, i)
         {
-            KDeltaw[i] =
+            kappas[i] = solid().kappa(0.0, (Tp[i] + nbrTp[i])/2.0);
-                solidPtr().kappa(0.0, (Tp[i] + nbrTw[i])/2.0)/thickness_[i];
         }
 
-        const scalarField q
+        scalarField KDeltaSolid(kappas/baffleThickness());
-        (
-            (Ti() - nbrTi)/(1.0/KDeltaw + 1.0/nbrh + 1.0/myh)
-        );
 
-        forAll(qDot, i)
+        scalarField alpha(KDeltaSolid - Qr/Tp);
-        {
-            if (Qs_[i] == 0)
-            {
-                this->refValue()[i] = Ti()[i] - q[i]/myh[i];
-                this->refGrad()[i] = 0.0;
-                this->valueFraction()[i] = 1.0;
-            }
-            else
-            {
-                if (q[i] > 0)
-                {
-                    this->refValue()[i] =
-                        nbrTp[i]
-                      - Q[i]*thickness_[i]/(2*KDeltaw[i]);
 
-                    this->refGrad()[i] = 0.0;
+        valueFraction() = alpha/(alpha + myKDelta);
-                    this->valueFraction()[i] =
+
-                        1.0
+        refValue() = (KDeltaSolid*nbrTp + Qs()/2.0)/alpha;
-                        /
-                        (
-                            1.0
-                          + patch().deltaCoeffs()[i]*kappaw[i]/KDeltaw[i]
-                        );
-                }
-                else if (q[i] < 0)
-                {
-                    this->refValue()[i] = 0.0;
-                    this->refGrad()[i] =
-                          (-nbrQDot[i] + Q[i]*thickness_[i])/kappaw[i];
-                    this->valueFraction()[i] = 0.0;
-                }
-                else
-                {
-                    scalar Qt = Q[i]*thickness_[i];
-                    this->refValue()[i] = 0.0;
-                    this->refGrad()[i] = Qt/2/kappaw[i];
-                    this->valueFraction()[i] = 0.0;
-                }
-            }
-        }
 
         if (debug)
         {
-            scalar Q = gSum(patch().magSf()*qDot);
+            scalar Q = gAverage(kappaw*snGrad());
             Info<< patch().boundaryMesh().mesh().name() << ':'
                 << patch().name() << ':'
                 << this->dimensionedInternalField().name() << " <- "
@@ -366,16 +396,21 @@ void thermalBaffle1DFvPatchScalarField<solidType>::updateCoeffs()
 }
 
 template<class solidType>
-void thermalBaffle1DFvPatchScalarField<solidType>:: write(Ostream& os) const
+void thermalBaffle1DFvPatchScalarField<solidType>::write(Ostream& os) const
 {
     mixedFvPatchScalarField::write(os);
-    os.writeKeyword("TName")
+    mappedPatchBase::write(os);
-        << TName_ << token::END_STATEMENT << nl;
+
-    thickness_.writeEntry("thickness", os);
+    if (this->owner())
-    os.writeKeyword("baffleActivated")
+    {
-        << baffleActivated_ << token::END_STATEMENT << nl;
+        baffleThickness().writeEntry("thickness", os);
-    Qs_.writeEntry("Qs", os);
+        Qs().writeEntry("Qs", os);
-    solidPtr().write(os);
+        solid().write(os);
+    }
+
+    os.writeKeyword("Qr")<< QrName_ << token::END_STATEMENT << nl;
+    os.writeKeyword("QrRelaxation")<< QrRelaxation_
+        << token::END_STATEMENT << nl;
 }
 
 

src/turbulenceModels/compressible/turbulenceModel/derivedFvPatchFields/thermalBaffle1D/thermalBaffle1DFvPatchScalarField.H

@@ -24,9 +24,68 @@ License
 Class
     Foam::thermalBaffle1DFvPatchScalarField
 
+Group
+    grpThermoBoundaryConditions
+
 Description
-    Boundary which solves the 1D steady state heat transfer equation
+
-    through a baffle.
+    This BC solves a steady 1D thermal baffle. The solid properties are
+    specify as dictionary. Optionaly radiative heat flux (Qr) can be
+    incorporated into the balance. Some under-relaxation might be needed on
+    Qr.
+    Baffle and solid properties need to be specified on the master side
+    of the baffle.
+
+    \heading Patch usage
+    Example of the boundary condition specification using constant
+    solid thermo :
+
+    \verbatim
+    myPatch_master
+    {
+        type   compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
+        samplePatch     myPatch_slave;
+
+        thickness       uniform 0.005;  // thickness [m]
+        Qs              uniform 100;    // heat flux [W/m2]
+
+        Qr              none;
+        relaxation      0;
+
+        // Solid thermo
+        specie
+        {
+            nMoles          1;
+            molWeight       20;
+        }
+        transport
+        {
+            kappa           1;
+        }
+        thermodynamics
+        {
+            Hf              0;
+            Cp              10;
+        }
+        equationOfState
+        {
+            rho             10;
+        }
+
+        value               uniform 300;
+    }
+
+    myPatch_slave
+    {
+        type   compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
+        samplePatch     myPatch_master_master;
+
+        Qr              none;
+        relaxation      0;
+    }
+    \endverbatim
+
+
 
 SourceFiles
     thermalBaffle1DFvPatchScalarField.C
@@ -38,7 +97,7 @@ SourceFiles
 
 #include "mixedFvPatchFields.H"
 #include "autoPtr.H"
-
+#include "mappedPatchBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,6 +113,7 @@ namespace compressible
 template<class solidType>
 class thermalBaffle1DFvPatchScalarField
 :
+    public mappedPatchBase,
     public mixedFvPatchScalarField
 {
     // Private data
@@ -65,10 +125,10 @@ class thermalBaffle1DFvPatchScalarField
         bool baffleActivated_;
 
         //- Baffle thickness [m]
-        scalarField thickness_;
+        mutable scalarField thickness_;
 
         //- Superficial heat source [W/m2]
-        scalarField Qs_;
+        mutable scalarField Qs_;
 
         //- Solid dictionary
         dictionary  solidDict_;
@@ -76,11 +136,29 @@ class thermalBaffle1DFvPatchScalarField
         //- Solid thermo
         mutable autoPtr<solidType> solidPtr_;
 
+        //- Chache Qr for relaxation
+        scalarField QrPrevious_;
+
+        //- Relaxation for Qr
+        scalar QrRelaxation_;
+
+        //- Name of the radiative heat flux in local region
+        const word QrName_;
+
 
     // Private members
 
-        //- Return non const solid thermo autoMap
+        //- Return const solid thermo
-        const solidType& solidPtr() const;
+        const solidType& solid() const;
+
+        //- Return Qs from master
+        const scalarField& Qs() const;
+
+        //- Return thickness from master
+        const scalarField& baffleThickness() const;
+
+        //- Is Owner
+        bool owner() const;
 
 
 public:
@@ -153,7 +231,6 @@ public:
 
     // Member functions
 
-
         // Mapping functions
 
             //- Map (and resize as needed) from self given a mapping object
@@ -170,8 +247,6 @@ public:
             );
 
 
-
-
         //- Update the coefficients associated with the patch field
         virtual void updateCoeffs();
 

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;
 }
 
 // ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/T

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 0 0 1 0 0 0 ];
+dimensions      [0 0 0 1 0 0 0];
 
 internalField   uniform 300;
 

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/U

@@ -46,6 +46,7 @@ boundaryField
     {
         type            empty;
     }
+
 }
 
 

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/alphat

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 1 -1 -1 0 0 0 0 ];
+dimensions      [1 -1 -1 0 0 0 0];
 
 internalField   uniform 0;
 
@@ -24,11 +24,13 @@ boundaryField
     floor
     {
         type            compressible::alphatWallFunction;
+        Prt             0.85;
         value           uniform 0;
     }
     ceiling
     {
         type            compressible::alphatWallFunction;
+        Prt             0.85;
         value           uniform 0;
     }
     inlet

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/baffle3DRegion/T

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [ 0 0 0 1 0 0 0 ];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    bottom
+    {
+        type                compressible::thermalBaffle;
+        Tnbr                T;
+        kappa               solidThermo;
+        kappaName           none;
+        QrNbr               none;
+        Qr                  none;
+        value               uniform 300;
+    }
+    side
+    {
+        type                zeroGradient;
+    }
+    top
+    {
+        type                compressible::thermalBaffle;
+        Tnbr                T;
+        kappa               solidThermo;
+        kappaName           none;
+        QrNbr               none;
+        Qr                  none;
+        value               uniform 300;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/epsilon

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -3 0 0 0 0 ];
+dimensions      [0 2 -3 0 0 0 0];
 
 internalField   uniform 0.01;
 
@@ -24,11 +24,17 @@ boundaryField
     floor
     {
         type            compressible::epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
         value           uniform 0.01;
     }
     ceiling
     {
         type            compressible::epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
         value           uniform 0.01;
     }
     inlet

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/1DBaffle/1DTemperatureMasterBafflePatches

@@ -5,27 +5,18 @@
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
-FoamFile
+
+T
 {
-    version     2.0;
+    type   compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
-    format      ascii;
+    samplePatch     baffle1DWall_slave;
-    class       dictionary;
-    location    "constant";
-    object      thermalBaffleProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermalBaffleModel  none;
+    thickness       uniform 0.005;  // thickness [m]
+    Qs              uniform 100;    // heat flux [W/m2]
 
-active          no;
+    # include "1DbaffleSolidThermo"
 
-regionName      none;
+    value               uniform 300;
-
-thermalBaffleCoeffs
-{
 }
 
-noThermoCoeffs
-{
-}
 // ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/1DBaffle/1DTemperatureSlaveBafflePatches

@@ -0,0 +1,16 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+T
+{
+   type     compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
+   samplePatch  baffle1DWall_master;
+   value        uniform 300;
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/1DBaffle/1DbaffleSolidThermo

@@ -5,23 +5,24 @@
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
-FoamFile
+
+specie
 {
-    version     2.0;
+    nMoles          1;
-    format      ascii;
+    molWeight       20;
-    class       volScalarField;
-    location    "0";
-    object      T;
 }
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+transport
-
-dimensions      [ 0 0 0 1 0 0 0 ];
-
-internalField   uniform 300;
-
-boundaryField
 {
+    kappa           1;
+}
+thermodynamics
+{
+    Hf              0;
+    Cp              10;
+}
+equationOfState
+{
+    rho             10;
 }
-
 
 // ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/3DBaffle/3DTemperatureMasterBafflePatches

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+T
+{
+        type                compressible::thermalBaffle;
+        Tnbr                T;
+        kappa               fluidThermo;
+        kappaName           none;
+        QrNbr               none;
+        Qr                  none;
+        value               uniform 300;
+
+        // Thermo baffle model
+        //thermalBaffleModel   thermalBaffle;
+        regionName           ${baffleRegionName};
+        active               yes;
+
+        # include "3DbaffleSolidThermo"
+
+        // New fvMesh (region) information
+        # include "extrudeModel"
+
+
+        // New mesh polyPatch information
+        bottomCoeffs
+        {
+            name        "bottom";
+            type        mappedWall;
+            sampleMode  nearestPatchFace;
+            samplePatch ${masterPatchName};
+            offsetMode  uniform;
+            offset      (0 0 0);
+        }
+
+        topCoeffs
+        {
+            name        "top";
+            type        mappedWall;
+            sampleMode  nearestPatchFace;
+            samplePatch ${slavePatchName};
+            offsetMode  uniform;
+            offset      (0 0 0);
+        }
+
+        sideCoeffs
+        {
+            name        "side";
+            type        patch;
+        }
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/3DBaffle/3DTemperatureSlaveBafflePatches

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+T
+{
+        type                compressible::thermalBaffle;
+        Tnbr                T;
+        kappa               fluidThermo;
+        kappaName           none;
+        QrNbr               none;
+        Qr                  none;
+        value               uniform 300;
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/3DBaffle/3DbaffleSolidThermo

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+ // Solid thermo
+thermoType
+{
+    type            heSolidThermo;
+    mixture         pureMixture;
+    transport       constIso;
+    thermo          hConst;
+    equationOfState rhoConst;
+    specie          specie;
+    energy          sensibleEnthalpy;
+}
+
+mixture
+{
+    specie
+    {
+        nMoles          1;
+        molWeight       20;
+    }
+    transport
+    {
+        kappa           ${Kappa};
+    }
+    thermodynamics
+    {
+        Hf              0;
+        Cp              ${Cp};
+    }
+    equationOfState
+    {
+        rho             ${rho};
+    }
+}
+
+radiation
+{
+    radiationModel  opaqueSolid;
+    absorptionEmissionModel none;
+    scatterModel    none;
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/3DBaffle/extrudeModel

@@ -5,34 +5,14 @@
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      extrudeToRegionMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-region          baffleRegion;
-
-faceZones       (baffleFaces2);
-
-oneD            false;
-
-extrudeModel    linearNormal;
-
-nLayers         10;
-
-expansionRatio  1;
-
-adaptMesh       yes; // apply directMapped to both regions
-
-sampleMode      nearestPatchFace;
 
+extrudeModel        linearNormal;
+nLayers             ${nLayers};
+expansionRatio      1;
+columnCells         ${oneD}; //3D
 linearNormalCoeffs
 {
-    thickness       0.02;
+    thickness           ${thickness};
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/baffle3DSetup

@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+baffleRegionName    baffle3DRegion;
+
+masterPatchName     baffle3DRegionMaster;
+
+slavePatchName      baffle3DRegionSlave;
+
+oneD                false;
+
+nLayers             50;
+
+thickness           0.02;
+
+Kappa               0.01;
+
+Cp                  15;
+
+rho                 80;
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/include/wallBafflePatches

@@ -0,0 +1,50 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+alphat
+{
+    type            compressible::alphatWallFunction;
+    value           uniform 0;
+}
+
+epsilon
+{
+    type            compressible::epsilonWallFunction;
+    value           uniform 0.01;
+}
+
+k
+{
+    type            compressible::kqRWallFunction;
+    value           uniform 0.01;
+}
+
+mut
+{
+    type            mutkWallFunction;
+    value           uniform 0;
+}
+
+p
+{
+    type            calculated;
+    value           uniform 101325;
+}
+
+p_rgh
+{
+    type            fixedFluxPressure;
+}
+
+U
+{
+    type            fixedValue;
+    value           uniform (0 0 0);
+}
+
+// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/k

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 0 2 -2 0 0 0 0 ];
+dimensions      [0 2 -2 0 0 0 0];
 
 internalField   uniform 0.1;
 
@@ -44,6 +44,7 @@ boundaryField
     {
         type            empty;
     }
+
 }
 
 

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/mut

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 1 -1 -1 0 0 0 0 ];
+dimensions      [1 -1 -1 0 0 0 0];
 
 internalField   uniform 0;
 
@@ -24,11 +24,17 @@ boundaryField
     floor
     {
         type            mutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
         value           uniform 0;
     }
     ceiling
     {
         type            mutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
         value           uniform 0;
     }
     inlet
@@ -45,6 +51,7 @@ boundaryField
     {
         type            empty;
     }
+
 }
 
 

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 1 -1 -2 0 0 0 0 ];
+dimensions      [1 -1 -2 0 0 0 0];
 
 internalField   uniform 101325;
 
@@ -24,27 +24,28 @@ boundaryField
     floor
     {
         type            calculated;
-        value           $internalField;
+        value           uniform 101325;
     }
     ceiling
     {
         type            calculated;
-        value           $internalField;
+        value           uniform 101325;
     }
     inlet
     {
         type            calculated;
-        value           $internalField;
+        value           uniform 101325;
     }
     outlet
     {
         type            calculated;
-        value           $internalField;
+        value           uniform 101325;
     }
     fixedWalls
     {
         type            empty;
     }
+
 }
 
 

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/0.org/p_rgh

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [ 1 -1 -2 0 0 0 0 ];
+dimensions      [1 -1 -2 0 0 0 0];
 
 internalField   uniform 101325;
 
@@ -24,27 +24,31 @@ boundaryField
     floor
     {
         type            fixedFluxPressure;
-        value           $internalField;
+        gradient        uniform 0;
+        value           uniform 101325;
     }
     ceiling
     {
         type            fixedFluxPressure;
-        value           $internalField;
+        gradient        uniform 0;
+        value           uniform 101325;
     }
     inlet
     {
         type            fixedFluxPressure;
-        value           $internalField;
+        gradient        uniform 0;
+        value           uniform 101325;
     }
     outlet
     {
         type            fixedValue;
-        value           $internalField;
+        value           uniform 101325;
     }
     fixedWalls
     {
         type            empty;
     }
+
 }
 
 

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/Allclean

@@ -5,7 +5,8 @@ cd ${0%/*} || exit 1    # run from this directory
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 cleanCase
-rm -rf constant/baffleRegion/polyMesh
+rm -rf constant/baffle3DRegion
-rm -rf sets 0
+rm -rf constant/polyMesh/boundary
+rm -rf 0
 
 # ----------------------------------------------------------------- end-of-file

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/Allrun

@@ -4,25 +4,11 @@
 # Get application name
 application=`getApplication`
 
-cp -r 0.org 0
 runApplication blockMesh
 
-runApplication topoSet
+cp -r 0.org 0
 
-unset FOAM_SETNAN
+# Create 1D and 3D baffles
-unset FOAM_SIGFPE
-
-# Create first baffle
 runApplication createBaffles -overwrite
 
-# Create region
-runApplication extrudeToRegionMesh -overwrite
-
-# Set Bc's for the region baffle
-runApplication changeDictionary -dict system/changeDictionaryDict.baffleRegion -literalRE
-rm log.changeDictionary
-
-# Reset proper values at the region
-runApplication changeDictionary -region baffleRegion -literalRE
-
 runApplication $application

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/polyMesh/blockMeshDict

@@ -22,10 +22,10 @@ vertices
     (10 0 0)
     (10 5 0)
     (0 5 0)
-    (0 0 10)
+    (0 0 1)
-    (10 0 10)
+    (10 0 1)
-    (10 5 10)
+    (10 5 1)
-    (0 5 10)
+    (0 5 1)
 );
 
 blocks

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/triSurface/baffle1D.stl

@@ -0,0 +1,16 @@
+solid ascii
+ facet normal -1 0 0
+  outer loop
+   vertex 0.3 0 0
+   vertex 0.3 0 0.1
+   vertex 0.3 0.2 0
+  endloop
+ endfacet
+ facet normal -1 0 0
+  outer loop
+   vertex 0.3 0.2 0.1
+   vertex 0.3 0.2 0
+   vertex 0.3 0 0.1
+  endloop
+ endfacet
+endsolid

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/constant/triSurface/baffle3D.stl

@@ -0,0 +1,58 @@
+solid ascii
+ facet normal -1 0 0
+  outer loop
+   vertex 0.59 0 0
+   vertex 0.59 0 0.05
+   vertex 0.59 0.1 0
+  endloop
+ endfacet
+ facet normal -1 0 0
+  outer loop
+   vertex 0.59 0.1 0.05
+   vertex 0.59 0.1 0
+   vertex 0.59 0 0.05
+  endloop
+ endfacet
+ facet normal -1 0 0
+  outer loop
+   vertex 0.59 0 0.05
+   vertex 0.59 0 0.1
+   vertex 0.59 0.1 0.05
+  endloop
+ endfacet
+ facet normal -1 0 0
+  outer loop
+   vertex 0.59 0.1 0.1
+   vertex 0.59 0.1 0.05
+   vertex 0.59 0 0.1
+  endloop
+ endfacet
+ facet normal -1 0 0
+  outer loop
+   vertex 0.59 0.1 0
+   vertex 0.59 0.1 0.05
+   vertex 0.59 0.2 0
+  endloop
+ endfacet
+ facet normal -1 0 0
+  outer loop
+   vertex 0.59 0.2 0.05
+   vertex 0.59 0.2 0
+   vertex 0.59 0.1 0.05
+  endloop
+ endfacet
+ facet normal -1 0 0
+  outer loop
+   vertex 0.59 0.1 0.05
+   vertex 0.59 0.1 0.1
+   vertex 0.59 0.2 0.05
+  endloop
+ endfacet
+ facet normal -1 0 0
+  outer loop
+   vertex 0.59 0.2 0.1
+   vertex 0.59 0.2 0.05
+   vertex 0.59 0.1 0.1
+  endloop
+ endfacet
+endsolid

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/baffleRegion/changeDictionaryDict

@@ -1,58 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      changeDictionaryDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dictionaryReplacement
-{
-    T
-    {
-        boundaryField
-        {
-            "region0_to.*"
-            {
-                type                compressible::thermalBaffle;
-                Tnbr                T;
-                kappa               solidThermo;
-                kappaName           none;
-                QrNbr               none;
-                Qr                  none;
-                value               uniform 300;
-            }
-            baffleFaces2_side
-            {
-                type            zeroGradient;
-            }
-            floor
-            {
-                type            fixedValue;
-                value           uniform 300;
-            }
-            fixedWalls
-            {
-                type            empty;
-            }
-        }
-    }
-
-    boundary
-    {
-        floor
-        {
-            type        patch;
-        }
-    }
-}
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict

@@ -1,125 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      changeDictionaryDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dictionaryReplacement
-{
-    alphat
-    {
-        boundaryField
-        {
-            "baffle.*"
-            {
-                type            compressible::alphatWallFunction;
-                value           uniform 0;
-            }
-        }
-    }
-    epsilon
-    {
-        boundaryField
-        {
-            "baffle.*"
-            {
-                type            compressible::epsilonWallFunction;
-                value           uniform 0.01;
-            }
-        }
-    }
-    k
-    {
-        boundaryField
-        {
-            "baffle.*"
-            {
-                type            compressible::kqRWallFunction;
-                value           uniform 0.01;
-            }
-        }
-    }
-    mut
-    {
-        boundaryField
-        {
-            "baffle.*"
-            {
-                type            mutkWallFunction;
-                value           uniform 0.0;
-            }
-        }
-    }
-    p
-    {
-        boundaryField
-        {
-             "baffle.*"
-            {
-                type            calculated;
-                value           uniform 101325;
-            }
-        }
-    }
-    p_rgh
-    {
-        boundaryField
-        {
-            "baffle.*"
-            {
-                type            fixedFluxPressure;
-                value           uniform 0;
-            }
-        }
-    }
-    T
-    {
-        boundaryField
-        {
-            "baffle.*"
-            {
-                type        compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
-                baffleActivated yes;
-                thickness       uniform 0.005;  // thickness [m]
-                Qs              uniform 100;    // heat flux [W/m2]
-                transport
-                {
-                    kappa           1.0;
-                }
-                thermodynamics
-                {
-                    Hf              0;
-                    Cp              0;
-                }
-                equationOfState
-                {
-                    rho             0;
-                }
-                value           uniform 300;
-            }
-        }
-    }
-    U
-    {
-        boundaryField
-        {
-            "baffle.*"
-            {
-                type            fixedValue;
-                value           uniform (0 0 0);
-            }
-        }
-    }
-}
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffle

@@ -1,130 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      changeDictionaryDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dictionaryReplacement
-{
-    alphat
-    {
-        boundaryField
-        {
-            "baffle1.*"
-            {
-                type            compressible::alphatWallFunction;
-                value           uniform 0;
-            }
-        }
-    }
-    epsilon
-    {
-        boundaryField
-        {
-            "baffle1.*"
-            {
-                type            compressible::epsilonWallFunction;
-                value           uniform 0.01;
-            }
-        }
-    }
-    k
-    {
-        boundaryField
-        {
-            "baffle1.*"
-            {
-                type            compressible::kqRWallFunction;
-                value           uniform 0.01;
-            }
-        }
-    }
-    mut
-    {
-        boundaryField
-        {
-            "baffle1.*"
-            {
-                type            mutkWallFunction;
-                value           uniform 0.0;
-            }
-        }
-    }
-    p
-    {
-        boundaryField
-        {
-            "baffle1.*"
-            {
-                type            calculated;
-                value           $internalField;
-            }
-        }
-    }
-    p_rgh
-    {
-        boundaryField
-        {
-            "baffle1.*"
-            {
-                type            fixedFluxPressure;
-                value           $internalField;
-            }
-        }
-    }
-    T
-    {
-        boundaryField
-        {
-            "baffle1Wall.*"
-            {
-                type   compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
-                baffleActivated yes;
-                thickness       uniform 0.005;  // thickness [m]
-                Qs              uniform 100;    // heat flux [W/m2]
-                specie
-                {
-                    nMoles          1;
-                    molWeight       20;
-                }
-                transport
-                {
-                    kappa           1;
-                }
-                thermodynamics
-                {
-                    Hf              0;
-                    Cp              10;
-                }
-                equationOfState
-                {
-                    rho             10;
-                }
-                value           uniform 300;
-            }
-        }
-    }
-    U
-    {
-        boundaryField
-        {
-            "baffle1.*"
-            {
-                type            fixedValue;
-                value           uniform (0 0 0);
-            }
-        }
-    }
-}
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/changeDictionaryDict.baffleRegion

@@ -1,171 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      changeDictionaryDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dictionaryReplacement
-{
-    alphat
-    {
-        boundaryField
-        {
-            "region0_to_.*"
-            {
-                type            compressible::alphatWallFunction;
-                value           uniform 0;
-            }
-        }
-    }
-    epsilon
-    {
-        boundaryField
-        {
-            "region0_to_.*"
-            {
-                type            compressible::epsilonWallFunction;
-                value           uniform 0.01;
-            }
-        }
-    }
-    k
-    {
-        boundaryField
-        {
-            "region0_to_.*"
-            {
-                type            compressible::kqRWallFunction;
-                value           uniform 0.01;
-            }
-
-        }
-    }
-    mut
-    {
-        boundaryField
-        {
-            "region0_to_.*"
-            {
-                type            mutkWallFunction;
-                value           uniform 0.0;
-            }
-        }
-    }
-    p
-    {
-        boundaryField
-        {
-            "region0_to_.*"
-            {
-                type            calculated;
-                value           $internalField;
-            }
-        }
-    }
-    p_rgh
-    {
-        boundaryField
-        {
-            "region0_to_.*"
-            {
-                type            fixedFluxPressure;
-                value           $internalField;
-            }
-
-        }
-    }
-    T
-    {
-        boundaryField
-        {
-
-            "region0_to.*"
-            {
-                type                compressible::thermalBaffle;
-
-                // Coupled BC.
-                Tnbr                T;
-                kappa               fluidThermo;
-                kappaName           none;
-                QrNbr               none;
-                Qr                  none;
-
-                // Thermo baffle model
-                thermalBaffleModel   thermalBaffle;
-                regionName          baffleRegion;
-                infoOutput          no;
-                active              yes;
-                thermalBaffleCoeffs
-                {
-                }
-
-                // Solid thermo
-                thermoType
-                {
-                    type            heSolidThermo;
-                    mixture         pureMixture;
-                    transport       constIso;
-                    thermo          hConst;
-                    equationOfState rhoConst;
-                    specie          specie;
-                    energy          sensibleEnthalpy;
-                }
-
-                mixture
-                {
-                    specie
-                    {
-                        nMoles          1;
-                        molWeight       20;
-                    }
-                    transport
-                    {
-                        kappa           0.01;
-                    }
-                    thermodynamics
-                    {
-                        Hf              0;
-                        Cp              15;
-                    }
-                    equationOfState
-                    {
-                        rho             80;
-                    }
-                }
-
-                radiation
-                {
-                    radiationModel  opaqueSolid;
-                    absorptionEmissionModel none;
-                    scatterModel    none;
-                }
-
-                value               uniform 300;
-            }
-        }
-    }
-    U
-    {
-        boundaryField
-        {
-            "region0_to_.*"
-            {
-                type            fixedValue;
-                value           uniform (0 0 0);
-            }
-
-        }
-    }
-}
-
-// ************************************************************************* //

tutorials/heatTransfer/buoyantSimpleFoam/circuitBoardCooling/system/createBafflesDict

@@ -19,109 +19,100 @@ FoamFile
 // faces.
 internalFacesOnly true;
 
-
 // Baffles to create.
 baffles
 {
     baffleFacesThermoBaffle1D
     {
         //- Use predefined faceZone to select faces and orientation.
-        type        faceZone;
+        type        searchableSurface;
-        zoneName    baffleFaces;
+        surface     triSurfaceMesh;
-
+        name        baffle1D.stl;
 
         patches
         {
             master
             {
                 //- Master side patch
-                name            baffle1Wall_0;
+                name            baffle1DWall_master;
-
-                type            mappedWall;
-                sampleMode      nearestPatchFace;
-                samplePatch     baffle1Wall_1;
-                offset          (0 0 0);
 
+                type            wall;
+                inGroups        (baffleWallGroup);
 
                 patchFields
                 {
-                    T
+                    #include "./0/include/wallBafflePatches"
-                    {
+                    #include "./0/include/1DBaffle/1DTemperatureMasterBafflePatches"
-                        type   compressible::thermalBaffle1D<hConstSolidThermoPhysics>;
-                        baffleActivated yes;
-                        thickness       uniform 0.005;  // thickness [m]
-                        Qs              uniform 100;    // heat flux [W/m2]
-                        specie
-                        {
-                            nMoles          1;
-                            molWeight       20;
-                        }
-                        transport
-                        {
-                            kappa           1;
-                        }
-                        thermodynamics
-                        {
-                            Hf              0;
-                            Cp              10;
-                        }
-                        equationOfState
-                        {
-                            rho             10;
-                        }
-
-                        value               uniform 300;
-                    }
-
-                    alphat
-                    {
-                        type            compressible::alphatWallFunction;
-                        value           uniform 0;
-                    }
-                    epsilon
-                    {
-                        type            compressible::epsilonWallFunction;
-                        value           uniform 0.01;
-                    }
-                    k
-                    {
-                        type            compressible::kqRWallFunction;
-                        value           uniform 0.01;
-                    }
-                    mut
-                    {
-                        type            mutkWallFunction;
-                        value           uniform 0;
-                    }
-                    p
-                    {
-                        type            calculated;
-                        value           uniform 101325;
-                    }
-                    p_rgh
-                    {
-                        type            fixedFluxPressure;
-                    }
-                    U
-                    {
-                        type            fixedValue;
-                        value           uniform (0 0 0);
-                    }
                 }
             }
             slave
             {
                 //- Slave side patch
-                name            baffle1Wall_1;
+                name            baffle1DWall_slave;
+
+                type            wall;
+                inGroups        (baffleWallGroup);
+
+                patchFields
+                {
+                    #include "./0/include/wallBafflePatches"
+                    #include "./0/include/1DBaffle/1DTemperatureSlaveBafflePatches"
+                }
+            }
+        }
+    }
+
+    #include "./0/include/baffle3DSetup"
+
+    baffleFacesThermoBaffle3D
+    {
+        type        searchableSurface;
+        surface     triSurfaceMesh;
+        name        baffle3D.stl;
+
+        patches
+        {
+            master
+            {
+                //- Master side patch
+                name            ${masterPatchName};
 
                 type            mappedWall;
+
+
+                type            interRegionMappedWallGenerator;
+                                inGroups        (baffleWallGroup);
+
                 sampleMode      nearestPatchFace;
-                samplePatch     baffle1Wall_0;
+                sampleRegion    ${baffleRegionName};
+                samplePatch     bottom;
+                offsetMode      uniform;
                 offset          (0 0 0);
 
                 patchFields
                 {
-                    ${...master.patchFields}
+                    #include "./0/include/wallBafflePatches"
+                    #include "./0/include/3DBaffle/3DTemperatureMasterBafflePatches"
+                }
+            }
+            slave
+            {
+                //- Slave side patch
+                name            ${slavePatchName};
+
+                type            mappedWall;
+                inGroups        (baffleWallGroup);
+
+                sampleMode      nearestPatchFace;
+                sampleRegion    ${baffleRegionName};
+                samplePatch     top;
+                offsetMode      uniform;
+                offset          (0 0 0);
+
+                patchFields
+                {
+                    #include "./0/include/wallBafflePatches"
+                    #include "./0/include/3DBaffle/3DTemperatureSlaveBafflePatches"
                 }
             }
         }