Block-Structured AMR Software Framework
Loading...
Searching...
No Matches
AMReX_MLCellLinOp.H
Go to the documentation of this file.
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
16namespace amrex {
17
25template <typename T>
26struct MLMGABCTag;
27
28template <typename MF>
29class MLCellLinOpT // NOLINT(cppcoreguidelines-virtual-class-destructor)
30 : public MLLinOpT<MF>
31{
32public:
33
34 using FAB = typename FabDataType<MF>::fab_type;
35 using RT = typename FabDataType<MF>::value_type;
36
37 using BCType = LinOpBCType;
38 using BCMode = typename MLLinOpT<MF>::BCMode;
39 using StateMode = typename MLLinOpT<MF>::StateMode;
40 using Location = typename MLLinOpT<MF>::Location;
41
42 MLCellLinOpT ();
43 ~MLCellLinOpT () override = default;
44
45 MLCellLinOpT (const MLCellLinOpT<MF>&) = delete;
46 MLCellLinOpT (MLCellLinOpT<MF>&&) = delete;
47 MLCellLinOpT<MF>& operator= (const MLCellLinOpT<MF>&) = delete;
48 MLCellLinOpT<MF>& operator= (MLCellLinOpT<MF>&&) = delete;
49
59 void define (const Vector<Geometry>& a_geom,
60 const Vector<BoxArray>& a_grids,
61 const Vector<DistributionMapping>& a_dmap,
62 const LPInfo& a_info = LPInfo(),
63 const Vector<FabFactory<FAB> const*>& a_factory = {});
64
78 void setLevelBC (int amrlev, const MF* levelbcdata,
79 const MF* robinbc_a = nullptr,
80 const MF* robinbc_b = nullptr,
81 const MF* robinbc_f = nullptr) final;
82
83 template <MultiFabLike AMF>
84 requires (!std::same_as<MF,AMF>)
85 void setLevelBC (int amrlev, const AMF* levelbcdata,
86 const AMF* robinbc_a = nullptr,
87 const AMF* robinbc_b = nullptr,
88 const AMF* robinbc_f = nullptr)
89 {
90 this->MLLinOpT<MF>::template setLevelBC<AMF>(amrlev, levelbcdata, robinbc_a, robinbc_b, robinbc_f);
91 }
92
93 bool needsUpdate () const override {
94 return MLLinOpT<MF>::needsUpdate();
95 }
96 void update () override;
97
104 void setGaussSeidel (bool flag) noexcept { m_use_gauss_seidel = flag; }
105
107 virtual bool isCrossStencil () const { return true; }
109 virtual bool isTensorOp () const { return false; }
110
117 void updateSolBC (int amrlev, const MF& crse_bcdata) const;
124 void updateCorBC (int amrlev, const MF& crse_bcdata) const;
125
137 virtual void applyBC (int amrlev, int mglev, MF& in, BCMode bc_mode, StateMode s_mode,
138 const MLMGBndryT<MF>* bndry=nullptr, bool skip_fillboundary=false) const;
139
141 BoxArray makeNGrids (int grid_size) const;
142
151 void restriction (int amrlev, int cmglev, MF& crse, MF& fine) const override;
152
162 void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const override;
163
173 void interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const override;
174
184 void interpolationAmr (int famrlev, MF& fine, const MF& crse,
185 IntVect const& nghost) const override;
186
196 void averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs,
197 const MF& fine_sol, const MF& fine_rhs) override;
198
210 void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode,
211 StateMode s_mode, const MLMGBndryT<MF>* bndry=nullptr) const override;
222 void smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
223 bool skip_fillboundary, int niter) const override;
224
234 void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b,
235 const MF* crse_bcdata=nullptr) override;
236
243 void prepareForFluxes (int amrlev, const MF* crse_bcdata = nullptr) override;
244
256 void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b,
257 BCMode bc_mode, const MF* crse_bcdata=nullptr) final;
258
272 void reflux (int crse_amrlev,
273 MF& res, const MF& crse_sol, const MF&,
274 MF&, MF& fine_sol, const MF&) const final;
283 void compFlux (int amrlev, const Array<MF*,AMREX_SPACEDIM>& fluxes,
284 MF& sol, Location loc) const override;
293 void compGrad (int amrlev, const Array<MF*,AMREX_SPACEDIM>& grad,
294 MF& sol, Location loc) const override;
295
303 void applyMetricTerm (int amrlev, int mglev, MF& rhs) const final;
311 void unapplyMetricTerm (int amrlev, int mglev, MF& rhs) const final;
320 Vector<RT> getSolvabilityOffset (int amrlev, int mglev,
321 MF const& rhs) const override;
330 void fixSolvabilityByOffset (int amrlev, int mglev, MF& rhs,
331 Vector<RT> const& offset) const override;
332
334 void prepareForSolve () override;
335
346 RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const final;
347
356 RT dotProductPrecond (Vector<MF const*> const& x, Vector<MF const*> const& y) const final;
357
365 RT norm2Precond (Vector<MF const*> const& x) const final;
366
367 virtual void Fapply (int amrlev, int mglev, MF& out, const MF& in) const = 0;
368 virtual void Fsmooth (int amrlev, int mglev, MF& sol, const MF& rhs, int redblack) const = 0;
369 virtual void FFlux (int amrlev, const MFIter& mfi,
370 const Array<FAB*,AMREX_SPACEDIM>& flux,
371 const FAB& sol, Location loc, int face_only=0) const = 0;
372
381 virtual void addInhomogNeumannFlux (int /*amrlev*/,
382 const Array<MF*,AMREX_SPACEDIM>& /*grad*/,
383 MF const& /*sol*/,
384 bool /*mult_bcoef*/) const {}
385
393 RT normInf (int amrlev, MF const& mf, bool local) const override;
394
400 void averageDownAndSync (Vector<MF>& sol) const override;
401
409 void avgDownResAmr (int clev, MF& cres, MF const& fres) const override;
410
412 void beginPrecondBC () override;
414 void endPrecondBC () override;
415
417 struct BCTL {
418 BoundCond type;
419 RT location;
420 };
422
423 Vector<std::unique_ptr<MF> > m_robin_bcval;
424
430 void setInterpBndryHalfWidth (int w) { m_interpbndry_halfwidth = w; }
431
432protected:
433
434 bool m_has_metric_term = false;
435
436 Vector<std::unique_ptr<MLMGBndryT<MF>> > m_bndry_sol;
437 Vector<std::unique_ptr<BndryRegisterT<MF>> > m_crse_sol_br;
438
439 Vector<std::unique_ptr<MLMGBndryT<MF>> > m_bndry_cor;
440 Vector<std::unique_ptr<BndryRegisterT<MF>> > m_crse_cor_br;
441
442 Vector<std::unique_ptr<MLMGBndryT<MF>> > m_bndry_sol_zero;
443
444 // In case of agglomeration, coarse MG grids on amr level 0 are
445 // not simply coarsened from fine MG grids. So we need to build
446 // bcond and bcloc for each MG level.
447 using RealTuple = Array<RT,2*BL_SPACEDIM>;
448 using BCTuple = Array<BoundCond,2*BL_SPACEDIM>;
450 class BndryCondLoc
451 {
452 public:
453 BndryCondLoc (const BoxArray& ba, const DistributionMapping& dm, int ncomp);
454
455 void setLOBndryConds (const Geometry& geom, const Real* dx,
456 const Vector<Array<BCType,AMREX_SPACEDIM> >& lobc,
457 const Vector<Array<BCType,AMREX_SPACEDIM> >& hibc,
458 IntVect const& ratio, const RealVect& interior_bloc,
459 const Array<Real,AMREX_SPACEDIM>& domain_bloc_lo,
460 const Array<Real,AMREX_SPACEDIM>& domain_bloc_hi,
461 LinOpBCType crse_fine_bc_type);
462
463 const Vector<BCTuple>& bndryConds (const MFIter& mfi) const noexcept {
464 return bcond[mfi];
465 }
466 const Vector<RealTuple>& bndryLocs (const MFIter& mfi) const noexcept {
467 return bcloc[mfi];
468 }
469 const BCTuple& bndryConds (const MFIter& mfi, int icomp) const noexcept {
470 return bcond[mfi][icomp];
471 }
472 const RealTuple& bndryLocs (const MFIter& mfi, int icomp) const noexcept {
473 return bcloc[mfi][icomp];
474 }
475 GpuArray<BCTL,2*AMREX_SPACEDIM> const* getBCTLPtr (const MFIter& mfi) const noexcept {
476 return bctl[mfi];
477 }
478 private:
479 LayoutData<Vector<BCTuple> > bcond;
480 LayoutData<Vector<RealTuple> > bcloc;
481 LayoutData<GpuArray<BCTL,2*AMREX_SPACEDIM>*> bctl;
482 Gpu::DeviceVector<GpuArray<BCTL,2*AMREX_SPACEDIM> > bctl_dv;
483 int m_ncomp;
484 };
486 Vector<Vector<std::unique_ptr<BndryCondLoc> > > m_bcondloc;
487
488 // used to save interpolation coefficients of the first interior cells
489 mutable Vector<Vector<BndryRegisterT<MF>> > m_undrrelxr;
490
491 // boundary cell flags for covered, not_covered, outside_domain
492 Vector<Vector<Array<MultiMask,2*AMREX_SPACEDIM> > > m_maskvals;
493
494 Vector<std::unique_ptr<iMultiFab> > m_norm_fine_mask;
495
496 mutable Vector<YAFluxRegisterT<MF>> m_fluxreg;
497
498 bool m_use_gauss_seidel = true; // use red-black Gauss-Seidel by default
499
500private:
501
502 void defineAuxData ();
503 void defineBC ();
504
505 void computeVolInv () const;
506 mutable Vector<Vector<RT> > m_volinv; // used by solvability fix
507
508 int m_interpbndry_halfwidth = 2;
509
510 mutable Vector<Vector<TagVector<MLMGABCTag<RT>>>> m_bc_tags;
511};
512
514
515template <typename T>
516struct MLMGABCTag {
517 Array4<int const> mask;
518 T bcloc;
519 Box bx;
520 BoundCond bctype;
521 int blen;
522 int comp;
523 Orientation face;
524 int local_index;
525
526 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
527 Box const& box() const noexcept { return bx; }
528};
529
530template <typename T>
531struct MLMGPSTag {
532 Array4<T> flo;
533 Array4<T> fhi;
534 Array4<int const> mlo;
535 Array4<int const> mhi;
536 T bcllo;
537 T bclhi;
538 Box bx;
539 BoundCond bctlo;
540 BoundCond bcthi;
541 int blen;
542 int comp;
543 int dir;
544
545 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
546 Box const& box() const noexcept { return bx; }
547};
548
549#ifdef AMREX_USE_EB
550template <typename T>
551struct MLMGPSEBTag {
552 Array4<T> flo;
553 Array4<T> fhi;
554 Array4<T const> ap;
555 Array4<int const> mlo;
556 Array4<int const> mhi;
557 T bcllo;
558 T bclhi;
559 Box bx;
560 BoundCond bctlo;
561 BoundCond bcthi;
562 int blen;
563 int comp;
564 int dir;
565
566 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
567 Box const& box() const noexcept { return bx; }
568};
569#endif
570
571template <typename MF>
572MLCellLinOpT<MF>::BndryCondLoc::BndryCondLoc (const BoxArray& ba,
573 const DistributionMapping& dm,
574 int ncomp)
575 : bcond(ba, dm),
576 bcloc(ba, dm),
577 bctl(ba, dm),
578 bctl_dv(bctl.local_size()*ncomp),
579 m_ncomp(ncomp)
580{
581 auto* dp = bctl_dv.data();
582 for (MFIter mfi(bcloc); mfi.isValid(); ++mfi) {
583 bcond[mfi].resize(ncomp);
584 bcloc[mfi].resize(ncomp);
585 bctl[mfi] = dp;
586 dp += ncomp;
587 }
588}
589
590template <typename MF>
591void
592MLCellLinOpT<MF>::BndryCondLoc::
593setLOBndryConds (const Geometry& geom, const Real* dx,
594 const Vector<Array<BCType,AMREX_SPACEDIM> >& lobc,
595 const Vector<Array<BCType,AMREX_SPACEDIM> >& hibc,
596 IntVect const& ratio, const RealVect& interior_bloc,
597 const Array<Real,AMREX_SPACEDIM>& domain_bloc_lo,
598 const Array<Real,AMREX_SPACEDIM>& domain_bloc_hi,
599 LinOpBCType crse_fine_bc_type)
600{
601 const Box& domain = geom.Domain();
602
603#ifdef AMREX_USE_OMP
604#pragma omp parallel
605#endif
606 for (MFIter mfi(bcloc); mfi.isValid(); ++mfi)
607 {
608 const Box& bx = mfi.validbox();
609 for (int icomp = 0; icomp < m_ncomp; ++icomp) {
610 RealTuple & bloc = bcloc[mfi][icomp];
611 BCTuple & bctag = bcond[mfi][icomp];
612 MLMGBndryT<MF>::setBoxBC(bloc, bctag, bx, domain,
613 lobc[icomp], hibc[icomp],
614 dx, ratio, interior_bloc,
615 domain_bloc_lo, domain_bloc_hi,
616 geom.isPeriodicArray(),
617 crse_fine_bc_type);
618 }
619 }
620
621 Gpu::PinnedVector<GpuArray<BCTL,2*AMREX_SPACEDIM> > hv;
622 hv.reserve(bctl_dv.size());
623 for (MFIter mfi(bctl); mfi.isValid(); ++mfi)
624 {
625 for (int icomp = 0; icomp < m_ncomp; ++icomp) {
626 GpuArray<BCTL,2*AMREX_SPACEDIM> tmp;
627 for (int m = 0; m < 2*AMREX_SPACEDIM; ++m) {
628 tmp[m].type = bcond[mfi][icomp][m];
629 tmp[m].location = bcloc[mfi][icomp][m];
630 }
631 hv.push_back(std::move(tmp));
632 }
633 }
634 Gpu::copyAsync(Gpu::hostToDevice, hv.begin(), hv.end(), bctl_dv.begin());
635 Gpu::streamSynchronize();
636}
637
639
640template <typename MF>
641MLCellLinOpT<MF>::MLCellLinOpT ()
642{
643 this->m_ixtype = IntVect::TheCellVector();
644}
645
646template <typename MF>
647void
648MLCellLinOpT<MF>::define (const Vector<Geometry>& a_geom,
649 const Vector<BoxArray>& a_grids,
650 const Vector<DistributionMapping>& a_dmap,
651 const LPInfo& a_info,
652 const Vector<FabFactory<FAB> const*>& a_factory)
653{
654 // defineAuxData needs all MG levels, which a deferred hierarchy lacks.
655 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!this->supportsAnisotropicCoarsening(),
656 "MLCellLinOp: anisotropic coarsening is not supported");
657 MLLinOpT<MF>::define(a_geom, a_grids, a_dmap, a_info, a_factory);
658 defineAuxData();
659 defineBC();
660}
661
662template <typename MF>
663void
664MLCellLinOpT<MF>::defineAuxData ()
665{
666 BL_PROFILE("MLCellLinOp::defineAuxData()");
667
668 m_undrrelxr.resize(this->m_num_amr_levels);
669 m_maskvals.resize(this->m_num_amr_levels);
670 m_fluxreg.resize(this->m_num_amr_levels-1);
671 m_norm_fine_mask.resize(this->m_num_amr_levels-1);
672 m_bc_tags.resize(this->m_num_amr_levels);
673
674 const int ncomp = this->getNComp();
675
676 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
677 {
678 m_undrrelxr[amrlev].resize(this->m_num_mg_levels[amrlev]);
679 m_bc_tags[amrlev].resize(this->m_num_mg_levels[amrlev]);
680 for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
681 {
682 m_undrrelxr[amrlev][mglev].define(this->m_grids[amrlev][mglev],
683 this->m_dmap[amrlev][mglev],
684 1, 0, 0, ncomp);
685 }
686 }
687
688 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
689 {
690 m_maskvals[amrlev].resize(this->m_num_mg_levels[amrlev]);
691 for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
692 {
693 for (OrientationIter oitr; oitr; ++oitr)
694 {
695 const Orientation face = oitr();
696 const int ngrow = 1;
697 const int extent = this->isCrossStencil() ? 0 : 1; // extend to corners
698 m_maskvals[amrlev][mglev][face].define(this->m_grids[amrlev][mglev],
699 this->m_dmap[amrlev][mglev],
700 this->m_geom[amrlev][mglev],
701 face, 0, ngrow, extent, 1, true);
702 }
703 }
704 }
705
706 for (int amrlev = 0; amrlev < this->m_num_amr_levels-1; ++amrlev)
707 {
708 const IntVect ratio{this->AMRRefRatioVect(amrlev)};
709 m_fluxreg[amrlev].define(this->m_grids[amrlev+1][0],
710 this->m_grids[amrlev][0],
711 this->m_dmap[amrlev+1][0],
712 this->m_dmap[amrlev][0],
713 this->m_geom[amrlev+1][0],
714 this->m_geom[amrlev][0],
715 ratio, amrlev+1, ncomp);
716 m_fluxreg[amrlev].setDeterministic(this->info.deterministic);
717 m_norm_fine_mask[amrlev] = std::make_unique<iMultiFab>
718 (makeFineMask(this->m_grids[amrlev][0], this->m_dmap[amrlev][0],
719 this->m_grids[amrlev+1][0],
720 ratio, 1, 0));
721 }
722
723#if (AMREX_SPACEDIM != 3)
724 m_has_metric_term = !this->m_geom[0][0].IsCartesian() && this->info.has_metric_term;
725#endif
726}
727
728template <typename MF>
729void
730MLCellLinOpT<MF>::defineBC ()
731{
732 BL_PROFILE("MLCellLinOp::defineBC()");
733
734 const int ncomp = this->getNComp();
735
736 m_bndry_sol.resize(this->m_num_amr_levels);
737 m_crse_sol_br.resize(this->m_num_amr_levels);
738
739 m_bndry_cor.resize(this->m_num_amr_levels);
740 m_crse_cor_br.resize(this->m_num_amr_levels);
741
742 m_robin_bcval.resize(this->m_num_amr_levels);
743
744 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
745 {
746 m_bndry_sol[amrlev] = std::make_unique<MLMGBndryT<MF>>(this->m_grids[amrlev][0],
747 this->m_dmap[amrlev][0],
748 ncomp,
749 this->m_geom[amrlev][0]);
750 }
751
752 for (int amrlev = 1; amrlev < this->m_num_amr_levels; ++amrlev)
753 {
754 const int in_rad = 0;
755 const int out_rad = 1;
756 const int extent_rad = 2;
757 const IntVect crse_ratio = this->AMRRefRatioVect(amrlev-1);
758 BoxArray cba = this->m_grids[amrlev][0];
759 cba.coarsen(crse_ratio);
760 m_crse_sol_br[amrlev] = std::make_unique<BndryRegisterT<MF>>
761 (cba, this->m_dmap[amrlev][0], in_rad, out_rad, extent_rad, ncomp);
762 }
763
764 for (int amrlev = 1; amrlev < this->m_num_amr_levels; ++amrlev)
765 {
766 const int in_rad = 0;
767 const int out_rad = 1;
768 const int extent_rad = 2;
769 const IntVect crse_ratio = this->AMRRefRatioVect(amrlev-1);
770 BoxArray cba = this->m_grids[amrlev][0];
771 cba.coarsen(crse_ratio);
772 m_crse_cor_br[amrlev] = std::make_unique<BndryRegisterT<MF>>
773 (cba, this->m_dmap[amrlev][0], in_rad, out_rad, extent_rad, ncomp);
774 m_crse_cor_br[amrlev]->setVal(RT(0.0));
775 }
776
777 // This has be to done after m_crse_cor_br is defined.
778 for (int amrlev = 1; amrlev < this->m_num_amr_levels; ++amrlev)
779 {
780 m_bndry_cor[amrlev] = std::make_unique<MLMGBndryT<MF>>
781 (this->m_grids[amrlev][0], this->m_dmap[amrlev][0], ncomp, this->m_geom[amrlev][0]);
782 MF bc_data(this->m_grids[amrlev][0], this->m_dmap[amrlev][0], ncomp, 1);
783 bc_data.setVal(0.0);
784
785 m_bndry_cor[amrlev]->setBndryValues(*m_crse_cor_br[amrlev], 0, bc_data, 0, 0, ncomp,
786 this->AMRRefRatioVect(amrlev-1),
787 InterpBndryDataT<MF>::IBD_max_order_DEF,
788 m_interpbndry_halfwidth);
789
790 Vector<Array<LinOpBCType,AMREX_SPACEDIM> > bclohi
791 (ncomp,Array<LinOpBCType,AMREX_SPACEDIM>{{AMREX_D_DECL(BCType::Dirichlet,
792 BCType::Dirichlet,
793 BCType::Dirichlet)}});
794 m_bndry_cor[amrlev]->setLOBndryConds(bclohi, bclohi, this->AMRRefRatioVect(amrlev-1), RealVect{});
795 }
796
797 m_bcondloc.resize(this->m_num_amr_levels);
798 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
799 {
800 m_bcondloc[amrlev].resize(this->m_num_mg_levels[amrlev]);
801 for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
802 {
803 m_bcondloc[amrlev][mglev] = std::make_unique<BndryCondLoc>(this->m_grids[amrlev][mglev],
804 this->m_dmap[amrlev][mglev],
805 ncomp);
806 }
807 }
808}
809
810template <typename MF>
811void
812MLCellLinOpT<MF>::setLevelBC (int amrlev, const MF* a_levelbcdata, const MF* robinbc_a,
813 const MF* robinbc_b, const MF* robinbc_f)
814{
815 BL_PROFILE("MLCellLinOp::setLevelBC()");
816
817 AMREX_ALWAYS_ASSERT(amrlev >= 0 && amrlev < this->m_num_amr_levels);
818
819 const int ncomp = this->getNComp();
820
821 MF zero;
822 IntVect ng(1);
823 if (this->hasHiddenDimension()) { ng[this->hiddenDirection()] = 0; }
824 if (a_levelbcdata == nullptr) {
825 zero.define(this->m_grids[amrlev][0], this->m_dmap[amrlev][0], ncomp, ng);
826 zero.setVal(RT(0.0));
827 } else {
828 AMREX_ALWAYS_ASSERT(a_levelbcdata->nGrowVect().allGE(ng));
829 }
830 const MF& bcdata = (a_levelbcdata == nullptr) ? zero : *a_levelbcdata;
831
832 IntVect br_ref_ratio(-1);
833
834 if (amrlev == 0)
835 {
836 if (this->needsCoarseDataForBC())
837 {
838 // AMREX_ALWAYS_ASSERT(!this->hasHiddenDimension());
839 if (this->hasHiddenDimension()) {
840 int hidden_dir = this->hiddenDirection();
841 AMREX_ALWAYS_ASSERT(this->m_coarse_data_crse_ratio[hidden_dir] == 1);
842 }
843 br_ref_ratio = this->m_coarse_data_crse_ratio.allGT(0) ? this->m_coarse_data_crse_ratio : IntVect(2);
844 if (this->m_crse_sol_br[amrlev] == nullptr && br_ref_ratio.allGT(0))
845 {
846 const int in_rad = 0;
847 const int out_rad = 1;
848 const int extent_rad = 2;
849 const IntVect crse_ratio = br_ref_ratio;
850 BoxArray cba = this->m_grids[amrlev][0];
851 cba.coarsen(crse_ratio);
852 this->m_crse_sol_br[amrlev] = std::make_unique<BndryRegisterT<MF>>
853 (cba, this->m_dmap[amrlev][0], in_rad, out_rad, extent_rad, ncomp);
854 }
855 if (this->m_coarse_data_for_bc != nullptr) {
856 AMREX_ALWAYS_ASSERT(this->m_coarse_data_crse_ratio.allGT(0));
857 const Box& cbx = amrex::coarsen(this->m_geom[0][0].Domain(), this->m_coarse_data_crse_ratio);
858 this->m_crse_sol_br[amrlev]->copyFrom(*(this->m_coarse_data_for_bc), 0, 0, 0, ncomp,
859 this->m_geom[0][0].periodicity(cbx));
860 } else {
861 this->m_crse_sol_br[amrlev]->setVal(RT(0.0));
862 }
863 this->m_bndry_sol[amrlev]->setBndryValues(*(this->m_crse_sol_br[amrlev]), 0,
864 bcdata, 0, 0, ncomp, br_ref_ratio,
865 InterpBndryDataT<MF>::IBD_max_order_DEF,
866 this->m_interpbndry_halfwidth);
867 br_ref_ratio = this->m_coarse_data_crse_ratio;
868 }
869 else
870 {
871 this->m_bndry_sol[amrlev]->setPhysBndryValues(bcdata,0,0,ncomp);
872 br_ref_ratio = IntVect(1);
873 }
874 }
875 else
876 {
877 this->m_bndry_sol[amrlev]->setPhysBndryValues(bcdata,0,0,ncomp);
878 br_ref_ratio = this->AMRRefRatioVect(amrlev-1);
879 }
880
881 auto crse_fine_bc_type = (amrlev == 0) ? this->m_coarse_fine_bc_type : LinOpBCType::Dirichlet;
882 this->m_bndry_sol[amrlev]->setLOBndryConds(this->m_lobc, this->m_hibc, br_ref_ratio,
883 this->m_coarse_bc_loc, crse_fine_bc_type);
884
885 const Real* dx = this->m_geom[amrlev][0].CellSize();
886 for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
887 {
888 this->m_bcondloc[amrlev][mglev]->setLOBndryConds(this->m_geom[amrlev][mglev], dx,
889 this->m_lobc, this->m_hibc,
890 br_ref_ratio, this->m_coarse_bc_loc,
891 this->m_domain_bloc_lo, this->m_domain_bloc_hi,
892 crse_fine_bc_type);
893 m_bc_tags[amrlev][mglev].undefine();
894 }
895
896 // In precond mode beginPrecondBC calls this with null Robin data; keep
897 // the user's a, b and f rather than dereferencing the null pointers.
898 if (this->hasRobinBC() && !this->m_precond_mode) {
899 AMREX_ASSERT(robinbc_a != nullptr && robinbc_b != nullptr && robinbc_f != nullptr);
900 this->m_robin_bcval[amrlev] = std::make_unique<MF>(this->m_grids[amrlev][0],
901 this->m_dmap[amrlev][0],
902 ncomp*3, 1);
903 const Box& domain = this->m_geom[amrlev][0].Domain();
904 MFItInfo mfi_info;
905 if (Gpu::notInLaunchRegion()) { mfi_info.SetDynamic(true); }
906#ifdef AMREX_USE_OMP
907#pragma omp parallel if (Gpu::notInLaunchRegion())
908#endif
909 for (MFIter mfi(*(this->m_robin_bcval[amrlev]), mfi_info); mfi.isValid(); ++mfi) {
910 Box const& vbx = mfi.validbox();
911 Array4<RT const> const& ra = robinbc_a->const_array(mfi);
912 Array4<RT const> const& rb = robinbc_b->const_array(mfi);
913 Array4<RT const> const& rf = robinbc_f->const_array(mfi);
914 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
915 const Box& blo = amrex::adjCellLo(vbx, idim);
916 const Box& bhi = amrex::adjCellHi(vbx, idim);
917 bool outside_domain_lo = !(domain.contains(blo));
918 bool outside_domain_hi = !(domain.contains(bhi));
919 if ((!outside_domain_lo) && (!outside_domain_hi)) { continue; }
920 for (int icomp = 0; icomp < ncomp; ++icomp) {
921 Array4<RT> const& rbc = (*(this->m_robin_bcval[amrlev]))[mfi].array(icomp*3);
922 if (this->m_lobc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_lo)
923 {
924 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(blo, i, j, k,
925 {
926 rbc(i,j,k,0) = ra(i,j,k,icomp);
927 rbc(i,j,k,1) = rb(i,j,k,icomp);
928 rbc(i,j,k,2) = rf(i,j,k,icomp);
929 });
930 }
931 if (this->m_hibc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_hi)
932 {
933 AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bhi, i, j, k,
934 {
935 rbc(i,j,k,0) = ra(i,j,k,icomp);
936 rbc(i,j,k,1) = rb(i,j,k,icomp);
937 rbc(i,j,k,2) = rf(i,j,k,icomp);
938 });
939 }
940 }
941 }
942 }
943 }
944}
945
946template <typename MF>
947void
952
953template <typename MF>
954void
955MLCellLinOpT<MF>::updateSolBC (int amrlev, const MF& crse_bcdata) const
956{
957 BL_PROFILE("MLCellLinOp::updateSolBC()");
958
959 AMREX_ALWAYS_ASSERT(amrlev > 0);
960 const int ncomp = this->getNComp();
961 m_crse_sol_br[amrlev]->copyFrom(crse_bcdata, 0, 0, 0, ncomp,
962 this->m_geom[amrlev-1][0].periodicity());
963 m_bndry_sol[amrlev]->updateBndryValues(*m_crse_sol_br[amrlev], 0, 0, ncomp,
964 this->AMRRefRatioVect(amrlev-1),
965 InterpBndryDataT<MF>::IBD_max_order_DEF,
966 m_interpbndry_halfwidth);
967}
968
969template <typename MF>
970void
971MLCellLinOpT<MF>::updateCorBC (int amrlev, const MF& crse_bcdata) const
972{
973 BL_PROFILE("MLCellLinOp::updateCorBC()");
974 AMREX_ALWAYS_ASSERT(amrlev > 0);
975 const int ncomp = this->getNComp();
976 m_crse_cor_br[amrlev]->copyFrom(crse_bcdata, 0, 0, 0, ncomp,
977 this->m_geom[amrlev-1][0].periodicity());
978 m_bndry_cor[amrlev]->updateBndryValues(*m_crse_cor_br[amrlev], 0, 0, ncomp,
979 this->AMRRefRatioVect(amrlev-1),
980 InterpBndryDataT<MF>::IBD_max_order_DEF,
981 m_interpbndry_halfwidth);
982}
983
984template <typename MF>
985void
986MLCellLinOpT<MF>::applyBC (int amrlev, int mglev, MF& in, BCMode bc_mode, StateMode,
987 const MLMGBndryT<MF>* bndry, bool skip_fillboundary) const
988{
989 BL_PROFILE("MLCellLinOp::applyBC()");
990 // No coarsened boundary values, cannot apply inhomog at mglev>0.
991 BL_ASSERT(mglev == 0 || bc_mode == BCMode::Homogeneous);
992 BL_ASSERT(bndry != nullptr || bc_mode == BCMode::Homogeneous);
993
994 const int ncomp = this->getNComp();
995 const int cross = isCrossStencil();
996 const int tensorop = isTensorOp();
997 if (!skip_fillboundary) {
998 in.FillBoundary(0, ncomp, this->m_geom[amrlev][mglev].periodicity(), cross);
999 }
1000
1001 int flagbc = bc_mode == BCMode::Inhomogeneous;
1002 const int imaxorder = this->maxorder;
1003
1004 const Real* dxinv = this->m_geom[amrlev][mglev].InvCellSize();
1005 const RT dxi = static_cast<RT>(dxinv[0]);
1006 const RT dyi = (AMREX_SPACEDIM >= 2) ? static_cast<RT>(dxinv[1]) : RT(1.0);
1007 const RT dzi = (AMREX_SPACEDIM == 3) ? static_cast<RT>(dxinv[2]) : RT(1.0);
1008
1009 const auto& maskvals = m_maskvals[amrlev][mglev];
1010 const auto& bcondloc = *m_bcondloc[amrlev][mglev];
1011
1012 FAB foofab(Box::TheUnitBox(),ncomp);
1013 const auto& foo = foofab.const_array();
1014
1015 MFItInfo mfi_info;
1016 if (Gpu::notInLaunchRegion()) { mfi_info.SetDynamic(true); }
1017
1018 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(cross || tensorop || Gpu::notInLaunchRegion(),
1019 "non-cross stencil not support for gpu");
1020
1021 const int hidden_direction = this->hiddenDirection();
1022
1023#ifdef AMREX_USE_GPU
1024 if ((cross || tensorop) && Gpu::inLaunchRegion())
1025 {
1026 if (! m_bc_tags[amrlev][mglev].is_defined()) {
1027 Vector<MLMGABCTag<RT>> tags;
1028 tags.reserve(in.local_size()*2*AMREX_SPACEDIM*ncomp);
1029 for (MFIter mfi(in); mfi.isValid(); ++mfi) {
1030 const Box& vbx = mfi.validbox();
1031 const auto & bdlv = bcondloc.bndryLocs(mfi);
1032 const auto & bdcv = bcondloc.bndryConds(mfi);
1033
1034 const int local_index = mfi.LocalIndex();
1035
1036 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1037 if (idim != hidden_direction) {
1038 const Orientation olo(idim,Orientation::low);
1039 const Orientation ohi(idim,Orientation::high);
1040 for (int icomp = 0; icomp < ncomp; ++icomp) {
1041 tags.emplace_back(MLMGABCTag<RT>{
1042 maskvals[olo].const_array(mfi),
1043 bdlv[icomp][olo],
1044 amrex::adjCell(vbx,olo),
1045 bdcv[icomp][olo], vbx.length(idim),
1046 icomp, olo, local_index
1047 });
1048 tags.emplace_back(MLMGABCTag<RT>{
1049 maskvals[ohi].const_array(mfi),
1050 bdlv[icomp][ohi],
1051 amrex::adjCell(vbx,ohi),
1052 bdcv[icomp][ohi], vbx.length(idim),
1053 icomp, ohi, local_index
1054 });
1055 }
1056 }
1057 }
1058 }
1059 m_bc_tags[amrlev][mglev].define(tags);
1060 }
1061
1062 MultiArray4<RT const> foo_ma;
1063 Array<MultiArray4<RT const>, 2*AMREX_SPACEDIM> bndry_arrays;
1064 for (OrientationIter oit; oit; ++oit) {
1065 const Orientation ori = oit();
1066 bndry_arrays[ori] = (bndry != nullptr) ?
1067 bndry->bndryValues(ori).arrays() : foo_ma;
1068 }
1069
1070 auto inma = in.arrays();
1071 ParallelFor(m_bc_tags[amrlev][mglev],
1072 [=] AMREX_GPU_DEVICE (int i, int j, int k, MLMGABCTag<RT> const& tag) noexcept
1073 {
1074 const auto& bcval = bndry_arrays[tag.face][tag.local_index];
1075 const int side = tag.face.faceDir();
1076 if (tag.face.coordDir() == 0) {
1077 mllinop_apply_bc_x(side, i, j, k, tag.blen, inma[tag.local_index],
1078 tag.mask, tag.bctype, tag.bcloc, bcval,
1079 imaxorder, dxi, flagbc, tag.comp);
1080 }
1081#if (AMREX_SPACEDIM > 1)
1082 else
1083#if (AMREX_SPACEDIM > 2)
1084 if (tag.face.coordDir() == 1)
1085#endif
1086 {
1087 mllinop_apply_bc_y(side, i, j, k, tag.blen, inma[tag.local_index],
1088 tag.mask, tag.bctype, tag.bcloc, bcval,
1089 imaxorder, dyi, flagbc, tag.comp);
1090 }
1091#if (AMREX_SPACEDIM > 2)
1092 else {
1093 mllinop_apply_bc_z(side, i, j, k, tag.blen, inma[tag.local_index],
1094 tag.mask, tag.bctype, tag.bcloc, bcval,
1095 imaxorder, dzi, flagbc, tag.comp);
1096 }
1097#endif
1098#endif
1099 });
1100 } else
1101#endif
1102 if (cross || tensorop)
1103 {
1104#ifdef AMREX_USE_OMP
1105#pragma omp parallel if (Gpu::notInLaunchRegion())
1106#endif
1107 for (MFIter mfi(in, mfi_info); mfi.isValid(); ++mfi)
1108 {
1109 const Box& vbx = mfi.validbox();
1110 const auto& iofab = in.array(mfi);
1111
1112 const auto & bdlv = bcondloc.bndryLocs(mfi);
1113 const auto & bdcv = bcondloc.bndryConds(mfi);
1114
1115 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
1116 {
1117 if (hidden_direction == idim) { continue; }
1118 const Orientation olo(idim,Orientation::low);
1119 const Orientation ohi(idim,Orientation::high);
1120 const Box blo = amrex::adjCellLo(vbx, idim);
1121 const Box bhi = amrex::adjCellHi(vbx, idim);
1122 const int blen = vbx.length(idim);
1123 const auto& mlo = maskvals[olo].array(mfi);
1124 const auto& mhi = maskvals[ohi].array(mfi);
1125 const auto& bvlo = (bndry != nullptr) ? bndry->bndryValues(olo).const_array(mfi) : foo;
1126 const auto& bvhi = (bndry != nullptr) ? bndry->bndryValues(ohi).const_array(mfi) : foo;
1127 for (int icomp = 0; icomp < ncomp; ++icomp) {
1128 const BoundCond bctlo = bdcv[icomp][olo];
1129 const BoundCond bcthi = bdcv[icomp][ohi];
1130 const RT bcllo = bdlv[icomp][olo];
1131 const RT bclhi = bdlv[icomp][ohi];
1132 if (idim == 0) {
1133 mllinop_apply_bc_x(0, blo, blen, iofab, mlo,
1134 bctlo, bcllo, bvlo,
1135 imaxorder, dxi, flagbc, icomp);
1136 mllinop_apply_bc_x(1, bhi, blen, iofab, mhi,
1137 bcthi, bclhi, bvhi,
1138 imaxorder, dxi, flagbc, icomp);
1139 } else if (idim == 1) {
1140 mllinop_apply_bc_y(0, blo, blen, iofab, mlo,
1141 bctlo, bcllo, bvlo,
1142 imaxorder, dyi, flagbc, icomp);
1143 mllinop_apply_bc_y(1, bhi, blen, iofab, mhi,
1144 bcthi, bclhi, bvhi,
1145 imaxorder, dyi, flagbc, icomp);
1146 } else {
1147 mllinop_apply_bc_z(0, blo, blen, iofab, mlo,
1148 bctlo, bcllo, bvlo,
1149 imaxorder, dzi, flagbc, icomp);
1150 mllinop_apply_bc_z(1, bhi, blen, iofab, mhi,
1151 bcthi, bclhi, bvhi,
1152 imaxorder, dzi, flagbc, icomp);
1153 }
1154 }
1155 }
1156 }
1157 }
1158 else
1159 {
1160#ifdef BL_NO_FORT
1161 amrex::Abort("amrex_mllinop_apply_bc not available when BL_NO_FORT=TRUE");
1162#else
1163 if constexpr (std::is_same_v<Real,RT>) {
1164#ifdef AMREX_USE_OMP
1165#pragma omp parallel
1166#endif
1167 for (MFIter mfi(in, mfi_info); mfi.isValid(); ++mfi)
1168 {
1169 const Box& vbx = mfi.validbox();
1170
1171 const auto & bdlv = bcondloc.bndryLocs(mfi);
1172 const auto & bdcv = bcondloc.bndryConds(mfi);
1173
1174 const RealTuple & bdl = bdlv[0];
1175 const BCTuple & bdc = bdcv[0];
1176
1177 for (OrientationIter oitr; oitr; ++oitr)
1178 {
1179 const Orientation ori = oitr();
1180
1181 int cdr = ori;
1182 RT bcl = bdl[ori];
1183 int bct = bdc[ori];
1184
1185 const auto& fsfab = (bndry != nullptr) ? bndry->bndryValues(ori)[mfi] : foofab;
1186
1187 const Mask& m = maskvals[ori][mfi];
1188
1189 amrex_mllinop_apply_bc(BL_TO_FORTRAN_BOX(vbx),
1190 BL_TO_FORTRAN_ANYD(in[mfi]),
1191 BL_TO_FORTRAN_ANYD(m),
1192 cdr, bct, bcl,
1193 BL_TO_FORTRAN_ANYD(fsfab),
1194 imaxorder, dxinv, flagbc, ncomp, cross);
1195 }
1196 }
1197 } else {
1198 amrex::Abort("Not supported");
1199 }
1200#endif
1201 }
1202}
1203
1204template <typename MF>
1205BoxArray
1206MLCellLinOpT<MF>::makeNGrids (int grid_size) const
1207{
1208 const Box& dombx = this->m_geom[0].back().Domain();
1209
1210 const BoxArray& old_ba = this->m_grids[0].back();
1211 const int N = old_ba.size();
1212 Vector<Box> bv;
1213 bv.reserve(N);
1214 for (int i = 0; i < N; ++i)
1215 {
1216 Box b = old_ba[i];
1217 b.coarsen(grid_size);
1218 b.refine(grid_size);
1219 IntVect sz = b.size();
1220 const IntVect nblks {AMREX_D_DECL(sz[0]/grid_size, sz[1]/grid_size, sz[2]/grid_size)};
1221
1222 IntVect big = b.smallEnd() + grid_size - 1;
1223 b.setBig(big);
1224
1225#if (AMREX_SPACEDIM == 3)
1226 for (int kk = 0; kk < nblks[2]; ++kk) {
1227#endif
1228#if (AMREX_SPACEDIM >= 2)
1229 for (int jj = 0; jj < nblks[1]; ++jj) {
1230#endif
1231 for (int ii = 0; ii < nblks[0]; ++ii)
1232 {
1233 IntVect shft{AMREX_D_DECL(ii*grid_size,jj*grid_size,kk*grid_size)};
1234 Box bb = amrex::shift(b,shft);
1235 bb &= dombx;
1236 bv.push_back(bb);
1237 }
1238#if (AMREX_SPACEDIM >= 2)
1239 }
1240#endif
1241#if (AMREX_SPACEDIM == 3)
1242 }
1243#endif
1244 }
1245
1246 std::sort(bv.begin(), bv.end());
1247 bv.erase(std::unique(bv.begin(), bv.end()), bv.end());
1248
1249 BoxList bl(std::move(bv));
1250
1251 return BoxArray{std::move(bl)};
1252}
1253
1254template <typename MF>
1255void
1256MLCellLinOpT<MF>::restriction (int amrlev, int cmglev, MF& crse, MF& fine) const
1257{
1258 const int ncomp = this->getNComp();
1259 IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[cmglev-1];
1260 // MG levels inside a fine AMR level are not coarsened in the hidden
1261 // direction either.
1262 if (this->hasHiddenDimension()) { ratio[this->hiddenDirection()] = 1; }
1263 amrex::average_down(fine, crse, 0, ncomp, ratio);
1264}
1265
1266template <typename MF>
1267void
1268MLCellLinOpT<MF>::interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const
1269{
1270 const int ncomp = this->getNComp();
1271
1272 Dim3 ratio3 = {.x = 2, .y = 2, .z = 2};
1273 IntVect ratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[fmglev];
1274 if (this->hasHiddenDimension()) { ratio[this->hiddenDirection()] = 1; }
1275 AMREX_D_TERM(ratio3.x = ratio[0];,
1276 ratio3.y = ratio[1];,
1277 ratio3.z = ratio[2];);
1278
1279#ifdef AMREX_USE_GPU
1280 if (Gpu::inLaunchRegion() && fine.isFusingCandidate()) {
1281 auto const& finema = fine.arrays();
1282 auto const& crsema = crse.const_arrays();
1283 ParallelFor(fine, IntVect(0), ncomp,
1284 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
1285 {
1286 int ic = amrex::coarsen(i,ratio3.x);
1287 int jc = amrex::coarsen(j,ratio3.y);
1288 int kc = amrex::coarsen(k,ratio3.z);
1289 finema[box_no](i,j,k,n) += crsema[box_no](ic,jc,kc,n);
1290 });
1291 if (!Gpu::inNoSyncRegion()) {
1292 Gpu::streamSynchronize();
1293 }
1294 } else
1295#endif
1296 {
1297#ifdef AMREX_USE_OMP
1298#pragma omp parallel if (Gpu::notInLaunchRegion())
1299#endif
1300 for (MFIter mfi(fine,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1301 {
1302 const Box& bx = mfi.tilebox();
1303 Array4<RT const> const& cfab = crse.const_array(mfi);
1304 Array4<RT> const& ffab = fine.array(mfi);
1305 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( bx, ncomp, i, j, k, n,
1306 {
1307 int ic = amrex::coarsen(i,ratio3.x);
1308 int jc = amrex::coarsen(j,ratio3.y);
1309 int kc = amrex::coarsen(k,ratio3.z);
1310 ffab(i,j,k,n) += cfab(ic,jc,kc,n);
1311 });
1312 }
1313 }
1314}
1315
1316template <typename MF>
1317void
1318MLCellLinOpT<MF>::interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const
1319{
1320 const int ncomp = this->getNComp();
1321
1322 const Geometry& crse_geom = this->Geom(amrlev,fmglev+1);
1323 IntVect refratio = (amrlev > 0) ? IntVect(2) : this->mg_coarsen_ratio_vec[fmglev];
1324 if (this->hasHiddenDimension()) { refratio[this->hiddenDirection()] = 1; }
1325 const IntVect ng = crse.nGrowVect();
1326
1327 MF cfine;
1328 const MF* cmf;
1329
1330 // Ghost cells not filled from valid cells are zero in both cases.
1331 if (amrex::isMFIterSafe(crse, fine))
1332 {
1333 crse.setBndry(RT(0.0));
1334 crse.FillBoundary(crse_geom.periodicity());
1335 cmf = &crse;
1336 }
1337 else
1338 {
1339 BoxArray cba = fine.boxArray();
1340 cba.coarsen(refratio);
1341 cfine.define(cba, fine.DistributionMap(), ncomp, ng, MFInfo().SetArena(The_Async_Arena()));
1342 cfine.setBndry(RT(0.0));
1343 cfine.ParallelCopy(crse, 0, 0, ncomp, IntVect(0), ng, crse_geom.periodicity());
1344 cmf = & cfine;
1345 }
1346
1347 bool isEB = fine.hasEBFabFactory();
1348 ignore_unused(isEB);
1349
1350#ifdef AMREX_USE_EB
1351 const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(&(fine.Factory()));
1352 const FabArray<EBCellFlagFab>* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr;
1353#endif
1354
1355 MFItInfo mfi_info;
1356 if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
1357#ifdef AMREX_USE_OMP
1358#pragma omp parallel if (Gpu::notInLaunchRegion())
1359#endif
1360 for (MFIter mfi(fine, mfi_info); mfi.isValid(); ++mfi)
1361 {
1362 const Box& bx = mfi.tilebox();
1363 const auto& ff = fine.array(mfi);
1364 const auto& cc = cmf->array(mfi);
1365#ifdef AMREX_USE_EB
1366 bool call_lincc;
1367 if (isEB)
1368 {
1369 const auto& flag = (*flags)[mfi];
1370 if (flag.getType(amrex::grow(bx,1)) == FabType::regular) {
1371 call_lincc = true;
1372 } else {
1373 Array4<EBCellFlag const> const& flg = flag.const_array();
1374 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx, tbx,
1375 {
1376 mlmg_eb_cc_interp_r<2>(tbx, ff, cc, flg, ncomp);
1377 });
1378
1379 call_lincc = false;
1380 }
1381 }
1382 else
1383 {
1384 call_lincc = true;
1385 }
1386#else
1387 const bool call_lincc = true;
1388#endif
1389 if (call_lincc)
1390 {
1391#if (AMREX_SPACEDIM == 3)
1392 if (this->hasHiddenDimension()) {
1393 Box const& bx_2d = this->compactify(bx);
1394 auto const& ff_2d = this->compactify(ff);
1395 auto const& cc_2d = this->compactify(cc);
1396 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx_2d, tbx,
1397 {
1398 TwoD::mlmg_lin_cc_interp_r2(tbx, ff_2d, cc_2d, ncomp);
1399 });
1400 } else
1401#endif
1402 {
1403 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx, tbx,
1404 {
1405 mlmg_lin_cc_interp_r2(tbx, ff, cc, ncomp);
1406 });
1407 }
1408 }
1409 }
1410}
1411
1412template <typename MF>
1413void
1414MLCellLinOpT<MF>::interpolationAmr (int famrlev, MF& fine, const MF& crse,
1415 IntVect const& /*nghost*/) const
1416{
1417 const int ncomp = this->getNComp();
1418 const IntVect refratioV = this->AMRRefRatioVect(famrlev-1);
1419 // AMRRefRatioVect sets hidden direction to 1 and active directions to the
1420 // scalar ratio; max() cleanly recovers the active ratio in all cases.
1421 const int refratio = refratioV.max();
1422
1423#ifdef AMREX_USE_EB
1424 const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(famrlev));
1425 const FabArray<EBCellFlagFab>* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr;
1426#endif
1427
1428 MFItInfo mfi_info;
1429 if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
1430#ifdef AMREX_USE_OMP
1431#pragma omp parallel if (Gpu::notInLaunchRegion())
1432#endif
1433 for (MFIter mfi(fine, mfi_info); mfi.isValid(); ++mfi)
1434 {
1435 const Box& bx = mfi.tilebox();
1436 auto const& ff = fine.array(mfi);
1437 auto const& cc = crse.const_array(mfi);
1438#ifdef AMREX_USE_EB
1439 bool call_lincc;
1440 if (factory)
1441 {
1442 const auto& flag = (*flags)[mfi];
1443 if (flag.getType(amrex::grow(bx,1)) == FabType::regular) {
1444 call_lincc = true;
1445 } else {
1446 Array4<EBCellFlag const> const& flg = flag.const_array();
1447 switch(refratio) {
1448 case 2:
1449 {
1450 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx, tbx,
1451 {
1452 mlmg_eb_cc_interp_r<2>(tbx, ff, cc, flg, ncomp);
1453 });
1454 break;
1455 }
1456 case 4:
1457 {
1458 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx, tbx,
1459 {
1460 mlmg_eb_cc_interp_r<4>(tbx, ff, cc, flg, ncomp);
1461 });
1462 break;
1463 }
1464 default:
1465 amrex::Abort("mlmg_eb_cc_interp: only refratio 2 and 4 are supported");
1466 }
1467
1468 call_lincc = false;
1469 }
1470 }
1471 else
1472 {
1473 call_lincc = true;
1474 }
1475#else
1476 const bool call_lincc = true;
1477#endif
1478 if (call_lincc)
1479 {
1480#if (AMREX_SPACEDIM == 3)
1481 if (this->hasHiddenDimension()) {
1482 Box const& bx_2d = this->compactify(bx);
1483 auto const& ff_2d = this->compactify(ff);
1484 auto const& cc_2d = this->compactify(cc);
1485 switch(refratio) {
1486 case 2:
1487 {
1488 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx_2d, tbx,
1489 {
1490 TwoD::mlmg_lin_cc_interp_r2(tbx, ff_2d, cc_2d, ncomp);
1491 });
1492 break;
1493 }
1494 case 4:
1495 {
1496 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx_2d, tbx,
1497 {
1498 TwoD::mlmg_lin_cc_interp_r4(tbx, ff_2d, cc_2d, ncomp);
1499 });
1500 break;
1501 }
1502 default:
1503 amrex::Abort("mlmg_lin_cc_interp: only refratio 2 and 4 are supported");
1504 }
1505 } else
1506#endif
1507 {
1508 switch(refratio) {
1509 case 2:
1510 {
1511 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx, tbx,
1512 {
1513 mlmg_lin_cc_interp_r2(tbx, ff, cc, ncomp);
1514 });
1515 break;
1516 }
1517 case 4:
1518 {
1519 AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE (bx, tbx,
1520 {
1521 mlmg_lin_cc_interp_r4(tbx, ff, cc, ncomp);
1522 });
1523 break;
1524 }
1525 default:
1526 amrex::Abort("mlmg_lin_cc_interp: only refratio 2 and 4 are supported");
1527 }
1528 }
1529 }
1530 }
1531}
1532
1533template <typename MF>
1534void
1535MLCellLinOpT<MF>::averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs,
1536 const MF& fine_sol, const MF& fine_rhs)
1537{
1538 const auto amrrr = this->AMRRefRatioVect(camrlev);
1539 const int ncomp = this->getNComp();
1540 amrex::average_down(fine_sol, crse_sol, 0, ncomp, amrrr);
1541 amrex::average_down(fine_rhs, crse_rhs, 0, ncomp, amrrr);
1542}
1543
1544template <typename MF>
1545void
1546MLCellLinOpT<MF>::apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode,
1547 StateMode s_mode, const MLMGBndryT<MF>* bndry) const
1548{
1549 BL_PROFILE("MLCellLinOp::apply()");
1550 applyBC(amrlev, mglev, in, bc_mode, s_mode, bndry);
1551 Fapply(amrlev, mglev, out, in);
1552}
1553
1554template <typename MF>
1555void
1556MLCellLinOpT<MF>::smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
1557 bool skip_fillboundary, int niter) const
1558{
1559 BL_PROFILE("MLCellLinOp::smooth()");
1560 for (int i = 0; i < niter; ++i) {
1561 for (int redblack = 0; redblack < 2; ++redblack)
1562 {
1563 applyBC(amrlev, mglev, sol, BCMode::Homogeneous, StateMode::Solution,
1564 nullptr, skip_fillboundary);
1565 Fsmooth(amrlev, mglev, sol, rhs, redblack);
1566 skip_fillboundary = false;
1567 }
1568 }
1569}
1570
1571template <typename MF>
1572void
1573MLCellLinOpT<MF>::solutionResidual (int amrlev, MF& resid, MF& x, const MF& b,
1574 const MF* crse_bcdata)
1575{
1576 BL_PROFILE("MLCellLinOp::solutionResidual()");
1577 const int ncomp = this->getNComp();
1578 if (crse_bcdata != nullptr) {
1579 updateSolBC(amrlev, *crse_bcdata);
1580 }
1581 const int mglev = 0;
1582 apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Solution,
1583 m_bndry_sol[amrlev].get());
1584
1585 AMREX_ASSERT(resid.nComp() == b.nComp());
1586 MF::Xpay(resid, RT(-1.0), b, 0, 0, ncomp, IntVect(0));
1587}
1588
1589template <typename MF>
1590void
1591MLCellLinOpT<MF>::prepareForFluxes (int amrlev, const MF* crse_bcdata)
1592{
1593 if (crse_bcdata != nullptr) {
1594 updateSolBC(amrlev, *crse_bcdata);
1595 }
1596}
1597
1598template <typename MF>
1599void
1600MLCellLinOpT<MF>::correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b,
1601 BCMode bc_mode, const MF* crse_bcdata)
1602{
1603 BL_PROFILE("MLCellLinOp::correctionResidual()");
1604 const int ncomp = this->getNComp();
1605 if (bc_mode == BCMode::Inhomogeneous)
1606 {
1607 if (crse_bcdata)
1608 {
1609 AMREX_ASSERT(mglev == 0 && amrlev > 0);
1610 updateCorBC(amrlev, *crse_bcdata);
1611 }
1612 apply(amrlev, mglev, resid, x, BCMode::Inhomogeneous, StateMode::Correction,
1613 m_bndry_cor[amrlev].get());
1614 }
1615 else
1616 {
1617 AMREX_ASSERT(crse_bcdata == nullptr);
1618 apply(amrlev, mglev, resid, x, BCMode::Homogeneous, StateMode::Correction, nullptr);
1619 }
1620
1621 MF::Xpay(resid, Real(-1.0), b, 0, 0, ncomp, IntVect(0));
1622}
1623
1624template <typename MF>
1625void
1626MLCellLinOpT<MF>::reflux (int crse_amrlev, MF& res, const MF& crse_sol, const MF&,
1627 MF&, MF& fine_sol, const MF&) const
1628{
1629 BL_PROFILE("MLCellLinOp::reflux()");
1630
1631 auto& fluxreg = m_fluxreg[crse_amrlev];
1632 fluxreg.reset();
1633
1634 const int ncomp = this->getNComp();
1635
1636 const int fine_amrlev = crse_amrlev+1;
1637
1638 Real dt = Real(1.0);
1639 const Real* crse_dx = this->m_geom[crse_amrlev][0].CellSize();
1640 const Real* fine_dx = this->m_geom[fine_amrlev][0].CellSize();
1641
1642 const int mglev = 0;
1643 applyBC(fine_amrlev, mglev, fine_sol, BCMode::Inhomogeneous, StateMode::Solution,
1644 m_bndry_sol[fine_amrlev].get());
1645
1646 MFItInfo mfi_info;
1647 if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
1648
1649#ifdef AMREX_USE_OMP
1650#pragma omp parallel if (Gpu::notInLaunchRegion())
1651#endif
1652 {
1653 Array<FAB,AMREX_SPACEDIM> flux;
1654 Array<FAB*,AMREX_SPACEDIM> pflux {{ AMREX_D_DECL(flux.data(), flux.data()+1, flux.data()+2) }};
1655 Array<FAB const*,AMREX_SPACEDIM> cpflux {{ AMREX_D_DECL(flux.data(), flux.data()+1, flux.data()+2) }};
1656
1657 for (MFIter mfi(crse_sol, mfi_info); mfi.isValid(); ++mfi)
1658 {
1659 if (fluxreg.CrseHasWork(mfi))
1660 {
1661 const Box& tbx = mfi.tilebox();
1662 AMREX_D_TERM(flux[0].resize(amrex::surroundingNodes(tbx,0),ncomp,The_Async_Arena());,
1663 flux[1].resize(amrex::surroundingNodes(tbx,1),ncomp,The_Async_Arena());,
1664 flux[2].resize(amrex::surroundingNodes(tbx,2),ncomp,The_Async_Arena()););
1665 FFlux(crse_amrlev, mfi, pflux, crse_sol[mfi], Location::FaceCentroid);
1666 fluxreg.CrseAdd(mfi, cpflux, crse_dx, dt, RunOn::Gpu);
1667 }
1668 }
1669
1670#ifdef AMREX_USE_OMP
1671#pragma omp barrier
1672#endif
1673
1674 for (MFIter mfi(fine_sol, mfi_info); mfi.isValid(); ++mfi)
1675 {
1676 if (fluxreg.FineHasWork(mfi))
1677 {
1678 const Box& tbx = mfi.tilebox();
1679 const int face_only = true;
1680 AMREX_D_TERM(flux[0].resize(amrex::surroundingNodes(tbx,0),ncomp,The_Async_Arena());,
1681 flux[1].resize(amrex::surroundingNodes(tbx,1),ncomp,The_Async_Arena());,
1682 flux[2].resize(amrex::surroundingNodes(tbx,2),ncomp,The_Async_Arena()););
1683 FFlux(fine_amrlev, mfi, pflux, fine_sol[mfi], Location::FaceCentroid, face_only);
1684 fluxreg.FineAdd(mfi, cpflux, fine_dx, dt, RunOn::Gpu);
1685 }
1686 }
1687 }
1688
1689 fluxreg.Reflux(res);
1690 this->applyOverset(crse_amrlev, res);
1691}
1692
1693template <typename MF>
1694void
1695MLCellLinOpT<MF>::compFlux (int amrlev, const Array<MF*,AMREX_SPACEDIM>& fluxes,
1696 MF& sol, Location loc) const
1697{
1698 BL_PROFILE("MLCellLinOp::compFlux()");
1699
1700 const int mglev = 0;
1701 const int ncomp = this->getNComp();
1702 applyBC(amrlev, mglev, sol, BCMode::Inhomogeneous, StateMode::Solution,
1703 m_bndry_sol[amrlev].get());
1704
1705 MFItInfo mfi_info;
1706 if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
1707
1708#ifdef AMREX_USE_OMP
1709#pragma omp parallel if (Gpu::notInLaunchRegion())
1710#endif
1711 {
1712 Array<FAB,AMREX_SPACEDIM> flux;
1713 Array<FAB*,AMREX_SPACEDIM> pflux {{ AMREX_D_DECL(flux.data(), flux.data()+1, flux.data()+2) }};
1714 for (MFIter mfi(sol, mfi_info); mfi.isValid(); ++mfi)
1715 {
1716 const Box& tbx = mfi.tilebox();
1717 AMREX_D_TERM(flux[0].resize(amrex::surroundingNodes(tbx,0),ncomp,The_Async_Arena());,
1718 flux[1].resize(amrex::surroundingNodes(tbx,1),ncomp,The_Async_Arena());,
1719 flux[2].resize(amrex::surroundingNodes(tbx,2),ncomp,The_Async_Arena()););
1720 FFlux(amrlev, mfi, pflux, sol[mfi], loc);
1721 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1722 const Box& nbx = mfi.nodaltilebox(idim);
1723 auto const& dst = fluxes[idim]->array(mfi);
1724 auto const& src = pflux[idim]->const_array();
1725 AMREX_HOST_DEVICE_PARALLEL_FOR_4D (nbx, ncomp, i, j, k, n,
1726 {
1727 dst(i,j,k,n) = src(i,j,k,n);
1728 });
1729 }
1730 }
1731 }
1732}
1733
1734template <typename MF>
1735void
1736MLCellLinOpT<MF>::compGrad (int amrlev, const Array<MF*,AMREX_SPACEDIM>& grad,
1737 MF& sol, Location /*loc*/) const
1738{
1739 BL_PROFILE("MLCellLinOp::compGrad()");
1740
1741 if (sol.nComp() > 1) {
1742 amrex::Abort("MLCellLinOp::compGrad called, but only works for single-component solves");
1743 }
1744
1745 const int mglev = 0;
1746 applyBC(amrlev, mglev, sol, BCMode::Inhomogeneous, StateMode::Solution,
1747 m_bndry_sol[amrlev].get());
1748
1749 const int ncomp = this->getNComp();
1750
1751 AMREX_D_TERM(const RT dxi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(0));,
1752 const RT dyi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(1));,
1753 const RT dzi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(2)););
1754#ifdef AMREX_USE_OMP
1755#pragma omp parallel if (Gpu::notInLaunchRegion())
1756#endif
1757 for (MFIter mfi(sol, TilingIfNotGPU()); mfi.isValid(); ++mfi)
1758 {
1759 AMREX_D_TERM(const Box& xbx = mfi.nodaltilebox(0);,
1760 const Box& ybx = mfi.nodaltilebox(1);,
1761 const Box& zbx = mfi.nodaltilebox(2););
1762 const auto& s = sol.array(mfi);
1763 AMREX_D_TERM(const auto& gx = grad[0]->array(mfi);,
1764 const auto& gy = grad[1]->array(mfi);,
1765 const auto& gz = grad[2]->array(mfi););
1766
1767 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( xbx, ncomp, i, j, k, n,
1768 {
1769 gx(i,j,k,n) = dxi*(s(i,j,k,n) - s(i-1,j,k,n));
1770 });
1771#if (AMREX_SPACEDIM >= 2)
1772 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( ybx, ncomp, i, j, k, n,
1773 {
1774 gy(i,j,k,n) = dyi*(s(i,j,k,n) - s(i,j-1,k,n));
1775 });
1776#endif
1777#if (AMREX_SPACEDIM == 3)
1778 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( zbx, ncomp, i, j, k, n,
1779 {
1780 gz(i,j,k,n) = dzi*(s(i,j,k,n) - s(i,j,k-1,n));
1781 });
1782#endif
1783 }
1784
1785 addInhomogNeumannFlux(amrlev, grad, sol, false);
1786}
1787
1788template <typename MF>
1789void
1790MLCellLinOpT<MF>::applyMetricTerm (int amrlev, int mglev, MF& rhs) const
1791{
1792 amrex::ignore_unused(amrlev,mglev,rhs);
1793#if (AMREX_SPACEDIM != 3)
1794 if (!m_has_metric_term) { return; }
1795
1796 const int ncomp = rhs.nComp();
1797
1798 bool cc = rhs.ixType().cellCentered(0);
1799
1800 const Geometry& geom = this->m_geom[amrlev][mglev];
1801 const RT dx = static_cast<RT>(geom.CellSize(0));
1802 const RT probxlo = static_cast<RT>(geom.ProbLo(0));
1803
1804#ifdef AMREX_USE_OMP
1805#pragma omp parallel if (Gpu::notInLaunchRegion())
1806#endif
1807 for (MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1808 {
1809 const Box& tbx = mfi.tilebox();
1810 auto const& rhsarr = rhs.array(mfi);
1811#if (AMREX_SPACEDIM == 1)
1812 if (cc) {
1813 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1814 {
1815 RT rc = probxlo + (RT(i)+RT(0.5))*dx;
1816 rhsarr(i,j,k,n) *= rc*rc;
1817 });
1818 } else {
1819 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1820 {
1821 RT re = probxlo + RT(i)*dx;
1822 rhsarr(i,j,k,n) *= re*re;
1823 });
1824 }
1825#elif (AMREX_SPACEDIM == 2)
1826 if (cc) {
1827 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1828 {
1829 RT rc = probxlo + (RT(i)+RT(0.5))*dx;
1830 rhsarr(i,j,k,n) *= rc;
1831 });
1832 } else {
1833 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1834 {
1835 RT re = probxlo + RT(i)*dx;
1836 rhsarr(i,j,k,n) *= re;
1837 });
1838 }
1839#endif
1840 }
1841#endif
1842}
1843
1844template <typename MF>
1845void
1846MLCellLinOpT<MF>::unapplyMetricTerm (int amrlev, int mglev, MF& rhs) const
1847{
1848 amrex::ignore_unused(amrlev,mglev,rhs);
1849#if (AMREX_SPACEDIM != 3)
1850 if (!m_has_metric_term) { return; }
1851
1852 const int ncomp = rhs.nComp();
1853
1854 bool cc = rhs.ixType().cellCentered(0);
1855
1856 const Geometry& geom = this->m_geom[amrlev][mglev];
1857 const RT dx = static_cast<RT>(geom.CellSize(0));
1858 const RT probxlo = static_cast<RT>(geom.ProbLo(0));
1859
1860#ifdef AMREX_USE_OMP
1861#pragma omp parallel if (Gpu::notInLaunchRegion())
1862#endif
1863 for (MFIter mfi(rhs,TilingIfNotGPU()); mfi.isValid(); ++mfi)
1864 {
1865 const Box& tbx = mfi.tilebox();
1866 auto const& rhsarr = rhs.array(mfi);
1867#if (AMREX_SPACEDIM == 1)
1868 if (cc) {
1869 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1870 {
1871 RT rcinv = RT(1.0)/(probxlo + (RT(i)+RT(0.5))*dx);
1872 rhsarr(i,j,k,n) *= rcinv*rcinv;
1873 });
1874 } else {
1875 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1876 {
1877 RT re = probxlo + RT(i)*dx;
1878 RT reinv = (re==RT(0.0)) ? RT(0.0) : RT(1.)/re;
1879 rhsarr(i,j,k,n) *= reinv*reinv;
1880 });
1881 }
1882#elif (AMREX_SPACEDIM == 2)
1883 if (cc) {
1884 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1885 {
1886 RT rcinv = RT(1.0)/(probxlo + (RT(i)+RT(0.5))*dx);
1887 rhsarr(i,j,k,n) *= rcinv;
1888 });
1889 } else {
1890 AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n,
1891 {
1892 RT re = probxlo + RT(i)*dx;
1893 RT reinv = (re==RT(0.0)) ? RT(0.0) : RT(1.)/re;
1894 rhsarr(i,j,k,n) *= reinv;
1895 });
1896 }
1897#endif
1898 }
1899#endif
1900}
1901
1902template <typename MF>
1903auto
1904MLCellLinOpT<MF>::getSolvabilityOffset (int amrlev, int mglev, MF const& rhs) const
1905 -> Vector<RT>
1906{
1907 computeVolInv();
1908
1909 const int ncomp = this->getNComp();
1910 Vector<RT> offset(ncomp);
1911
1912#ifdef AMREX_USE_EB
1913 const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(amrlev,mglev));
1914 if (factory && !factory->isAllRegular())
1915 {
1916 if constexpr (std::is_same<MF,MultiFab>()) {
1917 const MultiFab& vfrac = factory->getVolFrac();
1918 for (int c = 0; c < ncomp; ++c) {
1919 offset[c] = amrex::Dot(rhs, c, vfrac, 0, 1, IntVect(0), true)
1920 * m_volinv[amrlev][mglev];
1921 }
1922 } else {
1923 amrex::Abort("TODO: MLMG with EB only works with MultiFab");
1924 }
1925 }
1926 else
1927#endif
1928 {
1929 for (int c = 0; c < ncomp; ++c) {
1930 offset[c] = rhs.sum(c,IntVect(0),true) * m_volinv[amrlev][mglev];
1931 }
1932 }
1933
1934 ParallelAllReduce::Sum(offset.data(), ncomp, ParallelContext::CommunicatorSub());
1935
1936 return offset;
1937}
1938
1939template <typename MF>
1940void
1941MLCellLinOpT<MF>::fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, MF& rhs,
1942 Vector<RT> const& offset) const
1943{
1944 const int ncomp = this->getNComp();
1945 for (int c = 0; c < ncomp; ++c) {
1946 rhs.plus(-offset[c], c, 1);
1947 }
1948#ifdef AMREX_USE_EB
1949 if (!rhs.isAllRegular()) {
1950 if constexpr (std::is_same<MF,MultiFab>()) {
1951 amrex::EB_set_covered(rhs, 0, ncomp, 0, 0.0_rt);
1952 } else {
1953 amrex::Abort("amrex::EB_set_covered only works with MultiFab");
1954 }
1955 }
1956#endif
1957}
1958
1959template <typename MF>
1960void
1961MLCellLinOpT<MF>::prepareForSolve ()
1962{
1963 BL_PROFILE("MLCellLinOp::prepareForSolve()");
1964
1965 const int imaxorder = this->maxorder;
1966 const int ncomp = this->getNComp();
1967 const int hidden_direction = this->hiddenDirection();
1968 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev)
1969 {
1970 for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev)
1971 {
1972 const auto& bcondloc = *m_bcondloc[amrlev][mglev];
1973 const auto& maskvals = m_maskvals[amrlev][mglev];
1974
1975 const RT dxi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(0));
1976 const RT dyi = static_cast<RT>((AMREX_SPACEDIM >= 2) ? this->m_geom[amrlev][mglev].InvCellSize(1) : Real(1.0));
1977 const RT dzi = static_cast<RT>((AMREX_SPACEDIM == 3) ? this->m_geom[amrlev][mglev].InvCellSize(2) : Real(1.0));
1978
1979 auto& undrrelxr = this->m_undrrelxr[amrlev][mglev];
1980 MF foo(this->m_grids[amrlev][mglev], this->m_dmap[amrlev][mglev], ncomp, 0, MFInfo().SetAlloc(false));
1981
1982#ifdef AMREX_USE_EB
1983 const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->m_factory[amrlev][mglev].get());
1984 const FabArray<EBCellFlagFab>* flags =
1985 (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr;
1986 auto area = (factory) ? factory->getAreaFrac()
1987 : Array<const MultiCutFab*,AMREX_SPACEDIM>{AMREX_D_DECL(nullptr,nullptr,nullptr)};
1988 amrex::ignore_unused(area);
1989#endif
1990
1991#ifdef AMREX_USE_GPU
1992 if (Gpu::inLaunchRegion()) {
1993#ifdef AMREX_USE_EB
1994 if (factory && !factory->isAllRegular()) {
1995#if defined(AMREX_USE_CUDA) && defined(_WIN32)
1996 if (!std::is_same<MF,MultiFab>()) {
1997#else
1998 if constexpr (!std::is_same<MF,MultiFab>()) {
1999#endif
2000 amrex::Abort("MLCellLinOp with EB only works with MultiFab");
2001 } else {
2002 Vector<MLMGPSEBTag<RT>> tags;
2003 tags.reserve(foo.local_size()*AMREX_SPACEDIM*ncomp);
2004
2005 for (MFIter mfi(foo); mfi.isValid(); ++mfi)
2006 {
2007 const Box& vbx = mfi.validbox();
2008
2009 const auto & bdlv = bcondloc.bndryLocs(mfi);
2010 const auto & bdcv = bcondloc.bndryConds(mfi);
2011
2012 auto fabtyp = (flags) ? (*flags)[mfi].getType(vbx) : FabType::regular;
2013
2014 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
2015 {
2016 if (idim != hidden_direction && fabtyp != FabType::covered) {
2017 const Orientation olo(idim,Orientation::low);
2018 const Orientation ohi(idim,Orientation::high);
2019 auto const& ap = (fabtyp == FabType::singlevalued)
2020 ? area[idim]->const_array(mfi) : Array4<Real const>{};
2021 for (int icomp = 0; icomp < ncomp; ++icomp) {
2022 tags.emplace_back(MLMGPSEBTag<RT>{undrrelxr[olo].array(mfi),
2023 undrrelxr[ohi].array(mfi),
2024 ap,
2025 maskvals[olo].const_array(mfi),
2026 maskvals[ohi].const_array(mfi),
2027 bdlv[icomp][olo], bdlv[icomp][ohi],
2028 amrex::adjCell(vbx,olo),
2029 bdcv[icomp][olo], bdcv[icomp][ohi],
2030 vbx.length(idim), icomp, idim});
2031 }
2032 }
2033 }
2034 }
2035
2036 ParallelFor(tags,
2037 [=] AMREX_GPU_DEVICE (int i, int j, int k, MLMGPSEBTag<RT> const& tag) noexcept
2038 {
2039 if (tag.ap) {
2040 if (tag.dir == 0)
2041 {
2042 mllinop_comp_interp_coef0_x_eb
2043 (0, i , j, k, tag.blen, tag.flo, tag.mlo, tag.ap,
2044 tag.bctlo, tag.bcllo, imaxorder, dxi, tag.comp);
2045 mllinop_comp_interp_coef0_x_eb
2046 (1, i+tag.blen+1, j, k, tag.blen, tag.fhi, tag.mhi, tag.ap,
2047 tag.bcthi, tag.bclhi, imaxorder, dxi, tag.comp);
2048 }
2049#if (AMREX_SPACEDIM > 1)
2050 else
2051#if (AMREX_SPACEDIM > 2)
2052 if (tag.dir == 1)
2053#endif
2054 {
2055 mllinop_comp_interp_coef0_y_eb
2056 (0, i, j , k, tag.blen, tag.flo, tag.mlo, tag.ap,
2057 tag.bctlo, tag.bcllo, imaxorder, dyi, tag.comp);
2058 mllinop_comp_interp_coef0_y_eb
2059 (1, i, j+tag.blen+1, k, tag.blen, tag.fhi, tag.mhi, tag.ap,
2060 tag.bcthi, tag.bclhi, imaxorder, dyi, tag.comp);
2061 }
2062#if (AMREX_SPACEDIM > 2)
2063 else {
2064 mllinop_comp_interp_coef0_z_eb
2065 (0, i, j, k , tag.blen, tag.flo, tag.mlo, tag.ap,
2066 tag.bctlo, tag.bcllo, imaxorder, dzi, tag.comp);
2067 mllinop_comp_interp_coef0_z_eb
2068 (1, i, j, k+tag.blen+1, tag.blen, tag.fhi, tag.mhi, tag.ap,
2069 tag.bcthi, tag.bclhi, imaxorder, dzi, tag.comp);
2070 }
2071#endif
2072#endif
2073 } else {
2074 if (tag.dir == 0)
2075 {
2076 mllinop_comp_interp_coef0_x
2077 (0, i , j, k, tag.blen, tag.flo, tag.mlo,
2078 tag.bctlo, tag.bcllo, imaxorder, dxi, tag.comp);
2079 mllinop_comp_interp_coef0_x
2080 (1, i+tag.blen+1, j, k, tag.blen, tag.fhi, tag.mhi,
2081 tag.bcthi, tag.bclhi, imaxorder, dxi, tag.comp);
2082 }
2083#if (AMREX_SPACEDIM > 1)
2084 else
2085#if (AMREX_SPACEDIM > 2)
2086 if (tag.dir == 1)
2087#endif
2088 {
2089 mllinop_comp_interp_coef0_y
2090 (0, i, j , k, tag.blen, tag.flo, tag.mlo,
2091 tag.bctlo, tag.bcllo, imaxorder, dyi, tag.comp);
2092 mllinop_comp_interp_coef0_y
2093 (1, i, j+tag.blen+1, k, tag.blen, tag.fhi, tag.mhi,
2094 tag.bcthi, tag.bclhi, imaxorder, dyi, tag.comp);
2095 }
2096#if (AMREX_SPACEDIM > 2)
2097 else {
2098 mllinop_comp_interp_coef0_z
2099 (0, i, j, k , tag.blen, tag.flo, tag.mlo,
2100 tag.bctlo, tag.bcllo, imaxorder, dzi, tag.comp);
2101 mllinop_comp_interp_coef0_z
2102 (1, i, j, k+tag.blen+1, tag.blen, tag.fhi, tag.mhi,
2103 tag.bcthi, tag.bclhi, imaxorder, dzi, tag.comp);
2104 }
2105#endif
2106#endif
2107 }
2108 });
2109 }
2110 } else
2111#endif
2112 {
2113 Vector<MLMGPSTag<RT>> tags;
2114 tags.reserve(foo.local_size()*AMREX_SPACEDIM*ncomp);
2115
2116 for (MFIter mfi(foo); mfi.isValid(); ++mfi)
2117 {
2118 const Box& vbx = mfi.validbox();
2119
2120 const auto & bdlv = bcondloc.bndryLocs(mfi);
2121 const auto & bdcv = bcondloc.bndryConds(mfi);
2122
2123 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
2124 {
2125 if (idim != hidden_direction) {
2126 const Orientation olo(idim,Orientation::low);
2127 const Orientation ohi(idim,Orientation::high);
2128 for (int icomp = 0; icomp < ncomp; ++icomp) {
2129 tags.emplace_back(MLMGPSTag<RT>{undrrelxr[olo].array(mfi),
2130 undrrelxr[ohi].array(mfi),
2131 maskvals[olo].const_array(mfi),
2132 maskvals[ohi].const_array(mfi),
2133 bdlv[icomp][olo], bdlv[icomp][ohi],
2134 amrex::adjCell(vbx,olo),
2135 bdcv[icomp][olo], bdcv[icomp][ohi],
2136 vbx.length(idim), icomp, idim});
2137 }
2138 }
2139 }
2140 }
2141
2142 ParallelFor(tags,
2143 [=] AMREX_GPU_DEVICE (int i, int j, int k, MLMGPSTag<RT> const& tag) noexcept
2144 {
2145 if (tag.dir == 0)
2146 {
2147 mllinop_comp_interp_coef0_x
2148 (0, i , j, k, tag.blen, tag.flo, tag.mlo,
2149 tag.bctlo, tag.bcllo, imaxorder, dxi, tag.comp);
2150 mllinop_comp_interp_coef0_x
2151 (1, i+tag.blen+1, j, k, tag.blen, tag.fhi, tag.mhi,
2152 tag.bcthi, tag.bclhi, imaxorder, dxi, tag.comp);
2153 }
2154#if (AMREX_SPACEDIM > 1)
2155 else
2156#if (AMREX_SPACEDIM > 2)
2157 if (tag.dir == 1)
2158#endif
2159 {
2160 mllinop_comp_interp_coef0_y
2161 (0, i, j , k, tag.blen, tag.flo, tag.mlo,
2162 tag.bctlo, tag.bcllo, imaxorder, dyi, tag.comp);
2163 mllinop_comp_interp_coef0_y
2164 (1, i, j+tag.blen+1, k, tag.blen, tag.fhi, tag.mhi,
2165 tag.bcthi, tag.bclhi, imaxorder, dyi, tag.comp);
2166 }
2167#if (AMREX_SPACEDIM > 2)
2168 else {
2169 mllinop_comp_interp_coef0_z
2170 (0, i, j, k , tag.blen, tag.flo, tag.mlo,
2171 tag.bctlo, tag.bcllo, imaxorder, dzi, tag.comp);
2172 mllinop_comp_interp_coef0_z
2173 (1, i, j, k+tag.blen+1, tag.blen, tag.fhi, tag.mhi,
2174 tag.bcthi, tag.bclhi, imaxorder, dzi, tag.comp);
2175 }
2176#endif
2177#endif
2178 });
2179 }
2180 } else
2181#endif
2182 {
2183#ifdef AMREX_USE_OMP
2184#pragma omp parallel
2185#endif
2186 for (MFIter mfi(foo, MFItInfo{}.SetDynamic(true)); mfi.isValid(); ++mfi)
2187 {
2188 const Box& vbx = mfi.validbox();
2189
2190 const auto & bdlv = bcondloc.bndryLocs(mfi);
2191 const auto & bdcv = bcondloc.bndryConds(mfi);
2192
2193#ifdef AMREX_USE_EB
2194 auto fabtyp = (flags) ? (*flags)[mfi].getType(vbx) : FabType::regular;
2195#endif
2196 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
2197 {
2198 if (idim == hidden_direction) { continue; }
2199 const Orientation olo(idim,Orientation::low);
2200 const Orientation ohi(idim,Orientation::high);
2201 const Box blo = amrex::adjCellLo(vbx, idim);
2202 const Box bhi = amrex::adjCellHi(vbx, idim);
2203 const int blen = vbx.length(idim);
2204 const auto& mlo = maskvals[olo].array(mfi);
2205 const auto& mhi = maskvals[ohi].array(mfi);
2206 const auto& flo = undrrelxr[olo].array(mfi);
2207 const auto& fhi = undrrelxr[ohi].array(mfi);
2208 for (int icomp = 0; icomp < ncomp; ++icomp) {
2209 const BoundCond bctlo = bdcv[icomp][olo];
2210 const BoundCond bcthi = bdcv[icomp][ohi];
2211 const auto bcllo = bdlv[icomp][olo];
2212 const auto bclhi = bdlv[icomp][ohi];
2213#ifdef AMREX_USE_EB
2214 if (fabtyp == FabType::singlevalued) {
2215 if constexpr (!std::is_same<MF,MultiFab>()) {
2216 amrex::Abort("MLCellLinOp with EB only works with MultiFab");
2217 } else {
2218 auto const& ap = area[idim]->const_array(mfi);
2219 if (idim == 0) {
2220 mllinop_comp_interp_coef0_x_eb
2221 (0, blo, blen, flo, mlo, ap, bctlo, bcllo,
2222 imaxorder, dxi, icomp);
2223 mllinop_comp_interp_coef0_x_eb
2224 (1, bhi, blen, fhi, mhi, ap, bcthi, bclhi,
2225 imaxorder, dxi, icomp);
2226 } else if (idim == 1) {
2227 mllinop_comp_interp_coef0_y_eb
2228 (0, blo, blen, flo, mlo, ap, bctlo, bcllo,
2229 imaxorder, dyi, icomp);
2230 mllinop_comp_interp_coef0_y_eb
2231 (1, bhi, blen, fhi, mhi, ap, bcthi, bclhi,
2232 imaxorder, dyi, icomp);
2233 } else {
2234 mllinop_comp_interp_coef0_z_eb
2235 (0, blo, blen, flo, mlo, ap, bctlo, bcllo,
2236 imaxorder, dzi, icomp);
2237 mllinop_comp_interp_coef0_z_eb
2238 (1, bhi, blen, fhi, mhi, ap, bcthi, bclhi,
2239 imaxorder, dzi, icomp);
2240 }
2241 }
2242 } else if (fabtyp == FabType::regular)
2243#endif
2244 {
2245 if (idim == 0) {
2246 mllinop_comp_interp_coef0_x
2247 (0, blo, blen, flo, mlo, bctlo, bcllo,
2248 imaxorder, dxi, icomp);
2249 mllinop_comp_interp_coef0_x
2250 (1, bhi, blen, fhi, mhi, bcthi, bclhi,
2251 imaxorder, dxi, icomp);
2252 } else if (idim == 1) {
2253 mllinop_comp_interp_coef0_y
2254 (0, blo, blen, flo, mlo, bctlo, bcllo,
2255 imaxorder, dyi, icomp);
2256 mllinop_comp_interp_coef0_y
2257 (1, bhi, blen, fhi, mhi, bcthi, bclhi,
2258 imaxorder, dyi, icomp);
2259 } else {
2260 mllinop_comp_interp_coef0_z
2261 (0, blo, blen, flo, mlo, bctlo, bcllo,
2262 imaxorder, dzi, icomp);
2263 mllinop_comp_interp_coef0_z
2264 (1, bhi, blen, fhi, mhi, bcthi, bclhi,
2265 imaxorder, dzi, icomp);
2266 }
2267 }
2268 }
2269 }
2270 }
2271 }
2272 }
2273 }
2274}
2275
2276template <typename MF>
2277auto
2278MLCellLinOpT<MF>::xdoty (int /*amrlev*/, int /*mglev*/, const MF& x, const MF& y, bool local) const
2279 -> RT
2280{
2281 const int ncomp = this->getNComp();
2282 const IntVect nghost(0);
2283 RT result = amrex::Dot(x,0,y,0,ncomp,nghost,true);
2284 if (!local) {
2285 ParallelAllReduce::Sum(result, ParallelContext::CommunicatorSub());
2286 }
2287 return result;
2288}
2289
2290template <typename MF>
2291auto
2292MLCellLinOpT<MF>::dotProductPrecond (Vector<MF const*> const& x, Vector<MF const*> const& y) const -> RT
2293{
2294 const int ncomp = this->getNComp();
2295 const IntVect nghost(0);
2296 RT result = 0;
2297 for (int ilev = 0; ilev < this->NAMRLevels()-1; ++ilev) {
2298 result += amrex::Dot(*m_norm_fine_mask[ilev], *x[ilev], 0, *y[ilev], 0, ncomp, nghost, true);
2299 }
2300 result += amrex::Dot(*x[this->NAMRLevels()-1], 0,
2301 *y[this->NAMRLevels()-1], 0, ncomp, nghost, true);
2302 ParallelAllReduce::Sum(result, ParallelContext::CommunicatorSub());
2303 return result;
2304}
2305
2306template <typename MF>
2307auto
2308MLCellLinOpT<MF>::norm2Precond (Vector<MF const*> const& x) const -> RT
2309{
2310 const int ncomp = this->getNComp();
2311 const IntVect nghost(0);
2312 RT result = 0;
2313 for (int ilev = 0; ilev < this->NAMRLevels()-1; ++ilev) {
2314 result += amrex::Dot(*m_norm_fine_mask[ilev], *x[ilev], 0, ncomp, nghost, true);
2315 }
2316 result += amrex::Dot(*x[this->NAMRLevels()-1], 0, ncomp, nghost, true);
2317 ParallelAllReduce::Sum(result, ParallelContext::CommunicatorSub());
2318 return std::sqrt(result);
2319}
2320
2321template <typename MF>
2322void
2323MLCellLinOpT<MF>::computeVolInv () const
2324{
2325 if (!m_volinv.empty()) { return; }
2326
2327 m_volinv.resize(this->m_num_amr_levels);
2328 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
2329 m_volinv[amrlev].resize(this->NMGLevels(amrlev));
2330 }
2331
2332 // We don't need to compute for every level
2333
2334 auto f = [&] (int amrlev, int mglev) {
2335#ifdef AMREX_USE_EB
2336 const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(amrlev,mglev));
2337 if (factory && !factory->isAllRegular())
2338 {
2339 if constexpr (std::is_same<MF,MultiFab>()) {
2340 const auto& vfrac = factory->getVolFrac();
2341 m_volinv[amrlev][mglev] = vfrac.sum(0,true);
2342 } else {
2343 amrex::Abort("MLCellLinOp with EB only works with MultiFab");
2344 }
2345 }
2346 else
2347#endif
2348 {
2349 auto const npts = (this->m_coarse_fine_bc_type == LinOpBCType::Dirichlet)
2350 ? this->compactify(this->Geom(amrlev,mglev).Domain()).d_numPts()
2351 : this->m_grids[amrlev][mglev].d_numPts();
2352 AMREX_ASSERT(npts > 0.);
2353 m_volinv[amrlev][mglev] = RT(1.0 / npts);
2354 }
2355 };
2356
2357 // amrlev = 0, mglev = 0
2358 f(0,0);
2359
2360 int mgbottom = this->NMGLevels(0)-1;
2361 f(0,mgbottom);
2362
2363#ifdef AMREX_USE_EB
2364 RT temp1, temp2;
2365 const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(0,0));
2366 if (factory && !factory->isAllRegular())
2367 {
2368 ParallelAllReduce::Sum<RT>({m_volinv[0][0], m_volinv[0][mgbottom]},
2369 ParallelContext::CommunicatorSub());
2370 temp1 = RT(1.0)/m_volinv[0][0];
2371 temp2 = RT(1.0)/m_volinv[0][mgbottom];
2372 }
2373 else
2374 {
2375 temp1 = m_volinv[0][0];
2376 temp2 = m_volinv[0][mgbottom];
2377 }
2378 m_volinv[0][0] = temp1;
2379 m_volinv[0][mgbottom] = temp2;
2380#endif
2381}
2382
2383template <typename MF>
2384auto
2385MLCellLinOpT<MF>::normInf (int amrlev, MF const& mf, bool local) const -> RT
2386{
2387 const int ncomp = this->getNComp();
2388 const int finest_level = this->NAMRLevels() - 1;
2389 RT norm = RT(0.0);
2390#ifdef AMREX_USE_EB
2391 const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(this->Factory(amrlev));
2392 if (factory && !factory->isAllRegular()) {
2393#if defined(AMREX_USE_CUDA) && defined(_WIN32)
2394 if (!std::is_same<MF,MultiFab>()) {
2395#else
2396 if constexpr (!std::is_same<MF,MultiFab>()) {
2397#endif
2398 amrex::Abort("MLCellLinOpT with EB only works with MultiFab");
2399 } else {
2400 const MultiFab& vfrac = factory->getVolFrac();
2401 // |r*vfrac| with NaN as +inf; covered cells do not count.
2402 Math::detail::AbsNanToInf<Real> const absinf{};
2403 if (amrlev == finest_level) {
2404#ifdef AMREX_USE_GPU
2405 if (Gpu::inLaunchRegion()) {
2406 auto const& ma = mf.const_arrays();
2407 auto const& vfrac_ma = vfrac.const_arrays();
2408 norm = ParReduce(TypeList<ReduceOpMax>{}, TypeList<Real>{},
2409 mf, IntVect(0), ncomp,
2410 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n)
2411 -> GpuTuple<Real>
2412 {
2413 Real const vf = vfrac_ma[box_no](i,j,k);
2414 return (vf > Real(0)) ? absinf(ma[box_no](i,j,k,n))*vf : Real(0);
2415 });
2416 } else
2417#endif
2418 {
2419#ifdef AMREX_USE_OMP
2420#pragma omp parallel reduction(max:norm)
2421#endif
2422 for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) {
2423 Box const& bx = mfi.tilebox();
2424 auto const& fab = mf.const_array(mfi);
2425 auto const& v = vfrac.const_array(mfi);
2426 AMREX_LOOP_4D(bx, ncomp, i, j, k, n,
2427 {
2428 if (v(i,j,k) > Real(0)) {
2429 norm = std::max(norm, absinf(fab(i,j,k,n))*v(i,j,k));
2430 }
2431 });
2432 }
2433 }
2434 } else {
2435#ifdef AMREX_USE_GPU
2436 if (Gpu::inLaunchRegion()) {
2437 auto const& ma = mf.const_arrays();
2438 auto const& mask_ma = m_norm_fine_mask[amrlev]->const_arrays();
2439 auto const& vfrac_ma = vfrac.const_arrays();
2440 norm = ParReduce(TypeList<ReduceOpMax>{}, TypeList<Real>{},
2441 mf, IntVect(0), ncomp,
2442 [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n)
2443 -> GpuTuple<Real>
2444 {
2445 if (mask_ma[box_no](i,j,k)) {
2446 Real const vf = vfrac_ma[box_no](i,j,k);
2447 return (vf > Real(0)) ? absinf(ma[box_no](i,j,k,n))*vf : Real(0);
2448 } else {
2449 return Real(0.0);
2450 }
2451 });
2452 } else
2453#endif
2454 {
2455#ifdef AMREX_USE_OMP
2456#pragma omp parallel reduction(max:norm)
2457#endif
2458 for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) {
2459 Box const& bx = mfi.tilebox();
2460 auto const& fab = mf.const_array(mfi);
2461 auto const& mask = m_norm_fine_mask[amrlev]->const_array(mfi);
2462 auto const& v = vfrac.const_array(mfi);
2463 AMREX_LOOP_4D(bx, ncomp, i, j, k, n,
2464 {
2465 if (mask(i,j,k) && v(i,j,k) > Real(0)) {
2466 norm = std::max(norm, absinf(fab(i,j,k,n))*v(i,j,k));
2467 }
2468 });
2469 }
2470 }
2471 }
2472 }
2473 } else
2474#endif
2475 {
2476 if (amrlev == finest_level) {
2477 norm = mf.norminf(0, ncomp, IntVect(0), true);
2478 } else {
2479 norm = mf.norminf(*m_norm_fine_mask[amrlev], 0, ncomp, IntVect(0), true);
2480 }
2481 }
2482
2483 if (!local) { ParallelAllReduce::Max(norm, ParallelContext::CommunicatorSub()); }
2484 return norm;
2485}
2486
2487template <typename MF>
2488void
2489MLCellLinOpT<MF>::averageDownAndSync (Vector<MF>& sol) const
2490{
2491 int ncomp = this->getNComp();
2492 for (int falev = this->NAMRLevels()-1; falev > 0; --falev)
2493 {
2494#ifdef AMREX_USE_EB
2495 if (!sol[falev].isAllRegular()) {
2496 if constexpr (std::is_same<MF,MultiFab>()) {
2497 amrex::EB_average_down(sol[falev], sol[falev-1], 0, ncomp, this->AMRRefRatioVect(falev-1));
2498 } else {
2499 amrex::Abort("EB_average_down only works with MultiFab");
2500 }
2501 } else
2502#endif
2503 {
2504 amrex::average_down(sol[falev], sol[falev-1], 0, ncomp, this->AMRRefRatioVect(falev-1));
2505 }
2506 }
2507}
2508
2509template <typename MF>
2510void
2511MLCellLinOpT<MF>::avgDownResAmr (int clev, MF& cres, MF const& fres) const
2512{
2513#ifdef AMREX_USE_EB
2514 if (!fres.isAllRegular()) {
2515 if constexpr (std::is_same<MF,MultiFab>()) {
2516 amrex::EB_average_down(fres, cres, 0, this->getNComp(),
2517 this->AMRRefRatioVect(clev));
2518 } else {
2519 amrex::Abort("EB_average_down only works with MultiFab");
2520 }
2521 } else
2522#endif
2523 {
2524 amrex::average_down(fres, cres, 0, this->getNComp(),
2525 this->AMRRefRatioVect(clev));
2526 }
2527}
2528
2529template <typename MF>
2530void
2531MLCellLinOpT<MF>::beginPrecondBC ()
2532{
2533 this->m_precond_mode = true;
2534
2535 if (m_bndry_sol_zero.empty()) {
2536 m_bndry_sol_zero.resize(m_bndry_sol.size());
2537 const int ncomp = this->getNComp();
2538 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
2539 m_bndry_sol_zero[amrlev] = std::make_unique<MLMGBndryT<MF>>
2540 (this->m_grids[amrlev][0],
2541 this->m_dmap[amrlev][0],
2542 ncomp,
2543 this->m_geom[amrlev][0]);
2544 }
2545 std::swap(m_bndry_sol, m_bndry_sol_zero);
2546 MF const* coarse_data_for_bc_save = this->m_coarse_data_for_bc;
2547 this->m_coarse_data_for_bc = nullptr;
2548 for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
2549 this->setLevelBC(amrlev, nullptr);
2550 }
2551 this->m_coarse_data_for_bc = coarse_data_for_bc_save;
2552 } else {
2553 std::swap(m_bndry_sol, m_bndry_sol_zero);
2554 }
2555}
2556
2557template <typename MF>
2558void
2559MLCellLinOpT<MF>::endPrecondBC ()
2560{
2561 this->m_precond_mode = false;
2562 std::swap(m_bndry_sol, m_bndry_sol_zero);
2563}
2564
2565extern template class MLCellLinOpT<MultiFab>;
2566
2567using MLCellLinOp = MLCellLinOpT<MultiFab>;
2568
2569}
2570
2571#endif
BL_TO_FORTRAN_BOX
#define BL_TO_FORTRAN_BOX(x)
Definition AMReX_ArrayLim.H:51
BL_TO_FORTRAN_ANYD
#define BL_TO_FORTRAN_ANYD(x)
Definition AMReX_ArrayLim.H:44
BL_PROFILE
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:562
AMREX_ALWAYS_ASSERT_WITH_MESSAGE
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:49
BL_ASSERT
#define BL_ASSERT(EX)
Definition AMReX_BLassert.H:39
AMREX_ASSERT
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
AMREX_ALWAYS_ASSERT
#define AMREX_ALWAYS_ASSERT(EX)
Definition AMReX_BLassert.H:50
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:124
AMREX_HOST_DEVICE_PARALLEL_FOR_3D
#define AMREX_HOST_DEVICE_PARALLEL_FOR_3D(...)
Definition AMReX_GpuLaunchMacrosC.nolint.H:110
AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE
#define AMREX_GPU_LAUNCH_HOST_DEVICE_LAMBDA_RANGE(TN, TI, block)
Definition AMReX_GpuLaunchMacrosC.nolint.H:4
AMREX_HOST_DEVICE_PARALLEL_FOR_4D
#define AMREX_HOST_DEVICE_PARALLEL_FOR_4D(...)
Definition AMReX_GpuLaunchMacrosC.nolint.H:111
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
offset
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1139
cbx
Box cbx
Definition AMReX_HypreMLABecLap.cpp:1141
fine
Array4< Real > fine
Definition AMReX_InterpFaceRegister.cpp:90
mask
Array4< int const > mask
Definition AMReX_InterpFaceRegister.cpp:93
crse
Array4< Real const > crse
Definition AMReX_InterpFaceRegister.cpp:92
AMREX_LOOP_4D
#define AMREX_LOOP_4D(bx, ncomp, i, j, k, n, block)
Definition AMReX_Loop.nolint.H:16
s
GpuArray< MultiArray4< Real const >, 3 > s
Definition AMReX_MLEBNodeFDLaplacian.cpp:214
AMReX_MLLinOp.H
AMReX_MLLinOp_F.H
amrex_mllinop_apply_bc
void amrex_mllinop_apply_bc(const int *lo, const int *hi, amrex_real *phi, const int *philo, const int *phihi, const int *mask, const int *mlo, const int *mhi, int cdir, int bct, amrex_real bcl, const amrex_real *bcval, const int *blo, const int *bhi, int maxorder, const amrex_real *dxinv, int inhomog, int nc, int cross)
AMReX_MultiFabUtil.H
AMREX_D_TERM
#define AMREX_D_TERM(a, b, c)
Definition AMReX_SPACE.H:172
AMREX_D_DECL
#define AMREX_D_DECL(a, b, c)
Definition AMReX_SPACE.H:171
AMReX_YAFluxRegister.H
Yet-another flux register for refluxing.
AMReX_iMultiFab.H
amrex::BndryDataT::bndryValues
const FabSetT< MF > & bndryValues(Orientation face) const noexcept
Access the boundary values stored on orientation face.
Definition AMReX_BndryData.H:93
amrex::BoundCond
Maintain an identifier for boundary condition types.
Definition AMReX_BoundCond.H:25
amrex::BoxArray
Reference-counted collection of Boxes.
Definition AMReX_BoxArray.H:681
amrex::BoxArray::define
void define(const Box &bx)
Initialize the BoxArray from a single box.
Definition AMReX_BoxArray.cpp:352
amrex::BoxArray::coarsen
BoxArray & coarsen(int refinement_ratio)
Coarsen each Box in the BoxArray by refinement_ratio.
Definition AMReX_BoxArray.cpp:685
amrex::BoxArray::size
Long size() const noexcept
Return the number of boxes in the BoxArray.
Definition AMReX_BoxArray.H:758
amrex::BoxList
A list of Boxes sharing a common IndexType.
Definition AMReX_BoxList.H:109
amrex::BoxND::setBig
__host__ __device__ BoxND & setBig(const IntVectND< dim > &bg) noexcept
Redefine the big end of the BoxND.
Definition AMReX_Box.H:516
amrex::BoxND::length
__host__ __device__ IntVectND< dim > length() const noexcept
Return the length of the BoxND.
Definition AMReX_Box.H:167
amrex::BoxND::contains
__host__ __device__ bool contains(const IntVectND< dim > &p) const noexcept
Return true if argument is contained within BoxND.
Definition AMReX_Box.H:233
amrex::BoxND::size
__host__ __device__ IntVectND< dim > size() const noexcept
Return the length of the BoxND.
Definition AMReX_Box.H:160
amrex::BoxND::coarsen
__host__ __device__ BoxND & coarsen(int ref_ratio) noexcept
Coarsen BoxND by given (positive) refinement ratio. NOTE: if type(dir) = CELL centered: lo <- lo/rati...
Definition AMReX_Box.H:754
amrex::BoxND::refine
__host__ __device__ BoxND & refine(int ref_ratio) noexcept
Refine BoxND by given (positive) refinement ratio. NOTE: if type(dir) = CELL centered: lo <- lo*ratio...
Definition AMReX_Box.H:730
amrex::BoxND< 3 >::TheUnitBox
__host__ static __device__ BoxND TheUnitBox() noexcept
This static member function returns a constant reference to an object of type BoxND representing the ...
Definition AMReX_Box.H:783
amrex::BoxND::smallEnd
__host__ __device__ const IntVectND< dim > & smallEnd() const &noexcept
Return the inclusive lower bound of the box.
Definition AMReX_Box.H:124
amrex::CoordSys::CellSize
const Real * CellSize() const noexcept
Returns the cellsize for each coordinate direction.
Definition AMReX_CoordSys.H:79
amrex::DistributionMapping
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:51
amrex::EBFArrayBoxFactory
Definition AMReX_EBFabFactory.H:32
amrex::FabArray
An Array of FortranArrayBox(FAB)-like Objects.
Definition AMReX_FabArray.H:356
amrex::FabArray::const_arrays
MultiArray4< typename FabArray< FAB >::value_type const > const_arrays() const noexcept
Read-only convenience wrapper equivalent to arrays() const.
Definition AMReX_FabArray.H:731
amrex::FabArray::const_array
Array4< typename FabArray< FAB >::value_type const > const_array(const MFIter &mfi) const noexcept
Synonym for array(const MFIter&) that highlights read-only semantics.
Definition AMReX_FabArray.H:649
amrex::FabFactory
Abstract factory interface for creating, aliasing, and destroying FAB objects.
Definition AMReX_FabFactory.H:73
amrex::Geometry
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:85
amrex::Geometry::periodicity
Periodicity periodicity() const noexcept
Return the Periodicity based on the length of the domain.
Definition AMReX_Geometry.H:424
amrex::Geometry::ProbLo
const Real * ProbLo() const noexcept
Returns the lo end of the problem domain in each dimension.
Definition AMReX_Geometry.H:208
amrex::GpuTuple
GPU-compatible tuple.
Definition AMReX_Tuple.H:104
amrex::IntVectND::allGT
__host__ __device__ constexpr bool allGT(const IntVectND< dim > &rhs) const noexcept
Returns true if this is greater than argument for all components. NOTE: This is NOT a strict weak ord...
Definition AMReX_IntVect.H:517
amrex::IntVectND::setVal
__host__ __device__ IntVectND & setVal(int i, int val) noexcept
Set i'th coordinate of IntVectND to val.
Definition AMReX_IntVect.H:373
amrex::IntVectND< 3 >::TheCellVector
__host__ static __device__ constexpr IntVectND< dim > TheCellVector() noexcept
This static member function returns a reference to a constant IntVectND object, all of whose dim argu...
Definition AMReX_IntVect.H:832
amrex::IntVectND::max
__host__ __device__ constexpr int max() const noexcept
maximum (no absolute values) value
Definition AMReX_IntVect.H:313
amrex::InterpBndryDataT
An InterpBndryData object adds to a BndryData object the ability to manipulate and set the data store...
Definition AMReX_InterpBndryData.H:45
amrex::InterpBndryDataT::IBD_max_order_DEF
static constexpr int IBD_max_order_DEF
Definition AMReX_InterpBndryData.H:126
amrex::LayoutData
a one-thingy-per-box distributed object
Definition AMReX_LayoutData.H:13
amrex::MFIter
Iterator for looping ever tiles and boxes of amrex::FabArray based containers.
Definition AMReX_MFIter.H:88
amrex::MFIter::isValid
bool isValid() const noexcept
Is the iterator valid i.e. is it associated with a FAB?
Definition AMReX_MFIter.H:176
amrex::MLCellLinOpT
Definition AMReX_MLCellLinOp.H:31
amrex::MLCellLinOpT::Fsmooth
virtual void Fsmooth(int amrlev, int mglev, MF &sol, const MF &rhs, int redblack) const =0
amrex::MLCellLinOpT::getSolvabilityOffset
Vector< RT > getSolvabilityOffset(int amrlev, int mglev, MF const &rhs) const override
Compute the average offset needed to enforce solvability constraints.
Definition AMReX_MLCellLinOp.H:1904
amrex::MLCellLinOpT::averageDownSolutionRHS
void averageDownSolutionRHS(int camrlev, MF &crse_sol, MF &crse_rhs, const MF &fine_sol, const MF &fine_rhs) override
Average fine solution/RHS onto the next coarser AMR level.
Definition AMReX_MLCellLinOp.H:1535
amrex::MLCellLinOpT::averageDownAndSync
void averageDownAndSync(Vector< MF > &sol) const override
Average the solution hierarchy down (fine-to-coarse) and sync.
Definition AMReX_MLCellLinOp.H:2489
amrex::MLCellLinOpT::smooth
void smooth(int amrlev, int mglev, MF &sol, const MF &rhs, bool skip_fillboundary, int niter) const override
Perform niter smoothing iterations on the supplied residual equation.
Definition AMReX_MLCellLinOp.H:1556
amrex::MLCellLinOpT::makeNGrids
BoxArray makeNGrids(int grid_size) const
Helper that builds a BoxArray for NSolve with boxes no larger than the requested grid_size.
Definition AMReX_MLCellLinOp.H:1206
amrex::MLCellLinOpT::m_fluxreg
Vector< YAFluxRegisterT< MF > > m_fluxreg
Definition AMReX_MLCellLinOp.H:496
amrex::MLCellLinOpT::applyBC
virtual void applyBC(int amrlev, int mglev, MF &in, BCMode bc_mode, StateMode s_mode, const MLMGBndryT< MF > *bndry=nullptr, bool skip_fillboundary=false) const
Apply physical BCs (optionally skipping FillBoundary).
Definition AMReX_MLCellLinOp.H:986
amrex::MLCellLinOpT::updateSolBC
void updateSolBC(int amrlev, const MF &crse_bcdata) const
Refresh stored solution BC data from coarse inputs.
Definition AMReX_MLCellLinOp.H:955
amrex::MLCellLinOpT::m_bndry_sol
Vector< std::unique_ptr< MLMGBndryT< MF > > > m_bndry_sol
Definition AMReX_MLCellLinOp.H:436
amrex::MLCellLinOpT::operator=
MLCellLinOpT< MF > & operator=(const MLCellLinOpT< MF > &)=delete
amrex::MLCellLinOpT::compGrad
void compGrad(int amrlev, const Array< MF *, 3 > &grad, MF &sol, Location loc) const override
Compute directional gradients of the solution.
Definition AMReX_MLCellLinOp.H:1736
amrex::MLCellLinOpT::avgDownResAmr
void avgDownResAmr(int clev, MF &cres, MF const &fres) const override
Average a residual from a fine AMR level to its coarse parent.
Definition AMReX_MLCellLinOp.H:2511
amrex::MLCellLinOpT::reflux
void reflux(int crse_amrlev, MF &res, const MF &crse_sol, const MF &, MF &, MF &fine_sol, const MF &) const final
Reflux fine-level fluxes into the coarse residual.
Definition AMReX_MLCellLinOp.H:1626
amrex::MLCellLinOpT::dotProductPrecond
RT dotProductPrecond(Vector< MF const * > const &x, Vector< MF const * > const &y) const final
Dot product over the composite AMR hierarchy, excluding cells covered by finer levels (used when the ...
Definition AMReX_MLCellLinOp.H:2292
amrex::MLCellLinOpT::Fapply
virtual void Fapply(int amrlev, int mglev, MF &out, const MF &in) const =0
amrex::MLCellLinOpT::RT
typename FabDataType< MF >::value_type RT
Definition AMReX_MLCellLinOp.H:35
amrex::MLCellLinOpT::endPrecondBC
void endPrecondBC() override
Called when the operator stops being used as a preconditioner.
Definition AMReX_MLCellLinOp.H:2559
amrex::MLCellLinOpT::update
void update() override
Update for reuse.
Definition AMReX_MLCellLinOp.H:948
amrex::MLCellLinOpT::correctionResidual
void correctionResidual(int amrlev, int mglev, MF &resid, MF &x, const MF &b, BCMode bc_mode, const MF *crse_bcdata=nullptr) final
Compute the correction residual with optional coarse data.
Definition AMReX_MLCellLinOp.H:1600
amrex::MLCellLinOpT::MLCellLinOpT
MLCellLinOpT(const MLCellLinOpT< MF > &)=delete
amrex::MLCellLinOpT::BCMode
typename MLLinOpT< MF >::BCMode BCMode
Definition AMReX_MLCellLinOp.H:38
amrex::MLCellLinOpT::normInf
RT normInf(int amrlev, MF const &mf, bool local) const override
Infinity norm helper used by solvers and diagnostics.
Definition AMReX_MLCellLinOp.H:2385
amrex::MLCellLinOpT::isCrossStencil
virtual bool isCrossStencil() const
Whether the stencil is the cross shape.
Definition AMReX_MLCellLinOp.H:107
amrex::MLCellLinOpT::updateCorBC
void updateCorBC(int amrlev, const MF &crse_bcdata) const
Refresh stored correction BC data from coarse inputs.
Definition AMReX_MLCellLinOp.H:971
amrex::MLCellLinOpT::RealTuple
Array< RT, 2 *3 > RealTuple
Definition AMReX_MLCellLinOp.H:447
amrex::MLCellLinOpT::interpolation
void interpolation(int amrlev, int fmglev, MF &fine, const MF &crse) const override
Add the prolongation of coarse data onto the fine grid (fine += prolong(crse)).
Definition AMReX_MLCellLinOp.H:1268
amrex::MLCellLinOpT::xdoty
RT xdoty(int amrlev, int mglev, const MF &x, const MF &y, bool local) const final
Dot product helper.
Definition AMReX_MLCellLinOp.H:2278
amrex::MLCellLinOpT::BCTuple
Array< BoundCond, 2 *3 > BCTuple
Definition AMReX_MLCellLinOp.H:448
amrex::MLCellLinOpT::prepareForSolve
void prepareForSolve() override
Prepare multilevel metadata before MLMG iterates (coefficients, BC caches, etc.).
Definition AMReX_MLCellLinOp.H:1961
amrex::MLCellLinOpT::unapplyMetricTerm
void unapplyMetricTerm(int amrlev, int mglev, MF &rhs) const final
Remove metric scaling previously applied to the RHS.
Definition AMReX_MLCellLinOp.H:1846
amrex::MLCellLinOpT::apply
void apply(int amrlev, int mglev, MF &out, MF &in, BCMode bc_mode, StateMode s_mode, const MLMGBndryT< MF > *bndry=nullptr) const override
Apply the linear operator with boundary conditions.
Definition AMReX_MLCellLinOp.H:1546
amrex::MLCellLinOpT::applyMetricTerm
void applyMetricTerm(int amrlev, int mglev, MF &rhs) const final
Multiply the RHS by metric terms appropriate for curvilinear coordinates.
Definition AMReX_MLCellLinOp.H:1790
amrex::MLCellLinOpT::m_crse_cor_br
Vector< std::unique_ptr< BndryRegisterT< MF > > > m_crse_cor_br
Definition AMReX_MLCellLinOp.H:440
amrex::MLCellLinOpT::norm2Precond
RT norm2Precond(Vector< MF const * > const &x) const final
L2 norm over the composite AMR hierarchy, excluding cells covered by finer levels (used when the oper...
Definition AMReX_MLCellLinOp.H:2308
amrex::MLCellLinOpT::beginPrecondBC
void beginPrecondBC() override
Called when the operator starts being used as a preconditioner.
Definition AMReX_MLCellLinOp.H:2531
amrex::MLCellLinOpT::m_bcondloc
Vector< Vector< std::unique_ptr< BndryCondLoc > > > m_bcondloc
Definition AMReX_MLCellLinOp.H:486
amrex::MLCellLinOpT::restriction
void restriction(int amrlev, int cmglev, MF &crse, MF &fine) const override
Restrict a fine-grid field onto its coarse counterpart.
Definition AMReX_MLCellLinOp.H:1256
amrex::MLCellLinOpT::setGaussSeidel
void setGaussSeidel(bool flag) noexcept
Toggle Gauss–Seidel smoothing in place of Jacobi relaxation.
Definition AMReX_MLCellLinOp.H:104
amrex::MLCellLinOpT::m_robin_bcval
Vector< std::unique_ptr< MF > > m_robin_bcval
Definition AMReX_MLCellLinOp.H:423
amrex::MLCellLinOpT::define
void define(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const LPInfo &a_info=LPInfo(), const Vector< FabFactory< FAB > const * > &a_factory={})
Bind the operator to an AMR hierarchy.
Definition AMReX_MLCellLinOp.H:648
amrex::MLCellLinOpT::~MLCellLinOpT
~MLCellLinOpT() override=default
amrex::MLCellLinOpT::isTensorOp
virtual bool isTensorOp() const
Whether this operator is a tensor solve.
Definition AMReX_MLCellLinOp.H:109
amrex::MLCellLinOpT::setInterpBndryHalfWidth
void setInterpBndryHalfWidth(int w)
Control how many cells the interpolation boundary stencil spans.
Definition AMReX_MLCellLinOp.H:430
amrex::MLCellLinOpT::m_undrrelxr
Vector< Vector< BndryRegisterT< MF > > > m_undrrelxr
Definition AMReX_MLCellLinOp.H:489
amrex::MLCellLinOpT::MLCellLinOpT
MLCellLinOpT()
Definition AMReX_MLCellLinOp.H:641
amrex::MLCellLinOpT::MLCellLinOpT
MLCellLinOpT(MLCellLinOpT< MF > &&)=delete
amrex::MLCellLinOpT::FAB
typename FabDataType< MF >::fab_type FAB
Definition AMReX_MLCellLinOp.H:34
amrex::MLCellLinOpT::interpAssign
void interpAssign(int amrlev, int fmglev, MF &fine, MF &crse) const override
Overwrite fine data with the prolongation of coarse data (fine = prolong(crse)).
Definition AMReX_MLCellLinOp.H:1318
amrex::MLCellLinOpT::fixSolvabilityByOffset
void fixSolvabilityByOffset(int amrlev, int mglev, MF &rhs, Vector< RT > const &offset) const override
Apply solvability offsets to the RHS (subtracting the average).
Definition AMReX_MLCellLinOp.H:1941
amrex::MLCellLinOpT::m_bndry_cor
Vector< std::unique_ptr< MLMGBndryT< MF > > > m_bndry_cor
Definition AMReX_MLCellLinOp.H:439
amrex::MLCellLinOpT::compFlux
void compFlux(int amrlev, const Array< MF *, 3 > &fluxes, MF &sol, Location loc) const override
Compute face-centered fluxes from the supplied solution.
Definition AMReX_MLCellLinOp.H:1695
amrex::MLCellLinOpT::Location
typename MLLinOpT< MF >::Location Location
Definition AMReX_MLCellLinOp.H:40
amrex::MLCellLinOpT::prepareForFluxes
void prepareForFluxes(int amrlev, const MF *crse_bcdata=nullptr) override
Ensure BC caches are ready for flux computations (e.g., getFluxes).
Definition AMReX_MLCellLinOp.H:1591
amrex::MLCellLinOpT::solutionResidual
void solutionResidual(int amrlev, MF &resid, MF &x, const MF &b, const MF *crse_bcdata=nullptr) override
Compute the residual resid = b - A(x) using solution boundary data.
Definition AMReX_MLCellLinOp.H:1573
amrex::MLCellLinOpT::m_has_metric_term
bool m_has_metric_term
Definition AMReX_MLCellLinOp.H:434
amrex::MLCellLinOpT::m_use_gauss_seidel
bool m_use_gauss_seidel
Definition AMReX_MLCellLinOp.H:498
amrex::MLCellLinOpT::m_bndry_sol_zero
Vector< std::unique_ptr< MLMGBndryT< MF > > > m_bndry_sol_zero
Definition AMReX_MLCellLinOp.H:442
amrex::MLCellLinOpT::m_crse_sol_br
Vector< std::unique_ptr< BndryRegisterT< MF > > > m_crse_sol_br
Definition AMReX_MLCellLinOp.H:437
amrex::MLCellLinOpT::interpolationAmr
void interpolationAmr(int famrlev, MF &fine, const MF &crse, IntVect const &nghost) const override
Prolong AMR-level data during FMG initialization.
Definition AMReX_MLCellLinOp.H:1414
amrex::MLCellLinOpT::addInhomogNeumannFlux
virtual void addInhomogNeumannFlux(int, const Array< MF *, 3 > &, MF const &, bool) const
Optional hook for adding inhomogeneous Neumann contributions.
Definition AMReX_MLCellLinOp.H:381
amrex::MLCellLinOpT::m_norm_fine_mask
Vector< std::unique_ptr< iMultiFab > > m_norm_fine_mask
Definition AMReX_MLCellLinOp.H:494
amrex::MLCellLinOpT::StateMode
typename MLLinOpT< MF >::StateMode StateMode
Definition AMReX_MLCellLinOp.H:39
amrex::MLCellLinOpT::setLevelBC
void setLevelBC(int amrlev, const MF *levelbcdata, const MF *robinbc_a=nullptr, const MF *robinbc_b=nullptr, const MF *robinbc_f=nullptr) final
Provide per-level inhomogeneous boundary data.
Definition AMReX_MLCellLinOp.H:812
amrex::MLCellLinOpT::needsUpdate
bool needsUpdate() const override
Does it need update if it's reused?
Definition AMReX_MLCellLinOp.H:93
amrex::MLCellLinOpT::m_maskvals
Vector< Vector< Array< MultiMask, 2 *3 > > > m_maskvals
Definition AMReX_MLCellLinOp.H:492
amrex::MLCellLinOpT::FFlux
virtual void FFlux(int amrlev, const MFIter &mfi, const Array< FAB *, 3 > &flux, const FAB &sol, Location loc, int face_only=0) const =0
amrex::MLLinOpT
Abstract base class for multilevel linear operators used by MLMG and the bottom solvers.
Definition AMReX_MLLinOp.H:139
amrex::MLLinOpT::needsUpdate
virtual bool needsUpdate() const
Does it need update if it's reused?
Definition AMReX_MLLinOp.H:362
amrex::MLLinOpT::update
virtual void update()
Update for reuse.
Definition AMReX_MLLinOp.H:364
amrex::MLLinOpT::define
void define(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const LPInfo &a_info, const Vector< FabFactory< FAB > const * > &a_factory, bool eb_limit_coarsening=true)
Initialize the operator hierarchy on a set of AMR levels.
Definition AMReX_MLLinOp.H:1150
amrex::MLMGBndryT
Boundary helper for MLMG that manages coarse/fine and physical BC metadata.
Definition AMReX_MLMGBndry.H:20
amrex::MLMGBndryT::setBoxBC
static void setBoxBC(RealTuple &bloc, BCTuple &bctag, const Box &bx, const Box &domain, const Array< LinOpBCType, 3 > &lo, const Array< LinOpBCType, 3 > &hi, const Real *dx, IntVect const &ratio, const RealVect &interior_bloc, const Array< Real, 3 > &domain_bloc_lo, const Array< Real, 3 > &domain_bloc_hi, const GpuArray< int, 3 > &is_periodic, LinOpBCType a_crse_fine_bc_type)
Helper that sets up BC tuples for a single box.
Definition AMReX_MLMGBndry.H:163
amrex::Mask
Definition AMReX_Mask.H:33
amrex::MultiFab
A collection (stored as an array) of FArrayBox objects.
Definition AMReX_MultiFab.H:40
amrex::OrientationIter
An Iterator over the Orientation of Faces of a Box.
Definition AMReX_Orientation.H:135
amrex::Orientation
Encapsulation of the Orientation of the Faces of a Box.
Definition AMReX_Orientation.H:29
amrex::Orientation::low
@ low
Definition AMReX_Orientation.H:34
amrex::Orientation::high
@ high
Definition AMReX_Orientation.H:34
amrex::PODVector
Dynamically allocated vector for trivially copyable data.
Definition AMReX_PODVector.H:308
amrex::Vector
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:29
amrex::MultiFabLike
Checks if a type is amrex::MultiFab-like (i.e., FabArray<FAB>, where FAB is a BaseFabType).
Definition AMReX_Concepts.H:26
amrex::Real
amrex_real Real
Floating Point Type for Fields.
Definition AMReX_REAL.H:80
amrex::FabType::singlevalued
@ singlevalued
No multivalued cells; not entirely regular or covered.
amrex::FabType::covered
@ covered
Every cell in the region is covered.
amrex::FabType::regular
@ regular
Every cell in the region is regular.
amrex::adjCellHi
__host__ __device__ BoxND< dim > adjCellHi(const BoxND< dim > &b, int dir, int len=1) noexcept
Return the BoxND of length len adjacent to b on the high end along coordinate direction dir.
Definition AMReX_Box.H:1848
amrex::adjCellLo
__host__ __device__ BoxND< dim > adjCellLo(const BoxND< dim > &b, int dir, int len=1) noexcept
Return the BoxND of length len adjacent to b on the low end along coordinate direction dir.
Definition AMReX_Box.H:1817
amrex::surroundingNodes
__host__ __device__ BoxND< dim > surroundingNodes(const BoxND< dim > &b, int dir) noexcept
Return a BoxND with NODE based coordinates in direction dir that encloses BoxND b.
Definition AMReX_Box.H:1582
amrex::coarsen
__host__ __device__ BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Return a copy of b coarsened by the isotropic ratio ref_ratio.
Definition AMReX_Box.H:1469
amrex::adjCell
__host__ __device__ BoxND< dim > adjCell(const BoxND< dim > &b, Orientation face, int len=1) noexcept
Similar to adjCellLo and adjCellHi except that it operates on the given face of BoxND b.
Definition AMReX_Box.H:1880
amrex::shift
__host__ __device__ BoxND< dim > shift(const BoxND< dim > &b, int dir, int nzones) noexcept
Return a copy of b shifted by nzones cells in direction dir.
Definition AMReX_Box.H:1548
amrex::grow
__host__ __device__ BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Return a copy of b grown uniformly by i cells in every direction.
Definition AMReX_Box.H:1326
amrex::Array
std::array< T, N > Array
Definition AMReX_Array.H:31
amrex::The_Async_Arena
Arena * The_Async_Arena()
Definition AMReX_Arena.cpp:839
amrex::ParallelAllReduce::Sum
void Sum(Gpu::DeviceVector< T > &v, MPI_Comm comm)
Definition AMReX_GpuParallelReduce.H:37
amrex::ParallelAllReduce::Max
void Max(KeyValuePair< K, V > &vi, MPI_Comm comm)
Definition AMReX_ParallelReduce.H:133
amrex::Gpu::copyAsync
void copyAsync(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:228
amrex::Gpu::hostToDevice
static constexpr HostToDevice hostToDevice
Definition AMReX_GpuContainers.H:105
amrex::Gpu::streamSynchronize
void streamSynchronize() noexcept
Definition AMReX_GpuDevice.H:310
amrex::Gpu::inLaunchRegion
bool inLaunchRegion() noexcept
Definition AMReX_GpuControl.H:88
amrex::Gpu::notInLaunchRegion
bool notInLaunchRegion() noexcept
Definition AMReX_GpuControl.H:89
amrex::Gpu::inNoSyncRegion
bool inNoSyncRegion() noexcept
Definition AMReX_GpuControl.H:148
amrex::ParallelContext::CommunicatorSub
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition AMReX_ParallelContext.H:70
amrex
Definition AMReX_Amr.cpp:50
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
No-op helper that marks variables as intentionally unused.
Definition AMReX.H:273
amrex::norm
__host__ __device__ T norm(const GpuComplex< T > &a_z) noexcept
Return the norm (magnitude squared) of a complex number.
Definition AMReX_GpuComplex.H:349
amrex::ParReduce
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
amrex::average_down
void average_down(const MultiFab &S_fine, MultiFab &S_crse, const Geometry &fgeom, const Geometry &cgeom, int scomp, int ncomp, int rr)
Definition AMReX_MultiFabUtil.cpp:359
amrex::makeFineMask
iMultiFab makeFineMask(const BoxArray &cba, const DistributionMapping &cdm, const BoxArray &fba, const IntVect &ratio, int crse_value, int fine_value, MFInfo const &info)
Definition AMReX_MultiFabUtil.cpp:653
amrex::EB_set_covered
void EB_set_covered(MultiFab &mf, Real val)
Fill all covered cells with a single value val.
Definition AMReX_EBMultiFabUtil.cpp:21
amrex::Dot
FAB::value_type Dot(FabArray< FAB > const &x, int xcomp, FabArray< FAB > const &y, int ycomp, int ncomp, IntVect const &nghost, bool local=false)
Compute dot products of two FabArrays.
Definition AMReX_FabArrayUtility.H:1898
amrex::ParallelFor
void ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
Definition AMReX_CTOParallelForImpl.H:202
amrex::Box
BoxND< 3 > Box
Box is an alias for amrex::BoxND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:35
amrex::isMFIterSafe
bool isMFIterSafe(const FabArrayBase &x, const FabArrayBase &y)
Definition AMReX_MFIter.H:256
amrex::LinOpBCType
LinOpBCType
Definition AMReX_LO_BCTYPES.H:27
amrex::EB_average_down
void EB_average_down(const MultiFab &S_fine, MultiFab &S_crse, const MultiFab &vol_fine, const MultiFab &vfrac_fine, int scomp, int ncomp, const IntVect &ratio)
Volume-weighted average-down from fine to coarse using EB volume fractions.
Definition AMReX_EBMultiFabUtil.cpp:336
amrex::IntVect
IntVectND< 3 > IntVect
IntVect is an alias for amrex::IntVectND instantiated with AMREX_SPACEDIM.
Definition AMReX_BaseFwd.H:38
amrex::RealVect
RealVectND< 3 > RealVect
Definition AMReX_ParmParse.H:39
amrex::TilingIfNotGPU
bool TilingIfNotGPU() noexcept
Definition AMReX_MFIter.H:12
amrex::Abort
void Abort(const std::string &msg)
Print a fatal-error message to stderr and abort execution.
Definition AMReX.cpp:244
amrex::get
__host__ __device__ constexpr int get(IntVectND< dim > const &iv) noexcept
Get I'th element of IntVectND<dim>
Definition AMReX_IntVect.H:1360
amrex::ArrayND
A multidimensional array accessor.
Definition AMReX_Array4.H:289
amrex::Dim3
A simple struct holding 3 int values for a 3D index.
Definition AMReX_Dim3.H:24
amrex::Dim3::x
int x
Definition AMReX_Dim3.H:24
amrex::Dim3::z
int z
Definition AMReX_Dim3.H:24
amrex::Dim3::y
int y
Definition AMReX_Dim3.H:24
amrex::FabDataType
Type trait specialized for MultiFab-like types and their containers.
Definition AMReX_FabDataType.H:16
amrex::GpuArray
Fixed-size array that can be used on GPU.
Definition AMReX_Array.H:52
amrex::LPInfo
Configuration knobs for multilevel linear operators (grid agglomeration, metrics, etc....
Definition AMReX_MLLinOp.H:53
amrex::LinOpEnumType::StateMode
StateMode
Definition AMReX_MLLinOp.H:120
amrex::LinOpEnumType::BCMode
BCMode
Definition AMReX_MLLinOp.H:119
amrex::LinOpEnumType::Location
Location
Definition AMReX_MLLinOp.H:121
amrex::MFInfo
FabArray memory allocation information.
Definition AMReX_FabArray.H:73
amrex::MFItInfo
Definition AMReX_MFIter.H:20
amrex::MFItInfo::SetDynamic
MFItInfo & SetDynamic(bool f) noexcept
Definition AMReX_MFIter.H:43
amrex::MFItInfo::EnableTiling
MFItInfo & EnableTiling(const IntVect &ts=FabArrayBase::mfiter_tile_size) noexcept
Definition AMReX_MFIter.H:31
amrex::MLMGABCTag
Definition AMReX_MLCellLinOp.H:26
amrex::MultiArray4
Definition AMReX_FabArray.H:168
amrex::TypeList
Struct for holding types.
Definition AMReX_TypeList.H:13