Rosenbrock21: A-stable 2nd-order Rosenbrock ODE solver

src/ODE/Make/files

@@ -8,6 +8,7 @@ ODESolvers/Trapezoid/Trapezoid.C
 ODESolvers/RKF45/RKF45.C
 ODESolvers/RKCK45/RKCK45.C
 ODESolvers/RKDP45/RKDP45.C
+ODESolvers/Rosenbrock21/Rosenbrock21.C
 ODESolvers/Rosenbrock43/Rosenbrock43.C
 ODESolvers/rodas43/rodas43.C
 ODESolvers/SIBS/SIBS.C

src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C

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

src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::Rosenbrock21
+
+Description
+    Embedded Rosenbrock ODE solver of order (1)2.
+
+SourceFiles
+    Rosenbrock21.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Rosenbrock21_H
+#define Rosenbrock21_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Rosenbrock21 Declaration
+\*---------------------------------------------------------------------------*/
+
+class Rosenbrock21
+:
+    public ODESolver,
+    public adaptiveSolver
+{
+    // Private data
+
+        mutable scalarField k1_;
+        mutable scalarField k2_;
+        mutable scalarField err_;
+        mutable scalarField dydx_;
+        mutable scalarField dfdx_;
+        mutable scalarSquareMatrix dfdy_;
+        mutable scalarSquareMatrix a_;
+        mutable labelList pivotIndices_;
+
+        static const scalar
+            a21,
+            c21,
+            b1, b2,
+            gamma,
+            c2,
+            e1, e2,
+            d1, d2;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Rosenbrock21");
+
+
+    // Constructors
+
+        //- Construct from ODE
+        Rosenbrock21(const ODESystem& ode, const dictionary& dict);
+
+
+    // Member Functions
+
+        //- Solve a single step dx and return the error
+        scalar solve
+        (
+            const ODESystem& ode,
+            const scalar x0,
+            const scalarField& y0,
+            const scalarField& dydx0,
+            const scalar dx,
+            scalarField& y
+        ) const;
+
+        //- Solve the ODE system and the update the state
+        void solve
+        (
+            const ODESystem& ode,
+            scalar& x,
+            scalarField& y,
+            scalar& dxTry
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //