upwind versions of stencils

This commit is contained in:
mattijs
2009-06-12 16:41:13 +01:00
parent ff7be3fd03
commit 3efaa302a9
32 changed files with 915 additions and 78 deletions

View File

@ -95,6 +95,7 @@ DebugSwitches
Euler 0; Euler 0;
EulerImplicit 0; EulerImplicit 0;
EulerRotation 0; EulerRotation 0;
extendedCellToFaceStencil 0;
FDIC 0; FDIC 0;
FaceCellWave 0; FaceCellWave 0;
GAMG 0; GAMG 0;

View File

@ -273,6 +273,15 @@ Foam::mapDistribute::mapDistribute
} }
Foam::mapDistribute::mapDistribute(const mapDistribute& map)
:
constructSize_(map.constructSize_),
subMap_(map.subMap_),
constructMap_(map.constructMap_),
schedulePtr_()
{}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::mapDistribute::compact(const boolList& elemIsUsed) void Foam::mapDistribute::compact(const boolList& elemIsUsed)
@ -413,4 +422,24 @@ void Foam::mapDistribute::compact(const boolList& elemIsUsed)
} }
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void Foam::mapDistribute::operator=(const mapDistribute& rhs)
{
// Check for assignment to self
if (this == &rhs)
{
FatalErrorIn
(
"Foam::mapDistribute::operator=(const Foam::mapDistribute&)"
) << "Attempted assignment to self"
<< abort(FatalError);
}
constructSize_ = rhs.constructSize_;
subMap_ = rhs.subMap_;
constructMap_ = rhs.constructMap_;
schedulePtr_.clear();
}
// ************************************************************************* // // ************************************************************************* //

View File

@ -81,16 +81,6 @@ class mapDistribute
mutable autoPtr<List<labelPair> > schedulePtr_; mutable autoPtr<List<labelPair> > schedulePtr_;
// Private Member Functions
//- Disallow default bitwise copy construct
mapDistribute(const mapDistribute&);
//- Disallow default bitwise assignment
void operator=(const mapDistribute&);
public: public:
// Constructors // Constructors
@ -120,6 +110,9 @@ public:
const labelList& recvProcs const labelList& recvProcs
); );
//- Construct copy
mapDistribute(const mapDistribute&);
// Member Functions // Member Functions
@ -262,6 +255,11 @@ public:
"mapDistribute::updateMesh(const mapPolyMesh&)" "mapDistribute::updateMesh(const mapPolyMesh&)"
); );
} }
// Member Operators
void operator=(const mapDistribute&);
}; };

View File

@ -65,6 +65,7 @@ $(cellToFace)/MeshObjects/upwindCECCellToFaceStencilObject.C
$(cellToFace)/MeshObjects/upwindCFCCellToFaceStencilObject.C $(cellToFace)/MeshObjects/upwindCFCCellToFaceStencilObject.C
$(cellToFace)/MeshObjects/upwindCPCCellToFaceStencilObject.C $(cellToFace)/MeshObjects/upwindCPCCellToFaceStencilObject.C
$(cellToFace)/MeshObjects/upwindFECCellToFaceStencilObject.C $(cellToFace)/MeshObjects/upwindFECCellToFaceStencilObject.C
$(cellToFace)/MeshObjects/pureUpwindCFCCellToFaceStencilObject.C
faceToCell = $(extendedStencil)/faceToCell faceToCell = $(extendedStencil)/faceToCell
$(faceToCell)/fullStencils/faceToCellStencil.C $(faceToCell)/fullStencils/faceToCellStencil.C
@ -216,6 +217,10 @@ $(schemes)/quadraticFit/quadraticFit.C
$(schemes)/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C $(schemes)/quadraticLinearUpwindFit/quadraticLinearUpwindFit.C
$(schemes)/quadraticUpwindFit/quadraticUpwindFit.C $(schemes)/quadraticUpwindFit/quadraticUpwindFit.C
$(schemes)/cubicUpwindFit/cubicUpwindFit.C $(schemes)/cubicUpwindFit/cubicUpwindFit.C
/*
$(schemes)/quadraticLinearPureUpwindFit/quadraticLinearPureUpwindFit.C
*/
$(schemes)/linearPureUpwindFit/linearPureUpwindFit.C
limitedSchemes = $(surfaceInterpolation)/limitedSchemes limitedSchemes = $(surfaceInterpolation)/limitedSchemes
$(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C $(limitedSchemes)/limitedSurfaceInterpolationScheme/limitedSurfaceInterpolationSchemes.C

View File

@ -26,9 +26,6 @@ License
#include "CECCellToCellStencil.H" #include "CECCellToCellStencil.H"
#include "syncTools.H" #include "syncTools.H"
//#include "meshTools.H"
//#include "OFstream.H"
//#include "Time.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //

View File

@ -67,7 +67,14 @@ public:
: :
MeshObject<fvMesh, centredCECCellToFaceStencilObject>(mesh), MeshObject<fvMesh, centredCECCellToFaceStencilObject>(mesh),
extendedCentredCellToFaceStencil(CECCellToFaceStencil(mesh)) extendedCentredCellToFaceStencil(CECCellToFaceStencil(mesh))
{} {
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated centred stencil " << type()
<< nl << endl;
writeStencilStats(Info, stencil(), map());
}
}
// Destructor // Destructor

