Overset AMG support. Teseted serial version

This commit is contained in:
Hrvoje Jasak 2019-09-16 19:05:41 +01:00
parent 9537d07c25
commit 89be468228
16 changed files with 2043 additions and 42 deletions

View file

@ -31,6 +31,9 @@ oversetMesh/oversetMeshAddressing.C
oversetLduInterface/oversetLduInterface.C
oversetLduInterfaceField/oversetLduInterfaceField.C
oversetAMGInterface/oversetAMGInterface.C
oversetAMGInterfaceField/oversetAMGInterfaceField.C
oversetPolyPatch/oversetPolyPatch.C
oversetPointPatch/oversetPointPatch.C

View file

@ -0,0 +1,906 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "oversetAMGInterface.H"
#include "oversetMesh.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(oversetAMGInterface, 0);
addToRunTimeSelectionTable
(
AMGInterface,
oversetAMGInterface,
lduInterface
);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::oversetAMGInterface::initMap() const
{
if (mapPtr_)
{
FatalErrorInFunction
<< "map already calculated"
<< abort(FatalError);
}
if (!Pstream::parRun())
{
FatalErrorInFunction
<< "Requested calculation of send-receive addressing for a "
<< "serial run. This is not allowed"
<< abort(FatalError);
}
// Create send map: what I am sending to which processor
// There is no need to use a map: guess and resize
labelListList sendMap(Pstream::nProcs());
// Memory management
{
forAll (sendMap, procI)
{
sendMap[procI].setSize(donorCells_.size());
}
// Count number of entries
labelList nSendMap(Pstream::nProcs(), 0);
forAll (donorCellsProc_, dcpI)
{
// First index = proc
const label toProc = donorCellsProc_[dcpI];
// Second index = fill size
// Content = donor cell index
sendMap[toProc][nSendMap[toProc]] = dcpI;
nSendMap[toProc]++;
}
// Resize
forAll (sendMap, procI)
{
sendMap[procI].setSize(nSendMap[procI]);
}
}
// Create construct map: what I am receiving from each processor
labelListList constructMap(Pstream::nProcs());
// Memory management
{
forAll (constructMap, procI)
{
constructMap[procI].setSize(acceptorCells().size());
}
labelList nConstructMap(Pstream::nProcs(), 0);
forAll (donorProcForAcceptor_, dpaI)
{
// First index = proc
const label fromProc = donorProcForAcceptor_[dpaI];
// Second index = fill size
// Content = acceptor cell index
constructMap[fromProc][nConstructMap[fromProc]] = dpaI;
nConstructMap[fromProc]++;
}
// Resize
forAll (constructMap, procI)
{
constructMap[procI].setSize(nConstructMap[procI]);
}
}
mapPtr_ = new mapDistribute
(
acceptorCells().size(),
sendMap,
constructMap
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::oversetAMGInterface::oversetAMGInterface
(
const lduPrimitiveMesh& lduMesh,
const lduInterfacePtrsList& coarseInterfaces,
const lduInterface& fineInterface,
const labelField& localRestrictAddressing,
const labelField& neighbourRestrictAddressing
)
:
AMGInterface(lduMesh),
fineOversetInterface_(refCast<const oversetLduInterface>(fineInterface)),
mapPtr_(nullptr)
{
// Info<< "Constructing overset interface" << nl
// << "local restrict: " << localRestrictAddressing.size() << nl
// << "neighbour restrict: " << neighbourRestrictAddressing.size()
// << endl;
// Size of coarse interface is equal to max cluster size plus one
interfaceSize_ = max(localRestrictAddressing) + 1;
// Note: cannot deal with weighted interpolation because access
// to interpolation weights is coded under the name of the field being
// solved for, which is not available for AMG coarsening
// The coarsening across overset interface shall be done using only the
// master acceptor and a single weight
// HJ, 11/Sep/2019
// Initialise fine map
// Note: the mapDistribute contains a waitRequests call which cannot
// be limited to a comm. Therefore, all processors need to calculate
// the mapDistribute schedule before escaping the constructor,
// even if there are no overset cells available.
// HJ, 18/Oct/2016
if (Pstream::parRun())
{
fineOversetInterface_.map().schedule();
}
// Coarsen local donor cells and renumber. Record local donor index
// and communicate across
// Get fine donors
const labelList& fineDonors = fineOversetInterface_.donorCells();
// Get fine donors processor (to which donor is sent)
const labelList& fineDonorsProc = fineOversetInterface_.donorCellsProc();
// Memory management
{
// Into the hash table, using the coarse donor, record coarse processor
// from fine donor. They must be the same
HashTable<label, label, Hash<label> > coarseDonorProcMap;
// For all fine donors, find and collect restrict addressing.
// This will be the coarse donor
forAll (fineDonors, fdI)
{
coarseDonorProcMap.insert
(
localRestrictAddressing[fineDonors[fdI]],
fineDonorsProc[fdI]
);
}
// Collect coarse donor cells
donorCells_ = coarseDonorProcMap.sortedToc();
// Pick up coarse donorCellsProc from the fine
donorCellsProc_.setSize(donorCells_.size());
// This could be done with an iterator as well
// HJ, 13/Sep/2019
forAll (donorCells_, dcI)
{
donorCellsProc_[dcI] = coarseDonorProcMap.find(donorCells_[dcI])();
}
}
// Having established coarse donor list (in ascending order),
// for each fine donor find the position of the coarse donor in the
// list
// I will know coarse donor index and I need its location in the list
labelList coarseDonorIndex(max(donorCells_) + 1);
forAll (donorCells_, dcI)
{
coarseDonorIndex[donorCells_[dcI]] = dcI;
}
// Make a list across fine level, where each fine donor records
// coarse donor index
labelList acceptorCellDonorIndex
(
fineOversetInterface_.interfaceSize(),
-1
);
// Fill the list: for all fine donors, look up coarse donor,
// look up coarse donor index and put it into acceptorCellDonorIndex
forAll (fineDonors, fdI)
{
// Ger coarse donor index
const label cdi = localRestrictAddressing[fineDonors[fdI]];
// In the acceptor cell donor list, put coarse donor index
acceptorCellDonorIndex[fineDonors[fdI]] = coarseDonorIndex[cdi];
}
// Communicate coarse donor
acceptorCellDonorIndex =
fineOversetInterface_.internalFieldTransfer
(
Pstream::blocking,
acceptorCellDonorIndex
);
// Donor processor info needs to be communicated.
// Use fine level comms
labelList acceptorCellDonorProc
(
fineOversetInterface_.interfaceSize(),
Pstream::myProcNo()
);
// Communicate donor proc
acceptorCellDonorProc =
fineOversetInterface_.internalFieldTransfer
(
Pstream::blocking,
acceptorCellDonorProc
);
// Note: Guessing size of HashTable to fine interface size
// Coded neighbour index. Note: using long int to simplify encoding
// HJ, 1/Aug/2016
HashTable<DynamicList<long, 4>, long, Hash<long> > neighboursTable
(
Foam::max(128, fineOversetInterface_.interfaceSize()/4)
);
// Neighbour processor index
HashTable<DynamicList<label, 4>, label, Hash<label> > nbrsProcTable
(
Foam::max(128, fineOversetInterface_.interfaceSize()/4)
);
// Neighbour processor donor index
HashTable<DynamicList<label, 4>, label, Hash<label> > nbrsProcDonorIdTable
(
Foam::max(128, fineOversetInterface_.interfaceSize()/4)
);
// Neighbour face-faces addressing for a face with split neighbours
HashTable<DynamicList<DynamicList<label, 4>, 4>, label, Hash<label> >
faceFaceTable
(
Foam::max(128, fineOversetInterface_.interfaceSize()/4)
);
// Count the number of coarse faces
label nCoarseFaces = 0;
// Count the number of agglomeration pairs
label nAgglomPairs = 0;
// On the fine level, addressing is obtained from oversetMesh
if (fineOversetInterface_.fineLevel())
{
// Loop over all fringe faces
// Follow the pattern in oversetFvPatchField::initInterfaceMatrixUpdate
// Get mesh owner-neighbour addressing to visit cells around fringe
// faces
const unallocLabelList& own = lduMesh.lowerAddr();
const unallocLabelList& nei = lduMesh.upperAddr();
const labelList& fringeFaces =
fineOversetInterface_.overset().fringeFaces();
const labelList& fringeFaceCells =
fineOversetInterface_.overset().fringeFaceCells();
const boolList& fringeFaceFlips =
fineOversetInterface_.overset().fringeFaceFlips();
label curMasterProc, curSlaveProc, curSlaveProcDonorId;
long curMaster, curSlave;
// On the fine level, loop through fringe faces:
// Overset interfaces are between internal cells
forAll (fringeFaces, ffI)
{
const label& curFace = fringeFaces[ffI];
if (curFace < own.size())
{
curMaster = -1;
curMasterProc = -1;
curSlave = -1;
curSlaveProc = -1;
curSlaveProcDonorId = -1;
// Note. Signalling in global clustering requires
// me to recognise clustering from separate
// processors as separate. In the first phase,
// this will be used to recognise cluster from
// each processor as separate and in the second
// phase it will be used to filter local processor
// faces from the global patch. Currently, I am
// calculating unique cluster index as:
//
// id = cluster + procOffset*myProcID
// HJ, 1/Apr/2009
// Get live cell restrict
if (fringeFaceFlips[ffI])
{
// Face pointing into live cell: ie neighbour is live
curMaster = localRestrictAddressing[nei[curFace]];
}
else
{
// Face pointing out of live cell
curMaster = localRestrictAddressing[own[curFace]];
}
// Donor cell data is accessible from fringe
// fringeFaceCells gives acceptor index
curMasterProc = Pstream::myProcNo();
curSlave = neighbourRestrictAddressing[fringeFaceCells[ffI]];
curSlaveProc = acceptorCellDonorProc[fringeFaceCells[ffI]];
curSlaveProcDonorId =
acceptorCellDonorIndex[fringeFaceCells[ffI]];
// Info<< "A: " << curMaster << tab << curMasterProc << tab << curSlave << tab << curSlaveProc << tab << curSlaveProcDonorId << endl;
// Code in current master and slave
curMaster += procOffset*curMasterProc;
curSlave += procOffset*curSlaveProc;
// Look for the master cell. If it has already got a face,
// add the coefficient to the face. If not, create a new
// face
if (neighboursTable.found(curMaster))
{
// This master side face already exists
// Check all current neighbours to see if the current
// slave already exists. If so, add the coefficient.
DynamicList<long, 4>& curNbrs =
neighboursTable.find(curMaster)();
DynamicList<label, 4>& curNbrsProc =
nbrsProcTable.find(curMaster)();
DynamicList<label, 4>& curNbrsProcDonorId =
nbrsProcDonorIdTable.find(curMaster)();
DynamicList<DynamicList<label, 4>, 4>& curFaceFaces =
faceFaceTable.find(curMaster)();
// Search for coded neighbour
bool nbrFound = false;
forAll (curNbrs, curNbrI)
{
// Check neighbour slave
if (curNbrs[curNbrI] == curSlave)
{
nbrFound = true;
curFaceFaces[curNbrI].append(ffI);
// New agglomeration pair found in already
// existing pair
nAgglomPairs++;
break;
}
}
if (!nbrFound)
{
curNbrs.append(curSlave);
curNbrsProc.append(curSlaveProc);
curNbrsProcDonorId.append(curSlaveProcDonorId);
DynamicList<label, 4> newFF;
newFF.append(ffI);
curFaceFaces.append(newFF);
// New coarse face created for an existing master
nCoarseFaces++;
nAgglomPairs++;
}
}
else
{
// This master has got no neighbours yet.
// Add a neighbour, proc and a coefficient as a
// new list, thus creating a new face
DynamicList<long, 4> newNbrs;
newNbrs.append(curSlave);
neighboursTable.insert
(
curMaster,
newNbrs
);
DynamicList<label, 4> newNbrsProc;
newNbrsProc.append(curSlaveProc);
nbrsProcTable.insert
(
curMaster,
newNbrsProc
);
DynamicList<label, 4> newNbrsProcDonorId;
newNbrsProcDonorId.append(curSlaveProcDonorId);
nbrsProcDonorIdTable.insert
(
curMaster,
newNbrsProcDonorId
);
DynamicList<DynamicList<label, 4>, 4> newFF;
newFF.append(DynamicList<label, 4>());
newFF[0].append(ffI);
faceFaceTable.insert
(
curMaster,
newFF
);
// New coarse face created for a new master
nCoarseFaces++;
nAgglomPairs++;
}
} // end for all current neighbours
} // end for all fine faces
}
else
{
// Coarse level, addressing is stored in faceCells
// Perform analysis only for local faces
// HJ, 22/Jun/2016
label curMasterProc, curSlaveProc, curSlaveProcDonorId;
long curMaster, curSlave;
// Size check. Remove on debug
if
(
localRestrictAddressing.size()
!= fineOversetInterface_.interfaceSize()
)
{
FatalErrorInFunction
<< "Size check failed. localRestrictAddressing: "
<< localRestrictAddressing.size()
<< " faceCells: "
<< fineOversetInterface_.interfaceSize()
<< abort(FatalError);
}
// On the coarse level, acceptor cells are in faceCells_
// Donor data is served in the same order as acceptors: this is a
// matched interface
const labelList& fineAcceptorCells =
fineOversetInterface_.acceptorCells();
forAll (fineAcceptorCells, ffI)
{
curMaster = -1;
curMasterProc = -1;
curSlave = -1;
curSlaveProc = -1;
curSlaveProcDonorId = -1;
// Note. Signalling in global clustering requires
// me to recognise clustering from separate
// processors as separate. In the first phase,
// this will be used to recognise cluster from
// each processor as separate and in the second
// phase it will be used to filter local processor
// faces from the global patch. Currently, I am
// calculating unique cluster index as:
//
// id = cluster + procOffset*myProcID
// HJ, 1/Apr/2009
// Not sure if master-slave switching is needed here:
// each side keeps acceptors as local faceCells and they
// communicate not in pars but on a more complex pattern
// Master side
curMaster = localRestrictAddressing[fineAcceptorCells[ffI]];
curMasterProc = Pstream::myProcNo();
curSlave = neighbourRestrictAddressing[ffI];
curSlaveProc = acceptorCellDonorProc[ffI];
curSlaveProcDonorId = acceptorCellDonorIndex[ffI];
// Info<< "B: " << curMaster << tab << curMasterProc << tab << curSlave << tab << curSlaveProc << tab << curSlaveProcDonorId << endl;
// Code in current master and slave
curMaster += procOffset*curMasterProc;
curSlave += procOffset*curSlaveProc;
// Look for the master cell. If it has already got a face,
// add the coefficient to the face. If not, create a new face.
if (neighboursTable.found(curMaster))
{
// This master side face already exists
// Check all current neighbours to see if the current slave
// already exists and if so, add the fine face
// to the agglomeration.
DynamicList<long, 4>& curNbrs =
neighboursTable.find(curMaster)();
DynamicList<label, 4>& curNbrsProc =
nbrsProcTable.find(curMaster)();
DynamicList<label, 4>& curNbrsProcDonorId =
nbrsProcDonorIdTable.find(curMaster)();
DynamicList<DynamicList<label, 4>, 4>& curFaceFaces =
faceFaceTable.find(curMaster)();
// Search for coded neighbour
bool nbrFound = false;
forAll (curNbrs, curNbrI)
{
// Check neighbour slave
if (curNbrs[curNbrI] == curSlave)
{
nbrFound = true;
curFaceFaces[curNbrI].append(ffI);
// New agglomeration pair found in already
// existing pair
nAgglomPairs++;
break;
}
}
if (!nbrFound)
{
curNbrs.append(curSlave);
curNbrsProc.append(curSlaveProc);
curNbrsProcDonorId.append(curSlaveProcDonorId);
DynamicList<label, 4> newFF;
newFF.append(ffI);
curFaceFaces.append(newFF);
// New coarse face created for an existing master
nCoarseFaces++;
nAgglomPairs++;
}
}
else
{
// This master has got no neighbours yet. Add a neighbour
// and a coefficient, thus creating a new face
DynamicList<long, 4> newNbrs;
newNbrs.append(curSlave);
neighboursTable.insert
(
curMaster,
newNbrs
);
DynamicList<label, 4> newNbrsProc;
newNbrsProc.append(curSlaveProc);
nbrsProcTable.insert
(
curMaster,
newNbrsProc
);
DynamicList<label, 4> newNbrsProcDonorId;
newNbrsProcDonorId.append(curSlaveProcDonorId);
nbrsProcDonorIdTable.insert
(
curMaster,
newNbrsProcDonorId
);
DynamicList<DynamicList<label, 4>, 4> newFF;
newFF.append(DynamicList<label, 4>());
newFF[0].append(ffI);
faceFaceTable.insert
(
curMaster,
newFF
);
// New coarse face created for a new master
nCoarseFaces++;
nAgglomPairs++;
}
} // end for all fine faces
} // end of else in fine level (coarse level)
// Resize the lists
faceCells_.setSize(nCoarseFaces);
fineAddressing_.setSize(nAgglomPairs, -1);
restrictAddressing_.setSize(nAgglomPairs, -11);
// All weights are equal to 1: integral matching
restrictWeights_.setSize(nAgglomPairs, 1.0);
// Record donor processor for each acceptor (in faceCells_)
// to establish communication
donorCellForAcceptor_.setSize(nCoarseFaces);
donorProcForAcceptor_.setSize(nCoarseFaces);
donorProcIndexForAcceptor_.setSize(nCoarseFaces);
List<long> contents = neighboursTable.sortedToc();
// Note:
// Each patch has both donors and acceptors.
// Acceptors are master side: faceCells_ needs to address into the
// acceptors
// Donors are on the slave side: this is the data that will need to be
// sent accross to the donor
// Reset face counter for re-use
nCoarseFaces = 0;
nAgglomPairs = 0;
// On master side, the owner addressing is stored in table of contents
forAll (contents, masterI)
{
const label curMaster = contents[masterI];
const DynamicList<long, 4>& curNbrs =
neighboursTable.find(curMaster)();
const DynamicList<label, 4>& curNbrsProc =
nbrsProcTable.find(curMaster)();
const DynamicList<label, 4>& curNbrsProcDonorId =
nbrsProcDonorIdTable.find(curMaster)();
const DynamicList<DynamicList<label, 4>, 4>& curFaceFaces =
faceFaceTable.find(curMaster)();
forAll (curNbrs, curNbrI)
{
// Record new coarse acceptor
faceCells_[nCoarseFaces] = curMaster;
// Record donor for coarse acceptor (can be on different processor)
donorCellForAcceptor_[nCoarseFaces] = curNbrs[curNbrI];
// Record donor processor
donorProcForAcceptor_[nCoarseFaces] = curNbrsProc[curNbrI];
// Record donor processir index for acceptor
donorProcIndexForAcceptor_[nCoarseFaces] =
curNbrsProcDonorId[curNbrI];
// Get faces and weights
const DynamicList<label, 4>& facesIter = curFaceFaces[curNbrI];
forAll (facesIter, facesIterI)
{
fineAddressing_[nAgglomPairs] = facesIter[facesIterI];
restrictAddressing_[nAgglomPairs] = nCoarseFaces;
nAgglomPairs++;
}
nCoarseFaces++;
}
}
// Info<< "Finished level. Sizes " << nl
// << "faceCells_: " << faceCells_.size()
// << " fineAddressing_: " << fineAddressing_.size()
// << " restrictAddressing_: " << restrictAddressing_.size()
// << " donorCells_: " << donorCells_.size()
// << " donorCellForAcceptor_: " << donorCellForAcceptor_.size()
// << " donorProcIndexForAcceptor_: " << donorProcIndexForAcceptor_.size()
// << endl;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::oversetAMGInterface::~oversetAMGInterface()
{
deleteDemandDrivenData(mapPtr_);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::labelField> Foam::oversetAMGInterface::interfaceInternalField
(
const unallocLabelList& iF
) const
{
// Return complete field: needed for both donors and acceptors
// HJ, 16/Sep/2019
return tmp<labelField>(new labelField(iF));
}
Foam::tmp<Foam::scalarField> Foam::oversetAMGInterface::agglomerateCoeffs
(
const scalarField& fineCoeffs
) const
{
// Escape agglomerating internal coefficients: zero size
// HJ, 16/Sep/2019
if (fineCoeffs.empty())
{
return tmp<scalarField>(new scalarField());
}
tmp<scalarField> tcoarseCoeffs(new scalarField(size(), 0.0));
scalarField& coarseCoeffs = tcoarseCoeffs();
// Added weights to account for non-integral matching
forAll (restrictAddressing_, ffi)
{
coarseCoeffs[restrictAddressing_[ffi]] +=
restrictWeights_[ffi]*fineCoeffs[fineAddressing_[ffi]];
}
return tcoarseCoeffs;
}
bool Foam::oversetAMGInterface::master() const
{
return fineOversetInterface_.master();
}
bool Foam::oversetAMGInterface::fineLevel() const
{
return false;
}
Foam::label Foam::oversetAMGInterface::interfaceSize() const
{
return interfaceSize_;
}
const Foam::oversetMesh& Foam::oversetAMGInterface::overset() const
{
// Overset should not be accessed from coarse levels
FatalErrorInFunction
<< "Requested fine addressing at coarse level"
<< abort(FatalError);
// Dummy return
return fineOversetInterface_.overset();
}
const Foam::mapDistribute& Foam::oversetAMGInterface::map() const
{
if (!mapPtr_)
{
initMap();
}
return *mapPtr_;
}
void Foam::oversetAMGInterface::initTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& interfaceData
) const
{}
Foam::tmp<Foam::labelField> Foam::oversetAMGInterface::transfer
(
const Pstream::commsTypes,
const unallocLabelList& interfaceData
) const
{
return labelField::null();
}
void Foam::oversetAMGInterface::initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& iF
) const
{}
Foam::tmp<Foam::labelField> Foam::oversetAMGInterface::internalFieldTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& iF
) const
{
// Repackage donor data to acceptors
// Note: result is copied as internal field, re-scaled and passed across
if (Pstream::parRun())
{
tmp<labelField> tresult(new labelField(iF));
labelList& result = tresult();
map().distribute(result);
return tresult;
}
else
{
return tmp<labelField>(new labelField(iF, donorCellForAcceptor_));
}
}
void Foam::oversetAMGInterface::initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const scalarField& iF
) const
{}
Foam::tmp<Foam::scalarField> Foam::oversetAMGInterface::internalFieldTransfer
(
const Pstream::commsTypes,
const scalarField& iF
) const
{
// Repackage donor data to acceptors
// Note: result is copied as internal field, re-scaled and passed across
if (Pstream::parRun())
{
tmp<scalarField> tresult(new scalarField(iF));
scalarList& result = tresult();
map().distribute(result);
return tresult;
}
else
{
return tmp<scalarField>(new scalarField(iF, donorCellForAcceptor_));
}
}
// ************************************************************************* //

