From ea064d9cb00182dec984dfadd6039ce7834d1929 Mon Sep 17 00:00:00 2001 From: Henry Date: Tue, 5 Nov 2013 17:32:26 +0000 Subject: [PATCH 01/12] pimpleControl: Corrected bug introduced on last change --- .../general/solutionControl/pimpleControl/pimpleControlI.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H index 8ac86bdfa8..2a777908b8 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_); } From 3da28cb4bc3085598b799704b03f26c7fe49711e Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 6 Nov 2013 09:55:46 +0000 Subject: [PATCH 02/12] interDyMFoam: support relaxation of p_rgh --- applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H | 2 ++ 1 file changed, 2 insertions(+) diff --git a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H index 4fbd2666b4..c30b361e72 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); From 34c274761c74c5931470cd6bb460a25e9c2a700d Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 6 Nov 2013 09:56:23 +0000 Subject: [PATCH 03/12] solution: added set/unset to field relaxation diagnostics --- src/OpenFOAM/matrices/solution/solution.C | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/OpenFOAM/matrices/solution/solution.C b/src/OpenFOAM/matrices/solution/solution.C index 1a1c215444..ace5564e97 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"); From d179d9b6b99b22dc0be5aac2f3b1239185460e46 Mon Sep 17 00:00:00 2001 From: Henry Date: Wed, 6 Nov 2013 23:19:03 +0000 Subject: [PATCH 04/12] 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 05/12] 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 06/12] 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 07/12] 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 08/12] 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); } } From 8cae53f10b3072692f7f26c76bd0fe2f553cba0e Mon Sep 17 00:00:00 2001 From: Sergio Ferraris Date: Thu, 7 Nov 2013 14:25:03 +0000 Subject: [PATCH 09/12] ENH: Adding seulex ODE --- src/ODE/Make/files | 1 + src/ODE/ODESolvers/seulex/seulex.C | 501 +++++++++++++++++++++++++++++ src/ODE/ODESolvers/seulex/seulex.H | 160 +++++++++ 3 files changed, 662 insertions(+) create mode 100644 src/ODE/ODESolvers/seulex/seulex.C create mode 100644 src/ODE/ODESolvers/seulex/seulex.H diff --git a/src/ODE/Make/files b/src/ODE/Make/files index a22ec1d3be..34ef0d4cb1 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -14,5 +14,6 @@ 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/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C new file mode 100644 index 0000000000..69227f1963 --- /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 . + +\*---------------------------------------------------------------------------*/ + +#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 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 1.0) + { + return false; + } + } + + delTemp_ = yend; + LUBacksubstitute(a_, pivotIndices_, delTemp_); + del_ = delTemp_; + + if (dense_ && nn >= nstep-k-1) + { + ipt++; + for (label i=0;i 0; j--) + { + for (label i=0; i. + +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 + +// ************************************************************************* // From b2f8d205b687ab29220ceba37f59f5d738bbb9c3 Mon Sep 17 00:00:00 2001 From: Sergio Ferraris Date: Thu, 7 Nov 2013 14:27:11 +0000 Subject: [PATCH 10/12] BUG: Deleting Info statement --- .../chemistryModel/chemistryModel/chemistryModel.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel/chemistryModel.C index b02aaca39a..60c2550381 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::solve RR_[i][celli] = (c[i] - c0[i])*specieThermo_[i].W()/deltaT[celli]; } } - Info << "deltaTMin " << deltaTMin << endl; + return deltaTMin; } From bc2d925d61f85d4d373e526ec93ccf9d614f009d Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 7 Nov 2013 18:38:58 +0000 Subject: [PATCH 11/12] chemFoam tutorials: set ODE solver to seulex which performs best overall --- tutorials/combustion/chemFoam/README | 6 ++---- .../combustion/chemFoam/gri/constant/chemistryProperties | 2 +- .../combustion/chemFoam/h2/constant/chemistryProperties | 2 +- .../combustion/chemFoam/ic8h18/constant/chemistryProperties | 2 +- .../combustion/chemFoam/nc7h16/constant/chemistryProperties | 2 +- 5 files changed, 6 insertions(+), 8 deletions(-) diff --git a/tutorials/combustion/chemFoam/README b/tutorials/combustion/chemFoam/README index 50e94283bc..2cd8c89273 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 314f83cc2d..bf097e0ba3 100644 --- a/tutorials/combustion/chemFoam/gri/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/gri/constant/chemistryProperties @@ -33,7 +33,7 @@ EulerImplicitCoeffs odeCoeffs { - solver rodas32; + solver seulex; absTol 1e-12; relTol 1e-1; } diff --git a/tutorials/combustion/chemFoam/h2/constant/chemistryProperties b/tutorials/combustion/chemFoam/h2/constant/chemistryProperties index 314f83cc2d..b00cac494d 100644 --- a/tutorials/combustion/chemFoam/h2/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/h2/constant/chemistryProperties @@ -33,7 +33,7 @@ EulerImplicitCoeffs odeCoeffs { - solver rodas32; + solver seulex; //rodas32; absTol 1e-12; relTol 1e-1; } diff --git a/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties b/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties index 314f83cc2d..bf097e0ba3 100644 --- a/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/ic8h18/constant/chemistryProperties @@ -33,7 +33,7 @@ EulerImplicitCoeffs odeCoeffs { - solver rodas32; + solver seulex; absTol 1e-12; relTol 1e-1; } diff --git a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties index aee69a3239..a2a5861997 100644 --- a/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties +++ b/tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties @@ -33,7 +33,7 @@ EulerImplicitCoeffs odeCoeffs { - solver rodas32; + solver seulex; absTol 1e-14; relTol 1e-1; } From 3dec6f010126b98c562bade503b3768baedcfee4 Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 7 Nov 2013 18:59:48 +0000 Subject: [PATCH 12/12] pointPatchField: Correct mapping so that the patchType entry in the base-class is mapped --- .../basicSymmetry/basicSymmetryPointPatchField.C | 4 ++-- .../basic/calculated/calculatedPointPatchField.C | 4 ++-- .../basic/coupled/coupledPointPatchField.C | 4 ++-- .../basic/value/valuePointPatchField.C | 2 +- .../zeroGradient/zeroGradientPointPatchField.C | 4 ++-- .../constraint/empty/emptyPointPatchField.C | 4 ++-- .../constraint/wedge/wedgePointPatchField.C | 4 ++-- .../pointPatchField/pointPatchField.C | 16 ++++++++++++++++ .../pointPatchField/pointPatchField.H | 9 +++++++++ 9 files changed, 38 insertions(+), 13 deletions(-) diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/basicSymmetry/basicSymmetryPointPatchField.C index 3a47b78cac..925719ebcc 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::basicSymmetryPointPatchField const basicSymmetryPointPatchField& ptf, const pointPatch& p, const DimensionedField& iF, - const pointPatchFieldMapper& + const pointPatchFieldMapper& mapper ) : - pointPatchField(p, iF) + pointPatchField(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 a80e1421bb..05b0c876d2 100644 --- a/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C +++ b/src/OpenFOAM/fields/pointPatchFields/basic/calculated/calculatedPointPatchField.C @@ -68,10 +68,10 @@ calculatedPointPatchField::calculatedPointPatchField const calculatedPointPatchField& ptf, const pointPatch& p, const DimensionedField& iF, - const pointPatchFieldMapper& + const pointPatchFieldMapper& mapper ) : - pointPatchField(p, iF) + pointPatchField(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 dfae5f7d28..3ad9e51a79 100644 --- a/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C +++ b/src/OpenFOAM/fields/pointPatchFields/basic/coupled/coupledPointPatchField.C @@ -61,10 +61,10 @@ coupledPointPatchField::coupledPointPatchField const coupledPointPatchField& ptf, const pointPatch& p, const DimensionedField& iF, - const pointPatchFieldMapper& + const pointPatchFieldMapper& mapper ) : - pointPatchField(p, iF) + pointPatchField(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 d09bedf76b..d8fc9c43a3 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::valuePointPatchField const pointPatchFieldMapper& mapper ) : - pointPatchField(p, iF), + pointPatchField(ptf, p, iF, mapper), Field(ptf, mapper) {} diff --git a/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C index a10d879179..d6ce1e4cf1 100644 --- a/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C +++ b/src/OpenFOAM/fields/pointPatchFields/basic/zeroGradient/zeroGradientPointPatchField.C @@ -61,10 +61,10 @@ zeroGradientPointPatchField::zeroGradientPointPatchField const zeroGradientPointPatchField& ptf, const pointPatch& p, const DimensionedField& iF, - const pointPatchFieldMapper& + const pointPatchFieldMapper& mapper ) : - pointPatchField(p, iF) + pointPatchField(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 8fd4b5bc0b..8838f0fff5 100644 --- a/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C +++ b/src/OpenFOAM/fields/pointPatchFields/constraint/empty/emptyPointPatchField.C @@ -77,10 +77,10 @@ emptyPointPatchField::emptyPointPatchField const emptyPointPatchField& ptf, const pointPatch& p, const DimensionedField& iF, - const pointPatchFieldMapper& + const pointPatchFieldMapper& mapper ) : - pointPatchField(p, iF) + pointPatchField(ptf, p, iF, mapper) { if (!isType(this->patch())) { diff --git a/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/constraint/wedge/wedgePointPatchField.C index eac29f3251..6e8fa53184 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::wedgePointPatchField const wedgePointPatchField& ptf, const pointPatch& p, const DimensionedField& iF, - const pointPatchFieldMapper& + const pointPatchFieldMapper& mapper ) : - pointPatchField(p, iF) + pointPatchField(ptf, p, iF, mapper) { if (!isType(this->patch())) { diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C index 3e62d6a202..63eb039565 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.C @@ -63,6 +63,22 @@ pointPatchField::pointPatchField {} +template +Foam::pointPatchField::pointPatchField +( + const pointPatchField& ptf, + const pointPatch& p, + const DimensionedField& iF, + const pointPatchFieldMapper& +) +: + patch_(p), + internalField_(iF), + updated_(false), + patchType_(ptf.patchType_) +{} + + template pointPatchField::pointPatchField ( diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointPatchField.H index 439966188b..88e6bfecfe 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 onto a new patch + pointPatchField + ( + const pointPatchField&, + const pointPatch&, + const DimensionedField&, + const pointPatchFieldMapper& + ); + //- Construct as copy pointPatchField(const pointPatchField&);