View File

@ -67,7 +67,14 @@ public:
: :
MeshObject<fvMesh, centredCFCCellToFaceStencilObject>(mesh), MeshObject<fvMesh, centredCFCCellToFaceStencilObject>(mesh),
extendedCentredCellToFaceStencil(CFCCellToFaceStencil(mesh)) extendedCentredCellToFaceStencil(CFCCellToFaceStencil(mesh))
{} {
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated centred stencil " << type()
<< nl << endl;
writeStencilStats(Info, stencil(), map());
}
}
//- Destructor //- Destructor

View File

@ -67,7 +67,14 @@ public:
: :
MeshObject<fvMesh, centredCPCCellToFaceStencilObject>(mesh), MeshObject<fvMesh, centredCPCCellToFaceStencilObject>(mesh),
extendedCentredCellToFaceStencil(CPCCellToFaceStencil(mesh)) extendedCentredCellToFaceStencil(CPCCellToFaceStencil(mesh))
{} {
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated centred stencil " << type()
<< nl << endl;
writeStencilStats(Info, stencil(), map());
}
}
// Destructor // Destructor

View File

@ -67,7 +67,14 @@ public:
: :
MeshObject<fvMesh, centredFECCellToFaceStencilObject>(mesh), MeshObject<fvMesh, centredFECCellToFaceStencilObject>(mesh),
extendedCentredCellToFaceStencil(FECCellToFaceStencil(mesh)) extendedCentredCellToFaceStencil(FECCellToFaceStencil(mesh))
{} {
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated centred stencil " << type()
<< nl << endl;
writeStencilStats(Info, stencil(), map());
}
}
// Destructor // Destructor

View File

@ -0,0 +1,38 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 "pureUpwindCFCCellToFaceStencilObject.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(pureUpwindCFCCellToFaceStencilObject, 0);
}
// ************************************************************************* //

View File

@ -0,0 +1,95 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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::upwindCFCCellToFaceStencilObject
Description
SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef pureUpwindCFCCellToFaceStencilObject_H
#define pureUpwindCFCCellToFaceStencilObject_H
#include "extendedUpwindCellToFaceStencil.H"
#include "CFCCellToFaceStencil.H"
#include "MeshObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class pureUpwindCFCCellToFaceStencilObject Declaration
\*---------------------------------------------------------------------------*/
class pureUpwindCFCCellToFaceStencilObject
:
public MeshObject<fvMesh, pureUpwindCFCCellToFaceStencilObject>,
public extendedUpwindCellToFaceStencil
{
public:
TypeName("pureUpwindCFCCellToFaceStencil");
// Constructors
//- Construct from uncompacted face stencil
explicit pureUpwindCFCCellToFaceStencilObject
(
const fvMesh& mesh
)
:
MeshObject<fvMesh, pureUpwindCFCCellToFaceStencilObject>(mesh),
extendedUpwindCellToFaceStencil(CFCCellToFaceStencil(mesh))
{
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated pure upwind stencil " << type()
<< nl << endl;
writeStencilStats(Info, ownStencil(), ownMap());
}
}
// Destructor
virtual ~pureUpwindCFCCellToFaceStencilObject()
{}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -63,6 +63,7 @@ public:
explicit upwindCECCellToFaceStencilObject explicit upwindCECCellToFaceStencilObject
( (
const fvMesh& mesh, const fvMesh& mesh,
const bool pureUpwind,
const scalar minOpposedness const scalar minOpposedness
) )
: :
@ -70,9 +71,17 @@ public:
extendedUpwindCellToFaceStencil extendedUpwindCellToFaceStencil
( (
CECCellToFaceStencil(mesh), CECCellToFaceStencil(mesh),
pureUpwind,
minOpposedness minOpposedness
) )
{} {
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated off-centred stencil " << type()
<< nl << endl;
writeStencilStats(Info, ownStencil(), ownMap());
}
}
// Destructor // Destructor

View File

@ -35,7 +35,7 @@ SourceFiles
#define upwindCFCCellToFaceStencilObject_H #define upwindCFCCellToFaceStencilObject_H
#include "extendedUpwindCellToFaceStencil.H" #include "extendedUpwindCellToFaceStencil.H"
#include "FECCellToFaceStencil.H" #include "CFCCellToFaceStencil.H"
#include "MeshObject.H" #include "MeshObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -63,16 +63,25 @@ public:
explicit upwindCFCCellToFaceStencilObject explicit upwindCFCCellToFaceStencilObject
( (
const fvMesh& mesh, const fvMesh& mesh,
const bool pureUpwind,
const scalar minOpposedness const scalar minOpposedness
) )
: :
MeshObject<fvMesh, upwindCFCCellToFaceStencilObject>(mesh), MeshObject<fvMesh, upwindCFCCellToFaceStencilObject>(mesh),
extendedUpwindCellToFaceStencil extendedUpwindCellToFaceStencil
( (
FECCellToFaceStencil(mesh), CFCCellToFaceStencil(mesh),
pureUpwind,
minOpposedness minOpposedness
) )
{} {
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated off-centred stencil " << type()
<< nl << endl;
writeStencilStats(Info, ownStencil(), ownMap());
}
}
// Destructor // Destructor

