zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of github.com-OpenFOAM:OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
Henry Weller
2021-12-08 13:39:12 +00:00
parent 7dfb7146ea 053eed714d
commit bf042b39a7
19 changed files with 829 additions and 935 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
applications/utilities/postProcessing/miscellaneous/postChannel/Make/files
View File

@ -1,4 +0,0 @@
postChannel.C
channelIndex.C
EXE = $(FOAM_APPBIN)/postChannel

10
applications/utilities/postProcessing/miscellaneous/postChannel/Make/options
View File

@ -1,10 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lmeshTools \
-lfiniteVolume \
-lgenericPatchFields \
-lsampling

73
applications/utilities/postProcessing/miscellaneous/postChannel/calculateFields.H
View File

@ -1,73 +0,0 @@
/*
volTensorField gradU(fvc::grad(U));
volSymmTensorField D(symm(fvc::grad(U)));
volTensorField Dprim(symm(fvc::grad(U - UMean)));
volScalarField prod(-((U - UMean)*(U - UMean)) && D);
*/
/*
volScalarField txx
(
IOobject
(
"txx",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionSet(0, 1, -1, 0, 0)
);
txx =sqrt(Txx - (UMeanx*UMeanx));
txx.write();
volScalarField tyy
(
IOobject
(
"tyy",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionSet(0, 1, -1, 0, 0)
);
tyy = sqrt(Tyy - (UMeany*UMeany));
tyy.write();
volScalarField tzz
(
IOobject
(
"tzz",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionSet(0, 1, -1, 0, 0)
);
tzz = sqrt(Tzz - (UMeanz*UMeanz));
tzz.write();
volScalarField txy
(
IOobject
(
"txy",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionSet(0, 2, -2, 0, 0)
);
txy = Txy - (UMeanx*UMeany);
txy.write();
*/

293
applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
View File

