From a5d7374737bc8556b854e9360cc5b60a4fb24211 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Mon, 11 Jul 2016 17:27:04 +0100
Subject: [PATCH] ODESolvers: Add support for efficient ODE solver resizing

Note: this reuses the existing storage rather than costly reallocation
which requires the initial allocation to be sufficient for the largest
size the ODE system might have.  Attempt to set a size larger than the
initial size is a fatal error.
---
 src/ODE/ODESolvers/Euler/Euler.C              | 19 ++++-
 src/ODE/ODESolvers/Euler/Euler.H              | 16 +++--
 src/ODE/ODESolvers/EulerSI/EulerSI.C          | 22 ++++++
 src/ODE/ODESolvers/EulerSI/EulerSI.H          | 16 +++--
 src/ODE/ODESolvers/ODESolver/ODESolver.C      | 69 +++++++++++++++----
 src/ODE/ODESolvers/ODESolver/ODESolver.H      | 48 ++++++++-----
 src/ODE/ODESolvers/ODESolver/ODESolverI.H     | 67 ++++++++++++++++++
 src/ODE/ODESolvers/RKCK45/RKCK45.C            | 25 ++++++-
 src/ODE/ODESolvers/RKCK45/RKCK45.H            | 16 +++--
 src/ODE/ODESolvers/RKDP45/RKDP45.C            | 25 ++++++-
 src/ODE/ODESolvers/RKDP45/RKDP45.H            | 14 +++-
 src/ODE/ODESolvers/RKF45/RKF45.C              | 25 ++++++-
 src/ODE/ODESolvers/RKF45/RKF45.H              | 16 +++--
 .../ODESolvers/Rosenbrock12/Rosenbrock12.C    | 24 +++++++
 .../ODESolvers/Rosenbrock12/Rosenbrock12.H    | 16 +++--
 .../ODESolvers/Rosenbrock23/Rosenbrock23.C    | 25 +++++++
 .../ODESolvers/Rosenbrock23/Rosenbrock23.H    | 16 +++--
 .../ODESolvers/Rosenbrock34/Rosenbrock34.C    | 26 +++++++
 .../ODESolvers/Rosenbrock34/Rosenbrock34.H    | 16 +++--
 src/ODE/ODESolvers/SIBS/SIBS.C                | 20 ++++++
 src/ODE/ODESolvers/SIBS/SIBS.H                | 14 +++-
 src/ODE/ODESolvers/Trapezoid/Trapezoid.C      | 19 ++++-
 src/ODE/ODESolvers/Trapezoid/Trapezoid.H      | 16 +++--
 .../adaptiveSolver/adaptiveSolver.C           | 11 ++-
 .../adaptiveSolver/adaptiveSolver.H           |  7 +-
 src/ODE/ODESolvers/rodas23/rodas23.C          | 26 +++++++
 src/ODE/ODESolvers/rodas23/rodas23.H          | 16 +++--
 src/ODE/ODESolvers/rodas34/rodas34.C          | 28 ++++++++
 src/ODE/ODESolvers/rodas34/rodas34.H          | 16 +++--
 src/ODE/ODESolvers/seulex/seulex.C            | 24 +++++++
 src/ODE/ODESolvers/seulex/seulex.H            | 14 +++-
 src/OpenFOAM/matrices/Matrix/Matrix.H         |  3 +
 src/OpenFOAM/matrices/Matrix/MatrixI.H        |  8 +++
 .../matrices/SquareMatrix/SquareMatrix.H      |  3 +
 .../matrices/SquareMatrix/SquareMatrixI.H     |  7 ++
 .../chemistryModel/chemistrySolver/ode/ode.C  | 11 ++-
 .../chemistryModel/chemistrySolver/ode/ode.H  |  3 +-
 37 files changed, 657 insertions(+), 90 deletions(-)
 create mode 100644 src/ODE/ODESolvers/ODESolver/ODESolverI.H

