From d179d9b6b99b22dc0be5aac2f3b1239185460e46 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 6 Nov 2013 23:19:03 +0000
Subject: [PATCH] ODE solvers: Added Rosenbrock order 3(2) solvers
---
src/ODE/Make/files | 2 +
.../ODESolvers/Rosenbrock32/Rosenbrock32.C | 161 +++++++++++++++++
.../ODESolvers/Rosenbrock32/Rosenbrock32.H | 135 ++++++++++++++
src/ODE/ODESolvers/rodas32/rodas32.C | 166 ++++++++++++++++++
src/ODE/ODESolvers/rodas32/rodas32.H | 135 ++++++++++++++
5 files changed, 599 insertions(+)
create mode 100644 src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C
create mode 100644 src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H
create mode 100644 src/ODE/ODESolvers/rodas32/rodas32.C
create mode 100644 src/ODE/ODESolvers/rodas32/rodas32.H
diff --git a/src/ODE/Make/files b/src/ODE/Make/files
index a22ec1d3bed..c7e68837e02 100644
--- a/src/ODE/Make/files
+++ b/src/ODE/Make/files
@@ -9,7 +9,9 @@ ODESolvers/RKF45/RKF45.C
ODESolvers/RKCK45/RKCK45.C
ODESolvers/RKDP45/RKDP45.C
ODESolvers/Rosenbrock21/Rosenbrock21.C
+ODESolvers/Rosenbrock32/Rosenbrock32.C
ODESolvers/Rosenbrock43/Rosenbrock43.C
+ODESolvers/rodas32/rodas32.C
ODESolvers/rodas43/rodas43.C
ODESolvers/SIBS/SIBS.C
ODESolvers/SIBS/SIMPR.C
diff --git a/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C
new file mode 100644
index 00000000000..e471e3ffe18
--- /dev/null
+++ b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C
@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 "Rosenbrock32.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(Rosenbrock32, 0);
+ addToRunTimeSelectionTable(ODESolver, Rosenbrock32, dictionary);
+
+const scalar
+ Rosenbrock32::a21 = 1,
+ Rosenbrock32::a31 = 1,
+ Rosenbrock32::a32 = 0,
+
+ Rosenbrock32::c21 = -1.0156171083877702091975600115545,
+ Rosenbrock32::c31 = 4.0759956452537699824805835358067,
+ Rosenbrock32::c32 = 9.2076794298330791242156818474003,
+
+ Rosenbrock32::b1 = 1,
+ Rosenbrock32::b2 = 6.1697947043828245592553615689730,
+ Rosenbrock32::b3 = -0.4277225654321857332623837380651,
+
+ Rosenbrock32::e1 = 0.5,
+ Rosenbrock32::e2 = -2.9079558716805469821718236208017,
+ Rosenbrock32::e3 = 0.2235406989781156962736090927619,
+
+ Rosenbrock32::gamma = 0.43586652150845899941601945119356,
+ Rosenbrock32::c2 = 0.43586652150845899941601945119356,
+
+ Rosenbrock32::d1 = 0.43586652150845899941601945119356,
+ Rosenbrock32::d2 = 0.24291996454816804366592249683314,
+ Rosenbrock32::d3 = 2.1851380027664058511513169485832;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict)
+:
+ ODESolver(ode, dict),
+ adaptiveSolver(ode, dict),
+ k1_(n_),
+ k2_(n_),
+ k3_(n_),
+ err_(n_),
+ dydx_(n_),
+ dfdx_(n_),
+ dfdy_(n_, n_),
+ a_(n_, n_),
+ pivotIndices_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::scalar Foam::Rosenbrock32::solve
+(
+ const ODESystem& ode,
+ const scalar x0,
+ const scalarField& y0,
+ const scalarField& dydx0,
+ const scalar dx,
+ scalarField& y
+) const
+{
+ ode.jacobian(x0, y0, dfdx_, dfdy_);
+
+ for (register label i=0; i<n_; i++)
+ {
+ for (register label j=0; j<n_; j++)
+ {
+ a_[i][j] = -dfdy_[i][j];
+ }
+
+ a_[i][i] += 1.0/(gamma*dx);
+ }
+
+ LUDecompose(a_, pivotIndices_);
+
+ // Calculate k1:
+ forAll(k1_, i)
+ {
+ k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
+ }
+
+ LUBacksubstitute(a_, pivotIndices_, k1_);
+
+ // Calculate k2:
+ forAll(y, i)
+ {
+ y[i] = y0[i] + a21*k1_[i];
+ }
+
+ ode.derivatives(x0 + c2*dx, y, dydx_);
+
+ forAll(k2_, i)
+ {
+ k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
+ }
+
+ LUBacksubstitute(a_, pivotIndices_, k2_);
+
+ // Calculate k3:
+ forAll(k3_, i)
+ {
+ k3_[i] = dydx_[i] + dx*d3*dfdx_[i]
+ + (c31*k1_[i] + c32*k2_[i])/dx;
+ }
+
+ LUBacksubstitute(a_, pivotIndices_, k3_);
+
+ // Calculate error and update state:
+ forAll(y, i)
+ {
+ y[i] = y0[i] + b1*k1_[i] + b2*k2_[i] + b3*k3_[i];
+ err_[i] = e1*k1_[i] + e2*k2_[i] + e3*k3_[i];
+ }
+
+ return normalizeError(y0, y, err_);
+}
+
+
+void Foam::Rosenbrock32::solve
+(
+ const ODESystem& odes,
+ scalar& x,
+ scalarField& y,
+ scalar& dxTry
+) const
+{
+ adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //
diff --git a/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H
new file mode 100644
index 00000000000..8cd2e50dbd9
--- /dev/null
+++ b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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::Rosenbrock32
+
+Description
+ L-stable embedded Rosenbrock ODE solver of order (3)4.
+
+ References:
+ \verbatim
+ Sandu et al,
+ "Benchmarking stiff ODE solvers for atmospheric chemistry problems II
+ Rosenbrock solvers",
+ A. Sandu,
+ J.G. Verwer,
+ J.G. Blom,
+ E.J. Spee,
+ G.R. Carmichael,
+ F.A. Potra,
+ Atmospheric Environment, Volume 31, 1997, Issue 20, Pages 3459-3472
+ \endverbatim
+
+SourceFiles
+ Rosenbrock32.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Rosenbrock32_H
+#define Rosenbrock32_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class Rosenbrock32 Declaration
+\*---------------------------------------------------------------------------*/
+
+class Rosenbrock32
+:
+ public ODESolver,
+ public adaptiveSolver
+{
+ // Private data
+
+ mutable scalarField k1_;
+ mutable scalarField k2_;
+ mutable scalarField k3_;
+ mutable scalarField err_;
+ mutable scalarField dydx_;
+ mutable scalarField dfdx_;
+ mutable scalarSquareMatrix dfdy_;
+ mutable scalarSquareMatrix a_;
+ mutable labelList pivotIndices_;
+
+ static const scalar
+ a21, a31, a32,
+ c21, c31, c32,
+ b1, b2, b3,
+ e1, e2, e3,
+ gamma,
+ c2, c3,
+ d1, d2, d3;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("Rosenbrock32");
+
+
+ // Constructors
+
+ //- Construct from ODE
+ Rosenbrock32(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
+
+// ************************************************************************* //
diff --git a/src/ODE/ODESolvers/rodas32/rodas32.C b/src/ODE/ODESolvers/rodas32/rodas32.C
new file mode 100644
index 00000000000..6221788ea64
--- /dev/null
+++ b/src/ODE/ODESolvers/rodas32/rodas32.C
@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 "rodas32.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(rodas32, 0);
+ addToRunTimeSelectionTable(ODESolver, rodas32, dictionary);
+
+const scalar
+ rodas32::c3 = 1,
+ rodas32::d1 = 1.0/2.0,
+ rodas32::d2 = 3.0/2.0,
+ rodas32::a31 = 2,
+ rodas32::a41 = 2,
+ rodas32::c21 = 4,
+ rodas32::c31 = 1,
+ rodas32::c32 = -1,
+ rodas32::c41 = 1,
+ rodas32::c42 = -1,
+ rodas32::c43 = -8.0/3.0,
+ rodas32::gamma = 1.0/2.0;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict)
+:
+ ODESolver(ode, dict),
+ adaptiveSolver(ode, dict),
+ k1_(n_),
+ k2_(n_),
+ k3_(n_),
+ dy_(n_),
+ err_(n_),
+ dydx_(n_),
+ dfdx_(n_),
+ dfdy_(n_, n_),
+ a_(n_, n_),
+ pivotIndices_(n_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::scalar Foam::rodas32::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(k2_, i)
+ {
+ k2_[i] = dydx0[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
+ }
+
+ LUBacksubstitute(a_, pivotIndices_, k2_);
+
+ // Calculate k3:
+ forAll(y, i)
+ {
+ dy_[i] = a31*k1_[i];
+ y[i] = y0[i] + dy_[i];
+ }
+
+ ode.derivatives(x0 + dx, y, dydx_);
+
+ forAll(k3_, i)
+ {
+ k3_[i] = dydx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
+ }
+
+ LUBacksubstitute(a_, pivotIndices_, k3_);
+
+ // Calculate new state and error
+ forAll(y, i)
+ {
+ dy_[i] += k3_[i];
+ y[i] = y0[i] + dy_[i];
+ }
+
+ ode.derivatives(x0 + dx, y, dydx_);
+
+ forAll(err_, i)
+ {
+ err_[i] = dydx_[i] + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
+ }
+
+ LUBacksubstitute(a_, pivotIndices_, err_);
+
+ forAll(y, i)
+ {
+ y[i] = y0[i] + dy_[i] + err_[i];
+ }
+
+ return normalizeError(y0, y, err_);
+}
+
+
+void Foam::rodas32::solve
+(
+ const ODESystem& odes,
+ scalar& x,
+ scalarField& y,
+ scalar& dxTry
+) const
+{
+ adaptiveSolver::solve(odes, x, y, dxTry);
+}
+
+
+// ************************************************************************* //
diff --git a/src/ODE/ODESolvers/rodas32/rodas32.H b/src/ODE/ODESolvers/rodas32/rodas32.H
new file mode 100644
index 00000000000..5cd05e9a658
--- /dev/null
+++ b/src/ODE/ODESolvers/rodas32/rodas32.H
@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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::rodas32
+
+Description
+ L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (2)3.
+
+ References:
+ \verbatim
+ Sandu et al,
+ "Benchmarking stiff ODE solvers for atmospheric chemistry problems II
+ Rosenbrock solvers",
+ A. Sandu,
+ J.G. Verwer,
+ J.G. Blom,
+ E.J. Spee,
+ G.R. Carmichael,
+ F.A. Potra,
+ Atmospheric Environment, Volume 31, 1997, Issue 20, Pages 3459-3472
+ \endverbatim
+
+SourceFiles
+ rodas32.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rodas32_H
+#define rodas32_H
+
+#include "ODESolver.H"
+#include "adaptiveSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class rodas32 Declaration
+\*---------------------------------------------------------------------------*/
+
+class rodas32
+:
+ public ODESolver,
+ public adaptiveSolver
+{
+ // Private data
+
+ mutable scalarField k1_;
+ mutable scalarField k2_;
+ mutable scalarField k3_;
+ mutable scalarField dy_;
+ mutable scalarField err_;
+ mutable scalarField dydx_;
+ mutable scalarField dfdx_;
+ mutable scalarSquareMatrix dfdy_;
+ mutable scalarSquareMatrix a_;
+ mutable labelList pivotIndices_;
+
+ static const scalar
+ c3,
+ d1, d2,
+ a31,
+ a41,
+ c21, c31, c32,
+ c41, c42, c43,
+ gamma;
+
+public:
+
+ //- Runtime type information
+ TypeName("rodas32");
+
+
+ // Constructors
+
+ //- Construct from ODE
+ rodas32(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
+
+// ************************************************************************* //
--
GitLab