55 Foam::isoAdvection::isoAdvection
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())),
105 surf_(reconstructionSchemes::
New(alpha1_, phi_, U_, dict_)),
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)
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." 187 void Foam::isoAdvection::setProcessorPatches()
189 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
190 surfaceCellFacesOnProcPatches_.
clear();
191 surfaceCellFacesOnProcPatches_.resize(
patches.
size());
194 procPatchLabels_.clear();
199 isA<processorPolyPatch>(
patches[patchi])
203 procPatchLabels_.append(patchi);
209 void 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]);
243 void 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();
376 const label facei = bsFaces_[i];
377 const label patchi = boundaryMesh.patchID(facei);
379 if (!
phib[patchi].empty())
381 const label patchFacei = boundaryMesh[patchi].whichFace(facei);
383 const scalar phiP =
phib[patchi][patchFacei];
387 const scalar magSf = magSfb[patchi][patchFacei];
389 dVfb[patchi][patchFacei] = advectFace_.timeIntegratedFaceFlux
402 checkIfOnProcPatch(facei);
408 syncProcPatches(dVf_, phi_);
410 writeIsoFaces(isoFacePts);
412 DebugInfo <<
"Number of isoAdvector surface cells = " 417 void Foam::isoAdvection::setDownwindFaces
420 DynamicLabelList& downwindFaces
426 const labelList& own = mesh_.faceOwner();
428 const cell&
c =
cells[celli];
430 downwindFaces.
clear();
433 for (
const label facei:
c)
436 const scalar
phi = faceValue(phi_, facei);
438 if (own[facei] == celli)
442 downwindFaces.append(facei);
447 downwindFaces.append(facei);
451 downwindFaces.shrink();
455 Foam::scalar Foam::isoAdvection::netFlux
464 const cell&
c = mesh_.cells()[celli];
467 const labelList& own = mesh_.faceOwner();
469 for (
const label facei :
c)
471 const scalar dVff = faceValue(dVf, facei);
473 if (own[facei] == celli)
491 bool returnSyncedFaces
494 DynamicLabelList syncedFaces(0);
495 const polyBoundaryMesh&
patches = mesh_.boundaryMesh();
499 DynamicList<label> neighProcs;
503 for (
const label patchi : procPatchLabels_)
505 const processorPolyPatch& procPatch =
506 refCast<const processorPolyPatch>(
patches[patchi]);
507 const label nbrProci = procPatch.neighbProcNo();
510 neighProcs.push_uniq(nbrProci);
512 const scalarField& pFlux = dVf.boundaryField()[patchi];
513 const List<label>& surfCellFacesOnProcPatch =
514 surfaceCellFacesOnProcPatches_[patchi];
516 const UIndirectList<scalar> dVfPatch
519 surfCellFacesOnProcPatch
522 UOPstream toNbr(nbrProci, pBufs);
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();
538 List<scalar> nbrdVfs;
541 UIPstream fromNbr(nbrProci, pBufs);
542 fromNbr >> faceIDs >> nbrdVfs;
545 if (returnSyncedFaces)
547 List<label> syncedFaceI(faceIDs);
548 for (label& faceI : syncedFaceI)
550 faceI += procPatch.start();
552 syncedFaces.append(syncedFaceI);
557 Pout<<
"Received at time = " << mesh_.time().value()
558 <<
": surfCellFacesOnProcPatch = " << faceIDs <<
nl 559 <<
"Received at time = " << mesh_.time().value()
560 <<
": dVfPatch = " << nbrdVfs <<
endl;
564 scalarField& localFlux = dVf.boundaryFieldRef()[patchi];
568 const label facei = faceIDs[i];
569 localFlux[facei] = - nbrdVfs[i];
570 if (
debug &&
mag(localFlux[facei] + nbrdVfs[i]) > ROOTVSMALL)
572 Pout<<
"localFlux[facei] = " << localFlux[facei]
573 <<
" and nbrdVfs[i] = " << nbrdVfs[i]
574 <<
" for facei = " << facei <<
endl;
582 forAll(procPatchLabels_, patchLabeli)
584 const label patchi = procPatchLabels_[patchLabeli];
585 const scalarField& localFlux = dVf.boundaryField()[patchi];
586 Pout<<
"time = " << mesh_.time().value() <<
": localFlux = " 587 << localFlux <<
endl;
592 forAll(surfaceCellFacesOnProcPatches_, patchi)
594 surfaceCellFacesOnProcPatches_[patchi].clear();
602 void Foam::isoAdvection::checkIfOnProcPatch(
const label facei)
604 if (!mesh_.isInternalFace(facei))
606 const polyBoundaryMesh&
pbm = mesh_.boundaryMesh();
609 if (isA<processorPolyPatch>(
pbm[patchi]) && !
pbm[patchi].empty())
611 const label patchFacei =
pbm[patchi].whichFace(facei);
612 surfaceCellFacesOnProcPatches_[patchi].
append(patchFacei);
618 void Foam::isoAdvection::applyBruteForceBounding()
621 bool alpha1Changed =
false;
623 const scalar snapAlphaTol = dict_.getOrDefault<scalar>(
"snapTol", 0);
624 if (snapAlphaTol > 0)
628 *
pos0(alpha1_ - snapAlphaTol)
629 *
neg0(alpha1_ - (1.0 - snapAlphaTol))
630 +
pos0(alpha1_ - (1.0 - snapAlphaTol));
632 alpha1Changed =
true;
635 if (dict_.getOrDefault(
"clip",
true))
637 alpha1_.clamp_range(zero_one{});
638 alpha1Changed =
true;
643 alpha1_.correctBoundaryConditions();
650 if (!mesh_.time().writeTime())
return;
652 if (dict_.getOrDefault(
"writeSurfCells",
false))
659 mesh_.time().timeName(),
665 cSet.insert(surfCells_);
674 const DynamicList<List<point>>& faces
677 if (!writeIsoFacesToFile_ || !mesh_.time().writeTime())
return;
680 const fileName outputFile
682 mesh_.time().globalPath()
684 /
word::printf(
"isoFaces_%012d.obj", mesh_.time().timeIndex())
696 mkDir(outputFile.path());
697 OBJstream
os(outputFile);
698 Info<<
nl <<
"isoAdvection: writing iso faces to file: " 703 const DynamicList<List<point>>& procFacePts
707 for (
const List<point>& facePts : procFacePts)
716 mkDir(outputFile.path());
717 OBJstream
os(outputFile);
718 Info<<
nl <<
"isoAdvection: writing iso faces to file: " 721 for (
const List<point>& facePts : faces)
Ostream & writeFace(const UList< point > &points, const bool lines=true)
Write face loop points with lines/filled-polygon.
const polyBoundaryMesh & pbm
Surface integrate surfaceField creating a volField. Surface sum a surfaceField creating a volField...
const Internal::FieldType & primitiveField() const noexcept
Return a const-reference to the internal field values.
List< cell > cellList
List of cell.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Type gMin(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
const labelList & patchID() const
Per boundary face label the patch index.
void writeIsoFaces(const DynamicList< List< point >> &isoFacePts) const
Write isoface points to .obj file.
static bool & parRun() noexcept
Test if this a parallel run.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
const dimensionSet dimVol(dimVolume)
Older spelling for dimVolume.
const cellList & cells() const
GeometricBoundaryField< scalar, fvsPatchField, surfaceMesh > Boundary
Type of boundary fields.
Ignore writing from objectRegistry::writeObject()
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
const Time & time() const
Return the top-level database.
GeometricField< vector, fvPatchField, volMesh > volVectorField
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
conserve primitiveFieldRef()+
Macros for easy insertion into run-time selection tables.
#define forAll(list, i)
Loop across all elements in list.
void writeSurfaceCells() const
Return cellSet of surface cells.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
static void gatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
Gather data, but keep individual values separate. Uses the specified communication schedule...
const dimensionedScalar e
Elementary charge.
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Calculate the gradient of the given field.
void clear()
Clear the list, i.e. set size to zero.
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
Type * getObjectPtr(const word &name, const bool recursive=false) const
Return non-const pointer to the object of the given Type, using a const-cast to have it behave like a...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
#define DebugInFunction
Report an information message using Foam::Info.
label size() const noexcept
The number of entries in the list.
dimensionedScalar neg0(const dimensionedScalar &ds)
Reading is optional [identical to LAZY_READ].
const vectorField & cellCentres() const
#define DebugInfo
Report an information message using Foam::Info.
static word printf(const char *fmt, const PrimitiveType &val)
Use a printf-style formatter for a primitive.
const surfaceScalarField & magSf() const
Return cell face area magnitudes.
dimensionedScalar pos0(const dimensionedScalar &ds)
void clear()
Clear the patch list and all demand-driven data.
const word & constant() const noexcept
Return constant name.
int debug
Static debugging option.
Type gMax(const FieldField< Field, Type > &f)
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
IOdictionary porosityProperties(IOobject("porosityProperties", runTime.constant(), runTime, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE))
#define addProfilingInFunction(name)
Define profiling trigger with specified name and description corresponding to the compiler-defined fu...
virtual const fileName & name() const
Read/write access to the name of the stream.
const vectorField & faceCentres() const
void append(autoPtr< T > &ptr)
Move append an element to the end of the list.
vector point
Point is a vector.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const vectorField & faceAreas() const
const dimensionedScalar c
Speed of light in a vacuum.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
const polyBoundaryMesh & patches
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
"nonBlocking" : (MPI_Isend, MPI_Irecv)
messageStream Info
Information stream (stdout output on master, null elsewhere)
List< label > labelList
A List of labels.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
const labelListList & cellCells() const
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const labelListList & cellPoints() const
const scalarField & cellVolumes() const
static constexpr const zero Zero
Global zero (0)
const volScalarField & alpha1