exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Manually merged Henrik's DLB and some other stuff

Browse source
This commit is contained in:
Vuko Vukcevic 2019-01-23 13:03:44 +01:00
parent 75a6813d3a
commit 920ff7bcc4
58 changed files with 1225 additions and 209 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/conjugateHeatTransfer/fvPatchFields/chtRcThermalDiffusivity/chtRcThermalDiffusivityFvPatchScalarField.C
View file

@ -261,7 +261,7 @@ Foam::chtRcThermalDiffusivityFvPatchScalarField::calcThermalDiffusivity
const scalarField kOwn = fOwn/(1.0 - p.weights())/mld.magDelta(p.index());
const scalarField kNei = fNei/p.weights()/mld.magDelta(p.index());
tmp<scalarField> kTmp(new scalarField(p.size()));
scalarField& k = kTmp();

20
src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.C
View file

@ -34,6 +34,7 @@ License
#include "OSspecific.H"
#include "Map.H"
#include "DynamicList.H"
#include "globalMeshData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -53,6 +54,7 @@ Foam::domainDecomposition::domainDecomposition
nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
distributed_(false),
gfIndex_(mesh_),
gpIndex_(mesh_),
cellToProc_(mesh_.nCells()),
patchNbrCellToProc_(mesh_.boundaryMesh().size()),
procPointAddressing_(nProcs_),
@ -503,6 +505,24 @@ Foam::autoPtr<Foam::fvMesh> Foam::domainDecomposition::processorMesh
}
Foam::labelList Foam::domainDecomposition::globalPointIndex
(
const label procI
) const
{
const labelList& gppi = gpIndex_.globalLabel();
const labelList& ppAddr = procPointAddressing_[procI];
labelList globalPointIndex(ppAddr.size());
forAll (globalPointIndex, gpI)
{
globalPointIndex[gpI] = gppi[ppAddr[gpI]];
}
return globalPointIndex;
}
bool Foam::domainDecomposition::writeDecomposition()
{
Info<< "\nConstructing processor meshes" << endl;

7
src/decompositionMethods/decomposeReconstruct/decomposeTools/finiteVolume/domainDecomposition.H
View file

@ -43,6 +43,7 @@ SourceFiles
#include "PtrList.H"
#include "point.H"
#include "globalProcFaceIndex.H"
#include "globalProcPointIndex.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -72,6 +73,9 @@ class domainDecomposition
//- Global face index
globalProcFaceIndex gfIndex_;
//- Global face index
globalProcPointIndex gpIndex_;
//- Processor label for each cell
labelList cellToProc_;
@ -205,6 +209,9 @@ public:
const bool createPassiveProcPatches = false
) const;
//- Create and return global point index
labelList globalPointIndex(const label procI) const;
//- Return processor point addressing
const labelListList& procPointAddressing() const
{

8
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/fvFieldReconstructorReconstructFields.C
View file

@ -44,10 +44,10 @@ void Foam::fvFieldReconstructor::reconstructField
) const
{
// Get references to internal and boundary field
Field<Type>& iField = reconField.internalField();
Field<Type>& iField = reconField;
typename GeometricField<Type, fvPatchField, volMesh>::
GeometricBoundaryField& bouField = reconField.boundaryField();
GeometricBoundaryField& bouField = reconField.boundaryFieldNoStoreOldTimes();
forAll (procFields, procI)
{
@ -143,10 +143,10 @@ void Foam::fvFieldReconstructor::reconstructField
) const
{
// Create the internalField
Field<Type>& iField = reconField.internalField();
Field<Type>& iField = reconField;
typename GeometricField<Type, fvsPatchField, surfaceMesh>::
GeometricBoundaryField& bouField = reconField.boundaryField();
GeometricBoundaryField& bouField = reconField.boundaryFieldNoStoreOldTimes();
forAll (procFields, procI)
{

28
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesRebuild.C
View file

@ -652,10 +652,30 @@ Foam::processorMeshesReconstructor::reconstructMesh(const Time& db)
reconPatchSizes = 0;
// Prepare handling for globally shared points. This is equivalent
// to parallel processor points, but working on a PtrList of meshes
// on the same processor
sharedPoints sharedData(meshes_);
// HR 21.12.18 : A bit of a hack for the moment in order to support the new
// method (using global point addressing) and old method (syncing across
// processor BCs) of constructing the shared points. The old method is still
// needed since global point addressing does not exist when the sharedPoint
// object is constructed from reconstructParMesh.
//
// TODO: Unify the methods by constructing global point addressing from the
// mesh pieces when. Or even better, calculate procPointAddressing directly
// which will simplify the mesh merging immensely.
autoPtr<sharedPoints> sharedDataPtr;
if(globalPointIndex_.size())
{
// Determine globally shared points using global point
// adressing
sharedDataPtr.set(new sharedPoints(meshes_, globalPointIndex_));
}
else
{
// Prepare handling for globally shared points. This is equivalent
// to parallel processor points, but working on a PtrList of meshes
// on the same processor
sharedDataPtr.set(new sharedPoints(meshes_));
}
sharedPoints& sharedData = sharedDataPtr();
// Record global point index for shared points
labelList globalPointMapping(sharedData.nGlobalPoints(), -1);

2
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.C
View file

@ -71,6 +71,7 @@ Foam::processorMeshesReconstructor::processorMeshesReconstructor
:
meshName_(meshName),
meshes_(),
globalPointIndex_(),
pointProcAddressing_(),
faceProcAddressing_(),
cellProcAddressing_(),
@ -86,6 +87,7 @@ Foam::processorMeshesReconstructor::processorMeshesReconstructor
:
meshName_(meshName),
meshes_(databases.size()),
globalPointIndex_(),
pointProcAddressing_(),
faceProcAddressing_(),
cellProcAddressing_(),

9
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/processorMeshesReconstructor.H
View file

@ -72,6 +72,9 @@ class processorMeshesReconstructor
//- List of processor meshes
PtrList<fvMesh> meshes_;
//- global point index per sub-mesh
labelListList globalPointIndex_;
//- List of processor point addressing lists
PtrList<labelIOList> pointProcAddressing_;
@ -161,6 +164,12 @@ public:
return meshes_;
}
//- Return access to meshes
labelListList& globalPointIndex()
{
return globalPointIndex_;
}
//- Return point-processor addressing arrays
const PtrList<labelIOList>& pointProcAddressing() const;

171
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.C
View file

@ -451,9 +451,160 @@ void Foam::sharedPoints::calcSharedPoints()
}
void Foam::sharedPoints::calcSharedPoints
(
const labelListList& globalPointIndex
)
{
typedef HashTable<Pair<label>, label, Hash<label> > markTable;
markTable marker;
List<labelHashSet> procSharedPoints(meshes_.size());
// Mark up points for the first time
forAll (meshes_, meshI)
{
if (meshes_.set(meshI))
{
const labelList& curGlobalPointIndex = globalPointIndex[meshI];
labelHashSet& sharedPoints = procSharedPoints[meshI];
forAll (curGlobalPointIndex, pI)
{
const label gpi = curGlobalPointIndex[pI];
markTable::iterator iter = marker.find(gpi);
if (iter == Map<label>::end())
{
// Record meshI and pI
marker.insert(gpi, Pair<label>(-meshI - 1, pI));
}
else
{
if (iter().first() < 0)
{
if(iter().first() != -meshI - 1)
{
// Shared point detected. Add for both meshes
// Add for first mesh
const label firstMesh = -iter().first() - 1;
const label firstPointI = iter().second();
procSharedPoints[firstMesh].insert(firstPointI);
// Add for current mesh
sharedPoints.insert(pI);
// Count shared points and update bookkeeping
iter().first() = nGlobalPoints_;
nGlobalPoints_++;
}
}
else
{
// Existing shared point. Add for current mesh
sharedPoints.insert(pI);
}
}
}
}
}
forAll (meshes_, meshI)
{
if (meshes_.set(meshI))
{
const labelList& curGlobalPointIndex = globalPointIndex[meshI];
labelList& curSharedPoints = sharedPointLabels_[meshI];
curSharedPoints = procSharedPoints[meshI].toc();
labelList& curSharedAddr = sharedPointAddr_[meshI];
curSharedAddr.setSize(curSharedPoints.size());
forAll(curSharedPoints, i)
{
const label gpi = curSharedPoints[i];
curSharedAddr[i] = marker.find(curGlobalPointIndex[gpi])().first();
}
}
}
// debug
{
pointField gp(nGlobalPoints_);
boolList gpSet(nGlobalPoints_, false);
forAll (meshes_, meshI)
{
if (meshes_.set(meshI))
{
// Get points
const pointField& P = meshes_[meshI].points();
// Get list of local point labels that are globally shared
const labelList& curShared = sharedPointLabels_[meshI];
// Get index in global point list
const labelList& curAddr = sharedPointAddr_[meshI];
// Loop through all local points
forAll (curShared, i)
{
if (!gpSet[curAddr[i]])
{
// Point not set: set it
gp[curAddr[i]] = P[curShared[i]];
gpSet[curAddr[i]] = true;
}
else
{
// Point already set: check location
if (mag(gp[curAddr[i]] - P[curShared[i]]) > SMALL)
{
Info<< "MERGE MISMATCH: mesh" << meshI
<< " point: " << curShared[i]
<< " dist: " << gp[curAddr[i]] << " "
<< P[curShared[i]]
<< endl;
}
}
}
}
}
/*
// Grab marked points
OFstream ppp("points.vtk");
ppp << "# vtk DataFile Version 2.0" << nl
<< "points.vtk" << nl
<< "ASCII" << nl
<< "DATASET POLYDATA" << nl
<< "POINTS " << nGlobalPoints_ << " float" << nl;
Pout << "Global marking " << nGlobalPoints_ << endl;
forAll (gp, i)
{
ppp << float(gp[i].x()) << ' '
<< float(gp[i].y()) << ' '
<< float(gp[i].z())
<< nl;
Pout << gp[i] << endl;
}
*/
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sharedPoints::sharedPoints(const PtrList<fvMesh>& meshes)
Foam::sharedPoints::sharedPoints
(
const PtrList<fvMesh>& meshes
)
:
meshes_(meshes),
sharedPointAddr_(meshes_.size()),
@ -463,11 +614,19 @@ Foam::sharedPoints::sharedPoints(const PtrList<fvMesh>& meshes)
calcSharedPoints();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
// Foam::sharedPoints::~sharedPoints()
// {}
Foam::sharedPoints::sharedPoints
(
const PtrList<fvMesh>& meshes,
const labelListList& globalPointIndex
)
:
meshes_(meshes),
sharedPointAddr_(meshes_.size()),
sharedPointLabels_(meshes_.size()),
nGlobalPoints_(0)
{
calcSharedPoints(globalPointIndex);
}
// ************************************************************************* //

15
src/decompositionMethods/decomposeReconstruct/reconstructTools/finiteVolume/sharedPoints.H
View file

@ -117,13 +117,26 @@ class sharedPoints
//- Calculate shared points
void calcSharedPoints();
//- Calculate shared points
void calcSharedPoints(const labelListList& globalPointIndex);
public:
// Constructors
//- Construct from the list of meshes
sharedPoints(const PtrList<fvMesh>& meshes);
sharedPoints
(
const PtrList<fvMesh>& meshes
);
//- Construct from the list of meshes and global point addressing
sharedPoints
(
const PtrList<fvMesh>& meshes,
const labelListList& globalPointIndex
);
//- Destructor

14
src/dynamicMesh/topoChangerFvMesh/dynamicPolyRefinementFvMesh/dynamicPolyRefinementFvMesh.C
View file

@ -271,13 +271,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
(
refinementSelectionPtr_->refinementCellCandidates()()
);
Pout<< "Selected " << refCandidates.size()
Info<< "Selected " << refCandidates.size()
<< " refinement candidates."
<< endl;
}
else
{
Pout<< "Skipping refinement for this time-step..." << endl;
Info<< "Skipping refinement for this time-step..." << endl;
}
// Set cells to refine. Note: refinement needs to make sure that face
@ -297,13 +297,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
(
refinementSelectionPtr_->unrefinementPointCandidates()()
);
Pout<< "Selected " << unrefCandidates.size()
Info<< "Selected " << unrefCandidates.size()
<< " unrefinement candidates."
<< endl;
}
else
{
Pout<< "Skipping unrefinement for this time-step..." << endl;
Info<< "Skipping unrefinement for this time-step..." << endl;
}
// Set split points to unrefine around.
@ -344,13 +344,13 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
// some topo changes
if (sizeCellMap)
{
Pout<< "Successfully performed polyhedral refinement. "
Info<< "Successfully performed polyhedral refinement. "
<< "Changed from " << nOldCells << " to " << sizeCellMap
<< " cells." << endl;
}
else
{
Pout<< "Refinement/unrefinement not performed in this time step "
Info<< "Refinement/unrefinement not performed in this time step "
<< "since no cells were selected." << endl;
}
@ -363,7 +363,7 @@ bool Foam::dynamicPolyRefinementFvMesh::update()
curTimeIndex_ = time().timeIndex();
}
Pout<< "No refinement/unrefinement" << endl;
Info<< "No refinement/unrefinement" << endl;
// No refinement/unrefinement at this time step. Return false
return false;