diff --git a/src/ODE/ODESolvers/Euler/Euler.C b/src/ODE/ODESolvers/Euler/Euler.C
index 1496099ae89..fccecee001f 100644
--- a/src/ODE/ODESolvers/Euler/Euler.C
+++ b/src/ODE/ODESolvers/Euler/Euler.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,6 +47,23 @@ Foam::Euler::Euler(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::Euler::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(err_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::Euler::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/Euler/Euler.H b/src/ODE/ODESolvers/Euler/Euler.H
index fecfc57fa03..c1bf6b128d1 100644
--- a/src/ODE/ODESolvers/Euler/Euler.H
+++ b/src/ODE/ODESolvers/Euler/Euler.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,17 +75,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         Euler(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~Euler()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -95,7 +103,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.C b/src/ODE/ODESolvers/EulerSI/EulerSI.C
index 0648b266491..143187d5bfa 100644
--- a/src/ODE/ODESolvers/EulerSI/EulerSI.C
+++ b/src/ODE/ODESolvers/EulerSI/EulerSI.C
@@ -52,6 +52,28 @@ Foam::EulerSI::EulerSI(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::EulerSI::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(err_);
+        resizeField(dydx_);
+        resizeField(dfdx_);
+        resizeMatrix(dfdy_);
+        resizeMatrix(a_);
+        resizeField(pivotIndices_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::EulerSI::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.H b/src/ODE/ODESolvers/EulerSI/EulerSI.H
index a6703cc20c7..4ad1693fe8d 100644
--- a/src/ODE/ODESolvers/EulerSI/EulerSI.H
+++ b/src/ODE/ODESolvers/EulerSI/EulerSI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,17 +82,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         EulerSI(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~EulerSI()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -102,7 +110,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C
index 823f159097a..3dce5378e39 100644
--- a/src/ODE/ODESolvers/ODESolver/ODESolver.C
+++ b/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-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,11 +34,33 @@ namespace Foam
 }
 
 
+// * * * * * * * * * * * * * Private 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;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict)
 :
     odes_(ode),
+    maxN_(ode.nEqns()),
     n_(ode.nEqns()),
     absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
     relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-4)),
@@ -54,6 +76,7 @@ Foam::ODESolver::ODESolver
 )
 :
     odes_(ode),
+    maxN_(ode.nEqns()),
     n_(ode.nEqns()),
     absTol_(absTol),
     relTol_(relTol),
@@ -63,22 +86,42 @@ Foam::ODESolver::ODESolver
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::scalar Foam::ODESolver::normalizeError
-(
-    const scalarField& y0,
-    const scalarField& y,
-    const scalarField& err
-) const
+bool Foam::ODESolver::resize()
 {
-    // Calculate the maximum error
-    scalar maxErr = 0.0;
-    forAll(err, i)
+    if (odes_.nEqns() != n_)
     {
-        scalar tol = absTol_[i] + relTol_[i]*max(mag(y0[i]), mag(y[i]));
-        maxErr = max(maxErr, mag(err[i])/tol);
+        if (odes_.nEqns() > maxN_)
+        {
+            FatalErrorInFunction
+                << "Specified number of equations " << odes_.nEqns()
+                << " greater than maximum " << maxN_
+                << abort(FatalError);
+        }
+
+        n_ = odes_.nEqns();
+
+        resizeField(absTol_);
+        resizeField(relTol_);
+
+        return true;
+    }
+    else
+    {
+        return false;
     }
+}
 
-    return maxErr;
+
+void Foam::ODESolver::solve
+(
+    scalar& x,
+    scalarField& y,
+    scalar& dxTry
+) const
+{
+    stepState step(dxTry);
+    solve(x, y, step);
+    dxTry = step.dxTry;
 }
 
 
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H
index de4ef399f89..fe443187adb 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-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,8 +58,11 @@ protected:
         //- Reference to ODESystem
         const ODESystem& odes_;
 
-        //- Size of the ODESystem
-        label n_;
+        //- Maximum size of the ODESystem
+        const label maxN_;
+
+        //- Size of the ODESystem (adjustable)
+        mutable label n_;
 
         //- Absolute convergence tolerance per step
         scalarField absTol_;
@@ -90,6 +93,8 @@ protected:
 
 public:
 
+    friend class ODESystem;
+
     //- Runtime type information
     TypeName("ODESolver");
 
@@ -161,15 +166,25 @@ public:
 
     // Member Functions
 
-        scalarField& absTol()
-        {
-            return absTol_;
-        }
+        //- Return the number of equations to solve
+        inline label nEqns() const;
+
+        //- Return access to the absolute tolerance field
+        inline scalarField& absTol();
+
+        //- Return access to the relative tolerance field
+        inline scalarField& relTol();
+
+        //- Resize the ODE solver
+        virtual bool resize() = 0;
+
+        template<class Type>
+        static inline void resizeField(UList<Type>& f, const label n);
 
-        scalarField& relTol()
-        {
-            return relTol_;
-        }
+        template<class Type>
+        inline void resizeField(UList<Type>& f) const;
+
+        inline void resizeMatrix(scalarSquareMatrix& m) const;
 
         //- Solve the ODE system as far as possible upto dxTry
         //  adjusting the step as necessary to provide a solution within
@@ -180,12 +195,7 @@ public:
             scalar& x,
             scalarField& y,
             scalar& dxTry
-        ) const //= 0;
-        {
-            stepState step(dxTry);
-            solve(x, y, step);
-            dxTry = step.dxTry;
-        }
+        ) const;
 
         //- Solve the ODE system as far as possible upto dxTry
         //  adjusting the step as necessary to provide a solution within
