zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: dynamicRefineBalance: removed.

Browse Source
This commit is contained in:
mattijs
2018-11-15 15:03:19 +00:00
parent cbbafa7457
commit de2497461c
5 changed files with 17 additions and 1283 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/dynamicFvMesh/Make/files
View File

@ -5,7 +5,6 @@ dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
dynamicMultiMotionSolverFvMesh/dynamicMultiMotionSolverFvMesh.C
dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
dynamicRefineFvMesh/dynamicRefineFvMesh.C
dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
dynamicMotionSolverListFvMesh/dynamicMotionSolverListFvMesh.C
simplifiedDynamicFvMesh/simplifiedDynamicFvMeshes.C

843
src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.C
View File

@ -1,843 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 Tyler Voskuilen
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is a derivative work 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 "dynamicRefineBalancedFvMesh.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceInterpolate.H"
#include "volFields.H"
#include "polyTopoChange.H"
#include "surfaceFields.H"
#include "syncTools.H"
#include "pointFields.H"
#include "fvCFD.H"
#include "volPointInterpolation.H"
#include "pointMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(dynamicRefineBalancedFvMesh, 0);
addToRunTimeSelectionTable(dynamicFvMesh, dynamicRefineBalancedFvMesh, IOobject);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::dynamicRefineBalancedFvMesh::topParentID(label p)
{
label nextP = meshCutter().history().splitCells()[p].parent_;
if( nextP < 0 )
{
return p;
}
else
{
return topParentID(nextP);
}
}
Foam::List<Foam::scalar> Foam::dynamicRefineBalancedFvMesh::readRefinementPoints()
{
dictionary refineDict
(
IOdictionary
(
IOobject
(
"dynamicMeshDict",
time().constant(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).subDict("dynamicRefineFvMeshCoeffs")
);
List<scalar> refData(4, scalar(0));
refData[0] = readScalar(refineDict.lookup("unrefineLevel"));
refData[1] = readScalar(refineDict.lookup("lowerRefineLevel"));
refData[2] = readScalar(refineDict.lookup("refineInterval"));
refData[3] = readScalar(refineDict.lookup("maxRefinement"));
return refData;
}
void Foam::dynamicRefineBalancedFvMesh::updateRefinementField()
{
Info<< "Calculating internal refinement field" << endl;
volScalarField& intRefFld = *internalRefinementFieldPtr_;
volScalarField& targetFld = *targetLevelPtr_;
volScalarField& currFld = *isLevelPtr_;
// Set the internal refinement field to zero to start with
intRefFld = dimensionedScalar("zero",dimless,0.0);
// Get the cell level field from dynamicRefineFvMesh
const labelList& cellLevel = meshCutter().cellLevel();
// Read the points at which refinement and unrefinement occur from the
// dynamicMeshDict entries
List<scalar> refinePoints = readRefinementPoints();
// init list for refine/unrefine field (-1=unrefine, 1=refine, 0=do nothing)
labelList markRefineCells (this->nCells(), 0);
// init list for target refinement level per cell
labelList targetLevel (this->nCells(), 0);
// First fields
List<word> fieldNames = fields_.toc();
Field<scalar> refFld(nCells(),0.0);
forAll(fieldNames, i)
{
word fldName = fieldNames[i];
scalar minValue = fields_[fldName][0];
scalar maxValue = fields_[fldName][1];
label refineLevel = static_cast<label>(fields_[fldName][2]);
const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
// Limit the value of refFld based on its max level
forAll(fld, cellI)
{
if ((fld[cellI] >= minValue) && (fld[cellI] <= maxValue))
{
// increase targetLevel up to refineLevel
// BUT: do not decrease if cell already marked for higher refinement level by previous criterion
targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
}
}
}
// Then gradients
List<word> gradFieldNames = gradFields_.toc();
Field<scalar> cubeRtV = Foam::pow(this->V(),1.0/3.0);
forAll(gradFieldNames, i)
{
word fldName = gradFieldNames[i];
scalar minValue = gradFields_[fldName][0];
scalar maxValue = gradFields_[fldName][1];
label refineLevel = static_cast<label>(gradFields_[fldName][2]);
const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
refFld = mag(fvc::grad(fld)) * cubeRtV;
// Limit the value of refFld based on its max level
forAll(refFld, cellI)
{
if ((refFld[cellI] >= minValue) && (refFld[cellI] <= maxValue))
{
targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
}
}
}
// Then curls
List<word> curlFieldNames = curlFields_.toc();
forAll(curlFieldNames, i)
{
word fldName = curlFieldNames[i];
scalar minValue = curlFields_[fldName][0];
scalar maxValue = curlFields_[fldName][1];
label refineLevel = static_cast<label>(curlFields_[fldName][2]);
const volVectorField& fld = this->lookupObject<volVectorField>(fldName);
refFld = mag(fvc::curl(fld));
// Limit the value of refFld based on its max level
forAll(refFld, cellI)
{
if ((refFld[cellI] >= minValue) && (refFld[cellI] <= maxValue))
{
targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
}
}
}
// At the interface, assumed to be always the maximum refinement level
List<word> interfaceRefineField = interface_.toc();
forAll(interfaceRefineField, i)
{
word fldName = interfaceRefineField[i];
// read region of maximum refinement levels inside and outside of interface indicator field
// (does not need to be alpha, can also be a concentration field)
scalar innerRefLayers = readScalar(interface_[fldName].lookup("innerRefLayers"));
scalar outerRefLayers = readScalar(interface_[fldName].lookup("outerRefLayers"));
scalar nAddLayers(0);
if (interface_[fldName].found("nAddLayers"))
{
nAddLayers = readScalar(interface_[fldName].lookup("nAddLayers"));
}
label refineLevel = refinePoints[3];
if (interface_[fldName].found("maxRefineLevel"))
{
refineLevel = readScalar(interface_[fldName].lookup("maxRefineLevel"));
// to avoid user input mistakes, limit the value with the maximum allowed
refineLevel = min(refinePoints[3], refineLevel);
}
const volScalarField& fld = this->lookupObject<volScalarField>(fldName);
volScalarField isInterface(intRefFld * 0.0);
isInterface = dimensionedScalar("isInterface",dimless,0.0);
surfaceScalarField deltaAlpha = mag(fvc::snGrad(fld) / this->deltaCoeffs());
const labelUList& owner = this->owner();
const labelUList& neighbour = this->neighbour();
forAll(deltaAlpha, faceI)
{
//TODO: add a min max value of the specified field to be marked for refinement
// NOW: hard-coded method snGradAlpha: checks mag(alpha_N - alpha_P) > 0.1 ? 1:0
if (deltaAlpha[faceI] > 0.1) // currently a fixed prescribed value; should be read from díctionary
{
label own = owner[faceI];
label nei = neighbour[faceI];
// set isInterface field to one
isInterface[own] = 1.0;
isInterface[nei] = 1.0;
}
}
// assumed fld=0.5*(fldMax+fldMin) defines the interface
dimensionedScalar fldInterfaceValue(0.5*(max(fld)+min(fld)));
//-DD: old implementation based on face interpolation
// which results in slower transport in diagonal direction
// add inner refinement layers
//for(label i=0; i < innerRefLayers; i++)
//{
// isInterface += neg(- fvc::average(fvc::interpolate(isInterface)) * pos(fld - fldInterfaceValue));
// isInterface = neg(- isInterface);
//}
//
// add outer refinement layers
//for(label i=0; i < outerRefLayers; i++)
//{
// isInterface += neg(- fvc::average(fvc::interpolate(isInterface)) * pos(fldInterfaceValue - fld));
// isInterface = neg(- isInterface);
//}
//
//forAll(isInterface, cellI)
//{
// if (isInterface[cellI] > 0.5)
// {
// targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
// }
//}
//-DD: new version using volPointInterpolation (direction independent buffer layer)
const volPointInterpolation& pInterp = volPointInterpolation::New(*this);
const fvMesh& mesh = fld.mesh();
// add inner refinement layers
for(label i=0; i < innerRefLayers; i++)
{
volScalarField markInner(isInterface*pos(fld - fldInterfaceValue));
pointScalarField markLayerP(pInterp.interpolate(markInner));
forAll(mesh.C(), cellI)
{
scalar sum = 0.;
label nPoints = 0;
forAll(mesh.cellPoints()[cellI], pointI)
{
sum += markLayerP[mesh.cellPoints()[cellI][pointI]];
nPoints++;
}
if (nPoints > 0)
{
sum /= nPoints;
}
isInterface[cellI] += sum;
}
}
isInterface = pos(isInterface - SMALL);
// add outer refinement layers
for(label i=0; i < outerRefLayers; i++)
{
volScalarField markOuter(isInterface*pos(fldInterfaceValue - fld));
pointScalarField markLayerP(pInterp.interpolate(markOuter));
forAll(mesh.C(), cellI)
{
scalar sum = 0.;
label nPoints = 0;
forAll(mesh.cellPoints()[cellI], pointI)
{
sum += markLayerP[mesh.cellPoints()[cellI][pointI]];
nPoints++;
}
if (nPoints > 0)
{
sum /= nPoints;
}
isInterface[cellI] += sum;
}
}
isInterface = pos(isInterface - SMALL);
forAll(isInterface, cellI)
{
if (isInterface[cellI] > 0.5)
{
targetLevel[cellI] = max(targetLevel[cellI], refineLevel);
}
}
//-DD: old implementation based on face interpolation
// which results in slower transport in diagonal direction
// expand additional layers if specified:
if (nAddLayers > 0)
{
// loop over additional layers
for(label i=1; i < refineLevel; i++)
{
// #nAddLayers cells per additional level
for(label j=0; j < nAddLayers; j++)
{
isInterface += neg(- fvc::average(fvc::interpolate(isInterface)));
isInterface = neg(- isInterface);
}
forAll(isInterface, cellI)
{
if (isInterface[cellI] == 1.0)
{
targetLevel[cellI] = max(targetLevel[cellI], (refineLevel - i));
}
}
}
}
}
// The set refinement physical regions (force the mesh to stay refined
// near key features)
forAll(refinedRegions_, regionI)
{
const entry& region = refinedRegions_[regionI];
autoPtr<topoSetSource> source =
topoSetSource::New(region.keyword(), *this, region.dict());
cellSet selectedCellSet
(
*this,
"cellSet",
nCells()/10+1 //estimate
);
source->applyToSet
(
topoSetSource::NEW,
selectedCellSet
);
const labelList cells = selectedCellSet.toc();
label minLevel = readLabel(region.dict().lookup("minLevel"));
forAll(cells, i)
{
const label& cellI = cells[i];
targetLevel[cellI] = max(targetLevel[cellI], minLevel);
}
}
//-DD: add buffer layer based on targetLevel field to prevent 2-to-1 refinement
volScalarField blockedLevel = targetFld * 0.;
for (label currLayer=refinePoints[3]; currLayer>=1; currLayer--)
{
forAll (targetLevel, cellI)
{
if (targetLevel[cellI] >= currLayer)
{
blockedLevel[cellI] = targetLevel[cellI];
}
}
for (label i=0; i<nBufferLayers_; i++)
{
blockedLevel = max(blockedLevel, pos(fvc::average(fvc::interpolate(blockedLevel)) - SMALL)*(currLayer-1));
}
labelList blockLev(blockedLevel.internalField().size(), 0);
forAll (blockLev, cellI)
{
blockLev[cellI] = blockedLevel[cellI];
}
targetLevel = max(targetLevel, blockLev);
}
// mark cells to be refined/unrefined based on above criteria:
// simply check, if targetLevel lower or higher cellLevel
forAll (intRefFld.internalField(), cellI)
{
intRefFld.primitiveFieldRef()[cellI] = targetLevel[cellI] - cellLevel[cellI];
targetFld.primitiveFieldRef()[cellI] = targetLevel[cellI];
currFld.primitiveFieldRef()[cellI] = cellLevel[cellI];
}
intRefFld.correctBoundaryConditions();
Info<<"Min,max refinement field = " << Foam::min(intRefFld).value() << ", "
<< Foam::max(intRefFld).value() << endl;
}
void Foam::dynamicRefineBalancedFvMesh::readRefinementDict()
{
dictionary dynamicMeshDict
(
IOdictionary
(
IOobject
(
"dynamicMeshDict",
time().constant(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
)
);
const dictionary refineDict
(
IOdictionary
(
IOobject
(
"dynamicMeshDict",
time().constant(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).subDict("dynamicRefineFvMeshCoeffs")
);
{
wordList surfFlds;
if (refineDict.readIfPresent("mapSurfaceFields", surfFlds))
{
// Rework into hashtable.
mapSurfaceFields_.resize(surfFlds.size());
forAll(surfFlds, i)
{
mapSurfaceFields_.insert(surfFlds[i], surfFlds[i]);
}
}
}
nBufferLayers_ = readLabel(refineDict.lookup("nBufferLayers"));
if( dynamicMeshDict.isDict("refinementControls") )
{
dictionary refineControlDict =
dynamicMeshDict.subDict("refinementControls");
enableRefinementControl_ =
Switch(refineControlDict.lookup("enableRefinementControl"));
if( enableRefinementControl_ )
{
// Overwrite field name entry in dynamicRefineFvMeshCoeffs?
// For now you just have to be smart and enter
// 'internalRefinementField' for the field name manually
// Read HashTable of field-refinement scalars
if( refineControlDict.found("fields") )
{
fields_ = HashTable< List<scalar> >
(
refineControlDict.lookup("fields")
);
}
// Read HashTable of gradient-refinement scalars
if( refineControlDict.found("gradients") )
{
gradFields_ = HashTable< List<scalar> >
(
refineControlDict.lookup("gradients")
);
}
// Read HashTable of curl-refinement vectors
if( refineControlDict.found("curls") )
{
curlFields_ = HashTable< List<scalar> >
(
refineControlDict.lookup("curls")
);
}
// Read interface refinement data
if( refineControlDict.found("interface") )
{
interface_ = HashTable< dictionary >
(
refineControlDict.lookup("interface")
);
}
// Read refinement regions
if( refineControlDict.found("regions") )
{
refinedRegions_ = PtrList<entry>
(
refineControlDict.lookup("regions")
);
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::dynamicRefineBalancedFvMesh::dynamicRefineBalancedFvMesh
(
const IOobject& io
)
:
dynamicRefineFvMesh(io),
internalRefinementFieldPtr_(NULL),
targetLevelPtr_(NULL),
isLevelPtr_(NULL),
fields_(),
gradFields_(),
curlFields_(),
interface_(),
refinedRegions_(),
enableRefinementControl_(false),
nBufferLayers_(0),
rebalance_(false)
{
readRefinementDict();
if( enableRefinementControl_ )
{
internalRefinementFieldPtr_ = new volScalarField
(
IOobject
(
"internalRefinementField",
time().timeName(),
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
*this,
dimensionedScalar("zero", dimless, 0.0)
);
targetLevelPtr_ = new volScalarField
(
IOobject
(
"targetLevel",
time().timeName(),
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
*this,
dimensionedScalar("zero", dimless, 0.0)
);
isLevelPtr_ = new volScalarField
(
IOobject
(
"isLevel",
time().timeName(),
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
*this,
dimensionedScalar("zero", dimless, 0.0)
);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::dynamicRefineBalancedFvMesh::~dynamicRefineBalancedFvMesh()
{
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::dynamicRefineBalancedFvMesh::mapFields(const mapPolyMesh& mpm)
{
DebugVar(mpm.nOldCells());
dynamicRefineFvMesh::mapFields(mpm);
// Correct surface fields on introduced internal faces. These get
// created out-of-nothing so get an interpolated value.
mapNewInternalFaces<scalar>(mpm.faceMap());
mapNewInternalFaces<vector>(mpm.faceMap());
mapNewInternalFaces<sphericalTensor>(mpm.faceMap());
mapNewInternalFaces<symmTensor>(mpm.faceMap());
mapNewInternalFaces<tensor>(mpm.faceMap());
}
bool Foam::dynamicRefineBalancedFvMesh::update()
{
//Part 0 - Update internally calculated refinement field
readRefinementDict();
List<scalar> refinePoints = readRefinementPoints();
label refineInterval = refinePoints[2];
if( enableRefinementControl_ && time().timeIndex() > 0 && time().timeIndex() % refineInterval == 0 )
{
updateRefinementField();
}
//Part 1 - Call normal update from dynamicRefineFvMesh
bool hasChanged = dynamicRefineFvMesh::update();
if( Pstream::parRun() && hasChanged)
{
//Correct values on all coupled patches
correctBoundaries<scalar>();
correctBoundaries<vector>();
correctBoundaries<sphericalTensor>();
correctBoundaries<symmTensor>();
correctBoundaries<tensor>();
}
dictionary decomposeParDict
(
IOdictionary
(
IOobject
(
"decomposeParDict",
time().system(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
)
);
rebalance_ = false;
// Part 2 - Load Balancing
{
dictionary refineDict
(
IOdictionary
(
IOobject
(
"dynamicMeshDict",
time().constant(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE,
false
)
).subDict("dynamicRefineFvMeshCoeffs")
);
Switch enableBalancing = refineDict.lookup("enableBalancing");
if ( Pstream::parRun() && hasChanged )
{
const scalar allowableImbalance =
readScalar(refineDict.lookup("allowableImbalance"));
//First determine current level of imbalance - do this for all
// parallel runs with a changing mesh, even if balancing is disabled
label nGlobalCells = globalData().nTotalCells();
scalar idealNCells = scalar(nGlobalCells)/scalar(Pstream::nProcs());
scalar localImbalance = mag(scalar(nCells()) - idealNCells);
Foam::reduce(localImbalance, maxOp<scalar>());
scalar maxImbalance = localImbalance/idealNCells;
Info<<"Maximum imbalance = " << 100*maxImbalance << " %" << endl;
//If imbalanced, construct weighted coarse graph (level 0) with node
// weights equal to their number of subcells. This partitioning works
// as long as the number of level 0 cells is several times greater than
// the number of processors.
if( maxImbalance > allowableImbalance && enableBalancing)
{
Info << "\n**Solver hold for redistribution at time = " << time().timeName() << " s" << endl;
rebalance_ = true;
const labelIOList& cellLevel = meshCutter().cellLevel();
Map<label> coarseIDmap(100);
labelList uniqueIndex(nCells(),0);
label nCoarse = 0;
forAll(cells(), cellI)
{
if( cellLevel[cellI] > 0 )
{
uniqueIndex[cellI] = nCells() + topParentID
(
meshCutter().history().parentIndex(cellI)
);
}
else
{
uniqueIndex[cellI] = cellI;
}
if( coarseIDmap.insert(uniqueIndex[cellI], nCoarse) )
{
++nCoarse;
}
}
// Convert to local sequential indexing and calculate coarse
// points and weights
labelList localIndex(nCells(),0);
pointField coarsePoints(nCoarse,vector::zero);
scalarField coarseWeights(nCoarse,0.0);
forAll(uniqueIndex, cellI)
{
localIndex[cellI] = coarseIDmap[uniqueIndex[cellI]];
// If 2D refinement (quadtree) is ever implemented, this '3'
// should be set in general as the number of refinement
// dimensions.
label w = (1 << (3*cellLevel[cellI]));
coarseWeights[localIndex[cellI]] += 1.0;
coarsePoints[localIndex[cellI]] += C()[cellI]/w;
}
// Set up decomposer - a separate dictionary is used here so
// you can use a simple partitioning for decomposePar and
// ptscotch for the rebalancing (or any chosen algorithms)
autoPtr<decompositionMethod> decomposer
(
decompositionMethod::New
(
IOdictionary
(
IOobject
(
"balanceParDict",
time().system(),
*this,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
)
)
);
labelList finalDecomp = decomposer().decompose
(
*this,
localIndex,
coarsePoints,
coarseWeights
);
scalar tolDim = globalMeshData::matchTol_ * bounds().mag();
Info<< "Distributing the mesh ..." << endl;
fvMeshDistribute distributor(*this, tolDim);
Info<< "Mapping the fields ..." << endl;
autoPtr<mapDistributePolyMesh> map =
distributor.distribute(finalDecomp);
Info<< "Distribute the map ..." << endl;
meshCutter_.distribute(map);
Info << "Successfully distributed mesh" << endl;
scalarList procLoadNew (Pstream::nProcs(), 0.0);
procLoadNew[Pstream::myProcNo()] = this->nCells();
reduce(procLoadNew, sumOp<List<scalar> >());
scalar overallLoadNew = sum(procLoadNew);
scalar averageLoadNew = overallLoadNew/double(Pstream::nProcs());
Info << "New distribution: " << procLoadNew << endl;
Info << "Max deviation: " << max(Foam::mag(procLoadNew-averageLoadNew)/averageLoadNew)*100.0 << " %" << endl;
}
}
}
if( Pstream::parRun() && rebalance_)
{
//Correct values on all coupled patches
correctBoundaries<scalar>();
correctBoundaries<vector>();
correctBoundaries<sphericalTensor>();
correctBoundaries<symmTensor>();
correctBoundaries<tensor>();
}
return hasChanged;
}
// ************************************************************************* //

206
src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMesh.H
View File

@ -1,206 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 Tyler Voskuilen
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is a derivative work 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::dynamicRefineBalancedFvMesh
SourceFiles
dynamicRefineBalancedFvMesh.C
dynamicRefineBalancedFvMeshTemplates.C
Authors
T.G. Voskuilen ( https://github.com/tgvoskuilen/meshBalancing )
Daniel Deising <deising@mma.tu-darmstadt.de>
Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
All rights reserved.
Description
A fvMesh with built-in multi-criterion refinement
and run-time load balancing.
You may refer to this software as :
//- full bibliographic data to be provided
This code has been developed by :
Daniel Rettenmaier (main developer).
Method Development and Intellectual Property :
T.G. Voskuilen (Purdue University)
Timothée Pourpoint <timothee@purdue.edu> (Purdue University
Daniel Rettenmaier <rettenmaier@gsc.tu-darmstadt.de>
Daniel Deising <deising@mma.tu-darmstadt.de>
Holger Marschall <marschall@csi.tu-darmstadt.de>
Dieter Bothe <bothe@csi.tu-darmstadt.de>
Cameron Tropea <ctropea@sla.tu-darmstadt.de>
School of Aeronautics and Astronautics
Purdue University
Mathematical Modeling and Analysis
Institute for Fluid Mechanics and Aerodynamics
Center of Smart Interfaces
Technische Universitaet Darmstadt
If you use this software for your scientific work or your publications,
please don't forget to acknowledge explicitly the use of it.
\*---------------------------------------------------------------------------*/
#ifndef dynamicRefineBalancedFvMesh_H
#define dynamicRefineBalancedFvMesh_H
#include "dynamicFvMesh.H"
#include "dynamicRefineFvMesh.H"
#include "hexRef8.H"
#include "PackedBoolList.H"
#include "Switch.H"
#include "decompositionMethod.H"
#include "fvMeshDistribute.H"
#include "mapDistributePolyMesh.H"
#include "HashTable.H"
#include "topoSetSource.H"
#include "cellSet.H"
#include "PtrDictionary.H"
#include "dictionaryEntry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class dynamicRefineBalancedFvMesh Declaration
\*---------------------------------------------------------------------------*/
class dynamicRefineBalancedFvMesh
:
public dynamicRefineFvMesh
{
private:
//- (non-flux)SurfaceFields to map
HashTable<word> mapSurfaceFields_;
//-
volScalarField* internalRefinementFieldPtr_;
volScalarField* targetLevelPtr_;
volScalarField* isLevelPtr_;
//-
HashTable< List<scalar> > fields_;
//-
HashTable< List<scalar> > gradFields_;
//-
HashTable< List<scalar> > curlFields_;
//-
HashTable< dictionary > interface_;
//-
PtrList<entry> refinedRegions_;
//-
Switch enableRefinementControl_;
//-
void updateRefinementField();
//-
void readRefinementDict();
//-
List<scalar> readRefinementPoints();
//- Disallow default bitwise copy construct
dynamicRefineBalancedFvMesh(const dynamicRefineBalancedFvMesh&);
//- Disallow default bitwise assignment
void operator=(const dynamicRefineBalancedFvMesh&);
//-
label topParentID(label p);
//-
label nBufferLayers_;
//-
bool rebalance_;
//- Map non-flux surface<Type>Fields for new internal faces
// (from cell splitting)
template <class T>
void mapNewInternalFaces(const labelList& faceMap);
public:
//- Runtime type information
TypeName("dynamicRefineBalancedFvMesh");
// Constructors
//- Construct from IOobject
explicit dynamicRefineBalancedFvMesh(const IOobject& io);
//- Destructor
virtual ~dynamicRefineBalancedFvMesh();
// Member Functions
//- Update the mesh for both mesh motion and topology change
virtual bool update();
//- Map all fields in time using given map
virtual void mapFields(const mapPolyMesh& mpm);
//- Template to update all volField boundaries
template<class Type> void correctBoundaries();
//-
bool rebalance()
{
return rebalance_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "dynamicRefineBalancedFvMeshTemplates.C"
#endif
#endif
// ************************************************************************* //

217
src/dynamicFvMesh/dynamicRefineBalancedFvMesh/dynamicRefineBalancedFvMeshTemplates.C
View File

@ -1,217 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2014 Tyler Voskuilen
\\/ 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 "volMesh.H"
#include "fvPatchField.H"
#include "surfaceFields.H"
template <class T>
void Foam::dynamicRefineBalancedFvMesh::mapNewInternalFaces
(
const labelList& faceMap
)
{
typedef GeometricField<T, fvsPatchField, surfaceMesh> GeoField;
HashTable<GeoField*> sFlds(this->objectRegistry::lookupClass<GeoField>());
const labelUList& owner = this->faceOwner();
const labelUList& neighbour = this->faceNeighbour();
const dimensionedScalar deltaN = 1e-8 / pow(average(this->V()), 1.0 / 3.0);
forAllIter(typename HashTable<GeoField*>, sFlds, iter)
{
GeoField& sFld = *iter();
if (mapSurfaceFields_.found(iter.key()))
{
if (debug)
{
InfoInFunction << iter.key()<< endl;
}
// Create flat version of field
Field<T> tsFld(this->nFaces(), pTraits<T>::zero);
{
forAll(sFld.internalField(), iFace)
{
tsFld[iFace] = sFld.internalField()[iFace];
}
forAll(sFld.boundaryField(), iPatch)
{
const fvsPatchField<T>& pfld = sFld.boundaryField()[iPatch];
label start = pfld.patch().start();
forAll(pfld, faceI)
{
tsFld[faceI+start] = pfld[faceI];
}
}
}
// Loop over all faces
for (label facei = 0; facei < nInternalFaces(); facei++)
{
label oldFacei = faceMap[facei];
//- map surface field on newly generated faces
if (oldFacei == -1)
{
//- find owner and neighbour cell
const cell& ownFaces = this->cells()[owner[facei]];
const cell& neiFaces = this->cells()[neighbour[facei]];
//- loop over all owner/neighbour cell faces
//- and find already mapped ones (master-faces):
T tmpValue = pTraits<T>::zero;
scalar magFld = 0;
label counter = 0;
//- simple averaging of all neighbour master-faces
forAll(ownFaces, iFace)
{
if (faceMap[ownFaces[iFace]] != -1)
{
tmpValue += tsFld[ownFaces[iFace]];
magFld += mag(tsFld[ownFaces[iFace]]);
counter++;
}
}
forAll(neiFaces, iFace)
{
if (faceMap[neiFaces[iFace]] != -1)
{
tmpValue = tmpValue + tsFld[neiFaces[iFace]];
magFld += mag(tsFld[neiFaces[iFace]]);
counter++;
}
}
if(counter > 0)
{
if
(
GeometricField<T, fvsPatchField, surfaceMesh>::typeName
== "surfaceScalarField"
)
{
tmpValue /= counter;
}
else if
(
GeometricField<T, fvsPatchField, surfaceMesh>::typeName
== "surfaceVectorField"
)
{
magFld /= counter;
tmpValue *= magFld/(mag(tmpValue)+deltaN.value());
}
else
{
FatalErrorInFunction
<< "mapping implementation only valid for"
<< " scalar and vector fields! \n Field "
<< sFld.name() << " is of type: "
<< GeometricField<T, fvsPatchField, surfaceMesh>::typeName
<< abort(FatalError);
}
}
sFld[facei] = tmpValue;
}
}
}
}
}
template<class Type>
void Foam::dynamicRefineBalancedFvMesh::correctBoundaries()
{
typedef GeometricField<Type, fvPatchField, volMesh> GeoField;
HashTable<GeoField*> flds(this->objectRegistry::lookupClass<GeoField>());
forAllIter(typename HashTable<GeoField*>, flds, iter)
{
GeoField& fld = *iter();
//mimic "evaluate" but only for coupled patches (processor or cyclic)
// and only for blocking or nonBlocking comms (no scheduled comms)
if
(
Pstream::defaultCommsType == Pstream::commsTypes::blocking
|| Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
)
{
label nReq = Pstream::nRequests();
forAll(fld.boundaryField(), patchi)
{
if(fld.boundaryField()[patchi].coupled())
{
fld.boundaryFieldRef()[patchi].initEvaluate
(
Pstream::defaultCommsType
);
}
}
// Block for any outstanding requests
if
(
Pstream::parRun()
&& Pstream::defaultCommsType == Pstream::commsTypes::nonBlocking
)
{
Pstream::waitRequests(nReq);
}
forAll(fld.boundaryField(), patchi)
{
if(fld.boundaryField()[patchi].coupled())
{
fld.boundaryFieldRef()[patchi].evaluate
(
Pstream::defaultCommsType
);
}
}
}
else
{
//Scheduled patch updates not supported
FatalErrorIn
(
"dynamicRefineBalancedFvMeshTemplates::correctBoundaries"
) << "Unsuported communications type "
<< Pstream::commsTypeNames[Pstream::defaultCommsType]
<< exit(FatalError);
}
}
}
// ************************************************************************* //

33
src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
View File

@ -46,7 +46,7 @@ namespace Foam
Foam::label Foam::dynamicRefineFvMesh::count
(
const PackedBoolList& l,
const bitSet& l,
const unsigned int val
)
{
@ -72,7 +72,7 @@ Foam::label Foam::dynamicRefineFvMesh::count
void Foam::dynamicRefineFvMesh::calculateProtectedCells
(
PackedBoolList& unrefineableCell
bitSet& unrefineableCell
) const
{
if (protectedCell_.empty())
@ -470,7 +470,7 @@ Foam::dynamicRefineFvMesh::refine
// Update numbering of protectedCell_
if (protectedCell_.size())
{
PackedBoolList newProtectedCell(nCells());
bitSet newProtectedCell(nCells());
forAll(newProtectedCell, celli)
{
@ -648,7 +648,7 @@ Foam::dynamicRefineFvMesh::unrefine
// Update numbering of protectedCell_
if (protectedCell_.size())
{
PackedBoolList newProtectedCell(nCells());
bitSet newProtectedCell(nCells());
forAll(newProtectedCell, celli)
{
@ -751,7 +751,7 @@ void Foam::dynamicRefineFvMesh::selectRefineCandidates
const scalar lowerRefineLevel,
const scalar upperRefineLevel,
const scalarField& vFld,
PackedBoolList& candidateCell
bitSet& candidateCell
) const
{
// Get error per cell. Is -1 (not to be refined) to >0 (to be refined,
@ -784,7 +784,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
(
const label maxCells,
const label maxRefinement,
const PackedBoolList& candidateCell
const bitSet& candidateCell
) const
{
// Every refined cell causes 7 extra cells
@ -794,7 +794,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
// Mark cells that cannot be refined since they would trigger refinement
// of protected cells (since 2:1 cascade)
PackedBoolList unrefineableCell;
bitSet unrefineableCell;
calculateProtectedCells(unrefineableCell);
// Count current selection
@ -871,7 +871,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectRefineCells
Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
(
const scalar unrefineLevel,
const PackedBoolList& markedCell,
const bitSet& markedCell,
const scalarField& pFld
) const
{
@ -884,7 +884,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
// If we have any protected cells make sure they also are not being
// unrefined
PackedBoolList protectedPoint(nPoints());
bitSet protectedPoint(nPoints());
if (protectedCell_.size())
{
@ -976,7 +976,7 @@ Foam::labelList Foam::dynamicRefineFvMesh::selectUnrefinePoints
void Foam::dynamicRefineFvMesh::extendMarkedCells
(
PackedBoolList& markedCell
bitSet& markedCell
) const
{
// Mark faces using any marked cell
@ -1018,7 +1018,7 @@ void Foam::dynamicRefineFvMesh::extendMarkedCells
void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
(
PackedBoolList& protectedCell,
bitSet& protectedCell,
label& nProtected
) const
{
@ -1343,7 +1343,7 @@ bool Foam::dynamicRefineFvMesh::update()
refineDict.get<label>("nBufferLayers");
// Cells marked for refinement or otherwise protected from unrefinement.
PackedBoolList refineCell(nCells());
bitSet refineCell(nCells());
// Determine candidates for refinement (looking at field only)
selectRefineCandidates
@ -1384,7 +1384,7 @@ bool Foam::dynamicRefineFvMesh::update()
const labelList& cellMap = map().cellMap();
const labelList& reverseCellMap = map().reverseCellMap();
PackedBoolList newRefineCell(cellMap.size());
bitSet newRefineCell(cellMap.size());
forAll(cellMap, celli)
{
@ -1471,7 +1471,8 @@ bool Foam::dynamicRefineFvMesh::writeObject
(
IOstream::streamFormat fmt,
IOstream::versionNumber ver,
IOstream::compressionType cmp
IOstream::compressionType cmp,
const bool valid
) const
{
// Force refinement data to go to the current time directory.
@ -1479,8 +1480,8 @@ bool Foam::dynamicRefineFvMesh::writeObject
bool writeOk =
(
dynamicFvMesh::writeObject(fmt, ver, cmp)
&& meshCutter_.write()
dynamicFvMesh::writeObject(fmt, ver, cmp, valid)
&& meshCutter_.write(valid)
);
if (dumpLevel_)