58
src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalance.C
View file

@ -33,6 +33,7 @@ License
#include "addToRunTimeSelectionTable.H"
#include "cloud.H"
#include "cloudDistribute.H"
#include "meshObjectBase.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -138,10 +139,14 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
// Prepare receiving side
// Create the reconstructor
processorMeshesReconstructor meshRecon("reconstructed");
//
// HR 21.12.18 : Use empty domainname to avoid auto-created of
// fvSchemes/fvSolution
processorMeshesReconstructor meshRecon("");
PtrList<fvMesh>& procMeshes = meshRecon.meshes();
procMeshes.setSize(meshDecomp.nProcs());
meshRecon.globalPointIndex().setSize(meshDecomp.nProcs());
// Collect local fields for decomposition
clearOut();
@ -242,6 +247,29 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
}
}
// HR 14.12.18: Create dummy database pointing into the non-parallel case
// directory and disable the runTimeModifiable switch. dummyTime is used
// for decomposed, received and reconstructed fvMeshes (ie. before its data
// is moved into *this).
//
// The pointing into the non-parallel case directory is somewhat a hack to
// prevent auto-creation of fvSchemes and fvSolution (they exist in the root
// case). This is potentially dangerous if we write anything, but fvMeshes
// derived from this database are NO_WRITE so we should be fine. Making if
// non-runTimeModifiable prevents registration of fvSolution with the
// fileMonitor (addWatch). Not all processors will necessarily receive a
// mesh and watches will then cause dead-locks!
Time dummyTime
(
time().rootPath(),
time().globalCaseName(),
"system",
"constant",
false
);
const_cast<Switch&>(dummyTime.runTimeModifiable()) = false;
for (label procI = 0; procI < meshDecomp.nProcs(); procI++)
{
// Check if there is a mesh to send
@ -253,8 +281,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
autoPtr<fvMesh> procMeshPtr = meshDecomp.processorMesh
(
procI,
time(),
"processorPart" + Foam::name(procI),
dummyTime,
"", // HR 21.12.18 : Use empty domainname to avoid
// auto-created offvSchemes/fvSolution
true // Create passive processor patches
);
fvMesh& procMesh = procMeshPtr();
@ -286,6 +315,7 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
// Send the mesh and fields to target processor
toProc << procMesh << nl;
toProc << meshDecomp.globalPointIndex(procI) << nl;
// Send fields
sendFields(volScalarFields, fieldDecomposer, toProc);
@ -342,6 +372,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
procMeshPtr
);
meshRecon.globalPointIndex()[procI] =
meshDecomp.globalPointIndex(procI);
// Set local fields
// Note: first index is field index and second index is procI
insertFields
@ -411,6 +444,10 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
}
//HJ Insert clouds missing. HJ, 12/Oct/2018
//
// HR 18.11.18 - Not missing. Step is trivial and is treated in
// Cloud<ParticleType>::split which is called in the constructor
// of CloudDistibute
}
}
}
@ -445,9 +482,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
(
IOobject
(
"processorPart" + Foam::name(procI),
time().timeName(),
time(),
"", // "processorPart" + Foam::name(procI),
dummyTime.timeName(),
dummyTime,
IOobject::NO_READ,
IOobject::NO_WRITE
),
@ -456,6 +493,10 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
)
);
// Receive the global points addr
labelList gppi(fromProc);
meshRecon.globalPointIndex()[procI] = gppi;
// Receive fields
// Note: first index is field index and second index is procI
receiveFields
@ -535,7 +576,7 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
// Create the reconstructed mesh
autoPtr<fvMesh> reconstructedMeshPtr =
meshRecon.reconstructMesh(time());
meshRecon.reconstructMesh(dummyTime);
fvMesh& reconMesh = reconstructedMeshPtr();
Pout<< "Reconstructed mesh stats: "
@ -915,6 +956,9 @@ bool Foam::topoChangerFvMesh::loadBalance(const dictionary& decompDict)
);
}
// HR 13.12.18: Update the mesh objects
meshObjectBase::allUpdateTopology<polyMesh>(*this, meshMap);
// Debug: remove? HJ, 22/Oct/2018
// checkMesh(true);

