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_MultiFab.H>
6#include <AMReX_iMultiFab.H>
7#include <AMReX_LayoutData.H>
8#include <AMReX_MFIter.H>
9#include <AMReX_Array.H>
10#include <AMReX_Vector.H>
11#include <AMReX_MultiFabUtil_C.H>
12
13#include <AMReX_MultiFabUtilI.H>
14
15namespace amrex
16{
18 void average_node_to_cellcenter (MultiFab& cc, int dcomp,
19 const MultiFab& nd, int scomp,
20 int ncomp, int ngrow = 0);
21 void average_node_to_cellcenter (MultiFab& cc, int dcomp,
22 const MultiFab& nd, int scomp,
23 int ncomp, IntVect const& ng_vect);
24
31 void average_edge_to_cellcenter (MultiFab& cc, int dcomp,
32 const Vector<const MultiFab*>& edge,
33 int ngrow = 0);
34 void average_edge_to_cellcenter (MultiFab& cc, int dcomp,
35 const Vector<const MultiFab*>& edge,
36 IntVect const& ng_vect);
37
44 void average_face_to_cellcenter (MultiFab& cc, int dcomp,
45 const Vector<const MultiFab*>& fc,
46 int ngrow = 0);
47 void average_face_to_cellcenter (MultiFab& cc, int dcomp,
48 const Vector<const MultiFab*>& fc,
49 IntVect const& ng_vect);
50
52 template <typename CMF, typename FMF,
53 std::enable_if_t<IsFabArray_v<CMF> && IsFabArray_v<FMF>, int> = 0>
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 <typename CMF, typename FMF,
59 std::enable_if_t<IsFabArray_v<CMF> && IsFabArray_v<FMF>, int> = 0>
60 void average_face_to_cellcenter (CMF& cc, int dcomp,
61 const Array<const FMF*,AMREX_SPACEDIM>& fc,
62 IntVect const& ng_vect);
63
65 void average_face_to_cellcenter (MultiFab& cc,
66 const Vector<const MultiFab*>& fc,
67 const Geometry& geom);
69 void average_face_to_cellcenter (MultiFab& cc,
70 const Array<const MultiFab*,AMREX_SPACEDIM>& fc,
71 const Geometry& geom);
80 void average_cellcenter_to_face (const Vector<MultiFab*>& fc,
81 const MultiFab& cc,
82 const Geometry& geom,
83 int ncomp = 1,
84 bool use_harmonic_averaging = false,
85 int ngrow = 0);
87 void average_cellcenter_to_face (const Array<MultiFab*,AMREX_SPACEDIM>& fc,
88 const MultiFab& cc,
89 const Geometry& geom,
90 int ncomp = 1,
91 bool use_harmonic_averaging = false,
92 int ngrow = 0);
93
94 void average_cellcenter_to_face (const Array<MultiFab*,AMREX_SPACEDIM>& fc,
95 const MultiFab& cc,
96 const Geometry& geom,
97 int ncomp,
98 bool use_harmonic_averaging,
99 const Array<IntVect, AMREX_SPACEDIM>& ng_vects);
100
102 template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> = 0>
103 void average_down_faces (const Vector<const MF*>& fine,
104 const Vector<MF*>& crse,
105 const IntVect& ratio,
106 int ngcrse = 0);
108 template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> = 0>
109 void average_down_faces (const Vector<const MF*>& fine,
110 const Vector<MF*>& crse,
111 int ratio,
112 int ngcrse = 0);
114 template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> = 0>
115 void average_down_faces (const Array<const MF*,AMREX_SPACEDIM>& fine,
116 const Array<MF*,AMREX_SPACEDIM>& crse,
117 const IntVect& ratio,
118 int ngcrse = 0);
120 template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> = 0>
121 void average_down_faces (const Array<const MF*,AMREX_SPACEDIM>& fine,
122 const Array<MF*,AMREX_SPACEDIM>& crse,
123 int ratio,
124 int ngcrse = 0);
131 template <typename FAB>
132 void average_down_faces (const FabArray<FAB>& fine, FabArray<FAB>& crse,
133 const IntVect& ratio, int ngcrse=0);
134
135 // This version takes periodicity into account.
136 template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> = 0>
137 void average_down_faces (const Array<const MF*,AMREX_SPACEDIM>& fine,
138 const Array<MF*,AMREX_SPACEDIM>& crse,
139 const IntVect& ratio, const Geometry& crse_geom);
140 // This version takes periodicity into account.
141 template <typename FAB>
142 void average_down_faces (const FabArray<FAB>& fine, FabArray<FAB>& crse,
143 const IntVect& ratio, const Geometry& crse_geom);
144
146 void average_down_edges (const Vector<const MultiFab*>& fine,
147 const Vector<MultiFab*>& crse,
148 const IntVect& ratio,
149 int ngcrse = 0);
150 void average_down_edges (const Array<const MultiFab*,AMREX_SPACEDIM>& fine,
151 const Array<MultiFab*,AMREX_SPACEDIM>& crse,
152 const IntVect& ratio,
153 int ngcrse = 0);
157 void average_down_edges (const MultiFab& fine, MultiFab& crse,
158 const IntVect& ratio, int ngcrse=0);
159
161 template <typename FAB>
162 void average_down_nodal (const FabArray<FAB>& S_fine,
163 FabArray<FAB>& S_crse,
164 const IntVect& ratio,
165 int ngcrse = 0,
166 bool mfiter_is_definitely_safe=false);
167
174 void average_down (const MultiFab& S_fine, MultiFab& S_crse,
175 const Geometry& fgeom, const Geometry& cgeom,
176 int scomp, int ncomp, const IntVect& ratio);
177 void average_down (const MultiFab& S_fine, MultiFab& S_crse,
178 const Geometry& fgeom, const Geometry& cgeom,
179 int scomp, int ncomp, int rr);
180
184 template<typename FAB>
185 void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse,
186 int scomp, int ncomp, const IntVect& ratio);
187 template<typename FAB>
188 void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse,
189 int scomp, int ncomp, int rr);
190
193 void sum_fine_to_coarse (const MultiFab& S_Fine, MultiFab& S_crse,
194 int scomp, int ncomp,
195 const IntVect& ratio,
196 const Geometry& cgeom, const Geometry& fgeom);
197
199 void print_state (const MultiFab& mf, const IntVect& cell, int n=-1,
200 const IntVect& ng = IntVect::TheZeroVector());
201
203 void writeFabs (const MultiFab& mf, const std::string& name);
204 void writeFabs (const MultiFab& mf, int comp, int ncomp, const std::string& name);
205
208 std::unique_ptr<MultiFab> get_slice_data(int dir, Real coord,
209 const MultiFab& cc,
210 const Geometry& geom, int start_comp, int ncomp,
211 bool interpolate=false,
212 RealBox const& bnd_rbx = RealBox());
213
221 template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> FOO = 0>
222 Vector<typename MF::value_type> get_cell_data (MF const& mf, IntVect const& cell);
223
230 template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> FOO = 0>
231 MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box const& bnd_bx = Box());
232
236 template <typename FAB>
237 iMultiFab makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
238 int crse_value = 0, int fine_value = 1,
239 MFInfo const& info = MFInfo());
240 iMultiFab makeFineMask (const BoxArray& cba, const DistributionMapping& cdm,
241 const BoxArray& fba, const IntVect& ratio,
242 int crse_value = 0, int fine_value = 1,
243 MFInfo const& info = MFInfo());
244 template <typename FAB>
245 iMultiFab makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
246 Periodicity const& period, int crse_value, int fine_value,
247 MFInfo const& info = MFInfo());
248 iMultiFab makeFineMask (const BoxArray& cba, const DistributionMapping& cdm,
249 const IntVect& cnghost, const BoxArray& fba, const IntVect& ratio,
250 Periodicity const& period, int crse_value, int fine_value,
251 MFInfo const& info = MFInfo());
252 template <typename FAB>
253 iMultiFab makeFineMask (const FabArray<FAB>& cmf, const FabArray<FAB>& fmf,
254 const IntVect& cnghost, const IntVect& ratio,
255 Periodicity const& period, int crse_value, int fine_value,
256 MFInfo const& info = MFInfo());
257 template <typename FAB>
258 iMultiFab makeFineMask (const FabArray<FAB>& cmf, const FabArray<FAB>& fmf,
259 const IntVect& cnghost, const IntVect& ratio,
260 Periodicity const& period, int crse_value, int fine_value,
261 LayoutData<int>& has_cf,
262 MFInfo const& info = MFInfo());
263
264 MultiFab makeFineMask (const BoxArray& cba, const DistributionMapping& cdm,
265 const BoxArray& fba, const IntVect& ratio,
266 Real crse_value, Real fine_value,
267 MFInfo const& info = MFInfo());
268
270 void computeDivergence (MultiFab& divu, const Array<MultiFab const*,AMREX_SPACEDIM>& umac,
271 const Geometry& geom);
272
274 void computeGradient (MultiFab& grad, const Array<MultiFab const*,AMREX_SPACEDIM>& umac,
275 const Geometry& geom);
276
278 MultiFab ToMultiFab (const iMultiFab& imf);
280 FabArray<BaseFab<Long> > ToLongMultiFab (const iMultiFab& imf);
281
283 MultiFab periodicShift (MultiFab const& mf, IntVect const& offset,
284 Periodicity const& period);
285
287 template <typename T, typename U>
288 T cast (U const& mf_in)
289 {
290 T mf_out(mf_in.boxArray(), mf_in.DistributionMap(), mf_in.nComp(), mf_in.nGrowVect());
291
292#ifdef AMREX_USE_OMP
293#pragma omp parallel if (Gpu::notInLaunchRegion())
294#endif
295 for (MFIter mfi(mf_in); mfi.isValid(); ++mfi)
296 {
297 const Long n = mfi.fabbox().numPts() * mf_in.nComp();
298 auto * pdst = mf_out[mfi].dataPtr();
299 auto const* psrc = mf_in [mfi].dataPtr();
301 {
302 pdst[i] = static_cast<typename T::value_type>(psrc[i]); // NOLINT(bugprone-signed-char-misuse)
303 });
304 }
305 return mf_out;
306 }
307
368 template <typename Op, typename T, typename FAB, typename F,
369 std::enable_if_t<IsBaseFab<FAB>::value
370#ifndef AMREX_USE_CUDA
371 && IsCallableR<T,F,int,int,int,int>::value
372#endif
373 , int> FOO = 0>
374 BaseFab<T>
375 ReduceToPlane (int direction, Box const& domain, FabArray<FAB> const& mf, F const& f);
376
421 template <typename Op, typename FA, typename F,
422 std::enable_if_t<IsMultiFabLike_v<FA>
423#ifndef AMREX_USE_CUDA
424 && IsCallableR<typename FA::value_type,
425 F,int,int,int,int>::value
426#endif
427 , int> FOO = 0>
428 FA ReduceToPlaneMF (int direction, Box const& domain, FA const& mf, F const& f);
429
472 template <typename Op, typename FA, typename F,
473 std::enable_if_t<IsMultiFabLike_v<FA>
474#ifndef AMREX_USE_CUDA
475 && IsCallableR<typename FA::value_type,
476 F,int,int,int,int>::value
477#endif
478 , int> FOO = 0>
479 std::pair<FA,FA>
480 ReduceToPlaneMF2 (int direction, Box const& domain, FA const& mf, F const& f);
481
509 template <typename Op, typename FA, typename F,
510 std::enable_if_t<IsMultiFabLike_v<FA>
511#ifndef AMREX_USE_CUDA
512 && IsCallableR<typename FA::value_type,
513 F,int,int,int,int>::value
514#endif
515 , int> FOO = 0>
516 std::pair<FA,FA>
517 ReduceToPlaneMF2Patchy (int direction, Box const& domain, FA const& mf,
518 IntVect const& plane_max_grid_size, F const& f);
519
535 Gpu::HostVector<Real> sumToLine (MultiFab const& mf, int icomp, int ncomp,
536 Box const& domain, int direction, bool local = false);
537
545 Real volumeWeightedSum (Vector<MultiFab const*> const& mf, int icomp,
546 Vector<Geometry> const& geom,
547 Vector<IntVect> const& ratio,
548 bool local = false);
549
563 void FourthOrderInterpFromFineToCoarse (MultiFab& cmf, int scomp, int ncomp,
564 MultiFab const& fmf,
565 IntVect const& ratio);
566
577 void FillRandom (MultiFab& mf, int scomp, int ncomp);
578
590 void FillRandomNormal (MultiFab& mf, int scomp, int ncomp, Real mean, Real stddev);
591
603 [[nodiscard]] Vector<MultiFab> convexify (Vector<MultiFab const*> const& mf,
604 Vector<IntVect> const& refinement_ratio);
605}
606
607namespace amrex {
608
609template <typename FAB>
610iMultiFab
611makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
612 int crse_value, int fine_value, MFInfo const& info)
613{
614 return makeFineMask(cmf.boxArray(), cmf.DistributionMap(), cmf.nGrowVect(),
615 fba, ratio, Periodicity::NonPeriodic(), crse_value, fine_value,
616 info);
617}
618
619template <typename FAB>
620iMultiFab
621makeFineMask (const FabArray<FAB>& cmf, const BoxArray& fba, const IntVect& ratio,
622 Periodicity const& period, int crse_value, int fine_value,
623 MFInfo const& info)
624{
625 return makeFineMask(cmf.boxArray(), cmf.DistributionMap(), cmf.nGrowVect(),
626 fba, ratio, period, crse_value, fine_value, info);
627}
628
629template <typename FAB>
630iMultiFab
632 const IntVect& cnghost, const IntVect& ratio,
633 Periodicity const& period, int crse_value, int fine_value,
634 MFInfo const& info)
635{
636 iMultiFab mask(cmf.boxArray(), cmf.DistributionMap(), 1, cnghost, info);
637 mask.setVal(crse_value);
638
639 iMultiFab foo(amrex::coarsen(fmf.boxArray(),ratio), fmf.DistributionMap(),
640 1, 0, MFInfo().SetAlloc(false));
641 const FabArrayBase::CPC& cpc = mask.getCPC(cnghost,foo,IntVect::TheZeroVector(),period);
642 mask.setVal(fine_value, cpc, 0, 1);
643
644 return mask;
645}
646
647template <typename FAB>
648iMultiFab
650 const IntVect& cnghost, const IntVect& ratio,
651 Periodicity const& period, int crse_value, int fine_value,
652 LayoutData<int>& has_cf, MFInfo const& info)
653{
654 iMultiFab mask(cmf.boxArray(), cmf.DistributionMap(), 1, cnghost, info);
655 mask.setVal(crse_value);
656
657 iMultiFab foo(amrex::coarsen(fmf.boxArray(),ratio), fmf.DistributionMap(),
658 1, 0, MFInfo().SetAlloc(false));
659 const FabArrayBase::CPC& cpc = mask.getCPC(cnghost,foo,IntVect::TheZeroVector(),period);
660 mask.setVal(fine_value, cpc, 0, 1);
661
662 has_cf = mask.RecvLayoutMask(cpc);
663
664 return mask;
665}
666
669template <typename FAB>
671 const IntVect& ratio, int ngcrse, bool mfiter_is_definitely_safe)
672{
673 AMREX_ASSERT(fine.is_nodal());
674 AMREX_ASSERT(crse.is_nodal());
675 AMREX_ASSERT(crse.nComp() == fine.nComp());
676
677 int ncomp = crse.nComp();
678 using value_type = typename FAB::value_type;
679
680 if (mfiter_is_definitely_safe || isMFIterSafe(fine, crse))
681 {
682#ifdef AMREX_USE_OMP
683#pragma omp parallel if (Gpu::notInLaunchRegion())
684#endif
685 for (MFIter mfi(crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
686 {
687 const Box& bx = mfi.growntilebox(ngcrse);
688 Array4<value_type> const& crsearr = crse.array(mfi);
689 Array4<value_type const> const& finearr = fine.const_array(mfi);
690
692 {
693 amrex_avgdown_nodes(tbx,crsearr,finearr,0,0,ncomp,ratio);
694 });
695 }
696 }
697 else
698 {
699 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
700 ncomp, ngcrse);
701 average_down_nodal(fine, ctmp, ratio, ngcrse);
702 crse.ParallelCopy(ctmp,0,0,ncomp,ngcrse,ngcrse);
703 }
704}
705
706// *************************************************************************************************************
707
708// Average fine cell-based MultiFab onto crse cell-centered MultiFab.
709// We do NOT assume that the coarse layout is a coarsened version of the fine layout.
710// This version does NOT use volume-weighting
711template<typename FAB>
712void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse, int scomp, int ncomp, int rr)
713{
714 average_down(S_fine,S_crse,scomp,ncomp,rr*IntVect::TheUnitVector());
715}
716
717template<typename FAB>
718void average_down (const FabArray<FAB>& S_fine, FabArray<FAB>& S_crse,
719 int scomp, int ncomp, const IntVect& ratio)
720{
721 BL_PROFILE("amrex::average_down");
722 AMREX_ASSERT(S_crse.nComp() == S_fine.nComp());
723 AMREX_ASSERT((S_crse.is_cell_centered() && S_fine.is_cell_centered()) ||
724 (S_crse.is_nodal() && S_fine.is_nodal()));
725
726 using value_type = typename FAB::value_type;
727
728 bool is_cell_centered = S_crse.is_cell_centered();
729
730 //
731 // Coarsen() the fine stuff on processors owning the fine data.
732 //
733 BoxArray crse_S_fine_BA = S_fine.boxArray(); crse_S_fine_BA.coarsen(ratio);
734
735 if (crse_S_fine_BA == S_crse.boxArray() && S_fine.DistributionMap() == S_crse.DistributionMap())
736 {
737#ifdef AMREX_USE_GPU
738 if (Gpu::inLaunchRegion() && S_crse.isFusingCandidate()) {
739 auto const& crsema = S_crse.arrays();
740 auto const& finema = S_fine.const_arrays();
741 if (is_cell_centered) {
742 ParallelFor(S_crse, IntVect(0), ncomp,
743 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
744 {
745 amrex_avgdown(i,j,k,n,crsema[box_no],finema[box_no],scomp,scomp,ratio);
746 });
747 } else {
748 ParallelFor(S_crse, IntVect(0), ncomp,
749 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
750 {
751 amrex_avgdown_nodes(i,j,k,n,crsema[box_no],finema[box_no],scomp,scomp,ratio);
752 });
753 }
754 if (!Gpu::inNoSyncRegion()) {
756 }
757 } else
758#endif
759 {
760#ifdef AMREX_USE_OMP
761#pragma omp parallel if (Gpu::notInLaunchRegion())
762#endif
763 for (MFIter mfi(S_crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
764 {
765 // NOTE: The tilebox is defined at the coarse level.
766 const Box& bx = mfi.tilebox();
767 Array4<value_type> const& crsearr = S_crse.array(mfi);
768 Array4<value_type const> const& finearr = S_fine.const_array(mfi);
769
770 if (is_cell_centered) {
771 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
772 {
773 amrex_avgdown(i,j,k,n,crsearr,finearr,scomp,scomp,ratio);
774 });
775 } else {
776 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
777 {
778 amrex_avgdown_nodes(i,j,k,n,crsearr,finearr,scomp,scomp,ratio);
779 });
780 }
781 }
782 }
783 }
784 else
785 {
786 FabArray<FAB> crse_S_fine(crse_S_fine_BA, S_fine.DistributionMap(), ncomp, 0, MFInfo(),DefaultFabFactory<FAB>());
787
788#ifdef AMREX_USE_GPU
789 if (Gpu::inLaunchRegion() && crse_S_fine.isFusingCandidate()) {
790 auto const& crsema = crse_S_fine.arrays();
791 auto const& finema = S_fine.const_arrays();
792 if (is_cell_centered) {
793 ParallelFor(crse_S_fine, IntVect(0), ncomp,
794 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
795 {
796 amrex_avgdown(i,j,k,n,crsema[box_no],finema[box_no],0,scomp,ratio);
797 });
798 } else {
799 ParallelFor(crse_S_fine, IntVect(0), ncomp,
800 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
801 {
802 amrex_avgdown_nodes(i,j,k,n,crsema[box_no],finema[box_no],0,scomp,ratio);
803 });
804 }
805 if (!Gpu::inNoSyncRegion()) {
807 }
808 } else
809#endif
810 {
811#ifdef AMREX_USE_OMP
812#pragma omp parallel if (Gpu::notInLaunchRegion())
813#endif
814 for (MFIter mfi(crse_S_fine,TilingIfNotGPU()); mfi.isValid(); ++mfi)
815 {
816 // NOTE: The tilebox is defined at the coarse level.
817 const Box& bx = mfi.tilebox();
818 Array4<value_type> const& crsearr = crse_S_fine.array(mfi);
819 Array4<value_type const> const& finearr = S_fine.const_array(mfi);
820
821 // NOTE: We copy from component scomp of the fine fab into component 0 of the crse fab
822 // because the crse fab is a temporary which was made starting at comp 0, it is
823 // not part of the actual crse multifab which came in.
824
825 if (is_cell_centered) {
826 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
827 {
828 amrex_avgdown(i,j,k,n,crsearr,finearr,0,scomp,ratio);
829 });
830 } else {
831 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
832 {
833 amrex_avgdown_nodes(i,j,k,n,crsearr,finearr,0,scomp,ratio);
834 });
835 }
836 }
837 }
838
839 S_crse.ParallelCopy(crse_S_fine,0,scomp,ncomp);
840 }
841}
842
843
844
845
846
854template <typename F>
855Real
856NormHelper (const MultiFab& x, int xcomp,
857 const MultiFab& y, int ycomp,
858 F const& f,
859 int numcomp, IntVect nghost, bool local)
860{
861 BL_ASSERT(x.boxArray() == y.boxArray());
862 BL_ASSERT(x.DistributionMap() == y.DistributionMap());
863 BL_ASSERT(x.nGrowVect().allGE(nghost) && y.nGrowVect().allGE(nghost));
864
865 Real sm = Real(0.0);
866#ifdef AMREX_USE_GPU
867 if (Gpu::inLaunchRegion()) {
868 auto const& xma = x.const_arrays();
869 auto const& yma = y.const_arrays();
871 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple<Real>
872 {
873 Real t = Real(0.0);
874 auto const& xfab = xma[box_no];
875 auto const& yfab = yma[box_no];
876 for (int n = 0; n < numcomp; ++n) {
877 t += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
878 }
879 return t;
880 });
881 } else
882#endif
883 {
884#ifdef AMREX_USE_OMP
885#pragma omp parallel if (!system::regtest_reduction) reduction(+:sm)
886#endif
887 for (MFIter mfi(x,true); mfi.isValid(); ++mfi)
888 {
889 Box const& bx = mfi.growntilebox(nghost);
890 Array4<Real const> const& xfab = x.const_array(mfi);
891 Array4<Real const> const& yfab = y.const_array(mfi);
892 AMREX_LOOP_4D(bx, numcomp, i, j, k, n,
893 {
894 sm += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
895 });
896 }
897 }
898
899 if (!local) {
901 }
902
903 return sm;
904}
905
914template <typename MMF, typename Pred, typename F>
915Real
916NormHelper (const MMF& mask,
917 const MultiFab& x, int xcomp,
918 const MultiFab& y, int ycomp,
919 Pred const& pf,
920 F const& f,
921 int numcomp, IntVect nghost, bool local)
922{
923 BL_ASSERT(x.boxArray() == y.boxArray());
924 BL_ASSERT(x.boxArray() == mask.boxArray());
925 BL_ASSERT(x.DistributionMap() == y.DistributionMap());
926 BL_ASSERT(x.DistributionMap() == mask.DistributionMap());
927 BL_ASSERT(x.nGrowVect().allGE(nghost) && y.nGrowVect().allGE(nghost));
928 BL_ASSERT(mask.nGrowVect().allGE(nghost));
929
930 Real sm = Real(0.0);
931#ifdef AMREX_USE_GPU
932 if (Gpu::inLaunchRegion()) {
933 auto const& xma = x.const_arrays();
934 auto const& yma = y.const_arrays();
935 auto const& mma = mask.const_arrays();
937 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept -> GpuTuple<Real>
938 {
939 Real t = Real(0.0);
940 if (pf(mma[box_no](i,j,k))) {
941 auto const& xfab = xma[box_no];
942 auto const& yfab = yma[box_no];
943 for (int n = 0; n < numcomp; ++n) {
944 t += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
945 }
946 }
947 return t;
948 });
949 } else
950#endif
951 {
952#ifdef AMREX_USE_OMP
953#pragma omp parallel if (!system::regtest_reduction) reduction(+:sm)
954#endif
955 for (MFIter mfi(x,true); mfi.isValid(); ++mfi)
956 {
957 Box const& bx = mfi.growntilebox(nghost);
958 Array4<Real const> const& xfab = x.const_array(mfi);
959 Array4<Real const> const& yfab = y.const_array(mfi);
960 auto const& mfab = mask.const_array(mfi);
961 AMREX_LOOP_4D(bx, numcomp, i, j, k, n,
962 {
963 if (pf(mfab(i,j,k))) {
964 sm += f(xfab(i,j,k,xcomp+n) , yfab(i,j,k,ycomp+n));
965 }
966 });
967 }
968 }
969
970 if (!local) {
972 }
973
974 return sm;
975}
976
977template <typename CMF, typename FMF,
978 std::enable_if_t<IsFabArray_v<CMF> && IsFabArray_v<FMF>, int> FOO>
979void average_face_to_cellcenter (CMF& cc, int dcomp,
981 int ngrow)
982{
983 IntVect ng_vect(ngrow);
984 average_face_to_cellcenter(cc, dcomp, fc, ng_vect);
985}
986
987template <typename CMF, typename FMF,
988 std::enable_if_t<IsFabArray_v<CMF> && IsFabArray_v<FMF>, int> FOO>
989void average_face_to_cellcenter (CMF& cc, int dcomp,
991 IntVect const& ng_vect)
992{
993 AMREX_ASSERT(cc.nComp() >= dcomp + AMREX_SPACEDIM);
994 AMREX_ASSERT(fc[0]->nComp() == 1);
995
996#ifdef AMREX_USE_GPU
997 if (Gpu::inLaunchRegion() && cc.isFusingCandidate()) {
998 auto const& ccma = cc.arrays();
999 AMREX_D_TERM(auto const& fxma = fc[0]->const_arrays();,
1000 auto const& fyma = fc[1]->const_arrays();,
1001 auto const& fzma = fc[2]->const_arrays(););
1002 ParallelFor(cc, ng_vect,
1003 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
1004 {
1005#if (AMREX_SPACEDIM == 1)
1006 GeometryData gd{};
1007 gd.coord = 0;
1008#endif
1009 amrex_avg_fc_to_cc(i,j,k, ccma[box_no], AMREX_D_DECL(fxma[box_no],
1010 fyma[box_no],
1011 fzma[box_no]),
1012 dcomp
1013#if (AMREX_SPACEDIM == 1)
1014 , gd
1015#endif
1016 );
1017 });
1018 if (!Gpu::inNoSyncRegion()) {
1020 }
1021 } else
1022#endif
1023 {
1024#ifdef AMREX_USE_OMP
1025#pragma omp parallel if (Gpu::notInLaunchRegion())
1026#endif
1027 for (MFIter mfi(cc,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1028 {
1029 const Box bx = mfi.growntilebox(ng_vect);
1030 auto const& ccarr = cc.array(mfi);
1031 AMREX_D_TERM(auto const& fxarr = fc[0]->const_array(mfi);,
1032 auto const& fyarr = fc[1]->const_array(mfi);,
1033 auto const& fzarr = fc[2]->const_array(mfi););
1034
1035#if (AMREX_SPACEDIM == 1)
1037 {
1038 GeometryData gd;
1039 gd.coord = 0;
1040 amrex_avg_fc_to_cc(i,j,k, ccarr, fxarr, dcomp, gd);
1041 });
1042#else
1044 {
1045 amrex_avg_fc_to_cc(i,j,k, ccarr, AMREX_D_DECL(fxarr,fyarr,fzarr), dcomp);
1046 });
1047#endif
1048 }
1049 }
1050}
1051
1052template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int>>
1054 const Vector<MF*>& crse,
1055 const IntVect& ratio, int ngcrse)
1056{
1057 AMREX_ASSERT(fine.size() == AMREX_SPACEDIM && crse.size() == AMREX_SPACEDIM);
1059 {{AMREX_D_DECL(fine[0],fine[1],fine[2])}},
1061 {{AMREX_D_DECL(crse[0],crse[1],crse[2])}},
1062 ratio, ngcrse);
1063}
1064
1065template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int>>
1067 const Vector<MF*>& crse, int ratio, int ngcrse)
1068{
1069 average_down_faces(fine,crse,IntVect{ratio},ngcrse);
1070}
1071
1072template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int>>
1075 int ratio, int ngcrse)
1076{
1077 average_down_faces(fine,crse,IntVect{ratio},ngcrse);
1078}
1079
1080template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int>>
1083 const IntVect& ratio, int ngcrse)
1084{
1085 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1086 average_down_faces(*fine[idim], *crse[idim], ratio, ngcrse);
1087 }
1088}
1089
1090template <typename FAB>
1092 const IntVect& ratio, int ngcrse)
1093{
1094 BL_PROFILE("average_down_faces");
1095
1096 AMREX_ASSERT(crse.nComp() == fine.nComp());
1097 AMREX_ASSERT(fine.ixType() == crse.ixType());
1098 const auto type = fine.ixType();
1099 int dir;
1100 for (dir = 0; dir < AMREX_SPACEDIM; ++dir) {
1101 if (type.nodeCentered(dir)) { break; }
1102 }
1103 auto tmptype = type;
1104 tmptype.unset(dir);
1105 if (dir >= AMREX_SPACEDIM || !tmptype.cellCentered()) {
1106 amrex::Abort("average_down_faces: not face index type");
1107 }
1108 const int ncomp = crse.nComp();
1109 if (isMFIterSafe(fine, crse))
1110 {
1111#ifdef AMREX_USE_GPU
1112 if (Gpu::inLaunchRegion() && crse.isFusingCandidate()) {
1113 auto const& crsema = crse.arrays();
1114 auto const& finema = fine.const_arrays();
1115 ParallelFor(crse, IntVect(ngcrse), ncomp,
1116 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
1117 {
1118 amrex_avgdown_faces(i,j,k,n, crsema[box_no], finema[box_no], 0, 0, ratio, dir);
1119 });
1120 if (!Gpu::inNoSyncRegion()) {
1122 }
1123 } else
1124#endif
1125 {
1126#ifdef AMREX_USE_OMP
1127#pragma omp parallel if (Gpu::notInLaunchRegion())
1128#endif
1129 for (MFIter mfi(crse,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1130 {
1131 const Box& bx = mfi.growntilebox(ngcrse);
1132 auto const& crsearr = crse.array(mfi);
1133 auto const& finearr = fine.const_array(mfi);
1134 AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
1135 {
1136 amrex_avgdown_faces(i,j,k,n, crsearr, finearr, 0, 0, ratio, dir);
1137 });
1138 }
1139 }
1140 }
1141 else
1142 {
1143 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
1144 ncomp, ngcrse, MFInfo(), DefaultFabFactory<FAB>());
1145 average_down_faces(fine, ctmp, ratio, ngcrse);
1146 crse.ParallelCopy(ctmp,0,0,ncomp,ngcrse,ngcrse);
1147 }
1148}
1149
1150template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int>>
1153 const IntVect& ratio, const Geometry& crse_geom)
1154{
1155 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1156 average_down_faces(*fine[idim], *crse[idim], ratio, crse_geom);
1157 }
1158}
1159
1160template <typename FAB>
1162 const IntVect& ratio, const Geometry& crse_geom)
1163{
1164 FabArray<FAB> ctmp(amrex::coarsen(fine.boxArray(),ratio), fine.DistributionMap(),
1165 crse.nComp(), 0);
1166 average_down_faces(fine, ctmp, ratio, 0);
1167 crse.ParallelCopy(ctmp,0,0,crse.nComp(),0,0,crse_geom.periodicity());
1168}
1169
1170template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> FOO>
1172{
1173 using T = typename MF::value_type;
1174 const int ncomp = mf.nComp();
1175 Gpu::DeviceVector<T> dv(ncomp);
1176 auto* dp = dv.data();
1177 bool found = false;
1178 auto loc = cell.dim3();
1179 for (MFIter mfi(mf); mfi.isValid() && !found; ++mfi)
1180 {
1181 Box const& box = mfi.validbox();
1182 if (box.contains(cell)) {
1183 found = true;
1184 auto const& fab = mf.const_array(mfi);
1185 amrex::ParallelFor(1, [=] AMREX_GPU_DEVICE (int) noexcept
1186 {
1187 for (int n = 0; n < ncomp; ++n) {
1188 dp[n] = fab(loc.x,loc.y,loc.z,n);
1189 }
1190 });
1191 }
1192 }
1193 Vector<T> hv;
1194 if (found) {
1195 hv.resize(ncomp);
1196 Gpu::copy(Gpu::deviceToHost, dv.begin(), dv.end(), hv.begin());
1197 }
1198 return hv;
1199}
1200
1201template <typename MF, std::enable_if_t<IsFabArray<MF>::value,int> FOO>
1202MF get_line_data (MF const& mf, int dir, IntVect const& cell, Box const& bnd_bx)
1203{
1204 bool do_bnd = (!bnd_bx.isEmpty());
1205
1206 BoxArray const& ba = mf.boxArray();
1207 DistributionMapping const& dm = mf.DistributionMap();
1208 const auto nboxes = static_cast<int>(ba.size());
1209
1210 BoxList bl(ba.ixType());
1211 Vector<int> procmap;
1212 Vector<int> index_map;
1213 if (!do_bnd) {
1214 for (int i = 0; i < nboxes; ++i) {
1215 Box const& b = ba[i];
1216 IntVect lo = cell;
1217 lo[dir] = b.smallEnd(dir);
1218 if (b.contains(lo)) {
1219 IntVect hi = lo;
1220 hi[dir] = b.bigEnd(dir);
1221 Box b1d(lo,hi,b.ixType());
1222 bl.push_back(b1d);
1223 procmap.push_back(dm[i]);
1224 index_map.push_back(i);
1225 }
1226 }
1227 } else {
1228 for (int i = 0; i < nboxes; ++i) {
1229 Box const& b = ba[i];
1230 Box const& b1d = bnd_bx & b;
1231 if (b1d.ok()) {
1232 bl.push_back(b1d);
1233 procmap.push_back(dm[i]);
1234 index_map.push_back(i);
1235 }
1236 }
1237 }
1238
1239 if (bl.isEmpty()) {
1240 return MF();
1241 } else {
1242 BoxArray rba(std::move(bl));
1243 DistributionMapping rdm(std::move(procmap));
1244 MF rmf(rba, rdm, mf.nComp(), IntVect(0),
1245 MFInfo().SetArena(mf.arena()));
1246#ifdef AMREX_USE_OMP
1247#pragma omp parallel if (Gpu::notInLaunchRegion())
1248#endif
1249 for (MFIter mfi(rmf); mfi.isValid(); ++mfi) {
1250 Box const& b = mfi.validbox();
1251 auto const& dfab = rmf.array(mfi);
1252 auto const& sfab = mf.const_array(index_map[mfi.index()]);
1253 amrex::ParallelFor(b, mf.nComp(),
1254 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
1255 {
1256 dfab(i,j,k,n) = sfab(i,j,k,n);
1257 });
1258 }
1259 return rmf;
1260 }
1261}
1262
1264namespace detail {
1265template <typename Op, typename T, typename F>
1266void reduce_to_plane (Array4<T> const& ar, int direction, Box const& bx, int box_no,
1267 F const& f)
1268{
1269#if defined(AMREX_USE_GPU)
1270 Box b2d = bx;
1271 b2d.setRange(direction,0);
1272 const auto blo = amrex::lbound(bx);
1273 const auto len = amrex::length(bx);
1274 constexpr int nthreads = 128;
1275 auto nblocks = static_cast<int>(b2d.numPts());
1276#ifdef AMREX_USE_SYCL
1277 constexpr std::size_t shared_mem_bytes = sizeof(T)*Gpu::Device::warp_size;
1278 amrex::launch<nthreads>(nblocks, shared_mem_bytes, Gpu::gpuStream(),
1279 [=] AMREX_GPU_DEVICE (Gpu::Handler const& h)
1280 {
1281 int bid = h.blockIdx();
1282 int tid = h.threadIdx();
1283#else
1284 amrex::launch<nthreads>(nblocks, Gpu::gpuStream(),
1285 [=] AMREX_GPU_DEVICE ()
1286 {
1287 int bid = blockIdx.x;
1288 int tid = threadIdx.x;
1289#endif
1290 T tmp;
1291 Op().init(tmp);
1292 T* p;
1293 if (direction == 0) {
1294 int k = bid / len.y;
1295 int j = bid - k*len.y;
1296 k += blo.z;
1297 j += blo.y;
1298 for (int i = blo.x + tid; i < blo.x+len.x; i += nthreads) {
1299 Op().local_update(tmp, f(box_no,i,j,k));
1300 }
1301 p = ar.ptr(0,j,k);
1302 } else if (direction == 1) {
1303 int k = bid / len.x;
1304 int i = bid - k*len.x;
1305 k += blo.z;
1306 i += blo.x;
1307 for (int j = blo.y + tid; j < blo.y+len.y; j += nthreads) {
1308 Op().local_update(tmp, f(box_no,i,j,k));
1309 }
1310 p = ar.ptr(i,0,k);
1311 } else {
1312 int j = bid / len.x;
1313 int i = bid - j*len.x;
1314 j += blo.y;
1315 i += blo.x;
1316 for (int k = blo.z + tid; k < blo.z+len.z; k += nthreads) {
1317 Op().local_update(tmp, f(box_no,i,j,k));
1318 }
1319 p = ar.ptr(i,j,0);
1320 }
1321#ifdef AMREX_USE_SYCL
1322 Op().template parallel_update<T>(*p, tmp, h);
1323#else
1324 Op().template parallel_update<T,nthreads>(*p, tmp);
1325#endif
1326 });
1327#else
1328 // CPU
1329 if (direction == 0) {
1330 AMREX_LOOP_3D(bx, i, j, k,
1331 {
1332 Op().local_update(ar(0,j,k), f(box_no,i,j,k));
1333 });
1334 } else if (direction == 1) {
1335 AMREX_LOOP_3D(bx, i, j, k,
1336 {
1337 Op().local_update(ar(i,0,k), f(box_no,i,j,k));
1338 });
1339 } else {
1340 AMREX_LOOP_3D(bx, i, j, k,
1341 {
1342 Op().local_update(ar(i,j,0), f(box_no,i,j,k));
1343 });
1344 }
1345#endif
1346}
1347}
1349
1350template <typename Op, typename T, typename FAB, typename F,
1351 std::enable_if_t<IsBaseFab<FAB>::value
1352#ifndef AMREX_USE_CUDA
1353 && IsCallableR<T,F,int,int,int,int>::value
1354#endif
1355 , int> FOO>
1356BaseFab<T>
1357ReduceToPlane (int direction, Box const& a_domain, FabArray<FAB> const& mf, F const& f)
1358{
1359 Box const domain = amrex::convert(a_domain, mf.ixType());
1360
1361 Box domain2d = domain;
1362 domain2d.setRange(direction, 0);
1363
1364 T initval;
1365 Op().init(initval);
1366
1367 BaseFab<T> r(domain2d);
1368 r.template setVal<RunOn::Device>(initval);
1369 auto const& ar = r.array();
1370
1371 for (MFIter mfi(mf,MFItInfo().UseDefaultStream().DisableDeviceSync());
1372 mfi.isValid(); ++mfi)
1373 {
1374 Box bx = mfi.validbox() & domain;
1375 if (bx.ok()) {
1376 int box_no = mfi.LocalIndex();
1377 detail::reduce_to_plane<Op, T>(ar, direction, bx, box_no, f);
1378 }
1379 }
1380 if (!Gpu::inNoSyncRegion()) {
1382 }
1383
1384 return r;
1385}
1386
1388namespace detail {
1389template <typename Op, typename FA, typename F>
1390FA reduce_to_plane (int direction, Box const& domain, FA const& mf, F const& f)
1391{
1392 using T = typename FA::value_type;
1393
1394 Box const ndomain = amrex::convert(domain, mf.ixType());
1395
1396 auto npts = amrex::convert(mf.boxArray(),IntVect(0)).numPts();
1397 if (npts != amrex::convert(domain, amrex::IntVect(0)).numPts()) {
1398 amrex::Abort("ReduceToPlaneMF: mf's BoxArray must have a rectangular domain.");
1399 }
1400
1401 T initval;
1402 Op().init(initval);
1403
1404 BoxList bl = mf.boxArray().boxList();
1405 for (auto& b : bl) {
1406 b.setRange(direction, 0);
1407 }
1408 BoxArray ba(std::move(bl));
1409 FA tmpfa(ba, mf.DistributionMap(), 1, 0);
1410 tmpfa.setVal(initval);
1411
1412 for (MFIter mfi(mf); mfi.isValid(); ++mfi)
1413 {
1414 Box bx = mfi.validbox() & ndomain;
1415 if (bx.ok()) {
1416 int box_no = mfi.LocalIndex();
1417 detail::reduce_to_plane<Op, T>(tmpfa.array(mfi), direction, bx, box_no, f);
1418 }
1419 }
1420
1421 return tmpfa;
1422}
1423
1424}
1426
1427template <typename Op, typename FA, typename F,
1428 std::enable_if_t<IsMultiFabLike_v<FA>
1429#ifndef AMREX_USE_CUDA
1430 && IsCallableR<typename FA::value_type,
1431 F,int,int,int,int>::value
1432#endif
1433 , int> FOO>
1434FA ReduceToPlaneMF (int direction, Box const& domain, FA const& mf, F const& f)
1435{
1436 auto [fa3, fa2] = ReduceToPlaneMF2<Op>(direction, domain, mf, f);
1437 fa3.ParallelCopy(fa2);
1438 return std::move(fa3);
1439}
1440
1441template <typename Op, typename FA, typename F,
1442 std::enable_if_t<IsMultiFabLike_v<FA>
1443#ifndef AMREX_USE_CUDA
1444 && IsCallableR<typename FA::value_type,
1445 F,int,int,int,int>::value
1446#endif
1447 , int> FOO>
1448std::pair<FA,FA>
1449ReduceToPlaneMF2 (int direction, Box const& domain, FA const& mf, F const& f)
1450{
1451 using T = typename FA::value_type;
1452
1453 T initval;
1454 Op().init(initval);
1455
1456 auto tmpmf = detail::reduce_to_plane<Op>(direction, domain, mf, f);
1457
1458 BoxList bl2d(mf.ixType());
1459 Vector<int> procmap2d;
1460 auto const& ba3d = mf.boxArray();
1461 auto const& dm3d = mf.DistributionMap();
1462 int dlo = domain.smallEnd(direction);
1463 for (int i = 0, N = mf.size(); i < N; ++i) {
1464 Box b = ba3d[i];
1465 if (b.smallEnd(direction) <= dlo && dlo <= b.bigEnd(direction)) {
1466 b.setRange(direction, 0);
1467 bl2d.push_back(b);
1468 procmap2d.push_back(dm3d[i]);
1469 }
1470 }
1471
1472 BoxArray ba2d(std::move(bl2d));
1473 DistributionMapping dm2d(std::move(procmap2d));
1474
1475 FA mf2d(ba2d, dm2d, 1, 0);
1476 mf2d.setVal(initval);
1477
1478 static_assert(std::is_same_v<Op, ReduceOpSum>, "Currently only ReduceOpSum is supported.");
1479 mf2d.ParallelAdd(tmpmf);
1480
1481 return std::make_pair(std::move(tmpmf), std::move(mf2d));
1482}
1483
1484template <typename Op, typename FA, typename F,
1485 std::enable_if_t<IsMultiFabLike_v<FA>
1486#ifndef AMREX_USE_CUDA
1487 && IsCallableR<typename FA::value_type,
1488 F,int,int,int,int>::value
1489#endif
1490 , int> FOO>
1491std::pair<FA,FA>
1492ReduceToPlaneMF2Patchy (int direction, Box const& domain, FA const& mf,
1493 IntVect const& plane_max_grid_size, F const& f)
1494{
1495 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mf.is_cell_centered(),
1496 "ReduceToPlaneMF2Patchy: only cell-centered data is supported.");
1497
1498 Box const ndomain = amrex::convert(domain, mf.ixType());
1499 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ndomain.contains(mf.boxArray().minimalBox()),
1500 "ReduceToPlaneMF2Patchy: subdomains are not supported; domain must cover the full extent of mf.boxArray().");
1501
1502 using T = typename FA::value_type;
1503 T initval;
1504 Op().init(initval);
1505
1506 BoxList bl_patch2d = mf.boxArray().boxList();
1507 for (auto& b : bl_patch2d) {
1508 b.setRange(direction, 0);
1509 }
1510 BoxArray ba_patch2d(std::move(bl_patch2d));
1511 FA tmpmf(ba_patch2d, mf.DistributionMap(), 1, 0);
1512 tmpmf.setVal(initval);
1513
1514 for (MFIter mfi(mf); mfi.isValid(); ++mfi)
1515 {
1516 Box bx = mfi.validbox() & ndomain;
1517 if (bx.ok()) {
1518 int box_no = mfi.LocalIndex();
1519 detail::reduce_to_plane<Op, T>(tmpmf.array(mfi), direction, bx, box_no, f);
1520 }
1521 }
1522
1523 BoxList bl_unique = ba_patch2d.boxList();
1524 bl_unique.simplify(); // reduces work for .removeOverlap()
1525 BoxArray ba2d(std::move(bl_unique));
1526 ba2d.removeOverlap();
1527
1528 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(plane_max_grid_size.allGE(1),
1529 "ReduceToPlaneMF2Patchy: plane_max_grid_size must be >= 1 in every direction.");
1530 IntVect mgs2d = plane_max_grid_size;
1531 mgs2d[direction] = 1;
1532 ba2d.maxSize(mgs2d);
1533
1534 DistributionMapping dm2d(ba2d);
1535
1536 FA mf2d(ba2d, dm2d, 1, 0);
1537 mf2d.setVal(initval);
1538
1539 static_assert(std::is_same_v<Op, ReduceOpSum>, "Currently only ReduceOpSum is supported.");
1540 mf2d.ParallelAdd(tmpmf);
1541
1542 return std::make_pair(std::move(tmpmf), std::move(mf2d));
1543}
1544
1545}
1546
1547#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:190
Array4< T const > array() const noexcept
Definition AMReX_BaseFab.H:382
A collection of Boxes stored in an Array.
Definition AMReX_BoxArray.H:568
IndexType ixType() const noexcept
Return index type of this BoxArray.
Definition AMReX_BoxArray.H:858
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:615
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:654
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:123
__host__ __device__ bool isEmpty() const noexcept
Checks if it is an empty BoxND.
Definition AMReX_Box.H:204
__host__ __device__ bool contains(const IntVectND< dim > &p) const noexcept
Return true if argument is contained within BoxND.
Definition AMReX_Box.H:212
__host__ __device__ IndexTypeND< dim > ixType() const noexcept
Return the indexing type.
Definition AMReX_Box.H:135
__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:1108
__host__ __device__ bool ok() const noexcept
Checks if it is a proper BoxND (including a valid type).
Definition AMReX_Box.H:208
__host__ __device__ const IntVectND< dim > & smallEnd() const &noexcept
Return the inclusive lower bound of the box.
Definition AMReX_Box.H:111
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:2704
bool is_cell_centered() const noexcept
This tests on whether the FabArray is cell-centered.
Definition AMReX_FabArrayBase.cpp:2698
bool is_nodal() const noexcept
This tests on whether the FabArray is fully nodal.
Definition AMReX_FabArrayBase.cpp:2686
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:131
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:350
Array4< typename FabArray< FAB >::value_type const > const_array(const MFIter &mfi) const noexcept
Definition AMReX_FabArray.H:592
void ParallelCopy(const FabArray< FAB > &src, const Periodicity &period=Periodicity::NonPeriodic(), CpOp op=FabArrayBase::COPY)
Definition AMReX_FabArray.H:854
Array4< typename FabArray< FAB >::value_type const > array(const MFIter &mfi) const noexcept
Definition AMReX_FabArray.H:568
MultiArray4< typename FabArray< FAB >::value_type > arrays() noexcept
Definition AMReX_FabArray.H:640
MultiArray4< typename FabArray< FAB >::value_type const > const_arrays() const noexcept
Definition AMReX_FabArray.H:654
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:74
Periodicity periodicity() const noexcept
Definition AMReX_Geometry.H:356
GPU-compatible tuple.
Definition AMReX_Tuple.H:98
static constexpr int warp_size
Definition AMReX_GpuDevice.H:236
__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
__host__ __device__ constexpr Dim3 dim3() const noexcept
Definition AMReX_IntVect.H:265
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:28
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:1345
__host__ __device__ Dim3 lbound(Array4< T > const &a) noexcept
Return the inclusive lower bounds of an Array4 in Dim3 form.
Definition AMReX_Array4.H:1317
__host__ __device__ BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) coarsening ratio.
Definition AMReX_Box.H:1409
std::array< T, N > Array
Definition AMReX_Array.H:26
void Sum(T &v, MPI_Comm comm)
Definition AMReX_ParallelReduce.H:221
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:92
bool inNoSyncRegion() noexcept
Definition AMReX_GpuControl.H:152
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:49
__host__ __device__ BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Return a BoxND with different type.
Definition AMReX_Box.H:1558
void FillRandomNormal(MultiFab &mf, int scomp, int ncomp, Real mean, Real stddev)
Fill MultiFab with random numbers from normal distribution.
Definition AMReX_MultiFabUtil.cpp:1260
int nComp(FabArrayBase const &fa)
Definition AMReX_FabArrayBase.cpp:2851
void FillRandom(MultiFab &mf, int scomp, int ncomp)
Fill MultiFab with random numbers from uniform distribution.
Definition AMReX_MultiFabUtil.cpp:1247
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
Definition AMReX_CTOParallelForImpl.H:193
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:358
__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:596
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:1171
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_face_to_cellcenter(MultiFab &cc, int dcomp, const Vector< const MultiFab * > &fc, IntVect const &ng_vect)
Definition AMReX_MultiFabUtil.cpp:155
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:651
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:1357
std::pair< FA, FA > ReduceToPlaneMF2(int direction, Box const &domain, FA const &mf, F const &f)
Reduce FabArray/MultiFab data to plane FabArray.
Definition AMReX_MultiFabUtil.H:1449
Vector< MultiFab > convexify(Vector< MultiFab const * > const &mf, Vector< IntVect > const &refinement_ratio)
Convexify AMR data.
Definition AMReX_MultiFabUtil.cpp:1273
FA ReduceToPlaneMF(int direction, Box const &domain, FA const &mf, F const &f)
Reduce FabArray/MultiFab data to plane FabArray.
Definition AMReX_MultiFabUtil.H:1434
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:1053
MultiFab periodicShift(MultiFab const &mf, IntVect const &offset, Periodicity const &period)
Periodic shift MultiFab.
Definition AMReX_MultiFabUtil.cpp:840
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:1006
BoxND< 3 > Box
Box is an alias for amrex::BoxND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:30
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:586
void writeFabs(const MultiFab &mf, const std::string &name)
Write each fab individually.
Definition AMReX_MultiFabUtil.cpp:573
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:856
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:739
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:491
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:864
IntVectND< 3 > IntVect
IntVect is an alias for amrex::IntVectND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:33
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:568
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:437
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:240
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:62
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:1160
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:1492
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:240
FabArray< BaseFab< Long > > ToLongMultiFab(const iMultiFab &imf)
Convert iMultiFab to Long.
Definition AMReX_MultiFabUtil.cpp:591
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:1202
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:791
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:104
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:670
A multidimensional array accessor.
Definition AMReX_Array4.H:283
parallel copy or add
Definition AMReX_FabArrayBase.H:538
Definition AMReX_Geometry.H:25
int coord
Definition AMReX_Geometry.H:59
FabArray memory allocation information.
Definition AMReX_FabArray.H:66
Definition AMReX_MFIter.H:20
Struct for holding types.
Definition AMReX_TypeList.H:13