View file

@ -0,0 +1,281 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::oversetAMGInterface
Description
AMG agglomerated overset interface.
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved.
SourceFiles
oversetAMGInterface.C
\*---------------------------------------------------------------------------*/
#ifndef oversetAMGInterface_H
#define oversetAMGInterface_H
#include "AMGInterface.H"
#include "oversetLduInterface.H"
#include "mapDistribute.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class oversetAMGInterface Declaration
\*---------------------------------------------------------------------------*/
class oversetAMGInterface
:
public AMGInterface,
virtual public oversetLduInterface
{
// Private data
//- Reference tor the oversetLduInterface from which this is
// agglomerated
const oversetLduInterface& fineOversetInterface_;
//- Coarse interface size
label interfaceSize_;
// Local acceptor data
// Local acceptor index is stored in faceCells_
//- Donor cell index (can be on different processor)
// for local acceptor
labelList donorCellForAcceptor_;
//- Donor processor for local acceptor
labelList donorProcForAcceptor_;
//- Donor processor donor index for local acceptor
labelList donorProcIndexForAcceptor_;
//- Local donor cells
labelList donorCells_;
//- Donor cell processor (to which donor is sent)
labelList donorCellsProc_;
// Parallel communication
//- Map-distribute comms tool
mutable mapDistribute* mapPtr_;
// Private Member Functions
//- Disallow default bitwise copy construct
oversetAMGInterface(const oversetAMGInterface&);
//- Disallow default bitwise assignment
void operator=(const oversetAMGInterface&);
//- Init map
void initMap() const;
// Private static data
//- Processor cluster offset index
static const long procOffset = 12000000;
public:
//- Runtime type information
TypeName("overset");
// Constructors
//- Construct from fine level interface,
// local and neighbour restrict addressing
oversetAMGInterface
(
const lduPrimitiveMesh& lduMesh,
const lduInterfacePtrsList& coarseInterfaces,
const lduInterface& fineInterface,
const labelField& localRestrictAddressing,
const labelField& neighbourRestrictAddressing
);
//- Destructor
virtual ~oversetAMGInterface();
// Member Functions
// Access
//- Return true if interface is coupled
virtual bool coupled() const
{
return true;
}
// Communications support
//- Return the values of the given internal data adjacent to
// the interface as a field
virtual tmp<labelField> interfaceInternalField
(
const unallocLabelList& internalData
) const;
// Agglomeration
//- Agglomerating given fine-level coefficients
// In overset, special treatment is required for upper,
// which needs a virtual function. HJ, 16/Sep/2019
virtual tmp<scalarField> agglomerateCoeffs
(
const scalarField& fineCoeffs
) const;
// Interface transfer functions
//- Return access to overset mesh.
// Available only on the finest level
virtual const oversetMesh& overset() const;
//- Return acceptor cells
virtual const labelList& acceptorCells() const
{
return faceCells_;
}
//- Return donor cells for coarse-level addressing
virtual const labelList& donorCells() const
{
return donorCells_;
}
//- Return donor cells processor (to which donor is sent)
virtual const labelList& donorCellsProc() const
{
return donorCellsProc_;
}
//- Return mapDistribute
virtual const mapDistribute& map() const;
//- Initialise interface data transfer (face field). Dummy
virtual void initTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& interfaceData
) const;
//- Transfer and return neighbour field (face field). Dummy
virtual tmp<labelField> transfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& interfaceData
) const;
//- Initialise transfer of internal field adjacent to the interface
// For overset, this is fringe cell data
virtual void initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& iF
) const;
//- Transfer and return internal field adjacent to the interface
// For overset, this is fringe cell data
virtual tmp<labelField> internalFieldTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& iF
) const;
//- Initialise transfer of internal field adjacent to the interface
// For overset, this is fringe cell data
virtual void initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const scalarField& iF
) const;
//- Initialise transfer of internal field adjacent to the interface
// For overset, this is fringe cell data
virtual tmp<scalarField> internalFieldTransfer
(
const Pstream::commsTypes commsType,
const scalarField& iF
) const;
//- Overset interface functions
//- Does this side own the patch ?
virtual bool master() const;
//- Is this the fine level?
virtual bool fineLevel() const;
//- Return reference tor the oversetLduInterface from which this is
// agglomerated
const oversetLduInterface& fineOversetInterface() const
{
return fineOversetInterface_;
}
//- Return interface size. On the coarse level, this is the size
// of donor list
virtual label interfaceSize() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
// # include "oversetAMGInterfaceTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -0,0 +1,157 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "oversetAMGInterface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
tmp<Field<Type> > oversetAMGInterface::fastExpand(const UList<Type>& ff) const
{
// Rewrite, 1/Jun/2016
// To avoid creating zone-sized data and gather-scatter communication
// to the master, the optimised map-distribute call is implemented.
// The field is filled with local data which is then sent where needed
// through map-distribute.
// On return, the field is expanded to zone size but only filled with
// the data which is needed for the shadow
// HJ, 1/Jun/2016
if (ff.size() != this->size())
{
FatalErrorInFunction
<< "Wrong field size. ff: " << ff.size()
<< " interface: " << this->size()
<< abort(FatalError);
}
if (localParallel() || !Pstream::parRun())
{
// Field remains identical: no parallel communications required
tmp<Field<Type> > tresult(new Field<Type>(ff));
return tresult;
}
else
{
// Optimised mapDistribute
// Execute init reduce to calculate addressing if not already done
map();
// Prepare for distribute. Note: field will be expanded to zone size
// during the distribute operation
tmp<Field<Type> > tresult(new Field<Type>(ff));
List<Type>& expand = tresult();
map().distribute(expand);
return tresult;
}
}
template<class Type>
tmp<Field<Type> > oversetAMGInterface::fastReduce(const UList<Type>& ff) const
{
// Old algorithm: OBOSLETE
// Local processor contains faceCells part of the zone and requires
// zoneAddressing part.
// For fast communications, each processor will send the faceCells and
// zoneAddressing to the master. Master will assemble global zone
// and send off messages to all processors containing only
// the required data
// HJ, 24/Jun/2011
// Rewrite, 1/Jun/2016
// To avoid creating zone-sized data and gather-scatter communication
// to the master, the optimised map-distribute call is implemented.
// The field is filled with local data which is then sent where needed
// through map-distribute.
// On return, the field is expanded to zone size but only filled with
// the data which is needed for the shadow
// Having received the zone data, shadow data is extracted from the
// field size. Note: this works only on coarse levels, where one-on-one
// mapping applies
// HJ, 1/Jun/2016
if (ff.size() != this->size())
{
FatalErrorInFunction
<< "Wrong field size. ff: " << ff.size()
<< " interface: " << this->size()
<< abort(FatalError);
}
if (localParallel() || !Pstream::parRun())
{
// Field remains identical: no parallel communications required
tmp<Field<Type> > tresult(new Field<Type>(ff));
return tresult;
}
else
{
// Optimised mapDistribute
// Execute init reduce to calculate addressing if not already done
map();
// Prepare for distribute. Note: field will be expanded to zone size
// during the distribute operation
List<Type> expand = ff;
map().distribute(expand);
const labelList& shadowZa = shadowInterface().zoneAddressing();
// Prepare return field: zone size
tmp<Field<Type> > tresult
(
new Field<Type>(shadowZa.size())
);
Field<Type>& result = tresult();
// Filter from expanded field to zone size
forAll (shadowZa, shadowZaI)
{
result[shadowZaI] = expand[shadowZa[shadowZaI]];
}
return tresult;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

View file

@ -0,0 +1,113 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "oversetAMGInterfaceField.H"
#include "oversetFvPatchFields.H"
#include "addToRunTimeSelectionTable.H"
#include "lduMatrix.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(oversetAMGInterfaceField, 0);
addToRunTimeSelectionTable
(
AMGInterfaceField,
oversetAMGInterfaceField,
lduInterface
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::oversetAMGInterfaceField::oversetAMGInterfaceField
(
const AMGInterface& AMGCp,
const lduInterfaceField& fineInterface
)
:
AMGInterfaceField(AMGCp, fineInterface),
oversetInterface_(refCast<const oversetAMGInterface>(AMGCp))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::oversetAMGInterfaceField::initInterfaceMatrixUpdate
(
const scalarField& psiInternal,
scalarField& result,
const lduMatrix& m,
const scalarField& coeffs,
const direction cmpt,
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Get psi from donors
scalarField pnf =
oversetInterface_.internalFieldTransfer
(
commsType,
psiInternal
);
const labelList& acceptorCells = oversetInterface_.acceptorCells();
if (switchToLhs)
{
forAll (acceptorCells, elemI)
{
result[acceptorCells[elemI]] += coeffs[elemI]*pnf[elemI];
}
}
else
{
forAll (acceptorCells, elemI)
{
result[acceptorCells[elemI]] -= coeffs[elemI]*pnf[elemI];
}
}
}
void Foam::oversetAMGInterfaceField::updateInterfaceMatrix
(
const scalarField&,
scalarField&,
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType,
const bool switchToLhs
) const
{
// Nothing to do
}
// ************************************************************************* //

View file

@ -0,0 +1,147 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::oversetAMGInterfaceField
Description
AMG agglomerated overset interface field.
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved.
SourceFiles
oversetAMGInterfaceField.C
\*---------------------------------------------------------------------------*/
#ifndef oversetAMGInterfaceField_H
#define oversetAMGInterfaceField_H
#include "AMGInterfaceField.H"
#include "oversetAMGInterface.H"
#include "oversetLduInterfaceField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class oversetAMGInterfaceField Declaration
\*---------------------------------------------------------------------------*/
class oversetAMGInterfaceField
:
public AMGInterfaceField,
virtual public oversetLduInterfaceField
{
// Private Data Members
//- Local reference cast into the overset interface
const oversetAMGInterface& oversetInterface_;
// Private Member Functions
// Copy control
//- Delete default bitwise copy construct
oversetAMGInterfaceField(const oversetAMGInterfaceField&) = delete;
//- Delete default bitwise assignment
void operator=(const oversetAMGInterfaceField&) = delete;
public:
//- Runtime type information
TypeName("overset");
// Constructors
//- Construct from AMG interface and fine level interface field
oversetAMGInterfaceField
(
const AMGInterface& AMGCp,
const lduInterfaceField& fineInterfaceField
);
//- Destructor
virtual ~oversetAMGInterfaceField() = default;
// Member Functions
// Access functions
// Coupled interface matrix update
//- Transform neighbour field
virtual void transformCoupleField
(
scalarField& pnf,
const direction cmpt
) const
{
// Does nothing
}
//- Initialise neighbour matrix update
virtual void initInterfaceMatrixUpdate
(
const scalarField&,
scalarField&,
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
//- Update result field based on interface functionality
virtual void updateInterfaceMatrix
(
const scalarField&,
scalarField&,
const lduMatrix&,
const scalarField&,
const direction,
const Pstream::commsTypes commsType,
const bool switchToLhs
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View file

@ -75,14 +75,58 @@ const Foam::oversetPolyPatch& Foam::oversetFvPatch::oversetPatch() const
}
bool Foam::oversetFvPatch::master() const
{
return oversetPatch_.master();
}
bool Foam::oversetFvPatch::fineLevel() const
{
return true;
}
Foam::label Foam::oversetFvPatch::interfaceSize() const
{
return boundaryMesh().mesh().nCells();
}
const Foam::labelList& Foam::oversetFvPatch::acceptorCells() const
{
return overset().acceptorCells();
}
const Foam::labelList& Foam::oversetFvPatch::donorCells() const
{
return overset().donorCells();
}
const Foam::labelList& Foam::oversetFvPatch::donorCellsProc() const
{
return overset().donorCellsProc();
}
const Foam::mapDistribute& Foam::oversetFvPatch::map() const
{
return overset().map();
}
Foam::tmp<Foam::labelField> Foam::oversetFvPatch::interfaceInternalField
(
const unallocLabelList& internalData
) const
{
// Return the copy of the parameter since we need all the mapping for
// overset
return tmp<labelField>(new labelField(internalData));
// Return all internal values
return tmp<labelField>
(
new labelField(internalData)
);
}
@ -119,11 +163,12 @@ Foam::tmp<Foam::labelField> Foam::oversetFvPatch::internalFieldTransfer
) const
{
// Repackage donor data to acceptors
// Note: result is copied as internal field, re-scaled and passed across
tmp<labelField> tresult(new labelField(iF));
labelList& result = tresult();
overset().map().distribute(result);
result = overset().acceptorMasterData(result);
return tresult;
}

View file

@ -85,7 +85,7 @@ public:
oversetFvPatch(const polyPatch& patch, const fvBoundaryMesh& bm)
:
coupledFvPatch(patch, bm),
overPolyPatch_(refCast<const oversetPolyPatch>(patch))
oversetPatch_(refCast<const oversetPolyPatch>(patch))
{}
@ -102,7 +102,7 @@ public:
virtual tmp<vectorField> delta() const;
//- Return access to overset mesh
const oversetMesh& overset() const;
virtual const oversetMesh& overset() const;
//- Const access to overset poly patch
const oversetPolyPatch& oversetPatch() const;
@ -113,23 +113,46 @@ public:
return true;
}
// Interface transfer functions
//- Is this the master side?
virtual bool master() const;
//- Is this the fine level?
virtual bool fineLevel() const;
//- Return interface size. On the fine level, this is the size
// of internal field
virtual label interfaceSize() const;
//- Return acceptor cells
virtual const labelList& acceptorCells() const;
//- Return donor cells
virtual const labelList& donorCells() const;
//- Return donor cells processor (to which donor is sent)
virtual const labelList& donorCellsProc() const;
//- Return mapDistribute
virtual const mapDistribute& map() const;
//- Return the values of the given internal data adjacent to
// the interface as a field. For overset, this is acceptor data
// the interface as a field. For overset, this is fringe cell data
virtual tmp<labelField> interfaceInternalField
(
const unallocLabelList& internalData
) const;
//- Initialise interface data transfer
//- Initialise interface data transfer (face field). Dummy
virtual void initTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& interfaceData
) const;
//- Transfer and return neighbour field
//- Transfer and return neighbour field (face field). Dummy
virtual tmp<labelField> transfer
(
const Pstream::commsTypes commsType,
@ -137,6 +160,7 @@ public:
) const;
//- Initialise neighbour field transfer
// For overset, this is donor data
virtual void initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
@ -144,7 +168,7 @@ public:
) const;
//- Return neighbour field
// For overset, this is donor data for each acceptor
// For overset, this is donor data
virtual tmp<labelField> internalFieldTransfer
(
const Pstream::commsTypes commsType,

View file

@ -44,7 +44,7 @@ void oversetFvPatchField<Type>::setCoupledFringe
const bool coupled
)
{
forAll(psi.boundaryField(), patchI)
forAll (psi.boundaryField(), patchI)
{
fvPatchField<Type>& psip = psi.boundaryField()[patchI];
@ -66,7 +66,7 @@ void oversetFvPatchField<Type>::oversetInterpolate
)
{
// Loop through boundary field and find overset patch
forAll(psi.boundaryField(), patchI)
forAll (psi.boundaryField(), patchI)
{
fvPatchField<Type>& psip = psi.boundaryField()[patchI];
@ -570,6 +570,34 @@ void oversetFvPatchField<Type>::correctOffDiag
}
}
}
// Record the coefficients into internal/boundary coeffs for AMG
// coarsening
Field<Type>& bouCoeffs = eqn.boundaryCoeffs()[this->patch().index()];
bouCoeffs.setSize(fringeFaces.size());
// intCoeffs remain at zero size: one-sided multiplication
// Careful with sign. HJ, 16/Sep/2019
forAll (fringeFaceFlips, fringeI)
{
if (fringeFaceFlips[fringeI])
{
// Face pointing into live cell
// Add overset off-diagonal contribution to live cell
bouCoeffs[fringeI] =
-pTraits<Type>::one*fringeLowerCoeffs_[fringeI];
}
else
{
// Face pointing out of live cell
// Add overset off-diagonal contribution to live cell
bouCoeffs[fringeI] =
-pTraits<Type>::one*fringeUpperCoeffs_[fringeI];
}
}
}

