Line data Source code
1 :
2 : #ifndef AMREX_FABSET_H_
3 : #define AMREX_FABSET_H_
4 : #include <AMReX_Config.H>
5 :
6 : #include <AMReX_FabDataType.H>
7 : #include <AMReX_MultiFab.H>
8 : #include <AMReX_ParallelDescriptor.H>
9 : #include <AMReX_BLProfiler.H>
10 : #include <AMReX_VisMF.H>
11 :
12 : #ifdef AMREX_USE_OMP
13 : #include <omp.h>
14 : #endif
15 :
16 : #include <limits>
17 :
18 : namespace amrex {
19 :
20 : /**
21 : \brief A FabSet is a group of FArrayBox's. The grouping is designed
22 : specifically to represent regions along the boundary of Box's,
23 : and are used to implement boundary conditions to discretized
24 : partial differential equations.
25 :
26 : A FabSet is an array of pointers to FABs. The standard FAB operators,
27 : however, have been modified to be more useful for maintaining
28 : boundary conditions for partial differential equations discretized
29 : on boxes.
30 : Under normal circumstances, a FAB will be created for each face of a
31 : box. For a group of boxes, a FabSet will be the group of FABs at a
32 : particular orientation (ie. the lo-i side of each grid in a list).
33 :
34 : Since a FabSet FAB will likely be used to bound a grid box,
35 : FArrayBox::resize() operations are disallowed. Also, to preserve
36 : flexibility in applicable boundary scenarios, intersecting
37 : FABs in the FabSet are not guaranteed to contain identical data--thus
38 : copy operations from a FabSet to any FAB-like structure may be
39 : order-dependent.
40 :
41 : FabSets are used primarily as a data storage mechanism, and are
42 : manipulated by more sophisticated control classes.
43 : */
44 : template <typename MF>
45 : class FabSetT
46 : {
47 : friend class FabSetIter;
48 : friend class FluxRegister;
49 : public:
50 : using value_type = typename FabDataType<MF>::value_type;
51 : using FAB = typename FabDataType<MF>::fab_type;
52 :
53 : //
54 : //! The default constructor -- you must later call define().
55 : FabSetT () noexcept = default;
56 : //
57 : //! Construct a FabSetT<MF> of specified number of components on the grids.
58 : FabSetT (const BoxArray& grids, const DistributionMapping& dmap, int ncomp);
59 :
60 0 : ~FabSetT () = default;
61 :
62 : FabSetT (FabSetT<MF>&& rhs) noexcept = default;
63 :
64 : FabSetT (const FabSetT<MF>& rhs) = delete;
65 : FabSetT<MF>& operator= (const FabSetT<MF>& rhs) = delete;
66 : FabSetT<MF>& operator= (FabSetT<MF>&& rhs) = delete;
67 :
68 : //
69 : //! Define a FabSetT<MF> constructed via default constructor.
70 : void define (const BoxArray& grids, const DistributionMapping& dmap, int ncomp);
71 :
72 : FAB const& operator[] (const MFIter& mfi) const noexcept { return m_mf[mfi]; }
73 : FAB & operator[] (const MFIter& mfi) noexcept { return m_mf[mfi]; }
74 : FAB const& operator[] (int i) const noexcept { return m_mf[i]; }
75 : FAB & operator[] (int i) noexcept { return m_mf[i]; }
76 :
77 : Array4<value_type const> array (const MFIter& mfi) const noexcept { return m_mf.const_array(mfi); }
78 : Array4<value_type > array (const MFIter& mfi) noexcept { return m_mf.array(mfi); }
79 : Array4<value_type const> array (int i) const noexcept { return m_mf.const_array(i); }
80 : Array4<value_type > array (int i) noexcept { return m_mf.array(i); }
81 : Array4<value_type const> const_array (const MFIter& mfi) const noexcept { return m_mf.const_array(mfi); }
82 : Array4<value_type const> const_array (int i) const noexcept { return m_mf.const_array(i); }
83 :
84 : MultiArray4<value_type const> arrays () const noexcept { return m_mf.const_arrays(); }
85 : MultiArray4<value_type > arrays () noexcept { return m_mf.arrays(); }
86 : MultiArray4<value_type const> const_arrays () const noexcept { return m_mf.const_arrays(); }
87 :
88 : [[nodiscard]] Box fabbox (int K) const noexcept { return m_mf.fabbox(K); }
89 :
90 : [[nodiscard]] int size () const noexcept { return m_mf.size(); }
91 :
92 : [[nodiscard]] const BoxArray& boxArray () const noexcept { return m_mf.boxArray(); }
93 :
94 : [[nodiscard]] const DistributionMapping& DistributionMap () const noexcept
95 : { return m_mf.DistributionMap(); }
96 :
97 : [[nodiscard]] MF & multiFab () noexcept { return m_mf; }
98 : [[nodiscard]] MF const& multiFab () const noexcept { return m_mf; }
99 :
100 : [[nodiscard]] int nComp () const noexcept { return m_mf.nComp(); }
101 :
102 : void clear () { m_mf.clear(); }
103 :
104 : FabSetT<MF>& copyFrom (const FabSetT<MF>& src, int scomp, int dcomp, int ncomp);
105 :
106 : FabSetT<MF>& copyFrom (const MF& src, int ngrow, int scomp, int dcomp, int ncomp,
107 : const Periodicity& period = Periodicity::NonPeriodic());
108 :
109 : FabSetT<MF>& plusFrom (const FabSetT<MF>& src, int scomp, int dcomp, int ncomp);
110 :
111 : FabSetT<MF>& plusFrom (const MF& src, int ngrow, int scomp, int dcomp, int ncomp,
112 : const Periodicity& period = Periodicity::NonPeriodic());
113 :
114 : void copyTo (MF& dest, int ngrow, int scomp, int dcomp, int ncomp,
115 : const Periodicity& period = Periodicity::NonPeriodic()) const;
116 :
117 : void plusTo (MF& dest, int ngrow, int scomp, int dcomp, int ncomp,
118 : const Periodicity& period = Periodicity::NonPeriodic()) const;
119 :
120 : void setVal (value_type val);
121 :
122 : void setVal (value_type val, int comp, int num_comp);
123 :
124 : //!< Linear combination: this := a*this + b*src (FabSetT<MF>s must be commensurate).
125 : FabSetT<MF>& linComb (value_type a, value_type b, const FabSetT<MF>& src,
126 : int scomp, int dcomp, int ncomp);
127 :
128 : //!< Linear combination: this := a*mfa + b*mfb
129 : FabSetT<MF>& linComb (value_type a, const MF& mfa, int a_comp,
130 : value_type b, const MF& mfb, int b_comp,
131 : int dcomp, int ncomp, int ngrow);
132 :
133 : //
134 : //! Write (used for writing to checkpoint)
135 : void write (const std::string& name) const;
136 : //
137 : //! Read (used for reading from checkpoint)
138 : void read (const std::string& name);
139 :
140 : //!< Local copy function
141 : static void Copy (FabSetT<MF>& dst, const FabSetT<MF>& src);
142 :
143 : private:
144 : MF m_mf;
145 : };
146 :
147 : class FabSetIter
148 : : public MFIter
149 : {
150 : public:
151 : template <typename MF>
152 : explicit FabSetIter (const FabSetT<MF>& fs)
153 : : MFIter(fs.m_mf) { }
154 : };
155 :
156 : template <typename MF>
157 : FabSetT<MF>::FabSetT (const BoxArray& grids, const DistributionMapping& dmap, int ncomp)
158 : : m_mf(grids,dmap,ncomp,0)
159 : {}
160 :
161 : template <typename MF>
162 : void
163 : FabSetT<MF>::define (const BoxArray& grids, const DistributionMapping& dm, int ncomp)
164 : {
165 : m_mf.define(grids, dm, ncomp, 0);
166 : }
167 :
168 : template <typename MF>
169 : FabSetT<MF>&
170 : FabSetT<MF>::copyFrom (const FabSetT<MF>& src, int scomp, int dcomp, int ncomp)
171 : {
172 : if (boxArray() == src.boxArray() && DistributionMap() == src.DistributionMap()) {
173 : #ifdef AMREX_USE_OMP
174 : #pragma omp parallel if (Gpu::notInLaunchRegion())
175 : #endif
176 : for (FabSetIter fsi(*this); fsi.isValid(); ++fsi) {
177 : const Box& bx = fsi.validbox();
178 : auto const srcfab = src.array(fsi);
179 : auto dstfab = this->array(fsi);
180 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
181 : {
182 : dstfab(i,j,k,n+dcomp) = srcfab(i,j,k,n+scomp);
183 : });
184 : }
185 : } else {
186 : m_mf.ParallelCopy(src.m_mf,scomp,dcomp,ncomp);
187 : }
188 : return *this;
189 : }
190 :
191 : template <typename MF>
192 : FabSetT<MF>&
193 : FabSetT<MF>::copyFrom (const MF& src, int ngrow, int scomp, int dcomp, int ncomp,
194 : const Periodicity& period)
195 : {
196 : BL_ASSERT(boxArray() != src.boxArray());
197 : m_mf.ParallelCopy(src,scomp,dcomp,ncomp,ngrow,0,period);
198 : return *this;
199 : }
200 :
201 : template <typename MF>
202 : FabSetT<MF>&
203 : FabSetT<MF>::plusFrom (const FabSetT<MF>& src, int scomp, int dcomp, int ncomp)
204 : {
205 : if (boxArray() == src.boxArray() && DistributionMap() == src.DistributionMap()) {
206 : #ifdef AMREX_USE_OMP
207 : #pragma omp parallel if (Gpu::notInLaunchRegion())
208 : #endif
209 : for (FabSetIter fsi(*this); fsi.isValid(); ++fsi) {
210 : const Box& bx = fsi.validbox();
211 : auto const srcfab = src.array(fsi);
212 : auto dstfab = this->array(fsi);
213 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
214 : {
215 : dstfab(i,j,k,n+dcomp) += srcfab(i,j,k,n+scomp);
216 : });
217 : }
218 : } else {
219 : amrex::Abort("FabSetT<MF>::plusFrom: parallel plusFrom not supported");
220 : }
221 : return *this;
222 : }
223 :
224 : template <typename MF>
225 : FabSetT<MF>&
226 : FabSetT<MF>::plusFrom (const MF& src, int ngrow, int scomp, int dcomp, int ncomp,
227 : const Periodicity& period)
228 : {
229 : BL_ASSERT(boxArray() != src.boxArray());
230 : m_mf.ParallelCopy(src,scomp,dcomp,ncomp,ngrow,0,period,FabArrayBase::ADD);
231 : return *this;
232 : }
233 :
234 : template <typename MF>
235 : void
236 : FabSetT<MF>::copyTo (MF& dest, int ngrow, int scomp, int dcomp, int ncomp,
237 : const Periodicity& period) const
238 : {
239 : BL_ASSERT(boxArray() != dest.boxArray());
240 : dest.ParallelCopy(m_mf,scomp,dcomp,ncomp,0,ngrow,period);
241 : }
242 :
243 : template <typename MF>
244 : void
245 : FabSetT<MF>::plusTo (MF& dest, int ngrow, int scomp, int dcomp, int ncomp,
246 : const Periodicity& period) const
247 : {
248 : BL_ASSERT(boxArray() != dest.boxArray());
249 : dest.ParallelCopy(m_mf,scomp,dcomp,ncomp,0,ngrow,period,FabArrayBase::ADD);
250 : }
251 :
252 : template <typename MF>
253 : void
254 : FabSetT<MF>::setVal (value_type val)
255 : {
256 : const int ncomp = nComp();
257 : #ifdef AMREX_USE_OMP
258 : #pragma omp parallel if (Gpu::notInLaunchRegion())
259 : #endif
260 : for (FabSetIter fsi(*this); fsi.isValid(); ++fsi) {
261 : const Box& bx = fsi.validbox();
262 : auto fab = this->array(fsi);
263 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
264 : {
265 : fab(i,j,k,n) = val;
266 : });
267 : }
268 : }
269 :
270 : template <typename MF>
271 : void
272 : FabSetT<MF>::setVal (value_type val, int comp, int num_comp)
273 : {
274 : #ifdef AMREX_USE_OMP
275 : #pragma omp parallel if (Gpu::notInLaunchRegion())
276 : #endif
277 : for (FabSetIter fsi(*this); fsi.isValid(); ++fsi) {
278 : const Box& bx = fsi.validbox();
279 : auto fab = this->array(fsi);
280 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, num_comp, i, j, k, n,
281 : {
282 : fab(i,j,k,n+comp) = val;
283 : });
284 : }
285 : }
286 :
287 : // Linear combination this := a*this + b*src
288 : // Note: corresponding fabsets must be commensurate.
289 : template <typename MF>
290 : FabSetT<MF>&
291 : FabSetT<MF>::linComb (value_type a, value_type b, const FabSetT<MF>& src,
292 : int scomp, int dcomp, int ncomp)
293 : {
294 : BL_ASSERT(size() == src.size());
295 :
296 : #ifdef AMREX_USE_OMP
297 : #pragma omp parallel if (Gpu::notInLaunchRegion())
298 : #endif
299 : for (FabSetIter fsi(*this); fsi.isValid(); ++fsi)
300 : {
301 : const Box& bx = fsi.validbox();
302 : auto const srcfab = src.array(fsi);
303 : auto dstfab = this->array(fsi);
304 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
305 : {
306 : dstfab(i,j,k,n+dcomp) = a*dstfab(i,j,k,n+dcomp) + b*srcfab(i,j,k,n+scomp);
307 : });
308 : }
309 : return *this;
310 : }
311 :
312 : // Linear combination: this := a*mfa + b*mfb
313 : // CastroRadiation is the only code that uses this function.
314 : template <typename MF>
315 : FabSetT<MF>&
316 : FabSetT<MF>::linComb (value_type a, const MF& mfa, int a_comp,
317 : value_type b, const MF& mfb, int b_comp,
318 : int dcomp, int ncomp, int ngrow)
319 : {
320 : BL_PROFILE("FabSetT<MF>::linComb()");
321 : BL_ASSERT(mfa.nGrowVect().allGE(ngrow) && mfb.nGrowVect().allGE(ngrow));
322 : BL_ASSERT(mfa.boxArray() == mfb.boxArray());
323 : BL_ASSERT(boxArray() != mfa.boxArray());
324 :
325 : MF bdrya(boxArray(),DistributionMap(),ncomp,0,MFInfo());
326 : MF bdryb(boxArray(),DistributionMap(),ncomp,0,MFInfo());
327 :
328 : const auto huge = static_cast<value_type>(sizeof(value_type) == 8 ? 1.e200 : 1.e30);
329 :
330 : #ifdef AMREX_USE_OMP
331 : #pragma omp parallel if (Gpu::notInLaunchRegion())
332 : #endif
333 : for (MFIter mfi(bdrya); mfi.isValid(); ++mfi) // tiling is not safe for this BoxArray
334 : {
335 : const Box& bx = mfi.validbox();
336 : auto afab = bdrya.array(mfi);
337 : auto bfab = bdryb.array(mfi);
338 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
339 : {
340 : afab(i,j,k,n) = huge;
341 : bfab(i,j,k,n) = huge;
342 : });
343 : }
344 :
345 : bdrya.ParallelCopy(mfa,a_comp,0,ncomp,ngrow,0);
346 : bdryb.ParallelCopy(mfb,b_comp,0,ncomp,ngrow,0);
347 :
348 : #ifdef AMREX_USE_OMP
349 : #pragma omp parallel if (Gpu::notInLaunchRegion())
350 : #endif
351 : for (FabSetIter fsi(*this); fsi.isValid(); ++fsi)
352 : {
353 : const Box& bx = fsi.validbox();
354 : auto const afab = bdrya.array(fsi);
355 : auto const bfab = bdryb.array(fsi);
356 : auto dfab = this->array(fsi);
357 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
358 : {
359 : dfab(i,j,k,n+dcomp) = a*afab(i,j,k,n) + b*bfab(i,j,k,n);
360 : });
361 : }
362 :
363 : return *this;
364 : }
365 :
366 : template <typename MF>
367 : void
368 : FabSetT<MF>::write (const std::string& name) const
369 : {
370 : if (AsyncOut::UseAsyncOut()) {
371 : VisMF::AsyncWrite(m_mf,name);
372 : } else {
373 : VisMF::Write(m_mf,name);
374 : }
375 : }
376 :
377 : template <typename MF>
378 : void
379 : FabSetT<MF>::read (const std::string& name)
380 : {
381 : if (m_mf.empty()) {
382 : amrex::Abort("FabSetT<MF>::read: not predefined");
383 : }
384 : VisMF::Read(m_mf,name);
385 : }
386 :
387 : template <typename MF>
388 : void
389 : FabSetT<MF>::Copy (FabSetT<MF>& dst, const FabSetT<MF>& src)
390 : {
391 : BL_ASSERT(amrex::match(dst.boxArray(), src.boxArray()));
392 : BL_ASSERT(dst.DistributionMap() == src.DistributionMap());
393 : int ncomp = dst.nComp();
394 : #ifdef AMREX_USE_OMP
395 : #pragma omp parallel if (Gpu::notInLaunchRegion())
396 : #endif
397 : for (FabSetIter fsi(dst); fsi.isValid(); ++fsi) {
398 : const Box& bx = fsi.validbox();
399 : auto const srcfab = src.array(fsi);
400 : auto dstfab = dst.array(fsi);
401 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
402 : {
403 : dstfab(i,j,k,n) = srcfab(i,j,k,n);
404 : });
405 : }
406 : }
407 :
408 : using FabSet = FabSetT<MultiFab>;
409 : using fFabSet = FabSetT<fMultiFab>;
410 :
411 : }
412 :
413 : #endif /*_FABSET_H_*/
|