View File

@ -63,6 +63,7 @@ public:
explicit upwindCPCCellToFaceStencilObject explicit upwindCPCCellToFaceStencilObject
( (
const fvMesh& mesh, const fvMesh& mesh,
const bool pureUpwind,
const scalar minOpposedness const scalar minOpposedness
) )
: :
@ -70,9 +71,17 @@ public:
extendedUpwindCellToFaceStencil extendedUpwindCellToFaceStencil
( (
CPCCellToFaceStencil(mesh), CPCCellToFaceStencil(mesh),
pureUpwind,
minOpposedness minOpposedness
) )
{} {
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated off-centred stencil " << type()
<< nl << endl;
writeStencilStats(Info, ownStencil(), ownMap());
}
}
// Destructor // Destructor

View File

@ -63,6 +63,7 @@ public:
explicit upwindFECCellToFaceStencilObject explicit upwindFECCellToFaceStencilObject
( (
const fvMesh& mesh, const fvMesh& mesh,
const bool pureUpwind,
const scalar minOpposedness const scalar minOpposedness
) )
: :
@ -70,9 +71,17 @@ public:
extendedUpwindCellToFaceStencil extendedUpwindCellToFaceStencil
( (
FECCellToFaceStencil(mesh), FECCellToFaceStencil(mesh),
pureUpwind,
minOpposedness minOpposedness
) )
{} {
if (extendedCellToFaceStencil::debug)
{
Info<< "Generated off-centred stencil " << type()
<< nl << endl;
writeStencilStats(Info, ownStencil(), ownMap());
}
}
// Destructor // Destructor

View File

@ -29,8 +29,70 @@ License
#include "syncTools.H" #include "syncTools.H"
#include "SortableList.H" #include "SortableList.H"
/* * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * */
defineTypeNameAndDebug(Foam::extendedCellToFaceStencil, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::extendedCellToFaceStencil::writeStencilStats
(
Ostream& os,
const labelListList& stencil,
const mapDistribute& map
)
{
label sumSize = 0;
label nSum = 0;
label minSize = labelMax;
label maxSize = labelMin;
forAll(stencil, i)
{
const labelList& sCells = stencil[i];
if (sCells.size() > 0)
{
sumSize += sCells.size();
nSum++;
minSize = min(minSize, sCells.size());
maxSize = max(maxSize, sCells.size());
}
}
reduce(sumSize, sumOp<label>());
reduce(nSum, sumOp<label>());
reduce(minSize, minOp<label>());
reduce(maxSize, maxOp<label>());
os << "Stencil size :" << nl
<< " average : " << scalar(sumSize)/nSum << nl
<< " min : " << minSize << nl
<< " max : " << maxSize << nl
<< endl;
// Sum all sent data
label nSent = 0;
label nLocal = 0;
forAll(map.subMap(), procI)
{
if (procI != Pstream::myProcNo())
{
nSent += map.subMap()[procI].size();
}
else
{
nLocal += map.subMap()[procI].size();
}
}
os << "Local data size : " << returnReduce(nLocal, sumOp<label>()) << nl
<< "Sent data size : " << returnReduce(nSent, sumOp<label>()) << nl
<< endl;
}
Foam::autoPtr<Foam::mapDistribute> Foam::autoPtr<Foam::mapDistribute>
Foam::extendedCellToFaceStencil::calcDistributeMap Foam::extendedCellToFaceStencil::calcDistributeMap
( (
@ -211,8 +273,9 @@ Foam::extendedCellToFaceStencil::calcDistributeMap
} }
} }
// Constuct map for distribution of compact data. // Constuct map for distribution of compact data.
return autoPtr<mapDistribute> autoPtr<mapDistribute> mapPtr
( (
new mapDistribute new mapDistribute
( (
@ -222,6 +285,8 @@ Foam::extendedCellToFaceStencil::calcDistributeMap
true // reuse send/recv maps. true // reuse send/recv maps.
) )
); );
return mapPtr;
} }

View File

@ -88,8 +88,22 @@ private:
void operator=(const extendedCellToFaceStencil&); void operator=(const extendedCellToFaceStencil&);
protected:
//- Write some statistics about stencil
static void writeStencilStats
(
Ostream& os,
const labelListList& stencil,
const mapDistribute& map
);
public: public:
// Declare name of the class and its debug switch
ClassName("extendedCellToFaceStencil");
// Constructors // Constructors
//- Construct from mesh //- Construct from mesh

View File