View file

@ -46,8 +46,12 @@ SourceFiles
namespace Foam
{
// Forward declaration of classes
class oversetMesh;
class mapDistribute;
/*---------------------------------------------------------------------------*\
Class oversetLduInterface Declaration
Class oversetLduInterface Declaration
\*---------------------------------------------------------------------------*/
class oversetLduInterface
@ -58,15 +62,63 @@ public:
//- Runtime type information
TypeName("oversetLduInterface");
// Constructors
// Destructor
//- Construct null
oversetLduInterface()
{}
virtual ~oversetLduInterface();
//- Destructor
virtual ~oversetLduInterface();
// Member Functions
// Not currently needed. HJ, 14/Jan/2009
//- Is this the master side?
virtual bool master() const = 0;
//- Is this the fine level?
virtual bool fineLevel() const = 0;
//- Return interface size. Different on fine and coarse level
// due to different comms pattern
virtual label interfaceSize() const = 0;
//- Return oversetMesh for fine-level addressing
virtual const oversetMesh& overset() const = 0;
//- Return acceptor cells
virtual const labelList& acceptorCells() const = 0;
//- Return donor cells for coarse-level addressing
virtual const labelList& donorCells() const = 0;
//- Return donor cells processor (to which donor is sent)
virtual const labelList& donorCellsProc() const = 0;
//- Return mapDistribute
virtual const mapDistribute& map() const = 0;
// Interface transfer functions
//- Initialise transfer of internal field adjacent to the interface
// For overset, this is fringe cell data
virtual void initInternalFieldTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& iF
) const = 0;
//- Transfer and return internal field adjacent to the interface
// For overset, this is fringe cell data
virtual tmp<labelField> internalFieldTransfer
(
const Pstream::commsTypes commsType,
const unallocLabelList& iF
) const = 0;
};

