Block-Structured AMR Software Framework
Loading...
Searching...
No Matches
AMReX_MultiFabUtil.H
Go to the documentation of this file.
1#ifndef AMREX_MultiFabUtil_H_
2#define AMREX_MultiFabUtil_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_Concepts.H>
6#include <AMReX_MultiFab.H>
7#include <AMReX_iMultiFab.H>
8#include <AMReX_LayoutData.H>
9#include <AMReX_MFIter.H>
10#include <AMReX_Array.H>
11#include <AMReX_Vector.H>
12#include <AMReX_MultiFabUtil_C.H>
13
14#include <AMReX_MultiFabUtilI.H>
15
16namespace amrex
17{
19 void average_node_to_cellcenter (MultiFab& cc, int dcomp,
20 const MultiFab& nd, int scomp,
21 int ncomp, int ngrow = 0);
22 void average_node_to_cellcenter (MultiFab& cc, int dcomp,
23 const MultiFab& nd, int scomp,
24 int ncomp, IntVect const& ng_vect);
25
32 void average_edge_to_cellcenter (MultiFab& cc, int dcomp,
33 const Vector<const MultiFab*>& edge,
34 int ngrow = 0);
35 void average_edge_to_cellcenter (MultiFab& cc, int dcomp,
36 const Vector<const MultiFab*>& edge,
37 IntVect const& ng_vect);
38
45 void average_face_to_cellcenter (MultiFab& cc, int dcomp,
46 const Vector<const MultiFab*>& fc,
47 int ngrow = 0);
48 void average_face_to_cellcenter (MultiFab& cc, int dcomp,
49 const Vector<const MultiFab*>& fc,
50 IntVect const& ng_vect);
51
53 template <FabArrayType CMF, FabArrayType FMF>
54 void average_face_to_cellcenter (CMF& cc, int dcomp,
55 const Array<const FMF*,AMREX_SPACEDIM>& fc,
56 int ngrow = 0);
57
58 template <FabArrayType CMF, FabArrayType FMF>
59 void average_face_to_cellcenter (CMF& cc, int dcomp,
60 const Array<const FMF*,AMREX_SPACEDIM>& fc,
61 IntVect const& ng_vect);
62
64 void average_face_to_cellcenter (MultiFab& cc,
65 const Vector<const MultiFab*>& fc,
66 const Geometry& geom);
68 void average_face_to_cellcenter (MultiFab& cc,
69 const Array<const MultiFab*,AMREX_SPACEDIM>& fc,
70 const Geometry& geom);
79 void average_cellcenter_to_face (const Vector<MultiFab*>& fc,
80 const MultiFab& cc,
81 const Geometry& geom,
82 int ncomp = 1,
83 bool use_harmonic_averaging = false,
84 int ngrow = 0);
86 void average_cellcenter_to_face (const Array<MultiFab*,AMREX_SPACEDIM>& fc,
87 const MultiFab& cc,
88 const Geometry& geom,
89 int ncomp = 1,
90 bool use_harmonic_averaging = false,
91 int ngrow = 0);
92
93 void average_cellcenter_to_face (const Array<MultiFab*,AMREX_SPACEDIM>& fc,
94 const MultiFab& cc,
95 const Geometry& geom,
96 int ncomp,
97 bool use_harmonic_averaging,
98 const Array<IntVect, AMREX_SPACEDIM>& ng_vects);
99
101 template <FabArrayType MF>
102 void average_down_faces (const Vector<const MF*>& fine,
103 const Vector<MF*>& crse,
104 const IntVect& ratio,
105 int ngcrse = 0);
107 template <FabArrayType MF>
108 void average_down_faces (const Vector<const MF*>& fine,
109 const Vector<MF*>& crse,
110 int ratio,
111 int ngcrse = 0);
113 template <FabArrayType MF>
114 void average_down_faces (const Array<const MF*,AMREX_SPACEDIM>& fine,
115 const Array<MF*,AMREX_SPACEDIM>& crse,
116 const IntVect& ratio,
117 int ngcrse = 0);
119 template <FabArrayType MF>
120 void average_down_faces (const Array<const MF*,AMREX_SPACEDIM>& fine,
121 const Array<MF*,AMREX_SPACEDIM>& crse,
122 int ratio,
123 int ngcrse = 0);
130 template <BaseFabType FAB>
131 void average_down_faces (const FabArray<FAB>& fine, FabArray<FAB>& crse,
132 const IntVect& ratio, int ngcrse=0);
133
134 // This version takes periodicity into account.
135 template <FabArrayType MF>
136 void average_down_faces (const Array<const MF*,AMREX_SPACEDIM>& fine,
137 const Array<MF*,AMREX_SPACEDIM>& crse,
138 const IntVect& ratio, const Geometry& crse_geom);
139 // This version takes periodicity into account.
140 template <BaseFabType FAB>
141 void average_down_faces (const FabArray<FAB>& fine, FabArray<FAB>& crse,
142 const IntVect& ratio, const Geometry& crse_geom);
143
145 void average_down_edges (const Vector<const MultiFab*>& fine,
146 const Vector<MultiFab*>& crse,
147 const IntVect& ratio,
148 int ngcrse = 0);
149 void average_down_edges (const Array<const MultiFab*,AMREX_SPACEDIM>& fine,
150 const Array<MultiFab*,AMREX_SPACEDIM>& crse,
151 const IntVect& ratio,
152 int ngcrse = 0);
156 void average_down_edges (const MultiFab& fine, MultiFab& crse,
157 const IntVect& ratio, int ngcrse=0);
158
160 template <BaseFabType FAB>
161 void average_down_nodal (const FabArray<FAB>& S_fine,
162 FabArray<FAB>& S_crse,
163 const IntVect& ratio,
164 int ngcrse = 0,
165 bool mfiter_is_definitely_safe=false);
166
173 void average_down (const MultiFab& S_fine, MultiFab& S_crse,
174 const Geometry& fgeom, const Geometry& cgeom,
175 int scomp, int ncomp, const IntVect& ratio);
176 void average_down (const MultiFab& S_fine, MultiFab& S_crse,
177 const Geometry& fgeom, const Geometry& cgeom,
178 int scomp, int ncomp, int rr);
179
183 template <BaseFabType FAB>
184 void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse,
185 int scomp, int ncomp, const IntVect& ratio);
186 template <BaseFabType FAB>
187 void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse,
188 int scomp, int ncomp, int rr);
189
192 void sum_fine_to_coarse (const MultiFab& S_Fine, MultiFab& S_crse,
193 int scomp, int ncomp,
194 const IntVect& ratio,
195 const Geometry& cgeom, const Geometry& fgeom);
196
198 void print_state (const MultiFab& mf, const IntVect& cell, int n=-1,
199 const IntVect& ng = IntVect::TheZeroVector());
200
202 void writeFabs (const MultiFab& mf, const std::string& name);
203 void writeFabs (const MultiFab& mf, int comp, int ncomp, const std::string& name);
204
207 std::unique_ptr<MultiFab> get_slice_data(int dir, Real coord,
208 const MultiFab& cc,
209 const Geometry& geom, int start_comp, int ncomp,
210 bool interpolate=false,
211 RealBox const& bnd_rbx = RealBox());
212
220 template <FabArrayType MF>
221 Vector<typename MF::value_type> get_cell_data (MF const& mf, IntVect const& cell);
222
229 template <FabArrayType MF>
230 MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box const& bnd_bx = Box());
231
235 template <BaseFabType FAB>
236 iMultiFab makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
237 int crse_value = 0, int fine_value = 1,
238 MFInfo const& info = MFInfo());
239 iMultiFab makeFineMask (const BoxArray& cba, const DistributionMapping& cdm,
240 const BoxArray& fba, const IntVect& ratio,
241 int crse_value = 0, int fine_value = 1,
242 MFInfo const& info = MFInfo());
243 template <BaseFabType FAB>
244 iMultiFab makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
245 Periodicity const& period, int crse_value, int fine_value,
246 MFInfo const& info = MFInfo());
247 iMultiFab makeFineMask (const BoxArray& cba, const DistributionMapping& cdm,
248 const IntVect& cnghost, const BoxArray& fba, const IntVect& ratio,
249 Periodicity const& period, int crse_value, int fine_value,
250 MFInfo const& info = MFInfo());
251 template <BaseFabType FAB>
252 iMultiFab makeFineMask (const FabArray<FAB>& cmf, const FabArray<FAB>& fmf,
253 const IntVect& cnghost, const IntVect& ratio,
254 Periodicity const& period, int crse_value, int fine_value,
255 MFInfo const& info = MFInfo());
256 template <BaseFabType FAB>
257 iMultiFab makeFineMask (const FabArray<FAB>& cmf, const FabArray<FAB>& fmf,
258 const IntVect& cnghost, const IntVect& ratio,
259 Periodicity const& period, int crse_value, int fine_value,
260 LayoutData<int>& has_cf,
261 MFInfo const& info = MFInfo());
262
263 MultiFab makeFineMask (const BoxArray& cba, const DistributionMapping& cdm,
264 const BoxArray& fba, const IntVect& ratio,
265 Real crse_value, Real fine_value,
266 MFInfo const& info = MFInfo());
267
269 void computeDivergence (MultiFab& divu, const Array<MultiFab const*,AMREX_SPACEDIM>& umac,
270 const Geometry& geom);
271
273 void computeGradient (MultiFab& grad, const Array<MultiFab const*,AMREX_SPACEDIM>& umac,
274 const Geometry& geom);
275
277 MultiFab ToMultiFab (const iMultiFab& imf);
279 FabArray<BaseFab<Long> > ToLongMultiFab (const iMultiFab& imf);
280
282 MultiFab periodicShift (MultiFab const& mf, IntVect const& offset,
283 Periodicity const& period);
284
286 template <typename T, typename U>
287 T cast (U const& mf_in)
288 {
289 T mf_out(mf_in.boxArray(), mf_in.DistributionMap(), mf_in.nComp(), mf_in.nGrowVect());
290
291#ifdef AMREX_USE_OMP
292#pragma omp parallel if (Gpu::notInLaunchRegion())
293#endif
294 for (MFIter mfi(mf_in); mfi.isValid(); ++mfi)
295 {
296 const Long n = mfi.fabbox().numPts() * mf_in.nComp();
297 auto * pdst = mf_out[mfi].dataPtr();
298 auto const* psrc = mf_in [mfi].dataPtr();
300 {
301 pdst[i] = static_cast<typename T::value_type>(psrc[i]); // NOLINT(bugprone-signed-char-misuse)
302 });
303 }
304 return mf_out;
305 }
306
367 template <typename Op, typename T, BaseFabType FAB, typename F>
368#ifndef AMREX_USE_CUDA
369 requires (IsCallableR<T,F,int,int,int,int>::value)
370#endif
371 BaseFab<T>
372 ReduceToPlane (int direction, Box const& domain, FabArray<FAB> const& mf, F const& f);
373
418 template <typename Op, MultiFabLike FA, typename F>
419#ifndef AMREX_USE_CUDA
420 requires (IsCallableR<typename FA::value_type,
421 F,int,int,int,int>::value)
422#endif
423 FA ReduceToPlaneMF (int direction, Box const& domain, FA const& mf, F const& f);
424
479 template <typename Op, MultiFabLike FA, typename F>
480#ifndef AMREX_USE_CUDA
481 requires (IsCallableR<typename FA::value_type,
482 F,int,int,int,int>::value)
483#endif
484 std::pair<FA,FA>
485 ReduceToPlaneMF2 (int direction, Box const& domain, FA const& mf, F const& f,
486 IntVect const& nghost = IntVect(0),
487 Periodicity const& period = Periodicity::NonPeriodic());
488
516 template <typename Op, MultiFabLike FA, typename F>
517#ifndef AMREX_USE_CUDA
518 requires (IsCallableR<typename FA::value_type,
519 F,int,int,int,int>::value)
520#endif
521 std::pair<FA,FA>
522 ReduceToPlaneMF2Patchy (int direction, Box const& domain, FA const& mf,
523 IntVect const& plane_max_grid_size, F const& f);
524
542 Gpu::HostVector<Real> sumToLine (MultiFab const& mf, int icomp, int ncomp,
543 Box const& domain, int direction, bool local = false);
544
552 Real volumeWeightedSum (Vector<MultiFab const*> const& mf, int icomp,
553 Vector<Geometry> const& geom,
554 Vector<IntVect> const& ratio,
555 bool local = false);
556
570 void FourthOrderInterpFromFineToCoarse (MultiFab& cmf, int scomp, int ncomp,
571 MultiFab const& fmf,
572 IntVect const& ratio);
573
584 void FillRandom (MultiFab& mf, int scomp, int ncomp);
585
597 void FillRandomNormal (MultiFab& mf, int scomp, int ncomp, Real mean, Real stddev);
598
610 [[nodiscard]] Vector<MultiFab> convexify (Vector<MultiFab const*> const& mf,
611 Vector<IntVect> const& refinement_ratio);
612}
613
614namespace amrex {
615
616template <BaseFabType FAB>
617iMultiFab
618makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
619 int crse_value, int fine_value, MFInfo const& info)
620{
621 return makeFineMask(cmf.boxArray(), cmf.DistributionMap(), cmf.nGrowVect(),
622 fba, ratio, Periodicity::NonPeriodic(), crse_value, fine_value,
623 info);
624}
625
626template <BaseFabType FAB>
627iMultiFab
628makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
629 Periodicity const& period, int crse_value, int fine_value,
630 MFInfo const& info)
631{
632 return makeFineMask(cmf.boxArray(), cmf.DistributionMap(), cmf.nGrowVect(),
633 fba, ratio, period, crse_value, fine_value, info);
634}
635
636template <BaseFabType FAB>
637iMultiFab
639 const IntVect& cnghost, const IntVect& ratio,
640 Periodicity const& period, int crse_value, int fine_value,
641 MFInfo const& info)
642{
643 iMultiFab mask(cmf.boxArray(), cmf.DistributionMap(), 1, cnghost, info);
644 mask.setVal(crse_value);
645
646 iMultiFab foo(amrex::coarsen(fmf.boxArray(),ratio), fmf.DistributionMap(),
647 1, 0, MFInfo().SetAlloc(false));
648 const FabArrayBase::CPC& cpc = mask.getCPC(cnghost,foo,IntVect::TheZeroVector(),period);
649 mask.setVal(fine_value, cpc, 0, 1);
650
651 return mask;
652}
653
654template <BaseFabType FAB>
655iMultiFab
657 const IntVect& cnghost, const IntVect& ratio,
658 Periodicity const& period, int crse_value, int fine_value,
659 LayoutData<int>& has_cf, MFInfo const& info)
660{
661 iMultiFab mask(cmf.boxArray(), cmf.DistributionMap(), 1, cnghost, info);
662 mask.setVal(crse_value);
663
664 iMultiFab foo(amrex::coarsen(fmf.boxArray(),ratio), fmf.DistributionMap(),
665 1, 0, MFInfo().SetAlloc(false));
666 const FabArrayBase::CPC& cpc = mask.getCPC(cnghost,foo,IntVect::TheZeroVector(),period);
667 mask.setVal(fine_value, cpc, 0, 1);
668
669 has_cf = mask.RecvLayoutMask(cpc);
670
671 return mask;
672}
673
676template <BaseFabType FAB>
678 const IntVect& ratio, int ngcrse, bool mfiter_is_definitely_safe)
679{
680 AMREX_ASSERT(fine.is_nodal());
681 AMREX_ASSERT(crse.is_nodal());
682 AMREX_ASSERT(crse.nComp() == fine.nComp());
683
684 int ncomp = crse.nComp();
685 using value_type = typename FAB::value_type;
686
687 if (mfiter_is_definitely_safe || isMFIterSafe(fine, crse))
688 {
689#ifdef AMREX_USE_OMP
690#pragma omp parallel if (Gpu::notInLaunchRegion())
691#endif
692 for (MFIter mfi(crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
693 {
694 const Box& bx = mfi.growntilebox(ngcrse);
695 Array4<value_type> const& crsearr = crse.array(mfi);
696 Array4<value_type const> const& finearr = fine.const_array(mfi);
697
699 {
700 amrex_avgdown_nodes(tbx,crsearr,finearr,0,0,ncomp,ratio);
701 });
702 }
703 }
704 else
705 {
706 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
707 ncomp, ngcrse, MFInfo().SetArena(The_Async_Arena()));
708 average_down_nodal(fine, ctmp, ratio, ngcrse);
709 crse.ParallelCopy(ctmp,0,0,ncomp,ngcrse,ngcrse);
710 }
711}
712
713// *************************************************************************************************************
714
715// Average fine cell-based MultiFab onto crse cell-centered MultiFab.
716// We do NOT assume that the coarse layout is a coarsened version of the fine layout.
717// This version does NOT use volume-weighting
718template <BaseFabType FAB>
719void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse, int scomp, int ncomp, int rr)
720{
721 average_down(S_fine,S_crse,scomp,ncomp,rr*IntVect::TheUnitVector());
722}
723
724template <BaseFabType FAB>
725void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse,
726 int scomp, int ncomp, const IntVect& ratio)
727{
728 BL_PROFILE("amrex::average_down");
729 AMREX_ASSERT(S_crse.nComp() == S_fine.nComp());
730 AMREX_ASSERT((S_crse.is_cell_centered() && S_fine.is_cell_centered()) ||
731 (S_crse.is_nodal() && S_fine.is_nodal()));
732
733 using value_type = typename FAB::value_type;
734
735 bool is_cell_centered = S_crse.is_cell_centered();
736
737 //
738 // Coarsen() the fine stuff on processors owning the fine data.
739 //
740 BoxArray crse_S_fine_BA = S_fine.boxArray(); crse_S_fine_BA.coarsen(ratio);
741
742 if (crse_S_fine_BA == S_crse.boxArray() && S_fine.DistributionMap() == S_crse.DistributionMap())
743 {
744#ifdef AMREX_USE_GPU
745 if (Gpu::inLaunchRegion() && S_crse.isFusingCandidate()) {
746 auto const& crsema = S_crse.arrays();
747 auto const& finema = S_fine.const_arrays();
748 if (is_cell_centered) {
749 ParallelFor(S_crse, IntVect(0), ncomp,
750 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
751 {
752 amrex_avgdown(i,j,k,n,crsema[box_no],finema[box_no],scomp,scomp,ratio);
753 });
754 } else {
755 ParallelFor(S_crse, IntVect(0), ncomp,
756 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
757 {
758 amrex_avgdown_nodes(i,j,k,n,crsema[box_no],finema[box_no],scomp,scomp,ratio);
759 });
760 }
761 if (!Gpu::inNoSyncRegion()) {
763 }
764 } else
765#endif
766 {
767#ifdef AMREX_USE_OMP
768#pragma omp parallel if (Gpu::notInLaunchRegion())
769#endif
770 for (MFIter mfi(S_crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
771 {
772 // NOTE: The tilebox is defined at the coarse level.
773 const Box& bx = mfi.tilebox();
774 Array4<value_type> const& crsearr = S_crse.array(mfi);
775 Array4<value_type const> const& finearr = S_fine.const_array(mfi);
776
777 if (is_cell_centered) {
778 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
779 {
780 amrex_avgdown(i,j,k,n,crsearr,finearr,scomp,scomp,ratio);
781 });
782 } else {
783 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
784 {
785 amrex_avgdown_nodes(i,j,k,n,crsearr,finearr,scomp,scomp,ratio);
786 });
787 }
788 }
789 }
790 }
791 else
792 {
793 FabArray<FAB> crse_S_fine(crse_S_fine_BA, S_fine.DistributionMap(), ncomp, 0,
795
796#ifdef AMREX_USE_GPU
797 if (Gpu::inLaunchRegion() && crse_S_fine.isFusingCandidate()) {
798 auto const& crsema = crse_S_fine.arrays();
799 auto const& finema = S_fine.const_arrays();
800 if (is_cell_centered) {
801 ParallelFor(crse_S_fine, IntVect(0), ncomp,
802 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
803 {
804 amrex_avgdown(i,j,k,n,crsema[box_no],finema[box_no],0,scomp,ratio);
805 });
806 } else {
807 ParallelFor(crse_S_fine, IntVect(0), ncomp,
808 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
809 {
810 amrex_avgdown_nodes(i,j,k,n,crsema[box_no],finema[box_no],0,scomp,ratio);
811 });
812 }
813 if (!Gpu::inNoSyncRegion()) {
815 }
816 } else
817#endif
818 {
819#ifdef AMREX_USE_OMP
820#pragma omp parallel if (Gpu::notInLaunchRegion())
821#endif
822 for (MFIter mfi(crse_S_fine,TilingIfNotGPU()); mfi.isValid(); ++mfi)
823 {
824 // NOTE: The tilebox is defined at the coarse level.
825 const Box& bx = mfi.tilebox();
826 Array4<value_type> const& crsearr = crse_S_fine.array(mfi);
827 Array4<value_type const> const& finearr = S_fine.const_array(mfi);
828
829 // NOTE: We copy from component scomp of the fine fab into component 0 of the crse fab
830 // because the crse fab is a temporary which was made starting at comp 0, it is
831 // not part of the actual crse multifab which came in.
832
833 if (is_cell_centered) {
834 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
835 {
836 amrex_avgdown(i,j,k,n,crsearr,finearr,0,scomp,ratio);
837 });
838 } else {
839 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
840 {
841 amrex_avgdown_nodes(i,j,k,n,crsearr,finearr,0,scomp,ratio);
842 });
843 }
844 }
845 }
846
847 S_crse.ParallelCopy(crse_S_fine,0,scomp,ncomp);
848 }
849}
850
851
852
853
854
862template <typename F>
863Real
864NormHelper (const MultiFab& x, int xcomp,
865 const MultiFab& y, int ycomp,
866 F const& f,
867 int numcomp, IntVect nghost, bool local)
868{
869 BL_ASSERT(x.boxArray() == y.boxArray());
870 BL_ASSERT(x.DistributionMap() == y.DistributionMap());
871 BL_ASSERT(x.nGrowVect().allGE(nghost) && y.nGrowVect().allGE(nghost));
872
873 Real sm = Real(0.0);
874#ifdef AMREX_USE_GPU
875 if (Gpu::inLaunchRegion()) {
876 auto const& xma = x.const_arrays();
877 auto const& yma = y.const_arrays();
879 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple<Real>
880 {
881 Real t = Real(0.0);
882 auto const& xfab = xma[box_no];
883 auto const& yfab = yma[box_no];
884 for (int n = 0; n < numcomp; ++n) {
885 t += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
886 }
887 return t;
888 });
889 } else
890#endif
891 {
892#ifdef AMREX_USE_OMP
893#pragma omp parallel if (!system::regtest_reduction) reduction(+:sm)
894#endif
895 for (MFIter mfi(x,true); mfi.isValid(); ++mfi)
896 {
897 Box const& bx = mfi.growntilebox(nghost);
898 Array4<Real const> const& xfab = x.const_array(mfi);
899 Array4<Real const> const& yfab = y.const_array(mfi);
900 AMREX_LOOP_4D(bx, numcomp, i, j, k, n,
901 {
902 sm += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
903 });
904 }
905 }
906
907 if (!local) {
909 }
910
911 return sm;
912}
913
922template <typename MMF, typename Pred, typename F>
923Real
924NormHelper (const MMF& mask,
925 const MultiFab& x, int xcomp,
926 const MultiFab& y, int ycomp,
927 Pred const& pf,
928 F const& f,
929 int numcomp, IntVect nghost, bool local)
930{
931 BL_ASSERT(x.boxArray() == y.boxArray());
932 BL_ASSERT(x.boxArray() == mask.boxArray());
933 BL_ASSERT(x.DistributionMap() == y.DistributionMap());
934 BL_ASSERT(x.DistributionMap() == mask.DistributionMap());
935 BL_ASSERT(x.nGrowVect().allGE(nghost) && y.nGrowVect().allGE(nghost));
936 BL_ASSERT(mask.nGrowVect().allGE(nghost));
937
938 Real sm = Real(0.0);
939#ifdef AMREX_USE_GPU
940 if (Gpu::inLaunchRegion()) {
941 auto const& xma = x.const_arrays();
942 auto const& yma = y.const_arrays();
943 auto const& mma = mask.const_arrays();
945 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple<Real>
946 {
947 Real t = Real(0.0);
948 if (pf(mma[box_no](i,j,k))) {
949 auto const& xfab = xma[box_no];
950 auto const& yfab = yma[box_no];
951 for (int n = 0; n < numcomp; ++n) {
952 t += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
953 }
954 }
955 return t;
956 });
957 } else
958#endif
959 {
960#ifdef AMREX_USE_OMP
961#pragma omp parallel if (!system::regtest_reduction) reduction(+:sm)
962#endif
963 for (MFIter mfi(x,true); mfi.isValid(); ++mfi)
964 {
965 Box const& bx = mfi.growntilebox(nghost);
966 Array4<Real const> const& xfab = x.const_array(mfi);
967 Array4<Real const> const& yfab = y.const_array(mfi);
968 auto const& mfab = mask.const_array(mfi);
969 AMREX_LOOP_4D(bx, numcomp, i, j, k, n,
970 {
971 if (pf(mfab(i,j,k))) {
972 sm += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
973 }
974 });
975 }
976 }
977
978 if (!local) {
980 }
981
982 return sm;
983}
984
985template <FabArrayType CMF, FabArrayType FMF>
986void average_face_to_cellcenter (CMF& cc, int dcomp,
988 int ngrow)
989{
990 IntVect ng_vect(ngrow);
991 average_face_to_cellcenter(cc, dcomp, fc, ng_vect);
992}
993
994template <FabArrayType CMF, FabArrayType FMF>
995void average_face_to_cellcenter (CMF& cc, int dcomp,
997 IntVect const& ng_vect)
998{
999 AMREX_ASSERT(cc.nComp() >= dcomp + AMREX_SPACEDIM);
1000 AMREX_ASSERT(fc[0]->nComp() == 1);
1001
1002#ifdef AMREX_USE_GPU
1003 if (Gpu::inLaunchRegion() && cc.isFusingCandidate()) {
1004 auto const& ccma = cc.arrays();
1005 AMREX_D_TERM(auto const& fxma = fc[0]->const_arrays();,
1006 auto const& fyma = fc[1]->const_arrays();,
1007 auto const& fzma = fc[2]->const_arrays(););
1008 ParallelFor(cc, ng_vect,
1009 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
1010 {
1011#if (AMREX_SPACEDIM == 1)
1012 GeometryData gd{};
1013 gd.coord = 0;
1014#endif
1015 amrex_avg_fc_to_cc(i,j,k, ccma[box_no], AMREX_D_DECL(fxma[box_no],
1016 fyma[box_no],
1017 fzma[box_no]),
1018 dcomp
1019#if (AMREX_SPACEDIM == 1)
1020 , gd
1021#endif
1022 );
1023 });
1024 if (!Gpu::inNoSyncRegion()) {
1026 }
1027 } else
1028#endif
1029 {
1030#ifdef AMREX_USE_OMP
1031#pragma omp parallel if (Gpu::notInLaunchRegion())
1032#endif
1033 for (MFIter mfi(cc,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1034 {
1035 const Box bx = mfi.growntilebox(ng_vect);
1036 auto const& ccarr = cc.array(mfi);
1037 AMREX_D_TERM(auto const& fxarr = fc[0]->const_array(mfi);,
1038 auto const& fyarr = fc[1]->const_array(mfi);,
1039 auto const& fzarr = fc[2]->const_array(mfi););
1040
1041#if (AMREX_SPACEDIM == 1)
1043 {
1044 GeometryData gd;
1045 gd.coord = 0;
1046 amrex_avg_fc_to_cc(i,j,k, ccarr, fxarr, dcomp, gd);
1047 });
1048#else
1050 {
1051 amrex_avg_fc_to_cc(i,j,k, ccarr, AMREX_D_DECL(fxarr,fyarr,fzarr), dcomp);
1052 });
1053#endif
1054 }
1055 }
1056}
1057
1058template <FabArrayType MF>
1060 const Vector<MF*>& crse,
1061 const IntVect& ratio, int ngcrse)
1062{
1063 AMREX_ASSERT(fine.size() == AMREX_SPACEDIM && crse.size() == AMREX_SPACEDIM);
1065 {{AMREX_D_DECL(fine[0],fine[1],fine[2])}},
1067 {{AMREX_D_DECL(crse[0],crse[1],crse[2])}},
1068 ratio, ngcrse);
1069}
1070
1071template <FabArrayType MF>
1073 const Vector<MF*>& crse, int ratio, int ngcrse)
1074{
1075 average_down_faces(fine,crse,IntVect{ratio},ngcrse);
1076}
1077
1078template <FabArrayType MF>
1081 int ratio, int ngcrse)
1082{
1083 average_down_faces(fine,crse,IntVect{ratio},ngcrse);
1084}
1085
1086template <FabArrayType MF>
1089 const IntVect& ratio, int ngcrse)
1090{
1091 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1092 average_down_faces(*fine[idim], *crse[idim], ratio, ngcrse);
1093 }
1094}
1095
1096template <BaseFabType FAB>
1098 const IntVect& ratio, int ngcrse)
1099{
1100 BL_PROFILE("average_down_faces");
1101
1102 AMREX_ASSERT(crse.nComp() == fine.nComp());
1103 AMREX_ASSERT(fine.ixType() == crse.ixType());
1104 const auto type = fine.ixType();
1105 int dir;
1106 for (dir = 0; dir < AMREX_SPACEDIM; ++dir) {
1107 if (type.nodeCentered(dir)) { break; }
1108 }
1109 auto tmptype = type;
1110 tmptype.unset(dir);
1111 if (dir >= AMREX_SPACEDIM || !tmptype.cellCentered()) {
1112 amrex::Abort("average_down_faces: not face index type");
1113 }
1114 const int ncomp = crse.nComp();
1115 if (isMFIterSafe(fine, crse))
1116 {
1117#ifdef AMREX_USE_GPU
1118 if (Gpu::inLaunchRegion() && crse.isFusingCandidate()) {
1119 auto const& crsema = crse.arrays();
1120 auto const& finema = fine.const_arrays();
1121 ParallelFor(crse, IntVect(ngcrse), ncomp,
1122 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
1123 {
1124 amrex_avgdown_faces(i,j,k,n, crsema[box_no], finema[box_no], 0, 0, ratio, dir);
1125 });
1126 if (!Gpu::inNoSyncRegion()) {
1128 }
1129 } else
1130#endif
1131 {
1132#ifdef AMREX_USE_OMP
1133#pragma omp parallel if (Gpu::notInLaunchRegion())
1134#endif
1135 for (MFIter mfi(crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1136 {
1137 const Box& bx = mfi.growntilebox(ngcrse);
1138 auto const& crsearr = crse.array(mfi);
1139 auto const& finearr = fine.const_array(mfi);
1140 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
1141 {
1142 amrex_avgdown_faces(i,j,k,n, crsearr, finearr, 0, 0, ratio, dir);
1143 });
1144 }
1145 }
1146 }
1147 else
1148 {
1149 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
1150 ncomp, ngcrse, MFInfo().SetArena(The_Async_Arena()), DefaultFabFactory<FAB>());
1151 average_down_faces(fine, ctmp, ratio, ngcrse);
1152 crse.ParallelCopy(ctmp,0,0,ncomp,ngcrse,ngcrse);
1153 }
1154}
1155
1156template <FabArrayType MF>
1159 const IntVect& ratio, const Geometry& crse_geom)
1160{
1161 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1162 average_down_faces(*fine[idim], *crse[idim], ratio, crse_geom);
1163 }
1164}
1165
1166template <BaseFabType FAB>
1168 const IntVect& ratio, const Geometry& crse_geom)
1169{
1170 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
1171 crse.nComp(), 0, MFInfo().SetArena(The_Async_Arena()));
1172 average_down_faces(fine, ctmp, ratio, 0);
1173 crse.ParallelCopy(ctmp,0,0,crse.nComp(),0,0,crse_geom.periodicity());
1174}
1175
1176template <FabArrayType MF>
1178{
1179 using T = typename MF::value_type;
1180 const int ncomp = mf.nComp();
1181 Gpu::DeviceVector<T> dv(ncomp);
1182 auto* dp = dv.data();
1183 bool found = false;
1184 auto loc = cell.dim3();
1185 for (MFIter mfi(mf); mfi.isValid() && !found; ++mfi)
1186 {
1187 Box const& box = mfi.validbox();
1188 if (box.contains(cell)) {
1189 found = true;
1190 auto const& fab = mf.const_array(mfi);
1191 amrex::ParallelFor(1, [=] AMREX_GPU_DEVICE (int) noexcept
1192 {
1193 for (int n = 0; n < ncomp; ++n) {
1194 dp[n] = fab(loc.x,loc.y,loc.z,n);
1195 }
1196 });
1197 }
1198 }
1199 Vector<T> hv;
1200 if (found) {
1201 hv.resize(ncomp);
1202 Gpu::copy(Gpu::deviceToHost, dv.begin(), dv.end(), hv.begin());
1203 }
1204 return hv;
1205}
1206
1207template <FabArrayType MF>
1208MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box const& bnd_bx)
1209{
1210 bool do_bnd = (!bnd_bx.isEmpty());
1211
1212 BoxArray const& ba = mf.boxArray();
1213 DistributionMapping const& dm = mf.DistributionMap();
1214 const auto nboxes = static_cast<int>(ba.size());
1215
1216 BoxList bl(ba.ixType());
1217 Vector<int> procmap;
1218 Vector<int> index_map;
1219 if (!do_bnd) {
1220 for (int i = 0; i < nboxes; ++i) {
1221 Box const& b = ba[i];
1222 IntVect lo = cell;
1223 lo[dir] = b.smallEnd(dir);
1224 if (b.contains(lo)) {
1225 IntVect hi = lo;
1226 hi[dir] = b.bigEnd(dir);
1227 Box b1d(lo,hi,b.ixType());
1228 bl.push_back(b1d);
1229 procmap.push_back(dm[i]);
1230 index_map.push_back(i);
1231 }
1232 }
1233 } else {
1234 for (int i = 0; i < nboxes; ++i) {
1235 Box const& b = ba[i];
1236 Box const& b1d = bnd_bx & b;
1237 if (b1d.ok()) {
1238 bl.push_back(b1d);
1239 procmap.push_back(dm[i]);
1240 index_map.push_back(i);
1241 }
1242 }
1243 }
1244
1245 if (bl.isEmpty()) {
1246 return MF();
1247 } else {
1248 BoxArray rba(std::move(bl));
1249 DistributionMapping rdm(std::move(procmap));
1250 MF rmf(rba, rdm, mf.nComp(), IntVect(0),
1251 MFInfo().SetArena(mf.arena()));
1252#ifdef AMREX_USE_OMP
1253#pragma omp parallel if (Gpu::notInLaunchRegion())
1254#endif
1255 for (MFIter mfi(rmf); mfi.isValid(); ++mfi) {
1256 Box const& b = mfi.validbox();
1257 auto const& dfab = rmf.array(mfi);
1258 auto const& sfab = mf.const_array(index_map[mfi.index()]);
1259 amrex::ParallelFor(b, mf.nComp(),
1260 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
1261 {
1262 dfab(i,j,k,n) = sfab(i,j,k,n);
1263 });
1264 }
1265 return rmf;
1266 }
1267}
1268
1270namespace detail {
1271template <typename Op, typename T, typename F>
1272void reduce_to_plane (Array4<T> const& ar, int direction, Box const& bx, int box_no,
1273 F const& f)
1274{
1275#if defined(AMREX_USE_GPU)
1276 Box b2d = bx;
1277 b2d.setRange(direction,0);
1278 const auto blo = amrex::lbound(bx);
1279 const auto len = amrex::length(bx);
1280 constexpr int nthreads = 128;
1281 auto nblocks = static_cast<int>(b2d.numPts());
1282#ifdef AMREX_USE_SYCL
1283 constexpr std::size_t shared_mem_bytes = sizeof(T)*Gpu::Device::warp_size;
1284 amrex::launch<nthreads>(nblocks, shared_mem_bytes, Gpu::gpuStream(),
1285 [=] AMREX_GPU_DEVICE (Gpu::Handler const& h)
1286 {
1287 int bid = h.blockIdx();
1288 int tid = h.threadIdx();
1289#else
1290 amrex::launch<nthreads>(nblocks, Gpu::gpuStream(),
1291 [=] AMREX_GPU_DEVICE ()
1292 {
1293 int bid = blockIdx.x;
1294 int tid = threadIdx.x;
1295#endif
1296 T tmp;
1297 Op().init(tmp);
1298 T* p;
1299 if (direction == 0) {
1300 int k = bid / len.y;
1301 int j = bid - k*len.y;
1302 k += blo.z;
1303 j += blo.y;
1304 for (int i = blo.x + tid; i < blo.x+len.x; i += nthreads) {
1305 Op().local_update(tmp, f(box_no,i,j,k));
1306 }
1307 p = ar.ptr(0,j,k);
1308 } else if (direction == 1) {
1309 int k = bid / len.x;
1310 int i = bid - k*len.x;
1311 k += blo.z;
1312 i += blo.x;
1313 for (int j = blo.y + tid; j < blo.y+len.y; j += nthreads) {
1314 Op().local_update(tmp, f(box_no,i,j,k));
1315 }
1316 p = ar.ptr(i,0,k);
1317 } else {
1318 int j = bid / len.x;
1319 int i = bid - j*len.x;
1320 j += blo.y;
1321 i += blo.x;
1322 for (int k = blo.z + tid; k < blo.z+len.z; k += nthreads) {
1323 Op().local_update(tmp, f(box_no,i,j,k));
1324 }
1325 p = ar.ptr(i,j,0);
1326 }
1327#ifdef AMREX_USE_SYCL
1328 Op().template parallel_update<T>(*p, tmp, h);
1329#else
1330 Op().template parallel_update<T,nthreads>(*p, tmp);
1331#endif
1332 });
1333#else
1334 // CPU
1335 if (direction == 0) {
1336 AMREX_LOOP_3D(bx, i, j, k,
1337 {
1338 Op().local_update(ar(0,j,k), f(box_no,i,j,k));
1339 });
1340 } else if (direction == 1) {
1341 AMREX_LOOP_3D(bx, i, j, k,
1342 {
1343 Op().local_update(ar(i,0,k), f(box_no,i,j,k));
1344 });
1345 } else {
1346 AMREX_LOOP_3D(bx, i, j, k,
1347 {
1348 Op().local_update(ar(i,j,0), f(box_no,i,j,k));
1349 });
1350 }
1351#endif
1352}
1353}
1355
1356template <typename Op, typename T, BaseFabType FAB, typename F>
1357#ifndef AMREX_USE_CUDA
1358requires (IsCallableR<T,F,int,int,int,int>::value)
1359#endif
1360BaseFab<T>
1361ReduceToPlane (int direction, Box const& a_domain, FabArray<FAB> const& mf, F const& f)
1362{
1363 Box const domain = amrex::convert(a_domain, mf.ixType());
1364
1365 Box domain2d = domain;
1366 domain2d.setRange(direction, 0);
1367
1368 T initval;
1369 Op().init(initval);
1370
1371 BaseFab<T> r(domain2d);
1372 r.template setVal<RunOn::Device>(initval);
1373 auto const& ar = r.array();
1374
1375 for (MFIter mfi(mf,MFItInfo().UseDefaultStream().DisableDeviceSync());
1376 mfi.isValid(); ++mfi)
1377 {
1378 Box bx = mfi.validbox() & domain;
1379 if (bx.ok()) {
1380 int box_no = mfi.LocalIndex();
1381 detail::reduce_to_plane<Op, T>(ar, direction, bx, box_no, f);
1382 }
1383 }
1384 if (!Gpu::inNoSyncRegion()) {
1386 }
1387
1388 return r;
1389}
1390
1392namespace detail {
1393template <typename Op, MultiFabLike FA, typename F>
1394FA reduce_to_plane (int direction, Box const& domain, FA const& mf, F const& f,
1395 IntVect const& nghost = IntVect(0))
1396{
1397 using T = typename FA::value_type;
1398
1399 Box const ndomain = amrex::convert(domain, mf.ixType());
1400
1401 auto npts = amrex::convert(mf.boxArray(),IntVect(0)).numPts();
1402 if (npts != amrex::convert(domain, amrex::IntVect(0)).numPts()) {
1403 amrex::Abort("ReduceToPlaneMF: mf's BoxArray must have a rectangular domain.");
1404 }
1405
1406 // the reduction reads mf's guard cells, so it cannot include more guard
1407 // cells than mf actually has
1408 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nghost.allLE(mf.nGrowVect()),
1409 "ReduceToPlaneMF2: nghost must not exceed mf.nGrowVect().");
1410
1411 T initval;
1412 Op().init(initval);
1413
1414 // guard cells in the reduction direction are summed into the plane; guard
1415 // cells in the transverse directions are stored in the result's ghost cells
1416 // (and folded into the valid region by the caller's ParallelAdd)
1417 IntVect ng2d = nghost;
1418 ng2d[direction] = 0;
1419 Box const gdomain = amrex::grow(ndomain, nghost);
1420
1421 BoxList bl = mf.boxArray().boxList();
1422 for (auto& b : bl) {
1423 b.setRange(direction, 0);
1424 }
1425 BoxArray ba(std::move(bl));
1426 FA tmpfa(ba, mf.DistributionMap(), 1, ng2d);
1427 tmpfa.setVal(initval);
1428
1429 for (MFIter mfi(mf); mfi.isValid(); ++mfi)
1430 {
1431 Box bx = amrex::grow(mfi.validbox(), nghost) & gdomain;
1432 if (bx.ok()) {
1433 int box_no = mfi.LocalIndex();
1434 detail::reduce_to_plane<Op, T>(tmpfa.array(mfi), direction, bx, box_no, f);
1435 }
1436 }
1437
1438 return tmpfa;
1439}
1440
1441}
1443
1444template <typename Op, MultiFabLike FA, typename F>
1445#ifndef AMREX_USE_CUDA
1446requires (IsCallableR<typename FA::value_type,
1447 F,int,int,int,int>::value)
1448#endif
1449FA ReduceToPlaneMF (int direction, Box const& domain, FA const& mf, F const& f)
1450{
1451 auto [fa3, fa2] = ReduceToPlaneMF2<Op>(direction, domain, mf, f);
1452 fa3.ParallelCopy(fa2);
1453 return std::move(fa3);
1454}
1455
1456template <typename Op, MultiFabLike FA, typename F>
1457#ifndef AMREX_USE_CUDA
1458requires (IsCallableR<typename FA::value_type,
1459 F,int,int,int,int>::value)
1460#endif
1461std::pair<FA,FA>
1462ReduceToPlaneMF2 (int direction, Box const& domain, FA const& mf, F const& f,
1463 IntVect const& nghost, Periodicity const& period)
1464{
1465 using T = typename FA::value_type;
1466
1467 T initval;
1468 Op().init(initval);
1469
1470 auto tmpmf = detail::reduce_to_plane<Op>(direction, domain, mf, f, nghost);
1471
1472 BoxList bl2d(mf.ixType());
1473 Vector<int> procmap2d;
1474 auto const& ba3d = mf.boxArray();
1475 auto const& dm3d = mf.DistributionMap();
1476 int dlo = domain.smallEnd(direction);
1477 for (int i = 0, N = mf.size(); i < N; ++i) {
1478 Box b = ba3d[i];
1479 if (b.smallEnd(direction) <= dlo && dlo <= b.bigEnd(direction)) {
1480 b.setRange(direction, 0);
1481 bl2d.push_back(b);
1482 procmap2d.push_back(dm3d[i]);
1483 }
1484 }
1485
1486 BoxArray ba2d(std::move(bl2d));
1487 DistributionMapping dm2d(std::move(procmap2d));
1488
1489 FA mf2d(ba2d, dm2d, 1, 0);
1490 mf2d.setVal(initval);
1491
1492 static_assert(std::is_same_v<Op, ReduceOpSum>, "Currently only ReduceOpSum is supported.");
1493 // fold the transverse guard cells of the local results into the valid
1494 // region of the unique result (sums the quantity across box seams)
1495 IntVect ng2d = nghost;
1496 ng2d[direction] = 0;
1497 mf2d.ParallelAdd(tmpmf, 0, 0, 1, ng2d, IntVect(0), IntVect(0), period);
1498
1499 return std::make_pair(std::move(tmpmf), std::move(mf2d));
1500}
1501
1502template <typename Op, MultiFabLike FA, typename F>
1503#ifndef AMREX_USE_CUDA
1504requires (IsCallableR<typename FA::value_type,
1505 F,int,int,int,int>::value)
1506#endif
1507std::pair<FA,FA>
1508ReduceToPlaneMF2Patchy (int direction, Box const& domain, FA const& mf,
1509 IntVect const& plane_max_grid_size, F const& f)
1510{
1511 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mf.is_cell_centered(),
1512 "ReduceToPlaneMF2Patchy: only cell-centered data is supported.");
1513
1514 Box const ndomain = amrex::convert(domain, mf.ixType());
1515 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ndomain.contains(mf.boxArray().minimalBox()),
1516 "ReduceToPlaneMF2Patchy: subdomains are not supported; domain must cover the full extent of mf.boxArray().");
1517
1518 using T = typename FA::value_type;
1519 T initval;
1520 Op().init(initval);
1521
1522 BoxList bl_patch2d = mf.boxArray().boxList();
1523 for (auto& b : bl_patch2d) {
1524 b.setRange(direction, 0);
1525 }
1526 BoxArray ba_patch2d(std::move(bl_patch2d));
1527 FA tmpmf(ba_patch2d, mf.DistributionMap(), 1, 0);
1528 tmpmf.setVal(initval);
1529
1530 for (MFIter mfi(mf); mfi.isValid(); ++mfi)
1531 {
1532 Box bx = mfi.validbox() & ndomain;
1533 if (bx.ok()) {
1534 int box_no = mfi.LocalIndex();
1535 detail::reduce_to_plane<Op, T>(tmpmf.array(mfi), direction, bx, box_no, f);
1536 }
1537 }
1538
1539 BoxList bl_unique = ba_patch2d.boxList();
1540 bl_unique.simplify(); // reduces work for .removeOverlap()
1541 BoxArray ba2d(std::move(bl_unique));
1542 ba2d.removeOverlap();
1543
1544 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(plane_max_grid_size.allGE(1),
1545 "ReduceToPlaneMF2Patchy: plane_max_grid_size must be >= 1 in every direction.");
1546 IntVect mgs2d = plane_max_grid_size;
1547 mgs2d[direction] = 1;
1548 ba2d.maxSize(mgs2d);
1549
1550 DistributionMapping dm2d(ba2d);
1551
1552 FA mf2d(ba2d, dm2d, 1, 0);
1553 mf2d.setVal(initval);
1554
1555 static_assert(std::is_same_v<Op, ReduceOpSum>, "Currently only ReduceOpSum is supported.");
1556 mf2d.ParallelAdd(tmpmf);
1557
1558 return std::make_pair(std::move(tmpmf), std::move(mf2d));
1559}
1560
1561}
1562
1563#endif
Fixed-size array types for use on GPU and CPU.
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:562
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:49
#define BL_ASSERT(EX)
Definition AMReX_BLassert.H:39
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
#define AMREX_HOST_DEVICE_PARALLEL_FOR_1D(...)
Definition AMReX_GpuLaunchMacrosC.nolint.H:109
#define AMREX_HOST_DEVICE_PARALLEL_FOR_3D(...)
Definition AMReX_GpuLaunchMacrosC.nolint.H:110
#define AMREX_HOST_DEVICE_PARALLEL_FOR_4D(...)
Definition AMReX_GpuLaunchMacrosC.nolint.H:111
#define AMREX_LAUNCH_HOST_DEVICE_LAMBDA(...)
Definition AMReX_GpuLaunch.nolint.H:16
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1131
Real * pdst
Definition AMReX_HypreMLABecLap.cpp:1132
Array4< Real > fine
Definition AMReX_InterpFaceRegister.cpp:90
Array4< int const > mask
Definition AMReX_InterpFaceRegister.cpp:93
Array4< Real const > crse
Definition AMReX_InterpFaceRegister.cpp:92
#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
GpuArray< int, 3 > dlo
Definition AMReX_MLEBNodeFDLaplacian.cpp:1100
#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 FortranArrayBox(FAB)-like object.
Definition AMReX_BaseFab.H:222
Array4< T const > array() const noexcept
Create an Array4 view over all components.
Definition AMReX_BaseFab.H:475
Reference-counted collection of Boxes.
Definition AMReX_BoxArray.H:681
IndexType ixType() const noexcept
Return index type of this BoxArray.
Definition AMReX_BoxArray.H:1268
BoxList boxList() const
Create a BoxList from this BoxArray.
Definition AMReX_BoxArray.cpp:965
void removeOverlap(bool simplify=true)
Change the BoxArray to one with no overlap and then simplify it (see the simplify function in BoxList...
Definition AMReX_BoxArray.cpp:1474
BoxArray & coarsen(int refinement_ratio)
Coarsen each Box in the BoxArray by refinement_ratio.
Definition AMReX_BoxArray.cpp:685
Long size() const noexcept
Return the number of boxes in the BoxArray.
Definition AMReX_BoxArray.H:758
BoxArray & maxSize(int block_size)
Split any Box whose side length exceeds block_size.
Definition AMReX_BoxArray.cpp:553
A list of Boxes sharing a common IndexType.
Definition AMReX_BoxList.H:109
bool isEmpty() const noexcept
Is this BoxList empty?
Definition AMReX_BoxList.H:202
int simplify(bool best=false)
Merge adjacent Boxes in this BoxList. Return the number of Boxes merged. If "best" is specified we do...
Definition AMReX_BoxList.cpp:654
void push_back(const Box &bn)
Append a Box to this BoxList.
Definition AMReX_BoxList.H:151
__host__ __device__ const IntVectND< dim > & bigEnd() const &noexcept
Return the inclusive upper bound of the box.
Definition AMReX_Box.H:136
__host__ __device__ bool isEmpty() const noexcept
Checks if it is an empty BoxND.
Definition AMReX_Box.H:223
__host__ __device__ bool contains(const IntVectND< dim > &p) const noexcept
Return true if argument is contained within BoxND.
Definition AMReX_Box.H:233
__host__ __device__ IndexTypeND< dim > ixType() const noexcept
Return the indexing type.
Definition AMReX_Box.H:148
__host__ __device__ BoxND & setRange(int dir, int sm_index, int n_cells=1) noexcept
Set the entire range in a given direction, starting at sm_index with length n_cells....
Definition AMReX_Box.H:1143
__host__ __device__ bool ok() const noexcept
Return true if high bounds are >= low bounds and the index type is valid.
Definition AMReX_Box.H:229
__host__ __device__ const IntVectND< dim > & smallEnd() const &noexcept
Return the inclusive lower bound of the box.
Definition AMReX_Box.H:124
Default FabFactory that calls new / delete on the FAB type directly.
Definition AMReX_FabFactory.H:124
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:51
IntVect nGrowVect() const noexcept
Definition AMReX_FabArrayBase.H:85
bool isFusingCandidate() const noexcept
Is this a good candidate for kernel fusing?
Definition AMReX_FabArrayBase.cpp:2713
bool is_cell_centered() const noexcept
This tests on whether the FabArray is cell-centered.
Definition AMReX_FabArrayBase.cpp:2707
bool is_nodal() const noexcept
This tests on whether the FabArray is fully nodal.
Definition AMReX_FabArrayBase.cpp:2695
IndexType ixType() const noexcept
Return index type.
Definition AMReX_FabArrayBase.H:91
const DistributionMapping & DistributionMap() const noexcept
Return constant reference to associated DistributionMapping.
Definition AMReX_FabArrayBase.H:135
int nComp() const noexcept
Return number of variables (aka components) associated with each point.
Definition AMReX_FabArrayBase.H:88
const BoxArray & boxArray() const noexcept
Return a constant reference to the BoxArray that defines the valid region associated with this FabArr...
Definition AMReX_FabArrayBase.H:100
An Array of FortranArrayBox(FAB)-like Objects.
Definition AMReX_FabArray.H:356
void ParallelCopy(const FabArray< FAB > &src, const Periodicity &period=Periodicity::NonPeriodic(), CpOp op=FabArrayBase::COPY)
Definition AMReX_FabArray.H:971
MultiArray4< typename FabArray< FAB >::value_type > arrays() noexcept
Build, if needed, and return mutable Array4 views for local FABs.
Definition AMReX_FabArray.H:713
MultiArray4< typename FabArray< FAB >::value_type const > const_arrays() const noexcept
Read-only convenience wrapper equivalent to arrays() const.
Definition AMReX_FabArray.H:731
Array4< typename FabArray< FAB >::value_type const > array(const MFIter &mfi) const noexcept
Read-only Array4 view for the FAB referenced by iterator mfi.
Definition AMReX_FabArray.H:621
Array4< typename FabArray< FAB >::value_type const > const_array(const MFIter &mfi) const noexcept
Synonym for array(const MFIter&) that highlights read-only semantics.
Definition AMReX_FabArray.H:649
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:85
Periodicity periodicity() const noexcept
Return the Periodicity based on the length of the domain.
Definition AMReX_Geometry.H:424
GPU-compatible tuple.
Definition AMReX_Tuple.H:104
static constexpr int warp_size
Definition AMReX_GpuDevice.H:236
__host__ __device__ constexpr Dim3 dim3() const noexcept
Definition AMReX_IntVect.H:262
__host__ static __device__ constexpr IntVectND< dim > TheUnitVector() noexcept
This static member function returns a reference to a constant IntVectND object, all of whose dim argu...
Definition AMReX_IntVect.H:781
__host__ __device__ constexpr bool allGE(const IntVectND< dim > &rhs) const noexcept
Returns true if this is greater than or equal to argument for all components. NOTE: This is NOT a str...
Definition AMReX_IntVect.H:542
__host__ static __device__ constexpr IntVectND< dim > TheZeroVector() noexcept
This static member function returns a reference to a constant IntVectND object, all of whose dim argu...
Definition AMReX_IntVect.H:771
a one-thingy-per-box distributed object
Definition AMReX_LayoutData.H:13
Iterator for looping ever tiles and boxes of amrex::FabArray based containers.
Definition AMReX_MFIter.H:88
bool isValid() const noexcept
Is the iterator valid i.e. is it associated with a FAB?
Definition AMReX_MFIter.H:176
A collection (stored as an array) of FArrayBox objects.
Definition AMReX_MultiFab.H:40
Dynamically allocated vector for trivially copyable data.
Definition AMReX_PODVector.H:308
iterator begin() noexcept
Definition AMReX_PODVector.H:680
iterator end() noexcept
Definition AMReX_PODVector.H:684
T * data() noexcept
Definition AMReX_PODVector.H:672
This provides length of period for periodic domains. 0 means it is not periodic in that direction....
Definition AMReX_Periodicity.H:17
static const Periodicity & NonPeriodic() noexcept
Definition AMReX_Periodicity.cpp:52
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:29
A Collection of IArrayBoxes.
Definition AMReX_iMultiFab.H:34
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:80
amrex_long Long
Definition AMReX_INT.H:30
__host__ __device__ Dim3 length(Array4< T > const &a) noexcept
Return the spatial extents of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1379
__host__ __device__ Dim3 lbound(Array4< T > const &a) noexcept
Return the inclusive lower bounds of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1351
__host__ __device__ BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Return a copy of b converted to the nodal flags typ.
Definition AMReX_Box.H:1630
__host__ __device__ BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Return a copy of b coarsened by the isotropic ratio ref_ratio.
Definition AMReX_Box.H:1469
__host__ __device__ BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Return a copy of b grown uniformly by i cells in every direction.
Definition AMReX_Box.H:1326
std::array< T, N > Array
Definition AMReX_Array.H:31
Arena * The_Async_Arena()
Definition AMReX_Arena.cpp:825
void Sum(Gpu::DeviceVector< T > &v, MPI_Comm comm)
Definition AMReX_GpuParallelReduce.H:37
void copy(HostToDevice, InIter begin, InIter end, OutIter result) noexcept
A host-to-device copy routine. Note this is just a wrapper around memcpy, so it assumes contiguous st...
Definition AMReX_GpuContainers.H:128
static constexpr DeviceToHost deviceToHost
Definition AMReX_GpuContainers.H:106
void streamSynchronize() noexcept
Definition AMReX_GpuDevice.H:310
bool inLaunchRegion() noexcept
Definition AMReX_GpuControl.H:88
bool inNoSyncRegion() noexcept
Definition AMReX_GpuControl.H:148
gpuStream_t gpuStream() noexcept
Definition AMReX_GpuDevice.H:291
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition AMReX_ParallelContext.H:70
Definition AMReX_Amr.cpp:50
std::pair< FA, FA > ReduceToPlaneMF2Patchy(int direction, Box const &domain, FA const &mf, IntVect const &plane_max_grid_size, F const &f)
Reduce FabArray/MultiFab data to plane FabArray for patchy BoxArrays.
Definition AMReX_MultiFabUtil.H:1508
void FillRandomNormal(MultiFab &mf, int scomp, int ncomp, Real mean, Real stddev)
Fill MultiFab with random numbers from normal distribution.
Definition AMReX_MultiFabUtil.cpp:1246
void FillRandom(MultiFab &mf, int scomp, int ncomp)
Fill MultiFab with random numbers from uniform distribution.
Definition AMReX_MultiFabUtil.cpp:1233
ReduceData< Ts... >::Type ParReduce(TypeList< Ops... > operation_list, TypeList< Ts... > type_list, FabArray< FAB > const &fa, IntVect const &nghost, F &&f)
Parallel reduce for MultiFab/FabArray. The reduce result is local and it's the user's responsibility ...
Definition AMReX_ParReduce.H:48
void average_down(const MultiFab &S_fine, MultiFab &S_crse, const Geometry &fgeom, const Geometry &cgeom, int scomp, int ncomp, int rr)
Definition AMReX_MultiFabUtil.cpp:359
__host__ __device__ void cast(BaseFab< Tto > &tofab, BaseFab< Tfrom > const &fromfab, Box const &bx, SrcComp scomp, DestComp dcomp, NumComps ncomp) noexcept
Cast components from one BaseFab to another over a region.
Definition AMReX_BaseFabUtility.H:30
std::unique_ptr< MultiFab > get_slice_data(int dir, Real coord, const MultiFab &cc, const Geometry &geom, int start_comp, int ncomp, bool interpolate, RealBox const &bnd_rbx)
Definition AMReX_MultiFabUtil.cpp:598
void average_face_to_cellcenter(MultiFab &cc, int dcomp, const Vector< const MultiFab * > &fc, IntVect const &ng_vect)
Definition AMReX_MultiFabUtil.cpp:156
iMultiFab makeFineMask(const BoxArray &cba, const DistributionMapping &cdm, const BoxArray &fba, const IntVect &ratio, int crse_value, int fine_value, MFInfo const &info)
Definition AMReX_MultiFabUtil.cpp:653
Vector< MultiFab > convexify(Vector< MultiFab const * > const &mf, Vector< IntVect > const &refinement_ratio)
Convexify AMR data.
Definition AMReX_MultiFabUtil.cpp:1259
MF get_line_data(MF const &mf, int dir, IntVect const &cell, Box const &bnd_bx=Box())
Get data in a line of MultiFab/FabArray.
Definition AMReX_MultiFabUtil.H:1208
FA ReduceToPlaneMF(int direction, Box const &domain, FA const &mf, F const &f)
Reduce FabArray/MultiFab data to plane FabArray.
Definition AMReX_MultiFabUtil.H:1449
MultiFab periodicShift(MultiFab const &mf, IntVect const &offset, Periodicity const &period)
Periodic shift MultiFab.
Definition AMReX_MultiFabUtil.cpp:817
BaseFab< T > ReduceToPlane(int direction, Box const &domain, FabArray< FAB > const &mf, F const &f)
Reduce FabArray/MultiFab data to a plane Fab.
Definition AMReX_MultiFabUtil.H:1361
void ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
Definition AMReX_CTOParallelForImpl.H:202
Real volumeWeightedSum(Vector< MultiFab const * > const &mf, int icomp, Vector< Geometry > const &geom, Vector< IntVect > const &ratio, bool local)
Volume weighted sum for a vector of MultiFabs.
Definition AMReX_MultiFabUtil.cpp:988
BoxND< 3 > Box
Box is an alias for amrex::BoxND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:35
void average_down_faces(const Vector< const MF * > &fine, const Vector< MF * > &crse, const IntVect &ratio, int ngcrse=0)
Average fine face-based FabArray onto crse face-based FabArray.
Definition AMReX_MultiFabUtil.H:1059
bool isMFIterSafe(const FabArrayBase &x, const FabArrayBase &y)
Definition AMReX_MFIter.H:256
MultiFab ToMultiFab(const iMultiFab &imf)
Convert iMultiFab to MultiFab.
Definition AMReX_MultiFabUtil.cpp:588
Vector< typename MF::value_type > get_cell_data(MF const &mf, IntVect const &cell)
Get data in a cell of MultiFab/FabArray.
Definition AMReX_MultiFabUtil.H:1177
int nComp(FabArrayBase const &fa)
Convenience wrapper that forwards to fa.nComp().
Definition AMReX_FabArrayBase.cpp:2860
void writeFabs(const MultiFab &mf, const std::string &name)
Write each fab individually.
Definition AMReX_MultiFabUtil.cpp:575
Real NormHelper(const MultiFab &x, int xcomp, const MultiFab &y, int ycomp, F const &f, int numcomp, IntVect nghost, bool local)
Returns part of a norm based on two MultiFabs.
Definition AMReX_MultiFabUtil.H:864
void computeDivergence(MultiFab &divu, const Array< MultiFab const *, 3 > &umac, const Geometry &geom)
Computes divergence of face-data stored in the umac MultiFab.
Definition AMReX_MultiFabUtil.cpp:741
void average_down_edges(const Vector< const MultiFab * > &fine, const Vector< MultiFab * > &crse, const IntVect &ratio, int ngcrse)
Average fine edge-based MultiFab onto crse edge-based MultiFab.
Definition AMReX_MultiFabUtil.cpp:493
Gpu::HostVector< Real > sumToLine(MultiFab const &mf, int icomp, int ncomp, Box const &domain, int direction, bool local)
Sum MultiFab data to line.
Definition AMReX_MultiFabUtil.cpp:841
IntVectND< 3 > IntVect
IntVect is an alias for amrex::IntVectND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:38
std::pair< FA, FA > ReduceToPlaneMF2(int direction, Box const &domain, FA const &mf, F const &f, IntVect const &nghost=IntVect(0), Periodicity const &period=Periodicity::NonPeriodic())
Reduce FabArray/MultiFab data to plane FabArray.
Definition AMReX_MultiFabUtil.H:1462
void print_state(const MultiFab &mf, const IntVect &cell, const int n, const IntVect &ng)
Output state data for a single zone.
Definition AMReX_MultiFabUtil.cpp:570
bool TilingIfNotGPU() noexcept
Definition AMReX_MFIter.H:12
void sum_fine_to_coarse(const MultiFab &S_fine, MultiFab &S_crse, int scomp, int ncomp, const IntVect &ratio, const Geometry &cgeom, const Geometry &)
Definition AMReX_MultiFabUtil.cpp:438
void Abort(const std::string &msg)
Print a fatal-error message to stderr and abort execution.
Definition AMReX.cpp:242
const int[]
Definition AMReX_BLProfiler.cpp:1665
void average_node_to_cellcenter(MultiFab &cc, int dcomp, const MultiFab &nd, int scomp, int ncomp, int ngrow)
Average nodal-based MultiFab onto cell-centered MultiFab.
Definition AMReX_MultiFabUtil.cpp:63
void FourthOrderInterpFromFineToCoarse(MultiFab &cmf, int scomp, int ncomp, MultiFab const &fmf, IntVect const &ratio)
Fourth-order interpolation from fine to coarse level.
Definition AMReX_MultiFabUtil.cpp:1142
void average_cellcenter_to_face(const Vector< MultiFab * > &fc, const MultiFab &cc, const Geometry &geom, int ncomp, bool use_harmonic_averaging, int ngrow)
Average cell-centered MultiFab onto face-based MultiFab with geometric weighting.
Definition AMReX_MultiFabUtil.cpp:241
FabArray< BaseFab< Long > > ToLongMultiFab(const iMultiFab &imf)
Convert iMultiFab to Long.
Definition AMReX_MultiFabUtil.cpp:593
void computeGradient(MultiFab &grad, const Array< MultiFab const *, 3 > &umac, const Geometry &geom)
Computes gradient of face-data stored in the umac MultiFab.
Definition AMReX_MultiFabUtil.cpp:793
void average_edge_to_cellcenter(MultiFab &cc, int dcomp, const Vector< const MultiFab * > &edge, int ngrow)
Average edge-based MultiFab onto cell-centered MultiFab.
Definition AMReX_MultiFabUtil.cpp:105
void average_down_nodal(const FabArray< FAB > &S_fine, FabArray< FAB > &S_crse, const IntVect &ratio, int ngcrse=0, bool mfiter_is_definitely_safe=false)
Average fine nodal MultiFab onto coarse nodal MultiFab.
Definition AMReX_MultiFabUtil.H:677
A multidimensional array accessor.
Definition AMReX_Array4.H:289
parallel copy or add
Definition AMReX_FabArrayBase.H:630
GPU-safe snapshot of Geometry data suitable for capture in GPU kernels.
Definition AMReX_Geometry.H:36
int coord
Coordinate-system type: 0 for Cartesian, 1 for RZ, and 2 for spherical.
Definition AMReX_Geometry.H:70
FabArray memory allocation information.
Definition AMReX_FabArray.H:73
Definition AMReX_MFIter.H:20
Struct for holding types.
Definition AMReX_TypeList.H:13