24
src/dynamicMesh/topoChangerFvMesh/topoChangerFvMesh/topoChangerFvMeshLoadBalanceTemplates.C
View file

@ -126,7 +126,7 @@ void Foam::topoChangerFvMesh::receiveFields
receivedFields.setSize(nScalarFields);
fromProc.readBegin("topoChangerFvMeshReceiveFields");
forAll (receivedFields, fI)
{
word fieldName(fromProc);
@ -263,12 +263,12 @@ void Foam::topoChangerFvMesh::rebuildFields
// Resize internal field
if
(
masterField.internalField().size()
masterField.size()
!= GeoField::GeoMeshType::size(masterField.mesh())
)
{
Pout<< "Resizing internal field: old size = "
<< masterField.internalField().size()
<< masterField.size()
<< " new size = "
<< GeoField::GeoMeshType::size(masterField.mesh())
<< endl;
@ -280,7 +280,7 @@ void Foam::topoChangerFvMesh::rebuildFields
// Resize boundary (number of patches)
typename GeoField::GeometricBoundaryField& patchFields =
masterField.boundaryField();
masterField.boundaryFieldNoStoreOldTimes();
if (patchFields.size() != masterField.mesh().boundary().size())
{
@ -352,6 +352,22 @@ void Foam::topoChangerFvMesh::rebuildFields
fieldI++;
Pout<< "... done" << endl;
}
// HR 14.12.18: We create new processor boundary faces from internal
// faces. The values on these faces could be initialised by interpolation.
// Instead we choose to fix the values by evaluating the boundaries.
// I tried to execute evaluateCoupled() at the end of
// fvFieldReconstructor::reconstructField, but this fails in a strange way.
forAllConstIter
(
typename HashTable<const GeoField*>,
geoFields,
iter
)
{
GeoField& masterField = const_cast<GeoField&>(*iter());
masterField.boundaryField().evaluateCoupled();
}
}

21
src/finiteVolume/cfdTools/compressible/compressibleCourantNo.H
View file

@ -33,7 +33,12 @@ scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
scalar velMag = 0.0;
if (mesh.nInternalFaces())
// HR 26.06.18: A parallel run has at least two cells and therefore at least
// one internal face in the global mesh. It may be a processor boundary, but
// this is captured by max(mag(phi)).
// Old formulation hangs on parallel cases where one partition is degenerated
// to a single cell.
if (mesh.nInternalFaces() || Pstream::parRun())
{
surfaceScalarField phiOverRho = mag(phi)/fvc::interpolate(rho);
@ -47,6 +52,20 @@ if (mesh.nInternalFaces())
velMag = max(phiOverRho/mesh.magSf()).value();
}
else
{
// Single cell mesh: Co is still defined; use cell formulation
const scalar deltaT = runTime.deltaT().value();
const scalar deltaX = Foam::cbrt(mesh.V()[0]);
velMag = mag(U[0]);
CoNum = velMag*deltaT/deltaX;
meanCoNum = CoNum;
}
Info<< "Courant Number mean: " << meanCoNum
<< " max: " << CoNum

2
src/finiteVolume/cfdTools/incompressible/CourantNo.H
View file

@ -33,7 +33,7 @@ scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
scalar velMag = 0.0;
// HR 26.06.28: A parallel run has at least two cells and therefore at least
// HR 26.06.18: A parallel run has at least two cells and therefore at least
// one internal face in the global mesh. It may be a processor boundary, but
// this is captured by max(mag(phi)).
// Old formulation hangs on parallel cases where one partition is degenerated

8
src/finiteVolume/fields/fvPatchFields/basic/fixedGradient/fixedGradientFvPatchField.C
View file

@ -67,11 +67,11 @@ fixedGradientFvPatchField<Type>::fixedGradientFvPatchField
const dictionary& dict
)
:
fvPatchField<Type>(p, iF, dict),
// HR 15.12.18: Must not call evaluate during construction. Read value
// instead. This is needed for PLB.
fvPatchField<Type>(p, iF, dict, true),
gradient_("gradient", dict, p.size())
{
evaluate();
}
{}
template<class Type>

8
src/finiteVolume/fields/fvPatchFields/basic/mixed/mixedFvPatchField.C
View file

@ -54,13 +54,13 @@ mixedFvPatchField<Type>::mixedFvPatchField
const dictionary& dict
)
:
fvPatchField<Type>(p, iF, dict),
// HR 15.12.18: Must not call evaluate during construction. Read value
// instead. This is needed for PLB.
fvPatchField<Type>(p, iF, dict, true),
refValue_("refValue", dict, p.size()),
refGrad_("refGradient", dict, p.size()),
valueFraction_("valueFraction", dict, p.size())
{
evaluate();
}
{}
template<class Type>

2
src/finiteVolume/fields/fvPatchFields/derived/inletOutlet/inletOutletFvPatchField.C
View file

@ -62,7 +62,7 @@ inletOutletFvPatchField<Type>::inletOutletFvPatchField
{
// Read patch type
this->readPatchType(dict);
this->refValue() = Field<Type>("inletValue", dict, p.size());
if (dict.found("value"))

29
src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
View file

@ -70,7 +70,12 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
"fvSchemes",
obr.time().system(),
obr,
IOobject::READ_IF_PRESENT_IF_MODIFIED, // Allow default dictionary creation
(
obr.readOpt() == IOobject::MUST_READ
|| obr.readOpt() == IOobject::READ_IF_PRESENT
? IOobject::MUST_READ_IF_MODIFIED
: obr.readOpt()
),
IOobject::NO_WRITE
)
),
@ -168,21 +173,17 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
fluxRequired_(), // Do not read flux required option
defaultFluxRequired_(false)
{
if (!headerOk())
// HR 21.12.18 : Writing a default fvSchemes is not useful in PLB. Please
// specify MUST_READ on obr if you need this.
if
(
readOpt() == IOobject::MUST_READ
|| readOpt() == IOobject::MUST_READ_IF_MODIFIED
|| (readOpt() == IOobject::READ_IF_PRESENT && headerOk())
)
{
if (debug)
{
InfoIn
(
"fvSchemes::fvSchemes(const objectRegistry& obr)"
) << "fvSchemes dictionary not found. Creating default."
<< endl;
}
regIOobject::write();
read();
}
read();
}

3
src/finiteVolume/fvMesh/fvMesh.C
View file

@ -424,6 +424,9 @@ void Foam::fvMesh::resetFvPrimitives
// Reset fvPatches HJ, 16/Apr/2018
boundary_.resetFvPatches(resetFvPatchFlag);
// HR 14.12.18: Indicate that the mesh is changing to e.g. nearWallDist
polyMesh::changing(true);
// Clear all mesh data
clearOut();
}

2
src/finiteVolume/fvMesh/wallDist/nearWallDist.C
View file

@ -84,7 +84,7 @@ void Foam::nearWallDist::correct()
{
setSize(mesh_.boundary().size());
}
// Update size of GeometricBoundaryField
forAll (mesh_.boundary(), patchI)
{

1
src/foam/Make/files
View file

@ -447,6 +447,7 @@ $(globalMeshData)/globalMeshData.C
$(globalMeshData)/globalPoints.C
$(globalMeshData)/globalIndex.C
$(globalMeshData)/globalProcFaceIndex.C
$(globalMeshData)/globalProcPointIndex.C
$(polyMesh)/syncTools/syncTools.C

10
src/foam/fields/GeometricFields/GeometricField/GeometricField.C
View file

@ -728,6 +728,16 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryField()
}
// Return reference to GeometricBoundaryField
template<class Type, template<class> class PatchField, class GeoMesh>
typename
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField&
Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryFieldNoStoreOldTimes()
{
return boundaryField_;
}
// Store old-time field
template<class Type, template<class> class PatchField, class GeoMesh>
void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTimes() const