@@ -216,6 +226,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "ODESolverI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/ODE/ODESolvers/ODESolver/ODESolverI.H b/src/ODE/ODESolvers/ODESolver/ODESolverI.H
new file mode 100644
index 00000000000..e89f7bffaf1
--- /dev/null
+++ b/src/ODE/ODESolvers/ODESolver/ODESolverI.H
@@ -0,0 +1,67 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+inline Foam::label Foam::ODESolver::nEqns() const
+{
+    return n_;
+}
+
+
+inline Foam::scalarField& Foam::ODESolver::absTol()
+{
+    return absTol_;
+}
+
+
+inline Foam::scalarField& Foam::ODESolver::relTol()
+{
+    return relTol_;
+}
+
+
+template<class Type>
+inline void Foam::ODESolver::resizeField(UList<Type>& f, const label n)
+{
+    f.shallowCopy(UList<Type>(f.begin(), n));
+}
+
+
+template<class Type>
+inline void Foam::ODESolver::resizeField(UList<Type>& f) const
+{
+    resizeField(f, n_);
+}
+
+
+inline void Foam::ODESolver::resizeMatrix(scalarSquareMatrix& m) const
+{
+    m.shallowResize(n_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.C b/src/ODE/ODESolvers/RKCK45/RKCK45.C
index 30ae91de1a2..2e8b2685d82 100644
--- a/src/ODE/ODESolvers/RKCK45/RKCK45.C
+++ b/src/ODE/ODESolvers/RKCK45/RKCK45.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,6 +87,29 @@ Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::RKCK45::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(yTemp_);
+        resizeField(k2_);
+        resizeField(k3_);
+        resizeField(k4_);
+        resizeField(k5_);
+        resizeField(k6_);
+        resizeField(err_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::RKCK45::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.H b/src/ODE/ODESolvers/RKCK45/RKCK45.H
index c376a5f423f..356552761ea 100644
--- a/src/ODE/ODESolvers/RKCK45/RKCK45.H
+++ b/src/ODE/ODESolvers/RKCK45/RKCK45.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,17 +98,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         RKCK45(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~RKCK45()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -118,7 +126,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.C b/src/ODE/ODESolvers/RKDP45/RKDP45.C
index f9a3a78836e..479aecc0e9c 100644
--- a/src/ODE/ODESolvers/RKDP45/RKDP45.C
+++ b/src/ODE/ODESolvers/RKDP45/RKDP45.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,6 +92,29 @@ Foam::RKDP45::RKDP45(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::RKDP45::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(yTemp_);
+        resizeField(k2_);
+        resizeField(k3_);
+        resizeField(k4_);
+        resizeField(k5_);
+        resizeField(k6_);
+        resizeField(err_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::RKDP45::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.H b/src/ODE/ODESolvers/RKDP45/RKDP45.H
index 9998c15f6b0..75aa756f9cc 100644
--- a/src/ODE/ODESolvers/RKDP45/RKDP45.H
+++ b/src/ODE/ODESolvers/RKDP45/RKDP45.H
@@ -102,17 +102,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         RKDP45(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~RKDP45()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -122,7 +130,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/RKF45/RKF45.C b/src/ODE/ODESolvers/RKF45/RKF45.C
index 92804e06342..7ca9a232eaa 100644
--- a/src/ODE/ODESolvers/RKF45/RKF45.C
+++ b/src/ODE/ODESolvers/RKF45/RKF45.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -88,6 +88,29 @@ Foam::RKF45::RKF45(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::RKF45::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(yTemp_);
+        resizeField(k2_);
+        resizeField(k3_);
+        resizeField(k4_);
+        resizeField(k5_);
+        resizeField(k6_);
+        resizeField(err_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::RKF45::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/RKF45/RKF45.H b/src/ODE/ODESolvers/RKF45/RKF45.H
index 5198d5553f2..cc967715e1f 100644
--- a/src/ODE/ODESolvers/RKF45/RKF45.H
+++ b/src/ODE/ODESolvers/RKF45/RKF45.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,17 +102,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         RKF45(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~RKF45()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -122,7 +130,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C
index fa163aec0e9..5bb3cf7513b 100644
--- a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C
+++ b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C
@@ -66,6 +66,30 @@ Foam::Rosenbrock12::Rosenbrock12(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::Rosenbrock12::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(k1_);
+        resizeField(k2_);
+        resizeField(err_);
+        resizeField(dydx_);
+        resizeField(dfdx_);
+        resizeMatrix(dfdy_);
+        resizeMatrix(a_);
+        resizeField(pivotIndices_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::Rosenbrock12::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H
index b0bf1c79eff..44d86a8505c 100644
--- a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H
+++ b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,17 +92,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         Rosenbrock12(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~Rosenbrock12()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -112,7 +120,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C
index 3b6c8cf4635..f664726e83e 100644
--- a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C
+++ b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C
@@ -79,6 +79,31 @@ Foam::Rosenbrock23::Rosenbrock23(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::Rosenbrock23::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(k1_);
+        resizeField(k2_);
+        resizeField(k3_);
+        resizeField(err_);
+        resizeField(dydx_);
+        resizeField(dfdx_);
+        resizeMatrix(dfdy_);
+        resizeMatrix(a_);
+        resizeField(pivotIndices_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::Rosenbrock23::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H
index 38541a7c99c..571a9e7ff89 100644
--- a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H
+++ b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,17 +96,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         Rosenbrock23(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~Rosenbrock23()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -116,7 +124,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C
index 5b761b267fa..477de286165 100644
--- a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C
+++ b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C
@@ -124,6 +124,32 @@ Foam::Rosenbrock34::Rosenbrock34(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::Rosenbrock34::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(k1_);
+        resizeField(k2_);
+        resizeField(k3_);
+        resizeField(k4_);
+        resizeField(err_);
+        resizeField(dydx_);
+        resizeField(dfdx_);
+        resizeMatrix(dfdy_);
+        resizeMatrix(a_);
+        resizeField(pivotIndices_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::Rosenbrock34::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H
index 2af18e5d409..f72db63e7c0 100644
--- a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H
+++ b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,17 +104,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         Rosenbrock34(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~Rosenbrock34()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -124,7 +132,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C
index a095a7c5073..2bd63ece598 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.C
+++ b/src/ODE/ODESolvers/SIBS/SIBS.C
@@ -68,6 +68,26 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+bool Foam::SIBS::resize()
+{
+    if (ODESolver::resize())
+    {
+        resizeField(yTemp_);
+        resizeField(ySeq_);
+        resizeField(yErr_);
+        resizeField(dydx0_);
+        resizeField(dfdx_);
+        resizeMatrix(dfdy_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 void Foam::SIBS::solve
 (
     scalar& x,
diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H
index 0ef0a1d40cd..9f03753019c 100644
--- a/src/ODE/ODESolvers/SIBS/SIBS.H
+++ b/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-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -114,13 +114,21 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODE system
         SIBS(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~SIBS()
+    {}
+
+
     // Member Functions
 
-        void solve
+        //- Resize the ODE solver
+        virtual bool resize();
+
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/Trapezoid/Trapezoid.C b/src/ODE/ODESolvers/Trapezoid/Trapezoid.C
index 027308434d4..395ffc92d12 100644
--- a/src/ODE/ODESolvers/Trapezoid/Trapezoid.C
+++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,6 +47,23 @@ Foam::Trapezoid::Trapezoid(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::Trapezoid::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(err_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::Trapezoid::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/Trapezoid/Trapezoid.H b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H
index 667b6ee2adc..5dce76d5a8d 100644
--- a/src/ODE/ODESolvers/Trapezoid/Trapezoid.H
+++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,17 +65,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         Trapezoid(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~Trapezoid()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -85,7 +93,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C
index f868e921d09..330684261ce 100644
--- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C
+++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,6 +46,15 @@ Foam::adaptiveSolver::adaptiveSolver
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::adaptiveSolver::resize(const label n)
+{
+    ODESolver::resizeField(dydx0_, n);
+    ODESolver::resizeField(yTemp_, n);
+
+    return true;
+}
+
+
 void Foam::adaptiveSolver::solve
 (
     const ODESystem& odes,
diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H
index 0d7dcb035c2..eef90357313 100644
--- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H
+++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,6 +74,9 @@ public:
 
     // Member Functions
 
+        //- Resize the ODE solver
+        virtual bool resize(const label n);
+
         //- Solve a single step dx and return the error
         virtual scalar solve
         (
@@ -85,7 +88,7 @@ public:
         ) const = 0;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             const ODESystem& ode,
             scalar& x,
diff --git a/src/ODE/ODESolvers/rodas23/rodas23.C b/src/ODE/ODESolvers/rodas23/rodas23.C
index 31d5e8aceec..b11b508e0e5 100644
--- a/src/ODE/ODESolvers/rodas23/rodas23.C
+++ b/src/ODE/ODESolvers/rodas23/rodas23.C
@@ -70,6 +70,32 @@ Foam::rodas23::rodas23(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::rodas23::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(k1_);
+        resizeField(k2_);
+        resizeField(k3_);
+        resizeField(dy_);
+        resizeField(err_);
+        resizeField(dydx_);
+        resizeField(dfdx_);
+        resizeMatrix(dfdy_);
+        resizeMatrix(a_);
+        resizeField(pivotIndices_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::rodas23::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/rodas23/rodas23.H b/src/ODE/ODESolvers/rodas23/rodas23.H
index 75f0262813b..42839bfa8ff 100644
--- a/src/ODE/ODESolvers/rodas23/rodas23.H
+++ b/src/ODE/ODESolvers/rodas23/rodas23.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,17 +96,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         rodas23(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~rodas23()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -116,7 +124,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/rodas34/rodas34.C b/src/ODE/ODESolvers/rodas34/rodas34.C
index 5f403238d7f..2898b51c2d8 100644
--- a/src/ODE/ODESolvers/rodas34/rodas34.C
+++ b/src/ODE/ODESolvers/rodas34/rodas34.C
@@ -93,6 +93,34 @@ Foam::rodas34::rodas34(const ODESystem& ode, const dictionary& dict)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::rodas34::resize()
+{
+    if (ODESolver::resize())
+    {
+        adaptiveSolver::resize(n_);
+
+        resizeField(k1_);
+        resizeField(k2_);
+        resizeField(k3_);
+        resizeField(k4_);
+        resizeField(k5_);
+        resizeField(dy_);
+        resizeField(err_);
+        resizeField(dydx_);
+        resizeField(dfdx_);
+        resizeMatrix(dfdy_);
+        resizeMatrix(a_);
+        resizeField(pivotIndices_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 Foam::scalar Foam::rodas34::solve
 (
     const scalar x0,
diff --git a/src/ODE/ODESolvers/rodas34/rodas34.H b/src/ODE/ODESolvers/rodas34/rodas34.H
index 35cb468f6a4..9098a1e8b0b 100644
--- a/src/ODE/ODESolvers/rodas34/rodas34.H
+++ b/src/ODE/ODESolvers/rodas34/rodas34.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,17 +97,25 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         rodas34(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~rodas34()
+    {}
+
+
     // Member Functions
 
         //- Inherit solve from ODESolver
         using ODESolver::solve;
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve a single step dx and return the error
-        scalar solve
+        virtual scalar solve
         (
             const scalar x0,
             const scalarField& y0,
@@ -117,7 +125,7 @@ public:
         ) const;
 
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/ODE/ODESolvers/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C
index b07a9f49882..df163dd22d8 100644
--- a/src/ODE/ODESolvers/seulex/seulex.C
+++ b/src/ODE/ODESolvers/seulex/seulex.C
@@ -211,6 +211,30 @@ void Foam::seulex::extrapolate
 }
 
 
+bool Foam::seulex::resize()
+{
+    if (ODESolver::resize())
+    {
+        resizeField(dfdx_);
+        resizeMatrix(dfdy_);
+        resizeMatrix(a_);
+        resizeField(pivotIndices_);
+        resizeField(y0_);
+        resizeField(ySequence_);
+        resizeField(scale_);
+        resizeField(dy_);
+        resizeField(yTemp_);
+        resizeField(dydx_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
 void Foam::seulex::solve
 (
     scalar& x,
diff --git a/src/ODE/ODESolvers/seulex/seulex.H b/src/ODE/ODESolvers/seulex/seulex.H
index c3cd5401983..8739618ee7f 100644
--- a/src/ODE/ODESolvers/seulex/seulex.H
+++ b/src/ODE/ODESolvers/seulex/seulex.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -133,14 +133,22 @@ public:
 
     // Constructors
 
-        //- Construct from ODE
+        //- Construct from ODESystem
         seulex(const ODESystem& ode, const dictionary& dict);
 
 
+    //- Destructor
+    virtual ~seulex()
+    {}
+
+
     // Member Functions
 
+        //- Resize the ODE solver
+        virtual bool resize();
+
         //- Solve the ODE system and the update the state
-        void solve
+        virtual void solve
         (
             scalar& x,
             scalarField& y,
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index 22e289a689b..884979f3d85 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -249,6 +249,9 @@ public:
             //- Resize the matrix preserving the elements
             void setSize(const label m, const label n);
 
+            //- Resize the matrix without reallocating storage (unsafe)
+            inline void shallowResize(const label m, const label n);
+
 
         //- Return the transpose of the matrix
         Form T() const;
diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H
index 73992aac77b..2ba29e355db 100644
--- a/src/OpenFOAM/matrices/Matrix/MatrixI.H
+++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H
@@ -280,6 +280,14 @@ Foam::Matrix<Form, Type>::col
 }
 
 
+template<class Form, class Type>
+void Foam::Matrix<Form, Type>::shallowResize(const label m, const label n)
+{
+    mRows_ = m;
+    nCols_ = n;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Form, class Type>
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
index 443ac00bd2b..6baecc18e77 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H
@@ -113,6 +113,9 @@ public:
             //- Resize the matrix preserving the elements
             inline void setSize(const label m);
 
+            //- Resize the matrix without reallocating storage (unsafe)
+            inline void shallowResize(const label m);
+
 
     // Member operators
 
diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
index 12ecd2530f8..b85715ff07e 100644
--- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
+++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H
@@ -161,6 +161,13 @@ inline void Foam::SquareMatrix<Type>::setSize(const label m)
 }
 
 
+template<class Type>
+inline void Foam::SquareMatrix<Type>::shallowResize(const label m)
+{
+    Matrix<SquareMatrix<Type>, Type>::shallowResize(m, m);
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Type>
diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
index 2591f14f1ae..f9acbe96061 100644
--- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
+++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -61,7 +61,14 @@ void Foam::ode<ChemistryModel>::solve
     scalar& subDeltaT
 ) const
 {
-    label nSpecie = this->nSpecie();
+    // Reset the size of the ODE system to the simplified size when mechanism
+    // reduction is active
+    if (odeSolver_->resize())
+    {
+        odeSolver_->resizeField(cTp_);
+    }
+
+    const label nSpecie = this->nSpecie();
 
     // Copy the concentration, T and P to the total solve-vector
     for (int i=0; i<nSpecie; i++)
diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H
index 1b301dca4e6..0d1509f0c55 100644
--- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H
+++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H
@@ -55,7 +55,8 @@ class ode
     // Private data
 
         dictionary coeffsDict_;
-        autoPtr<ODESolver> odeSolver_;
+
+        mutable autoPtr<ODESolver> odeSolver_;
 
         // Solver data
         mutable scalarField cTp_;
-- 
GitLab