63 return (idx.
x() >= 0 && idx.
y() >= 0 && idx.
z() >= 0);
76 {
"Aw",
"surface area per unit volume" },
83 {
"Bv",
"area blockage" },
84 {
"B",
"area blockage" },
86 {
"Blong",
"longitudinal blockage" },
108static constexpr unsigned outputPrecision = 8;
112 double brs,
double brr,
bool blocked,
113 double surr_br,
double surr_dr,
114 scalar* drags_p, scalar* dragr_p,
121void write_scalarField
124 const scalar& deflt,
const scalarMinMax& limits,
const char *wall_bc,
130void write_uniformField
132 const word& fieldName,
const scalar& deflt,
const char *wall_bc,
145void write_symmTensorField
154void write_symmTensorFieldV
163void write_blocked_face_list
172 double limit_par,
const fileName& casepath
175void write_blockedCellsSet
186static inline scalar averageSurrounding
195 mat(i,j) + mat(i,j+1) + mat(i,j+2)
196 + mat(i+1,j) + mat(i+1,j+2)
197 + mat(i+2,j) + mat(i+2,j+1) + mat(i+2,j+2)
207 <<
word(
"uniform") << token::SPACE
208 << val << token::END_STATEMENT <<
nl;
213static void make_header
221 string note = fieldNotes(
object);
223 IOobject::writeBanner(
os);
225 os <<
"FoamFile\n{\n"
227 <<
" format ascii;\n"
228 <<
" class " << clsName <<
";\n";
232 os <<
" note " << note <<
";\n";
235 if (!location.empty())
237 os <<
" location " << location <<
";\n";
240 os <<
" object " <<
object <<
";\n"
243 IOobject::writeDivider(
os) <<
nl;
258 <<
"No PDRblock set" <<
nl
267 const int xdim = faceDims.
x();
268 const int ydim = faceDims.
y();
269 const int zdim = faceDims.
z();
294 Info<<
"Calculating fields" <<
nl;
310 const label maxDim =
cmptMax(surrDims);
330 for (label iz = 0; iz <= zdim; ++iz)
333 (iz == 0 ? 0 : iz == zdim ? zdim - 2 : iz - 1);
348 for (label iz = 0; iz < surrDims.
z(); ++iz)
353 neiBlock(0, iz) = neiBlock(cellDims.
y(), iz);
354 neiDrag(0, iz) = neiDrag(cellDims.
y(), iz);
355 neiBlock(ydim, iz) = neiBlock(1, iz);
356 neiDrag(ydim, iz) = neiDrag(1, iz);
360 neiBlock(0, iz) = neiBlock(1, iz);
361 neiDrag(0, iz) = neiDrag(1, iz);
362 neiBlock(ydim, iz) = neiBlock(cellDims.
y(), iz);
363 neiDrag(ydim, iz) = neiDrag(cellDims.
y(), iz);
371 const scalar cell_vol = pdrBlock.
V(ix,iy,iz);
373 const scalar surr_br = averageSurrounding(neiBlock, iy, iz);
374 const scalar surr_dr = averageSurrounding(neiDrag, iy, iz);
385 &(cbdi(ix,iy,iz).
x()),
402 for (label ix = 0; ix <= xdim; ++ix)
405 (ix == 0 ? 0 : ix == xdim ? xdim - 2 : ix - 1);
419 for (label ix = 0; ix < surrDims.
x(); ++ix)
421 neiBlock(0, ix) = neiBlock(1, ix);
422 neiDrag(0, ix) = neiDrag(1, ix);
423 neiBlock(zdim, ix) = neiBlock(cellDims.
z(), ix);
424 neiDrag(zdim, ix) = neiDrag(cellDims.
z(), ix);
431 const scalar cell_vol = pdrBlock.
V(ix,iy,iz);
433 const scalar surr_br = averageSurrounding(neiBlock, iz, ix);
434 const scalar surr_dr = averageSurrounding(neiDrag, iz, ix);
445 &(cbdi(ix,iy,iz).
y()),
462 for (label iy = 0; iy <= ydim; ++iy)
468 iyy = (iy == 0 ? ydim - 2 : iy == ydim ? 0 : iy - 1);
472 iyy = (iy == 0 ? 0 : iy == ydim ? ydim - 2 : iy - 1);
487 for (label iy = 0; iy < surrDims.
y(); ++iy)
489 neiBlock(0, iy) = neiBlock(1, iy);
490 neiDrag(0, iy) = neiDrag(1, iy);
491 neiBlock(xdim, iy) = neiBlock(cellDims.
x(), iy);
492 neiDrag(xdim, iy) = neiDrag(cellDims.
x(), iy);
499 const scalar cell_vol = pdrBlock.
V(ix,iy,iz);
501 const scalar surr_br = averageSurrounding(neiBlock, ix, iy);
502 const scalar surr_dr = averageSurrounding(neiDrag, ix, iy);
513 &(cbdi(ix,iy,iz).z()),
532 const scalar dx = pdrBlock.
dx(ix);
533 const scalar dy = pdrBlock.
dy(iy);
534 const scalar dz = pdrBlock.
dz(iz);
535 const scalar cell_vol = pdrBlock.
V(ix, iy, iz);
536 const scalar cell_size = pdrBlock.
width(ix, iy, iz);
546 if (cbdi(ix,iy,iz).
x() > maxCT ) { cbdi(ix,iy,iz).x() = maxCT; } ;
547 if (cbdi(ix,iy,iz).
y() > maxCT ) { cbdi(ix,iy,iz).y() = maxCT; } ;
548 if (cbdi(ix,iy,iz).z() > maxCT ) { cbdi(ix,iy,iz).z() = maxCT; } ;
550 surf(ix,iy,iz) /= cell_vol;
556 const scalar vb =
v_block(ix,iy,iz);
567 || (
surf(ix,iy,iz) <= 0.0 )
577 double nn, lobs, lobsMax;
579 if ( nn < 1.0 ) { nn = 1.0; }
582 if ( nn < 1.0 ) { nn = 1.0; }
584 if ( lobs > lobsMax )
590 if ( nn < 1.0 ) { nn = 1.0; }
592 if ( lobs > lobsMax )
611 else if (
v_block(ix,iy,iz) > 1)
649 write_blocked_face_list
656 write_blockedCellsSet
660 casepath,
"blockedCellsSet"
665 "betav", arr.
v_block, 1, {0, pars.cong_max_betav},
"zeroGradient",
676 const scalar cell_vol = pdrBlock.
V(ix, iy, iz);
699 const scalar cell_23 =
::pow(cell_vol, 2.0/3.0);
709 Info<<
"Writing field files" <<
nl;
721 "Aw", arr.
surf, 0, {0, 1000},
"zeroGradient",
725 write_symmTensorField
730 write_symmTensorFieldV
732 "CT", cbdi,
Zero,
"zeroGradient",
741 "Nv", arr.
obs_count, 0, {0, 1000},
"zeroGradient",
745 write_symmTensorFieldV
756 "N", arr.
obs_count, 0, {0, 1000},
"zeroGradient",
760 write_symmTensorFieldV
787 eff_block /= pdrBlock.
width(ix, iy, iz);
793 if (
dirn_block(ix,iy,iz)[cmpt] || eff_block[cmpt] < 0)
801 drag_s(ix,iy,iz).xx() = eff_block.x();
802 drag_s(ix,iy,iz).yy() = eff_block.y();
803 drag_s(ix,iy,iz).zz() = eff_block.z();
805 cbdi(ix,iy,iz).x() = 1.0 / (
betai_inv1(ix,iy,iz).x() + 1.0);
806 cbdi(ix,iy,iz).y() = 1.0 / (
betai_inv1(ix,iy,iz).y() + 1.0);
807 cbdi(ix,iy,iz).z() = 1.0 / (
betai_inv1(ix,iy,iz).z() + 1.0);
809 if (cbdi(ix,iy,iz).z() < 0 || cbdi(ix,iy,iz).z() > 1.0)
812 <<
"beta_i problem. z-betai_inv1=" <<
betai_inv1(ix,iy,iz).z()
813 <<
" beta_i=" << cbdi(ix,iy,iz).z()
831 write_symmTensorField
840 write_symmTensorField
849 write_symmTensorFieldV
857 write_symmTensorFieldV
873 "b", 1.0,
"zeroGradient",
893 "ft", 0,
"zeroGradient",
917 "Xi", 1,
"zeroGradient",
923 "Xp", 1,
"zeroGradient",
929 "GRxp", 0,
"zeroGradient",
935 "GRep", 0,
"zeroGradient",
941 "RPers", 0,
"zeroGradient",
945 write_pU_fields(meshIndexing,
patches, casepath);
966 "combustFlag", 1,
"zeroGradient",
984 "H2OPS", 0,
"zeroGradient",
990 "AIR", 0,
"zeroGradient",
996 "Ydefault", 0,
"zeroGradient",
1002 "eRatio", 1,
"zeroGradient",
1008 "sprayFlag", 1,
"zeroGradient",
1024 calculateAndWrite(*
this, meshIndexing, casepath,
patches);
1030 double brs,
double brr,
bool blocked,
1031 double surr_br,
double surr_dr,
1032 scalar* drags_p, scalar* dragr_p,
1039 scalar br = brr + brs;
1043 if (count > 1.0) { br /= std::sqrt(count); }
1045 const scalar
alpha =
1048 ? (1.0 - 0.5 * br) / (1.0 - br) / (1.0 - br)
1056 const scalar
expon =
1059 ?
min(
max((surr_br / br - 0.25) * 4.0 / 3.0, scalar(0)), scalar(1))
1066 *dragr_p *= alpha_r;
1067 *drags_p *=
::pow(alpha_s, 1.09);
1069 if ( *cbdi_p < 0.0 ) { *cbdi_p = 0.0; }
1073 if ( *drags_p < 0.0 ) { *drags_p = 0.0; }
1077 if ( (surr_dr * 0.25) > *drags_p )
1079 *drags_p = surr_dr * 0.25;
1083 if ( blocked ) { *cbdi_p = 0.0; *drags_p = 0.0; *dragr_p = 0.0; }
1093 <<
"No blockage information - PDRblock is not set" <<
nl;
1100 scalar totCount = 0;
1101 scalar totVolBlock = 0;
1114 totVolBlock += v_block(ijk) * pdrBlock.
V(ijk);
1115 totArea += surf(ijk);
1117 totCount +=
max(0, obs_count(ijk));
1119 totDrag.x() +=
max(0, drag_s(ijk).xx());
1120 totDrag.y() +=
max(0, drag_s(ijk).yy());
1121 totDrag.z() +=
max(0, drag_s(ijk).zz());
1125 totBlock[cmpt] +=
max(0, area_block_s(ijk)[cmpt]);
1126 totBlock[cmpt] +=
max(0, area_block_r(ijk)[cmpt]);
1133 <<
"Volume blockage: " << totVolBlock <<
nl
1134 <<
"Total drag: " << totDrag <<
nl
1135 <<
"Total count: " << totCount <<
nl
1136 <<
"Total area blockage: " << totBlock <<
nl
1137 <<
"Total surface area: " << totArea <<
nl;
1145static void tail_field
1149 const char* wall_bc,
1157 putUniform(
os,
"value", deflt);
1163 const word& patchName =
patches[patchi].patchName;
1174 putUniform(
os,
"value", deflt);
1183 else if (
patches[patchi].patchType == 0)
1188 putUniform(
os,
"value", deflt);
1196 putUniform(
os,
"value", deflt);
1250void write_scalarField
1253 const scalar& deflt,
const scalarMinMax& limits,
const char *wall_bc,
1269 <<
"nonuniform List<scalar>" <<
nl
1272 for (label celli=0; celli < meshIndexing.
nCells(); ++celli)
1276 if (!isGoodIndex(cellIdx))
1296 putUniform(
os,
"inletValue", deflt);
1297 putUniform(
os,
"value", deflt);
1302 tail_field(
os, deflt, wall_bc,
patches);
1312void write_uniformField
1314 const word& fieldName,
const scalar& deflt,
const char *wall_bc,
1328 putUniform(
os,
"internalField", deflt);
1337 if (fieldName ==
"alphat" || fieldName ==
"nut")
1345 putUniform(
os,
"inletValue", deflt);
1348 putUniform(
os,
"value", deflt);
1352 tail_field(
os, deflt, wall_bc,
patches);
1400 for (label patchi = 0; patchi < 3; ++patchi)
1408 for (label patchi = 3; patchi <
patches.
size(); ++patchi)
1411 const word& patchName =
p.patchName;
1413 if (
p.patchType == 0)
1418 os.
writeEntry(
"fileName",
"<case>" / (patchName +
".dat"));
1426 os.
writeEntry(
"type",
"activePressureForceBaffleVelocity");
1489 const char *wall_bc =
"zeroGradient;\n\trho\trho";
1499 putUniform(
os,
"internalField", deflt);
1512 putUniform(
os,
"value", deflt);
1516 tail_field(
os, deflt, wall_bc,
patches);
1527void write_symmTensorField
1529 const word& fieldName,
1531 const symmTensor& deflt,
const char *wall_bc,
1546 <<
"nonuniform List<symmTensor>" <<
nl
1549 for (label celli=0; celli < meshIndexing.
nCells(); ++celli)
1553 if (!isGoodIndex(cellIdx))
1571 putUniform(
os,
"inletValue", deflt);
1572 putUniform(
os,
"value", deflt);
1577 tail_field(
os, deflt, wall_bc,
patches);
1587void write_symmTensorFieldV
1589 const word& fieldName,
1591 const symmTensor& deflt,
const char *wall_bc,
1606 <<
"nonuniform List<symmTensor>" <<
nl
1611 for (label celli=0; celli < meshIndexing.
nCells(); ++celli)
1615 if (!isGoodIndex(cellIdx))
1639 putUniform(
os,
"inletValue", deflt);
1640 putUniform(
os,
"value", deflt);
1645 tail_field(
os, deflt, wall_bc,
patches);
1655void write_blocked_face_list
1663 double limit_par,
const fileName& casepath
1676 for (label facei=0; facei < meshIndexing.
nFaces(); ++facei)
1681 if (!isGoodIndex(faceIdx))
1686 const label ix = faceIdx.
x();
1687 const label iy = faceIdx.
y();
1688 const label iz = faceIdx.
z();
1705 val = face_block(faceIdx).x();
1706 patchId = face_patch(faceIdx).x();
1720 ++n_blocked_faces(ix,iy,iz).x();
1721 sub_count(ix,iy,iz).x() = obs_count(ix,iy,iz);
1727 ++n_blocked_faces(ix-1,iy,iz).x();
1728 sub_count(ix-1,iy,iz).x() = obs_count(ix-1,iy,iz);
1736 val = face_block(faceIdx).y();
1737 patchId = face_patch(faceIdx).y();
1751 ++n_blocked_faces(ix,iy,iz).y();
1752 sub_count(ix,iy,iz).y() = obs_count(ix,iy,iz);
1758 ++n_blocked_faces(ix,iy-1,iz).y();
1759 sub_count(ix,iy-1,iz).y() = obs_count(ix,iy-1,iz);
1767 val = face_block(faceIdx).z();
1768 patchId = face_patch(faceIdx).z();
1782 ++n_blocked_faces(ix,iy,iz).z();
1783 sub_count(ix,iy,iz).z() = obs_count(ix,iy,iz);
1789 ++n_blocked_faces(ix,iy,iz-1).z();
1790 sub_count(ix,iy,iz-1).z() = obs_count(ix,iy,iz-1);
1802 else if (val > limit_par)
1820 if (!
isDir(setsDir))
1830 const word& patchName =
patches[patchi].patchName;
1834 make_header(
os,
"polyMesh/sets",
"faceSet", patchName);
1837 const auto& fnd = usedFaces.
cfind(patchi);
1839 if (fnd.good() && (*fnd).any())
1841 os <<
nl << (*fnd).toc() <<
nl;
1856 OFstream os(casepath /
"system/PDRMeshDict");
1858 make_header(
os,
"system",
"dictionary",
"PDRMeshDict");
1865 const word& patchName =
p.patchName;
1866 const word setName = patchName +
"Set";
1868 if (
p.patchType == 0)
1875 else if (
p.patchType > 0)
1877 panelNames.
append(setName);
1886 const word& patchName =
p.patchName;
1887 const word setName = patchName +
"Set";
1889 if (
p.patchType > 0)
1910void write_blockedCellsSet
1917 const word& listName
1920 if (listName.empty())
1943 for (label celli=0; celli < meshIndexing.
nCells(); ++celli)
1947 if (!isGoodIndex(cellIdx))
1952 if (
fld(cellIdx) < limit_par)
1954 blockedCell.
set(celli);
1960 const label n_bfaces =
cmptSum(blocked);
1968 if (blocked[cmpt] > 1) ++n_bpairs;
1973 Info<<
"block " << celli <<
" from "
1974 << blocked <<
" -> ("
1975 << n_bfaces <<
' ' << n_bpairs
1986 blockedCell.
set(celli);
1992 make_header(
os,
"constant/polyMesh/sets",
"cellSet", listName);
1994 if (blockedCell.
any())
1996 os << blockedCell.
toc();
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
Preparation of fields for PDRFoam.
#define MIN_COUNT_FOR_SIZE
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void append(const T &val)
Copy append an element to the end of this list.
A HashTable similar to std::unordered_map.
bool set(const Key &key, const T &obj)
Copy assign a new entry, overwriting existing entries.
const_iterator cfind(const Key &key) const
Find and return an const_iterator set at the hashed entry.
static Ostream & writeEndDivider(Ostream &os)
Write the standard end file divider.
Generic templated field type with i-j-k addressing.
A HashTable to objects of type <T> with a label key.
const T & clip(const T &val) const
Output to file stream, using an OSstream.
virtual int precision() const
Get precision of output field.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
virtual Ostream & endBlock()
Write end block group.
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
virtual Ostream & writeKeyword(const keyType &kw)
Write the keyword followed by an appropriate indentation.
virtual Ostream & beginBlock(const keyType &kw)
Write begin block group with the given name.
Work array definitions for PDR fields.
void blockageSummary() const
Summary of the blockages.
IjkField< vector > drag_r
Directional drag from round obstacles.
IjkField< vector > betai_inv1
IjkField< vector > along_block
IjkField< scalar > v_block
Volume blockage.
IjkField< vector > area_block_r
Summed area blockage (directional) from round obstacles.
IjkField< vector > sub_count
Number of obstacles parallel to specified direction.
IjkField< vector > area_block_s
Summed area blockage (directional) from sharp obstacles.
const PDRblock & block() const
Reference to PDRblock.
IjkField< scalar > obs_size
Obstacle size in cell.
IjkField< Vector< bool > > dirn_block
A total directional blockage in the cell.
IjkField< scalar > obs_count
Number of obstacles in cell.
IjkField< vector > grating_count
IjkField< symmTensor > drag_s
Tensorial drag from sharp obstacles.
IjkField< vector > face_block
IjkField< labelVector > face_patch
Face field for (directional) for patch Id.
IjkField< scalar > surf
Surface area in cell.
A single block x-y-z rectilinear mesh addressable as i,j,k with simplified creation....
scalar V(const label i, const label j, const label k) const
Cell volume at i,j,k position.
scalar dx(const label i) const
Cell size in x-direction at i position.
scalar dy(const label j) const
Cell size in y-direction at j position.
scalar dz(const label k) const
Cell size in z-direction at k position.
scalar width(const label i, const label j, const label k) const
Characteristic cell size at i,j,k position.
OpenFOAM/PDRblock addressing information.
label nCells() const
The number of cells.
List< labelVector > cellIndex
For each cell, the corresponding i-j-k address.
List< direction > faceOrient
For each face, the x/y/z orientation.
labelVector cellDims
The cell i-j-k addressing range.
List< labelVector > faceIndex
For each face, the corresponding i-j-k address.
label nFaces() const
The number of faces.
labelVector faceDims
The face i-j-k addressing range.
string UPatchBc
"fixedValue;value uniform (0 0 0)"
scalar blockedFacePar
Faces with area blockage greater than this are blocked.
scalar empty_lobs_fac
Lobs in empty cell is this * cube root of cell volume.
scalar maxCR
Upper limit on CR (CT also gets limited)
word groundPatchName
The name for the "ground" patch.
scalar blockedCellPoros
Cells with porosity less than this are blocked.
scalar outerCombFac
Value for outer region.
word outerPatchName
The name for the "outer" patch.
Bookkeeping for patch definitions.
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
static const SymmTensor I
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
label size() const noexcept
The number of elements in the list.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
void set(const bitSet &bitset)
Set specified bits from another bitset.
labelList toc() const
The indices of the on bits as a sorted labelList.
bool any() const
True if any bits in this bitset are set.
Creates a single block of cells from point coordinates, numbers of cells in each direction and an exp...
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
A class for handling file names.
static std::string path(const std::string &str)
Return directory path name (part before last /)
void calculateAndWrite()
Calculate the fields and write.
label size() const
Return the total i*j*k size.
static constexpr direction nComponents
Number of components in bool is 1.
@ END_STATEMENT
End entry [isseparator].
@ BEGIN_LIST
Begin list [isseparator].
@ END_LIST
End list [isseparator].
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionSet dimPressure
const dimensionSet dimViscosity
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
const dimensionSet dimless
Dimensionless.
List< label > labelList
A List of labels.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
const dimensionSet dimVelocity
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
bool isNull(const T *ptr)
True if ptr is a pointer (of type T) to the nullObject.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
bool equal(const T &s1, const T &s2)
Compare two values for equality.
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0)
Foam::PDRparams pars
Globals for program parameters (ugly hack)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
errorManipArg< error, int > exit(error &err, const int errNo=1)
bool isDir(const fileName &name, const bool followLink=true)
Does the name exist as a DIRECTORY in the file system?
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.
static const char *const typeName
The type name used in ensight case files.
A non-counting (dummy) refCount.