diff --git a/applications/utilities/preProcessing/PDRsetFields/Make/files b/applications/utilities/preProcessing/PDRsetFields/Make/files
new file mode 100644
index 0000000000..af1c28d4ca
--- /dev/null
+++ b/applications/utilities/preProcessing/PDRsetFields/Make/files
@@ -0,0 +1,19 @@
+PDRsetFields.C
+
+PDRarrays.C
+PDRarraysAnalyse.C
+PDRarraysCalc.C
+PDRmeshArrays.C
+PDRparams.C
+PDRpatchDef.C
+PDRlegacyMeshSpec.C
+PDRutilsIntersect.C
+PDRutilsOverlap.C
+
+obstacles/PDRobstacle.C
+obstacles/PDRobstacleIO.C
+obstacles/PDRobstacleTypes.C
+obstacles/PDRobstacleLegacyIO.C
+obstacles/PDRobstacleLegacyRead.C
+
+EXE = $(FOAM_APPBIN)/PDRsetFields
diff --git a/applications/utilities/preProcessing/PDRsetFields/Make/options b/applications/utilities/preProcessing/PDRsetFields/Make/options
new file mode 100644
index 0000000000..eb9d10ceb0
--- /dev/null
+++ b/applications/utilities/preProcessing/PDRsetFields/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+ -Iobstacles \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude \
+ -I$(LIB_SRC)/surfMesh/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/mesh/blockMesh/lnInclude
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -lfileFormats \
+ -lsurfMesh \
+ -lmeshTools \
+ -lblockMesh
diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRarrays.C b/applications/utilities/preProcessing/PDRsetFields/PDRarrays.C
new file mode 100644
index 0000000000..f15cd52fd6
--- /dev/null
+++ b/applications/utilities/preProcessing/PDRsetFields/PDRarrays.C
@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016 Shell Research Ltd.
+ Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRarrays.H"
+#include "PDRblock.H"
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Smaller helper to resize matrix and assign to Zero
+template
+inline void resizeMatrix(SquareMatrix& mat, const label n)
+{
+ mat.setSize(n);
+ mat = Zero;
+}
+
+
+// Smaller helper to resize i-j-k field and assign to uniform value,
+// normally Zero
+template
+inline void resizeField
+(
+ IjkField& fld,
+ const labelVector& ijk,
+ const T& val = T(Zero)
+)
+{
+ fld.resize(ijk);
+ fld = val;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::PDRarrays::PDRarrays()
+:
+ pdrBlock_(std::cref(PDRblock::null()))
+{}
+
+
+Foam::PDRarrays::PDRarrays(const PDRblock& pdrBlock)
+:
+ PDRarrays()
+{
+ reset(pdrBlock);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::PDRarrays::reset(const PDRblock& pdrBlock)
+{
+ pdrBlock_ = std::cref(pdrBlock);
+
+ // Resize all the major arrays, which are grouped in the structure arrp
+ // All the relevant dimensions are in PDRblock
+
+ // Cell-based addressing
+ const labelVector cellDims = pdrBlock.sizes();
+
+ // Face or point-based addressing
+ const labelVector faceDims(cellDims + labelVector::one);
+
+ // Max addressing dimensions for 2D arrays, with some extra space
+ // These will be used for any combination of x,y,z,
+ // so need to be dimensioned to the maximum size in both directions
+ const label maxDim = cmptMax(pdrBlock.sizes()) + 2;
+
+ resizeField(v_block, cellDims);
+ resizeField(surf, cellDims);
+
+ resizeField(area_block_s, cellDims);
+ resizeField(area_block_r, cellDims);
+ resizeField(dirn_block, cellDims);
+
+ resizeField(face_block, faceDims);
+
+ resizeField(along_block, cellDims);
+
+ resizeField(betai_inv1, cellDims);
+
+ resizeField(obs_count, cellDims);
+ resizeField(sub_count, cellDims);
+ resizeField(grating_count, cellDims);
+
+ resizeField(drag_s, cellDims);
+ resizeField(drag_r, cellDims);
+
+ resizeField(obs_size, cellDims);
+
+ for (auto& list : overlap_1d)
+ {
+ list.resize(maxDim);
+ list = Zero;
+ }
+
+ resizeMatrix(aboverlap, maxDim);
+ resizeMatrix(abperim, maxDim);
+ resizeMatrix(a_lblock, maxDim);
+ resizeMatrix(b_lblock, maxDim);
+ resizeMatrix(ac_lblock, maxDim);
+ resizeMatrix(bc_lblock, maxDim);
+ resizeMatrix(c_count, maxDim);
+ resizeMatrix(c_drag, maxDim);
+
+ resizeField(face_patch, faceDims, labelVector::uniform(-1));
+ resizeField(hole_in_face, faceDims);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRarrays.H b/applications/utilities/preProcessing/PDRsetFields/PDRarrays.H
new file mode 100644
index 0000000000..908b45a37f
--- /dev/null
+++ b/applications/utilities/preProcessing/PDRsetFields/PDRarrays.H
@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016 Shell Research Ltd.
+ Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Class
+ Foam::PDRarrays
+
+Description
+ Work array definitions for PDR fields
+
+SourceFiles
+ PDRarrays.C
+ PDRarraysCalc.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRarrays_H
+#define PDRarrays_H
+
+#include "symmTensor.H"
+#include "symmTensor2D.H"
+#include "SquareMatrix.H"
+#include "IjkField.H"
+#include
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class PDRblock;
+class PDRmeshArrays;
+class PDRobstacle;
+class PDRpatchDef;
+
+
+/*---------------------------------------------------------------------------*\
+ Class PDRarrays Declaration
+\*---------------------------------------------------------------------------*/
+
+class PDRarrays
+{
+ //- Reference to PDRblock
+ std::reference_wrapper pdrBlock_;
+
+public:
+
+ // Data Members
+ // Entries used for analysis and when writing fields
+
+ //- Volume blockage
+ IjkField v_block;
+
+ //- Surface area in cell
+ IjkField surf;
+
+ //- Obstacle size in cell
+ IjkField obs_size;
+
+ //- Summed area blockage (directional) from sharp obstacles
+ IjkField area_block_s;
+
+ //- Summed area blockage (directional) from round obstacles
+ IjkField area_block_r;
+
+ //- A total directional blockage in the cell
+ IjkField> dirn_block;
+
+ //- Face area blockage for face,
+ //- summed from cell centre-plane to cell centre-plane
+ IjkField face_block;
+
+ //- Longitudinal area blockage from obstacles that extend all the way
+ //- through the cell in a given direction.
+ IjkField along_block;
+
+ IjkField betai_inv1;
+
+ //- Number of obstacles in cell.
+ // Can be non-integer if an obstacle does not pass all way through cell
+ IjkField obs_count;
+
+ //- Number of obstacles parallel to specified direction
+ IjkField sub_count;
+
+ //- Addition to count to account for grating comprises many bars
+ //- (to get Lobs right)
+ IjkField grating_count;
+
+ //- Tensorial drag from sharp obstacles
+ IjkField drag_s;
+
+ //- Directional drag from round obstacles
+ IjkField drag_r;
+
+
+ // Next arrays are for 2D calculations of intersection
+
+ // One-dimensional scratch areas for cell overlaps
+ Vector> overlap_1d;
+
+ // In two dimensions, area of cell covered by circle
+ SquareMatrix aboverlap;
+
+ // In two dimensions, length of perimeter of circle witthin cell
+ SquareMatrix abperim;
+
+ // For offset cells, i.e. face blockage
+ SquareMatrix a_lblock, b_lblock;
+
+ // For centred cells
+ SquareMatrix ac_lblock, bc_lblock;
+
+ // The count in the cells
+ SquareMatrix c_count;
+
+ //- Cell-centred drag
+ SquareMatrix c_drag;
+
+ //- Face field for (directional) for patch Id
+ IjkField face_patch;
+
+ //- Face field for (directional) hole in face
+ IjkField> hole_in_face;
+
+
+ // Constructors
+
+ //- Construct null
+ PDRarrays();
+
+ //- Construct and reset
+ explicit PDRarrays(const PDRblock& pdrBlock);
+
+
+ //- Destructor
+ ~PDRarrays() = default;
+
+
+ // Member Functions
+
+ //- Reset PDRblock reference, resize and zero arrays
+ void reset(const PDRblock& pdrBlock);
+
+ //- Reference to PDRblock
+ const PDRblock& block() const
+ {
+ return pdrBlock_.get();
+ }
+
+ //- Summary of the blockages
+ // For diagnostics and general overview
+ void blockageSummary() const;
+
+ //- Add cylinder blockage
+ void addCylinder(const PDRobstacle& obs);
+
+ //- Add general (non-cylinder) blockage
+ void addBlockage
+ (
+ const PDRobstacle& obs,
+ DynamicList& patches,
+ const int volumeSign
+ );
+
+
+ static void calculateAndWrite
+ (
+ PDRarrays& arr,
+ const PDRmeshArrays& meshIndexing,
+ const fileName& casepath,
+ const UList& patches
+ );
+
+ void calculateAndWrite
+ (
+ const fileName& casepath,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches
+ );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRarraysAnalyse.C b/applications/utilities/preProcessing/PDRsetFields/PDRarraysAnalyse.C
new file mode 100644
index 0000000000..a72c769ab1
--- /dev/null
+++ b/applications/utilities/preProcessing/PDRsetFields/PDRarraysAnalyse.C
@@ -0,0 +1,1099 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016 Shell Research Ltd.
+ Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRsetFields.H"
+#include "PDRobstacle.H"
+#include "PDRpatchDef.H"
+#include "PDRutils.H"
+#include "PDRutilsInternal.H"
+#include "ListOps.H"
+
+#include
+#include
+#include
+
+#ifndef FULLDEBUG
+#define NDEBUG
+#endif
+#include
+
+using namespace Foam;
+using namespace Foam::PDRutils;
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+// Cell blockage.
+//
+// Use simple sum, because this will eventually be divided by a notional
+// number of rows to give a per-row blockage ratio
+//
+// b is the (fractional) area blockage. f is 1 or 0.5 to split between ends
+// Thus if b isclose to 1, the obstacle is totally blocking the cell in this direction,
+// and we could modify the behaviour if we wish.
+inline static void add_blockage_c
+(
+ scalar& a,
+ bool& blocked,
+ const scalar b,
+ const scalar f = 1.0
+)
+{
+ a += b * f;
+ if (b > pars.blockageNoCT)
+ {
+ blocked = true;
+ }
+}
+
+
+// Face blockage
+//
+// Adds more area blockage to existing amount by assuming partial overlap,
+// i.e. multiplying porosities.
+//
+// Simple addition if the existing amount is negative, because negative
+// blocks (summed first) should just cancel out part of positive blocks.
+inline static void add_blockage_f
+(
+ scalar& a,
+ const scalar b,
+ bool isHole
+)
+{
+ if (a > 0.0)
+ {
+ // Both positive
+ a = 1.0 - (1.0 - a) * (1.0 - b);
+ }
+ else if (b < pars.blockedFacePar || isHole)
+ {
+ // Add until it eventually becomes positive
+ a += b;
+ }
+ else
+ {
+ // If one obstacle blocks face, face is blocked, regardless of
+ // overlap calculations, unless an input negative obstacle makes a
+ // hole in it
+ a = b;
+ }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::PDRarrays::addCylinder(const PDRobstacle& obs)
+{
+ if (equal(obs.vbkge, 0))
+ {
+ return;
+ }
+
+ if (isNull(block()))
+ {
+ FatalErrorInFunction
+ << "No PDRblock set" << nl
+ << exit(FatalError);
+ }
+
+ const PDRblock& pdrBlock = block();
+ const PDRblock::location& xgrid = pdrBlock.grid().x();
+ const PDRblock::location& ygrid = pdrBlock.grid().y();
+ const PDRblock::location& zgrid = pdrBlock.grid().z();
+
+ scalarList& xoverlap = overlap_1d.x();
+ scalarList& yoverlap = overlap_1d.y();
+ scalarList& zoverlap = overlap_1d.z();
+
+ int cxmin, cxmax, cymin, cymax, czmin, czmax;
+ int cfxmin, cfxmax, cfymin, cfymax, cfzmin, cfzmax;
+
+ scalar rad_a, rad_b;
+ vector area_block(Zero);
+
+ if (obs.typeId == PDRobstacle::CYLINDER)
+ {
+ rad_a = rad_b = 0.5*obs.dia();
+ }
+ else
+ {
+ rad_a = 0.5*(obs.wa * ::cos(obs.theta()) + obs.wb * ::sin(obs.theta()));
+ rad_b = 0.5*(obs.wa * ::sin(obs.theta()) + obs.wb * ::cos(obs.theta()));
+ }
+
+ switch (obs.orient)
+ {
+ case vector::Z:
+ {
+ // Determine the part of the grid (potentially) covered by this obstacle.
+ one_d_overlap
+ (
+ obs.x() - rad_a, obs.x() + rad_a,
+ pdrBlock.grid().x(),
+ xoverlap, &cxmin, &cxmax, &cfxmin, &cfxmax
+ ); assert(cxmax >=0);
+
+ one_d_overlap
+ (
+ obs.y() - rad_b, obs.y() + rad_b,
+ pdrBlock.grid().y(),
+ yoverlap, &cymin, &cymax, &cfymin, &cfymax
+ ); assert(cymax >=0);
+
+ one_d_overlap
+ (
+ obs.z(), obs.z() + obs.len(),
+ pdrBlock.grid().z(),
+ zoverlap, &czmin, &czmax, &cfzmin, &cfzmax
+ ); assert(czmax >=0);
+
+ // The area of intersection with each 2D cell in an x-y plane.
+ // a corresponds to x, and b to y.
+ circle_overlap
+ (
+ obs.x(), obs.y(), obs.dia(),
+ obs.theta(), obs.wa, obs.wb,
+ pdrBlock.grid().x(), cxmin, cfxmax,
+ pdrBlock.grid().y(), cymin, cfymax,
+ aboverlap, abperim, a_lblock, ac_lblock, c_count, c_drag,
+ b_lblock, bc_lblock
+ );
+
+
+ for (label ix = cxmin; ix <= cfxmax; ix++)
+ {
+ for (label iy = cymin; iy <= cfymax; iy++)
+ {
+ for (label iz = czmin; iz <= cfzmax; iz++)
+ {
+ const scalar vol_block = aboverlap(ix,iy) * zoverlap[iz];
+ v_block(ix,iy,iz) += vol_block;
+ surf(ix,iy,iz) += abperim(ix,iy) * zoverlap[iz] * zgrid.width(iz);
+
+ // In the 2D case, the ends of the cylinder appear to
+ // be in the cell even when not, making surf and
+ // obs_size wrong. So leave out ends.
+
+ if (!pars.two_d && (iz == czmin || iz == czmax))
+ {
+ // End cells
+ const scalar both_ends_fac = (czmin == czmax ? 2.0 : 1.0);
+
+ surf(ix,iy,iz) += aboverlap(ix,iy)
+ * xgrid.width(ix) * ygrid.width(iy) * both_ends_fac;
+ }
+
+ const scalar temp = c_count(ix,iy) * zoverlap[iz];
+
+ obs_count(ix,iy,iz) += temp;
+ sub_count(ix,iy,iz).z() += temp;
+
+ if (!pars.two_d && (iz == czmin || iz == czmax))
+ {
+ // End faces
+ const scalar both_ends_fac = (czmin == czmax ? 2.0 : 1.0);
+
+ sub_count(ix,iy,iz).z() -= temp * both_ends_fac / 2.0;
+ }
+
+ // Keep the blockage and drag of round obst separate
+ // from the sharp for the moment because different
+ // blockage ratio corrections will be needed later.
+ //
+ // Only the relevant diagonal elements of drag
+ // are stored; other are zero.
+
+ area_block.x() = ac_lblock(ix,iy) * zoverlap[iz];
+ area_block.y() = bc_lblock(ix,iy) * zoverlap[iz];
+
+ // Do not want blockage and drag across the end
+ // except at perimeter.
+ if (aboverlap(ix,iy) < pars.blockedFacePar)
+ {
+ if (obs.typeId == PDRobstacle::CYLINDER)
+ {
+ drag_r(ix,iy,iz).x() += c_drag(ix,iy).xx() * zoverlap[iz] / xgrid.width(ix) / ygrid.width(iy);
+ drag_r(ix,iy,iz).y() += c_drag(ix,iy).yy() * zoverlap[iz] / xgrid.width(ix) / ygrid.width(iy);
+
+ add_blockage_c(area_block_r(ix,iy,iz).x(), dirn_block(ix,iy,iz).x(), area_block.x(), 1.0);
+ add_blockage_c(area_block_r(ix,iy,iz).y(), dirn_block(ix,iy,iz).y(), area_block.y(), 1.0);
+ }
+ else
+ {
+ drag_s(ix,iy,iz).xx() += c_drag(ix,iy).xx() * zoverlap[iz] / xgrid.width(ix) / ygrid.width(iy);
+ drag_s(ix,iy,iz).yy() += c_drag(ix,iy).yy() * zoverlap[iz] / xgrid.width(ix) / ygrid.width(iy);
+
+ add_blockage_c(area_block_s(ix,iy,iz).x(), dirn_block(ix,iy,iz).x(), area_block.x(), 1.0);
+ add_blockage_c(area_block_s(ix,iy,iz).y(), dirn_block(ix,iy,iz).y(), area_block.y(), 1.0);
+ }
+ }
+ // Here we accumulate 1/betai - 1.
+ betai_inv1(ix,iy,iz).x() += vol_block / (1.0 - area_block.x() + floatSMALL);
+ betai_inv1(ix,iy,iz).y() += vol_block / (1.0 - area_block.y() + floatSMALL);
+ betai_inv1(ix,iy,iz).z() += vol_block / (1.0 - aboverlap(ix,iy) + floatSMALL);
+
+ // The off-diagonal elements of drag are stored in "drag" for BOTH round and sharp ostacles
+ drag_s(ix,iy,iz).xy() += c_drag(ix,iy).xy() * zoverlap[iz] / xgrid.width(ix) / ygrid.width(iy);
+
+ add_blockage_f(face_block(ix,iy,iz).x(), a_lblock(ix,iy) * zoverlap[iz], hole_in_face(ix,iy,iz).x());
+ add_blockage_f(face_block(ix,iy,iz).y(), b_lblock(ix,iy) * zoverlap[iz], hole_in_face(ix,iy,iz).y());
+ }
+
+ // Face blockage in the z direction
+ if (czmin == czmax)
+ {
+ // Does not span cell. Block nearest face.
+ add_blockage_f(face_block(ix,iy,cfzmax).z(), aboverlap(ix,iy), hole_in_face(ix,iy,cfzmax).z());
+ }
+ else
+ {
+ // In at least two cells.
+ // Block first and last faces overlapped
+ add_blockage_f(face_block(ix,iy,czmin+1).z(), aboverlap(ix,iy), hole_in_face(ix,iy,czmin+1).z());
+ if (czmax > czmin+1)
+ {
+ add_blockage_f(face_block(ix,iy,czmax).z(), aboverlap(ix,iy), hole_in_face(ix,iy,czmax).z());
+ }
+ }
+
+ // z_block is used to work out the blockage ratio for each
+ // "row" of sub-grid obstacles so this cylinder should
+ // not have any eeffct in cells that it completely spans.
+ // Hence statement commented out below and replaced by
+ // two lines after this loop. That longitudinal clockage
+ // goes into new array z_aling_block
+
+ for (label iz = czmin+1; iz < czmax; ++iz)
+ {
+ // Internal only
+ along_block(ix,iy,iz).z() += aboverlap(ix,iy);
+ }
+
+ // Longitudinal drag only applies at ends of cylinder.
+ // If cylinder spans more than 1 cell, apply half at each
+ // end.
+
+ drag_s(ix,iy,czmin).zz() += aboverlap(ix,iy) * xgrid.width(ix) * ygrid.width(iy) / 2.0;
+ drag_s(ix,iy,czmax).zz() += aboverlap(ix,iy) * xgrid.width(ix) * ygrid.width(iy) / 2.0;
+
+ add_blockage_c(area_block_s(ix,iy,czmin).z(), dirn_block(ix,iy,czmin).z(), aboverlap(ix,iy), 0.5);
+ add_blockage_c(area_block_s(ix,iy,czmax).z(), dirn_block(ix,iy,czmax).z(), aboverlap(ix,iy), 0.5);
+ }
+ }
+
+ break;
+ }
+
+ case vector::X: // orientation
+ {
+ // x-direction cylinder. a,b are y,z.
+ one_d_overlap
+ (
+ obs.x(), obs.x() + obs.len(),
+ pdrBlock.grid().x(),
+ xoverlap, &cxmin, &cxmax, &cfxmin, &cfxmax
+ ); assert(cxmax >=0);
+
+ one_d_overlap
+ (
+ obs.y() - rad_a, obs.y() + rad_a,
+ pdrBlock.grid().y(),
+ yoverlap, &cymin, &cymax, &cfymin, &cfymax
+ ); assert(cymax >=0);
+
+ one_d_overlap
+ (
+ obs.z() - rad_b, obs.z() + rad_b,
+ pdrBlock.grid().z(),
+ zoverlap, &czmin, &czmax, &cfzmin, &cfzmax
+ ); assert(czmax >=0);
+
+ circle_overlap
+ (
+ obs.y(), obs.z(), obs.dia(),
+ obs.theta(), obs.wa, obs.wb,
+ pdrBlock.grid().y(), cymin, cfymax,
+ pdrBlock.grid().z(), czmin, cfzmax,
+ aboverlap, abperim, a_lblock, ac_lblock, c_count, c_drag,
+ b_lblock, bc_lblock
+ );
+
+
+ for (label iy = cymin; iy <= cfymax; iy++)
+ {
+ for (label iz = czmin; iz <= cfzmax; iz++)
+ {
+ for (label ix = cxmin; ix <= cxmax; ix++)
+ {
+ const scalar vol_block = aboverlap(iy,iz) * xoverlap[ix];
+ v_block(ix,iy,iz) += vol_block;
+ surf(ix,iy,iz) += abperim(iy,iz) * xoverlap[ix] * xgrid.width(ix);
+
+ if (ix == cxmin || ix == cxmax)
+ {
+ // End cells
+ const scalar both_ends_fac = (cxmin == cxmax ? 2.0 : 1.0);
+
+ surf(ix,iy,iz) += aboverlap(iy,iz)
+ * ygrid.width(iy) * zgrid.width(iz) * both_ends_fac;
+ }
+
+ const scalar temp = c_count(iy,iz) * xoverlap[ix];
+ obs_count(ix,iy,iz) += temp;
+ sub_count(ix,iy,iz).x() += temp;
+
+ if (ix == cfxmin || ix == cfxmax)
+ {
+ // End faces
+ const scalar both_ends_fac = (cfxmin == cfxmax ? 2.0 : 1.0);
+
+ sub_count(ix,iy,iz).x() -= temp * both_ends_fac / 2.0;
+ }
+
+ area_block.y() = ac_lblock(iy,iz) * xoverlap[ix];
+ area_block.z() = bc_lblock(iy,iz) * xoverlap[ix];
+
+ // Do not want blockage and drag across the end
+ // except at perimeter.
+ if (aboverlap(iy,iz) < pars.blockedFacePar)
+ {
+ if (obs.typeId == PDRobstacle::CYLINDER)
+ {
+ drag_r(ix,iy,iz).y() += c_drag(iy,iz).xx() * xoverlap[ix] / ygrid.width(iy) / zgrid.width(iz);
+ drag_r(ix,iy,iz).z() += c_drag(iy,iz).yy() * xoverlap[ix] / ygrid.width(iy) / zgrid.width(iz);
+
+ add_blockage_c(area_block_r(ix,iy,iz).y(), dirn_block(ix,iy,iz).y(), area_block.y(), 1.0);
+ add_blockage_c(area_block_r(ix,iy,iz).z(), dirn_block(ix,iy,iz).z(), area_block.z(), 1.0);
+ }
+ else
+ {
+ drag_s(ix,iy,iz).yy() += c_drag(iy,iz).xx() * xoverlap[ix] / ygrid.width(iy) / zgrid.width(iz);
+ drag_s(ix,iy,iz).zz() += c_drag(iy,iz).yy() * xoverlap[ix] / ygrid.width(iy) / zgrid.width(iz);
+
+ add_blockage_c(area_block_s(ix,iy,iz).y(), dirn_block(ix,iy,iz).y(), area_block.y(), 1.0);
+ add_blockage_c(area_block_s(ix,iy,iz).z(), dirn_block(ix,iy,iz).z(), area_block.z(), 1.0);
+ }
+ }
+ betai_inv1(ix,iy,iz).y() += vol_block / (1.0 - area_block.y() + floatSMALL);
+ betai_inv1(ix,iy,iz).z() += vol_block / (1.0 - area_block.z() + floatSMALL);
+ betai_inv1(ix,iy,iz).x() += vol_block / (1.0 - aboverlap(iy,iz) + floatSMALL);
+
+ // The off-diagonal elements of drag are stored in "drag" for BOTH round and sharp ostacles
+ drag_s(ix,iy,iz).yz() += c_drag(iy,iz).xy() * xoverlap[ix] / ygrid.width(iy) / zgrid.width(iz);
+
+ add_blockage_f(face_block(ix,iy,iz).y(), a_lblock(iy,iz) * xoverlap[ix], hole_in_face(ix,iy,iz).y());
+ add_blockage_f(face_block(ix,iy,iz).z(), b_lblock(iy,iz) * xoverlap[ix], hole_in_face(ix,iy,iz).z());
+ }
+ if (cxmin == cxmax)
+ {
+ // Does not span cell. Block nearest face.
+ add_blockage_f(face_block(cfxmax,iy,iz).x(), aboverlap(iy,iz), hole_in_face(cfxmax,iy,iz).x());
+ }
+ else
+ {
+ // In at least two cells.
+ // Block first and last faces overlapped
+ add_blockage_f(face_block(cxmin+1,iy,iz).x(), aboverlap(iy,iz), hole_in_face(cxmin+1,iy,iz).x());
+ if (cxmax > cxmin+1)
+ {
+ add_blockage_f(face_block(cxmax,iy,iz).x(), aboverlap(iy,iz), hole_in_face(cxmax,iy,iz).x());
+ }
+ }
+
+ for (label ix = cxmin+1; ix < cxmax; ++ix)
+ {
+ // Internal only
+ along_block(ix,iy,iz).x() += aboverlap(iy,iz);
+ }
+ drag_s(cxmin,iy,iz).xx() += aboverlap(iy,iz) * ygrid.width(iy) * zgrid.width(iz) / 2.0;
+ drag_s(cxmax,iy,iz).xx() += aboverlap(iy,iz) * ygrid.width(iy) * zgrid.width(iz) / 2.0;
+
+ add_blockage_c(area_block_s(cxmin,iy,iz).x(), dirn_block(cxmin,iy,iz).x(), aboverlap(iy,iz), 0.5);
+ add_blockage_c(area_block_s(cxmax,iy,iz).x(), dirn_block(cxmax,iy,iz).x(), aboverlap(iy,iz), 0.5);
+ }
+ }
+
+ break;
+ }
+
+ case vector::Y: // orientation
+ {
+ // y-direction cylinder. a,b are z,x.
+ one_d_overlap
+ (
+ obs.x() - rad_b, obs.x() + rad_b,
+ pdrBlock.grid().x(),
+ xoverlap, &cxmin, &cxmax, &cfxmin, &cfxmax
+ ); assert(cxmax >=0);
+
+ one_d_overlap
+ (
+ obs.y(), obs.y() + obs.len(),
+ pdrBlock.grid().y(),
+ yoverlap, &cymin, &cymax, &cfymin, &cfymax
+ ); assert(cymax >=0);
+
+ one_d_overlap
+ (
+ obs.z() - rad_a, obs.z() + rad_a,
+ pdrBlock.grid().z(),
+ zoverlap, &czmin, &czmax, &cfzmin, &cfzmax
+ ); assert(czmax >=0);
+
+ circle_overlap
+ (
+ obs.z(), obs.x(), obs.dia(),
+ obs.theta(), obs.wa, obs.wb,
+ pdrBlock.grid().z(), czmin, cfzmax,
+ pdrBlock.grid().x(), cxmin, cfxmax,
+ aboverlap, abperim, a_lblock, ac_lblock, c_count, c_drag,
+ b_lblock, bc_lblock
+ );
+
+
+ for (label iz = czmin; iz <= cfzmax; iz++)
+ {
+ for (label ix = cxmin; ix <= cfxmax; ix++)
+ {
+ for (label iy = cymin; iy <= cymax; iy++)
+ {
+ const scalar vol_block = aboverlap(iz,ix) * yoverlap[iy];
+ v_block(ix,iy,iz) += vol_block;
+ surf(ix,iy,iz) += abperim(iz,ix) * yoverlap[iy] * ygrid.width(iy);
+
+ if (iy == cymin || iy == cymax)
+ {
+ // End cells
+ const scalar both_ends_fac = (cymin == cymax ? 2.0 : 1.0);
+
+ surf(ix,iy,iz) += aboverlap(iz,ix)
+ * zgrid.width(iz) * xgrid.width(ix) * both_ends_fac;
+ }
+
+ const scalar temp = c_count(iz,ix) * yoverlap[iy];
+
+ obs_count(ix,iy,iz) += temp;
+ sub_count(ix,iy,iz).y() += temp;
+
+ if (iy == cfymin || iy == cfymax)
+ {
+ // End faces
+ const scalar both_ends_fac = (cfymin == cfymax ? 2.0 : 1.0);
+
+ sub_count(ix,iy,iz).y() -= temp * both_ends_fac / 2.0;
+ }
+
+ area_block.z() = ac_lblock(iz,ix) * yoverlap[iy];
+ area_block.x() = bc_lblock(iz,ix) * yoverlap[iy];
+
+ // Do not want blockage and drag across the end
+ // except at perimeter.
+ if (aboverlap(iz,ix) < pars.blockedFacePar)
+ {
+ if (obs.typeId == PDRobstacle::CYLINDER)
+ {
+ drag_r(ix,iy,iz).z() += c_drag(iz,ix).xx() * yoverlap[iy] / zgrid.width(iz) / xgrid.width(ix);
+ drag_r(ix,iy,iz).x() += c_drag(iz,ix).yy() * yoverlap[iy] / zgrid.width(iz) / xgrid.width(ix);
+
+ add_blockage_c(area_block_r(ix,iy,iz).z(), dirn_block(ix,iy,iz).z(), area_block.z(), 1.0);
+ add_blockage_c(area_block_r(ix,iy,iz).x(), dirn_block(ix,iy,iz).x(), area_block.x(), 1.0);
+ }
+ else
+ {
+ drag_s(ix,iy,iz).zz() += c_drag(iz,ix).xx() * yoverlap[iy] / zgrid.width(iz) / xgrid.width(ix);
+ drag_s(ix,iy,iz).xx() += c_drag(iz,ix).yy() * yoverlap[iy] / zgrid.width(iz) / xgrid.width(ix);
+
+ add_blockage_c(area_block_s(ix,iy,iz).z(), dirn_block(ix,iy,iz).z(), area_block.z(), 1.0);
+ add_blockage_c(area_block_s(ix,iy,iz).x(), dirn_block(ix,iy,iz).x(), area_block.x(), 1.0);
+ }
+ }
+ betai_inv1(ix,iy,iz).z() += vol_block / (1.0 - area_block.z() + floatSMALL);
+ betai_inv1(ix,iy,iz).x() += vol_block / (1.0 - area_block.x() + floatSMALL);
+ betai_inv1(ix,iy,iz).y() += vol_block / (1.0 - aboverlap(iz,ix) + floatSMALL);
+
+ // The off-diagonal elements of drag are stored in "drag" for BOTH round and sharp obstacles
+ drag_s(ix,iy,iz).xz() += c_drag(iz,ix).xy() * yoverlap[iy] / zgrid.width(iz) / xgrid.width(ix);
+
+ add_blockage_f(face_block(ix,iy,iz).z(), a_lblock(iz,ix) * yoverlap[iy], hole_in_face(ix,iy,iz).z());
+ add_blockage_f(face_block(ix,iy,iz).x(), b_lblock(iz,ix) * yoverlap[iy], hole_in_face(ix,iy,iz).x());
+ }
+ if (cymin == cymax)
+ {
+ // Does not span cell. Block nearest face.
+ add_blockage_f(face_block(ix,cfymax,iz).y(), aboverlap(iz,ix), hole_in_face(ix,cfymax,iz).y());
+ }
+ else
+ {
+ // In at least two cells.
+ // Block first and last faces overlapped
+ add_blockage_f(face_block(ix,cymin+1,iz).y(), aboverlap(iz,ix), hole_in_face(ix,cymin+1,iz).y());
+ if (cymax > cymin+1)
+ {
+ add_blockage_f(face_block(ix,cymax,iz).y(), aboverlap(iz,ix), hole_in_face(ix,cymax,iz).y());
+ }
+ }
+
+ for (label iy = cymin+1; iy < cymax; ++iy)
+ {
+ // Internal only
+ along_block(ix,iy,iz).y() += aboverlap(iz,ix);
+ }
+
+ drag_s(ix,cymin,iz).yy() += aboverlap(iz,ix) * zgrid.width(iz) * xgrid.width(ix) / 2.0;
+ drag_s(ix,cymax,iz).yy() += aboverlap(iz,ix) * zgrid.width(iz) * xgrid.width(ix) / 2.0;
+
+ add_blockage_c(area_block_s(ix,cymin,iz).y(), dirn_block(ix,cymin,iz).y(), aboverlap(iz,ix), 0.5);
+ add_blockage_c(area_block_s(ix,cymax,iz).y(), dirn_block(ix,cymax,iz).y(), aboverlap(iz,ix), 0.5);
+ }
+ }
+
+ break;
+ }
+
+ default: // orientation
+ {
+ Info<< "Unexpected orientation " << int(obs.orient) << nl;
+ break;
+ }
+ }
+}
+
+
+void Foam::PDRarrays::addBlockage
+(
+ const PDRobstacle& obs,
+ DynamicList& patches,
+ const int volumeSign
+)
+{
+ // The volumeSign indicates whether this pass is for negative or positive
+ // obstacles. Both if 0.
+ if
+ (
+ equal(obs.vbkge, 0)
+ || (volumeSign < 0 && obs.vbkge >= 0)
+ || (volumeSign > 0 && obs.vbkge < 0)
+ )
+ {
+ return;
+ }
+
+ if (isNull(block()))
+ {
+ FatalErrorInFunction
+ << "No PDRblock set" << nl
+ << exit(FatalError);
+ }
+
+ const PDRblock& pdrBlock = block();
+ const PDRblock::location& xgrid = pdrBlock.grid().x();
+ const PDRblock::location& ygrid = pdrBlock.grid().y();
+ const PDRblock::location& zgrid = pdrBlock.grid().z();
+
+ scalarList& xoverlap = overlap_1d.x();
+ scalarList& yoverlap = overlap_1d.y();
+ scalarList& zoverlap = overlap_1d.z();
+
+
+ // 0 will be used later for faces found to be blocked.
+ // 2 is used for wall-function faces.
+ label patchNum = PDRpatchDef::LAST_PREDEFINED;
+
+ // Only the part of the panel that covers full cell faces will be used
+ // so later should keep the panel in the list plus a hole (-ve obstacle) for
+ // part that will become blowoff b.c.
+
+ int indir = 0;
+
+ // Panel or patch
+ const bool isPatch =
+ (
+ (obs.typeId == PDRobstacle::LOUVRE_BLOWOFF && obs.blowoff_type > 0)
+ || (obs.typeId == PDRobstacle::RECT_PATCH)
+ );
+
+ if (isPatch)
+ {
+ const auto& identifier = obs.identifier;
+
+ const auto spc = identifier.find_first_of(" \t\n\v\f\r");
+
+ const word patchName = word::validate(identifier.substr(0, spc));
+
+ patchNum = ListOps::find
+ (
+ patches,
+ [=](const PDRpatchDef& p){ return patchName == p.patchName; },
+ 1 // skip 0 (blocked face)
+ );
+
+ if (patchNum < 1)
+ {
+ // The patch name was not already in the list
+ patchNum = patches.size();
+
+ patches.append(PDRpatchDef(patchName));
+ }
+
+
+ PDRpatchDef& p = patches[patchNum];
+
+ if (obs.typeId == PDRobstacle::RECT_PATCH)
+ {
+ indir = obs.inlet_dirn;
+ p.patchType = 0;
+ }
+ else
+ {
+ p.patchType = obs.blowoff_type;
+ p.blowoffPress = obs.blowoff_press;
+ p.blowoffTime = obs.blowoff_time;
+ if (obs.span.x() < 1e-5)
+ {
+ indir = 1;
+ }
+ else if (obs.span.y() < 1e-5)
+ {
+ indir = 2;
+ }
+ else if (obs.span.z() < 1e-5)
+ {
+ indir = 3;
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Blowoff panel should have zero thickness" << nl
+ << exit(FatalError);
+ }
+ }
+ }
+
+ int cxmin, cxmax, cfxmin, cfxmax;
+ one_d_overlap
+ (
+ obs.x(), obs.x() + obs.span.x(),
+ pdrBlock.grid().x(),
+ xoverlap, &cxmin, &cxmax, &cfxmin, &cfxmax
+ ); assert(cxmax >=0);
+
+ int cymin, cymax, cfymin, cfymax;
+ one_d_overlap
+ (
+ obs.y(), obs.y() + obs.span.y(),
+ pdrBlock.grid().y(),
+ yoverlap, &cymin, &cymax, &cfymin, &cfymax
+ ); assert(cymax >=0);
+
+ int czmin, czmax, cfzmin, cfzmax;
+ one_d_overlap
+ (
+ obs.z(), obs.z() + obs.span.z(),
+ pdrBlock.grid().z(),
+ zoverlap, &czmin, &czmax, &cfzmin, &cfzmax
+ ); assert(czmax >=0);
+
+ two_d_overlap(xoverlap, cxmin, cxmax, yoverlap, cymin, cymax, aboverlap);
+
+
+ const scalar vbkge = obs.vbkge;
+ const scalar xbkge = obs.xbkge;
+ const scalar ybkge = obs.ybkge;
+ const scalar zbkge = obs.zbkge;
+
+ // Compensate for double-counting of drag if two edges in same cell
+ const vector double_f
+ (
+ ((cxmin == cxmax) ? 0.5 : 1.0),
+ ((cymin == cymax) ? 0.5 : 1.0),
+ ((czmin == czmax) ? 0.5 : 1.0)
+ );
+
+ for (label ix = cxmin; ix <= cfxmax; ix++)
+ {
+ const scalar xov = xoverlap[ix];
+
+ scalar area, cell_area, temp;
+
+ for (label iy = cymin; iy <= cfymax; iy++)
+ {
+ const scalar yov = yoverlap[iy];
+
+ for (label iz = czmin; iz <= cfzmax; iz++)
+ {
+ const scalar zov = zoverlap[iz];
+
+ if
+ (
+ isPatch
+ &&
+ (
+ (indir == -1 && ix == cfxmin)
+ || (indir == 1 && ix == cfxmax)
+ || (indir == -2 && iy == cfymin)
+ || (indir == 2 && iy == cfymax)
+ || (indir == -3 && iz == cfzmin)
+ || (indir == 3 && iz == cfzmax)
+ )
+ )
+ {
+ /* Type RECT_PATCH (16) exists to set all faces it covers to be in a particular patch
+ usually an inlet or outlet.
+ ?? If the face not on a cell boundary, this will move it to the lower-cordinate
+ face of the relevant cell. It should be at the face of teh volume blocked by
+ the obstacle. But, if that is not at a cell boundary, the obstacle will be
+ putting blockage in front of teh vent, so we should be checking that it is
+ at a cell boundary. */
+
+ switch (indir) // Face orientation
+ {
+ // X
+ case -1:
+ case 1:
+ if (yov * zov > pars.blockedFacePar)
+ {
+ face_patch(ix,iy,iz).x() = patchNum;
+ }
+ break;
+
+ // Y
+ case -2:
+ case 2:
+ if (zov * xov > pars.blockedFacePar)
+ {
+ face_patch(ix,iy,iz).y() = patchNum;
+ }
+ break;
+
+ // Z
+ case -3:
+ case 3:
+ if (xov * yov > pars.blockedFacePar)
+ {
+ face_patch(ix,iy,iz).z() = patchNum;
+ }
+ break;
+ }
+ } // End of code for patch
+
+
+ const scalar vol_block = aboverlap(ix,iy) * zov * vbkge;
+ v_block(ix,iy,iz) += vol_block;
+
+ // These are the face blockages
+ if ((ix > cxmin && ix <= cxmax) || (cxmin == cxmax && ix == cfxmax))
+ {
+ temp = yov * zov * xbkge;
+
+ // Has -ve volumeSign only when processing user-defined
+ // -ve blocks
+ if (volumeSign < 0)
+ {
+ hole_in_face(ix,iy,iz).x() = true;
+ }
+ add_blockage_f(face_block(ix,iy,iz).x(), temp, hole_in_face(ix,iy,iz).x());
+ if (temp > pars.blockedFacePar && ! hole_in_face(ix,iy,iz).x() && !isPatch)
+ {
+ // Put faces of block in wall patch
+ face_patch(ix,iy,iz).x() = PDRpatchDef::WALL_PATCH;
+ }
+ }
+ if ((iy > cymin && iy <= cymax) || (cymin == cymax && iy == cfymax))
+ {
+ temp = zov * xov * ybkge;
+ if (volumeSign < 0)
+ {
+ hole_in_face(ix,iy,iz).y() = true;
+ }
+ add_blockage_f(face_block(ix,iy,iz).y(), temp, hole_in_face(ix,iy,iz).y());
+ if (temp > pars.blockedFacePar && ! hole_in_face(ix,iy,iz).y() && !isPatch)
+ {
+ face_patch(ix,iy,iz).y() = PDRpatchDef::WALL_PATCH;
+ }
+ }
+ if ((iz > czmin && iz <= czmax) || (czmin == czmax && iz == cfzmax))
+ {
+ temp = xov * yov * zbkge;
+ if (volumeSign < 0)
+ {
+ hole_in_face(ix,iy,iz).z() = true;
+ }
+ add_blockage_f(face_block(ix,iy,iz).z(), temp, hole_in_face(ix,iy,iz).z());
+ if (temp > pars.blockedFacePar && ! hole_in_face(ix,iy,iz).z() && !isPatch)
+ {
+ face_patch(ix,iy,iz).z() = PDRpatchDef::WALL_PATCH;
+ }
+ }
+
+ // These are the interior blockages
+ /* As for cylinders, blockage that extends longitudinally all the way through the cell
+ should not be in x_block etc., but it does go into new arrays along_block etc. */
+ area = yov * zov * xbkge; // Note this is fraction of cell c-s area
+ if (ix < cxmin || ix > cxmax)
+ {}
+ else if (ix > cxmin && ix < cxmax && xbkge >= 1.0)
+ {
+ along_block(ix,iy,iz).x() += area;
+ betai_inv1(ix,iy,iz).x() += vol_block / (1.0 - area + floatSMALL);
+ }
+ else if (ix == cxmin || ix == cxmax)
+ {
+ // If front and back of the obstacle are not in the
+ // same cell, put half in each
+ const scalar double_f = (cxmin == cxmax ? 1.0 : 0.5);
+
+ add_blockage_c(area_block_s(ix,iy,iz).x(), dirn_block(ix,iy,iz).x(), area, double_f);
+ cell_area = (ygrid.width(iy) * zgrid.width(iz));
+ surf(ix,iy,iz) += double_f * area * cell_area;
+ betai_inv1(ix,iy,iz).x() += vol_block / (1.0 - area + floatSMALL);
+
+ // To get Lobs right for a grating, allow for the fact that it is series of bars by having additional count
+ if (obs.typeId == PDRobstacle::GRATING && obs.orient == vector::X)
+ {
+ // * cell_area to make dimensional, then / std::sqrt(cell_area) to get width
+ temp = area * std::sqrt(cell_area) / obs.slat_width - 1.0;
+ if (temp > 0.0) { grating_count(ix,iy,iz).x() += temp; }
+ }
+ }
+
+ /* As for cylinders, blockage that extends longitudinally all the way through the cell
+ should not be in x_block etc., but it does go into new arrays along_block etc. */
+ area = zov * xov * ybkge;
+ if (iy < cymin || iy > cymax)
+ {}
+ else if (iy > cymin && iy < cymax && ybkge >= 1.0)
+ {
+ along_block(ix,iy,iz).y() += area;
+ betai_inv1(ix,iy,iz).y() += vol_block / (1.0 - area + floatSMALL);
+ }
+ else if (iy == cymin || iy == cymax)
+ {
+ // If front and back of the obstacle are not in the
+ // same cell, put half in each
+ const scalar double_f = (cymin == cymax ? 1.0 : 0.5);
+
+ add_blockage_c(area_block_s(ix,iy,iz).y(), dirn_block(ix,iy,iz).y(), area, double_f);
+ cell_area = (zgrid.width(iz) * xgrid.width(ix));
+ surf(ix,iy,iz) += double_f * area * cell_area;
+ betai_inv1(ix,iy,iz).y() += vol_block / (1.0 - area + floatSMALL);
+
+ if (obs.typeId == PDRobstacle::GRATING && obs.orient == vector::Y)
+ {
+ // * cell_area to make dimensional, then / std::sqrt(cell_area) to get width
+ temp = area * std::sqrt(cell_area) / obs.slat_width - 1.0;
+ if (temp > 0.0) { grating_count(ix,iy,iz).y() += temp; }
+ }
+ }
+
+ area = xov * yov * zbkge;
+ if (iz < czmin || iz > czmax)
+ {}
+ else if (iz > czmin && iz < czmax && zbkge >= 1.0)
+ {
+ along_block(ix,iy,iz).z() += area;
+ betai_inv1(ix,iy,iz).z() += vol_block / (1.0 - area + floatSMALL);
+ }
+ else if (iz == czmin || iz == czmax)
+ {
+ // If front and back of the obstacle are not in the
+ // same cell, put half in each
+ const scalar double_f = (czmin == czmax ? 1.0 : 0.5);
+
+ add_blockage_c(area_block_s(ix,iy,iz).z(), dirn_block(ix,iy,iz).z(), area, double_f);
+ cell_area = (xgrid.width(ix) * ygrid.width(iy));
+ surf(ix,iy,iz) += double_f * area * cell_area;
+ betai_inv1(ix,iy,iz).z() += vol_block / (1.0 - area + floatSMALL);
+
+ if (obs.typeId == PDRobstacle::GRATING && obs.orient == vector::Z)
+ {
+ // * cell_area to make dimensional, then / std::sqrt(cell_area) to get width
+ temp = area * std::sqrt(cell_area) / obs.slat_width - 1.0;
+ if (temp > 0.0) { grating_count(ix,iy,iz).z() += temp; }
+ }
+ }
+ }
+ }
+ }
+
+ if (obs.typeId == PDRobstacle::RECT_PATCH)
+ {
+ // Was only needed to set face_patch values
+ return;
+ }
+
+
+ /* A narrow obstacle completely crossing the cell adds 1 to the count for the transverse directions
+ If all four edges are not in the relevant cell, we take 1/4 for each edge that is in.
+ If it does not totally cross the cell, then the count is reduced proportionately
+ If it is porous, then the count is reduced proportionately.
+ ?? Should it be? At least this is consistent with effect going smoothly to zero as porosity approaches 1 ??
+
+ Note that more than 1 can be added for one obstacle, if sides and ends are in the cell.
+
+ We do all the x-aligned edges first, both for y-blockage and z-blockage. */
+
+ /* Intersection obstacles can have more edges than the intersecting obstacles that
+ generated them, so not good to incorporate these in N and drag. */
+
+ for (label ix = cxmin; ix <= cxmax; ix++)
+ {
+ // Factor of 0.25 because 4 edges to be done
+ scalar olap25 = 0.25 * xoverlap[ix];
+
+ const scalar temp =
+ ((pars.noIntersectN && vbkge < 0.334) ? 0 : (olap25 * vbkge));
+
+ obs_count(ix,cymin,czmin) += temp;
+ obs_count(ix,cymax,czmin) += temp;
+ obs_count(ix,cymin,czmax) += temp;
+ obs_count(ix,cymax,czmax) += temp;
+
+ sub_count(ix,cymin,czmin).x() += temp;
+ sub_count(ix,cymax,czmin).x() += temp;
+ sub_count(ix,cymin,czmax).x() += temp;
+ sub_count(ix,cymax,czmax).x() += temp;
+
+ // The 0.25 becomes 0.5 to allow for front/back faces in drag direction
+ olap25 *= 2.0;
+
+ drag_s(ix,cymin,czmin).yy() += zoverlap[czmin] * double_f.z() * olap25 * ybkge / ygrid.width(cymin);
+ drag_s(ix,cymax,czmin).yy() += zoverlap[czmin] * double_f.z() * olap25 * ybkge / ygrid.width(cymax);
+ drag_s(ix,cymin,czmax).yy() += zoverlap[czmax] * double_f.z() * olap25 * ybkge / ygrid.width(cymin);
+ drag_s(ix,cymax,czmax).yy() += zoverlap[czmax] * double_f.z() * olap25 * ybkge / ygrid.width(cymax);
+
+ drag_s(ix,cymin,czmin).zz() += yoverlap[cymin] * double_f.y() * olap25 * zbkge / zgrid.width(czmin);
+ drag_s(ix,cymax,czmin).zz() += yoverlap[cymax] * double_f.y() * olap25 * zbkge / zgrid.width(czmin);
+ drag_s(ix,cymin,czmax).zz() += yoverlap[cymin] * double_f.y() * olap25 * zbkge / zgrid.width(czmax);
+ drag_s(ix,cymax,czmax).zz() += yoverlap[cymax] * double_f.y() * olap25 * zbkge / zgrid.width(czmax);
+
+ // Porous obstacles do not only have drag at edges
+ if (xbkge < 1.0)
+ {
+ for (label iy = cymin+1; iy < cymax; iy++)
+ {
+ for (label iz = czmin+1; iz < czmax; iz++)
+ {
+ // Internal only
+ drag_s(ix,iy,iz).xx() = xbkge / xgrid.width(ix);
+ }
+ }
+ }
+ }
+
+ for (label iy = cymin; iy <= cymax; iy++)
+ {
+ scalar olap25 = 0.25 * yoverlap[iy];
+ const scalar temp =
+ ((pars.noIntersectN && vbkge < 0.334) ? 0 : (olap25 * vbkge));
+
+ obs_count(cxmin,iy,czmin) += temp;
+ obs_count(cxmax,iy,czmin) += temp;
+ obs_count(cxmin,iy,czmax) += temp;
+ obs_count(cxmax,iy,czmax) += temp;
+
+ sub_count(cxmin,iy,czmin).y() += temp;
+ sub_count(cxmax,iy,czmin).y() += temp;
+ sub_count(cxmin,iy,czmax).y() += temp;
+ sub_count(cxmax,iy,czmax).y() += temp;
+
+ olap25 *= 2.0;
+
+ if (iy > cymin && iy < cymax) // Avoid re-doing corners already done above
+ {
+ drag_s(cxmin,iy,czmin).zz() += xoverlap[cxmin] * double_f.x() * olap25 * zbkge / zgrid.width(czmin);
+ drag_s(cxmax,iy,czmin).zz() += xoverlap[cxmin] * double_f.x() * olap25 * zbkge / zgrid.width(czmin);
+ drag_s(cxmin,iy,czmax).zz() += xoverlap[cxmax] * double_f.x() * olap25 * zbkge / zgrid.width(czmax);
+ drag_s(cxmax,iy,czmax).zz() += xoverlap[cxmax] * double_f.x() * olap25 * zbkge / zgrid.width(czmax);
+ }
+ drag_s(cxmin,iy,czmin).xx() += zoverlap[czmin] * double_f.z() * olap25 * xbkge / xgrid.width(cxmin);
+ drag_s(cxmax,iy,czmin).xx() += zoverlap[czmax] * double_f.z() * olap25 * xbkge / xgrid.width(cxmax);
+ drag_s(cxmin,iy,czmax).xx() += zoverlap[czmin] * double_f.z() * olap25 * xbkge / xgrid.width(cxmin);
+ drag_s(cxmax,iy,czmax).xx() += zoverlap[czmax] * double_f.z() * olap25 * xbkge / xgrid.width(cxmax);
+
+ // Porous obstacles do not only have drag at edges
+ if (ybkge < 1.0)
+ {
+ for (label iz = czmin+1; iz < czmax; iz++)
+ {
+ for (label ix = cxmin+1; ix < cxmax; ix++)
+ {
+ // Internal only
+ drag_s(ix,iy,iz).yy() = ybkge / ygrid.width(iy);
+ }
+ }
+ }
+ }
+
+ for (label iz = czmin; iz <= czmax; iz++)
+ {
+ scalar olap25 = 0.25 * zoverlap[iz];
+ const scalar temp =
+ ((pars.noIntersectN && vbkge < 0.334) ? 0 : (olap25 * vbkge));
+
+ obs_count(cxmin,cymin,iz) += temp;
+ obs_count(cxmin,cymax,iz) += temp;
+ obs_count(cxmax,cymin,iz) += temp;
+ obs_count(cxmax,cymax,iz) += temp;
+
+ sub_count(cxmin,cymin,iz).z() += temp;
+ sub_count(cxmin,cymax,iz).z() += temp;
+ sub_count(cxmax,cymin,iz).z() += temp;
+ sub_count(cxmax,cymax,iz).z() += temp;
+
+ olap25 *= 2.0;
+
+ if (iz > czmin && iz < czmax) // Avoid re-doing corners already done above
+ {
+ drag_s(cxmin,cymin,iz).xx() += yoverlap[cymin] * double_f.y() * olap25 * xbkge / xgrid.width(cxmin);
+ drag_s(cxmax,cymin,iz).xx() += yoverlap[cymin] * double_f.y() * olap25 * xbkge / xgrid.width(cxmax);
+ drag_s(cxmin,cymax,iz).xx() += yoverlap[cymax] * double_f.y() * olap25 * xbkge / xgrid.width(cxmin);
+ drag_s(cxmax,cymax,iz).xx() += yoverlap[cymax] * double_f.y() * olap25 * xbkge / xgrid.width(cxmax);
+
+ drag_s(cxmin,cymin,iz).yy() += xoverlap[cxmin] * double_f.x() * olap25 * ybkge / ygrid.width(cymin);
+ drag_s(cxmax,cymin,iz).yy() += xoverlap[cxmax] * double_f.x() * olap25 * ybkge / ygrid.width(cymin);
+ drag_s(cxmin,cymax,iz).yy() += xoverlap[cxmin] * double_f.x() * olap25 * ybkge / ygrid.width(cymax);
+ drag_s(cxmax,cymax,iz).yy() += xoverlap[cxmax] * double_f.x() * olap25 * ybkge / ygrid.width(cymax);
+ }
+
+ // Porous obstacles do not only have drag at edges
+ if (zbkge < 1.0)
+ {
+ for (label ix = cxmin+1; ix < cxmax; ix++)
+ {
+ for (label iy = cymin+1; iy < cymax; iy++)
+ {
+ // Internal only
+ drag_s(ix,iy,iz).zz() = zbkge / zgrid.width(iz);
+ }
+ }
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRarraysCalc.C b/applications/utilities/preProcessing/PDRsetFields/PDRarraysCalc.C
new file mode 100644
index 0000000000..e74ff71d7d
--- /dev/null
+++ b/applications/utilities/preProcessing/PDRsetFields/PDRarraysCalc.C
@@ -0,0 +1,1984 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016 Shell Research Ltd.
+ Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+\*---------------------------------------------------------------------------*/
+
+#include "PDRarrays.H"
+#include "PDRblock.H"
+#include "PDRpatchDef.H"
+#include "PDRmeshArrays.H"
+#include "PDRparams.H"
+
+#include "PDRsetFields.H"
+
+#include "bitSet.H"
+#include "DynamicList.H"
+#include "dimensionSet.H"
+#include "symmTensor.H"
+#include "SquareMatrix.H"
+#include "IjkField.H"
+#include "MinMax.H"
+#include "volFields.H"
+#include "OFstream.H"
+#include "OSspecific.H"
+
+#ifndef FULLDEBUG
+#define NDEBUG
+#endif
+#include
+
+using namespace Foam;
+
+HashTable fieldNotes
+({
+ { "Lobs", "" },
+ { "Aw", "surface area per unit volume" },
+ { "CR", "" },
+ { "CT", "" },
+ { "N", "" },
+ { "ns", "" },
+ { "Nv", "" },
+ { "nsv", "" },
+ { "Bv", "area blockage" },
+ { "B", "area blockage" },
+ { "betai", "" },
+ { "Blong", "longitudinal blockage" },
+ { "Ep", "1/Lobs" },
+});
+
+
+// calc_fields
+
+
+// Local Functions
+/*
+// calc_drag_etc
+make_header
+tail_field
+write_scalarField
+write_uniformField
+write_symmTensorField
+write_pU_fields
+write_blocked_face_list
+write_blockedCellsSet
+*/
+
+// Somewhat similar to what the C-fprintf would have had
+static constexpr unsigned outputPrecision = 8;
+
+void calc_drag_etc
+(
+ double brs, double brr, bool blocked,
+ double surr_br, double surr_dr,
+ scalar* drags_p, scalar* dragr_p,
+ double count,
+ scalar* cbdi_p,
+ double cell_vol
+);
+
+
+void write_scalarField
+(
+ const word& fieldName, const IjkField& fld,
+ const scalar& deflt, const scalarMinMax& limits, const char *wall_bc,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const dimensionSet& dims, const fileName& casepath
+);
+
+void write_uniformField
+(
+ const word& fieldName, const scalar& deflt, const char *wall_bc,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const dimensionSet& dims, const fileName& casepath
+);
+
+void write_pU_fields
+(
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const fileName& casepath
+);
+
+void write_symmTensorField
+(
+ const word& fieldName, const IjkField& fld,
+ const symmTensor& deflt, const char *wall_bc,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const dimensionSet& dims, const fileName& casepath
+);
+
+void write_symmTensorFieldV
+(
+ const word& fieldName, const IjkField& fld,
+ const symmTensor& deflt, const char *wall_bc,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const dimensionSet& dims, const fileName& casepath
+);
+
+void write_blocked_face_list
+(
+ const IjkField& face_block,
+ const IjkField& face_patch,
+ const IjkField& obs_count,
+ IjkField& sub_count,
+ IjkField>& n_blocked_faces,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ double limit_par, const fileName& casepath
+);
+
+void write_blockedCellsSet
+(
+ const IjkField& fld,
+ const PDRmeshArrays& meshIndexing, double limit_par,
+ const IjkField>& n_blocked_faces,
+ const fileName& casepath,
+ const word& listName
+);
+
+
+// The average values of surrounding an array position
+static inline scalar averageSurrounding
+(
+ const SquareMatrix& mat,
+ const label i,
+ const label j
+)
+{
+ return
+ (
+ mat(i,j) + mat(i,j+1) + mat(i,j+2)
+ + mat(i+1,j) /* centre */ + mat(i+1,j+2)
+ + mat(i+2,j) + mat(i+2,j+1) + mat(i+2,j+2)
+ ) / 8.0; // Average
+}
+
+
+// Helper
+template
+static inline Ostream& putUniform(Ostream& os, const word& key, const Type& val)
+{
+ os.writeKeyword(key) << "uniform " << val << token::END_STATEMENT << nl;
+ return os;
+}
+
+
+static void make_header
+(
+ Ostream& os,
+ const fileName& location,
+ const word& clsName,
+ const word& object
+)
+{
+ string note = fieldNotes(object);
+
+ IOobject::writeBanner(os);
+
+ os << "FoamFile\n{\n"
+ << " version 2.0;\n"
+ << " format ascii;\n"
+ << " class " << clsName << ";\n";
+
+ if (!note.empty())
+ {
+ os << " note " << note << ";\n";
+ }
+
+ if (!location.empty())
+ {
+ os << " location " << location << ";\n";
+ }
+
+ os << " object " << object << ";\n"
+ << "}\n";
+
+ IOobject::writeDivider(os) << nl;
+}
+
+
+void Foam::PDRarrays::calculateAndWrite
+(
+ PDRarrays& arr,
+ const PDRmeshArrays& meshIndexing,
+ const fileName& casepath,
+ const UList& patches
+)
+{
+ if (isNull(arr.block()))
+ {
+ FatalErrorInFunction
+ << "No PDRblock set" << nl
+ << exit(FatalError);
+ }
+
+ const PDRblock& pdrBlock = arr.block();
+
+ const labelVector& cellDims = meshIndexing.cellDims;
+ const labelVector& faceDims = meshIndexing.faceDims;
+
+ const int xdim = faceDims.x();
+ const int ydim = faceDims.y();
+ const int zdim = faceDims.z();
+ const scalar maxCT = pars.maxCR * pars.cb_r;
+
+
+ // Later used to store the total effective blockage ratio per cell/direction
+ IjkField& drag_s = arr.drag_s;
+
+ IjkField& drag_r = arr.drag_r;
+
+ const IjkField& area_block_s = arr.area_block_s;
+ const IjkField& area_block_r = arr.area_block_r;
+ const IjkField>& dirn_block = arr.dirn_block;
+
+ const IjkField& betai_inv1 = arr.betai_inv1;
+
+ IjkField& obs_count = arr.obs_count;
+ IjkField& sub_count = arr.sub_count; // ns. Later used to hold longitudinal blockage
+ const IjkField& grating_count = arr.grating_count;
+
+ IjkField& v_block = arr.v_block;
+ IjkField& surf = arr.surf;
+
+ // Lobs. Later used for initial Ep
+ IjkField& obs_size = arr.obs_size;
+
+ Info<< "Calculating fields" << nl;
+
+ // Local scratch arrays
+
+ // The turbulance generation field CT.
+ // Later used to to hold the beta_i in tensor form
+ IjkField cbdi(cellDims, Zero);
+
+
+ // For 2D addressing it is convenient to just use the max dimension
+ // and avoid resizing when handling each direction.
+
+ // Dimension of the cells and a layer of surrounding halo cells
+ const labelVector surrDims = (faceDims + labelVector::uniform(2));
+
+ // Max addressing dimensions
+ const label maxDim = cmptMax(surrDims);
+
+ // Blockage-ratio correction to the drag
+ //
+ // neiBlock:
+ // 2-D for averaging the blockage ratio of neighbouring cells.
+ // It extends one cell outside the domain in each direction,
+ // so the indices are offset by 1.
+ // neiDrag:
+ // 2-D array for averaging the drag ratio of neighbouring cells
+
+ SquareMatrix neiBlock(maxDim, Zero);
+ SquareMatrix neiDrag(maxDim, Zero);
+
+ // X blockage, drag
+
+ for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
+ {
+ for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
+ {
+ for (label iz = 0; iz <= zdim; ++iz)
+ {
+ const label izz =
+ (iz == 0 ? 0 : iz == zdim ? zdim - 2 : iz - 1);
+
+ neiBlock(iy+1, iz) =
+ (
+ area_block_s(ix,iy,izz).x()
+ + area_block_r(ix,iy,izz).x()
+ );
+
+ neiDrag(iy+1, iz) =
+ (
+ drag_s(ix,iy,izz).xx() * pars.cd_s
+ + drag_r(ix,iy,izz).x() * pars.cd_r
+ );
+ }
+ }
+ for (label iz = 0; iz < surrDims.z(); ++iz)
+ {
+ if (pars.yCyclic)
+ {
+ // Cyclic in y
+ neiBlock(0, iz) = neiBlock(cellDims.y(), iz);
+ neiDrag(0, iz) = neiDrag(cellDims.y(), iz);
+ neiBlock(ydim, iz) = neiBlock(1, iz);
+ neiDrag(ydim, iz) = neiDrag(1, iz);
+ }
+ else
+ {
+ neiBlock(0, iz) = neiBlock(1, iz);
+ neiDrag(0, iz) = neiDrag(1, iz);
+ neiBlock(ydim, iz) = neiBlock(cellDims.y(), iz);
+ neiDrag(ydim, iz) = neiDrag(cellDims.y(), iz);
+ }
+ }
+
+ for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
+ {
+ for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
+ {
+ const scalar cell_vol = pdrBlock.V(ix,iy,iz);
+
+ const scalar surr_br = averageSurrounding(neiBlock, iy, iz);
+ const scalar surr_dr = averageSurrounding(neiDrag, iy, iz);
+
+ calc_drag_etc
+ (
+ area_block_s(ix,iy,iz).x(),
+ area_block_r(ix,iy,iz).x(),
+ dirn_block(ix,iy,iz).x(),
+ surr_br, surr_dr,
+ &(drag_s(ix,iy,iz).xx()),
+ &(drag_r(ix,iy,iz).x()),
+ obs_count(ix,iy,iz),
+ &(cbdi(ix,iy,iz).x()),
+ cell_vol
+ );
+ }
+ }
+ }
+
+
+ // Y blockage, drag
+
+ neiBlock = Zero;
+ neiDrag = Zero;
+
+ for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
+ {
+ for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
+ {
+ for (label ix = 0; ix <= xdim; ++ix)
+ {
+ const label ixx =
+ (ix == 0 ? 0 : ix == xdim ? xdim - 2 : ix - 1);
+
+ neiBlock(iz+1, ix) =
+ (
+ area_block_s(ixx,iy,iz).y()
+ + area_block_r(ixx,iy,iz).y()
+ );
+ neiDrag(iz+1, ix) =
+ (
+ drag_s(ixx,iy,iz).yy() * pars.cd_s
+ + drag_r(ixx,iy,iz).y() * pars.cd_r
+ );
+ }
+ }
+ for (label ix = 0; ix < surrDims.x(); ++ix)
+ {
+ neiBlock(0, ix) = neiBlock(1, ix);
+ neiDrag(0, ix) = neiDrag(1, ix);
+ neiBlock(zdim, ix) = neiBlock(cellDims.z(), ix);
+ neiDrag(zdim, ix) = neiDrag(cellDims.z(), ix);
+ }
+
+ for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
+ {
+ for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
+ {
+ const scalar cell_vol = pdrBlock.V(ix,iy,iz);
+
+ const scalar surr_br = averageSurrounding(neiBlock, iz, ix);
+ const scalar surr_dr = averageSurrounding(neiDrag, iz, ix);
+
+ calc_drag_etc
+ (
+ area_block_s(ix,iy,iz).y(),
+ area_block_r(ix,iy,iz).y(),
+ dirn_block(ix,iy,iz).y(),
+ surr_br, surr_dr,
+ &(drag_s(ix,iy,iz).yy()),
+ &(drag_r(ix,iy,iz).y()),
+ obs_count(ix,iy,iz),
+ &(cbdi(ix,iy,iz).y()),
+ cell_vol
+ );
+ }
+ }
+ }
+
+
+ // Z blockage, drag
+
+ neiBlock = Zero;
+ neiDrag = Zero;
+
+ for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
+ {
+ for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
+ {
+ for (label iy = 0; iy <= ydim; ++iy)
+ {
+ label iyy;
+
+ if (pars.yCyclic)
+ {
+ iyy = (iy == 0 ? ydim - 2 : iy == ydim ? 0 : iy - 1);
+ }
+ else
+ {
+ iyy = (iy == 0 ? 0 : iy == ydim ? ydim - 2 : iy - 1);
+ }
+
+ neiBlock(ix+1, iy) =
+ (
+ area_block_s(ix,iyy,iz).z()
+ + area_block_r(ix,iyy,iz).z()
+ );
+ neiDrag(ix+1, iy) =
+ (
+ drag_s(ix,iyy,iz).zz() * pars.cd_s
+ + drag_r(ix,iyy,iz).z() * pars.cd_r
+ );
+ }
+ }
+ for (label iy = 0; iy < surrDims.y(); ++iy)
+ {
+ neiBlock(0, iy) = neiBlock(1, iy);
+ neiDrag(0, iy) = neiDrag(1, iy);
+ neiBlock(xdim, iy) = neiBlock(cellDims.x(), iy);
+ neiDrag(xdim, iy) = neiDrag(cellDims.x(), iy);
+ }
+
+ for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
+ {
+ for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
+ {
+ const scalar cell_vol = pdrBlock.V(ix,iy,iz);
+
+ const scalar surr_br = averageSurrounding(neiBlock, ix, iy);
+ const scalar surr_dr = averageSurrounding(neiDrag, ix, iy);
+
+ calc_drag_etc
+ (
+ area_block_s(ix,iy,iz).z(),
+ area_block_r(ix,iy,iz).z(),
+ dirn_block(ix,iy,iz).z(),
+ surr_br, surr_dr,
+ &(drag_s(ix,iy,iz).zz()),
+ &(drag_r(ix,iy,iz).z()),
+ obs_count(ix,iy,iz),
+ &(cbdi(ix,iy,iz).z()),
+ cell_vol
+ );
+ }
+ }
+ }
+
+ neiBlock.clear();
+ neiDrag.clear();
+
+
+ // Calculate other parameters
+
+ for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
+ {
+ for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
+ {
+ for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
+ {
+ const scalar dx = pdrBlock.dx(ix);
+ const scalar dy = pdrBlock.dy(iy);
+ const scalar dz = pdrBlock.dz(iz);
+ const scalar cell_vol = pdrBlock.V(ix, iy, iz);
+ const scalar cell_size = pdrBlock.width(ix, iy, iz);
+
+ drag_s(ix,iy,iz).xy() *= pars.cd_s;
+ drag_s(ix,iy,iz).xz() *= pars.cd_s;
+ drag_s(ix,iy,iz).yz() *= pars.cd_s;
+
+ if (drag_s(ix,iy,iz).xx() > pars.maxCR) { drag_s(ix,iy,iz).xx() = pars.maxCR; } ;
+ if (drag_s(ix,iy,iz).yy() > pars.maxCR) { drag_s(ix,iy,iz).yy() = pars.maxCR; } ;
+ if (drag_s(ix,iy,iz).zz() > pars.maxCR) { drag_s(ix,iy,iz).zz() = pars.maxCR; } ;
+
+ if (cbdi(ix,iy,iz).x() > maxCT ) { cbdi(ix,iy,iz).x() = maxCT; } ;
+ if (cbdi(ix,iy,iz).y() > maxCT ) { cbdi(ix,iy,iz).y() = maxCT; } ;
+ if (cbdi(ix,iy,iz).z() > maxCT ) { cbdi(ix,iy,iz).z() = maxCT; } ;
+
+ surf(ix,iy,iz) /= cell_vol;
+
+ /* Calculate length scale of obstacles in each cell
+ Result is stored in surf. */
+
+ {
+ const scalar vb = v_block(ix,iy,iz);
+
+ if
+ (
+ (
+ ((area_block_s(ix,iy,iz).x() + area_block_r(ix,iy,iz).x()) < MIN_AB_FOR_SIZE)
+ && ((area_block_s(ix,iy,iz).y() + area_block_r(ix,iy,iz).y()) < MIN_AB_FOR_SIZE)
+ && ((area_block_s(ix,iy,iz).z() + area_block_r(ix,iy,iz).z()) < MIN_AB_FOR_SIZE)
+ )
+ || ( vb > MAX_VB_FOR_SIZE )
+ || ((obs_count(ix,iy,iz) + cmptSum(grating_count(ix,iy,iz))) < MIN_COUNT_FOR_SIZE)
+ || ( surf(ix,iy,iz) <= 0.0 )
+ )
+ {
+ obs_size(ix,iy,iz) = cell_size * pars.empty_lobs_fac;
+ }
+ else
+ {
+ /* A small sliver of a large cylinder ina cell can give large surface area
+ but low volume, hence snall "size". Therefore the vol/area formulation
+ is only fully implemented when count is at least COUNT_FOR_SIZE.*/
+ double nn, lobs, lobsMax;
+ nn = obs_count(ix,iy,iz) - sub_count(ix,iy,iz).x() + grating_count(ix,iy,iz).x();
+ if ( nn < 1.0 ) { nn = 1.0; }
+ lobsMax = (area_block_s(ix,iy,iz).x() + area_block_r(ix,iy,iz).x()) / nn * std::sqrt( dy * dz );
+ nn = obs_count(ix,iy,iz) - sub_count(ix,iy,iz).y() + grating_count(ix,iy,iz).y();
+ if ( nn < 1.0 ) { nn = 1.0; }
+ lobs = (area_block_s(ix,iy,iz).y() + area_block_r(ix,iy,iz).y()) / nn * std::sqrt( dz * dx );
+ if ( lobs > lobsMax )
+ {
+ lobsMax = lobs;
+ }
+
+ nn = obs_count(ix,iy,iz) - sub_count(ix,iy,iz).z() + grating_count(ix,iy,iz).z();
+ if ( nn < 1.0 ) { nn = 1.0; }
+ lobs = (area_block_s(ix,iy,iz).z() + area_block_r(ix,iy,iz).z()) / nn * std::sqrt( dx * dy );
+ if ( lobs > lobsMax )
+ {
+ lobsMax = lobs;
+ }
+
+ obs_size(ix,iy,iz) = lobsMax;
+ }
+ }
+
+ /* The formulation correctly deals with triple intersections. For quadruple intersections
+ and worse, there are very many second level overlaps and the resulting volume can be large
+ positive. However, many or all of these may be eliminated because of the minimum volume of
+ overlap blocks. Then the result can be negative volume - constrain accordingly
+ */
+
+ if (v_block(ix,iy,iz) < 0)
+ {
+ v_block(ix,iy,iz) = 0;
+ }
+ else if (v_block(ix,iy,iz) > 1)
+ {
+ v_block(ix,iy,iz) = 1;
+ }
+
+ /* We can get -ve sub_count (ns) if two pipes/bars intersect and the dominat direction
+ of the (-ve) intersection block is not the same as either of the intersecting obstacles.
+ Also, if we have two hirizontal abrs intersecting, the overlap block can have vertical
+ edges in a cell where the original bars do not. This can give -ve N and ns.
+ Negative N is removed by write_scalar. */
+
+ for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+ {
+ if (sub_count(ix,iy,iz)[cmpt] < 0)
+ {
+ sub_count(ix,iy,iz)[cmpt] = 0;
+ }
+ }
+
+ v_block(ix,iy,iz) = 1.0 - v_block(ix,iy,iz); // Now porosity
+ }
+ }
+ }
+
+
+//*** Now we start writing the fields *********//
+
+ /* v_block is now porosity
+ The maximum value does not override the default value placed in the external cells,
+ so pars.cong_max_betav can be set just below 1 to mark the congested-region cells
+ for use by the adaptive mesh refinement. */
+
+ IjkField> n_blocked_faces
+ (
+ faceDims,
+ Vector::uniform(0)
+ );
+
+ write_blocked_face_list
+ (
+ arr.face_block, arr.face_patch,
+ obs_count, sub_count, n_blocked_faces,
+ meshIndexing, patches,
+ pars.blockedFacePar, casepath
+ );
+ write_blockedCellsSet
+ (
+ arr.v_block,
+ meshIndexing, pars.blockedCellPoros, n_blocked_faces,
+ casepath, "blockedCellsSet"
+ );
+
+ write_scalarField
+ (
+ "betav", arr.v_block, 1, {0, pars.cong_max_betav}, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+
+ for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
+ {
+ for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
+ {
+ for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
+ {
+ const scalar cell_vol = pdrBlock.V(ix, iy, iz);
+
+ /* After the correction to set the number of obstacles normal to a blocked face
+ to be zero, we can have N and all the components of ns the same. Then there
+ are no obstacles in the cell as the number in each direction is n minus ns component),
+ but N is not zero. This can cause problems. We reduce all four numbers by the same amount,
+ which is OK as only the difference is used except when N is checked to se if there are
+ any obstacles in then cell. */
+
+ scalar nmin = cmptMin(sub_count(ix,iy,iz));
+
+ sub_count(ix,iy,iz).x() -= nmin;
+ sub_count(ix,iy,iz).y() -= nmin;
+ sub_count(ix,iy,iz).z() -= nmin;
+
+ obs_count(ix,iy,iz) -= nmin;
+
+ assert(obs_count(ix,iy,iz) > -1);
+ if ( pars.new_fields )
+ {
+ /* New fields Nv and nsv are intensive quantities that stay unchanged as a cell is subdivided
+ We do not divide by cell volume because we assume that typical obstacle
+ is a cylinder passing through the cell */
+ const scalar cell_23 = ::pow(cell_vol, 2.0/3.0);
+ obs_count(ix,iy,iz) /= cell_23;
+ sub_count(ix,iy,iz) /= cell_23;
+ }
+ }
+ }
+ }
+
+
+ {
+ Info<< "Writing field files" << nl;
+
+ // obs_size is now the integral scale of the generated turbulence
+ write_scalarField
+ (
+ "Lobs", arr.obs_size, DEFAULT_LOBS, {0, 10}, "zeroGradient",
+ meshIndexing, patches,
+ dimLength, casepath
+ );
+ // surf is now surface area per unit volume
+ write_scalarField
+ (
+ "Aw", arr.surf, 0, {0, 1000}, "zeroGradient",
+ meshIndexing, patches,
+ inv(dimLength), casepath
+ );
+ write_symmTensorField
+ (
+ "CR", arr.drag_s, Zero, "zeroGradient",
+ meshIndexing, patches, inv(dimLength), casepath
+ );
+ write_symmTensorFieldV
+ (
+ "CT", cbdi, Zero, "zeroGradient",
+ meshIndexing, patches,
+ inv(dimLength), casepath
+ );
+ if ( pars.new_fields )
+ {
+ // These have been divided by cell volume ^ (2/3)
+ write_scalarField
+ (
+ "Nv", arr.obs_count, 0, {0, 1000}, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_symmTensorFieldV
+ (
+ "nsv", arr.sub_count, Zero, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ }
+ else
+ {
+ write_scalarField
+ (
+ "N", arr.obs_count, 0, {0, 1000}, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_symmTensorFieldV
+ (
+ "ns", arr.sub_count, Zero, "zeroGradient",
+ meshIndexing, patches, dimless, casepath
+ );
+ }
+
+ // Compute some further variables; store in already used arrays
+ // Re-use the drag array
+ drag_s = Zero;
+
+ for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
+ {
+ for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
+ {
+ for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
+ {
+ // Effective blockage ratio per cell/direction
+ vector eff_block =
+ (
+ area_block_s(ix,iy,iz) * pars.cd_s/pars.cd_r
+ + area_block_r(ix,iy,iz)
+ );
+
+ // Convert from B to Bv
+ if (pars.new_fields)
+ {
+ eff_block /= pdrBlock.width(ix, iy, iz);
+ }
+
+ // Effective blockage is zero when faces are blocked
+ for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+ {
+ if (dirn_block(ix,iy,iz)[cmpt] || eff_block[cmpt] < 0)
+ {
+ eff_block[cmpt] = 0;
+ }
+ }
+
+ // Use the drag array to store the total effective blockage ratio per cell/direction
+ // - off-diagonal already zeroed
+ drag_s(ix,iy,iz).xx() = eff_block.x();
+ drag_s(ix,iy,iz).yy() = eff_block.y();
+ drag_s(ix,iy,iz).zz() = eff_block.z();
+
+ cbdi(ix,iy,iz).x() = 1.0 / (betai_inv1(ix,iy,iz).x() + 1.0);
+ cbdi(ix,iy,iz).y() = 1.0 / (betai_inv1(ix,iy,iz).y() + 1.0);
+ cbdi(ix,iy,iz).z() = 1.0 / (betai_inv1(ix,iy,iz).z() + 1.0);
+
+ if (cbdi(ix,iy,iz).z() < 0 || cbdi(ix,iy,iz).z() > 1.0)
+ {
+ WarningInFunction
+ << "beta_i problem. z-betai_inv1=" << betai_inv1(ix,iy,iz).z()
+ << " beta_i=" << cbdi(ix,iy,iz).z()
+ << nl;
+ }
+
+ //Use the obs_size array to store Ep
+ //We use Ep/(Xp-0.999) as length scale to avoid divide by zero,
+ // so this is OK for initial Xp=1.
+ obs_size(ix,iy,iz) = 0.001 / obs_size(ix,iy,iz);
+
+ // Use the count array to store the combustion flag ( --1 everywhere in rectangular cells).
+ obs_count(ix,iy,iz) = 1.0;
+ }
+ }
+ }
+
+ // drag array holds area blockage
+ if ( pars.new_fields )
+ {
+ write_symmTensorField
+ (
+ "Bv", arr.drag_s, Zero, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ }
+ else
+ {
+ write_symmTensorField
+ (
+ "B", arr.drag_s, Zero, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ }
+
+ // cbdi array holds beta_i
+ write_symmTensorFieldV
+ (
+ "betai", cbdi, symmTensor::I, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+
+ // The longitudinal blockage
+ write_symmTensorFieldV
+ (
+ "Blong", arr.along_block, Zero, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+
+ // obs_size array now contains 1/Lobs
+ write_scalarField
+ (
+ "Ep", arr.obs_size, DEFAULT_EP, {0, 10}, "zeroGradient",
+ meshIndexing, patches,
+ inv(dimLength), casepath
+ );
+ write_uniformField
+ (
+ "b", 1.0, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_uniformField
+ (
+ "k", DEFAULT_K, K_WALL_FN,
+ meshIndexing, patches,
+ sqr(dimVelocity),
+ casepath
+ );
+
+ write_uniformField
+ (
+ "epsilon", DEFAULT_EPS, EPS_WALL_FN,
+ meshIndexing, patches,
+ sqr(dimVelocity)/dimTime, casepath
+ );
+ write_uniformField
+ (
+ "ft", 0, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_uniformField
+ (
+ "Su", DEFAULT_SU, "zeroGradient",
+ meshIndexing, patches,
+ dimVelocity, casepath
+ );
+ write_uniformField
+ (
+ "T", DEFAULT_T, "zeroGradient",
+ meshIndexing, patches,
+ dimTemperature, casepath
+ );
+ write_uniformField
+ (
+ "Tu", DEFAULT_T, "zeroGradient",
+ meshIndexing, patches,
+ dimTemperature, casepath
+ );
+ write_uniformField
+ (
+ "Xi", 1, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_uniformField
+ (
+ "Xp", 1, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_uniformField
+ (
+ "GRxp", 0, "zeroGradient",
+ meshIndexing, patches,
+ inv(dimTime), casepath
+ );
+ write_uniformField
+ (
+ "GRep", 0, "zeroGradient",
+ meshIndexing, patches,
+ inv(dimLength*dimTime), casepath
+ );
+ write_uniformField
+ (
+ "RPers", 0, "zeroGradient",
+ meshIndexing, patches,
+ inv(dimTime), casepath
+ );
+ write_pU_fields(meshIndexing, patches, casepath);
+
+ write_uniformField
+ (
+ "alphat", 0, ALPHAT_WALL,
+ meshIndexing, patches,
+ dimMass/(dimLength*dimTime),
+ casepath
+ );
+ write_uniformField
+ (
+ "mut", 0, MUT_WALL_FN,
+ meshIndexing, patches,
+ dimDynamicViscosity, casepath
+ );
+ // combustFlag is 1 in rectangular region, 0 or 1 elsewhere
+ // (although user could set it to another value)
+ if (equal(pars.outerCombFac, 1))
+ {
+ write_uniformField
+ (
+ "combustFlag", 1, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ }
+ else
+ {
+ write_scalarField
+ (
+ "combustFlag", arr.obs_count, pars.outerCombFac, {0, 1}, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ }
+ if ( pars.deluge )
+ {
+ write_uniformField
+ (
+ "H2OPS", 0, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_uniformField
+ (
+ "AIR", 0, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_uniformField
+ (
+ "Ydefault", 0, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_uniformField
+ (
+ "eRatio", 1, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ write_uniformField
+ (
+ "sprayFlag", 1, "zeroGradient",
+ meshIndexing, patches,
+ dimless, casepath
+ );
+ }
+ }
+}
+
+
+void Foam::PDRarrays::calculateAndWrite
+(
+ const fileName& casepath,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches
+)
+{
+ calculateAndWrite(*this, meshIndexing, casepath, patches);
+}
+
+
+void calc_drag_etc
+(
+ double brs, double brr, bool blocked,
+ double surr_br, double surr_dr,
+ scalar* drags_p, scalar* dragr_p,
+ double count,
+ scalar* cbdi_p,
+ double cell_vol
+)
+{
+ // Total blockage ratio
+ scalar br = brr + brs;
+
+ // Idealise obstacle arrangement as sqrt(count) rows.
+ // Make br the blockage ratio for each row.
+ if (count > 1.0) { br /= std::sqrt(count); }
+
+ const scalar alpha =
+ (
+ br < 0.99
+ ? (1.0 - 0.5 * br) / (1.0 - br) / (1.0 - br)
+ : GREAT
+ );
+
+ // For the moment keep separate the two contributions to the blockage-corrected drag
+ /* An isolated long obstcale will have two of the surronding eight cells with the same blockage,
+ so surr_br would be br/4. In this case no correction. Rising to full correction when
+ all surrounding cells have the same blockage. */
+ const scalar expon =
+ (
+ br > 0.0
+ ? min(max((surr_br / br - 0.25) * 4.0 / 3.0, scalar(0)), scalar(1))
+ : 0.0
+ );
+
+ const scalar alpha_r = ::pow(alpha, 0.5 + 0.5 * expon);
+ const scalar alpha_s = ::pow(alpha, expon);
+
+ *dragr_p *= alpha_r;
+ *drags_p *= ::pow(alpha_s, 1.09);
+ *cbdi_p = ( pars.cb_r * pars.cd_r * *dragr_p + pars.cb_s * pars.cd_s * *drags_p );
+ if ( *cbdi_p < 0.0 ) { *cbdi_p = 0.0; }
+
+ // Finally sum the drag.
+ *drags_p = ( *drags_p * pars.cd_s + *dragr_p * pars.cd_r );
+ if ( *drags_p < 0.0 ) { *drags_p = 0.0; }
+ /* If well-blocked cells are surrounded by empty cells, the flow just goes round
+ and the drag parameters have little effect. So, for any cells much more empty
+ than the surrounding cells, we put some CR in there as well. */
+ if ( (surr_dr * 0.25) > *drags_p )
+ {
+ *drags_p = surr_dr * 0.25;
+ *cbdi_p = *drags_p * (pars.cb_r + pars.cb_s ) * 0.5;
+ // Don't know whether surr. stuff was round or sharp; use average of cb factors
+ }
+ if ( blocked ) { *cbdi_p = 0.0; *drags_p = 0.0; *dragr_p = 0.0; }
+}
+
+
+void Foam::PDRarrays::blockageSummary() const
+{
+ if (isNull(block()))
+ {
+ WarningInFunction
+ << nl
+ << "No blockage information - PDRblock is not set" << nl;
+ return;
+ }
+
+ const PDRblock& pdrBlock = block();
+
+ scalar totArea = 0;
+ scalar totCount = 0;
+ scalar totVolBlock = 0;
+
+ vector totBlock(Zero);
+ vector totDrag(Zero);
+
+ for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
+ {
+ for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
+ {
+ for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
+ {
+ const labelVector ijk(ix,iy,iz);
+
+ totVolBlock += v_block(ijk) * pdrBlock.V(ijk);
+ totArea += surf(ijk);
+
+ totCount += max(0, obs_count(ijk));
+
+ totDrag.x() += max(0, drag_s(ijk).xx());
+ totDrag.y() += max(0, drag_s(ijk).yy());
+ totDrag.z() += max(0, drag_s(ijk).zz());
+
+ for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+ {
+ totBlock[cmpt] += max(0, area_block_s(ijk)[cmpt]);
+ totBlock[cmpt] += max(0, area_block_r(ijk)[cmpt]);
+ }
+ }
+ }
+ }
+
+ Info<< nl
+ << "Volume blockage: " << totVolBlock << nl
+ << "Total drag: " << totDrag << nl
+ << "Total count: " << totCount << nl
+ << "Total area blockage: " << totBlock << nl
+ << "Total surface area: " << totArea << nl;
+}
+
+
+// ------------------------------------------------------------------------- //
+
+// Another temporary measure
+template
+static void tail_field
+(
+ Ostream& os,
+ const Type& deflt,
+ const char* wall_bc,
+ const UList& patches
+)
+{
+ // seaGround
+ {
+ os.beginBlock("seaGround");
+ os.writeKeyword("type") << wall_bc << token::END_STATEMENT << nl;
+ putUniform(os, "value", deflt);
+ os.endBlock();
+ }
+
+ forAll(patches, patchi)
+ {
+ const word& patchName = patches[patchi].patchName;
+
+ if (PDRpatchDef::BLOCKED_FACE == patchi)
+ {
+ // blockedFaces
+ os.beginBlock(patchName);
+
+ // No wall functions for blockedFaces patch unless selected
+ if (pars.blockedFacesWallFn)
+ {
+ os.writeKeyword("type") << wall_bc << token::END_STATEMENT << nl;
+ putUniform(os, "value", deflt);
+ }
+ else
+ {
+ os.writeEntry("type", "zeroGradient");
+ }
+
+ os.endBlock();
+ }
+ else if (patches[patchi].patchType == 0)
+ {
+ os.beginBlock(patchName);
+
+ os.writeKeyword("type") << wall_bc << token::END_STATEMENT << nl;
+ putUniform(os, "value", deflt);
+
+ os.endBlock();
+ }
+ else
+ {
+ os.beginBlock(word(patchName + "Wall"));
+ os.writeKeyword("type") << wall_bc << token::END_STATEMENT << nl;
+ putUniform(os, "value", deflt);
+ os.endBlock();
+
+ os.beginBlock(word(patchName + "Cyclic_half0"));
+ os.writeEntry("type", "cyclic");
+ os.endBlock();
+
+ os.beginBlock(word(patchName + "Cyclic_half1"));
+ os.writeEntry("type", "cyclic");
+ os.endBlock();
+ }
+ }
+
+ if (pars.yCyclic)
+ {
+ os.beginBlock("Cyclic_half0");
+ os.writeEntry("type", "cyclic");
+ os.endBlock();
+
+ os.beginBlock("Cyclic_half1");
+ os.writeEntry("type", "cyclic");
+ os.endBlock();
+ }
+ else
+ {
+ os.beginBlock("ySymmetry");
+ os.writeEntry("type", "symmetryPlane");
+ os.endBlock();
+ }
+
+ if ( pars.two_d )
+ {
+ os.beginBlock("z_boundaries");
+ os.writeEntry("type", "empty");
+ os.endBlock();
+ }
+ if ( pars.outer_orthog )
+ {
+ os.beginBlock("outer_inner");
+ os.writeEntry("type", "cyclicAMI");
+ os.writeEntry("neighbourPatch", "inner_outer");
+ os.endBlock();
+
+ os.beginBlock("inner_outer");
+ os.writeEntry("type", "cyclicAMI");
+ os.writeEntry("neighbourPatch", "outer_inner");
+ os.endBlock();
+ }
+}
+
+
+// ------------------------------------------------------------------------- //
+
+void write_scalarField
+(
+ const word& fieldName, const IjkField& fld,
+ const scalar& deflt, const scalarMinMax& limits, const char *wall_bc,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const dimensionSet& dims, const fileName& casepath
+)
+{
+ fileName path = (casepath / pars.timeName / fieldName);
+ OFstream os(path);
+ os.precision(outputPrecision);
+
+ make_header(os, "", volScalarField::typeName, fieldName);
+
+ os.writeEntry("dimensions", dims);
+
+ os << nl;
+ os.writeKeyword("internalField")
+ << "nonuniform List" << nl
+ << meshIndexing.nCells() << nl << token::BEGIN_LIST << nl;
+
+ for (label celli=0; celli < meshIndexing.nCells(); ++celli)
+ {
+ const labelVector& cellIdx = meshIndexing.cellIndex[celli];
+
+ if (cmptMin(cellIdx) < 0)
+ {
+ os << deflt << nl;
+ continue;
+ }
+
+ os << limits.clip(fld(cellIdx)) << nl;
+ }
+
+ os << token::END_LIST << token::END_STATEMENT << nl;
+
+ os << nl;
+ os.beginBlock("boundaryField");
+
+ // outer
+ {
+ os.beginBlock("outer");
+
+ os.writeEntry("type", "inletOutlet");
+ putUniform(os, "inletValue", deflt);
+ putUniform(os, "value", deflt);
+
+ os.endBlock();
+ }
+
+ tail_field(os, deflt, wall_bc, patches);
+
+ os.endBlock(); // boundaryField
+
+ IOobject::writeEndDivider(os);
+}
+
+
+// ------------------------------------------------------------------------- //
+
+void write_uniformField
+(
+ const word& fieldName, const scalar& deflt, const char *wall_bc,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const dimensionSet& dims, const fileName& casepath
+)
+{
+ OFstream os(casepath / pars.timeName / fieldName);
+ os.precision(outputPrecision);
+
+ make_header(os, "", volScalarField::typeName, fieldName);
+
+ os.writeEntry("dimensions", dims);
+
+ os << nl;
+ putUniform(os, "internalField", deflt);
+
+ os << nl;
+ os.beginBlock("boundaryField");
+
+ // outer
+ {
+ os.beginBlock("outer");
+ if (fieldName == "alphat" || fieldName == "mut")
+ {
+ // Different b.c. for alphat & mut
+ os.writeEntry("type", "calculated");
+ }
+ else
+ {
+ os.writeEntry("type", "inletOutlet");
+ putUniform(os, "inletValue", deflt);
+ }
+
+ putUniform(os, "value", deflt);
+ os.endBlock();
+ }
+
+ tail_field(os, deflt, wall_bc, patches);
+
+ os.endBlock(); // boundaryField
+
+ IOobject::writeEndDivider(os);
+}
+
+
+// ------------------------------------------------------------------------- //
+
+void write_pU_fields
+(
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const fileName& casepath
+)
+{
+ // Velocity field
+ {
+ OFstream os(casepath / pars.timeName / "Ubet");
+ os.precision(outputPrecision);
+
+ make_header(os, "", volVectorField::typeName, "Ubet");
+
+ os.writeEntry("dimensions", dimVelocity);
+
+ os << nl;
+ putUniform(os, "internalField", vector::zero);
+
+ os << nl;
+ os.beginBlock("boundaryField");
+
+ // "outer"
+ {
+ os.beginBlock("outer");
+ os.writeEntry("type", "inletOutlet");
+ putUniform(os, "inletValue", vector::zero);
+ os.endBlock();
+ }
+
+ // seaGround
+ {
+ os.beginBlock("seaGround");
+ os.writeEntry("type", "zeroGradient");
+ os.endBlock();
+ }
+
+ // Patch 0 is the blocked faces' and 1 is mergingFaces for ignition cell
+ for (label patchi = 0; patchi < 3; ++patchi)
+ {
+ os.beginBlock(patches[patchi].patchName);
+ os.writeKeyword("type") << pars.UPatchBc.c_str()
+ << token::END_STATEMENT << nl;
+ os.endBlock();
+ }
+
+ for (label patchi = 3; patchi < patches.size(); ++patchi)
+ {
+ const PDRpatchDef& p = patches[patchi];
+ const word& patchName = p.patchName;
+
+ if (p.patchType == 0)
+ {
+ os.beginBlock(patchName);
+
+ os.writeEntry("type", "timeVaryingMappedFixedValue");
+ os.writeEntry("fileName", "" / (patchName + ".dat"));
+ os.writeEntry("outOfBounds", "clamp");
+ putUniform(os, "value", vector::zero);
+ os.endBlock();
+ }
+ else
+ {
+ os.beginBlock(word(patchName + "Wall"));
+ os.writeEntry("type", "activePressureForceBaffleVelocity");
+
+ os.writeEntry("cyclicPatch", word(patchName + "Cyclic_half0"));
+ os.writeEntry("openFraction", 0); // closed
+ os.writeEntry("openingTime", p.blowoffTime);
+ os.writeEntry("maxOpenFractionDelta", 0.1);
+ os.writeEntry("forceBased", "false");
+ os.writeEntry("opening", "true");
+
+ putUniform(os, "value", vector::zero);
+ os.endBlock();
+
+ os.beginBlock(word(patchName + "Cyclic_half0"));
+ os.writeEntry("type", "cyclic");
+ putUniform(os, "value", vector::zero);
+ os.endBlock();
+
+ os.beginBlock(word(patchName + "Cyclic_half1"));
+ os.writeEntry("type", "cyclic");
+ putUniform(os, "value", vector::zero);
+ os.endBlock();
+ }
+ }
+
+ if (pars.yCyclic)
+ {
+ os.beginBlock("yCyclic_half0");
+ os.writeEntry("type", "cyclic");
+ os.endBlock();
+
+ os.beginBlock("yCyclic_half1");
+ os.writeEntry("type", "cyclic");
+ os.endBlock();
+ }
+ else
+ {
+ os.beginBlock("ySymmetry");
+ os.writeEntry("type", "symmetryPlane");
+ os.endBlock();
+ }
+
+ if ( pars.outer_orthog )
+ {
+ os.beginBlock("outer_inner");
+ os.writeEntry("type", "cyclicAMI");
+ os.writeEntry("neighbourPatch", "inner_outer");
+ os.endBlock();
+
+ os.beginBlock("inner_outer");
+ os.writeEntry("type", "cyclicAMI");
+ os.writeEntry("neighbourPatch", "outer_inner");
+ }
+
+ os.endBlock(); // boundaryField
+
+ IOobject::writeEndDivider(os);
+ }
+
+
+ // Pressure field
+ {
+ const scalar deflt = DEFAULT_P;
+ const char *wall_bc = "zeroGradient;\n\trho\trho";
+
+ OFstream os(casepath / pars.timeName / "p");
+ os.precision(outputPrecision);
+
+ make_header(os, "", volScalarField::typeName, "p");
+
+ os.writeEntry("dimensions", dimPressure);
+
+ os << nl;
+ putUniform(os, "internalField", deflt);
+
+ os << nl;
+ os.beginBlock("boundaryField");
+
+ // "outer"
+ {
+ os.beginBlock("outer");
+ os.writeEntry("type", "waveTransmissive");
+ os.writeEntry("gamma", 1.3);
+ os.writeEntry("fieldInf", deflt);
+ os.writeEntry("lInf", 5);
+ putUniform(os, "value", deflt);
+ os.endBlock();
+ }
+
+ tail_field(os, deflt, wall_bc, patches);
+
+ os.endBlock(); // boundaryField
+
+ IOobject::writeEndDivider(os);
+ }
+}
+
+
+// ------------------------------------------------------------------------- //
+
+void write_symmTensorField
+(
+ const word& fieldName,
+ const IjkField& fld,
+ const symmTensor& deflt, const char *wall_bc,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const dimensionSet& dims, const fileName& casepath
+)
+{
+ OFstream os(casepath / pars.timeName / fieldName);
+ os.precision(outputPrecision);
+
+ make_header(os, "", volSymmTensorField::typeName, fieldName);
+
+ os.writeEntry("dimensions", dims);
+
+ os << nl;
+ os.writeKeyword("internalField")
+ << "nonuniform List" << nl
+ << meshIndexing.nCells() << nl << token::BEGIN_LIST << nl;
+
+ for (label celli=0; celli < meshIndexing.nCells(); ++celli)
+ {
+ const labelVector& cellIdx = meshIndexing.cellIndex[celli];
+
+ if (cmptMin(cellIdx) < 0)
+ {
+ os << deflt << nl;
+ continue;
+ }
+
+ os << fld(cellIdx) << nl;
+ }
+ os << token::END_LIST << token::END_STATEMENT << nl;
+
+ os << nl;
+ os.beginBlock("boundaryField");
+
+ // outer
+ {
+ os.beginBlock("outer");
+
+ os.writeEntry("type", "inletOutlet");
+ putUniform(os, "inletValue", deflt);
+ putUniform(os, "value", deflt);
+
+ os.endBlock();
+ }
+
+ tail_field(os, deflt, wall_bc, patches);
+
+ os.endBlock(); // boundaryField
+
+ IOobject::writeEndDivider(os);
+}
+
+
+// Write a volSymmTensorField but with vectors as input.
+// The off-diagonals are zero.
+void write_symmTensorFieldV
+(
+ const word& fieldName,
+ const IjkField& fld,
+ const symmTensor& deflt, const char *wall_bc,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ const dimensionSet& dims, const fileName& casepath
+)
+{
+ OFstream os(casepath / pars.timeName / fieldName);
+ os.precision(outputPrecision);
+
+ make_header(os, "", volSymmTensorField::typeName, fieldName);
+
+ os.writeEntry("dimensions", dims);
+
+ os << nl;
+ os.writeKeyword("internalField")
+ << "nonuniform List" << nl
+ << meshIndexing.nCells() << nl << token::BEGIN_LIST << nl;
+
+ symmTensor val(symmTensor::zero);
+
+ for (label celli=0; celli < meshIndexing.nCells(); ++celli)
+ {
+ const labelVector& cellIdx = meshIndexing.cellIndex[celli];
+
+ if (cmptMin(cellIdx) < 0)
+ {
+ os << deflt << nl;
+ continue;
+ }
+
+ const vector& vec = fld(cellIdx);
+
+ val.xx() = vec.x();
+ val.yy() = vec.y();
+ val.zz() = vec.z();
+
+ os << val << nl;
+ }
+ os << token::END_LIST << token::END_STATEMENT << nl;
+
+ os << nl;
+ os.beginBlock("boundaryField");
+
+ // outer
+ {
+ os.beginBlock("outer");
+
+ os.writeEntry("type", "inletOutlet");
+ putUniform(os, "inletValue", deflt);
+ putUniform(os, "value", deflt);
+
+ os.endBlock();
+ }
+
+ tail_field(os, deflt, wall_bc, patches);
+
+ os.endBlock(); // boundaryField
+
+ IOobject::writeEndDivider(os);
+}
+
+
+// ------------------------------------------------------------------------- //
+
+void write_blocked_face_list
+(
+ const IjkField& face_block,
+ const IjkField& face_patch,
+ const IjkField& obs_count, IjkField& sub_count,
+ IjkField>& n_blocked_faces,
+ const PDRmeshArrays& meshIndexing,
+ const UList& patches,
+ double limit_par, const fileName& casepath
+)
+{
+ /* Create the lists of face numbers for faces that have already been defined as
+ belonging to (inlet) patches), and others that are found to be blocked.
+ Then write these out to set files, */
+
+ const labelVector& cellDims = meshIndexing.cellDims;
+
+ Map usedFaces;
+
+ Info<< "Number of patches: " << patches.size() << nl;
+
+ for (label facei=0; facei < meshIndexing.nFaces(); ++facei)
+ {
+ // The related i-j-k face index for the mesh face
+ const labelVector& faceIdx = meshIndexing.faceIndex[facei];
+
+ if (cmptMin(faceIdx) < 0)
+ {
+ continue;
+ }
+
+ const label ix = faceIdx.x();
+ const label iy = faceIdx.y();
+ const label iz = faceIdx.z();
+ const direction orient = meshIndexing.faceOrient[facei];
+
+ label patchId = -1;
+ scalar val(Zero);
+
+ /* A bit messy to be changing sub_count here. but there is a problem of generation
+ of subgrid flame area Xp when the flame approaches a blocked wall. the fix is to make
+ the normal component of "n" zero in the cells adjacent to the blocked face. That component
+ of n is zero when that component of sub_count i.e. ns) equals count (i.e. N). */
+ {
+ switch (orient)
+ {
+ case vector::X:
+ {
+ // face_block is the face blockage;
+ // face_patch is the patch number on the face (if any)
+ val = face_block(faceIdx).x();
+ patchId = face_patch(faceIdx).x();
+
+ if
+ (
+ val > limit_par
+ && iy < cellDims[vector::Y]
+ && iz < cellDims[vector::Z]
+ )
+ {
+ // n_blocked_faces:
+ // count of x-faces blocked for this cell
+
+ if (ix < cellDims[vector::X])
+ {
+ ++n_blocked_faces(ix,iy,iz).x();
+ sub_count(ix,iy,iz).x() = obs_count(ix,iy,iz);
+ }
+
+ if (ix > 0)
+ {
+ // And the neighbouring cell
+ ++n_blocked_faces(ix-1,iy,iz).x();
+ sub_count(ix-1,iy,iz).x() = obs_count(ix-1,iy,iz);
+ }
+ }
+ }
+ break;
+
+ case vector::Y:
+ {
+ val = face_block(faceIdx).y();
+ patchId = face_patch(faceIdx).y();
+
+ if
+ (
+ val > limit_par
+ && iz < cellDims[vector::Z]
+ && ix < cellDims[vector::X]
+ )
+ {
+ // n_blocked_faces:
+ // count of y-faces blocked for this cell
+
+ if (iy < cellDims[vector::Y])
+ {
+ ++n_blocked_faces(ix,iy,iz).y();
+ sub_count(ix,iy,iz).y() = obs_count(ix,iy,iz);
+ }
+
+ if (iy > 0)
+ {
+ // And the neighbouring cell
+ ++n_blocked_faces(ix,iy-1,iz).y();
+ sub_count(ix,iy-1,iz).y() = obs_count(ix,iy-1,iz);
+ }
+ }
+ }
+ break;
+
+ case vector::Z:
+ {
+ val = face_block(faceIdx).z();
+ patchId = face_patch(faceIdx).z();
+
+ if
+ (
+ val > limit_par
+ && ix < cellDims[vector::X]
+ && iy < cellDims[vector::Y]
+ )
+ {
+ // n_blocked_faces:
+ // count of z-faces blocked for this cell
+
+ if (iz < cellDims[vector::Z])
+ {
+ ++n_blocked_faces(ix,iy,iz).z();
+ sub_count(ix,iy,iz).z() = obs_count(ix,iy,iz);
+ }
+
+ if (iz > 0)
+ {
+ // And the neighbouring cell
+ ++n_blocked_faces(ix,iy,iz-1).z();
+ sub_count(ix,iy,iz-1).z() = obs_count(ix,iy,iz-1);
+ }
+ }
+ }
+ break;
+ }
+
+ if (patchId > 0)
+ {
+ // If this face is on a defined patch add to list
+ usedFaces(patchId).set(facei);
+ }
+ else if (val > limit_par)
+ {
+ // Add to blocked faces list
+ usedFaces(PDRpatchDef::BLOCKED_FACE).set(facei);
+ }
+ }
+ }
+
+ // Write in time or constant dir
+ const bool hasPolyMeshTimeDir = isDir(casepath/pars.timeName/"polyMesh");
+
+ const fileName setsDir =
+ (
+ casepath
+ / (hasPolyMeshTimeDir ? pars.timeName : word("constant"))
+ / fileName("polyMesh/sets")
+ );
+
+ if (!isDir(setsDir))
+ {
+ mkDir(setsDir);
+ }
+
+
+ // Create as blockedFaces Set file for each patch, including
+ // basic blocked faces
+ forAll(patches, patchi)
+ {
+ const word& patchName = patches[patchi].patchName;
+
+ OFstream os(setsDir / (patchName + "Set"));
+
+ make_header(os, "polyMesh/sets", "faceSet", patchName);
+
+ // Check for blocked faces
+ const auto& fnd = usedFaces.cfind(patchi);
+
+ if (fnd.good() && (*fnd).any())
+ {
+ os << nl << (*fnd).toc() << nl;
+ }
+ else
+ {
+ os << nl << labelList() << nl;
+ }
+
+ IOobject::writeEndDivider(os);
+ }
+
+ // Create the PDRMeshDict, listing the blocked faces sets and their patch names
+
+ {
+ DynamicList panelNames;
+
+ OFstream os(casepath / "system/PDRMeshDict");
+
+ make_header(os, "system", "dictionary", "PDRMeshDict");
+
+ os.writeEntry("blockedCells", "blockedCellsSet");
+ os << nl << "blockedFaces" << nl << token::BEGIN_LIST << nl;
+
+ for (const PDRpatchDef& p : patches)
+ {
+ const word& patchName = p.patchName;
+ const word setName = patchName + "Set";
+
+ if (p.patchType == 0) // Patch
+ {
+ os << " " << token::BEGIN_LIST
+ << setName << token::SPACE
+ << patchName << token::END_LIST
+ << nl;
+ }
+ else if (p.patchType > 0) // Panel
+ {
+ panelNames.append(setName);
+ }
+ }
+
+ os << token::END_LIST << token::END_STATEMENT << nl << nl;
+ os.beginBlock("coupledFaces");
+
+ for (const PDRpatchDef& p : patches)
+ {
+ const word& patchName = p.patchName;
+ const word setName = patchName + "Set";
+
+ if (p.patchType > 0) // Panel
+ {
+ os.beginBlock(setName);
+ os.writeEntry("wallPatchName", word(patchName + "Wall"));
+ os.writeEntry("cyclicMasterPatchName", patchName);
+ os.endBlock();
+ }
+ }
+ os.endBlock() << nl;
+
+ os.writeEntry("defaultPatch", "blockedFaces");
+
+ IOobject::writeEndDivider(os);
+
+ // Write panelList
+ OFstream(casepath / "panelList")()
+ << panelNames << token::END_STATEMENT << nl;
+ }
+}
+
+
+void write_blockedCellsSet
+(
+ const IjkField& fld,
+ const PDRmeshArrays& meshIndexing,
+ double limit_par,
+ const IjkField>& n_blocked_faces,
+ const fileName& casepath,
+ const word& listName
+)
+{
+ if (listName.empty())
+ {
+ return;
+ }
+
+ // Write in time or constant dir
+ const bool hasPolyMeshTimeDir = isDir(casepath/pars.timeName/"polyMesh");
+
+ const fileName path =
+ (
+ casepath
+ / (hasPolyMeshTimeDir ? pars.timeName : word("constant"))
+ / fileName("polyMesh/sets")
+ / listName
+ );
+
+ if (!isDir(path.path()))
+ {
+ mkDir(path.path());
+ }
+
+ bitSet blockedCell;
+
+ for (label celli=0; celli < meshIndexing.nCells(); ++celli)
+ {
+ const labelVector& cellIdx = meshIndexing.cellIndex[celli];
+
+ if (cmptMin(cellIdx) < 0)
+ {
+ continue;
+ }
+
+ if (fld(cellIdx) < limit_par)
+ {
+ blockedCell.set(celli);
+ continue;
+ }
+
+ const Vector& blocked = n_blocked_faces(cellIdx);
+
+ const label n_bfaces = cmptSum(blocked);
+
+ label n_bpairs = 0;
+
+ if (n_bfaces > 1)
+ {
+ for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
+ {
+ if (blocked[cmpt] > 1) ++n_bpairs;
+ }
+
+ #if 0
+ // Extra debugging
+ Info<<"block " << celli << " from "
+ << blocked << " -> ("
+ << n_bfaces << ' ' << n_bpairs
+ << ')' << nl;
+ #endif
+ }
+
+ if
+ (
+ n_bfaces >= pars.nFacesToBlockC
+ || n_bpairs >= pars.nPairsToBlockC
+ )
+ {
+ blockedCell.set(celli);
+ }
+ }
+
+
+ OFstream os(path);
+ make_header(os, "constant/polyMesh/sets", "cellSet", listName);
+
+ if (blockedCell.any())
+ {
+ os << blockedCell.toc();
+ }
+ else
+ {
+ os << labelList();
+ }
+
+ os << token::END_STATEMENT << nl;
+
+ IOobject::writeEndDivider(os);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRlegacy.H b/applications/utilities/preProcessing/PDRsetFields/PDRlegacy.H
new file mode 100644
index 0000000000..bbb660f04a
--- /dev/null
+++ b/applications/utilities/preProcessing/PDRsetFields/PDRlegacy.H
@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016 Shell Research Ltd.
+ Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+Namespace
+ Foam::PDRlegacy
+
+Description
+ Legacy and transitional routines
+
+SourceFiles
+ PDRlegacyMeshSpec.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PDRlegacy_H
+#define PDRlegacy_H
+
+#include "PDRblock.H"
+#include "fileName.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class ISstream;
+class PDRblock;
+
+/*---------------------------------------------------------------------------*\
+ Namespace PDRlegacy
+\*---------------------------------------------------------------------------*/
+
+namespace PDRlegacy
+{
+ void print_info(const PDRblock& block);
+
+ void read_mesh_spec(const fileName& casepath, PDRblock& pdrBlock);
+
+ void read_mesh_spec(ISstream& is, PDRblock& pdrBlock);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRlegacyMeshSpec.C b/applications/utilities/preProcessing/PDRsetFields/PDRlegacyMeshSpec.C
new file mode 100644
index 0000000000..720ec86359
--- /dev/null
+++ b/applications/utilities/preProcessing/PDRsetFields/PDRlegacyMeshSpec.C
@@ -0,0 +1,340 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016 Shell Research Ltd.
+ Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+\*---------------------------------------------------------------------------*/
+
+// Read spec that resemble this:
+//
+// xmesh
+// (
+// ( -8.18 14 0.83 )
+// ( -0.69 11 1.00 )
+// ( 0.69 14 1.20 )
+// ( 8.18 0 )
+// )
+//
+// ymesh
+// (
+// ...
+// )
+
+#include "PDRlegacy.H"
+
+// OpenFOAM includes
+#include "error.H"
+#include "IFstream.H"
+#include "stringOps.H"
+#include "OSspecific.H"
+
+#define XMESH_TAG "xmesh"
+#define YMESH_TAG "ymesh"
+#define ZMESH_TAG "zmesh"
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace PDRlegacy
+{
+namespace Detail
+{
+
+//- Simple structure for processing a mesh spec entry
+// Handles an entry with (float), (float int), or (float int float)
+//
+// This is for handling a sub-spec that resembles this:
+//
+// (
+// ( -8.18 14 0.83 )
+// ( -0.69 11 1.00 )
+// ( 0.69 14 1.20 )
+// ( 8.18 0 )
+// )
+struct pdrMeshSpecLine
+{
+ scalar knot;
+ label ndiv;
+ scalar factor;
+
+ pdrMeshSpecLine() : knot(0), ndiv(0), factor(0) {}
+
+ //- Cheap means to avoid uniform list output
+ bool operator!=(const pdrMeshSpecLine&) const
+ {
+ return false;
+ }
+};
+
+
+// Read mesh-spec entry
+Istream& operator>>(Istream& is, pdrMeshSpecLine& spec)
+{
+ spec.knot = 0;
+ spec.ndiv = 0;
+ spec.factor = 0;
+
+ is.readBegin("pdrMeshSpecLine");
+
+ // Must have a point
+ is >> spec.knot;
+
+ token tok(is);
+ if (tok.isLabel())
+ {
+ spec.ndiv = tok.labelToken();
+
+ if (spec.ndiv)
+ {
+ is >> spec.factor;
+ }
+ }
+ else
+ {
+ is.putBack(tok);
+ }
+
+ is.readEnd("pdrMeshSpecLine");
+
+ is.check(FUNCTION_NAME);
+ return is;
+}
+
+#ifdef FULLDEBUG
+// Write mesh-spec entry
+Ostream& operator<<(Ostream& os, const pdrMeshSpecLine& spec)
+{
+ os << token::BEGIN_LIST << spec.knot;
+
+ if (spec.ndiv)
+ {
+ os << token::SPACE << spec.ndiv
+ << token::SPACE << spec.factor;
+ }
+
+ os << token::END_LIST;
+
+ return os;
+}
+#endif
+
+
+void read_spec(ISstream& is, const direction cmpt, List& gridPoint)
+{
+ if (!gridPoint.empty())
+ {
+ FatalErrorInFunction
+ << "Duplicate specification of "
+ << vector::componentNames[cmpt]
+ << " grid"
+ << exit(FatalError);
+ }
+
+ List specs(is);
+
+ if (specs.size() < 2)
+ {
+ FatalErrorInFunction
+ << "Grid specification for " << vector::componentNames[cmpt]
+ << " is too small. Need at least two points!" << nl
+ << exit(FatalError);
+ }
+
+ specs.last().ndiv = 0; // safety
+
+
+ DynamicList knots;
+ DynamicList