Line data Source code
1 : #ifndef AMREX_ML_CELL_LINOP_H_
2 : #define AMREX_ML_CELL_LINOP_H_
3 : #include <AMReX_Config.H>
4 :
5 : #include <AMReX_MLLinOp.H>
6 : #include <AMReX_iMultiFab.H>
7 : #include <AMReX_MultiFabUtil.H>
8 : #include <AMReX_YAFluxRegister.H>
9 : #include <AMReX_MLLinOp_K.H>
10 : #include <AMReX_MLMG_K.H>
11 :
12 : #ifndef BL_NO_FORT
13 : #include <AMReX_MLLinOp_F.H>
14 : #endif
15 :
16 : namespace amrex {
17 :
18 : template <typename MF>
19 : class MLCellLinOpT // NOLINT(cppcoreguidelines-virtual-class-destructor)
20 : : public MLLinOpT<MF>
21 : {
22 : public:
23 :
24 : using FAB = typename MF::fab_type;
25 : using RT = typename MF::value_type;
26 :
27 : using BCType = LinOpBCType;
28 : using BCMode = typename MLLinOpT<MF>::BCMode;
29 : using StateMode = typename MLLinOpT<MF>::StateMode;
30 : using Location = typename MLLinOpT<MF>::Location;
31 :
32 : MLCellLinOpT ();
33 0 : ~MLCellLinOpT () override = default;
34 :
35 : MLCellLinOpT (const MLCellLinOpT<MF>&) = delete;
36 : MLCellLinOpT (MLCellLinOpT<MF>&&) = delete;
37 : MLCellLinOpT<MF>& operator= (const MLCellLinOpT<MF>&) = delete;
38 : MLCellLinOpT<MF>& operator= (MLCellLinOpT<MF>&&) = delete;
39 :
40 : void define (const Vector<Geometry>& a_geom,
41 : const Vector<BoxArray>& a_grids,
42 : const Vector<DistributionMapping>& a_dmap,
43 : const LPInfo& a_info = LPInfo(),
44 : const Vector<FabFactory<FAB> const*>& a_factory = {});
45 :
46 : void setLevelBC (int amrlev, const MF* levelbcdata,
47 : const MF* robinbc_a = nullptr,
48 : const MF* robinbc_b = nullptr,
49 : const MF* robinbc_f = nullptr) final;
50 :
51 : using MLLinOpT<MF>::setLevelBC;
52 :
53 : bool needsUpdate () const override {
54 : return MLLinOpT<MF>::needsUpdate();
55 : }
56 : void update () override;
57 :
58 : virtual bool isCrossStencil () const { return true; }
59 : virtual bool isTensorOp () const { return false; }
60 :
61 : void updateSolBC (int amrlev, const MF& crse_bcdata) const;
62 : void updateCorBC (int amrlev, const MF& crse_bcdata) const;
63 :
64 : virtual void applyBC (int amrlev, int mglev, MF& in, BCMode bc_mode, StateMode s_mode,
65 : const MLMGBndryT<MF>* bndry=nullptr, bool skip_fillboundary=false) const;
66 :
67 : BoxArray makeNGrids (int grid_size) const;
68 :
69 : void restriction (int, int, MF& crse, MF& fine) const override;
70 :
71 : void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const override;
72 :
73 : void interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const override;
74 :
75 : void interpolationAmr (int famrlev, MF& fine, const MF& crse,
76 : IntVect const& nghost) const override;
77 :
78 : void averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs,
79 : const MF& fine_sol, const MF& fine_rhs) override;
80 :
81 : void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode,
82 : StateMode s_mode, const MLMGBndryT<MF>* bndry=nullptr) const override;
83 : void smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
84 : bool skip_fillboundary=false) const final;
85 :
86 : void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b,
87 : const MF* crse_bcdata=nullptr) override;
88 :
89 : void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b,
90 : BCMode bc_mode, const MF* crse_bcdata=nullptr) final;
91 :
92 : // The assumption is crse_sol's boundary has been filled, but not fine_sol.
93 : void reflux (int crse_amrlev,
94 : MF& res, const MF& crse_sol, const MF&,
95 : MF&, MF& fine_sol, const MF&) const final;
96 : void compFlux (int amrlev, const Array<MF*,AMREX_SPACEDIM>& fluxes,
97 : MF& sol, Location loc) const override;
98 : void compGrad (int amrlev, const Array<MF*,AMREX_SPACEDIM>& grad,
99 : MF& sol, Location loc) const override;
100 :
101 : void applyMetricTerm (int amrlev, int mglev, MF& rhs) const final;
102 : void unapplyMetricTerm (int amrlev, int mglev, MF& rhs) const final;
103 : Vector<RT> getSolvabilityOffset (int amrlev, int mglev,
104 : MF const& rhs) const override;
105 : void fixSolvabilityByOffset (int amrlev, int mglev, MF& rhs,
106 : Vector<RT> const& offset) const override;
107 :
108 : void prepareForSolve () override;
109 :
110 : RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const final;
111 :
112 : virtual void Fapply (int amrlev, int mglev, MF& out, const MF& in) const = 0;
113 : virtual void Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const = 0;
114 : virtual void FFlux (int amrlev, const MFIter& mfi,
115 : const Array<FAB*,AMREX_SPACEDIM>& flux,
116 : const FAB& sol, Location loc, int face_only=0) const = 0;
117 :
118 : virtual void addInhomogNeumannFlux (int /*amrlev*/,
119 : const Array<MF*,AMREX_SPACEDIM>& /*grad*/,
120 : MF const& /*sol*/,
121 : bool /*mult_bcoef*/) const {}
122 :
123 : RT normInf (int amrlev, MF const& mf, bool local) const override;
124 :
125 : void averageDownAndSync (Vector<MF>& sol) const override;
126 :
127 : void avgDownResAmr (int clev, MF& cres, MF const& fres) const override;
128 :
129 : struct BCTL {
130 : BoundCond type;
131 : RT location;
132 : };
133 :
134 : Vector<std::unique_ptr<MF> > m_robin_bcval;
135 :
136 : protected:
137 :
138 : bool m_has_metric_term = false;
139 :
140 : Vector<std::unique_ptr<MLMGBndryT<MF>> > m_bndry_sol;
141 : Vector<std::unique_ptr<BndryRegisterT<MF>> > m_crse_sol_br;
142 :
143 : Vector<std::unique_ptr<MLMGBndryT<MF>> > m_bndry_cor;
144 : Vector<std::unique_ptr<BndryRegisterT<MF>> > m_crse_cor_br;
145 :
146 : // In case of agglomeration, coarse MG grids on amr level 0 are
147 : // not simply coarsened from fine MG grids. So we need to build
148 : // bcond and bcloc for each MG level.
149 : using RealTuple = Array<RT,2*BL_SPACEDIM>;
150 : using BCTuple = Array<BoundCond,2*BL_SPACEDIM>;
151 : class BndryCondLoc
152 : {
153 : public:
154 : BndryCondLoc (const BoxArray& ba, const DistributionMapping& dm, int ncomp);
155 :
156 : void setLOBndryConds (const Geometry& geom, const Real* dx,
157 : const Vector<Array<BCType,AMREX_SPACEDIM> >& lobc,
158 : const Vector<Array<BCType,AMREX_SPACEDIM> >& hibc,
159 : IntVect const& ratio, const RealVect& interior_bloc,
160 : const Array<Real,AMREX_SPACEDIM>& domain_bloc_lo,
161 : const Array<Real,AMREX_SPACEDIM>& domain_bloc_hi,
162 : LinOpBCType crse_fine_bc_type);
163 :
164 : const Vector<BCTuple>& bndryConds (const MFIter& mfi) const noexcept {
165 : return bcond[mfi];
166 : }
167 : const Vector<RealTuple>& bndryLocs (const MFIter& mfi) const noexcept {
168 : return bcloc[mfi];
169 : }
170 : const BCTuple& bndryConds (const MFIter& mfi, int icomp) const noexcept {
171 : return bcond[mfi][icomp];
172 : }
173 : const RealTuple& bndryLocs (const MFIter& mfi, int icomp) const noexcept {
174 : return bcloc[mfi][icomp];
175 : }
176 : GpuArray<BCTL,2*AMREX_SPACEDIM> const* getBCTLPtr (const MFIter& mfi) const noexcept {
177 : return bctl[mfi];
178 : }
179 : private:
180 : LayoutData<Vector<BCTuple> > bcond;
181 : LayoutData<Vector<RealTuple> > bcloc;
182 : LayoutData<GpuArray<BCTL,2*AMREX_SPACEDIM>*> bctl;
183 : Gpu::DeviceVector<GpuArray<BCTL,2*AMREX_SPACEDIM> > bctl_dv;
184 : int m_ncomp;
185 : };
186 : Vector<Vector<std::unique_ptr<BndryCondLoc> > > m_bcondloc;
187 :
188 : // used to save interpolation coefficients of the first interior cells
189 : mutable Vector<Vector<BndryRegisterT<MF>> > m_undrrelxr;
190 :
191 : // boundary cell flags for covered, not_covered, outside_domain
192 : Vector<Vector<Array<MultiMask,2*AMREX_SPACEDIM> > > m_maskvals;
193 :
194 : Vector<std::unique_ptr<iMultiFab> > m_norm_fine_mask;
195 :
196 : mutable Vector<YAFluxRegisterT<MF>> m_fluxreg;
197 :
198 : private:
199 :
200 : void defineAuxData ();
201 : void defineBC ();
202 :
203 : void computeVolInv () const;
204 : mutable Vector<Vector<RT> > m_volinv; // used by solvability fix
205 : };
206 :
207 : template <typename T>
208 : struct MLMGABCTag {
209 : Array4<T> fab;
210 : Array4<T const> bcval_lo;
211 : Array4<T const> bcval_hi;
212 : Array4<int const> mask_lo;
213 : Array4<int const> mask_hi;
214 : T bcloc_lo;
215 : T bcloc_hi;
216 : Box bx;
217 : BoundCond bctype_lo;
218 : BoundCond bctype_hi;
219 : int blen;
220 : int comp;
221 : int dir;
222 :
223 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
224 : Box const& box() const noexcept { return bx; }
225 : };
226 :
227 : template <typename T>
228 : struct MLMGPSTag {
229 : Array4<T> flo;
230 : Array4<T> fhi;
231 : Array4<int const> mlo;
232 : Array4<int const> mhi;
233 : T bcllo;
234 : T bclhi;
235 : Box bx;
236 : BoundCond bctlo;
237 : BoundCond bcthi;
238 : int blen;
239 : int comp;
240 : int dir;
241 :
242 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
243 : Box const& box() const noexcept { return bx; }
244 : };
245 :
246 : #ifdef AMREX_USE_EB
247 : template <typename T>
248 : struct MLMGPSEBTag {
249 : Array4<T> flo;
250 : Array4<T> fhi;
251 : Array4<T const> ap;
252 : Array4<int const> mlo;
253 : Array4<int const> mhi;
254 : T bcllo;
255 : T bclhi;
256 : Box bx;
257 : BoundCond bctlo;
258 : BoundCond bcthi;
259 : int blen;
260 : int comp;
261 : int dir;
262 :
263 : [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
264 : Box const& box() const noexcept { return bx; }
265 : };
266 : #endif
267 :
268 : template <typename MF>
269 : MLCellLinOpT<MF>::BndryCondLoc::BndryCondLoc (const BoxArray& ba,
270 : const DistributionMapping& dm,
271 : int ncomp)
272 : : bcond(ba, dm),
273 : bcloc(ba, dm),
274 : bctl(ba, dm),
275 : bctl_dv(bctl.local_size()*ncomp),
276 : m_ncomp(ncomp)
277 : {
278 : auto* dp = bctl_dv.data();
279 : for (MFIter mfi(bcloc); mfi.isValid(); ++mfi) {
280 : bcond[mfi].resize(ncomp);
281 : bcloc[mfi].resize(ncomp);
282 : bctl[mfi] = dp;
283 : dp += ncomp;
284 : }
285 : }
286 :
287 : template <typename MF>
288 : void
289 : MLCellLinOpT<MF>::BndryCondLoc::
290 : setLOBndryConds (const Geometry& geom, const Real* dx,
291 : const Vector<Array<BCType,AMREX_SPACEDIM> >& lobc,
292 : const Vector<Array<BCType,AMREX_SPACEDIM> >& hibc,
293 : IntVect const& ratio, const RealVect& interior_bloc,
294 : const Array<Real,AMREX_SPACEDIM>& domain_bloc_lo,
295 : const Array<Real,AMREX_SPACEDIM>& domain_bloc_hi,
296 : LinOpBCType crse_fine_bc_type)
297 : {
298 : const Box& domain = geom.Domain();
299 :
300 : #ifdef AMREX_USE_OMP
301 : #pragma omp parallel
302 : #endif
303 : for (MFIter mfi(bcloc); mfi.isValid(); ++mfi)
304 : {
305 : const Box& bx = mfi.validbox();
306 : for (int icomp = 0; icomp < m_ncomp; ++icomp) {
307 : RealTuple & bloc = bcloc[mfi][icomp];
308 : BCTuple & bctag = bcond[mfi][icomp];
309 : MLMGBndryT<MF>::setBoxBC(bloc, bctag, bx, domain,
310 : lobc[icomp], hibc[icomp],
311 : dx, ratio, interior_bloc,
312 : domain_bloc_lo, domain_bloc_hi,
313 : geom.isPeriodicArray(),
314 : crse_fine_bc_type);
315 : }
316 : }
317 :
318 : Gpu::PinnedVector<GpuArray<BCTL,2*AMREX_SPACEDIM> > hv;
319 : hv.reserve(bctl_dv.size());
320 : for (MFIter mfi(bctl); mfi.isValid(); ++mfi)
321 : {
322 : for (int icomp = 0; icomp < m_ncomp; ++icomp) {
323 : GpuArray<BCTL,2*AMREX_SPACEDIM> tmp;
324 : for (int m = 0; m < 2*AMREX_SPACEDIM; ++m) {
325 : tmp[m].type = bcond[mfi][icomp][m];
326 : tmp[m].location = bcloc[mfi][icomp][m];
327 : }
328 : hv.push_back(std::move(tmp));
329 : }
330 : }
331 : Gpu::copyAsync(Gpu::hostToDevice, hv.begin(), hv.end(), bctl_dv.begin());
332 : Gpu::streamSynchronize();
333 : }
334 :
335 : template <typename MF>
336 : MLCellLinOpT<MF>::MLCellLinOpT ()
337 : {
338 : this->m_ixtype = IntVect::TheCellVector();
339 : }
340 :
341 : template <typename MF>
342 : void
343 : MLCellLinOpT<MF>::define (const Vector<Geometry>& a_geom,
344 : const Vector<BoxArray>& a_grids,
345 : const Vector<DistributionMapping>& a_dmap,
346 : const LPInfo& a_info,
347 : const Vector<FabFactory<FAB> const*>& a_factory)
348 : {
349 : MLLinOpT<MF>::define(a_geom, a_grids, a_dmap, a_info, a_factory);
350 : defineAuxData();
351 : defineBC();
352 : }
353 :
354 : template <typename MF>
355 : void
356 : MLCellLinOpT<MF>::defineAuxData ()
357 : {
358 : BL_PROFILE("MLCellLinOp::defineAuxData()");
359 :
360 : m_undrrelxr.resize(this->m_num_amr_levels);
361 : m_maskvals.resize(this->m_num_amr_levels);
362 : m_fluxreg.resize(this->m_num_amr_levels-1);
363 : m_norm_fine_mask.resize(this->m_num_amr_levels-1);
364 :
365 : const int ncomp = this->getNComp();
366 :
367 : for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
368 : {
369 : m_undrrelxr[amrlev].resize(this->m_num_mg_levels[amrlev]);
370 : for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
371 : {
372 : m_undrrelxr[amrlev][mglev].define(this->m_grids[amrlev][mglev],
373 : this->m_dmap[amrlev][mglev],
374 : 1, 0, 0, ncomp);
375 : }
376 : }
377 :
378 : for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
379 : {
380 : m_maskvals[amrlev].resize(this->m_num_mg_levels[amrlev]);
381 : for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
382 : {
383 : for (OrientationIter oitr; oitr; ++oitr)
384 : {
385 : const Orientation face = oitr();
386 : const int ngrow = 1;
387 : const int extent = this->isCrossStencil() ? 0 : 1; // extend to corners
388 : m_maskvals[amrlev][mglev][face].define(this->m_grids[amrlev][mglev],
389 : this->m_dmap[amrlev][mglev],
390 : this->m_geom[amrlev][mglev],
391 : face, 0, ngrow, extent, 1, true);
392 : }
393 : }
394 : }
395 :
396 : for (int amrlev = 0; amrlev < this->m_num_amr_levels-1; ++amrlev)
397 : {
398 : const IntVect ratio{this->m_amr_ref_ratio[amrlev]};
399 : m_fluxreg[amrlev].define(this->m_grids[amrlev+1][0],
400 : this->m_grids[amrlev][0],
401 : this->m_dmap[amrlev+1][0],
402 : this->m_dmap[amrlev][0],
403 : this->m_geom[amrlev+1][0],
404 : this->m_geom[amrlev][0],
405 : ratio, amrlev+1, ncomp);
406 : m_norm_fine_mask[amrlev] = std::make_unique<iMultiFab>
407 : (makeFineMask(this->m_grids[amrlev][0], this->m_dmap[amrlev][0],
408 : this->m_grids[amrlev+1][0],
409 : ratio, 1, 0));
410 : }
411 :
412 : #if (AMREX_SPACEDIM != 3)
413 : m_has_metric_term = !this->m_geom[0][0].IsCartesian() && this->info.has_metric_term;
414 : #endif
415 : }
416 :
417 : template <typename MF>
418 : void
419 : MLCellLinOpT<MF>::defineBC ()
420 : {
421 : BL_PROFILE("MLCellLinOp::defineBC()");
422 :
423 : const int ncomp = this->getNComp();
424 :
425 : m_bndry_sol.resize(this->m_num_amr_levels);
426 : m_crse_sol_br.resize(this->m_num_amr_levels);
427 :
428 : m_bndry_cor.resize(this->m_num_amr_levels);
429 : m_crse_cor_br.resize(this->m_num_amr_levels);
430 :
431 : m_robin_bcval.resize(this->m_num_amr_levels);
432 :
433 : for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
434 : {
435 : m_bndry_sol[amrlev] = std::make_unique<MLMGBndryT<MF>>(this->m_grids[amrlev][0],
436 : this->m_dmap[amrlev][0],
437 : ncomp,
438 : this->m_geom[amrlev][0]);
439 : }
440 :
441 : for (int amrlev = 1; amrlev < this->m_num_amr_levels; ++amrlev)
442 : {
443 : const int in_rad = 0;
444 : const int out_rad = 1;
445 : const int extent_rad = 2;
446 : const int crse_ratio = this->m_amr_ref_ratio[amrlev-1];
447 : BoxArray cba = this->m_grids[amrlev][0];
448 : cba.coarsen(crse_ratio);
449 : m_crse_sol_br[amrlev] = std::make_unique<BndryRegisterT<MF>>
450 : (cba, this->m_dmap[amrlev][0], in_rad, out_rad, extent_rad, ncomp);
451 : }
452 :
453 : for (int amrlev = 1; amrlev < this->m_num_amr_levels; ++amrlev)
454 : {
455 : const int in_rad = 0;
456 : const int out_rad = 1;
457 : const int extent_rad = 2;
458 : const int crse_ratio = this->m_amr_ref_ratio[amrlev-1];
459 : BoxArray cba = this->m_grids[amrlev][0];
460 : cba.coarsen(crse_ratio);
461 : m_crse_cor_br[amrlev] = std::make_unique<BndryRegisterT<MF>>
462 : (cba, this->m_dmap[amrlev][0], in_rad, out_rad, extent_rad, ncomp);
463 : m_crse_cor_br[amrlev]->setVal(RT(0.0));
464 : }
465 :
466 : // This has be to done after m_crse_cor_br is defined.
467 : for (int amrlev = 1; amrlev < this->m_num_amr_levels; ++amrlev)
468 : {
469 : m_bndry_cor[amrlev] = std::make_unique<MLMGBndryT<MF>>
470 : (this->m_grids[amrlev][0], this->m_dmap[amrlev][0], ncomp, this->m_geom[amrlev][0]);
471 : MF bc_data(this->m_grids[amrlev][0], this->m_dmap[amrlev][0], ncomp, 1);
472 : bc_data.setVal(0.0);
473 :
474 : m_bndry_cor[amrlev]->setBndryValues(*m_crse_cor_br[amrlev], 0, bc_data, 0, 0, ncomp,
475 : IntVect(this->m_amr_ref_ratio[amrlev-1]));
476 :
477 : Vector<Array<LinOpBCType,AMREX_SPACEDIM> > bclohi
478 : (ncomp,Array<LinOpBCType,AMREX_SPACEDIM>{{AMREX_D_DECL(BCType::Dirichlet,
479 : BCType::Dirichlet,
480 : BCType::Dirichlet)}});
481 : m_bndry_cor[amrlev]->setLOBndryConds(bclohi, bclohi, this->m_amr_ref_ratio[amrlev-1], RealVect{});
482 : }
483 :
484 : m_bcondloc.resize(this->m_num_amr_levels);
485 : for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
486 : {
487 : m_bcondloc[amrlev].resize(this->m_num_mg_levels[amrlev]);
488 : for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
489 : {
490 : m_bcondloc[amrlev][mglev] = std::make_unique<BndryCondLoc>(this->m_grids[amrlev][mglev],
491 : this->m_dmap[amrlev][mglev],
492 : ncomp);
493 : }
494 : }
495 : }
496 :
497 : template <typename MF>
498 : void
499 : MLCellLinOpT<MF>::setLevelBC (int amrlev, const MF* a_levelbcdata, const MF* robinbc_a,
500 : const MF* robinbc_b, const MF* robinbc_f)
501 : {
502 : BL_PROFILE("MLCellLinOp::setLevelBC()");
503 :
504 : AMREX_ALWAYS_ASSERT(amrlev >= 0 && amrlev < this->m_num_amr_levels);
505 :
506 : const int ncomp = this->getNComp();
507 :
508 : MF zero;
509 : IntVect ng(1);
510 : if (this->hasHiddenDimension()) { ng[this->hiddenDirection()] = 0; }
511 : if (a_levelbcdata == nullptr) {
512 : zero.define(this->m_grids[amrlev][0], this->m_dmap[amrlev][0], ncomp, ng);
513 : zero.setVal(RT(0.0));
514 : } else {
515 : AMREX_ALWAYS_ASSERT(a_levelbcdata->nGrowVect().allGE(ng));
516 : }
517 : const MF& bcdata = (a_levelbcdata == nullptr) ? zero : *a_levelbcdata;
518 :
519 : IntVect br_ref_ratio(-1);
520 :
521 : if (amrlev == 0)
522 : {
523 : if (this->needsCoarseDataForBC())
524 : {
525 : AMREX_ALWAYS_ASSERT(!this->hasHiddenDimension());
526 : br_ref_ratio = this->m_coarse_data_crse_ratio.allGT(0) ? this->m_coarse_data_crse_ratio : IntVect(2);
527 : if (m_crse_sol_br[amrlev] == nullptr && br_ref_ratio.allGT(0))
528 : {
529 : const int in_rad = 0;
530 : const int out_rad = 1;
531 : const int extent_rad = 2;
532 : const IntVect crse_ratio = br_ref_ratio;
533 : BoxArray cba = this->m_grids[amrlev][0];
534 : cba.coarsen(crse_ratio);
535 : m_crse_sol_br[amrlev] = std::make_unique<BndryRegisterT<MF>>
536 : (cba, this->m_dmap[amrlev][0], in_rad, out_rad, extent_rad, ncomp);
537 : }
538 : if (this->m_coarse_data_for_bc != nullptr) {
539 : AMREX_ALWAYS_ASSERT(this->m_coarse_data_crse_ratio.allGT(0));
540 : const Box& cbx = amrex::coarsen(this->m_geom[0][0].Domain(), this->m_coarse_data_crse_ratio);
541 : m_crse_sol_br[amrlev]->copyFrom(*this->m_coarse_data_for_bc, 0, 0, 0, ncomp,
542 : this->m_geom[0][0].periodicity(cbx));
543 : } else {
544 : m_crse_sol_br[amrlev]->setVal(RT(0.0));
545 : }
546 : m_bndry_sol[amrlev]->setBndryValues(*m_crse_sol_br[amrlev], 0,
547 : bcdata, 0, 0, ncomp, br_ref_ratio);
548 : br_ref_ratio = this->m_coarse_data_crse_ratio;
549 : }
550 : else
551 : {
552 : m_bndry_sol[amrlev]->setPhysBndryValues(bcdata,0,0,ncomp);
553 : br_ref_ratio = IntVect(1);
554 : }
555 : }
556 : else
557 : {
558 : m_bndry_sol[amrlev]->setPhysBndryValues(bcdata,0,0,ncomp);
559 : br_ref_ratio = IntVect(this->m_amr_ref_ratio[amrlev-1]);
560 : }
561 :
562 : auto crse_fine_bc_type = (amrlev == 0) ? this->m_coarse_fine_bc_type : LinOpBCType::Dirichlet;
563 : m_bndry_sol[amrlev]->setLOBndryConds(this->m_lobc, this->m_hibc, br_ref_ratio,
564 : this->m_coarse_bc_loc, crse_fine_bc_type);
565 :
566 : const Real* dx = this->m_geom[amrlev][0].CellSize();
567 : for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
568 : {
569 : m_bcondloc[amrlev][mglev]->setLOBndryConds(this->m_geom[amrlev][mglev], dx,
570 : this->m_lobc, this->m_hibc,
571 : br_ref_ratio, this->m_coarse_bc_loc,
572 : this->m_domain_bloc_lo, this->m_domain_bloc_hi,
573 : crse_fine_bc_type);
574 : }
575 :
576 : if (this->hasRobinBC()) {
577 : AMREX_ASSERT(robinbc_a != nullptr && robinbc_b != nullptr && robinbc_f != nullptr);
578 : m_robin_bcval[amrlev] = std::make_unique<MF>(this->m_grids[amrlev][0],
579 : this->m_dmap[amrlev][0],
580 : ncomp*3, 1);
581 : const Box& domain = this->m_geom[amrlev][0].Domain();
582 : MFItInfo mfi_info;
583 : if (Gpu::notInLaunchRegion()) { mfi_info.SetDynamic(true); }
584 : #ifdef AMREX_USE_OMP
585 : #pragma omp parallel if (Gpu::notInLaunchRegion())
586 : #endif
587 : for (MFIter mfi(*m_robin_bcval[amrlev], mfi_info); mfi.isValid(); ++mfi) {
588 : Box const& vbx = mfi.validbox();
589 : Array4<RT const> const& ra = robinbc_a->const_array(mfi);
590 : Array4<RT const> const& rb = robinbc_b->const_array(mfi);
591 : Array4<RT const> const& rf = robinbc_f->const_array(mfi);
592 : for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
593 : const Box& blo = amrex::adjCellLo(vbx, idim);
594 : const Box& bhi = amrex::adjCellHi(vbx, idim);
595 : bool outside_domain_lo = !(domain.contains(blo));
596 : bool outside_domain_hi = !(domain.contains(bhi));
597 : if ((!outside_domain_lo) && (!outside_domain_hi)) { continue; }
598 : for (int icomp = 0; icomp < ncomp; ++icomp) {
599 : Array4<RT> const& rbc = (*m_robin_bcval[amrlev])[mfi].array(icomp*3);
600 : if (this->m_lobc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_lo)
601 : {
602 : AMREX_HOST_DEVICE_PARALLEL_FOR_3D(blo, i, j, k,
603 : {
604 : rbc(i,j,k,0) = ra(i,j,k,icomp);
605 : rbc(i,j,k,1) = rb(i,j,k,icomp);
606 : rbc(i,j,k,2) = rf(i,j,k,icomp);
607 : });
608 : }
609 : if (this->m_hibc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_hi)
610 : {
611 : AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bhi, i, j, k,
612 : {
613 : rbc(i,j,k,0) = ra(i,j,k,icomp);
614 : rbc(i,j,k,1) = rb(i,j,k,icomp);
615 : rbc(i,j,k,2) = rf(i,j,k,icomp);
616 : });
617 : }
618 : }
619 : }
620 : }
621 : }
622 : }
623 :
624 : template <typename MF>
625 : void
626 : MLCellLinOpT<MF>::update ()
627 : {
628 : if (MLLinOpT<MF>::needsUpdate()) { MLLinOpT<MF>::update(); }
629 : }
630 :
631 : template <typename MF>
632 : void
633 : MLCellLinOpT<MF>::updateSolBC (int amrlev, const MF& crse_bcdata) const
634 : {
635 : BL_PROFILE("MLCellLinOp::updateSolBC()");
636 :
637 : AMREX_ALWAYS_ASSERT(amrlev > 0);
638 : const int ncomp = this->getNComp();
639 : m_crse_sol_br[amrlev]->copyFrom(crse_bcdata, 0, 0, 0, ncomp,
640 : this->m_geom[amrlev-1][0].periodicity());
641 : m_bndry_sol[amrlev]->updateBndryValues(*m_crse_sol_br[amrlev], 0, 0, ncomp,
642 : IntVect(this->m_amr_ref_ratio[amrlev-1]));
643 : }
644 :
645 : template <typename MF>
646 : void
647 : MLCellLinOpT<MF>::updateCorBC (int amrlev, const MF& crse_bcdata) const
648 : {
649 : BL_PROFILE("MLCellLinOp::updateCorBC()");
650 : AMREX_ALWAYS_ASSERT(amrlev > 0);
651 : const int ncomp = this->getNComp();
652 : m_crse_cor_br[amrlev]->copyFrom(crse_bcdata, 0, 0, 0, ncomp,
653 : this->m_geom[amrlev-1][0].periodicity());
654 : m_bndry_cor[amrlev]->updateBndryValues(*m_crse_cor_br[amrlev], 0, 0, ncomp,
655 : IntVect(this->m_amr_ref_ratio[amrlev-1]));
656 : }
657 :
658 : template <typename MF>
659 : void
660 : MLCellLinOpT<MF>::applyBC (int amrlev, int mglev, MF& in, BCMode bc_mode, StateMode,
661 : const MLMGBndryT<MF>* bndry, bool skip_fillboundary) const
662 : {
663 : BL_PROFILE("MLCellLinOp::applyBC()");
664 : // No coarsened boundary values, cannot apply inhomog at mglev>0.
665 : BL_ASSERT(mglev == 0 || bc_mode == BCMode::Homogeneous);
666 : BL_ASSERT(bndry != nullptr || bc_mode == BCMode::Homogeneous);
667 :
668 : const int ncomp = this->getNComp();
669 : const int cross = isCrossStencil();
670 : const int tensorop = isTensorOp();
671 : if (!skip_fillboundary) {
672 : in.FillBoundary(0, ncomp, this->m_geom[amrlev][mglev].periodicity(), cross);
673 : }
674 :
675 : int flagbc = bc_mode == BCMode::Inhomogeneous;
676 : const int imaxorder = this->maxorder;
677 :
678 : const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
679 : const RT dxi = static_cast<RT>(dxinv[0]);
680 : const RT dyi = (AMREX_SPACEDIM >= 2) ? static_cast<RT>(dxinv[1]) : RT(1.0);
681 : const RT dzi = (AMREX_SPACEDIM == 3) ? static_cast<RT>(dxinv[2]) : RT(1.0);
682 :
683 : const auto& maskvals = m_maskvals[amrlev][mglev];
684 : const auto& bcondloc = *m_bcondloc[amrlev][mglev];
685 :
686 : FAB foofab(Box::TheUnitBox(),ncomp);
687 : const auto& foo = foofab.const_array();
688 :
689 : MFItInfo mfi_info;
690 : if (Gpu::notInLaunchRegion()) { mfi_info.SetDynamic(true); }
691 :
692 : AMREX_ALWAYS_ASSERT_WITH_MESSAGE(cross || tensorop || Gpu::notInLaunchRegion(),
693 : "non-cross stencil not support for gpu");
694 :
695 : const int hidden_direction = this->hiddenDirection();
696 :
697 : #ifdef AMREX_USE_GPU
698 : if ((cross || tensorop) && Gpu::inLaunchRegion())
699 : {
700 : Vector<MLMGABCTag<RT>> tags;
701 : tags.reserve(in.local_size()*AMREX_SPACEDIM*ncomp);
702 :
703 : for (MFIter mfi(in); mfi.isValid(); ++mfi) {
704 : const Box& vbx = mfi.validbox();
705 : const auto& iofab = in.array(mfi);
706 : const auto & bdlv = bcondloc.bndryLocs(mfi);
707 : const auto & bdcv = bcondloc.bndryConds(mfi);
708 :
709 : for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
710 : if (idim != hidden_direction) {
711 : const Orientation olo(idim,Orientation::low);
712 : const Orientation ohi(idim,Orientation::high);
713 : const auto& bvlo = (bndry != nullptr) ?
714 : bndry->bndryValues(olo).const_array(mfi) : foo;
715 : const auto& bvhi = (bndry != nullptr) ?
716 : bndry->bndryValues(ohi).const_array(mfi) : foo;
717 : for (int icomp = 0; icomp < ncomp; ++icomp) {
718 : tags.emplace_back(MLMGABCTag<RT>{iofab, bvlo, bvhi,
719 : maskvals[olo].const_array(mfi),
720 : maskvals[ohi].const_array(mfi),
721 : bdlv[icomp][olo], bdlv[icomp][ohi],
722 : amrex::adjCell(vbx,olo),
723 : bdcv[icomp][olo], bdcv[icomp][ohi],
724 : vbx.length(idim), icomp, idim});
725 : }
726 : }
727 : }
728 : }
729 :
730 : ParallelFor(tags,
731 : [=] AMREX_GPU_DEVICE (int i, int j, int k, MLMGABCTag<RT> const& tag) noexcept
732 : {
733 : if (tag.dir == 0)
734 : {
735 : mllinop_apply_bc_x(0, i, j, k, tag.blen, tag.fab,
736 : tag.mask_lo, tag.bctype_lo, tag.bcloc_lo, tag.bcval_lo,
737 : imaxorder, dxi, flagbc, tag.comp);
738 : mllinop_apply_bc_x(1, i+tag.blen+1, j, k, tag.blen, tag.fab,
739 : tag.mask_hi, tag.bctype_hi, tag.bcloc_hi, tag.bcval_hi,
740 : imaxorder, dxi, flagbc, tag.comp);
741 : }
742 : #if (AMREX_SPACEDIM > 1)
743 : else
744 : #if (AMREX_SPACEDIM > 2)
745 : if (tag.dir == 1)
746 : #endif
747 : {
748 : mllinop_apply_bc_y(0, i, j, k, tag.blen, tag.fab,
749 : tag.mask_lo, tag.bctype_lo, tag.bcloc_lo, tag.bcval_lo,
750 : imaxorder, dyi, flagbc, tag.comp);
751 : mllinop_apply_bc_y(1, i, j+tag.blen+1, k, tag.blen, tag.fab,
752 : tag.mask_hi, tag.bctype_hi, tag.bcloc_hi, tag.bcval_hi,
753 : imaxorder, dyi, flagbc, tag.comp);
754 : }
755 : #if (AMREX_SPACEDIM > 2)
756 : else {
757 : mllinop_apply_bc_z(0, i, j, k, tag.blen, tag.fab,
758 : tag.mask_lo, tag.bctype_lo, tag.bcloc_lo, tag.bcval_lo,
759 : imaxorder, dzi, flagbc, tag.comp);
760 : mllinop_apply_bc_z(1, i, j, k+tag.blen+1, tag.blen, tag.fab,
761 : tag.mask_hi, tag.bctype_hi, tag.bcloc_hi, tag.bcval_hi,
762 : imaxorder, dzi, flagbc, tag.comp);
763 : }
764 : #endif
765 : #endif
766 : });
767 : } else
768 : #endif
769 : if (cross || tensorop)
770 : {
771 : #ifdef AMREX_USE_OMP
772 : #pragma omp parallel if (Gpu::notInLaunchRegion())
773 : #endif
774 : for (MFIter mfi(in, mfi_info); mfi.isValid(); ++mfi)
775 : {
776 : const Box& vbx = mfi.validbox();
777 : const auto& iofab = in.array(mfi);
778 :
779 : const auto & bdlv = bcondloc.bndryLocs(mfi);
780 : const auto & bdcv = bcondloc.bndryConds(mfi);
781 :
782 : for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
783 : {
784 : if (hidden_direction == idim) { continue; }
785 : const Orientation olo(idim,Orientation::low);
786 : const Orientation ohi(idim,Orientation::high);
787 : const Box blo = amrex::adjCellLo(vbx, idim);
788 : const Box bhi = amrex::adjCellHi(vbx, idim);
789 : const int blen = vbx.length(idim);
790 : const auto& mlo = maskvals[olo].array(mfi);
791 : const auto& mhi = maskvals[ohi].array(mfi);
792 : const auto& bvlo = (bndry != nullptr) ? bndry->bndryValues(olo).const_array(mfi) : foo;
793 : const auto& bvhi = (bndry != nullptr) ? bndry->bndryValues(ohi).const_array(mfi) : foo;
794 : for (int icomp = 0; icomp < ncomp; ++icomp) {
795 : const BoundCond bctlo = bdcv[icomp][olo];
796 : const BoundCond bcthi = bdcv[icomp][ohi];
797 : const RT bcllo = bdlv[icomp][olo];
798 : const RT bclhi = bdlv[icomp][ohi];
799 : if (idim == 0) {
800 : mllinop_apply_bc_x(0, blo, blen, iofab, mlo,
801 : bctlo, bcllo, bvlo,
802 : imaxorder, dxi, flagbc, icomp);
803 : mllinop_apply_bc_x(1, bhi, blen, iofab, mhi,
804 : bcthi, bclhi, bvhi,
805 : imaxorder, dxi, flagbc, icomp);
806 : } else if (idim == 1) {
807 : mllinop_apply_bc_y(0, blo, blen, iofab, mlo,
808 : bctlo, bcllo, bvlo,
809 : imaxorder, dyi, flagbc, icomp);
810 : mllinop_apply_bc_y(1, bhi, blen, iofab, mhi,
811 : bcthi, bclhi, bvhi,
812 : imaxorder, dyi, flagbc, icomp);
813 : } else {
814 : mllinop_apply_bc_z(0, blo, blen, iofab, mlo,
815 : bctlo, bcllo, bvlo,
816 : imaxorder, dzi, flagbc, icomp);
817 : mllinop_apply_bc_z(1, bhi, blen, iofab, mhi,
818 : bcthi, bclhi, bvhi,
819 : imaxorder, dzi, flagbc, icomp);
820 : }
821 : }
822 : }
823 : }
824 : }
825 : else
826 : {
827 : #ifdef BL_NO_FORT
828 : amrex::Abort("amrex_mllinop_apply_bc not available when BL_NO_FORT=TRUE");
829 : #else
830 : if constexpr (std::is_same_v<Real,RT>) {
831 : #ifdef AMREX_USE_OMP
832 : #pragma omp parallel
833 : #endif
834 : for (MFIter mfi(in, mfi_info); mfi.isValid(); ++mfi)
835 : {
836 : const Box& vbx = mfi.validbox();
837 :
838 : const auto & bdlv = bcondloc.bndryLocs(mfi);
839 : const auto & bdcv = bcondloc.bndryConds(mfi);
840 :
841 : const RealTuple & bdl = bdlv[0];
842 : const BCTuple & bdc = bdcv[0];
843 :
844 : for (OrientationIter oitr; oitr; ++oitr)
845 : {
846 : const Orientation ori = oitr();
847 :
848 : int cdr = ori;
849 : RT bcl = bdl[ori];
850 : int bct = bdc[ori];
851 :
852 : const auto& fsfab = (bndry != nullptr) ? bndry->bndryValues(ori)[mfi] : foofab;
853 :
854 : const Mask& m = maskvals[ori][mfi];
855 :
856 : amrex_mllinop_apply_bc(BL_TO_FORTRAN_BOX(vbx),
857 : BL_TO_FORTRAN_ANYD(in[mfi]),
858 : BL_TO_FORTRAN_ANYD(m),
859 : cdr, bct, bcl,
860 : BL_TO_FORTRAN_ANYD(fsfab),
861 : imaxorder, dxinv, flagbc, ncomp, cross);
862 : }
863 : }
864 : } else {
865 : amrex::Abort("Not supported");
866 : }
867 : #endif
868 : }
869 : }
870 :
871 : template <typename MF>
872 : BoxArray
873 : MLCellLinOpT<MF>::makeNGrids (int grid_size) const
874 : {
875 : const Box& dombx = this->m_geom[0].back().Domain();
876 :
877 : const BoxArray& old_ba = this->m_grids[0].back();
878 : const int N = old_ba.size();
879 : Vector<Box> bv;
880 : bv.reserve(N);
881 : for (int i = 0; i < N; ++i)
882 : {
883 : Box b = old_ba[i];
884 : b.coarsen(grid_size);
885 : b.refine(grid_size);
886 : IntVect sz = b.size();
887 : const IntVect nblks {AMREX_D_DECL(sz[0]/grid_size, sz[1]/grid_size, sz[2]/grid_size)};
888 :
889 : IntVect big = b.smallEnd() + grid_size - 1;
890 : b.setBig(big);
891 :
892 : #if (AMREX_SPACEDIM == 3)
893 : for (int kk = 0; kk < nblks[2]; ++kk) {
894 : #endif
895 : #if (AMREX_SPACEDIM >= 2)
896 : for (int jj = 0; jj < nblks[1]; ++jj) {
897 : #endif
898 : for (int ii = 0; ii < nblks[0]; ++ii)
899 : {
900 : IntVect shft{AMREX_D_DECL(ii*grid_size,jj*grid_size,kk*grid_size)};
901 : Box bb = amrex::shift(b,shft);
902 : bb &= dombx;
903 : bv.push_back(bb);
904 : }
905 : #if (AMREX_SPACEDIM >= 2)
906 : }
907 : #endif
908 : #if (AMREX_SPACEDIM == 3)
909 : }
910 : #endif
911 : }
912 :
913 : std::sort(bv.begin(), bv.end());
914 : bv.erase(std::unique(bv.begin(), bv.end()), bv.end());
915 :
916 : BoxList bl(std::move(bv));
917 :
918 : return BoxArray{std::move(bl)};
919 : }
920 :
921 : template <typename MF>
922 : void
923 : MLCellLinOpT<MF>::restriction (int amrlev, int cmglev, MF& crse, MF& fine) const
924 : {
925 : const int ncomp = this->getNComp();
926 : IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[cmglev-1];
927 : amrex::average_down(fine, crse, 0, ncomp, ratio);
928 : }
929 :
930 : template <typename MF>
931 : void
932 : MLCellLinOpT<MF>::interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const
933 : {
934 : const int ncomp = this->getNComp();
935 :
936 : Dim3 ratio3 = {2,2,2};
937 : IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[fmglev];
938 : AMREX_D_TERM(ratio3.x = ratio[0];,
939 : ratio3.y = ratio[1];,
940 : ratio3.z = ratio[2];);
941 :
942 : #ifdef AMREX_USE_GPU
943 : if (Gpu::inLaunchRegion() && fine.isFusingCandidate()) {
944 : auto const& finema = fine.arrays();
945 : auto const& crsema = crse.const_arrays();
946 : ParallelFor(fine, IntVect(0), ncomp,
947 : [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
948 : {
949 : int ic = amrex::coarsen(i,ratio3.x);
950 : int jc = amrex::coarsen(j,ratio3.y);
951 : int kc = amrex::coarsen(k,ratio3.z);
952 : finema[box_no](i,j,k,n) += crsema[box_no](ic,jc,kc,n);
953 : });
954 : Gpu::streamSynchronize();
955 : } else
956 : #endif
957 : {
958 : #ifdef AMREX_USE_OMP
959 : #pragma omp parallel if (Gpu::notInLaunchRegion())
960 : #endif
961 : for (MFIter mfi(fine,TilingIfNotGPU()); mfi.isValid(); ++mfi)
962 : {
963 : const Box& bx = mfi.tilebox();
964 : Array4<RT const> const& cfab = crse.const_array(mfi);
965 : Array4<RT> const& ffab = fine.array(mfi);
966 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
967 : {
968 : int ic = amrex::coarsen(i,ratio3.x);
969 : int jc = amrex::coarsen(j,ratio3.y);
970 : int kc = amrex::coarsen(k,ratio3.z);
971 : ffab(i,j,k,n) += cfab(ic,jc,kc,n);
972 : });
973 : }
974 : }
975 : }
976 :
977 : template <typename MF>
978 : void
979 : MLCellLinOpT<MF>::interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const
980 : {
981 : const int ncomp = this->getNComp();
982 :
983 : const Geometry& crse_geom = this->Geom(amrlev,fmglev+1);
984 : const IntVect refratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[fmglev];
985 : const IntVect ng = crse.nGrowVect();
986 :
987 : MF cfine;
988 : const MF* cmf;
989 :
990 : if (amrex::isMFIterSafe(crse, fine))
991 : {
992 : crse.FillBoundary(crse_geom.periodicity());
993 : cmf = &crse;
994 : }
995 : else
996 : {
997 : BoxArray cba = fine.boxArray();
998 : cba.coarsen(refratio);
999 : cfine.define(cba, fine.DistributionMap(), ncomp, ng);
1000 : cfine.setVal(RT(0.0));
1001 : cfine.ParallelCopy(crse, 0, 0, ncomp, IntVect(0), ng, crse_geom.periodicity());
1002 : cmf = & cfine;
1003 : }
1004 :
1005 : bool isEB = fine.hasEBFabFactory();
1006 : ignore_unused(isEB);
1007 :
1008 : #ifdef AMREX_USE_EB
1009 : const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(&(fine.Factory()));
1010 : const FabArray<EBCellFlagFab>* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr;
1011 : #endif
1012 :
1013 : MFItInfo mfi_info;
1014 : if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
1015 : #ifdef AMREX_USE_OMP
1016 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1017 : #endif
1018 : for (MFIter mfi(fine, mfi_info); mfi.isValid(); ++mfi)
1019 : {
1020 : const Box& bx = mfi.tilebox();
1021 : const auto& ff = fine.array(mfi);
1022 : const auto& cc = cmf->array(mfi);
1023 : #ifdef AMREX_USE_EB
1024 : bool call_lincc;
1025 : if (isEB)
1026 : {
1027 : const auto& flag = (*flags)[mfi];
1028 : if (flag.getType(amrex::grow(bx,1)) == FabType::regular) {
1029 : call_lincc = true;
1030 : } else {
1031 : Array4<EBCellFlag const> const& flg = flag.const_array();
1032 : AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx,
1033 : {
1034 : mlmg_eb_cc_interp_r<2>(tbx, ff, cc, flg, ncomp);
1035 : });
1036 :
1037 : call_lincc = false;
1038 : }
1039 : }
1040 : else
1041 : {
1042 : call_lincc = true;
1043 : }
1044 : #else
1045 : const bool call_lincc = true;
1046 : #endif
1047 : if (call_lincc)
1048 : {
1049 : #if (AMREX_SPACEDIM == 3)
1050 : if (this->hasHiddenDimension()) {
1051 : Box const& bx_2d = this->compactify(bx);
1052 : auto const& ff_2d = this->compactify(ff);
1053 : auto const& cc_2d = this->compactify(cc);
1054 : AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx_2d, tbx,
1055 : {
1056 : TwoD::mlmg_lin_cc_interp_r2(tbx, ff_2d, cc_2d, ncomp);
1057 : });
1058 : } else
1059 : #endif
1060 : {
1061 : AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx,
1062 : {
1063 : mlmg_lin_cc_interp_r2(tbx, ff, cc, ncomp);
1064 : });
1065 : }
1066 : }
1067 : }
1068 : }
1069 :
1070 : template <typename MF>
1071 : void
1072 : MLCellLinOpT<MF>::interpolationAmr (int famrlev, MF& fine, const MF& crse,
1073 : IntVect const& /*nghost*/) const
1074 : {
1075 : const int ncomp = this->getNComp();
1076 : const int refratio = this->AMRRefRatio(famrlev-1);
1077 :
1078 : #ifdef AMREX_USE_EB
1079 : const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(famrlev));
1080 : const FabArray<EBCellFlagFab>* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr;
1081 : #endif
1082 :
1083 : MFItInfo mfi_info;
1084 : if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
1085 : #ifdef AMREX_USE_OMP
1086 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1087 : #endif
1088 : for (MFIter mfi(fine, mfi_info); mfi.isValid(); ++mfi)
1089 : {
1090 : const Box& bx = mfi.tilebox();
1091 : auto const& ff = fine.array(mfi);
1092 : auto const& cc = crse.const_array(mfi);
1093 : #ifdef AMREX_USE_EB
1094 : bool call_lincc;
1095 : if (factory)
1096 : {
1097 : const auto& flag = (*flags)[mfi];
1098 : if (flag.getType(amrex::grow(bx,1)) == FabType::regular) {
1099 : call_lincc = true;
1100 : } else {
1101 : Array4<EBCellFlag const> const& flg = flag.const_array();
1102 : switch(refratio) {
1103 : case 2:
1104 : {
1105 : AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx,
1106 : {
1107 : mlmg_eb_cc_interp_r<2>(tbx, ff, cc, flg, ncomp);
1108 : });
1109 : break;
1110 : }
1111 : case 4:
1112 : {
1113 : AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx,
1114 : {
1115 : mlmg_eb_cc_interp_r<4>(tbx, ff, cc, flg, ncomp);
1116 : });
1117 : break;
1118 : }
1119 : default:
1120 : amrex::Abort("mlmg_eb_cc_interp: only refratio 2 and 4 are supported");
1121 : }
1122 :
1123 : call_lincc = false;
1124 : }
1125 : }
1126 : else
1127 : {
1128 : call_lincc = true;
1129 : }
1130 : #else
1131 : const bool call_lincc = true;
1132 : #endif
1133 : if (call_lincc)
1134 : {
1135 : switch(refratio) {
1136 : case 2:
1137 : {
1138 : AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx,
1139 : {
1140 : mlmg_lin_cc_interp_r2(tbx, ff, cc, ncomp);
1141 : });
1142 : break;
1143 : }
1144 : case 4:
1145 : {
1146 : AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx,
1147 : {
1148 : mlmg_lin_cc_interp_r4(tbx, ff, cc, ncomp);
1149 : });
1150 : break;
1151 : }
1152 : default:
1153 : amrex::Abort("mlmg_lin_cc_interp: only refratio 2 and 4 are supported");
1154 : }
1155 : }
1156 : }
1157 : }
1158 :
1159 : template <typename MF>
1160 : void
1161 : MLCellLinOpT<MF>::averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs,
1162 : const MF& fine_sol, const MF& fine_rhs)
1163 : {
1164 : const auto amrrr = this->AMRRefRatio(camrlev);
1165 : const int ncomp = this->getNComp();
1166 : amrex::average_down(fine_sol, crse_sol, 0, ncomp, amrrr);
1167 : amrex::average_down(fine_rhs, crse_rhs, 0, ncomp, amrrr);
1168 : }
1169 :
1170 : template <typename MF>
1171 : void
1172 : MLCellLinOpT<MF>::apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode,
1173 : StateMode s_mode, const MLMGBndryT<MF>* bndry) const
1174 : {
1175 : BL_PROFILE("MLCellLinOp::apply()");
1176 : applyBC(amrlev, mglev, in, bc_mode, s_mode, bndry);
1177 : Fapply(amrlev, mglev, out, in);
1178 : }
1179 :
1180 : template <typename MF>
1181 : void
1182 : MLCellLinOpT<MF>::smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
1183 : bool skip_fillboundary) const
1184 : {
1185 : BL_PROFILE("MLCellLinOp::smooth()");
1186 : for (int redblack = 0; redblack < 2; ++redblack)
1187 : {
1188 : applyBC(amrlev, mglev, sol, BCMode::Homogeneous, StateMode::Solution,
1189 : nullptr, skip_fillboundary);
1190 : Fsmooth(amrlev, mglev, sol, rhs, redblack);
1191 : skip_fillboundary = false;
1192 : }
1193 : }
1194 :
1195 : template <typename MF>
1196 : void
1197 : MLCellLinOpT<MF>::solutionResidual (int amrlev, MF& resid, MF& x, const MF& b,
1198 : const MF* crse_bcdata)
1199 : {
1200 : BL_PROFILE("MLCellLinOp::solutionResidual()");
1201 : const int ncomp = this->getNComp();
1202 : if (crse_bcdata != nullptr) {
1203 : updateSolBC(amrlev, *crse_bcdata);
1204 : }
1205 : const int mglev = 0;
1206 : apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution,
1207 : m_bndry_sol[amrlev].get());
1208 :
1209 : AMREX_ASSERT(resid.nComp() == b.nComp());
1210 : MF::Xpay(resid, RT(-1.0), b, 0, 0, ncomp, IntVect(0));
1211 : }
1212 :
1213 : template <typename MF>
1214 : void
1215 : MLCellLinOpT<MF>::correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b,
1216 : BCMode bc_mode, const MF* crse_bcdata)
1217 : {
1218 : BL_PROFILE("MLCellLinOp::correctionResidual()");
1219 : const int ncomp = this->getNComp();
1220 : if (bc_mode == BCMode::Inhomogeneous)
1221 : {
1222 : if (crse_bcdata)
1223 : {
1224 : AMREX_ASSERT(mglev == 0 && amrlev > 0);
1225 : updateCorBC(amrlev, *crse_bcdata);
1226 : }
1227 : apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Correction,
1228 : m_bndry_cor[amrlev].get());
1229 : }
1230 : else
1231 : {
1232 : AMREX_ASSERT(crse_bcdata == nullptr);
1233 : apply(amrlev, mglev, resid, x, BCMode::Homogeneous, StateMode::Correction, nullptr);
1234 : }
1235 :
1236 : MF::Xpay(resid, Real(-1.0), b, 0, 0, ncomp, IntVect(0));
1237 : }
1238 :
1239 : template <typename MF>
1240 : void
1241 : MLCellLinOpT<MF>::reflux (int crse_amrlev, MF& res, const MF& crse_sol, const MF&,
1242 : MF&, MF& fine_sol, const MF&) const
1243 : {
1244 : BL_PROFILE("MLCellLinOp::reflux()");
1245 :
1246 : auto& fluxreg = m_fluxreg[crse_amrlev];
1247 : fluxreg.reset();
1248 :
1249 : const int ncomp = this->getNComp();
1250 :
1251 : const int fine_amrlev = crse_amrlev+1;
1252 :
1253 : Real dt = Real(1.0);
1254 : const Real* crse_dx = this->m_geom[crse_amrlev][0].CellSize();
1255 : const Real* fine_dx = this->m_geom[fine_amrlev][0].CellSize();
1256 :
1257 : const int mglev = 0;
1258 : applyBC(fine_amrlev, mglev, fine_sol, BCMode::Inhomogeneous, StateMode::Solution,
1259 : m_bndry_sol[fine_amrlev].get());
1260 :
1261 : MFItInfo mfi_info;
1262 : if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
1263 :
1264 : #ifdef AMREX_USE_OMP
1265 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1266 : #endif
1267 : {
1268 : Array<FAB,AMREX_SPACEDIM> flux;
1269 : Array<FAB*,AMREX_SPACEDIM> pflux {{ AMREX_D_DECL(flux.data(), flux.data()+1, flux.data()+2) }};
1270 : Array<FAB const*,AMREX_SPACEDIM> cpflux {{ AMREX_D_DECL(flux.data(), flux.data()+1, flux.data()+2) }};
1271 :
1272 : for (MFIter mfi(crse_sol, mfi_info); mfi.isValid(); ++mfi)
1273 : {
1274 : if (fluxreg.CrseHasWork(mfi))
1275 : {
1276 : const Box& tbx = mfi.tilebox();
1277 : AMREX_D_TERM(flux[0].resize(amrex::surroundingNodes(tbx,0),ncomp);,
1278 : flux[1].resize(amrex::surroundingNodes(tbx,1),ncomp);,
1279 : flux[2].resize(amrex::surroundingNodes(tbx,2),ncomp););
1280 : AMREX_D_TERM(Elixir elifx = flux[0].elixir();,
1281 : Elixir elify = flux[1].elixir();,
1282 : Elixir elifz = flux[2].elixir(););
1283 : FFlux(crse_amrlev, mfi, pflux, crse_sol[mfi], Location::FaceCentroid);
1284 : fluxreg.CrseAdd(mfi, cpflux, crse_dx, dt, RunOn::Gpu);
1285 : }
1286 : }
1287 :
1288 : #ifdef AMREX_USE_OMP
1289 : #pragma omp barrier
1290 : #endif
1291 :
1292 : for (MFIter mfi(fine_sol, mfi_info); mfi.isValid(); ++mfi)
1293 : {
1294 : if (fluxreg.FineHasWork(mfi))
1295 : {
1296 : const Box& tbx = mfi.tilebox();
1297 : const int face_only = true;
1298 : AMREX_D_TERM(flux[0].resize(amrex::surroundingNodes(tbx,0),ncomp);,
1299 : flux[1].resize(amrex::surroundingNodes(tbx,1),ncomp);,
1300 : flux[2].resize(amrex::surroundingNodes(tbx,2),ncomp););
1301 : AMREX_D_TERM(Elixir elifx = flux[0].elixir();,
1302 : Elixir elify = flux[1].elixir();,
1303 : Elixir elifz = flux[2].elixir(););
1304 : FFlux(fine_amrlev, mfi, pflux, fine_sol[mfi], Location::FaceCentroid, face_only);
1305 : fluxreg.FineAdd(mfi, cpflux, fine_dx, dt, RunOn::Gpu);
1306 : }
1307 : }
1308 : }
1309 :
1310 : fluxreg.Reflux(res);
1311 : }
1312 :
1313 : template <typename MF>
1314 : void
1315 : MLCellLinOpT<MF>::compFlux (int amrlev, const Array<MF*,AMREX_SPACEDIM>& fluxes,
1316 : MF& sol, Location loc) const
1317 : {
1318 : BL_PROFILE("MLCellLinOp::compFlux()");
1319 :
1320 : const int mglev = 0;
1321 : const int ncomp = this->getNComp();
1322 : applyBC(amrlev, mglev, sol, BCMode::Inhomogeneous, StateMode::Solution,
1323 : m_bndry_sol[amrlev].get());
1324 :
1325 : MFItInfo mfi_info;
1326 : if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
1327 :
1328 : #ifdef AMREX_USE_OMP
1329 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1330 : #endif
1331 : {
1332 : Array<FAB,AMREX_SPACEDIM> flux;
1333 : Array<FAB*,AMREX_SPACEDIM> pflux {{ AMREX_D_DECL(flux.data(), flux.data()+1, flux.data()+2) }};
1334 : for (MFIter mfi(sol, mfi_info); mfi.isValid(); ++mfi)
1335 : {
1336 : const Box& tbx = mfi.tilebox();
1337 : AMREX_D_TERM(flux[0].resize(amrex::surroundingNodes(tbx,0),ncomp);,
1338 : flux[1].resize(amrex::surroundingNodes(tbx,1),ncomp);,
1339 : flux[2].resize(amrex::surroundingNodes(tbx,2),ncomp););
1340 : AMREX_D_TERM(Elixir elifx = flux[0].elixir();,
1341 : Elixir elify = flux[1].elixir();,
1342 : Elixir elifz = flux[2].elixir(););
1343 : FFlux(amrlev, mfi, pflux, sol[mfi], loc);
1344 : for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1345 : const Box& nbx = mfi.nodaltilebox(idim);
1346 : auto const& dst = fluxes[idim]->array(mfi);
1347 : auto const& src = pflux[idim]->const_array();
1348 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D (nbx, ncomp, i, j, k, n,
1349 : {
1350 : dst(i,j,k,n) = src(i,j,k,n);
1351 : });
1352 : }
1353 : }
1354 : }
1355 : }
1356 :
1357 : template <typename MF>
1358 : void
1359 : MLCellLinOpT<MF>::compGrad (int amrlev, const Array<MF*,AMREX_SPACEDIM>& grad,
1360 : MF& sol, Location /*loc*/) const
1361 : {
1362 : BL_PROFILE("MLCellLinOp::compGrad()");
1363 :
1364 : if (sol.nComp() > 1) {
1365 : amrex::Abort("MLCellLinOp::compGrad called, but only works for single-component solves");
1366 : }
1367 :
1368 : const int mglev = 0;
1369 : applyBC(amrlev, mglev, sol, BCMode::Inhomogeneous, StateMode::Solution,
1370 : m_bndry_sol[amrlev].get());
1371 :
1372 : const int ncomp = this->getNComp();
1373 :
1374 : AMREX_D_TERM(const RT dxi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(0));,
1375 : const RT dyi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(1));,
1376 : const RT dzi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(2)););
1377 : #ifdef AMREX_USE_OMP
1378 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1379 : #endif
1380 : for (MFIter mfi(sol, TilingIfNotGPU()); mfi.isValid(); ++mfi)
1381 : {
1382 : AMREX_D_TERM(const Box& xbx = mfi.nodaltilebox(0);,
1383 : const Box& ybx = mfi.nodaltilebox(1);,
1384 : const Box& zbx = mfi.nodaltilebox(2););
1385 : const auto& s = sol.array(mfi);
1386 : AMREX_D_TERM(const auto& gx = grad[0]->array(mfi);,
1387 : const auto& gy = grad[1]->array(mfi);,
1388 : const auto& gz = grad[2]->array(mfi););
1389 :
1390 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( xbx, ncomp, i, j, k, n,
1391 : {
1392 : gx(i,j,k,n) = dxi*(s(i,j,k,n) - s(i-1,j,k,n));
1393 : });
1394 : #if (AMREX_SPACEDIM >= 2)
1395 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( ybx, ncomp, i, j, k, n,
1396 : {
1397 : gy(i,j,k,n) = dyi*(s(i,j,k,n) - s(i,j-1,k,n));
1398 : });
1399 : #endif
1400 : #if (AMREX_SPACEDIM == 3)
1401 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( zbx, ncomp, i, j, k, n,
1402 : {
1403 : gz(i,j,k,n) = dzi*(s(i,j,k,n) - s(i,j,k-1,n));
1404 : });
1405 : #endif
1406 : }
1407 :
1408 : addInhomogNeumannFlux(amrlev, grad, sol, false);
1409 : }
1410 :
1411 : template <typename MF>
1412 : void
1413 : MLCellLinOpT<MF>::applyMetricTerm (int amrlev, int mglev, MF& rhs) const
1414 : {
1415 : amrex::ignore_unused(amrlev,mglev,rhs);
1416 : #if (AMREX_SPACEDIM != 3)
1417 : if (!m_has_metric_term) { return; }
1418 :
1419 : const int ncomp = rhs.nComp();
1420 :
1421 : bool cc = rhs.ixType().cellCentered(0);
1422 :
1423 : const Geometry& geom = this->m_geom[amrlev][mglev];
1424 : const RT dx = static_cast<RT>(geom.CellSize(0));
1425 : const RT probxlo = static_cast<RT>(geom.ProbLo(0));
1426 :
1427 : #ifdef AMREX_USE_OMP
1428 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1429 : #endif
1430 : for (MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1431 : {
1432 : const Box& tbx = mfi.tilebox();
1433 : auto const& rhsarr = rhs.array(mfi);
1434 : #if (AMREX_SPACEDIM == 1)
1435 : if (cc) {
1436 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1437 : {
1438 : RT rc = probxlo + (i+RT(0.5))*dx;
1439 : rhsarr(i,j,k,n) *= rc*rc;
1440 : });
1441 : } else {
1442 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1443 : {
1444 : RT re = probxlo + i*dx;
1445 : rhsarr(i,j,k,n) *= re*re;
1446 : });
1447 : }
1448 : #elif (AMREX_SPACEDIM == 2)
1449 : if (cc) {
1450 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1451 : {
1452 : RT rc = probxlo + (i+RT(0.5))*dx;
1453 : rhsarr(i,j,k,n) *= rc;
1454 : });
1455 : } else {
1456 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1457 : {
1458 : RT re = probxlo + i*dx;
1459 : rhsarr(i,j,k,n) *= re;
1460 : });
1461 : }
1462 : #endif
1463 : }
1464 : #endif
1465 : }
1466 :
1467 : template <typename MF>
1468 : void
1469 : MLCellLinOpT<MF>::unapplyMetricTerm (int amrlev, int mglev, MF& rhs) const
1470 : {
1471 : amrex::ignore_unused(amrlev,mglev,rhs);
1472 : #if (AMREX_SPACEDIM != 3)
1473 : if (!m_has_metric_term) { return; }
1474 :
1475 : const int ncomp = rhs.nComp();
1476 :
1477 : bool cc = rhs.ixType().cellCentered(0);
1478 :
1479 : const Geometry& geom = this->m_geom[amrlev][mglev];
1480 : const RT dx = static_cast<RT>(geom.CellSize(0));
1481 : const RT probxlo = static_cast<RT>(geom.ProbLo(0));
1482 :
1483 : #ifdef AMREX_USE_OMP
1484 : #pragma omp parallel if (Gpu::notInLaunchRegion())
1485 : #endif
1486 : for (MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1487 : {
1488 : const Box& tbx = mfi.tilebox();
1489 : auto const& rhsarr = rhs.array(mfi);
1490 : #if (AMREX_SPACEDIM == 1)
1491 : if (cc) {
1492 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1493 : {
1494 : RT rcinv = RT(1.0)/(probxlo + (i+RT(0.5))*dx);
1495 : rhsarr(i,j,k,n) *= rcinv*rcinv;
1496 : });
1497 : } else {
1498 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1499 : {
1500 : RT re = probxlo + i*dx;
1501 : RT reinv = (re==RT(0.0)) ? RT(0.0) : RT(1.)/re;
1502 : rhsarr(i,j,k,n) *= reinv*reinv;
1503 : });
1504 : }
1505 : #elif (AMREX_SPACEDIM == 2)
1506 : if (cc) {
1507 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1508 : {
1509 : RT rcinv = RT(1.0)/(probxlo + (i+RT(0.5))*dx);
1510 : rhsarr(i,j,k,n) *= rcinv;
1511 : });
1512 : } else {
1513 : AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1514 : {
1515 : RT re = probxlo + i*dx;
1516 : RT reinv = (re==RT(0.0)) ? RT(0.0) : RT(1.)/re;
1517 : rhsarr(i,j,k,n) *= reinv;
1518 : });
1519 : }
1520 : #endif
1521 : }
1522 : #endif
1523 : }
1524 :
1525 : template <typename MF>
1526 : auto
1527 : MLCellLinOpT<MF>::getSolvabilityOffset (int amrlev, int mglev, MF const& rhs) const
1528 : -> Vector<RT>
1529 : {
1530 : computeVolInv();
1531 :
1532 : const int ncomp = this->getNComp();
1533 : Vector<RT> offset(ncomp);
1534 :
1535 : #ifdef AMREX_USE_EB
1536 : const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(amrlev,mglev));
1537 : if (factory && !factory->isAllRegular())
1538 : {
1539 : if constexpr (std::is_same<MF,MultiFab>()) {
1540 : const MultiFab& vfrac = factory->getVolFrac();
1541 : for (int c = 0; c < ncomp; ++c) {
1542 : offset[c] = amrex::Dot(rhs, c, vfrac, 0, 1, IntVect(0), true)
1543 : * m_volinv[amrlev][mglev];
1544 : }
1545 : } else {
1546 : amrex::Abort("TODO: MLMG with EB only works with MultiFab");
1547 : }
1548 : }
1549 : else
1550 : #endif
1551 : {
1552 : for (int c = 0; c < ncomp; ++c) {
1553 : offset[c] = rhs.sum(c,IntVect(0),true) * m_volinv[amrlev][mglev];
1554 : }
1555 : }
1556 :
1557 : ParallelAllReduce::Sum(offset.data(), ncomp, ParallelContext::CommunicatorSub());
1558 :
1559 : return offset;
1560 : }
1561 :
1562 : template <typename MF>
1563 : void
1564 : MLCellLinOpT<MF>::fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, MF& rhs,
1565 : Vector<RT> const& offset) const
1566 : {
1567 : const int ncomp = this->getNComp();
1568 : for (int c = 0; c < ncomp; ++c) {
1569 : rhs.plus(-offset[c], c, 1);
1570 : }
1571 : #ifdef AMREX_USE_EB
1572 : if (!rhs.isAllRegular()) {
1573 : if constexpr (std::is_same<MF,MultiFab>()) {
1574 : amrex::EB_set_covered(rhs, 0, ncomp, 0, 0.0_rt);
1575 : } else {
1576 : amrex::Abort("amrex::EB_set_covered only works with MultiFab");
1577 : }
1578 : }
1579 : #endif
1580 : }
1581 :
1582 : template <typename MF>
1583 : void
1584 : MLCellLinOpT<MF>::prepareForSolve ()
1585 : {
1586 : BL_PROFILE("MLCellLinOp::prepareForSolve()");
1587 :
1588 : const int imaxorder = this->maxorder;
1589 : const int ncomp = this->getNComp();
1590 : const int hidden_direction = this->hiddenDirection();
1591 : for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
1592 : {
1593 : for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
1594 : {
1595 : const auto& bcondloc = *m_bcondloc[amrlev][mglev];
1596 : const auto& maskvals = m_maskvals[amrlev][mglev];
1597 :
1598 : const RT dxi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(0));
1599 : const RT dyi = static_cast<RT>((AMREX_SPACEDIM >= 2) ? this->m_geom[amrlev][mglev].InvCellSize(1) : Real(1.0));
1600 : const RT dzi = static_cast<RT>((AMREX_SPACEDIM == 3) ? this->m_geom[amrlev][mglev].InvCellSize(2) : Real(1.0));
1601 :
1602 : auto& undrrelxr = this->m_undrrelxr[amrlev][mglev];
1603 : MF foo(this->m_grids[amrlev][mglev], this->m_dmap[amrlev][mglev], ncomp, 0, MFInfo().SetAlloc(false));
1604 :
1605 : #ifdef AMREX_USE_EB
1606 : const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->m_factory[amrlev][mglev].get());
1607 : const FabArray<EBCellFlagFab>* flags =
1608 : (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr;
1609 : auto area = (factory) ? factory->getAreaFrac()
1610 : : Array<const MultiCutFab*,AMREX_SPACEDIM>{AMREX_D_DECL(nullptr,nullptr,nullptr)};
1611 : amrex::ignore_unused(area);
1612 : #endif
1613 :
1614 : #ifdef AMREX_USE_GPU
1615 : if (Gpu::inLaunchRegion()) {
1616 : #ifdef AMREX_USE_EB
1617 : if (factory && !factory->isAllRegular()) {
1618 : if constexpr (!std::is_same<MF,MultiFab>()) {
1619 : amrex::Abort("MLCellLinOp with EB only works with MultiFab");
1620 : } else {
1621 : Vector<MLMGPSEBTag<RT>> tags;
1622 : tags.reserve(foo.local_size()*AMREX_SPACEDIM*ncomp);
1623 :
1624 : for (MFIter mfi(foo); mfi.isValid(); ++mfi)
1625 : {
1626 : const Box& vbx = mfi.validbox();
1627 :
1628 : const auto & bdlv = bcondloc.bndryLocs(mfi);
1629 : const auto & bdcv = bcondloc.bndryConds(mfi);
1630 :
1631 : auto fabtyp = (flags) ? (*flags)[mfi].getType(vbx) : FabType::regular;
1632 :
1633 : for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
1634 : {
1635 : if (idim != hidden_direction && fabtyp != FabType::covered) {
1636 : const Orientation olo(idim,Orientation::low);
1637 : const Orientation ohi(idim,Orientation::high);
1638 : auto const& ap = (fabtyp == FabType::singlevalued)
1639 : ? area[idim]->const_array(mfi) : Array4<Real const>{};
1640 : for (int icomp = 0; icomp < ncomp; ++icomp) {
1641 : tags.emplace_back(MLMGPSEBTag<RT>{undrrelxr[olo].array(mfi),
1642 : undrrelxr[ohi].array(mfi),
1643 : ap,
1644 : maskvals[olo].const_array(mfi),
1645 : maskvals[ohi].const_array(mfi),
1646 : bdlv[icomp][olo], bdlv[icomp][ohi],
1647 : amrex::adjCell(vbx,olo),
1648 : bdcv[icomp][olo], bdcv[icomp][ohi],
1649 : vbx.length(idim), icomp, idim});
1650 : }
1651 : }
1652 : }
1653 : }
1654 :
1655 : ParallelFor(tags,
1656 : [=] AMREX_GPU_DEVICE (int i, int j, int k, MLMGPSEBTag<RT> const& tag) noexcept
1657 : {
1658 : if (tag.ap) {
1659 : if (tag.dir == 0)
1660 : {
1661 : mllinop_comp_interp_coef0_x_eb
1662 : (0, i , j, k, tag.blen, tag.flo, tag.mlo, tag.ap,
1663 : tag.bctlo, tag.bcllo, imaxorder, dxi, tag.comp);
1664 : mllinop_comp_interp_coef0_x_eb
1665 : (1, i+tag.blen+1, j, k, tag.blen, tag.fhi, tag.mhi, tag.ap,
1666 : tag.bcthi, tag.bclhi, imaxorder, dxi, tag.comp);
1667 : }
1668 : #if (AMREX_SPACEDIM > 1)
1669 : else
1670 : #if (AMREX_SPACEDIM > 2)
1671 : if (tag.dir == 1)
1672 : #endif
1673 : {
1674 : mllinop_comp_interp_coef0_y_eb
1675 : (0, i, j , k, tag.blen, tag.flo, tag.mlo, tag.ap,
1676 : tag.bctlo, tag.bcllo, imaxorder, dyi, tag.comp);
1677 : mllinop_comp_interp_coef0_y_eb
1678 : (1, i, j+tag.blen+1, k, tag.blen, tag.fhi, tag.mhi, tag.ap,
1679 : tag.bcthi, tag.bclhi, imaxorder, dyi, tag.comp);
1680 : }
1681 : #if (AMREX_SPACEDIM > 2)
1682 : else {
1683 : mllinop_comp_interp_coef0_z_eb
1684 : (0, i, j, k , tag.blen, tag.flo, tag.mlo, tag.ap,
1685 : tag.bctlo, tag.bcllo, imaxorder, dzi, tag.comp);
1686 : mllinop_comp_interp_coef0_z_eb
1687 : (1, i, j, k+tag.blen+1, tag.blen, tag.fhi, tag.mhi, tag.ap,
1688 : tag.bcthi, tag.bclhi, imaxorder, dzi, tag.comp);
1689 : }
1690 : #endif
1691 : #endif
1692 : } else {
1693 : if (tag.dir == 0)
1694 : {
1695 : mllinop_comp_interp_coef0_x
1696 : (0, i , j, k, tag.blen, tag.flo, tag.mlo,
1697 : tag.bctlo, tag.bcllo, imaxorder, dxi, tag.comp);
1698 : mllinop_comp_interp_coef0_x
1699 : (1, i+tag.blen+1, j, k, tag.blen, tag.fhi, tag.mhi,
1700 : tag.bcthi, tag.bclhi, imaxorder, dxi, tag.comp);
1701 : }
1702 : #if (AMREX_SPACEDIM > 1)
1703 : else
1704 : #if (AMREX_SPACEDIM > 2)
1705 : if (tag.dir == 1)
1706 : #endif
1707 : {
1708 : mllinop_comp_interp_coef0_y
1709 : (0, i, j , k, tag.blen, tag.flo, tag.mlo,
1710 : tag.bctlo, tag.bcllo, imaxorder, dyi, tag.comp);
1711 : mllinop_comp_interp_coef0_y
1712 : (1, i, j+tag.blen+1, k, tag.blen, tag.fhi, tag.mhi,
1713 : tag.bcthi, tag.bclhi, imaxorder, dyi, tag.comp);
1714 : }
1715 : #if (AMREX_SPACEDIM > 2)
1716 : else {
1717 : mllinop_comp_interp_coef0_z
1718 : (0, i, j, k , tag.blen, tag.flo, tag.mlo,
1719 : tag.bctlo, tag.bcllo, imaxorder, dzi, tag.comp);
1720 : mllinop_comp_interp_coef0_z
1721 : (1, i, j, k+tag.blen+1, tag.blen, tag.fhi, tag.mhi,
1722 : tag.bcthi, tag.bclhi, imaxorder, dzi, tag.comp);
1723 : }
1724 : #endif
1725 : #endif
1726 : }
1727 : });
1728 : }
1729 : } else
1730 : #endif
1731 : {
1732 : Vector<MLMGPSTag<RT>> tags;
1733 : tags.reserve(foo.local_size()*AMREX_SPACEDIM*ncomp);
1734 :
1735 : for (MFIter mfi(foo); mfi.isValid(); ++mfi)
1736 : {
1737 : const Box& vbx = mfi.validbox();
1738 :
1739 : const auto & bdlv = bcondloc.bndryLocs(mfi);
1740 : const auto & bdcv = bcondloc.bndryConds(mfi);
1741 :
1742 : for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
1743 : {
1744 : if (idim != hidden_direction) {
1745 : const Orientation olo(idim,Orientation::low);
1746 : const Orientation ohi(idim,Orientation::high);
1747 : for (int icomp = 0; icomp < ncomp; ++icomp) {
1748 : tags.emplace_back(MLMGPSTag<RT>{undrrelxr[olo].array(mfi),
1749 : undrrelxr[ohi].array(mfi),
1750 : maskvals[olo].const_array(mfi),
1751 : maskvals[ohi].const_array(mfi),
1752 : bdlv[icomp][olo], bdlv[icomp][ohi],
1753 : amrex::adjCell(vbx,olo),
1754 : bdcv[icomp][olo], bdcv[icomp][ohi],
1755 : vbx.length(idim), icomp, idim});
1756 : }
1757 : }
1758 : }
1759 : }
1760 :
1761 : ParallelFor(tags,
1762 : [=] AMREX_GPU_DEVICE (int i, int j, int k, MLMGPSTag<RT> const& tag) noexcept
1763 : {
1764 : if (tag.dir == 0)
1765 : {
1766 : mllinop_comp_interp_coef0_x
1767 : (0, i , j, k, tag.blen, tag.flo, tag.mlo,
1768 : tag.bctlo, tag.bcllo, imaxorder, dxi, tag.comp);
1769 : mllinop_comp_interp_coef0_x
1770 : (1, i+tag.blen+1, j, k, tag.blen, tag.fhi, tag.mhi,
1771 : tag.bcthi, tag.bclhi, imaxorder, dxi, tag.comp);
1772 : }
1773 : #if (AMREX_SPACEDIM > 1)
1774 : else
1775 : #if (AMREX_SPACEDIM > 2)
1776 : if (tag.dir == 1)
1777 : #endif
1778 : {
1779 : mllinop_comp_interp_coef0_y
1780 : (0, i, j , k, tag.blen, tag.flo, tag.mlo,
1781 : tag.bctlo, tag.bcllo, imaxorder, dyi, tag.comp);
1782 : mllinop_comp_interp_coef0_y
1783 : (1, i, j+tag.blen+1, k, tag.blen, tag.fhi, tag.mhi,
1784 : tag.bcthi, tag.bclhi, imaxorder, dyi, tag.comp);
1785 : }
1786 : #if (AMREX_SPACEDIM > 2)
1787 : else {
1788 : mllinop_comp_interp_coef0_z
1789 : (0, i, j, k , tag.blen, tag.flo, tag.mlo,
1790 : tag.bctlo, tag.bcllo, imaxorder, dzi, tag.comp);
1791 : mllinop_comp_interp_coef0_z
1792 : (1, i, j, k+tag.blen+1, tag.blen, tag.fhi, tag.mhi,
1793 : tag.bcthi, tag.bclhi, imaxorder, dzi, tag.comp);
1794 : }
1795 : #endif
1796 : #endif
1797 : });
1798 : }
1799 : } else
1800 : #endif
1801 : {
1802 : #ifdef AMREX_USE_OMP
1803 : #pragma omp parallel
1804 : #endif
1805 : for (MFIter mfi(foo, MFItInfo{}.SetDynamic(true)); mfi.isValid(); ++mfi)
1806 : {
1807 : const Box& vbx = mfi.validbox();
1808 :
1809 : const auto & bdlv = bcondloc.bndryLocs(mfi);
1810 : const auto & bdcv = bcondloc.bndryConds(mfi);
1811 :
1812 : #ifdef AMREX_USE_EB
1813 : auto fabtyp = (flags) ? (*flags)[mfi].getType(vbx) : FabType::regular;
1814 : #endif
1815 : for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
1816 : {
1817 : if (idim == hidden_direction) { continue; }
1818 : const Orientation olo(idim,Orientation::low);
1819 : const Orientation ohi(idim,Orientation::high);
1820 : const Box blo = amrex::adjCellLo(vbx, idim);
1821 : const Box bhi = amrex::adjCellHi(vbx, idim);
1822 : const int blen = vbx.length(idim);
1823 : const auto& mlo = maskvals[olo].array(mfi);
1824 : const auto& mhi = maskvals[ohi].array(mfi);
1825 : const auto& flo = undrrelxr[olo].array(mfi);
1826 : const auto& fhi = undrrelxr[ohi].array(mfi);
1827 : for (int icomp = 0; icomp < ncomp; ++icomp) {
1828 : const BoundCond bctlo = bdcv[icomp][olo];
1829 : const BoundCond bcthi = bdcv[icomp][ohi];
1830 : const auto bcllo = bdlv[icomp][olo];
1831 : const auto bclhi = bdlv[icomp][ohi];
1832 : #ifdef AMREX_USE_EB
1833 : if (fabtyp == FabType::singlevalued) {
1834 : if constexpr (!std::is_same<MF,MultiFab>()) {
1835 : amrex::Abort("MLCellLinOp with EB only works with MultiFab");
1836 : } else {
1837 : auto const& ap = area[idim]->const_array(mfi);
1838 : if (idim == 0) {
1839 : mllinop_comp_interp_coef0_x_eb
1840 : (0, blo, blen, flo, mlo, ap, bctlo, bcllo,
1841 : imaxorder, dxi, icomp);
1842 : mllinop_comp_interp_coef0_x_eb
1843 : (1, bhi, blen, fhi, mhi, ap, bcthi, bclhi,
1844 : imaxorder, dxi, icomp);
1845 : } else if (idim == 1) {
1846 : mllinop_comp_interp_coef0_y_eb
1847 : (0, blo, blen, flo, mlo, ap, bctlo, bcllo,
1848 : imaxorder, dyi, icomp);
1849 : mllinop_comp_interp_coef0_y_eb
1850 : (1, bhi, blen, fhi, mhi, ap, bcthi, bclhi,
1851 : imaxorder, dyi, icomp);
1852 : } else {
1853 : mllinop_comp_interp_coef0_z_eb
1854 : (0, blo, blen, flo, mlo, ap, bctlo, bcllo,
1855 : imaxorder, dzi, icomp);
1856 : mllinop_comp_interp_coef0_z_eb
1857 : (1, bhi, blen, fhi, mhi, ap, bcthi, bclhi,
1858 : imaxorder, dzi, icomp);
1859 : }
1860 : }
1861 : } else if (fabtyp == FabType::regular)
1862 : #endif
1863 : {
1864 : if (idim == 0) {
1865 : mllinop_comp_interp_coef0_x
1866 : (0, blo, blen, flo, mlo, bctlo, bcllo,
1867 : imaxorder, dxi, icomp);
1868 : mllinop_comp_interp_coef0_x
1869 : (1, bhi, blen, fhi, mhi, bcthi, bclhi,
1870 : imaxorder, dxi, icomp);
1871 : } else if (idim == 1) {
1872 : mllinop_comp_interp_coef0_y
1873 : (0, blo, blen, flo, mlo, bctlo, bcllo,
1874 : imaxorder, dyi, icomp);
1875 : mllinop_comp_interp_coef0_y
1876 : (1, bhi, blen, fhi, mhi, bcthi, bclhi,
1877 : imaxorder, dyi, icomp);
1878 : } else {
1879 : mllinop_comp_interp_coef0_z
1880 : (0, blo, blen, flo, mlo, bctlo, bcllo,
1881 : imaxorder, dzi, icomp);
1882 : mllinop_comp_interp_coef0_z
1883 : (1, bhi, blen, fhi, mhi, bcthi, bclhi,
1884 : imaxorder, dzi, icomp);
1885 : }
1886 : }
1887 : }
1888 : }
1889 : }
1890 : }
1891 : }
1892 : }
1893 : }
1894 :
1895 : template <typename MF>
1896 : auto
1897 : MLCellLinOpT<MF>::xdoty (int /*amrlev*/, int /*mglev*/, const MF& x, const MF& y, bool local) const
1898 : -> RT
1899 : {
1900 : const int ncomp = this->getNComp();
1901 : const IntVect nghost(0);
1902 : RT result = amrex::Dot(x,0,y,0,ncomp,nghost,true);
1903 : if (!local) {
1904 : ParallelAllReduce::Sum(result, ParallelContext::CommunicatorSub());
1905 : }
1906 : return result;
1907 : }
1908 :
1909 : template <typename MF>
1910 : void
1911 : MLCellLinOpT<MF>::computeVolInv () const
1912 : {
1913 : if (!m_volinv.empty()) { return; }
1914 :
1915 : m_volinv.resize(this->m_num_amr_levels);
1916 : for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
1917 : m_volinv[amrlev].resize(this->NMGLevels(amrlev));
1918 : }
1919 :
1920 : AMREX_ASSERT(this->m_coarse_fine_bc_type == LinOpBCType::Dirichlet || ! this->hasHiddenDimension());
1921 :
1922 : // We don't need to compute for every level
1923 :
1924 : auto f = [&] (int amrlev, int mglev) {
1925 : #ifdef AMREX_USE_EB
1926 : const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(amrlev,mglev));
1927 : if (factory && !factory->isAllRegular())
1928 : {
1929 : if constexpr (std::is_same<MF,MultiFab>()) {
1930 : const auto& vfrac = factory->getVolFrac();
1931 : m_volinv[amrlev][mglev] = vfrac.sum(0,true);
1932 : } else {
1933 : amrex::Abort("MLCellLinOp with EB only works with MultiFab");
1934 : }
1935 : }
1936 : else
1937 : #endif
1938 : {
1939 : if (this->m_coarse_fine_bc_type == LinOpBCType::Dirichlet) {
1940 : m_volinv[amrlev][mglev]
1941 : = RT(1.0 / this->compactify(this->Geom(amrlev,mglev).Domain()).d_numPts());
1942 : } else {
1943 : m_volinv[amrlev][mglev]
1944 : = RT(1.0 / this->m_grids[amrlev][mglev].d_numPts());
1945 : }
1946 : }
1947 : };
1948 :
1949 : // amrlev = 0, mglev = 0
1950 : f(0,0);
1951 :
1952 : int mgbottom = this->NMGLevels(0)-1;
1953 : f(0,mgbottom);
1954 :
1955 : #ifdef AMREX_USE_EB
1956 : RT temp1, temp2;
1957 : const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(0,0));
1958 : if (factory && !factory->isAllRegular())
1959 : {
1960 : ParallelAllReduce::Sum<RT>({m_volinv[0][0], m_volinv[0][mgbottom]},
1961 : ParallelContext::CommunicatorSub());
1962 : temp1 = RT(1.0)/m_volinv[0][0];
1963 : temp2 = RT(1.0)/m_volinv[0][mgbottom];
1964 : }
1965 : else
1966 : {
1967 : temp1 = m_volinv[0][0];
1968 : temp2 = m_volinv[0][mgbottom];
1969 : }
1970 : m_volinv[0][0] = temp1;
1971 : m_volinv[0][mgbottom] = temp2;
1972 : #endif
1973 : }
1974 :
1975 : template <typename MF>
1976 : auto
1977 : MLCellLinOpT<MF>::normInf (int amrlev, MF const& mf, bool local) const -> RT
1978 : {
1979 : const int ncomp = this->getNComp();
1980 : const int finest_level = this->NAMRLevels() - 1;
1981 : RT norm = RT(0.0);
1982 : #ifdef AMREX_USE_EB
1983 : if (! mf.isAllRegular()) {
1984 : if constexpr (!std::is_same<MF,MultiFab>()) {
1985 : amrex::Abort("MLCellLinOpT with EB only works with MultiFab");
1986 : } else {
1987 : const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(amrlev));
1988 : const MultiFab& vfrac = factory->getVolFrac();
1989 : if (amrlev == finest_level) {
1990 : #ifdef AMREX_USE_GPU
1991 : if (Gpu::inLaunchRegion()) {
1992 : auto const& ma = mf.const_arrays();
1993 : auto const& vfrac_ma = vfrac.const_arrays();
1994 : norm = ParReduce(TypeList<ReduceOpMax>{}, TypeList<Real>{},
1995 : mf, IntVect(0), ncomp,
1996 : [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n)
1997 : -> GpuTuple<Real>
1998 : {
1999 : return std::abs(ma[box_no](i,j,k,n)
2000 : *vfrac_ma[box_no](i,j,k));
2001 : });
2002 : } else
2003 : #endif
2004 : {
2005 : #ifdef AMREX_USE_OMP
2006 : #pragma omp parallel reduction(max:norm)
2007 : #endif
2008 : for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) {
2009 : Box const& bx = mfi.tilebox();
2010 : auto const& fab = mf.const_array(mfi);
2011 : auto const& v = vfrac.const_array(mfi);
2012 : AMREX_LOOP_4D(bx, ncomp, i, j, k, n,
2013 : {
2014 : norm = std::max(norm, std::abs(fab(i,j,k,n)*v(i,j,k)));
2015 : });
2016 : }
2017 : }
2018 : } else {
2019 : #ifdef AMREX_USE_GPU
2020 : if (Gpu::inLaunchRegion()) {
2021 : auto const& ma = mf.const_arrays();
2022 : auto const& mask_ma = m_norm_fine_mask[amrlev]->const_arrays();
2023 : auto const& vfrac_ma = vfrac.const_arrays();
2024 : norm = ParReduce(TypeList<ReduceOpMax>{}, TypeList<Real>{},
2025 : mf, IntVect(0), ncomp,
2026 : [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n)
2027 : -> GpuTuple<Real>
2028 : {
2029 : if (mask_ma[box_no](i,j,k)) {
2030 : return std::abs(ma[box_no](i,j,k,n)
2031 : *vfrac_ma[box_no](i,j,k));
2032 : } else {
2033 : return Real(0.0);
2034 : }
2035 : });
2036 : } else
2037 : #endif
2038 : {
2039 : #ifdef AMREX_USE_OMP
2040 : #pragma omp parallel reduction(max:norm)
2041 : #endif
2042 : for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) {
2043 : Box const& bx = mfi.tilebox();
2044 : auto const& fab = mf.const_array(mfi);
2045 : auto const& mask = m_norm_fine_mask[amrlev]->const_array(mfi);
2046 : auto const& v = vfrac.const_array(mfi);
2047 : AMREX_LOOP_4D(bx, ncomp, i, j, k, n,
2048 : {
2049 : if (mask(i,j,k)) {
2050 : norm = std::max(norm, std::abs(fab(i,j,k,n)*v(i,j,k)));
2051 : }
2052 : });
2053 : }
2054 : }
2055 : }
2056 : }
2057 : } else
2058 : #endif
2059 : {
2060 : if (amrlev == finest_level) {
2061 : norm = mf.norminf(0, ncomp, IntVect(0), true);
2062 : } else {
2063 : norm = mf.norminf(*m_norm_fine_mask[amrlev], 0, ncomp, IntVect(0), true);
2064 : }
2065 : }
2066 :
2067 : if (!local) { ParallelAllReduce::Max(norm, ParallelContext::CommunicatorSub()); }
2068 : return norm;
2069 : }
2070 :
2071 : template <typename MF>
2072 : void
2073 : MLCellLinOpT<MF>::averageDownAndSync (Vector<MF>& sol) const
2074 : {
2075 : int ncomp = this->getNComp();
2076 : for (int falev = this->NAMRLevels()-1; falev > 0; --falev)
2077 : {
2078 : #ifdef AMREX_USE_EB
2079 : if (!sol[falev].isAllRegular()) {
2080 : if constexpr (std::is_same<MF,MultiFab>()) {
2081 : amrex::EB_average_down(sol[falev], sol[falev-1], 0, ncomp, this->AMRRefRatio(falev-1));
2082 : } else {
2083 : amrex::Abort("EB_average_down only works with MultiFab");
2084 : }
2085 : } else
2086 : #endif
2087 : {
2088 : amrex::average_down(sol[falev], sol[falev-1], 0, ncomp, this->AMRRefRatio(falev-1));
2089 : }
2090 : }
2091 : }
2092 :
2093 : template <typename MF>
2094 : void
2095 : MLCellLinOpT<MF>::avgDownResAmr (int clev, MF& cres, MF const& fres) const
2096 : {
2097 : #ifdef AMREX_USE_EB
2098 : if (!fres.isAllRegular()) {
2099 : if constexpr (std::is_same<MF,MultiFab>()) {
2100 : amrex::EB_average_down(fres, cres, 0, this->getNComp(),
2101 : this->AMRRefRatio(clev));
2102 : } else {
2103 : amrex::Abort("EB_average_down only works with MultiFab");
2104 : }
2105 : } else
2106 : #endif
2107 : {
2108 : amrex::average_down(fres, cres, 0, this->getNComp(),
2109 : this->AMRRefRatio(clev));
2110 : }
2111 : }
2112 :
2113 : extern template class MLCellLinOpT<MultiFab>;
2114 :
2115 : using MLCellLinOp = MLCellLinOpT<MultiFab>;
2116 :
2117 : }
2118 :
2119 : #endif
|
