songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: geoVoF module has the capability to run AMR with load balancing

Browse Source 
code style and quality improvements
renamed recon::centre to interfaceCentre.{groupName}
ranmed recon::normal to interfaceNormal.{groupName}
centre and normal field are not written by default
This commit is contained in:
HenningScheufler 2021-05-12 14:32:53 +02:00 committed by Andrew Heather
parent 7e4cd55882
commit 9dada5f3f2
33 changed files with 481 additions and 484 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
View File

@ -91,10 +91,11 @@ Foam::zoneDistribute::coupledFacesPatch() const
Foam::zoneDistribute::zoneDistribute(const fvMesh& mesh)
:
MeshObject<fvMesh, Foam::UpdateableMeshObject, zoneDistribute>(mesh),
stencil_(mesh),
MeshObject<fvMesh, Foam::TopologicalMeshObject, zoneDistribute>(mesh),
coupledBoundaryPoints_(coupledFacesPatch()().meshPoints()),
send_(Pstream::nProcs())
send_(Pstream::nProcs()),
stencil_(zoneCPCStencil::New(mesh)),
gblIdx_(stencil_.globalNumbering())
{
}
@ -103,15 +104,11 @@ Foam::zoneDistribute::zoneDistribute(const fvMesh& mesh)
Foam::zoneDistribute& Foam::zoneDistribute::New(const fvMesh& mesh)
{
zoneDistribute* ptr = mesh.thisDb().getObjectPtr<zoneDistribute>
(
zoneDistribute::typeName
);
auto* ptr = mesh.thisDb().getObjectPtr<zoneDistribute>("zoneDistribute");
if (!ptr)
{
ptr = new zoneDistribute(mesh);
regIOobject::store(ptr);
}
@ -123,24 +120,22 @@ Foam::zoneDistribute& Foam::zoneDistribute::New(const fvMesh& mesh)
void Foam::zoneDistribute::updateStencil(const boolList& zone)
{
stencil_.updateStencil(zone);
zoneCPCStencil::New(mesh_).updateStencil(zone);
}
void Foam::zoneDistribute::setUpCommforZone
(
const boolList& zone,
bool updateStencil
)
void Foam::zoneDistribute::setUpCommforZone(const boolList& zone,bool updateStencil)
{
zoneCPCStencil& stencil = zoneCPCStencil::New(mesh_);
if (updateStencil)
{
stencil_.updateStencil(zone);
stencil.updateStencil(zone);
}
const labelHashSet comms = stencil_.needsComm();
const labelHashSet comms = stencil.needsComm();
List<labelHashSet> needed_(Pstream::nProcs());
List<labelHashSet> needed(Pstream::nProcs());
if (Pstream::parRun())
{
@ -150,11 +145,10 @@ void Foam::zoneDistribute::setUpCommforZone
{
for (const label gblIdx : stencil_[celli])
{
if (!globalNumbering().isLocal(gblIdx))
if (!gblIdx_.isLocal(gblIdx))
{
const label procID =
globalNumbering().whichProcID(gblIdx);
needed_[procID].insert(gblIdx);
const label procID = gblIdx_.whichProcID (gblIdx);
needed[procID].insert(gblIdx);
}
}
}
@ -170,7 +164,7 @@ void Foam::zoneDistribute::setUpCommforZone
// Put data into send buffer
UOPstream toDomain(domain, pBufs);
toDomain << needed_[domain];
toDomain << needed[domain];
}
}
@ -193,13 +187,6 @@ void Foam::zoneDistribute::setUpCommforZone
}
void Foam::zoneDistribute::updateMesh(const mapPolyMesh& mpm)
{
if (mesh_.topoChanging())
{
coupledBoundaryPoints_ = coupledFacesPatch()().meshPoints();
}
}
// ************************************************************************* //

25
src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
View File

@ -64,12 +64,10 @@ namespace Foam
class zoneDistribute
:
public MeshObject<fvMesh, UpdateableMeshObject, zoneDistribute>
public MeshObject<fvMesh, TopologicalMeshObject, zoneDistribute>
{
// Private Data
//- cell-point-cell stencil elements are in global addressing
zoneCPCStencil stencil_;
//- labels of the points on coupled patches
labelList coupledBoundaryPoints_;
@ -80,6 +78,10 @@ class zoneDistribute
//- Return patch of all coupled faces.
autoPtr<indirectPrimitivePatch> coupledFacesPatch() const;
zoneCPCStencil& stencil_;
const globalIndex& gblIdx_;
//- Gives patchNumber and patchFaceNumber for a given
//- Geometric volume field
@ -118,6 +120,11 @@ public:
static zoneDistribute& New(const fvMesh&);
//- Destructor
virtual ~zoneDistribute() = default;
// Member Functions
//- Update stencil with boolList the size has to match mesh nCells
@ -127,15 +134,15 @@ public:
void updateStencil(const boolList& zone);
//- Stencil reference
const labelListList& getStencil()
const labelListList& getStencil() noexcept
{
return stencil_;
}
//- Addressing reference
const globalIndex& globalNumbering() const
const globalIndex& globalNumbering() const noexcept
{
return stencil_.globalNumbering();
return gblIdx_;
}
//- Gives patchNumber and patchFaceNumber for a given
@ -164,13 +171,7 @@ public:
const GeometricField<Type, fvPatchField, volMesh>& phi
);
virtual void updateMesh(const mapPolyMesh& mpm);
virtual bool movePoints()
{
// do nothing
return false;
}
};

13
src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
View File

@ -43,11 +43,8 @@ Type Foam::zoneDistribute::getLocalValue
{
return phi[localIdx];
}
else
{
return faceValue(phi,localIdx);
}
return faceValue(phi,localIdx);
}
@ -88,9 +85,9 @@ Type Foam::zoneDistribute::getValue
const label gblIdx
) const
{
if (globalNumbering().isLocal(gblIdx))
if (gblIdx_.isLocal(gblIdx))
{
const label idx = globalNumbering().toLocal(gblIdx);
const label idx = gblIdx_.toLocal(gblIdx);
return getLocalValue(phi,idx);
}
else // from other proc
@ -115,7 +112,6 @@ Foam::Map<Foam::Field<Type>> Foam::zoneDistribute::getFields
<< "do not match. Did the mesh change?"
<< exit(FatalError);
return Map<Field<Type>>();
}
@ -160,7 +156,6 @@ Foam::Map<Type> Foam::zoneDistribute::getDatafromOtherProc
<< "do not match. Did the mesh change?"
<< exit(FatalError);
return Map<Type>();
}
@ -177,7 +172,7 @@ Foam::Map<Type> Foam::zoneDistribute::getDatafromOtherProc
sendValues[domaini].insert
(
sendIdx,
getLocalValue(phi,globalNumbering().toLocal(sendIdx))
getLocalValue(phi,gblIdx_.toLocal(sendIdx))
);
}
}

29
src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.C
View File

@ -123,6 +123,7 @@ void Foam::zoneCPCStencil::calcPointBoundaryData
Foam::zoneCPCStencil::zoneCPCStencil(const fvMesh& mesh)
:
MeshObject<fvMesh, Foam::TopologicalMeshObject, zoneCPCStencil>(mesh),
zoneCellStencils(mesh),
nonEmptyBoundaryPoints_(nonEmptyFacesPatch()().meshPoints()),
uptodate_(mesh.nCells(), false)
@ -131,6 +132,19 @@ Foam::zoneCPCStencil::zoneCPCStencil(const fvMesh& mesh)
validBoundaryFaces(isValidBFace_);
}
Foam::zoneCPCStencil& Foam::zoneCPCStencil::New(const fvMesh& mesh)
{
auto* ptr = mesh.thisDb().getObjectPtr<zoneCPCStencil>("zoneCPCStencil");
if (!ptr)
{
ptr = new zoneCPCStencil(mesh);
regIOobject::store(ptr);
}
return *ptr;
}
void Foam::zoneCPCStencil::calculateStencil
(
@ -224,19 +238,4 @@ void Foam::zoneCPCStencil::calculateStencil
}
void Foam::zoneCPCStencil::updateMesh(const mapPolyMesh& mpm)
{
if (mesh_.topoChanging())
{
// resize map and globalIndex
zoneCellStencils::updateMesh(mpm);
nonEmptyBoundaryPoints_ = nonEmptyFacesPatch()().meshPoints();
uptodate_.resize(mesh_.nCells());
uptodate_ = false;
validBoundaryFaces(isValidBFace_);
}
}
// ************************************************************************* //

