Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H

@@ -54,6 +54,8 @@
         {
             phi = phiHbyA - p_rghEqn.flux();
 
+            p_rgh.relax();
+
             U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);

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

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;
 }
 
 

src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H

@@ -25,7 +25,18 @@ Class
     Foam::Rosenbrock21
 
 Description
-    Embedded Rosenbrock ODE solver of order (1)2.
+    L-stable embedded Rosenbrock ODE solver of order (1)2.
+
+    References:
+    \verbatim
+        "A second-order Rosenbrock method applied to
+         photochemical dispersion problems",
+        J. G. Verwer,
+        E. J. Spee,
+        J. G. Blom,
+        W. Hundsdorfer,
+        Siam Journal on Scientific Computing 01/1999; 20(4):1456-1480.
+    \endverbatim
 
 SourceFiles
     Rosenbrock21.C

src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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);
+}
+
+
+// ************************************************************************* //

src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+// ************************************************************************* //

src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C

@@ -34,12 +34,46 @@ namespace Foam
     addToRunTimeSelectionTable(ODESolver, Rosenbrock43, dictionary);
 
 const scalar
-
-    Rosenbrock43::a21 = 2.0,
+    // L-Stable constants from Hairer et. al.
+    Rosenbrock43::a21 = 2,
+    Rosenbrock43::a31 = 1.867943637803922,
+    Rosenbrock43::a32 = 0.2344449711399156,
+
+    Rosenbrock43::c21 = -7.137615036412310,
+    Rosenbrock43::c31 = 2.580708087951457,
+    Rosenbrock43::c32 = 0.6515950076447975,
+    Rosenbrock43::c41 = -2.137148994382534,
+    Rosenbrock43::c42 = -0.3214669691237626,
+    Rosenbrock43::c43 = -0.6949742501781779,
+
+    Rosenbrock43::b1 = 2.255570073418735,
+    Rosenbrock43::b2 = 0.2870493262186792,
+    Rosenbrock43::b3 = 0.435317943184018,
+    Rosenbrock43::b4 = 1.093502252409163,
+
+    Rosenbrock43::e1 = -0.2815431932141155,
+    Rosenbrock43::e2 = -0.0727619912493892,
+    Rosenbrock43::e3 = -0.1082196201495311,
+    Rosenbrock43::e4 = -1.093502252409163,
+
+    Rosenbrock43::gamma = 0.57282,
+    Rosenbrock43::c2 = 1.14564,
+    Rosenbrock43::c3 = 0.65521686381559,
+
+    Rosenbrock43::d1 = 0.57282,
+    Rosenbrock43::d2 = -1.769193891319233,
+    Rosenbrock43::d3 = 0.7592633437920482,
+    Rosenbrock43::d4 = -0.1049021087100450;
+
+    // Constants by Shampine
+    // More accurate than the L-Stable coefficients for small step-size
+    // but less stable for large step-size
+    /*
+    Rosenbrock43::a21 = 2,
     Rosenbrock43::a31 = 48.0/25.0,
     Rosenbrock43::a32 = 6.0/25.0,
 
-    Rosenbrock43::c21 = -8.0,
+    Rosenbrock43::c21 = -8,
     Rosenbrock43::c31 = 372.0/25.0,
     Rosenbrock43::c32 = 12.0/5.0,
 
@@ -54,17 +88,18 @@ const scalar
 
     Rosenbrock43::e1 = 34.0/108.0,
     Rosenbrock43::e2 = 7.0/36.0,
-    Rosenbrock43::e3 = 0.0,
+    Rosenbrock43::e3 = 0,
     Rosenbrock43::e4 = 125.0/108.0,
 
-    Rosenbrock43::gamma = 0.5,
-    Rosenbrock43::c2 = 1.0,
+    Rosenbrock43::gamma = 1.0/2.0,
+    Rosenbrock43::c2 = 1,
     Rosenbrock43::c3  = 3.0/5.0,
 
     Rosenbrock43::d1 = 1.0/2.0,
     Rosenbrock43::d2 = -3.0/2.0,
     Rosenbrock43::d3 = 605.0/250.0,
     Rosenbrock43::d4 = 29.0/250.0;
+    */
 }
 
 

src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.H

@@ -25,7 +25,7 @@ Class
     Foam::Rosenbrock43
 
 Description
-    Embedded Rosenbrock ODE solver of order (3)4.
+    L-stable embedded Rosenbrock ODE solver of order (3)4.
 
     Based on code from:
     \verbatim
@@ -37,12 +37,14 @@ Description
         Springer-Verlag, Berlin. 1996.
     \endverbatim
 
-    Constants from:
+    Also provided are the constants from:
     \verbatim
         "Implementation of Rosenbrock Methods"
         Shampine, L. F.,
         ACM Transactions on Mathematical Software, vol. 8, 1982, pp. 93–113.
     \endverbatim
+    with which the scheme is more accurate than with the L-Stable coefficients
+    for small step-size but less stable for large step-size.
 
 SourceFiles
     Rosenbrock43.C

src/ODE/ODESolvers/rodas32/rodas32.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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);
+}