4
src/foam/fields/GeometricFields/GeometricField/GeometricField.H
View file

@ -447,6 +447,10 @@ public:
//- Return reference to GeometricBoundaryField for const field
inline const GeometricBoundaryField& boundaryField() const;
//- Return reference to GeometricBoundaryField without storing old
//- times
GeometricBoundaryField& boundaryFieldNoStoreOldTimes();
//- Return the time index of the field
inline label timeIndex() const;

7
src/foam/fields/cloud/cloud.C
View file

@ -64,6 +64,13 @@ Foam::autoPtr<Foam::cloudDistribute> Foam::cloud::cloudDist
}
Foam::labelList Foam::cloud::nParticlesPerCell() const
{
NotImplemented;
return labelList(0);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::cloud::~cloud()

4
src/foam/fields/cloud/cloud.H
View file

@ -101,7 +101,9 @@ public:
//- Return size of the cloud
virtual label size() const;
//- Count and return number of particles per cell
virtual labelList nParticlesPerCell() const;
// Edit
//- Remap the cells of particles corresponding to the

32
src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
View file

@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "pointBoundaryMesh.H"
#include "polyMesh.H"
#include "polyBoundaryMesh.H"
#include "facePointPatch.H"
#include "globalPointPatch.H"
@ -105,10 +106,39 @@ void Foam::pointBoundaryMesh::movePoints()
}
void Foam::pointBoundaryMesh::updateMesh()
void Foam::pointBoundaryMesh::updateMesh
(
const polyMesh& pMesh
)
{
const polyBoundaryMesh& pBoundary = pMesh.boundaryMesh();
pointPatchList& patches = *this;
// 21.1.19 HR : Patches need to be recreated in PLB
patches.clear();
patches.setSize(pBoundary.size());
forAll(patches, patchI)
{
patches.set
(
patchI,
facePointPatch::New(pBoundary[patchI], *this).ptr()
);
}
// Add the globalPointPatch
if(pMesh.globalData().nGlobalPoints())
{
patches.setSize(pBoundary.size() + 1);
patches.set
(
patches.size() - 1,
new globalPointPatch(*this, patches.size() - 1)
);
}
forAll(patches, patchi)
{
patches[patchi].initUpdateMesh();

6
src/foam/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H
View file

@ -46,6 +46,7 @@ namespace Foam
class pointMesh;
class polyBoundaryMesh;
class globalPointPatch;
class polyMesh;
/*---------------------------------------------------------------------------*\
Class pointBoundaryMesh Declaration
@ -105,7 +106,10 @@ public:
void movePoints();
//- Correct polyBoundaryMesh after topology update
void updateMesh();
void updateMesh
(
const polyMesh& pMesh
);
};

5
src/foam/meshes/pointMesh/pointMesh.C
View file

@ -122,7 +122,10 @@ bool Foam::pointMesh::updateMesh(const mapPolyMesh& mpm) const
{
// Casting const-ness to answer the interface of meshObject
// HJ, 30/Aug/2010
const_cast<pointBoundaryMesh&>(boundary_).updateMesh();
const_cast<pointBoundaryMesh&>(boundary_).updateMesh
(
GeoMesh<polyMesh>::mesh_
);
// Map all registered point fields
mapFields(mpm);

6
src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.C
View file

@ -118,7 +118,8 @@ void Foam::globalProcFaceIndex::calcFaceIndex()
OPstream toProc
(
Pstream::nonBlocking,
// HR 12.12.18: nonBlocking fails on PLB of Aachen bomb
Pstream::blocking,
procPatch.neighbProcNo()
);
toProc<< curFaceLabels;
@ -154,7 +155,8 @@ void Foam::globalProcFaceIndex::calcFaceIndex()
// Receive the data from master and insert into the list
IPstream fromProc
(
Pstream::nonBlocking,
// HR 12.12.18: nonBlocking fails on PLB of Aachen bomb
Pstream::blocking,
procPatch.neighbProcNo()
);

2
src/foam/meshes/polyMesh/globalMeshData/globalProcFaceIndex.H
View file

@ -120,7 +120,7 @@ public:
return procFaceOffset_;
}
//- Return lobal face label for all faces of current mesh
//- Return global face label for all faces of current mesh
// Sized to number of live faces in the mesh
const labelList globalLabel() const
{

319
src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.C Normal file
View file

@ -0,0 +1,319 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "globalProcPointIndex.H"
#include "processorPolyPatch.H"
#include "globalMeshData.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::globalProcPointIndex::calcPointIndex()
{
// Count and mark points in the following order
// 1. global points
// 2. processor patches
// 3. regular patches
// 4. anything else ("internal points")
//
// The marking is a follows
// -4 : Not visited yet
// -3 : Found in a regular patch
// -2 : Found in slave processor patch
// -1 : Found in a master processor patch
// 0 - nGlobalPoints-1 : Global point ID
// 1. Mark the global points
const labelList& spl = mesh_.globalData().sharedPointLabels();
const labelList& spa = mesh_.globalData().sharedPointAddr();
forAll(spl, spI)
{
globalLabel_[spl[spI]] = spa[spI];
}
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
labelList procPointCounter(Pstream::nProcs(), 0);
labelList patchPointCounter(patches.size(), 0);
label& nProcPoints = procPointCounter[Pstream::myProcNo()];
// Master processor patches - first pass
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& patchMeshPoints = pp.meshPoints();
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
label& nPatchPoints = patchPointCounter[patchI];
if (procPatch.master())
{
forAll (patchMeshPoints, mpI)
{
label& pl = globalLabel_[patchMeshPoints[mpI]];
if (pl < 0)
{
nProcPoints++;
nPatchPoints++;
pl = -1;
}
}
}
}
}
// Slave processor patches - first pass
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& patchMeshPoints = pp.meshPoints();
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
if (!procPatch.master())
{
forAll (patchMeshPoints, mpI)
{
label& pl = globalLabel_[patchMeshPoints[mpI]];
if (pl == -4)
{
pl = -2;
}
}
}
}
}
// Regular patches - first pass
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& patchMeshPoints = pp.meshPoints();
if (!isA<processorPolyPatch>(pp))
{
label& nPatchPoints = patchPointCounter[patchI];
forAll (patchMeshPoints, mpI)
{
label& pl = globalLabel_[patchMeshPoints[mpI]];
if (pl == -4)
{
nProcPoints++;
nPatchPoints++;
pl = -3;
}
}
}
}
// "Internal" points - first pass
label nInternalPoints = 0;
forAll (globalLabel_, pI)
{
if (globalLabel_[pI] == -4)
{
nInternalPoints++;
nProcPoints++;
}
}
// Gather-Scatter counter data
Pstream::gatherList(procPointCounter);
Pstream::scatterList(procPointCounter);
// Convert counters to offsets
procPointOffset_[0] = mesh_.globalData().nGlobalPoints();
for (label procI = 1; procI < procPointOffset_.size(); procI++)
{
procPointOffset_[procI] =
procPointCounter[procI-1] + procPointOffset_[procI-1];
}
patchPointOffset_[0] =
procPointOffset_[Pstream::myProcNo()] + nInternalPoints;
for (label patchI = 1; patchI < patchPointOffset_.size(); patchI++)
{
patchPointOffset_[patchI] =
patchPointCounter[patchI-1] + patchPointOffset_[patchI-1];
}
const pointField& points = mesh_.points();
// Master processor patches - second pass
// Send patch offset to slave and assign global labels to points on
// master processor patches and regular patches
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
const labelList& patchMeshPoints = pp.meshPoints();
const label patchPO = patchPointOffset_[patchI];
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
if (procPatch.master())
{
// Master processor patch
// Send the offset to the slave and mark the points through a
// face loop to establish an order that can be untangled on
// slave side
OPstream toProc
(
Pstream::nonBlocking,
procPatch.neighbProcNo()
);
toProc<< patchPO;
pointField pointLocs(patchMeshPoints.size());
label glPointLabelsI = 0;
forAll (pp, faceI)
{
const face& curFace = pp[faceI];
forAll (curFace, fpI)
{
const label pointI = curFace[fpI];
label& pl = globalLabel_[pointI];
if (pl == -1)
{
pl = patchPO + glPointLabelsI;
pointLocs[glPointLabelsI] = points[pointI];
glPointLabelsI++;
}
}
}
}
}
}
// Slave processor and regular patches - second pass
// Receive data on slave and assign global labels to points on
// master processor patches and regular patches
forAll (patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (isA<processorPolyPatch>(pp))
{
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(pp);
if (!procPatch.master())
{
// Slave processor patch - second pass
// Receive the offset from master and mark the points by taking
// into account that the points in the faces are in reverse
// order
IPstream fromProc
(
Pstream::nonBlocking,
procPatch.neighbProcNo()
);
label masterPatchPO(readLabel(fromProc));
label glPointLabelsI = 0;
forAll (pp, faceI)
{
const face curFace = pp[faceI].reverseFace();
forAll (curFace, fpI)
{
const label pointI = curFace[fpI];
label& pl = globalLabel_[pointI];
if (pl == -2)
{
pl = masterPatchPO + glPointLabelsI;
glPointLabelsI++;
}
}
}
}
}
else
{
// Regular patch - second pass
const labelList& patchMeshPoints = pp.meshPoints();
const label patchPO = patchPointOffset_[patchI];
label glPointLabelsI = 0;
forAll (patchMeshPoints, mpI)
{
label& pl = globalLabel_[patchMeshPoints[mpI]];
if (pl == -3)
{
pl = patchPO + glPointLabelsI;
glPointLabelsI++;
}
}
}
}
// "Internal" points - second pass
nInternalPoints = 0;
forAll (globalLabel_, pI)
{
label& pl = globalLabel_[pI];
if (pl == -4)
{
pl = procPointOffset_[Pstream::myProcNo()] + nInternalPoints;
nInternalPoints++;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::globalProcPointIndex::globalProcPointIndex(const polyMesh& mesh)
:
mesh_(mesh),
procPointOffset_(Pstream::nProcs()+1),
patchPointOffset_(mesh_.boundaryMesh().size()+1),
globalLabel_(mesh_.nPoints(), -4)
{
calcPointIndex();
}
// ************************************************************************* //

135
src/foam/meshes/polyMesh/globalMeshData/globalProcPointIndex.H Normal file
View file

@ -0,0 +1,135 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::globalProcPointIndex
Description
The class creates a unique global point index for each processor points
(pair) in the mesh. Master and slave processor points carry the same index.
Point offsets counts a number of unique points on each processor,
excluding slave processor patch points and global points.
Currently, faces are ordered with internal faces first, followed by patch
faces in patch order, excluding slave processor patches.
If needed, this can be changed to group processor faces with internal faces
to facilitate parallel I/O. HJ, 4/May/2018
Author
Henrik Rusche, Wikki GmbH
SourceFiles
globalProcPointIndex.C
\*---------------------------------------------------------------------------*/
#ifndef globalProcPointIndex_H
#define globalProcPointIndex_H
#include "polyMesh.H"
#include "labelList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class globalProcPointIndex Declaration
\*---------------------------------------------------------------------------*/
class globalProcPointIndex
{
// Private data
//- Mesh reference
const polyMesh& mesh_;
//- Processor point index offset
labelList procPointOffset_;
//- point index offset for patches in local mesh
labelList patchPointOffset_;
//- Global poins label for all points of current mesh
// Sized to number of live points in the mesh
labelList globalLabel_;
// Private Member Functions
//- Disallow default bitwise copy construct
globalProcPointIndex(const globalProcPointIndex&);
//- Disallow default bitwise assignment
void operator=(const globalProcPointIndex&);
//- Calculate point index
void calcPointIndex();
public:
// Constructors
//- Construct from mesh
globalProcPointIndex(const polyMesh&);
//- Destructor - default
// Member Functions
//- Return point index offset per processor
inline const labelList& procPointOffset() const
{
return procPointOffset_;
}
inline const labelList& patchPointOffset() const
{
return patchPointOffset_;
}
//- Return global point label for all points of current mesh
// Sized to number of live points in the mesh
const labelList globalLabel() const
{
return globalLabel_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/immersedBoundary/immersedBoundaryDynamicMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh/immersedBoundaryDynamicRefineSolidBodyMotionFvMesh.C
View file

@ -127,7 +127,7 @@ bool immersedBoundaryDynamicRefineSolidBodyMotionFvMesh::update()
{
hasChanged = loadBalance(refinementDict());
}
// Execute dummy mesh motion for the background mesh
const pointField oldPoints = allPoints();
fvMesh::movePoints(oldPoints);

2
src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryEpsilonWallFunction/immersedBoundaryEpsilonWallFunctionFvPatchScalarField.C
View file

@ -68,7 +68,7 @@ immersedBoundaryEpsilonWallFunctionFvPatchScalarField
p,
Switch(dict.lookup("setDeadValue")),
readScalar(dict.lookup("deadValue"))
)
)
{}

2
src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryKqRWallFunction/immersedBoundaryKqRWallFunctionFvPatchField.C
View file

@ -150,7 +150,7 @@ immersedBoundaryKqRWallFunctionFvPatchField
(
ptf.ibPatch(),
ptf.setDeadValue(),
ptf.deadValue()
ptf.deadValue()
)
{
this->checkType();

2
src/immersedBoundary/immersedBoundaryTurbulence/wallFunctions/immersedBoundaryOmegaWallFunctions/immersedBoundaryOmegaWallFunctionFvPatchScalarField.C
View file

@ -66,7 +66,7 @@ immersedBoundaryOmegaWallFunctionFvPatchScalarField
p,
Switch(dict.lookup("setDeadValue")),
readScalar(dict.lookup("deadValue"))
)
)
{}

2
src/lagrangian/basic/CloudDistributeTemplate/CloudDistributeTemplate.H
View file

@ -47,7 +47,7 @@ class Cloud;
class polyMesh;
class cloud;
/*---------------------------------------------------------------------------*\
Class Cloud Declaration
\*---------------------------------------------------------------------------*/

34
src/lagrangian/basic/CloudTemplate/CloudTemplate.C
View file

@ -506,6 +506,40 @@ void Foam::Cloud<ParticleType>::rebuild
}
template<class ParticleType>
Foam::labelList Foam::Cloud<ParticleType>::nParticlesPerCell() const
{
Pout << "In nParticlesPerCell" << endl;
labelList nppc (pMesh().nCells(), 0);
forAllConstIter(typename Cloud<ParticleType>, *this, pIter)
{
const ParticleType& p = pIter();
const label celli = p.cell();
// Check
if (celli < 0 || celli >= pMesh().nCells())
{
FatalErrorIn
(
"Foam::Cloud<ParticleType>::nParticlesPerCell()"
)
<< "Illegal cell number " << celli
<< " at position " << p.position() << nl
<< "Cell number should be between 0 and "
<< pMesh().nCells()-1 << nl
<< "On this mesh the particle should be in cell "
<< pMesh().findCell(p.position())
<< exit(FatalError);
}
nppc[celli]++;
}
return nppc;
}
template<class ParticleType>
void Foam::Cloud<ParticleType>::writePositions() const
{

3
src/lagrangian/basic/CloudTemplate/CloudTemplate.H
View file

@ -290,6 +290,9 @@ public:
const PtrList<labelIOList>& faceProcAddressing
);
//- Count and return number of particles per cell
virtual labelList nParticlesPerCell() const;
// Read

