From d179d9b6b99b22dc0be5aac2f3b1239185460e46 Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 6 Nov 2013 23:19:03 +0000 Subject: [PATCH 1/5] 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 a22ec1d3be..c7e68837e0 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 0000000000..e471e3ffe1 --- /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 . + +\*---------------------------------------------------------------------------*/ + +#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. + +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 0000000000..6221788ea6 --- /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 . + +\*---------------------------------------------------------------------------*/ + +#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. + +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 + +// ************************************************************************* // From 5023472727ff68597033eb3a34970a3ad2f5fb5e Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 6 Nov 2013 23:20:08 +0000 Subject: [PATCH 2/5] Cleanup --- src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C | 4 ++-- src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C | 2 +- src/ODE/ODESolvers/rodas43/rodas43.C | 3 --- src/ODE/ODESolvers/rodas43/rodas43.H | 3 +-- 4 files changed, 4 insertions(+), 8 deletions(-) diff --git a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C index de6b968348..13371a8e90 100644 --- a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C +++ b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C @@ -42,8 +42,8 @@ const scalar Rosenbrock21::b2 = (1.0/2.0)/gamma, Rosenbrock21::e1 = b1 - 1.0/gamma, Rosenbrock21::e2 = b2, - Rosenbrock21::d1 = 1, - Rosenbrock21::d2 = -1; + Rosenbrock21::d1 = gamma, + Rosenbrock21::d2 = -gamma; } diff --git a/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C index 04b55dbbe6..061390fdbc 100644 --- a/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C +++ b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C @@ -91,7 +91,7 @@ const scalar Rosenbrock43::e3 = 0, Rosenbrock43::e4 = 125.0/108.0, - Rosenbrock43::gamma = 0.5, + Rosenbrock43::gamma = 1.0/2.0, Rosenbrock43::c2 = 1, Rosenbrock43::c3 = 3.0/5.0, diff --git a/src/ODE/ODESolvers/rodas43/rodas43.C b/src/ODE/ODESolvers/rodas43/rodas43.C index 74a3da1d80..92b78a7ab4 100644 --- a/src/ODE/ODESolvers/rodas43/rodas43.C +++ b/src/ODE/ODESolvers/rodas43/rodas43.C @@ -37,9 +37,6 @@ const scalar rodas43::c2 = 0.386, rodas43::c3 = 0.21, rodas43::c4 = 0.63, - rodas43::bet2p = 0.0317, - rodas43::bet3p = 0.0635, - rodas43::bet4p = 0.3438, rodas43::d1 = 0.25, rodas43::d2 = -0.1043, rodas43::d3 = 0.1035, diff --git a/src/ODE/ODESolvers/rodas43/rodas43.H b/src/ODE/ODESolvers/rodas43/rodas43.H index a9b8ed7fc2..e77c0b423b 100644 --- a/src/ODE/ODESolvers/rodas43/rodas43.H +++ b/src/ODE/ODESolvers/rodas43/rodas43.H @@ -25,7 +25,7 @@ Class Foam::rodas43 Description - Stiffly-stable embedded Rosenbrock ODE solver of order (3)4. + L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (3)4. References: \verbatim @@ -79,7 +79,6 @@ class rodas43 static const scalar c2, c3, c4, - bet2p, bet3p, bet4p, d1, d2, d3, d4, a21, a31, a32, a41, a42, a43, From e6108fa20595e37498ab9baf45ac54888e2a3d29 Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 6 Nov 2013 23:20:23 +0000 Subject: [PATCH 3/5] forces: switch-off logging by default --- src/postProcessing/functionObjects/forces/forces/forces.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C index c24fa72729..74b892aa1e 100644 --- a/src/postProcessing/functionObjects/forces/forces/forces.C +++ b/src/postProcessing/functionObjects/forces/forces/forces.C @@ -616,7 +616,7 @@ void Foam::forces::read(const dictionary& dict) { initialised_ = false; - log_ = dict.lookupOrDefault("log", true); + log_ = dict.lookupOrDefault("log", false); if (log_) { From 69ffc9051f08262360252d6298e6bd4329d5db28 Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 7 Nov 2013 09:38:20 +0000 Subject: [PATCH 4/5] chemFoam tutorials: Changed ODE solver to rodas32 which is fastest and robust --- .../chemFoam/gri/constant/chemistryProperties | 6 +++--- .../chemFoam/h2/constant/chemistryProperties | 12 +++++++++--- .../chemFoam/ic8h18/constant/chemistryProperties | 12 +++++++++--- .../chemFoam/nc7h16/constant/chemistryProperties | 12 +++++++++--- 4 files changed, 30 insertions(+), 12 deletions(-) diff --git a/tutorials/combustion/chemFoam/gri/constant/chemistryProperties b/tutorials/combustion/chemFoam/gri/constant/chemistryProperties index e65509f396..314f83cc2d 100644 --- a/tutorials/combustion/chemFoam/gri/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/gri/constant/chemistryProperties @@ -23,7 +23,7 @@ chemistryType chemistry on; -initialChemicalTimeStep 1e-6; +initialChemicalTimeStep 1e-7; EulerImplicitCoeffs { @@ -33,9 +33,9 @@ EulerImplicitCoeffs odeCoeffs { - solver Rosenbrock43; + solver rodas32; absTol 1e-12; - relTol 1e-2; + relTol 1e-1; } diff --git a/tutorials/combustion/chemFoam/h2/constant/chemistryProperties b/tutorials/combustion/chemFoam/h2/constant/chemistryProperties index 2c0182864e..314f83cc2d 100644 --- a/tutorials/combustion/chemFoam/h2/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/h2/constant/chemistryProperties @@ -23,13 +23,19 @@ chemistryType chemistry on; -initialChemicalTimeStep 1e-10; +initialChemicalTimeStep 1e-7; + +EulerImplicitCoeffs +{ + cTauChem 1; + equilibriumRateLimiter off; +} odeCoeffs { - solver SIBS; + solver rodas32; absTol 1e-12; - relTol 1e-2; + relTol 1e-1; } diff --git a/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties b/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties index 7528b62231..314f83cc2d 100644 --- a/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties @@ -23,13 +23,19 @@ chemistryType chemistry on; -initialChemicalTimeStep 1e-10; +initialChemicalTimeStep 1e-7; + +EulerImplicitCoeffs +{ + cTauChem 1; + equilibriumRateLimiter off; +} odeCoeffs { - solver SIBS; + solver rodas32; absTol 1e-12; - relTol 1e-3; + relTol 1e-1; } diff --git a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties index d7e58ddc14..aee69a3239 100644 --- a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties @@ -23,13 +23,19 @@ chemistryType chemistry on; -initialChemicalTimeStep 1e-10; +initialChemicalTimeStep 1e-7; + +EulerImplicitCoeffs +{ + cTauChem 10; + equilibriumRateLimiter no; +} odeCoeffs { - solver SIBS; + solver rodas32; absTol 1e-14; - relTol 1e-4; + relTol 1e-1; } From 3e715dfb706e7379934f152d031b835271fb0468 Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 7 Nov 2013 12:23:17 +0000 Subject: [PATCH 5/5] floatingObject interDyMFoam/ras/floatingObject: Relax the acceleration --- .../interDyMFoam/ras/floatingObject/0.org/pointDisplacement | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement index 33d0400057..daff376a9d 100644 --- a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement +++ b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement @@ -10,7 +10,6 @@ FoamFile version 2.0; format ascii; class pointVectorField; - location "0.01"; object pointDisplacement; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -37,8 +36,9 @@ boundaryField centreOfMass (0.5 0.5 0.5); momentOfInertia (0.08622222 0.08622222 0.144); mass 9.6; - rhoInf 1; // needed only for solvers solving for kinematic pressure + rhoInf 1; report on; + accelerationRelaxation 0.3; value uniform (0 0 0); } }