diff --git a/applications/solvers/multiphase/interFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/alphaEqn.H
index f6e57df771d897333bcc161ab381f66f4767f4b5..f0e9dab4036dd0cde2e72f377e37290158bd008f 100644
--- a/applications/solvers/multiphase/interFoam/alphaEqn.H
+++ b/applications/solvers/multiphase/interFoam/alphaEqn.H
@@ -6,7 +6,8 @@
     phic = min(interface.cAlpha()*phic, max(phic));
     surfaceScalarField phir(phic*interface.nHatf());
 
-    if (pimple.firstIter() && MULESCorr)
+    //***HGW if (pimple.firstIter() && MULESCorr)
+    if (MULESCorr)
     {
         fvScalarMatrix alpha1Eqn
         (
diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C
index c5dec296615f95d8fb5c8a6d64f4e7c0cb46461f..e5d30268285bdd0ce8af6e89c58eaa3e69dc7b48 100644
--- a/applications/test/ODE/Test-ODE.C
+++ b/applications/test/ODE/Test-ODE.C
@@ -27,7 +27,7 @@ Description
 
 #include "argList.H"
 #include "IOmanip.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "ODESolver.H"
 #include "RK.H"
 
@@ -37,7 +37,7 @@ using namespace Foam;
 
 class testODE
 :
-    public ODE
+    public ODESystem
 {
 
 public:
@@ -107,9 +107,13 @@ int main(int argc, char *argv[])
     argList::validArgs.append("ODESolver");
     argList args(argc, argv);
 
+    // Create the ODE system
     testODE ode;
+
+    // Create the selected ODE system solver
     autoPtr<ODESolver> odeSolver = ODESolver::New(args[1], ode);
 
+    // Initialise the ODE system fields
     scalar xStart = 1.0;
     scalarField yStart(ode.nEqns());
     yStart[0] = ::Foam::j0(xStart);
@@ -117,6 +121,7 @@ int main(int argc, char *argv[])
     yStart[2] = ::Foam::jn(2, xStart);
     yStart[3] = ::Foam::jn(3, xStart);
 
+    // Print the evolution of the solution and the time-step
     scalarField dyStart(ode.nEqns());
     ode.derivatives(xStart, yStart, dyStart);
 
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C
index 9246605b51ded2be40c90280f81eb3538bded533..e9218641a70d3e16a565d20b42706ec0cb157c15 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.C
+++ b/src/ODE/ODESolvers/KRR4/KRR4.C
@@ -55,7 +55,7 @@ const scalar
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::KRR4::KRR4(const ODE& ode)
+Foam::KRR4::KRR4(const ODESystem& ode)
 :
     ODESolver(ode),
     yTemp_(n_, 0.0),
@@ -76,7 +76,7 @@ Foam::KRR4::KRR4(const ODE& ode)
 
 void Foam::KRR4::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
@@ -168,7 +168,7 @@ void Foam::KRR4::solve
             (
                 "void Foam::KRR4::solve"
                 "("
-                    "const ODE&, "
+                    "const ODESystem&, "
                     "scalar&, "
                     "scalarField&, "
                     "scalarField&, "
@@ -206,7 +206,7 @@ void Foam::KRR4::solve
     (
         "void Foam::KRR4::solve"
         "("
-            "const ODE&, "
+            "const ODESystem&, "
             "scalar&, "
             "scalarField&, "
             "scalarField&, "
diff --git a/src/ODE/ODESolvers/KRR4/KRR4.H b/src/ODE/ODESolvers/KRR4/KRR4.H
index 17e9a8caa04b1ef7ff023e0c962b4431171d8294..1eb434a0ba87501a4bebe85d7908875968288805 100644
--- a/src/ODE/ODESolvers/KRR4/KRR4.H
+++ b/src/ODE/ODESolvers/KRR4/KRR4.H
@@ -88,14 +88,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        KRR4(const ODE& ode);
+        KRR4(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C
index 7a981dfd1611a249fb35395a37760ef84fe41282..5f7188a079251aebd852bff60cfe045c468206df 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.C
+++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C
@@ -36,7 +36,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::ODESolver::ODESolver(const ODE& ode)
+Foam::ODESolver::ODESolver(const ODESystem& ode)
 :
     n_(ode.nEqns()),
     yScale_(n_),
@@ -48,7 +48,7 @@ Foam::ODESolver::ODESolver(const ODE& ode)
 
 void Foam::ODESolver::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
@@ -102,7 +102,7 @@ void Foam::ODESolver::solve
     FatalErrorIn
     (
         "ODESolver::solve"
-        "(const ODE& ode, const scalar xStart, const scalar xEnd,"
+        "(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
         "scalarField& yStart, const scalar eps, scalar& hEst) const"
     )   << "Too many integration steps"
         << exit(FatalError);
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H
index 79c3fc5d1b7698f3d50d100ab99471ce10008ac1..04315ea7dccb73abe7a142a6f9dac5f433406bcf 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.H
+++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,7 +25,7 @@ Class
     Foam::ODESolver
 
 Description
-    Selection for ODE solver
+    Abstract base-class for ODE system solvers
 
 SourceFiles
     ODESolver.C
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef ODESolver_H
 #define ODESolver_H
 
-#include "ODE.H"
+#include "ODESystem.H"
 #include "typeInfo.H"
 #include "autoPtr.H"
 
@@ -82,7 +82,7 @@ public:
             autoPtr,
             ODESolver,
             ODE,
-            (const ODE& ode),
+            (const ODESystem& ode),
             (ode)
         );
 
@@ -90,7 +90,7 @@ public:
     // Constructors
 
         //- Construct for given ODE
-        ODESolver(const ODE& ode);
+        ODESolver(const ODESystem& ode);
 
 
     // Selectors
@@ -99,7 +99,7 @@ public:
         static autoPtr<ODESolver> New
         (
             const word& ODESolverTypeName,
-            const ODE& ode
+            const ODESystem& ode
         );
 
 
@@ -112,7 +112,7 @@ public:
 
         virtual void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
@@ -126,7 +126,7 @@ public:
 
         virtual void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar xStart,
             const scalar xEnd,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C b/src/ODE/ODESolvers/ODESolver/ODESolverNew.C
index 186e5577262833b9cf1d15aed16b0c9917d10ce1..aa8ed538cc980933fbfd81c31b6c9a32042c4025 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C
+++ b/src/ODE/ODESolvers/ODESolver/ODESolverNew.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
@@ -30,7 +30,7 @@ License
 Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
 (
     const Foam::word& ODESolverTypeName,
-    const Foam::ODE& ode
+    const Foam::ODESystem& ode
 )
 {
     Info<< "Selecting ODE solver " << ODESolverTypeName << endl;
@@ -42,7 +42,8 @@ Foam::autoPtr<Foam::ODESolver> Foam::ODESolver::New
     {
         FatalErrorIn
         (
-            "ODESolver::New(const word& ODESolverTypeName, const ODE& ode)"
+            "ODESolver::New"
+            "(const word& ODESolverTypeName, const ODESystem& ode)"
         )   << "Unknown ODESolver type "
             << ODESolverTypeName << nl << nl
             << "Valid ODESolvers are : " << endl
diff --git a/src/ODE/ODESolvers/RK/RK.C b/src/ODE/ODESolvers/RK/RK.C
index 6ae5036078219974bb8eb2ebf7ba597e7c029259..06c5e6ae811a2b12175bb276900b531a2fe2443b 100644
--- a/src/ODE/ODESolvers/RK/RK.C
+++ b/src/ODE/ODESolvers/RK/RK.C
@@ -54,7 +54,7 @@ const scalar
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::RK::RK(const ODE& ode)
+Foam::RK::RK(const ODESystem& ode)
 :
     ODESolver(ode),
     yTemp_(n_, 0.0),
@@ -72,7 +72,7 @@ Foam::RK::RK(const ODE& ode)
 
 void Foam::RK::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar x,
     const scalarField& y,
     const scalarField& dydx,
@@ -142,7 +142,7 @@ void Foam::RK::solve
 
 void Foam::RK::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
diff --git a/src/ODE/ODESolvers/RK/RK.H b/src/ODE/ODESolvers/RK/RK.H
index 5ab15fa0d19a0530d211dcc27c1d38cd327a6cbd..20fa3aa789b302839b342f0cf9fe70436d985cb0 100644
--- a/src/ODE/ODESolvers/RK/RK.H
+++ b/src/ODE/ODESolvers/RK/RK.H
@@ -81,14 +81,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        RK(const ODE& ode);
+        RK(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar x,
             const scalarField& y,
             const scalarField& dydx,
@@ -100,7 +100,7 @@ public:
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C
index d86f60d5c9c4d920059ec1a73fd918e6a5f32f85..95ab44776ce87d8a49401ee25206843add960b6e 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.C
+++ b/src/ODE/ODESolvers/SIBS/SIBS.C
@@ -47,7 +47,7 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::SIBS::SIBS(const ODE& ode)
+Foam::SIBS::SIBS(const ODESystem& ode)
 :
     ODESolver(ode),
     a_(iMaxX_, 0.0),
@@ -70,7 +70,7 @@ Foam::SIBS::SIBS(const ODE& ode)
 
 void Foam::SIBS::solve
 (
-    const ODE& ode,
+    const ODESystem& ode,
     scalar& x,
     scalarField& y,
     scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H
index 7c1d9c6fbeeecb9a6304d892506a5afec00d8061..a031b1c83b1e9db531126bbbc6a0ea1b4a949f07 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.H
+++ b/src/ODE/ODESolvers/SIBS/SIBS.H
@@ -78,7 +78,7 @@ class SIBS
 
         void SIMPR
         (
-            const ODE& ode,
+            const ODESystem& ode,
             const scalar xStart,
             const scalarField& y,
             const scalarField& dydx,
@@ -110,14 +110,14 @@ public:
     // Constructors
 
         //- Construct from ODE
-        SIBS(const ODE& ode);
+        SIBS(const ODESystem& ode);
 
 
     // Member Functions
 
         void solve
         (
-            const ODE& ode,
+            const ODESystem& ode,
             scalar& x,
             scalarField& y,
             scalarField& dydx,
diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C
index 9c67ce11588e4ec01fce0d1ff5760c4b9a147840..36402a202d2604675fbfb785bf9844e3dc5b5067 100644
--- a/src/ODE/ODESolvers/SIBS/SIMPR.C
+++ b/src/ODE/ODESolvers/SIBS/SIMPR.C
@@ -29,7 +29,7 @@ License
 
 void Foam::SIBS::SIMPR
 (
-    const ODE& ode,
+    const ODESystem& ode,
     const scalar xStart,
     const scalarField& y,
     const scalarField& dydx,
diff --git a/src/ODE/ODE/ODE.H b/src/ODE/ODESystem/ODESystem.H
similarity index 80%
rename from src/ODE/ODE/ODE.H
rename to src/ODE/ODESystem/ODESystem.H
index 4c98d8dc4cc6782a50723cef95d5da9648774716..90da37c427c7d6f2a313767558e9787217682005 100644
--- a/src/ODE/ODE/ODE.H
+++ b/src/ODE/ODESystem/ODESystem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,15 +22,15 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::ODE
+    Foam::ODESystem
 
 Description
-    Abstract base class for the ODE solvers.
+    Abstract base class for the systems of ordinary differential equations.
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ODE_H
-#define ODE_H
+#ifndef ODESystem_H
+#define ODESystem_H
 
 #include "scalarField.H"
 #include "scalarMatrices.H"
@@ -41,30 +41,32 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class ODE Declaration
+                       Class ODESystem Declaration
 \*---------------------------------------------------------------------------*/
 
-class ODE
+class ODESystem
 {
 
 public:
 
     // Constructors
 
-        //- Construct null
-        ODE()
+        //- Construct null8
+        ODESystem()
         {}
 
 
     //- Destructor
-    virtual ~ODE()
+    virtual ~ODESystem()
     {}
 
 
     // Member Functions
 
+        //- Return the number of equations in the system
         virtual label nEqns() const = 0;
 
+        //- Calculate the derivatives in dydx
         virtual void derivatives
         (
             const scalar x,
@@ -72,7 +74,8 @@ public:
             scalarField& dydx
         ) const = 0;
 
-
+        //- Calculate the Jacobian of the system
+        //  Need by the stiff-system solvers
         virtual void jacobian
         (
             const scalar x,
diff --git a/src/OpenFOAM/primitives/Tuple2/Tuple2.H b/src/OpenFOAM/primitives/Tuple2/Tuple2.H
index 26ca000575bba9c812ea540b3cf6b5fd76592ca4..4b7b45a1f84cb4a2eb673b91ea12411aa4cb69cc 100644
--- a/src/OpenFOAM/primitives/Tuple2/Tuple2.H
+++ b/src/OpenFOAM/primitives/Tuple2/Tuple2.H
@@ -25,7 +25,7 @@ Class
     Foam::Tuple2
 
 Description
-    A 2-tuple.
+    A 2-tuple for storing two objects of different types.
 
 SeeAlso
     Foam::Pair for storing two objects of identical types.
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
index 59ff974315322ef3ac6bcad5e39631ee047aa9e0..3045282447067284e583c0ac88b00d27c18a7b49 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.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
@@ -53,7 +53,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_("rho"),
     lookupGravity_(-1),
     g_(vector::zero),
-    relaxationFactor_(1)
+    curTimeIndex_(-1)
 {}
 
 
@@ -71,7 +71,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(dict.lookupOrDefault<word>("rhoName", "rho")),
     lookupGravity_(-1),
     g_(vector::zero),
-    relaxationFactor_(dict.lookupOrDefault<scalar>("relaxationFactor", 1))
+    curTimeIndex_(-1)
 {
     if (rhoName_ == "rhoInf")
     {
@@ -115,7 +115,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(ptf.rhoName_),
     lookupGravity_(ptf.lookupGravity_),
     g_(ptf.g_),
-    relaxationFactor_(ptf.relaxationFactor_)
+    curTimeIndex_(-1)
 {}
 
 
@@ -133,7 +133,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
     rhoName_(ptf.rhoName_),
     lookupGravity_(ptf.lookupGravity_),
     g_(ptf.g_),
-    relaxationFactor_(ptf.relaxationFactor_)
+    curTimeIndex_(-1)
 {}
 
 
@@ -203,6 +203,13 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     const Time& t = mesh.time();
     const pointPatch& ptPatch = this->patch();
 
+    // Store the motion state at the beginning of the time-step
+    if (curTimeIndex_ != this->db().time().timeIndex())
+    {
+        motion_.newTime();
+        curTimeIndex_ = this->db().time().timeIndex();
+    }
+
     // Patch force data is valid for the current positions, so
     // calculate the forces on the motion object from this data, then
     // update the positions
@@ -231,7 +238,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
         g_ = g.value();
     }
 
-    motion_.updateForce
+    motion_.updateAcceleration
     (
         f.forceEff() + g_*motion_.mass(),
         f.momentEff(),
@@ -240,11 +247,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
 
     Field<vector>::operator=
     (
-        relaxationFactor_
-       *(
-           motion_.currentPosition(initialPoints_)
-         - initialPoints_
-       )
+        motion_.currentPosition(initialPoints_) - initialPoints_
     );
 
     fixedValuePointPatchField<vector>::updateCoeffs();
@@ -267,9 +270,6 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::write(Ostream& os) const
         os.writeKeyword("g") << g_ << token::END_STATEMENT << nl;
     }
 
-    os.writeKeyword("relaxationFactor")
-        << relaxationFactor_ << token::END_STATEMENT << nl;
-
     motion_.write(os);
 
     initialPoints_.writeEntry("initialPoints", os);
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
index db17ef3fae7444467023b5d9a7f055de18e451b4..a4d9865a3637d9cc5742cdd1b93ddf64ec30e937 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.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
@@ -84,8 +84,8 @@ class sixDoFRigidBodyDisplacementPointPatchVectorField
         //- Gravity vector to store when not available from the db
         vector g_;
 
-        //- Optional under-relaxation factor for the motion
-        scalar relaxationFactor_;
+        //- Current time index (used for updating)
+        label curTimeIndex_;
 
 
 public:
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index cb6d4cbe63e20616c8017fd7e3f59e9095bd2576..4e218fc0d316944711bffbe75cf7cdb2663fdec7 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
@@ -164,6 +164,7 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
 Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
 :
     motionState_(),
+    motionState0_(),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -173,8 +174,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
     initialQ_(I),
     momentOfInertia_(diagTensor::one*VSMALL),
     mass_(VSMALL),
-    cDamp_(0.0),
-    aLim_(VGREAT),
+    aRelax_(1.0),
     report_(false)
 {}
 
@@ -191,8 +191,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     const point& initialCentreOfMass,
     const tensor& initialQ,
     const diagTensor& momentOfInertia,
-    scalar cDamp,
-    scalar aLim,
+    scalar aRelax,
     bool report
 )
 :
@@ -205,6 +204,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
         pi,
         tau
     ),
+    motionState0_(motionState_),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -214,8 +214,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     initialQ_(initialQ),
     momentOfInertia_(momentOfInertia),
     mass_(mass),
-    cDamp_(cDamp),
-    aLim_(aLim),
+    aRelax_(aRelax),
     report_(report)
 {}
 
@@ -223,6 +222,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
 Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
 :
     motionState_(dict),
+    motionState0_(motionState_),
     restraints_(),
     restraintNames_(),
     constraints_(),
@@ -238,8 +238,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const dictionary& dict)
     ),
     momentOfInertia_(dict.lookup("momentOfInertia")),
     mass_(readScalar(dict.lookup("mass"))),
-    cDamp_(dict.lookupOrDefault<scalar>("accelerationDampingCoeff", 0.0)),
-    aLim_(dict.lookupOrDefault<scalar>("accelerationLimit", VGREAT)),
+    aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
     report_(dict.lookupOrDefault<Switch>("report", false))
 {
     addRestraints(dict);
@@ -254,6 +253,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
 )
 :
     motionState_(sDoFRBM.motionState_),
+    motionState0_(sDoFRBM.motionState0_),
     restraints_(sDoFRBM.restraints_),
     restraintNames_(sDoFRBM.restraintNames_),
     constraints_(sDoFRBM.constraints_),
@@ -263,8 +263,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
     initialQ_(sDoFRBM.initialQ_),
     momentOfInertia_(sDoFRBM.momentOfInertia_),
     mass_(sDoFRBM.mass_),
-    cDamp_(sDoFRBM.cDamp_),
-    aLim_(sDoFRBM.aLim_),
+    aRelax_(sDoFRBM.aRelax_),
     report_(sDoFRBM.report_)
 {}
 
@@ -376,91 +375,79 @@ void Foam::sixDoFRigidBodyMotion::updatePosition
 
     if (Pstream::master())
     {
-        vector aClip = a();
-        scalar aMag = mag(aClip);
-
-        if (aMag > SMALL)
-        {
-            aClip /= aMag;
-        }
-
-        if (aMag > aLim_)
-        {
-            WarningIn
-            (
-                "void Foam::sixDoFRigidBodyMotion::updatePosition"
-                "("
-                    "scalar, "
-                    "scalar"
-                ")"
-            )
-                << "Limited acceleration " << a()
-                << " to " << aClip*aLim_
-                << endl;
-        }
-
-        v() += 0.5*(1 - cDamp_)*deltaT0*aClip*min(aMag, aLim_);
-
-        pi() += 0.5*(1 - cDamp_)*deltaT0*tau();
+        v() = v0() + 0.5*deltaT0*a();
+        pi() = pi0() + 0.5*deltaT0*tau();
 
         // Leapfrog move part
-        centreOfMass() += deltaT*v();
+        centreOfMass() = centreOfMass0() + deltaT*v();
 
         // Leapfrog orientation adjustment
-        rotate(Q(), pi(), deltaT);
+        Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
+        Q() = Qpi.first();
+        pi() = Qpi.second();
     }
 
     Pstream::scatter(motionState_);
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updateForce
+void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
     const vector& fGlobal,
     const vector& tauGlobal,
     scalar deltaT
 )
 {
+    static bool first = true;
+
     // Second leapfrog velocity adjust part, required after motion and
-    // force calculation
+    // acceleration calculation
 
     if (Pstream::master())
     {
-        a() = fGlobal/mass_;
+        // Save the previous iteration accelerations for relaxation
+        vector aPrevIter = a();
+        vector tauPrevIter = tau();
 
+        // Calculate new accelerations
+        a() = fGlobal/mass_;
         tau() = (Q().T() & tauGlobal);
-
         applyRestraints();
 
+        // Relax accelerations on all but first iteration
+        if (!first)
+        {
+            a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
+            tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
+        }
+        first = false;
+
+        // Apply constraints after relaxation
         applyConstraints(deltaT);
 
-        vector aClip = a();
-        scalar aMag = mag(aClip);
+        // Correct velocities
+        v() += 0.5*deltaT*a();
+        pi() += 0.5*deltaT*tau();
 
-        if (aMag > SMALL)
+        if (report_)
         {
-            aClip /= aMag;
+            status();
         }
+    }
+
+    Pstream::scatter(motionState_);
+}
 
-        if (aMag > aLim_)
-        {
-            WarningIn
-            (
-                "void Foam::sixDoFRigidBodyMotion::updateForce"
-                "("
-                    "const vector&, "
-                    "const vector&, "
-                    "scalar"
-                ")"
-            )
-                << "Limited acceleration " << a()
-                << " to " << aClip*aLim_
-                << endl;
-        }
 
-        v() += 0.5*(1 - cDamp_)*deltaT*aClip*min(aMag, aLim_);
+void Foam::sixDoFRigidBodyMotion::updateVelocity(scalar deltaT)
+{
+    // Second leapfrog velocity adjust part, required after motion and
+    // acceleration calculation
 
-        pi() += 0.5*(1 - cDamp_)*deltaT*tau();
+    if (Pstream::master())
+    {
+        v() += 0.5*deltaT*a();
+        pi() += 0.5*deltaT*tau();
 
         if (report_)
         {
@@ -472,7 +459,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
 }
 
 
-void Foam::sixDoFRigidBodyMotion::updateForce
+void Foam::sixDoFRigidBodyMotion::updateAcceleration
 (
     const pointField& positions,
     const vectorField& forces,
@@ -493,7 +480,7 @@ void Foam::sixDoFRigidBodyMotion::updateForce
         }
     }
 
-    updateForce(fGlobal, tauGlobal, deltaT);
+    updateAcceleration(fGlobal, tauGlobal, deltaT);
 }
 
 
@@ -506,19 +493,14 @@ Foam::point Foam::sixDoFRigidBodyMotion::predictedPosition
 ) const
 {
     vector vTemp = v() + deltaT*(a() + deltaForce/mass_);
-
     vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
-
-    point centreOfMassTemp = centreOfMass() + deltaT*vTemp;
-
-    tensor QTemp = Q();
-
-    rotate(QTemp, piTemp, deltaT);
+    point centreOfMassTemp = centreOfMass0() + deltaT*vTemp;
+    Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
 
     return
     (
         centreOfMassTemp
-      + (QTemp & initialQ_.T() & (pInitial - initialCentreOfMass_))
+      + (QpiTemp.first() & initialQ_.T() & (pInitial - initialCentreOfMass_))
     );
 }
 
@@ -531,12 +513,9 @@ Foam::vector Foam::sixDoFRigidBodyMotion::predictedOrientation
 ) const
 {
     vector piTemp = pi() + deltaT*(tau() + (Q().T() & deltaMoment));
+    Tuple2<tensor, vector> QpiTemp = rotate(Q0(), piTemp, deltaT);
 
-    tensor QTemp = Q();
-
-    rotate(QTemp, piTemp, deltaT);
-
-    return (QTemp & initialQ_.T() & vInitial);
+    return (QpiTemp.first() & initialQ_.T() & vInitial);
 }
 
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
index 4ead010504bfe5a22ec2bfdee2a4b1b6eb1481a3..dcdc9f441bbd26f8d07ff9515cfb78aa10c5af2f 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.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
@@ -59,7 +59,7 @@ SourceFiles
 #include "pointField.H"
 #include "sixDoFRigidBodyMotionRestraint.H"
 #include "sixDoFRigidBodyMotionConstraint.H"
-
+#include "Tuple2.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -87,6 +87,9 @@ class sixDoFRigidBodyMotion
         //- Motion state data object
         sixDoFRigidBodyMotionState motionState_;
 
+        //- Motion state data object for previous time-step
+        sixDoFRigidBodyMotionState motionState0_;
+
         //- Restraints on the motion
         PtrList<sixDoFRigidBodyMotionRestraint> restraints_;
 
@@ -116,14 +119,8 @@ class sixDoFRigidBodyMotion
         //- Mass of the body
         scalar mass_;
 
-        //- Acceleration damping coefficient.  Modify applied acceleration:
-        //  v1 = v0 + a*dt - cDamp*a*dt
-        //     = v0 + dt*f*(1 - cDamp)/m
-        //  Increases effective mass by 1/(1 - cDamp).
-        scalar cDamp_;
-
-        //- Acceleration magnitude limit - clips large accelerations
-        scalar aLim_;
+        //- Acceleration relaxation coefficient
+        scalar aRelax_;
 
         //- Switch to turn reporting of motion data on and off
         Switch report_;
@@ -143,8 +140,14 @@ class sixDoFRigidBodyMotion
         //  frame z-axis by the given angle
         inline tensor rotationTensorZ(scalar deltaT) const;
 
-        //- Apply rotation tensors to Q for the given torque (pi) and deltaT
-        inline void rotate(tensor& Q, vector& pi, scalar deltaT) const;
+        //- Apply rotation tensors to Q0 for the given torque (pi) and deltaT
+        //  and return the rotated Q and pi as a tuple
+        inline Tuple2<tensor, vector> rotate
+        (
+            const tensor& Q0,
+            const vector& pi,
+            const scalar deltaT
+        ) const;
 
         //- Apply the restraints to the object
         void applyRestraints();
@@ -152,6 +155,7 @@ class sixDoFRigidBodyMotion
         //- Apply the constraints to the object
         void applyConstraints(scalar deltaT);
 
+
         // Access functions retained as private because of the risk of
         // confusion over what is a body local frame vector and what is global
 
@@ -199,6 +203,21 @@ class sixDoFRigidBodyMotion
             //- Return access to torque
             inline const vector& tau() const;
 
+            //- Return access to the orientation at previous time-step
+            inline const tensor& Q0() const;
+
+            //- Return access to velocity at previous time-step
+            inline const vector& v0() const;
+
+            //- Return access to acceleration at previous time-step
+            inline const vector& a0() const;
+
+            //- Return access to angular momentum at previous time-step
+            inline const vector& pi0() const;
+
+            //- Return access to torque at previous time-step
+            inline const vector& tau0() const;
+
 
         // Edit
 
@@ -244,8 +263,7 @@ public:
             const point& initialCentreOfMass,
             const tensor& initialQ,
             const diagTensor& momentOfInertia,
-            scalar cDamp = 0.0,
-            scalar aLim = VGREAT,
+            scalar aRelax = 1.0,
             bool report = false
         );
 
@@ -279,18 +297,22 @@ public:
             scalar deltaT0
         );
 
-        //- Second leapfrog velocity adjust part, required after motion and
-        // force calculation
-        void updateForce
+        //- Second leapfrog velocity adjust part
+        //  required after motion and force calculation
+        void updateAcceleration
         (
             const vector& fGlobal,
             const vector& tauGlobal,
             scalar deltaT
         );
 
+        //- Second leapfrog velocity adjust part
+        //  required after motion and force calculation
+        void updateVelocity(scalar deltaT);
+
         //- Global forces supplied at locations, calculating net force
         //  and moment
-        void updateForce
+        void updateAcceleration
         (
             const pointField& positions,
             const vectorField& forces,
@@ -354,6 +376,9 @@ public:
             //- Return const access to the centre of mass
             inline const point& centreOfMass() const;
 
+            //- Return const access to the centre of mass at previous time-step
+            inline const point& centreOfMass0() const;
+
             //- Return access to the inertia tensor
             inline const diagTensor& momentOfInertia() const;
 
@@ -366,6 +391,9 @@ public:
 
         // Edit
 
+            //- Store the motion state at the beginning of the time-step
+            inline void newTime();
+
             //- Return non-const access to the centre of mass
             inline point& centreOfMass();
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
index 19f6da1711868eb0c20573e63eb836ba718569b4..5c1eb5e7ebe6d91c0698ab69305f3b1eeb9c076d 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.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,16 +61,19 @@ Foam::sixDoFRigidBodyMotion::rotationTensorZ(scalar phi) const
 }
 
 
-inline void Foam::sixDoFRigidBodyMotion::rotate
+inline Foam::Tuple2<Foam::tensor, Foam::vector>
+Foam::sixDoFRigidBodyMotion::rotate
 (
-    tensor& Q,
-    vector& pi,
-    scalar deltaT
+    const tensor& Q0,
+    const vector& pi0,
+    const scalar deltaT
 ) const
 {
-    tensor R;
+    Tuple2<tensor, vector> Qpi(Q0, pi0);
+    tensor& Q = Qpi.first();
+    vector& pi = Qpi.second();
 
-    R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
+    tensor R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
     pi = pi & R;
     Q = Q & R;
 
@@ -89,6 +92,8 @@ inline void Foam::sixDoFRigidBodyMotion::rotate
     R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
     pi = pi & R;
     Q = Q & R;
+
+    return Qpi;
 }
 
 
@@ -176,6 +181,36 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
 }
 
 
+inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q0() const
+{
+    return motionState0_.Q();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::v0() const
+{
+    return motionState0_.v();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::a0() const
+{
+    return motionState0_.a();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::pi0() const
+{
+    return motionState0_.pi();
+}
+
+
+inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau0() const
+{
+    return motionState0_.tau();
+}
+
+
 inline Foam::point& Foam::sixDoFRigidBodyMotion::initialCentreOfMass()
 {
     return initialCentreOfMass_;
@@ -281,6 +316,12 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass() const
 }
 
 
+inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass0() const
+{
+    return motionState0_.centreOfMass();
+}
+
+
 inline const Foam::diagTensor&
 Foam::sixDoFRigidBodyMotion::momentOfInertia() const
 {
@@ -300,6 +341,12 @@ inline bool Foam::sixDoFRigidBodyMotion::report() const
 }
 
 
+inline void Foam::sixDoFRigidBodyMotion::newTime()
+{
+    motionState0_ = motionState_;
+}
+
+
 inline Foam::point& Foam::sixDoFRigidBodyMotion::centreOfMass()
 {
     return motionState_.centreOfMass();
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
index 5cc32dbe5431a83b927e28e9a52f78ed1f82db14..e16fd6d5d51a5d32cf416f54c3b719d0cd33d246 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.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
@@ -40,10 +40,8 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
         << momentOfInertia_ << token::END_STATEMENT << nl;
     os.writeKeyword("mass")
         << mass_ << token::END_STATEMENT << nl;
-    os.writeKeyword("accelerationDampingCoeff")
-        << cDamp_ << token::END_STATEMENT << nl;
-    os.writeKeyword("accelerationLimit")
-        << aLim_ << token::END_STATEMENT << nl;
+    os.writeKeyword("accelerationRelaxation")
+        << aRelax_ << token::END_STATEMENT << nl;
     os.writeKeyword("report")
         << report_ << token::END_STATEMENT << nl;
 
diff --git a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
index af8a089907a9f5561994c3f7ffe195425f696c65..4d717b0f8b8ab214a6d6ff1d0ba8487fe17dc756 100644
--- a/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/postProcessing/functionObjects/forces/pointPatchFields/derived/uncoupledSixDoFRigidBodyDisplacement/uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField.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
@@ -160,7 +160,11 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
     }
 
     // Do not modify the accelerations, except with gravity, where available
-    motion_.updateForce(gravity*motion_.mass(), vector::zero, t.deltaTValue());
+    motion_.updateAcceleration
+    (
+        gravity*motion_.mass(),
+        vector::zero, t.deltaTValue()
+    );
 
     Field<vector>::operator=
     (
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index 6d5ad3d562f7ae175e8138479d9606530c9ec1e5..8357df6dfff01a0669be17a99168a6b8830fb526 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -36,7 +36,7 @@ Foam::chemistryModel<CompType, ThermoType>::chemistryModel
 )
 :
     CompType(mesh),
-    ODE(),
+    ODESystem(),
     Y_(this->thermo().composition().Y()),
     reactions_
     (
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
index 6da538d3d9a6c5891412956121652c5bfd1f0a94..87bba8c2bb24e09fa89ec02f5290b5d1afe8d210 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.H
@@ -39,7 +39,7 @@ SourceFiles
 #define chemistryModel_H
 
 #include "Reaction.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "volFieldsFwd.H"
 #include "simpleMatrix.H"
 #include "DimensionedField.H"
@@ -60,7 +60,7 @@ template<class CompType, class ThermoType>
 class chemistryModel
 :
     public CompType,
-    public ODE
+    public ODESystem
 {
     // Private Member Functions
 
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
index b60b5c7eb54696cf3b67974a759d957d9a30da44..802a8ab87d5b9647070a225f1219dfbbde10a0e3 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.C
@@ -37,7 +37,7 @@ solidChemistryModel
 )
 :
     CompType(mesh),
-    ODE(),
+    ODESystem(),
     Ys_(this->solidThermo().composition().Y()),
     reactions_
     (
diff --git a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
index 09b04d06a9a4c2d20d8cdf2df89b19ec9d8e79d0..23cbbde6b7182e1690fdfc8b0b98625e5916f6f8 100644
--- a/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
+++ b/src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H
@@ -40,7 +40,7 @@ SourceFiles
 #define solidChemistryModel_H
 
 #include "Reaction.H"
-#include "ODE.H"
+#include "ODESystem.H"
 #include "volFieldsFwd.H"
 #include "DimensionedField.H"
 #include "simpleMatrix.H"
@@ -61,7 +61,7 @@ template<class CompType, class SolidThermo>
 class solidChemistryModel
 :
     public CompType,
-    public ODE
+    public ODESystem
 {
     // Private Member Functions
 
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
index efed1dd047d74766c45ce412b36d29ca7ddfdd16..2a0a90efa67a634198a39d9b79ec7f9606a5be1a 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/controlDict
@@ -17,7 +17,7 @@ FoamFile
 
 application     interDyMFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
 startTime       0;
 
@@ -29,7 +29,7 @@ deltaT          0.01;
 
 writeControl    adjustableRunTime;
 
-writeInterval   0.1;
+writeInterval   0.2;
 
 purgeWrite      0;
 
@@ -47,9 +47,9 @@ runTimeModifiable yes;
 
 adjustTimeStep  yes;
 
-maxCo           0.5;
-maxAlphaCo      0.5;
-maxDeltaT       0.01;
+maxCo           5;
+maxAlphaCo      5;
+maxDeltaT       1;
 
 libs
 (
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
index 852f5e0c74425f84b7f2d62fe4dda0a2ab1fec4a..df7685da654a4a4fa7b159dc88722b7c7a521f5a 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSchemes
@@ -27,7 +27,7 @@ gradSchemes
 
 divSchemes
 {
-    div(rho*phi,U)  Gauss vanLeerV;
+    div(rhoPhi,U)  Gauss vanLeerV;
     div(phi,alpha)  Gauss vanLeer;
     div(phirb,alpha) Gauss interfaceCompression;
     div(phi,k)      Gauss upwind;
@@ -55,7 +55,7 @@ fluxRequired
     default         no;
     p_rgh;
     pcorr;
-    alpha;
+    alpha.water;
 }
 
 
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
index d3c114226bcbad9cfd73fbe56db2662662c3a4b9..c539caec767bd1663cc77089af78b605b21446b0 100644
--- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
+++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    cellDisplacement
+    "cellDisplacement.*"
     {
         solver          GAMG;
         tolerance       1e-5;
@@ -29,14 +29,24 @@ solvers
         mergeLevels     1;
     }
 
-    alpha.water
+    "alpha.water.*"
     {
         nAlphaCorr      1;
-        nAlphaSubCycles 3;
+        nAlphaSubCycles 1;
         cAlpha          1;
+
+        alphaOuterCorrectors  yes;
+
+        MULESCorr       yes;
+        nLimiterIter    5;
+
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-8;
+        relTol          0;
     }
 
-    pcorr
+    "pcorr.*"
     {
         solver          PCG;
         preconditioner
@@ -120,9 +130,11 @@ solvers
 PIMPLE
 {
     momentumPredictor no;
-    nCorrectors     2;
+    nOuterCorrectors 6;
+    nCorrectors      1;
     nNonOrthogonalCorrectors 0;
     correctPhi      yes;
+    moveMeshOuterCorrectors yes;
 }
 
 relaxationFactors
@@ -132,7 +144,7 @@ relaxationFactors
     }
     equations
     {
-        "(U|k|epsilon|omega|R|nuTilda).*" 1;
+        ".*" 1;
     }
 }