2
src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C
View file

@ -370,7 +370,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
// constrainToMeshCentre is inside of an if-statement
mesh.geometricD();
mesh.bounds();
// Prepare for next time step
label parcelsAdded = 0;
scalar massAdded = 0.0;

3
src/overset/oversetMesh/donorAcceptor/donorAcceptor.H
View file

@ -134,7 +134,8 @@ public:
donorProcNo_(-1),
donorPoint_(vector::zero),
extendedDonorCells_(),
extendedDonorPoints_()
extendedDonorPoints_(),
withinBB_(false)
{}
//- Construct from acceptor data

125
src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C
View file

@ -67,54 +67,44 @@ void Foam::adaptiveOverlapFringe::suitabilityFractionSlope
scalar& beta
) const
{
// Linear regression coefficients equations implemented down there:
// beta = (n*sum(x_i*y_i)-sum(x_i)*sum(y_i)) / (n*sum(x_i^2)-(sum(x_i))^2)
// alpha = sum(y_i)/n - beta* (sum(x_i)/n
// Linear regression coefficients y = alpha + beta*x
// beta = sum((x_i - x_mean)*(y_i - y_mean))/sum(x_i - x_mean)^2
// alpha = y_mean - beta*x_mean
// x -> iteration number
// y -> suitability
// Iteration ordinal numbers sum (iteration ordinal number represents x
// value in x-y coordinate system)
scalar iterSum = 0;
// Calculate mean values first
scalar iterMean = 0;
scalar suitabilityMean = 0;
// Sum of squared iteration ordinal numbers
scalar iterSquaredSum = 0;
// Average donor/acceptor suitability fraction values sum (average
// suitability represents y value in x-y coordinate system)
scalar suitabilitySum = 0;
// Sum of iteration ordinal number multiplied by donor/acceptor suitability
// fraction value from one iteration
scalar iterXSuitabilitySum = 0;
// Loop through all iteration data objects and calculate values needed for
// linear regression calculation
forAllConstIter(FIFOStack<iterationData>, iterHist, it)
{
iterSum += it().iteration();
iterSquaredSum += pow(it().iteration(), 2);
suitabilitySum += it().suitability();
iterXSuitabilitySum += it().iteration()*it().suitability();
iterMean += it().iteration();
suitabilityMean += it().suitability();
}
// Number of iterations is equal to iterationDataHistory_ size
const label n = iterationDataHistory_.size();
reduce(iterMean, sumOp<scalar>());
iterMean /= iterHist.size();
reduce(suitabilityMean, sumOp<scalar>());
suitabilityMean /= iterHist.size();
// Calculate slope (i.e. if y = alpha + beta*x, slope is beta)
beta =
(n*iterXSuitabilitySum - iterSum*suitabilitySum)/
(n*iterSquaredSum - pow(iterSum, 2));
// Calculate denominator and numerator for beta
scalar n = 0;
scalar d = 0;
// Calculate mean iteration number - actually we don't need to calculate
// it beacuse we are only interested in trend, i.e we need only slope
// (beta). Maybe I will need it for diagrams in my master thesis, so
// calculate it for now.
const scalar meanIter = iterSum/n;
forAllConstIter(FIFOStack<iterationData>, iterHist, it)
{
n += (it().iteration() - iterMean)*
(it().suitability() - suitabilityMean);
d += Foam::sqr((it().iteration() - iterMean));
}
// Calculate mean suitability from all iterations
const scalar meanSuitability = suitabilitySum/n;
reduce(n, sumOp<scalar>());
reduce(d, sumOp<scalar>());
// Calculate intercept
alpha = meanSuitability - beta*meanIter;
// Calculate regression coefficients
beta = n/d;
alpha = suitabilityMean - beta*iterMean;
}
@ -179,20 +169,30 @@ void Foam::adaptiveOverlapFringe::calcAddressing() const
// holes)
boolList eligibleAcceptors(mesh.nCells(), true);
// Read user specified holes into allHoles list. Note: if the cell set
// is not found, the list will be empty
cellSet allHoles
(
mesh,
holesSetName_,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
);
// Read user specified holes into allHoles list. Note: use cellZone rather
// than cellSet to have correct behaviour on dynamic mesh simulations
// We will silently proceed if the zone is not found since this option is
// not mandatory but is useful in certain cases
// Get zone index
const label zoneID = mesh.cellZones().findZoneID(holesZoneName_);
// Create a hash set for allHoles
labelHashSet allHoles;
if (zoneID > -1)
{
// Get the zone for holes and append them to set
const labelList& specifiedHoles = mesh.cellZones()[zoneID];
allHoles.insert(specifiedHoles);
}
// else silently proceed without user-specified holes
// Extend allHoles with cutHoles
forAll (cutHoles, chI)
{
// Note: cellSet is a hashSet so it automatically removes duplicates
// Note: duplicated are removed because we're using hash set
allHoles.insert(cutHoles[chI]);
}
@ -393,7 +393,7 @@ Foam::adaptiveOverlapFringe::adaptiveOverlapFringe
acceptorsPtr_(NULL),
finalDonorAcceptorsPtr_(NULL),
holesSetName_(dict.lookupOrDefault<word>("holes", word())),
holesZoneName_(dict.lookupOrDefault<word>("holes", word())),
initPatchNames_
(
dict.lookupOrDefault<wordList>("initPatchNames", wordList())
@ -408,6 +408,14 @@ Foam::adaptiveOverlapFringe::adaptiveOverlapFringe
(
dict.lookupOrDefault<label>("specifiedIterationsNumber", 5)
),
minSuitabilityRate_
(
dict.lookupOrDefault<scalar>("minSuitabilityRate", 0.0)
),
maxIter_
(
dict.lookupOrDefault<label>("maximumIterations", 100)
),
relativeCounter_(specifiedIterationsNumber_),
additionalIterations_
(
@ -667,6 +675,7 @@ bool Foam::adaptiveOverlapFringe::updateIteration
(additionalIterations_)
&& (fringeIter_ == relativeCounter_)
&& (specifiedIterationsNumber_ > 1)
&& (fringeIter_ < maxIter_)
)
{
// Slope (i.e. if y = alpha + beta*x, slope is beta)
@ -695,7 +704,7 @@ bool Foam::adaptiveOverlapFringe::updateIteration
// If the donor/acceptor suitability gradient is positive,
// make 1 additional iteration
if (beta > 0)
if (beta > minSuitabilityRate_)
{
// Pop the bottom element of the stack
iterationDataHistory_.pop();
@ -710,7 +719,7 @@ bool Foam::adaptiveOverlapFringe::updateIteration
<< " will be performed." << endl;
}
}
}
}
// Go through unsuitable donor/acceptor pairs from previous iteration and
// find a new batch of acceptors and holes for the next iteration if:
@ -1076,16 +1085,26 @@ bool Foam::adaptiveOverlapFringe::updateIteration
reduce(nAccToHoles, sumOp<label>());
Info<< "Converted " << nAccToHoles << " acceptors to holes." << endl;
// Bugfix: Although we have found suitable overlap, we need to update
// acceptors as well because eligible donors for acceptors of other
// regions are calculated based on these acceptors (and holes)
labelList& acceptors = *acceptorsPtr_;
acceptors.setSize(finalDAPairs.size());
forAll (acceptors, aI)
{
acceptors[aI] = finalDAPairs[aI].acceptorCell();
}
// Construct final donor/acceptor list
finalDonorAcceptorsPtr_ = new donorAcceptorList
(
finalDAPairs
finalDAPairs.xfer()
);
// Construct final fringe holes list
fringeHolesPtr_ = new labelList
(
fringeHoles
fringeHoles.xfer()
);
// Print information

