From a60fe9c7b029877f0f23267b623ddf4f97ebe68a Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Mon, 16 Dec 2019 21:50:47 +0100 Subject: [PATCH] ENH: PDRsetFields utility (#1216) - the PDRsetFields utility processes a set of geometrical obstructions to determine the equivalent blockage effects. These fields are necessary inputs for PDRFoam calculations. After setting up the geometries, the -dry-run option can be used to generate a VTK file for diagnosis and post-processing purposes. - this is an initial release, with improvements slated for the future. NOTE - the field results may be less than fully reliable when run in single-precision. This howver does not represent a realistic restriction since the prepared fields target a combustion application which will invariably be double-precision. --- .../preProcessing/PDRsetFields/Make/files | 19 + .../preProcessing/PDRsetFields/Make/options | 14 + .../preProcessing/PDRsetFields/PDRarrays.C | 141 ++ .../preProcessing/PDRsetFields/PDRarrays.H | 215 ++ .../PDRsetFields/PDRarraysAnalyse.C | 1099 +++++++++ .../PDRsetFields/PDRarraysCalc.C | 1984 +++++++++++++++++ .../preProcessing/PDRsetFields/PDRlegacy.H | 75 + .../PDRsetFields/PDRlegacyMeshSpec.C | 340 +++ .../PDRsetFields/PDRmeshArrays.C | 262 +++ .../PDRsetFields/PDRmeshArrays.H | 131 ++ .../preProcessing/PDRsetFields/PDRparams.C | 111 + .../preProcessing/PDRsetFields/PDRparams.H | 165 ++ .../preProcessing/PDRsetFields/PDRpatchDef.C | 52 + .../preProcessing/PDRsetFields/PDRpatchDef.H | 123 + .../preProcessing/PDRsetFields/PDRsetFields.C | 351 +++ .../preProcessing/PDRsetFields/PDRsetFields.H | 142 ++ .../preProcessing/PDRsetFields/PDRutils.H | 62 + .../PDRsetFields/PDRutilsInternal.H | 221 ++ .../PDRsetFields/PDRutilsIntersect.C | 712 ++++++ .../PDRsetFields/PDRutilsOverlap.C | 767 +++++++ .../PDRsetFields/obstacles/ObstaclesDict | 179 ++ .../PDRsetFields/obstacles/PDRobstacle.C | 737 ++++++ .../PDRsetFields/obstacles/PDRobstacle.H | 477 ++++ .../PDRsetFields/obstacles/PDRobstacleI.H | 46 + .../PDRsetFields/obstacles/PDRobstacleIO.C | 353 +++ .../obstacles/PDRobstacleLegacyIO.C | 475 ++++ .../obstacles/PDRobstacleLegacyRead.C | 567 +++++ .../PDRsetFields/obstacles/PDRobstacleTypes.C | 512 +++++ .../PDRsetFields/obstacles/PDRobstacleTypes.H | 300 +++ 29 files changed, 10632 insertions(+) create mode 100644 applications/utilities/preProcessing/PDRsetFields/Make/files create mode 100644 applications/utilities/preProcessing/PDRsetFields/Make/options create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRarrays.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRarrays.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRarraysAnalyse.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRarraysCalc.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRlegacy.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRlegacyMeshSpec.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRmeshArrays.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRmeshArrays.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRparams.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRparams.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRpatchDef.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRpatchDef.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRsetFields.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRsetFields.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRutils.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRutilsInternal.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRutilsIntersect.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/PDRutilsOverlap.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/ObstaclesDict create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacle.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacle.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleI.H create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleIO.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleLegacyIO.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleLegacyRead.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleTypes.C create mode 100644 applications/utilities/preProcessing/PDRsetFields/obstacles/PDRobstacleTypes.H 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