powerLawLopesdaCosta.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2018-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 #include "powerLawLopesdaCosta.H"
28 #include "geometricOneField.H"
29 #include "fvMatrices.H"
30 #include "triSurfaceMesh.H"
31 #include "triSurfaceTools.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  namespace porosityModels
38  {
41  }
42 }
43 
44 
45 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
46 
48 (
49  const word& name,
50  const word& modelType,
51  const fvMesh& mesh,
52  const dictionary& dict
53 )
54 :
55  zoneName_(name + ":porousZone")
56 {
57  dictionary coeffs(dict.optionalSubDict(modelType + "Coeffs"));
58 
59  // Vertical direction within porous region
60  vector zDir(coeffs.lookup("zDir"));
61 
62  // Span of the search for the cells in the porous region
63  scalar searchSpan(coeffs.lookup<scalar>("searchSpan"));
64 
65  // Top surface file name defining the extent of the porous region
66  word topSurfaceFileName(coeffs.lookup("topSurface"));
67 
68  // List of the ground patches defining the lower surface
69  // of the porous region
70  List<wordRe> groundPatches(coeffs.lookup("groundPatches"));
71 
72  // Functional form of the porosity surface area per unit volume as a
73  // function of the normalised vertical position
75  (
76  Function1<scalar>::New("Sigma", dict)
77  );
78 
79  // Searchable triSurface for the top of the porous region
80  triSurfaceMesh searchSurf
81  (
82  IOobject
83  (
84  topSurfaceFileName,
85  mesh.time().constant(),
87  mesh.time()
88  )
89  );
90 
91  // Limit the porous cell search to those within the lateral and vertical
92  // extent of the top surface
93 
94  boundBox surfBounds(searchSurf.points());
95  boundBox searchBounds
96  (
97  surfBounds.min() - searchSpan*zDir, surfBounds.max()
98  );
99 
100  const pointField& C = mesh.cellCentres();
101 
102  // List of cells within the porous region
103  labelList porousCells(C.size());
104  label porousCelli = 0;
105 
106  forAll(C, celli)
107  {
108  if (searchBounds.contains(C[celli]))
109  {
110  porousCells[porousCelli++] = celli;
111  }
112  }
113 
114  porousCells.setSize(porousCelli);
115 
116  pointField start(porousCelli);
117  pointField end(porousCelli);
118 
119  forAll(porousCells, porousCelli)
120  {
121  start[porousCelli] = C[porousCells[porousCelli]];
122  end[porousCelli] = start[porousCelli] + searchSpan*zDir;
123  }
124 
125  // Field of distance between the cell centre and the corresponding point
126  // on the porous region top surface
127  scalarField zTop(porousCelli);
128 
129  // Search the porous cells for those with a corresponding point on the
130  // porous region top surface
131  List<pointIndexHit> hitInfo;
132  searchSurf.findLine(start, end, hitInfo);
133 
134  porousCelli = 0;
135 
136  forAll(porousCells, celli)
137  {
138  const pointIndexHit& hit = hitInfo[celli];
139 
140  if (hit.hit())
141  {
142  porousCells[porousCelli] = porousCells[celli];
143 
144  zTop[porousCelli] =
145  (hit.hitPoint() - C[porousCells[porousCelli]]) & zDir;
146 
147  porousCelli++;
148  }
149  }
150 
151  // Resize the porous cells list and height field
152  porousCells.setSize(porousCelli);
153  zTop.setSize(porousCelli);
154 
155  // Create a triSurface for the ground patch(s)
156  triSurface groundSurface
157  (
159  (
160  mesh.boundaryMesh(),
161  mesh.boundaryMesh().patchSet(groundPatches),
162  searchBounds
163  )
164  );
165 
166  // Combine the ground triSurfaces across all processors
167  if (Pstream::parRun())
168  {
169  List<List<labelledTri>> groundSurfaceProcTris(Pstream::nProcs());
170  List<pointField> groundSurfaceProcPoints(Pstream::nProcs());
171 
172  groundSurfaceProcTris[Pstream::myProcNo()] = groundSurface;
173  groundSurfaceProcPoints[Pstream::myProcNo()] = groundSurface.points();
174 
175  Pstream::gatherList(groundSurfaceProcTris);
176  Pstream::scatterList(groundSurfaceProcTris);
177  Pstream::gatherList(groundSurfaceProcPoints);
178  Pstream::scatterList(groundSurfaceProcPoints);
179 
180  label nTris = 0;
181  forAll(groundSurfaceProcTris, i)
182  {
183  nTris += groundSurfaceProcTris[i].size();
184  }
185 
186  List<labelledTri> groundSurfaceTris(nTris);
187  label trii = 0;
188  label offset = 0;
189  forAll(groundSurfaceProcTris, i)
190  {
191  forAll(groundSurfaceProcTris[i], j)
192  {
193  groundSurfaceTris[trii] = groundSurfaceProcTris[i][j];
194  groundSurfaceTris[trii][0] += offset;
195  groundSurfaceTris[trii][1] += offset;
196  groundSurfaceTris[trii][2] += offset;
197  trii++;
198  }
199  offset += groundSurfaceProcPoints[i].size();
200  }
201 
202  label nPoints = 0;
203  forAll(groundSurfaceProcPoints, i)
204  {
205  nPoints += groundSurfaceProcPoints[i].size();
206  }
207 
208  pointField groundSurfacePoints(nPoints);
209 
210  label pointi = 0;
211  forAll(groundSurfaceProcPoints, i)
212  {
213  forAll(groundSurfaceProcPoints[i], j)
214  {
215  groundSurfacePoints[pointi++] = groundSurfaceProcPoints[i][j];
216  }
217  }
218 
219  groundSurface = triSurface(groundSurfaceTris, groundSurfacePoints);
220  }
221 
222  // Create a searchable triSurface for the ground
223  triSurfaceSearch groundSearch(groundSurface);
224 
225  // Search the porous cells for the corresponding point on the ground
226 
227  start.setSize(porousCelli);
228  end.setSize(porousCelli);
229 
230  forAll(porousCells, porousCelli)
231  {
232  start[porousCelli] = C[porousCells[porousCelli]];
233  end[porousCelli] = start[porousCelli] - searchSpan*zDir;
234  }
235 
236  groundSearch.findLine(start, end, hitInfo);
237 
238  scalarField zBottom(porousCelli);
239 
240  forAll(porousCells, porousCelli)
241  {
242  const pointIndexHit& hit = hitInfo[porousCelli];
243 
244  if (hit.hit())
245  {
246  zBottom[porousCelli] =
247  (C[porousCells[porousCelli]] - hit.hitPoint()) & zDir;
248  }
249  }
250 
251  // Create the normalised height field
252  scalarField zNorm(zBottom/(zBottom + zTop));
253 
254  // Create the porosity surface area per unit volume zone field
255  Sigma_ = SigmaFunc->value(zNorm);
256 
257  // Create the porous region cellZone and add to the mesh cellZones
258 
259  meshCellZones& cellZones = const_cast<meshCellZones&>(mesh.cellZones());
260 
261  label zoneID = cellZones.findZoneID(zoneName_);
262 
263  if (zoneID == -1)
264  {
265  zoneID = cellZones.size();
266  cellZones.setSize(zoneID + 1);
267 
268  cellZones.set
269  (
270  zoneID,
271  new cellZone
272  (
273  zoneName_,
274  porousCells,
275  zoneID,
276  cellZones
277  )
278  );
279  }
280  else
281  {
283  << "Unable to create porous cellZone " << zoneName_
284  << ": zone already exists"
285  << abort(FatalError);
286  }
287 }
288 
289 
291 (
292  const word& name,
293  const word& modelType,
294  const fvMesh& mesh,
295  const dictionary& dict,
296  const word& dummyCellZoneName
297 )
298 :
299  powerLawLopesdaCostaZone(name, modelType, mesh, dict),
301  (
302  name,
303  modelType,
304  mesh,
305  dict,
307  ),
308  Cd_(coeffs_.lookup<scalar>("Cd")),
309  C1_(coeffs_.lookup<scalar>("C1")),
310  rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho"))
311 {}
312 
313 
314 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
315 
317 {}
318 
319 
320 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
321 
322 const Foam::scalarField&
324 {
325  return Sigma_;
326 }
327 
328 
330 {}
331 
332 
334 (
335  const volVectorField& U,
336  const volScalarField& rho,
337  const volScalarField& mu,
338  vectorField& force
339 ) const
340 {
341  scalarField Udiag(U.size(), 0.0);
342  const scalarField& V = mesh_.V();
343 
344  apply(Udiag, V, rho, U);
345 
346  force = Udiag*U;
347 }
348 
349 
351 (
352  fvVectorMatrix& UEqn
353 ) const
354 {
355  const vectorField& U = UEqn.psi();
356  const scalarField& V = mesh_.V();
357  scalarField& Udiag = UEqn.diag();
358 
359  if (UEqn.dimensions() == dimForce)
360  {
361  const volScalarField& rho =
362  mesh_.lookupObject<volScalarField>(rhoName_);
363 
364  apply(Udiag, V, rho, U);
365  }
366  else
367  {
368  apply(Udiag, V, geometricOneField(), U);
369  }
370 }
371 
372 
374 (
375  fvVectorMatrix& UEqn,
376  const volScalarField& rho,
377  const volScalarField& mu
378 ) const
379 {
380  const vectorField& U = UEqn.psi();
381  const scalarField& V = mesh_.V();
382  scalarField& Udiag = UEqn.diag();
383 
384  apply(Udiag, V, rho, U);
385 }
386 
387 
389 (
390  const fvVectorMatrix& UEqn,
391  volTensorField& AU
392 ) const
393 {
394  const vectorField& U = UEqn.psi();
395 
396  if (UEqn.dimensions() == dimForce)
397  {
398  const volScalarField& rho =
399  mesh_.lookupObject<volScalarField>(rhoName_);
400 
401  apply(AU, rho, U);
402  }
403  else
404  {
405  apply(AU, geometricOneField(), U);
406  }
407 }
408 
409 
411 {
412  os << indent << name_ << endl;
413  dict_.write(os);
414 
415  return true;
416 }
417 
418 
419 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find first intersection on segment from start to end.
Run-time selectable general function of one variable.
Definition: Function1.H:52
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
Graphite solid properties.
Definition: C.H:48
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: MeshZones.C:341
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual Ostream & write(const char)=0
Write character.
defineTypeNameAndDebug(powerLawLopesdaCosta, 0)
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
Variant of the power law porosity model with spatially varying drag coefficient.
virtual void correct(fvVectorMatrix &UEqn) const
Add resistance.
U
Definition: pEqn.H:72
const meshCellZones & cellZones() const
Return cell zones.
Definition: polyMesh.H:501
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
const scalarField & Sigma() const
Return the porosity surface area per unit volume zone field.
addToRunTimeSelectionTable(porosityModel, powerLawLopesdaCosta, mesh)
powerLawLopesdaCosta(const word &name, const word &modelType, const fvMesh &mesh, const dictionary &dict, const word &cellZoneName)
virtual tmp< pointField > points() const
Get the points that define the surface.
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:290
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:372
labelHashSet patchSet(const UList< wordRe > &patchNames, const bool warnNotFound=true, const bool usePatchGroups=true) const
Return the set of patch IDs corresponding to the given names.
Macros for easy insertion into run-time selection tables.
powerLawLopesdaCostaZone(const word &name, const word &modelType, const fvMesh &mesh, const dictionary &dict)
Constructor.
A surface geometry formed of discrete facets, e.g. triangles and/or quadrilaterals, defined in a file using formats such as Wavefront OBJ, or stereolithography STL.
Helper class to search on triSurface.
const Point & hitPoint() const
Return hit point.
virtual void calcForce(const volVectorField &U, const volScalarField &rho, const volScalarField &mu, vectorField &force) const
Calculate the porosity force.
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &info) const
const dictionary & optionalSubDict(const word &) const
Find and return a sub-dictionary if found.
Definition: dictionary.C:1060
static const word & geometryDir()
Return the geometry directory name.
A class for handling words, derived from string.
Definition: word.H:59
label nPoints
const word & constant() const
Return constant name.
Definition: TimePaths.H:123
bool hit() const
Is there a hit.
const Field< PointType > & points() const
Return reference to global points.
virtual void calcTransformModelData()
Transform the model data wrt mesh changes.
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
const dimensionSet dimForce
const vectorField & cellCentres() const
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
A subset of mesh cells.
Definition: cellZone.H:61
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void setSize(const label)
Reset size of List.
Definition: List.C:281
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:95
A special matrix type and solver, designed for finite volume solutions of scalar equations.
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
static triSurface triangulate(const polyBoundaryMesh &mBesh, const labelHashSet &includePatches, const bool verbose=false)
Simple triangulation of (selected patches of) boundaryMesh. Needs.
scalarField & diag()
Definition: lduMatrix.C:186
Triangulated surface description with patch information.
Definition: triSurface.H:66
const word zoneName_
Automatically generated zone name for this porous zone.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
Top level model for porosity models.
Definition: porosityModel.H:54
Namespace for OpenFOAM.
const dimensionSet & dimensions() const
Definition: fvMatrix.H:295