10
src/overset/oversetMesh/oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.H
View file

@ -105,7 +105,7 @@ class adaptiveOverlapFringe
//- Name of the cell set defining initial holes (empty by default).
// Useful when the resolution of the background mesh is much
// coarser than the front mesh and no hole is found
const word holesSetName_;
const word holesZoneName_;
//- Optional list of patches to start the iterative fringe assembly
// process (empty list by default). Useful when we actually have
@ -125,6 +125,14 @@ class adaptiveOverlapFringe
// be made
const label specifiedIterationsNumber_;
//- User specified suitability trend rate. If the trend is larger
// than the specified trend rate, additional iterations are
// performed
const scalar minSuitabilityRate_;
//- Maximum number of iterations
const label maxIter_;
//- Relative iteration counter
mutable label relativeCounter_;

15
src/overset/oversetMesh/oversetFringe/overlapFringe/donorSuitability/donorSuitability/donorSuitability.H
View file

@ -224,6 +224,21 @@ public:
<< " on processor: " << daPair.acceptorProcNo()
<< " did not find donor candidate."
<< nl
<< "Additional information: "
<< nl
<< (
daPair.withinBB()
? " Within bounding box"
: " Not within bounding box"
)
<< nl << tab
<< "Donor index: " << daPair.donorCell()
<< nl << tab
<< "Donor processor index: " << daPair.donorProcNo()
<< nl << tab
<< "Number of extended donors: "
<< daPair.extendedDonorCells().size()
<< nl << nl
<< "Please review your fringe assembly settings"
<< " (or try using adaptiveOverlap fringe algorithm)."
<< abort(FatalError);