View file

@ -66,9 +66,8 @@ public:
{}
// Destructor
virtual ~oversetLduInterfaceField();
//- Destructor
virtual ~oversetLduInterfaceField();
// Member Functions

View file

@ -64,6 +64,7 @@ Foam::oversetMesh::oversetMesh(const fvMesh& mesh)
acceptorCellsPtr_(nullptr),
donorCellsPtr_(nullptr),
donorCellsProcPtr_(nullptr),
holeCellsPtr_(nullptr),
oversetTypesPtr_(nullptr),

View file

@ -112,6 +112,9 @@ private:
//- Donor cell labels
mutable labelList* donorCellsPtr_;
//- Donor cell processor (to which donor is sent)
mutable labelList* donorCellsProcPtr_;
//- Hole cell labels
mutable labelList* holeCellsPtr_;
@ -276,6 +279,9 @@ public:
//- Return donor cells
const labelList& donorCells() const;
//- Return donor cells processor (to which donor is sent)
const labelList& donorCellsProc() const;
//- Return hole cells
const labelList& holeCells() const;
@ -351,20 +357,52 @@ public:
const word& fieldName
) const;
//- Interpolate to acceptors
//- Expand from master donor to acceptor: injection
template<class Type>
void acceptorMasterData
(
UList<Type>& accF,
const UList<Type>& cellField
) const;
//- Expand from master donor to acceptor: injection
template<class Type>
tmp<Field<Type> > acceptorMasterData
(
const UList<Type>& cellField
) const;
//- Expand from all donors to acceptor: raw donor data
template<class Type>
void acceptorAllData
(
FieldField<Field, Type>& accF,
const UList<Type>& cellField
) const;
//- Expand from master donor to acceptor: injection
template<class Type>
tmp<FieldField<Field, Type> > acceptorAllData
(
const UList<Type>& cellField
) const;
//- Interpolate to acceptors. Supports interpolation weight
// and extended addressing on a per-field basis
template<class Type>
void interpolate
(
Field<Type>& accF,
const Field<Type>& cellField,
UList<Type>& accF,
const UList<Type>& cellField,
const word& fieldName
) const;
//- Interpolate to acceptors
//- Interpolate to acceptors. Supports interpolation weight
// and extended addressing on a per-field basis
template<class Type>
tmp<Field<Type> > interpolate
(
const Field<Type>& cellField,
const UList<Type>& cellField,
const word& fieldName
) const;

