PDRarraysCalc.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2016 Shell Research Ltd.
9 Copyright (C) 2019-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "PDRarrays.H"
30#include "PDRblock.H"
31#include "PDRpatchDef.H"
32#include "PDRmeshArrays.H"
33#include "PDRparams.H"
34
35#include "PDRsetFields.H"
36
37#include "bitSet.H"
38#include "DynamicList.H"
39#include "dimensionSet.H"
40#include "symmTensor.H"
41#include "SquareMatrix.H"
42#include "IjkField.H"
43#include "MinMax.H"
44#include "volFields.H"
45#include "OFstream.H"
46#include "OSspecific.H"
47
48#ifndef FULLDEBUG
49#ifndef NDEBUG
50#define NDEBUG
51#endif
52#endif
53#include <cassert>
54
55// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
56
57namespace
58{
59
60// A good ijk index has all components >= 0
61static inline bool isGoodIndex(const Foam::labelVector& idx)
62{
63 return (idx.x() >= 0 && idx.y() >= 0 && idx.z() >= 0);
64}
65
66} // End anonymous namespace
67
68
69// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
70
71using namespace Foam;
72
73static Foam::HashTable<Foam::string> fieldNotes
74({
75 { "Lobs", "" },
76 { "Aw", "surface area per unit volume" },
77 { "CR", "" },
78 { "CT", "" },
79 { "N", "" },
80 { "ns", "" },
81 { "Nv", "" },
82 { "nsv", "" },
83 { "Bv", "area blockage" },
84 { "B", "area blockage" },
85 { "betai", "" },
86 { "Blong", "longitudinal blockage" },
87 { "Ep", "1/Lobs" },
88});
89
90
91// calc_fields
92
93
94// Local Functions
95/*
96// calc_drag_etc
97make_header
98tail_field
99write_scalarField
100write_uniformField
101write_symmTensorField
102write_pU_fields
103write_blocked_face_list
104write_blockedCellsSet
105*/
106
107// Somewhat similar to what the C-fprintf would have had
108static constexpr unsigned outputPrecision = 8;
109
110void calc_drag_etc
111(
112 double brs, double brr, bool blocked,
113 double surr_br, double surr_dr,
114 scalar* drags_p, scalar* dragr_p,
115 double count,
116 scalar* cbdi_p,
117 double cell_vol
118);
119
120
121void write_scalarField
122(
123 const word& fieldName, const IjkField<scalar>& fld,
124 const scalar& deflt, const scalarMinMax& limits, const char *wall_bc,
125 const PDRmeshArrays& meshIndexing,
127 const dimensionSet& dims, const fileName& casepath
128);
129
130void write_uniformField
131(
132 const word& fieldName, const scalar& deflt, const char *wall_bc,
133 const PDRmeshArrays& meshIndexing,
135 const dimensionSet& dims, const fileName& casepath
136);
137
138void write_pU_fields
139(
140 const PDRmeshArrays& meshIndexing,
142 const fileName& casepath
143);
144
145void write_symmTensorField
146(
147 const word& fieldName, const IjkField<symmTensor>& fld,
148 const symmTensor& deflt, const char *wall_bc,
149 const PDRmeshArrays& meshIndexing,
151 const dimensionSet& dims, const fileName& casepath
152);
153
154void write_symmTensorFieldV
155(
156 const word& fieldName, const IjkField<vector>& fld,
157 const symmTensor& deflt, const char *wall_bc,
158 const PDRmeshArrays& meshIndexing,
160 const dimensionSet& dims, const fileName& casepath
161);
162
163void write_blocked_face_list
164(
165 const IjkField<vector>& face_block,
166 const IjkField<labelVector>& face_patch,
167 const IjkField<scalar>& obs_count,
168 IjkField<vector>& sub_count,
169 IjkField<Vector<direction>>& n_blocked_faces,
170 const PDRmeshArrays& meshIndexing,
172 double limit_par, const fileName& casepath
173);
174
175void write_blockedCellsSet
176(
177 const IjkField<scalar>& fld,
178 const PDRmeshArrays& meshIndexing, double limit_par,
179 const IjkField<Vector<direction>>& n_blocked_faces,
180 const fileName& casepath,
181 const word& listName
182);
183
184
185// The average values of surrounding an array position
186static inline scalar averageSurrounding
187(
188 const SquareMatrix<scalar>& mat,
189 const label i,
190 const label j
191)
192{
193 return
194 (
195 mat(i,j) + mat(i,j+1) + mat(i,j+2)
196 + mat(i+1,j) /* centre */ + mat(i+1,j+2)
197 + mat(i+2,j) + mat(i+2,j+1) + mat(i+2,j+2)
198 ) / 8.0; // Average
199}
200
201
202// Helper
203template<class Type>
204static inline Ostream& putUniform(Ostream& os, const word& key, const Type& val)
205{
206 os.writeKeyword(key)
207 << word("uniform") << token::SPACE
208 << val << token::END_STATEMENT << nl;
209 return os;
210}
211
212
213static void make_header
214(
215 Ostream& os,
216 const fileName& location,
217 const word& clsName,
218 const word& object
219)
220{
221 string note = fieldNotes(object);
222
223 IOobject::writeBanner(os);
224
225 os << "FoamFile\n{\n"
226 << " version 2.0;\n"
227 << " format ascii;\n"
228 << " class " << clsName << ";\n";
229
230 if (!note.empty())
231 {
232 os << " note " << note << ";\n";
233 }
234
235 if (!location.empty())
236 {
237 os << " location " << location << ";\n";
238 }
239
240 os << " object " << object << ";\n"
241 << "}\n";
242
243 IOobject::writeDivider(os) << nl;
244}
245
246
248(
249 PDRarrays& arr,
250 const PDRmeshArrays& meshIndexing,
251 const fileName& casepath,
253)
254{
255 if (isNull(arr.block()))
256 {
258 << "No PDRblock set" << nl
259 << exit(FatalError);
260 }
261
262 const PDRblock& pdrBlock = arr.block();
263
264 const labelVector& cellDims = meshIndexing.cellDims;
265 const labelVector& faceDims = meshIndexing.faceDims;
266
267 const int xdim = faceDims.x();
268 const int ydim = faceDims.y();
269 const int zdim = faceDims.z();
270 const scalar maxCT = pars.maxCR * pars.cb_r;
271
272
273 // Later used to store the total effective blockage ratio per cell/direction
275
277
281
283
285 IjkField<vector>& sub_count = arr.sub_count; // ns. Later used to hold longitudinal blockage
287
290
291 // Lobs. Later used for initial Ep
293
294 Info<< "Calculating fields" << nl;
295
296 // Local scratch arrays
297
298 // The turbulance generation field CT.
299 // Later used to to hold the beta_i in tensor form
300 IjkField<vector> cbdi(cellDims, Zero);
301
302
303 // For 2D addressing it is convenient to just use the max dimension
304 // and avoid resizing when handling each direction.
305
306 // Dimension of the cells and a layer of surrounding halo cells
307 const labelVector surrDims = (faceDims + labelVector::uniform(2));
308
309 // Max addressing dimensions
310 const label maxDim = cmptMax(surrDims);
311
312 // Blockage-ratio correction to the drag
313 //
314 // neiBlock:
315 // 2-D for averaging the blockage ratio of neighbouring cells.
316 // It extends one cell outside the domain in each direction,
317 // so the indices are offset by 1.
318 // neiDrag:
319 // 2-D array for averaging the drag ratio of neighbouring cells
320
321 SquareMatrix<scalar> neiBlock(maxDim, Zero);
322 SquareMatrix<scalar> neiDrag(maxDim, Zero);
323
324 // X blockage, drag
325
326 for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
327 {
328 for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
329 {
330 for (label iz = 0; iz <= zdim; ++iz)
331 {
332 const label izz =
333 (iz == 0 ? 0 : iz == zdim ? zdim - 2 : iz - 1);
334
335 neiBlock(iy+1, iz) =
336 (
337 area_block_s(ix,iy,izz).x()
338 + area_block_r(ix,iy,izz).x()
339 );
340
341 neiDrag(iy+1, iz) =
342 (
343 drag_s(ix,iy,izz).xx() * pars.cd_s
344 + drag_r(ix,iy,izz).x() * pars.cd_r
345 );
346 }
347 }
348 for (label iz = 0; iz < surrDims.z(); ++iz)
349 {
350 if (pars.yCyclic)
351 {
352 // Cyclic in y
353 neiBlock(0, iz) = neiBlock(cellDims.y(), iz);
354 neiDrag(0, iz) = neiDrag(cellDims.y(), iz);
355 neiBlock(ydim, iz) = neiBlock(1, iz);
356 neiDrag(ydim, iz) = neiDrag(1, iz);
357 }
358 else
359 {
360 neiBlock(0, iz) = neiBlock(1, iz);
361 neiDrag(0, iz) = neiDrag(1, iz);
362 neiBlock(ydim, iz) = neiBlock(cellDims.y(), iz);
363 neiDrag(ydim, iz) = neiDrag(cellDims.y(), iz);
364 }
365 }
366
367 for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
368 {
369 for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
370 {
371 const scalar cell_vol = pdrBlock.V(ix,iy,iz);
372
373 const scalar surr_br = averageSurrounding(neiBlock, iy, iz);
374 const scalar surr_dr = averageSurrounding(neiDrag, iy, iz);
375
376 calc_drag_etc
377 (
378 area_block_s(ix,iy,iz).x(),
379 area_block_r(ix,iy,iz).x(),
380 dirn_block(ix,iy,iz).x(),
381 surr_br, surr_dr,
382 &(drag_s(ix,iy,iz).xx()),
383 &(drag_r(ix,iy,iz).x()),
384 obs_count(ix,iy,iz),
385 &(cbdi(ix,iy,iz).x()),
386 cell_vol
387 );
388 }
389 }
390 }
391
392
393 // Y blockage, drag
394
395 neiBlock = Zero;
396 neiDrag = Zero;
397
398 for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
399 {
400 for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
401 {
402 for (label ix = 0; ix <= xdim; ++ix)
403 {
404 const label ixx =
405 (ix == 0 ? 0 : ix == xdim ? xdim - 2 : ix - 1);
406
407 neiBlock(iz+1, ix) =
408 (
409 area_block_s(ixx,iy,iz).y()
410 + area_block_r(ixx,iy,iz).y()
411 );
412 neiDrag(iz+1, ix) =
413 (
414 drag_s(ixx,iy,iz).yy() * pars.cd_s
415 + drag_r(ixx,iy,iz).y() * pars.cd_r
416 );
417 }
418 }
419 for (label ix = 0; ix < surrDims.x(); ++ix)
420 {
421 neiBlock(0, ix) = neiBlock(1, ix);
422 neiDrag(0, ix) = neiDrag(1, ix);
423 neiBlock(zdim, ix) = neiBlock(cellDims.z(), ix);
424 neiDrag(zdim, ix) = neiDrag(cellDims.z(), ix);
425 }
426
427 for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
428 {
429 for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
430 {
431 const scalar cell_vol = pdrBlock.V(ix,iy,iz);
432
433 const scalar surr_br = averageSurrounding(neiBlock, iz, ix);
434 const scalar surr_dr = averageSurrounding(neiDrag, iz, ix);
435
436 calc_drag_etc
437 (
438 area_block_s(ix,iy,iz).y(),
439 area_block_r(ix,iy,iz).y(),
440 dirn_block(ix,iy,iz).y(),
441 surr_br, surr_dr,
442 &(drag_s(ix,iy,iz).yy()),
443 &(drag_r(ix,iy,iz).y()),
444 obs_count(ix,iy,iz),
445 &(cbdi(ix,iy,iz).y()),
446 cell_vol
447 );
448 }
449 }
450 }
451
452
453 // Z blockage, drag
454
455 neiBlock = Zero;
456 neiDrag = Zero;
457
458 for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
459 {
460 for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
461 {
462 for (label iy = 0; iy <= ydim; ++iy)
463 {
464 label iyy;
465
466 if (pars.yCyclic)
467 {
468 iyy = (iy == 0 ? ydim - 2 : iy == ydim ? 0 : iy - 1);
469 }
470 else
471 {
472 iyy = (iy == 0 ? 0 : iy == ydim ? ydim - 2 : iy - 1);
473 }
474
475 neiBlock(ix+1, iy) =
476 (
477 area_block_s(ix,iyy,iz).z()
478 + area_block_r(ix,iyy,iz).z()
479 );
480 neiDrag(ix+1, iy) =
481 (
482 drag_s(ix,iyy,iz).zz() * pars.cd_s
483 + drag_r(ix,iyy,iz).z() * pars.cd_r
484 );
485 }
486 }
487 for (label iy = 0; iy < surrDims.y(); ++iy)
488 {
489 neiBlock(0, iy) = neiBlock(1, iy);
490 neiDrag(0, iy) = neiDrag(1, iy);
491 neiBlock(xdim, iy) = neiBlock(cellDims.x(), iy);
492 neiDrag(xdim, iy) = neiDrag(cellDims.x(), iy);
493 }
494
495 for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
496 {
497 for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
498 {
499 const scalar cell_vol = pdrBlock.V(ix,iy,iz);
500
501 const scalar surr_br = averageSurrounding(neiBlock, ix, iy);
502 const scalar surr_dr = averageSurrounding(neiDrag, ix, iy);
503
504 calc_drag_etc
505 (
506 area_block_s(ix,iy,iz).z(),
507 area_block_r(ix,iy,iz).z(),
508 dirn_block(ix,iy,iz).z(),
509 surr_br, surr_dr,
510 &(drag_s(ix,iy,iz).zz()),
511 &(drag_r(ix,iy,iz).z()),
512 obs_count(ix,iy,iz),
513 &(cbdi(ix,iy,iz).z()),
514 cell_vol
515 );
516 }
517 }
518 }
519
520 neiBlock.clear();
521 neiDrag.clear();
522
523
524 // Calculate other parameters
525
526 for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
527 {
528 for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
529 {
530 for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
531 {
532 const scalar dx = pdrBlock.dx(ix);
533 const scalar dy = pdrBlock.dy(iy);
534 const scalar dz = pdrBlock.dz(iz);
535 const scalar cell_vol = pdrBlock.V(ix, iy, iz);
536 const scalar cell_size = pdrBlock.width(ix, iy, iz);
537
538 drag_s(ix,iy,iz).xy() *= pars.cd_s;
539 drag_s(ix,iy,iz).xz() *= pars.cd_s;
540 drag_s(ix,iy,iz).yz() *= pars.cd_s;
541
542 if (drag_s(ix,iy,iz).xx() > pars.maxCR) { drag_s(ix,iy,iz).xx() = pars.maxCR; } ;
543 if (drag_s(ix,iy,iz).yy() > pars.maxCR) { drag_s(ix,iy,iz).yy() = pars.maxCR; } ;
544 if (drag_s(ix,iy,iz).zz() > pars.maxCR) { drag_s(ix,iy,iz).zz() = pars.maxCR; } ;
545
546 if (cbdi(ix,iy,iz).x() > maxCT ) { cbdi(ix,iy,iz).x() = maxCT; } ;
547 if (cbdi(ix,iy,iz).y() > maxCT ) { cbdi(ix,iy,iz).y() = maxCT; } ;
548 if (cbdi(ix,iy,iz).z() > maxCT ) { cbdi(ix,iy,iz).z() = maxCT; } ;
549
550 surf(ix,iy,iz) /= cell_vol;
551
552 /* Calculate length scale of obstacles in each cell
553 Result is stored in surf. */
554
555 {
556 const scalar vb = v_block(ix,iy,iz);
557
558 if
559 (
560 (
561 ((area_block_s(ix,iy,iz).x() + area_block_r(ix,iy,iz).x()) < MIN_AB_FOR_SIZE)
562 && ((area_block_s(ix,iy,iz).y() + area_block_r(ix,iy,iz).y()) < MIN_AB_FOR_SIZE)
563 && ((area_block_s(ix,iy,iz).z() + area_block_r(ix,iy,iz).z()) < MIN_AB_FOR_SIZE)
564 )
565 || ( vb > MAX_VB_FOR_SIZE )
566 || ((obs_count(ix,iy,iz) + cmptSum(grating_count(ix,iy,iz))) < MIN_COUNT_FOR_SIZE)
567 || ( surf(ix,iy,iz) <= 0.0 )
568 )
569 {
570 obs_size(ix,iy,iz) = cell_size * pars.empty_lobs_fac;
571 }
572 else
573 {
574 /* A small sliver of a large cylinder ina cell can give large surface area
575 but low volume, hence snall "size". Therefore the vol/area formulation
576 is only fully implemented when count is at least COUNT_FOR_SIZE.*/
577 double nn, lobs, lobsMax;
578 nn = obs_count(ix,iy,iz) - sub_count(ix,iy,iz).x() + grating_count(ix,iy,iz).x();
579 if ( nn < 1.0 ) { nn = 1.0; }
580 lobsMax = (area_block_s(ix,iy,iz).x() + area_block_r(ix,iy,iz).x()) / nn * std::sqrt( dy * dz );
581 nn = obs_count(ix,iy,iz) - sub_count(ix,iy,iz).y() + grating_count(ix,iy,iz).y();
582 if ( nn < 1.0 ) { nn = 1.0; }
583 lobs = (area_block_s(ix,iy,iz).y() + area_block_r(ix,iy,iz).y()) / nn * std::sqrt( dz * dx );
584 if ( lobs > lobsMax )
585 {
586 lobsMax = lobs;
587 }
588
589 nn = obs_count(ix,iy,iz) - sub_count(ix,iy,iz).z() + grating_count(ix,iy,iz).z();
590 if ( nn < 1.0 ) { nn = 1.0; }
591 lobs = (area_block_s(ix,iy,iz).z() + area_block_r(ix,iy,iz).z()) / nn * std::sqrt( dx * dy );
592 if ( lobs > lobsMax )
593 {
594 lobsMax = lobs;
595 }
596
597 obs_size(ix,iy,iz) = lobsMax;
598 }
599 }
600
601 /* The formulation correctly deals with triple intersections. For quadruple intersections
602 and worse, there are very many second level overlaps and the resulting volume can be large
603 positive. However, many or all of these may be eliminated because of the minimum volume of
604 overlap blocks. Then the result can be negative volume - constrain accordingly
605 */
606
607 if (v_block(ix,iy,iz) < 0)
608 {
609 v_block(ix,iy,iz) = 0;
610 }
611 else if (v_block(ix,iy,iz) > 1)
612 {
613 v_block(ix,iy,iz) = 1;
614 }
615
616 /* We can get -ve sub_count (ns) if two pipes/bars intersect and the dominat direction
617 of the (-ve) intersection block is not the same as either of the intersecting obstacles.
618 Also, if we have two hirizontal abrs intersecting, the overlap block can have vertical
619 edges in a cell where the original bars do not. This can give -ve N and ns.
620 Negative N is removed by write_scalar. */
621
622 for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
623 {
624 if (sub_count(ix,iy,iz)[cmpt] < 0)
625 {
626 sub_count(ix,iy,iz)[cmpt] = 0;
627 }
628 }
629
630 v_block(ix,iy,iz) = 1.0 - v_block(ix,iy,iz); // Now porosity
631 }
632 }
633 }
634
635
636//*** Now we start writing the fields *********//
637
638 /* v_block is now porosity
639 The maximum value does not override the default value placed in the external cells,
640 so pars.cong_max_betav can be set just below 1 to mark the congested-region cells
641 for use by the adaptive mesh refinement. */
642
643 IjkField<Vector<direction>> n_blocked_faces
644 (
645 faceDims,
647 );
648
649 write_blocked_face_list
650 (
651 arr.face_block, arr.face_patch,
652 obs_count, sub_count, n_blocked_faces,
653 meshIndexing, patches,
654 pars.blockedFacePar, casepath
655 );
656 write_blockedCellsSet
657 (
658 arr.v_block,
659 meshIndexing, pars.blockedCellPoros, n_blocked_faces,
660 casepath, "blockedCellsSet"
661 );
662
663 write_scalarField
664 (
665 "betav", arr.v_block, 1, {0, pars.cong_max_betav}, "zeroGradient",
666 meshIndexing, patches,
667 dimless, casepath
668 );
669
670 for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
671 {
672 for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
673 {
674 for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
675 {
676 const scalar cell_vol = pdrBlock.V(ix, iy, iz);
677
678 /* After the correction to set the number of obstacles normal to a blocked face
679 to be zero, we can have N and all the components of ns the same. Then there
680 are no obstacles in the cell as the number in each direction is n minus ns component),
681 but N is not zero. This can cause problems. We reduce all four numbers by the same amount,
682 which is OK as only the difference is used except when N is checked to se if there are
683 any obstacles in then cell. */
684
685 scalar nmin = cmptMin(sub_count(ix,iy,iz));
686
687 sub_count(ix,iy,iz).x() -= nmin;
688 sub_count(ix,iy,iz).y() -= nmin;
689 sub_count(ix,iy,iz).z() -= nmin;
690
691 obs_count(ix,iy,iz) -= nmin;
692
693 assert(obs_count(ix,iy,iz) > -1);
694 if ( pars.new_fields )
695 {
696 /* New fields Nv and nsv are intensive quantities that stay unchanged as a cell is subdivided
697 We do not divide by cell volume because we assume that typical obstacle
698 is a cylinder passing through the cell */
699 const scalar cell_23 = ::pow(cell_vol, 2.0/3.0);
700 obs_count(ix,iy,iz) /= cell_23;
701 sub_count(ix,iy,iz) /= cell_23;
702 }
703 }
704 }
705 }
706
707
708 {
709 Info<< "Writing field files" << nl;
710
711 // obs_size is now the integral scale of the generated turbulence
712 write_scalarField
713 (
714 "Lobs", arr.obs_size, DEFAULT_LOBS, {0, 10}, "zeroGradient",
715 meshIndexing, patches,
716 dimLength, casepath
717 );
718 // surf is now surface area per unit volume
719 write_scalarField
720 (
721 "Aw", arr.surf, 0, {0, 1000}, "zeroGradient",
722 meshIndexing, patches,
723 inv(dimLength), casepath
724 );
725 write_symmTensorField
726 (
727 "CR", arr.drag_s, Zero, "zeroGradient",
728 meshIndexing, patches, inv(dimLength), casepath
729 );
730 write_symmTensorFieldV
731 (
732 "CT", cbdi, Zero, "zeroGradient",
733 meshIndexing, patches,
734 inv(dimLength), casepath
735 );
736 if ( pars.new_fields )
737 {
738 // These have been divided by cell volume ^ (2/3)
739 write_scalarField
740 (
741 "Nv", arr.obs_count, 0, {0, 1000}, "zeroGradient",
742 meshIndexing, patches,
743 dimless, casepath
744 );
745 write_symmTensorFieldV
746 (
747 "nsv", arr.sub_count, Zero, "zeroGradient",
748 meshIndexing, patches,
749 dimless, casepath
750 );
751 }
752 else
753 {
754 write_scalarField
755 (
756 "N", arr.obs_count, 0, {0, 1000}, "zeroGradient",
757 meshIndexing, patches,
758 dimless, casepath
759 );
760 write_symmTensorFieldV
761 (
762 "ns", arr.sub_count, Zero, "zeroGradient",
763 meshIndexing, patches, dimless, casepath
764 );
765 }
766
767 // Compute some further variables; store in already used arrays
768 // Re-use the drag array
769 drag_s = Zero;
770
771 for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
772 {
773 for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
774 {
775 for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
776 {
777 // Effective blockage ratio per cell/direction
778 vector eff_block =
779 (
780 area_block_s(ix,iy,iz) * pars.cd_s/pars.cd_r
781 + area_block_r(ix,iy,iz)
782 );
783
784 // Convert from B to Bv
785 if (pars.new_fields)
786 {
787 eff_block /= pdrBlock.width(ix, iy, iz);
788 }
789
790 // Effective blockage is zero when faces are blocked
791 for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
792 {
793 if (dirn_block(ix,iy,iz)[cmpt] || eff_block[cmpt] < 0)
794 {
795 eff_block[cmpt] = 0;
796 }
797 }
798
799 // Use the drag array to store the total effective blockage ratio per cell/direction
800 // - off-diagonal already zeroed
801 drag_s(ix,iy,iz).xx() = eff_block.x();
802 drag_s(ix,iy,iz).yy() = eff_block.y();
803 drag_s(ix,iy,iz).zz() = eff_block.z();
804
805 cbdi(ix,iy,iz).x() = 1.0 / (betai_inv1(ix,iy,iz).x() + 1.0);
806 cbdi(ix,iy,iz).y() = 1.0 / (betai_inv1(ix,iy,iz).y() + 1.0);
807 cbdi(ix,iy,iz).z() = 1.0 / (betai_inv1(ix,iy,iz).z() + 1.0);
808
809 if (cbdi(ix,iy,iz).z() < 0 || cbdi(ix,iy,iz).z() > 1.0)
810 {
812 << "beta_i problem. z-betai_inv1=" << betai_inv1(ix,iy,iz).z()
813 << " beta_i=" << cbdi(ix,iy,iz).z()
814 << nl;
815 }
816
817 //Use the obs_size array to store Ep
818 //We use Ep/(Xp-0.999) as length scale to avoid divide by zero,
819 // so this is OK for initial Xp=1.
820 obs_size(ix,iy,iz) = 0.001 / obs_size(ix,iy,iz);
821
822 // Use the count array to store the combustion flag ( --1 everywhere in rectangular cells).
823 obs_count(ix,iy,iz) = 1.0;
824 }
825 }
826 }
827
828 // drag array holds area blockage
829 if ( pars.new_fields )
830 {
831 write_symmTensorField
832 (
833 "Bv", arr.drag_s, Zero, "zeroGradient",
834 meshIndexing, patches,
835 dimless, casepath
836 );
837 }
838 else
839 {
840 write_symmTensorField
841 (
842 "B", arr.drag_s, Zero, "zeroGradient",
843 meshIndexing, patches,
844 dimless, casepath
845 );
846 }
847
848 // cbdi array holds beta_i
849 write_symmTensorFieldV
850 (
851 "betai", cbdi, symmTensor::I, "zeroGradient",
852 meshIndexing, patches,
853 dimless, casepath
854 );
855
856 // The longitudinal blockage
857 write_symmTensorFieldV
858 (
859 "Blong", arr.along_block, Zero, "zeroGradient",
860 meshIndexing, patches,
861 dimless, casepath
862 );
863
864 // obs_size array now contains 1/Lobs
865 write_scalarField
866 (
867 "Ep", arr.obs_size, DEFAULT_EP, {0, 10}, "zeroGradient",
868 meshIndexing, patches,
869 inv(dimLength), casepath
870 );
871 write_uniformField
872 (
873 "b", 1.0, "zeroGradient",
874 meshIndexing, patches,
875 dimless, casepath
876 );
877 write_uniformField
878 (
879 "k", DEFAULT_K, K_WALL_FN,
880 meshIndexing, patches,
882 casepath
883 );
884
885 write_uniformField
886 (
887 "epsilon", DEFAULT_EPS, EPS_WALL_FN,
888 meshIndexing, patches,
889 sqr(dimVelocity)/dimTime, casepath
890 );
891 write_uniformField
892 (
893 "ft", 0, "zeroGradient",
894 meshIndexing, patches,
895 dimless, casepath
896 );
897 write_uniformField
898 (
899 "Su", DEFAULT_SU, "zeroGradient",
900 meshIndexing, patches,
901 dimVelocity, casepath
902 );
903 write_uniformField
904 (
905 "T", DEFAULT_T, "zeroGradient",
906 meshIndexing, patches,
907 dimTemperature, casepath
908 );
909 write_uniformField
910 (
911 "Tu", DEFAULT_T, "zeroGradient",
912 meshIndexing, patches,
913 dimTemperature, casepath
914 );
915 write_uniformField
916 (
917 "Xi", 1, "zeroGradient",
918 meshIndexing, patches,
919 dimless, casepath
920 );
921 write_uniformField
922 (
923 "Xp", 1, "zeroGradient",
924 meshIndexing, patches,
925 dimless, casepath
926 );
927 write_uniformField
928 (
929 "GRxp", 0, "zeroGradient",
930 meshIndexing, patches,
931 inv(dimTime), casepath
932 );
933 write_uniformField
934 (
935 "GRep", 0, "zeroGradient",
936 meshIndexing, patches,
937 inv(dimLength*dimTime), casepath
938 );
939 write_uniformField
940 (
941 "RPers", 0, "zeroGradient",
942 meshIndexing, patches,
943 inv(dimTime), casepath
944 );
945 write_pU_fields(meshIndexing, patches, casepath);
946
947 write_uniformField
948 (
949 "alphat", 0, ALPHAT_WALL,
950 meshIndexing, patches,
952 casepath
953 );
954 write_uniformField
955 (
956 "nut", 0, NUT_WALL_FN,
957 meshIndexing, patches,
958 dimViscosity, casepath
959 );
960 // combustFlag is 1 in rectangular region, 0 or 1 elsewhere
961 // (although user could set it to another value)
962 if (equal(pars.outerCombFac, 1))
963 {
964 write_uniformField
965 (
966 "combustFlag", 1, "zeroGradient",
967 meshIndexing, patches,
968 dimless, casepath
969 );
970 }
971 else
972 {
973 write_scalarField
974 (
975 "combustFlag", arr.obs_count, pars.outerCombFac, {0, 1}, "zeroGradient",
976 meshIndexing, patches,
977 dimless, casepath
978 );
979 }
980 if ( pars.deluge )
981 {
982 write_uniformField
983 (
984 "H2OPS", 0, "zeroGradient",
985 meshIndexing, patches,
986 dimless, casepath
987 );
988 write_uniformField
989 (
990 "AIR", 0, "zeroGradient",
991 meshIndexing, patches,
992 dimless, casepath
993 );
994 write_uniformField
995 (
996 "Ydefault", 0, "zeroGradient",
997 meshIndexing, patches,
998 dimless, casepath
999 );
1000 write_uniformField
1001 (
1002 "eRatio", 1, "zeroGradient",
1003 meshIndexing, patches,
1004 dimless, casepath
1005 );
1006 write_uniformField
1007 (
1008 "sprayFlag", 1, "zeroGradient",
1009 meshIndexing, patches,
1010 dimless, casepath
1011 );
1012 }
1013 }
1014}
1015
1016
1018(
1019 const fileName& casepath,
1020 const PDRmeshArrays& meshIndexing,
1022)
1023{
1024 calculateAndWrite(*this, meshIndexing, casepath, patches);
1025}
1026
1027
1028void calc_drag_etc
1029(
1030 double brs, double brr, bool blocked,
1031 double surr_br, double surr_dr,
1032 scalar* drags_p, scalar* dragr_p,
1033 double count,
1034 scalar* cbdi_p,
1035 double cell_vol
1036)
1037{
1038 // Total blockage ratio
1039 scalar br = brr + brs;
1040
1041 // Idealise obstacle arrangement as sqrt(count) rows.
1042 // Make br the blockage ratio for each row.
1043 if (count > 1.0) { br /= std::sqrt(count); }
1044
1045 const scalar alpha =
1046 (
1047 br < 0.99
1048 ? (1.0 - 0.5 * br) / (1.0 - br) / (1.0 - br)
1049 : GREAT
1050 );
1051
1052 // For the moment keep separate the two contributions to the blockage-corrected drag
1053 /* An isolated long obstcale will have two of the surronding eight cells with the same blockage,
1054 so surr_br would be br/4. In this case no correction. Rising to full correction when
1055 all surrounding cells have the same blockage. */
1056 const scalar expon =
1057 (
1058 br > 0.0
1059 ? min(max((surr_br / br - 0.25) * 4.0 / 3.0, scalar(0)), scalar(1))
1060 : 0.0
1061 );
1062
1063 const scalar alpha_r = ::pow(alpha, 0.5 + 0.5 * expon);
1064 const scalar alpha_s = ::pow(alpha, expon);
1065
1066 *dragr_p *= alpha_r;
1067 *drags_p *= ::pow(alpha_s, 1.09);
1068 *cbdi_p = ( pars.cb_r * pars.cd_r * *dragr_p + pars.cb_s * pars.cd_s * *drags_p );
1069 if ( *cbdi_p < 0.0 ) { *cbdi_p = 0.0; }
1070
1071 // Finally sum the drag.
1072 *drags_p = ( *drags_p * pars.cd_s + *dragr_p * pars.cd_r );
1073 if ( *drags_p < 0.0 ) { *drags_p = 0.0; }
1074 /* If well-blocked cells are surrounded by empty cells, the flow just goes round
1075 and the drag parameters have little effect. So, for any cells much more empty
1076 than the surrounding cells, we put some CR in there as well. */
1077 if ( (surr_dr * 0.25) > *drags_p )
1078 {
1079 *drags_p = surr_dr * 0.25;
1080 *cbdi_p = *drags_p * (pars.cb_r + pars.cb_s ) * 0.5;
1081 // Don't know whether surr. stuff was round or sharp; use average of cb factors
1082 }
1083 if ( blocked ) { *cbdi_p = 0.0; *drags_p = 0.0; *dragr_p = 0.0; }
1084}
1085
1086
1088{
1089 if (isNull(block()))
1090 {
1092 << nl
1093 << "No blockage information - PDRblock is not set" << nl;
1094 return;
1095 }
1096
1097 const PDRblock& pdrBlock = block();
1098
1099 scalar totArea = 0;
1100 scalar totCount = 0;
1101 scalar totVolBlock = 0;
1102
1103 vector totBlock(Zero);
1104 vector totDrag(Zero);
1105
1106 for (label iz = 0; iz < pdrBlock.size(vector::Z); ++iz)
1107 {
1108 for (label iy = 0; iy < pdrBlock.size(vector::Y); ++iy)
1109 {
1110 for (label ix = 0; ix < pdrBlock.size(vector::X); ++ix)
1111 {
1112 const labelVector ijk(ix,iy,iz);
1113
1114 totVolBlock += v_block(ijk) * pdrBlock.V(ijk);
1115 totArea += surf(ijk);
1116
1117 totCount += max(0, obs_count(ijk));
1118
1119 totDrag.x() += max(0, drag_s(ijk).xx());
1120 totDrag.y() += max(0, drag_s(ijk).yy());
1121 totDrag.z() += max(0, drag_s(ijk).zz());
1122
1123 for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
1124 {
1125 totBlock[cmpt] += max(0, area_block_s(ijk)[cmpt]);
1126 totBlock[cmpt] += max(0, area_block_r(ijk)[cmpt]);
1127 }
1128 }
1129 }
1130 }
1131
1132 Info<< nl
1133 << "Volume blockage: " << totVolBlock << nl
1134 << "Total drag: " << totDrag << nl
1135 << "Total count: " << totCount << nl
1136 << "Total area blockage: " << totBlock << nl
1137 << "Total surface area: " << totArea << nl;
1138}
1139
1140
1141// ------------------------------------------------------------------------- //
1142
1143// Another temporary measure
1144template<class Type>
1145static void tail_field
1146(
1147 Ostream& os,
1148 const Type& deflt,
1149 const char* wall_bc,
1151)
1152{
1153 // ground
1154 {
1156 os.writeKeyword("type") << wall_bc << token::END_STATEMENT << nl;
1157 putUniform(os, "value", deflt);
1158 os.endBlock();
1159 }
1160
1161 forAll(patches, patchi)
1162 {
1163 const word& patchName = patches[patchi].patchName;
1164
1165 if (PDRpatchDef::BLOCKED_FACE == patchi)
1166 {
1167 // blockedFaces
1168 os.beginBlock(patchName);
1169
1170 // No wall functions for blockedFaces patch unless selected
1172 {
1173 os.writeKeyword("type") << wall_bc << token::END_STATEMENT << nl;
1174 putUniform(os, "value", deflt);
1175 }
1176 else
1177 {
1178 os.writeEntry("type", "zeroGradient");
1179 }
1180
1181 os.endBlock();
1182 }
1183 else if (patches[patchi].patchType == 0)
1184 {
1185 os.beginBlock(patchName);
1186
1187 os.writeKeyword("type") << wall_bc << token::END_STATEMENT << nl;
1188 putUniform(os, "value", deflt);
1189
1190 os.endBlock();
1191 }
1192 else
1193 {
1194 os.beginBlock(word(patchName + "Wall"));
1195 os.writeKeyword("type") << wall_bc << token::END_STATEMENT << nl;
1196 putUniform(os, "value", deflt);
1197 os.endBlock();
1198
1199 os.beginBlock(word(patchName + "Cyclic_half0"));
1200 os.writeEntry("type", "cyclic");
1201 os.endBlock();
1202
1203 os.beginBlock(word(patchName + "Cyclic_half1"));
1204 os.writeEntry("type", "cyclic");
1205 os.endBlock();
1206 }
1207 }
1208
1209 if (pars.yCyclic)
1210 {
1211 os.beginBlock("Cyclic_half0");
1212 os.writeEntry("type", "cyclic");
1213 os.endBlock();
1214
1215 os.beginBlock("Cyclic_half1");
1216 os.writeEntry("type", "cyclic");
1217 os.endBlock();
1218 }
1219 else
1220 {
1221 os.beginBlock("ySymmetry");
1222 os.writeEntry("type", "symmetryPlane");
1223 os.endBlock();
1224 }
1225
1226 if (pars.two_d)
1227 {
1228 os.beginBlock("z_boundaries");
1229 os.writeEntry("type", "empty");
1230 os.endBlock();
1231 }
1232
1233 if (pars.outer_orthog)
1234 {
1235 os.beginBlock("outer_inner");
1236 os.writeEntry("type", "cyclicAMI");
1237 os.writeEntry("neighbourPatch", "inner_outer");
1238 os.endBlock();
1239
1240 os.beginBlock("inner_outer");
1241 os.writeEntry("type", "cyclicAMI");
1242 os.writeEntry("neighbourPatch", "outer_inner");
1243 os.endBlock();
1244 }
1245}
1246
1247
1248// ------------------------------------------------------------------------- //
1249
1250void write_scalarField
1251(
1252 const word& fieldName, const IjkField<scalar>& fld,
1253 const scalar& deflt, const scalarMinMax& limits, const char *wall_bc,
1254 const PDRmeshArrays& meshIndexing,
1256 const dimensionSet& dims, const fileName& casepath
1257)
1258{
1259 fileName path = (casepath / pars.timeName / fieldName);
1260 OFstream os(path);
1261 os.precision(outputPrecision);
1262
1263 make_header(os, "", volScalarField::typeName, fieldName);
1264
1265 os.writeEntry("dimensions", dims);
1266
1267 os << nl;
1268 os.writeKeyword("internalField")
1269 << "nonuniform List<scalar>" << nl
1270 << meshIndexing.nCells() << nl << token::BEGIN_LIST << nl;
1271
1272 for (label celli=0; celli < meshIndexing.nCells(); ++celli)
1273 {
1274 const labelVector& cellIdx = meshIndexing.cellIndex[celli];
1275
1276 if (!isGoodIndex(cellIdx))
1277 {
1278 os << deflt << nl;
1279 continue;
1280 }
1281
1282 os << limits.clip(fld(cellIdx)) << nl;
1283 }
1284
1286
1287 os << nl;
1288 os.beginBlock("boundaryField");
1289
1290
1291 // outer
1292 {
1294
1295 os.writeEntry("type", "inletOutlet");
1296 putUniform(os, "inletValue", deflt);
1297 putUniform(os, "value", deflt);
1298
1299 os.endBlock();
1300 }
1301
1302 tail_field(os, deflt, wall_bc, patches);
1303
1304 os.endBlock(); // boundaryField
1305
1307}
1308
1309
1310// ------------------------------------------------------------------------- //
1311
1312void write_uniformField
1313(
1314 const word& fieldName, const scalar& deflt, const char *wall_bc,
1315 const PDRmeshArrays& meshIndexing,
1317 const dimensionSet& dims, const fileName& casepath
1318)
1319{
1320 OFstream os(casepath / pars.timeName / fieldName);
1321 os.precision(outputPrecision);
1322
1323 make_header(os, "", volScalarField::typeName, fieldName);
1324
1325 os.writeEntry("dimensions", dims);
1326
1327 os << nl;
1328 putUniform(os, "internalField", deflt);
1329
1330 os << nl;
1331 os.beginBlock("boundaryField");
1332
1333 // outer
1334 {
1336
1337 if (fieldName == "alphat" || fieldName == "nut")
1338 {
1339 // Different b.c. for alphat & nut
1340 os.writeEntry("type", "calculated");
1341 }
1342 else
1343 {
1344 os.writeEntry("type", "inletOutlet");
1345 putUniform(os, "inletValue", deflt);
1346 }
1347
1348 putUniform(os, "value", deflt);
1349 os.endBlock();
1350 }
1351
1352 tail_field(os, deflt, wall_bc, patches);
1353
1354 os.endBlock(); // boundaryField
1355
1357}
1358
1359
1360// ------------------------------------------------------------------------- //
1361
1362void write_pU_fields
1363(
1364 const PDRmeshArrays& meshIndexing,
1366 const fileName& casepath
1367)
1368{
1369 // Velocity field
1370 {
1371 OFstream os(casepath / pars.timeName / "U");
1372 os.precision(outputPrecision);
1373
1374 make_header(os, "", volVectorField::typeName, "U");
1375
1376 os.writeEntry("dimensions", dimVelocity);
1377
1378 os << nl;
1379 putUniform(os, "internalField", vector::zero);
1380
1381 os << nl;
1382 os.beginBlock("boundaryField");
1383
1384 // outer
1385 {
1387 os.writeEntry("type", "inletOutlet");
1388 putUniform(os, "inletValue", vector::zero);
1389 os.endBlock();
1390 }
1391
1392 // ground
1393 {
1395 os.writeEntry("type", "zeroGradient");
1396 os.endBlock();
1397 }
1398
1399 // Patch 0 is the blocked faces' and 1 is mergingFaces for ignition cell
1400 for (label patchi = 0; patchi < 3; ++patchi)
1401 {
1402 os.beginBlock(patches[patchi].patchName);
1403 os.writeKeyword("type") << pars.UPatchBc.c_str()
1405 os.endBlock();
1406 }
1407
1408 for (label patchi = 3; patchi < patches.size(); ++patchi)
1409 {
1410 const PDRpatchDef& p = patches[patchi];
1411 const word& patchName = p.patchName;
1412
1413 if (p.patchType == 0)
1414 {
1415 os.beginBlock(patchName);
1416
1417 os.writeEntry("type", "timeVaryingMappedFixedValue");
1418 os.writeEntry("fileName", "<case>" / (patchName + ".dat"));
1419 os.writeEntry("outOfBounds", "clamp");
1420 putUniform(os, "value", vector::zero);
1421 os.endBlock();
1422 }
1423 else
1424 {
1425 os.beginBlock(word(patchName + "Wall"));
1426 os.writeEntry("type", "activePressureForceBaffleVelocity");
1427
1428 os.writeEntry("cyclicPatch", word(patchName + "Cyclic_half0"));
1429 os.writeEntry("openFraction", 0); // closed
1430 os.writeEntry("openingTime", p.blowoffTime);
1431 os.writeEntry("minThresholdValue", p.blowoffPress);
1432 os.writeEntry("maxOpenFractionDelta", 0.1);
1433 os.writeEntry("forceBased", "false");
1434 os.writeEntry("opening", "true");
1435
1436 putUniform(os, "value", vector::zero);
1437 os.endBlock();
1438
1439 os.beginBlock(word(patchName + "Cyclic_half0"));
1440 os.writeEntry("type", "cyclic");
1441 putUniform(os, "value", vector::zero);
1442 os.endBlock();
1443
1444 os.beginBlock(word(patchName + "Cyclic_half1"));
1445 os.writeEntry("type", "cyclic");
1446 putUniform(os, "value", vector::zero);
1447 os.endBlock();
1448 }
1449 }
1450
1451 if (pars.yCyclic)
1452 {
1453 os.beginBlock("yCyclic_half0");
1454 os.writeEntry("type", "cyclic");
1455 os.endBlock();
1456
1457 os.beginBlock("yCyclic_half1");
1458 os.writeEntry("type", "cyclic");
1459 os.endBlock();
1460 }
1461 else
1462 {
1463 os.beginBlock("ySymmetry");
1464 os.writeEntry("type", "symmetryPlane");
1465 os.endBlock();
1466 }
1467
1468 if ( pars.outer_orthog )
1469 {
1470 os.beginBlock("outer_inner");
1471 os.writeEntry("type", "cyclicAMI");
1472 os.writeEntry("neighbourPatch", "inner_outer");
1473 os.endBlock();
1474
1475 os.beginBlock("inner_outer");
1476 os.writeEntry("type", "cyclicAMI");
1477 os.writeEntry("neighbourPatch", "outer_inner");
1478 }
1479
1480 os.endBlock(); // boundaryField
1481
1483 }
1484
1485
1486 // Pressure field
1487 {
1488 const scalar deflt = DEFAULT_P;
1489 const char *wall_bc = "zeroGradient;\n\trho\trho";
1490
1491 OFstream os(casepath / pars.timeName / "p");
1492 os.precision(outputPrecision);
1493
1494 make_header(os, "", volScalarField::typeName, "p");
1495
1496 os.writeEntry("dimensions", dimPressure);
1497
1498 os << nl;
1499 putUniform(os, "internalField", deflt);
1500
1501 os << nl;
1502 os.beginBlock("boundaryField");
1503
1504 // outer
1505 {
1507
1508 os.writeEntry("type", "waveTransmissive");
1509 os.writeEntry("gamma", 1.3);
1510 os.writeEntry("fieldInf", deflt);
1511 os.writeEntry("lInf", 5);
1512 putUniform(os, "value", deflt);
1513 os.endBlock();
1514 }
1515
1516 tail_field(os, deflt, wall_bc, patches);
1517
1518 os.endBlock(); // boundaryField
1519
1521 }
1522}
1523
1524
1525// ------------------------------------------------------------------------- //
1526
1527void write_symmTensorField
1528(
1529 const word& fieldName,
1531 const symmTensor& deflt, const char *wall_bc,
1532 const PDRmeshArrays& meshIndexing,
1534 const dimensionSet& dims, const fileName& casepath
1535)
1536{
1537 OFstream os(casepath / pars.timeName / fieldName);
1538 os.precision(outputPrecision);
1539
1540 make_header(os, "", volSymmTensorField::typeName, fieldName);
1541
1542 os.writeEntry("dimensions", dims);
1543
1544 os << nl;
1545 os.writeKeyword("internalField")
1546 << "nonuniform List<symmTensor>" << nl
1547 << meshIndexing.nCells() << nl << token::BEGIN_LIST << nl;
1548
1549 for (label celli=0; celli < meshIndexing.nCells(); ++celli)
1550 {
1551 const labelVector& cellIdx = meshIndexing.cellIndex[celli];
1552
1553 if (!isGoodIndex(cellIdx))
1554 {
1555 os << deflt << nl;
1556 continue;
1557 }
1558
1559 os << fld(cellIdx) << nl;
1560 }
1562
1563 os << nl;
1564 os.beginBlock("boundaryField");
1565
1566 // outer
1567 {
1569
1570 os.writeEntry("type", "inletOutlet");
1571 putUniform(os, "inletValue", deflt);
1572 putUniform(os, "value", deflt);
1573
1574 os.endBlock();
1575 }
1576
1577 tail_field(os, deflt, wall_bc, patches);
1578
1579 os.endBlock(); // boundaryField
1580
1582}
1583
1584
1585// Write a volSymmTensorField but with vectors as input.
1586// The off-diagonals are zero.
1587void write_symmTensorFieldV
1588(
1589 const word& fieldName,
1590 const IjkField<vector>& fld,
1591 const symmTensor& deflt, const char *wall_bc,
1592 const PDRmeshArrays& meshIndexing,
1594 const dimensionSet& dims, const fileName& casepath
1595)
1596{
1597 OFstream os(casepath / pars.timeName / fieldName);
1598 os.precision(outputPrecision);
1599
1600 make_header(os, "", volSymmTensorField::typeName, fieldName);
1601
1602 os.writeEntry("dimensions", dims);
1603
1604 os << nl;
1605 os.writeKeyword("internalField")
1606 << "nonuniform List<symmTensor>" << nl
1607 << meshIndexing.nCells() << nl << token::BEGIN_LIST << nl;
1608
1610
1611 for (label celli=0; celli < meshIndexing.nCells(); ++celli)
1612 {
1613 const labelVector& cellIdx = meshIndexing.cellIndex[celli];
1614
1615 if (!isGoodIndex(cellIdx))
1616 {
1617 os << deflt << nl;
1618 continue;
1619 }
1620
1621 const vector& vec = fld(cellIdx);
1622
1623 val.xx() = vec.x();
1624 val.yy() = vec.y();
1625 val.zz() = vec.z();
1626
1627 os << val << nl;
1628 }
1630
1631 os << nl;
1632 os.beginBlock("boundaryField");
1633
1634 // outer
1635 {
1637
1638 os.writeEntry("type", "inletOutlet");
1639 putUniform(os, "inletValue", deflt);
1640 putUniform(os, "value", deflt);
1641
1642 os.endBlock();
1643 }
1644
1645 tail_field(os, deflt, wall_bc, patches);
1646
1647 os.endBlock(); // boundaryField
1648
1650}
1651
1652
1653// ------------------------------------------------------------------------- //
1654
1655void write_blocked_face_list
1656(
1657 const IjkField<vector>& face_block,
1658 const IjkField<labelVector>& face_patch,
1659 const IjkField<scalar>& obs_count, IjkField<vector>& sub_count,
1660 IjkField<Vector<direction>>& n_blocked_faces,
1661 const PDRmeshArrays& meshIndexing,
1663 double limit_par, const fileName& casepath
1664)
1665{
1666 /* Create the lists of face numbers for faces that have already been defined as
1667 belonging to (inlet) patches), and others that are found to be blocked.
1668 Then write these out to set files, */
1669
1670 const labelVector& cellDims = meshIndexing.cellDims;
1671
1672 Map<bitSet> usedFaces;
1673
1674 Info<< "Number of patches: " << patches.size() << nl;
1675
1676 for (label facei=0; facei < meshIndexing.nFaces(); ++facei)
1677 {
1678 // The related i-j-k face index for the mesh face
1679 const labelVector& faceIdx = meshIndexing.faceIndex[facei];
1680
1681 if (!isGoodIndex(faceIdx))
1682 {
1683 continue;
1684 }
1685
1686 const label ix = faceIdx.x();
1687 const label iy = faceIdx.y();
1688 const label iz = faceIdx.z();
1689 const direction orient = meshIndexing.faceOrient[facei];
1690
1691 label patchId = -1;
1692 scalar val(Zero);
1693
1694 /* A bit messy to be changing sub_count here. but there is a problem of generation
1695 of subgrid flame area Xp when the flame approaches a blocked wall. the fix is to make
1696 the normal component of "n" zero in the cells adjacent to the blocked face. That component
1697 of n is zero when that component of sub_count i.e. ns) equals count (i.e. N). */
1698 {
1699 switch (orient)
1700 {
1701 case vector::X:
1702 {
1703 // face_block is the face blockage;
1704 // face_patch is the patch number on the face (if any)
1705 val = face_block(faceIdx).x();
1706 patchId = face_patch(faceIdx).x();
1707
1708 if
1709 (
1710 val > limit_par
1711 && iy < cellDims[vector::Y]
1712 && iz < cellDims[vector::Z]
1713 )
1714 {
1715 // n_blocked_faces:
1716 // count of x-faces blocked for this cell
1717
1718 if (ix < cellDims[vector::X])
1719 {
1720 ++n_blocked_faces(ix,iy,iz).x();
1721 sub_count(ix,iy,iz).x() = obs_count(ix,iy,iz);
1722 }
1723
1724 if (ix > 0)
1725 {
1726 // And the neighbouring cell
1727 ++n_blocked_faces(ix-1,iy,iz).x();
1728 sub_count(ix-1,iy,iz).x() = obs_count(ix-1,iy,iz);
1729 }
1730 }
1731 }
1732 break;
1733
1734 case vector::Y:
1735 {
1736 val = face_block(faceIdx).y();
1737 patchId = face_patch(faceIdx).y();
1738
1739 if
1740 (
1741 val > limit_par
1742 && iz < cellDims[vector::Z]
1743 && ix < cellDims[vector::X]
1744 )
1745 {
1746 // n_blocked_faces:
1747 // count of y-faces blocked for this cell
1748
1749 if (iy < cellDims[vector::Y])
1750 {
1751 ++n_blocked_faces(ix,iy,iz).y();
1752 sub_count(ix,iy,iz).y() = obs_count(ix,iy,iz);
1753 }
1754
1755 if (iy > 0)
1756 {
1757 // And the neighbouring cell
1758 ++n_blocked_faces(ix,iy-1,iz).y();
1759 sub_count(ix,iy-1,iz).y() = obs_count(ix,iy-1,iz);
1760 }
1761 }
1762 }
1763 break;
1764
1765 case vector::Z:
1766 {
1767 val = face_block(faceIdx).z();
1768 patchId = face_patch(faceIdx).z();
1769
1770 if
1771 (
1772 val > limit_par
1773 && ix < cellDims[vector::X]
1774 && iy < cellDims[vector::Y]
1775 )
1776 {
1777 // n_blocked_faces:
1778 // count of z-faces blocked for this cell
1779
1780 if (iz < cellDims[vector::Z])
1781 {
1782 ++n_blocked_faces(ix,iy,iz).z();
1783 sub_count(ix,iy,iz).z() = obs_count(ix,iy,iz);
1784 }
1785
1786 if (iz > 0)
1787 {
1788 // And the neighbouring cell
1789 ++n_blocked_faces(ix,iy,iz-1).z();
1790 sub_count(ix,iy,iz-1).z() = obs_count(ix,iy,iz-1);
1791 }
1792 }
1793 }
1794 break;
1795 }
1796
1797 if (patchId > 0)
1798 {
1799 // If this face is on a defined patch add to list
1800 usedFaces(patchId).set(facei);
1801 }
1802 else if (val > limit_par)
1803 {
1804 // Add to blocked faces list
1805 usedFaces(PDRpatchDef::BLOCKED_FACE).set(facei);
1806 }
1807 }
1808 }
1809
1810 // Write in time or constant dir
1811 const bool hasPolyMeshTimeDir = isDir(casepath/pars.timeName/"polyMesh");
1812
1813 const fileName setsDir =
1814 (
1815 casepath
1816 / (hasPolyMeshTimeDir ? pars.timeName : word("constant"))
1817 / fileName("polyMesh/sets")
1818 );
1819
1820 if (!isDir(setsDir))
1821 {
1822 mkDir(setsDir);
1823 }
1824
1825
1826 // Create as blockedFaces Set file for each patch, including
1827 // basic blocked faces
1828 forAll(patches, patchi)
1829 {
1830 const word& patchName = patches[patchi].patchName;
1831
1832 OFstream os(setsDir / (patchName + "Set"));
1833
1834 make_header(os, "polyMesh/sets", "faceSet", patchName);
1835
1836 // Check for blocked faces
1837 const auto& fnd = usedFaces.cfind(patchi);
1838
1839 if (fnd.good() && (*fnd).any())
1840 {
1841 os << nl << (*fnd).toc() << nl;
1842 }
1843 else
1844 {
1845 os << nl << labelList() << nl;
1846 }
1847
1849 }
1850
1851 // Create the PDRMeshDict, listing the blocked faces sets and their patch names
1852
1853 {
1854 DynamicList<word> panelNames;
1855
1856 OFstream os(casepath / "system/PDRMeshDict");
1857
1858 make_header(os, "system", "dictionary", "PDRMeshDict");
1859
1860 os.writeEntry("blockedCells", "blockedCellsSet");
1861 os << nl << "blockedFaces" << nl << token::BEGIN_LIST << nl;
1862
1863 for (const PDRpatchDef& p : patches)
1864 {
1865 const word& patchName = p.patchName;
1866 const word setName = patchName + "Set";
1867
1868 if (p.patchType == 0) // Patch
1869 {
1870 os << " " << token::BEGIN_LIST
1871 << setName << token::SPACE
1872 << patchName << token::END_LIST
1873 << nl;
1874 }
1875 else if (p.patchType > 0) // Panel
1876 {
1877 panelNames.append(setName);
1878 }
1879 }
1880
1882 os.beginBlock("coupledFaces");
1883
1884 for (const PDRpatchDef& p : patches)
1885 {
1886 const word& patchName = p.patchName;
1887 const word setName = patchName + "Set";
1888
1889 if (p.patchType > 0) // Panel
1890 {
1891 os.beginBlock(setName);
1892 os.writeEntry("wallPatch", word(patchName + "Wall"));
1893 os.writeEntry("cyclicMasterPatch", word(patchName + "Cyclic_half0"));
1894 os.endBlock();
1895 }
1896 }
1897 os.endBlock() << nl;
1898
1899 os.writeEntry("defaultPatch", "blockedFaces");
1900
1902
1903 // Write panelList
1904 OFstream(casepath / "panelList")()
1905 << panelNames << token::END_STATEMENT << nl;
1906 }
1907}
1908
1909
1910void write_blockedCellsSet
1911(
1912 const IjkField<scalar>& fld,
1913 const PDRmeshArrays& meshIndexing,
1914 double limit_par,
1915 const IjkField<Vector<direction>>& n_blocked_faces,
1916 const fileName& casepath,
1917 const word& listName
1918)
1919{
1920 if (listName.empty())
1921 {
1922 return;
1923 }
1924
1925 // Write in time or constant dir
1926 const bool hasPolyMeshTimeDir = isDir(casepath/pars.timeName/"polyMesh");
1927
1928 const fileName path =
1929 (
1930 casepath
1931 / (hasPolyMeshTimeDir ? pars.timeName : word("constant"))
1932 / fileName("polyMesh/sets")
1933 / listName
1934 );
1935
1936 if (!isDir(path.path()))
1937 {
1938 mkDir(path.path());
1939 }
1940
1941 bitSet blockedCell;
1942
1943 for (label celli=0; celli < meshIndexing.nCells(); ++celli)
1944 {
1945 const labelVector& cellIdx = meshIndexing.cellIndex[celli];
1946
1947 if (!isGoodIndex(cellIdx))
1948 {
1949 continue;
1950 }
1951
1952 if (fld(cellIdx) < limit_par)
1953 {
1954 blockedCell.set(celli);
1955 continue;
1956 }
1957
1958 const Vector<direction>& blocked = n_blocked_faces(cellIdx);
1959
1960 const label n_bfaces = cmptSum(blocked);
1961
1962 label n_bpairs = 0;
1963
1964 if (n_bfaces > 1)
1965 {
1966 for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
1967 {
1968 if (blocked[cmpt] > 1) ++n_bpairs;
1969 }
1970
1971 #if 0
1972 // Extra debugging
1973 Info<<"block " << celli << " from "
1974 << blocked << " -> ("
1975 << n_bfaces << ' ' << n_bpairs
1976 << ')' << nl;
1977 #endif
1978 }
1979
1980 if
1981 (
1982 n_bfaces >= pars.nFacesToBlockC
1983 || n_bpairs >= pars.nPairsToBlockC
1984 )
1985 {
1986 blockedCell.set(celli);
1987 }
1988 }
1989
1990
1991 OFstream os(path);
1992 make_header(os, "constant/polyMesh/sets", "cellSet", listName);
1993
1994 if (blockedCell.any())
1995 {
1996 os << blockedCell.toc();
1997 }
1998 else
1999 {
2000 os << labelList();
2001 }
2002
2004
2006}
2007
2008
2009// ************************************************************************* //
scalar y
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
Preparation of fields for PDRFoam.
#define EPS_WALL_FN
Definition: PDRsetFields.H:70
#define DEFAULT_T
Definition: PDRsetFields.H:61
#define MIN_COUNT_FOR_SIZE
Definition: PDRsetFields.H:85
#define NUT_WALL_FN
Definition: PDRsetFields.H:73
#define DEFAULT_P
Definition: PDRsetFields.H:62
#define ALPHAT_WALL
Definition: PDRsetFields.H:71
#define MAX_VB_FOR_SIZE
Definition: PDRsetFields.H:83
#define DEFAULT_LOBS
Definition: PDRsetFields.H:64
#define DEFAULT_SU
Definition: PDRsetFields.H:63
#define K_WALL_FN
Definition: PDRsetFields.H:69
#define MIN_AB_FOR_SIZE
Definition: PDRsetFields.H:82
#define DEFAULT_K
Definition: PDRsetFields.H:59
#define DEFAULT_EPS
Definition: PDRsetFields.H:60
#define DEFAULT_EP
Definition: PDRsetFields.H:66
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
A HashTable similar to std::unordered_map.
Definition: HashTable.H:123
bool set(const Key &key, const T &obj)
Copy assign a new entry, overwriting existing entries.
Definition: HashTableI.H:202
const_iterator cfind(const Key &key) const
Find and return an const_iterator set at the hashed entry.
Definition: HashTableI.H:141
static Ostream & writeEndDivider(Ostream &os)
Write the standard end file divider.
Generic templated field type with i-j-k addressing.
Definition: IjkField.H:56
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
const T & clip(const T &val) const
Definition: MinMaxI.H:223
Output to file stream, using an OSstream.
Definition: OFstream.H:57
virtual int precision() const
Get precision of output field.
Definition: OSstream.C:326
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
virtual Ostream & endBlock()
Write end block group.
Definition: Ostream.C:105
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
Definition: Ostream.H:239
virtual Ostream & writeKeyword(const keyType &kw)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:57
virtual Ostream & beginBlock(const keyType &kw)
Write begin block group with the given name.
Definition: Ostream.C:87
Work array definitions for PDR fields.
Definition: PDRarrays.H:64
void blockageSummary() const
Summary of the blockages.
IjkField< vector > drag_r
Directional drag from round obstacles.
Definition: PDRarrays.H:116
IjkField< vector > betai_inv1
Definition: PDRarrays.H:99
IjkField< vector > along_block
Definition: PDRarrays.H:97
IjkField< scalar > v_block
Volume blockage.
Definition: PDRarrays.H:74
IjkField< vector > area_block_r
Summed area blockage (directional) from round obstacles.
Definition: PDRarrays.H:86
IjkField< vector > sub_count
Number of obstacles parallel to specified direction.
Definition: PDRarrays.H:106
IjkField< vector > area_block_s
Summed area blockage (directional) from sharp obstacles.
Definition: PDRarrays.H:83
const PDRblock & block() const
Reference to PDRblock.
Definition: PDRarrays.H:168
IjkField< scalar > obs_size
Obstacle size in cell.
Definition: PDRarrays.H:80
IjkField< Vector< bool > > dirn_block
A total directional blockage in the cell.
Definition: PDRarrays.H:89
IjkField< scalar > obs_count
Number of obstacles in cell.
Definition: PDRarrays.H:103
IjkField< vector > grating_count
Definition: PDRarrays.H:110
IjkField< symmTensor > drag_s
Tensorial drag from sharp obstacles.
Definition: PDRarrays.H:113
IjkField< vector > face_block
Definition: PDRarrays.H:93
IjkField< labelVector > face_patch
Face field for (directional) for patch Id.
Definition: PDRarrays.H:143
IjkField< scalar > surf
Surface area in cell.
Definition: PDRarrays.H:77
A single block x-y-z rectilinear mesh addressable as i,j,k with simplified creation....
Definition: PDRblock.H:156
scalar V(const label i, const label j, const label k) const
Cell volume at i,j,k position.
Definition: PDRblockI.H:266
scalar dx(const label i) const
Cell size in x-direction at i position.
Definition: PDRblockI.H:166
scalar dy(const label j) const
Cell size in y-direction at j position.
Definition: PDRblockI.H:178
scalar dz(const label k) const
Cell size in z-direction at k position.
Definition: PDRblockI.H:190
scalar width(const label i, const label j, const label k) const
Characteristic cell size at i,j,k position.
Definition: PDRblockI.H:283
OpenFOAM/PDRblock addressing information.
Definition: PDRmeshArrays.H:66
label nCells() const
The number of cells.
List< labelVector > cellIndex
For each cell, the corresponding i-j-k address.
Definition: PDRmeshArrays.H:79
List< direction > faceOrient
For each face, the x/y/z orientation.
Definition: PDRmeshArrays.H:85
labelVector cellDims
The cell i-j-k addressing range.
Definition: PDRmeshArrays.H:73
List< labelVector > faceIndex
For each face, the corresponding i-j-k address.
Definition: PDRmeshArrays.H:82
label nFaces() const
The number of faces.
labelVector faceDims
The face i-j-k addressing range.
Definition: PDRmeshArrays.H:76
string UPatchBc
"fixedValue;value uniform (0 0 0)"
Definition: PDRparams.H:71
scalar blockedFacePar
Faces with area blockage greater than this are blocked.
Definition: PDRparams.H:130
scalar empty_lobs_fac
Lobs in empty cell is this * cube root of cell volume.
Definition: PDRparams.H:116
scalar maxCR
Upper limit on CR (CT also gets limited)
Definition: PDRparams.H:133
bool outer_orthog
Definition: PDRparams.H:84
word groundPatchName
The name for the "ground" patch.
Definition: PDRparams.H:66
bool blockedFacesWallFn
Definition: PDRparams.H:82
scalar blockedCellPoros
Cells with porosity less than this are blocked.
Definition: PDRparams.H:127
int nPairsToBlockC
Definition: PDRparams.H:94
scalar outerCombFac
Value for outer region.
Definition: PDRparams.H:119
int nFacesToBlockC
Definition: PDRparams.H:90
word outerPatchName
The name for the "outer" patch.
Definition: PDRparams.H:69
Bookkeeping for patch definitions.
Definition: PDRpatchDef.H:54
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
Definition: SquareMatrix.H:66
static const SymmTensor I
Definition: SymmTensor.H:74
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
labelList toc() const
The indices of the on bits as a sorted labelList.
Definition: bitSet.C:474
bool any() const
True if any bits in this bitset are set.
Definition: bitSetI.H:469
Creates a single block of cells from point coordinates, numbers of cells in each direction and an exp...
Definition: block.H:61
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
Definition: dimensionSet.H:109
Particle-size distribution model wherein random samples are drawn from the doubly-truncated uniform p...
Definition: uniform.H:164
A class for handling file names.
Definition: fileName.H:76
static std::string path(const std::string &str)
Return directory path name (part before last /)
Definition: fileNameI.H:176
void calculateAndWrite()
Calculate the fields and write.
label size() const
Return the total i*j*k size.
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
@ END_STATEMENT
End entry [isseparator].
Definition: token.H:154
@ BEGIN_LIST
Begin list [isseparator].
Definition: token.H:155
@ END_LIST
End list [isseparator].
Definition: token.H:156
@ SPACE
Space [isspace].
Definition: token.H:125
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
volScalarField & p
const polyBoundaryMesh & patches
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
label patchId(-1)
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
const dimensionSet dimPressure
const dimensionSet dimViscosity
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
const dimensionSet dimless
Dimensionless.
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
bool mkDir(const fileName &pathName, mode_t mode=0777)
Make a directory and return an error if it could not be created.
Definition: MSwindows.C:515
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
const dimensionSet dimVelocity
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
bool isNull(const T *ptr)
True if ptr is a pointer (of type T) to the nullObject.
Definition: nullObject.H:192
messageStream Info
Information stream (stdout output on master, null elsewhere)
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
bool equal(const T &s1, const T &s2)
Compare two values for equality.
Definition: doubleFloat.H:46
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0)
Definition: dimensionSets.H:54
Foam::PDRparams pars
Globals for program parameters (ugly hack)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
uint8_t direction
Definition: direction.H:56
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
bool isDir(const fileName &name, const bool followLink=true)
Does the name exist as a DIRECTORY in the file system?
Definition: MSwindows.C:651
const double expon
Definition: doubleFloat.H:63
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:51
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & alpha
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
static const char *const typeName
The type name used in ensight case files.
A non-counting (dummy) refCount.
Definition: refCount.H:59