63
src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.C
View file

@ -188,26 +188,37 @@ void Foam::overlapFringe::calcAddressing() const
// holes)
boolList eligibleAcceptors(mesh.nCells(), true);
// Read user specified holes into allHoles list. Note: if the cell set
// is not found, the list will be empty
labelList allHoles
(
cellSet
(
mesh,
holesSetName_,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
).toc()
);
// Read user specified holes into allHoles list. Note: use cellZone rather
// than cellSet to have correct behaviour on dynamic mesh simulations
// We will silently proceed if the zone is not found since this option is
// not mandatory but is useful in certain cases
// Get zone index
const label zoneID = mesh.cellZones().findZoneID(holesZoneName_);
// Create a hash set for allHoles
labelHashSet allHoles;
if (zoneID > -1)
{
// Get the zone for holes and append them to set
const labelList& specifiedHoles = mesh.cellZones()[zoneID];
allHoles.insert(specifiedHoles);
}
// else silently proceed without user-specified holes
// Extend allHoles with cutHoles
allHoles.append(cutHoles);
forAll (cutHoles, chI)
{
// Note: duplicated are removed because we're using hash set
allHoles.insert(cutHoles[chI]);
}
// Mark all holes
forAll (allHoles, hI)
forAllConstIter (labelHashSet, allHoles, iter)
{
const label& holeCellI = allHoles[hI];
const label& holeCellI = iter.key();
// Mask eligible acceptors
eligibleAcceptors[holeCellI] = false;
@ -225,10 +236,10 @@ void Foam::overlapFringe::calcAddressing() const
dynamicLabelList candidateAcceptors(mesh.nCells());
// Loop through all holes and find acceptor candidates
forAll (allHoles, hI)
forAllConstIter (labelHashSet, allHoles, iter)
{
// Get neighbours of this hole cell
const labelList& hNbrs = cc[allHoles[hI]];
const labelList& hNbrs = cc[iter.key()];
// Loop through neighbours of this hole cell
forAll (hNbrs, nbrI)
@ -370,7 +381,7 @@ void Foam::overlapFringe::calcAddressing() const
// Transfer the acceptor list and allocate empty fringeHoles list, which
// may be populated in updateIteration member function
acceptorsPtr_ = new labelList(candidateAcceptors.xfer());
fringeHolesPtr_ = new labelList(cutHoles);
fringeHolesPtr_ = new labelList(allHoles.toc().xfer());
}
@ -398,7 +409,7 @@ Foam::overlapFringe::overlapFringe
acceptorsPtr_(NULL),
finalDonorAcceptorsPtr_(NULL),
holesSetName_(dict.lookupOrDefault<word>("holes", word())),
holesZoneName_(dict.lookupOrDefault<word>("holes", word())),
initPatchNames_
(
dict.lookupOrDefault<wordList>("initPatchNames", wordList())
@ -687,8 +698,15 @@ bool Foam::overlapFringe::updateIteration
}
}
// Tranfer back the allFringeHoles dynamic list into member data
fringeHolesPtr_->transfer(allFringeHoles);
// Bugfix: Although we have found suitable overlap, we need to update
// acceptors as well because eligible donors for acceptors of other
// regions are calculated based on these acceptors (and holes)
labelList& acceptors = *acceptorsPtr_;
acceptors.setSize(finalDAPairs.size());
forAll (acceptors, aI)
{
acceptors[aI] = finalDAPairs[aI].acceptorCell();
}
// Transfer ownership of the final donor/acceptor list to the
// finalDonorAcceptorsPtr_
@ -697,6 +715,9 @@ bool Foam::overlapFringe::updateIteration
finalDAPairs.xfer()
);
// Tranfer back the allFringeHoles dynamic list into member data
fringeHolesPtr_->transfer(allFringeHoles);
// At least 100*minGlobalFraction_ % of suitable donor/acceptor pairs
// have been found.
Info<< "Converted " << nAccToHoles << " acceptors to holes."

2
src/overset/oversetMesh/oversetFringe/overlapFringe/overlapFringe/overlapFringe.H
View file

@ -96,7 +96,7 @@ class overlapFringe
//- Name of the cell set defining initial holes (empty by default).
// Useful when the resolution of the background mesh is much
// coarser than the front mesh and no hole is found
const word holesSetName_;
const word holesZoneName_;
//- Optional list of patches to start the iterative fringe assembly
// process (empty list by default). Useful when we actually have

21
src/overset/oversetMesh/oversetMesh/oversetMeshAddressing.C
View file

@ -55,18 +55,33 @@ void Foam::oversetMesh::calcCellClassification() const
nHoleCells += regions_[regionI].holes().size();
}
Pout<< "Number of acceptor cells: " << nAcceptorCells << endl;
acceptorCellsPtr_ = new labelList(nAcceptorCells);
labelList& acceptor = *acceptorCellsPtr_;
Pout<< "Number of donor cells: " << nDonorCells << endl;
donorCellsPtr_ = new labelList(nDonorCells);
labelList& donor = *donorCellsPtr_;
Pout<< "Number of hole cells: " << nHoleCells << endl;
holeCellsPtr_ = new labelList(nHoleCells);
labelList& hole = *holeCellsPtr_;
// Print out processor specific information in debug
if (oversetMesh::debug)
{
Pout<< "Number of acceptor cells: " << nAcceptorCells << endl;
Pout<< "Number of donor cells: " << nDonorCells << endl;
Pout<< "Number of hole cells: " << nHoleCells << endl;
}
// Else print global information
else
{
Info<< "Number of acceptor cells: "
<< returnReduce(nAcceptorCells, sumOp<label>()) << endl;
Info<< "Number of donor cells: "
<< returnReduce(nDonorCells, sumOp<label>()) << endl;
Info<< "Number of hole cells: "
<< returnReduce(nHoleCells, sumOp<label>()) << endl;
}
// Reset counters
nAcceptorCells = 0;
nDonorCells = 0;

92
src/overset/oversetMesh/oversetRegion/oversetRegion.C
View file

@ -190,18 +190,9 @@ void Foam::oversetRegion::calcDonorAcceptorCells() const
// Update global flag
foundGlobalOverlap &= regionFoundSuitableOverlap;
// If the overlap has not been found for this region, we need to
// reset:
// - holeCells (depend on fringe holes)
// - eligibleDonors (depend on fringe holes and acceptors),
// - cellSearch (depends on eligible donors).
if (!regionFoundSuitableOverlap)
{
deleteDemandDrivenData(curRegion.cellSearchPtr_);
}
deleteDemandDrivenData(curRegion.holeCellsPtr_);
deleteDemandDrivenData(curRegion.eligibleDonorCellsPtr_);
// Deletion of demand driven data relocated to the end of
// updateDonorAcceptors() member function. Makes more sense to have
// it there. VV, 13/Jan/2019.
}
} while (!foundGlobalOverlap);
@ -1317,6 +1308,20 @@ bool Foam::oversetRegion::updateDonorAcceptors() const
// STAGE 11: Filter possibly multiple remote donors
// Sanity check before moving on. For each initial acceptor, we must
// receive at least 1 corresponding donor (even if it is not found,
// indicating an orphan cell) from different processors.
if (a.size() > completeDonorAcceptorList.size())
{
FatalErrorIn("void oversetRegion::updateDonorAcceptors() const")
<< "Size of initial acceptor set: " << a.size()
<< " is larger than the size of the distributed "
<< " donor/acceptor list: " << completeDonorAcceptorList.size()
<< nl
<< "This should not have happened..."
<< abort(FatalError);
}
// Create a masking field indicating that a certain acceptor has been
// visited
boolList isVisited(a.size(), false);
@ -1379,35 +1384,16 @@ bool Foam::oversetRegion::updateDonorAcceptors() const
curDA.donorPoint(),
curDA.withinBB()
);
// Bugfix: also need to reset extended donors since a better
// candidate has been found. VV, 1/Jan/2019
curDACombined.setExtendedDonors(curDA);
}
}
}
// Update withinBB flag if the donor is within bounding box of acceptor
// (previously we checked whether the acceptor is within bounding box of
// donor)
forAll (combinedDonorAcceptorList, daI)
{
donorAcceptor& curDA = combinedDonorAcceptorList[daI];
// If the acceptor is not within bounding box of donor, set the flag
// other way around
if (!curDA.withinBB())
{
curDA.setWithinBB
(
mesh_.pointInCellBB
(
curDA.donorPoint(),
curDA.acceptorCell()
)
);
}
}
// Check whether all acceptors have been visited. Used for testing/debugging
// parallel comms
if (oversetMesh::debug)
// Check whether all acceptors have been visited. Useful check if in
// no-debug mode
{
bool allVisited = true;
@ -1435,12 +1421,35 @@ bool Foam::oversetRegion::updateDonorAcceptors() const
<< nl
<< "Try switching off useLocalBoundingBoxes for all regions"
<< nl
<< "(this optimisation is switched on by default)."
<< "(this optimisation is switched off by default)."
<< abort(FatalError);
}
}
}
// Update withinBB flag if the donor is within bounding box of acceptor
// (previously we checked whether the acceptor is within bounding box of
// donor)
forAll (combinedDonorAcceptorList, daI)
{
donorAcceptor& curDA = combinedDonorAcceptorList[daI];
// If the acceptor is not within bounding box of donor, set the flag
// other way around
if (!curDA.withinBB())
{
curDA.setWithinBB
(
mesh_.pointInCellBB
(
curDA.donorPoint(),
curDA.acceptorCell()
)
);
}
}
// STAGE 12: Finish the iteration by updating the fringe, which will
// actually hold final and some intermediate steps for donor/acceptor
// assembly process
@ -1450,6 +1459,13 @@ bool Foam::oversetRegion::updateDonorAcceptors() const
bool suitableOverlapFound =
fringePtr_->updateIteration(combinedDonorAcceptorList);
// Delete all necessary demand driven data for this region since we have
// just updated the iteration. Therefore, cellSearch, eligibleDonors and
// holes need to be updated
deleteDemandDrivenData(cellSearchPtr_);
deleteDemandDrivenData(eligibleDonorCellsPtr_);
deleteDemandDrivenData(holeCellsPtr_);
return suitableOverlapFound;
}