View file

@ -30,14 +30,11 @@ License
#include "oversetFvPatchFields.H"
#include "demandDrivenData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::oversetMesh::calcCellClassification() const
{
if (acceptorCellsPtr_ || donorCellsPtr_ || holeCellsPtr_)
if (acceptorCellsPtr_ || donorCellsPtr_ || donorCellsPtr_ || holeCellsPtr_)
{
FatalErrorIn("void oversetMesh::calcCellClassification() const")
<< "Cell classification already calculated"
@ -62,6 +59,9 @@ void Foam::oversetMesh::calcCellClassification() const
donorCellsPtr_ = new labelList(nDonorCells);
labelList& donor = *donorCellsPtr_;
donorCellsProcPtr_ = new labelList(nDonorCells);
labelList& donorProc = *donorCellsProcPtr_;
holeCellsPtr_ = new labelList(nHoleCells);
labelList& hole = *holeCellsPtr_;
@ -105,6 +105,7 @@ void Foam::oversetMesh::calcCellClassification() const
forAll (curDonors, dI)
{
donor[nDonorCells] = curDonors[dI].donorCell();
donorProc[nDonorCells] = curDonors[dI].acceptorProcNo();
nDonorCells++;
}
@ -1126,6 +1127,7 @@ void Foam::oversetMesh::clearOut() const
{
deleteDemandDrivenData(acceptorCellsPtr_);
deleteDemandDrivenData(donorCellsPtr_);
deleteDemandDrivenData(donorCellsProcPtr_);
deleteDemandDrivenData(holeCellsPtr_);
deleteDemandDrivenData(oversetTypesPtr_);
@ -1174,6 +1176,17 @@ const Foam::labelList& Foam::oversetMesh::donorCells() const
}
const Foam::labelList& Foam::oversetMesh::donorCellsProc() const
{
if (!donorCellsProcPtr_)
{
calcCellClassification();
}
return *donorCellsProcPtr_;
}
const Foam::labelList& Foam::oversetMesh::holeCells() const
{
if (!holeCellsPtr_)

View file

@ -27,11 +27,210 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::oversetMesh::acceptorMasterData
(
UList<Type>& accF,
const UList<Type>& cellF
) const
{
// Check sizes
if
(
accF.size() != this->acceptorCells().size()
|| cellF.size() != this->mesh().nCells()
)
{
FatalErrorInFunction
<< "Size of fields does not correspond to interpolation" << nl
<< "Source field size = " << cellF.size()
<< " mesh size = " << this->mesh().nCells()
<< " target field size = " << accF.size()
<< " acceptor list size = " << this->acceptorCells().size()
<< abort(FatalError);
}
// Create a copy of the field to interpolate and distribute its donor data
// across processors.
List<Type> donorField(cellF);
// Note: mapDistribute::distribute changes the size of the field it operates
// on to correspond to the size of received donor data
map().distribute(donorField);
// Get remote donor to local acceptor map. Allows easy indexing of donor
// that came from a (possibly) remote processor via its processor number and
// cell number.
const List<labelField>& remoteDAAddressing =
remoteDonorToLocalAcceptorAddr();
// Note: accF field indexed by number of acceptors (region-wise
// incremental, see oversetMesh::calcCellClassification() for details)
label aI = 0;
// Loop through all regions to account for all acceptors
forAll (regions_, regionI)
{
// Get acceptors for this region
const donorAcceptorList& curAcceptors = regions_[regionI].acceptors();
// Loop through all acceptors of this region
forAll (curAcceptors, regionAccI)
{
// Get necessary acceptor information
const donorAcceptor& curDA = curAcceptors[regionAccI];
// Get necessary donor information
const label donorCellI = curDA.donorCell();
// Get remote donor to local acceptor addressing for this donor
// processor. Note: it is assumed that all donors for an acceptor
// come from the same processor (though it can be remote processor)
const labelField& donorProcAddr =
remoteDAAddressing[curDA.donorProcNo()];
// Master donor contribution
accF[aI] = donorField[donorProcAddr[donorCellI]];
// Increment acceptor index
++aI;
}
}
}
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::oversetMesh::acceptorMasterData
(
const UList<Type>& cellF
) const
{
tmp<Field<Type> > tresult(new Field<Type>(this->acceptorCells().size()));
UList<Type>& result = tresult();
acceptorMasterData(result, cellF);
return tresult;
}
template<class Type>
void Foam::oversetMesh::acceptorAllData
(
FieldField<Field, Type>& accF,
const UList<Type>& cellF
) const
{
// Check sizes
if
(
accF.size() != this->acceptorCells().size()
|| cellF.size() != this->mesh().nCells()
)
{
FatalErrorInFunction
<< "Size of fields does not correspond to interpolation" << nl
<< "Source field size = " << cellF.size()
<< " mesh size = " << this->mesh().nCells()
<< " target field size = " << accF.size()
<< " acceptor list size = " << this->acceptorCells().size()
<< abort(FatalError);
}
// Create a copy of the field to interpolate and distribute its donor data
// across processors.
List<Type> donorField(cellF);
// Note: mapDistribute::distribute changes the size of the field it operates
// on to correspond to the size of received donor data
map().distribute(donorField);
// Get remote donor to local acceptor map. Allows easy indexing of donor
// that came from a (possibly) remote processor via its processor number and
// cell number.
const List<labelField>& remoteDAAddressing =
remoteDonorToLocalAcceptorAddr();
// Note: accF field indexed by number of acceptors (region-wise
// incremental, see oversetMesh::calcCellClassification() for details)
label aI = 0;
// Loop through all regions to account for all acceptors
forAll (regions_, regionI)
{
// Get acceptors for this region
const donorAcceptorList& curAcceptors = regions_[regionI].acceptors();
// Loop through all acceptors of this region
forAll (curAcceptors, regionAccI)
{
// Get necessary acceptor information
const donorAcceptor& curDA = curAcceptors[regionAccI];
// Get necessary donor information
const label donorCellI = curDA.donorCell();
// Neighbouring donor contributions
const donorAcceptor::DynamicLabelList& nbrDonors =
curDA.extendedDonorCells();
// Get remote donor to local acceptor addressing for this donor
// processor. Note: it is assumed that all donors for an acceptor
// come from the same processor (though it can be remote processor)
const labelField& donorProcAddr =
remoteDAAddressing[curDA.donorProcNo()];
// Collect donor contributions for this acceptor
// Note that accF is addressed by aI instead of acceptor cell
// indices. See oversetMesh::calcCellClassification for details
accF[aI].setSize(nbrDonors.size() + 1);
// Master donor contribution (first entry corresponds to
// the weight for master donor)
accF[aI][0] = donorField[donorProcAddr[donorCellI]];
forAll (nbrDonors, nbrDonorI)
{
// Get extended neighbour cell label
const label nbrDonorCellI = nbrDonors[nbrDonorI];
// Note indexing with + 1 offset since the first entry
// is for the master donor and is already accounted for
accF[aI][nbrDonorI + 1] =
donorField[donorProcAddr[nbrDonorCellI]];
}
// Increment acceptor index
++aI;
}
}
}
template<class Type>
Foam::tmp<Foam::FieldField<Foam::Field, Type> >
Foam::oversetMesh::acceptorAllData
(
const UList<Type>& cellF
) const
{
tmp<Field<Type> > tresult
(
new FieldField<Field, Type>(this->acceptorCells().size())
);
FieldField<Field, Type>& result = tresult();
acceptorAllData(result, cellF);
return tresult;
}
template<class Type>
void Foam::oversetMesh::interpolate
(
Field<Type>& accF,
const Field<Type>& cellF,
UList<Type>& accF,
const UList<Type>& cellF,
const word& fieldName
) const
{
@ -42,14 +241,8 @@ void Foam::oversetMesh::interpolate
|| cellF.size() != this->mesh().nCells()
)
{
FatalErrorIn
(
"void oversetMesh::donorToAcceptor\n"
"(\n"
" Field<Type>& accF,\n"
" const Field<Type>& ,cellF\n"
") const"
) << "Size of fields does not correspond to interpolation" << nl
FatalErrorInFunction
<< "Size of fields does not correspond to interpolation" << nl
<< "Source field size = " << cellF.size()
<< " mesh size = " << this->mesh().nCells()
<< " target field size = " << accF.size()
@ -59,7 +252,7 @@ void Foam::oversetMesh::interpolate
// Create a copy of the field to interpolate and distribute its donor data
// across processors.
Field<Type> donorField = cellF;
List<Type> donorField(cellF);
// Note: mapDistribute::distribute changes the size of the field it operates
// on to correspond to the size of received donor data
@ -145,12 +338,12 @@ void Foam::oversetMesh::interpolate
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::oversetMesh::interpolate
(
const Field<Type>& cellF,
const UList<Type>& cellF,
const word& fieldName
) const
{
tmp<Field<Type> > tresult(new Field<Type>(this->acceptorCells().size()));
Field<Type>& result = tresult();
UList<Type>& result = tresult();
interpolate(result, cellF, fieldName);

View file

@ -63,7 +63,8 @@ void Foam::oversetRegion::calcDonorRegions() const
(
"void oversetRegion::calcDonorRegions() const"
) << "Region " << name_ << " specified as the donor "
<< "of itself. List of donor regions: " << donorRegionNames_ << nl
<< "of itself. List of donor regions: "
<< donorRegionNames_ << nl
<< "This is not allowed: check oversetMesh definition"
<< abort(FatalError);
}