diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
index 4fbd2666b4ffd108ba263561c6076e4884fb69fa..c30b361e72f9badfc98172483c455fc469552275 100644
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/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);
diff --git a/src/ODE/Make/files b/src/ODE/Make/files
index a22ec1d3bedce5f75bb0a02275d85c48e0b5f857..7c361626eda507cc5a56ed4e2945c5420cab5dc2 100644
--- a/src/ODE/Make/files
+++ b/src/ODE/Make/files
@@ -9,10 +9,13 @@ 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
ODESolvers/SIBS/polyExtrapolate.C
+ODESolvers/seulex/seulex.C
LIB = $(FOAM_LIBBIN)/libODE
diff --git a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C
index de6b968348ac0b5811d0d99bf44349e91070f4b3..13371a8e90a2bb909aa71107087a1d38606b14e0 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/Rosenbrock32/Rosenbrock32.C b/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C
new file mode 100644
index 0000000000000000000000000000000000000000..e471e3ffe18da6b1ed605462c5dde0976260bbb2
--- /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 0000000000000000000000000000000000000000..8cd2e50dbd9818aabcf0cb30b9455dd777c83ff5
--- /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/Rosenbrock43/Rosenbrock43.C b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C
index 04b55dbbe61db7663fd3a300028e8b8a5f234be5..061390fdbc9e343251ac7d04da11001fdcdbedc4 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/rodas32/rodas32.C b/src/ODE/ODESolvers/rodas32/rodas32.C
new file mode 100644
index 0000000000000000000000000000000000000000..6221788ea6424d0ec3d0d3d1953ebd20af598a7b
--- /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 0000000000000000000000000000000000000000..5cd05e9a658053bf9dc113a457a9e8efce7d8e02
--- /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
+
+// ************************************************************************* //
diff --git a/src/ODE/ODESolvers/rodas43/rodas43.C b/src/ODE/ODESolvers/rodas43/rodas43.C
index 74a3da1d80b26df6c9f2c700907c64560450dea2..92b78a7ab46aa2b21e24ac0802bcb5fce77750c2 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 a9b8ed7fc2c82b23e82ead011c4075f17ec022c1..e77c0b423b285550779b1f19d45309cd3ed0e840 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,
diff --git a/src/ODE/ODESolvers/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C
new file mode 100644
index 0000000000000000000000000000000000000000..69227f19638adb94d6c44481069d585acaa74f06
--- /dev/null
+++ b/src/ODE/ODESolvers/seulex/seulex.C
@@ -0,0 +1,501 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 "seulex.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(seulex, 0);
+
+ addToRunTimeSelectionTable(ODESolver, seulex, dictionary);
+
+ const scalar
+ seulex::EPS = VSMALL,
+ seulex::STEPFAC1 = 0.6,
+ seulex::STEPFAC2 = 0.93,
+ seulex::STEPFAC3 = 0.1,
+ seulex::STEPFAC4 = 4.0,
+ seulex::STEPFAC5 = 0.5,
+ seulex::KFAC1 = 0.7,
+ seulex::KFAC2 = 0.9;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
+:
+ ODESolver(ode, dict),
+ ode_(ode),
+ nSeq_(IMAXX),
+ cost_(IMAXX),
+ dfdx_(n_),
+ factrl_(IMAXX),
+ table_(KMAXX,n_),
+ fSave_((IMAXX-1)*(IMAXX+1)/2 + 2,n_),
+ dfdy_(n_),
+ calcJac_(false),
+ dense_(false),
+ a_(n_),
+ coeff_(IMAXX,IMAXX),
+ pivotIndices_(n_, 0.0),
+ hOpt_(IMAXX),
+ work_(IMAXX),
+ ySaved_(n_),
+ ySequence_(n_),
+ scale_(n_),
+ del_(n_),
+ yTemp_(n_),
+ dyTemp_(n_),
+ delTemp_(n_)
+{
+ const scalar costfunc = 1.0, costjac = 5.0, costlu = 1.0, costsolve = 1.0;
+
+ jacRedo_ = min(1.0e-4, min(relTol_));
+ theta_ = 2.0*jacRedo_;
+ nSeq_[0] = 2;
+ nSeq_[1] = 3;
+
+ for (label i = 2; i < IMAXX; i++)
+ {
+ nSeq_[i] = 2*nSeq_[i-2];
+ }
+ cost_[0] = costjac + costlu + nSeq_[0]*(costfunc + costsolve);
+
+ for (label k=0;k<KMAXX;k++)
+ {
+ cost_[k+1] = cost_[k] + (nSeq_[k+1]-1)*(costfunc + costsolve) + costlu;
+ }
+
+ hnext_ =- VGREAT;
+
+ //NOTE: the first element of relTol_ and absTol_ are used here.
+ scalar logfact =- log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
+
+ kTarg_ = min(1, min(KMAXX - 1, label(logfact)));
+
+ for (label k = 0; k < IMAXX; k++)
+ {
+ for (label l = 0; l < k; l++)
+ {
+ scalar ratio = scalar(nSeq_[k])/nSeq_[l];
+ coeff_[k][l] = 1.0/(ratio - 1.0);
+ }
+ }
+
+ factrl_[0] = 1.0;
+ for (label k = 0; k < IMAXX - 1; k++)
+ {
+ factrl_[k+1] = (k+1)*factrl_[k];
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+void Foam::seulex::solve
+(
+ const ODESystem& ode,
+ scalar& x,
+ scalarField& y,
+ scalar& dxTry
+) const
+{
+ step(dxTry, x, y);
+ dxTry = hnext_;
+}
+
+
+void Foam::seulex::step(const scalar& htry, scalar& x, scalarField& y) const
+{
+ static bool first_step = true, last_step = false;
+ static bool forward, reject = false, prev_reject = false;
+ static scalar errold;
+ label i, k;
+ scalar fac, h, hnew, err;
+ bool firstk;
+
+ work_[0] = 1.e30;
+ h = htry;
+ forward = h > 0 ? true : false;
+
+ ySaved_ = y;
+
+ if (h != hnext_ && !first_step)
+ {
+ last_step = true;
+ }
+
+ if (reject)
+ {
+ prev_reject = true;
+ last_step = false;
+ theta_=2.0*jacRedo_;
+ }
+
+ for (i = 0; i < n_; i++)
+ {
+ scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
+ }
+
+ reject = false;
+ firstk = true;
+ hnew = mag(h);
+
+ if (theta_ > jacRedo_ && !calcJac_)
+ {
+ ode_.jacobian(x, y, dfdx_, dfdy_);
+ calcJac_ = true;
+ }
+
+ while (firstk || reject)
+ {
+ h = forward ? hnew : -hnew;
+ firstk = false;
+ reject = false;
+
+ if (mag(h) <= mag(x)*EPS)
+ {
+ WarningIn("seulex::step(const scalar")
+ << "step size underflow in step :" << h << endl;
+ }
+ label ipt=-1;
+
+ for (k = 0; k <= kTarg_ + 1; k++)
+ {
+ bool success = dy(x, ySaved_, h, k, ySequence_, ipt, scale_);
+
+ if (!success)
+ {
+ reject = true;
+ hnew = mag(h)*STEPFAC5;
+ break;
+ }
+
+ if (k == 0)
+ {
+ y = ySequence_;
+ }
+ else
+ {
+ for (i = 0; i < n_; i++)
+ table_[k-1][i] = ySequence_[i];
+ }
+ if (k != 0)
+ {
+ polyextr(k, table_, y);
+ err = 0.0;
+ for (i = 0; i < n_; i++)
+ {
+ scale_[i] = absTol_[i] + relTol_[i]*mag(ySaved_[i]);
+ err += sqr((y[i] - table_[0][i])/scale_[i]);
+ }
+ err = sqrt(err/n_);
+ if (err > 1.0/EPS || (k > 1 && err >= errold))
+ {
+ reject = true;
+ hnew = mag(h)*STEPFAC5;
+ break;
+ }
+ errold = min(4.0*err, 1.0);
+ scalar expo = 1.0/(k + 1);
+ scalar facmin = pow(STEPFAC3, expo);
+ if (err == 0.0)
+ {
+ fac = 1.0/facmin;
+ }
+ else
+ {
+ fac = STEPFAC2/pow(err/STEPFAC1, expo);
+ fac = max(facmin/STEPFAC4, min(1.0/facmin, fac));
+ }
+ hOpt_[k] = mag(h*fac);
+ work_[k] = cost_[k]/hOpt_[k];
+
+ if ((first_step || last_step) && err <= 1.0)
+ {
+ break;
+ }
+ if (k == kTarg_-1 && !prev_reject && !first_step && !last_step)
+ {
+ if (err <= 1.0)
+ {
+ break;
+ }
+ else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0)
+ {
+ reject = true;
+ kTarg_ = k;
+ if (kTarg_>1 && work_[k-1] < KFAC1*work_[k])
+ {
+ kTarg_--;
+ }
+ hnew = hOpt_[kTarg_];
+ break;
+ }
+ }
+ if (k == kTarg_)
+ {
+ if (err <= 1.0)
+ {
+ break;
+ }
+ else if (err > nSeq_[k + 1]*2.0)
+ {
+ reject = true;
+ if (kTarg_>1 && work_[k-1] < KFAC1*work_[k])
+ {
+ kTarg_--;
+ }
+ hnew = hOpt_[kTarg_];
+ break;
+ }
+ }
+ if (k == kTarg_+1)
+ {
+ if (err > 1.0)
+ {
+ reject = true;
+ if (kTarg_ > 1 && work_[kTarg_-1] < KFAC1*work_[kTarg_])
+ {
+ kTarg_--;
+ }
+ hnew = hOpt_[kTarg_];
+ }
+ break;
+ }
+ }
+ }
+ if (reject)
+ {
+ prev_reject = true;
+ if (!calcJac_)
+ {
+ theta_ = 2.0*jacRedo_;
+
+ if (theta_ > jacRedo_ && !calcJac_)
+ {
+ ode_.jacobian(x, y, dfdx_, dfdy_);
+ calcJac_ = true;
+ }
+ }
+ }
+ }
+
+ calcJac_ = false;
+
+ x += h;
+ hdid_ = h;
+ first_step = false;
+
+ label kopt;
+ if (k == 1)
+ {
+ kopt = 2;
+ }
+ else if (k <= kTarg_)
+ {
+ kopt=k;
+ if (work_[k-1] < KFAC1*work_[k])
+ {
+ kopt = k - 1;
+ }
+ else if (work_[k] < KFAC2*work_[k - 1])
+ {
+ kopt = min(k + 1, KMAXX - 1);
+ }
+ }
+ else
+ {
+ kopt = k - 1;
+ if (k > 2 && work_[k-2] < KFAC1*work_[k - 1])
+ {
+ kopt = k - 2;
+ }
+ if (work_[k] < KFAC2*work_[kopt])
+ {
+ kopt = min(k, KMAXX - 1);
+ }
+ }
+
+ if (prev_reject)
+ {
+ kTarg_ = min(kopt, k);
+ hnew = min(mag(h), hOpt_[kTarg_]);
+ prev_reject = false;
+ }
+ else
+ {
+ if (kopt <= k)
+ {
+ hnew = hOpt_[kopt];
+ }
+ else
+ {
+ if (k < kTarg_ && work_[k] < KFAC2*work_[k - 1])
+ {
+ hnew = hOpt_[k]*cost_[kopt + 1]/cost_[k];
+ }
+ else
+ {
+ hnew = hOpt_[k]*cost_[kopt]/cost_[k];
+ }
+ }
+ kTarg_ = kopt;
+ }
+
+ if (forward)
+ {
+ hnext_ = hnew;
+ }
+ else
+ {
+ hnext_ =- hnew;
+ }
+}
+
+
+bool Foam::seulex::dy
+(
+ const scalar& x,
+ scalarField& y,
+ const scalar htot,
+ const label k,
+ scalarField& yend,
+ label& ipt,
+ scalarField& scale
+) const
+{
+ label nstep = nSeq_[k];
+ scalar h = htot/nstep;
+ for (label i = 0; i < n_; i++)
+ {
+ for (label j = 0; j < n_; j++) a_[i][j] = -dfdy_[i][j];
+ a_[i][i] += 1.0/h;
+ }
+ LUDecompose(a_, pivotIndices_);
+ scalar xnew = x + h;
+
+ ode_.derivatives(xnew, y, del_);
+
+ yTemp_ = y;
+
+ LUBacksubstitute(a_, pivotIndices_, del_);
+
+ if (dense_ && nstep == k + 1)
+ {
+ ipt++;
+ for (label i = 0; i < n_; i++)
+ {
+ fSave_[ipt][i] = del_[i];
+ }
+ }
+ for (label nn = 1; nn < nstep; nn++)
+ {
+ for (label i=0;i<n_;i++)
+ {
+ yTemp_[i] += del_[i];
+ }
+ xnew += h;
+
+ ode_.derivatives(xnew, yTemp_, yend);
+
+ if (nn == 1 && k<=1)
+ {
+ scalar del1=0.0;
+ for (label i = 0; i < n_; i++)
+ {
+ del1 += sqr(del_[i]/scale[i]);
+ }
+ del1 = sqrt(del1);
+
+ ode_.derivatives(x+h, yTemp_, dyTemp_);
+ for (label i=0;i<n_;i++)
+ {
+ del_[i] = dyTemp_[i] - del_[i]/h;
+ }
+
+ LUBacksubstitute(a_, pivotIndices_, del_);
+
+ scalar del2 = 0.0;
+ for (label i = 0; i <n_ ; i++)
+ {
+ del2 += sqr(del_[i]/scale[i]);
+ }
+ del2 = sqrt(del2);
+ theta_ = del2 / min(1.0, del1 + SMALL);
+
+ if (theta_ > 1.0)
+ {
+ return false;
+ }
+ }
+
+ delTemp_ = yend;
+ LUBacksubstitute(a_, pivotIndices_, delTemp_);
+ del_ = delTemp_;
+
+ if (dense_ && nn >= nstep-k-1)
+ {
+ ipt++;
+ for (label i=0;i<n_;i++)
+ {
+ fSave_[ipt][i]=del_[i];
+ }
+ }
+ }
+
+ yend = yTemp_ + del_;
+ return true;
+}
+
+
+void Foam::seulex::polyextr
+(
+ const label k,
+ scalarRectangularMatrix& table,
+ scalarField& last
+) const
+{
+ label l=last.size();
+ for (label j = k - 1; j > 0; j--)
+ {
+ for (label i=0; i<l; i++)
+ {
+ table[j-1][i] =
+ table[j][i] + coeff_[k][j]*(table[j][i] - table[j-1][i]);
+ }
+ }
+ for (label i = 0; i < l; i++)
+ {
+ last[i] = table[0][i] + coeff_[k][0]*(table[0][i] - last[i]);
+ }
+}
+
+// ************************************************************************* //
diff --git a/src/ODE/ODESolvers/seulex/seulex.H b/src/ODE/ODESolvers/seulex/seulex.H
new file mode 100644
index 0000000000000000000000000000000000000000..71be1b3b3f144d49c3c21c943c66cd8796060df3
--- /dev/null
+++ b/src/ODE/ODESolvers/seulex/seulex.H
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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::seulex
+
+Description
+ An extrapolation-algorithm, based on the linearly implicit Euler method
+ with step size control and order selection.
+
+ The implementation is based on the SEULEX code in
+ \verbatim
+ "Solving Ordinary Differential Equations II: Stiff
+ and Differential-Algebraic Problems, second edition",
+ Hairer, E.,
+ Nørsett, S.,
+ Wanner, G.,
+ Springer-Verlag, Berlin. 1996.
+ \endverbatim
+
+SourceFiles
+ seulex.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef seulex_H
+#define seulex_H
+
+#include "ODESolver.H"
+
+#include "scalarFieldField.H"
+#include "scalarMatrices.H"
+#include "labelField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class seulex Declaration
+\*---------------------------------------------------------------------------*/
+
+class seulex
+:
+ public ODESolver
+{
+
+ // Private data
+
+ //- seulex constants
+ static const scalar EPS;
+
+ static const label KMAXX=12, IMAXX = KMAXX + 1;
+
+ static const scalar
+ STEPFAC1, STEPFAC2, STEPFAC3, STEPFAC4, STEPFAC5, KFAC1, KFAC2;
+
+
+ //- Reference to ODESystem
+ const ODESystem& ode_;
+
+ // Temporary fields
+ mutable label kTarg_;
+ mutable labelField nSeq_;
+ mutable scalarField cost_, dfdx_, factrl_;
+ mutable scalarRectangularMatrix table_, fSave_;
+ mutable scalarSquareMatrix dfdy_;
+ mutable scalar jacRedo_, theta_;
+ mutable bool calcJac_, dense_;
+ mutable scalarSquareMatrix a_, coeff_;
+
+
+ mutable scalar hnext_, hdid_;
+
+ mutable labelList pivotIndices_;
+
+ // Fields space for "step" function
+ mutable scalarField hOpt_ ,work_, ySaved_, ySequence_, scale_;
+
+ // Fields used in "dy" function
+ mutable scalarField del_, yTemp_, dyTemp_, delTemp_;
+
+
+ // Private Member Functions
+
+ void step(const scalar& htry, scalar& x, scalarField& y) const;
+
+ bool dy
+ (
+ const scalar& x,
+ scalarField& y,
+ const scalar htot,
+ const label k,
+ scalarField& yend,
+ label& ipt,
+ scalarField& scale
+ ) const;
+
+
+ void polyextr
+ (
+ const label k,
+ scalarRectangularMatrix& table,
+ scalarField& last
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("seulex");
+
+
+ // Constructors
+
+ //- Construct from ODE
+ seulex(const ODESystem& ode, const dictionary& dict);
+
+
+ // Member Functions
+ void solve
+ (
+ const ODESystem& ode,
+ scalar& x,
+ scalarField& y,
+ scalar& dxTry
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
index 3a47b78cacd0cb01d69195c2f6821bb7459182c1..925719ebcc318ef5e4b1c0edc2ed491ddd613f74 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C
@@ -59,10 +59,10 @@ Foam::basicSymmetryPointPatchField<Type>::basicSymmetryPointPatchField
const basicSymmetryPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
- const pointPatchFieldMapper&
+ const pointPatchFieldMapper& mapper
)
:
- pointPatchField<Type>(p, iF)
+ pointPatchField<Type>(ptf, p, iF, mapper)
{}
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
index a80e1421bbc0a9a2454972b3a62fe713689d116c..05b0c876d26361fc71d338d860f87e7ae5874e4c 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C
@@ -68,10 +68,10 @@ calculatedPointPatchField<Type>::calculatedPointPatchField
const calculatedPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
- const pointPatchFieldMapper&
+ const pointPatchFieldMapper& mapper
)
:
- pointPatchField<Type>(p, iF)
+ pointPatchField<Type>(ptf, p, iF, mapper)
{}
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
index dfae5f7d2839c30b2b54015eb91eb739675d32a4..3ad9e51a79cb9f228bac1e366e0cefc077535875 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C
@@ -61,10 +61,10 @@ coupledPointPatchField<Type>::coupledPointPatchField
const coupledPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
- const pointPatchFieldMapper&
+ const pointPatchFieldMapper& mapper
)
:
- pointPatchField<Type>(p, iF)
+ pointPatchField<Type>(ptf, p, iF, mapper)
{}
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
index d09bedf76bd0408b11973cb0179058b7b87a0089..d8fc9c43a38dd274cc63b909385ebac11dff39ff 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/value/valuePointPatchField.C
@@ -108,7 +108,7 @@ Foam::valuePointPatchField<Type>::valuePointPatchField
const pointPatchFieldMapper& mapper
)
:
- pointPatchField<Type>(p, iF),
+ pointPatchField<Type>(ptf, p, iF, mapper),
Field<Type>(ptf, mapper)
{}
diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
index a10d879179faba2f4ec1fb3a41819128b5f33c82..d6ce1e4cf1c53632caad45e7afcb2c511fed0991 100644
--- a/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C
@@ -61,10 +61,10 @@ zeroGradientPointPatchField<Type>::zeroGradientPointPatchField
const zeroGradientPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
- const pointPatchFieldMapper&
+ const pointPatchFieldMapper& mapper
)
:
- pointPatchField<Type>(p, iF)
+ pointPatchField<Type>(ptf, p, iF, mapper)
{}
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
index 8fd4b5bc0b3d83f6cf5c642cd78380005b304d7d..8838f0fff5b90fab5408eeda2738ded66e010b21 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C
@@ -77,10 +77,10 @@ emptyPointPatchField<Type>::emptyPointPatchField
const emptyPointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
- const pointPatchFieldMapper&
+ const pointPatchFieldMapper& mapper
)
:
- pointPatchField<Type>(p, iF)
+ pointPatchField<Type>(ptf, p, iF, mapper)
{
if (!isType<emptyPointPatch>(this->patch()))
{
diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
index eac29f3251cc74a4329d6499c48d7222441d8ba1..6e8fa53184f183c3a7bdca0a73bde90a6698c301 100644
--- a/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C
@@ -74,10 +74,10 @@ Foam::wedgePointPatchField<Type>::wedgePointPatchField
const wedgePointPatchField<Type>& ptf,
const pointPatch& p,
const DimensionedField<Type, pointMesh>& iF,
- const pointPatchFieldMapper&
+ const pointPatchFieldMapper& mapper
)
:
- pointPatchField<Type>(p, iF)
+ pointPatchField<Type>(ptf, p, iF, mapper)
{
if (!isType<wedgePointPatch>(this->patch()))
{
diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C
index 3e62d6a202093d98fa413b5ed6953efdef875f00..63eb0395652752665ac429cd18aa35bac90a268d 100644
--- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C
+++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C
@@ -63,6 +63,22 @@ pointPatchField<Type>::pointPatchField
{}
+template<class Type>
+Foam::pointPatchField<Type>::pointPatchField
+(
+ const pointPatchField<Type>& ptf,
+ const pointPatch& p,
+ const DimensionedField<Type, pointMesh>& iF,
+ const pointPatchFieldMapper&
+)
+:
+ patch_(p),
+ internalField_(iF),
+ updated_(false),
+ patchType_(ptf.patchType_)
+{}
+
+
template<class Type>
pointPatchField<Type>::pointPatchField
(
diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
index 439966188b40b8da4a9b9b2c90970f4a4d18d6f0..88e6bfecfe652d19a5d89e1cc6501facdf2bb8f3 100644
--- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
+++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H
@@ -167,6 +167,15 @@ public:
const dictionary&
);
+ //- Construct by mapping given patchField<Type> onto a new patch
+ pointPatchField
+ (
+ const pointPatchField<Type>&,
+ const pointPatch&,
+ const DimensionedField<Type, pointMesh>&,
+ const pointPatchFieldMapper&
+ );
+
//- Construct as copy
pointPatchField(const pointPatchField<Type>&);
diff --git a/src/OpenFOAM/matrices/solution/solution.C b/src/OpenFOAM/matrices/solution/solution.C
index 1a1c215444f78e026d8aa5eb40034d58562c582d..ace5564e970ec643df761a1727cdba94ee827e10 100644
--- a/src/OpenFOAM/matrices/solution/solution.C
+++ b/src/OpenFOAM/matrices/solution/solution.C
@@ -262,7 +262,9 @@ bool Foam::solution::relaxField(const word& name) const
{
if (debug)
{
- Info<< "Find variable relaxation factor for " << name << endl;
+ Info<< "Field relaxation factor for " << name
+ << " is " << (fieldRelaxDict_.found(name) ? "set" : "unset")
+ << endl;
}
return fieldRelaxDict_.found(name) || fieldRelaxDict_.found("default");
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
index 8ac86bdfa886f1c37929c6c2586931eea2b3951d..2a777908b829697e691c76aaac2567316c7b2743 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
@@ -79,7 +79,7 @@ inline bool Foam::pimpleControl::firstIter() const
inline bool Foam::pimpleControl::finalIter() const
{
- return converged_ || (corr_ == corrPISO_);
+ return converged_ || (corr_ == nCorrPIMPLE_);
}
diff --git a/src/postProcessing/functionObjects/forces/forces/forces.C b/src/postProcessing/functionObjects/forces/forces/forces.C
index c24fa72729d8bdae8a5c2ca3f1734188a034caad..74b892aa1e9627b5c744bf48977176d162251f12 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<Switch>("log", true);
+ log_ = dict.lookupOrDefault<Switch>("log", false);
if (log_)
{
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
index b02aaca39a021c439b9299497491af88758b3d95..60c2550381984d3ff06a47394de3211ecd088dc4 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C
@@ -740,7 +740,7 @@ Foam::scalar Foam::chemistryModel<CompType, ThermoType>::solve
RR_[i][celli] = (c[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
}
}
- Info << "deltaTMin " << deltaTMin << endl;
+
return deltaTMin;
}
diff --git a/tutorials/combustion/chemFoam/README b/tutorials/combustion/chemFoam/README
index 50e94283bcefb7958faaaf29660e35b9882b669b..2cd8c89273293ba04ca4114a55f86b1ac3c62feb 100644
--- a/tutorials/combustion/chemFoam/README
+++ b/tutorials/combustion/chemFoam/README
@@ -1,13 +1,11 @@
Stiff chemistry solver validation test cases
gri : GRI-Mech 3.0. CH4 combustion, 53 species, 325 reactions
-
h2 : H2 combustion, 10 species, 27 reactions
-
nc7h16 : n-Heptane combustion, 544 species, 2446 reactions
-
ic8h18 : iso-Octane combustion, 874 species, 3796 reactions
Results interpreted in 'validation' sub-folder, where OpenFOAM results
-are compared against those predicted by Senkin (CHEMKIN II)
+are compared against those predicted by CHEMKIN II.
+Overall the best performing ODE solver is seulex followed closely by rodas32.
diff --git a/tutorials/combustion/chemFoam/gri/constant/chemistryProperties b/tutorials/combustion/chemFoam/gri/constant/chemistryProperties
index e65509f3962c8ec4c430ea2fab6bf99729d409fe..bf097e0ba39679d2ec60ded2d59bb8ffe51fa227 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 seulex;
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 2c0182864e6b26ac74a17767dff98e0d97028087..b00cac494d4cf528afe6947d96595c7a8058ce49 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 seulex; //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 7528b62231ac1def9a2febfeae547ff4a9919274..bf097e0ba39679d2ec60ded2d59bb8ffe51fa227 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 seulex;
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 d7e58ddc143c14f35dbd33969ef4acdb6eba393d..a2a586199750a45e70f0b577a05c075a57b213e0 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 seulex;
absTol 1e-14;
- relTol 1e-4;
+ relTol 1e-1;
}
diff --git a/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement b/tutorials/multiphase/interDyMFoam/ras/floatingObject/0.org/pointDisplacement
index 33d0400057f51c5f7bf1a7e67ea576815ba76078..daff376a9d323b0fc015d085f6f15ccf2e51e110 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);
}
}