mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-12-28 03:37:59 +00:00
912009c458
All remote contributions to interpolation stencils now get added as 'processor' type lduInterfaces. This guarantees a consistent matrix, e.g. initial residual is normalised to 1. Second change is the normalisation of the interpolation discretisation which uses the diagonal from the unmodified equation. This helps GAMG.
689 lines
20 KiB
C
689 lines
20 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2014-2019 OpenCFD Ltd.
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify i
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "dynamicOversetFvMesh.H"
|
|
#include "addToRunTimeSelectionTable.H"
|
|
#include "cellCellStencilObject.H"
|
|
#include "zeroGradientFvPatchFields.H"
|
|
#include "lduPrimitiveProcessorInterface.H"
|
|
#include "globalIndex.H"
|
|
|
|
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
|
|
|
namespace Foam
|
|
{
|
|
defineTypeNameAndDebug(dynamicOversetFvMesh, 0);
|
|
addToRunTimeSelectionTable(dynamicFvMesh, dynamicOversetFvMesh, IOobject);
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
|
|
|
|
bool Foam::dynamicOversetFvMesh::updateAddressing() const
|
|
{
|
|
const cellCellStencilObject& overlap = Stencil::New(*this);
|
|
|
|
// The (processor-local part of the) stencil determines the local
|
|
// faces to add to the matrix. tbd: parallel
|
|
const labelListList& stencil = overlap.cellStencil();
|
|
|
|
// Get the base addressing
|
|
const lduAddressing& baseAddr = dynamicMotionSolverFvMesh::lduAddr();
|
|
|
|
// Add to the base addressing
|
|
labelList lowerAddr;
|
|
labelList upperAddr;
|
|
label nExtraFaces;
|
|
|
|
const globalIndex globalNumbering(baseAddr.size());
|
|
labelListList localFaceCells;
|
|
labelListList remoteFaceCells;
|
|
|
|
labelList globalCellIDs(overlap.cellInterpolationMap().constructSize());
|
|
forAll(baseAddr, cellI)
|
|
{
|
|
globalCellIDs[cellI] = globalNumbering.toGlobal(cellI);
|
|
}
|
|
overlap.cellInterpolationMap().distribute(globalCellIDs);
|
|
|
|
|
|
reverseFaceMap_ = fvMeshPrimitiveLduAddressing::addAddressing
|
|
(
|
|
baseAddr,
|
|
stencil,
|
|
nExtraFaces,
|
|
lowerAddr,
|
|
upperAddr,
|
|
stencilFaces_,
|
|
globalNumbering,
|
|
globalCellIDs,
|
|
localFaceCells,
|
|
remoteFaceCells
|
|
);
|
|
|
|
if (debug)
|
|
{
|
|
Pout<< "dynamicOversetFvMesh::update() : extended addressing from"
|
|
<< " nFaces:" << baseAddr.lowerAddr().size()
|
|
<< " to nFaces:" << lowerAddr.size()
|
|
<< " nExtraFaces:" << nExtraFaces << endl;
|
|
}
|
|
|
|
|
|
// Now for the tricky bits. We want to hand out processor faces according
|
|
// to the localFaceCells/remoteFaceCells. Ultimately we need
|
|
// per entry in stencil:
|
|
// - the patch (or -1 for internal faces)
|
|
// - the face (is either an internal face index or a patch face index)
|
|
|
|
stencilPatches_.setSize(stencilFaces_.size());
|
|
|
|
// Per processor to owner (local)/neighbour (remote)
|
|
List<DynamicList<label>> procOwner(Pstream::nProcs());
|
|
List<DynamicList<label>> dynProcNeighbour(Pstream::nProcs());
|
|
forAll(stencil, celli)
|
|
{
|
|
const labelList& nbrs = stencil[celli];
|
|
stencilPatches_[celli].setSize(nbrs.size(), -1);
|
|
|
|
forAll(nbrs, nbri)
|
|
{
|
|
const label nbrCelli = nbrs[nbri];
|
|
if (stencilFaces_[celli][nbri] == -1)
|
|
{
|
|
label globalNbr = globalCellIDs[nbrCelli];
|
|
label proci = globalNumbering.whichProcID(globalNbr);
|
|
label remoteCelli = globalNumbering.toLocal(proci, globalNbr);
|
|
|
|
// Overwrite the face to be a patch face
|
|
stencilFaces_[celli][nbri] = procOwner[proci].size();
|
|
stencilPatches_[celli][nbri] = proci;
|
|
procOwner[proci].append(celli);
|
|
dynProcNeighbour[proci].append(remoteCelli);
|
|
|
|
//Pout<< "From neighbour proc:" << proci
|
|
// << " allocating patchFace:" << stencilFaces_[celli][nbri]
|
|
// << " to get remote cell " << remoteCelli
|
|
// << " onto local cell " << celli << endl;
|
|
}
|
|
}
|
|
}
|
|
labelListList procNeighbour(dynProcNeighbour.size());
|
|
forAll(procNeighbour, i)
|
|
{
|
|
procNeighbour[i] = std::move(dynProcNeighbour[i]);
|
|
}
|
|
labelListList mySendCells;
|
|
Pstream::exchange<labelList, label>(procNeighbour, mySendCells);
|
|
|
|
label nbri = 0;
|
|
forAll(procOwner, proci)
|
|
{
|
|
if (procOwner[proci].size())
|
|
{
|
|
nbri++;
|
|
}
|
|
if (mySendCells[proci].size())
|
|
{
|
|
nbri++;
|
|
}
|
|
}
|
|
remoteStencilInterfaces_.setSize(nbri);
|
|
nbri = 0;
|
|
|
|
// E.g. if proc1 needs some data from proc2 and proc2 needs some data from
|
|
// proc1. We first want the pair : proc1 receive and proc2 send
|
|
// and then the pair : proc1 send, proc2 receive
|
|
|
|
|
|
labelList procToInterface(Pstream::nProcs(), -1);
|
|
|
|
forAll(procOwner, proci)
|
|
{
|
|
if (proci < Pstream::myProcNo() && procOwner[proci].size())
|
|
{
|
|
//Pout<< "Adding interface " << nbri
|
|
// << " to receive my " << procOwner[proci]
|
|
// << " from " << proci << endl;
|
|
|
|
procToInterface[proci] = nbri;
|
|
remoteStencilInterfaces_.set
|
|
(
|
|
nbri++,
|
|
new lduPrimitiveProcessorInterface
|
|
(
|
|
procOwner[proci],
|
|
Pstream::myProcNo(),
|
|
proci,
|
|
tensorField(0),
|
|
Pstream::msgType()+2
|
|
)
|
|
);
|
|
}
|
|
else if (proci > Pstream::myProcNo() && mySendCells[proci].size())
|
|
{
|
|
//Pout<< "Adding interface " << nbri
|
|
// << " to send my " << mySendCells[proci]
|
|
// << " to " << proci << endl;
|
|
remoteStencilInterfaces_.set
|
|
(
|
|
nbri++,
|
|
new lduPrimitiveProcessorInterface
|
|
(
|
|
mySendCells[proci],
|
|
Pstream::myProcNo(),
|
|
proci,
|
|
tensorField(0),
|
|
Pstream::msgType()+2
|
|
)
|
|
);
|
|
}
|
|
}
|
|
forAll(procOwner, proci)
|
|
{
|
|
if (proci > Pstream::myProcNo() && procOwner[proci].size())
|
|
{
|
|
//Pout<< "Adding interface " << nbri
|
|
// << " to receive my " << procOwner[proci]
|
|
// << " from " << proci << endl;
|
|
procToInterface[proci] = nbri;
|
|
remoteStencilInterfaces_.set
|
|
(
|
|
nbri++,
|
|
new lduPrimitiveProcessorInterface
|
|
(
|
|
procOwner[proci],
|
|
Pstream::myProcNo(),
|
|
proci,
|
|
tensorField(0),
|
|
Pstream::msgType()+2
|
|
)
|
|
);
|
|
}
|
|
else if (proci < Pstream::myProcNo() && mySendCells[proci].size())
|
|
{
|
|
//Pout<< "Adding interface " << nbri
|
|
// << " to send my " << mySendCells[proci]
|
|
// << " to " << proci << endl;
|
|
remoteStencilInterfaces_.set
|
|
(
|
|
nbri++,
|
|
new lduPrimitiveProcessorInterface
|
|
(
|
|
mySendCells[proci],
|
|
Pstream::myProcNo(),
|
|
proci,
|
|
tensorField(0),
|
|
Pstream::msgType()+2
|
|
)
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// Rewrite stencilPatches now we know the actual interface (procToInterface)
|
|
for (auto& patches : stencilPatches_)
|
|
{
|
|
for (auto& interface : patches)
|
|
{
|
|
if (interface != -1)
|
|
{
|
|
interface = procToInterface[interface]+boundary().size();
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Get addressing and interfaces of all interfaces
|
|
|
|
|
|
List<const labelUList*> patchAddr;
|
|
{
|
|
const fvBoundaryMesh& fvp = boundary();
|
|
|
|
patchAddr.setSize(fvp.size() + remoteStencilInterfaces_.size());
|
|
|
|
allInterfaces_ = dynamicMotionSolverFvMesh::interfaces();
|
|
allInterfaces_.setSize(patchAddr.size());
|
|
|
|
forAll(fvp, patchI)
|
|
{
|
|
patchAddr[patchI] = &fvp[patchI].faceCells();
|
|
}
|
|
forAll(remoteStencilInterfaces_, i)
|
|
{
|
|
label patchI = fvp.size()+i;
|
|
const lduPrimitiveProcessorInterface& pp =
|
|
remoteStencilInterfaces_[i];
|
|
|
|
//Pout<< "at patch:" << patchI
|
|
// << " have procPatch:" << pp.type()
|
|
// << " from:" << pp.myProcNo()
|
|
// << " to:" << pp.neighbProcNo()
|
|
// << " with fc:" << pp.faceCells().size() << endl;
|
|
|
|
patchAddr[patchI] = &pp.faceCells();
|
|
allInterfaces_.set(patchI, &pp);
|
|
}
|
|
}
|
|
const lduSchedule ps
|
|
(
|
|
lduPrimitiveMesh::nonBlockingSchedule<processorLduInterface>
|
|
(
|
|
allInterfaces_
|
|
)
|
|
);
|
|
|
|
lduPtr_.reset
|
|
(
|
|
new fvMeshPrimitiveLduAddressing
|
|
(
|
|
nCells(),
|
|
std::move(lowerAddr),
|
|
std::move(upperAddr),
|
|
patchAddr,
|
|
ps
|
|
)
|
|
);
|
|
|
|
|
|
// Check
|
|
if (debug)
|
|
{
|
|
const lduAddressing& addr = lduPtr_(); //this->lduAddr();
|
|
|
|
Pout<< "Adapted addressing:"
|
|
<< " lower:" << addr.lowerAddr().size()
|
|
<< " upper:" << addr.upperAddr().size() << endl;
|
|
|
|
lduInterfacePtrsList iFaces = this->interfaces();
|
|
// Using lduAddressing::patch
|
|
forAll(patchAddr, patchI)
|
|
{
|
|
Pout<< " " << patchI << "\tpatchAddr:"
|
|
<< addr.patchAddr(patchI).size()
|
|
<< endl;
|
|
}
|
|
|
|
// Using interfaces
|
|
Pout<< "iFaces:" << iFaces.size() << endl;
|
|
forAll(iFaces, patchI)
|
|
{
|
|
if (iFaces.set(patchI))
|
|
{
|
|
Pout<< " " << patchI << "\tiFace:" << iFaces[patchI].type()
|
|
<< endl;
|
|
}
|
|
}
|
|
|
|
Pout<< "end of printing." << endl;
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
void Foam::dynamicOversetFvMesh::writeAgglomeration
|
|
(
|
|
const GAMGAgglomeration& agglom
|
|
) const
|
|
{
|
|
labelList cellToCoarse(identity(nCells()));
|
|
labelListList coarseToCell(invertOneToMany(nCells(), cellToCoarse));
|
|
|
|
// Write initial agglomeration
|
|
{
|
|
volScalarField scalarAgglomeration
|
|
(
|
|
IOobject
|
|
(
|
|
"agglomeration",
|
|
this->time().timeName(),
|
|
*this,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
*this,
|
|
dimensionedScalar(dimless, Zero)
|
|
);
|
|
scalarField& fld = scalarAgglomeration.primitiveFieldRef();
|
|
forAll(fld, celli)
|
|
{
|
|
fld[celli] = cellToCoarse[celli];
|
|
}
|
|
fld /= max(fld);
|
|
correctBoundaryConditions
|
|
<
|
|
volScalarField,
|
|
oversetFvPatchField<scalar>
|
|
>(scalarAgglomeration.boundaryFieldRef(), false);
|
|
scalarAgglomeration.write();
|
|
|
|
Info<< "Writing initial cell distribution to "
|
|
<< this->time().timeName() << endl;
|
|
}
|
|
|
|
|
|
for (label level = 0; level < agglom.size(); level++)
|
|
{
|
|
const labelList& addr = agglom.restrictAddressing(level);
|
|
label coarseSize = max(addr)+1;
|
|
|
|
Info<< "Level : " << level << endl
|
|
<< returnReduce(addr.size(), sumOp<label>()) << endl
|
|
<< " current size : "
|
|
<< returnReduce(addr.size(), sumOp<label>()) << endl
|
|
<< " agglomerated size : "
|
|
<< returnReduce(coarseSize, sumOp<label>()) << endl;
|
|
|
|
forAll(addr, fineI)
|
|
{
|
|
const labelList& cellLabels = coarseToCell[fineI];
|
|
forAll(cellLabels, i)
|
|
{
|
|
cellToCoarse[cellLabels[i]] = addr[fineI];
|
|
}
|
|
}
|
|
coarseToCell = invertOneToMany(coarseSize, cellToCoarse);
|
|
|
|
// Write agglomeration
|
|
{
|
|
volScalarField scalarAgglomeration
|
|
(
|
|
IOobject
|
|
(
|
|
"agglomeration_" + Foam::name(level),
|
|
this->time().timeName(),
|
|
*this,
|
|
IOobject::NO_READ,
|
|
IOobject::AUTO_WRITE
|
|
),
|
|
*this,
|
|
dimensionedScalar(dimless, Zero)
|
|
);
|
|
scalarField& fld = scalarAgglomeration.primitiveFieldRef();
|
|
forAll(fld, celli)
|
|
{
|
|
fld[celli] = cellToCoarse[celli];
|
|
}
|
|
//if (normalise)
|
|
//{
|
|
// fld /= max(fld);
|
|
//}
|
|
correctBoundaryConditions
|
|
<
|
|
volScalarField,
|
|
oversetFvPatchField<scalar>
|
|
>(scalarAgglomeration.boundaryFieldRef(), false);
|
|
scalarAgglomeration.write();
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::dynamicOversetFvMesh::dynamicOversetFvMesh(const IOobject& io)
|
|
:
|
|
dynamicMotionSolverFvMesh(io),
|
|
active_(false)
|
|
{
|
|
// Load stencil (but do not update)
|
|
(void)Stencil::New(*this, false);
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
|
|
|
Foam::dynamicOversetFvMesh::~dynamicOversetFvMesh()
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
const Foam::lduAddressing& Foam::dynamicOversetFvMesh::lduAddr() const
|
|
{
|
|
if (!active_)
|
|
{
|
|
return dynamicMotionSolverFvMesh::lduAddr();
|
|
}
|
|
if (lduPtr_.empty())
|
|
{
|
|
// Build extended addressing
|
|
updateAddressing();
|
|
}
|
|
return *lduPtr_;
|
|
}
|
|
|
|
|
|
Foam::lduInterfacePtrsList Foam::dynamicOversetFvMesh::interfaces() const
|
|
{
|
|
if (!active_)
|
|
{
|
|
return dynamicMotionSolverFvMesh::interfaces();
|
|
}
|
|
if (lduPtr_.empty())
|
|
{
|
|
// Build extended addressing
|
|
updateAddressing();
|
|
}
|
|
return allInterfaces_;
|
|
}
|
|
|
|
|
|
const Foam::fvMeshPrimitiveLduAddressing&
|
|
Foam::dynamicOversetFvMesh::primitiveLduAddr() const
|
|
{
|
|
if (lduPtr_.empty())
|
|
{
|
|
FatalErrorInFunction
|
|
<< "Extended addressing not allocated" << abort(FatalError);
|
|
}
|
|
|
|
return *lduPtr_;
|
|
}
|
|
|
|
|
|
bool Foam::dynamicOversetFvMesh::update()
|
|
{
|
|
if (dynamicMotionSolverFvMesh::update())
|
|
{
|
|
// Calculate the local extra faces for the interpolation. Note: could
|
|
// let demand-driven lduAddr() trigger it but just to make sure.
|
|
updateAddressing();
|
|
|
|
// Addressing and/or weights have changed. Make interpolated cells
|
|
// up to date with donors
|
|
interpolateFields();
|
|
|
|
return true;
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
|
|
Foam::word Foam::dynamicOversetFvMesh::baseName(const word& name)
|
|
{
|
|
if (name.endsWith("_0"))
|
|
{
|
|
return baseName(name.substr(0, name.size()-2));
|
|
}
|
|
|
|
return name;
|
|
}
|
|
|
|
|
|
bool Foam::dynamicOversetFvMesh::interpolateFields()
|
|
{
|
|
// Add the stencil suppression list
|
|
wordHashSet suppressed(Stencil::New(*this).nonInterpolatedFields());
|
|
|
|
// Use whatever the solver has set up as suppression list
|
|
const dictionary* dictPtr
|
|
(
|
|
this->schemesDict().findDict("oversetInterpolationSuppressed")
|
|
);
|
|
if (dictPtr)
|
|
{
|
|
suppressed.insert(dictPtr->toc());
|
|
}
|
|
|
|
interpolate<volScalarField>(suppressed);
|
|
interpolate<volVectorField>(suppressed);
|
|
interpolate<volSphericalTensorField>(suppressed);
|
|
interpolate<volSymmTensorField>(suppressed);
|
|
interpolate<volTensorField>(suppressed);
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
|
|
bool Foam::dynamicOversetFvMesh::writeObject
|
|
(
|
|
IOstream::streamFormat fmt,
|
|
IOstream::versionNumber ver,
|
|
IOstream::compressionType cmp,
|
|
const bool valid
|
|
) const
|
|
{
|
|
bool ok = dynamicMotionSolverFvMesh::writeObject(fmt, ver, cmp, valid);
|
|
|
|
// For postprocessing : write cellTypes and zoneID
|
|
{
|
|
const cellCellStencilObject& overlap = Stencil::New(*this);
|
|
|
|
const labelUList& cellTypes = overlap.cellTypes();
|
|
|
|
volScalarField volTypes
|
|
(
|
|
IOobject
|
|
(
|
|
"cellTypes",
|
|
this->time().timeName(),
|
|
*this,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
*this,
|
|
dimensionedScalar(dimless, Zero),
|
|
zeroGradientFvPatchScalarField::typeName
|
|
);
|
|
|
|
forAll(volTypes.internalField(), cellI)
|
|
{
|
|
volTypes[cellI] = cellTypes[cellI];
|
|
}
|
|
volTypes.correctBoundaryConditions();
|
|
volTypes.writeObject(fmt, ver, cmp, valid);
|
|
}
|
|
{
|
|
volScalarField volZoneID
|
|
(
|
|
IOobject
|
|
(
|
|
"zoneID",
|
|
this->time().timeName(),
|
|
*this,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
*this,
|
|
dimensionedScalar(dimless, Zero),
|
|
zeroGradientFvPatchScalarField::typeName
|
|
);
|
|
|
|
const cellCellStencilObject& overlap = Stencil::New(*this);
|
|
const labelIOList& zoneID = overlap.zoneID();
|
|
|
|
forAll(zoneID, cellI)
|
|
{
|
|
volZoneID[cellI] = zoneID[cellI];
|
|
}
|
|
volZoneID.correctBoundaryConditions();
|
|
volZoneID.writeObject(fmt, ver, cmp, valid);
|
|
}
|
|
if (debug)
|
|
{
|
|
const cellCellStencilObject& overlap = Stencil::New(*this);
|
|
const labelIOList& zoneID = overlap.zoneID();
|
|
const labelListList& cellStencil = overlap.cellStencil();
|
|
|
|
labelList donorZoneID(zoneID);
|
|
overlap.cellInterpolationMap().distribute(donorZoneID);
|
|
|
|
forAll(cellStencil, cellI)
|
|
{
|
|
const labelList& stencil = cellStencil[cellI];
|
|
if (stencil.size())
|
|
{
|
|
donorZoneID[cellI] = zoneID[stencil[0]];
|
|
for (label i = 1; i < stencil.size(); i++)
|
|
{
|
|
if (zoneID[stencil[i]] != donorZoneID[cellI])
|
|
{
|
|
WarningInFunction << "Mixed donor meshes for cell "
|
|
<< cellI << " at " << C()[cellI]
|
|
<< " donors:" << UIndirectList<point>(C(), stencil)
|
|
<< endl;
|
|
donorZoneID[cellI] = -2;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
volScalarField volDonorZoneID
|
|
(
|
|
IOobject
|
|
(
|
|
"donorZoneID",
|
|
this->time().timeName(),
|
|
*this,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE,
|
|
false
|
|
),
|
|
*this,
|
|
dimensionedScalar("minOne", dimless, scalar(-1)),
|
|
zeroGradientFvPatchScalarField::typeName
|
|
);
|
|
forAll(donorZoneID, celli)
|
|
{
|
|
volDonorZoneID[celli] = donorZoneID[celli];
|
|
}
|
|
//- Do not correctBoundaryConditions since re-interpolates!
|
|
//volDonorZoneID.correctBoundaryConditions();
|
|
volDonorZoneID.writeObject(fmt, ver, cmp, valid);
|
|
}
|
|
|
|
return ok;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|