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 this->arena()->streamOrderedFree(this->dptr, alloc_stream);
1710#else
1711 this->free(this->dptr);
1712#endif
1713
1714 if (this->nvar > 1) {
1715 amrex::update_fab_stats(-this->truesize/this->nvar, -this->truesize, sizeof(T));
1716 } else {
1717 amrex::update_fab_stats(0, -this->truesize, sizeof(T));
1718 }
1719 }
1720
1721 this->dptr = nullptr;
1722 this->truesize = 0;
1723 }
1724}
1725
1726template <class T>
1727std::unique_ptr<T,DataDeleter>
1729{
1730 std::unique_ptr<T,DataDeleter> r(nullptr, DataDeleter{this->arena()});
1731 if (this->dptr && this->ptr_owner) {
1732 r.reset(this->dptr);
1733 this->ptr_owner = false;
1734 if (this->nvar > 1) {
1735 amrex::update_fab_stats(-this->truesize/this->nvar, -this->truesize, sizeof(T));
1736 } else {
1737 amrex::update_fab_stats(0, -this->truesize, sizeof(T));
1738 }
1739 }
1740 return r;
1741}
1742
1743template <class T>
1744template <RunOn run_on>
1745std::size_t
1747 int srccomp,
1748 int numcomp,
1749 void* dst) const noexcept
1750{
1751 BL_ASSERT(box().contains(srcbox));
1752 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= nComp());
1753
1754 if (srcbox.ok())
1755 {
1756 Array4<T> d(static_cast<T*>(dst),amrex::begin(srcbox),amrex::end(srcbox),numcomp);
1757 Array4<T const> const& s = this->const_array();
1758 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, srcbox, numcomp, i, j, k, n,
1759 {
1760 d(i,j,k,n) = s(i,j,k,n+srccomp);
1761 });
1762 return sizeof(T)*d.size();
1763 }
1764 else
1765 {
1766 return 0;
1767 }
1768}
1769
1770template <class T>
1771template <RunOn run_on, typename BUF>
1772std::size_t
1774 int dstcomp,
1775 int numcomp,
1776 const void* src) noexcept
1777{
1778 BL_ASSERT(box().contains(dstbox));
1779 BL_ASSERT(dstcomp >= 0 && dstcomp+numcomp <= nComp());
1780
1781 if (dstbox.ok())
1782 {
1783 Array4<BUF const> s(static_cast<BUF const*>(src), amrex::begin(dstbox),
1784 amrex::end(dstbox), numcomp);
1785 Array4<T> const& d = this->array();
1786 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, dstbox, numcomp, i, j, k, n,
1787 {
1788 d(i,j,k,n+dstcomp) = static_cast<T>(s(i,j,k,n));
1789 });
1790 return sizeof(BUF)*s.size();
1791 }
1792 else
1793 {
1794 return 0;
1795 }
1796}
1797
1798template <class T>
1799template <RunOn run_on, typename BUF>
1800std::size_t
1802 int dstcomp,
1803 int numcomp,
1804 const void* src) noexcept
1805{
1806 BL_ASSERT(box().contains(dstbox));
1807 BL_ASSERT(dstcomp >= 0 && dstcomp+numcomp <= nComp());
1808
1809 if (dstbox.ok())
1810 {
1811 Array4<BUF const> s(static_cast<BUF const*>(src), amrex::begin(dstbox),
1812 amrex::end(dstbox), numcomp);
1813 Array4<T> const& d = this->array();
1814 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, dstbox, numcomp, i, j, k, n,
1815 {
1816 d(i,j,k,n+dstcomp) += static_cast<T>(s(i,j,k,n));
1817 });
1818 return sizeof(BUF)*s.size();
1819 }
1820 else
1821 {
1822 return 0;
1823 }
1824}
1825
1826template <class T>
1827template <RunOn run_on>
1828void
1829BaseFab<T>::setComplement (T const& x, const Box& b, int ns, int num) noexcept
1830{
1831 this->setComplement<run_on>(x, b, DestComp{ns}, NumComps{num});
1832}
1833
1834template <class T>
1835template <RunOn run_on>
1836void
1838{
1839 this->abs<run_on>(this->domain,0,this->nvar);
1840}
1841
1842template <class T>
1843template <RunOn run_on>
1844void
1845BaseFab<T>::abs (int comp, int numcomp) noexcept
1846{
1847 this->abs<run_on>(this->domain,comp,numcomp);
1848}
1849
1850template <class T>
1851template <RunOn run_on>
1852void
1853BaseFab<T>::abs (const Box& subbox, int comp, int numcomp) noexcept
1854{
1855 Array4<T> const& a = this->array();
1856 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, subbox, numcomp, i, j, k, n,
1857 {
1858 a(i,j,k,n+comp) = std::abs(a(i,j,k,n+comp));
1859 });
1860}
1861
1862template <class T>
1863template <RunOn run_on>
1864Real
1866 int scomp, int ncomp) const noexcept
1867{
1868 BL_ASSERT(this->domain.contains(subbox));
1869 BL_ASSERT(scomp >= 0 && scomp + ncomp <= this->nvar);
1870
1871 Array4<T const> const& a = this->const_array();
1872 Array4<int const> const& m = mask.const_array();
1873 Real r = 0.0;
1874#ifdef AMREX_USE_GPU
1875 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
1876 ReduceOps<ReduceOpMax> reduce_op;
1877 ReduceData<Real> reduce_data(reduce_op);
1878 using ReduceTuple = ReduceData<Real>::Type;
1879 reduce_op.eval(subbox, reduce_data,
1880 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
1881 {
1882 Real t = 0.0;
1883 if (m(i,j,k)) {
1884 for (int n = 0; n < ncomp; ++n) {
1885 t = amrex::max(t,static_cast<Real>(std::abs(a(i,j,k,n+scomp))));
1886 }
1887 }
1888 return {t};
1889 });
1890 ReduceTuple hv = reduce_data.value(reduce_op);
1891 r = amrex::get<0>(hv);
1892 } else
1893#endif
1894 {
1895 amrex::LoopOnCpu(subbox, ncomp, [=,&r] (int i, int j, int k, int n) noexcept
1896 {
1897 if (m(i,j,k)) {
1898 Real t = static_cast<Real>(std::abs(a(i,j,k,n+scomp)));
1899 r = amrex::max(r,t);
1900 }
1901 });
1902 }
1903 return r;
1904}
1905
1906template <class T>
1907template <RunOn run_on>
1908Real
1909BaseFab<T>::norm (int p, int comp, int numcomp) const noexcept
1910{
1911 return norm<run_on>(this->domain,p,comp,numcomp);
1912}
1913
1914template <class T>
1915template <RunOn run_on>
1916Real
1917BaseFab<T>::norm (const Box& subbox, int p, int comp, int numcomp) const noexcept
1918{
1919 BL_ASSERT(this->domain.contains(subbox));
1920 BL_ASSERT(comp >= 0 && comp + numcomp <= this->nvar);
1921
1922 Array4<T const> const& a = this->const_array();
1923 Real nrm = 0.;
1924#ifdef AMREX_USE_GPU
1925 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
1926 if (p == 0) {
1927 ReduceOps<ReduceOpMax> reduce_op;
1928 ReduceData<Real> reduce_data(reduce_op);
1929 using ReduceTuple = ReduceData<Real>::Type;
1930 reduce_op.eval(subbox, reduce_data,
1931 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
1932 {
1933 Real t = 0.0;
1934 for (int n = 0; n < numcomp; ++n) {
1935 t = amrex::max(t, static_cast<Real>(std::abs(a(i,j,k,n+comp))));
1936 }
1937 return {t};
1938 });
1939 ReduceTuple hv = reduce_data.value(reduce_op);
1940 nrm = amrex::get<0>(hv);
1941 } else if (p == 1) {
1942 ReduceOps<ReduceOpSum> reduce_op;
1943 ReduceData<Real> reduce_data(reduce_op);
1944 using ReduceTuple = ReduceData<Real>::Type;
1945 reduce_op.eval(subbox, reduce_data,
1946 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
1947 {
1948 Real t = 0.0;
1949 for (int n = 0; n < numcomp; ++n) {
1950 t += static_cast<Real>(std::abs(a(i,j,k,n+comp)));
1951 }
1952 return {t};
1953 });
1954 ReduceTuple hv = reduce_data.value(reduce_op);
1955 nrm = amrex::get<0>(hv);
1956 } else if (p == 2) {
1957 ReduceOps<ReduceOpSum> reduce_op;
1958 ReduceData<Real> reduce_data(reduce_op);
1959 using ReduceTuple = ReduceData<Real>::Type;
1960 reduce_op.eval(subbox, reduce_data,
1961 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
1962 {
1963 Real t = 0.0;
1964 for (int n = 0; n < numcomp; ++n) {
1965 t += static_cast<Real>(a(i,j,k,n+comp)*a(i,j,k,n+comp));
1966 }
1967 return {t};
1968 });
1969 ReduceTuple hv = reduce_data.value(reduce_op);
1970 nrm = amrex::get<0>(hv);
1971 } else {
1972 amrex::Error("BaseFab<T>::norm: wrong p");
1973 }
1974 } else
1975#endif
1976 {
1977 if (p == 0) {
1978 amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
1979 {
1980 Real t = static_cast<Real>(std::abs(a(i,j,k,n+comp)));
1981 nrm = amrex::max(nrm,t);
1982 });
1983 } else if (p == 1) {
1984 amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
1985 {
1986 nrm += std::abs(a(i,j,k,n+comp));
1987 });
1988 } else if (p == 2) {
1989 amrex::LoopOnCpu(subbox, numcomp, [=,&nrm] (int i, int j, int k, int n) noexcept
1990 {
1991 nrm += a(i,j,k,n+comp)*a(i,j,k,n+comp);
1992 });
1993 } else {
1994 amrex::Error("BaseFab<T>::norm: wrong p");
1995 }
1996 }
1997
1998 return nrm;
1999}
2000
2001template <class T>
2002template <RunOn run_on>
2003T
2004BaseFab<T>::min (int comp) const noexcept
2005{
2006 return this->min<run_on>(this->domain,comp);
2007}
2008
2009template <class T>
2010template <RunOn run_on>
2011T
2012BaseFab<T>::min (const Box& subbox, int comp) const noexcept
2013{
2014 Array4<T const> const& a = this->const_array(comp);
2015#ifdef AMREX_USE_GPU
2016 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2017 ReduceOps<ReduceOpMin> reduce_op;
2018 ReduceData<T> reduce_data(reduce_op);
2019 using ReduceTuple = typename decltype(reduce_data)::Type;
2020 reduce_op.eval(subbox, reduce_data,
2021 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2022 {
2023 return { a(i,j,k) };
2024 });
2025 ReduceTuple hv = reduce_data.value(reduce_op);
2026 return amrex::get<0>(hv);
2027 } else
2028#endif
2029 {
2030 T r = std::numeric_limits<T>::max();
2031 amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
2032 {
2033 r = amrex::min(r, a(i,j,k));
2034 });
2035 return r;
2036 }
2037}
2038
2039template <class T>
2040template <RunOn run_on>
2041T
2042BaseFab<T>::max (int comp) const noexcept
2043{
2044 return this->max<run_on>(this->domain,comp);
2045}
2046
2047template <class T>
2048template <RunOn run_on>
2049T
2050BaseFab<T>::max (const Box& subbox, int comp) const noexcept
2051{
2052 Array4<T const> const& a = this->const_array(comp);
2053#ifdef AMREX_USE_GPU
2054 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2055 ReduceOps<ReduceOpMax> reduce_op;
2056 ReduceData<T> reduce_data(reduce_op);
2057 using ReduceTuple = typename decltype(reduce_data)::Type;
2058 reduce_op.eval(subbox, reduce_data,
2059 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2060 {
2061 return { a(i,j,k) };
2062 });
2063 ReduceTuple hv = reduce_data.value(reduce_op);
2064 return amrex::get<0>(hv);
2065 } else
2066#endif
2067 {
2068 T r = std::numeric_limits<T>::lowest();
2069 amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
2070 {
2071 r = amrex::max(r, a(i,j,k));
2072 });
2073 return r;
2074 }
2075}
2076
2077template <class T>
2078template <RunOn run_on>
2079std::pair<T,T>
2080BaseFab<T>::minmax (int comp) const noexcept
2081{
2082 return this->minmax<run_on>(this->domain,comp);
2083}
2084
2085template <class T>
2086template <RunOn run_on>
2087std::pair<T,T>
2088BaseFab<T>::minmax (const Box& subbox, int comp) const noexcept
2089{
2090 Array4<T const> const& a = this->const_array(comp);
2091#ifdef AMREX_USE_GPU
2092 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2094 ReduceData<T,T> reduce_data(reduce_op);
2095 using ReduceTuple = typename decltype(reduce_data)::Type;
2096 reduce_op.eval(subbox, reduce_data,
2097 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2098 {
2099 auto const x = a(i,j,k);
2100 return { x, x };
2101 });
2102 ReduceTuple hv = reduce_data.value(reduce_op);
2103 return std::make_pair(amrex::get<0>(hv), amrex::get<1>(hv));
2104 } else
2105#endif
2106 {
2107 T rmax = std::numeric_limits<T>::lowest();
2108 T rmin = std::numeric_limits<T>::max();
2109 amrex::LoopOnCpu(subbox, [=,&rmin,&rmax] (int i, int j, int k) noexcept
2110 {
2111 auto const x = a(i,j,k);
2112 rmin = amrex::min(rmin, x);
2113 rmax = amrex::max(rmax, x);
2114 });
2115 return std::make_pair(rmin,rmax);
2116 }
2117}
2118
2119template <class T>
2120template <RunOn run_on>
2121T
2122BaseFab<T>::maxabs (int comp) const noexcept
2123{
2124 return this->maxabs<run_on>(this->domain,comp);
2125}
2126
2127template <class T>
2128template <RunOn run_on>
2129T
2130BaseFab<T>::maxabs (const Box& subbox, int comp) const noexcept
2131{
2132 Array4<T const> const& a = this->const_array(comp);
2133#ifdef AMREX_USE_GPU
2134 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2135 ReduceOps<ReduceOpMax> reduce_op;
2136 ReduceData<T> reduce_data(reduce_op);
2137 using ReduceTuple = typename decltype(reduce_data)::Type;
2138 reduce_op.eval(subbox, reduce_data,
2139 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2140 {
2141 return { std::abs(a(i,j,k)) };
2142 });
2143 ReduceTuple hv = reduce_data.value(reduce_op);
2144 return amrex::get<0>(hv);
2145 } else
2146#endif
2147 {
2148 T r = 0;
2149 amrex::LoopOnCpu(subbox, [=,&r] (int i, int j, int k) noexcept
2150 {
2151 r = amrex::max(r, std::abs(a(i,j,k)));
2152 });
2153 return r;
2154 }
2155}
2156
2157
2158template <class T>
2159template <RunOn run_on>
2160IntVect
2161BaseFab<T>::indexFromValue (Box const& subbox, int comp, T const& value) const noexcept
2162{
2163 Array4<T const> const& a = this->const_array(comp);
2164#ifdef AMREX_USE_GPU
2165 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2166 Array<int,1+AMREX_SPACEDIM> ha{0,AMREX_D_DECL(std::numeric_limits<int>::lowest(),
2167 std::numeric_limits<int>::lowest(),
2168 std::numeric_limits<int>::lowest())};
2169 Gpu::DeviceVector<int> dv(1+AMREX_SPACEDIM);
2170 int* p = dv.data();
2171 Gpu::htod_memcpy_async(p, ha.data(), sizeof(int)*ha.size());
2172 amrex::ParallelFor(subbox, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
2173 {
2174 int* flag = p;
2175 if ((*flag == 0) && (a(i,j,k) == value)) {
2176 if (Gpu::Atomic::Exch(flag,1) == 0) {
2177 AMREX_D_TERM(p[1] = i;,
2178 p[2] = j;,
2179 p[3] = k;);
2180 }
2181 }
2182 });
2183 Gpu::dtoh_memcpy_async(ha.data(), p, sizeof(int)*ha.size());
2185 return IntVect(AMREX_D_DECL(ha[1],ha[2],ha[3]));
2186 } else
2187#endif
2188 {
2189 AMREX_LOOP_3D(subbox, i, j, k,
2190 {
2191 if (a(i,j,k) == value) { return IntVect(AMREX_D_DECL(i,j,k)); }
2192 });
2193 return IntVect::TheMinVector();
2194 }
2195}
2196
2197template <class T>
2198template <RunOn run_on>
2199IntVect
2200BaseFab<T>::minIndex (int comp) const noexcept
2201{
2202 return this->minIndex<run_on>(this->domain,comp);
2203}
2204
2205template <class T>
2206template <RunOn run_on>
2207IntVect
2208BaseFab<T>::minIndex (const Box& subbox, int comp) const noexcept
2209{
2210 T min_val = this->min<run_on>(subbox, comp);
2211 return this->indexFromValue<run_on>(subbox, comp, min_val);
2212}
2213
2214template <class T>
2215template <RunOn run_on>
2216void
2217BaseFab<T>::minIndex (const Box& subbox, Real& min_val, IntVect& min_idx, int comp) const noexcept
2218{
2219 min_val = this->min<run_on>(subbox, comp);
2220 min_idx = this->indexFromValue<run_on>(subbox, comp, min_val);
2221}
2222
2223template <class T>
2224template <RunOn run_on>
2225IntVect
2226BaseFab<T>::maxIndex (int comp) const noexcept
2227{
2228 return this->maxIndex<run_on>(this->domain,comp);
2229}
2230
2231template <class T>
2232template <RunOn run_on>
2233IntVect
2234BaseFab<T>::maxIndex (const Box& subbox, int comp) const noexcept
2235{
2236 T max_val = this->max<run_on>(subbox, comp);
2237 return this->indexFromValue<run_on>(subbox, comp, max_val);
2238}
2239
2240template <class T>
2241template <RunOn run_on>
2242void
2243BaseFab<T>::maxIndex (const Box& subbox, Real& max_val, IntVect& max_idx, int comp) const noexcept
2244{
2245 max_val = this->max<run_on>(subbox, comp);
2246 max_idx = this->indexFromValue<run_on>(subbox, comp, max_val);
2247}
2248
2249template <class T>
2250template <RunOn run_on>
2251int
2252BaseFab<T>::maskLT (BaseFab<int>& mask, T const& val, int comp) const noexcept
2253{
2254 mask.resize(this->domain,1);
2255 int cnt = 0;
2256 Array4<int> const& m = mask.array();
2257 Array4<T const> const& a = this->const_array(comp);
2258#ifdef AMREX_USE_GPU
2259 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2260 ReduceOps<ReduceOpSum> reduce_op;
2261 ReduceData<int> reduce_data(reduce_op);
2262 using ReduceTuple = typename decltype(reduce_data)::Type;
2263 reduce_op.eval(this->domain, reduce_data,
2264 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2265 {
2266 int t;
2267 if (a(i,j,k) < val) {
2268 m(i,j,k) = 1;
2269 t = 1;
2270 } else {
2271 m(i,j,k) = 0;
2272 t = 0;
2273 }
2274 return {t};
2275 });
2276 ReduceTuple hv = reduce_data.value(reduce_op);
2277 cnt = amrex::get<0>(hv);
2278 } else
2279#endif
2280 {
2281 AMREX_LOOP_3D(this->domain, i, j, k,
2282 {
2283 if (a(i,j,k) < val) {
2284 m(i,j,k) = 1;
2285 ++cnt;
2286 } else {
2287 m(i,j,k) = 0;
2288 }
2289 });
2290 }
2291
2292 return cnt;
2293}
2294
2295template <class T>
2296template <RunOn run_on>
2297int
2298BaseFab<T>::maskLE (BaseFab<int>& mask, T const& val, int comp) const noexcept
2299{
2300 mask.resize(this->domain,1);
2301 int cnt = 0;
2302 Array4<int> const& m = mask.array();
2303 Array4<T const> const& a = this->const_array(comp);
2304#ifdef AMREX_USE_GPU
2305 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2306 ReduceOps<ReduceOpSum> reduce_op;
2307 ReduceData<int> reduce_data(reduce_op);
2308 using ReduceTuple = typename decltype(reduce_data)::Type;
2309 reduce_op.eval(this->domain, reduce_data,
2310 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2311 {
2312 int t;
2313 if (a(i,j,k) <= val) {
2314 m(i,j,k) = 1;
2315 t = 1;
2316 } else {
2317 m(i,j,k) = 0;
2318 t = 0;
2319 }
2320 return {t};
2321 });
2322 ReduceTuple hv = reduce_data.value(reduce_op);
2323 cnt = amrex::get<0>(hv);
2324 } else
2325#endif
2326 {
2327 AMREX_LOOP_3D(this->domain, i, j, k,
2328 {
2329 if (a(i,j,k) <= val) {
2330 m(i,j,k) = 1;
2331 ++cnt;
2332 } else {
2333 m(i,j,k) = 0;
2334 }
2335 });
2336 }
2337
2338 return cnt;
2339}
2340
2341template <class T>
2342template <RunOn run_on>
2343int
2344BaseFab<T>::maskEQ (BaseFab<int>& mask, T const& val, int comp) const noexcept
2345{
2346 mask.resize(this->domain,1);
2347 int cnt = 0;
2348 Array4<int> const& m = mask.array();
2349 Array4<T const> const& a = this->const_array(comp);
2350#ifdef AMREX_USE_GPU
2351 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2352 ReduceOps<ReduceOpSum> reduce_op;
2353 ReduceData<int> reduce_data(reduce_op);
2354 using ReduceTuple = typename decltype(reduce_data)::Type;
2355 reduce_op.eval(this->domain, reduce_data,
2356 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2357 {
2358 int t;
2359 if (a(i,j,k) == val) {
2360 m(i,j,k) = 1;
2361 t = 1;
2362 } else {
2363 m(i,j,k) = 0;
2364 t = 0;
2365 }
2366 return {t};
2367 });
2368 ReduceTuple hv = reduce_data.value(reduce_op);
2369 cnt = amrex::get<0>(hv);
2370 } else
2371#endif
2372 {
2373 AMREX_LOOP_3D(this->domain, i, j, k,
2374 {
2375 if (a(i,j,k) == val) {
2376 m(i,j,k) = 1;
2377 ++cnt;
2378 } else {
2379 m(i,j,k) = 0;
2380 }
2381 });
2382 }
2383
2384 return cnt;
2385}
2386
2387template <class T>
2388template <RunOn run_on>
2389int
2390BaseFab<T>::maskGT (BaseFab<int>& mask, T const& val, int comp) const noexcept
2391{
2392 mask.resize(this->domain,1);
2393 int cnt = 0;
2394 Array4<int> const& m = mask.array();
2395 Array4<T const> const& a = this->const_array(comp);
2396#ifdef AMREX_USE_GPU
2397 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2398 ReduceOps<ReduceOpSum> reduce_op;
2399 ReduceData<int> reduce_data(reduce_op);
2400 using ReduceTuple = typename decltype(reduce_data)::Type;
2401 reduce_op.eval(this->domain, reduce_data,
2402 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2403 {
2404 int t;
2405 if (a(i,j,k) > val) {
2406 m(i,j,k) = 1;
2407 t = 1;
2408 } else {
2409 m(i,j,k) = 0;
2410 t = 0;
2411 }
2412 return {t};
2413 });
2414 ReduceTuple hv = reduce_data.value(reduce_op);
2415 cnt = amrex::get<0>(hv);
2416 } else
2417#endif
2418 {
2419 AMREX_LOOP_3D(this->domain, i, j, k,
2420 {
2421 if (a(i,j,k) > val) {
2422 m(i,j,k) = 1;
2423 ++cnt;
2424 } else {
2425 m(i,j,k) = 0;
2426 }
2427 });
2428 }
2429
2430 return cnt;
2431}
2432
2433template <class T>
2434template <RunOn run_on>
2435int
2436BaseFab<T>::maskGE (BaseFab<int>& mask, T const& val, int comp) const noexcept
2437{
2438 mask.resize(this->domain,1);
2439 int cnt = 0;
2440 Array4<int> const& m = mask.array();
2441 Array4<T const> const& a = this->const_array(comp);
2442#ifdef AMREX_USE_GPU
2443 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2444 ReduceOps<ReduceOpSum> reduce_op;
2445 ReduceData<int> reduce_data(reduce_op);
2446 using ReduceTuple = typename decltype(reduce_data)::Type;
2447 reduce_op.eval(this->domain, reduce_data,
2448 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2449 {
2450 int t;
2451 if (a(i,j,k) >= val) {
2452 m(i,j,k) = 1;
2453 t = 1;
2454 } else {
2455 m(i,j,k) = 0;
2456 t = 0;
2457 }
2458 return {t};
2459 });
2460 ReduceTuple hv = reduce_data.value(reduce_op);
2461 cnt = amrex::get<0>(hv);
2462 } else
2463#endif
2464 {
2465 AMREX_LOOP_3D(this->domain, i, j, k,
2466 {
2467 if (a(i,j,k) >= val) {
2468 m(i,j,k) = 1;
2469 ++cnt;
2470 } else {
2471 m(i,j,k) = 0;
2472 }
2473 });
2474 }
2475
2476 return cnt;
2477}
2478
2479template <class T>
2480template <RunOn run_on>
2483{
2484 Box ovlp(this->domain);
2485 ovlp &= x.domain;
2486 return ovlp.ok() ? this->atomicAdd<run_on>(x,ovlp,ovlp,0,0,this->nvar) : *this;
2487}
2488
2489template <class T>
2490template <RunOn run_on>
2492BaseFab<T>::saxpy (T a, const BaseFab<T>& x, const Box& srcbox, const Box& destbox,
2493 int srccomp, int destcomp, int numcomp) noexcept
2494{
2495 BL_ASSERT(srcbox.ok());
2496 BL_ASSERT(x.box().contains(srcbox));
2497 BL_ASSERT(destbox.ok());
2498 BL_ASSERT(box().contains(destbox));
2499 BL_ASSERT(destbox.sameSize(srcbox));
2500 BL_ASSERT( srccomp >= 0 && srccomp+numcomp <= x.nComp());
2501 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2502
2503 Array4<T> const& d = this->array();
2504 Array4<T const> const& s = x.const_array();
2505 const auto dlo = amrex::lbound(destbox);
2506 const auto slo = amrex::lbound(srcbox);
2507 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
2508 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2509 {
2510 d(i,j,k,n+destcomp) += a * s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
2511 });
2512
2513 return *this;
2514}
2515
2516template <class T>
2517template <RunOn run_on>
2519BaseFab<T>::saxpy (T a, const BaseFab<T>& x) noexcept
2520{
2521 Box ovlp(this->domain);
2522 ovlp &= x.domain;
2523 return ovlp.ok() ? saxpy<run_on>(a,x,ovlp,ovlp,0,0,this->nvar) : *this;
2524}
2525
2526template <class T>
2527template <RunOn run_on>
2530 const Box& srcbox, const Box& destbox,
2531 int srccomp, int destcomp, int numcomp) noexcept
2532{
2533 BL_ASSERT(srcbox.ok());
2534 BL_ASSERT(x.box().contains(srcbox));
2535 BL_ASSERT(destbox.ok());
2536 BL_ASSERT(box().contains(destbox));
2537 BL_ASSERT(destbox.sameSize(srcbox));
2538 BL_ASSERT( srccomp >= 0 && srccomp+numcomp <= x.nComp());
2539 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2540
2541 Array4<T> const& d = this->array();
2542 Array4<T const> const& s = x.const_array();
2543 const auto dlo = amrex::lbound(destbox);
2544 const auto slo = amrex::lbound(srcbox);
2545 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
2546 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2547 {
2548 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);
2549 });
2550
2551 return *this;
2552}
2553
2554template <class T>
2555template <RunOn run_on>
2557BaseFab<T>::addproduct (const Box& destbox, int destcomp, int numcomp,
2558 const BaseFab<T>& src1, int comp1,
2559 const BaseFab<T>& src2, int comp2) noexcept
2560{
2561 BL_ASSERT(destbox.ok());
2562 BL_ASSERT(box().contains(destbox));
2563 BL_ASSERT( comp1 >= 0 && comp1+numcomp <= src1.nComp());
2564 BL_ASSERT( comp2 >= 0 && comp2+numcomp <= src2.nComp());
2565 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2566
2567 Array4<T> const& d = this->array();
2568 Array4<T const> const& s1 = src1.const_array();
2569 Array4<T const> const& s2 = src2.const_array();
2570 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2571 {
2572 d(i,j,k,n+destcomp) += s1(i,j,k,n+comp1) * s2(i,j,k,n+comp2);
2573 });
2574
2575 return *this;
2576}
2577
2578template <class T>
2579template <RunOn run_on>
2581BaseFab<T>::linComb (const BaseFab<T>& f1, const Box& b1, int comp1,
2582 const BaseFab<T>& f2, const Box& b2, int comp2,
2583 Real alpha, Real beta, const Box& b,
2584 int comp, int numcomp) noexcept
2585{
2586 BL_ASSERT(b1.ok());
2587 BL_ASSERT(f1.box().contains(b1));
2588 BL_ASSERT(b2.ok());
2589 BL_ASSERT(f2.box().contains(b2));
2590 BL_ASSERT(b.ok());
2591 BL_ASSERT(box().contains(b));
2592 BL_ASSERT(b.sameSize(b1));
2593 BL_ASSERT(b.sameSize(b2));
2594 BL_ASSERT(comp1 >= 0 && comp1+numcomp <= f1.nComp());
2595 BL_ASSERT(comp2 >= 0 && comp2+numcomp <= f2.nComp());
2596 BL_ASSERT(comp >= 0 && comp +numcomp <= nComp());
2597
2598 Array4<T> const& d = this->array();
2599 Array4<T const> const& s1 = f1.const_array();
2600 Array4<T const> const& s2 = f2.const_array();
2601 const auto dlo = amrex::lbound(b);
2602 const auto slo1 = amrex::lbound(b1);
2603 const auto slo2 = amrex::lbound(b2);
2604 const Dim3 off1{slo1.x-dlo.x,slo1.y-dlo.y,slo1.z-dlo.z};
2605 const Dim3 off2{slo2.x-dlo.x,slo2.y-dlo.y,slo2.z-dlo.z};
2606
2607 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, b, numcomp, i, j, k, n,
2608 {
2609 d(i,j,k,n+comp) = alpha*s1(i+off1.x,j+off1.y,k+off1.z,n+comp1)
2610 + beta*s2(i+off2.x,j+off2.y,k+off2.z,n+comp2);
2611 });
2612 return *this;
2613}
2614
2615template <class T>
2616template <RunOn run_on>
2617T
2618BaseFab<T>::dot (const Box& xbx, int xcomp,
2619 const BaseFab<T>& y, const Box& ybx, int ycomp,
2620 int numcomp) const noexcept
2621{
2622 BL_ASSERT(xbx.ok());
2623 BL_ASSERT(box().contains(xbx));
2624 BL_ASSERT(y.box().contains(ybx));
2625 BL_ASSERT(xbx.sameSize(ybx));
2626 BL_ASSERT(xcomp >= 0 && xcomp+numcomp <= nComp());
2627 BL_ASSERT(ycomp >= 0 && ycomp+numcomp <= y.nComp());
2628
2629 T r = 0;
2630
2631 const auto xlo = amrex::lbound(xbx);
2632 const auto ylo = amrex::lbound(ybx);
2633 const Dim3 offset{ylo.x-xlo.x,ylo.y-xlo.y,ylo.z-xlo.z};
2634 Array4<T const> const& xa = this->const_array();
2635 Array4<T const> const& ya = y.const_array();
2636
2637#ifdef AMREX_USE_GPU
2638 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2639 ReduceOps<ReduceOpSum> reduce_op;
2640 ReduceData<T> reduce_data(reduce_op);
2641 using ReduceTuple = typename decltype(reduce_data)::Type;
2642 reduce_op.eval(xbx, reduce_data,
2643 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2644 {
2645 T t = 0;
2646 for (int n = 0; n < numcomp; ++n) {
2647 t += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp);
2648 }
2649 return {t};
2650 });
2651 ReduceTuple hv = reduce_data.value(reduce_op);
2652 r = amrex::get<0>(hv);
2653 } else
2654#endif
2655 {
2656 AMREX_LOOP_4D(xbx, numcomp, i, j, k, n,
2657 {
2658 r += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp);
2659 });
2660 }
2661
2662 return r;
2663}
2664
2665template <class T>
2666template <RunOn run_on>
2667T
2668BaseFab<T>::dotmask (const BaseFab<int>& mask, const Box& xbx, int xcomp,
2669 const BaseFab<T>& y, const Box& ybx, int ycomp,
2670 int numcomp) const noexcept
2671{
2672 BL_ASSERT(xbx.ok());
2673 BL_ASSERT(box().contains(xbx));
2674 BL_ASSERT(y.box().contains(ybx));
2675 BL_ASSERT(xbx.sameSize(ybx));
2676 BL_ASSERT(xcomp >= 0 && xcomp+numcomp <= nComp());
2677 BL_ASSERT(ycomp >= 0 && ycomp+numcomp <= y.nComp());
2678
2679 T r = 0;
2680
2681 const auto xlo = amrex::lbound(xbx);
2682 const auto ylo = amrex::lbound(ybx);
2683 const Dim3 offset{ylo.x-xlo.x,ylo.y-xlo.y,ylo.z-xlo.z};
2684
2685 Array4<T const> const& xa = this->const_array();
2686 Array4<T const> const& ya = y.const_array();
2687 Array4<int const> const& ma = mask.const_array();
2688
2689#ifdef AMREX_USE_GPU
2690 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
2691 ReduceOps<ReduceOpSum> reduce_op;
2692 ReduceData<T> reduce_data(reduce_op);
2693 using ReduceTuple = typename decltype(reduce_data)::Type;
2694 reduce_op.eval(xbx, reduce_data,
2695 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
2696 {
2697 int m = static_cast<int>(static_cast<bool>(ma(i,j,k)));
2698 T t = 0;
2699 for (int n = 0; n < numcomp; ++n) {
2700 t += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp) * m;
2701 }
2702 return {t};
2703 });
2704 ReduceTuple hv = reduce_data.value(reduce_op);
2705 r = amrex::get<0>(hv);
2706 } else
2707#endif
2708 {
2709 AMREX_LOOP_4D(xbx, numcomp, i, j, k, n,
2710 {
2711 int m = static_cast<int>(static_cast<bool>(ma(i,j,k)));
2712 r += xa(i,j,k,n+xcomp) * ya(i+offset.x,j+offset.y,k+offset.z,n+ycomp) * m;
2713 });
2714 }
2715
2716 return r;
2717}
2718
2719template <class T>
2720template <RunOn run_on>
2721T
2722BaseFab<T>::sum (int comp, int numcomp) const noexcept
2723{
2724 return this->sum<run_on>(this->domain, DestComp{comp}, NumComps{numcomp});
2725}
2726
2727template <class T>
2728template <RunOn run_on>
2729T
2730BaseFab<T>::sum (const Box& subbox, int comp, int numcomp) const noexcept
2731{
2732 return this->sum<run_on>(subbox, DestComp{comp}, NumComps{numcomp});
2733}
2734
2735template <class T>
2736template <RunOn run_on>
2738BaseFab<T>::negate (int comp, int numcomp) noexcept
2739{
2740 return this->negate<run_on>(this->domain, DestComp{comp}, NumComps{numcomp});
2741}
2742
2743template <class T>
2744template <RunOn run_on>
2746BaseFab<T>::negate (const Box& b, int comp, int numcomp) noexcept
2747{
2748 return this->negate<run_on>(b, DestComp{comp}, NumComps{numcomp});
2749}
2750
2751template <class T>
2752template <RunOn run_on>
2754BaseFab<T>::invert (T const& r, int comp, int numcomp) noexcept
2755{
2756 return this->invert<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
2757}
2758
2759template <class T>
2760template <RunOn run_on>
2762BaseFab<T>::invert (T const& r, const Box& b, int comp, int numcomp) noexcept
2763{
2764 return this->invert<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
2765}
2766
2767template <class T>
2768template <RunOn run_on>
2770BaseFab<T>::plus (T const& r, int comp, int numcomp) noexcept
2771{
2772 return this->plus<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
2773}
2774
2775template <class T>
2776template <RunOn run_on>
2778BaseFab<T>::plus (T const& r, const Box& b, int comp, int numcomp) noexcept
2779{
2780 return this->plus<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
2781}
2782
2783template <class T>
2784template <RunOn run_on>
2786BaseFab<T>::plus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
2787{
2788 return this->plus<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
2789}
2790
2791template <class T>
2792template <RunOn run_on>
2794BaseFab<T>::atomicAdd (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
2795{
2796 Box ovlp(this->domain);
2797 ovlp &= src.domain;
2798 return ovlp.ok() ? this->atomicAdd<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
2799}
2800
2801template <class T>
2802template <RunOn run_on>
2804BaseFab<T>::plus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
2805 int numcomp) noexcept
2806{
2807 return this->plus<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
2808}
2809
2810template <class T>
2811template <RunOn run_on>
2813BaseFab<T>::atomicAdd (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
2814 int numcomp) noexcept
2815{
2816 Box ovlp(this->domain);
2817 ovlp &= src.domain;
2818 ovlp &= subbox;
2819 return ovlp.ok() ? this->atomicAdd<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
2820}
2821
2822template <class T>
2823template <RunOn run_on>
2825BaseFab<T>::plus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
2826 int srccomp, int destcomp, int numcomp) noexcept
2827{
2828 BL_ASSERT(destbox.ok());
2829 BL_ASSERT(src.box().contains(srcbox));
2830 BL_ASSERT(box().contains(destbox));
2831 BL_ASSERT(destbox.sameSize(srcbox));
2832 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
2833 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2834
2835 Array4<T> const& d = this->array();
2836 Array4<T const> const& s = src.const_array();
2837 const auto dlo = amrex::lbound(destbox);
2838 const auto slo = amrex::lbound(srcbox);
2839 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
2840 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2841 {
2842 d(i,j,k,n+destcomp) += s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
2843 });
2844
2845 return *this;
2846}
2847
2848namespace detail
2849{
2850
2851template <RunOn run_on, typename T,
2852 std::enable_if_t<HasAtomicAdd<T>::value,int> FOO = 0>
2853void basefab_atomic_add (BaseFab<T>& dfab, const BaseFab<T>& sfab,
2854 const Box& srcbox, const Box& destbox,
2855 int srccomp, int destcomp, int numcomp) noexcept
2856{
2857 Array4<T> const& d = dfab.array();
2858 Array4<T const> const& s = sfab.const_array();
2859 const auto dlo = amrex::lbound(destbox);
2860 const auto slo = amrex::lbound(srcbox);
2861 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
2862 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
2863 {
2864 T* p = d.ptr(i,j,k,n+destcomp);
2865 HostDevice::Atomic::Add(p, s(i+offset.x,j+offset.y,k+offset.z,n+srccomp));
2866 });
2867}
2868
2869template <RunOn run_on, typename T,
2870 std::enable_if_t<! HasAtomicAdd<T>::value,int> FOO = 0>
2871void basefab_atomic_add (BaseFab<T>& dfab, const BaseFab<T>& sfab,
2872 const Box& srcbox, const Box& destbox,
2873 int srccomp, int destcomp, int numcomp) noexcept
2874{
2875 amrex::ignore_unused(dfab, sfab, srcbox, destbox, srccomp, destcomp, numcomp);
2876 amrex::Abort("BaseFab: atomicAdd not supported");
2877}
2878
2879}
2880
2881template <class T>
2882template <RunOn run_on>
2883BaseFab<T>&
2884BaseFab<T>::atomicAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
2885 int srccomp, int destcomp, int numcomp) noexcept
2886{
2887 BL_ASSERT(destbox.ok());
2888 BL_ASSERT(src.box().contains(srcbox));
2889 BL_ASSERT(box().contains(destbox));
2890 BL_ASSERT(destbox.sameSize(srcbox));
2891 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
2892 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2893
2894 detail::basefab_atomic_add<run_on>(*this, src, srcbox, destbox,
2895 srccomp, destcomp, numcomp);
2896
2897 return *this;
2898}
2899
2900template <class T>
2901template <RunOn run_on>
2903BaseFab<T>::lockAdd (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
2904 int srccomp, int destcomp, int numcomp) noexcept
2905{
2906#if defined(AMREX_USE_OMP) && (AMREX_SPACEDIM > 1)
2907#if defined(AMREX_USE_GPU)
2908 if (run_on == RunOn::Host || Gpu::notInLaunchRegion()) {
2909#endif
2910 BL_ASSERT(destbox.ok());
2911 BL_ASSERT(src.box().contains(srcbox));
2912 BL_ASSERT(box().contains(destbox));
2913 BL_ASSERT(destbox.sameSize(srcbox));
2914 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
2915 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
2916
2917 Array4<T> const& d = this->array();
2918 Array4<T const> const& s = src.const_array();
2919 auto const& dlo = amrex::lbound(destbox);
2920 auto const& dhi = amrex::ubound(destbox);
2921 auto const& len = amrex::length(destbox);
2922 auto const& slo = amrex::lbound(srcbox);
2923 Dim3 const offset{slo.x-dlo.x, slo.y-dlo.y, slo.z-dlo.z};
2924
2925 int planedim;
2926 int nplanes;
2927 int plo;
2928 if (len.z == 1) {
2929 planedim = 1;
2930 nplanes = len.y;
2931 plo = dlo.y;
2932 } else {
2933 planedim = 2;
2934 nplanes = len.z;
2935 plo = dlo.z;
2936 }
2937
2938 auto* mask = (bool*) amrex_mempool_alloc(sizeof(bool)*nplanes);
2939 for (int ip = 0; ip < nplanes; ++ip) {
2940 mask[ip] = false;
2941 }
2942
2943 int mm = 0;
2944 int planes_left = nplanes;
2945 while (planes_left > 0) {
2946 AMREX_ASSERT(mm < nplanes);
2947 auto const m = mm + plo;
2948 auto* lock = OpenMP::get_lock(m);
2949 if (omp_test_lock(lock))
2950 {
2951 auto lo = dlo;
2952 auto hi = dhi;
2953 if (planedim == 1) {
2954 lo.y = m;
2955 hi.y = m;
2956 } else {
2957 lo.z = m;
2958 hi.z = m;
2959 }
2960
2961 for (int n = 0; n < numcomp; ++n) {
2962 for (int k = lo.z; k <= hi.z; ++k) {
2963 for (int j = lo.y; j <= hi.y; ++j) {
2964 auto * pdst = d.ptr(dlo.x,j ,k ,n+destcomp);
2965 auto const* psrc = s.ptr(slo.x,j+offset.y,k+offset.z,n+ srccomp);
2966#pragma omp simd
2967 for (int ii = 0; ii < len.x; ++ii) {
2968 pdst[ii] += psrc[ii];
2969 }
2970 }
2971 }
2972 }
2973
2974 mask[mm] = true;
2975 --planes_left;
2976 omp_unset_lock(lock);
2977 if (planes_left == 0) { break; }
2978 }
2979
2980 ++mm;
2981 for (int ip = 0; ip < nplanes; ++ip) {
2982 int new_mm = (mm+ip) % nplanes;
2983 if ( ! mask[new_mm] ) {
2984 mm = new_mm;
2985 break;
2986 }
2987 }
2988 }
2989
2991
2992 return *this;
2993
2994#if defined(AMREX_USE_GPU)
2995 } else {
2996 return this->template atomicAdd<run_on>(src, srcbox, destbox, srccomp, destcomp, numcomp);
2997 }
2998#endif
2999#else
3000 return this->template atomicAdd<run_on>(src, srcbox, destbox, srccomp, destcomp, numcomp);
3001#endif
3002}
3003
3004template <class T>
3005template <RunOn run_on>
3007BaseFab<T>::minus (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
3008{
3009 return this->minus<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3010}
3011
3012template <class T>
3013template <RunOn run_on>
3015BaseFab<T>::minus (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
3016{
3017 return this->minus<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3018}
3019
3020template <class T>
3021template <RunOn run_on>
3023BaseFab<T>::minus (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
3024 int srccomp, int destcomp, int numcomp) noexcept
3025{
3026 BL_ASSERT(destbox.ok());
3027 BL_ASSERT(src.box().contains(srcbox));
3028 BL_ASSERT(box().contains(destbox));
3029 BL_ASSERT(destbox.sameSize(srcbox));
3030 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
3031 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
3032
3033 Array4<T> const& d = this->array();
3034 Array4<T const> const& s = src.const_array();
3035 const auto dlo = amrex::lbound(destbox);
3036 const auto slo = amrex::lbound(srcbox);
3037 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
3038 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
3039 {
3040 d(i,j,k,n+destcomp) -= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
3041 });
3042
3043 return *this;
3044}
3045
3046template <class T>
3047template <RunOn run_on>
3049BaseFab<T>::mult (T const& r, int comp, int numcomp) noexcept
3050{
3051 return this->mult<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
3052}
3053
3054template <class T>
3055template <RunOn run_on>
3057BaseFab<T>::mult (T const& r, const Box& b, int comp, int numcomp) noexcept
3058{
3059 return this->mult<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
3060}
3061
3062template <class T>
3063template <RunOn run_on>
3065BaseFab<T>::mult (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
3066{
3067 return this->mult<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3068}
3069
3070template <class T>
3071template <RunOn run_on>
3073BaseFab<T>::mult (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
3074{
3075 return this->mult<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3076}
3077
3078template <class T>
3079template <RunOn run_on>
3081BaseFab<T>::mult (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
3082 int srccomp, int destcomp, int numcomp) noexcept
3083{
3084 BL_ASSERT(destbox.ok());
3085 BL_ASSERT(src.box().contains(srcbox));
3086 BL_ASSERT(box().contains(destbox));
3087 BL_ASSERT(destbox.sameSize(srcbox));
3088 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
3089 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
3090
3091 Array4<T> const& d = this->array();
3092 Array4<T const> const& s = src.const_array();
3093 const auto dlo = amrex::lbound(destbox);
3094 const auto slo = amrex::lbound(srcbox);
3095 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
3096 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
3097 {
3098 d(i,j,k,n+destcomp) *= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
3099 });
3100
3101 return *this;
3102}
3103
3104template <class T>
3105template <RunOn run_on>
3107BaseFab<T>::divide (T const& r, int comp, int numcomp) noexcept
3108{
3109 return this->divide<run_on>(r, this->domain, DestComp{comp}, NumComps{numcomp});
3110}
3111
3112template <class T>
3113template <RunOn run_on>
3115BaseFab<T>::divide (T const& r, const Box& b, int comp, int numcomp) noexcept
3116{
3117 return this->divide<run_on>(r, b, DestComp{comp}, NumComps{numcomp});
3118}
3119
3120template <class T>
3121template <RunOn run_on>
3123BaseFab<T>::divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
3124{
3125 return this->divide<run_on>(src, this->domain, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3126}
3127
3128template <class T>
3129template <RunOn run_on>
3131BaseFab<T>::divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp, int numcomp) noexcept
3132{
3133 return this->divide<run_on>(src, subbox, SrcComp{srccomp}, DestComp{destcomp}, NumComps{numcomp});
3134}
3135
3136template <class T>
3137template <RunOn run_on>
3139BaseFab<T>::divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
3140 int srccomp, int destcomp, int numcomp) noexcept
3141{
3142 BL_ASSERT(destbox.ok());
3143 BL_ASSERT(src.box().contains(srcbox));
3144 BL_ASSERT(box().contains(destbox));
3145 BL_ASSERT(destbox.sameSize(srcbox));
3146 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
3147 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
3148
3149 Array4<T> const& d = this->array();
3150 Array4<T const> const& s = src.const_array();
3151 const auto dlo = amrex::lbound(destbox);
3152 const auto slo = amrex::lbound(srcbox);
3153 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
3154 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
3155 {
3156 d(i,j,k,n+destcomp) /= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
3157 });
3158
3159 return *this;
3160}
3161
3162template <class T>
3163template <RunOn run_on>
3166{
3167 Box ovlp(this->domain);
3168 ovlp &= src.domain;
3169 return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,0,0,this->nvar) : *this;
3170}
3171
3172template <class T>
3173template <RunOn run_on>
3175BaseFab<T>::protected_divide (const BaseFab<T>& src, int srccomp, int destcomp, int numcomp) noexcept
3176{
3177 Box ovlp(this->domain);
3178 ovlp &= src.domain;
3179 return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
3180}
3181
3182template <class T>
3183template <RunOn run_on>
3185BaseFab<T>::protected_divide (const BaseFab<T>& src, const Box& subbox, int srccomp, int destcomp,
3186 int numcomp) noexcept
3187{
3188 Box ovlp(this->domain);
3189 ovlp &= src.domain;
3190 ovlp &= subbox;
3191 return ovlp.ok() ? this->protected_divide<run_on>(src,ovlp,ovlp,srccomp,destcomp,numcomp) : *this;
3192}
3193
3194template <class T>
3195template <RunOn run_on>
3197BaseFab<T>::protected_divide (const BaseFab<T>& src, const Box& srcbox, const Box& destbox,
3198 int srccomp, int destcomp, int numcomp) noexcept
3199{
3200 BL_ASSERT(destbox.ok());
3201 BL_ASSERT(src.box().contains(srcbox));
3202 BL_ASSERT(box().contains(destbox));
3203 BL_ASSERT(destbox.sameSize(srcbox));
3204 BL_ASSERT(srccomp >= 0 && srccomp+numcomp <= src.nComp());
3205 BL_ASSERT(destcomp >= 0 && destcomp+numcomp <= nComp());
3206
3207 Array4<T> const& d = this->array();
3208 Array4<T const> const& s = src.const_array();
3209 const auto dlo = amrex::lbound(destbox);
3210 const auto slo = amrex::lbound(srcbox);
3211 const Dim3 offset{slo.x-dlo.x,slo.y-dlo.y,slo.z-dlo.z};
3212 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, destbox, numcomp, i, j, k, n,
3213 {
3214 if (s(i+offset.x,j+offset.y,k+offset.z,n+srccomp) != 0) {
3215 d(i,j,k,n+destcomp) /= s(i+offset.x,j+offset.y,k+offset.z,n+srccomp);
3216 }
3217 });
3218
3219 return *this;
3220}
3221
3232template <class T>
3233template <RunOn run_on>
3235BaseFab<T>::linInterp (const BaseFab<T>& f1, const Box& b1, int comp1,
3236 const BaseFab<T>& f2, const Box& b2, int comp2,
3237 Real t1, Real t2, Real t,
3238 const Box& b, int comp, int numcomp) noexcept
3239{
3240 if (amrex::almostEqual(t1,t2) || amrex::almostEqual(t1,t)) {
3241 return copy<run_on>(f1,b1,comp1,b,comp,numcomp);
3242 } else if (amrex::almostEqual(t2,t)) {
3243 return copy<run_on>(f2,b2,comp2,b,comp,numcomp);
3244 } else {
3245 Real alpha = (t2-t)/(t2-t1);
3246 Real beta = (t-t1)/(t2-t1);
3247 return linComb<run_on>(f1,b1,comp1,f2,b2,comp2,alpha,beta,b,comp,numcomp);
3248 }
3249}
3250
3251template <class T>
3252template <RunOn run_on>
3254BaseFab<T>::linInterp (const BaseFab<T>& f1, int comp1,
3255 const BaseFab<T>& f2, int comp2,
3256 Real t1, Real t2, Real t,
3257 const Box& b, int comp, int numcomp) noexcept
3258{
3259 if (amrex::almostEqual(t1,t2) || amrex::almostEqual(t1,t)) {
3260 return copy<run_on>(f1,b,comp1,b,comp,numcomp);
3261 } else if (amrex::almostEqual(t2,t)) {
3262 return copy<run_on>(f2,b,comp2,b,comp,numcomp);
3263 } else {
3264 Real alpha = (t2-t)/(t2-t1);
3265 Real beta = (t-t1)/(t2-t1);
3266 return linComb<run_on>(f1,b,comp1,f2,b,comp2,alpha,beta,b,comp,numcomp);
3267 }
3268}
3269
3270//
3271// New interfaces
3272//
3273
3274template <class T>
3275template <RunOn run_on>
3276void
3277BaseFab<T>::setVal (T const& val) noexcept
3278{
3279 this->setVal<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3280}
3281
3282template <class T>
3283template <RunOn run_on>
3284void
3285BaseFab<T>::setVal (T const& x, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3286{
3287 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3288 Array4<T> const& a = this->array();
3289 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
3290 {
3291 a(i,j,k,n+dcomp.i) = x;
3292 });
3293}
3294
3295template <class T>
3296template <RunOn run_on>
3297void
3298BaseFab<T>::setValIf (T const& val, const BaseFab<int>& mask) noexcept
3299{
3300 this->setValIf<run_on>(val, this->domain, mask, DestComp{0}, NumComps{this->nvar});
3301}
3302
3303template <class T>
3304template <RunOn run_on>
3305void
3306BaseFab<T>::setValIf (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept
3307{
3308 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3309 Array4<T> const& a = this->array();
3310 Array4<int const> const& m = mask.const_array();
3311 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
3312 {
3313 if (m(i,j,k)) { a(i,j,k,n+dcomp.i) = val; }
3314 });
3315}
3316
3317template <class T>
3318template <RunOn run_on>
3319void
3320BaseFab<T>::setValIfNot (T const& val, const BaseFab<int>& mask) noexcept
3321{
3322 this->setValIfNot<run_on>(val, this->domain, mask, DestComp{0}, NumComps{this->nvar});
3323}
3324
3325template <class T>
3326template <RunOn run_on>
3327void
3328BaseFab<T>::setValIfNot (T const& val, Box const& bx, const BaseFab<int>& mask, DestComp dcomp, NumComps ncomp) noexcept
3329{
3330 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3331 Array4<T> const& a = this->array();
3332 Array4<int const> const& m = mask.const_array();
3333 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG (run_on, bx, ncomp.n, i, j, k, n,
3334 {
3335 if (!m(i,j,k)) { a(i,j,k,n+dcomp.i) = val; }
3336 });
3337}
3338
3339template <class T>
3340template <RunOn run_on>
3341void
3342BaseFab<T>::setComplement (T const& x, const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
3343{
3344#ifdef AMREX_USE_GPU
3345 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3346 Array4<T> const& a = this->array();
3347 amrex::ParallelFor(this->domain, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
3348 {
3349 if (! bx.contains(IntVect(AMREX_D_DECL(i,j,k)))) {
3350 for (int n = dcomp.i; n < ncomp.n+dcomp.i; ++n) {
3351 a(i,j,k,n) = x;
3352 }
3353 }
3354 });
3355 } else
3356#endif
3357 {
3358 const BoxList b_lst = amrex::boxDiff(this->domain,bx);
3359 for (auto const& b : b_lst) {
3360 this->setVal<RunOn::Host>(x, b, dcomp, ncomp);
3361 }
3362 }
3363}
3364
3365template <class T>
3366template <RunOn run_on>
3368BaseFab<T>::copy (const BaseFab<T>& src) noexcept
3369{
3370 this->copy<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3371 return *this;
3372}
3373
3374template <class T>
3375template <RunOn run_on>
3378 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3379{
3380 AMREX_ASSERT(this->domain.sameType(src.domain));
3381 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3382 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3383
3384 bx &= src.domain;
3385
3386 Array4<T> const& d = this->array();
3387 Array4<T const> const& s = src.const_array();
3388 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3389 {
3390 d(i,j,k,n+dcomp.i) = s(i,j,k,n+scomp.i);
3391 });
3392
3393 return *this;
3394}
3395
3396template <class T>
3397template <RunOn run_on>
3399BaseFab<T>::plus (T const& val) noexcept
3400{
3401 return this->plus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3402}
3403
3404template <class T>
3405template <RunOn run_on>
3407BaseFab<T>::operator+= (T const& val) noexcept
3408{
3409 return this->plus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3410}
3411
3412template <class T>
3413template <RunOn run_on>
3415BaseFab<T>::plus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3416{
3417 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3418
3419 Array4<T> const& a = this->array();
3420 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3421 {
3422 a(i,j,k,n+dcomp.i) += val;
3423 });
3424
3425 return *this;
3426}
3427
3428template <class T>
3429template <RunOn run_on>
3431BaseFab<T>::plus (const BaseFab<T>& src) noexcept
3432{
3433 return this->plus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3434}
3435
3436template <class T>
3437template <RunOn run_on>
3440{
3441 return this->plus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3442}
3443
3444template <class T>
3445template <RunOn run_on>
3448 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3449{
3450 AMREX_ASSERT(this->domain.sameType(src.domain));
3451 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3452 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3453
3454 bx &= src.domain;
3455
3456 Array4<T> const& d = this->array();
3457 Array4<T const> const& s = src.const_array();
3458 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3459 {
3460 d(i,j,k,n+dcomp.i) += s(i,j,k,n+scomp.i);
3461 });
3462
3463 return *this;
3464}
3465
3466template <class T>
3467template <RunOn run_on>
3469BaseFab<T>::minus (T const& val) noexcept
3470{
3471 return this->minus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3472}
3473
3474template <class T>
3475template <RunOn run_on>
3477BaseFab<T>::operator-= (T const& val) noexcept
3478{
3479 return this->minus<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3480}
3481
3482template <class T>
3483template <RunOn run_on>
3485BaseFab<T>::minus (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3486{
3487 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3488
3489 Array4<T> const& a = this->array();
3490 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3491 {
3492 a(i,j,k,n+dcomp.i) -= val;
3493 });
3494
3495 return *this;
3496}
3497
3498template <class T>
3499template <RunOn run_on>
3501BaseFab<T>::minus (const BaseFab<T>& src) noexcept
3502{
3503 return this->minus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3504}
3505
3506template <class T>
3507template <RunOn run_on>
3510{
3511 return this->minus<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3512}
3513
3514template <class T>
3515template <RunOn run_on>
3518 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3519{
3520 AMREX_ASSERT(this->domain.sameType(src.domain));
3521 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3522 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3523
3524 bx &= src.domain;
3525
3526 Array4<T> const& d = this->array();
3527 Array4<T const> const& s = src.const_array();
3528 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3529 {
3530 d(i,j,k,n+dcomp.i) -= s(i,j,k,n+scomp.i);
3531 });
3532
3533 return *this;
3534}
3535
3536template <class T>
3537template <RunOn run_on>
3539BaseFab<T>::mult (T const& val) noexcept
3540{
3541 return this->mult<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3542}
3543
3544template <class T>
3545template <RunOn run_on>
3547BaseFab<T>::operator*= (T const& val) noexcept
3548{
3549 return this->mult<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3550}
3551
3552template <class T>
3553template <RunOn run_on>
3555BaseFab<T>::mult (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3556{
3557 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3558
3559 Array4<T> const& a = this->array();
3560 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3561 {
3562 a(i,j,k,n+dcomp.i) *= val;
3563 });
3564
3565 return *this;
3566}
3567
3568template <class T>
3569template <RunOn run_on>
3571BaseFab<T>::mult (const BaseFab<T>& src) noexcept
3572{
3573 return this->mult<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3574}
3575
3576template <class T>
3577template <RunOn run_on>
3580{
3581 return this->mult<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3582}
3583
3584template <class T>
3585template <RunOn run_on>
3588 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3589{
3590 AMREX_ASSERT(this->domain.sameType(src.domain));
3591 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3592 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3593
3594 bx &= src.domain;
3595
3596 Array4<T> const& d = this->array();
3597 Array4<T const> const& s = src.const_array();
3598 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3599 {
3600 d(i,j,k,n+dcomp.i) *= s(i,j,k,n+scomp.i);
3601 });
3602
3603 return *this;
3604}
3605
3606template <class T>
3607template <RunOn run_on>
3609BaseFab<T>::divide (T const& val) noexcept
3610{
3611 return this->divide<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3612}
3613
3614template <class T>
3615template <RunOn run_on>
3617BaseFab<T>::operator/= (T const& val) noexcept
3618{
3619 return this->divide<run_on>(val, this->domain, DestComp{0}, NumComps{this->nvar});
3620}
3621
3622template <class T>
3623template <RunOn run_on>
3625BaseFab<T>::divide (T const& val, Box const& bx, DestComp dcomp, NumComps ncomp) noexcept
3626{
3627 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3628
3629 Array4<T> const& a = this->array();
3630 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3631 {
3632 a(i,j,k,n+dcomp.i) /= val;
3633 });
3634
3635 return *this;
3636}
3637
3638template <class T>
3639template <RunOn run_on>
3641BaseFab<T>::divide (const BaseFab<T>& src) noexcept
3642{
3643 return this->divide<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3644}
3645
3646template <class T>
3647template <RunOn run_on>
3650{
3651 return this->divide<run_on>(src, this->domain, SrcComp{0}, DestComp{0}, NumComps{this->nvar});
3652}
3653
3654template <class T>
3655template <RunOn run_on>
3658 SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
3659{
3660 AMREX_ASSERT(this->domain.sameType(src.domain));
3661 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3662 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3663
3664 bx &= src.domain;
3665
3666 Array4<T> const& d = this->array();
3667 Array4<T const> const& s = src.const_array();
3668 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3669 {
3670 d(i,j,k,n+dcomp.i) /= s(i,j,k,n+scomp.i);
3671 });
3672
3673 return *this;
3674}
3675
3676template <class T>
3677template <RunOn run_on>
3680{
3681 return this->negate<run_on>(this->domain, DestComp{0}, NumComps{this->nvar});
3682}
3683
3684template <class T>
3685template <RunOn run_on>
3687BaseFab<T>::negate (const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
3688{
3689 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3690
3691 Array4<T> const& a = this->array();
3692 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3693 {
3694 a(i,j,k,n+dcomp.i) = -a(i,j,k,n+dcomp.i);
3695 });
3696
3697 return *this;
3698}
3699
3700template <class T>
3701template <RunOn run_on>
3703BaseFab<T>::invert (T const& r) noexcept
3704{
3705 return this->invert<run_on>(r, this->domain, DestComp{0}, NumComps{this->nvar});
3706}
3707
3708template <class T>
3709template <RunOn run_on>
3711BaseFab<T>::invert (T const& r, const Box& bx, DestComp dcomp, NumComps ncomp) noexcept
3712{
3713 BL_ASSERT(dcomp.i >= 0 && dcomp.i + ncomp.n <= this->nvar);
3714
3715 Array4<T> const& a = this->array();
3716 AMREX_HOST_DEVICE_PARALLEL_FOR_4D_FLAG(run_on, bx, ncomp.n, i, j, k, n,
3717 {
3718 a(i,j,k,n+dcomp.i) = r / a(i,j,k,n+dcomp.i);
3719 });
3720
3721 return *this;
3722}
3723
3724template <class T>
3725template <RunOn run_on>
3726T
3727BaseFab<T>::sum (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept
3728{
3729 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3730
3731 T r = 0;
3732 Array4<T const> const& a = this->const_array();
3733#ifdef AMREX_USE_GPU
3734 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3735 ReduceOps<ReduceOpSum> reduce_op;
3736 ReduceData<T> reduce_data(reduce_op);
3737 using ReduceTuple = typename decltype(reduce_data)::Type;
3738 reduce_op.eval(bx, reduce_data,
3739 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
3740 {
3741 T t = 0;
3742 for (int n = 0; n < ncomp.n; ++n) {
3743 t += a(i,j,k,n+dcomp.i);
3744 }
3745 return { t };
3746 });
3747 ReduceTuple hv = reduce_data.value(reduce_op);
3748 r = amrex::get<0>(hv);
3749 } else
3750#endif
3751 {
3752 amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
3753 {
3754 r += a(i,j,k,n+dcomp.i);
3755 });
3756 }
3757
3758 return r;
3759}
3760
3761template <class T>
3762template <RunOn run_on>
3763T
3764BaseFab<T>::dot (const BaseFab<T>& src, const Box& bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept
3765{
3766 AMREX_ASSERT(this->domain.sameType(src.domain));
3767 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3768 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3769
3770 T r = 0;
3771 Array4<T const> const& d = this->const_array();
3772 Array4<T const> const& s = src.const_array();
3773#ifdef AMREX_USE_GPU
3774 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3775 ReduceOps<ReduceOpSum> reduce_op;
3776 ReduceData<T> reduce_data(reduce_op);
3777 using ReduceTuple = typename decltype(reduce_data)::Type;
3778 reduce_op.eval(bx, reduce_data,
3779 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
3780 {
3781 T t = 0;
3782 for (int n = 0; n < ncomp.n; ++n) {
3783 t += d(i,j,k,n+dcomp.i) * s(i,j,k,n+scomp.i);
3784 }
3785 return { t };
3786 });
3787 ReduceTuple hv = reduce_data.value(reduce_op);
3788 r = amrex::get<0>(hv);
3789 } else
3790#endif
3791 {
3792 amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
3793 {
3794 r += d(i,j,k,n+dcomp.i) * s(i,j,k,n+scomp.i);
3795 });
3796 }
3797
3798 return r;
3799}
3800
3801template <class T>
3802template <RunOn run_on>
3803T
3804BaseFab<T>::dot (const Box& bx, int destcomp, int numcomp) const noexcept
3805{
3806 return dot<run_on>(bx, DestComp{destcomp}, NumComps{numcomp});
3807}
3808
3809
3810template <class T>
3811template <RunOn run_on>
3812T
3813BaseFab<T>::dot (const Box& bx, DestComp dcomp, NumComps ncomp) const noexcept
3814{
3815 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3816
3817 T r = 0;
3818 Array4<T const> const& a = this->const_array();
3819#ifdef AMREX_USE_GPU
3820 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3821 ReduceOps<ReduceOpSum> reduce_op;
3822 ReduceData<T> reduce_data(reduce_op);
3823 using ReduceTuple = typename decltype(reduce_data)::Type;
3824 reduce_op.eval(bx, reduce_data,
3825 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
3826 {
3827 T t = 0;
3828 for (int n = 0; n < ncomp.n; ++n) {
3829 t += a(i,j,k,n+dcomp.i)*a(i,j,k,n+dcomp.i);
3830 }
3831 return { t };
3832 });
3833 ReduceTuple hv = reduce_data.value(reduce_op);
3834 r = amrex::get<0>(hv);
3835 } else
3836#endif
3837 {
3838 amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
3839 {
3840 r += a(i,j,k,n+dcomp.i)*a(i,j,k,n+dcomp.i);
3841 });
3842 }
3843
3844 return r;
3845}
3846
3847template <class T>
3848template <RunOn run_on>
3849T
3850BaseFab<T>::dotmask (const BaseFab<T>& src, const Box& bx, const BaseFab<int>& mask,
3851 SrcComp scomp, DestComp dcomp, NumComps ncomp) const noexcept
3852{
3853 AMREX_ASSERT(this->domain.sameType(src.domain));
3854 AMREX_ASSERT(this->domain.sameType(mask.domain));
3855 AMREX_ASSERT(scomp.i >= 0 && scomp.i+ncomp.n <= src.nvar);
3856 AMREX_ASSERT(dcomp.i >= 0 && dcomp.i+ncomp.n <= this->nvar);
3857
3858 T r = 0;
3859 Array4<T const> const& d = this->const_array();
3860 Array4<T const> const& s = src.const_array();
3861 Array4<int const> const& m = mask.const_array();
3862#ifdef AMREX_USE_GPU
3863 if (run_on == RunOn::Device && Gpu::inLaunchRegion()) {
3864 ReduceOps<ReduceOpSum> reduce_op;
3865 ReduceData<T> reduce_data(reduce_op);
3866 using ReduceTuple = typename decltype(reduce_data)::Type;
3867 reduce_op.eval(bx, reduce_data,
3868 [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
3869 {
3870 T t = 0;
3871 T mi = static_cast<T>(static_cast<int>(static_cast<bool>(m(i,j,k))));
3872 for (int n = 0; n < ncomp.n; ++n) {
3873 t += d(i,j,k,n+dcomp.i)*s(i,j,k,n+scomp.i)*mi;
3874 }
3875 return { t };
3876 });
3877 ReduceTuple hv = reduce_data.value(reduce_op);
3878 r = amrex::get<0>(hv);
3879 } else
3880#endif
3881 {
3882 amrex::LoopOnCpu(bx, ncomp.n, [=,&r] (int i, int j, int k, int n) noexcept
3883 {
3884 int mi = static_cast<int>(static_cast<bool>(m(i,j,k)));
3885 r += d(i,j,k,n+dcomp.i)*s(i,j,k,n+scomp.i)*mi;
3886 });
3887 }
3888
3889 return r;
3890}
3891
3892}
3893
3894#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:1143
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1139
Real * pdst
Definition AMReX_HypreMLABecLap.cpp:1140
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:127
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:2436
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:2492
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:3377
T sum(int comp, int numcomp=1) const noexcept
Returns sum of given component of FAB state vector.
Definition AMReX_BaseFab.H:2722
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:1865
BaseFab< T > & divide(T const &val) noexcept
Scalar division on the whole domain and all components.
Definition AMReX_BaseFab.H:3609
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:3517
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:2903
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:1746
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:1801
std::size_t nBytesOwned() const noexcept
Definition AMReX_BaseFab.H:271
BaseFab< T > & copy(const BaseFab< T > &src) noexcept
Definition AMReX_BaseFab.H:3368
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:2557
BaseFab< T > & minus(T const &val) noexcept
Scalar subtraction on the whole domain and all components.
Definition AMReX_BaseFab.H:3469
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:2252
BaseFab< T > & plus(T const &val) noexcept
Scalar addition on the whole domain and all components.
Definition AMReX_BaseFab.H:3399
BaseFab< T > & mult(T const &val, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3555
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:3587
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:2581
void define()
Allocates memory for the BaseFab<T>.
Definition AMReX_BaseFab.H:1458
BaseFab< T > & operator*=(T const &val) noexcept
Definition AMReX_BaseFab.H:3547
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:3049
BaseFab< T > & atomicAdd(const BaseFab< T > &x) noexcept
Atomic FAB addition (a[i] <- a[i] + b[i]).
Definition AMReX_BaseFab.H:2482
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:2344
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:1728
void setVal(T const &x, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3285
BaseFab< T > & operator-=(T const &val) noexcept
Definition AMReX_BaseFab.H:3477
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:3306
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:2161
BaseFab< T > & mult(const BaseFab< T > &src) noexcept
Definition AMReX_BaseFab.H:3571
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:2298
void setValIf(T const &val, const BaseFab< int > &mask) noexcept
Definition AMReX_BaseFab.H:3298
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:3447
void setValIfNot(T const &val, const BaseFab< int > &mask) noexcept
Definition AMReX_BaseFab.H:3320
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:2529
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:1773
BaseFab< T > & negate() noexcept
on the whole domain and all components
Definition AMReX_BaseFab.H:3679
BaseFab< T > & minus(const BaseFab< T > &src) noexcept
Definition AMReX_BaseFab.H:3501
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:3007
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:2122
BaseFab< T > & operator+=(T const &val) noexcept
Definition AMReX_BaseFab.H:3407
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:3342
BaseFab< T > & minus(T const &val, Box const &bx, DestComp dcomp, NumComps ncomp) noexcept
Do nothing if bx is empty.
Definition AMReX_BaseFab.H:3485
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:3277
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:3625
T max(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2042
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:2618
BaseFab< T > & divide(T const &r, int comp, int numcomp=1) noexcept
As above except specify which components.
Definition AMReX_BaseFab.H:3107
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:1909
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:3617
std::pair< T, T > minmax(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2080
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:2778
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:3328
BaseFab< T > & mult(T const &val) noexcept
Scalar multiplication on the whole domain and all components.
Definition AMReX_BaseFab.H:3539
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:3657
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:2226
BaseFab< T > & protected_divide(const BaseFab< T > &src) noexcept
Divide wherever "src" is "true" or "non-zero".
Definition AMReX_BaseFab.H:3165
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:2762
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:2390
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:2668
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:3431
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:3641
void prefetchToHost() const noexcept
Definition AMReX_BaseFab.H:1203
T min(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2004
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:1829
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:3235
IntVect minIndex(int comp=0) const noexcept
Definition AMReX_BaseFab.H:2200
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:1837
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:3415
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 int deviceId() noexcept
Definition AMReX_GpuDevice.cpp:691
static int devicePropMajor() noexcept
Definition AMReX_GpuDevice.H:203
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:488
Definition AMReX_Reduce.H:612
std::enable_if_t< IsFabArray< MF >::value > eval(MF const &mf, IntVect const &nghost, D &reduce_data, F &&f)
Definition AMReX_Reduce.H:746
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:310
void dtoh_memcpy_async(void *p_h, const void *p_d, const std::size_t sz) noexcept
Definition AMReX_GpuDevice.H:435
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:421
gpuStream_t gpuStream() noexcept
Definition AMReX_GpuDevice.H:291
__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:183
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:234
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:240
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