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

Immersed Boundary mapDistribute comms update

Browse source
This commit is contained in:
Hrvoje Jasak 2017-11-09 13:24:16 +00:00
parent 9fe436ac41
commit 51a5c5e0ab
6 changed files with 690 additions and 497 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

813
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.C
View file

File diff suppressed because it is too large Load diff

96
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatch.H
View file

@ -27,9 +27,15 @@ Class
Description
Immersed boundary FV patch
Notes
Rewrite of the immersed boundary method for increased efficiency
and robustness. Main changes
- using octree for fast search of eligible cells using indexedOctree
- parallel communication using mapDistribute
Author
Zeljko Tukovic
Reorganisation by Hrvoje Jasak
Reorganisation and optimisation by Hrvoje Jasak
SourceFiles
immersedBoundaryFvPatch.C
@ -48,11 +54,13 @@ SourceFiles
#include "volFieldsFwd.H"
#include "surfaceFieldsFwd.H"
#include "dynamicLabelList.H"
#include "labelPair.H"
#include "scalarList.H"
#include "vectorList.H"
#include "FieldFields.H"
#include "scalarMatrices.H"
#include "indexedOctree.H"
#include "treeDataCell.H"
#include "mapDistribute.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -84,14 +92,14 @@ class immersedBoundaryFvPatch
// Static data
//- Fitting angle rejection factor (deg)
// Cells within the the radius will be used for the fitting function
static const debug::tolerancesSwitch angleFactor_;
//- Fitting radius factor
// Cells within the the radius will be used for the fitting function
// Cells within the radius will be used for the fitting function
static const debug::tolerancesSwitch radiusFactor_;
//- Fitting angle rejection factor (deg)
// Cells outside of the angle will not be used for the fitting function
static const debug::tolerancesSwitch angleFactor_;
//- Maximum number of rows in cell-cell search
static const debug::optimisationSwitch maxCellCellRows_;
@ -166,17 +174,20 @@ class immersedBoundaryFvPatch
// (extended stencil)
mutable labelListList* ibCellCellsPtr_;
//- List of cells needed by neighbour processors
//- List of cells needed by other processors
mutable labelListList* ibProcCellsPtr_;
//- Centres of cells from neighbour processors
mutable vectorListList* ibProcCentresPtr_;
//- Processor mapDistribute
mutable mapDistribute* ibMapPtr_;
//- Gamma of cells from neighbour processors
mutable scalarListList* ibProcGammaPtr_;
//- Centres of cells from other processors
mutable vectorList* ibProcCentresPtr_;
//- Cell-proc-cell addressing
mutable List<List<labelPair> >* ibCellProcCellsPtr_;
//- Gamma of cells from other processors (needed for sampling)
mutable scalarList* ibProcGammaPtr_;
//- Neighbour IB cell-proc-cell addressing
mutable labelListList* ibCellProcCellsPtr_;
//- Dead cells list
mutable labelList* deadCellsPtr_;
@ -221,6 +232,13 @@ class immersedBoundaryFvPatch
mutable dynamicLabelList triFacesInMesh_;
// Cell surface search
//- Cell octree for search in local region
// Used in donor cell search
mutable indexedOctree<treeDataCell>* cellSearchPtr_;
// Private Member Functions
@ -229,6 +247,9 @@ class immersedBoundaryFvPatch
//- Clear all demand-driven data
void clearOut();
//- Return cell set size estimate for hash sets
label csEst() const;
// Make demand-driven data
@ -297,10 +318,11 @@ class immersedBoundaryFvPatch
//- Make triangular surface face area vectors
void makeTriSf() const;
//- Find nearest cell
label findNearestCell(const point& location) const;
//- Calculate cell octree
void calcCellSearch() const;
//- Return extended cell-cell addressing
//- Prepare a list of cell candidates for cell-cell addressing
// The list is sorted in increasing distance from pt
void findCellCells
(
const vector& pt,
@ -446,24 +468,27 @@ public:
//- Return IB cell extended stencil
const labelListList& ibCellCells() const;
//- Return neighbour proc centres
// These are centres from neighbouring processors the
// local processor needs to receive
const vectorListList& ibProcCentres() const;
//- Return neighbour proc gamma
// These are gamma values from neighbouring processors the
// local processor needs to receive
const scalarListList& ibProcGamma() const;
//- Return neighbour cell addressing
const List<List<labelPair> >& ibCellProcCells() const;
//- Return neighbour proc cells
// These are cell data that other processors need from
// local processor
const labelListList& ibProcCells() const;
//- Processor mapDistribute
const mapDistribute& ibMap() const;
//- Return neighbour proc centres
// These are centres from neighbouring processors the
// local processor needs to receive
const vectorList& ibProcCentres() const;
//- Return neighbour cell addressing
const labelListList& ibCellProcCells() const;
//- Return neighbour proc gamma
// These are gamma values from neighbouring processors the
// local processor needs to receive
const scalarList& ibProcGamma() const;
//- Return dead cells
const labelList& deadCells() const;
@ -510,14 +535,11 @@ public:
const dynamicLabelList& triFacesInMesh() const;
// Parallel communication
// Donor cell search functionality
//- Send and receive
template<class Type>
tmp<FieldField<Field, Type> > sendAndReceive
(
const Field<Type>& psi
) const;
//- Return cell octree for search in local region
// Only eligible cells are in the tree
const indexedOctree<treeDataCell>& cellSearch() const;
// Interpolation functions to and from triangular surface

31
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatchLeastSquaresFit.C
View file

@ -52,7 +52,7 @@ void Foam::immersedBoundaryFvPatch::makeInvDirichletMatrices() const
// Get addressing
const labelList& ibc = ibCells();
const labelListList& ibcc = ibCellCells();
const List<List<labelPair> >& ibcProcC = ibCellProcCells();
const labelListList& ibcProcC = ibCellProcCells();
const vectorField& ibp = ibPoints();
invDirichletMatricesPtr_ =
@ -66,7 +66,7 @@ void Foam::immersedBoundaryFvPatch::makeInvDirichletMatrices() const
// Initialize maxRowSum for debug
scalar maxRowSum = 0.0;
const vectorListList& procC = ibProcCentres();
const vectorList& procC = ibProcCentres();
label nCoeffs = 5;
@ -78,7 +78,7 @@ void Foam::immersedBoundaryFvPatch::makeInvDirichletMatrices() const
forAll (idm, cellI)
{
const labelList& interpCells = ibcc[cellI];
const List<labelPair>& interpProcCells = ibcProcC[cellI];
const labelList& interpProcCells = ibcProcC[cellI];
vectorField allPoints
(
@ -107,14 +107,7 @@ void Foam::immersedBoundaryFvPatch::makeInvDirichletMatrices() const
// Processor cells
for (label i = 0; i < interpProcCells.size(); i++)
{
allPoints[pointID++] =
procC
[
interpProcCells[i].first()
]
[
interpProcCells[i].second()
];
allPoints[pointID++] = procC[interpProcCells[i]];
}
// Weights calculation
@ -286,7 +279,7 @@ void Foam::immersedBoundaryFvPatch::makeInvNeumannMatrices() const
// Get addressing
const labelList& ibc = ibCells();
const labelListList& ibcc = ibCellCells();
const List<List<labelPair> >& ibcProcC = ibCellProcCells();
const labelListList& ibcProcC = ibCellProcCells();
const vectorField& ibp = ibPoints();
// Note: the algorithm is originally written with inward-facing normals
@ -305,7 +298,7 @@ void Foam::immersedBoundaryFvPatch::makeInvNeumannMatrices() const
// Initialize maxRowSum for debug
scalar maxRowSum = 0.0;
const vectorListList& procC = ibProcCentres();
const vectorList& procC = ibProcCentres();
label nCoeffs = 6;
@ -317,7 +310,7 @@ void Foam::immersedBoundaryFvPatch::makeInvNeumannMatrices() const
forAll (inm, cellI)
{
const labelList& interpCells = ibcc[cellI];
const List<labelPair>& interpProcCells = ibcProcC[cellI];
const labelList& interpProcCells = ibcProcC[cellI];
vectorField allPoints
(
@ -340,14 +333,7 @@ void Foam::immersedBoundaryFvPatch::makeInvNeumannMatrices() const
// Processor cells
for (label i = 0; i < interpProcCells.size(); i++)
{
allPoints[pointID++] =
procC
[
interpProcCells[i].first()
]
[
interpProcCells[i].second()
];
allPoints[pointID++] = procC[interpProcCells[i]];
}
// Weights calculation
@ -359,7 +345,6 @@ void Foam::immersedBoundaryFvPatch::makeInvNeumannMatrices() const
// Calculate weights
scalarField W = 1 - curR/(1.1*max(curR));
inm.set
(
cellI,

27
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatchSamplingWeights.C
View file

@ -50,7 +50,7 @@ void Foam::immersedBoundaryFvPatch::makeIbSamplingWeights() const
// Get addressing
const labelList& ibc = ibCells();
const labelListList& ibcc = ibCellCells();
const List<List<labelPair> >& ibcProcC = ibCellProcCells();
const labelListList& ibcProcC = ibCellProcCells();
// Initialise the weights
ibSamplingWeightsPtr_ = new scalarListList(ibc.size());
@ -74,8 +74,8 @@ void Foam::immersedBoundaryFvPatch::makeIbSamplingWeights() const
const scalarField& gammaIn = gamma().internalField();
const vectorField& CIn = mesh_.C().internalField();
const scalarListList& gammaProc = ibProcGamma();
const vectorListList& CProc = ibProcCentres();
const scalarList& gammaProc = ibProcGamma();
const vectorList& CProc = ibProcCentres();
// Go through all cellCells and calculate inverse distance for
// all live points
@ -105,7 +105,7 @@ void Foam::immersedBoundaryFvPatch::makeIbSamplingWeights() const
}
// Processor weights
const List<labelPair>& interpProcCells = ibcProcC[cellI];
const labelList& interpProcCells = ibcProcC[cellI];
scalarList& curProcCW = cellProcWeights[cellI];
@ -113,26 +113,11 @@ void Foam::immersedBoundaryFvPatch::makeIbSamplingWeights() const
{
if
(
gammaProc
[
interpProcCells[cProcI].first()
]
[
interpProcCells[cProcI].second()
] > SMALL
gammaProc[interpProcCells[cProcI]] > SMALL
)
{
curProcCW[cProcI] =
1/mag
(
CProc
[
interpProcCells[cProcI].first()
]
[
interpProcCells[cProcI].second()
] - curP
);
1/mag(CProc[interpProcCells[cProcI]]);
sumW += curProcCW[cProcI];
}

115
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatch/immersedBoundaryFvPatchTemplates.C
View file

@ -28,86 +28,6 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
Foam::tmp<Foam::FieldField<Foam::Field, Type> >
Foam::immersedBoundaryFvPatch::sendAndReceive
(
const Field<Type>& psi
) const
{
tmp<FieldField<Field, Type> > tprocPsi
(
new FieldField<Field, Type>(Pstream::nProcs())
);
FieldField<Field, Type>& procPsi = tprocPsi();
// This requires a rewrite useng mapDistribute
// HJ, 11/Aug/2016
forAll (procPsi, procI)
{
procPsi.set
(
procI,
new Field<Type>
(
ibProcCentres()[procI].size(),
pTraits<Type>::zero
)
);
}
if (Pstream::parRun())
{
// Send
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
if (procI != Pstream::myProcNo())
{
// Do not send empty lists
if (!ibProcCells()[procI].empty())
{
Field<Type> curPsi(psi, ibProcCells()[procI]);
// Parallel data exchange
OPstream toProc
(
Pstream::blocking,
procI,
curPsi.size()*sizeof(Type)
);
toProc << curPsi;
}
}
}
// Receive
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
if (procI != Pstream::myProcNo())
{
// Do not receive empty lists
if (!procPsi[procI].empty())
{
// Parallel data exchange
IPstream fromProc
(
Pstream::blocking,
procI,
procPsi[procI].size()*sizeof(Type)
);
fromProc >> procPsi[procI];
}
}
}
}
return tprocPsi;
}
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::immersedBoundaryFvPatch::toIbPoints
@ -206,8 +126,8 @@ Foam::immersedBoundaryFvPatch::toTriFaces
<< abort(FatalError);
}
const labelListList& ctfAddr = cellsToTriAddr();
const scalarListList& ctfWeights = cellsToTriWeights();
const labelListList& ctfAddr = this->cellsToTriAddr();
const scalarListList& ctfWeights = this->cellsToTriWeights();
tmp<Field<Type> > tIbPsi
(
@ -270,13 +190,13 @@ Foam::immersedBoundaryFvPatch::toSamplingPoints
}
// Get addressing
const labelList& ibc = ibCells();
const labelListList& ibcc = ibCellCells();
const List<List<labelPair> >& ibcProcC = ibCellProcCells();
const labelList& ibc = this->ibCells();
const labelListList& ibcc = this->ibCellCells();
const labelListList& ibcProcC = this->ibCellProcCells();
// Get weights
const scalarListList& cellWeights = ibSamplingWeights();
const scalarListList& cellProcWeights = ibSamplingProcWeights();
const scalarListList& cellWeights = this->ibSamplingWeights();
const scalarListList& cellProcWeights = this->ibSamplingProcWeights();
tmp<Field<Type> > tIbPsi
(
@ -296,26 +216,25 @@ Foam::immersedBoundaryFvPatch::toSamplingPoints
}
}
// Parallel communication for psi
FieldField<Field, Type> procCellValues = sendAndReceive(cellValues);
Field<Type> procCellValues;
// Parallel communication for psi, using mapDistribute
if (Pstream::parRun())
{
procCellValues = cellValues;
ibMap().distribute(procCellValues);
}
// Do interpolation, cell data from other processors
forAll (ibc, cellI)
{
const List<labelPair>& curProcCells = ibcProcC[cellI];
const labelList& curProcCells = ibcProcC[cellI];
const scalarList& curProcWeights = cellProcWeights[cellI];
forAll (curProcCells, cpcI)
{
ibPsi[cellI] +=
curProcWeights[cpcI]*
procCellValues
[
curProcCells[cpcI].first()
]
[
curProcCells[cpcI].second()
];
curProcWeights[cpcI]*procCellValues[curProcCells[cpcI]];
}
}

49
src/immersedBoundary/immersedBoundary/immersedBoundaryFvPatchField/immersedBoundaryFvPatchField.C
View file

@ -92,7 +92,7 @@ immersedBoundaryFvPatchField<Type>::imposeDirichletCondition() const
// Get addressing
const labelList& ibc = ibPatch_.ibCells();
const labelListList& ibCellCells = ibPatch_.ibCellCells();
const List<List<labelPair> >& ibCellProcCells = ibPatch_.ibCellProcCells();
const labelListList& ibCellProcCells = ibPatch_.ibCellProcCells();
const PtrList<scalarRectangularMatrix>& invMat =
ibPatch_.invDirichletMatrices();
@ -140,7 +140,17 @@ immersedBoundaryFvPatchField<Type>::imposeDirichletCondition() const
counter++;
// Parallel communication for psi
FieldField<Field, Type> procPsi = ibPatch_.sendAndReceive(psiI);
// Parallel communication for psi
List<Type> procPsi;
if (Pstream::parRun())
{
// Get reference to mapDistribute for global comms
const mapDistribute& ibm = ibPatch_.ibMap();
procPsi = psiI;
ibm.distribute(procPsi);
}
// Prepare error normalisation
scalar maxMagPolyPsi = 0;
@ -151,7 +161,7 @@ immersedBoundaryFvPatchField<Type>::imposeDirichletCondition() const
const labelList& curCells = ibCellCells[cellI];
const List<labelPair>& curProcCells = ibCellProcCells[cellI];
const labelList& curProcCells = ibCellProcCells[cellI];
const scalarRectangularMatrix& curInvMatrix = invMat[cellI];
@ -171,14 +181,7 @@ immersedBoundaryFvPatchField<Type>::imposeDirichletCondition() const
for (label i = 0; i < curProcCells.size(); i++)
{
source[pointID++] =
procPsi
[
curProcCells[i].first()
]
[
curProcCells[i].second()
]
- ibValue_[cellI];
procPsi[curProcCells[i]] - ibValue_[cellI];
}
for (label i = 0; i < nCoeffs; i++)
@ -288,7 +291,7 @@ immersedBoundaryFvPatchField<Type>::imposeNeumannCondition() const
const labelListList& ibCellCells = ibPatch_.ibCellCells();
const List<List<labelPair> >& ibCellProcCells = ibPatch_.ibCellProcCells();
const labelListList& ibCellProcCells = ibPatch_.ibCellProcCells();
const PtrList<scalarRectangularMatrix>& invMat =
ibPatch_.invNeumannMatrices();
@ -336,7 +339,16 @@ immersedBoundaryFvPatchField<Type>::imposeNeumannCondition() const
counter++;
// Parallel communication for psi
FieldField<Field, Type> procPsi = ibPatch_.sendAndReceive(psiI);
List<Type> procPsi;
if (Pstream::parRun())
{
// Get reference to mapDistribute for global comms
const mapDistribute& ibm = ibPatch_.ibMap();
procPsi = psiI;
ibm.distribute(procPsi);
}
// Prepare error normalisation
scalar maxMagPolyPsi = 0;
@ -346,7 +358,7 @@ immersedBoundaryFvPatchField<Type>::imposeNeumannCondition() const
label curCell = ibc[cellI];
const labelList& curCells = ibCellCells[cellI];
const List<labelPair>& curProcCells = ibCellProcCells[cellI];
const labelList& curProcCells = ibCellProcCells[cellI];
const scalarRectangularMatrix& curInvMatrix = invMat[cellI];
@ -367,14 +379,7 @@ immersedBoundaryFvPatchField<Type>::imposeNeumannCondition() const
for (label i = 0; i < curProcCells.size(); i++)
{
source[pointID++] =
procPsi
[
curProcCells[i].first()
]
[
curProcCells[i].second()
];
source[pointID++] = procPsi[curProcCells[i]];
}
for (label i = 0; i < nCoeffs; i++)