@ -101,7 +101,10 @@ Foam::extendedCellToFaceStencil::weightedSum
( (
fld.name(), fld.name(),
mesh.time().timeName(), mesh.time().timeName(),
mesh mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
), ),
mesh, mesh,
dimensioned<Type> dimensioned<Type>

View File

@ -375,10 +375,12 @@ void Foam::extendedUpwindCellToFaceStencil::transportStencils
Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil
( (
const cellToFaceStencil& stencil, const cellToFaceStencil& stencil,
const bool pureUpwind,
const scalar minOpposedness const scalar minOpposedness
) )
: :
extendedCellToFaceStencil(stencil.mesh()) extendedCellToFaceStencil(stencil.mesh()),
pureUpwind_(pureUpwind)
{ {
//forAll(stencil, faceI) //forAll(stencil, faceI)
//{ //{
@ -430,6 +432,134 @@ Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil
stencil.globalNumbering(), stencil.globalNumbering(),
neiStencil_ neiStencil_
); );
// stencil now in compact form
if (pureUpwind_)
{
const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh());
List<List<point> > stencilPoints(ownStencil_.size());
// Owner stencil
// ~~~~~~~~~~~~~
collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints);
// Mask off all stencil points on wrong side of face
forAll(stencilPoints, faceI)
{
const point& fc = mesh.faceCentres()[faceI];
const vector& fArea = mesh.faceAreas()[faceI];
const List<point>& points = stencilPoints[faceI];
const labelList& stencil = ownStencil_[faceI];
DynamicList<label> newStencil(stencil.size());
forAll(points, i)
{
if (((points[i]-fc) & fArea) < 0)
{
newStencil.append(stencil[i]);
}
}
if (newStencil.size() != stencil.size())
{
ownStencil_[faceI].transfer(newStencil);
}
}
// Neighbour stencil
// ~~~~~~~~~~~~~~~~~
collectData(neiMapPtr_(), neiStencil_, mesh.C(), stencilPoints);
// Mask off all stencil points on wrong side of face
forAll(stencilPoints, faceI)
{
const point& fc = mesh.faceCentres()[faceI];
const vector& fArea = mesh.faceAreas()[faceI];
const List<point>& points = stencilPoints[faceI];
const labelList& stencil = neiStencil_[faceI];
DynamicList<label> newStencil(stencil.size());
forAll(points, i)
{
if (((points[i]-fc) & fArea) > 0)
{
newStencil.append(stencil[i]);
}
}
if (newStencil.size() != stencil.size())
{
neiStencil_[faceI].transfer(newStencil);
}
}
// Note: could compact schedule as well. for if cells are not needed
// across any boundary anymore. However relatively rare.
}
}
Foam::extendedUpwindCellToFaceStencil::extendedUpwindCellToFaceStencil
(
const cellToFaceStencil& stencil
)
:
extendedCellToFaceStencil(stencil.mesh()),
pureUpwind_(true)
{
// Calculate stencil points with full stencil
ownStencil_ = stencil;
ownMapPtr_ = calcDistributeMap
(
stencil.mesh(),
stencil.globalNumbering(),
ownStencil_
);
const fvMesh& mesh = dynamic_cast<const fvMesh&>(stencil.mesh());
List<List<point> > stencilPoints(ownStencil_.size());
collectData(ownMapPtr_(), ownStencil_, mesh.C(), stencilPoints);
// Split stencil into owner and neighbour
neiStencil_.setSize(ownStencil_.size());
forAll(stencilPoints, faceI)
{
const point& fc = mesh.faceCentres()[faceI];
const vector& fArea = mesh.faceAreas()[faceI];
const List<point>& points = stencilPoints[faceI];
const labelList& stencil = ownStencil_[faceI];
DynamicList<label> newOwnStencil(stencil.size());
DynamicList<label> newNeiStencil(stencil.size());
forAll(points, i)
{
if (((points[i]-fc) & fArea) > 0)
{
newNeiStencil.append(stencil[i]);
}
else
{
newOwnStencil.append(stencil[i]);
}
}
if (newNeiStencil.size() > 0)
{
ownStencil_[faceI].transfer(newOwnStencil);
neiStencil_[faceI].transfer(newNeiStencil);
}
}
// Should compact schedule. Or have both return the same schedule.
neiMapPtr_.reset(new mapDistribute(ownMapPtr_()));
} }

View File