18
src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCPCStencil.H
View File

@ -58,6 +58,12 @@ namespace Foam
class zoneCPCStencil
:
public MeshObject
<
fvMesh,
TopologicalMeshObject,
zoneCPCStencil
>,
public zoneCellStencils
{
// Private Data
@ -97,6 +103,11 @@ class zoneCPCStencil
labelListList& globalCellCells
);
//- No copy construct
zoneCPCStencil(const zoneCPCStencil&) = delete;
//- No copy assignment
void operator=(const zoneCPCStencil&) = delete;
public:
@ -109,13 +120,10 @@ public:
//- Construct from all cells and boundary faces
explicit zoneCPCStencil(const fvMesh&);
// Selectors
// Member Functions
static zoneCPCStencil& New(const fvMesh&);
//- Calculates per cell the neighbour data
//- (= cell or boundary in global numbering).
//- First element is always cell itself!
virtual void updateMesh(const mapPolyMesh& mpm);
};

54
src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.C
View File

@ -43,7 +43,7 @@ namespace Foam
Foam::autoPtr<Foam::indirectPrimitivePatch>
Foam::zoneCellStencils::nonEmptyFacesPatch() const
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const polyBoundaryMesh& patches = meshRef_.boundaryMesh();
label nNonEmpty = 0;
@ -74,10 +74,10 @@ Foam::zoneCellStencils::nonEmptyFacesPatch() const
(
IndirectList<face>
(
mesh_.faces(),
meshRef_.faces(),
nonEmptyFaces
),
mesh_.points()
meshRef_.points()
);
}
@ -85,7 +85,7 @@ Foam::zoneCellStencils::nonEmptyFacesPatch() const
Foam::autoPtr<Foam::indirectPrimitivePatch>
Foam::zoneCellStencils::allCoupledFacesPatch() const
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const polyBoundaryMesh& patches = meshRef_.boundaryMesh();
label nCoupled = 0;
@ -116,19 +116,19 @@ Foam::zoneCellStencils::allCoupledFacesPatch() const
(
IndirectList<face>
(
mesh_.faces(),
meshRef_.faces(),
coupledFaces
),
mesh_.points()
meshRef_.points()
);
}
void Foam::zoneCellStencils::validBoundaryFaces(boolList& isValidBFace) const
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
const polyBoundaryMesh& patches = meshRef_.boundaryMesh();
isValidBFace.setSize(mesh_.nBoundaryFaces());
isValidBFace.setSize(meshRef_.nBoundaryFaces());
isValidBFace = true;
@ -136,7 +136,7 @@ void Foam::zoneCellStencils::validBoundaryFaces(boolList& isValidBFace) const
{
if (pp.coupled() || isA<emptyPolyPatch>(pp))
{
label bFacei = pp.start()-mesh_.nInternalFaces();
label bFacei = pp.start() - meshRef_.nInternalFaces();
forAll(pp, i)
{
isValidBFace[bFacei++] = false;
@ -190,22 +190,20 @@ void Foam::zoneCellStencils::insertFaceCells
labelHashSet& globals
) const
{
const labelList& own = mesh_.faceOwner();
const labelList& nei = mesh_.faceNeighbour();
const labelList& own = meshRef_.faceOwner();
const labelList& nei = meshRef_.faceNeighbour();
forAll(faceLabels, i)
for (const label facei : faceLabels)
{
label facei = faceLabels[i];
label globalOwn = globalNumbering().toGlobal(own[facei]);
const label globalOwn = globalNumbering().toGlobal(own[facei]);
if (globalOwn != exclude0 && globalOwn != exclude1)
{
globals.insert(globalOwn);
}
if (mesh_.isInternalFace(facei))
if (meshRef_.isInternalFace(facei))
{
label globalNei = globalNumbering().toGlobal(nei[facei]);
const label globalNei = globalNumbering().toGlobal(nei[facei]);
if (globalNei != exclude0 && globalNei != exclude1)
{
globals.insert(globalNei);
@ -213,13 +211,14 @@ void Foam::zoneCellStencils::insertFaceCells
}
else
{
const label bFacei = facei-mesh_.nInternalFaces();
const label bFacei = facei - meshRef_.nInternalFaces();
if (isValidBFace[bFacei])
{
label globalI = globalNumbering().toGlobal
(
mesh_.nCells() + bFacei
meshRef_.nCells()
+ bFacei
);
if (globalI != exclude0 && globalI != exclude1)
@ -258,25 +257,12 @@ Foam::labelList Foam::zoneCellStencils::calcFaceCells
Foam::zoneCellStencils::zoneCellStencils(const fvMesh& mesh)
:
MeshObject<fvMesh, Foam::UpdateableMeshObject, zoneCellStencils>(mesh),
labelListList(mesh.nCells()),
meshRef_(mesh),
needComm_(),
globalNumbering_(mesh_.nCells() + mesh_.nBoundaryFaces())
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::zoneCellStencils::updateMesh(const mapPolyMesh& mpm)
globalNumbering_(meshRef_.nCells() + meshRef_.nBoundaryFaces())
{
if (mesh_.topoChanging())
{
globalNumbering_ =
globalIndex(mesh_.nCells() + mesh_.nBoundaryFaces());
static_cast<labelListList&>(*this) = labelListList(mesh_.nCells());
needComm_.clear();
}
}

38
src/finiteVolume/fvMesh/zoneDistribute/zoneStencils/zoneCellStencils.H
View File

@ -58,15 +58,19 @@ namespace Foam
class zoneCellStencils
:
public MeshObject<fvMesh, UpdateableMeshObject, zoneCellStencils>,
public labelListList
{
protected:
// Protected Members
// Protected Data
//- const reference to fvMesh
const fvMesh& meshRef_;
//- cells requiring processor communciation
labelHashSet needComm_;
//- Global numbering for cells and boundary faces
globalIndex globalNumbering_;
@ -114,6 +118,12 @@ protected:
) = 0;
//- No copy construct
zoneCellStencils(const zoneCellStencils&) = delete;
//- No copy assignment
void operator=(const zoneCellStencils&) = delete;
public:
// Declare name of the class and its debug switch
@ -125,12 +135,15 @@ public:
//- Construct from all cells and boundary faces
explicit zoneCellStencils(const fvMesh&);
// Member Functions
//- Calculates per cell the neighbour data
// (= cell or boundary in global numbering).
// First element is always cell itself!
//- Destructor
virtual ~zoneCellStencils() = default;
// Member Functions
void updateStencil
(
const boolList& zone
@ -139,23 +152,20 @@ public:
calculateStencil(zone,*this);
}
const labelHashSet& needsComm()
const labelHashSet& needsComm() noexcept
{
return needComm_;
}
//- Global numbering for cells and boundary faces
const globalIndex& globalNumbering() const
const polyMesh& mesh() const noexcept
{
return globalNumbering_;
return meshRef_;
}
virtual void updateMesh(const mapPolyMesh& mpm);
virtual bool movePoints()
//- Global numbering for cells and boundary faces
const globalIndex& globalNumbering() const noexcept
{
// Do nothing
return false;
return globalNumbering_;
}
};

60
src/functionObjects/field/setFlow/setFlow.C
View File

@ -291,48 +291,42 @@ bool Foam::functionObjects::setFlow::execute()
Uc.replace(vector::Z, sin(2*pi*x)*sqr(sin(pi*z)));
U = U & R_;
U.correctBoundaryConditions();
// Calculating incompressible flux based on stream function
// Calculating phi
// Note: R_ rotation not implemented in phi calculation
const scalarField xp(mesh_.points().component(0) - origin_[0]);
const scalarField yp(mesh_.points().component(1) - origin_[1]);
const scalarField zp(mesh_.points().component(2) - origin_[2]);
const scalarField psi((1.0/pi)*sqr(sin(pi*xp))*sqr(sin(pi*zp)));
const vectorField Cf(mesh_.Cf().primitiveField() - origin_);
const scalarField Xf(Cf.component(vector::X));
const scalarField Yf(Cf.component(vector::Y));
const scalarField Zf(Cf.component(vector::Z));
vectorField Uf(Xf.size());
Uf.replace(vector::X, -sin(2*pi*Zf)*sqr(sin(pi*Xf)));
Uf.replace(vector::Y, scalar(0));
Uf.replace(vector::Z, sin(2*pi*Xf)*sqr(sin(pi*Zf)));
scalarField& phic = phi.primitiveFieldRef();
forAll(phic, fi)
{
phic[fi] = 0;
const face& f = mesh_.faces()[fi];
const label nPoints = f.size();
forAll(f, fpi)
{
const label p1 = f[fpi];
const label p2 = f[(fpi + 1) % nPoints];
phic[fi] += 0.5*(psi[p1] + psi[p2])*(yp[p2] - yp[p1]);
}
}
const vectorField& Sfc = mesh_.Sf().primitiveField();
phic = Uf & Sfc;
surfaceScalarField::Boundary& phibf = phi.boundaryFieldRef();
const surfaceVectorField::Boundary& Sfbf =
mesh_.Sf().boundaryField();
const surfaceVectorField::Boundary& Cfbf =
mesh_.Cf().boundaryField();
forAll(phibf, patchi)
{
scalarField& phif = phibf[patchi];
const label start = mesh_.boundaryMesh()[patchi].start();
forAll(phif, fi)
{
phif[fi] = 0;
const face& f = mesh_.faces()[start + fi];
const label nPoints = f.size();
forAll(f, fpi)
{
const label p1 = f[fpi];
const label p2 = f[(fpi + 1) % nPoints];
phif[fi] += 0.5*(psi[p1] + psi[p2])*(yp[p2] - yp[p1]);
}
}
const vectorField& Sff = Sfbf[patchi];
const vectorField& Cff = Cfbf[patchi];
const scalarField xf(Cff.component(vector::X));
const scalarField yf(Cff.component(vector::Y));
const scalarField zf(Cff.component(vector::Z));
vectorField Ufb(xf.size());
Ufb.replace(vector::X, -sin(2*pi*zf)*sqr(sin(pi*xf)));
Ufb.replace(vector::Y, scalar(0));
Ufb.replace(vector::Z, sin(2*pi*xf)*sqr(sin(pi*zf)));
phif = Ufb & Sff;
}
break;

115
src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C
View File

@ -38,6 +38,9 @@ License
#include "cellSet.H"
#include "meshTools.H"
#include "OBJstream.H"
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -88,6 +91,7 @@ Foam::isoAdvection::isoAdvection
dimensionedScalar(dimVol/dimTime, Zero)
),
advectionTime_(0),
timeIndex_(-1),
// Tolerances and solution controls
nAlphaBounds_(dict_.getOrDefault<label>("nAlphaBounds", 3)),
@ -126,28 +130,69 @@ Foam::isoAdvection::isoAdvection
mesh_.cells();
// Get boundary mesh and resize the list for parallel comms
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
surfaceCellFacesOnProcPatches_.resize(patches.size());
// Append all processor patch labels to the list
forAll(patches, patchi)
{
if
(
isA<processorPolyPatch>(patches[patchi])
&& patches[patchi].size() > 0
)
{
procPatchLabels_.append(patchi);
}
}
setProcessorPatches();
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::isoAdvection::setProcessorPatches()
{
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
surfaceCellFacesOnProcPatches_.clear();
surfaceCellFacesOnProcPatches_.resize(patches.size());
// Append all processor patch labels to the list
procPatchLabels_.clear();
forAll(patches, patchi)
{
if
(
isA<processorPolyPatch>(patches[patchi])
&& !patches[patchi].empty()
)
{
procPatchLabels_.append(patchi);
}
}
}
void Foam::isoAdvection::extendMarkedCells
(
bitSet& markedCell
) const
{
// Mark faces using any marked cell
bitSet markedFace(mesh_.nFaces());
for (const label celli : markedCell)
{
markedFace.set(mesh_.cells()[celli]); // set multiple faces
}
syncTools::syncFaceList(mesh_, markedFace, orEqOp<unsigned int>());
// Update cells using any markedFace
for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
{
if (markedFace.test(facei))
{
markedCell.set(mesh_.faceOwner()[facei]);
markedCell.set(mesh_.faceNeighbour()[facei]);
}
}
for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
{
if (markedFace.test(facei))
{
markedCell.set(mesh_.faceOwner()[facei]);
}
}
}
void Foam::isoAdvection::timeIntegratedFlux()
{
// Get time step
@ -196,8 +241,7 @@ void Foam::isoAdvection::timeIntegratedFlux()
// Cell is cut
const point x0 = surf_->centre()[celli];
vector n0 = -surf_->normal()[celli];
n0 /= (mag(n0));
const vector n0(normalised(-surf_->normal()[celli]));
// Get the speed of the isoface by interpolating velocity and
// dotting it with isoface unit normal
@ -212,10 +256,8 @@ void Foam::isoAdvection::timeIntegratedFlux()
// Note: looping over all cell faces - in reduced-D, some of
// these faces will be on empty patches
const cell& celliFaces = cellFaces[celli];
forAll(celliFaces, fi)
for (const label facei : celliFaces)
{
const label facei = celliFaces[fi];
if (mesh_.isInternalFace(facei))
{
bool isDownwindFace = false;
@ -279,7 +321,7 @@ void Foam::isoAdvection::timeIntegratedFlux()
const label patchi = boundaryMesh.patchID()[facei - nInternalFaces];
const label start = boundaryMesh[patchi].start();
if (phib[patchi].size())
if (!phib[patchi].empty())
{
const label patchFacei = facei - start;
const scalar phiP = phib[patchi][patchFacei];
@ -332,10 +374,9 @@ void Foam::isoAdvection::setDownwindFaces
downwindFaces.clear();
// Check all faces of the cell
forAll(c, fi)
for (const label facei: c)
{
// Get face and corresponding flux
const label facei = c[fi];
const scalar phi = faceValue(phi_, facei);
if (own[facei] == celli)
@ -369,9 +410,8 @@ Foam::scalar Foam::isoAdvection::netFlux
// Get mesh data
const labelList& own = mesh_.faceOwner();
forAll(c, fi)
for (const label facei : c)
{
const label facei = c[fi];
const scalar dVff = faceValue(dVf, facei);
if (own[facei] == celli)
@ -403,10 +443,8 @@ Foam::DynamicList<Foam::label> Foam::isoAdvection::syncProcPatches
PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
// Send
forAll(procPatchLabels_, i)
for (const label patchi : procPatchLabels_)
{
const label patchi = procPatchLabels_[i];
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(patches[patchi]);
@ -429,10 +467,8 @@ Foam::DynamicList<Foam::label> Foam::isoAdvection::syncProcPatches
// Receive and combine
forAll(procPatchLabels_, patchLabeli)
for (const label patchi : procPatchLabels_)
{
const label patchi = procPatchLabels_[patchLabeli];
const processorPolyPatch& procPatch =
refCast<const processorPolyPatch>(patches[patchi]);
@ -466,7 +502,7 @@ Foam::DynamicList<Foam::label> Foam::isoAdvection::syncProcPatches
{
const label facei = faceIDs[i];
localFlux[facei] = - nbrdVfs[i];
if (debug && mag(localFlux[facei] + nbrdVfs[i]) > 10*SMALL)
if (debug && mag(localFlux[facei] + nbrdVfs[i]) > ROOTVSMALL)
{
Pout<< "localFlux[facei] = " << localFlux[facei]
<< " and nbrdVfs[i] = " << nbrdVfs[i]
@ -505,7 +541,7 @@ void Foam::isoAdvection::checkIfOnProcPatch(const label facei)
const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
if (isA<processorPolyPatch>(pbm[patchi]) && pbm[patchi].size())
if (isA<processorPolyPatch>(pbm[patchi]) && !pbm[patchi].empty())
{
const label patchFacei = pbm[patchi].whichFace(facei);
surfaceCellFacesOnProcPatches_[patchi].append(patchFacei);
@ -514,8 +550,6 @@ void Foam::isoAdvection::checkIfOnProcPatch(const label facei)
}
void Foam::isoAdvection::applyBruteForceBounding()
{
bool alpha1Changed = false;
@ -568,6 +602,7 @@ void Foam::isoAdvection::writeSurfaceCells() const
}
}
void Foam::isoAdvection::writeIsoFaces
(
const DynamicList<List<point>>& faces
@ -603,10 +638,8 @@ void Foam::isoAdvection::writeIsoFaces
const DynamicList<List<point>>& procFacePts =
allProcFaces[proci];
forAll(procFacePts, i)
for (const List<point>& facePts : procFacePts)
{
const List<point>& facePts = procFacePts[i];
if (facePts.size() != f.size())
{
f = face(identity(facePts.size()));
@ -625,10 +658,8 @@ void Foam::isoAdvection::writeIsoFaces
<< os.name() << nl << endl;
face f;
forAll(faces, i)
for (const List<point>& facePts : faces)
{
const List<point>& facePts = faces[i];
if (facePts.size() != f.size())
{
f = face(identity(facePts.size()));

24
src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.H
View File

@ -115,6 +115,8 @@ class isoAdvection
//- Time spent performing interface advection
scalar advectionTime_;
//- store timeIndex
label timeIndex_;
// Switches and tolerances. Tolerances need to go into toleranceSwitches
@ -179,6 +181,11 @@ class isoAdvection
// Advection functions
//- Extend markedCell with cell-face-cell.
void extendMarkedCells(bitSet& markedCell) const;
void setProcessorPatches();
//- For each face calculate volumetric face transport during dt
void timeIntegratedFlux();
@ -202,7 +209,7 @@ class isoAdvection
template < class SpType, class SuType >
void boundFlux
(
const scalarField& alpha1,
const bitSet& nextToInterface,
surfaceScalarField& dVfcorrectionValues,
DynamicLabelList& correctedFaces,
const SpType& Sp,
@ -274,6 +281,8 @@ class isoAdvection
// list of surface cell faces on processor patches
void checkIfOnProcPatch(const label facei);
//- Apply the bounding based on user inputs
void applyBruteForceBounding();
public:
@ -303,25 +312,22 @@ public:
template < class SpType, class SuType >
void advect(const SpType& Sp, const SuType& Su);
//- Apply the bounding based on user inputs
void applyBruteForceBounding();
// Access functions
//- Return reconstructionSchemes
reconstructionSchemes& surf()
reconstructionSchemes& surf() noexcept
{
return surf_();
}
//- Return alpha field
const volScalarField& alpha() const
const volScalarField& alpha() const noexcept
{
return alpha1_;
}
//- Return the controls dictionary
const dictionary& dict() const
const dictionary& dict() const noexcept
{
return dict_;
}
@ -365,13 +371,13 @@ public:
}
//- reference to alphaPhi
const surfaceScalarField& alphaPhi() const
const surfaceScalarField& alphaPhi() const noexcept
{
return alphaPhi_;
}
//- time spend in the advection step
scalar advectionTime() const
scalar advectionTime() const noexcept
{
return advectionTime_;
}

60
src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvectionTemplates.C
View File

@ -120,7 +120,7 @@ void Foam::isoAdvection::limitFluxes
{
DebugInFunction << endl;
const scalar aTol = 1.0e-12; // Note: tolerances
const scalar aTol = 100*SMALL; // Note: tolerances
scalar maxAlphaMinus1 = gMax(alpha1_) - 1; // max(alphaNew - 1);
scalar minAlpha = gMin(alpha1_); // min(alphaNew);
const label nOvershoots = 20; // sum(pos0(alphaNew - 1 - aTol));
@ -134,6 +134,10 @@ void Foam::isoAdvection::limitFluxes
surfaceScalarField dVfcorrectionValues("dVfcorrectionValues", dVf_*0.0);
bitSet needBounding(mesh_.nCells(),false);
needBounding.set(surfCells_);
extendMarkedCells(needBounding);
// Loop number of bounding steps
for (label n = 0; n < nAlphaBounds_; n++)
@ -143,8 +147,8 @@ void Foam::isoAdvection::limitFluxes
DebugInfo << "boundAlpha... " << endl;
DynamicList<label> correctedFaces(3*nOvershoots);
dVfcorrectionValues = dimensionedScalar(dimVolume, Zero);
boundFlux(alpha1In_, dVfcorrectionValues, correctedFaces,Sp,Su);
dVfcorrectionValues = dimensionedScalar("0",dimVolume,0.0);
boundFlux(needBounding, dVfcorrectionValues, correctedFaces,Sp,Su);
correctedFaces.append
(
@ -152,25 +156,24 @@ void Foam::isoAdvection::limitFluxes
);
labelHashSet alreadyUpdated;
forAll(correctedFaces, fi)
for (const label facei : correctedFaces)
{
label facei = correctedFaces[fi];
if (alreadyUpdated.insert(facei))
{
checkIfOnProcPatch(facei);
const label own = owner[facei];
alpha1_[own] +=
-faceValue(dVfcorrectionValues, facei)/mesh_.V()[own];
alpha1_[own] -=
faceValue(dVfcorrectionValues, facei)/mesh_.V()[own];
if (mesh_.isInternalFace(facei))
{
const label nei = neighbour[facei];
alpha1_[nei] -=
-faceValue(dVfcorrectionValues, facei)/mesh_.V()[nei];
alpha1_[nei] +=
faceValue(dVfcorrectionValues, facei)/mesh_.V()[nei];
}
// Change to treat boundaries consistently
scalar corrVf =
const scalar corrVf =
faceValue(dVf_, facei)
+ faceValue(dVfcorrectionValues, facei);
@ -212,7 +215,7 @@ void Foam::isoAdvection::limitFluxes
template<class SpType, class SuType>
void Foam::isoAdvection::boundFlux
(
const scalarField& alpha1,
const bitSet& nextToInterface,
surfaceScalarField& dVfcorrectionValues,
DynamicList<label>& correctedFaces,
const SpType& Sp,
@ -223,7 +226,7 @@ void Foam::isoAdvection::boundFlux
scalar rDeltaT = 1/mesh_.time().deltaTValue();
correctedFaces.clear();
scalar aTol = 10*SMALL; // Note: tolerances
const scalar aTol = 100*SMALL; // Note: tolerances
const scalarField& meshV = mesh_.cellVolumes();
const scalar dt = mesh_.time().deltaTValue();
@ -235,13 +238,15 @@ void Foam::isoAdvection::boundFlux
const volScalarField& alphaOld = alpha1_.oldTime();
// Loop through alpha cell centred field
forAll(alpha1, celli)
for (label celli: nextToInterface)
{
if (alpha1[celli] < -aTol || alpha1[celli] > 1 + aTol)
if (alpha1_[celli] < -aTol || alpha1_[celli] > 1 + aTol)
{
const scalar Vi = meshV[celli];
scalar alphaOvershoot = pos0(alpha1[celli]-1)*(alpha1[celli]-1)
+ neg0(alpha1[celli])*alpha1[celli];
scalar alphaOvershoot =
pos0(alpha1_[celli] - 1)*(alpha1_[celli] - 1)
+ neg0(alpha1_[celli])*alpha1_[celli];
scalar fluidToPassOn = alphaOvershoot*Vi;
label nFacesToPassFluidThrough = 1;
@ -251,7 +256,7 @@ void Foam::isoAdvection::boundFlux
// not filled and to which dVf < phi*dt
for (label iter=0; iter<10; iter++)
{
if(mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
if (mag(alphaOvershoot) < aTol || nFacesToPassFluidThrough == 0)
{
break;
}
@ -271,9 +276,8 @@ void Foam::isoAdvection::boundFlux
scalar dVftot = 0;
nFacesToPassFluidThrough = 0;
forAll(downwindFaces, fi)
for (const label facei : downwindFaces)
{
const label facei = downwindFaces[fi];
const scalar phif = faceValue(phi_, facei);
const scalar dVff =
@ -343,7 +347,7 @@ void Foam::isoAdvection::boundFlux
scalar alpha1New =
(
alphaOld[celli]*rDeltaT + Su[celli]
alphaOld[celli]*rDeltaT + Su[celli]
- netFlux(dVf_, celli)/Vi*rDeltaT
- netFlux(dVfcorrectionValues, celli)/Vi*rDeltaT
)
@ -351,7 +355,7 @@ void Foam::isoAdvection::boundFlux
(rDeltaT - Sp[celli]);
alphaOvershoot =
pos0(alpha1New-1)*(alpha1New-1)
pos0(alpha1New - 1)*(alpha1New - 1)
+ neg0(alpha1New)*alpha1New;
fluidToPassOn = alphaOvershoot*Vi;
@ -372,13 +376,25 @@ void Foam::isoAdvection::advect(const SpType& Sp, const SuType& Su)
{
DebugInFunction << endl;
if (mesh_.topoChanging())
{
setProcessorPatches();
}
scalar advectionStartTime = mesh_.time().elapsedCpuTime();
scalar rDeltaT = 1/mesh_.time().deltaTValue();
const scalar rDeltaT = 1/mesh_.time().deltaTValue();
// reconstruct the interface
surf_->reconstruct();
if (timeIndex_ < mesh_.time().timeIndex())
{
timeIndex_= mesh_.time().timeIndex();
surf_->normal().oldTime() = surf_->normal();
surf_->centre().oldTime() = surf_->centre();
}
// Initialising dVf with upwind values
// i.e. phi[facei]*alpha1[upwindCell[facei]]*dt
dVf_ = upwind<scalar>(mesh_, phi_).flux(alpha1_)*mesh_.time().deltaT();

32
src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.C
View File

@ -37,8 +37,17 @@ int Foam::cutCell::debug = 0;
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cutCell::cutCell(const fvMesh&)
{}
Foam::cutCell::cutCell
(
const fvMesh& mesh
)
{
// required as otherwise setAlphaFields might not work in parallel
mesh.C();
mesh.V();
mesh.Cf();
mesh.magSf();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -47,7 +56,8 @@ void Foam::cutCell::calcCellData
(
const DynamicList<point>& cutFaceCentres,
const DynamicList<vector>& cutFaceAreas,
vector& subCellCentre, scalar& subCellVolume
vector& subCellCentre,
scalar& subCellVolume
)
{
// Clear the fields for accumulation
@ -159,11 +169,27 @@ void Foam::cutCell::calcIsoFacePointsFromEdges
DynamicList<point>& facePoints
)
{
if (mag(faceArea) < VSMALL)
{
facePoints.clear();
return;
}
const vector zhat = normalised(faceArea);
vector xhat = faceEdges[0][0] - faceCentre;
xhat = (xhat - (xhat & zhat)*zhat);
xhat.normalise();
if (mag(xhat) == 0)
{
facePoints.clear();
return;
}
vector yhat = normalised(zhat ^ xhat);
if (mag(yhat) == 0)
{
facePoints.clear();
return;
}
yhat.normalise();
// Calculating all intersection points
DynamicList<point> unsortedFacePoints(3 * faceEdges.size());

2
src/transportModels/geometricVoF/cellCuts/cutCell/cutCell.H
View File

@ -99,7 +99,7 @@ public:
// Constructors
//- Construct from fvMesh
explicit cutCell(const fvMesh& unused);
explicit cutCell(const fvMesh& mesh);
};

46
src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.C
View File

@ -115,7 +115,7 @@ Foam::label Foam::cutCellIso::calcSubCell
// In the rare but occuring cases where a cell is only touched at a
// point or a line the isoFaceArea_ will have zero length and here the
// cell should be treated as either completely empty or full.
if (mag(faceArea_) < 10*SMALL)
if (mag(faceArea_) < ROOTVSMALL)
{
if (nFaceBelowInterface == 0)
{
@ -170,21 +170,9 @@ Foam::label Foam::cutCellIso::calcSubCell
}
const Foam::point& Foam::cutCellIso::subCellCentre() const
{
return subCellCentre_;
}
Foam::scalar Foam::cutCellIso::subCellVolume() const
{
return subCellVolume_;
}
const Foam::DynamicList<Foam::point>& Foam::cutCellIso::facePoints()
{
if (facePoints_.size() == 0)
if (facePoints_.empty())
{
// get face points in sorted order
calcIsoFacePointsFromEdges
@ -200,36 +188,6 @@ const Foam::DynamicList<Foam::point>& Foam::cutCellIso::facePoints()
}
const Foam::point& Foam::cutCellIso::faceCentre() const
{
return faceCentre_;
}
const Foam::vector& Foam::cutCellIso::faceArea() const
{
return faceArea_;
}
Foam::label Foam::cutCellIso::cellStatus() const
{
return cellStatus_;
}
Foam::scalar Foam::cutCellIso::VolumeOfFluid() const
{
return VOF_;
}
Foam::scalar Foam::cutCellIso::cutValue() const
{
return cutValue_;
}
void Foam::cutCellIso::clearStorage()
{
cellI_ = -1;

41
src/transportModels/geometricVoF/cellCuts/cutCell/cutCellIso.H
View File

@ -148,28 +148,49 @@ class cutCellIso
label calcSubCell(const label cellI, const scalar cutValue);
//- Returns subCellCentre
const point& subCellCentre() const;
const point& subCellCentre() const noexcept
{
return subCellCentre_;
}
//- Returns subCellVolume
scalar subCellVolume() const;
scalar subCellVolume() const noexcept
{
return subCellVolume_;
}
//- Returns the points of the cutting isoface
//- Returns the points of the cutting PLICface
const DynamicList<point>& facePoints();
//- Returns the centre of the cutting isoface
const point& faceCentre() const;
//- Returns the centre of the cutting PLICface
const point& faceCentre() const noexcept
{
return faceCentre_;
}
//- Returns the area normal vector of the cutting isoface
const vector& faceArea() const;
//- Returns the area normal vector of the cutting PLICface
const vector& faceArea() const noexcept
{
return faceArea_;
}
//- Returns cellStatus
label cellStatus() const;
label cellStatus() const noexcept
{
return cellStatus_;
}
//- Returns volume of fluid value
scalar VolumeOfFluid() const;
scalar VolumeOfFluid() const noexcept
{
return VOF_;
}
//- Returns cutValue
scalar cutValue() const;
scalar cutValue() const noexcept
{
return cutValue_;
}
//- Resets internal values
void clearStorage();

46
src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.C
View File

@ -117,7 +117,7 @@ Foam::label Foam::cutCellPLIC::calcSubCell
// In the rare but occuring cases where a cell is only touched at a
// point or a line the isoFaceArea_ will have zero length and here the
// cell should be treated as either completely empty or full.
if (mag(faceArea_) < 10*SMALL)
if (mag(faceArea_) < ROOTVSMALL)
{
if (nFaceBelowInterface == 0)
{
@ -172,21 +172,9 @@ Foam::label Foam::cutCellPLIC::calcSubCell
}
const Foam::point& Foam::cutCellPLIC::subCellCentre() const
{
return subCellCentre_;
}
Foam::scalar Foam::cutCellPLIC::subCellVolume() const
{
return subCellVolume_;
}
const Foam::DynamicList<Foam::point>& Foam::cutCellPLIC::facePoints()
{
if (facePoints_.size() == 0)
if (facePoints_.empty())
{
// get face points in sorted order
calcIsoFacePointsFromEdges
@ -202,36 +190,6 @@ const Foam::DynamicList<Foam::point>& Foam::cutCellPLIC::facePoints()
}
const Foam::point& Foam::cutCellPLIC::faceCentre() const
{
return faceCentre_;
}
const Foam::vector& Foam::cutCellPLIC::faceArea() const
{
return faceArea_;
}
Foam::scalar Foam::cutCellPLIC::VolumeOfFluid() const
{
return VOF_;
}
Foam::label Foam::cutCellPLIC::cellStatus() const
{
return cellStatus_;
}
Foam::scalar Foam::cutCellPLIC::cutValue() const
{
return cutValue_;
}
void Foam::cutCellPLIC::clearStorage()
{
cellI_ = -1;

35
src/transportModels/geometricVoF/cellCuts/cutCell/cutCellPLIC.H
View File

@ -147,28 +147,49 @@ class cutCellPLIC
);
//- Returns subCellCentre
const point& subCellCentre() const;
const point& subCellCentre() const noexcept
{
return subCellCentre_;
}
//- Returns subCellVolume
scalar subCellVolume() const;
scalar subCellVolume() const noexcept
{
return subCellVolume_;
}
//- Returns the points of the cutting PLICface
const DynamicList<point>& facePoints();
//- Returns the centre of the cutting PLICface
const point& faceCentre() const;
const point& faceCentre() const noexcept
{
return faceCentre_;
}
//- Returns the area normal vector of the cutting PLICface
const vector& faceArea() const;
const vector& faceArea() const noexcept
{
return faceArea_;
}
//- Returns cellStatus
label cellStatus() const;
label cellStatus() const noexcept
{
return cellStatus_;
}
//- Returns volume of fluid value
scalar VolumeOfFluid() const;
scalar VolumeOfFluid() const noexcept
{
return VOF_;
}
//- Returns cutValue
scalar cutValue() const;
scalar cutValue() const noexcept
{
return cutValue_;
}
//- Resets internal values
void clearStorage();

17
src/transportModels/geometricVoF/cellCuts/cutFace/cutFace.C
View File

@ -83,8 +83,8 @@ void Foam::cutFace::calcSubFace
if
(
(pointStatus[idx] < 0 && pointStatus[nextIdx] > 0) ||
(pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
(pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
|| (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
) // cut on edge cut Value is zero
{
label nextP = f.nextLabel(idx);
@ -161,8 +161,8 @@ void Foam::cutFace::calcSubFace
if
(
(pointStatus[idx] < 0 && pointStatus[nextIdx] > 0) ||
(pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
(pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
|| (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
) // cut on edge cut Value is zero
{
label nextP = f.nextLabel(idx);
@ -232,8 +232,8 @@ void Foam::cutFace::calcSubFace
if
(
(pointStatus[idx] < 0 && pointStatus[nextIdx] > 0) ||
(pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
(pointStatus[idx] < 0 && pointStatus[nextIdx] > 0)
|| (pointStatus[idx] > 0 && pointStatus[nextIdx] < 0)
)
{
label nextP = f.nextLabel(idx);
@ -328,7 +328,10 @@ void Foam::cutFace::calcSubFaceCentreAndArea
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cutFace::cutFace(const fvMesh& mesh)
Foam::cutFace::cutFace
(
const fvMesh& mesh
)
:
mesh_(mesh)
{}

28
src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.C
View File

@ -246,7 +246,7 @@ Foam::scalar Foam::cutFaceAdvect::timeIntegratedFaceFlux
// Check if pTimes changes direction more than twice when looping face
label nShifts = 0;
forAll(pTimes_, pi) // i have no clue what this is
forAll(pTimes_, pi)
{
const label oldEdgeSign =
sign(pTimes_[(pi + 1) % nPoints] - pTimes_[pi]);
@ -356,7 +356,7 @@ Foam::scalar Foam::cutFaceAdvect::timeIntegratedFaceFlux
// Check if pTimes changes direction more than twice when looping face
label nShifts = 0;
forAll(pTimes_, pi) // i have no clue what this is
forAll(pTimes_, pi)
{
const label oldEdgeSign =
sign(pTimes_[(pi + 1) % nPoints] - pTimes_[pi]);
@ -984,30 +984,6 @@ void Foam::cutFaceAdvect::cutPoints
}
const Foam::point& Foam::cutFaceAdvect::subFaceCentre() const
{
return subFaceCentre_;
}
const Foam::vector& Foam::cutFaceAdvect::subFaceArea() const
{
return subFaceArea_;
}
const Foam::DynamicList<Foam::point>& Foam::cutFaceAdvect::subFacePoints() const
{
return subFacePoints_;
}
const Foam::DynamicList<Foam::point>& Foam::cutFaceAdvect::surfacePoints() const
{
return surfacePoints_;
}
void Foam::cutFaceAdvect::clearStorage()
{
subFaceCentre_ = Zero;

23
src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceAdvect.H
View File

@ -232,16 +232,29 @@ class cutFaceAdvect
);
//- Returns centre of cutted face
const point& subFaceCentre() const;
const point& subFaceCentre() const noexcept
{
return subFaceCentre_;
}
//- Returns area vector of cutted face
const vector& subFaceArea() const;
const vector& subFaceArea() const noexcept
{
return subFaceArea_;
}
//- Returns the cut edge of the cutted face
const DynamicList<point>& subFacePoints() const;
const DynamicList<point>& subFacePoints() const noexcept
{
return subFacePoints_;
}
//- Returns point of the face in sorted of cutted face
const DynamicList<point>& surfacePoints() const noexcept
{
return surfacePoints_;
}
//- Returns point of the cutted face in sorted order
const DynamicList<point>& surfacePoints() const;
//- Resets internal variables
void clearStorage();

24
src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.C
View File

@ -111,30 +111,6 @@ Foam::label Foam::cutFaceIso::calcSubFace
}
const Foam::point& Foam::cutFaceIso::subFaceCentre() const
{
return subFaceCentre_;
}
const Foam::vector& Foam::cutFaceIso::subFaceArea() const
{
return subFaceArea_;
}
const Foam::DynamicList<Foam::point>& Foam::cutFaceIso::subFacePoints() const
{
return subFacePoints_;
}
const Foam::DynamicList<Foam::point>& Foam::cutFaceIso::surfacePoints() const
{
return surfacePoints_;
}
void Foam::cutFaceIso::clearStorage()
{
subFaceCentre_ = Zero;

20
src/transportModels/geometricVoF/cellCuts/cutFace/cutFaceIso.H
View File

@ -124,16 +124,28 @@ public:
);
//- Returns centre of cutted face
const point& subFaceCentre() const;
const point& subFaceCentre() const noexcept
{
return subFaceCentre_;
}
//- Returns area vector of cutted face
const vector& subFaceArea() const;
const vector& subFaceArea() const noexcept
{
return subFaceArea_;
}
//- Returns the cut edge of the cutted face
const DynamicList<point>& subFacePoints() const;
const DynamicList<point>& subFacePoints() const noexcept
{
return subFacePoints_;
}
//- Returns point of the face in sorted of cutted face
const DynamicList<point>& surfacePoints() const;
const DynamicList<point>& surfacePoints() const noexcept
{
return surfacePoints_;
}
//- Resets internal variables
void clearStorage();

24
src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.C
View File

@ -115,30 +115,6 @@ Foam::label Foam::cutFacePLIC::calcSubFace
}
const Foam::point& Foam::cutFacePLIC::subFaceCentre() const
{
return subFaceCentre_;
}
const Foam::vector& Foam::cutFacePLIC::subFaceArea() const
{
return subFaceArea_;
}
const Foam::DynamicList<Foam::point>& Foam::cutFacePLIC::subFacePoints() const
{
return subFacePoints_;
}
const Foam::DynamicList<Foam::point>& Foam::cutFacePLIC::surfacePoints() const
{
return surfacePoints_;
}
void Foam::cutFacePLIC::clearStorage()
{
subFaceCentre_ = Zero;

20
src/transportModels/geometricVoF/cellCuts/cutFace/cutFacePLIC.H
View File

@ -123,16 +123,28 @@ public:
);
//- Returns centre of cutted face
const point& subFaceCentre() const;
const point& subFaceCentre() const noexcept
{
return subFaceCentre_;
}
//- Returns area vector of cutted face
const vector& subFaceArea() const;
const vector& subFaceArea() const noexcept
{
return subFaceArea_;
}
//- Returns the cut edge of the cutted face
const DynamicList<point>& subFacePoints() const;
const DynamicList<point>& subFacePoints() const noexcept
{
return subFacePoints_;
}
//- Returns point of the face in sorted of cutted face
const DynamicList<point>& surfacePoints() const;
const DynamicList<point>& surfacePoints() const noexcept
{
return surfacePoints_;
}
//- Resets internal variables
void clearStorage();

5
src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
View File

@ -104,7 +104,7 @@ void Foam::reconstructedDistanceFunction::markCellsNearSurf
// do coupled face first
Map<bool> syncMap;
for (int level=0;level<=neiRingLevel;level++)
for (label level=0;level<=neiRingLevel;level++)
{
// parallel
if (level > 0)
@ -280,9 +280,8 @@ const Foam::volScalarField& Foam::reconstructedDistanceFunction::constructRDF
scalar avgWeight = 0;
const point p = mesh_.C()[celli];
forAll(stencil[celli], i)
for (const label gblIdx : stencil[celli])
{
const label gblIdx = stencil[celli][i];
vector n = -distribute.getValue(normal, mapNormal, gblIdx);
if (mag(n) != 0)
{

4
src/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.H
View File

@ -108,12 +108,12 @@ public:
surfaceVectorField::Boundary& nHatb
);
const volScalarField& cellDistLevel() const
const volScalarField& cellDistLevel() const noexcept
{
return cellDistLevel_;
}
const boolList& nextToInterface() const
const boolList& nextToInterface() const noexcept
{
return nextToInterface_;
}

15
src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.C
View File

@ -48,21 +48,23 @@ void Foam::reconstruction::gradAlpha::gradSurf(const volScalarField& phi)
{
leastSquareGrad<scalar> lsGrad("polyDegree1",mesh_.geometricD());
exchangeFields_.setUpCommforZone(interfaceCell_,true);
zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
exchangeFields.setUpCommforZone(interfaceCell_,true);
Map<vector> mapCC
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, mesh_.C())
exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
);
Map<scalar> mapPhi
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, phi)
exchangeFields.getDatafromOtherProc(interfaceCell_, phi)
);
DynamicField<vector> cellCentre(100);
DynamicField<scalar> phiValues(100);
const labelListList& stencil = exchangeFields_.getStencil();
const labelListList& stencil = exchangeFields.getStencil();
forAll(interfaceLabels_, i)
{
@ -75,11 +77,11 @@ void Foam::reconstruction::gradAlpha::gradSurf(const volScalarField& phi)
{
cellCentre.append
(
exchangeFields_.getValue(mesh_.C(), mapCC, gblIdx)
exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
);
phiValues.append
(
exchangeFields_.getValue(phi, mapPhi, gblIdx)
exchangeFields.getValue(phi, mapPhi, gblIdx)
);
}
@ -111,7 +113,6 @@ Foam::reconstruction::gradAlpha::gradAlpha
interfaceNormal_(fvc::grad(alpha1)),
isoFaceTol_(modelDict().getOrDefault<scalar>("isoFaceTol", 1e-8)),
surfCellTol_(modelDict().getOrDefault<scalar>("surfCellTol", 1e-8)),
exchangeFields_(zoneDistribute::New(mesh_)),
sIterPLIC_(mesh_,surfCellTol_)
{
reconstruct();

4
src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/gradAlpha/gradAlpha.H
View File

@ -83,13 +83,9 @@ class gradAlpha
// Those with surfCellTol_ < alpha1 < 1 - surfCellTol_
scalar surfCellTol_;
//- Provides stencil and map
zoneDistribute& exchangeFields_;
//- SurfaceIterator finds the plane centre for specified VOF value
surfaceIteratorPLIC sIterPLIC_;
// Private Member Functions
//- Compute gradient at the surfaces

58
src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C
View File

@ -48,27 +48,28 @@ namespace reconstruction
void Foam::reconstruction::plicRDF::interpolateNormal()
{
scalar dt = mesh_.time().deltaTValue();
zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
leastSquareGrad<scalar> lsGrad("polyDegree1",mesh_.geometricD());
exchangeFields_.setUpCommforZone(interfaceCell_,false);
exchangeFields.setUpCommforZone(interfaceCell_,false);
Map<vector> mapCentre
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, centre_)
exchangeFields.getDatafromOtherProc(interfaceCell_, centre_)
);
Map<vector> mapNormal
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, normal_)
exchangeFields.getDatafromOtherProc(interfaceCell_, normal_)
);
Map<vector> mapCC
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, mesh_.C())
exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
);
Map<scalar> mapAlpha
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, alpha1_)
exchangeFields.getDatafromOtherProc(interfaceCell_, alpha1_)
);
DynamicField<vector > cellCentre(100);
@ -76,7 +77,7 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
DynamicList<vector> foundNormals(30);
const labelListList& stencil = exchangeFields_.getStencil();
const labelListList& stencil = exchangeFields.getStencil();
forAll(interfaceLabels_, i)
{
@ -84,17 +85,16 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
vector estimatedNormal{Zero};
scalar weight{0};
foundNormals.clear();
forAll(stencil[celli], i)
for (const label gblIdx : stencil[celli])
{
const label gblIdx = stencil[celli][i];
vector n =
exchangeFields_.getValue(normal_, mapNormal, gblIdx);
exchangeFields.getValue(normal_, mapNormal, gblIdx);
point p = mesh_.C()[celli]-U_[celli]*dt;
if (mag(n) != 0)
{
n /= mag(n);
vector centre =
exchangeFields_.getValue(centre_, mapCentre, gblIdx);
exchangeFields.getValue(centre_, mapCentre, gblIdx);
vector distanceToIntSeg = (tensor::I- n*n) & (p - centre);
estimatedNormal += n /max(mag(distanceToIntSeg), SMALL);
weight += 1/max(mag(distanceToIntSeg), SMALL);
@ -143,11 +143,11 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
const label gblIdx = stencil[celli][i];
cellCentre.append
(
exchangeFields_.getValue(mesh_.C(), mapCC, gblIdx)
exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
);
alphaValues.append
(
exchangeFields_.getValue(alpha1_, mapAlpha, gblIdx)
exchangeFields.getValue(alpha1_, mapAlpha, gblIdx)
);
}
cellCentre -= mesh_.C()[celli];
@ -160,22 +160,23 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
{
leastSquareGrad<scalar> lsGrad("polyDegree1", mesh_.geometricD());
zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
exchangeFields_.setUpCommforZone(interfaceCell_, false);
exchangeFields.setUpCommforZone(interfaceCell_, false);
Map<vector> mapCC
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, mesh_.C())
exchangeFields.getDatafromOtherProc(interfaceCell_, mesh_.C())
);
Map<scalar> mapPhi
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, phi)
exchangeFields.getDatafromOtherProc(interfaceCell_, phi)
);
DynamicField<vector> cellCentre(100);
DynamicField<scalar> phiValues(100);
const labelListList& stencil = exchangeFields_.getStencil();
const labelListList& stencil = exchangeFields.getStencil();
forAll(interfaceLabels_, i)
{
@ -188,11 +189,11 @@ void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
{
cellCentre.append
(
exchangeFields_.getValue(mesh_.C(), mapCC, gblIdx)
exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
);
phiValues.append
(
exchangeFields_.getValue(phi, mapPhi, gblIdx)
exchangeFields.getValue(phi, mapPhi, gblIdx)
);
}
@ -204,6 +205,8 @@ void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
void Foam::reconstruction::plicRDF::setInitNormals(bool interpolate)
{
zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
interfaceLabels_.clear();
forAll(alpha1_, celli)
@ -218,7 +221,7 @@ void Foam::reconstruction::plicRDF::setInitNormals(bool interpolate)
RDF_.markCellsNearSurf(interfaceCell_, 1);
const boolList& nextToInterface_ = RDF_.nextToInterface();
exchangeFields_.updateStencil(nextToInterface_);
exchangeFields.updateStencil(nextToInterface_);
if (interpolate)
{
@ -236,15 +239,15 @@ void Foam::reconstruction::plicRDF::calcResidual
List<normalRes>& normalResidual
)
{
exchangeFields_.setUpCommforZone(interfaceCell_,false);
zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
exchangeFields.setUpCommforZone(interfaceCell_,false);
Map<vector> mapNormal
(
exchangeFields_.getDatafromOtherProc(interfaceCell_, normal_)
exchangeFields.getDatafromOtherProc(interfaceCell_, normal_)
);
const labelListList& stencil = exchangeFields_.getStencil();
const labelListList& stencil = exchangeFields.getStencil();
forAll(interfaceLabels_, i)
{
@ -265,7 +268,7 @@ void Foam::reconstruction::plicRDF::calcResidual
{
const label gblIdx = stencil[celli][j];
vector normal =
exchangeFields_.getValue(normal_, mapNormal, gblIdx);
exchangeFields.getValue(normal_, mapNormal, gblIdx);
if (mag(normal) != 0 && j != 0)
{
@ -328,7 +331,6 @@ Foam::reconstruction::plicRDF::plicRDF
iteration_(modelDict().getOrDefault<label>("iterations", 5)),
interpolateNormal_(modelDict().getOrDefault("interpolateNormal", true)),
RDF_(mesh_),
exchangeFields_(zoneDistribute::New(mesh_)),
sIterPLIC_(mesh_,surfCellTol_)
{
setInitNormals(false);
@ -375,6 +377,7 @@ Foam::reconstruction::plicRDF::plicRDF
void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
{
zoneDistribute& exchangeFields = zoneDistribute::New(mesh_);
const bool uptodate = alreadyReconstructed(forceUpdate);
if (uptodate && !forceUpdate)
@ -403,7 +406,7 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
bitSet tooCoarse(mesh_.nCells(),false);
for (int iter=0; iter<iteration_; ++iter)
for (label iter=0; iter<iteration_; ++iter)
{
forAll(interfaceLabels_, i)
{
@ -452,7 +455,7 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
nextToInterface_,
centre_,
normal_,
exchangeFields_,
exchangeFields,
false
);
RDF_.updateContactAngle(alpha1_, U_, nHatb);
@ -492,7 +495,6 @@ void Foam::reconstruction::plicRDF::reconstruct(bool forceUpdate)
resCounter++;
}
}
reduce(avgRes,sumOp<scalar>());

3
src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.H
View File

@ -118,9 +118,6 @@ class plicRDF
//- Calculates the RDF function
reconstructedDistanceFunction RDF_;
//- Provides stencil and map
zoneDistribute& exchangeFields_;
//- surfaceIterator finds the plane centre for specified VOF value
surfaceIteratorPLIC sIterPLIC_;

12
src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.C
View File

@ -106,11 +106,13 @@ Foam::reconstructionSchemes::reconstructionSchemes
(
IOobject
(
"recon::normal",
IOobject::groupName("interfaceNormal", alpha1.group()),
alpha1_.mesh().time().timeName(),
alpha1_.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
dict.getOrDefault("writeFields",false)
? IOobject::AUTO_WRITE
: IOobject::NO_WRITE
),
alpha1_.mesh(),
dimensionedVector(dimArea, Zero)
@ -119,11 +121,13 @@ Foam::reconstructionSchemes::reconstructionSchemes
(
IOobject
(
"recon::centre",
IOobject::groupName("interfaceCentre", alpha1.group()),
alpha1_.mesh().time().timeName(),
alpha1_.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
dict.getOrDefault("writeFields",false)
? IOobject::AUTO_WRITE
: IOobject::NO_WRITE
),
alpha1_.mesh(),
dimensionedVector(dimLength, Zero)

24
src/transportModels/geometricVoF/reconstructionSchemes/reconstructionSchemes.H
View File

@ -209,13 +209,13 @@ public:
// Member Functions
//- Access to the model dictionary
dictionary& modelDict()
dictionary& modelDict() noexcept
{
return reconstructionSchemesCoeffs_;
}
//- Const access to the model dictionary
const dictionary& modelDict() const
const dictionary& modelDict() const noexcept
{
return reconstructionSchemesCoeffs_;
}
@ -223,26 +223,38 @@ public:
//- Reconstruct the interface
virtual void reconstruct(bool forceUpdate = true) = 0;
//- const-Reference to interface normals
const volVectorField& normal() const noexcept
{
return normal_;
}
//- const-Reference to interface centres
const volVectorField& centre() const noexcept
{
return centre_;
}
//- Reference to interface normals
const volVectorField& normal() const
volVectorField& normal() noexcept
{
return normal_;
}
//- Reference to interface centres
const volVectorField& centre() const
volVectorField& centre() noexcept
{
return centre_;
}
//- Does the cell contain interface
const boolList& interfaceCell() const
const boolList& interfaceCell() const noexcept
{
return interfaceCell_;
}
//- List of cells with an interface
const DynamicField<label>& interfaceLabels() const
const DynamicField<label>& interfaceLabels() const noexcept
{
return interfaceLabels_;
}

1
tutorials/multiphase/interIsoFoam/damBreak/system/fvSolution
View File

@ -25,6 +25,7 @@ solvers
snapTol 1e-12;
clip true;
reconstructionScheme isoAlpha;
writeFields true;
nAlphaSubCycles 1;
cAlpha 1; // Note: cAlpha is not used by isoAdvector but must