@ -1,293 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "channelIndex.H"
#include "boolList.H"
#include "syncTools.H"
#include "OFstream.H"
#include "meshTools.H"
#include "Time.H"
#include "SortableList.H"
// * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * * //
namespace Foam
{
template<>
const char* Foam::NamedEnum
<
Foam::vector::components,
3
>::names[] =
{
"x",
"y",
"z"
};
}
const Foam::NamedEnum<Foam::vector::components, 3>
Foam::channelIndex::vectorComponentsNames_;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Determines face blocking
void Foam::channelIndex::walkOppositeFaces
(
const polyMesh& mesh,
const labelList& startFaces,
boolList& blockedFace
)
{
const cellList& cells = mesh.cells();
const faceList& faces = mesh.faces();
label nBnd = mesh.nFaces() - mesh.nInternalFaces();
DynamicList<label> frontFaces(startFaces);
forAll(frontFaces, i)
{
label facei = frontFaces[i];
blockedFace[facei] = true;
}
while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
{
// Transfer across.
boolList isFrontBndFace(nBnd, false);
forAll(frontFaces, i)
{
label facei = frontFaces[i];
if (!mesh.isInternalFace(facei))
{
isFrontBndFace[facei-mesh.nInternalFaces()] = true;
}
}
syncTools::swapBoundaryFaceList(mesh, isFrontBndFace);
// Add
forAll(isFrontBndFace, i)
{
label facei = mesh.nInternalFaces()+i;
if (isFrontBndFace[i] && !blockedFace[facei])
{
blockedFace[facei] = true;
frontFaces.append(facei);
}
}
// Transfer across cells
DynamicList<label> newFrontFaces(frontFaces.size());
forAll(frontFaces, i)
{
label facei = frontFaces[i];
{
const cell& ownCell = cells[mesh.faceOwner()[facei]];
label oppositeFacei = ownCell.opposingFaceLabel(facei, faces);
if (oppositeFacei == -1)
{
FatalErrorInFunction
<< "Face:" << facei << " owner cell:" << ownCell
<< " is not a hex?" << abort(FatalError);
}
else
{
if (!blockedFace[oppositeFacei])
{
blockedFace[oppositeFacei] = true;
newFrontFaces.append(oppositeFacei);
}
}
}
if (mesh.isInternalFace(facei))
{
const cell& neiCell = mesh.cells()[mesh.faceNeighbour()[facei]];
label oppositeFacei = neiCell.opposingFaceLabel(facei, faces);
if (oppositeFacei == -1)
{
FatalErrorInFunction
<< "Face:" << facei << " neighbour cell:" << neiCell
<< " is not a hex?" << abort(FatalError);
}
else
{
if (!blockedFace[oppositeFacei])
{
blockedFace[oppositeFacei] = true;
newFrontFaces.append(oppositeFacei);
}
}
}
}
frontFaces.transfer(newFrontFaces);
}
}
// Calculate regions.
void Foam::channelIndex::calcLayeredRegions
(
const polyMesh& mesh,
const labelList& startFaces
)
{
boolList blockedFace(mesh.nFaces(), false);
walkOppositeFaces
(
mesh,
startFaces,
blockedFace
);
if (false)
{
OFstream str(mesh.time().path()/"blockedFaces.obj");
label vertI = 0;
forAll(blockedFace, facei)
{
if (blockedFace[facei])
{
const face& f = mesh.faces()[facei];
forAll(f, fp)
{
meshTools::writeOBJ(str, mesh.points()[f[fp]]);
}
str<< 'f';
forAll(f, fp)
{
str << ' ' << vertI+fp+1;
}
str << nl;
vertI += f.size();
}
}
}
// Do analysis for connected regions
cellRegion_.reset(new regionSplit(mesh, blockedFace));
Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
// Sum number of entries per region
regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
// Average cell centres to determine ordering.
pointField regionCc
(
regionSum(mesh.cellCentres())
/ regionCount_
);
SortableList<scalar> sortComponent(regionCc.component(dir_));
sortMap_ = sortComponent.indices();
y_ = sortComponent;
if (symmetric_)
{
y_.setSize(cellRegion_().nRegions()/2);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::channelIndex::channelIndex
(
const polyMesh& mesh,
const dictionary& dict
)
:
symmetric_(readBool(dict.lookup("symmetric"))),
dir_(vectorComponentsNames_.read(dict.lookup("component")))
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
const wordList patchNames(dict.lookup("patches"));
label nFaces = 0;
forAll(patchNames, i)
{
const label patchi = patches.findPatchID(patchNames[i]);
if (patchi == -1)
{
FatalErrorInFunction
<< "Illegal patch " << patchNames[i]
<< ". Valid patches are " << patches.name()
<< exit(FatalError);
}
nFaces += patches[patchi].size();
}
labelList startFaces(nFaces);
nFaces = 0;
forAll(patchNames, i)
{
const polyPatch& pp = patches[patchNames[i]];
forAll(pp, j)
{
startFaces[nFaces++] = pp.start()+j;
}
}
// Calculate regions.
calcLayeredRegions(mesh, startFaces);
}
Foam::channelIndex::channelIndex
(
const polyMesh& mesh,
const labelList& startFaces,
const bool symmetric,
const direction dir
)
:
symmetric_(symmetric),
dir_(dir)
{
// Calculate regions.
calcLayeredRegions(mesh, startFaces);
}
// ************************************************************************* //

160
applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.H
View File

@ -1,160 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::channelIndex
Description
Does averaging of fields over layers of cells. Assumes layered mesh.
SourceFiles
channelIndex.C
\*---------------------------------------------------------------------------*/
#ifndef channelIndex_H
#define channelIndex_H
#include "regionSplit.H"
#include "direction.H"
#include "scalarField.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class channelIndex Declaration
\*---------------------------------------------------------------------------*/
class channelIndex
{
// Private Data
static const NamedEnum<vector::components, 3> vectorComponentsNames_;
//- Is mesh symmetric
const bool symmetric_;
//- Direction to sort
const direction dir_;
//- Per cell the global region
autoPtr<regionSplit> cellRegion_;
//- Per global region the number of cells (scalarField so we can use
// field algebra)
scalarField regionCount_;
//- From sorted region back to unsorted global region
labelList sortMap_;
//- Sorted component of cell centres
scalarField y_;
// Private Member Functions
void walkOppositeFaces
(
const polyMesh& mesh,
const labelList& startFaces,
boolList& blockedFace
);
void calcLayeredRegions
(
const polyMesh& mesh,
const labelList& startFaces
);
public:
// Constructors
//- Construct from dictionary
channelIndex(const polyMesh&, const dictionary&);
//- Construct from supplied starting faces
channelIndex
(
const polyMesh& mesh,
const labelList& startFaces,
const bool symmetric,
const direction dir
);
//- Disallow default bitwise copy construction
channelIndex(const channelIndex&) = delete;
// Member Functions
// Access
//- Sum field per region
template<class T>
Field<T> regionSum(const Field<T>& cellField) const;
//- Collapse a field to a line
template<class T>
Field<T> collapse
(
const Field<T>& vsf,
const bool asymmetric=false
) const;
//- Return the field of Y locations from the cell centres
const scalarField& y() const
{
return y_;
}
// Member Operators
//- Disallow default bitwise assignment
void operator=(const channelIndex&);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "channelIndexTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

102
applications/utilities/postProcessing/miscellaneous/postChannel/channelIndexTemplates.C
View File

@ -1,102 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "channelIndex.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class T>
Foam::Field<T> Foam::channelIndex::regionSum(const Field<T>& cellField) const
{
Field<T> regionField(cellRegion_().nRegions(), Zero);
forAll(cellRegion_(), celli)
{
regionField[cellRegion_()[celli]] += cellField[celli];
}
// Global sum
Pstream::listCombineGather(regionField, plusEqOp<T>());
Pstream::listCombineScatter(regionField);
return regionField;
}
template<class T>
Foam::Field<T> Foam::channelIndex::collapse
(
const Field<T>& cellField,
const bool asymmetric
) const
{
// Average and order
const Field<T> summedField(regionSum(cellField));
Field<T> regionField
(
summedField
/ regionCount_,
sortMap_
);
// Symmetry?
if (symmetric_)
{
label nlb2 = cellRegion_().nRegions()/2;
if (asymmetric)
{
for (label j=0; j<nlb2; j++)
{
regionField[j] =
0.5
* (
regionField[j]
- regionField[cellRegion_().nRegions() - j - 1]
);
}
}
else
{
for (label j=0; j<nlb2; j++)
{
regionField[j] =
0.5
* (
regionField[j]
+ regionField[cellRegion_().nRegions() - j - 1]
);
}
}
regionField.setSize(nlb2);
}
return regionField;
}
// ************************************************************************* //

69
applications/utilities/postProcessing/miscellaneous/postChannel/collapse.H
View File

@ -1,69 +0,0 @@
scalarField UMeanXvalues
(
channelIndexing.collapse(UMean.component(vector::X)())
);
scalarField UMeanYvalues
(
channelIndexing.collapse(UMean.component(vector::Y)())
);
scalarField UMeanZvalues
(
channelIndexing.collapse(UMean.component(vector::Z)())
);
scalarField RxxValues(channelIndexing.collapse(Rxx));
scalarField RyyValues(channelIndexing.collapse(Ryy));
scalarField RzzValues(channelIndexing.collapse(Rzz));
scalarField RxyValues(channelIndexing.collapse(Rxy, true));
scalarField pPrime2MeanValues(channelIndexing.collapse(pPrime2Mean));
/*
scalarField epsilonValues(channelIndexing.collapse(epsilonMean));
scalarField nuMeanValues(channelIndexing.collapse(nuMean));
scalarField nuPrimeValues(channelIndexing.collapse(nuPrime));
scalarField gammaDotMeanValues(channelIndexing.collapse(gammaDotMean));
scalarField gammaDotPrimeValues(channelIndexing.collapse(gammaDotPrime));
*/
scalarField urmsValues(sqrt(mag(RxxValues)));
scalarField vrmsValues(sqrt(mag(RyyValues)));
scalarField wrmsValues(sqrt(mag(RzzValues)));
scalarField kValues
(
0.5*(sqr(urmsValues) + sqr(vrmsValues) + sqr(wrmsValues))
);
setWriter::New(runTime.graphFormat())->write
(
runTime.globalPath()
/functionObjects::writeFile::outputPrefix
/args.executable()
/runTime.timeName(),
args.executable(),
coordSet(true, "y", channelIndexing.y()),
"Uf", UMeanXvalues,
"u", urmsValues,
"v", vrmsValues,
"w", wrmsValues,
"uv", RxyValues,
"k", kValues,
"pPrime2Mean", pPrime2MeanValues
// "epsilon", epsilonValues,
// "nu", nuMeanValues,
// "nuPrime", nuPrimeValues,
// "gammaDot", gammaDotMeanValues,
// "gammaDotPrime", gammaDotPrimeValues
);

97
applications/utilities/postProcessing/miscellaneous/postChannel/postChannel.C
View File

@ -1,97 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
postChannel
Description
Post-processes data from channel flow calculations.
For each time: calculate: txx, txy,tyy, txy,
eps, prod, vorticity, enstrophy and helicity. Assuming that the mesh
is periodic in the x and z directions, collapse Umeanx, Umeany, txx,
txy and tyy to a line and print them as standard output.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "timeSelector.H"
#include "volFields.H"
#include "channelIndex.H"
#include "setWriter.H"
#include "writeFile.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::noParallel();
timeSelector::addOptions();
#include "setRootCase.H"
#include "createTime.H"
// Get times list
instantList timeDirs = timeSelector::select0(runTime, args);
#include "createMeshNoChangers.H"
#include "readPhysicalProperties.H"
// Setup channel indexing for averaging over channel down to a line
IOdictionary channelDict
(
IOobject
(
"postChannelDict",
mesh.time().system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
channelIndex channelIndexing(mesh, channelDict);
// For each time step read all fields
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
Info<< "Collapsing fields for time " << runTime.timeName() << endl;
#include "readFields.H"
#include "calculateFields.H"
// Average fields over channel down to a line
#include "collapse.H"
}
Info<< "\nEnd\n" << endl;
return 0;
}
// ************************************************************************* //

112
applications/utilities/postProcessing/miscellaneous/postChannel/readFields.H
View File

@ -1,112 +0,0 @@
typeIOobject<volVectorField> UMeanHeader
(
"UMean",
runTime.timeName(),
mesh,
IOobject::MUST_READ
);
if (!UMeanHeader.headerOk())
{
Info<< " No UMean field" << endl;
continue;
}
volVectorField UMean
(
UMeanHeader,
mesh
);
volSymmTensorField UPrime2Mean
(
IOobject
(
"UPrime2Mean",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField Rxx(UPrime2Mean.component(symmTensor::XX));
volScalarField Ryy(UPrime2Mean.component(symmTensor::YY));
volScalarField Rzz(UPrime2Mean.component(symmTensor::ZZ));
volScalarField Rxy(UPrime2Mean.component(symmTensor::XY));
volScalarField pPrime2Mean
(
IOobject
(
"pPrime2Mean",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
/*
volScalarField epsilonMean
(
IOobject
(
"epsilonMean",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField nuMean
(
IOobject
(
"nuMean",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField gammaDotMean
(
IOobject
(
"gammaDotMean",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField nuPrime2
(
IOobject
(
"nuPrime",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField nuPrime = sqrt(mag(nuPrime2 - sqr(nuMean)));
volScalarField gammaDotPrime2
(
IOobject
(
"gammaDotPrime",
runTime.timeName(),
mesh,
IOobject::MUST_READ
),
mesh
);
volScalarField gammaDotPrime = sqrt(mag(gammaDotPrime2 -sqr(gammaDotMean)));
*/

13
applications/utilities/postProcessing/miscellaneous/postChannel/readPhysicalProperties.H
View File

@ -1,13 +0,0 @@
Info<< nl << "Reading physicalProperties" << endl;
IOdictionary transportProperties
(
IOobject
(
"physicalProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);

26
etc/caseDicts/postProcessing/graphs/graphLayerAverage Normal file
View File

@ -0,0 +1,26 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: dev
\\/ M anipulation |
-------------------------------------------------------------------------------
Description
Generates plots of fields averaged over the layers in the mesh
\*---------------------------------------------------------------------------*/
patches (<patchNames>); // Patches from which layers extrude
zones (); // Zones from which layers extrude
axis <axisName>; // The independent variable of the graph. Can
// be "x", "y" or "z".
symmetric false; // Are the layers symmetric about the centre?
fields (<fieldsNames>); // Fields to plot
#includeEtc "caseDicts/postProcessing/graphs/graphLayerAverage.cfg"
// ************************************************************************* //

17
etc/caseDicts/postProcessing/graphs/graphLayerAverage.cfg Normal file
View File

@ -0,0 +1,17 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: dev
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
type layerAverage;
libs ("libfieldFunctionObjects.so");
writeControl writeTime;
writeInterval 1;
setFormat raw;
// ************************************************************************* //

2
src/functionObjects/field/Make/files
View File

@ -81,4 +81,6 @@ cylindrical/cylindricalFunctionObject.C
uniform/uniform.C
layerAverage/layerAverage.C
LIB = $(FOAM_LIBBIN)/libfieldFunctionObjects

420
src/functionObjects/field/layerAverage/layerAverage.C Normal file
View File

@ -0,0 +1,420 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "layerAverage.H"
#include "regionSplit.H"
#include "syncTools.H"
#include "volFields.H"
#include "writeFile.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace functionObjects
{
defineTypeNameAndDebug(layerAverage, 0);
addToRunTimeSelectionTable(functionObject, layerAverage, dictionary);
}
template<>
const char*
Foam::NamedEnum<Foam::vector::components, 3>::names[] =
{"x", "y", "z"};
}
const Foam::NamedEnum<Foam::vector::components, 3>
Foam::functionObjects::layerAverage::axisNames_;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::functionObjects::layerAverage::walkOppositeFaces
(
const labelList& startFaces,
const boolList& startFaceIntoOwners,
boolList& blockedFace,
boolList& cellIsLayer
)
{
// Initialise the front
DynamicList<label> frontFaces(startFaces);
DynamicList<bool> frontFaceIntoOwners(startFaceIntoOwners);
DynamicList<label> newFrontFaces(frontFaces.size());
DynamicList<bool> newFrontFaceIntoOwners(frontFaceIntoOwners.size());
// Set the start and end faces as blocked
UIndirectList<bool>(blockedFace, startFaces) = true;
// Iterate until the front is empty
while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
{
// Transfer front faces across couplings
boolList bndFaceIsFront(mesh_.nFaces() - mesh_.nInternalFaces(), false);
forAll(frontFaces, i)
{
const label facei = frontFaces[i];
const label intoOwner = frontFaceIntoOwners[i];
if (!mesh_.isInternalFace(facei) && !intoOwner)
{
bndFaceIsFront[facei - mesh_.nInternalFaces()] = true;
}
}
syncTools::swapBoundaryFaceList(mesh_, bndFaceIsFront);
// Set new front faces
forAll(bndFaceIsFront, i)
{
const label facei = mesh_.nInternalFaces() + i;
if (bndFaceIsFront[i] && !blockedFace[facei])
{
frontFaces.append(facei);
frontFaceIntoOwners.append(true);
}
}
// Transfer across cells
newFrontFaces.clear();
newFrontFaceIntoOwners.clear();
forAll(frontFaces, i)
{
const label facei = frontFaces[i];
const label intoOwner = frontFaceIntoOwners[i];
const label celli =
intoOwner
? mesh_.faceOwner()[facei]
: mesh_.isInternalFace(facei) ? mesh_.faceNeighbour()[facei] : -1;
if (celli != -1)
{
cellIsLayer[celli] = true;
const label oppositeFacei =
mesh_.cells()[celli]
.opposingFaceLabel(facei, mesh_.faces());
const bool oppositeIntoOwner =
mesh_.faceOwner()[oppositeFacei] != celli;
if (oppositeFacei == -1)
{
FatalErrorInFunction
<< "Cannot find face on cell " << mesh_.cells()[celli]
<< " opposing face " << mesh_.faces()[facei]
<< ". Mesh is not layered." << exit(FatalError);
}
else if (!blockedFace[oppositeFacei])
{
blockedFace[oppositeFacei] = true;
newFrontFaces.append(oppositeFacei);
newFrontFaceIntoOwners.append(oppositeIntoOwner);
}
}
}
frontFaces.transfer(newFrontFaces);
frontFaceIntoOwners.transfer(newFrontFaceIntoOwners);
}
}
void Foam::functionObjects::layerAverage::calcLayers()
{
DynamicList<label> startFaces;
DynamicList<bool> startFaceIntoOwners;
forAll(patchIDs_, i)
{
const polyPatch& pp = mesh_.boundaryMesh()[patchIDs_[i]];
startFaces.append(identity(pp.size()) + pp.start());
startFaceIntoOwners.append(boolList(pp.size(), true));
}
forAll(zoneIDs_, i)
{
const faceZone& fz = mesh_.faceZones()[zoneIDs_[i]];
startFaces.append(fz);
startFaceIntoOwners.append(fz.flipMap());
}
// Identify faces that separate the layers
boolList blockedFace(mesh_.nFaces(), false);
boolList cellIsLayer(mesh_.nCells(), false);
walkOppositeFaces
(
startFaces,
startFaceIntoOwners,
blockedFace,
cellIsLayer
);
// Do analysis for connected layers
regionSplit rs(mesh_, blockedFace);
nLayers_ = rs.nRegions();
cellLayer_.transfer(rs);
// Get rid of regions that are not layers
label layeri0 = labelMax, layeri1 = -labelMax;
forAll(cellLayer_, celli)
{
if (cellIsLayer[celli])
{
layeri0 = min(layeri0, cellLayer_[celli]);
layeri1 = max(layeri1, cellLayer_[celli]);
}
else
{
cellLayer_[celli] = -1;
}
}
reduce(layeri0, minOp<label>());
reduce(layeri1, maxOp<label>());
nLayers_ = layeri0 != labelMax ? 1 + layeri1 - layeri0 : 0;
forAll(cellLayer_, celli)
{
if (cellLayer_[celli] != -1)
{
cellLayer_[celli] -= layeri0;
}
}
// Report
if (nLayers_ != 0)
{
Info<< " Detected " << nLayers_ << " layers" << nl << endl;
}
else
{
WarningInFunction
<< "No layers detected" << endl;
}
// Sum number of entries per layer
layerCount_ = sum(scalarField(mesh_.nCells(), 1));
// Average the cell centres
const pointField layerCentres(sum(mesh_.cellCentres())/layerCount_);
// Sort by the direction component
x_ = layerCentres.component(axis_);
sortMap_ = identity(layerCentres.size());
stableSort(sortMap_, UList<scalar>::less(x_));
inplaceReorder(sortMap_, x_);
// If symmetric, keep only half of the coordinates
if (symmetric_)
{
x_.setSize(nLayers_/2);
}
}
template<>
Foam::vector
Foam::functionObjects::layerAverage::symmetricCoeff<Foam::vector>() const
{
vector result = vector::one;
result[axis_] = -1;
return result;
}
template<>
Foam::symmTensor
Foam::functionObjects::layerAverage::symmetricCoeff<Foam::symmTensor>() const
{
return sqr(symmetricCoeff<vector>());
}
template<>
Foam::tensor
Foam::functionObjects::layerAverage::symmetricCoeff<Foam::tensor>() const
{
return symmetricCoeff<symmTensor>();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::functionObjects::layerAverage::layerAverage
(
const word& name,
const Time& runTime,
const dictionary& dict
)
:
fvMeshFunctionObject(name, runTime, dict)
{
read(dict);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::functionObjects::layerAverage::~layerAverage()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::functionObjects::layerAverage::read(const dictionary& dict)
{
Info<< type() << " " << name() << ":" << nl;
patchIDs_ =
mesh_.boundaryMesh().patchSet
(
dict.lookupOrDefault<wordReList>("patches", wordReList())
).toc();
zoneIDs_ =
findStrings
(
dict.lookupOrDefault<wordReList>("zones", wordReList()),
mesh_.faceZones().names()
);
if (patchIDs_.empty() && zoneIDs_.empty())
{
WarningInFunction
<< "No patches or zones specified" << endl;
}
symmetric_ = dict.lookupOrDefault<bool>("symmetric", false);
axis_ = axisNames_.read(dict.lookup("axis"));
fields_ = dict.lookup<wordList>("fields");
formatter_ = setWriter::New(dict.lookup("setFormat"), dict);
calcLayers();
return true;
}
Foam::wordList Foam::functionObjects::layerAverage::fields() const
{
return fields_;
}
bool Foam::functionObjects::layerAverage::execute()
{
return true;
}
bool Foam::functionObjects::layerAverage::write()
{
// Create list of available fields
wordList fieldNames;
forAll(fields_, fieldi)
{
if
(
false
#define FoundTypeField(Type, nullArg) \
|| foundObject<VolField<Type>>(fields_[fieldi])
FOR_ALL_FIELD_TYPES(FoundTypeField)
#undef FoundTypeField
)
{
fieldNames.append(fields_[fieldi]);
}
else
{
cannotFindObject(fields_[fieldi]);
}
}
// Collapse the fields
#define DeclareTypeValueSets(Type, nullArg) \
PtrList<Field<Type>> Type##ValueSets(fieldNames.size());
FOR_ALL_FIELD_TYPES(DeclareTypeValueSets);
#undef DeclareTypeValueSets
forAll(fieldNames, fieldi)
{
#define CollapseTypeFields(Type, nullArg) \
if (mesh_.foundObject<VolField<Type>>(fieldNames[fieldi])) \
{ \
const VolField<Type>& field = \
mesh_.lookupObject<VolField<Type>>(fieldNames[fieldi]); \
\
Type##ValueSets.set \
( \
fieldi, \
average(field.primitiveField()) \
); \
}
FOR_ALL_FIELD_TYPES(CollapseTypeFields);
#undef CollapseTypeFields
}
// Output directory
fileName outputPath =
mesh_.time().globalPath()/writeFile::outputPrefix/name();
if (mesh_.name() != fvMesh::defaultRegion)
{
outputPath = outputPath/mesh_.name();
}
// Write
formatter_->write
(
outputPath/mesh_.time().timeName(),
typeName,
coordSet(true, axisNames_[axis_], x_),
fieldNames
#define TypeValueSetsParameter(Type, nullArg) , Type##ValueSets
FOR_ALL_FIELD_TYPES(TypeValueSetsParameter)
#undef TypeValueSetsParameter
);
return true;
}
void Foam::functionObjects::layerAverage::updateMesh(const mapPolyMesh&)
{
Info<< type() << " " << name() << ":" << nl;
calcLayers();
}
void Foam::functionObjects::layerAverage::movePoints(const polyMesh&)
{
Info<< type() << " " << name() << ":" << nl;
calcLayers();
}
// ************************************************************************* //

243
src/functionObjects/field/layerAverage/layerAverage.H Normal file
View File

@ -0,0 +1,243 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::layerAverage
Description
Generates plots of fields averaged over the layers in the mesh
Example of function object specification:
\verbatim
layerAverage1
{
type layerAverage;
libs ("libfieldFunctionObjects.so");
writeControl writeTime;
setFormat raw;
patches (bottom);
zones (quarterPlane threeQuartersPlane);
component y;
symmetric true;
fields (pMean pPrime2Mean UMean UPrime2Mean k);
}
\endverbatim
Usage
\table
Property | Description | Required | Default value
type | Type name: layerAverage | yes |
setFormat | Output plotting format | yes |
patches | Patches that layers extrude from | no | ()
zones | Face zones that the layers extrude from | no | ()
component | Component of position to plot against | yes |
symmetric | Is the geometry symmetric around the centre layer? \
| no | false
field | Fields to average and plot | yes |
\endtable
SourceFiles
layerAverage.C
layerAverageTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef layerAverage_H
#define layerAverage_H
#include "fvMeshFunctionObject.H"
#include "setWriter.H"
#include "boolList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace functionObjects
{
/*---------------------------------------------------------------------------*\
Class layerAverage Declaration
\*---------------------------------------------------------------------------*/
class layerAverage
:
public fvMeshFunctionObject
{
// Private Static Data
//- Names of vector components
static const NamedEnum<vector::components, 3> axisNames_;
// Private Data
//- Patches which form the start of the layers
labelList patchIDs_;
//- Zones which form the start of the layers
labelList zoneIDs_;
//- Zones on which the layers are considered to end
labelList endZoneIDs_;
//- Is the case symmetric?
bool symmetric_;
//- The direction over which to plot the results
vector::components axis_;
//- Per cell the global layer
label nLayers_;
//- Per cell the global layer
labelList cellLayer_;
//- Per global layer the number of cells
scalarField layerCount_;
//- From sorted layer back to unsorted global layer
labelList sortMap_;
//- Sorted component of cell centres
scalarField x_;
//- Fields to sample
wordList fields_;
//- Set formatter
autoPtr<setWriter> formatter_;
// Private Member Functions
//- Walk through layers marking faces that separate layers and cells
// that are within layers
void walkOppositeFaces
(
const labelList& startFaces,
const boolList& startFaceIntoOwners,
boolList& blockedFace,
boolList& cellIsLayer
);
//- Create the layer information, the sort map, and the scalar axis
void calcLayers();
//- Return the coefficient to multiply onto symmetric values
template<class T>
T symmetricCoeff() const;
//- Sum field per layer
template<class T>
Field<T> sum(const Field<T>& cellField) const;
//- Average a field per layer
template<class T>
Field<T> average(const Field<T>& cellField) const;
public:
//- Runtime type information
TypeName("layerAverage");
// Constructors
//- Construct from Time and dictionary
layerAverage
(
const word& name,
const Time& runTime,
const dictionary& dict
);
//- Disallow default bitwise copy construction
layerAverage(const layerAverage&) = delete;
//- Destructor
virtual ~layerAverage();
// Member Functions
//- Read the field average data
virtual bool read(const dictionary&);
//- Return the list of fields required
virtual wordList fields() const;
//- Do nothing
virtual bool execute();
//- Calculate and write the graphs
virtual bool write();
//- Update for changes of mesh
virtual void updateMesh(const mapPolyMesh&);
//- Update for mesh point-motion
virtual void movePoints(const polyMesh&);
// Member Operators
//- Disallow default bitwise assignment
void operator=(const layerAverage&) = delete;
};
template<>
vector layerAverage::symmetricCoeff<vector>() const;
template<>
symmTensor layerAverage::symmetricCoeff<symmTensor>() const;
template<>
tensor layerAverage::symmetricCoeff<tensor>() const;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace functionObjects
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "layerAverageTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

89
src/functionObjects/field/layerAverage/layerAverageTemplates.C Normal file
View File

@ -0,0 +1,89 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "layerAverage.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class T>
T Foam::functionObjects::layerAverage::symmetricCoeff() const
{
return pTraits<T>::one;
}
template<class T>
Foam::Field<T> Foam::functionObjects::layerAverage::sum
(
const Field<T>& cellField
) const
{
Field<T> layerField(nLayers_, Zero);
forAll(cellLayer_, celli)
{
if (cellLayer_[celli] != -1)
{
layerField[cellLayer_[celli]] += cellField[celli];
}
}
Pstream::listCombineGather(layerField, plusEqOp<T>());
Pstream::listCombineScatter(layerField);
return layerField;
}
template<class T>
Foam::Field<T> Foam::functionObjects::layerAverage::average
(
const Field<T>& cellField
) const
{
// Sum, average and sort
Field<T> layerField(sum(cellField)/layerCount_, sortMap_);
// Handle symmetry
if (symmetric_)
{
const T coeff = symmetricCoeff<T>();
for (label i=0; i<nLayers_/2; i++)
{
const label j = nLayers_ - i - 1;
layerField[i] =
(layerField[i] + cmptMultiply(coeff, layerField[j]))/2;
}
layerField.setSize(nLayers_/2);
}
return layerField;
}
// ************************************************************************* //

7
test/postProcessing/channel/system/controlDict
View File

@ -248,6 +248,13 @@ graphsFunctions
end=(0.500001 0.500001 0),
fields=(p U)
)
#includeFunc graphLayerAverage
(
patches=(inlet),
axis=x,
fields=(p U)
)
}
lagrangianFunctions

2
tutorials/incompressible/pimpleFoam/LES/channel395/Allrun
View File

@ -17,6 +17,4 @@ runApplication decomposePar -cellDist
runParallel $application
runApplication reconstructPar
runApplication postChannel
#------------------------------------------------------------------------------

25
tutorials/incompressible/pimpleFoam/LES/channel395/system/controlDict
View File

@ -47,6 +47,31 @@ runTimeModifiable true;
functions
{
#includeFunc fieldAverage(U, p, prime2Mean = yes)
#includeFunc uniform
(
fieldType = volScalarField,
name = half,
dimensions = [0 0 0 0 0 0 0],
value = 0.5
)
#includeFunc mag(UPrime2Mean)
#includeFunc multiply(half, mag(UPrime2Mean), result = k)
#includeFunc graphLayerAverage
(
funcName = layerAverage,
patches = (bottomWall),
axis = y,
symmetric = yes,
pMean,
pPrime2Mean,
UMean,
UPrime2Mean,
k
);
}
// ************************************************************************* //