@ -27,7 +27,14 @@ Class
Description Description
Creates upwind stencil by shifting a centred stencil to upwind and downwind Creates upwind stencil by shifting a centred stencil to upwind and downwind
faces. faces and optionally removing all non-(up/down)wind faces ('pureUpwind').
Note: the minOpposedness parameter is to decide which upwind and
downwind faces to combine the stencils from. If myArea is the
local area and upwindArea
the area of the possible upwind candidate it will be included if
(upwindArea & myArea)/magSqr(myArea) > minOpposedness
so this includes both cosine and area. WIP.
SourceFiles SourceFiles
extendedUpwindCellToFaceStencil.C extendedUpwindCellToFaceStencil.C
@ -57,6 +64,9 @@ class extendedUpwindCellToFaceStencil
{ {
// Private data // Private data
//- Does stencil contain upwind points only
const bool pureUpwind_;
//- Swap map for getting neigbouring data //- Swap map for getting neigbouring data
autoPtr<mapDistribute> ownMapPtr_; autoPtr<mapDistribute> ownMapPtr_;
autoPtr<mapDistribute> neiMapPtr_; autoPtr<mapDistribute> neiMapPtr_;
@ -115,16 +125,31 @@ public:
// Constructors // Constructors
//- Construct from mesh and uncompacted face stencil //- Construct from mesh and uncompacted centred face stencil.
// Transports facestencil to create owner and neighbour versions.
// pureUpwind to remove any remaining downwind cells.
extendedUpwindCellToFaceStencil extendedUpwindCellToFaceStencil
( (
const cellToFaceStencil&, const cellToFaceStencil&,
const bool pureUpwind,
const scalar minOpposedness const scalar minOpposedness
); );
//- Construct from mesh and uncompacted centred face stencil. Splits
// stencil into owner and neighbour (so always pure upwind)
extendedUpwindCellToFaceStencil
(
const cellToFaceStencil&
);
// Member Functions // Member Functions
bool pureUpwind() const
{
return pureUpwind_;
}
//- Return reference to the parallel distribution map //- Return reference to the parallel distribution map
const mapDistribute& ownMap() const const mapDistribute& ownMap() const
{ {

View File

@ -54,7 +54,10 @@ Foam::extendedUpwindCellToFaceStencil::weightedSum
( (
fld.name(), fld.name(),
mesh.time().timeName(), mesh.time().timeName(),
mesh mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
), ),
mesh, mesh,
dimensioned<Type> dimensioned<Type>

View File

@ -26,6 +26,8 @@ Class
Foam::extendedFaceToCellStencil Foam::extendedFaceToCellStencil
Description Description
Note: transformations on coupled patches not supported. Problem is the
positions of cells reachable through these patches.
SourceFiles SourceFiles
extendedFaceToCellStencil.C extendedFaceToCellStencil.C

View File

@ -49,7 +49,7 @@ Foam::CentredFitData<Polynomial>::CentredFitData
Polynomial Polynomial
> >
( (
mesh, stencil, linearLimitFactor, centralWeight mesh, stencil, true, linearLimitFactor, centralWeight
), ),
coeffs_(mesh.nFaces()) coeffs_(mesh.nFaces())
{ {

View File

@ -36,12 +36,14 @@ Foam::FitData<Form, ExtendedStencil, Polynomial>::FitData
( (
const fvMesh& mesh, const fvMesh& mesh,
const ExtendedStencil& stencil, const ExtendedStencil& stencil,
const bool linearCorrection,
const scalar linearLimitFactor, const scalar linearLimitFactor,
const scalar centralWeight const scalar centralWeight
) )
: :
MeshObject<fvMesh, Form>(mesh), MeshObject<fvMesh, Form>(mesh),
stencil_(stencil), stencil_(stencil),
linearCorrection_(linearCorrection),
linearLimitFactor_(linearLimitFactor), linearLimitFactor_(linearLimitFactor),
centralWeight_(centralWeight), centralWeight_(centralWeight),
# ifdef SPHERICAL_GEOMETRY # ifdef SPHERICAL_GEOMETRY
@ -143,7 +145,10 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
// Setup the point weights // Setup the point weights
scalarList wts(C.size(), scalar(1)); scalarList wts(C.size(), scalar(1));
wts[0] = centralWeight_; wts[0] = centralWeight_;
if (linearCorrection_)
{
wts[1] = centralWeight_; wts[1] = centralWeight_;
}
// Reference point // Reference point
point p0 = this->mesh().faceCentres()[facei]; point p0 = this->mesh().faceCentres()[facei];
@ -218,10 +223,20 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
} }
} }
if (linearCorrection_)
{
goodFit = goodFit =
(mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin) (mag(coeffsi[0] - wLin) < linearLimitFactor_*wLin)
&& (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin)) && (mag(coeffsi[1] - (1 - wLin)) < linearLimitFactor_*(1 - wLin))
&& maxCoeffi <= 1; && maxCoeffi <= 1;
}
else
{
// Upwind: weight on face is 1.
goodFit =
(mag(coeffsi[0] - 1.0) < linearLimitFactor_*1.0)
&& maxCoeffi <= 1;
}
// if (goodFit && iIt > 0) // if (goodFit && iIt > 0)
// { // {
@ -251,7 +266,10 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
// } // }
wts[0] *= 10; wts[0] *= 10;
if (linearCorrection_)
{
wts[1] *= 10; wts[1] *= 10;
}
for(label j = 0; j < B.m(); j++) for(label j = 0; j < B.m(); j++)
{ {
@ -269,11 +287,19 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
if (goodFit) if (goodFit)
{ {
// Remove the uncorrected linear ocefficients if (linearCorrection_)
{
// Remove the uncorrected linear coefficients
coeffsi[0] -= wLin; coeffsi[0] -= wLin;
coeffsi[1] -= 1 - wLin; coeffsi[1] -= 1 - wLin;
} }
else else
{
// Remove the uncorrected upwind coefficients
coeffsi[0] -= 1.0;
}
}
else
{ {
// if (debug) // if (debug)
// { // {

View File

@ -26,7 +26,11 @@ Class
Foam::FitData Foam::FitData
Description Description
Data for the upwinded and centred polynomial fit interpolation schemes Data for the upwinded and centred polynomial fit interpolation schemes.
The linearCorrection_ determines whether the fit is for a corrected
linear scheme (first two coefficients are corrections for owner and
neighbour) or a pure upwind scheme (first coefficient is correction for
owner ; weight on face taken as 1).
SourceFiles SourceFiles
FitData.C FitData.C
@ -58,7 +62,11 @@ class FitData
//- The stencil the fit is based on //- The stencil the fit is based on
const ExtendedStencil& stencil_; const ExtendedStencil& stencil_;
//- Factor the fit is allowed to deviate from linear. //- Is scheme correction on linear (true) or on upwind (false)
const bool linearCorrection_;
//- Factor the fit is allowed to deviate from the base scheme
// (linear or pure upwind)
// This limits the amount of high-order correction and increases // This limits the amount of high-order correction and increases
// stability on bad meshes // stability on bad meshes
const scalar linearLimitFactor_; const scalar linearLimitFactor_;
@ -84,11 +92,6 @@ class FitData
const label faci const label faci
); );
//- Calculate the fit for the all the mesh faces
// and set the coefficients
// virtual void calcFit();
public: public:
//TypeName("FitData"); //TypeName("FitData");
@ -101,6 +104,7 @@ public:
( (
const fvMesh& mesh, const fvMesh& mesh,
const ExtendedStencil& stencil, const ExtendedStencil& stencil,
const bool linearCorrection,
const scalar linearLimitFactor, const scalar linearLimitFactor,
const scalar centralWeight const scalar centralWeight
); );
@ -119,12 +123,18 @@ public:
return stencil_; return stencil_;
} }
bool linearCorrection() const
{
return linearCorrection_;
}
//- Calculate the fit for the specified face and set the coefficients //- Calculate the fit for the specified face and set the coefficients
void calcFit void calcFit
( (
scalarList& coeffsi, // coefficients to be set scalarList& coeffsi, // coefficients to be set
const List<point>&, // Stencil points const List<point>&, // Stencil points
const scalar wLin, // Linear weight const scalar wLin, // Weight for linear approximation (weights
// nearest neighbours)
const label faci // Current face index const label faci // Current face index
); );
@ -132,7 +142,7 @@ public:
virtual void calcFit() = 0; virtual void calcFit() = 0;
//- Delete the data when the mesh moves not implemented //- Recalculate weights (but not stencil) when the mesh moves
bool movePoints(); bool movePoints();
}; };

