Block-Structured AMR Software Framework
Loading...
Searching...
No Matches
AMReX_BaseFab.H
Go to the documentation of this file.
1#ifndef AMREX_BASEFAB_H_
2#define AMREX_BASEFAB_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_Algorithm.H>
6#include <AMReX_Extension.H>
7#include <AMReX_BLassert.H>
8#include <AMReX_Array.H>
9#include <AMReX_Box.H>
10#include <AMReX_Loop.H>
11#include <AMReX_BoxList.H>
12#include <AMReX_BArena.H>
13#include <AMReX_CArena.H>
14#include <AMReX_DataAllocator.H>
15#include <AMReX_REAL.H>
16#include <AMReX_BLProfiler.H>
17#include <AMReX_BoxIterator.H>
18#include <AMReX_MakeType.H>
19#include <AMReX_Utility.H>
20#include <AMReX_Reduce.H>
21#include <AMReX_Gpu.H>
22#include <AMReX_Scan.H>
23#include <AMReX_Math.H>
24#include <AMReX_OpenMP.H>
25#include <AMReX_MemPool.H>
26#include <AMReX_TypeTraits.H>
27
28#include <cmath>
29#include <cstdlib>
30#include <algorithm>
31#include <limits>
32#include <climits>
33#include <array>
34#include <type_traits>
35#include <memory>
36#include <atomic>
37#include <utility>
38
39namespace amrex
40{
41
43
44extern std::atomic<Long> atomic_total_bytes_allocated_in_fabs;
45extern std::atomic<Long> atomic_total_bytes_allocated_in_fabs_hwm;
46extern std::atomic<Long> atomic_total_cells_allocated_in_fabs;
47extern std::atomic<Long> atomic_total_cells_allocated_in_fabs_hwm;
48extern Long private_total_bytes_allocated_in_fabs;
49extern Long private_total_bytes_allocated_in_fabs_hwm;
50extern Long private_total_cells_allocated_in_fabs;
51extern Long private_total_cells_allocated_in_fabs_hwm;
52#ifdef AMREX_USE_OMP
53#pragma omp threadprivate(private_total_bytes_allocated_in_fabs)
54#pragma omp threadprivate(private_total_bytes_allocated_in_fabs_hwm)
55#pragma omp threadprivate(private_total_cells_allocated_in_fabs)
56#pragma omp threadprivate(private_total_cells_allocated_in_fabs_hwm)
57#endif
58
60
66void update_fab_stats (Long n, Long s, std::size_t szt) noexcept;
67
68void BaseFab_Initialize ();
69void BaseFab_Finalize ();
70
71struct SrcComp {
73 explicit SrcComp (int ai) noexcept : i(ai) {}
74 int i;
75};
76
77struct DestComp {
79 explicit DestComp (int ai) noexcept : i(ai) {}
80 int i;
81};
82
83struct NumComps {
85 explicit NumComps (int an) noexcept : n(an) {}
86 int n;
87};
88
89template <typename T>
92makeArray4 (T* p, Box const& bx, int ncomp) noexcept
93{
94 return Array4<T>{p, amrex::begin(bx), amrex::end(bx), ncomp};
95}
96
97template <typename T>
98std::enable_if_t<std::is_arithmetic_v<T>>
99placementNew (T* const /*ptr*/, Long /*n*/)
100{}
101
102template <typename T>
103std::enable_if_t<std::is_trivially_default_constructible_v<T>
104 && !std::is_arithmetic_v<T>>
105placementNew (T* const ptr, Long n)
106{
107 for (Long i = 0; i < n; ++i) {
108 new (ptr+i) T;
109 }
110}
111
112template <typename T>
113std::enable_if_t<!std::is_trivially_default_constructible_v<T>>
114placementNew (T* const ptr, Long n)
115{
117 {
118 new (ptr+i) T;
119 });
120}
121
122template <typename T>
123std::enable_if_t<std::is_trivially_destructible_v<T>>
124placementDelete (T* const /*ptr*/, Long /*n*/)
125{}
126
127template <typename T>
128std::enable_if_t<!std::is_trivially_destructible_v<T>>
129placementDelete (T* const ptr, Long n)
130{
132 {
133 (ptr+i)->~T();
134 });
135}
136
187template <class T>
189 : public DataAllocator
190{
191public:
192
193 template <class U> friend class BaseFab;
194
195 using value_type = T;
196
198 BaseFab () noexcept = default;
199
200 explicit BaseFab (Arena* ar) noexcept;
201
202 BaseFab (const Box& bx, int n, Arena* ar);
203
205 explicit BaseFab (const Box& bx, int n = 1, bool alloc = true,
206 bool shared = false, Arena* ar = nullptr);
207
208 BaseFab (const BaseFab<T>& rhs, MakeType make_type, int scomp, int ncomp);
209
215 BaseFab (const Box& bx, int ncomp, T* p);
216 BaseFab (const Box& bx, int ncomp, T const* p);
217
218 explicit BaseFab (Array4<T> const& a) noexcept;
219
220 explicit BaseFab (Array4<T> const& a, IndexType t) noexcept;
221
222 explicit BaseFab (Array4<T const> const& a) noexcept;
223
224 explicit BaseFab (Array4<T const> const& a, IndexType t) noexcept;
225
227 virtual ~BaseFab () noexcept;
228
229 BaseFab (const BaseFab<T>& rhs) = delete;
230 BaseFab<T>& operator= (const BaseFab<T>& rhs) = delete;
231
232 BaseFab (BaseFab<T>&& rhs) noexcept;
233 BaseFab<T>& operator= (BaseFab<T>&& rhs) noexcept;
234
235 template <RunOn run_on AMREX_DEFAULT_RUNON>
236 BaseFab& operator= (T const&) noexcept;
237
238 static void Initialize();
239 static void Finalize();
240
254 void resize (const Box& b, int N = 1, Arena* ar = nullptr);
255
256 template <class U=T, std::enable_if_t<std::is_trivially_destructible_v<U>,int> = 0>
257 [[nodiscard]] Elixir elixir () noexcept;
258
263 void clear () noexcept;
264
266 [[nodiscard]] std::unique_ptr<T,DataDeleter> release () noexcept;
267
269 [[nodiscard]] std::size_t nBytes () const noexcept { return this->truesize*sizeof(T); }
270
271 [[nodiscard]] std::size_t nBytesOwned () const noexcept {
272 return (this->ptr_owner) ? nBytes() : 0;
273 }
274
276 [[nodiscard]] std::size_t nBytes (const Box& bx, int ncomps) const noexcept
277 { return bx.numPts() * sizeof(T) * ncomps; }
278
280 [[nodiscard]] int nComp () const noexcept { return this->nvar; }
281
283 [[nodiscard]] const int* nCompPtr() const noexcept {
284 return &(this->nvar);
285 }
286
288 [[nodiscard]] Long numPts () const noexcept { return this->domain.numPts(); }
289
291 [[nodiscard]] Long size () const noexcept { return this->nvar*this->domain.numPts(); }
292
294 [[nodiscard]] const Box& box () const noexcept { return this->domain; }
295
300 [[nodiscard]] IntVect length () const noexcept { return this->domain.length(); }
301
306 [[nodiscard]] const IntVect& smallEnd () const noexcept { return this->domain.smallEnd(); }
307
309 [[nodiscard]] const IntVect& bigEnd () const noexcept { return this->domain.bigEnd(); }
310
319 [[nodiscard]] const int* loVect () const noexcept { return this->domain.loVect(); }
320
329 [[nodiscard]] const int* hiVect () const noexcept { return this->domain.hiVect(); }
330
335 [[nodiscard]] bool contains (const BaseFab<T>& fab) const noexcept
336 {
337 return box().contains(fab.box()) && this->nvar <= fab.nvar;
338 }
339
344 [[nodiscard]] bool contains (const Box& bx) const noexcept { return box().contains(bx); }
345
355 [[nodiscard]] T* dataPtr (int n = 0) noexcept {
356 if (this->dptr) {
357 return &(this->dptr[n*this->domain.numPts()]);
358 } else {
359 return nullptr;
360 }
361 }
362
364 [[nodiscard]] const T* dataPtr (int n = 0) const noexcept {
365 if (this->dptr) {
366 return &(this->dptr[n*this->domain.numPts()]);
367 } else {
368 return nullptr;
369 }
370 }
371
372 [[nodiscard]] T* dataPtr (const IntVect& p, int n = 0) noexcept;
373
374 [[nodiscard]] const T* dataPtr (const IntVect& p, int n = 0) const noexcept;
375
376 void setPtr (T* p, Long sz) noexcept { AMREX_ASSERT(this->dptr == nullptr && this->truesize == 0); this->dptr = p; this->truesize = sz; }
377
378 void prefetchToHost () const noexcept;
379 void prefetchToDevice () const noexcept;
380
381 [[nodiscard]] AMREX_FORCE_INLINE
382 Array4<T const> array () const noexcept
383 {
384 return makeArray4<T const>(this->dptr, this->domain, this->nvar);
385 }
386
387 [[nodiscard]] AMREX_FORCE_INLINE
388 Array4<T const> array (int start_comp) const noexcept
389 {
390 return Array4<T const>(makeArray4<T const>(this->dptr, this->domain, this->nvar),start_comp);
391 }
392
393 [[nodiscard]] AMREX_FORCE_INLINE
394 Array4<T const> array (int start_comp, int num_comps) const noexcept
395 {
396 return Array4<T const>(makeArray4<T const>(this->dptr, this->domain, this->nvar), start_comp, num_comps);
397 }
398
399 [[nodiscard]] AMREX_FORCE_INLINE
400 Array4<T> array () noexcept
401 {
402 return makeArray4<T>(this->dptr, this->domain, this->nvar);
403 }
404
405 [[nodiscard]] AMREX_FORCE_INLINE
406 Array4<T> array (int start_comp) noexcept
407 {
408 return Array4<T>(makeArray4<T>(this->dptr, this->domain, this->nvar),start_comp);
409 }
410
411 [[nodiscard]] AMREX_FORCE_INLINE
412 Array4<T> array (int start_comp, int num_comps) noexcept
413 {
414 return Array4<T>(makeArray4<T>(this->dptr, this->domain, this->nvar), start_comp, num_comps);
415 }
416
417 [[nodiscard]] AMREX_FORCE_INLINE
418 Array4<T const> const_array () const noexcept
419 {
420 return makeArray4<T const>(this->dptr, this->domain, this->nvar);
421 }
422
423 [[nodiscard]] AMREX_FORCE_INLINE
424 Array4<T const> const_array (int start_comp) const noexcept
425 {
426 return Array4<T const>(makeArray4<T const>(this->dptr, this->domain, this->nvar),start_comp);
427 }
428
429 [[nodiscard]] AMREX_FORCE_INLINE
430 Array4<T const> const_array (int start_comp, int num_comps) const noexcept
431 {
432 return Array4<T const>(makeArray4<T const>(this->dptr, this->domain, this->nvar), start_comp, num_comps);
433 }
434
436 [[nodiscard]] bool isAllocated () const noexcept { return this->dptr != nullptr; }
437
444 [[nodiscard]] T& operator() (const IntVect& p, int N) noexcept;
445
447 [[nodiscard]] T& operator() (const IntVect& p) noexcept;
448
450 [[nodiscard]] const T& operator() (const IntVect& p, int N) const noexcept;
451
453 [[nodiscard]] const T& operator() (const IntVect& p) const noexcept;
454
460 void getVal (T* data, const IntVect& pos, int N, int numcomp) const noexcept;
462 void getVal (T* data, const IntVect& pos) const noexcept;
463
464 template <RunOn run_on AMREX_DEFAULT_RUNON,
465 class U=T, std::enable_if_t<std::is_same_v<U,float> || std::is_same_v<U,double>,int> FOO = 0>
466 void fill_snan () noexcept;
467
474 template <RunOn run_on AMREX_DEFAULT_RUNON>
475 void setVal (T const& x, const Box& bx, int dcomp, int ncomp) noexcept;
477 template <RunOn run_on AMREX_DEFAULT_RUNON>
478 void setVal (T const& x, const Box& bx, int N = 0) noexcept;
480 template <RunOn run_on AMREX_DEFAULT_RUNON>
481 void setVal (T const& x, int N) noexcept;
482
483 template <RunOn run_on AMREX_DEFAULT_RUNON>
484 void setValIfNot (T const& val, const Box& bx, const BaseFab<int>& mask, int nstart, int num) noexcept;
485
491 template <RunOn run_on AMREX_DEFAULT_RUNON>
492 void setComplement (T const& x, const Box& b, int ns, int num) noexcept;
493
510 template <RunOn run_on AMREX_DEFAULT_RUNON>
511 BaseFab<T>& copy (const BaseFab<T>& src, const Box& srcbox, int srccomp,
512 const Box& destbox, int destcomp, int numcomp) noexcept;
513
520 template <RunOn run_on AMREX_DEFAULT_RUNON>
521 BaseFab<T>& copy (const BaseFab<T>& src, int srccomp, int destcomp,
522 int numcomp = 1) noexcept;
529 template <RunOn run_on AMREX_DEFAULT_RUNON>
530 BaseFab<T>& copy (const BaseFab<T>& src, const Box& destbox) noexcept;
531
533 template <RunOn run_on AMREX_DEFAULT_RUNON>
534 std::size_t copyToMem (const Box& srcbox, int srccomp,
535 int numcomp, void* dst) const noexcept;
536
538 template <RunOn run_on AMREX_DEFAULT_RUNON, typename BUF = T>
539 std::size_t copyFromMem (const Box& dstbox, int dstcomp,
540 int numcomp, const void* src) noexcept;
541
543 template <RunOn run_on AMREX_DEFAULT_RUNON, typename BUF = T>
544 std::size_t addFromMem (const Box& dstbox, int dstcomp,
545 int numcomp, const void* src) noexcept;
546
552 BaseFab<T>& shift (const IntVect& v) noexcept;
558 BaseFab<T>& shift (int idir, int n_cell) noexcept;
564 BaseFab<T>& shiftHalf (int dir, int n_cell) noexcept;
570 BaseFab<T>& shiftHalf (const IntVect& v) noexcept;
571
572 template <RunOn run_on AMREX_DEFAULT_RUNON>
573 [[nodiscard]] Real norminfmask (const Box& subbox, const BaseFab<int>& mask, int scomp=0, int ncomp=1) const noexcept;
574
581 template <RunOn run_on AMREX_DEFAULT_RUNON>
582 [[nodiscard]] Real norm (int p, int scomp = 0, int numcomp = 1) const noexcept;
583
585 template <RunOn run_on AMREX_DEFAULT_RUNON>
586 [[nodiscard]] Real norm (const Box& subbox, int p, int scomp = 0, int numcomp = 1) const noexcept;
588 template <RunOn run_on AMREX_DEFAULT_RUNON>
589 void abs () noexcept;
591 template <RunOn run_on AMREX_DEFAULT_RUNON>
592 void abs (int comp, int numcomp=1) noexcept;
596 template <RunOn run_on AMREX_DEFAULT_RUNON>
597 void abs (const Box& subbox, int comp = 0, int numcomp=1) noexcept;
601 template <RunOn run_on AMREX_DEFAULT_RUNON>
602 [[nodiscard]] T min (int comp = 0) const noexcept;
606 template <RunOn run_on AMREX_DEFAULT_RUNON>
607 [[nodiscard]] T min (const Box& subbox, int comp = 0) const noexcept;
611 template <RunOn run_on AMREX_DEFAULT_RUNON>
612 [[nodiscard]] T max (int comp = 0) const noexcept;
616 template <RunOn run_on AMREX_DEFAULT_RUNON>
617 [[nodiscard]] T max (const Box& subbox, int comp = 0) const noexcept;
621 template <RunOn run_on AMREX_DEFAULT_RUNON>
622 [[nodiscard]] std::pair<T,T> minmax (int comp = 0) const noexcept;
626 template <RunOn run_on AMREX_DEFAULT_RUNON>
627 [[nodiscard]] std::pair<T,T> minmax (const Box& subbox, int comp = 0) const noexcept;
631 template <RunOn run_on AMREX_DEFAULT_RUNON>
632 [[nodiscard]] T maxabs (int comp = 0) const noexcept;
636 template <RunOn run_on AMREX_DEFAULT_RUNON>
637 [[nodiscard]] T maxabs (const Box& subbox, int comp = 0) const noexcept;
638
639 /*(
640 * \return location of given value
641 */
642 template <RunOn run_on AMREX_DEFAULT_RUNON>
643 [[nodiscard]] IntVect indexFromValue (const Box& subbox, int comp, T const& value) const noexcept;
644
648 template <RunOn run_on AMREX_DEFAULT_RUNON>
649 [[nodiscard]] IntVect minIndex (int comp = 0) const noexcept;
654 template <RunOn run_on AMREX_DEFAULT_RUNON>
655 [[nodiscard]] IntVect minIndex (const Box& subbox, int comp = 0) const noexcept;
660 template <RunOn run_on AMREX_DEFAULT_RUNON>
661 void minIndex (const Box& subbox, Real& min_val, IntVect& min_idx, int comp = 0) const noexcept;
662
666 template <RunOn run_on AMREX_DEFAULT_RUNON>
667 [[nodiscard]] IntVect maxIndex (int comp = 0) const noexcept;
672 template <RunOn run_on AMREX_DEFAULT_RUNON>
673 [[nodiscard]] IntVect maxIndex (const Box& subbox, int comp = 0) const noexcept;
678 template <RunOn run_on AMREX_DEFAULT_RUNON>
679 void maxIndex (const Box& subbox, Real& max_value, IntVect& max_idx, int comp = 0) const noexcept;
680
687 template <RunOn run_on AMREX_DEFAULT_RUNON>
688 int maskLT (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
690 template <RunOn run_on AMREX_DEFAULT_RUNON>
691 int maskLE (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
692
694 template <RunOn run_on AMREX_DEFAULT_RUNON>
695 int maskEQ (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
697 template <RunOn run_on AMREX_DEFAULT_RUNON>
698 int maskGT (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
700 template <RunOn run_on AMREX_DEFAULT_RUNON>
701 int maskGE (BaseFab<int>& mask, T const& val, int comp = 0) const noexcept;
702
704 template <RunOn run_on AMREX_DEFAULT_RUNON>
705 [[nodiscard]] T sum (int comp, int numcomp = 1) const noexcept;
707 template <RunOn run_on AMREX_DEFAULT_RUNON>
708 [[nodiscard]] T sum (const Box& subbox, int comp, int numcomp = 1) const noexcept;
709
711 template <RunOn run_on AMREX_DEFAULT_RUNON>
712 BaseFab<T>& invert (T const& r, const Box& b, int comp=0, int numcomp=1) noexcept;
714 template <RunOn run_on AMREX_DEFAULT_RUNON>
715 BaseFab<T>& invert (T const& r, int comp, int numcomp=1) noexcept;
716
718 template <RunOn run_on AMREX_DEFAULT_RUNON>
719 BaseFab<T>& negate (const Box& b, int comp=0, int numcomp=1) noexcept;
721 template <RunOn run_on AMREX_DEFAULT_RUNON>
722 BaseFab<T>& negate (int comp, int numcomp=1) noexcept;
723
725 template <RunOn run_on AMREX_DEFAULT_RUNON>
726 BaseFab<T>& plus (T const& r, const Box& b, int comp=0, int numcomp=1) noexcept;
727
729 template <RunOn run_on AMREX_DEFAULT_RUNON>
730 BaseFab<T>& plus (T const& r, int comp, int numcomp=1) noexcept;
731
737 template <RunOn run_on AMREX_DEFAULT_RUNON>
738 BaseFab<T>& plus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
744 template <RunOn run_on AMREX_DEFAULT_RUNON>
745 BaseFab<T>& plus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp=1) noexcept;
750 template <RunOn run_on AMREX_DEFAULT_RUNON>
751 BaseFab<T>& plus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
752 int srccomp, int destcomp, int numcomp=1) noexcept;
753
755 template <RunOn run_on AMREX_DEFAULT_RUNON>
756 BaseFab<T>& atomicAdd (const BaseFab<T>& x) noexcept;
757
763 template <RunOn run_on AMREX_DEFAULT_RUNON>
764 BaseFab<T>& atomicAdd (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
770 template <RunOn run_on AMREX_DEFAULT_RUNON>
771 BaseFab<T>& atomicAdd (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
772 int numcomp=1) noexcept;
777 template <RunOn run_on AMREX_DEFAULT_RUNON>
778 BaseFab<T>& atomicAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
779 int srccomp, int destcomp, int numcomp=1) noexcept;
780
786 template <RunOn run_on AMREX_DEFAULT_RUNON>
787 BaseFab<T>& lockAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
788 int srccomp, int destcomp, int numcomp) noexcept;
789
791 template <RunOn run_on AMREX_DEFAULT_RUNON>
792 BaseFab<T>& saxpy (T a, const BaseFab<T>& x, const Box& srcbox, const Box& destbox,
793 int srccomp, int destcomp, int numcomp=1) noexcept;
795 template <RunOn run_on AMREX_DEFAULT_RUNON>
796 BaseFab<T>& saxpy (T a, const BaseFab<T>& x) noexcept;
797
799 template <RunOn run_on AMREX_DEFAULT_RUNON>
800 BaseFab<T>& xpay (T a, const BaseFab<T>& x, const Box& srcbox, const Box& destbox,
801 int srccomp, int destcomp, int numcomp=1) noexcept;
802
804 template <RunOn run_on AMREX_DEFAULT_RUNON>
805 BaseFab<T>& addproduct (const Box& destbox, int destcomp, int numcomp,
806 const BaseFab<T>& src1, int comp1,
807 const BaseFab<T>& src2, int comp2) noexcept;
808
814 template <RunOn run_on AMREX_DEFAULT_RUNON>
815 BaseFab<T>& minus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
821 template <RunOn run_on AMREX_DEFAULT_RUNON>
822 BaseFab<T>& minus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
823 int numcomp=1) noexcept;
828 template <RunOn run_on AMREX_DEFAULT_RUNON>
829 BaseFab<T>& minus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
830 int srccomp, int destcomp, int numcomp=1) noexcept;
831
833 template <RunOn run_on AMREX_DEFAULT_RUNON>
834 BaseFab<T>& mult (T const& r, int comp, int numcomp=1) noexcept;
838 template <RunOn run_on AMREX_DEFAULT_RUNON>
839 BaseFab<T>& mult (T const& r, const Box& b, int comp=0, int numcomp=1) noexcept;
840
846 template <RunOn run_on AMREX_DEFAULT_RUNON>
847 BaseFab<T>& mult (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
848
854 template <RunOn run_on AMREX_DEFAULT_RUNON>
855 BaseFab<T>& mult (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
856 int numcomp=1) noexcept;
857
862 template <RunOn run_on AMREX_DEFAULT_RUNON>
863 BaseFab<T>& mult (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
864 int srccomp, int destcomp, int numcomp=1) noexcept;
865
867 template <RunOn run_on AMREX_DEFAULT_RUNON>
868 BaseFab<T>& divide (T const& r, int comp, int numcomp=1) noexcept;
869
871 template <RunOn run_on AMREX_DEFAULT_RUNON>
872 BaseFab<T>& divide (T const& r, const Box& b, int comp=0, int numcomp=1) noexcept;
873
880 template <RunOn run_on AMREX_DEFAULT_RUNON>
881 BaseFab<T>& divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
887 template <RunOn run_on AMREX_DEFAULT_RUNON>
888 BaseFab<T>& divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
889 int numcomp=1) noexcept;
894 template <RunOn run_on AMREX_DEFAULT_RUNON>
895 BaseFab<T>& divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
896 int srccomp, int destcomp, int numcomp=1) noexcept;
900 template <RunOn run_on AMREX_DEFAULT_RUNON>
901 BaseFab<T>& protected_divide (const BaseFab<T>& src) noexcept;
902
910 template <RunOn run_on AMREX_DEFAULT_RUNON>
911 BaseFab<T>& protected_divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp=1) noexcept;
912
919 template <RunOn run_on AMREX_DEFAULT_RUNON>
920 BaseFab<T>& protected_divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
921 int numcomp=1) noexcept;
922
928 template <RunOn run_on AMREX_DEFAULT_RUNON>
929 BaseFab<T>& protected_divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
930 int srccomp, int destcomp, int numcomp=1) noexcept;
931
942 template <RunOn run_on AMREX_DEFAULT_RUNON>
943 BaseFab<T>& linInterp (const BaseFab<T>& f1, const Box& b1, int comp1,
944 const BaseFab<T>& f2, const Box& b2, int comp2,
945 Real t1, Real t2, Real t,
946 const Box& b, int comp, int numcomp = 1) noexcept;
947
949 template <RunOn run_on AMREX_DEFAULT_RUNON>
950 BaseFab<T>& linInterp (const BaseFab<T>& f1, int comp1,
951 const BaseFab<T>& f2, int comp2,
952 Real t1, Real t2, Real t,
953 const Box& b, int comp, int numcomp = 1) noexcept;
954
964 template <RunOn run_on AMREX_DEFAULT_RUNON>
965 BaseFab<T>& linComb (const BaseFab<T>& f1, const Box& b1, int comp1,
966 const BaseFab<T>& f2, const Box& b2, int comp2,
967 Real alpha, Real beta, const Box& b,
968 int comp, int numcomp = 1) noexcept;
969
971 template <RunOn run_on AMREX_DEFAULT_RUNON>
972 [[nodiscard]] T dot (const Box& xbx, int xcomp, const BaseFab<T>& y, const Box& ybx, int ycomp,
973 int numcomp = 1) const noexcept;
974
975 template <RunOn run_on AMREX_DEFAULT_RUNON>
976 [[nodiscard]] T dotmask (const BaseFab<int>& mask, const Box& xbx, int xcomp,
977 const BaseFab<T>& y, const Box& ybx, int ycomp,
978 int numcomp) const noexcept;
979
981 void SetBoxType (const IndexType& typ) noexcept { this->domain.setType(typ); }
982
983 //
984 // New interfaces
985 //
986
988 template <RunOn run_on AMREX_DEFAULT_RUNON>
989 void setVal (T const& val) noexcept;
990 //
992 template <RunOn run_on AMREX_DEFAULT_RUNON>
993 void setVal (T const& x, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
994
995 template <RunOn run_on AMREX_DEFAULT_RUNON>
996 void setValIf (T const& val, const BaseFab<int>& mask) noexcept;
997 //
999 template <RunOn run_on AMREX_DEFAULT_RUNON>
1000 void setValIf (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept;
1001
1002 template <RunOn run_on AMREX_DEFAULT_RUNON>
1003 void setValIfNot (T const& val, const BaseFab<int>& mask) noexcept;
1004 //
1006 template <RunOn run_on AMREX_DEFAULT_RUNON>
1007 void setValIfNot (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept;
1008
1010 template <RunOn run_on AMREX_DEFAULT_RUNON>
1011 void setComplement (T const& x, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
1012
1018 template <RunOn run_on AMREX_DEFAULT_RUNON>
1019 BaseFab<T>& copy (const BaseFab<T>& src) noexcept;
1020 //
1022 template <RunOn run_on AMREX_DEFAULT_RUNON>
1023 BaseFab<T>& copy (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
1024
1026 template <RunOn run_on AMREX_DEFAULT_RUNON>
1027 BaseFab<T>& plus (T const& val) noexcept;
1028 //
1029 template <RunOn run_on AMREX_DEFAULT_RUNON>
1030 BaseFab<T>& operator+= (T const& val) noexcept;
1031 //
1033 template <RunOn run_on AMREX_DEFAULT_RUNON>
1034 BaseFab<T>& plus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
1040 template <RunOn run_on AMREX_DEFAULT_RUNON>
1041 BaseFab<T>& plus (const BaseFab<T>& src) noexcept;
1042 //
1043 template <RunOn run_on AMREX_DEFAULT_RUNON>
1044 BaseFab<T>& operator+= (const BaseFab<T>& src) noexcept;
1045 //
1047 template <RunOn run_on AMREX_DEFAULT_RUNON>
1048 BaseFab<T>& plus (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
1049
1051 template <RunOn run_on AMREX_DEFAULT_RUNON>
1052 BaseFab<T>& minus (T const& val) noexcept;
1053 //
1054 template <RunOn run_on AMREX_DEFAULT_RUNON>
1055 BaseFab<T>& operator-= (T const& val) noexcept;
1056 //
1058 template <RunOn run_on AMREX_DEFAULT_RUNON>
1059 BaseFab<T>& minus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
1065 template <RunOn run_on AMREX_DEFAULT_RUNON>
1066 BaseFab<T>& minus (const BaseFab<T>& src) noexcept;
1067 //
1068 template <RunOn run_on AMREX_DEFAULT_RUNON>
1069 BaseFab<T>& operator-= (const BaseFab<T>& src) noexcept;
1070 //
1072 template <RunOn run_on AMREX_DEFAULT_RUNON>
1073 BaseFab<T>& minus (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
1074
1076 template <RunOn run_on AMREX_DEFAULT_RUNON>
1077 BaseFab<T>& mult (T const& val) noexcept;
1078 //
1079 template <RunOn run_on AMREX_DEFAULT_RUNON>
1080 BaseFab<T>& operator*= (T const& val) noexcept;
1081 //
1083 template <RunOn run_on AMREX_DEFAULT_RUNON>
1084 BaseFab<T>& mult (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
1090 template <RunOn run_on AMREX_DEFAULT_RUNON>
1091 BaseFab<T>& mult (const BaseFab<T>& src) noexcept;
1092 //
1093 template <RunOn run_on AMREX_DEFAULT_RUNON>
1094 BaseFab<T>& operator*= (const BaseFab<T>& src) noexcept;
1095 //
1097 template <RunOn run_on AMREX_DEFAULT_RUNON>
1098 BaseFab<T>& mult (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
1099
1101 template <RunOn run_on AMREX_DEFAULT_RUNON>
1102 BaseFab<T>& divide (T const& val) noexcept;
1103 //
1104 template <RunOn run_on AMREX_DEFAULT_RUNON>
1105 BaseFab<T>& operator/= (T const& val) noexcept;
1106 //
1108 template <RunOn run_on AMREX_DEFAULT_RUNON>
1109 BaseFab<T>& divide (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept;
1115 template <RunOn run_on AMREX_DEFAULT_RUNON>
1116 BaseFab<T>& divide (const BaseFab<T>& src) noexcept;
1117 //
1118 template <RunOn run_on AMREX_DEFAULT_RUNON>
1119 BaseFab<T>& operator/= (const BaseFab<T>& src) noexcept;
1120 //
1122 template <RunOn run_on AMREX_DEFAULT_RUNON>
1123 BaseFab<T>& divide (const BaseFab<T>& src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept;
1124
1126 template <RunOn run_on AMREX_DEFAULT_RUNON>
1127 BaseFab<T>& negate () noexcept;
1128 //
1129 template <RunOn run_on AMREX_DEFAULT_RUNON>
1130 BaseFab<T>& negate (const Box& bx, DestComp dcomp, NumComps ncomp) noexcept;
1131
1133 template <RunOn run_on AMREX_DEFAULT_RUNON>
1134 BaseFab<T>& invert (T const& r) noexcept;
1135 //
1136 template <RunOn run_on AMREX_DEFAULT_RUNON>
1137 BaseFab<T>& invert (T const& r, const Box& bx, DestComp dcomp, NumComps ncomp) noexcept;
1138
1140 template <RunOn run_on AMREX_DEFAULT_RUNON>
1141 [[nodiscard]] T sum (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept;
1142
1144 template <RunOn run_on AMREX_DEFAULT_RUNON>
1145 [[nodiscard]] T dot (const BaseFab<T>& src, const Box& bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept;
1146
1148 template <RunOn run_on AMREX_DEFAULT_RUNON>
1149 [[nodiscard]] T dot (const Box& bx, int destcomp, int numcomp) const noexcept;
1150
1152 template <RunOn run_on AMREX_DEFAULT_RUNON>
1153 [[nodiscard]] T dot (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept;
1154
1156 template <RunOn run_on AMREX_DEFAULT_RUNON>
1157 [[nodiscard]] T dotmask (const BaseFab<T>& src, const Box& bx, const BaseFab<int>& mask,
1158 SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept;
1159
1160protected:
1162 void define ();
1163
1164 T* dptr = nullptr;
1166 int nvar = 0;
1168 bool ptr_owner = false;
1169 bool shared_memory = false;
1170#ifdef AMREX_USE_GPU
1172#endif
1173};
1174
1175template <class T>
1177T*
1178BaseFab<T>::dataPtr (const IntVect& p, int n) noexcept
1179{
1180 AMREX_ASSERT(n >= 0);
1181 AMREX_ASSERT(n < this->nvar);
1182 AMREX_ASSERT(!(this->dptr == nullptr));
1183 AMREX_ASSERT(this->domain.contains(p));
1184
1185 return this->dptr + (this->domain.index(p)+n*this->domain.numPts());
1186}
1187
1188template <class T>
1190const T*
1191BaseFab<T>::dataPtr (const IntVect& p, int n) const noexcept
1192{
1193 AMREX_ASSERT(n >= 0);
1194 AMREX_ASSERT(n < this->nvar);
1195 AMREX_ASSERT(!(this->dptr == nullptr));
1196 AMREX_ASSERT(this->domain.contains(p));
1197
1198 return this->dptr + (this->domain.index(p)+n*this->domain.numPts());
1199}
1200
1201template <class T>
1202void
1204{
1205#ifdef AMREX_USE_GPU
1206 if (this->arena()->isManaged()) {
1207#if defined(AMREX_USE_SYCL)
1208 // xxxxx SYCL todo: prefetchToHost
1209 // std::size_t s = sizeof(T)*this->nvar*this->domain.numPts();
1210 // auto& q = Gpu::Device::streamQueue();
1211 // q.submit([&] (sycl::handler& h) { h.prefetch(this->dptr, s); });
1212#elif defined(AMREX_USE_CUDA) && !defined(_WIN32)
1213 if (Gpu::Device::devicePropMajor() >= 6) {
1214 std::size_t s = sizeof(T)*this->nvar*this->domain.numPts();
1215#if defined(__CUDACC__) && (__CUDACC_VER_MAJOR__ >= 13)
1216 cudaMemLocation location = {};
1217 location.type = cudaMemLocationTypeDevice;
1218 location.id = cudaCpuDeviceId;
1219 AMREX_CUDA_SAFE_CALL(cudaMemPrefetchAsync(this->dptr, s, location, 0,
1220 Gpu::gpuStream()));
1221#else
1222 AMREX_CUDA_SAFE_CALL(cudaMemPrefetchAsync(this->dptr, s,
1223 cudaCpuDeviceId,
1224 Gpu::gpuStream()));
1225#endif
1226 }
1227#elif defined(AMREX_USE_HIP)
1228 // xxxxx HIP FIX HERE after managed memory is supported
1229#endif
1230 }
1231#endif
1232}
1233
1234template <class T>
1235void
1237{
1238#ifdef AMREX_USE_GPU
1239 if (this->arena()->isManaged()) {
1240#if defined(AMREX_USE_SYCL)
1241 std::size_t s = sizeof(T)*this->nvar*this->domain.numPts();
1242 auto& q = Gpu::Device::streamQueue();
1243 q.submit([&] (sycl::handler& h) { h.prefetch(this->dptr, s); });
1244#elif defined(AMREX_USE_CUDA) && !defined(_WIN32)
1245 if (Gpu::Device::devicePropMajor() >= 6) {
1246 std::size_t s = sizeof(T)*this->nvar*this->domain.numPts();
1247#if defined(__CUDACC__) && (__CUDACC_VER_MAJOR__ >= 13)
1248 cudaMemLocation location = {};
1249 location.type = cudaMemLocationTypeDevice;
1250 location.id = Gpu::Device::deviceId();
1251 AMREX_CUDA_SAFE_CALL(cudaMemPrefetchAsync(this->dptr, s, location, 0,
1252 Gpu::gpuStream()));
1253#else
1254 AMREX_CUDA_SAFE_CALL(cudaMemPrefetchAsync(this->dptr, s,
1256 Gpu::gpuStream()));
1257#endif
1258 }
1259#elif defined(AMREX_USE_HIP)
1260 // xxxxx HIP FIX HERE after managed memory is supported
1261#endif
1262 }
1263#endif
1264}
1265
1266template <class T>
1268T&
1269BaseFab<T>::operator() (const IntVect& p, int n) noexcept
1270{
1271 AMREX_ASSERT(n >= 0);
1272 AMREX_ASSERT(n < this->nvar);
1273 AMREX_ASSERT(!(this->dptr == nullptr));
1274 AMREX_ASSERT(this->domain.contains(p));
1275
1276 return this->dptr[this->domain.index(p)+n*this->domain.numPts()];
1277}
1278
1279template <class T>
1281T&
1283{
1284 AMREX_ASSERT(!(this->dptr == nullptr));
1285 AMREX_ASSERT(this->domain.contains(p));
1286
1287 return this->dptr[this->domain.index(p)];
1288}
1289
1290template <class T>
1292const T&
1293BaseFab<T>::operator() (const IntVect& p, int n) const noexcept
1294{
1295 AMREX_ASSERT(n >= 0);
1296 AMREX_ASSERT(n < this->nvar);
1297 AMREX_ASSERT(!(this->dptr == nullptr));
1298 AMREX_ASSERT(this->domain.contains(p));
1299
1300 return this->dptr[this->domain.index(p)+n*this->domain.numPts()];
1301}
1302
1303template <class T>
1305const T&
1306BaseFab<T>::operator() (const IntVect& p) const noexcept
1307{
1308 AMREX_ASSERT(!(this->dptr == nullptr));
1309 AMREX_ASSERT(this->domain.contains(p));
1310
1311 return this->dptr[this->domain.index(p)];
1312}
1313
1314template <class T>
1315void
1317 const IntVect& pos,
1318 int n,
1319 int numcomp) const noexcept
1320{
1321 const int loc = this->domain.index(pos);
1322 const Long sz = this->domain.numPts();
1323
1324 AMREX_ASSERT(!(this->dptr == nullptr));
1325 AMREX_ASSERT(n >= 0 && n + numcomp <= this->nvar);
1326
1327 for (int k = 0; k < numcomp; k++) {
1328 data[k] = this->dptr[loc+(n+k)*sz];
1329 }
1330}
1331
1332template <class T>
1333void
1335 const IntVect& pos) const noexcept
1336{
1337 getVal(data,pos,0,this->nvar);
1338}
1339
1340template <class T>
1342BaseFab<T>::shift (const IntVect& v) noexcept
1343{
1344 this->domain += v;
1345 return *this;
1346}
1347
1348template <class T>
1350BaseFab<T>::shift (int idir, int n_cell) noexcept
1351{
1352 this->domain.shift(idir,n_cell);
1353 return *this;
1354}
1355
1356template <class T>
1357BaseFab<T> &
1359{
1360 this->domain.shiftHalf(v);
1361 return *this;
1362}
1363
1364template <class T>
1365BaseFab<T> &
1366BaseFab<T>::shiftHalf (int idir, int n_cell) noexcept
1367{
1368 this->domain.shiftHalf(idir,n_cell);
1369 return *this;
1370}
1371
1372template <class T>
1373template <RunOn run_on, class U,
1374 std::enable_if_t<std::is_same_v<U,float> || std::is_same_v<U,double>, int> FOO>
1375void
1377{
1378 amrex::fill_snan<run_on>(this->dptr, this->truesize);
1379}
1380
1381template <class T>
1382template <RunOn run_on>
1383void
1384BaseFab<T>::setVal (T const& x, const Box& bx, int n) noexcept
1385{
1386 this->setVal<run_on>(x, bx, DestComp{n}, NumComps{1});
1387}
1388
1389template <class T>
1390template <RunOn run_on>
1391void
1392BaseFab<T>::setVal (T const& x, int n) noexcept
1393{
1394 this->setVal<run_on>(x, this->domain, DestComp{n}, NumComps{1});
1395}
1396
1397template <class T>
1398template <RunOn run_on>
1399void
1400BaseFab<T>::setVal (T const& x, const Box& bx, int dcomp, int ncomp) noexcept
1401{
1402 this->setVal<run_on>(x, bx, DestComp{dcomp}, NumComps{ncomp});
1403}
1404
1405template <class T>
1406template <RunOn run_on>
1407void
1408BaseFab<T>::setValIfNot (T const& val, const Box& bx, const BaseFab<int>& mask, int ns, int num) noexcept
1409{
1410 this->setValIfNot<run_on>(val, bx, mask, DestComp{ns}, NumComps{num});
1411}
1412
1413template <class T>
1414template <RunOn run_on>
1416BaseFab<T>::copy (const BaseFab<T>& src, const Box& srcbox, int srccomp,
1417 const Box& destbox, int destcomp, int numcomp) noexcept
1418{
1419 AMREX_ASSERT(destbox.ok());
1420 AMREX_ASSERT(srcbox.sameSize(destbox));
1421 AMREX_ASSERT(src.box().contains(srcbox));
1422 AMREX_ASSERT(this->domain.contains(destbox));
1423 AMREX_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
1424 AMREX_ASSERT(destcomp >= 0 && destcomp+numcomp <= this->nvar);
1425
1426 Array4<T> const& d = this->array();
1427 Array4<T const> const& s = src.const_array();
1428 const auto dlo = amrex::lbound(destbox);
1429 const auto slo = amrex::lbound(srcbox);
1430 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
1431
1432 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
1433 {
1434 d(i,j,k,n+destcomp) = s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
1435 });
1436
1437 return *this;
1438}
1439
1440template <class T>
1441template <RunOn run_on>
1443BaseFab<T>::copy (const BaseFab<T>& src, const Box& destbox) noexcept
1444{
1445 return this->copy<run_on>(src, destbox, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
1446}
1447
1448template <class T>
1449template <RunOn run_on>
1451BaseFab<T>::copy (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
1452{
1453 return copy<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
1454}
1455
1456template <class T>
1457void
1459{
1460 AMREX_ASSERT(this->dptr == nullptr);
1461 AMREX_ASSERT(this->domain.numPts() > 0);
1462 AMREX_ASSERT(this->nvar >= 0);
1463 if (this->nvar == 0) { return; }
1464 AMREX_ASSERT(std::numeric_limits<Long>::max()/this->nvar > this->domain.numPts());
1465
1466 this->truesize = this->nvar*this->domain.numPts();
1467 this->ptr_owner = true;
1468 this->dptr = static_cast<T*>(this->alloc(this->truesize*sizeof(T)));
1469#ifdef AMREX_USE_GPU
1470 this->alloc_stream = Gpu::gpuStream();
1471#endif
1472
1473 placementNew(this->dptr, this->truesize);
1474
1475 amrex::update_fab_stats(this->domain.numPts(), this->truesize, sizeof(T));
1476
1477 if constexpr (std::is_same_v<T,float> || std::is_same_v<T,double>) {
1478 if (amrex::InitSNaN() && this->truesize > 0) {
1479#ifdef AMREX_USE_GPU
1480 if (Gpu::inLaunchRegion() && arena()->isDeviceAccessible()) {
1481 this->template fill_snan<RunOn::Device>();
1483 } else
1484#endif
1485 {
1486 this->template fill_snan<RunOn::Host>();
1487 }
1488 }
1489 }
1490}
1491
1492template <class T>
1494 : DataAllocator{ar}
1495{}
1496
1497template <class T>
1498BaseFab<T>::BaseFab (const Box& bx, int n, Arena* ar)
1499 : DataAllocator{ar}, domain(bx), nvar(n)
1500{
1501 define();
1502}
1503
1504template <class T>
1505BaseFab<T>::BaseFab (const Box& bx, int n, bool alloc, bool shared, Arena* ar)
1506 : DataAllocator{ar}, domain(bx), nvar(n), shared_memory(shared)
1507{
1508 if (!this->shared_memory && alloc) { define(); }
1509}
1510
1511template <class T>
1512BaseFab<T>::BaseFab (const BaseFab<T>& rhs, MakeType make_type, int scomp, int ncomp)
1513 : DataAllocator{rhs.arena()},
1514 dptr(const_cast<T*>(rhs.dataPtr(scomp))),
1515 domain(rhs.domain), nvar(ncomp),
1516 truesize(ncomp*rhs.domain.numPts())
1517{
1518 AMREX_ASSERT(scomp+ncomp <= rhs.nComp());
1519 if (make_type == amrex::make_deep_copy)
1520 {
1521 this->dptr = nullptr;
1522 define();
1523 this->copy<RunOn::Device>(rhs, this->domain, scomp, this->domain, 0, ncomp);
1524 } else if (make_type == amrex::make_alias) {
1525 ; // nothing to do
1526 } else {
1527 amrex::Abort("BaseFab: unknown MakeType");
1528 }
1529}
1530
1531template<class T>
1532BaseFab<T>::BaseFab (const Box& bx, int ncomp, T* p)
1533 : dptr(p), domain(bx), nvar(ncomp), truesize(bx.numPts()*ncomp)
1534{
1535}
1536
1537template<class T>
1538BaseFab<T>::BaseFab (const Box& bx, int ncomp, T const* p)
1539 : dptr(const_cast<T*>(p)), domain(bx), nvar(ncomp), truesize(bx.numPts()*ncomp)
1540{
1541}
1542
1543template<class T>
1545 : dptr(a.p),
1546 domain(IntVect(AMREX_D_DECL(a.begin[0],a.begin[1],a.begin[2])),
1547 IntVect(AMREX_D_DECL(a.end[0]-1,a.end[1]-1,a.end[2]-1))),
1548 nvar(a.nComp()), truesize(a.size())
1549{}
1550
1551template<class T>
1553 : dptr(a.p),
1554 domain(IntVect(AMREX_D_DECL(a.begin[0],a.begin[1],a.begin[2])),
1555 IntVect(AMREX_D_DECL(a.end[0]-1,a.end[1]-1,a.end[2]-1)), t),
1556 nvar(a.nComp()), truesize(a.size())
1557{}
1558
1559template<class T>
1561 : dptr(const_cast<T*>(a.p)),
1562 domain(IntVect(AMREX_D_DECL(a.begin[0],a.begin[1],a.begin[2])),
1563 IntVect(AMREX_D_DECL(a.end[0]-1,a.end[1]-1,a.end[2]-1))),
1564 nvar(a.nComp()), truesize(a.size())
1565{}
1566
1567template<class T>
1569 : dptr(const_cast<T*>(a.p)),
1570 domain(IntVect(AMREX_D_DECL(a.begin[0],a.begin[1],a.begin[2])),
1571 IntVect(AMREX_D_DECL(a.end[0]-1,a.end[1]-1,a.end[2]-1)), t),
1572 nvar(a.nComp()), truesize(a.size())
1573{}
1574
1575template <class T>
1577{
1578 clear();
1579}
1580
1581template <class T>
1583 : DataAllocator{rhs.arena()},
1584 dptr(rhs.dptr), domain(rhs.domain),
1585 nvar(rhs.nvar), truesize(rhs.truesize),
1586 ptr_owner(rhs.ptr_owner), shared_memory(rhs.shared_memory)
1587#ifdef AMREX_USE_GPU
1588 , alloc_stream(rhs.alloc_stream)
1589#endif
1590{
1591 rhs.dptr = nullptr;
1592 rhs.ptr_owner = false;
1593}
1594
1595template <class T>
1596BaseFab<T>&
1598{
1599 if (this != &rhs) {
1600 clear();
1601 DataAllocator::operator=(rhs);
1602 dptr = rhs.dptr;
1603 domain = rhs.domain;
1604 nvar = rhs.nvar;
1605 truesize = rhs.truesize;
1606 ptr_owner = rhs.ptr_owner;
1607 shared_memory = rhs.shared_memory;
1608#ifdef AMREX_USE_GPU
1609 alloc_stream = rhs.alloc_stream;
1610#endif
1611
1612 rhs.dptr = nullptr;
1613 rhs.ptr_owner = false;
1614 }
1615 return *this;
1616}
1617
1618template <class T>
1619template <RunOn run_on>
1621BaseFab<T>::operator= (T const& t) noexcept
1622{
1623 setVal<run_on>(t);
1624 return *this;
1625}
1626
1627template <class T>
1628void
1629BaseFab<T>::resize (const Box& b, int n, Arena* ar)
1630{
1631 this->nvar = n;
1632 this->domain = b;
1633
1634 if (ar == nullptr) {
1635 ar = m_arena;
1636 }
1637
1638 if (arena() != DataAllocator(ar).arena()) {
1639 clear();
1640 m_arena = ar;
1641 define();
1642 }
1643 else if (this->dptr == nullptr || !this->ptr_owner)
1644 {
1645 if (this->shared_memory) {
1646 amrex::Abort("BaseFab::resize: BaseFab in shared memory cannot increase size");
1647 }
1648
1649 this->dptr = nullptr;
1650 define();
1651 }
1652 else if (this->nvar*this->domain.numPts() > this->truesize
1653#ifdef AMREX_USE_GPU
1654 || (arena()->isStreamOrderedArena() && alloc_stream != Gpu::gpuStream())
1655#endif
1656 )
1657 {
1658 if (this->shared_memory) {
1659 amrex::Abort("BaseFab::resize: BaseFab in shared memory cannot increase size");
1660 }
1661
1662 clear();
1663
1664 define();
1665 }
1666}
1667
1668template <class T>
1669template <class U, std::enable_if_t<std::is_trivially_destructible_v<U>,int>>
1670Elixir
1672{
1673 bool o;
1674 if (Gpu::inLaunchRegion()) {
1675 o = this->ptr_owner;
1676 this->ptr_owner = false;
1677 if (o && this->dptr) {
1678 if (this->nvar > 1) {
1679 amrex::update_fab_stats(-this->truesize/this->nvar, -this->truesize, sizeof(T));
1680 } else {
1681 amrex::update_fab_stats(0, -this->truesize, sizeof(T));
1682 }
1683 }
1684 } else {
1685 o = false;
1686 }
1687 return Elixir((o ? this->dptr : nullptr), this->arena());
1688}
1689
1690template <class T>
1691void
1693{
1694 if (this->dptr)
1695 {
1696 //
1697 // Call T::~T() on the to-be-destroyed memory.
1698 //
1699 if (this->ptr_owner)
1700 {
1701 if (this->shared_memory)
1702 {
1703 amrex::Abort("BaseFab::clear: BaseFab cannot be owner of shared memory");
1704 }
1705
1706 placementDelete(this->dptr, this->truesize);
1707
1708#ifdef AMREX_USE_GPU
1709 auto current_stream = Gpu::Device::gpuStream();
1710 Gpu::Device::setStream(alloc_stream);
1711#endif
1712 this->free(this->dptr);
1713#ifdef AMREX_USE_GPU
1714 Gpu::Device::setStream(current_stream);
1715#endif
1716
1717 if (this->nvar > 1) {
1718 amrex::update_fab_stats(-this->truesize/this->nvar, -this->truesize, sizeof(T));
1719 } else {
1720 amrex::update_fab_stats(0, -this->truesize, sizeof(T));
1721 }
1722 }
1723
1724 this->dptr = nullptr;
1725 this->truesize = 0;
1726 }
1727}
1728
1729template <class T>
1730std::unique_ptr<T,DataDeleter>
1732{
1733 std::unique_ptr<T,DataDeleter> r(nullptr, DataDeleter{this->arena()});
1734 if (this->dptr && this->ptr_owner) {
1735 r.reset(this->dptr);
1736 this->ptr_owner = false;
1737 if (this->nvar > 1) {
1738 amrex::update_fab_stats(-this->truesize/this->nvar, -this->truesize, sizeof(T));
1739 } else {
1740 amrex::update_fab_stats(0, -this->truesize, sizeof(T));
1741 }
1742 }
1743 return r;
1744}
1745
1746template <class T>
1747template <RunOn run_on>
1748std::size_t
1750 int srccomp,
1751 int numcomp,
1752 void* dst) const noexcept
1753{
1754 BL_ASSERT(box().contains(srcbox));
1755 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= nComp());
1756
1757 if (srcbox.ok())
1758 {
1759 Array4<T> d(static_cast<T*>(dst),amrex::begin(srcbox),amrex::end(srcbox),numcomp);
1760 Array4<T const> const& s = this->const_array();
1761 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, srcbox, numcomp, i, j, k, n,
1762 {
1763 d(i,j,k,n) = s(i,j,k,n+srccomp);
1764 });
1765 return sizeof(T)*d.size();
1766 }
1767 else
1768 {
1769 return 0;
1770 }
1771}
1772
1773template <class T>
1774template <RunOn run_on, typename BUF>
1775std::size_t
1777 int dstcomp,
1778 int numcomp,
1779 const void* src) noexcept
1780{
1781 BL_ASSERT(box().contains(dstbox));
1782 BL_ASSERT(dstcomp >= 0 && dstcomp+numcomp <= nComp());
1783
1784 if (dstbox.ok())
1785 {
1786 Array4<BUF const> s(static_cast<BUF const*>(src), amrex::begin(dstbox),
1787 amrex::end(dstbox), numcomp);
1788 Array4<T> const& d = this->array();
1789 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, dstbox, numcomp, i, j, k, n,
1790 {
1791 d(i,j,k,n+dstcomp) = static_cast<T>(s(i,j,k,n));
1792 });
1793 return sizeof(BUF)*s.size();
1794 }
1795 else
1796 {
1797 return 0;
1798 }
1799}
1800
1801template <class T>
1802template <RunOn run_on, typename BUF>
1803std::size_t
1805 int dstcomp,
1806 int numcomp,
1807 const void* src) noexcept
1808{
1809 BL_ASSERT(box().contains(dstbox));
1810 BL_ASSERT(dstcomp >= 0 && dstcomp+numcomp <= nComp());
1811
1812 if (dstbox.ok())
1813 {
1814 Array4<BUF const> s(static_cast<BUF const*>(src), amrex::begin(dstbox),
1815 amrex::end(dstbox), numcomp);
1816 Array4<T> const& d = this->array();
1817 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, dstbox, numcomp, i, j, k, n,
1818 {
1819 d(i,j,k,n+dstcomp) += static_cast<T>(s(i,j,k,n));
1820 });
1821 return sizeof(BUF)*s.size();
1822 }
1823 else
1824 {
1825 return 0;
1826 }
1827}
1828
1829template <class T>
1830template <RunOn run_on>
1831void
1832BaseFab<T>::setComplement (T const& x, const Box& b, int ns, int num) noexcept
1833{
1834 this->setComplement<run_on>(x, b, DestComp{ns}, NumComps{num});
1835}
1836
1837template <class T>
1838template <RunOn run_on>
1839void
1841{
1842 this->abs<run_on>(this->domain,0,this->nvar);
1843}
1844
1845template <class T>
1846template <RunOn run_on>
1847void
1848BaseFab<T>::abs (int comp, int numcomp) noexcept
1849{
1850 this->abs<run_on>(this->domain,comp,numcomp);
1851}
1852
1853template <class T>
1854template <RunOn run_on>
1855void
1856BaseFab<T>::abs (const Box& subbox, int comp, int numcomp) noexcept
1857{
1858 Array4<T> const& a = this->array();
1859 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, subbox, numcomp, i, j, k, n,
1860 {
1861 a(i,j,k,n+comp) = std::abs(a(i,j,k,n+comp));
1862 });
1863}
1864
1865template <class T>
1866template <RunOn run_on>
1867Real
1869 int scomp, int ncomp) const noexcept
1870{
1871 BL_ASSERT(this->domain.contains(subbox));
1872 BL_ASSERT(scomp >= 0 && scomp + ncomp <= this->nvar);
1873
1874 Array4<T const> const& a = this->const_array();
1875 Array4<int const> const& m = mask.const_array();
1876 Real r = 0.0;
1877#ifdef AMREX_USE_GPU
1878 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
1879 ReduceOps<ReduceOpMax> reduce_op;
1880 ReduceData<Real> reduce_data(reduce_op);
1881 using ReduceTuple = ReduceData<Real>::Type;
1882 reduce_op.eval(subbox, reduce_data,
1883 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
1884 {
1885 Real t = 0.0;
1886 if (m(i,j,k)) {
1887 for (int n = 0; n < ncomp; ++n) {
1888 t = amrex::max(t,static_cast<Real>(std::abs(a(i,j,k,n+scomp))));
1889 }
1890 }
1891 return {t};
1892 });
1893 ReduceTuple hv = reduce_data.value(reduce_op);
1894 r = amrex::get<0>(hv);
1895 } else
1896#endif
1897 {
1898 amrex::LoopOnCpu(subbox, ncomp, [=,&r] (int i, int j, int k, int n) noexcept
1899 {
1900 if (m(i,j,k)) {
1901 Real t = static_cast<Real>(std::abs(a(i,j,k,n+scomp)));
1902 r = amrex::max(r,t);
1903 }
1904 });
1905 }
1906 return r;
1907}
1908
1909template <class T>
1910template <RunOn run_on>
1911Real
1912BaseFab<T>::norm (int p, int comp, int numcomp) const noexcept
1913{
1914 return norm<run_on>(this->domain,p,comp,numcomp);
1915}
1916
1917template <class T>
1918template <RunOn run_on>
1919Real
1920BaseFab<T>::norm (const Box& subbox, int p, int comp, int numcomp) const noexcept
1921{
1922 BL_ASSERT(this->domain.contains(subbox));
1923 BL_ASSERT(comp >= 0 && comp + numcomp <= this->nvar);
1924
1925 Array4<T const> const& a = this->const_array();
1926 Real nrm = 0.;
1927#ifdef AMREX_USE_GPU
1928 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
1929 if (p == 0) {
1930 ReduceOps<ReduceOpMax> reduce_op;
1931 ReduceData<Real> reduce_data(reduce_op);
1932 using ReduceTuple = ReduceData<Real>::Type;
1933 reduce_op.eval(subbox, reduce_data,
1934 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
1935 {
1936 Real t = 0.0;
1937 for (int n = 0; n < numcomp; ++n) {
1938 t = amrex::max(t, static_cast<Real>(std::abs(a(i,j,k,n+comp))));
1939 }
1940 return {t};
1941 });
1942 ReduceTuple hv = reduce_data.value(reduce_op);
1943 nrm = amrex::get<0>(hv);
1944 } else if (p == 1) {
1945 ReduceOps<ReduceOpSum> reduce_op;
1946 ReduceData<Real> reduce_data(reduce_op);
1947 using ReduceTuple = ReduceData<Real>::Type;
1948 reduce_op.eval(subbox, reduce_data,
1949 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
1950 {
1951 Real t = 0.0;
1952 for (int n = 0; n < numcomp; ++n) {
1953 t += static_cast<Real>(std::abs(a(i,j,k,n+comp)));
1954 }
1955 return {t};
1956 });
1957 ReduceTuple hv = reduce_data.value(reduce_op);
1958 nrm = amrex::get<0>(hv);
1959 } else if (p == 2) {
1960 ReduceOps<ReduceOpSum> reduce_op;
1961 ReduceData<Real> reduce_data(reduce_op);
1962 using ReduceTuple = ReduceData<Real>::Type;
1963 reduce_op.eval(subbox, reduce_data,
1964 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
1965 {
1966 Real t = 0.0;
1967 for (int n = 0; n < numcomp; ++n) {
1968 t += static_cast<Real>(a(i,j,k,n+comp)*a(i,j,k,n+comp));
1969 }
1970 return {t};
1971 });
1972 ReduceTuple hv = reduce_data.value(reduce_op);
1973 nrm = amrex::get<0>(hv);
1974 } else {
1975 amrex::Error("BaseFab<T>::norm: wrong p");
1976 }
1977 } else
1978#endif
1979 {
1980 if (p == 0) {
1981 amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
1982 {
1983 Real t = static_cast<Real>(std::abs(a(i,j,k,n+comp)));
1984 nrm = amrex::max(nrm,t);
1985 });
1986 } else if (p == 1) {
1987 amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
1988 {
1989 nrm += std::abs(a(i,j,k,n+comp));
1990 });
1991 } else if (p == 2) {
1992 amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
1993 {
1994 nrm += a(i,j,k,n+comp)*a(i,j,k,n+comp);
1995 });
1996 } else {
1997 amrex::Error("BaseFab<T>::norm: wrong p");
1998 }
1999 }
2000
2001 return nrm;
2002}
2003
2004template <class T>
2005template <RunOn run_on>
2006T
2007BaseFab<T>::min (int comp) const noexcept
2008{
2009 return this->min<run_on>(this->domain,comp);
2010}
2011
2012template <class T>
2013template <RunOn run_on>
2014T
2015BaseFab<T>::min (const Box& subbox, int comp) const noexcept
2016{
2017 Array4<T const> const& a = this->const_array(comp);
2018#ifdef AMREX_USE_GPU
2019 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2020 ReduceOps<ReduceOpMin> reduce_op;
2021 ReduceData<T> reduce_data(reduce_op);
2022 using ReduceTuple = typename decltype(reduce_data)::Type;
2023 reduce_op.eval(subbox, reduce_data,
2024 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2025 {
2026 return { a(i,j,k) };
2027 });
2028 ReduceTuple hv = reduce_data.value(reduce_op);
2029 return amrex::get<0>(hv);
2030 } else
2031#endif
2032 {
2033 T r = std::numeric_limits<T>::max();
2034 amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
2035 {
2036 r = amrex::min(r, a(i,j,k));
2037 });
2038 return r;
2039 }
2040}
2041
2042template <class T>
2043template <RunOn run_on>
2044T
2045BaseFab<T>::max (int comp) const noexcept
2046{
2047 return this->max<run_on>(this->domain,comp);
2048}
2049
2050template <class T>
2051template <RunOn run_on>
2052T
2053BaseFab<T>::max (const Box& subbox, int comp) const noexcept
2054{
2055 Array4<T const> const& a = this->const_array(comp);
2056#ifdef AMREX_USE_GPU
2057 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2058 ReduceOps<ReduceOpMax> reduce_op;
2059 ReduceData<T> reduce_data(reduce_op);
2060 using ReduceTuple = typename decltype(reduce_data)::Type;
2061 reduce_op.eval(subbox, reduce_data,
2062 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2063 {
2064 return { a(i,j,k) };
2065 });
2066 ReduceTuple hv = reduce_data.value(reduce_op);
2067 return amrex::get<0>(hv);
2068 } else
2069#endif
2070 {
2071 T r = std::numeric_limits<T>::lowest();
2072 amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
2073 {
2074 r = amrex::max(r, a(i,j,k));
2075 });
2076 return r;
2077 }
2078}
2079
2080template <class T>
2081template <RunOn run_on>
2082std::pair<T,T>
2083BaseFab<T>::minmax (int comp) const noexcept
2084{
2085 return this->minmax<run_on>(this->domain,comp);
2086}
2087
2088template <class T>
2089template <RunOn run_on>
2090std::pair<T,T>
2091BaseFab<T>::minmax (const Box& subbox, int comp) const noexcept
2092{
2093 Array4<T const> const& a = this->const_array(comp);
2094#ifdef AMREX_USE_GPU
2095 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2097 ReduceData<T,T> reduce_data(reduce_op);
2098 using ReduceTuple = typename decltype(reduce_data)::Type;
2099 reduce_op.eval(subbox, reduce_data,
2100 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2101 {
2102 auto const x = a(i,j,k);
2103 return { x, x };
2104 });
2105 ReduceTuple hv = reduce_data.value(reduce_op);
2106 return std::make_pair(amrex::get<0>(hv), amrex::get<1>(hv));
2107 } else
2108#endif
2109 {
2110 T rmax = std::numeric_limits<T>::lowest();
2111 T rmin = std::numeric_limits<T>::max();
2112 amrex::LoopOnCpu(subbox, [=,&rmin,&rmax] (int i, int j, int k) noexcept
2113 {
2114 auto const x = a(i,j,k);
2115 rmin = amrex::min(rmin, x);
2116 rmax = amrex::max(rmax, x);
2117 });
2118 return std::make_pair(rmin,rmax);
2119 }
2120}
2121
2122template <class T>
2123template <RunOn run_on>
2124T
2125BaseFab<T>::maxabs (int comp) const noexcept
2126{
2127 return this->maxabs<run_on>(this->domain,comp);
2128}
2129
2130template <class T>
2131template <RunOn run_on>
2132T
2133BaseFab<T>::maxabs (const Box& subbox, int comp) const noexcept
2134{
2135 Array4<T const> const& a = this->const_array(comp);
2136#ifdef AMREX_USE_GPU
2137 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2138 ReduceOps<ReduceOpMax> reduce_op;
2139 ReduceData<T> reduce_data(reduce_op);
2140 using ReduceTuple = typename decltype(reduce_data)::Type;
2141 reduce_op.eval(subbox, reduce_data,
2142 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2143 {
2144 return { std::abs(a(i,j,k)) };
2145 });
2146 ReduceTuple hv = reduce_data.value(reduce_op);
2147 return amrex::get<0>(hv);
2148 } else
2149#endif
2150 {
2151 T r = 0;
2152 amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
2153 {
2154 r = amrex::max(r, std::abs(a(i,j,k)));
2155 });
2156 return r;
2157 }
2158}
2159
2160
2161template <class T>
2162template <RunOn run_on>
2163IntVect
2164BaseFab<T>::indexFromValue (Box const& subbox, int comp, T const& value) const noexcept
2165{
2166 Array4<T const> const& a = this->const_array(comp);
2167#ifdef AMREX_USE_GPU
2168 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2169 Array<int,1+AMREX_SPACEDIM> ha{0,AMREX_D_DECL(std::numeric_limits<int>::lowest(),
2170 std::numeric_limits<int>::lowest(),
2171 std::numeric_limits<int>::lowest())};
2172 Gpu::DeviceVector<int> dv(1+AMREX_SPACEDIM);
2173 int* p = dv.data();
2174 Gpu::htod_memcpy_async(p, ha.data(), sizeof(int)*ha.size());
2175 amrex::ParallelFor(subbox, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
2176 {
2177 int* flag = p;
2178 if ((*flag == 0) && (a(i,j,k) == value)) {
2179 if (Gpu::Atomic::Exch(flag,1) == 0) {
2180 AMREX_D_TERM(p[1] = i;,
2181 p[2] = j;,
2182 p[3] = k;);
2183 }
2184 }
2185 });
2186 Gpu::dtoh_memcpy_async(ha.data(), p, sizeof(int)*ha.size());
2188 return IntVect(AMREX_D_DECL(ha[1],ha[2],ha[2]));
2189 } else
2190#endif
2191 {
2192 AMREX_LOOP_3D(subbox, i, j, k,
2193 {
2194 if (a(i,j,k) == value) { return IntVect(AMREX_D_DECL(i,j,k)); }
2195 });
2196 return IntVect::TheMinVector();
2197 }
2198}
2199
2200template <class T>
2201template <RunOn run_on>
2202IntVect
2203BaseFab<T>::minIndex (int comp) const noexcept
2204{
2205 return this->minIndex<run_on>(this->domain,comp);
2206}
2207
2208template <class T>
2209template <RunOn run_on>
2210IntVect
2211BaseFab<T>::minIndex (const Box& subbox, int comp) const noexcept
2212{
2213 T min_val = this->min<run_on>(subbox, comp);
2214 return this->indexFromValue<run_on>(subbox, comp, min_val);
2215}
2216
2217template <class T>
2218template <RunOn run_on>
2219void
2220BaseFab<T>::minIndex (const Box& subbox, Real& min_val, IntVect& min_idx, int comp) const noexcept
2221{
2222 min_val = this->min<run_on>(subbox, comp);
2223 min_idx = this->indexFromValue<run_on>(subbox, comp, min_val);
2224}
2225
2226template <class T>
2227template <RunOn run_on>
2228IntVect
2229BaseFab<T>::maxIndex (int comp) const noexcept
2230{
2231 return this->maxIndex<run_on>(this->domain,comp);
2232}
2233
2234template <class T>
2235template <RunOn run_on>
2236IntVect
2237BaseFab<T>::maxIndex (const Box& subbox, int comp) const noexcept
2238{
2239 T max_val = this->max<run_on>(subbox, comp);
2240 return this->indexFromValue<run_on>(subbox, comp, max_val);
2241}
2242
2243template <class T>
2244template <RunOn run_on>
2245void
2246BaseFab<T>::maxIndex (const Box& subbox, Real& max_val, IntVect& max_idx, int comp) const noexcept
2247{
2248 max_val = this->max<run_on>(subbox, comp);
2249 max_idx = this->indexFromValue<run_on>(subbox, comp, max_val);
2250}
2251
2252template <class T>
2253template <RunOn run_on>
2254int
2255BaseFab<T>::maskLT (BaseFab<int>& mask, T const& val, int comp) const noexcept
2256{
2257 mask.resize(this->domain,1);
2258 int cnt = 0;
2259 Array4<int> const& m = mask.array();
2260 Array4<T const> const& a = this->const_array(comp);
2261#ifdef AMREX_USE_GPU
2262 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2263 ReduceOps<ReduceOpSum> reduce_op;
2264 ReduceData<int> reduce_data(reduce_op);
2265 using ReduceTuple = typename decltype(reduce_data)::Type;
2266 reduce_op.eval(this->domain, reduce_data,
2267 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2268 {
2269 int t;
2270 if (a(i,j,k) < val) {
2271 m(i,j,k) = 1;
2272 t = 1;
2273 } else {
2274 m(i,j,k) = 0;
2275 t = 0;
2276 }
2277 return {t};
2278 });
2279 ReduceTuple hv = reduce_data.value(reduce_op);
2280 cnt = amrex::get<0>(hv);
2281 } else
2282#endif
2283 {
2284 AMREX_LOOP_3D(this->domain, i, j, k,
2285 {
2286 if (a(i,j,k) < val) {
2287 m(i,j,k) = 1;
2288 ++cnt;
2289 } else {
2290 m(i,j,k) = 0;
2291 }
2292 });
2293 }
2294
2295 return cnt;
2296}
2297
2298template <class T>
2299template <RunOn run_on>
2300int
2301BaseFab<T>::maskLE (BaseFab<int>& mask, T const& val, int comp) const noexcept
2302{
2303 mask.resize(this->domain,1);
2304 int cnt = 0;
2305 Array4<int> const& m = mask.array();
2306 Array4<T const> const& a = this->const_array(comp);
2307#ifdef AMREX_USE_GPU
2308 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2309 ReduceOps<ReduceOpSum> reduce_op;
2310 ReduceData<int> reduce_data(reduce_op);
2311 using ReduceTuple = typename decltype(reduce_data)::Type;
2312 reduce_op.eval(this->domain, reduce_data,
2313 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2314 {
2315 int t;
2316 if (a(i,j,k) <= val) {
2317 m(i,j,k) = 1;
2318 t = 1;
2319 } else {
2320 m(i,j,k) = 0;
2321 t = 0;
2322 }
2323 return {t};
2324 });
2325 ReduceTuple hv = reduce_data.value(reduce_op);
2326 cnt = amrex::get<0>(hv);
2327 } else
2328#endif
2329 {
2330 AMREX_LOOP_3D(this->domain, i, j, k,
2331 {
2332 if (a(i,j,k) <= val) {
2333 m(i,j,k) = 1;
2334 ++cnt;
2335 } else {
2336 m(i,j,k) = 0;
2337 }
2338 });
2339 }
2340
2341 return cnt;
2342}
2343
2344template <class T>
2345template <RunOn run_on>
2346int
2347BaseFab<T>::maskEQ (BaseFab<int>& mask, T const& val, int comp) const noexcept
2348{
2349 mask.resize(this->domain,1);
2350 int cnt = 0;
2351 Array4<int> const& m = mask.array();
2352 Array4<T const> const& a = this->const_array(comp);
2353#ifdef AMREX_USE_GPU
2354 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2355 ReduceOps<ReduceOpSum> reduce_op;
2356 ReduceData<int> reduce_data(reduce_op);
2357 using ReduceTuple = typename decltype(reduce_data)::Type;
2358 reduce_op.eval(this->domain, reduce_data,
2359 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2360 {
2361 int t;
2362 if (a(i,j,k) == val) {
2363 m(i,j,k) = 1;
2364 t = 1;
2365 } else {
2366 m(i,j,k) = 0;
2367 t = 0;
2368 }
2369 return {t};
2370 });
2371 ReduceTuple hv = reduce_data.value(reduce_op);
2372 cnt = amrex::get<0>(hv);
2373 } else
2374#endif
2375 {
2376 AMREX_LOOP_3D(this->domain, i, j, k,
2377 {
2378 if (a(i,j,k) == val) {
2379 m(i,j,k) = 1;
2380 ++cnt;
2381 } else {
2382 m(i,j,k) = 0;
2383 }
2384 });
2385 }
2386
2387 return cnt;
2388}
2389
2390template <class T>
2391template <RunOn run_on>
2392int
2393BaseFab<T>::maskGT (BaseFab<int>& mask, T const& val, int comp) const noexcept
2394{
2395 mask.resize(this->domain,1);
2396 int cnt = 0;
2397 Array4<int> const& m = mask.array();
2398 Array4<T const> const& a = this->const_array(comp);
2399#ifdef AMREX_USE_GPU
2400 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2401 ReduceOps<ReduceOpSum> reduce_op;
2402 ReduceData<int> reduce_data(reduce_op);
2403 using ReduceTuple = typename decltype(reduce_data)::Type;
2404 reduce_op.eval(this->domain, reduce_data,
2405 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2406 {
2407 int t;
2408 if (a(i,j,k) > val) {
2409 m(i,j,k) = 1;
2410 t = 1;
2411 } else {
2412 m(i,j,k) = 0;
2413 t = 0;
2414 }
2415 return {t};
2416 });
2417 ReduceTuple hv = reduce_data.value(reduce_op);
2418 cnt = amrex::get<0>(hv);
2419 } else
2420#endif
2421 {
2422 AMREX_LOOP_3D(this->domain, i, j, k,
2423 {
2424 if (a(i,j,k) > val) {
2425 m(i,j,k) = 1;
2426 ++cnt;
2427 } else {
2428 m(i,j,k) = 0;
2429 }
2430 });
2431 }
2432
2433 return cnt;
2434}
2435
2436template <class T>
2437template <RunOn run_on>
2438int
2439BaseFab<T>::maskGE (BaseFab<int>& mask, T const& val, int comp) const noexcept
2440{
2441 mask.resize(this->domain,1);
2442 int cnt = 0;
2443 Array4<int> const& m = mask.array();
2444 Array4<T const> const& a = this->const_array(comp);
2445#ifdef AMREX_USE_GPU
2446 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2447 ReduceOps<ReduceOpSum> reduce_op;
2448 ReduceData<int> reduce_data(reduce_op);
2449 using ReduceTuple = typename decltype(reduce_data)::Type;
2450 reduce_op.eval(this->domain, reduce_data,
2451 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2452 {
2453 int t;
2454 if (a(i,j,k) >= val) {
2455 m(i,j,k) = 1;
2456 t = 1;
2457 } else {
2458 m(i,j,k) = 0;
2459 t = 0;
2460 }
2461 return {t};
2462 });
2463 ReduceTuple hv = reduce_data.value(reduce_op);
2464 cnt = amrex::get<0>(hv);
2465 } else
2466#endif
2467 {
2468 AMREX_LOOP_3D(this->domain, i, j, k,
2469 {
2470 if (a(i,j,k) >= val) {
2471 m(i,j,k) = 1;
2472 ++cnt;
2473 } else {
2474 m(i,j,k) = 0;
2475 }
2476 });
2477 }
2478
2479 return cnt;
2480}
2481
2482template <class T>
2483template <RunOn run_on>
2486{
2487 Box ovlp(this->domain);
2488 ovlp &= x.domain;
2489 return ovlp.ok() ? this->atomicAdd<run_on>(x,ovlp,ovlp,0,0,this->nvar) : *this;
2490}
2491
2492template <class T>
2493template <RunOn run_on>
2495BaseFab<T>::saxpy (T a, const BaseFab<T>& x, const Box& srcbox, const Box& destbox,
2496 int srccomp, int destcomp, int numcomp) noexcept
2497{
2498 BL_ASSERT(srcbox.ok());
2499 BL_ASSERT(x.box().contains(srcbox));
2500 BL_ASSERT(destbox.ok());
2501 BL_ASSERT(box().contains(destbox));
2502 BL_ASSERT(destbox.sameSize(srcbox));
2503 BL_ASSERT( srccomp >= 0 && srccomp+numcomp <= x.nComp());
2504 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2505
2506 Array4<T> const& d = this->array();
2507 Array4<T const> const& s = x.const_array();
2508 const auto dlo = amrex::lbound(destbox);
2509 const auto slo = amrex::lbound(srcbox);
2510 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
2511 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2512 {
2513 d(i,j,k,n+destcomp) += a * s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
2514 });
2515
2516 return *this;
2517}
2518
2519template <class T>
2520template <RunOn run_on>
2522BaseFab<T>::saxpy (T a, const BaseFab<T>& x) noexcept
2523{
2524 Box ovlp(this->domain);
2525 ovlp &= x.domain;
2526 return ovlp.ok() ? saxpy<run_on>(a,x,ovlp,ovlp,0,0,this->nvar) : *this;
2527}
2528
2529template <class T>
2530template <RunOn run_on>
2533 const Box& srcbox, const Box& destbox,
2534 int srccomp, int destcomp, int numcomp) noexcept
2535{
2536 BL_ASSERT(srcbox.ok());
2537 BL_ASSERT(x.box().contains(srcbox));
2538 BL_ASSERT(destbox.ok());
2539 BL_ASSERT(box().contains(destbox));
2540 BL_ASSERT(destbox.sameSize(srcbox));
2541 BL_ASSERT( srccomp >= 0 && srccomp+numcomp <= x.nComp());
2542 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2543
2544 Array4<T> const& d = this->array();
2545 Array4<T const> const& s = x.const_array();
2546 const auto dlo = amrex::lbound(destbox);
2547 const auto slo = amrex::lbound(srcbox);
2548 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
2549 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2550 {
2551 d(i,j,k,n+destcomp) = s(i+offset.x,j+offset.y,k+offset.z,n+srccomp) + a*d(i,j,k,n+destcomp);
2552 });
2553
2554 return *this;
2555}
2556
2557template <class T>
2558template <RunOn run_on>
2560BaseFab<T>::addproduct (const Box& destbox, int destcomp, int numcomp,
2561 const BaseFab<T>& src1, int comp1,
2562 const BaseFab<T>& src2, int comp2) noexcept
2563{
2564 BL_ASSERT(destbox.ok());
2565 BL_ASSERT(box().contains(destbox));
2566 BL_ASSERT( comp1 >= 0 && comp1+numcomp <= src1.nComp());
2567 BL_ASSERT( comp2 >= 0 && comp2+numcomp <= src2.nComp());
2568 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2569
2570 Array4<T> const& d = this->array();
2571 Array4<T const> const& s1 = src1.const_array();
2572 Array4<T const> const& s2 = src2.const_array();
2573 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2574 {
2575 d(i,j,k,n+destcomp) += s1(i,j,k,n+comp1) * s2(i,j,k,n+comp2);
2576 });
2577
2578 return *this;
2579}
2580
2581template <class T>
2582template <RunOn run_on>
2584BaseFab<T>::linComb (const BaseFab<T>& f1, const Box& b1, int comp1,
2585 const BaseFab<T>& f2, const Box& b2, int comp2,
2586 Real alpha, Real beta, const Box& b,
2587 int comp, int numcomp) noexcept
2588{
2589 BL_ASSERT(b1.ok());
2590 BL_ASSERT(f1.box().contains(b1));
2591 BL_ASSERT(b2.ok());
2592 BL_ASSERT(f2.box().contains(b2));
2593 BL_ASSERT(b.ok());
2594 BL_ASSERT(box().contains(b));
2595 BL_ASSERT(b.sameSize(b1));
2596 BL_ASSERT(b.sameSize(b2));
2597 BL_ASSERT(comp1 >= 0 && comp1+numcomp <= f1.nComp());
2598 BL_ASSERT(comp2 >= 0 && comp2+numcomp <= f2.nComp());
2599 BL_ASSERT(comp >= 0 && comp +numcomp <= nComp());
2600
2601 Array4<T> const& d = this->array();
2602 Array4<T const> const& s1 = f1.const_array();
2603 Array4<T const> const& s2 = f2.const_array();
2604 const auto dlo = amrex::lbound(b);
2605 const auto slo1 = amrex::lbound(b1);
2606 const auto slo2 = amrex::lbound(b2);
2607 const Dim3 off1{slo1.x-dlo.x,slo1.y-dlo.y,slo1.z-dlo.z};
2608 const Dim3 off2{slo2.x-dlo.x,slo2.y-dlo.y,slo2.z-dlo.z};
2609
2610 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, b, numcomp, i, j, k, n,
2611 {
2612 d(i,j,k,n+comp) = alpha*s1(i+off1.x,j+off1.y,k+off1.z,n+comp1)
2613 + beta*s2(i+off2.x,j+off2.y,k+off2.z,n+comp2);
2614 });
2615 return *this;
2616}
2617
2618template <class T>
2619template <RunOn run_on>
2620T
2621BaseFab<T>::dot (const Box& xbx, int xcomp,
2622 const BaseFab<T>& y, const Box& ybx, int ycomp,
2623 int numcomp) const noexcept
2624{
2625 BL_ASSERT(xbx.ok());
2626 BL_ASSERT(box().contains(xbx));
2627 BL_ASSERT(y.box().contains(ybx));
2628 BL_ASSERT(xbx.sameSize(ybx));
2629 BL_ASSERT(xcomp >= 0 && xcomp+numcomp <= nComp());
2630 BL_ASSERT(ycomp >= 0 && ycomp+numcomp <= y.nComp());
2631
2632 T r = 0;
2633
2634 const auto xlo = amrex::lbound(xbx);
2635 const auto ylo = amrex::lbound(ybx);
2636 const Dim3 offset{ylo.x-xlo.x,ylo.y-xlo.y,ylo.z-xlo.z};
2637 Array4<T const> const& xa = this->const_array();
2638 Array4<T const> const& ya = y.const_array();
2639
2640#ifdef AMREX_USE_GPU
2641 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2642 ReduceOps<ReduceOpSum> reduce_op;
2643 ReduceData<T> reduce_data(reduce_op);
2644 using ReduceTuple = typename decltype(reduce_data)::Type;
2645 reduce_op.eval(xbx, reduce_data,
2646 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2647 {
2648 T t = 0;
2649 for (int n = 0; n < numcomp; ++n) {
2650 t += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp);
2651 }
2652 return {t};
2653 });
2654 ReduceTuple hv = reduce_data.value(reduce_op);
2655 r = amrex::get<0>(hv);
2656 } else
2657#endif
2658 {
2659 AMREX_LOOP_4D(xbx, numcomp, i, j, k, n,
2660 {
2661 r += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp);
2662 });
2663 }
2664
2665 return r;
2666}
2667
2668template <class T>
2669template <RunOn run_on>
2670T
2671BaseFab<T>::dotmask (const BaseFab<int>& mask, const Box& xbx, int xcomp,
2672 const BaseFab<T>& y, const Box& ybx, int ycomp,
2673 int numcomp) const noexcept
2674{
2675 BL_ASSERT(xbx.ok());
2676 BL_ASSERT(box().contains(xbx));
2677 BL_ASSERT(y.box().contains(ybx));
2678 BL_ASSERT(xbx.sameSize(ybx));
2679 BL_ASSERT(xcomp >= 0 && xcomp+numcomp <= nComp());
2680 BL_ASSERT(ycomp >= 0 && ycomp+numcomp <= y.nComp());
2681
2682 T r = 0;
2683
2684 const auto xlo = amrex::lbound(xbx);
2685 const auto ylo = amrex::lbound(ybx);
2686 const Dim3 offset{ylo.x-xlo.x,ylo.y-xlo.y,ylo.z-xlo.z};
2687
2688 Array4<T const> const& xa = this->const_array();
2689 Array4<T const> const& ya = y.const_array();
2690 Array4<int const> const& ma = mask.const_array();
2691
2692#ifdef AMREX_USE_GPU
2693 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2694 ReduceOps<ReduceOpSum> reduce_op;
2695 ReduceData<T> reduce_data(reduce_op);
2696 using ReduceTuple = typename decltype(reduce_data)::Type;
2697 reduce_op.eval(xbx, reduce_data,
2698 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2699 {
2700 int m = static_cast<int>(static_cast<bool>(ma(i,j,k)));
2701 T t = 0;
2702 for (int n = 0; n < numcomp; ++n) {
2703 t += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp) * m;
2704 }
2705 return {t};
2706 });
2707 ReduceTuple hv = reduce_data.value(reduce_op);
2708 r = amrex::get<0>(hv);
2709 } else
2710#endif
2711 {
2712 AMREX_LOOP_4D(xbx, numcomp, i, j, k, n,
2713 {
2714 int m = static_cast<int>(static_cast<bool>(ma(i,j,k)));
2715 r += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp) * m;
2716 });
2717 }
2718
2719 return r;
2720}
2721
2722template <class T>
2723template <RunOn run_on>
2724T
2725BaseFab<T>::sum (int comp, int numcomp) const noexcept
2726{
2727 return this->sum<run_on>(this->domain, DestComp{comp}, NumComps{numcomp});
2728}
2729
2730template <class T>
2731template <RunOn run_on>
2732T
2733BaseFab<T>::sum (const Box& subbox, int comp, int numcomp) const noexcept
2734{
2735 return this->sum<run_on>(subbox, DestComp{comp}, NumComps{numcomp});
2736}
2737
2738template <class T>
2739template <RunOn run_on>
2741BaseFab<T>::negate (int comp, int numcomp) noexcept
2742{
2743 return this->negate<run_on>(this->domain, DestComp{comp}, NumComps{numcomp});
2744}
2745
2746template <class T>
2747template <RunOn run_on>
2749BaseFab<T>::negate (const Box& b, int comp, int numcomp) noexcept
2750{
2751 return this->negate<run_on>(b, DestComp{comp}, NumComps{numcomp});
2752}
2753
2754template <class T>
2755template <RunOn run_on>
2757BaseFab<T>::invert (T const& r, int comp, int numcomp) noexcept
2758{
2759 return this->invert<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
2760}
2761
2762template <class T>
2763template <RunOn run_on>
2765BaseFab<T>::invert (T const& r, const Box& b, int comp, int numcomp) noexcept
2766{
2767 return this->invert<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
2768}
2769
2770template <class T>
2771template <RunOn run_on>
2773BaseFab<T>::plus (T const& r, int comp, int numcomp) noexcept
2774{
2775 return this->plus<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
2776}
2777
2778template <class T>
2779template <RunOn run_on>
2781BaseFab<T>::plus (T const& r, const Box& b, int comp, int numcomp) noexcept
2782{
2783 return this->plus<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
2784}
2785
2786template <class T>
2787template <RunOn run_on>
2789BaseFab<T>::plus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
2790{
2791 return this->plus<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
2792}
2793
2794template <class T>
2795template <RunOn run_on>
2797BaseFab<T>::atomicAdd (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
2798{
2799 Box ovlp(this->domain);
2800 ovlp &= src.domain;
2801 return ovlp.ok() ? this->atomicAdd<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
2802}
2803
2804template <class T>
2805template <RunOn run_on>
2807BaseFab<T>::plus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
2808 int numcomp) noexcept
2809{
2810 return this->plus<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
2811}
2812
2813template <class T>
2814template <RunOn run_on>
2816BaseFab<T>::atomicAdd (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
2817 int numcomp) noexcept
2818{
2819 Box ovlp(this->domain);
2820 ovlp &= src.domain;
2821 ovlp &= subbox;
2822 return ovlp.ok() ? this->atomicAdd<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
2823}
2824
2825template <class T>
2826template <RunOn run_on>
2828BaseFab<T>::plus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
2829 int srccomp, int destcomp, int numcomp) noexcept
2830{
2831 BL_ASSERT(destbox.ok());
2832 BL_ASSERT(src.box().contains(srcbox));
2833 BL_ASSERT(box().contains(destbox));
2834 BL_ASSERT(destbox.sameSize(srcbox));
2835 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
2836 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2837
2838 Array4<T> const& d = this->array();
2839 Array4<T const> const& s = src.const_array();
2840 const auto dlo = amrex::lbound(destbox);
2841 const auto slo = amrex::lbound(srcbox);
2842 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
2843 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2844 {
2845 d(i,j,k,n+destcomp) += s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
2846 });
2847
2848 return *this;
2849}
2850
2851namespace detail
2852{
2853
2854template <RunOn run_on, typename T,
2855 std::enable_if_t<HasAtomicAdd<T>::value,int> FOO = 0>
2856void basefab_atomic_add (BaseFab<T>& dfab, const BaseFab<T>& sfab,
2857 const Box& srcbox, const Box& destbox,
2858 int srccomp, int destcomp, int numcomp) noexcept
2859{
2860 Array4<T> const& d = dfab.array();
2861 Array4<T const> const& s = sfab.const_array();
2862 const auto dlo = amrex::lbound(destbox);
2863 const auto slo = amrex::lbound(srcbox);
2864 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
2865 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2866 {
2867 T* p = d.ptr(i,j,k,n+destcomp);
2868 HostDevice::Atomic::Add(p, s(i+offset.x,j+offset.y,k+offset.z,n+srccomp));
2869 });
2870}
2871
2872template <RunOn run_on, typename T,
2873 std::enable_if_t<! HasAtomicAdd<T>::value,int> FOO = 0>
2874void basefab_atomic_add (BaseFab<T>& dfab, const BaseFab<T>& sfab,
2875 const Box& srcbox, const Box& destbox,
2876 int srccomp, int destcomp, int numcomp) noexcept
2877{
2878 amrex::ignore_unused(dfab, sfab, srcbox, destbox, srccomp, destcomp, numcomp);
2879 amrex::Abort("BaseFab: atomicAdd not supported");
2880}
2881
2882}
2883
2884template <class T>
2885template <RunOn run_on>
2886BaseFab<T>&
2887BaseFab<T>::atomicAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
2888 int srccomp, int destcomp, int numcomp) noexcept
2889{
2890 BL_ASSERT(destbox.ok());
2891 BL_ASSERT(src.box().contains(srcbox));
2892 BL_ASSERT(box().contains(destbox));
2893 BL_ASSERT(destbox.sameSize(srcbox));
2894 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
2895 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2896
2897 detail::basefab_atomic_add<run_on>(*this, src, srcbox, destbox,
2898 srccomp, destcomp, numcomp);
2899
2900 return *this;
2901}
2902
2903template <class T>
2904template <RunOn run_on>
2906BaseFab<T>::lockAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
2907 int srccomp, int destcomp, int numcomp) noexcept
2908{
2909#if defined(AMREX_USE_OMP) && (AMREX_SPACEDIM > 1)
2910#if defined(AMREX_USE_GPU)
2911 if (run_on == RunOn::Host || Gpu::notInLaunchRegion()) {
2912#endif
2913 BL_ASSERT(destbox.ok());
2914 BL_ASSERT(src.box().contains(srcbox));
2915 BL_ASSERT(box().contains(destbox));
2916 BL_ASSERT(destbox.sameSize(srcbox));
2917 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
2918 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2919
2920 Array4<T> const& d = this->array();
2921 Array4<T const> const& s = src.const_array();
2922 auto const& dlo = amrex::lbound(destbox);
2923 auto const& dhi = amrex::ubound(destbox);
2924 auto const& len = amrex::length(destbox);
2925 auto const& slo = amrex::lbound(srcbox);
2926 Dim3 const offset{slo.x-dlo.x, slo.y-dlo.y, slo.z-dlo.z};
2927
2928 int planedim;
2929 int nplanes;
2930 int plo;
2931 if (len.z == 1) {
2932 planedim = 1;
2933 nplanes = len.y;
2934 plo = dlo.y;
2935 } else {
2936 planedim = 2;
2937 nplanes = len.z;
2938 plo = dlo.z;
2939 }
2940
2941 auto* mask = (bool*) amrex_mempool_alloc(sizeof(bool)*nplanes);
2942 for (int ip = 0; ip < nplanes; ++ip) {
2943 mask[ip] = false;
2944 }
2945
2946 int mm = 0;
2947 int planes_left = nplanes;
2948 while (planes_left > 0) {
2949 AMREX_ASSERT(mm < nplanes);
2950 auto const m = mm + plo;
2951 auto* lock = OpenMP::get_lock(m);
2952 if (omp_test_lock(lock))
2953 {
2954 auto lo = dlo;
2955 auto hi = dhi;
2956 if (planedim == 1) {
2957 lo.y = m;
2958 hi.y = m;
2959 } else {
2960 lo.z = m;
2961 hi.z = m;
2962 }
2963
2964 for (int n = 0; n < numcomp; ++n) {
2965 for (int k = lo.z; k <= hi.z; ++k) {
2966 for (int j = lo.y; j <= hi.y; ++j) {
2967 auto * pdst = d.ptr(dlo.x,j ,k ,n+destcomp);
2968 auto const* psrc = s.ptr(slo.x,j+offset.y,k+offset.z,n+ srccomp);
2969#pragma omp simd
2970 for (int ii = 0; ii < len.x; ++ii) {
2971 pdst[ii] += psrc[ii];
2972 }
2973 }
2974 }
2975 }
2976
2977 mask[mm] = true;
2978 --planes_left;
2979 omp_unset_lock(lock);
2980 if (planes_left == 0) { break; }
2981 }
2982
2983 ++mm;
2984 for (int ip = 0; ip < nplanes; ++ip) {
2985 int new_mm = (mm+ip) % nplanes;
2986 if ( ! mask[new_mm] ) {
2987 mm = new_mm;
2988 break;
2989 }
2990 }
2991 }
2992
2994
2995 return *this;
2996
2997#if defined(AMREX_USE_GPU)
2998 } else {
2999 return this->template atomicAdd<run_on>(src, srcbox, destbox, srccomp, destcomp, numcomp);
3000 }
3001#endif
3002#else
3003 return this->template atomicAdd<run_on>(src, srcbox, destbox, srccomp, destcomp, numcomp);
3004#endif
3005}
3006
3007template <class T>
3008template <RunOn run_on>
3010BaseFab<T>::minus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
3011{
3012 return this->minus<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3013}
3014
3015template <class T>
3016template <RunOn run_on>
3018BaseFab<T>::minus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
3019{
3020 return this->minus<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3021}
3022
3023template <class T>
3024template <RunOn run_on>
3026BaseFab<T>::minus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
3027 int srccomp, int destcomp, int numcomp) noexcept
3028{
3029 BL_ASSERT(destbox.ok());
3030 BL_ASSERT(src.box().contains(srcbox));
3031 BL_ASSERT(box().contains(destbox));
3032 BL_ASSERT(destbox.sameSize(srcbox));
3033 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
3034 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
3035
3036 Array4<T> const& d = this->array();
3037 Array4<T const> const& s = src.const_array();
3038 const auto dlo = amrex::lbound(destbox);
3039 const auto slo = amrex::lbound(srcbox);
3040 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
3041 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
3042 {
3043 d(i,j,k,n+destcomp) -= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
3044 });
3045
3046 return *this;
3047}
3048
3049template <class T>
3050template <RunOn run_on>
3052BaseFab<T>::mult (T const& r, int comp, int numcomp) noexcept
3053{
3054 return this->mult<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
3055}
3056
3057template <class T>
3058template <RunOn run_on>
3060BaseFab<T>::mult (T const& r, const Box& b, int comp, int numcomp) noexcept
3061{
3062 return this->mult<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
3063}
3064
3065template <class T>
3066template <RunOn run_on>
3068BaseFab<T>::mult (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
3069{
3070 return this->mult<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3071}
3072
3073template <class T>
3074template <RunOn run_on>
3076BaseFab<T>::mult (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
3077{
3078 return this->mult<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3079}
3080
3081template <class T>
3082template <RunOn run_on>
3084BaseFab<T>::mult (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
3085 int srccomp, int destcomp, int numcomp) noexcept
3086{
3087 BL_ASSERT(destbox.ok());
3088 BL_ASSERT(src.box().contains(srcbox));
3089 BL_ASSERT(box().contains(destbox));
3090 BL_ASSERT(destbox.sameSize(srcbox));
3091 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
3092 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
3093
3094 Array4<T> const& d = this->array();
3095 Array4<T const> const& s = src.const_array();
3096 const auto dlo = amrex::lbound(destbox);
3097 const auto slo = amrex::lbound(srcbox);
3098 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
3099 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
3100 {
3101 d(i,j,k,n+destcomp) *= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
3102 });
3103
3104 return *this;
3105}
3106
3107template <class T>
3108template <RunOn run_on>
3110BaseFab<T>::divide (T const& r, int comp, int numcomp) noexcept
3111{
3112 return this->divide<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
3113}
3114
3115template <class T>
3116template <RunOn run_on>
3118BaseFab<T>::divide (T const& r, const Box& b, int comp, int numcomp) noexcept
3119{
3120 return this->divide<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
3121}
3122
3123template <class T>
3124template <RunOn run_on>
3126BaseFab<T>::divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
3127{
3128 return this->divide<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3129}
3130
3131template <class T>
3132template <RunOn run_on>
3134BaseFab<T>::divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
3135{
3136 return this->divide<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3137}
3138
3139template <class T>
3140template <RunOn run_on>
3142BaseFab<T>::divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
3143 int srccomp, int destcomp, int numcomp) noexcept
3144{
3145 BL_ASSERT(destbox.ok());
3146 BL_ASSERT(src.box().contains(srcbox));
3147 BL_ASSERT(box().contains(destbox));
3148 BL_ASSERT(destbox.sameSize(srcbox));
3149 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
3150 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
3151
3152 Array4<T> const& d = this->array();
3153 Array4<T const> const& s = src.const_array();
3154 const auto dlo = amrex::lbound(destbox);
3155 const auto slo = amrex::lbound(srcbox);
3156 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
3157 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
3158 {
3159 d(i,j,k,n+destcomp) /= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
3160 });
3161
3162 return *this;
3163}
3164
3165template <class T>
3166template <RunOn run_on>
3169{
3170 Box ovlp(this->domain);
3171 ovlp &= src.domain;
3172 return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,0,0,this->nvar) : *this;
3173}
3174
3175template <class T>
3176template <RunOn run_on>
3178BaseFab<T>::protected_divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
3179{
3180 Box ovlp(this->domain);
3181 ovlp &= src.domain;
3182 return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
3183}
3184
3185template <class T>
3186template <RunOn run_on>
3188BaseFab<T>::protected_divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
3189 int numcomp) noexcept
3190{
3191 Box ovlp(this->domain);
3192 ovlp &= src.domain;
3193 ovlp &= subbox;
3194 return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
3195}
3196
3197template <class T>
3198template <RunOn run_on>
3200BaseFab<T>::protected_divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
3201 int srccomp, int destcomp, int numcomp) noexcept
3202{
3203 BL_ASSERT(destbox.ok());
3204 BL_ASSERT(src.box().contains(srcbox));
3205 BL_ASSERT(box().contains(destbox));
3206 BL_ASSERT(destbox.sameSize(srcbox));
3207 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
3208 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
3209
3210 Array4<T> const& d = this->array();
3211 Array4<T const> const& s = src.const_array();
3212 const auto dlo = amrex::lbound(destbox);
3213 const auto slo = amrex::lbound(srcbox);
3214 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
3215 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
3216 {
3217 if (s(i+offset.x,j+offset.y,k+offset.z,n+srccomp) != 0) {
3218 d(i,j,k,n+destcomp) /= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
3219 }
3220 });
3221
3222 return *this;
3223}
3224
3235template <class T>
3236template <RunOn run_on>
3238BaseFab<T>::linInterp (const BaseFab<T>& f1, const Box& b1, int comp1,
3239 const BaseFab<T>& f2, const Box& b2, int comp2,
3240 Real t1, Real t2, Real t,
3241 const Box& b, int comp, int numcomp) noexcept
3242{
3243 if (amrex::almostEqual(t1,t2) || amrex::almostEqual(t1,t)) {
3244 return copy<run_on>(f1,b1,comp1,b,comp,numcomp);
3245 } else if (amrex::almostEqual(t2,t)) {
3246 return copy<run_on>(f2,b2,comp2,b,comp,numcomp);
3247 } else {
3248 Real alpha = (t2-t)/(t2-t1);
3249 Real beta = (t-t1)/(t2-t1);
3250 return linComb<run_on>(f1,b1,comp1,f2,b2,comp2,alpha,beta,b,comp,numcomp);
3251 }
3252}
3253
3254template <class T>
3255template <RunOn run_on>
3257BaseFab<T>::linInterp (const BaseFab<T>& f1, int comp1,
3258 const BaseFab<T>& f2, int comp2,
3259 Real t1, Real t2, Real t,
3260 const Box& b, int comp, int numcomp) noexcept
3261{
3262 if (amrex::almostEqual(t1,t2) || amrex::almostEqual(t1,t)) {
3263 return copy<run_on>(f1,b,comp1,b,comp,numcomp);
3264 } else if (amrex::almostEqual(t2,t)) {
3265 return copy<run_on>(f2,b,comp2,b,comp,numcomp);
3266 } else {
3267 Real alpha = (t2-t)/(t2-t1);
3268 Real beta = (t-t1)/(t2-t1);
3269 return linComb<run_on>(f1,b,comp1,f2,b,comp2,alpha,beta,b,comp,numcomp);
3270 }
3271}
3272
3273//
3274// New interfaces
3275//
3276
3277template <class T>
3278template <RunOn run_on>
3279void
3280BaseFab<T>::setVal (T const& val) noexcept
3281{
3282 this->setVal<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3283}
3284
3285template <class T>
3286template <RunOn run_on>
3287void
3288BaseFab<T>::setVal (T const& x, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3289{
3290 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3291 Array4<T> const& a = this->array();
3292 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
3293 {
3294 a(i,j,k,n+dcomp.i) = x;
3295 });
3296}
3297
3298template <class T>
3299template <RunOn run_on>
3300void
3301BaseFab<T>::setValIf (T const& val, const BaseFab<int>& mask) noexcept
3302{
3303 this->setValIf<run_on>(val, this->domain, mask, DestComp{0}, NumComps{this->nvar});
3304}
3305
3306template <class T>
3307template <RunOn run_on>
3308void
3309BaseFab<T>::setValIf (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept
3310{
3311 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3312 Array4<T> const& a = this->array();
3313 Array4<int const> const& m = mask.const_array();
3314 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
3315 {
3316 if (m(i,j,k)) { a(i,j,k,n+dcomp.i) = val; }
3317 });
3318}
3319
3320template <class T>
3321template <RunOn run_on>
3322void
3323BaseFab<T>::setValIfNot (T const& val, const BaseFab<int>& mask) noexcept
3324{
3325 this->setValIfNot<run_on>(val, this->domain, mask, DestComp{0}, NumComps{this->nvar});
3326}
3327
3328template <class T>
3329template <RunOn run_on>
3330void
3331BaseFab<T>::setValIfNot (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept
3332{
3333 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3334 Array4<T> const& a = this->array();
3335 Array4<int const> const& m = mask.const_array();
3336 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
3337 {
3338 if (!m(i,j,k)) { a(i,j,k,n+dcomp.i) = val; }
3339 });
3340}
3341
3342template <class T>
3343template <RunOn run_on>
3344void
3345BaseFab<T>::setComplement (T const& x, const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
3346{
3347#ifdef AMREX_USE_GPU
3348 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3349 Array4<T> const& a = this->array();
3350 amrex::ParallelFor(this->domain, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
3351 {
3352 if (! bx.contains(IntVect(AMREX_D_DECL(i,j,k)))) {
3353 for (int n = dcomp.i; n < ncomp.n+dcomp.i; ++n) {
3354 a(i,j,k,n) = x;
3355 }
3356 }
3357 });
3358 } else
3359#endif
3360 {
3361 const BoxList b_lst = amrex::boxDiff(this->domain,bx);
3362 for (auto const& b : b_lst) {
3363 this->setVal<RunOn::Host>(x, b, dcomp, ncomp);
3364 }
3365 }
3366}
3367
3368template <class T>
3369template <RunOn run_on>
3371BaseFab<T>::copy (const BaseFab<T>& src) noexcept
3372{
3373 this->copy<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3374 return *this;
3375}
3376
3377template <class T>
3378template <RunOn run_on>
3381 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3382{
3383 AMREX_ASSERT(this->domain.sameType(src.domain));
3384 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3385 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3386
3387 bx &= src.domain;
3388
3389 Array4<T> const& d = this->array();
3390 Array4<T const> const& s = src.const_array();
3391 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3392 {
3393 d(i,j,k,n+dcomp.i) = s(i,j,k,n+scomp.i);
3394 });
3395
3396 return *this;
3397}
3398
3399template <class T>
3400template <RunOn run_on>
3402BaseFab<T>::plus (T const& val) noexcept
3403{
3404 return this->plus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3405}
3406
3407template <class T>
3408template <RunOn run_on>
3410BaseFab<T>::operator+= (T const& val) noexcept
3411{
3412 return this->plus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3413}
3414
3415template <class T>
3416template <RunOn run_on>
3418BaseFab<T>::plus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3419{
3420 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3421
3422 Array4<T> const& a = this->array();
3423 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3424 {
3425 a(i,j,k,n+dcomp.i) += val;
3426 });
3427
3428 return *this;
3429}
3430
3431template <class T>
3432template <RunOn run_on>
3434BaseFab<T>::plus (const BaseFab<T>& src) noexcept
3435{
3436 return this->plus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3437}
3438
3439template <class T>
3440template <RunOn run_on>
3443{
3444 return this->plus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3445}
3446
3447template <class T>
3448template <RunOn run_on>
3451 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3452{
3453 AMREX_ASSERT(this->domain.sameType(src.domain));
3454 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3455 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3456
3457 bx &= src.domain;
3458
3459 Array4<T> const& d = this->array();
3460 Array4<T const> const& s = src.const_array();
3461 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3462 {
3463 d(i,j,k,n+dcomp.i) += s(i,j,k,n+scomp.i);
3464 });
3465
3466 return *this;
3467}
3468
3469template <class T>
3470template <RunOn run_on>
3472BaseFab<T>::minus (T const& val) noexcept
3473{
3474 return this->minus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3475}
3476
3477template <class T>
3478template <RunOn run_on>
3480BaseFab<T>::operator-= (T const& val) noexcept
3481{
3482 return this->minus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3483}
3484
3485template <class T>
3486template <RunOn run_on>
3488BaseFab<T>::minus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3489{
3490 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3491
3492 Array4<T> const& a = this->array();
3493 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3494 {
3495 a(i,j,k,n+dcomp.i) -= val;
3496 });
3497
3498 return *this;
3499}
3500
3501template <class T>
3502template <RunOn run_on>
3504BaseFab<T>::minus (const BaseFab<T>& src) noexcept
3505{
3506 return this->minus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3507}
3508
3509template <class T>
3510template <RunOn run_on>
3513{
3514 return this->minus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3515}
3516
3517template <class T>
3518template <RunOn run_on>
3521 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3522{
3523 AMREX_ASSERT(this->domain.sameType(src.domain));
3524 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3525 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3526
3527 bx &= src.domain;
3528
3529 Array4<T> const& d = this->array();
3530 Array4<T const> const& s = src.const_array();
3531 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3532 {
3533 d(i,j,k,n+dcomp.i) -= s(i,j,k,n+scomp.i);
3534 });
3535
3536 return *this;
3537}
3538
3539template <class T>
3540template <RunOn run_on>
3542BaseFab<T>::mult (T const& val) noexcept
3543{
3544 return this->mult<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3545}
3546
3547template <class T>
3548template <RunOn run_on>
3550BaseFab<T>::operator*= (T const& val) noexcept
3551{
3552 return this->mult<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3553}
3554
3555template <class T>
3556template <RunOn run_on>
3558BaseFab<T>::mult (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3559{
3560 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3561
3562 Array4<T> const& a = this->array();
3563 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3564 {
3565 a(i,j,k,n+dcomp.i) *= val;
3566 });
3567
3568 return *this;
3569}
3570
3571template <class T>
3572template <RunOn run_on>
3574BaseFab<T>::mult (const BaseFab<T>& src) noexcept
3575{
3576 return this->mult<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3577}
3578
3579template <class T>
3580template <RunOn run_on>
3583{
3584 return this->mult<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3585}
3586
3587template <class T>
3588template <RunOn run_on>
3591 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3592{
3593 AMREX_ASSERT(this->domain.sameType(src.domain));
3594 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3595 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3596
3597 bx &= src.domain;
3598
3599 Array4<T> const& d = this->array();
3600 Array4<T const> const& s = src.const_array();
3601 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3602 {
3603 d(i,j,k,n+dcomp.i) *= s(i,j,k,n+scomp.i);
3604 });
3605
3606 return *this;
3607}
3608
3609template <class T>
3610template <RunOn run_on>
3612BaseFab<T>::divide (T const& val) noexcept
3613{
3614 return this->divide<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3615}
3616
3617template <class T>
3618template <RunOn run_on>
3620BaseFab<T>::operator/= (T const& val) noexcept
3621{
3622 return this->divide<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3623}
3624
3625template <class T>
3626template <RunOn run_on>
3628BaseFab<T>::divide (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3629{
3630 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3631
3632 Array4<T> const& a = this->array();
3633 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3634 {
3635 a(i,j,k,n+dcomp.i) /= val;
3636 });
3637
3638 return *this;
3639}
3640
3641template <class T>
3642template <RunOn run_on>
3644BaseFab<T>::divide (const BaseFab<T>& src) noexcept
3645{
3646 return this->divide<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3647}
3648
3649template <class T>
3650template <RunOn run_on>
3653{
3654 return this->divide<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3655}
3656
3657template <class T>
3658template <RunOn run_on>
3661 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3662{
3663 AMREX_ASSERT(this->domain.sameType(src.domain));
3664 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3665 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3666
3667 bx &= src.domain;
3668
3669 Array4<T> const& d = this->array();
3670 Array4<T const> const& s = src.const_array();
3671 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3672 {
3673 d(i,j,k,n+dcomp.i) /= s(i,j,k,n+scomp.i);
3674 });
3675
3676 return *this;
3677}
3678
3679template <class T>
3680template <RunOn run_on>
3683{
3684 return this->negate<run_on>(this->domain, DestComp{0}, NumComps{this->nvar});
3685}
3686
3687template <class T>
3688template <RunOn run_on>
3690BaseFab<T>::negate (const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
3691{
3692 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3693
3694 Array4<T> const& a = this->array();
3695 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3696 {
3697 a(i,j,k,n+dcomp.i) = -a(i,j,k,n+dcomp.i);
3698 });
3699
3700 return *this;
3701}
3702
3703template <class T>
3704template <RunOn run_on>
3706BaseFab<T>::invert (T const& r) noexcept
3707{
3708 return this->invert<run_on>(r, this->domain, DestComp{0}, NumComps{this->nvar});
3709}
3710
3711template <class T>
3712template <RunOn run_on>
3714BaseFab<T>::invert (T const& r, const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
3715{
3716 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3717
3718 Array4<T> const& a = this->array();
3719 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3720 {
3721 a(i,j,k,n+dcomp.i) = r / a(i,j,k,n+dcomp.i);
3722 });
3723
3724 return *this;
3725}
3726
3727template <class T>
3728template <RunOn run_on>
3729T
3730BaseFab<T>::sum (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept
3731{
3732 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3733
3734 T r = 0;
3735 Array4<T const> const& a = this->const_array();
3736#ifdef AMREX_USE_GPU
3737 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3738 ReduceOps<ReduceOpSum> reduce_op;
3739 ReduceData<T> reduce_data(reduce_op);
3740 using ReduceTuple = typename decltype(reduce_data)::Type;
3741 reduce_op.eval(bx, reduce_data,
3742 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
3743 {
3744 T t = 0;
3745 for (int n = 0; n < ncomp.n; ++n) {
3746 t += a(i,j,k,n+dcomp.i);
3747 }
3748 return { t };
3749 });
3750 ReduceTuple hv = reduce_data.value(reduce_op);
3751 r = amrex::get<0>(hv);
3752 } else
3753#endif
3754 {
3755 amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
3756 {
3757 r += a(i,j,k,n+dcomp.i);
3758 });
3759 }
3760
3761 return r;
3762}
3763
3764template <class T>
3765template <RunOn run_on>
3766T
3767BaseFab<T>::dot (const BaseFab<T>& src, const Box& bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept
3768{
3769 AMREX_ASSERT(this->domain.sameType(src.domain));
3770 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3771 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3772
3773 T r = 0;
3774 Array4<T const> const& d = this->const_array();
3775 Array4<T const> const& s = src.const_array();
3776#ifdef AMREX_USE_GPU
3777 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3778 ReduceOps<ReduceOpSum> reduce_op;
3779 ReduceData<T> reduce_data(reduce_op);
3780 using ReduceTuple = typename decltype(reduce_data)::Type;
3781 reduce_op.eval(bx, reduce_data,
3782 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
3783 {
3784 T t = 0;
3785 for (int n = 0; n < ncomp.n; ++n) {
3786 t += d(i,j,k,n+dcomp.i) * s(i,j,k,n+scomp.i);
3787 }
3788 return { t };
3789 });
3790 ReduceTuple hv = reduce_data.value(reduce_op);
3791 r = amrex::get<0>(hv);
3792 } else
3793#endif
3794 {
3795 amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
3796 {
3797 r += d(i,j,k,n+dcomp.i) * s(i,j,k,n+scomp.i);
3798 });
3799 }
3800
3801 return r;
3802}
3803
3804template <class T>
3805template <RunOn run_on>
3806T
3807BaseFab<T>::dot (const Box& bx, int destcomp, int numcomp) const noexcept
3808{
3809 return dot<run_on>(bx, DestComp{destcomp}, NumComps{numcomp});
3810}
3811
3812
3813template <class T>
3814template <RunOn run_on>
3815T
3816BaseFab<T>::dot (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept
3817{
3818 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3819
3820 T r = 0;
3821 Array4<T const> const& a = this->const_array();
3822#ifdef AMREX_USE_GPU
3823 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3824 ReduceOps<ReduceOpSum> reduce_op;
3825 ReduceData<T> reduce_data(reduce_op);
3826 using ReduceTuple = typename decltype(reduce_data)::Type;
3827 reduce_op.eval(bx, reduce_data,
3828 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
3829 {
3830 T t = 0;
3831 for (int n = 0; n < ncomp.n; ++n) {
3832 t += a(i,j,k,n+dcomp.i)*a(i,j,k,n+dcomp.i);
3833 }
3834 return { t };
3835 });
3836 ReduceTuple hv = reduce_data.value(reduce_op);
3837 r = amrex::get<0>(hv);
3838 } else
3839#endif
3840 {
3841 amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
3842 {
3843 r += a(i,j,k,n+dcomp.i)*a(i,j,k,n+dcomp.i);
3844 });
3845 }
3846
3847 return r;
3848}
3849
3850template <class T>
3851template <RunOn run_on>
3852T
3853BaseFab<T>::dotmask (const BaseFab<T>& src, const Box& bx, const BaseFab<int>& mask,
3854 SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept
3855{
3856 AMREX_ASSERT(this->domain.sameType(src.domain));
3857 AMREX_ASSERT(this->domain.sameType(mask.domain));
3858 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3859 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3860
3861 T r = 0;
3862 Array4<T const> const& d = this->const_array();
3863 Array4<T const> const& s = src.const_array();
3864 Array4<int const> const& m = mask.const_array();
3865#ifdef AMREX_USE_GPU
3866 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3867 ReduceOps<ReduceOpSum> reduce_op;
3868 ReduceData<T> reduce_data(reduce_op);
3869 using ReduceTuple = typename decltype(reduce_data)::Type;
3870 reduce_op.eval(bx, reduce_data,
3871 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
3872 {
3873 T t = 0;
3874 T mi = static_cast<T>(static_cast<int>(static_cast<bool>(m(i,j,k))));
3875 for (int n = 0; n < ncomp.n; ++n) {
3876 t += d(i,j,k,n+dcomp.i)*s(i,j,k,n+scomp.i)*mi;
3877 }
3878 return { t };
3879 });
3880 ReduceTuple hv = reduce_data.value(reduce_op);
3881 r = amrex::get<0>(hv);
3882 } else
3883#endif
3884 {
3885 amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
3886 {
3887 int mi = static_cast<int>(static_cast<bool>(m(i,j,k)));
3888 r += d(i,j,k,n+dcomp.i)*s(i,j,k,n+scomp.i)*mi;
3889 });
3890 }
3891
3892 return r;
3893}
3894
3895}
3896
3897#endif /*BL_BASEFAB_H*/
#define BL_ASSERT(EX)
Definition AMReX_BLassert.H:39
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_DEFAULT_RUNON
Definition AMReX_GpuControl.H:73
#define AMREX_CUDA_SAFE_CALL(call)
Definition AMReX_GpuError.H:73
#define AMREX_HOST_DEVICE_FOR_1D(...)
Definition AMReX_GpuLaunchMacrosC.nolint.H:105
#define AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(where_to_run, box, nc, i, j, k, n, block)
Definition AMReX_GpuLaunch.nolint.H:73
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
int idir
Definition AMReX_HypreMLABecLap.cpp:1093
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
Real * pdst
Definition AMReX_HypreMLABecLap.cpp:1090
Array4< int const > mask
Definition AMReX_InterpFaceRegister.cpp:93
#define AMREX_LOOP_3D(bx, i, j, k, block)
Definition AMReX_Loop.nolint.H:4
#define AMREX_LOOP_4D(bx, ncomp, i, j, k, n, block)
Definition AMReX_Loop.nolint.H:16
void amrex_mempool_free(void *p)
Definition AMReX_MemPool.cpp:80
void * amrex_mempool_alloc(size_t nbytes)
Definition AMReX_MemPool.cpp:74
#define AMREX_D_TERM(a, b, c)
Definition AMReX_SPACE.H:172
#define AMREX_D_DECL(a, b, c)
Definition AMReX_SPACE.H:171
A virtual base class for objects that manage their own dynamic memory allocation.
Definition AMReX_Arena.H:124
A FortranArrayBox(FAB)-like object.
Definition AMReX_BaseFab.H:190
int maskGE(BaseFab< int > &mask, T const &val, int comp=0) const noexcept
Same as above except mark cells with value greater than or equal to val.
Definition AMReX_BaseFab.H:2439
BaseFab< T > & saxpy(T a, const BaseFab< T > &x, const Box &srcbox, const Box &destbox, int srccomp, int destcomp, int numcomp=1) noexcept
FAB SAXPY (y[i] <- y[i] + a * x[i]), in place.
Definition AMReX_BaseFab.H:2495
Array4< T const > const_array() const noexcept
Definition AMReX_BaseFab.H:418
BaseFab< T > & copy(const BaseFab< T > &src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx does not intersect with src fab.
Definition AMReX_BaseFab.H:3380
T sum(int comp, int numcomp=1) const noexcept
Returns sum of given component of FAB state vector.
Definition AMReX_BaseFab.H:2725
gpuStream_t alloc_stream
Definition AMReX_BaseFab.H:1171
Real norminfmask(const Box &subbox, const BaseFab< int > &mask, int scomp=0, int ncomp=1) const noexcept
Definition AMReX_BaseFab.H:1868
BaseFab< T > & divide(T const &val) noexcept
Scalar division on the whole domain and all components.
Definition AMReX_BaseFab.H:3612
const int * hiVect() const noexcept
Returns the upper corner of the domain.
Definition AMReX_BaseFab.H:329
BaseFab< T > & minus(const BaseFab< T > &src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx does not intersect with src fab.
Definition AMReX_BaseFab.H:3520
BaseFab< T > & lockAdd(const BaseFab< T > &src, const Box &srcbox, const Box &destbox, int srccomp, int destcomp, int numcomp) noexcept
Atomically add srcbox region of src FAB to destbox region of this FAB. The srcbox and destbox must be...
Definition AMReX_BaseFab.H:2906
static void Initialize()
std::size_t copyToMem(const Box &srcbox, int srccomp, int numcomp, void *dst) const noexcept
Copy from the srcbox of this Fab to raw memory and return the number of bytes copied.
Definition AMReX_BaseFab.H:1749
std::size_t addFromMem(const Box &dstbox, int dstcomp, int numcomp, const void *src) noexcept
Add from raw memory to the dstbox of this Fab and return the number of bytes copied.
Definition AMReX_BaseFab.H:1804
std::size_t nBytesOwned() const noexcept
Definition AMReX_BaseFab.H:271
BaseFab< T > & copy(const BaseFab< T > &src) noexcept
Definition AMReX_BaseFab.H:3371
BaseFab< T > & addproduct(const Box &destbox, int destcomp, int numcomp, const BaseFab< T > &src1, int comp1, const BaseFab< T > &src2, int comp2) noexcept
y[i] <- y[i] + x1[i] * x2[i])
Definition AMReX_BaseFab.H:2560
BaseFab< T > & minus(T const &val) noexcept
Scalar subtraction on the whole domain and all components.
Definition AMReX_BaseFab.H:3472
int maskLT(BaseFab< int > &mask, T const &val, int comp=0) const noexcept
Compute mask array with value of 1 in cells where BaseFab has value less than val,...
Definition AMReX_BaseFab.H:2255
BaseFab< T > & plus(T const &val) noexcept
Scalar addition on the whole domain and all components.
Definition AMReX_BaseFab.H:3402
BaseFab< T > & mult(T const &val, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3558
BaseFab< T > & mult(const BaseFab< T > &src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx does not intersect with src fab.
Definition AMReX_BaseFab.H:3590
std::size_t nBytes(const Box &bx, int ncomps) const noexcept
Returns bytes used in the Box for those components.
Definition AMReX_BaseFab.H:276
void setPtr(T *p, Long sz) noexcept
Definition AMReX_BaseFab.H:376
BaseFab< T > & linComb(const BaseFab< T > &f1, const Box &b1, int comp1, const BaseFab< T > &f2, const Box &b2, int comp2, Real alpha, Real beta, const Box &b, int comp, int numcomp=1) noexcept
Linear combination. Result is alpha*f1 + beta*f2. Data is taken from b1 region of f1,...
Definition AMReX_BaseFab.H:2584
void define()
Allocates memory for the BaseFab<T>.
Definition AMReX_BaseFab.H:1458
BaseFab< T > & operator*=(T const &val) noexcept
Definition AMReX_BaseFab.H:3550
void resize(const Box &b, int N=1, Arena *ar=nullptr)
This function resizes a BaseFab so it covers the Box b with N components.
Definition AMReX_BaseFab.H:1629
const IntVect & smallEnd() const noexcept
Returns the lower corner of the domain See class Box for analogue.
Definition AMReX_BaseFab.H:306
BaseFab< T > & mult(T const &r, int comp, int numcomp=1) noexcept
Scalar multiplication, except control which components are multiplied.
Definition AMReX_BaseFab.H:3052
BaseFab< T > & atomicAdd(const BaseFab< T > &x) noexcept
Atomic FAB addition (a[i] <- a[i] + b[i]).
Definition AMReX_BaseFab.H:2485
int maskEQ(BaseFab< int > &mask, T const &val, int comp=0) const noexcept
Same as above except mark cells with value equal to val.
Definition AMReX_BaseFab.H:2347
const int * loVect() const noexcept
Returns the lower corner of the domain.
Definition AMReX_BaseFab.H:319
bool contains(const BaseFab< T > &fab) const noexcept
Returns true if the domain of fab is totally contained within the domain of this BaseFab.
Definition AMReX_BaseFab.H:335
bool isAllocated() const noexcept
Returns true if the data for the FAB has been allocated.
Definition AMReX_BaseFab.H:436
std::unique_ptr< T, DataDeleter > release() noexcept
Release ownership of memory.
Definition AMReX_BaseFab.H:1731
void setVal(T const &x, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3288
BaseFab< T > & operator-=(T const &val) noexcept
Definition AMReX_BaseFab.H:3480
const Box & box() const noexcept
Returns the domain (box) where the array is defined.
Definition AMReX_BaseFab.H:294
void setValIf(T const &val, Box const &bx, const BaseFab< int > &mask, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3309
static void Finalize()
Array4< T > array() noexcept
Definition AMReX_BaseFab.H:400
IntVect indexFromValue(const Box &subbox, int comp, T const &value) const noexcept
Definition AMReX_BaseFab.H:2164
BaseFab< T > & mult(const BaseFab< T > &src) noexcept
Definition AMReX_BaseFab.H:3574
bool shared_memory
Is the memory allocated in shared memory?
Definition AMReX_BaseFab.H:1169
int maskLE(BaseFab< int > &mask, T const &val, int comp=0) const noexcept
Same as above except mark cells with value less than or equal to val.
Definition AMReX_BaseFab.H:2301
void setValIf(T const &val, const BaseFab< int > &mask) noexcept
Definition AMReX_BaseFab.H:3301
BaseFab< T > & plus(const BaseFab< T > &src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx does not intersect with src fab.
Definition AMReX_BaseFab.H:3450
void setValIfNot(T const &val, const BaseFab< int > &mask) noexcept
Definition AMReX_BaseFab.H:3323
BaseFab< T > & xpay(T a, const BaseFab< T > &x, const Box &srcbox, const Box &destbox, int srccomp, int destcomp, int numcomp=1) noexcept
FAB XPAY (y[i] <- x[i] + a * y[i])
Definition AMReX_BaseFab.H:2532
T & operator()(const IntVect &p, int N) noexcept
Returns a reference to the Nth component value defined at position p in the domain....
Definition AMReX_BaseFab.H:1269
std::size_t nBytes() const noexcept
Returns how many bytes used.
Definition AMReX_BaseFab.H:269
std::size_t copyFromMem(const Box &dstbox, int dstcomp, int numcomp, const void *src) noexcept
Copy from raw memory to the dstbox of this Fab and return the number of bytes copied.
Definition AMReX_BaseFab.H:1776
BaseFab< T > & negate() noexcept
on the whole domain and all components
Definition AMReX_BaseFab.H:3682
BaseFab< T > & minus(const BaseFab< T > &src) noexcept
Definition AMReX_BaseFab.H:3504
BaseFab< T > & minus(const BaseFab< T > &src, int srccomp, int destcomp, int numcomp=1) noexcept
Subtract src components (srccomp:srccomp+numcomp-1) to this FABs components (destcomp:destcomp+numcom...
Definition AMReX_BaseFab.H:3010
T value_type
Definition AMReX_BaseFab.H:195
void SetBoxType(const IndexType &typ) noexcept
Change the Box type without change the length.
Definition AMReX_BaseFab.H:981
Array4< T const > array() const noexcept
Definition AMReX_BaseFab.H:382
T maxabs(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2125
BaseFab< T > & operator+=(T const &val) noexcept
Definition AMReX_BaseFab.H:3410
void setComplement(T const &x, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
setVal on the complement of bx in the fab's domain
Definition AMReX_BaseFab.H:3345
BaseFab< T > & minus(T const &val, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3488
Long truesize
nvar*numpts that was allocated on heap.
Definition AMReX_BaseFab.H:1167
void setVal(T const &val) noexcept
Set value on the whole domain and all components.
Definition AMReX_BaseFab.H:3280
const int * nCompPtr() const noexcept
for calls to fortran.
Definition AMReX_BaseFab.H:283
Array4< T const > const_array(int start_comp, int num_comps) const noexcept
Definition AMReX_BaseFab.H:430
Box domain
My index space.
Definition AMReX_BaseFab.H:1165
bool contains(const Box &bx) const noexcept
Returns true if bx is totally contained within the domain of this BaseFab.
Definition AMReX_BaseFab.H:344
T * dptr
The data pointer.
Definition AMReX_BaseFab.H:1164
BaseFab< T > & shift(const IntVect &v) noexcept
Perform shifts upon the domain of the BaseFab. They are completely analogous to the corresponding Box...
Definition AMReX_BaseFab.H:1342
BaseFab< T > & divide(T const &val, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3628
T max(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2045
BaseFab< T > & copy(const BaseFab< T > &src, const Box &srcbox, int srccomp, const Box &destbox, int destcomp, int numcomp) noexcept
The copy functions copy the contents of one BaseFab into another. The destination BaseFab is always t...
Definition AMReX_BaseFab.H:1416
Array4< T > array(int start_comp, int num_comps) noexcept
Definition AMReX_BaseFab.H:412
int nvar
Number components.
Definition AMReX_BaseFab.H:1166
T dot(const Box &xbx, int xcomp, const BaseFab< T > &y, const Box &ybx, int ycomp, int numcomp=1) const noexcept
Dot product of x (i.e.,this) and y.
Definition AMReX_BaseFab.H:2621
BaseFab< T > & divide(T const &r, int comp, int numcomp=1) noexcept
As above except specify which components.
Definition AMReX_BaseFab.H:3110
Real norm(int p, int scomp=0, int numcomp=1) const noexcept
Compute the Lp-norm of this FAB using components (scomp : scomp+ncomp-1). p < 0 -> ERROR p = 0 -> inf...
Definition AMReX_BaseFab.H:1912
Array4< T const > array(int start_comp) const noexcept
Definition AMReX_BaseFab.H:388
BaseFab< T > & operator/=(T const &val) noexcept
Definition AMReX_BaseFab.H:3620
std::pair< T, T > minmax(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2083
Array4< T const > const_array(int start_comp) const noexcept
Definition AMReX_BaseFab.H:424
void fill_snan() noexcept
Definition AMReX_BaseFab.H:1376
void setVal(T const &x, const Box &bx, int dcomp, int ncomp) noexcept
The setVal functions set sub-regions in the BaseFab to a constant value. This most general form speci...
Definition AMReX_BaseFab.H:1400
Long size() const noexcept
Returns the total number of points of all components.
Definition AMReX_BaseFab.H:291
BaseFab< T > & plus(T const &r, const Box &b, int comp=0, int numcomp=1) noexcept
Scalar addition (a[i] <- a[i] + r), most general.
Definition AMReX_BaseFab.H:2781
BaseFab< T > & operator=(const BaseFab< T > &rhs)=delete
void getVal(T *data, const IntVect &pos, int N, int numcomp) const noexcept
This function puts numcomp component values, starting at component N, from position pos in the domain...
Definition AMReX_BaseFab.H:1316
const IntVect & bigEnd() const noexcept
Returns the upper corner of the domain. See class Box for analogue.
Definition AMReX_BaseFab.H:309
Array4< T > array(int start_comp) noexcept
Definition AMReX_BaseFab.H:406
Elixir elixir() noexcept
Definition AMReX_BaseFab.H:1671
Long numPts() const noexcept
Returns the number of points.
Definition AMReX_BaseFab.H:288
const T * dataPtr(int n=0) const noexcept
Same as above except works on const FABs.
Definition AMReX_BaseFab.H:364
void setValIfNot(T const &val, Box const &bx, const BaseFab< int > &mask, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3331
BaseFab< T > & mult(T const &val) noexcept
Scalar multiplication on the whole domain and all components.
Definition AMReX_BaseFab.H:3542
BaseFab< T > & divide(const BaseFab< T > &src, Box bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx does not intersect with src fab.
Definition AMReX_BaseFab.H:3660
void setValIfNot(T const &val, const Box &bx, const BaseFab< int > &mask, int nstart, int num) noexcept
Definition AMReX_BaseFab.H:1408
BaseFab< T > & shiftHalf(int dir, int n_cell) noexcept
Perform shifts upon the domain of the BaseFab. They are completely analogous to the corresponding Box...
Definition AMReX_BaseFab.H:1366
void prefetchToDevice() const noexcept
Definition AMReX_BaseFab.H:1236
T * dataPtr(int n=0) noexcept
Returns a pointer to an object of type T that is the value of the Nth component associated with the c...
Definition AMReX_BaseFab.H:355
bool ptr_owner
Owner of T*?
Definition AMReX_BaseFab.H:1168
virtual ~BaseFab() noexcept
The destructor deletes the array memory.
Definition AMReX_BaseFab.H:1576
IntVect length() const noexcept
Returns a pointer to an array of SPACEDIM integers giving the length of the domain in each direction.
Definition AMReX_BaseFab.H:300
IntVect maxIndex(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2229
BaseFab< T > & protected_divide(const BaseFab< T > &src) noexcept
Divide wherever "src" is "true" or "non-zero".
Definition AMReX_BaseFab.H:3168
friend class BaseFab
Definition AMReX_BaseFab.H:193
BaseFab< T > & invert(T const &r, const Box &b, int comp=0, int numcomp=1) noexcept
Most general version, specify subbox and which components.
Definition AMReX_BaseFab.H:2765
Array4< T const > array(int start_comp, int num_comps) const noexcept
Definition AMReX_BaseFab.H:394
int nComp() const noexcept
Returns the number of components.
Definition AMReX_BaseFab.H:280
int maskGT(BaseFab< int > &mask, T const &val, int comp=0) const noexcept
Same as above except mark cells with value greater than val.
Definition AMReX_BaseFab.H:2393
T dotmask(const BaseFab< int > &mask, const Box &xbx, int xcomp, const BaseFab< T > &y, const Box &ybx, int ycomp, int numcomp) const noexcept
Definition AMReX_BaseFab.H:2671
void clear() noexcept
The function returns the BaseFab to the invalid state. The memory is freed.
Definition AMReX_BaseFab.H:1692
BaseFab< T > & plus(const BaseFab< T > &src) noexcept
Definition AMReX_BaseFab.H:3434
BaseFab() noexcept=default
Construct an empty BaseFab, which must be resized (see BaseFab::resize) before use.
BaseFab< T > & divide(const BaseFab< T > &src) noexcept
Definition AMReX_BaseFab.H:3644
void prefetchToHost() const noexcept
Definition AMReX_BaseFab.H:1203
T min(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2007
void setComplement(T const &x, const Box &b, int ns, int num) noexcept
This function is analogous to the fourth form of setVal above, except that instead of setting values ...
Definition AMReX_BaseFab.H:1832
BaseFab< T > & linInterp(const BaseFab< T > &f1, const Box &b1, int comp1, const BaseFab< T > &f2, const Box &b2, int comp2, Real t1, Real t2, Real t, const Box &b, int comp, int numcomp=1) noexcept
Linear interpolation / extrapolation. Result is (t2-t)/(t2-t1)*f1 + (t-t1)/(t2-t1)*f2 Data is taken f...
Definition AMReX_BaseFab.H:3238
IntVect minIndex(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2203
T * dataPtr(const IntVect &p, int n=0) noexcept
Definition AMReX_BaseFab.H:1178
void abs() noexcept
Compute absolute value for all components of this FAB.
Definition AMReX_BaseFab.H:1840
void getVal(T *data, const IntVect &pos) const noexcept
Same as above, except that starts at component 0 and copies all comps.
Definition AMReX_BaseFab.H:1334
BaseFab< T > & plus(T const &val, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3418
A class for managing a List of Boxes that share a common IndexType. This class implements operations ...
Definition AMReX_BoxList.H:52
__host__ __device__ const IntVectND< dim > & bigEnd() const &noexcept
Return the inclusive upper bound of the box.
Definition AMReX_Box.H:123
__host__ __device__ const int * hiVect() const &noexcept
Return a constant pointer the array of high end coordinates. Useful for calls to FORTRAN.
Definition AMReX_Box.H:191
__host__ __device__ Long numPts() const noexcept
Return the number of points contained in the BoxND.
Definition AMReX_Box.H:356
__host__ __device__ IntVectND< dim > length() const noexcept
Return the length of the BoxND.
Definition AMReX_Box.H:154
__host__ __device__ bool contains(const IntVectND< dim > &p) const noexcept
Return true if argument is contained within BoxND.
Definition AMReX_Box.H:212
__host__ __device__ const int * loVect() const &noexcept
Return a constant pointer the array of low end coordinates. Useful for calls to FORTRAN.
Definition AMReX_Box.H:186
__host__ __device__ BoxND & setType(const IndexTypeND< dim > &t) noexcept
Set indexing type.
Definition AMReX_Box.H:505
__host__ __device__ bool ok() const noexcept
Checks if it is a proper BoxND (including a valid type).
Definition AMReX_Box.H:208
__host__ __device__ const IntVectND< dim > & smallEnd() const &noexcept
Return the inclusive lower bound of the box.
Definition AMReX_Box.H:111
GPU-compatible tuple.
Definition AMReX_Tuple.H:98
static gpuStream_t setStream(gpuStream_t s) noexcept
Definition AMReX_GpuDevice.cpp:735
static gpuStream_t gpuStream() noexcept
Definition AMReX_GpuDevice.H:78
static int deviceId() noexcept
Definition AMReX_GpuDevice.cpp:683
static int devicePropMajor() noexcept
Definition AMReX_GpuDevice.H:166
Definition AMReX_GpuElixir.H:13
__host__ static __device__ constexpr IntVectND< dim > TheMinVector() noexcept
Definition AMReX_IntVect.H:727
Dynamically allocated vector for trivially copyable data.
Definition AMReX_PODVector.H:308
T * data() noexcept
Definition AMReX_PODVector.H:666
Definition AMReX_Reduce.H:453
Type value()
Definition AMReX_Reduce.H:485
Definition AMReX_Reduce.H:582
std::enable_if_t< IsFabArray< MF >::value > eval(MF const &mf, IntVect const &nghost, D &reduce_data, F &&f)
Definition AMReX_Reduce.H:716
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:79
amrex_long Long
Definition AMReX_INT.H:30
__host__ __device__ Dim3 ubound(Array4< T > const &a) noexcept
Return the inclusive upper bounds of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1331
__host__ __device__ Dim3 length(Array4< T > const &a) noexcept
Return the spatial extents of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1345
__host__ __device__ Dim3 lbound(Array4< T > const &a) noexcept
Return the inclusive lower bounds of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1317
std::array< T, N > Array
Definition AMReX_Array.H:26
__host__ __device__ AMREX_FORCE_INLINE T Exch(T *address, T val) noexcept
Definition AMReX_GpuAtomic.H:488
void streamSynchronize() noexcept
Definition AMReX_GpuDevice.H:263
void dtoh_memcpy_async(void *p_h, const void *p_d, const std::size_t sz) noexcept
Definition AMReX_GpuDevice.H:315
bool inLaunchRegion() noexcept
Definition AMReX_GpuControl.H:92
bool notInLaunchRegion() noexcept
Definition AMReX_GpuControl.H:93
void htod_memcpy_async(void *p_d, const void *p_h, const std::size_t sz) noexcept
Definition AMReX_GpuDevice.H:301
gpuStream_t gpuStream() noexcept
Definition AMReX_GpuDevice.H:244
__host__ __device__ AMREX_FORCE_INLINE void Add(T *const sum, T const value) noexcept
Definition AMReX_GpuAtomic.H:622
Definition AMReX_Amr.cpp:49
MakeType
Definition AMReX_MakeType.H:7
@ make_deep_copy
Definition AMReX_MakeType.H:7
@ make_alias
Definition AMReX_MakeType.H:7
__host__ __device__ void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:139
int nComp(FabArrayBase const &fa)
Definition AMReX_FabArrayBase.cpp:2851
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
Definition AMReX_CTOParallelForImpl.H:193
__host__ __device__ Array4< T > makeArray4(T *p, Box const &bx, int ncomp) noexcept
Definition AMReX_BaseFab.H:92
RunOn
Definition AMReX_GpuControl.H:69
std::enable_if_t< std::is_arithmetic_v< T > > placementNew(T *const, Long)
Definition AMReX_BaseFab.H:99
cudaStream_t gpuStream_t
Definition AMReX_GpuControl.H:83
bool InitSNaN() noexcept
Definition AMReX.cpp:173
Long TotalBytesAllocatedInFabs() noexcept
Definition AMReX_BaseFab.cpp:66
BoxND< 3 > Box
Box is an alias for amrex::BoxND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:30
void BaseFab_Initialize()
Definition AMReX_BaseFab.cpp:30
__host__ __device__ constexpr const T & min(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:24
void BaseFab_Finalize()
Definition AMReX_BaseFab.cpp:59
__host__ __device__ Dim3 begin(BoxND< dim > const &box) noexcept
Definition AMReX_Box.H:2006
void ResetTotalBytesAllocatedInFabsHWM() noexcept
Definition AMReX_BaseFab.cpp:134
BoxList boxDiff(const Box &b1in, const Box &b2)
Returns BoxList defining the compliment of b2 in b1in.
Definition AMReX_BoxList.cpp:599
__host__ __device__ std::enable_if_t< std::is_floating_point_v< T >, bool > almostEqual(T x, T y, int ulp=2)
Definition AMReX_Algorithm.H:122
IntVectND< 3 > IntVect
IntVect is an alias for amrex::IntVectND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:33
Long TotalBytesAllocatedInFabsHWM() noexcept
Definition AMReX_BaseFab.cpp:83
__host__ __device__ constexpr const T & max(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:44
Long TotalCellsAllocatedInFabsHWM() noexcept
Definition AMReX_BaseFab.cpp:117
Long TotalCellsAllocatedInFabs() noexcept
Definition AMReX_BaseFab.cpp:100
void Error(const std::string &msg)
Print out message to cerr and exit via amrex::Abort().
Definition AMReX.cpp:224
std::enable_if_t< std::is_trivially_destructible_v< T > > placementDelete(T *const, Long)
Definition AMReX_BaseFab.H:124
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:230
void LoopOnCpu(Dim3 lo, Dim3 hi, F const &f) noexcept
Definition AMReX_Loop.H:365
void update_fab_stats(Long n, Long s, size_t szt) noexcept
Definition AMReX_BaseFab.cpp:146
__host__ __device__ Dim3 end(BoxND< dim > const &box) noexcept
Definition AMReX_Box.H:2015
A multidimensional array accessor.
Definition AMReX_Array4.H:283
__host__ __device__ constexpr std::size_t size() const noexcept
Total number of elements in the ArrayND's index region.
Definition AMReX_Array4.H:638
__host__ __device__ T * ptr(idx... i) const noexcept
Multi-index ptr() for accessing pointer to element.
Definition AMReX_Array4.H:556
Definition AMReX_DataAllocator.H:9
void * alloc(std::size_t sz) const noexcept
Definition AMReX_DataAllocator.H:16
Arena * arena() const noexcept
Definition AMReX_DataAllocator.H:24
Definition AMReX_DataAllocator.H:29
Definition AMReX_BaseFab.H:77
int i
Definition AMReX_BaseFab.H:80
__host__ __device__ DestComp(int ai) noexcept
Definition AMReX_BaseFab.H:79
Definition AMReX_Dim3.H:12
int x
Definition AMReX_Dim3.H:12
Definition AMReX_BaseFab.H:83
__host__ __device__ NumComps(int an) noexcept
Definition AMReX_BaseFab.H:85
int n
Definition AMReX_BaseFab.H:86
Definition AMReX_BaseFab.H:71
__host__ __device__ SrcComp(int ai) noexcept
Definition AMReX_BaseFab.H:73
int i
Definition AMReX_BaseFab.H:74