a60fe9c7b0
- 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.
568 lines
16 KiB
C
568 lines
16 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / 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 <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "PDRsetFields.H"
|
|
#include "PDRobstacle.H"
|
|
#include "volumeType.H"
|
|
|
|
using namespace Foam;
|
|
using namespace Foam::constant;
|
|
|
|
#undef USE_ZERO_INSTANCE_GROUPS
|
|
// #define USE_ZERO_INSTANCE_GROUPS
|
|
|
|
|
|
// Counting
|
|
Foam::labelPair Foam::PDRlegacy::readObstacleFiles
|
|
(
|
|
const fileName& obsFileDir,
|
|
const wordList& obsFileNames,
|
|
Map<obstacleGrouping>& groups
|
|
)
|
|
{
|
|
// Default group with single instance and position (0,0,0)
|
|
groups(0).clear();
|
|
groups(0).append(point::zero);
|
|
|
|
string buffer;
|
|
|
|
if (!obsFileNames.empty())
|
|
{
|
|
Info<< "Counting groups in obstacle files" << nl;
|
|
}
|
|
for (const word& inputFile : obsFileNames)
|
|
{
|
|
Info<< " file: " << inputFile << nl;
|
|
|
|
fileName path = (obsFileDir / inputFile);
|
|
|
|
IFstream is(path);
|
|
|
|
while (is.good())
|
|
{
|
|
// Process each line of obstacle files
|
|
is.getLine(buffer);
|
|
|
|
const auto firstCh = buffer.find_first_not_of(" \t\n\v\f\r");
|
|
|
|
if
|
|
(
|
|
firstCh == std::string::npos
|
|
|| buffer[firstCh] == '#'
|
|
)
|
|
{
|
|
// Empty line or comment line
|
|
continue;
|
|
}
|
|
|
|
int typeId;
|
|
double x, y, z; // Double (not scalar) to match sscanf spec
|
|
|
|
if
|
|
(
|
|
sscanf(buffer.c_str(), "%d %lf %lf %lf", &typeId, &x, &y, &z)<4
|
|
|| typeId == 0
|
|
|| typeId == PDRobstacle::MESH_PLANE
|
|
)
|
|
{
|
|
continue;
|
|
}
|
|
|
|
x *= pars.scale;
|
|
y *= pars.scale;
|
|
z *= pars.scale;
|
|
|
|
const label groupId = typeId / 100;
|
|
typeId %= 100;
|
|
|
|
if (typeId == PDRobstacle::OLD_INLET)
|
|
{
|
|
Info<< "Ignored old-inlet type" << nl;
|
|
continue;
|
|
}
|
|
else if (typeId == PDRobstacle::GRATING && pars.ignoreGratings)
|
|
{
|
|
Info<< "Ignored grating" << nl;
|
|
continue;
|
|
}
|
|
|
|
if (typeId == 0)
|
|
{
|
|
// Defining a group location
|
|
groups(groupId).append(x, y, z);
|
|
}
|
|
else if (PDRobstacle::isCylinder(typeId))
|
|
{
|
|
// Increment cylinder count for the group
|
|
groups(groupId).addCylinder();
|
|
}
|
|
else
|
|
{
|
|
// Increment obstacle count for the group
|
|
groups(groupId).addObstacle();
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
label nTotalObs = 0;
|
|
label nTotalCyl = 0;
|
|
|
|
label nMissedObs = 0;
|
|
label nMissedCyl = 0;
|
|
|
|
forAllConstIters(groups, iter)
|
|
{
|
|
const auto& group = iter.val();
|
|
|
|
nTotalObs += group.nTotalObstacle();
|
|
nTotalCyl += group.nTotalCylinder();
|
|
|
|
if (group.empty())
|
|
{
|
|
nMissedObs += group.nObstacle();
|
|
nMissedCyl += group.nCylinder();
|
|
}
|
|
}
|
|
|
|
for (const label groupId : groups.sortedToc())
|
|
{
|
|
const auto& group = groups[groupId];
|
|
|
|
if (groupId)
|
|
{
|
|
if (group.size())
|
|
{
|
|
Info<< "Found " << group.size()
|
|
<< " instances of group " << groupId << " ("
|
|
<< group.nObstacle() << " obstacles "
|
|
<< group.nCylinder() << " cylinders)"
|
|
<< nl;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// The group 0 is for ungrouped obstacles
|
|
Info<< "Found "
|
|
<< group.nObstacle() << " obstacles "
|
|
<< group.nCylinder() << " cylinders not in groups" << nl;
|
|
}
|
|
}
|
|
|
|
Info<< "Number of obstacles: "
|
|
<< (nTotalObs + nTotalCyl) << " ("
|
|
<< nTotalCyl << " cylinders)" << nl;
|
|
|
|
if (nMissedObs + nMissedCyl)
|
|
{
|
|
#ifdef USE_ZERO_INSTANCE_GROUPS
|
|
|
|
nTotalObs += nMissedObs;
|
|
nTotalCyl += nMissedCyl;
|
|
Info<< "Adding " << (nMissedObs + nMissedCyl)
|
|
<< " obstacles in groups without instances to default group" << nl;
|
|
|
|
#else
|
|
|
|
Warning
|
|
<< nl << "Found " << (nMissedObs + nMissedCyl)
|
|
<< " obstacles in groups without instances" << nl << nl;
|
|
|
|
if (pars.debugLevel > 1)
|
|
{
|
|
for (const label groupId : groups.sortedToc())
|
|
{
|
|
const auto& group = groups[groupId];
|
|
|
|
if
|
|
(
|
|
groupId && group.empty()
|
|
&& (group.nObstacle() || group.nCylinder())
|
|
)
|
|
{
|
|
Info<< " Group " << groupId << " ("
|
|
<< group.nObstacle() << " obstacles "
|
|
<< group.nCylinder() << " cylinders)"
|
|
<< nl;
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
|
|
return labelPair(nTotalObs, nTotalCyl);
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::PDRlegacy::readObstacleFiles
|
|
(
|
|
const fileName& obsFileDir,
|
|
const wordList& obsFileNames,
|
|
const Map<obstacleGrouping>& groups,
|
|
const boundBox& meshBb,
|
|
|
|
DynamicList<PDRobstacle>& blocks,
|
|
DynamicList<PDRobstacle>& cylinders
|
|
)
|
|
{
|
|
// Catch programming errors
|
|
if (!groups.found(0))
|
|
{
|
|
FatalErrorInFunction
|
|
<< "No default group 0 defined!" << nl
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
scalar totVolume = 0;
|
|
label nOutside = 0;
|
|
label nProtruding = 0;
|
|
|
|
scalar shift = pars.obs_expand;
|
|
|
|
string buffer;
|
|
|
|
if (!obsFileNames.empty())
|
|
{
|
|
Info<< "Reading obstacle files" << nl;
|
|
}
|
|
|
|
for (const word& inputFile : obsFileNames)
|
|
{
|
|
Info<< " file: " << inputFile << nl;
|
|
|
|
fileName path = (obsFileDir / inputFile);
|
|
|
|
IFstream is(path);
|
|
|
|
label lineNo = 0;
|
|
while (is.good())
|
|
{
|
|
// Process each line of obstacle files
|
|
++lineNo;
|
|
is.getLine(buffer);
|
|
|
|
const auto firstCh = buffer.find_first_not_of(" \t\n\v\f\r");
|
|
|
|
if
|
|
(
|
|
firstCh == std::string::npos
|
|
|| buffer[firstCh] == '#'
|
|
)
|
|
{
|
|
// Empty line or comment line
|
|
continue;
|
|
}
|
|
|
|
// Quick reject
|
|
|
|
int typeId; // Int (not label) to match sscanf spec
|
|
double x, y, z; // Double (not scalar) to match sscanf spec
|
|
|
|
if
|
|
(
|
|
sscanf(buffer.c_str(), "%d %lf %lf %lf", &typeId, &x, &y, &z) < 4
|
|
|| typeId == 0
|
|
|| typeId == PDRobstacle::MESH_PLANE
|
|
)
|
|
{
|
|
continue;
|
|
}
|
|
|
|
int groupId = typeId / 100;
|
|
typeId %= 100;
|
|
|
|
if
|
|
(
|
|
typeId == PDRobstacle::OLD_INLET
|
|
|| (typeId == PDRobstacle::GRATING && pars.ignoreGratings)
|
|
)
|
|
{
|
|
// Silent - already warned during counting
|
|
continue;
|
|
}
|
|
|
|
if (typeId == 0)
|
|
{
|
|
// Group location - not an obstacle
|
|
continue;
|
|
}
|
|
|
|
if (!groups.found(groupId))
|
|
{
|
|
// Catch programming errors.
|
|
// - group should be there after the previous read
|
|
Warning
|
|
<< "Encountered undefined group: " << groupId << nl;
|
|
continue;
|
|
}
|
|
|
|
#ifdef USE_ZERO_INSTANCE_GROUPS
|
|
const obstacleGrouping& group =
|
|
(
|
|
groups[groups[groupId].size() ? groupId : 0]
|
|
);
|
|
#else
|
|
const obstacleGrouping& group = groups[groupId];
|
|
#endif
|
|
|
|
// Add the obstacle to the list with different position
|
|
// offsets according to its group. Non-group obstacles
|
|
// are treated as group 0, which has a single instance
|
|
// with position (0,0,0) and are added only once.
|
|
|
|
PDRobstacle scanObs;
|
|
|
|
if
|
|
(
|
|
!scanObs.setFromLegacy
|
|
(
|
|
(groupId * 100) + typeId,
|
|
buffer,
|
|
lineNo,
|
|
inputFile
|
|
)
|
|
)
|
|
{
|
|
continue;
|
|
}
|
|
|
|
scanObs.scale(pars.scale);
|
|
|
|
// Ignore anything below minimum width
|
|
if (scanObs.tooSmall(pars.min_width))
|
|
{
|
|
continue;
|
|
}
|
|
|
|
|
|
for (const point& origin : group)
|
|
{
|
|
// A different (very small) shift for each obstacle
|
|
// so that faces cannot be coincident
|
|
shift += floatSMALL;
|
|
const scalar shift2 = shift * 2.0;
|
|
|
|
switch (typeId)
|
|
{
|
|
case PDRobstacle::CYLINDER:
|
|
{
|
|
// Make a copy
|
|
PDRobstacle obs(scanObs);
|
|
|
|
// Offset for the position of the group
|
|
obs.pt += origin;
|
|
|
|
// Shift the end outwards so, if exactly on
|
|
// cell boundary, now overlap cell.
|
|
// So included in Aw.
|
|
obs.pt -= point::uniform(shift);
|
|
obs.len() += shift2;
|
|
obs.dia() -= floatSMALL;
|
|
|
|
|
|
// Trim against the mesh bounds
|
|
// - ignore if it doesn't overlap
|
|
const volumeType vt = obs.trim(meshBb);
|
|
|
|
switch (vt)
|
|
{
|
|
case volumeType::OUTSIDE:
|
|
++nOutside;
|
|
continue; // Can ignore the rest
|
|
break;
|
|
|
|
case volumeType::MIXED:
|
|
++nProtruding;
|
|
break;
|
|
|
|
default:
|
|
break;
|
|
}
|
|
|
|
// Later for position sorting
|
|
switch (obs.orient)
|
|
{
|
|
case vector::X:
|
|
obs.sortBias = obs.len();
|
|
break;
|
|
case vector::Y:
|
|
obs.sortBias = 0.5*obs.dia();
|
|
break;
|
|
case vector::Z:
|
|
obs.sortBias = 0.5*obs.dia();
|
|
break;
|
|
}
|
|
|
|
totVolume += obs.volume();
|
|
cylinders.append(obs);
|
|
|
|
break;
|
|
}
|
|
|
|
case PDRobstacle::DIAG_BEAM:
|
|
{
|
|
// Make a copy
|
|
PDRobstacle obs(scanObs);
|
|
|
|
// Offset for the position of the group
|
|
obs.pt += origin;
|
|
|
|
// Shift the end outwards so, if exactly on
|
|
// cell boundary, now overlap cell.
|
|
// So included in Aw.
|
|
obs.pt -= point::uniform(shift);
|
|
obs.len() += shift2;
|
|
obs.wa += shift2;
|
|
obs.wb += shift2;
|
|
|
|
totVolume += obs.volume();
|
|
cylinders.append(obs);
|
|
|
|
break;
|
|
}
|
|
|
|
case PDRobstacle::CUBOID_1: // Cuboid "Type 1"
|
|
case PDRobstacle::LOUVRE_BLOWOFF: // Louvred wall or blow-off panel
|
|
case PDRobstacle::CUBOID: // Cuboid
|
|
case PDRobstacle::WALL_BEAM: // Beam against wall (treated here as normal cuboid)
|
|
case PDRobstacle::GRATING: // Grating
|
|
case PDRobstacle::RECT_PATCH: // Inlet, outlet or ather b.c. (rectangular)
|
|
{
|
|
// Make a copy
|
|
PDRobstacle obs(scanObs);
|
|
|
|
// Offset for the position of the group
|
|
obs.pt += origin;
|
|
|
|
if (typeId == PDRobstacle::GRATING)
|
|
{
|
|
if (obs.slat_width <= 0)
|
|
{
|
|
obs.slat_width = pars.def_grating_slat_w;
|
|
}
|
|
}
|
|
|
|
// Shift the end outwards so, if exactly on
|
|
// cell boundary, now overlap cell.
|
|
// So included in Aw.
|
|
obs.pt -= point::uniform(shift);
|
|
obs.span += point::uniform(shift2);
|
|
|
|
|
|
// Trim against the mesh bounds
|
|
// - ignore if it doesn't overlap
|
|
const volumeType vt = obs.trim(meshBb);
|
|
|
|
switch (vt)
|
|
{
|
|
case volumeType::OUTSIDE:
|
|
++nOutside;
|
|
continue; // Can ignore the rest
|
|
break;
|
|
|
|
case volumeType::MIXED:
|
|
++nProtruding;
|
|
break;
|
|
|
|
default:
|
|
break;
|
|
}
|
|
|
|
totVolume += obs.volume();
|
|
|
|
blocks.append(obs);
|
|
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (nOutside || nProtruding)
|
|
{
|
|
Info<< "Warning: " << nOutside << " obstacles outside domain, "
|
|
<< nProtruding << " obstacles partly outside domain" << nl;
|
|
}
|
|
}
|
|
|
|
|
|
// #ifdef FULLDEBUG
|
|
// Info<< blocks << nl << cylinders << nl;
|
|
// #endif
|
|
|
|
return totVolume;
|
|
}
|
|
|
|
|
|
Foam::scalar Foam::PDRobstacle::legacyReadFiles
|
|
(
|
|
const fileName& obsFileDir,
|
|
const wordList& obsFileNames,
|
|
const boundBox& meshBb,
|
|
DynamicList<PDRobstacle>& blocks,
|
|
DynamicList<PDRobstacle>& cylinders
|
|
)
|
|
{
|
|
// Still just with legacy reading
|
|
|
|
// Count the obstacles and get the group locations
|
|
Map<PDRlegacy::obstacleGrouping> groups;
|
|
|
|
const labelPair obsCounts =
|
|
PDRlegacy::readObstacleFiles(obsFileDir, obsFileNames, groups);
|
|
|
|
const label nObstacle = obsCounts.first();
|
|
const label nCylinder = obsCounts.second();
|
|
|
|
// Info<< "grouping: " << groups << endl;
|
|
|
|
if (!nObstacle && !nCylinder)
|
|
{
|
|
FatalErrorInFunction
|
|
<< "No obstacles in domain" << nl
|
|
<< exit(FatalError);
|
|
}
|
|
|
|
blocks.clear();
|
|
blocks.reserve(4 * max(4, nObstacle));
|
|
|
|
cylinders.clear();
|
|
cylinders.reserve(4 * max(4, nCylinder));
|
|
|
|
return PDRlegacy::readObstacleFiles
|
|
(
|
|
obsFileDir, obsFileNames, groups,
|
|
meshBb,
|
|
blocks,
|
|
cylinders
|
|
);
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|