View File

@ -0,0 +1,186 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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::PureUpwindFitScheme
Description
Upwind biased fit surface interpolation scheme that applies an explicit
correction to upwind.
\*---------------------------------------------------------------------------*/
#ifndef PureUpwindFitScheme_H
#define PureUpwindFitScheme_H
#include "UpwindFitData.H"
#include "upwind.H"
#include "Switch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class PureUpwindFitScheme Declaration
\*---------------------------------------------------------------------------*/
template<class Type, class Polynomial, class Stencil>
class PureUpwindFitScheme
:
public upwind<Type>
{
// Private Data
//- Factor the fit is allowed to deviate from linear.
// This limits the amount of high-order correction and increases
// stability on bad meshes
const scalar linearLimitFactor_;
//- Weights for central stencil
const scalar centralWeight_;
// Private Member Functions
//- Disallow default bitwise copy construct
PureUpwindFitScheme(const PureUpwindFitScheme&);
//- Disallow default bitwise assignment
void operator=(const PureUpwindFitScheme&);
public:
//- Runtime type information
TypeName("PureUpwindFitScheme");
// Constructors
//- Construct from mesh and Istream
// The name of the flux field is read from the Istream and looked-up
// from the mesh objectRegistry
PureUpwindFitScheme(const fvMesh& mesh, Istream& is)
:
upwind<Type>
(
mesh,
mesh.lookupObject<surfaceScalarField>(word(is))
),
linearLimitFactor_(readScalar(is)),
centralWeight_(1000)
{}
//- Construct from mesh, faceFlux and Istream
PureUpwindFitScheme
(
const fvMesh& mesh,
const surfaceScalarField& faceFlux,
Istream& is
)
:
upwind<Type>(mesh, faceFlux),
linearLimitFactor_(readScalar(is)),
centralWeight_(1000)
{}
// Member Functions
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return true;
}
//- Return the explicit correction to the face-interpolate
virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
correction
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
const fvMesh& mesh = this->mesh();
// Use the owner/neighbour splitting constructor
const extendedUpwindCellToFaceStencil& stencil = Stencil::New(mesh);
const UpwindFitData<Polynomial>& ufd =
UpwindFitData<Polynomial>::New
(
mesh,
stencil,
false, //offset to upwind
linearLimitFactor_,
centralWeight_
);
const List<scalarList>& fo = ufd.owncoeffs();
const List<scalarList>& fn = ufd.neicoeffs();
return stencil.weightedSum(this->faceFlux_, vf, fo, fn);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Add the patch constructor functions to the hash tables
#define makePureUpwindFitSurfaceInterpolationTypeScheme(SS, POLYNOMIAL, STENCIL, TYPE) \
\
typedef PureUpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> \
PureUpwindFitScheme##TYPE##POLYNOMIAL##STENCIL##_; \
defineTemplateTypeNameAndDebugWithName \
(PureUpwindFitScheme##TYPE##POLYNOMIAL##STENCIL##_, #SS, 0); \
\
surfaceInterpolationScheme<TYPE>::addMeshConstructorToTable \
<PureUpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> > \
add##SS##STENCIL##TYPE##MeshConstructorToTable_; \
\
surfaceInterpolationScheme<TYPE>::addMeshFluxConstructorToTable \
<PureUpwindFitScheme<TYPE, POLYNOMIAL, STENCIL> > \
add##SS##STENCIL##TYPE##MeshFluxConstructorToTable_;
#define makePureUpwindFitSurfaceInterpolationScheme(SS, POLYNOMIAL, STENCIL) \
\
makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,scalar) \
makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,vector) \
makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,sphericalTensor) \
makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,symmTensor) \
makePureUpwindFitSurfaceInterpolationTypeScheme(SS,POLYNOMIAL,STENCIL,tensor)
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -38,6 +38,7 @@ Foam::UpwindFitData<Polynomial>::UpwindFitData
( (
const fvMesh& mesh, const fvMesh& mesh,
const extendedUpwindCellToFaceStencil& stencil, const extendedUpwindCellToFaceStencil& stencil,
const bool linearCorrection,
const scalar linearLimitFactor, const scalar linearLimitFactor,
const scalar centralWeight const scalar centralWeight
) )
@ -49,7 +50,7 @@ Foam::UpwindFitData<Polynomial>::UpwindFitData
Polynomial Polynomial
> >
( (
mesh, stencil, linearLimitFactor, centralWeight mesh, stencil, linearCorrection, linearLimitFactor, centralWeight
), ),
owncoeffs_(mesh.nFaces()), owncoeffs_(mesh.nFaces()),
neicoeffs_(mesh.nFaces()) neicoeffs_(mesh.nFaces())
@ -77,30 +78,25 @@ void Foam::UpwindFitData<Polynomial>::calcFit()
{ {
const fvMesh& mesh = this->mesh(); const fvMesh& mesh = this->mesh();
const surfaceScalarField& w = mesh.surfaceInterpolation::weights();
const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField();
// Owner stencil weights
// ~~~~~~~~~~~~~~~~~~~~~
// Get the cell/face centres in stencil order. // Get the cell/face centres in stencil order.
// Upwind face stencils no good for triangles or tets. List<List<point> > stencilPoints(mesh.nFaces());
// Need bigger stencils
List<List<point> > ownStencilPoints(mesh.nFaces());
this->stencil().collectData this->stencil().collectData
( (
this->stencil().ownMap(), this->stencil().ownMap(),
this->stencil().ownStencil(), this->stencil().ownStencil(),
mesh.C(), mesh.C(),
ownStencilPoints stencilPoints
);
List<List<point> > neiStencilPoints(mesh.nFaces());
this->stencil().collectData
(
this->stencil().neiMap(),
this->stencil().neiStencil(),
mesh.C(),
neiStencilPoints
); );
// find the fit coefficients for every owner and neighbour of ever face // find the fit coefficients for every owner
const surfaceScalarField& w = mesh.surfaceInterpolation::weights();
Pout<< "-- Owner --" << endl;
for(label facei = 0; facei < mesh.nInternalFaces(); facei++) for(label facei = 0; facei < mesh.nInternalFaces(); facei++)
{ {
FitData FitData
@ -108,16 +104,17 @@ void Foam::UpwindFitData<Polynomial>::calcFit()
UpwindFitData<Polynomial>, UpwindFitData<Polynomial>,
extendedUpwindCellToFaceStencil, extendedUpwindCellToFaceStencil,
Polynomial Polynomial
>::calcFit(owncoeffs_[facei], ownStencilPoints[facei], w[facei], facei); >::calcFit(owncoeffs_[facei], stencilPoints[facei], w[facei], facei);
FitData
<
UpwindFitData<Polynomial>,
extendedUpwindCellToFaceStencil,
Polynomial
>::calcFit(neicoeffs_[facei], neiStencilPoints[facei], w[facei], facei);
}
const surfaceScalarField::GeometricBoundaryField& bw = w.boundaryField(); Pout<< " facei:" << facei
<< " at:" << mesh.faceCentres()[facei] << endl;
forAll(owncoeffs_[facei], i)
{
Pout<< " point:" << stencilPoints[facei][i]
<< "\tweight:" << owncoeffs_[facei][i]
<< endl;
}
}
forAll(bw, patchi) forAll(bw, patchi)
{ {
@ -136,8 +133,58 @@ void Foam::UpwindFitData<Polynomial>::calcFit()
Polynomial Polynomial
>::calcFit >::calcFit
( (
owncoeffs_[facei], ownStencilPoints[facei], pw[i], facei owncoeffs_[facei], stencilPoints[facei], pw[i], facei
); );
facei++;
}
}
}
// Neighbour stencil weights
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Note:reuse stencilPoints since is major storage
this->stencil().collectData
(
this->stencil().neiMap(),
this->stencil().neiStencil(),
mesh.C(),
stencilPoints
);
// find the fit coefficients for every neighbour
Pout<< "-- Neighbour --" << endl;
for(label facei = 0; facei < mesh.nInternalFaces(); facei++)
{
FitData
<
UpwindFitData<Polynomial>,
extendedUpwindCellToFaceStencil,
Polynomial
>::calcFit(neicoeffs_[facei], stencilPoints[facei], w[facei], facei);
Pout<< " facei:" << facei
<< " at:" << mesh.faceCentres()[facei] << endl;
forAll(neicoeffs_[facei], i)
{
Pout<< " point:" << stencilPoints[facei][i]
<< "\tweight:" << neicoeffs_[facei][i]
<< endl;
}
}
forAll(bw, patchi)
{
const fvsPatchScalarField& pw = bw[patchi];
if (pw.coupled())
{
label facei = pw.patch().patch().start();
forAll(pw, i)
{
FitData FitData
< <
UpwindFitData<Polynomial>, UpwindFitData<Polynomial>,
@ -145,7 +192,7 @@ void Foam::UpwindFitData<Polynomial>::calcFit()
Polynomial Polynomial
>::calcFit >::calcFit
( (
neicoeffs_[facei], neiStencilPoints[facei], pw[i], facei neicoeffs_[facei], stencilPoints[facei], pw[i], facei
); );
facei++; facei++;
} }

View File

@ -27,7 +27,9 @@ Class
Description Description
Data for the quadratic fit correction interpolation scheme to be used with Data for the quadratic fit correction interpolation scheme to be used with
upwind biased stencil upwind biased stencil.
- linearCorrection = true : fit calculated for corrected linear scheme
- linearCorrection = false : fit calculated for corrected upwind scheme
SourceFiles SourceFiles
UpwindFitData.C UpwindFitData.C
@ -90,6 +92,7 @@ public:
( (
const fvMesh& mesh, const fvMesh& mesh,
const extendedUpwindCellToFaceStencil& stencil, const extendedUpwindCellToFaceStencil& stencil,
const bool linearCorrection,
const scalar linearLimitFactor, const scalar linearLimitFactor,
const scalar centralWeight const scalar centralWeight
); );

View File

@ -129,6 +129,7 @@ public:
const extendedUpwindCellToFaceStencil& stencil = Stencil::New const extendedUpwindCellToFaceStencil& stencil = Stencil::New
( (
mesh, mesh,
false, //pureUpwind
scalar(0.5) scalar(0.5)
); );
@ -137,6 +138,7 @@ public:
( (
mesh, mesh,
stencil, stencil,
true, //calculate as offset to linear
linearLimitFactor_, linearLimitFactor_,
centralWeight_ centralWeight_
); );

View File

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 "PureUpwindFitScheme.H"
#include "linearFitPolynomial.H"
#include "pureUpwindCFCCellToFaceStencilObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTemplateTypeNameAndDebug
(
UpwindFitData<linearFitPolynomial>,
0
);
makePureUpwindFitSurfaceInterpolationScheme
(
linearPureUpwindFit,
linearFitPolynomial,
pureUpwindCFCCellToFaceStencilObject
);
}
// ************************************************************************* //

View File

@ -0,0 +1,45 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 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 "PureUpwindFitScheme.H"
#include "quadraticLinearUpwindFitPolynomial.H"
#include "upwindCFCCellToFaceStencilObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Use stencil with three upwind cells:
// upwindCFCCellToFaceStencilObject + pureUpwind
makePureUpwindFitSurfaceInterpolationScheme
(
quadraticLinearPureUpwindFit,
quadraticLinearUpwindFitPolynomial,
upwindCFCCellToFaceStencilObject
);
}
// ************************************************************************* //

View File

@ -26,7 +26,7 @@ License
#include "UpwindFitScheme.H" #include "UpwindFitScheme.H"
#include "quadraticLinearUpwindFitPolynomial.H" #include "quadraticLinearUpwindFitPolynomial.H"
#include "upwindCFCCellToFaceStencilObject.H" #include "upwindFECCellToFaceStencilObject.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -42,7 +42,7 @@ namespace Foam
( (
quadraticLinearUpwindFit, quadraticLinearUpwindFit,
quadraticLinearUpwindFitPolynomial, quadraticLinearUpwindFitPolynomial,
upwindCFCCellToFaceStencilObject upwindFECCellToFaceStencilObject
); );
} }