98 nAlphaBounds_(dict_.getOrDefault<label>(
"nAlphaBounds", 3)),
99 isoFaceTol_(dict_.getOrDefault<scalar>(
"isoFaceTol", 1
e-8)),
100 surfCellTol_(dict_.getOrDefault<scalar>(
"surfCellTol", 1
e-8)),
101 writeIsoFacesToFile_(dict_.getOrDefault(
"writeIsoFaces", false)),
104 surfCells_(label(0.2*mesh_.nCells())),
107 bsFaces_(label(0.2*mesh_.nBoundaryFaces())),
108 bsx0_(bsFaces_.size()),
109 bsn0_(bsFaces_.size()),
110 bsUn0_(bsFaces_.size()),
113 porosityEnabled_(dict_.getOrDefault<
bool>(
"porosityEnabled", false)),
114 porosityPtr_(nullptr),
117 procPatchLabels_(mesh_.
boundary().size()),
118 surfaceCellFacesOnProcPatches_(0)
120 cutFaceAdvect::debug = debug;
138 setProcessorPatches();
146 "porosityProperties",
157 if (porosityEnabled_)
170 <<
"Porosity field has values <= 0 or > 1"
177 <<
"Porosity enabled in constant/porosityProperties "
178 <<
"but no porosity field is found in object registry."
187void Foam::isoAdvection::setProcessorPatches()
190 surfaceCellFacesOnProcPatches_.
clear();
191 surfaceCellFacesOnProcPatches_.resize(
patches.
size());
194 procPatchLabels_.clear();
199 isA<processorPolyPatch>(
patches[patchi])
203 procPatchLabels_.append(patchi);
209void Foam::isoAdvection::extendMarkedCells
215 bitSet markedFace(mesh_.nFaces());
217 for (
const label celli : markedCell)
219 markedFace.set(mesh_.cells()[celli]);
225 for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
227 if (markedFace.test(facei))
229 markedCell.set(mesh_.faceOwner()[facei]);
230 markedCell.set(mesh_.faceNeighbour()[facei]);
233 for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); ++facei)
235 if (markedFace.test(facei))
237 markedCell.set(mesh_.faceOwner()[facei]);
243void Foam::isoAdvection::timeIntegratedFlux()
247 const scalar dt = mesh_.time().deltaTValue();
250 interpolationCellPoint<vector> UInterp(U_);
253 label nSurfaceCells = 0;
261 const scalarField& magSfIn = mesh_.magSf().primitiveField();
265 const cellList& cellFaces = mesh_.cells();
266 const labelList& own = mesh_.faceOwner();
270 DynamicList<List<point>> isoFacePts;
271 const DynamicField<label>& interfaceLabels = surf_->interfaceLabels();
277 const label celli = interfaceLabels[i];
278 const point x0(surf_->centre()[celli]);
280 Un0[i] = UInterp.interpolate(x0, celli) & n0;
284 if (porosityEnabled_)
288 const label celli = interfaceLabels[i];
289 Un0[i] /= porosityPtr_->primitiveField()[celli];
294 forAll(interfaceLabels, i)
296 const label celli = interfaceLabels[i];
297 if (
mag(surf_->normal()[celli]) != 0)
302 surfCells_.append(celli);
303 const point x0(surf_->centre()[celli]);
307 <<
"\n------------ Cell " << celli <<
" with alpha1 = "
308 << alpha1In_[celli] <<
" and 1-alpha1 = "
309 << 1.0 - alpha1In_[celli] <<
" ------------"
315 const cell& celliFaces = cellFaces[celli];
316 for (
const label facei : celliFaces)
318 if (mesh_.isInternalFace(facei))
320 bool isDownwindFace =
false;
322 if (celli == own[facei])
324 if (phiIn[facei] >= 0)
326 isDownwindFace =
true;
331 if (phiIn[facei] < 0)
333 isDownwindFace =
true;
339 dVfIn[facei] = advectFace_.timeIntegratedFaceFlux
354 bsFaces_.append(facei);
357 bsUn0_.append(Un0[i]);
367 const polyBoundaryMesh& boundaryMesh = mesh_.boundaryMesh();
371 const label nInternalFaces = mesh_.nInternalFaces();
377 const label facei = bsFaces_[i];
378 const label patchi = boundaryMesh.patchID()[facei - nInternalFaces];
379 const label start = boundaryMesh[patchi].start();
381 if (!
phib[patchi].empty())
383 const label patchFacei = facei - start;
384 const scalar phiP =
phib[patchi][patchFacei];
388 const scalar magSf = magSfb[patchi][patchFacei];
390 dVfb[patchi][patchFacei] = advectFace_.timeIntegratedFaceFlux
403 checkIfOnProcPatch(facei);
409 syncProcPatches(dVf_, phi_);
411 writeIsoFaces(isoFacePts);
413 DebugInfo <<
"Number of isoAdvector surface cells = "
418void Foam::isoAdvection::setDownwindFaces
421 DynamicLabelList& downwindFaces
427 const labelList& own = mesh_.faceOwner();
429 const cell&
c =
cells[celli];
431 downwindFaces.
clear();
434 for (
const label facei: c)
437 const scalar
phi = faceValue(phi_, facei);
439 if (own[facei] == celli)
443 downwindFaces.append(facei);
448 downwindFaces.append(facei);
452 downwindFaces.shrink();
456Foam::scalar Foam::isoAdvection::netFlux
465 const cell&
c = mesh_.cells()[celli];
468 const labelList& own = mesh_.faceOwner();
470 for (
const label facei : c)
472 const scalar dVff = faceValue(dVf, facei);
474 if (own[facei] == celli)
492 bool returnSyncedFaces
495 DynamicLabelList syncedFaces(0);
496 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
500 DynamicList<label> neighProcs;
504 for (
const label patchi : procPatchLabels_)
506 const processorPolyPatch& procPatch =
507 refCast<const processorPolyPatch>(
patches[patchi]);
508 const label nbrProci = procPatch.neighbProcNo();
510 neighProcs.append(nbrProci);
511 UOPstream toNbr(nbrProci, pBufs);
513 const scalarField& pFlux = dVf.boundaryField()[patchi];
514 const List<label>& surfCellFacesOnProcPatch =
515 surfaceCellFacesOnProcPatches_[patchi];
517 const UIndirectList<scalar> dVfPatch
520 surfCellFacesOnProcPatch
523 toNbr << surfCellFacesOnProcPatch << dVfPatch;
527 pBufs.finishedNeighbourSends(neighProcs);
531 for (
const label patchi : procPatchLabels_)
533 const processorPolyPatch& procPatch =
534 refCast<const processorPolyPatch>(
patches[patchi]);
535 const label nbrProci = procPatch.neighbProcNo();
537 UIPstream fromNeighb(nbrProci, pBufs);
539 List<scalar> nbrdVfs;
541 fromNeighb >> faceIDs >> nbrdVfs;
542 if (returnSyncedFaces)
544 List<label> syncedFaceI(faceIDs);
545 for (label& faceI : syncedFaceI)
547 faceI += procPatch.start();
549 syncedFaces.append(syncedFaceI);
554 Pout<<
"Received at time = " << mesh_.time().value()
555 <<
": surfCellFacesOnProcPatch = " << faceIDs <<
nl
556 <<
"Received at time = " << mesh_.time().value()
557 <<
": dVfPatch = " << nbrdVfs <<
endl;
561 scalarField& localFlux = dVf.boundaryFieldRef()[patchi];
565 const label facei = faceIDs[i];
566 localFlux[facei] = - nbrdVfs[i];
567 if (debug &&
mag(localFlux[facei] + nbrdVfs[i]) > ROOTVSMALL)
569 Pout<<
"localFlux[facei] = " << localFlux[facei]
570 <<
" and nbrdVfs[i] = " << nbrdVfs[i]
571 <<
" for facei = " << facei <<
endl;
579 forAll(procPatchLabels_, patchLabeli)
581 const label patchi = procPatchLabels_[patchLabeli];
582 const scalarField& localFlux = dVf.boundaryField()[patchi];
583 Pout<<
"time = " << mesh_.time().value() <<
": localFlux = "
584 << localFlux <<
endl;
589 forAll(surfaceCellFacesOnProcPatches_, patchi)
591 surfaceCellFacesOnProcPatches_[patchi].clear();
599void Foam::isoAdvection::checkIfOnProcPatch(
const label facei)
601 if (!mesh_.isInternalFace(facei))
603 const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
604 const label patchi = pbm.patchID()[facei - mesh_.nInternalFaces()];
606 if (isA<processorPolyPatch>(pbm[patchi]) && !pbm[patchi].empty())
608 const label patchFacei = pbm[patchi].whichFace(facei);
609 surfaceCellFacesOnProcPatches_[patchi].append(patchFacei);
615void Foam::isoAdvection::applyBruteForceBounding()
618 bool alpha1Changed =
false;
620 const scalar snapAlphaTol = dict_.getOrDefault<scalar>(
"snapTol", 0);
621 if (snapAlphaTol > 0)
625 *
pos0(alpha1_ - snapAlphaTol)
626 *
neg0(alpha1_ - (1.0 - snapAlphaTol))
627 +
pos0(alpha1_ - (1.0 - snapAlphaTol));
629 alpha1Changed =
true;
632 if (dict_.getOrDefault(
"clip",
true))
634 alpha1_ =
min(scalar(1),
max(scalar(0), alpha1_));
635 alpha1Changed =
true;
640 alpha1_.correctBoundaryConditions();
647 if (!mesh_.time().writeTime())
return;
649 if (dict_.getOrDefault(
"writeSurfCells",
false))
656 mesh_.time().timeName(),
674 if (!writeIsoFacesToFile_ || !mesh_.time().writeTime())
return;
679 mesh_.time().globalPath()
681 /
word::printf(
"isoFaces_%012d.obj", mesh_.time().timeIndex())
695 Info<<
nl <<
"isoAdvection: writing iso faces to file: "
699 forAll(allProcFaces, proci)
706 if (facePts.size() !=
f.
size())
720 Info<<
nl <<
"isoAdvection: writing iso faces to file: "
726 if (facePts.size() !=
f.
size())
Macros for easy insertion into run-time selection tables.
const volScalarField & alpha1
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
GeometricBoundaryField< scalar, fvsPatchField, surfaceMesh > Boundary
Type of boundary fields.
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Defines the attributes of an object for which implicit objectRegistry management is supported,...
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void clear()
Clear the list, i.e. set size to zero.
OFstream that keeps track of vertices.
virtual Ostream & write(const char c)
Write character.
virtual const fileName & name() const
Read/write access to the name of the stream.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void gatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
void clear()
Clear the PtrList. Delete allocated entries and set size to zero.
const word & constant() const
Return constant name.
void size(const label n)
Older name for setAddressableSize.
@ nonBlocking
"nonBlocking"
static bool & parRun() noexcept
Test if this a parallel run.
label size() const noexcept
The number of elements in the list.
A collection of cell labels.
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
A face is a list of labels corresponding to mesh vertices.
A class for handling file names.
static std::string path(const std::string &str)
Return directory path name (part before last /)
const Time & time() const
Return the top-level database.
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
An implementation of the isoAdvector geometric Volume-of-Fluid method advancing the provided volume f...
void writeIsoFaces(const DynamicList< List< point > > &isoFacePts) const
Write isoface points to .obj file.
void writeSurfaceCells() const
Return cellSet of surface cells.
bool foundObject(const word &name, const bool recursive=false) const
Is the named Type found?
Type * getObjectPtr(const word &name, const bool recursive=false) const
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
const labelList & patchID() const
Per boundary face label the patch index.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
const vectorField & faceCentres() const
const scalarField & cellVolumes() const
const vectorField & cellCentres() const
const labelListList & cellCells() const
const labelListList & cellPoints() const
const vectorField & faceAreas() const
const cellList & cells() const
int myProcNo() const noexcept
Return processor number.
Original code supplied by Henning Scheufler, DLR (2019)
virtual bool write(const bool valid=true) const
Write using setting from DB.
splitCell * master() const
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
static word printf(const char *fmt, const PrimitiveType &val)
Use a printf-style formatter for a primitive.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const polyBoundaryMesh & patches
IOdictionary porosityProperties(IOobject("porosityProperties", runTime.constant(), runTime, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE))
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
Calculate the gradient of the given field.
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField.
#define DebugInfo
Report an information message using Foam::Info.
#define DebugInFunction
Report an information message using Foam::Info.
const dimensionedScalar c
Speed of light in a vacuum.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
dimensionedScalar pos0(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
List< cell > cellList
A List of cells.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
messageStream Info
Information stream (stdout output on master, null elsewhere)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
vector point
Point is a vector.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
dimensionedScalar neg0(const dimensionedScalar &ds)
Type gMin(const FieldField< Field, Type > &f)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Type gMax(const FieldField< Field, Type > &f)
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
constexpr char nl
The newline '\n' character (0x0a)
#define addProfilingInFunction(name)
#define forAll(list, i)
Loop across all elements in list.