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
540 Gpu::HostVector<Real> sumToLine (MultiFab const& mf, int icomp, int ncomp,
541 Box const& domain, int direction, bool local = false);
542
550 Real volumeWeightedSum (Vector<MultiFab const*> const& mf, int icomp,
551 Vector<Geometry> const& geom,
552 Vector<IntVect> const& ratio,
553 bool local = false);
554
568 void FourthOrderInterpFromFineToCoarse (MultiFab& cmf, int scomp, int ncomp,
569 MultiFab const& fmf,
570 IntVect const& ratio);
571
582 void FillRandom (MultiFab& mf, int scomp, int ncomp);
583
595 void FillRandomNormal (MultiFab& mf, int scomp, int ncomp, Real mean, Real stddev);
596
608 [[nodiscard]] Vector<MultiFab> convexify (Vector<MultiFab const*> const& mf,
609 Vector<IntVect> const& refinement_ratio);
610}
611
612namespace amrex {
613
614template <BaseFabType FAB>
615iMultiFab
616makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
617 int crse_value, int fine_value, MFInfo const& info)
618{
619 return makeFineMask(cmf.boxArray(), cmf.DistributionMap(), cmf.nGrowVect(),
620 fba, ratio, Periodicity::NonPeriodic(), crse_value, fine_value,
621 info);
622}
623
624template <BaseFabType FAB>
625iMultiFab
626makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
627 Periodicity const& period, int crse_value, int fine_value,
628 MFInfo const& info)
629{
630 return makeFineMask(cmf.boxArray(), cmf.DistributionMap(), cmf.nGrowVect(),
631 fba, ratio, period, crse_value, fine_value, info);
632}
633
634template <BaseFabType FAB>
635iMultiFab
637 const IntVect& cnghost, const IntVect& ratio,
638 Periodicity const& period, int crse_value, int fine_value,
639 MFInfo const& info)
640{
641 iMultiFab mask(cmf.boxArray(), cmf.DistributionMap(), 1, cnghost, info);
642 mask.setVal(crse_value);
643
644 iMultiFab foo(amrex::coarsen(fmf.boxArray(),ratio), fmf.DistributionMap(),
645 1, 0, MFInfo().SetAlloc(false));
646 const FabArrayBase::CPC& cpc = mask.getCPC(cnghost,foo,IntVect::TheZeroVector(),period);
647 mask.setVal(fine_value, cpc, 0, 1);
648
649 return mask;
650}
651
652template <BaseFabType FAB>
653iMultiFab
655 const IntVect& cnghost, const IntVect& ratio,
656 Periodicity const& period, int crse_value, int fine_value,
657 LayoutData<int>& has_cf, MFInfo const& info)
658{
659 iMultiFab mask(cmf.boxArray(), cmf.DistributionMap(), 1, cnghost, info);
660 mask.setVal(crse_value);
661
662 iMultiFab foo(amrex::coarsen(fmf.boxArray(),ratio), fmf.DistributionMap(),
663 1, 0, MFInfo().SetAlloc(false));
664 const FabArrayBase::CPC& cpc = mask.getCPC(cnghost,foo,IntVect::TheZeroVector(),period);
665 mask.setVal(fine_value, cpc, 0, 1);
666
667 has_cf = mask.RecvLayoutMask(cpc);
668
669 return mask;
670}
671
674template <BaseFabType FAB>
676 const IntVect& ratio, int ngcrse, bool mfiter_is_definitely_safe)
677{
678 AMREX_ASSERT(fine.is_nodal());
679 AMREX_ASSERT(crse.is_nodal());
680 AMREX_ASSERT(crse.nComp() == fine.nComp());
681
682 int ncomp = crse.nComp();
683 using value_type = typename FAB::value_type;
684
685 if (mfiter_is_definitely_safe || isMFIterSafe(fine, crse))
686 {
687#ifdef AMREX_USE_OMP
688#pragma omp parallel if (Gpu::notInLaunchRegion())
689#endif
690 for (MFIter mfi(crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
691 {
692 const Box& bx = mfi.growntilebox(ngcrse);
693 Array4<value_type> const& crsearr = crse.array(mfi);
694 Array4<value_type const> const& finearr = fine.const_array(mfi);
695
697 {
698 amrex_avgdown_nodes(tbx,crsearr,finearr,0,0,ncomp,ratio);
699 });
700 }
701 }
702 else
703 {
704 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
705 ncomp, ngcrse);
706 average_down_nodal(fine, ctmp, ratio, ngcrse);
707 crse.ParallelCopy(ctmp,0,0,ncomp,ngcrse,ngcrse);
708 }
709}
710
711// *************************************************************************************************************
712
713// Average fine cell-based MultiFab onto crse cell-centered MultiFab.
714// We do NOT assume that the coarse layout is a coarsened version of the fine layout.
715// This version does NOT use volume-weighting
716template <BaseFabType FAB>
717void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse, int scomp, int ncomp, int rr)
718{
719 average_down(S_fine,S_crse,scomp,ncomp,rr*IntVect::TheUnitVector());
720}
721
722template <BaseFabType FAB>
723void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse,
724 int scomp, int ncomp, const IntVect& ratio)
725{
726 BL_PROFILE("amrex::average_down");
727 AMREX_ASSERT(S_crse.nComp() == S_fine.nComp());
728 AMREX_ASSERT((S_crse.is_cell_centered() && S_fine.is_cell_centered()) ||
729 (S_crse.is_nodal() && S_fine.is_nodal()));
730
731 using value_type = typename FAB::value_type;
732
733 bool is_cell_centered = S_crse.is_cell_centered();
734
735 //
736 // Coarsen() the fine stuff on processors owning the fine data.
737 //
738 BoxArray crse_S_fine_BA = S_fine.boxArray(); crse_S_fine_BA.coarsen(ratio);
739
740 if (crse_S_fine_BA == S_crse.boxArray() && S_fine.DistributionMap() == S_crse.DistributionMap())
741 {
742#ifdef AMREX_USE_GPU
743 if (Gpu::inLaunchRegion() && S_crse.isFusingCandidate()) {
744 auto const& crsema = S_crse.arrays();
745 auto const& finema = S_fine.const_arrays();
746 if (is_cell_centered) {
747 ParallelFor(S_crse, IntVect(0), ncomp,
748 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
749 {
750 amrex_avgdown(i,j,k,n,crsema[box_no],finema[box_no],scomp,scomp,ratio);
751 });
752 } else {
753 ParallelFor(S_crse, IntVect(0), ncomp,
754 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
755 {
756 amrex_avgdown_nodes(i,j,k,n,crsema[box_no],finema[box_no],scomp,scomp,ratio);
757 });
758 }
759 if (!Gpu::inNoSyncRegion()) {
761 }
762 } else
763#endif
764 {
765#ifdef AMREX_USE_OMP
766#pragma omp parallel if (Gpu::notInLaunchRegion())
767#endif
768 for (MFIter mfi(S_crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
769 {
770 // NOTE: The tilebox is defined at the coarse level.
771 const Box& bx = mfi.tilebox();
772 Array4<value_type> const& crsearr = S_crse.array(mfi);
773 Array4<value_type const> const& finearr = S_fine.const_array(mfi);
774
775 if (is_cell_centered) {
776 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
777 {
778 amrex_avgdown(i,j,k,n,crsearr,finearr,scomp,scomp,ratio);
779 });
780 } else {
781 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
782 {
783 amrex_avgdown_nodes(i,j,k,n,crsearr,finearr,scomp,scomp,ratio);
784 });
785 }
786 }
787 }
788 }
789 else
790 {
791 FabArray<FAB> crse_S_fine(crse_S_fine_BA, S_fine.DistributionMap(), ncomp, 0, MFInfo(),DefaultFabFactory<FAB>());
792
793#ifdef AMREX_USE_GPU
794 if (Gpu::inLaunchRegion() && crse_S_fine.isFusingCandidate()) {
795 auto const& crsema = crse_S_fine.arrays();
796 auto const& finema = S_fine.const_arrays();
797 if (is_cell_centered) {
798 ParallelFor(crse_S_fine, IntVect(0), ncomp,
799 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
800 {
801 amrex_avgdown(i,j,k,n,crsema[box_no],finema[box_no],0,scomp,ratio);
802 });
803 } else {
804 ParallelFor(crse_S_fine, IntVect(0), ncomp,
805 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
806 {
807 amrex_avgdown_nodes(i,j,k,n,crsema[box_no],finema[box_no],0,scomp,ratio);
808 });
809 }
810 if (!Gpu::inNoSyncRegion()) {
812 }
813 } else
814#endif
815 {
816#ifdef AMREX_USE_OMP
817#pragma omp parallel if (Gpu::notInLaunchRegion())
818#endif
819 for (MFIter mfi(crse_S_fine,TilingIfNotGPU()); mfi.isValid(); ++mfi)
820 {
821 // NOTE: The tilebox is defined at the coarse level.
822 const Box& bx = mfi.tilebox();
823 Array4<value_type> const& crsearr = crse_S_fine.array(mfi);
824 Array4<value_type const> const& finearr = S_fine.const_array(mfi);
825
826 // NOTE: We copy from component scomp of the fine fab into component 0 of the crse fab
827 // because the crse fab is a temporary which was made starting at comp 0, it is
828 // not part of the actual crse multifab which came in.
829
830 if (is_cell_centered) {
831 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
832 {
833 amrex_avgdown(i,j,k,n,crsearr,finearr,0,scomp,ratio);
834 });
835 } else {
836 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
837 {
838 amrex_avgdown_nodes(i,j,k,n,crsearr,finearr,0,scomp,ratio);
839 });
840 }
841 }
842 }
843
844 S_crse.ParallelCopy(crse_S_fine,0,scomp,ncomp);
845 }
846}
847
848
849
850
851
859template <typename F>
860Real
861NormHelper (const MultiFab& x, int xcomp,
862 const MultiFab& y, int ycomp,
863 F const& f,
864 int numcomp, IntVect nghost, bool local)
865{
866 BL_ASSERT(x.boxArray() == y.boxArray());
867 BL_ASSERT(x.DistributionMap() == y.DistributionMap());
868 BL_ASSERT(x.nGrowVect().allGE(nghost) && y.nGrowVect().allGE(nghost));
869
870 Real sm = Real(0.0);
871#ifdef AMREX_USE_GPU
872 if (Gpu::inLaunchRegion()) {
873 auto const& xma = x.const_arrays();
874 auto const& yma = y.const_arrays();
876 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple<Real>
877 {
878 Real t = Real(0.0);
879 auto const& xfab = xma[box_no];
880 auto const& yfab = yma[box_no];
881 for (int n = 0; n < numcomp; ++n) {
882 t += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
883 }
884 return t;
885 });
886 } else
887#endif
888 {
889#ifdef AMREX_USE_OMP
890#pragma omp parallel if (!system::regtest_reduction) reduction(+:sm)
891#endif
892 for (MFIter mfi(x,true); mfi.isValid(); ++mfi)
893 {
894 Box const& bx = mfi.growntilebox(nghost);
895 Array4<Real const> const& xfab = x.const_array(mfi);
896 Array4<Real const> const& yfab = y.const_array(mfi);
897 AMREX_LOOP_4D(bx, numcomp, i, j, k, n,
898 {
899 sm += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
900 });
901 }
902 }
903
904 if (!local) {
906 }
907
908 return sm;
909}
910
919template <typename MMF, typename Pred, typename F>
920Real
921NormHelper (const MMF& mask,
922 const MultiFab& x, int xcomp,
923 const MultiFab& y, int ycomp,
924 Pred const& pf,
925 F const& f,
926 int numcomp, IntVect nghost, bool local)
927{
928 BL_ASSERT(x.boxArray() == y.boxArray());
929 BL_ASSERT(x.boxArray() == mask.boxArray());
930 BL_ASSERT(x.DistributionMap() == y.DistributionMap());
931 BL_ASSERT(x.DistributionMap() == mask.DistributionMap());
932 BL_ASSERT(x.nGrowVect().allGE(nghost) && y.nGrowVect().allGE(nghost));
933 BL_ASSERT(mask.nGrowVect().allGE(nghost));
934
935 Real sm = Real(0.0);
936#ifdef AMREX_USE_GPU
937 if (Gpu::inLaunchRegion()) {
938 auto const& xma = x.const_arrays();
939 auto const& yma = y.const_arrays();
940 auto const& mma = mask.const_arrays();
942 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple<Real>
943 {
944 Real t = Real(0.0);
945 if (pf(mma[box_no](i,j,k))) {
946 auto const& xfab = xma[box_no];
947 auto const& yfab = yma[box_no];
948 for (int n = 0; n < numcomp; ++n) {
949 t += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
950 }
951 }
952 return t;
953 });
954 } else
955#endif
956 {
957#ifdef AMREX_USE_OMP
958#pragma omp parallel if (!system::regtest_reduction) reduction(+:sm)
959#endif
960 for (MFIter mfi(x,true); mfi.isValid(); ++mfi)
961 {
962 Box const& bx = mfi.growntilebox(nghost);
963 Array4<Real const> const& xfab = x.const_array(mfi);
964 Array4<Real const> const& yfab = y.const_array(mfi);
965 auto const& mfab = mask.const_array(mfi);
966 AMREX_LOOP_4D(bx, numcomp, i, j, k, n,
967 {
968 if (pf(mfab(i,j,k))) {
969 sm += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
970 }
971 });
972 }
973 }
974
975 if (!local) {
977 }
978
979 return sm;
980}
981
982template <FabArrayType CMF, FabArrayType FMF>
983void average_face_to_cellcenter (CMF& cc, int dcomp,
985 int ngrow)
986{
987 IntVect ng_vect(ngrow);
988 average_face_to_cellcenter(cc, dcomp, fc, ng_vect);
989}
990
991template <FabArrayType CMF, FabArrayType FMF>
992void average_face_to_cellcenter (CMF& cc, int dcomp,
994 IntVect const& ng_vect)
995{
996 AMREX_ASSERT(cc.nComp() >= dcomp + AMREX_SPACEDIM);
997 AMREX_ASSERT(fc[0]->nComp() == 1);
998
999#ifdef AMREX_USE_GPU
1000 if (Gpu::inLaunchRegion() && cc.isFusingCandidate()) {
1001 auto const& ccma = cc.arrays();
1002 AMREX_D_TERM(auto const& fxma = fc[0]->const_arrays();,
1003 auto const& fyma = fc[1]->const_arrays();,
1004 auto const& fzma = fc[2]->const_arrays(););
1005 ParallelFor(cc, ng_vect,
1006 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
1007 {
1008#if (AMREX_SPACEDIM == 1)
1009 GeometryData gd{};
1010 gd.coord = 0;
1011#endif
1012 amrex_avg_fc_to_cc(i,j,k, ccma[box_no], AMREX_D_DECL(fxma[box_no],
1013 fyma[box_no],
1014 fzma[box_no]),
1015 dcomp
1016#if (AMREX_SPACEDIM == 1)
1017 , gd
1018#endif
1019 );
1020 });
1021 if (!Gpu::inNoSyncRegion()) {
1023 }
1024 } else
1025#endif
1026 {
1027#ifdef AMREX_USE_OMP
1028#pragma omp parallel if (Gpu::notInLaunchRegion())
1029#endif
1030 for (MFIter mfi(cc,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1031 {
1032 const Box bx = mfi.growntilebox(ng_vect);
1033 auto const& ccarr = cc.array(mfi);
1034 AMREX_D_TERM(auto const& fxarr = fc[0]->const_array(mfi);,
1035 auto const& fyarr = fc[1]->const_array(mfi);,
1036 auto const& fzarr = fc[2]->const_array(mfi););
1037
1038#if (AMREX_SPACEDIM == 1)
1040 {
1041 GeometryData gd;
1042 gd.coord = 0;
1043 amrex_avg_fc_to_cc(i,j,k, ccarr, fxarr, dcomp, gd);
1044 });
1045#else
1047 {
1048 amrex_avg_fc_to_cc(i,j,k, ccarr, AMREX_D_DECL(fxarr,fyarr,fzarr), dcomp);
1049 });
1050#endif
1051 }
1052 }
1053}
1054
1055template <FabArrayType MF>
1057 const Vector<MF*>& crse,
1058 const IntVect& ratio, int ngcrse)
1059{
1060 AMREX_ASSERT(fine.size() == AMREX_SPACEDIM && crse.size() == AMREX_SPACEDIM);
1062 {{AMREX_D_DECL(fine[0],fine[1],fine[2])}},
1064 {{AMREX_D_DECL(crse[0],crse[1],crse[2])}},
1065 ratio, ngcrse);
1066}
1067
1068template <FabArrayType MF>
1070 const Vector<MF*>& crse, int ratio, int ngcrse)
1071{
1072 average_down_faces(fine,crse,IntVect{ratio},ngcrse);
1073}
1074
1075template <FabArrayType MF>
1078 int ratio, int ngcrse)
1079{
1080 average_down_faces(fine,crse,IntVect{ratio},ngcrse);
1081}
1082
1083template <FabArrayType MF>
1086 const IntVect& ratio, int ngcrse)
1087{
1088 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1089 average_down_faces(*fine[idim], *crse[idim], ratio, ngcrse);
1090 }
1091}
1092
1093template <BaseFabType FAB>
1095 const IntVect& ratio, int ngcrse)
1096{
1097 BL_PROFILE("average_down_faces");
1098
1099 AMREX_ASSERT(crse.nComp() == fine.nComp());
1100 AMREX_ASSERT(fine.ixType() == crse.ixType());
1101 const auto type = fine.ixType();
1102 int dir;
1103 for (dir = 0; dir < AMREX_SPACEDIM; ++dir) {
1104 if (type.nodeCentered(dir)) { break; }
1105 }
1106 auto tmptype = type;
1107 tmptype.unset(dir);
1108 if (dir >= AMREX_SPACEDIM || !tmptype.cellCentered()) {
1109 amrex::Abort("average_down_faces: not face index type");
1110 }
1111 const int ncomp = crse.nComp();
1112 if (isMFIterSafe(fine, crse))
1113 {
1114#ifdef AMREX_USE_GPU
1115 if (Gpu::inLaunchRegion() && crse.isFusingCandidate()) {
1116 auto const& crsema = crse.arrays();
1117 auto const& finema = fine.const_arrays();
1118 ParallelFor(crse, IntVect(ngcrse), ncomp,
1119 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
1120 {
1121 amrex_avgdown_faces(i,j,k,n, crsema[box_no], finema[box_no], 0, 0, ratio, dir);
1122 });
1123 if (!Gpu::inNoSyncRegion()) {
1125 }
1126 } else
1127#endif
1128 {
1129#ifdef AMREX_USE_OMP
1130#pragma omp parallel if (Gpu::notInLaunchRegion())
1131#endif
1132 for (MFIter mfi(crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1133 {
1134 const Box& bx = mfi.growntilebox(ngcrse);
1135 auto const& crsearr = crse.array(mfi);
1136 auto const& finearr = fine.const_array(mfi);
1137 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
1138 {
1139 amrex_avgdown_faces(i,j,k,n, crsearr, finearr, 0, 0, ratio, dir);
1140 });
1141 }
1142 }
1143 }
1144 else
1145 {
1146 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
1147 ncomp, ngcrse, MFInfo(), DefaultFabFactory<FAB>());
1148 average_down_faces(fine, ctmp, ratio, ngcrse);
1149 crse.ParallelCopy(ctmp,0,0,ncomp,ngcrse,ngcrse);
1150 }
1151}
1152
1153template <FabArrayType MF>
1156 const IntVect& ratio, const Geometry& crse_geom)
1157{
1158 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1159 average_down_faces(*fine[idim], *crse[idim], ratio, crse_geom);
1160 }
1161}
1162
1163template <BaseFabType FAB>
1165 const IntVect& ratio, const Geometry& crse_geom)
1166{
1167 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
1168 crse.nComp(), 0);
1169 average_down_faces(fine, ctmp, ratio, 0);
1170 crse.ParallelCopy(ctmp,0,0,crse.nComp(),0,0,crse_geom.periodicity());
1171}
1172
1173template <FabArrayType MF>
1175{
1176 using T = typename MF::value_type;
1177 const int ncomp = mf.nComp();
1178 Gpu::DeviceVector<T> dv(ncomp);
1179 auto* dp = dv.data();
1180 bool found = false;
1181 auto loc = cell.dim3();
1182 for (MFIter mfi(mf); mfi.isValid() && !found; ++mfi)
1183 {
1184 Box const& box = mfi.validbox();
1185 if (box.contains(cell)) {
1186 found = true;
1187 auto const& fab = mf.const_array(mfi);
1188 amrex::ParallelFor(1, [=] AMREX_GPU_DEVICE (int) noexcept
1189 {
1190 for (int n = 0; n < ncomp; ++n) {
1191 dp[n] = fab(loc.x,loc.y,loc.z,n);
1192 }
1193 });
1194 }
1195 }
1196 Vector<T> hv;
1197 if (found) {
1198 hv.resize(ncomp);
1199 Gpu::copy(Gpu::deviceToHost, dv.begin(), dv.end(), hv.begin());
1200 }
1201 return hv;
1202}
1203
1204template <FabArrayType MF>
1205MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box const& bnd_bx)
1206{
1207 bool do_bnd = (!bnd_bx.isEmpty());
1208
1209 BoxArray const& ba = mf.boxArray();
1210 DistributionMapping const& dm = mf.DistributionMap();
1211 const auto nboxes = static_cast<int>(ba.size());
1212
1213 BoxList bl(ba.ixType());
1214 Vector<int> procmap;
1215 Vector<int> index_map;
1216 if (!do_bnd) {
1217 for (int i = 0; i < nboxes; ++i) {
1218 Box const& b = ba[i];
1219 IntVect lo = cell;
1220 lo[dir] = b.smallEnd(dir);
1221 if (b.contains(lo)) {
1222 IntVect hi = lo;
1223 hi[dir] = b.bigEnd(dir);
1224 Box b1d(lo,hi,b.ixType());
1225 bl.push_back(b1d);
1226 procmap.push_back(dm[i]);
1227 index_map.push_back(i);
1228 }
1229 }
1230 } else {
1231 for (int i = 0; i < nboxes; ++i) {
1232 Box const& b = ba[i];
1233 Box const& b1d = bnd_bx & b;
1234 if (b1d.ok()) {
1235 bl.push_back(b1d);
1236 procmap.push_back(dm[i]);
1237 index_map.push_back(i);
1238 }
1239 }
1240 }
1241
1242 if (bl.isEmpty()) {
1243 return MF();
1244 } else {
1245 BoxArray rba(std::move(bl));
1246 DistributionMapping rdm(std::move(procmap));
1247 MF rmf(rba, rdm, mf.nComp(), IntVect(0),
1248 MFInfo().SetArena(mf.arena()));
1249#ifdef AMREX_USE_OMP
1250#pragma omp parallel if (Gpu::notInLaunchRegion())
1251#endif
1252 for (MFIter mfi(rmf); mfi.isValid(); ++mfi) {
1253 Box const& b = mfi.validbox();
1254 auto const& dfab = rmf.array(mfi);
1255 auto const& sfab = mf.const_array(index_map[mfi.index()]);
1256 amrex::ParallelFor(b, mf.nComp(),
1257 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
1258 {
1259 dfab(i,j,k,n) = sfab(i,j,k,n);
1260 });
1261 }
1262 return rmf;
1263 }
1264}
1265
1267namespace detail {
1268template <typename Op, typename T, typename F>
1269void reduce_to_plane (Array4<T> const& ar, int direction, Box const& bx, int box_no,
1270 F const& f)
1271{
1272#if defined(AMREX_USE_GPU)
1273 Box b2d = bx;
1274 b2d.setRange(direction,0);
1275 const auto blo = amrex::lbound(bx);
1276 const auto len = amrex::length(bx);
1277 constexpr int nthreads = 128;
1278 auto nblocks = static_cast<int>(b2d.numPts());
1279#ifdef AMREX_USE_SYCL
1280 constexpr std::size_t shared_mem_bytes = sizeof(T)*Gpu::Device::warp_size;
1281 amrex::launch<nthreads>(nblocks, shared_mem_bytes, Gpu::gpuStream(),
1282 [=] AMREX_GPU_DEVICE (Gpu::Handler const& h)
1283 {
1284 int bid = h.blockIdx();
1285 int tid = h.threadIdx();
1286#else
1287 amrex::launch<nthreads>(nblocks, Gpu::gpuStream(),
1288 [=] AMREX_GPU_DEVICE ()
1289 {
1290 int bid = blockIdx.x;
1291 int tid = threadIdx.x;
1292#endif
1293 T tmp;
1294 Op().init(tmp);
1295 T* p;
1296 if (direction == 0) {
1297 int k = bid / len.y;
1298 int j = bid - k*len.y;
1299 k += blo.z;
1300 j += blo.y;
1301 for (int i = blo.x + tid; i < blo.x+len.x; i += nthreads) {
1302 Op().local_update(tmp, f(box_no,i,j,k));
1303 }
1304 p = ar.ptr(0,j,k);
1305 } else if (direction == 1) {
1306 int k = bid / len.x;
1307 int i = bid - k*len.x;
1308 k += blo.z;
1309 i += blo.x;
1310 for (int j = blo.y + tid; j < blo.y+len.y; j += nthreads) {
1311 Op().local_update(tmp, f(box_no,i,j,k));
1312 }
1313 p = ar.ptr(i,0,k);
1314 } else {
1315 int j = bid / len.x;
1316 int i = bid - j*len.x;
1317 j += blo.y;
1318 i += blo.x;
1319 for (int k = blo.z + tid; k < blo.z+len.z; k += nthreads) {
1320 Op().local_update(tmp, f(box_no,i,j,k));
1321 }
1322 p = ar.ptr(i,j,0);
1323 }
1324#ifdef AMREX_USE_SYCL
1325 Op().template parallel_update<T>(*p, tmp, h);
1326#else
1327 Op().template parallel_update<T,nthreads>(*p, tmp);
1328#endif
1329 });
1330#else
1331 // CPU
1332 if (direction == 0) {
1333 AMREX_LOOP_3D(bx, i, j, k,
1334 {
1335 Op().local_update(ar(0,j,k), f(box_no,i,j,k));
1336 });
1337 } else if (direction == 1) {
1338 AMREX_LOOP_3D(bx, i, j, k,
1339 {
1340 Op().local_update(ar(i,0,k), f(box_no,i,j,k));
1341 });
1342 } else {
1343 AMREX_LOOP_3D(bx, i, j, k,
1344 {
1345 Op().local_update(ar(i,j,0), f(box_no,i,j,k));
1346 });
1347 }
1348#endif
1349}
1350}
1352
1353template <typename Op, typename T, BaseFabType FAB, typename F>
1354#ifndef AMREX_USE_CUDA
1355requires (IsCallableR<T,F,int,int,int,int>::value)
1356#endif
1357BaseFab<T>
1358ReduceToPlane (int direction, Box const& a_domain, FabArray<FAB> const& mf, F const& f)
1359{
1360 Box const domain = amrex::convert(a_domain, mf.ixType());
1361
1362 Box domain2d = domain;
1363 domain2d.setRange(direction, 0);
1364
1365 T initval;
1366 Op().init(initval);
1367
1368 BaseFab<T> r(domain2d);
1369 r.template setVal<RunOn::Device>(initval);
1370 auto const& ar = r.array();
1371
1372 for (MFIter mfi(mf,MFItInfo().UseDefaultStream().DisableDeviceSync());
1373 mfi.isValid(); ++mfi)
1374 {
1375 Box bx = mfi.validbox() & domain;
1376 if (bx.ok()) {
1377 int box_no = mfi.LocalIndex();
1378 detail::reduce_to_plane<Op, T>(ar, direction, bx, box_no, f);
1379 }
1380 }
1381 if (!Gpu::inNoSyncRegion()) {
1383 }
1384
1385 return r;
1386}
1387
1389namespace detail {
1390template <typename Op, MultiFabLike FA, typename F>
1391FA reduce_to_plane (int direction, Box const& domain, FA const& mf, F const& f,
1392 IntVect const& nghost = IntVect(0))
1393{
1394 using T = typename FA::value_type;
1395
1396 Box const ndomain = amrex::convert(domain, mf.ixType());
1397
1398 auto npts = amrex::convert(mf.boxArray(),IntVect(0)).numPts();
1399 if (npts != amrex::convert(domain, amrex::IntVect(0)).numPts()) {
1400 amrex::Abort("ReduceToPlaneMF: mf's BoxArray must have a rectangular domain.");
1401 }
1402
1403 // the reduction reads mf's guard cells, so it cannot include more guard
1404 // cells than mf actually has
1405 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(nghost.allLE(mf.nGrowVect()),
1406 "ReduceToPlaneMF2: nghost must not exceed mf.nGrowVect().");
1407
1408 T initval;
1409 Op().init(initval);
1410
1411 // guard cells in the reduction direction are summed into the plane; guard
1412 // cells in the transverse directions are stored in the result's ghost cells
1413 // (and folded into the valid region by the caller's ParallelAdd)
1414 IntVect ng2d = nghost;
1415 ng2d[direction] = 0;
1416 Box const gdomain = amrex::grow(ndomain, nghost);
1417
1418 BoxList bl = mf.boxArray().boxList();
1419 for (auto& b : bl) {
1420 b.setRange(direction, 0);
1421 }
1422 BoxArray ba(std::move(bl));
1423 FA tmpfa(ba, mf.DistributionMap(), 1, ng2d);
1424 tmpfa.setVal(initval);
1425
1426 for (MFIter mfi(mf); mfi.isValid(); ++mfi)
1427 {
1428 Box bx = amrex::grow(mfi.validbox(), nghost) & gdomain;
1429 if (bx.ok()) {
1430 int box_no = mfi.LocalIndex();
1431 detail::reduce_to_plane<Op, T>(tmpfa.array(mfi), direction, bx, box_no, f);
1432 }
1433 }
1434
1435 return tmpfa;
1436}
1437
1438}
1440
1441template <typename Op, MultiFabLike FA, typename F>
1442#ifndef AMREX_USE_CUDA
1443requires (IsCallableR<typename FA::value_type,
1444 F,int,int,int,int>::value)
1445#endif
1446FA ReduceToPlaneMF (int direction, Box const& domain, FA const& mf, F const& f)
1447{
1448 auto [fa3, fa2] = ReduceToPlaneMF2<Op>(direction, domain, mf, f);
1449 fa3.ParallelCopy(fa2);
1450 return std::move(fa3);
1451}
1452
1453template <typename Op, MultiFabLike FA, typename F>
1454#ifndef AMREX_USE_CUDA
1455requires (IsCallableR<typename FA::value_type,
1456 F,int,int,int,int>::value)
1457#endif
1458std::pair<FA,FA>
1459ReduceToPlaneMF2 (int direction, Box const& domain, FA const& mf, F const& f,
1460 IntVect const& nghost, Periodicity const& period)
1461{
1462 using T = typename FA::value_type;
1463
1464 T initval;
1465 Op().init(initval);
1466
1467 auto tmpmf = detail::reduce_to_plane<Op>(direction, domain, mf, f, nghost);
1468
1469 BoxList bl2d(mf.ixType());
1470 Vector<int> procmap2d;
1471 auto const& ba3d = mf.boxArray();
1472 auto const& dm3d = mf.DistributionMap();
1473 int dlo = domain.smallEnd(direction);
1474 for (int i = 0, N = mf.size(); i < N; ++i) {
1475 Box b = ba3d[i];
1476 if (b.smallEnd(direction) <= dlo && dlo <= b.bigEnd(direction)) {
1477 b.setRange(direction, 0);
1478 bl2d.push_back(b);
1479 procmap2d.push_back(dm3d[i]);
1480 }
1481 }
1482
1483 BoxArray ba2d(std::move(bl2d));
1484 DistributionMapping dm2d(std::move(procmap2d));
1485
1486 FA mf2d(ba2d, dm2d, 1, 0);
1487 mf2d.setVal(initval);
1488
1489 static_assert(std::is_same_v<Op, ReduceOpSum>, "Currently only ReduceOpSum is supported.");
1490 // fold the transverse guard cells of the local results into the valid
1491 // region of the unique result (sums the quantity across box seams)
1492 IntVect ng2d = nghost;
1493 ng2d[direction] = 0;
1494 mf2d.ParallelAdd(tmpmf, 0, 0, 1, ng2d, IntVect(0), IntVect(0), period);
1495
1496 return std::make_pair(std::move(tmpmf), std::move(mf2d));
1497}
1498
1499template <typename Op, MultiFabLike FA, typename F>
1500#ifndef AMREX_USE_CUDA
1501requires (IsCallableR<typename FA::value_type,
1502 F,int,int,int,int>::value)
1503#endif
1504std::pair<FA,FA>
1505ReduceToPlaneMF2Patchy (int direction, Box const& domain, FA const& mf,
1506 IntVect const& plane_max_grid_size, F const& f)
1507{
1508 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mf.is_cell_centered(),
1509 "ReduceToPlaneMF2Patchy: only cell-centered data is supported.");
1510
1511 Box const ndomain = amrex::convert(domain, mf.ixType());
1512 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ndomain.contains(mf.boxArray().minimalBox()),
1513 "ReduceToPlaneMF2Patchy: subdomains are not supported; domain must cover the full extent of mf.boxArray().");
1514
1515 using T = typename FA::value_type;
1516 T initval;
1517 Op().init(initval);
1518
1519 BoxList bl_patch2d = mf.boxArray().boxList();
1520 for (auto& b : bl_patch2d) {
1521 b.setRange(direction, 0);
1522 }
1523 BoxArray ba_patch2d(std::move(bl_patch2d));
1524 FA tmpmf(ba_patch2d, mf.DistributionMap(), 1, 0);
1525 tmpmf.setVal(initval);
1526
1527 for (MFIter mfi(mf); mfi.isValid(); ++mfi)
1528 {
1529 Box bx = mfi.validbox() & ndomain;
1530 if (bx.ok()) {
1531 int box_no = mfi.LocalIndex();
1532 detail::reduce_to_plane<Op, T>(tmpmf.array(mfi), direction, bx, box_no, f);
1533 }
1534 }
1535
1536 BoxList bl_unique = ba_patch2d.boxList();
1537 bl_unique.simplify(); // reduces work for .removeOverlap()
1538 BoxArray ba2d(std::move(bl_unique));
1539 ba2d.removeOverlap();
1540
1541 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(plane_max_grid_size.allGE(1),
1542 "ReduceToPlaneMF2Patchy: plane_max_grid_size must be >= 1 in every direction.");
1543 IntVect mgs2d = plane_max_grid_size;
1544 mgs2d[direction] = 1;
1545 ba2d.maxSize(mgs2d);
1546
1547 DistributionMapping dm2d(ba2d);
1548
1549 FA mf2d(ba2d, dm2d, 1, 0);
1550 mf2d.setVal(initval);
1551
1552 static_assert(std::is_same_v<Op, ReduceOpSum>, "Currently only ReduceOpSum is supported.");
1553 mf2d.ParallelAdd(tmpmf);
1554
1555 return std::make_pair(std::move(tmpmf), std::move(mf2d));
1556}
1557
1558}
1559
1560#endif
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#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:1139
Real * pdst
Definition AMReX_HypreMLABecLap.cpp:1140
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
#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:194
Array4< T const > array() const noexcept
Definition AMReX_BaseFab.H:387
A collection of Boxes stored in an Array.
Definition AMReX_BoxArray.H:564
IndexType ixType() const noexcept
Return index type of this BoxArray.
Definition AMReX_BoxArray.H:854
BoxList boxList() const
Create a BoxList from this BoxArray.
Definition AMReX_BoxArray.cpp:949
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:1456
BoxArray & coarsen(int refinement_ratio)
Coarsen each Box in the BoxArray to the specified ratio.
Definition AMReX_BoxArray.cpp:672
Long size() const noexcept
Return the number of boxes in the BoxArray.
Definition AMReX_BoxArray.H:611
BoxArray & maxSize(int block_size)
Forces each Box in BoxArray to have sides <= block_size.
Definition AMReX_BoxArray.cpp:553
A class for managing a List of Boxes that share a common IndexType. This class implements operations ...
Definition AMReX_BoxList.H:52
bool isEmpty() const noexcept
Is this BoxList empty?
Definition AMReX_BoxList.H:135
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:652
void push_back(const Box &bn)
Append a Box to this BoxList.
Definition AMReX_BoxList.H:93
__host__ __device__ const IntVectND< dim > & bigEnd() const &noexcept
Return the inclusive upper bound of the box.
Definition AMReX_Box.H:124
__host__ __device__ bool isEmpty() const noexcept
Checks if it is an empty BoxND.
Definition AMReX_Box.H:208
__host__ __device__ bool contains(const IntVectND< dim > &p) const noexcept
Return true if argument is contained within BoxND.
Definition AMReX_Box.H:216
__host__ __device__ IndexTypeND< dim > ixType() const noexcept
Return the indexing type.
Definition AMReX_Box.H:136
__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:1117
__host__ __device__ bool ok() const noexcept
Checks if it is a proper BoxND (including a valid type).
Definition AMReX_Box.H:212
__host__ __device__ const IntVectND< dim > & smallEnd() const &noexcept
Return the inclusive lower bound of the box.
Definition AMReX_Box.H:112
Definition AMReX_FabFactory.H:76
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:43
IntVect nGrowVect() const noexcept
Definition AMReX_FabArrayBase.H:80
bool isFusingCandidate() const noexcept
Is this a good candidate for kernel fusing?
Definition AMReX_FabArrayBase.cpp:2705
bool is_cell_centered() const noexcept
This tests on whether the FabArray is cell-centered.
Definition AMReX_FabArrayBase.cpp:2699
bool is_nodal() const noexcept
This tests on whether the FabArray is fully nodal.
Definition AMReX_FabArrayBase.cpp:2687
IndexType ixType() const noexcept
Return index type.
Definition AMReX_FabArrayBase.H:86
const DistributionMapping & DistributionMap() const noexcept
Return constant reference to associated DistributionMapping.
Definition AMReX_FabArrayBase.H:130
int nComp() const noexcept
Return number of variables (aka components) associated with each point.
Definition AMReX_FabArrayBase.H:83
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:95
An Array of FortranArrayBox(FAB)-like Objects.
Definition AMReX_FabArray.H:344
void ParallelCopy(const FabArray< FAB > &src, const Periodicity &period=Periodicity::NonPeriodic(), CpOp op=FabArrayBase::COPY)
Definition AMReX_FabArray.H:873
MultiArray4< typename FabArray< FAB >::value_type > arrays() noexcept
Definition AMReX_FabArray.H:633
MultiArray4< typename FabArray< FAB >::value_type const > const_arrays() const noexcept
Definition AMReX_FabArray.H:647
Array4< typename FabArray< FAB >::value_type const > array(const MFIter &mfi) const noexcept
Definition AMReX_FabArray.H:561
Array4< typename FabArray< FAB >::value_type const > const_array(const MFIter &mfi) const noexcept
Definition AMReX_FabArray.H:585
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:75
Periodicity periodicity() const noexcept
Definition AMReX_Geometry.H:361
GPU-compatible tuple.
Definition AMReX_Tuple.H:98
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:172
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:674
iterator end() noexcept
Definition AMReX_PODVector.H:678
T * data() noexcept
Definition AMReX_PODVector.H:666
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:79
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:1373
__host__ __device__ Dim3 lbound(Array4< T > const &a) noexcept
Return the inclusive lower bounds of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1345
__host__ __device__ BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) coarsening ratio.
Definition AMReX_Box.H:1418
__host__ __device__ BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Grow BoxND in all directions by given amount.
Definition AMReX_Box.H:1289
std::array< T, N > Array
Definition AMReX_Array.H:26
void Sum(Gpu::DeviceVector< T > &v, MPI_Comm comm)
Definition AMReX_GpuParallelReduce.H:34
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:1505
__host__ __device__ BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Return a BoxND with different type.
Definition AMReX_Box.H:1567
void FillRandomNormal(MultiFab &mf, int scomp, int ncomp, Real mean, Real stddev)
Fill MultiFab with random numbers from normal distribution.
Definition AMReX_MultiFabUtil.cpp:1261
int nComp(FabArrayBase const &fa)
Definition AMReX_FabArrayBase.cpp:2852
void FillRandom(MultiFab &mf, int scomp, int ncomp)
Fill MultiFab with random numbers from uniform distribution.
Definition AMReX_MultiFabUtil.cpp:1248
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
Definition AMReX_BaseFabUtility.H:13
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:597
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:652
Vector< MultiFab > convexify(Vector< MultiFab const * > const &mf, Vector< IntVect > const &refinement_ratio)
Convexify AMR data.
Definition AMReX_MultiFabUtil.cpp:1274
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:1205
FA ReduceToPlaneMF(int direction, Box const &domain, FA const &mf, F const &f)
Reduce FabArray/MultiFab data to plane FabArray.
Definition AMReX_MultiFabUtil.H:1446
MultiFab periodicShift(MultiFab const &mf, IntVect const &offset, Periodicity const &period)
Periodic shift MultiFab.
Definition AMReX_MultiFabUtil.cpp:841
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:1358
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:1007
BoxND< 3 > Box
Box is an alias for amrex::BoxND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:30
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:1056
bool isMFIterSafe(const FabArrayBase &x, const FabArrayBase &y)
Definition AMReX_MFIter.H:252
MultiFab ToMultiFab(const iMultiFab &imf)
Convert iMultiFab to MultiFab.
Definition AMReX_MultiFabUtil.cpp:587
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:1174
void writeFabs(const MultiFab &mf, const std::string &name)
Write each fab individually.
Definition AMReX_MultiFabUtil.cpp:574
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:861
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:740
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:492
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:865
IntVectND< 3 > IntVect
IntVect is an alias for amrex::IntVectND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:33
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:1459
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:569
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 out message to cerr and exit via abort().
Definition AMReX.cpp:241
const int[]
Definition AMReX_BLProfiler.cpp:1664
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:1161
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:592
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:792
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 node-based MultiFab onto crse node-centered MultiFab.
Definition AMReX_MultiFabUtil.H:675
A multidimensional array accessor.
Definition AMReX_Array4.H:285
parallel copy or add
Definition AMReX_FabArrayBase.H:537
Definition AMReX_Geometry.H:26
int coord
Definition AMReX_Geometry.H:60
FabArray memory allocation information.
Definition AMReX_FabArray.H:68
Definition AMReX_MFIter.H:20
Struct for holding types.
Definition AMReX_TypeList.H:13