4
src/sampling/sampledSet/curve/curveSet.C
View file

@ -65,7 +65,7 @@ bool Foam::curveSet::trackToBoundary
const vector offset =
sampleCoords_[sampleI + 1] - sampleCoords_[sampleI];
const scalar smallDist = mag(tol*offset);
const scalar smallDist = mag(tol_()*offset);
point oldPos = trackPt;
label facei = -1;
@ -174,7 +174,7 @@ void Foam::curveSet::calcSamples
{
const vector offset =
sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
const scalar smallDist = mag(tol*offset);
const scalar smallDist = mag(tol_()*offset);
// Get all boundary intersections

16
src/sampling/sampledSet/face/faceOnlySet.C
View file

@ -43,6 +43,16 @@ namespace Foam
}
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::debug::optimisationSwitch
Foam::faceOnlySet::maxNSteps_
(
"maximumFaceOnlySetStepNumber",
100
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -58,7 +68,7 @@ bool Foam::faceOnlySet::trackToBoundary
{
// distance vector between sampling points
const vector offset = end_ - start_;
const vector smallVec = tol*offset;
const vector smallVec = tol_()*offset;
const scalar smallDist = mag(smallVec);
// Alias
@ -115,7 +125,7 @@ void Foam::faceOnlySet::calcSamples
const vector offset = (end_ - start_);
const vector normOffset = offset/mag(offset);
const vector smallVec = tol*offset;
const vector smallVec = tol_()*offset;
const scalar smallDist = mag(smallVec);
@ -198,7 +208,7 @@ void Foam::faceOnlySet::calcSamples
// index in bHits; current boundary intersection
label bHitI = 1;
while(true)
while(segmentI <= maxNSteps_())
{
if (trackFaceI != -1)
{

7
src/sampling/sampledSet/face/faceOnlySet.H
View file

@ -63,6 +63,13 @@ class faceOnlySet
point end_;
// Static data
//- Maximum number of steps for intersection search. Introduced to avoid
// inifinite loops (IG 4/Dec/2018)
static const debug::optimisationSwitch maxNSteps_;
// Private Member Functions
//- Samples from startTrackPt/CellI. Updates particle/samplePt/sampleI

24
src/sampling/sampledSet/sampledSet/sampledSet.C
View file

@ -33,13 +33,29 @@ License
namespace Foam
{
const scalar sampledSet::tol = 1e-6;
defineTypeNameAndDebug(sampledSet, 0);
defineRunTimeSelectionTable(sampledSet, word);
}
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::debug::tolerancesSwitch
Foam::sampledSet::tol_
(
"sampleSetTolerance",
1e-5
);
const Foam::debug::tolerancesSwitch
Foam::sampledSet::pushTol_
(
"sampleSetPushTolerance",
0.1
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::sampledSet::getBoundaryCell(const label faceI) const
@ -189,7 +205,7 @@ Foam::point Foam::sampledSet::pushIn
const point& cellCtr = mesh().cellCentres()[cellI];
point newSample =
facePt + tol*(cellCtr - facePt);
facePt + pushTol_()*(cellCtr - facePt);
if (!searchEngine().pointInCell(newSample, cellI))
{
@ -221,7 +237,7 @@ bool Foam::sampledSet::getTrackingPoint
label& trackFaceI
) const
{
const scalar smallDist = mag(tol*offset);
const scalar smallDist = mag(tol_()*offset);
bool isGoodSample = false;

8
src/sampling/sampledSet/sampledSet/sampledSet.H
View file

@ -48,6 +48,8 @@ SourceFiles
#include "runTimeSelectionTables.H"
#include "autoPtr.H"
#include "coordSet.H"
#include "debug.H"
#include "tolerancesSwitch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -201,7 +203,11 @@ public:
//- Tolerance when comparing points. Usually relative to difference
// between start_ and end_
const static scalar tol;
static const debug::tolerancesSwitch tol_;
//- Tolerance used to push points out of the face. It is relative to the
// distance between the cell face centre and cell centre.
static const debug::tolerancesSwitch pushTol_;
// Constructors

4
src/sampling/sampledSet/uniform/uniformSet.C
View file

@ -96,7 +96,7 @@ bool Foam::uniformSet::trackToBoundary
{
// distance vector between sampling points
const vector offset = (end_ - start_)/(nPoints_ - 1);
const vector smallVec = tol*offset;
const vector smallVec = tol_()*offset;
const scalar smallDist = mag(smallVec);
// Alias
@ -233,7 +233,7 @@ void Foam::uniformSet::calcSamples
const vector offset = (end_ - start_)/(nPoints_ - 1);
const vector normOffset = offset/mag(offset);
const vector smallVec = tol*offset;
const vector smallVec = tol_()*offset;
const scalar smallDist = mag(smallVec);
// Get all boundary intersections

17
src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.C
View file

@ -57,8 +57,21 @@ Foam::basicChemistryModel::basicChemistryModel
chemistry_(lookup("chemistry")),
deltaTChem_
(
mesh.nCells(),
readScalar(lookup("initialChemicalTimeStep"))
IOobject
(
"deltaTChem",
mesh.time().constant(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar
(
"deltaTChem0",
dimTime,
readScalar(lookup("initialChemicalTimeStep"))
)
)
{}

6
src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H
View file

@ -38,7 +38,7 @@ SourceFiles
#include "IOdictionary.H"
#include "Switch.H"
#include "scalarField.H"
#include "volFields.H"
#include "volFieldsFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -77,7 +77,9 @@ protected:
Switch chemistry_;
//- Latest estimation of integration step
scalarField deltaTChem_;
// HR 14.12.18: Upgraded to volScalarField in order to enable automatic
// mapping/resizing after a topo change. Also needed for PLB.
volScalarField deltaTChem_;
// Protected member functions

2
src/turbulenceModels/incompressible/RAS/RASModel/RASModel.H
View file

@ -298,7 +298,7 @@ public:
) const;
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct() = 0;
virtual void correct();
//- Read RASProperties dictionary
virtual bool read() = 0;

2
src/turbulenceModels/incompressible/turbulenceModel/turbulenceModel.H
View file

@ -207,7 +207,7 @@ public:
virtual tmp<fvVectorMatrix> divDevReff() const = 0;
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct() = 0;
virtual void correct();
//- Read turbulenceProperties dictionary
virtual bool read() = 0;