From de2ac5afd21ea52829fa922a69366e0d5751f308 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 6 Oct 2008 15:25:57 +0100 Subject: [PATCH 01/11] new systemCall function object, and added -lsampling to options files --- src/postProcessing/Allwmake | 1 + src/postProcessing/forces/Make/options | 1 + src/postProcessing/minMaxFields/Make/options | 3 +- src/postProcessing/systemCall/IOsystemCall.H | 50 ++++++ src/postProcessing/systemCall/Make/files | 4 + src/postProcessing/systemCall/Make/options | 9 ++ src/postProcessing/systemCall/systemCall.C | 91 +++++++++++ src/postProcessing/systemCall/systemCall.H | 147 ++++++++++++++++++ .../systemCall/systemCallFunctionObject.C | 43 +++++ .../systemCall/systemCallFunctionObject.H | 55 +++++++ 10 files changed, 403 insertions(+), 1 deletion(-) create mode 100644 src/postProcessing/systemCall/IOsystemCall.H create mode 100644 src/postProcessing/systemCall/Make/files create mode 100644 src/postProcessing/systemCall/Make/options create mode 100644 src/postProcessing/systemCall/systemCall.C create mode 100644 src/postProcessing/systemCall/systemCall.H create mode 100644 src/postProcessing/systemCall/systemCallFunctionObject.C create mode 100644 src/postProcessing/systemCall/systemCallFunctionObject.H diff --git a/src/postProcessing/Allwmake b/src/postProcessing/Allwmake index 1d889e8b76..e1f7cc56d8 100755 --- a/src/postProcessing/Allwmake +++ b/src/postProcessing/Allwmake @@ -7,5 +7,6 @@ wmake libso forces wmake libso fieldAverage wmake libso foamCalcFunctions wmake libso minMaxFields +wmake libso systemCall # ----------------------------------------------------------------- end-of-file diff --git a/src/postProcessing/forces/Make/options b/src/postProcessing/forces/Make/options index 929d4c25d1..3ae6adad4e 100644 --- a/src/postProcessing/forces/Make/options +++ b/src/postProcessing/forces/Make/options @@ -11,6 +11,7 @@ EXE_INC = \ LIB_LIBS = \ -lfiniteVolume \ -lmeshTools \ + -lsampling \ -lincompressibleTransportModels \ -lincompressibleRASModels \ -lincompressibleLESModels \ diff --git a/src/postProcessing/minMaxFields/Make/options b/src/postProcessing/minMaxFields/Make/options index 0954a90689..5166bcc9e3 100644 --- a/src/postProcessing/minMaxFields/Make/options +++ b/src/postProcessing/minMaxFields/Make/options @@ -5,4 +5,5 @@ EXE_INC = \ LIB_LIBS = \ -lfiniteVolume \ - -lmeshTools + -lmeshTools \ + -lsampling diff --git a/src/postProcessing/systemCall/IOsystemCall.H b/src/postProcessing/systemCall/IOsystemCall.H new file mode 100644 index 0000000000..fff6e7dbdb --- /dev/null +++ b/src/postProcessing/systemCall/IOsystemCall.H @@ -0,0 +1,50 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Typedef + Foam::IOsystemCall + +Description + Instance of the generic IOOutputFilter for systemCall. + +\*---------------------------------------------------------------------------*/ + +#ifndef IOsystemCall_H +#define IOsystemCall_H + +#include "systemCall.H" +#include "IOOutputFilter.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + typedef IOOutputFilter IOsystemCall; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/systemCall/Make/files b/src/postProcessing/systemCall/Make/files new file mode 100644 index 0000000000..20bafa8dfa --- /dev/null +++ b/src/postProcessing/systemCall/Make/files @@ -0,0 +1,4 @@ +systemCall.C +systemCallFunctionObject.C + +LIB = $(FOAM_LIBBIN)/libsystemCall diff --git a/src/postProcessing/systemCall/Make/options b/src/postProcessing/systemCall/Make/options new file mode 100644 index 0000000000..5166bcc9e3 --- /dev/null +++ b/src/postProcessing/systemCall/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +LIB_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling diff --git a/src/postProcessing/systemCall/systemCall.C b/src/postProcessing/systemCall/systemCall.C new file mode 100644 index 0000000000..b1756776c5 --- /dev/null +++ b/src/postProcessing/systemCall/systemCall.C @@ -0,0 +1,91 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "systemCall.H" +#include "dictionary.H" +#include "Time.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(systemCall, 0); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::systemCall::systemCall +( + const word& name, + const objectRegistry& obr, + const dictionary& dict, + const bool loadFromFiles +) +: + name_(name), + obr_(obr), + active_(true), + executeCalls_(), + writeCalls_() +{ + read(dict); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::systemCall::~systemCall() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::systemCall::read(const dictionary& dict) +{ + dict.lookup("executeCalls") >> executeCalls_; + dict.lookup("writeCalls") >> writeCalls_; +} + + +void Foam::systemCall::execute() +{ + forAll(executeCalls_, callI) + { + ::system(executeCalls_[callI].c_str()); + } +} + +void Foam::systemCall::write() +{ + forAll(writeCalls_, callI) + { + ::system(writeCalls_[callI].c_str()); + } +} + + +// ************************************************************************* // diff --git a/src/postProcessing/systemCall/systemCall.H b/src/postProcessing/systemCall/systemCall.H new file mode 100644 index 0000000000..1d070cd6a9 --- /dev/null +++ b/src/postProcessing/systemCall/systemCall.H @@ -0,0 +1,147 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::systemCall + +Description + Executes system calls, entered in the form of a string list + +SourceFiles + systemCall.C + IOsystemCall.H + +\*---------------------------------------------------------------------------*/ + +#ifndef systemCall_H +#define systemCall_H + +#include "stringList.H" +#include "pointFieldFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class objectRegistry; +class dictionary; +class mapPolyMesh; + +/*---------------------------------------------------------------------------*\ + Class systemCall Declaration +\*---------------------------------------------------------------------------*/ + +class systemCall +{ +protected: + + // Private data + + //- Name of this set of forces, + // Also used as the name of the probes directory. + word name_; + + const objectRegistry& obr_; + + //- on/off switch + bool active_; + + //- List of calls to execute - every step + stringList executeCalls_; + + //- List of calls to execute - write steps + stringList writeCalls_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + systemCall(const systemCall&); + + //- Disallow default bitwise assignment + void operator=(const systemCall&); + + +public: + + //- Runtime type information + TypeName("systemCall"); + + + // Constructors + + //- Construct for given objectRegistry and dictionary. + // Allow the possibility to load fields from files + systemCall + ( + const word& name, + const objectRegistry&, + const dictionary&, + const bool loadFromFiles = false + ); + + + // Destructor + + virtual ~systemCall(); + + + // Member Functions + + //- Return name of the system call set + virtual const word& name() const + { + return name_; + } + + //- Read the system calls + virtual void read(const dictionary&); + + //- Execute + virtual void execute(); + + //- Write + virtual void write(); + + //- Update for changes of mesh + virtual void updateMesh(const mapPolyMesh&) + {} + + //- Update for changes of mesh + virtual void movePoints(const pointField&) + {} +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/postProcessing/systemCall/systemCallFunctionObject.C b/src/postProcessing/systemCall/systemCallFunctionObject.C new file mode 100644 index 0000000000..6e6b2e646b --- /dev/null +++ b/src/postProcessing/systemCall/systemCallFunctionObject.C @@ -0,0 +1,43 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "systemCallFunctionObject.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineNamedTemplateTypeNameAndDebug(systemCallFunctionObject, 0); + + addToRunTimeSelectionTable + ( + functionObject, + systemCallFunctionObject, + dictionary + ); +} + +// ************************************************************************* // diff --git a/src/postProcessing/systemCall/systemCallFunctionObject.H b/src/postProcessing/systemCall/systemCallFunctionObject.H new file mode 100644 index 0000000000..6f8b6ba635 --- /dev/null +++ b/src/postProcessing/systemCall/systemCallFunctionObject.H @@ -0,0 +1,55 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Typedef + Foam::systemCallFunctionObject + +Description + FunctionObject wrapper around systemCall to allow them to be created via + the functions list within controlDict. + +SourceFiles + systemCallFunctionObject.C + +\*---------------------------------------------------------------------------*/ + +#ifndef systemCallFunctionObject_H +#define systemCallFunctionObject_H + +#include "systemCall.H" +#include "OutputFilterFunctionObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + typedef OutputFilterFunctionObject + systemCallFunctionObject; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From 45422cda59e62450cbf6b4d0ce3c61ba1cd2d230 Mon Sep 17 00:00:00 2001 From: henry Date: Thu, 9 Oct 2008 11:23:05 +0100 Subject: [PATCH 02/11] Changed the access to derived properties to ensure the old-time level values are stored if required. --- .../basic/hThermo/hThermo.C | 54 ++++++------ .../hMixtureThermo/hMixtureThermo.C | 52 ++++++------ .../hhuMixtureThermo/hhuMixtureThermo.C | 82 +++++++++++-------- 3 files changed, 101 insertions(+), 87 deletions(-) diff --git a/src/thermophysicalModels/basic/hThermo/hThermo.C b/src/thermophysicalModels/basic/hThermo/hThermo.C index 19e43d2410..7e32f7a063 100644 --- a/src/thermophysicalModels/basic/hThermo/hThermo.C +++ b/src/thermophysicalModels/basic/hThermo/hThermo.C @@ -28,15 +28,10 @@ License #include "fvMesh.H" #include "fixedValueFvPatchFields.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -hThermo::hThermo(const fvMesh& mesh) +Foam::hThermo::hThermo(const fvMesh& mesh) : basicThermo(mesh), MixtureType(*this, mesh), @@ -56,9 +51,12 @@ hThermo::hThermo(const fvMesh& mesh) hBoundaryTypes() ) { - forAll(h_, celli) + scalarField& hCells = h_.internalField(); + const scalarField& TCells = T_.internalField(); + + forAll(hCells, celli) { - h_[celli] = this->cellMixture(celli).H(T_[celli]); + hCells[celli] = this->cellMixture(celli).H(TCells[celli]); } forAll(h_.boundaryField(), patchi) @@ -76,25 +74,33 @@ hThermo::hThermo(const fvMesh& mesh) // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template -hThermo::~hThermo() +Foam::hThermo::~hThermo() {} // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template -void hThermo::calculate() +void Foam::hThermo::calculate() { - forAll(T_, celli) + const scalarField& hCells = h_.internalField(); + const scalarField& pCells = p_.internalField(); + + scalarField& TCells = T_.internalField(); + scalarField& psiCells = psi_.internalField(); + scalarField& muCells = mu_.internalField(); + scalarField& alphaCells = alpha_.internalField(); + + forAll(TCells, celli) { const typename MixtureType::thermoType& mixture_ = this->cellMixture(celli); - T_[celli] = mixture_.TH(h_[celli], T_[celli]); - psi_[celli] = mixture_.psi(p_[celli], T_[celli]); + TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]); + psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]); - mu_[celli] = mixture_.mu(T_[celli]); - alpha_[celli] = mixture_.alpha(T_[celli]); + muCells[celli] = mixture_.mu(TCells[celli]); + alphaCells[celli] = mixture_.alpha(TCells[celli]); } forAll(T_.boundaryField(), patchi) @@ -143,7 +149,7 @@ void hThermo::calculate() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -void hThermo::correct() +void Foam::hThermo::correct() { if (debug) { @@ -163,7 +169,7 @@ void hThermo::correct() template -tmp hThermo::h +Foam::tmp Foam::hThermo::h ( const scalarField& T, const labelList& cells @@ -182,7 +188,7 @@ tmp hThermo::h template -tmp hThermo::h +Foam::tmp Foam::hThermo::h ( const scalarField& T, const label patchi @@ -201,7 +207,7 @@ tmp hThermo::h template -tmp hThermo::Cp +Foam::tmp Foam::hThermo::Cp ( const scalarField& T, const label patchi @@ -220,7 +226,7 @@ tmp hThermo::Cp template -tmp hThermo::Cp() const +Foam::tmp Foam::hThermo::Cp() const { const fvMesh& mesh = T_.mesh(); @@ -258,7 +264,7 @@ tmp hThermo::Cp() const template -tmp hThermo::Cv() const +Foam::tmp Foam::hThermo::Cv() const { const fvMesh& mesh = T_.mesh(); @@ -303,7 +309,7 @@ tmp hThermo::Cv() const template -bool hThermo::read() +bool Foam::hThermo::read() { if (basicThermo::read()) { @@ -317,8 +323,4 @@ bool hThermo::read() } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - // ************************************************************************* // diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C index 2d72a4f688..064e7c3eba 100644 --- a/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C +++ b/src/thermophysicalModels/combustion/mixtureThermos/hMixtureThermo/hMixtureThermo.C @@ -28,22 +28,20 @@ License #include "fvMesh.H" #include "fixedValueFvPatchFields.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -hMixtureThermo::hMixtureThermo(const fvMesh& mesh) +Foam::hMixtureThermo::hMixtureThermo(const fvMesh& mesh) : hCombustionThermo(mesh), MixtureType(*this, mesh) { - forAll(h_, celli) + scalarField& hCells = h_.internalField(); + const scalarField& TCells = T_.internalField(); + + forAll(hCells, celli) { - h_[celli] = this->cellMixture(celli).H(T_[celli]); + hCells[celli] = this->cellMixture(celli).H(TCells[celli]); } forAll(h_.boundaryField(), patchi) @@ -61,25 +59,33 @@ hMixtureThermo::hMixtureThermo(const fvMesh& mesh) // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template -hMixtureThermo::~hMixtureThermo() +Foam::hMixtureThermo::~hMixtureThermo() {} // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template -void hMixtureThermo::calculate() +void Foam::hMixtureThermo::calculate() { - forAll(T_, celli) + const scalarField& hCells = h_.internalField(); + const scalarField& pCells = p_.internalField(); + + scalarField& TCells = T_.internalField(); + scalarField& psiCells = psi_.internalField(); + scalarField& muCells = mu_.internalField(); + scalarField& alphaCells = alpha_.internalField(); + + forAll(TCells, celli) { const typename MixtureType::thermoType& mixture_ = this->cellMixture(celli); - T_[celli] = mixture_.TH(h_[celli], T_[celli]); - psi_[celli] = mixture_.psi(p_[celli], T_[celli]); + TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]); + psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]); - mu_[celli] = mixture_.mu(T_[celli]); - alpha_[celli] = mixture_.alpha(T_[celli]); + muCells[celli] = mixture_.mu(TCells[celli]); + alphaCells[celli] = mixture_.alpha(TCells[celli]); } forAll(T_.boundaryField(), patchi) @@ -128,7 +134,7 @@ void hMixtureThermo::calculate() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -tmp hMixtureThermo::h +Foam::tmp Foam::hMixtureThermo::h ( const scalarField& T, const labelList& cells @@ -147,7 +153,7 @@ tmp hMixtureThermo::h template -tmp hMixtureThermo::h +Foam::tmp Foam::hMixtureThermo::h ( const scalarField& T, const label patchi @@ -166,7 +172,7 @@ tmp hMixtureThermo::h template -tmp hMixtureThermo::Cp +Foam::tmp Foam::hMixtureThermo::Cp ( const scalarField& T, const label patchi @@ -186,7 +192,7 @@ tmp hMixtureThermo::Cp template -tmp hMixtureThermo::Cp() const +Foam::tmp Foam::hMixtureThermo::Cp() const { const fvMesh& mesh = T_.mesh(); @@ -224,7 +230,7 @@ tmp hMixtureThermo::Cp() const template -void hMixtureThermo::correct() +void Foam::hMixtureThermo::correct() { if (debug) { @@ -244,7 +250,7 @@ void hMixtureThermo::correct() template -bool hMixtureThermo::read() +bool Foam::hMixtureThermo::read() { if (hCombustionThermo::read()) { @@ -258,8 +264,4 @@ bool hMixtureThermo::read() } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - // ************************************************************************* // diff --git a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C index 04bf08e330..0086e8cd14 100644 --- a/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C +++ b/src/thermophysicalModels/combustion/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C @@ -28,23 +28,23 @@ License #include "fvMesh.H" #include "fixedValueFvPatchFields.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -hhuMixtureThermo::hhuMixtureThermo(const fvMesh& mesh) +Foam::hhuMixtureThermo::hhuMixtureThermo(const fvMesh& mesh) : hhuCombustionThermo(mesh), MixtureType(*this, mesh) { - forAll(h_, celli) + scalarField& hCells = h_.internalField(); + scalarField& huCells = hu_.internalField(); + const scalarField& TCells = T_.internalField(); + const scalarField& TuCells = Tu_.internalField(); + + forAll(hCells, celli) { - h_[celli] = this->cellMixture(celli).H(T_[celli]); - hu_[celli] = this->cellReactants(celli).H(Tu_[celli]); + hCells[celli] = this->cellMixture(celli).H(TCells[celli]); + huCells[celli] = this->cellReactants(celli).H(TuCells[celli]); } forAll(h_.boundaryField(), patchi) @@ -71,27 +71,38 @@ hhuMixtureThermo::hhuMixtureThermo(const fvMesh& mesh) // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template -hhuMixtureThermo::~hhuMixtureThermo() +Foam::hhuMixtureThermo::~hhuMixtureThermo() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -void hhuMixtureThermo::calculate() +void Foam::hhuMixtureThermo::calculate() { - forAll(T_, celli) + const scalarField& hCells = h_.internalField(); + const scalarField& huCells = hu_.internalField(); + const scalarField& pCells = p_.internalField(); + + scalarField& TCells = T_.internalField(); + scalarField& TuCells = Tu_.internalField(); + scalarField& psiCells = psi_.internalField(); + scalarField& muCells = mu_.internalField(); + scalarField& alphaCells = alpha_.internalField(); + + forAll(TCells, celli) { - const typename MixtureType::thermoType& mixture_ = + const typename MixtureType::thermoType& mixture_ = this->cellMixture(celli); - T_[celli] = mixture_.TH(h_[celli], T_[celli]); - psi_[celli] = mixture_.psi(p_[celli], T_[celli]); + TCells[celli] = mixture_.TH(hCells[celli], TCells[celli]); + psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]); - mu_[celli] = mixture_.mu(T_[celli]); - alpha_[celli] = mixture_.alpha(T_[celli]); + muCells[celli] = mixture_.mu(TCells[celli]); + alphaCells[celli] = mixture_.alpha(TCells[celli]); - Tu_[celli] = this->cellReactants(celli).TH(hu_[celli], Tu_[celli]); + TuCells[celli] = + this->cellReactants(celli).TH(huCells[celli], TuCells[celli]); } forAll(T_.boundaryField(), patchi) @@ -144,7 +155,7 @@ void hhuMixtureThermo::calculate() template -void hhuMixtureThermo::correct() +void Foam::hhuMixtureThermo::correct() { if (debug) { @@ -163,7 +174,7 @@ void hhuMixtureThermo::correct() } template -tmp hhuMixtureThermo::h +Foam::tmp Foam::hhuMixtureThermo::h ( const scalarField& T, const labelList& cells @@ -182,7 +193,7 @@ tmp hhuMixtureThermo::h template -tmp hhuMixtureThermo::h +Foam::tmp Foam::hhuMixtureThermo::h ( const scalarField& T, const label patchi @@ -201,7 +212,7 @@ tmp hhuMixtureThermo::h template -tmp hhuMixtureThermo::Cp +Foam::tmp Foam::hhuMixtureThermo::Cp ( const scalarField& T, const label patchi @@ -221,7 +232,7 @@ tmp hhuMixtureThermo::Cp template -tmp hhuMixtureThermo::Cp() const +Foam::tmp Foam::hhuMixtureThermo::Cp() const { const fvMesh& mesh = T_.mesh(); @@ -260,7 +271,7 @@ tmp hhuMixtureThermo::Cp() const template -tmp hhuMixtureThermo::hu +Foam::tmp Foam::hhuMixtureThermo::hu ( const scalarField& Tu, const labelList& cells @@ -279,7 +290,7 @@ tmp hhuMixtureThermo::hu template -tmp hhuMixtureThermo::hu +Foam::tmp Foam::hhuMixtureThermo::hu ( const scalarField& Tu, const label patchi @@ -298,7 +309,7 @@ tmp hhuMixtureThermo::hu template -tmp hhuMixtureThermo::Tb() const +Foam::tmp Foam::hhuMixtureThermo::Tb() const { tmp tTb ( @@ -342,7 +353,8 @@ tmp hhuMixtureThermo::Tb() const template -tmp hhuMixtureThermo::psiu() const +Foam::tmp +Foam::hhuMixtureThermo::psiu() const { tmp tpsiu ( @@ -388,7 +400,8 @@ tmp hhuMixtureThermo::psiu() const template -tmp hhuMixtureThermo::psib() const +Foam::tmp +Foam::hhuMixtureThermo::psib() const { tmp tpsib ( @@ -426,7 +439,8 @@ tmp hhuMixtureThermo::psib() const forAll(ppsib, facei) { ppsib[facei] = - this->patchFaceReactants(patchi, facei).psi(pp[facei], pTb[facei]); + this->patchFaceReactants + (patchi, facei).psi(pp[facei], pTb[facei]); } } @@ -435,7 +449,7 @@ tmp hhuMixtureThermo::psib() const template -tmp hhuMixtureThermo::muu() const +Foam::tmp Foam::hhuMixtureThermo::muu() const { tmp tmuu ( @@ -478,7 +492,7 @@ tmp hhuMixtureThermo::muu() const template -tmp hhuMixtureThermo::mub() const +Foam::tmp Foam::hhuMixtureThermo::mub() const { tmp tmub ( @@ -521,7 +535,7 @@ tmp hhuMixtureThermo::mub() const template -bool hhuMixtureThermo::read() +bool Foam::hhuMixtureThermo::read() { if (hhuCombustionThermo::read()) { @@ -535,8 +549,4 @@ bool hhuMixtureThermo::read() } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - // ************************************************************************* // From b4263cca93beac0a5996dc6cfc7b8dba723581cd Mon Sep 17 00:00:00 2001 From: henry Date: Thu, 9 Oct 2008 14:19:24 +0100 Subject: [PATCH 03/11] Added more aggressive limiting of the quadratic correction. --- .../schemes/quadraticFit/quadraticFitData.C | 76 ++++++++++++++----- 1 file changed, 57 insertions(+), 19 deletions(-) diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C index 7271bad0ff..79f9107e97 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/quadraticFit/quadraticFitData.C @@ -113,7 +113,9 @@ Foam::quadraticFitData::quadraticFitData { interpPolySize[faci] = calcFit(stencilPoints[faci], faci); } + interpPolySize.write(); + if (debug) { Info<< "quadraticFitData::quadraticFitData() :" @@ -206,7 +208,7 @@ Foam::label Foam::quadraticFitData::calcFit // calculate the matrix of the polynomial components scalarRectangularMatrix B(C.size(), minSize_, scalar(0)); - for(label ip = 0; ip < C.size(); ip++) + for(label ip = 0; ip < C.size(); ip++) { const point& p = C[ip]; @@ -228,6 +230,7 @@ Foam::label Foam::quadraticFitData::calcFit pz /= scale; label is = 0; + B[ip][is++] = wts[0]*wts[ip]; B[ip][is++] = wts[0]*wts[ip]*px; B[ip][is++] = wts[ip]*sqr(px); @@ -253,21 +256,36 @@ Foam::label Foam::quadraticFitData::calcFit scalarList singVals(minSize_); label nSVDzeros = 0; + const GeometricField& w = + mesh().surfaceInterpolation::weights(); + bool goodFit = false; - for(int iIt = 0; iIt < 35 && !goodFit; iIt++) + for(int iIt = 0; iIt < 10 && !goodFit; iIt++) { SVD svd(B, SMALL); scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[0][0]; scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[0][1]; - goodFit = sign(fit0) == sign(fit1) && fit0 < 1 && fit1 < 1; + //goodFit = (fit0 > 0 && fit1 > 0); + + goodFit = + (mag(fit0 - w[faci])/w[faci] < 0.5) + && (mag(fit1 - (1 - w[faci]))/(1 - w[faci]) < 0.5); + + //scalar w0Err = fit0/w[faci]; + //scalar w1Err = fit1/(1 - w[faci]); + + //goodFit = + // (w0Err > 0.5 && w0Err < 1.5) + // && (w1Err > 0.5 && w1Err < 1.5); if (goodFit) { fit_[faci][0] = fit0; fit_[faci][1] = fit1; - for(label i = 2; i < stencilSize; i++) + + for(label i=2; i& w = - mesh().surfaceInterpolation::weights(); - - // remove the uncorrected linear coefficients - - fit_[faci][0] -= w[faci]; - fit_[faci][1] -= 1 - w[faci]; return minSize_ - nSVDzeros; } From b8a5d2ef595571ba7c0f8d63b7dba2f01c60fe18 Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 10 Oct 2008 14:55:46 +0100 Subject: [PATCH 04/11] minor cosmetics change --- .../parcels/Templates/ReactingParcel/ReactingParcel.C | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C index 5fe3aa1633..d7b2f2b969 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C @@ -57,6 +57,7 @@ void Foam::ReactingParcel::calcCoupled // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ const vector U0 = this->U_; const scalar mass0 = this->mass(); + const scalar cp0 = this->cp_; const scalar np0 = this->nParticle_; const scalar T0 = this->T_; @@ -116,11 +117,11 @@ void Foam::ReactingParcel::calcCoupled T1 -= td.constProps().Ldevol() *sum(dMassMT) - /(0.5*(mass0 + mass1)*this->cp_); + /(0.5*(mass0 + mass1)*cp0); // Add retained enthalpy from surface reaction to particle and remove // from gas - T1 += dhRet/(0.5*(mass0 + mass1)*this->cp_); + T1 += dhRet/(0.5*(mass0 + mass1)*cp0); dhTrans -= dhRet; // Correct dhTrans to account for enthalpy of evolved volatiles From 8293b56e179c806ca188b5f36161dd23968b7a2a Mon Sep 17 00:00:00 2001 From: andy Date: Fri, 10 Oct 2008 16:52:10 +0100 Subject: [PATCH 05/11] consistency update --- src/postProcessing/forces/forces/forces.C | 18 +++++++++--------- src/postProcessing/forces/forces/forces.H | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/postProcessing/forces/forces/forces.C b/src/postProcessing/forces/forces/forces.C index 17254aec3d..d3db7ed569 100644 --- a/src/postProcessing/forces/forces/forces.C +++ b/src/postProcessing/forces/forces/forces.C @@ -81,7 +81,7 @@ Foam::tmp Foam::forces::devRhoReff() const const basicThermo& thermo = obr_.lookupObject("thermophysicalProperties"); - const volVectorField& U = obr_.lookupObject(Uname_); + const volVectorField& U = obr_.lookupObject(UName_); return -thermo.mu()*dev(twoSymm(fvc::grad(U))); } @@ -94,7 +94,7 @@ Foam::tmp Foam::forces::devRhoReff() const obr_.lookupObject ("transportProperties"); - const volVectorField& U = obr_.lookupObject(Uname_); + const volVectorField& U = obr_.lookupObject(UName_); return -rhoRef_*laminarT.nu()*dev(twoSymm(fvc::grad(U))); } @@ -105,7 +105,7 @@ Foam::tmp Foam::forces::devRhoReff() const dimensionedScalar nu(transportProperties.lookup("nu")); - const volVectorField& U = obr_.lookupObject(Uname_); + const volVectorField& U = obr_.lookupObject(UName_); return -rhoRef_*nu*dev(twoSymm(fvc::grad(U))); } @@ -149,7 +149,7 @@ Foam::forces::forces log_(false), patchSet_(), pName_(""), - Uname_(""), + UName_(""), rhoRef_(0), CofR_(vector::zero), forcesFilePtr_(NULL) @@ -190,18 +190,18 @@ void Foam::forces::read(const dictionary& dict) // Optional entries U and p pName_ = dict.lookupOrDefault("pName", "p"); - Uname_ = dict.lookupOrDefault("Uname", "U"); + UName_ = dict.lookupOrDefault("UName", "U"); - // Check whether Uname and pName exists, if not deactivate forces + // Check whether UName and pName exists, if not deactivate forces if ( - !obr_.foundObject(Uname_) + !obr_.foundObject(UName_) || !obr_.foundObject(pName_) ) { active_ = false; WarningIn("void forces::read(const dictionary& dict)") - << "Could not find " << Uname_ << " or " + << "Could not find " << UName_ << " or " << pName_ << " in database." << nl << " De-activating forces." << endl; @@ -299,7 +299,7 @@ void Foam::forces::write() Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const { - const volVectorField& U = obr_.lookupObject(Uname_); + const volVectorField& U = obr_.lookupObject(UName_); const volScalarField& p = obr_.lookupObject(pName_); const fvMesh& mesh = U.mesh(); diff --git a/src/postProcessing/forces/forces/forces.H b/src/postProcessing/forces/forces/forces.H index 285a0acea1..35b1bdf034 100644 --- a/src/postProcessing/forces/forces/forces.H +++ b/src/postProcessing/forces/forces/forces.H @@ -130,7 +130,7 @@ protected: word pName_; //- Name of velocity field - word Uname_; + word UName_; //- Reference density needed for incompressible calculations scalar rhoRef_; From 4e4f9e0a0a1518fa4124ea8955fe0201f6b8e7ab Mon Sep 17 00:00:00 2001 From: henry Date: Sun, 12 Oct 2008 11:42:20 +0100 Subject: [PATCH 06/11] Added quadraticFit snGrad. --- .../quadraticFitSnGrad/quadraticFitSnGrad.H | 161 +++++++++ .../quadraticFitSnGradData.C | 325 ++++++++++++++++++ .../quadraticFitSnGradData.H | 131 +++++++ .../quadraticFitSnGrad/quadraticFitSnGrads.C | 44 +++ 4 files changed, 661 insertions(+) create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H create mode 100644 src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H new file mode 100644 index 0000000000..b74c3d2f20 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrad.H @@ -0,0 +1,161 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + quadraticFitSnGrad + +Description + Simple central-difference snGrad scheme with quadratic fit correction from + a larger stencil. + +SourceFiles + quadraticFitSnGrad.C + +\*---------------------------------------------------------------------------*/ + +#ifndef quadraticFitSnGrad_H +#define quadraticFitSnGrad_H + +#include "snGradScheme.H" +#include "quadraticFitSnGradData.H" +#include "extendedStencil.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fv +{ + +/*---------------------------------------------------------------------------*\ + Class quadraticFitSnGrad Declaration +\*---------------------------------------------------------------------------*/ + +template +class quadraticFitSnGrad +: + public snGradScheme +{ + // Private Data + //- weights for central stencil + const scalar centralWeight_; + + // Private Member Functions + + //- Disallow default bitwise assignment + void operator=(const quadraticFitSnGrad&); + + +public: + + //- Runtime type information + TypeName("quadraticFit"); + + + // Constructors + + //- Construct from mesh and scheme data + quadraticFitSnGrad + ( + const fvMesh& mesh, + const scalar centralWeight + ) + : + snGradScheme(mesh), + centralWeight_(centralWeight) + {} + + + //- Construct from mesh and data stream + quadraticFitSnGrad(const fvMesh& mesh, Istream& is) + : + snGradScheme(mesh), + centralWeight_(readScalar(is)) + {} + + + // Destructor + + virtual ~quadraticFitSnGrad() {} + + + // Member Functions + + //- Return the interpolation weighting factors for the given field + virtual tmp deltaCoeffs + ( + const GeometricField& + ) const + { + return this->mesh().deltaCoeffs(); + } + + //- Return true if this scheme uses an explicit correction + virtual bool corrected() const + { + return true; + } + + //- Return the explicit correction to the quadraticFitSnGrad + // for the given field + virtual tmp > + correction(const GeometricField& vf) const + { + const fvMesh& mesh = this->mesh(); + + const quadraticFitSnGradData& qfd = quadraticFitSnGradData::New + ( + mesh, + centralWeight_ + ); + + const extendedStencil& stencil = qfd.stencil(); + const List& f = qfd.fit(); + + tmp > sft + = stencil.weightedSum(vf, f); + + sft().dimensions() /= dimLength; + + return sft; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fv + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C new file mode 100644 index 0000000000..133f1e7933 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.C @@ -0,0 +1,325 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "quadraticFitSnGradData.H" +#include "surfaceFields.H" +#include "volFields.H" +#include "SVD.H" +#include "syncTools.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(quadraticFitSnGradData, 0); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::quadraticFitSnGradData::quadraticFitSnGradData +( + const fvMesh& mesh, + const scalar cWeight +) +: + MeshObject(mesh), + centralWeight_(cWeight), + #ifdef SPHERICAL_GEOMETRY + dim_(2), + #else + dim_(mesh.nGeometricD()), + #endif + minSize_ + ( + dim_ == 1 ? 3 : + dim_ == 2 ? 6 : + dim_ == 3 ? 9 : 0 + ), + stencil_(mesh), + fit_(mesh.nInternalFaces()) +{ + if (debug) + { + Info << "Contructing quadraticFitSnGradData" << endl; + } + + // check input + if (centralWeight_ < 1 - SMALL) + { + FatalErrorIn("quadraticFitSnGradData::quadraticFitSnGradData") + << "centralWeight requested = " << centralWeight_ + << " should not be less than one" + << exit(FatalError); + } + + if (minSize_ == 0) + { + FatalErrorIn("quadraticFitSnGradData") + << " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError); + } + + // store the polynomial size for each face to write out + surfaceScalarField snGradPolySize + ( + IOobject + ( + "quadraticFitSnGradPolySize", + "constant", + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("quadraticFitSnGradPolySize", dimless, scalar(0)) + ); + + // Get the cell/face centres in stencil order. + // Centred face stencils no good for triangles of tets. Need bigger stencils + List > stencilPoints(stencil_.stencil().size()); + stencil_.collectData + ( + mesh.C(), + stencilPoints + ); + + // find the fit coefficients for every face in the mesh + + for(label faci = 0; faci < mesh.nInternalFaces(); faci++) + { + snGradPolySize[faci] = calcFit(stencilPoints[faci], faci); + } + + if (debug) + { + snGradPolySize.write(); + Info<< "quadraticFitSnGradData::quadraticFitSnGradData() :" + << "Finished constructing polynomialFit data" + << endl; + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::quadraticFitSnGradData::findFaceDirs +( + vector& idir, // value changed in return + vector& jdir, // value changed in return + vector& kdir, // value changed in return + const fvMesh& mesh, + const label faci +) +{ + idir = mesh.Sf()[faci]; + idir /= mag(idir); + + #ifndef SPHERICAL_GEOMETRY + if (mesh.nGeometricD() <= 2) // find the normal direcion + { + if (mesh.directions()[0] == -1) + { + kdir = vector(1, 0, 0); + } + else if (mesh.directions()[1] == -1) + { + kdir = vector(0, 1, 0); + } + else + { + kdir = vector(0, 0, 1); + } + } + else // 3D so find a direction in the place of the face + { + const face& f = mesh.faces()[faci]; + kdir = mesh.points()[f[0]] - mesh.points()[f[1]]; + } + #else + // Spherical geometry so kdir is the radial direction + kdir = mesh.Cf()[faci]; + #endif + + if (mesh.nGeometricD() == 3) + { + // Remove the idir component from kdir and normalise + kdir -= (idir & kdir)*idir; + + scalar magk = mag(kdir); + + if (magk < SMALL) + { + FatalErrorIn("findFaceDirs") << " calculated kdir = zero" + << exit(FatalError); + } + else + { + kdir /= magk; + } + } + + jdir = kdir ^ idir; +} + + +Foam::label Foam::quadraticFitSnGradData::calcFit +( + const List& C, + const label faci +) +{ + vector idir(1,0,0); + vector jdir(0,1,0); + vector kdir(0,0,1); + findFaceDirs(idir, jdir, kdir, mesh(), faci); + + scalarList wts(C.size(), scalar(1)); + wts[0] = centralWeight_; + wts[1] = centralWeight_; + + point p0 = mesh().faceCentres()[faci]; + scalar scale = 0; + + // calculate the matrix of the polynomial components + scalarRectangularMatrix B(C.size(), minSize_, scalar(0)); + + for(label ip = 0; ip < C.size(); ip++) + { + const point& p = C[ip]; + + scalar px = (p - p0)&idir; + scalar py = (p - p0)&jdir; + #ifdef SPHERICAL_GEOMETRY + scalar pz = mag(p) - mag(p0); + #else + scalar pz = (p - p0)&kdir; + #endif + + if (ip == 0) scale = max(max(mag(px), mag(py)), mag(pz)); + + px /= scale; + py /= scale; + pz /= scale; + + label is = 0; + + B[ip][is++] = wts[0]*wts[ip]; + B[ip][is++] = wts[0]*wts[ip]*px; + B[ip][is++] = wts[ip]*sqr(px); + + if (dim_ >= 2) + { + B[ip][is++] = wts[ip]*py; + B[ip][is++] = wts[ip]*px*py; + B[ip][is++] = wts[ip]*sqr(py); + } + if (dim_ == 3) + { + B[ip][is++] = wts[ip]*pz; + B[ip][is++] = wts[ip]*px*pz; + //B[ip][is++] = wts[ip]*py*pz; + B[ip][is++] = wts[ip]*sqr(pz); + } + } + + // Set the fit + label stencilSize = C.size(); + fit_[faci].setSize(stencilSize); + scalarList singVals(minSize_); + label nSVDzeros = 0; + + const scalar& deltaCoeff = mesh().deltaCoeffs()[faci]; + + bool goodFit = false; + for(int iIt = 0; iIt < 10 && !goodFit; iIt++) + { + SVD svd(B, SMALL); + + scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[1][0]/scale; + scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[1][1]/scale; + + goodFit = + fit0 < 0 && fit1 > 0 + && mag(fit0 + deltaCoeff) < 0.5*deltaCoeff + && mag(fit1 - deltaCoeff) < 0.5*deltaCoeff; + + if (goodFit) + { + fit_[faci][0] = fit0; + fit_[faci][1] = fit1; + for(label i = 2; i < stencilSize; i++) + { + fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[1][i]/scale; + } + singVals = svd.S(); + nSVDzeros = svd.nZeros(); + } + else // (not good fit so increase weight in the centre and for linear) + { + wts[0] *= 10; + wts[1] *= 10; + + for(label i = 0; i < B.n(); i++) + { + B[i][0] *= 10; + B[i][1] *= 10; + } + + for(label j = 0; j < B.m(); j++) + { + B[0][j] *= 10; + B[1][j] *= 10; + } + } + } + + if (goodFit) + { + // remove the uncorrected snGradScheme coefficients + fit_[faci][0] += deltaCoeff; + fit_[faci][1] -= deltaCoeff; + } + else + { + Pout<< "quadratifFitSnGradData could not fit face " << faci + << " fit_[faci][0] = " << fit_[faci][0] + << " fit_[faci][1] = " << fit_[faci][1] + << " deltaCoeff = " << deltaCoeff << endl; + fit_[faci] = 0; + } + + return minSize_ - nSVDzeros; +} + +bool Foam::quadraticFitSnGradData::movePoints() +{ + notImplemented("quadraticFitSnGradData::movePoints()"); + + return true; +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H new file mode 100644 index 0000000000..0f5fda5ce6 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGradData.H @@ -0,0 +1,131 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + quadraticFitSnGradData + +Description + Data for the quadratic fit correction snGrad scheme + +SourceFiles + quadraticFitSnGradData.C + +\*---------------------------------------------------------------------------*/ + +#ifndef quadraticFitSnGradData_H +#define quadraticFitSnGradData_H + +#include "MeshObject.H" +#include "fvMesh.H" +#include "extendedStencil.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class quadraticFitSnGradData Declaration +\*---------------------------------------------------------------------------*/ + +class quadraticFitSnGradData +: + public MeshObject +{ + // Private data + + const scalar centralWeight_; + const label dim_; + + //- minimum stencil size + const label minSize_; + + //- Extended stencil addressing + extendedStencil stencil_; + + //- For each cell in the mesh store the values which multiply the + // values of the stencil to obtain the gradient for each direction + List fit_; + + + // Private member functions + + //- Find the normal direction and i, j and k directions for face faci + static void findFaceDirs + ( + vector& idir, // value changed in return + vector& jdir, // value changed in return + vector& kdir, // value changed in return + const fvMesh& mesh, + const label faci + ); + + label calcFit(const List&, const label faci); + + +public: + + TypeName("quadraticFitSnGradData"); + + + // Constructors + + explicit quadraticFitSnGradData + ( + const fvMesh& mesh, + const scalar cWeight + ); + + + // Destructor + + virtual ~quadraticFitSnGradData() + {}; + + + // Member functions + + //- Return reference to the stencil + const extendedStencil& stencil() const + { + return stencil_; + } + + //- Return reference to fit coefficients + const List& fit() const { return fit_; } + + //- Delete the data when the mesh moves not implemented + virtual bool movePoints(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C new file mode 100644 index 0000000000..aee4af29e3 --- /dev/null +++ b/src/finiteVolume/finiteVolume/snGradSchemes/quadraticFitSnGrad/quadraticFitSnGrads.C @@ -0,0 +1,44 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + Simple central-difference snGrad scheme with quadratic fit correction from + a larger stencil. + +\*---------------------------------------------------------------------------*/ + +#include "quadraticFitSnGrad.H" +#include "fvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fv +{ + makeSnGradScheme(quadraticFitSnGrad) +} +} + +// ************************************************************************* // From 8342283907e145a0c72ec746ca083ab51fde7f25 Mon Sep 17 00:00:00 2001 From: henry Date: Sun, 12 Oct 2008 11:43:08 +0100 Subject: [PATCH 07/11] Added DDES and IDDES SGS models for incompressible LES. --- .../SpalartAllmarasDDES/SpalartAllmarasDDES.C | 120 ++++++++++ .../SpalartAllmarasDDES/SpalartAllmarasDDES.H | 118 ++++++++++ .../IDDESDelta/IDDESDelta.C | 140 ++++++++++++ .../IDDESDelta/IDDESDelta.H | 113 +++++++++ .../SpalartAllmarasIDDES.C | 214 ++++++++++++++++++ .../SpalartAllmarasIDDES.H | 132 +++++++++++ 6 files changed, 837 insertions(+) create mode 100644 src/turbulenceModels/LES/incompressible/SpalartAllmarasDDES/SpalartAllmarasDDES.C create mode 100644 src/turbulenceModels/LES/incompressible/SpalartAllmarasDDES/SpalartAllmarasDDES.H create mode 100644 src/turbulenceModels/LES/incompressible/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C create mode 100644 src/turbulenceModels/LES/incompressible/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.H create mode 100644 src/turbulenceModels/LES/incompressible/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C create mode 100644 src/turbulenceModels/LES/incompressible/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H diff --git a/src/turbulenceModels/LES/incompressible/SpalartAllmarasDDES/SpalartAllmarasDDES.C b/src/turbulenceModels/LES/incompressible/SpalartAllmarasDDES/SpalartAllmarasDDES.C new file mode 100644 index 0000000000..6f5435a268 --- /dev/null +++ b/src/turbulenceModels/LES/incompressible/SpalartAllmarasDDES/SpalartAllmarasDDES.C @@ -0,0 +1,120 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "SpalartAllmarasDDES.H" +#include "addToRunTimeSelectionTable.H" +#include "wallDist.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressible +{ +namespace LESModels +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(SpalartAllmarasDDES, 0); +addToRunTimeSelectionTable(LESModel, SpalartAllmarasDDES, dictionary); + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +tmp SpalartAllmarasDDES::rd +( + const volScalarField& visc, + const volScalarField& S +) const +{ + volScalarField d = wallDist(mesh_).y(); + + tmp trd + ( + new volScalarField + ( + min + ( + visc + /( + max + ( + S, + dimensionedScalar("SMALL", S.dimensions(), SMALL) + )*sqr(kappa_*d) + + dimensionedScalar + ( + "ROOTVSMALL", + dimensionSet(0, 2 , -1, 0, 0), + ROOTVSMALL + ) + ), scalar(10.0) + ) + ) + ); + + return trd; +} + + +tmp SpalartAllmarasDDES::fd(const volScalarField& S) +{ + return 1.0 - tanh(pow3(8.0*rd(nuSgs_ + transport().nu(), S))); +} + + +void SpalartAllmarasDDES::dTildaUpdate(const volScalarField& S) +{ + dTilda_ = + wallDist(mesh_).y() + - fd(S)*max + ( + dimensionedScalar("zero", dimLength, 0.0), + wallDist(mesh_).y() - CDES_*delta() + ); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +SpalartAllmarasDDES::SpalartAllmarasDDES +( + const volVectorField& U, + const surfaceScalarField& phi, + transportModel& transport +) +: + SpalartAllmaras(U, phi, transport, typeName) +{} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace LESModels +} // End namespace incompressible +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/turbulenceModels/LES/incompressible/SpalartAllmarasDDES/SpalartAllmarasDDES.H b/src/turbulenceModels/LES/incompressible/SpalartAllmarasDDES/SpalartAllmarasDDES.H new file mode 100644 index 0000000000..d905c34710 --- /dev/null +++ b/src/turbulenceModels/LES/incompressible/SpalartAllmarasDDES/SpalartAllmarasDDES.H @@ -0,0 +1,118 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::LESmodels::SpalartAllmarasDDES + +Description + SpalartAllmaras DDES LES turbulence model for incompressible flows + + Reference: + P.R. Spalart, S. Deck, S., M.L.Shur, K.D. Squires, M.Kh Strelets, and + A. Travin. `A new version of detached-eddy simulation, resistant to + ambiguous grid densities'. Theor. Comp. Fluid Dyn., 20:181-195, 2006. + +SourceFiles + SpalartAllmarasDDES.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SpalartAllmarasDDES_H +#define SpalartAllmarasDDES_H + +#include "SpalartAllmaras.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressible +{ +namespace LESModels +{ + +/*---------------------------------------------------------------------------*\ + Class SpalartAllmarasDDES Declaration +\*---------------------------------------------------------------------------*/ + +class SpalartAllmarasDDES +: + public SpalartAllmaras +{ + // Private member functions + + // Disallow default bitwise copy construct and assignment + SpalartAllmarasDDES(const SpalartAllmarasDDES&); + SpalartAllmarasDDES& operator=(const SpalartAllmarasDDES&); + + +protected: + + // Protected member functions + + tmp fd(const volScalarField& S); + tmp rd + ( + const volScalarField& visc, + const volScalarField& S + ) const; + + //- Length scale calculation + virtual void dTildaUpdate(const volScalarField& S); + + +public: + + //- Runtime type information + TypeName("SpalartAllmarasDDES"); + + + // Constructors + + //- Constructor from components + SpalartAllmarasDDES + ( + const volVectorField& U, + const surfaceScalarField& phi, + transportModel& transport + ); + + + //- Destructor + virtual ~SpalartAllmarasDDES() + {} +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace LESModels +} // End namespace incompressible +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/turbulenceModels/LES/incompressible/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C b/src/turbulenceModels/LES/incompressible/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C new file mode 100644 index 0000000000..36e6a46a4e --- /dev/null +++ b/src/turbulenceModels/LES/incompressible/SpalartAllmarasIDDES/IDDESDelta/IDDESDelta.C @@ -0,0 +1,140 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "IDDESDelta.H" +#include "addToRunTimeSelectionTable.H" +#include "wallDist.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(IDDESDelta, 0); + addToRunTimeSelectionTable(LESdelta, IDDESDelta, dictionary); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::IDDESDelta::calcDelta() +{ + const Vector