Line data Source code
1 : #ifndef AMREX_GEOMETRY_H_
2 : #define AMREX_GEOMETRY_H_
3 : #include <AMReX_Config.H>
4 :
5 : #include <AMReX_Array.H>
6 : #include <AMReX_CoordSys.H>
7 : #include <AMReX_ParallelDescriptor.H>
8 : #include <AMReX_RealBox.H>
9 : #include <AMReX_Periodicity.H>
10 :
11 : #ifdef AMREX_USE_OMP
12 : #include <omp.h>
13 : #endif
14 :
15 : #include <iosfwd>
16 : #include <map>
17 :
18 : namespace amrex {
19 :
20 : class MultiFab;
21 : class DistributionMapping;
22 : class BoxArray;
23 :
24 : struct GeometryData
25 : {
26 : //! Returns the cellsize for each coordinate direction.
27 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
28 : const Real* CellSize () const noexcept { return dx; }
29 : //
30 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
31 : Real CellSize (int dir) const noexcept { return dx[dir]; }
32 : //! Returns the lo end of the problem domain in each dimension.
33 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
34 : const Real* ProbLo () const noexcept { return prob_domain.lo(); }
35 : //
36 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
37 : Real ProbLo (int dir) const noexcept { return prob_domain.lo(dir); }
38 : //! Returns the hi end of the problem domain in each dimension.
39 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
40 : const Real* ProbHi () const noexcept { return prob_domain.hi(); }
41 : //
42 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
43 : Real ProbHi (int dir) const noexcept { return prob_domain.hi(dir); }
44 : //! Returns our rectangular domain.
45 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
46 : const Box& Domain () const noexcept { return domain; }
47 : //! Returns whether the domain is periodic in the given direction.
48 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
49 : int isPeriodic (const int i) const noexcept { return is_periodic[i]; }
50 : //! Coordinates type
51 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
52 : int Coord () const noexcept { return coord; }
53 :
54 : public:
55 : RealBox prob_domain;
56 : Box domain;
57 : Real dx[AMREX_SPACEDIM];
58 : int is_periodic[AMREX_SPACEDIM]; /**< For each dimension, 0 if the domain is non-periodic and 1 if it is. */
59 : int coord;
60 : };
61 :
62 : /**
63 : * \class Geometry
64 : * \brief Rectangular problem domain geometry.
65 : *
66 : * This class describes problem domain and coordinate system for
67 : * RECTANGULAR problem domains. Since the problem domain is RECTANGULAR,
68 : * periodicity is meaningful.
69 : */
70 : class Geometry
71 : :
72 : public CoordSys
73 : {
74 : public:
75 : /** \brief The default constructor.
76 : *
77 : * Leaves object in an unusable state. A "define" method must be called before use.
78 : */
79 : Geometry () noexcept;
80 :
81 : /**
82 : * Constructs a fully-defined geometry object that describes the problem domain and coordinate system.
83 : *
84 : * \param dom cell-centered `Box` specifying the index space bounds of the domain.
85 : * \param rb A pointer to a `RealBox` specifying the real space bounds of the domain. Optional. If not specified,
86 : * will be ParmParsed from "geometry.prob_lo" and "geometry.prob_hi".
87 : * \param coord Specifies the coordinate system type. Optional. Can be one of:
88 : * CoordSys::cartesian
89 : * CoordSys::RZ
90 : * CoordSys::SPHERICAL
91 : * Default is cartesian.
92 : * \param is_per pointer to memory space specifying periodicity in each coordinate direction.
93 : * Optional. Default is non-periodic.
94 : */
95 : explicit Geometry (const Box& dom,
96 : const RealBox* rb = nullptr,
97 : int coord = -1,
98 : int const* is_per = nullptr) noexcept;
99 : /**
100 : * Constructs a fully-defined geometry object that describes the problem domain and coordinate system.
101 : *
102 : * \param dom cell-centered `Box` specifying the index space bounds of the domain.
103 : * \param rb A pointer to a `RealBox` specifying the real space bounds of the domain. Optional. If not specified,
104 : * will be ParmParsed from "geometry.prob_lo" and "geometry.prob_hi".
105 : * \param coord Specifies the coordinate system type. Optional. Can be one of:
106 : * CoordSys::cartesian
107 : * CoordSys::RZ
108 : * CoordSys::SPHERICAL
109 : * Default is cartesian.
110 : * \param is_per amrex::Array specifying periodicity in each coordinate direction. Optional.
111 : * Default is non-periodic.
112 : */
113 : Geometry (const Box& dom, const RealBox& rb, int coord,
114 : Array<int,AMREX_SPACEDIM> const& is_per) noexcept;
115 :
116 : ~Geometry () = default;
117 : Geometry (const Geometry& rhs) = default;
118 : Geometry (Geometry&& rhs) noexcept = default;
119 : Geometry& operator= (const Geometry& rhs) = default;
120 : Geometry& operator= (Geometry&& rhs) noexcept = default;
121 :
122 : //! Returns non-static copy of geometry's stored data.
123 : [[nodiscard]] GeometryData data() const noexcept {
124 : return { prob_domain, domain, {AMREX_D_DECL(dx[0],dx[1],dx[2])},
125 : {AMREX_D_DECL(is_periodic[0], is_periodic[1], is_periodic[2])},
126 : static_cast<int>(c_sys) };
127 : }
128 :
129 : //! Read static values from ParmParse database.
130 : static void Setup (const RealBox* rb = nullptr, int coord = -1, int const* isper = nullptr) noexcept;
131 : static void ResetDefaultProbDomain (const RealBox& rb) noexcept;
132 : static void ResetDefaultPeriodicity (const Array<int,AMREX_SPACEDIM>& is_per) noexcept;
133 : static void ResetDefaultCoord (int coord) noexcept;
134 :
135 : /**
136 : * Defines a geometry object that describes the problem domain and coordinate system for rectangular problem
137 : * domains. Useful for when the Geometry object has been constructed with the default constructor.
138 : *
139 : * \param dom cell-centered `Box` specifying the index space bounds of the domain.
140 : * \param rb A pointer to a `RealBox` specifying the real space bounds of the domain. Optional. If not specified,
141 : * will be ParmParsed from "geometry.prob_lo" and "geometry.prob_hi".
142 : * \param coord Specifies the coordinate system type. Optional. Can be one of:
143 : * CoordSys::cartesian
144 : * CoordSys::RZ
145 : * CoordSys::SPHERICAL
146 : * Default is cartesian.
147 : * \param is_per pointer to memory space specifying periodicity in each coordinate direction.
148 : * Optional. Default is non-periodic.
149 : */
150 : void define (const Box& dom, const RealBox* rb = nullptr, int coord = -1, int const* is_per = nullptr) noexcept;
151 :
152 : /**
153 : * Defines a geometry object that describes the problem domain and coordinate system for rectangular problem
154 : * domains. Useful for when the Geometry object has been constructed with the default constructor.
155 : *
156 : * \param dom cell-centered `Box` specifying the index space bounds of the domain.
157 : * \param rb A pointer to a `RealBox` specifying the real space bounds of the domain. Optional. If not specified,
158 : * will be ParmParsed from "geometry.prob_lo" and "geometry.prob_hi".
159 : * \param coord Specifies the coordinate system type. Optional. Can be one of:
160 : * CoordSys::cartesian
161 : * CoordSys::RZ
162 : * CoordSys::SPHERICAL
163 : * Default is cartesian.
164 : * \param is_per amrex::Array specifying periodicity in each coordinate direction. Optional.
165 : * Default is non-periodic.
166 : */
167 : void define (const Box& dom, const RealBox& rb, int coord, Array<int,AMREX_SPACEDIM> const& is_per) noexcept;
168 :
169 : //! Returns the problem domain.
170 : [[nodiscard]] const RealBox& ProbDomain () const noexcept { return prob_domain; }
171 : //! Sets the problem domain.
172 : void ProbDomain (const RealBox& rb) noexcept
173 : {
174 : prob_domain = rb;
175 : computeRoundoffDomain();
176 : }
177 : //! Returns the lo end of the problem domain in each dimension.
178 1057510 : [[nodiscard]] const Real* ProbLo () const noexcept { return prob_domain.lo(); }
179 : //! Returns the hi end of the problem domain in each dimension.
180 100 : [[nodiscard]] const Real* ProbHi () const noexcept { return prob_domain.hi(); }
181 : //! Returns the lo end of the problem domain in specified direction.
182 80 : [[nodiscard]] Real ProbLo (int dir) const noexcept { return prob_domain.lo(dir); }
183 : //! Returns the hi end of the problem domain in specified direction.
184 40 : [[nodiscard]] Real ProbHi (int dir) const noexcept { return prob_domain.hi(dir); }
185 :
186 : [[nodiscard]] GpuArray<Real,AMREX_SPACEDIM> ProbLoArray () const noexcept {
187 : return {{AMREX_D_DECL(prob_domain.lo(0),prob_domain.lo(1),prob_domain.lo(2))}};
188 : }
189 :
190 : [[nodiscard]] GpuArray<Real,AMREX_SPACEDIM> ProbHiArray () const noexcept {
191 : return {{AMREX_D_DECL(prob_domain.hi(0),prob_domain.hi(1),prob_domain.hi(2))}};
192 : }
193 :
194 : [[nodiscard]] GpuArray<ParticleReal,AMREX_SPACEDIM> ProbLoArrayInParticleReal () const noexcept {
195 : return roundoff_lo;
196 : }
197 :
198 : [[nodiscard]] GpuArray<ParticleReal,AMREX_SPACEDIM> ProbHiArrayInParticleReal () const noexcept {
199 : return roundoff_hi;
200 : }
201 :
202 : //! Returns the overall size of the domain by multiplying the ProbLength's together
203 : [[nodiscard]] Real ProbSize () const noexcept
204 : {
205 : return AMREX_D_TERM(prob_domain.length(0),*prob_domain.length(1),*prob_domain.length(2));
206 : }
207 : //! Returns length of problem domain in specified dimension.
208 : [[nodiscard]] Real ProbLength (int dir) const noexcept { return prob_domain.length(dir); }
209 : //! Returns our rectangular domain.
210 4071290 : [[nodiscard]] const Box& Domain () const noexcept { return domain; }
211 : //! Sets our rectangular domain.
212 : void Domain (const Box& bx) noexcept
213 : {
214 : AMREX_ASSERT(bx.cellCentered());
215 : domain = bx;
216 : computeRoundoffDomain();
217 : }
218 : //! Define a multifab of areas and volumes with given grow factor.
219 : void GetVolume (MultiFab& vol,
220 : const BoxArray& grds,
221 : const DistributionMapping& dm,
222 : int grow) const;
223 : //! Fill the pre-built multifab with volume
224 : void GetVolume (MultiFab& vol) const;
225 :
226 : void GetVolume (FArrayBox& vol,
227 : const BoxArray& grds,
228 : int idx,
229 : int grow) const;
230 :
231 : //! Return the volume of the specified cell.
232 : AMREX_GPU_HOST_DEVICE AMREX_INLINE
233 : static Real Volume (const IntVect& point, const GeometryData& geomdata)
234 : {
235 : const auto *dx = geomdata.CellSize();
236 :
237 : Real vol;
238 :
239 : #if AMREX_SPACEDIM == 1
240 :
241 : auto coord = geomdata.Coord();
242 :
243 : if (coord == CoordSys::cartesian) {
244 : // Cartesian
245 :
246 : vol = dx[0];
247 : }
248 : else {
249 : // Spherical
250 :
251 : Real rl = geomdata.ProbLo()[0] + static_cast<Real>(point[0]) * dx[0];
252 : Real rr = rl + dx[0];
253 :
254 : constexpr Real pi = Real(3.1415926535897932);
255 : vol = (4.0_rt / 3.0_rt) * pi * dx[0] * (rl * rl + rl * rr + rr * rr);
256 : }
257 :
258 : #elif AMREX_SPACEDIM == 2
259 :
260 : auto coord = geomdata.Coord();
261 :
262 : if (coord == CoordSys::cartesian) {
263 : // Cartesian
264 :
265 : vol = dx[0] * dx[1];
266 : }
267 : else {
268 : // Cylindrical
269 :
270 : Real r_l = geomdata.ProbLo()[0] + static_cast<Real>(point[0]) * dx[0];
271 : Real r_r = geomdata.ProbLo()[0] + static_cast<Real>(point[0]+1) * dx[0];
272 :
273 : constexpr Real pi = Real(3.1415926535897932);
274 : vol = pi * (r_l + r_r) * dx[0] * dx[1];
275 : }
276 :
277 : #else
278 :
279 : amrex::ignore_unused(point);
280 :
281 : // Cartesian
282 :
283 : vol = dx[0] * dx[1] * dx[2];
284 :
285 : #endif
286 :
287 : return vol;
288 : }
289 :
290 : /**
291 : * \brief Compute d(log(A))/dr at cell centers in given region and
292 : * stuff the results into the passed MultiFab.
293 : */
294 : void GetDLogA (MultiFab& dloga,
295 : const BoxArray& grds,
296 : const DistributionMapping& dm,
297 : int dir,
298 : int grow) const;
299 : /**
300 : * \brief Compute area of cell faces in given region and stuff
301 : * stuff the results into the passed MultiFab.
302 : */
303 : void GetFaceArea (MultiFab& area,
304 : const BoxArray& grds,
305 : const DistributionMapping& dm,
306 : int dir,
307 : int grow) const;
308 : //! Fill the pre-built multifab with area
309 : void GetFaceArea (MultiFab& area,
310 : int dir) const;
311 :
312 : void GetFaceArea (FArrayBox& area,
313 : const BoxArray& grds,
314 : int idx,
315 : int dir,
316 : int grow) const;
317 : //! Is the domain periodic in the specified direction?
318 4732180 : [[nodiscard]] bool isPeriodic (int dir) const noexcept { return is_periodic[dir]; }
319 : //! Is domain periodic in any direction?
320 126 : [[nodiscard]] bool isAnyPeriodic () const noexcept
321 : {
322 126 : return AMREX_D_TERM(isPeriodic(0),||isPeriodic(1),||isPeriodic(2));
323 : }
324 : //! Is domain periodic in all directions?
325 : [[nodiscard]] bool isAllPeriodic () const noexcept
326 : {
327 : return AMREX_D_TERM(isPeriodic(0),&&isPeriodic(1),&&isPeriodic(2));
328 : }
329 : [[nodiscard]] Array<int,AMREX_SPACEDIM> isPeriodic () const noexcept {
330 : return {{AMREX_D_DECL(static_cast<int>(is_periodic[0]),
331 : static_cast<int>(is_periodic[1]),
332 : static_cast<int>(is_periodic[2]))}};
333 : }
334 : [[nodiscard]] GpuArray<int,AMREX_SPACEDIM> isPeriodicArray () const noexcept {
335 : return {{AMREX_D_DECL(static_cast<int>(is_periodic[0]),
336 : static_cast<int>(is_periodic[1]),
337 : static_cast<int>(is_periodic[2]))}};
338 : }
339 : //! What's period in specified direction?
340 : [[nodiscard]] int period (int dir) const noexcept { BL_ASSERT(is_periodic[dir]); return domain.length(dir); }
341 :
342 869672 : [[nodiscard]] Periodicity periodicity () const noexcept {
343 869672 : return Periodicity(IntVect(AMREX_D_DECL(domain.length(0) * is_periodic[0],
344 : domain.length(1) * is_periodic[1],
345 869672 : domain.length(2) * is_periodic[2])));
346 : }
347 :
348 0 : [[nodiscard]] Periodicity periodicity (const Box& b) const noexcept {
349 : AMREX_ASSERT(b.cellCentered());
350 0 : return Periodicity(IntVect(AMREX_D_DECL(b.length(0) * is_periodic[0],
351 : b.length(1) * is_periodic[1],
352 0 : b.length(2) * is_periodic[2])));
353 : }
354 :
355 : /**
356 : * \brief Compute Array of shifts which will translate src so that it will
357 : * intersect target with non-zero intersection. the array will be
358 : * resized internally, so anything previously there will be gone
359 : * DO NOT return non-periodic shifts, even if the box's do
360 : * intersect without shifting. The logic is that you will only do
361 : * this as a special case if there is some periodicity.
362 : */
363 : void periodicShift (const Box& target,
364 : const Box& src,
365 : Vector<IntVect>& out) const noexcept;
366 :
367 : //! Return domain box with non-periodic directions grown by ngrow.
368 : [[nodiscard]] Box growNonPeriodicDomain (IntVect const& ngrow) const noexcept;
369 : //! Return domain box with non-periodic directions grown by ngrow.
370 : [[nodiscard]] Box growNonPeriodicDomain (int ngrow) const noexcept;
371 : //! Return domain box with periodic directions grown by ngrow.
372 : [[nodiscard]] Box growPeriodicDomain (IntVect const& ngrow) const noexcept;
373 : //! Return domain box with periodic directions grown by ngrow.
374 : [[nodiscard]] Box growPeriodicDomain (int ngrow) const noexcept;
375 :
376 : //! Set periodicity flags and return the old flags.
377 : //! Note that, unlike Periodicity class, the flags are just boolean.
378 : //!
379 : Array<int,AMREX_SPACEDIM>
380 : setPeriodicity (Array<int,AMREX_SPACEDIM> const& period) noexcept {
381 : Array<int,AMREX_SPACEDIM> r{{AMREX_D_DECL(is_periodic[0],
382 : is_periodic[1],
383 : is_periodic[2])}};
384 : AMREX_D_TERM(is_periodic[0] = period[0];,
385 : is_periodic[1] = period[1];,
386 : is_periodic[2] = period[2];);
387 : return r;
388 : }
389 :
390 : void coarsen (IntVect const& rr) {
391 : domain.coarsen(rr);
392 : for (int i = 0; i < AMREX_SPACEDIM; ++i) {
393 : dx[i] = (ProbHi(i) - ProbLo(i)) / static_cast<Real>(domain.length(i));
394 : inv_dx[i] = Real(1.)/dx[i];
395 : }
396 : }
397 :
398 : void refine (IntVect const& rr) {
399 : domain.refine(rr);
400 : for (int i = 0; i < AMREX_SPACEDIM; ++i) {
401 : dx[i] = (ProbHi(i) - ProbLo(i)) / static_cast<Real>(domain.length(i));
402 : inv_dx[i] = Real(1.)/dx[i];
403 : }
404 : }
405 :
406 : /**
407 : * \brief Returns true if a point is outside the roundoff domain.
408 : * All particles with positions inside the roundoff domain
409 : * are sure to be mapped to cells inside the Domain() box. Note that
410 : * the same need not be true for all points inside ProbDomain().
411 : */
412 : [[nodiscard]] bool outsideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const;
413 :
414 : /**
415 : * \brief Returns true if a point is inside the roundoff domain.
416 : * All particles with positions inside the roundoff domain
417 : * are sure to be mapped to cells inside the Domain() box. Note that
418 : * the same need not be true for all points inside ProbDomain().
419 : */
420 : [[nodiscard]] bool insideRoundoffDomain (AMREX_D_DECL(ParticleReal x, ParticleReal y, ParticleReal z)) const;
421 :
422 : /**
423 : * \brief Compute the roundoff domain. Public because it contains an
424 : * extended host / device lambda.
425 : */
426 : void computeRoundoffDomain ();
427 :
428 : private:
429 : void read_params ();
430 :
431 : // is_periodic and RealBox used to be static
432 : bool is_periodic[AMREX_SPACEDIM] = {AMREX_D_DECL(false,false,false)};
433 : RealBox prob_domain;
434 :
435 : // Due to round-off errors, not all floating point numbers for which plo >= x < phi
436 : // will map to a cell that is inside "domain". "roundoff_{lo,hi}" store
437 : // positions that are very close to that in prob_domain, and for which all doubles and floats
438 : // between them will map to a cell inside domain.
439 : GpuArray<ParticleReal , AMREX_SPACEDIM> roundoff_lo, roundoff_hi;
440 :
441 : //
442 : Box domain;
443 :
444 : friend std::istream& operator>> (std::istream&, Geometry&);
445 : };
446 :
447 :
448 : //! Nice ASCII output.
449 : std::ostream& operator<< (std::ostream&, const Geometry&);
450 : //! Nice ASCII input.
451 : std::istream& operator>> (std::istream&, Geometry&);
452 :
453 : inline
454 : Geometry
455 : coarsen (Geometry const& fine, IntVect const& rr) {
456 : Geometry r{fine};
457 : r.coarsen(rr);
458 : return r;
459 : }
460 :
461 : inline
462 : Geometry
463 : coarsen (Geometry const& fine, int rr) { return coarsen(fine, IntVect(rr)); }
464 :
465 : inline
466 : Geometry
467 : refine (Geometry const& crse, IntVect const& rr) {
468 : Geometry r{crse};
469 : r.refine(rr);
470 : return r;
471 : }
472 :
473 : inline
474 : Geometry
475 : refine (Geometry const& crse, int rr) { return refine(crse, IntVect(rr)); }
476 :
477 : inline
478 : const Geometry&
479 : DefaultGeometry () {
480 : return *(AMReX::top()->getDefaultGeometry());
481 : }
482 :
483 : }
484 :
485 : #endif /*_GEOMETRY_H_*/
|
