Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
AMReX_MLLinOp.H
Go to the documentation of this file.
1#ifndef AMREX_ML_LINOP_H_
2#define AMREX_ML_LINOP_H_
3#include <AMReX_Config.H>
4
5#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
6#include <AMReX_Hypre.H>
8#endif
9
10#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
11#include <AMReX_PETSc.H>
12#endif
13
14#ifdef AMREX_USE_EB
16#include <AMReX_MultiCutFab.H>
17#endif
18
19#include <AMReX_Any.H>
20#include <AMReX_BndryRegister.H>
21#include <AMReX_FabDataType.H>
22#include <AMReX_MLMGBndry.H>
23#include <AMReX_MultiFabUtil.H>
24
25#include <algorithm>
26#include <string>
27
28namespace amrex {
29
30enum class BottomSolver : int {
32};
33
34struct LPInfo
35{
36 bool do_agglomeration = true;
37 bool do_consolidation = true;
38 bool do_semicoarsening = false;
39 int agg_grid_size = -1;
40 int con_grid_size = -1;
41 int con_ratio = 2;
42 int con_strategy = 3;
43 bool has_metric_term = true;
48
49 LPInfo& setAgglomeration (bool x) noexcept { do_agglomeration = x; return *this; }
50 LPInfo& setConsolidation (bool x) noexcept { do_consolidation = x; return *this; }
51 LPInfo& setSemicoarsening (bool x) noexcept { do_semicoarsening = x; return *this; }
52 LPInfo& setAgglomerationGridSize (int x) noexcept { agg_grid_size = x; return *this; }
53 LPInfo& setConsolidationGridSize (int x) noexcept { con_grid_size = x; return *this; }
54 LPInfo& setConsolidationRatio (int x) noexcept { con_ratio = x; return *this; }
55 LPInfo& setConsolidationStrategy (int x) noexcept { con_strategy = x; return *this; }
56 LPInfo& setMetricTerm (bool x) noexcept { has_metric_term = x; return *this; }
57 LPInfo& setMaxCoarseningLevel (int n) noexcept { max_coarsening_level = n; return *this; }
58 LPInfo& setMaxSemicoarseningLevel (int n) noexcept { max_semicoarsening_level = n; return *this; }
59 LPInfo& setSemicoarseningDirection (int n) noexcept { semicoarsening_direction = n; return *this; }
60 LPInfo& setHiddenDirection (int n) noexcept { hidden_direction = n; return *this; }
61
62 [[nodiscard]] bool hasHiddenDimension () const noexcept {
63 return hidden_direction >=0 && hidden_direction < AMREX_SPACEDIM;
64 }
65
66 static constexpr int getDefaultAgglomerationGridSize () {
67#ifdef AMREX_USE_GPU
68 return 32;
69#else
70 return AMREX_D_PICK(32, 16, 8);
71#endif
72 }
73
74 static constexpr int getDefaultConsolidationGridSize () {
75#ifdef AMREX_USE_GPU
76 return 32;
77#else
78 return AMREX_D_PICK(32, 16, 8);
79#endif
80 }
81};
82
89
90template <typename T> class MLMGT;
91template <typename T> class MLCGSolverT;
92template <typename T> class MLPoissonT;
93template <typename T> class MLABecLaplacianT;
94template <typename T> class GMRESMLMGT;
95
96template <typename MF>
98{
99public:
100
101 template <typename T> friend class MLMGT;
102 template <typename T> friend class MLCGSolverT;
103 template <typename T> friend class MLPoissonT;
104 template <typename T> friend class MLABecLaplacianT;
105 template <typename T> friend class GMRESMLMGT;
106
107 using MFType = MF;
110
111 using BCType = LinOpBCType;
115
116 MLLinOpT () = default;
117 virtual ~MLLinOpT () = default;
118
119 MLLinOpT (const MLLinOpT<MF>&) = delete;
120 MLLinOpT (MLLinOpT<MF>&&) = delete;
123
124 void define (const Vector<Geometry>& a_geom,
125 const Vector<BoxArray>& a_grids,
126 const Vector<DistributionMapping>& a_dmap,
127 const LPInfo& a_info,
128 const Vector<FabFactory<FAB> const*>& a_factory,
129 bool eb_limit_coarsening = true);
130
131 [[nodiscard]] virtual std::string name () const { return std::string("Unspecified"); }
132
144 const Array<BCType,AMREX_SPACEDIM>& hibc) noexcept;
145
156 const Vector<Array<BCType,AMREX_SPACEDIM> >& hibc) noexcept;
157
168 const Array<Real,AMREX_SPACEDIM>& hi_bcloc) noexcept;
169
177 [[nodiscard]] bool needsCoarseDataForBC () const noexcept { return m_needs_coarse_data_for_bc; }
178
200 void setCoarseFineBC (const MF* crse, int crse_ratio,
201 LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;
202
203 void setCoarseFineBC (const MF* crse, IntVect const& crse_ratio,
204 LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;
205
206 template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int> = 0>
207 void setCoarseFineBC (const AMF* crse, int crse_ratio,
208 LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;
209
210 template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int> = 0>
211 void setCoarseFineBC (const AMF* crse, IntVect const& crse_ratio,
212 LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;
213
232 virtual void setLevelBC (int /*amrlev*/, const MF* /*levelbcdata*/,
233 const MF* /*robinbc_a*/ = nullptr,
234 const MF* /*robinbc_b*/ = nullptr,
235 const MF* /*robinbc_f*/ = nullptr) = 0;
236
237 template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int> = 0>
238 void setLevelBC (int amrlev, const AMF* levelbcdata,
239 const AMF* robinbc_a = nullptr,
240 const AMF* robinbc_b = nullptr,
241 const AMF* robinbc_f = nullptr);
242
244 void setVerbose (int v) noexcept { verbose = v; }
245
247 void setMaxOrder (int o) noexcept { maxorder = o; }
249 [[nodiscard]] int getMaxOrder () const noexcept { return maxorder; }
250
256 [[nodiscard]] bool getEnforceSingularSolvable () const noexcept { return enforceSingularSolvable; }
257
258 [[nodiscard]] virtual BottomSolver getDefaultBottomSolver () const { return BottomSolver::bicgstab; }
259
261 [[nodiscard]] virtual int getNComp () const { return 1; }
262
263 [[nodiscard]] virtual int getNGrow (int /*a_lev*/ = 0, int /*mg_lev*/ = 0) const { return 0; }
264
266 [[nodiscard]] virtual bool needsUpdate () const { return false; }
268 virtual void update () {}
269
279 virtual void restriction (int amrlev, int cmglev, MF& crse, MF& fine) const = 0;
280
289 virtual void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const = 0;
290
299 virtual void interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const
300 {
301 amrex::ignore_unused(amrlev, fmglev, fine, crse);
302 amrex::Abort("MLLinOpT::interpAssign: Must be implemented for FMG cycle");
303 }
304
313 virtual void interpolationAmr (int famrlev, MF& fine, const MF& crse,
314 IntVect const& nghost) const
315 {
316 amrex::ignore_unused(famrlev, fine, crse, nghost);
317 amrex::Abort("MLLinOpT::interpolationAmr: Must be implemented for composite solves across multiple AMR levels");
318 }
319
329 virtual void averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs,
330 const MF& fine_sol, const MF& fine_rhs)
331 {
332 amrex::ignore_unused(camrlev, crse_sol, crse_rhs, fine_sol, fine_rhs);
333 amrex::Abort("MLLinOpT::averageDownSolutionRHS: Must be implemented for composite solves across multiple AMR levels");
334 }
335
347 virtual void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode,
348 StateMode s_mode, const MLMGBndryT<MF>* bndry=nullptr) const = 0;
349
359 virtual void smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
360 bool skip_fillboundary=false) const = 0;
361
363 virtual void normalize (int amrlev, int mglev, MF& mf) const {
364 amrex::ignore_unused(amrlev, mglev, mf);
365 }
366
376 virtual void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b,
377 const MF* crse_bcdata=nullptr) = 0;
378
379 virtual void prepareForFluxes (int /*amrlev*/, const MF* /*crse_bcdata*/ = nullptr) {}
380
392 virtual void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b,
393 BCMode bc_mode, const MF* crse_bcdata=nullptr) = 0;
394
406 virtual void reflux (int crse_amrlev,
407 MF& res, const MF& crse_sol, const MF& crse_rhs,
408 MF& fine_res, MF& fine_sol, const MF& fine_rhs) const
409 {
410 amrex::ignore_unused(crse_amrlev, res, crse_sol, crse_rhs, fine_res,
411 fine_sol, fine_rhs);
412 amrex::Abort("MLLinOpT::reflux: Must be implemented for composite solves across multiple AMR levels");
413 }
414
423 virtual void compFlux (int /*amrlev*/, const Array<MF*,AMREX_SPACEDIM>& /*fluxes*/,
424 MF& /*sol*/, Location /*loc*/) const
425 {
426 amrex::Abort("AMReX_MLLinOp::compFlux::How did we get here?");
427 }
428
437 virtual void compGrad (int /*amrlev*/, const Array<MF*,AMREX_SPACEDIM>& /*grad*/,
438 MF& /*sol*/, Location /*loc*/) const
439 {
440 amrex::Abort("AMReX_MLLinOp::compGrad::How did we get here?");
441 }
442
444 virtual void applyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const {}
446 virtual void unapplyMetricTerm (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/) const {}
447
449 virtual void unimposeNeumannBC (int /*amrlev*/, MF& /*rhs*/) const {}
450
452 virtual void applyInhomogNeumannTerm (int /*amrlev*/, MF& /*rhs*/) const {}
453
455 virtual void applyOverset (int /*amlev*/, MF& /*rhs*/) const {}
456
458 [[nodiscard]] virtual bool scaleRHS (int /*amrlev*/, MF* /*rhs*/) const {
459 return false;
460 }
461
463 virtual Vector<RT> getSolvabilityOffset (int /*amrlev*/, int /*mglev*/,
464 MF const& /*rhs*/) const { return {}; }
465
467 virtual void fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, MF& /*rhs*/,
468 Vector<RT> const& /*offset*/) const {}
469
470 virtual void prepareForSolve () = 0;
471
472 virtual void preparePrecond () {}
473
476 virtual void setDirichletNodesToZero (int /*amrlev*/, int /*mglev*/,
477 MF& /*mf*/) const
478 {
479 amrex::Warning("This function might need to be implemented for GMRES to work with this LinOp.");
480 }
481
483 [[nodiscard]] virtual bool isSingular (int amrlev) const = 0;
485 [[nodiscard]] virtual bool isBottomSingular () const = 0;
486
488 virtual RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const = 0;
489
490 virtual RT dotProductPrecond (Vector<MF const*> const& x, Vector<MF const*> const& y) const;
491
492 virtual RT norm2Precond (Vector<MF const*> const& x) const;
493
494 virtual std::unique_ptr<MLLinOpT<MF>> makeNLinOp (int /*grid_size*/) const
495 {
496 amrex::Abort("MLLinOp::makeNLinOp: N-Solve not supported");
497 return nullptr;
498 }
499
500 virtual void getFluxes (const Vector<Array<MF*,AMREX_SPACEDIM> >& /*a_flux*/,
501 const Vector<MF*>& /*a_sol*/,
502 Location /*a_loc*/) const {
503 amrex::Abort("MLLinOp::getFluxes: How did we get here?");
504 }
505 virtual void getFluxes (const Vector<MF*>& /*a_flux*/,
506 const Vector<MF*>& /*a_sol*/) const {
507 amrex::Abort("MLLinOp::getFluxes: How did we get here?");
508 }
509
510#ifdef AMREX_USE_EB
511 virtual void getEBFluxes (const Vector<MF*>& /*a_flux*/,
512 const Vector<MF*>& /*a_sol*/) const {
513 amrex::Abort("MLLinOp::getEBFluxes: How did we get here?");
514 }
515#endif
516
517#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
518 [[nodiscard]] virtual std::unique_ptr<Hypre> makeHypre (Hypre::Interface /*hypre_interface*/) const {
519 amrex::Abort("MLLinOp::makeHypre: How did we get here?");
520 return {nullptr};
521 }
522 [[nodiscard]] virtual std::unique_ptr<HypreNodeLap> makeHypreNodeLap(
523 int /*bottom_verbose*/,
524 const std::string& /* options_namespace */) const
525 {
526 amrex::Abort("MLLinOp::makeHypreNodeLap: How did we get here?");
527 return {nullptr};
528 }
529#endif
530
531#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
532 [[nodiscard]] virtual std::unique_ptr<PETScABecLap> makePETSc () const {
533 amrex::Abort("MLLinOp::makePETSc: How did we get here?");
534 return {nullptr};
535 }
536#endif
537
538 [[nodiscard]] virtual bool supportNSolve () const { return false; }
539
540 virtual void copyNSolveSolution (MF&, MF const&) const {}
541
542 virtual void postSolve (Vector<MF*> const& /*sol*/) const {}
543
544 [[nodiscard]] virtual RT normInf (int amrlev, MF const& mf, bool local) const = 0;
545
546 virtual void averageDownAndSync (Vector<MF>& sol) const = 0;
547
548 virtual void avgDownResAmr (int clev, MF& cres, MF const& fres) const
549 {
550 amrex::ignore_unused(clev, cres, fres);
551 amrex::Abort("MLLinOpT::avgDownResAmr: Must be implemented for composite solves across multiple AMR levels");
552 }
553
554 // This function is needed for FMG cycle, but not V-cycle.
555 virtual void avgDownResMG (int clev, MF& cres, MF const& fres) const;
556
557 virtual void beginPrecondBC () { m_precond_mode = true; }
558 virtual void endPrecondBC () { m_precond_mode = false; }
559
560 [[nodiscard]] bool isMFIterSafe (int amrlev, int mglev1, int mglev2) const;
561
563 [[nodiscard]] int NAMRLevels () const noexcept { return m_num_amr_levels; }
564
566 [[nodiscard]] int NMGLevels (int amrlev) const noexcept { return m_num_mg_levels[amrlev]; }
567
568 [[nodiscard]] const Geometry& Geom (int amr_lev, int mglev=0) const noexcept { return m_geom[amr_lev][mglev]; }
569
570 // BC
573 // Need to save the original copy because we change the BC type to
574 // Neumann for inhomogeneous Neumann and Robin.
577
578protected:
579
580 static constexpr int mg_coarsen_ratio = 2;
581 static constexpr int mg_box_min_width = 2;
583
585
586 int verbose = 0;
587
588 int maxorder = 3;
589
591
594
596 const MLLinOpT<MF>* m_parent = nullptr;
597
599
600 bool m_do_agglomeration = false;
601 bool m_do_consolidation = false;
602
605
612
617 CommContainer (MPI_Comm m) noexcept : comm(m) {}
618 CommContainer (const CommContainer&) = delete;
620 void operator= (const CommContainer&) = delete;
621 void operator= (CommContainer&&) = delete;
622 ~CommContainer () { // NOLINT(modernize-use-equals-default)
623#ifdef BL_USE_MPI
624 if (comm != MPI_COMM_NULL) { MPI_Comm_free(&comm); }
625#endif
626 }
627 };
628 std::unique_ptr<CommContainer> m_raii_comm;
629
632
634 LinOpBCType m_coarse_fine_bc_type = LinOpBCType::Dirichlet;
637 const MF* m_coarse_data_for_bc = nullptr;
639
640 bool m_precond_mode = false;
641
643 [[nodiscard]] const Vector<int>& AMRRefRatio () const noexcept { return m_amr_ref_ratio; }
644
646 [[nodiscard]] int AMRRefRatio (int amr_lev) const noexcept { return m_amr_ref_ratio[amr_lev]; }
647
648 [[nodiscard]] FabFactory<FAB> const* Factory (int amr_lev, int mglev=0) const noexcept {
649 return m_factory[amr_lev][mglev].get();
650 }
651
652 [[nodiscard]] GpuArray<BCType,AMREX_SPACEDIM> LoBC (int icomp = 0) const noexcept {
654 m_lobc[icomp][1],
655 m_lobc[icomp][2])}};
656 }
657 [[nodiscard]] GpuArray<BCType,AMREX_SPACEDIM> HiBC (int icomp = 0) const noexcept {
659 m_hibc[icomp][1],
660 m_hibc[icomp][2])}};
661 }
662
663 [[nodiscard]] bool hasBC (BCType bct) const noexcept;
664 [[nodiscard]] bool hasInhomogNeumannBC () const noexcept;
665 [[nodiscard]] bool hasRobinBC () const noexcept;
666
667 [[nodiscard]] virtual bool supportRobinBC () const noexcept { return false; }
668 [[nodiscard]] virtual bool supportInhomogNeumannBC () const noexcept { return false; }
669
670#ifdef BL_USE_MPI
671 [[nodiscard]] bool isBottomActive () const noexcept { return m_bottom_comm != MPI_COMM_NULL; }
672#else
673 [[nodiscard]] bool isBottomActive () const noexcept { return true; }
674#endif
675 [[nodiscard]] MPI_Comm BottomCommunicator () const noexcept { return m_bottom_comm; }
676 [[nodiscard]] MPI_Comm Communicator () const noexcept { return m_default_comm; }
677
678 void setCoarseFineBCLocation (const RealVect& cloc) noexcept { m_coarse_bc_loc = cloc; }
679
680 [[nodiscard]] bool doAgglomeration () const noexcept { return m_do_agglomeration; }
681 [[nodiscard]] bool doConsolidation () const noexcept { return m_do_consolidation; }
682 [[nodiscard]] bool doSemicoarsening () const noexcept { return m_do_semicoarsening; }
683
684 [[nodiscard]] bool isCellCentered () const noexcept { return m_ixtype == 0; }
685
686 [[nodiscard]] virtual IntVect getNGrowVectRestriction () const {
687 return isCellCentered() ? IntVect(0) : IntVect(1);
688 }
689
690 virtual void make (Vector<Vector<MF> >& mf, IntVect const& ng) const;
691
692 [[nodiscard]] virtual MF make (int amrlev, int mglev, IntVect const& ng) const;
693
694 [[nodiscard]] virtual MF makeAlias (MF const& mf) const;
695
696 [[nodiscard]] virtual MF makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const;
697
698 [[nodiscard]] virtual MF makeCoarseAmr (int famrlev, IntVect const& ng) const;
699
700 [[nodiscard]] virtual std::unique_ptr<FabFactory<FAB> > makeFactory (int /*amrlev*/, int /*mglev*/) const {
701 return std::make_unique<DefaultFabFactory<FAB>>();
702 }
703
704 virtual void resizeMultiGrid (int new_size);
705
706 [[nodiscard]] bool hasHiddenDimension () const noexcept { return info.hasHiddenDimension(); }
707 [[nodiscard]] int hiddenDirection () const noexcept { return info.hidden_direction; }
708 [[nodiscard]] Box compactify (Box const& b) const noexcept;
709
710 template <typename T>
711 [[nodiscard]] Array4<T> compactify (Array4<T> const& a) const noexcept
712 {
713 if (info.hidden_direction == 0) {
714 return Array4<T>(a.dataPtr(), {a.begin.y,a.begin.z,0}, {a.end.y,a.end.z,1}, a.nComp());
715 } else if (info.hidden_direction == 1) {
716 return Array4<T>(a.dataPtr(), {a.begin.x,a.begin.z,0}, {a.end.x,a.end.z,1}, a.nComp());
717 } else if (info.hidden_direction == 2) {
718 return Array4<T>(a.dataPtr(), {a.begin.x,a.begin.y,0}, {a.end.x,a.end.y,1}, a.nComp());
719 } else {
720 return a;
721 }
722 }
723
724 template <typename T>
725 [[nodiscard]] T get_d0 (T const& dx, T const& dy, T const&) const noexcept
726 {
727 if (info.hidden_direction == 0) {
728 return dy;
729 } else {
730 return dx;
731 }
732 }
733
734 template <typename T>
735 [[nodiscard]] T get_d1 (T const&, T const& dy, T const& dz) const noexcept
736 {
737 if (info.hidden_direction == 0 || info.hidden_direction == 1) {
738 return dz;
739 } else {
740 return dy;
741 }
742 }
743
744private:
745
746 void defineGrids (const Vector<Geometry>& a_geom,
747 const Vector<BoxArray>& a_grids,
748 const Vector<DistributionMapping>& a_dmap,
749 const Vector<FabFactory<FAB> const*>& a_factory);
750 void defineBC ();
753 int ratio, int strategy);
755
756 virtual void checkPoint (std::string const& /*file_name*/) const {
757 amrex::Abort("MLLinOp:checkPoint: not implemented");
758 }
759
764};
765
766template <typename MF>
767void
769 const Vector<BoxArray>& a_grids,
770 const Vector<DistributionMapping>& a_dmap,
771 const LPInfo& a_info,
772 const Vector<FabFactory<FAB> const*>& a_factory,
773 [[maybe_unused]] bool eb_limit_coarsening)
774{
775 BL_PROFILE("MLLinOp::define()");
776
777 info = a_info;
778#ifdef AMREX_USE_GPU
780 {
781 if (info.agg_grid_size <= 0) { info.agg_grid_size = AMREX_D_PICK(32, 16, 8); }
782 if (info.con_grid_size <= 0) { info.con_grid_size = AMREX_D_PICK(32, 16, 8); }
783 }
784 else
785#endif
786 {
787 if (info.agg_grid_size <= 0) { info.agg_grid_size = LPInfo::getDefaultAgglomerationGridSize(); }
788 if (info.con_grid_size <= 0) { info.con_grid_size = LPInfo::getDefaultConsolidationGridSize(); }
789 }
790
791#ifdef AMREX_USE_EB
792 if (!a_factory.empty() && eb_limit_coarsening) {
793 const auto *f = dynamic_cast<EBFArrayBoxFactory const*>(a_factory[0]);
794 if (f) {
795 info.max_coarsening_level = std::min(info.max_coarsening_level,
796 f->maxCoarseningLevel());
797 }
798 }
799#endif
800 defineGrids(a_geom, a_grids, a_dmap, a_factory);
801 defineBC();
802}
803
804template <typename MF>
805void
807 const Vector<BoxArray>& a_grids,
808 const Vector<DistributionMapping>& a_dmap,
809 const Vector<FabFactory<FAB> const*>& a_factory)
810{
811 BL_PROFILE("MLLinOp::defineGrids()");
812
813#ifdef AMREX_USE_EB
814 if ( ! a_factory.empty() ) {
815 auto const* ebf = dynamic_cast<EBFArrayBoxFactory const*>(a_factory[0]);
816 if (ebf && !(ebf->isAllRegular())) { // Has non-trivial EB
817 mg_domain_min_width = 4;
818 }
819 }
820#endif
821
822 m_num_amr_levels = 0;
823 for (int amrlev = 0; amrlev < a_geom.size(); amrlev++) {
824 if (!a_grids[amrlev].empty()) {
825 m_num_amr_levels++;
826 }
827 }
828
829 m_amr_ref_ratio.resize(m_num_amr_levels);
830 m_num_mg_levels.resize(m_num_amr_levels);
831
832 m_geom.resize(m_num_amr_levels);
833 m_grids.resize(m_num_amr_levels);
834 m_dmap.resize(m_num_amr_levels);
835 m_factory.resize(m_num_amr_levels);
836
837 m_default_comm = ParallelContext::CommunicatorSub();
838
839 const RealBox& rb = a_geom[0].ProbDomain();
840 const int coord = a_geom[0].Coord();
841 const Array<int,AMREX_SPACEDIM>& is_per = a_geom[0].isPeriodic();
842
843 IntVect mg_coarsen_ratio_v(mg_coarsen_ratio);
844 IntVect mg_box_min_width_v(mg_box_min_width);
845 IntVect mg_domain_min_width_v(mg_domain_min_width);
846 if (hasHiddenDimension()) {
847 AMREX_ASSERT_WITH_MESSAGE(AMREX_SPACEDIM == 3 && m_num_amr_levels == 1,
848 "Hidden direction only supported for 3d level solve");
849 mg_coarsen_ratio_v[info.hidden_direction] = 1;
850 mg_box_min_width_v[info.hidden_direction] = 0;
851 mg_domain_min_width_v[info.hidden_direction] = 0;
852 }
853
854 // fine amr levels
855 for (int amrlev = m_num_amr_levels-1; amrlev > 0; --amrlev)
856 {
857 m_num_mg_levels[amrlev] = 1;
858 m_geom[amrlev].push_back(a_geom[amrlev]);
859 m_grids[amrlev].push_back(a_grids[amrlev]);
860 m_dmap[amrlev].push_back(a_dmap[amrlev]);
861 if (amrlev < a_factory.size()) {
862 m_factory[amrlev].emplace_back(a_factory[amrlev]->clone());
863 } else {
864 m_factory[amrlev].push_back(std::make_unique<DefaultFabFactory<FAB>>());
865 }
866
867 IntVect rr = mg_coarsen_ratio_v;
868 const Box& dom = a_geom[amrlev].Domain();
869 for (int i = 0; i < 2; ++i)
870 {
871 if (!dom.coarsenable(rr)) { amrex::Abort("MLLinOp: Uncoarsenable domain"); }
872
873 const Box& cdom = amrex::coarsen(dom,rr);
874 if (cdom == a_geom[amrlev-1].Domain()) { break; }
875
876 ++(m_num_mg_levels[amrlev]);
877
878 m_geom[amrlev].emplace_back(cdom, rb, coord, is_per);
879
880 m_grids[amrlev].push_back(a_grids[amrlev]);
881 AMREX_ASSERT(m_grids[amrlev].back().coarsenable(rr));
882 m_grids[amrlev].back().coarsen(rr);
883
884 m_dmap[amrlev].push_back(a_dmap[amrlev]);
885
886 rr *= mg_coarsen_ratio_v;
887 }
888
889#if (AMREX_SPACEDIM > 1)
890 if (hasHiddenDimension()) {
891 m_amr_ref_ratio[amrlev-1] = rr[AMREX_SPACEDIM-info.hidden_direction];
892 } else
893#endif
894 {
895 m_amr_ref_ratio[amrlev-1] = rr[0];
896 }
897 }
898
899 // coarsest amr level
900 m_num_mg_levels[0] = 1;
901 m_geom[0].push_back(a_geom[0]);
902 m_grids[0].push_back(a_grids[0]);
903 m_dmap[0].push_back(a_dmap[0]);
904 if (!a_factory.empty()) {
905 m_factory[0].emplace_back(a_factory[0]->clone());
906 } else {
907 m_factory[0].push_back(std::make_unique<DefaultFabFactory<FAB>>());
908 }
909
910 m_domain_covered.resize(m_num_amr_levels, false);
911 auto npts0 = m_grids[0][0].numPts();
912 m_domain_covered[0] = (npts0 == compactify(m_geom[0][0].Domain()).numPts());
913 for (int amrlev = 1; amrlev < m_num_amr_levels; ++amrlev)
914 {
915 if (!m_domain_covered[amrlev-1]) { break; }
916 m_domain_covered[amrlev] = (m_grids[amrlev][0].numPts() ==
917 compactify(m_geom[amrlev][0].Domain()).numPts());
918 }
919
920 Box aggbox;
921 bool aggable = false;
922
923 if (m_grids[0][0].size() > 1 && info.do_agglomeration)
924 {
925 if (m_domain_covered[0])
926 {
927 aggbox = m_geom[0][0].Domain();
928 if (hasHiddenDimension()) {
929 aggbox.makeSlab(hiddenDirection(), m_grids[0][0][0].smallEnd(hiddenDirection()));
930 }
931 aggable = true;
932 }
933 else
934 {
935 aggbox = m_grids[0][0].minimalBox();
936 aggable = (aggbox.numPts() == npts0);
937 }
938 }
939
940 bool agged = false;
941 bool coned = false;
942 int agg_lev = 0, con_lev = 0;
943
944 AMREX_ALWAYS_ASSERT( ! (info.do_semicoarsening && info.hasHiddenDimension())
945 && info.semicoarsening_direction >= -1
946 && info.semicoarsening_direction < AMREX_SPACEDIM );
947
948 if (info.do_agglomeration && aggable)
949 {
950 Box dbx = m_geom[0][0].Domain();
951 Box bbx = aggbox;
952 Real const nbxs = static_cast<Real>(m_grids[0][0].size());
953 Real const threshold_npts = static_cast<Real>(AMREX_D_TERM(info.agg_grid_size,
954 *info.agg_grid_size,
955 *info.agg_grid_size));
956 Vector<Box> domainboxes{dbx};
957 Vector<Box> boundboxes{bbx};
958 Vector<int> agg_flag{false};
959 Vector<IntVect> accum_coarsen_ratio{IntVect(1)};
960 int numsclevs = 0;
961
962 for (int lev = 0; lev < info.max_coarsening_level; ++lev)
963 {
964 IntVect rr_level = mg_coarsen_ratio_v;
965 bool const do_semicoarsening_level = info.do_semicoarsening
966 && numsclevs < info.max_semicoarsening_level;
967 if (do_semicoarsening_level
968 && info.semicoarsening_direction != -1)
969 {
970 rr_level[info.semicoarsening_direction] = 1;
971 }
972 IntVect is_coarsenable;
973 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
974 IntVect rr_dir(1);
975 rr_dir[idim] = rr_level[idim];
976 is_coarsenable[idim] = dbx.coarsenable(rr_dir, mg_domain_min_width_v)
977 && bbx.coarsenable(rr_dir, mg_box_min_width_v);
978 if (!is_coarsenable[idim] && do_semicoarsening_level
979 && info.semicoarsening_direction == -1)
980 {
981 is_coarsenable[idim] = true;
982 rr_level[idim] = 1;
983 }
984 }
985 if (is_coarsenable != IntVect(1) || rr_level == IntVect(1)) {
986 break;
987 }
988 if (do_semicoarsening_level && info.semicoarsening_direction == -1) {
989 // make sure there is at most one direction that is not coarsened
990 int n_ones = AMREX_D_TERM( static_cast<int>(rr_level[0] == 1),
991 + static_cast<int>(rr_level[1] == 1),
992 + static_cast<int>(rr_level[2] == 1));
993 if (n_ones > 1) { break; }
994 }
995 if (rr_level != mg_coarsen_ratio_v) {
996 ++numsclevs;
997 }
998
999 accum_coarsen_ratio.push_back(accum_coarsen_ratio.back()*rr_level);
1000 domainboxes.push_back(dbx.coarsen(rr_level));
1001 boundboxes.push_back(bbx.coarsen(rr_level));
1002 bool to_agg = (bbx.d_numPts() / nbxs) < 0.999*threshold_npts;
1003 agg_flag.push_back(to_agg);
1004 }
1005
1006 for (int lev = 1, nlevs = static_cast<int>(domainboxes.size()); lev < nlevs; ++lev) {
1007 if (!agged && !agg_flag[lev] &&
1008 a_grids[0].coarsenable(accum_coarsen_ratio[lev], mg_box_min_width_v))
1009 {
1010 m_grids[0].push_back(amrex::coarsen(a_grids[0], accum_coarsen_ratio[lev]));
1011 m_dmap[0].push_back(a_dmap[0]);
1012 } else {
1013 IntVect cr = domainboxes[lev-1].length() / domainboxes[lev].length();
1014 if (!m_grids[0].back().coarsenable(cr)) {
1015 break; // average_down would fail if fine boxarray is not coarsenable.
1016 }
1017 m_grids[0].emplace_back(boundboxes[lev]);
1018 IntVect max_grid_size(info.agg_grid_size);
1019 if (info.do_semicoarsening && info.max_semicoarsening_level >= lev
1020 && info.semicoarsening_direction != -1)
1021 {
1022 IntVect blen = amrex::enclosedCells(boundboxes[lev]).size();
1023 AMREX_D_TERM(int mgs_0 = (max_grid_size[0]+blen[0]-1) / blen[0];,
1024 int mgs_1 = (max_grid_size[1]+blen[1]-1) / blen[1];,
1025 int mgs_2 = (max_grid_size[2]+blen[2]-1) / blen[2]);
1026 max_grid_size[info.semicoarsening_direction]
1027 *= AMREX_D_TERM(mgs_0, *mgs_1, *mgs_2);
1028 }
1029 m_grids[0].back().maxSize(max_grid_size);
1030 m_dmap[0].push_back(DistributionMapping());
1031 if (!agged) {
1032 agged = true;
1033 agg_lev = lev;
1034 }
1035 }
1036 m_geom[0].emplace_back(domainboxes[lev],rb,coord,is_per);
1037 }
1038 }
1039 else
1040 {
1041 Long consolidation_threshold = 0;
1042 Real avg_npts = 0.0;
1043 if (info.do_consolidation) {
1044 avg_npts = static_cast<Real>(a_grids[0].d_numPts()) / static_cast<Real>(ParallelContext::NProcsSub());
1045 consolidation_threshold = AMREX_D_TERM(Long(info.con_grid_size),
1046 *info.con_grid_size,
1047 *info.con_grid_size);
1048 }
1049
1050 Box const& dom0 = a_geom[0].Domain();
1051 IntVect rr_vec(1);
1052 int numsclevs = 0;
1053 for (int lev = 0; lev < info.max_coarsening_level; ++lev)
1054 {
1055 IntVect rr_level = mg_coarsen_ratio_v;
1056 bool do_semicoarsening_level = info.do_semicoarsening
1057 && numsclevs < info.max_semicoarsening_level;
1058 if (do_semicoarsening_level
1059 && info.semicoarsening_direction != -1)
1060 {
1061 rr_level[info.semicoarsening_direction] = 1;
1062 }
1063 IntVect is_coarsenable;
1064 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1065 IntVect rr_dir(1);
1066 rr_dir[idim] = rr_vec[idim] * rr_level[idim];
1067 is_coarsenable[idim] = dom0.coarsenable(rr_dir, mg_domain_min_width_v)
1068 && a_grids[0].coarsenable(rr_dir, mg_box_min_width_v);
1069 if (!is_coarsenable[idim] && do_semicoarsening_level
1070 && info.semicoarsening_direction == -1)
1071 {
1072 is_coarsenable[idim] = true;
1073 rr_level[idim] = 1;
1074 }
1075 }
1076 if (is_coarsenable != IntVect(1) || rr_level == IntVect(1)) {
1077 break;
1078 }
1079 if (do_semicoarsening_level && info.semicoarsening_direction == -1) {
1080 // make sure there is at most one direction that is not coarsened
1081 int n_ones = AMREX_D_TERM( static_cast<int>(rr_level[0] == 1),
1082 + static_cast<int>(rr_level[1] == 1),
1083 + static_cast<int>(rr_level[2] == 1));
1084 if (n_ones > 1) { break; }
1085 }
1086 if (rr_level != mg_coarsen_ratio_v) {
1087 ++numsclevs;
1088 }
1089 rr_vec *= rr_level;
1090
1091 m_geom[0].emplace_back(amrex::coarsen(dom0, rr_vec), rb, coord, is_per);
1092 m_grids[0].push_back(amrex::coarsen(a_grids[0], rr_vec));
1093
1094 if (info.do_consolidation)
1095 {
1096 if (avg_npts/static_cast<Real>(AMREX_D_TERM(rr_vec[0], *rr_vec[1], *rr_vec[2]))
1097 < Real(0.999)*static_cast<Real>(consolidation_threshold))
1098 {
1099 coned = true;
1100 con_lev = m_dmap[0].size();
1101 m_dmap[0].push_back(DistributionMapping());
1102 }
1103 else
1104 {
1105 m_dmap[0].push_back(m_dmap[0].back());
1106 }
1107 }
1108 else
1109 {
1110 m_dmap[0].push_back(a_dmap[0]);
1111 }
1112 }
1113 }
1114
1115 m_num_mg_levels[0] = m_grids[0].size();
1116
1117 for (int mglev = 0; mglev < m_num_mg_levels[0] - 1; mglev++){
1118 const Box& fine_domain = m_geom[0][mglev].Domain();
1119 const Box& crse_domain = m_geom[0][mglev+1].Domain();
1120 mg_coarsen_ratio_vec.push_back(fine_domain.length()/crse_domain.length());
1121 }
1122
1123 for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) {
1124 if (AMRRefRatio(amrlev) == 4 && mg_coarsen_ratio_vec.empty()) {
1125 mg_coarsen_ratio_vec.push_back(IntVect(2));
1126 }
1127 }
1128
1129 if (agged)
1130 {
1131 makeAgglomeratedDMap(m_grids[0], m_dmap[0]);
1132 }
1133 else if (coned)
1134 {
1135 makeConsolidatedDMap(m_grids[0], m_dmap[0], info.con_ratio, info.con_strategy);
1136 }
1137
1138 if (agged || coned)
1139 {
1140 m_bottom_comm = makeSubCommunicator(m_dmap[0].back());
1141 }
1142 else
1143 {
1144 m_bottom_comm = m_default_comm;
1145 }
1146
1147 m_do_agglomeration = agged;
1148 m_do_consolidation = coned;
1149
1150 if (verbose > 1) {
1151 if (agged) {
1152 Print() << "MLLinOp::defineGrids(): agglomerated AMR level 0 starting at MG level "
1153 << agg_lev << " of " << m_num_mg_levels[0] << "\n";
1154 } else if (coned) {
1155 Print() << "MLLinOp::defineGrids(): consolidated AMR level 0 starting at MG level "
1156 << con_lev << " of " << m_num_mg_levels[0]
1157 << " (ratio = " << info.con_ratio << ")" << "\n";
1158 } else {
1159 Print() << "MLLinOp::defineGrids(): no agglomeration or consolidation of AMR level 0\n";
1160 }
1161 }
1162
1163 for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
1164 {
1165 for (int mglev = 1; mglev < m_num_mg_levels[amrlev]; ++mglev)
1166 {
1167 m_factory[amrlev].emplace_back(makeFactory(amrlev,mglev));
1168 }
1169 }
1170
1171 for (int amrlev = 1; amrlev < m_num_amr_levels; ++amrlev)
1172 {
1173 AMREX_ASSERT_WITH_MESSAGE(m_grids[amrlev][0].coarsenable(m_amr_ref_ratio[amrlev-1]),
1174 "MLLinOp: grids not coarsenable between AMR levels");
1175 }
1176}
1177
1178template <typename MF>
1179void
1181{
1182 m_needs_coarse_data_for_bc = !m_domain_covered[0];
1183
1184 levelbc_raii.resize(m_num_amr_levels);
1185 robin_a_raii.resize(m_num_amr_levels);
1186 robin_b_raii.resize(m_num_amr_levels);
1187 robin_f_raii.resize(m_num_amr_levels);
1188}
1189
1190template <typename MF>
1191void
1193 const Array<BCType,AMREX_SPACEDIM>& a_hibc) noexcept
1194{
1195 const int ncomp = getNComp();
1196 setDomainBC(Vector<Array<BCType,AMREX_SPACEDIM> >(ncomp,a_lobc),
1197 Vector<Array<BCType,AMREX_SPACEDIM> >(ncomp,a_hibc));
1198}
1199
1200template <typename MF>
1201void
1203 const Vector<Array<BCType,AMREX_SPACEDIM> >& a_hibc) noexcept
1204{
1205 const int ncomp = getNComp();
1206 AMREX_ASSERT_WITH_MESSAGE(ncomp == a_lobc.size() && ncomp == a_hibc.size(),
1207 "MLLinOp::setDomainBC: wrong size");
1208 m_lobc = a_lobc;
1209 m_hibc = a_hibc;
1210 m_lobc_orig = m_lobc;
1211 m_hibc_orig = m_hibc;
1212 for (int icomp = 0; icomp < ncomp; ++icomp) {
1213 for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
1214 if (m_geom[0][0].isPeriodic(idim)) {
1215 AMREX_ALWAYS_ASSERT(m_lobc[icomp][idim] == BCType::Periodic &&
1216 m_hibc[icomp][idim] == BCType::Periodic);
1217 } else {
1218 AMREX_ALWAYS_ASSERT(m_lobc[icomp][idim] != BCType::Periodic &&
1219 m_hibc[icomp][idim] != BCType::Periodic);
1220 }
1221
1222 if (m_lobc[icomp][idim] == LinOpBCType::inhomogNeumann ||
1223 m_lobc[icomp][idim] == LinOpBCType::Robin)
1224 {
1225 m_lobc[icomp][idim] = LinOpBCType::Neumann;
1226 }
1227
1228 if (m_hibc[icomp][idim] == LinOpBCType::inhomogNeumann ||
1229 m_hibc[icomp][idim] == LinOpBCType::Robin)
1230 {
1231 m_hibc[icomp][idim] = LinOpBCType::Neumann;
1232 }
1233 }
1234 }
1235
1236 if (hasHiddenDimension()) {
1237 const int hd = hiddenDirection();
1238 for (int n = 0; n < ncomp; ++n) {
1239 m_lobc[n][hd] = LinOpBCType::Neumann;
1240 m_hibc[n][hd] = LinOpBCType::Neumann;
1241 }
1242 }
1243
1244 if (hasInhomogNeumannBC() && !supportInhomogNeumannBC()) {
1245 amrex::Abort("Inhomogeneous Neumann BC not supported");
1246 }
1247 if (hasRobinBC() && !supportRobinBC()) {
1248 amrex::Abort("Robin BC not supported");
1249 }
1250}
1251
1252template <typename MF>
1253bool
1254MLLinOpT<MF>::hasBC (BCType bct) const noexcept
1255{
1256 int ncomp = m_lobc_orig.size();
1257 for (int n = 0; n < ncomp; ++n) {
1258 for (int idim = 0; idim <AMREX_SPACEDIM; ++idim) {
1259 if (m_lobc_orig[n][idim] == bct || m_hibc_orig[n][idim] == bct) {
1260 return true;
1261 }
1262 }
1263 }
1264 return false;
1265}
1266
1267template <typename MF>
1268bool
1270{
1271 return hasBC(BCType::inhomogNeumann);
1272}
1273
1274template <typename MF>
1275bool
1277{
1278 return hasBC(BCType::Robin);
1279}
1280
1281template <typename MF>
1282Box
1283MLLinOpT<MF>::compactify (Box const& b) const noexcept
1284{
1285#if (AMREX_SPACEDIM == 3)
1286 if (info.hasHiddenDimension()) {
1287 const auto& lo = b.smallEnd();
1288 const auto& hi = b.bigEnd();
1289 if (info.hidden_direction == 0) {
1290 return Box(IntVect(lo[1],lo[2],0), IntVect(hi[1],hi[2],0), b.ixType());
1291 } else if (info.hidden_direction == 1) {
1292 return Box(IntVect(lo[0],lo[2],0), IntVect(hi[0],hi[2],0), b.ixType());
1293 } else {
1294 return Box(IntVect(lo[0],lo[1],0), IntVect(hi[0],hi[1],0), b.ixType());
1295 }
1296 } else
1297#endif
1298 {
1299 return b;
1300 }
1301}
1302
1303template <typename MF>
1304void
1307{
1308 BL_PROFILE("MLLinOp::makeAgglomeratedDMap");
1309
1310 BL_ASSERT(!dm[0].empty());
1311 for (int i = 1, N=static_cast<int>(ba.size()); i < N; ++i)
1312 {
1313 if (dm[i].empty())
1314 {
1315 const std::vector< std::vector<int> >& sfc = DistributionMapping::makeSFC(ba[i]);
1316
1317 const int nprocs = ParallelContext::NProcsSub();
1318 AMREX_ASSERT(static_cast<int>(sfc.size()) == nprocs);
1319
1320 Vector<int> pmap(ba[i].size());
1321 for (int iproc = 0; iproc < nprocs; ++iproc) {
1322 int grank = ParallelContext::local_to_global_rank(iproc);
1323 for (int ibox : sfc[iproc]) {
1324 pmap[ibox] = grank;
1325 }
1326 }
1327 dm[i].define(std::move(pmap));
1328 }
1329 }
1330}
1331
1332template <typename MF>
1333void
1336 int ratio, int strategy)
1337{
1338 BL_PROFILE("MLLinOp::makeConsolidatedDMap()");
1339
1340 int factor = 1;
1341 BL_ASSERT(!dm[0].empty());
1342 for (int i = 1, N=static_cast<int>(ba.size()); i < N; ++i)
1343 {
1344 if (dm[i].empty())
1345 {
1346 factor *= ratio;
1347
1348 const int nprocs = ParallelContext::NProcsSub();
1349 const auto& pmap_fine = dm[i-1].ProcessorMap();
1350 Vector<int> pmap(pmap_fine.size());
1351 ParallelContext::global_to_local_rank(pmap.data(), pmap_fine.data(), static_cast<int>(pmap.size()));
1352 if (strategy == 1) {
1353 for (auto& x: pmap) {
1354 x /= ratio;
1355 }
1356 } else if (strategy == 2) {
1357 int nprocs_con = static_cast<int>(std::ceil(static_cast<Real>(nprocs)
1358 / static_cast<Real>(factor)));
1359 for (auto& x: pmap) {
1360 auto d = std::div(x,nprocs_con);
1361 x = d.rem;
1362 }
1363 } else if (strategy == 3) {
1364 if (factor == ratio) {
1365 const std::vector< std::vector<int> >& sfc = DistributionMapping::makeSFC(ba[i]);
1366 for (int iproc = 0; iproc < nprocs; ++iproc) {
1367 for (int ibox : sfc[iproc]) {
1368 pmap[ibox] = iproc;
1369 }
1370 }
1371 }
1372 for (auto& x: pmap) {
1373 x /= ratio;
1374 }
1375 }
1376
1378 dm[i].define(std::move(pmap));
1379 } else {
1380 Vector<int> pmap_g(pmap.size());
1381 ParallelContext::local_to_global_rank(pmap_g.data(), pmap.data(), static_cast<int>(pmap.size()));
1382 dm[i].define(std::move(pmap_g));
1383 }
1384 }
1385 }
1386}
1387
1388template <typename MF>
1391{
1392 BL_PROFILE("MLLinOp::makeSubCommunicator()");
1393
1394#ifdef BL_USE_MPI
1395
1396 Vector<int> newgrp_ranks = dm.ProcessorMap();
1397 std::sort(newgrp_ranks.begin(), newgrp_ranks.end());
1398 auto last = std::unique(newgrp_ranks.begin(), newgrp_ranks.end());
1399 newgrp_ranks.erase(last, newgrp_ranks.end());
1400
1401 MPI_Comm newcomm;
1402 MPI_Group defgrp, newgrp;
1403 MPI_Comm_group(m_default_comm, &defgrp);
1405 MPI_Group_incl(defgrp, static_cast<int>(newgrp_ranks.size()), newgrp_ranks.data(), &newgrp);
1406 } else {
1407 Vector<int> local_newgrp_ranks(newgrp_ranks.size());
1408 ParallelContext::global_to_local_rank(local_newgrp_ranks.data(),
1409 newgrp_ranks.data(), static_cast<int>(newgrp_ranks.size()));
1410 MPI_Group_incl(defgrp, static_cast<int>(local_newgrp_ranks.size()), local_newgrp_ranks.data(), &newgrp);
1411 }
1412
1413 MPI_Comm_create(m_default_comm, newgrp, &newcomm);
1414
1415 m_raii_comm = std::make_unique<CommContainer>(newcomm);
1416
1417 MPI_Group_free(&defgrp);
1418 MPI_Group_free(&newgrp);
1419
1420 return newcomm;
1421#else
1423 return m_default_comm;
1424#endif
1425}
1426
1427template <typename MF>
1428void
1430 const Array<Real,AMREX_SPACEDIM>& hi_bcloc) noexcept
1431{
1432 m_domain_bloc_lo = lo_bcloc;
1433 m_domain_bloc_hi = hi_bcloc;
1434}
1435
1436template <typename MF>
1437void
1438MLLinOpT<MF>::setCoarseFineBC (const MF* crse, int crse_ratio,
1439 LinOpBCType bc_type) noexcept
1440{
1441 setCoarseFineBC(crse, IntVect(crse_ratio), bc_type);
1442}
1443
1444template <typename MF>
1445void
1446MLLinOpT<MF>::setCoarseFineBC (const MF* crse, IntVect const& crse_ratio,
1447 LinOpBCType bc_type) noexcept
1448{
1449 m_coarse_data_for_bc = crse;
1450 m_coarse_data_crse_ratio = crse_ratio;
1451 m_coarse_fine_bc_type = bc_type;
1452}
1453
1454template <typename MF>
1455template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int>>
1456void
1457MLLinOpT<MF>::setCoarseFineBC (const AMF* crse, int crse_ratio,
1458 LinOpBCType bc_type) noexcept
1459{
1460 setCoarseFineBC(crse, IntVect(crse_ratio), bc_type);
1461}
1462
1463template <typename MF>
1464template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int>>
1465void
1466MLLinOpT<MF>::setCoarseFineBC (const AMF* crse, IntVect const& crse_ratio,
1467 LinOpBCType bc_type) noexcept
1468{
1469 m_coarse_data_for_bc_raii = MF(crse->boxArray(), crse->DistributionMap(),
1470 crse->nComp(), crse->nGrowVect());
1471 m_coarse_data_for_bc_raii.LocalCopy(*crse, 0, 0, crse->nComp(),
1472 crse->nGrowVect());
1473 m_coarse_data_for_bc = &m_coarse_data_for_bc_raii;
1474 m_coarse_data_crse_ratio = crse_ratio;
1475 m_coarse_fine_bc_type = bc_type;
1476}
1477
1478template <typename MF>
1479void
1481{
1482 mf.clear();
1483 mf.resize(m_num_amr_levels);
1484 for (int alev = 0; alev < m_num_amr_levels; ++alev) {
1485 mf[alev].resize(m_num_mg_levels[alev]);
1486 for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) {
1487 mf[alev][mlev] = make(alev, mlev, ng);
1488 }
1489 }
1490}
1491
1492template <typename MF>
1493MF
1494MLLinOpT<MF>::make (int amrlev, int mglev, IntVect const& ng) const
1495{
1496 if constexpr (IsMultiFabLike_v<MF>) {
1497 return MF(amrex::convert(m_grids[amrlev][mglev], m_ixtype),
1498 m_dmap[amrlev][mglev], getNComp(), ng, MFInfo(),
1499 *m_factory[amrlev][mglev]);
1500 } else {
1501 amrex::ignore_unused(amrlev, mglev, ng);
1502 amrex::Abort("MLLinOpT::make: how did we get here?");
1503 return {};
1504 }
1505}
1506
1507template <typename MF>
1508MF
1509MLLinOpT<MF>::makeAlias (MF const& mf) const
1510{
1511 if constexpr (IsMultiFabLike_v<MF>) {
1512 return MF(mf, amrex::make_alias, 0, mf.nComp());
1513 } else {
1515 amrex::Abort("MLLinOpT::makeAlias: how did we get here?");
1516 return {};
1517 }
1518}
1519
1520template <typename MF>
1521MF
1522MLLinOpT<MF>::makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const
1523{
1524 if constexpr (IsMultiFabLike_v<MF>) {
1525 BoxArray cba = m_grids[amrlev][mglev];
1526 IntVect ratio = (amrlev > 0) ? IntVect(2) : mg_coarsen_ratio_vec[mglev];
1527 cba.coarsen(ratio);
1528 cba.convert(m_ixtype);
1529 return MF(cba, m_dmap[amrlev][mglev], getNComp(), ng);
1530 } else {
1531 amrex::ignore_unused(amrlev, mglev, ng);
1532 amrex::Abort("MLLinOpT::makeCoarseMG: how did we get here?");
1533 return {};
1534 }
1535}
1536
1537template <typename MF>
1538MF
1539MLLinOpT<MF>::makeCoarseAmr (int famrlev, IntVect const& ng) const
1540{
1541 if constexpr (IsMultiFabLike_v<MF>) {
1542 BoxArray cba = m_grids[famrlev][0];
1543 IntVect ratio(AMRRefRatio(famrlev-1));
1544 cba.coarsen(ratio);
1545 cba.convert(m_ixtype);
1546 return MF(cba, m_dmap[famrlev][0], getNComp(), ng);
1547 } else {
1548 amrex::ignore_unused(famrlev, ng);
1549 amrex::Abort("MLLinOpT::makeCoarseAmr: how did we get here?");
1550 return {};
1551 }
1552}
1553
1554template <typename MF>
1555void
1557{
1558 if (new_size <= 0 || new_size >= m_num_mg_levels[0]) { return; }
1559
1560 m_num_mg_levels[0] = new_size;
1561
1562 m_geom[0].resize(new_size);
1563 m_grids[0].resize(new_size);
1564 m_dmap[0].resize(new_size);
1565 m_factory[0].resize(new_size);
1566
1567 if (m_bottom_comm != m_default_comm) {
1568 m_bottom_comm = makeSubCommunicator(m_dmap[0].back());
1569 }
1570}
1571
1572template <typename MF>
1573void
1574MLLinOpT<MF>::avgDownResMG (int clev, MF& cres, MF const& fres) const
1575{
1576 amrex::ignore_unused(clev, cres, fres);
1577 if constexpr (amrex::IsFabArray<MF>::value) {
1578 const int ncomp = this->getNComp();
1579#ifdef AMREX_USE_EB
1580 if (!fres.isAllRegular()) {
1581 if constexpr (std::is_same<MF,MultiFab>()) {
1582 amrex::EB_average_down(fres, cres, 0, ncomp,
1583 mg_coarsen_ratio_vec[clev-1]);
1584 } else {
1585 amrex::Abort("EB_average_down only works with MultiFab");
1586 }
1587 } else
1588#endif
1589 {
1590 amrex::average_down(fres, cres, 0, ncomp, mg_coarsen_ratio_vec[clev-1]);
1591 }
1592 } else {
1593 amrex::Abort("For non-FabArray, MLLinOpT<MF>::avgDownResMG should be overridden.");
1594 }
1595}
1596
1597template <typename MF>
1598bool
1599MLLinOpT<MF>::isMFIterSafe (int amrlev, int mglev1, int mglev2) const
1600{
1601 return m_dmap[amrlev][mglev1] == m_dmap[amrlev][mglev2]
1602 && BoxArray::SameRefs(m_grids[amrlev][mglev1], m_grids[amrlev][mglev2]);
1603}
1604
1605template <typename MF>
1606template <typename AMF, std::enable_if_t<!std::is_same_v<MF,AMF>,int>>
1607void
1608MLLinOpT<MF>::setLevelBC (int amrlev, const AMF* levelbcdata,
1609 const AMF* robinbc_a, const AMF* robinbc_b,
1610 const AMF* robinbc_f)
1611{
1612 const int ncomp = this->getNComp();
1613 if (levelbcdata) {
1614 levelbc_raii[amrlev] = std::make_unique<MF>(levelbcdata->boxArray(),
1615 levelbcdata->DistributionMap(),
1616 ncomp, levelbcdata->nGrowVect());
1617 levelbc_raii[amrlev]->LocalCopy(*levelbcdata, 0, 0, ncomp,
1618 levelbcdata->nGrowVect());
1619 } else {
1620 levelbc_raii[amrlev].reset();
1621 }
1622
1623 if (robinbc_a) {
1624 robin_a_raii[amrlev] = std::make_unique<MF>(robinbc_a->boxArray(),
1625 robinbc_a->DistributionMap(),
1626 ncomp, robinbc_a->nGrowVect());
1627 robin_a_raii[amrlev]->LocalCopy(*robinbc_a, 0, 0, ncomp,
1628 robinbc_a->nGrowVect());
1629 } else {
1630 robin_a_raii[amrlev].reset();
1631 }
1632
1633 if (robinbc_b) {
1634 robin_b_raii[amrlev] = std::make_unique<MF>(robinbc_b->boxArray(),
1635 robinbc_b->DistributionMap(),
1636 ncomp, robinbc_b->nGrowVect());
1637 robin_b_raii[amrlev]->LocalCopy(*robinbc_b, 0, 0, ncomp,
1638 robinbc_b->nGrowVect());
1639 } else {
1640 robin_b_raii[amrlev].reset();
1641 }
1642
1643 if (robinbc_f) {
1644 robin_f_raii[amrlev] = std::make_unique<MF>(robinbc_f->boxArray(),
1645 robinbc_f->DistributionMap(),
1646 ncomp, robinbc_f->nGrowVect());
1647 robin_f_raii[amrlev]->LocalCopy(*robinbc_f, 0, 0, ncomp,
1648 robinbc_f->nGrowVect());
1649 } else {
1650 robin_f_raii[amrlev].reset();
1651 }
1652
1653 this->setLevelBC(amrlev, levelbc_raii[amrlev].get(), robin_a_raii[amrlev].get(),
1654 robin_b_raii[amrlev].get(), robin_f_raii[amrlev].get());
1655}
1656
1657template <typename MF>
1658auto
1660{
1661 AMREX_ALWAYS_ASSERT(NAMRLevels() == 1);
1662 return xdoty(0,0,*x[0],*y[0],false);
1663}
1664
1665template <typename MF>
1666auto
1668{
1669 AMREX_ALWAYS_ASSERT(NAMRLevels() == 1);
1670 auto r = xdoty(0,0,*x[0],*x[0],false);
1671 return std::sqrt(r);
1672}
1673
1674extern template class MLLinOpT<MultiFab>;
1675
1677
1678}
1679
1680#endif
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define BL_ASSERT(EX)
Definition AMReX_BLassert.H:39
#define AMREX_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:37
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
#define AMREX_ALWAYS_ASSERT(EX)
Definition AMReX_BLassert.H:50
Array4< Real > fine
Definition AMReX_InterpFaceRegister.cpp:90
Array4< Real const > crse
Definition AMReX_InterpFaceRegister.cpp:92
#define AMREX_D_TERM(a, b, c)
Definition AMReX_SPACE.H:129
#define AMREX_D_PICK(a, b, c)
Definition AMReX_SPACE.H:151
int MPI_Comm
Definition AMReX_ccse-mpi.H:47
int MPI_Group
Definition AMReX_ccse-mpi.H:48
static constexpr int MPI_COMM_NULL
Definition AMReX_ccse-mpi.H:55
A collection of Boxes stored in an Array.
Definition AMReX_BoxArray.H:550
static bool SameRefs(const BoxArray &lhs, const BoxArray &rhs)
whether two BoxArrays share the same data
Definition AMReX_BoxArray.H:820
BoxArray & coarsen(int refinement_ratio)
Coarsen each Box in the BoxArray to the specified ratio.
BoxArray & convert(IndexType typ)
Apply Box::convert(IndexType) to each Box in the BoxArray.
AMREX_GPU_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:708
AMREX_GPU_HOST_DEVICE BoxND & makeSlab(int direction, int slab_index) noexcept
Definition AMReX_Box.H:779
AMREX_GPU_HOST_DEVICE double d_numPts() const noexcept
Returns the number of points contained in the BoxND. This is intended for use only in diagnostic mess...
Definition AMReX_Box.H:366
AMREX_GPU_HOST_DEVICE IntVectND< dim > length() const noexcept
Return the length of the BoxND.
Definition AMReX_Box.H:146
AMREX_GPU_HOST_DEVICE bool coarsenable(const IntVectND< dim > &refrat, const IntVectND< dim > &min_width) const noexcept
Definition AMReX_Box.H:745
AMREX_GPU_HOST_DEVICE const IntVectND< dim > & smallEnd() const &noexcept
Get the smallend of the BoxND.
Definition AMReX_Box.H:105
AMREX_GPU_HOST_DEVICE Long numPts() const noexcept
Returns the number of points contained in the BoxND.
Definition AMReX_Box.H:346
Definition AMReX_FabFactory.H:76
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:41
static DistributionMapping makeSFC(const MultiFab &weight, bool sort=true)
Definition AMReX_DistributionMapping.cpp:1764
const Vector< int > & ProcessorMap() const noexcept
Returns a constant reference to the mapping of boxes in the underlying BoxArray to the CPU that holds...
Definition AMReX_DistributionMapping.cpp:47
Definition AMReX_EBFabFactory.H:24
Definition AMReX_FabFactory.H:50
Solve using GMRES with multigrid as preconditioner.
Definition AMReX_GMRES_MLMG.H:21
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:73
Interface
Definition AMReX_Hypre.H:21
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr std::size_t size() noexcept
Definition AMReX_IntVect.H:725
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE IntVectND< new_dim > resize(int fill_extra=0) const noexcept
Returns a new IntVectND of size new_dim by either shrinking or expanding this IntVectND.
Definition AMReX_IntVect.H:768
Definition AMReX_MLABecLaplacian.H:15
Definition AMReX_MLCGSolver.H:12
Definition AMReX_MLLinOp.H:98
const MF * m_coarse_data_for_bc
Definition AMReX_MLLinOp.H:637
virtual void smooth(int amrlev, int mglev, MF &sol, const MF &rhs, bool skip_fillboundary=false) const =0
Smooth.
void setDomainBC(const Array< BCType, AMREX_SPACEDIM > &lobc, const Array< BCType, AMREX_SPACEDIM > &hibc) noexcept
Boundary of the whole domain.
Definition AMReX_MLLinOp.H:1192
Vector< Vector< std::unique_ptr< FabFactory< FAB > > > > m_factory
Definition AMReX_MLLinOp.H:610
virtual void avgDownResMG(int clev, MF &cres, MF const &fres) const
Definition AMReX_MLLinOp.H:1574
int NAMRLevels() const noexcept
Return the number of AMR levels.
Definition AMReX_MLLinOp.H:563
bool m_do_consolidation
Definition AMReX_MLLinOp.H:601
virtual void compGrad(int, const Array< MF *, AMREX_SPACEDIM > &, MF &, Location) const
Compute gradients of the solution.
Definition AMReX_MLLinOp.H:437
bool isCellCentered() const noexcept
Definition AMReX_MLLinOp.H:684
IntVect m_ixtype
Definition AMReX_MLLinOp.H:598
void setVerbose(int v) noexcept
Set verbosity.
Definition AMReX_MLLinOp.H:244
Array< Real, AMREX_SPACEDIM > m_domain_bloc_hi
Definition AMReX_MLLinOp.H:631
bool isMFIterSafe(int amrlev, int mglev1, int mglev2) const
Definition AMReX_MLLinOp.H:1599
RealVect m_coarse_bc_loc
Definition AMReX_MLLinOp.H:636
virtual bool needsUpdate() const
Does it need update if it's reused?
Definition AMReX_MLLinOp.H:266
virtual void interpolation(int amrlev, int fmglev, MF &fine, const MF &crse) const =0
Add interpolated coarse MG level data to fine MG level data.
virtual void setLevelBC(int, const MF *, const MF *=nullptr, const MF *=nullptr, const MF *=nullptr)=0
Set boundary conditions for given level. For cell-centered solves only.
virtual MF make(int amrlev, int mglev, IntVect const &ng) const
Definition AMReX_MLLinOp.H:1494
FabFactory< FAB > const * Factory(int amr_lev, int mglev=0) const noexcept
Definition AMReX_MLLinOp.H:648
T get_d0(T const &dx, T const &dy, T const &) const noexcept
Definition AMReX_MLLinOp.H:725
MPI_Comm BottomCommunicator() const noexcept
Definition AMReX_MLLinOp.H:675
void setEnforceSingularSolvable(bool o) noexcept
Definition AMReX_MLLinOp.H:253
MPI_Comm Communicator() const noexcept
Definition AMReX_MLLinOp.H:676
virtual void getFluxes(const Vector< Array< MF *, AMREX_SPACEDIM > > &, const Vector< MF * > &, Location) const
Definition AMReX_MLLinOp.H:500
int mg_domain_min_width
Definition AMReX_MLLinOp.H:582
void setMaxOrder(int o) noexcept
Set order of interpolation at coarse/fine boundary.
Definition AMReX_MLLinOp.H:247
Vector< Array< BCType, AMREX_SPACEDIM > > m_hibc_orig
Definition AMReX_MLLinOp.H:576
virtual void interpAssign(int amrlev, int fmglev, MF &fine, MF &crse) const
Overwrite fine MG level data with interpolated coarse data.
Definition AMReX_MLLinOp.H:299
virtual std::string name() const
Definition AMReX_MLLinOp.H:131
void setCoarseFineBC(const AMF *crse, IntVect const &crse_ratio, LinOpBCType bc_type=LinOpBCType::Dirichlet) noexcept
Definition AMReX_MLLinOp.H:1466
bool doAgglomeration() const noexcept
Definition AMReX_MLLinOp.H:680
Vector< std::unique_ptr< MF > > levelbc_raii
Definition AMReX_MLLinOp.H:760
MF m_coarse_data_for_bc_raii
Definition AMReX_MLLinOp.H:638
virtual void setDirichletNodesToZero(int, int, MF &) const
Definition AMReX_MLLinOp.H:476
MLLinOpT< MF > & operator=(const MLLinOpT< MF > &)=delete
std::unique_ptr< CommContainer > m_raii_comm
Definition AMReX_MLLinOp.H:628
bool m_do_semicoarsening
Definition AMReX_MLLinOp.H:603
bool hasRobinBC() const noexcept
Definition AMReX_MLLinOp.H:1276
Array< Real, AMREX_SPACEDIM > m_domain_bloc_lo
Definition AMReX_MLLinOp.H:630
virtual void resizeMultiGrid(int new_size)
Definition AMReX_MLLinOp.H:1556
Vector< Vector< BoxArray > > m_grids
Definition AMReX_MLLinOp.H:608
virtual MF makeCoarseAmr(int famrlev, IntVect const &ng) const
Definition AMReX_MLLinOp.H:1539
Vector< std::unique_ptr< MF > > robin_f_raii
Definition AMReX_MLLinOp.H:763
bool m_do_agglomeration
Definition AMReX_MLLinOp.H:600
void setCoarseFineBC(const AMF *crse, int crse_ratio, LinOpBCType bc_type=LinOpBCType::Dirichlet) noexcept
Definition AMReX_MLLinOp.H:1457
virtual MF makeAlias(MF const &mf) const
Definition AMReX_MLLinOp.H:1509
void setDomainBC(const Vector< Array< BCType, AMREX_SPACEDIM > > &lobc, const Vector< Array< BCType, AMREX_SPACEDIM > > &hibc) noexcept
Boundary of the whole domain.
Definition AMReX_MLLinOp.H:1202
Array4< T > compactify(Array4< T > const &a) const noexcept
Definition AMReX_MLLinOp.H:711
static constexpr int mg_coarsen_ratio
Definition AMReX_MLLinOp.H:580
virtual void copyNSolveSolution(MF &, MF const &) const
Definition AMReX_MLLinOp.H:540
virtual void solutionResidual(int amrlev, MF &resid, MF &x, const MF &b, const MF *crse_bcdata=nullptr)=0
Compute residual for solution.
int getMaxOrder() const noexcept
Get order of interpolation at coarse/fine boundary.
Definition AMReX_MLLinOp.H:249
virtual int getNComp() const
Return number of components.
Definition AMReX_MLLinOp.H:261
virtual void getFluxes(const Vector< MF * > &, const Vector< MF * > &) const
Definition AMReX_MLLinOp.H:505
void defineBC()
Definition AMReX_MLLinOp.H:1180
void setCoarseFineBCLocation(const RealVect &cloc) noexcept
Definition AMReX_MLLinOp.H:678
Vector< int > m_amr_ref_ratio
Definition AMReX_MLLinOp.H:593
MPI_Comm m_default_comm
Definition AMReX_MLLinOp.H:613
static void makeAgglomeratedDMap(const Vector< BoxArray > &ba, Vector< DistributionMapping > &dm)
Definition AMReX_MLLinOp.H:1305
bool isBottomActive() const noexcept
Definition AMReX_MLLinOp.H:673
virtual Vector< RT > getSolvabilityOffset(int, int, MF const &) const
get offset for fixing solvability
Definition AMReX_MLLinOp.H:463
void defineGrids(const Vector< Geometry > &a_geom, const Vector< BoxArray > &a_grids, const Vector< DistributionMapping > &a_dmap, const Vector< FabFactory< FAB > const * > &a_factory)
Definition AMReX_MLLinOp.H:806
virtual BottomSolver getDefaultBottomSolver() const
Definition AMReX_MLLinOp.H:258
typename FabDataType< MF >::fab_type FAB
Definition AMReX_MLLinOp.H:108
virtual RT normInf(int amrlev, MF const &mf, bool local) const =0
Vector< int > m_num_mg_levels
Definition AMReX_MLLinOp.H:595
bool hasBC(BCType bct) const noexcept
Definition AMReX_MLLinOp.H:1254
static void makeConsolidatedDMap(const Vector< BoxArray > &ba, Vector< DistributionMapping > &dm, int ratio, int strategy)
Definition AMReX_MLLinOp.H:1334
Vector< Vector< DistributionMapping > > m_dmap
Definition AMReX_MLLinOp.H:609
int verbose
Definition AMReX_MLLinOp.H:586
virtual void checkPoint(std::string const &) const
Definition AMReX_MLLinOp.H:756
IntVect m_coarse_data_crse_ratio
Definition AMReX_MLLinOp.H:635
LinOpBCType BCType
Definition AMReX_MLLinOp.H:111
virtual void correctionResidual(int amrlev, int mglev, MF &resid, MF &x, const MF &b, BCMode bc_mode, const MF *crse_bcdata=nullptr)=0
Compute residual for the residual-correction form, resid = b - L(x)
virtual void compFlux(int, const Array< MF *, AMREX_SPACEDIM > &, MF &, Location) const
Compute fluxes.
Definition AMReX_MLLinOp.H:423
const Vector< int > & AMRRefRatio() const noexcept
Return AMR refinement ratios.
Definition AMReX_MLLinOp.H:643
MLLinOpT(MLLinOpT< MF > &&)=delete
Vector< Array< BCType, AMREX_SPACEDIM > > m_lobc_orig
Definition AMReX_MLLinOp.H:575
void setCoarseFineBC(const MF *crse, int crse_ratio, LinOpBCType bc_type=LinOpBCType::Dirichlet) noexcept
Set coarse/fine boundary conditions. For cell-centered solves only.
Definition AMReX_MLLinOp.H:1438
virtual void apply(int amrlev, int mglev, MF &out, MF &in, BCMode bc_mode, StateMode s_mode, const MLMGBndryT< MF > *bndry=nullptr) const =0
Apply the linear operator, out = L(in)
GpuArray< BCType, AMREX_SPACEDIM > HiBC(int icomp=0) const noexcept
Definition AMReX_MLLinOp.H:657
void setDomainBCLoc(const Array< Real, AMREX_SPACEDIM > &lo_bcloc, const Array< Real, AMREX_SPACEDIM > &hi_bcloc) noexcept
Set location of domain boundaries.
Definition AMReX_MLLinOp.H:1429
bool needsCoarseDataForBC() const noexcept
Needs coarse data for bc?
Definition AMReX_MLLinOp.H:177
typename FabDataType< MF >::value_type RT
Definition AMReX_MLLinOp.H:109
virtual void update()
Update for reuse.
Definition AMReX_MLLinOp.H:268
bool m_precond_mode
Definition AMReX_MLLinOp.H:640
virtual std::unique_ptr< FabFactory< FAB > > makeFactory(int, int) const
Definition AMReX_MLLinOp.H:700
virtual void applyMetricTerm(int, int, MF &) const
apply metric terms if there are any
Definition AMReX_MLLinOp.H:444
virtual void unapplyMetricTerm(int, int, MF &) const
unapply metric terms if there are any
Definition AMReX_MLLinOp.H:446
bool hasHiddenDimension() const noexcept
Definition AMReX_MLLinOp.H:706
virtual bool isBottomSingular() const =0
Is the bottom of MG singular?
virtual void reflux(int crse_amrlev, MF &res, const MF &crse_sol, const MF &crse_rhs, MF &fine_res, MF &fine_sol, const MF &fine_rhs) const
Reflux at AMR coarse/fine boundary.
Definition AMReX_MLLinOp.H:406
virtual IntVect getNGrowVectRestriction() const
Definition AMReX_MLLinOp.H:686
Vector< std::unique_ptr< MF > > robin_b_raii
Definition AMReX_MLLinOp.H:762
static constexpr int mg_box_min_width
Definition AMReX_MLLinOp.H:581
virtual RT dotProductPrecond(Vector< MF const * > const &x, Vector< MF const * > const &y) const
Definition AMReX_MLLinOp.H:1659
virtual MF makeCoarseMG(int amrlev, int mglev, IntVect const &ng) const
Definition AMReX_MLLinOp.H:1522
virtual void fixSolvabilityByOffset(int, int, MF &, Vector< RT > const &) const
fix solvability by subtracting offset from RHS
Definition AMReX_MLLinOp.H:467
const Geometry & Geom(int amr_lev, int mglev=0) const noexcept
Definition AMReX_MLLinOp.H:568
virtual void endPrecondBC()
Definition AMReX_MLLinOp.H:558
Vector< Array< BCType, AMREX_SPACEDIM > > m_lobc
Definition AMReX_MLLinOp.H:571
int hiddenDirection() const noexcept
Definition AMReX_MLLinOp.H:707
Vector< int > m_domain_covered
Definition AMReX_MLLinOp.H:611
const MLLinOpT< MF > * m_parent
Definition AMReX_MLLinOp.H:596
bool doSemicoarsening() const noexcept
Definition AMReX_MLLinOp.H:682
virtual bool supportNSolve() const
Definition AMReX_MLLinOp.H:538
virtual bool supportRobinBC() const noexcept
Definition AMReX_MLLinOp.H:667
virtual void normalize(int amrlev, int mglev, MF &mf) const
Divide mf by the diagonal component of the operator. Used by bicgstab.
Definition AMReX_MLLinOp.H:363
virtual void applyInhomogNeumannTerm(int, MF &) const
Extra terms introduced when we treat inhomogeneous Nuemann BC as homogeneous.
Definition AMReX_MLLinOp.H:452
virtual void avgDownResAmr(int clev, MF &cres, MF const &fres) const
Definition AMReX_MLLinOp.H:548
Vector< Vector< Geometry > > m_geom
first Vector is for amr level and second is mg level
Definition AMReX_MLLinOp.H:607
virtual RT norm2Precond(Vector< MF const * > const &x) const
Definition AMReX_MLLinOp.H:1667
MLLinOpT(const MLLinOpT< MF > &)=delete
virtual bool scaleRHS(int, MF *) const
scale RHS to fix solvability
Definition AMReX_MLLinOp.H:458
virtual void averageDownAndSync(Vector< MF > &sol) const =0
void setCoarseFineBC(const MF *crse, IntVect const &crse_ratio, LinOpBCType bc_type=LinOpBCType::Dirichlet) noexcept
Definition AMReX_MLLinOp.H:1446
GpuArray< BCType, AMREX_SPACEDIM > LoBC(int icomp=0) const noexcept
Definition AMReX_MLLinOp.H:652
MLLinOpT()=default
Box compactify(Box const &b) const noexcept
Definition AMReX_MLLinOp.H:1283
bool m_needs_coarse_data_for_bc
Definition AMReX_MLLinOp.H:633
int maxorder
Definition AMReX_MLLinOp.H:588
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)
Definition AMReX_MLLinOp.H:768
Vector< IntVect > mg_coarsen_ratio_vec
Definition AMReX_MLLinOp.H:604
MF MFType
Definition AMReX_MLLinOp.H:107
virtual void preparePrecond()
Definition AMReX_MLLinOp.H:472
virtual ~MLLinOpT()=default
virtual void averageDownSolutionRHS(int camrlev, MF &crse_sol, MF &crse_rhs, const MF &fine_sol, const MF &fine_rhs)
Average-down data from fine AMR level to coarse AMR level.
Definition AMReX_MLLinOp.H:329
LPInfo info
Definition AMReX_MLLinOp.H:584
MPI_Comm makeSubCommunicator(const DistributionMapping &dm)
Definition AMReX_MLLinOp.H:1390
Vector< Array< BCType, AMREX_SPACEDIM > > m_hibc
Definition AMReX_MLLinOp.H:572
int NMGLevels(int amrlev) const noexcept
Return the number of MG levels at given AMR level.
Definition AMReX_MLLinOp.H:566
T get_d1(T const &, T const &dy, T const &dz) const noexcept
Definition AMReX_MLLinOp.H:735
bool enforceSingularSolvable
Definition AMReX_MLLinOp.H:590
void setLevelBC(int amrlev, const AMF *levelbcdata, const AMF *robinbc_a=nullptr, const AMF *robinbc_b=nullptr, const AMF *robinbc_f=nullptr)
Definition AMReX_MLLinOp.H:1608
virtual RT xdoty(int amrlev, int mglev, const MF &x, const MF &y, bool local) const =0
x dot y, used by the bottom solver
virtual void interpolationAmr(int famrlev, MF &fine, const MF &crse, IntVect const &nghost) const
Interpolation between AMR levels.
Definition AMReX_MLLinOp.H:313
bool doConsolidation() const noexcept
Definition AMReX_MLLinOp.H:681
virtual std::unique_ptr< MLLinOpT< MF > > makeNLinOp(int) const
Definition AMReX_MLLinOp.H:494
virtual bool supportInhomogNeumannBC() const noexcept
Definition AMReX_MLLinOp.H:668
virtual void prepareForFluxes(int, const MF *=nullptr)
Definition AMReX_MLLinOp.H:379
LinOpBCType m_coarse_fine_bc_type
Definition AMReX_MLLinOp.H:634
virtual bool isSingular(int amrlev) const =0
Is it singular on given AMR level?
virtual void postSolve(Vector< MF * > const &) const
Definition AMReX_MLLinOp.H:542
int AMRRefRatio(int amr_lev) const noexcept
Return AMR refinement ratio at given AMR level.
Definition AMReX_MLLinOp.H:646
MPI_Comm m_bottom_comm
Definition AMReX_MLLinOp.H:614
virtual void restriction(int amrlev, int cmglev, MF &crse, MF &fine) const =0
Restriction onto coarse MG level.
virtual void beginPrecondBC()
Definition AMReX_MLLinOp.H:557
virtual void applyOverset(int, MF &) const
for overset solver only
Definition AMReX_MLLinOp.H:455
virtual void unimposeNeumannBC(int, MF &) const
This is needed for our nodal projection solver.
Definition AMReX_MLLinOp.H:449
bool getEnforceSingularSolvable() const noexcept
Definition AMReX_MLLinOp.H:256
virtual void prepareForSolve()=0
virtual void make(Vector< Vector< MF > > &mf, IntVect const &ng) const
Definition AMReX_MLLinOp.H:1480
bool hasInhomogNeumannBC() const noexcept
Definition AMReX_MLLinOp.H:1269
Vector< std::unique_ptr< MF > > robin_a_raii
Definition AMReX_MLLinOp.H:761
virtual int getNGrow(int=0, int=0) const
Definition AMReX_MLLinOp.H:263
int m_num_amr_levels
Definition AMReX_MLLinOp.H:592
Definition AMReX_MLMGBndry.H:12
Definition AMReX_MLMG.H:12
Definition AMReX_MLPoisson.H:16
This class provides the user with a few print options.
Definition AMReX_Print.H:35
A Box with real dimensions. A RealBox is OK iff volume >= 0.
Definition AMReX_RealBox.H:21
A Real vector in SpaceDim-dimensional space.
Definition AMReX_RealVect.H:32
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:27
Long size() const noexcept
Definition AMReX_Vector.H:50
bool notInLaunchRegion() noexcept
Definition AMReX_GpuControl.H:87
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition AMReX_ParallelContext.H:70
int local_to_global_rank(int rank) noexcept
translate between local rank and global rank
Definition AMReX_ParallelContext.H:95
int global_to_local_rank(int rank) noexcept
Definition AMReX_ParallelContext.H:98
int NProcsSub() noexcept
number of ranks in current frame
Definition AMReX_ParallelContext.H:74
MPI_Comm Communicator() noexcept
Definition AMReX_ParallelDescriptor.H:210
Definition AMReX_Amr.cpp:49
@ make_alias
Definition AMReX_MakeType.H:7
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:314
BoxND< AMREX_SPACEDIM > Box
Definition AMReX_BaseFwd.H:27
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > enclosedCells(const BoxND< dim > &b, int dir) noexcept
Returns a BoxND with CELL based coordinates in direction dir that is enclosed by b....
Definition AMReX_Box.H:1463
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Returns a BoxND with different type.
Definition AMReX_Box.H:1435
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > coarsen(const BoxND< dim > &b, int ref_ratio) noexcept
Coarsen BoxND by given (positive) refinement ratio. NOTE: if type(dir) = CELL centered: lo <- lo/rati...
Definition AMReX_Box.H:1304
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)
Definition AMReX_EBMultiFabUtil.cpp:336
IntVectND< AMREX_SPACEDIM > IntVect
Definition AMReX_BaseFwd.H:30
AMREX_GPU_HOST_DEVICE constexpr GpuTupleElement< I, GpuTuple< Ts... > >::type & get(GpuTuple< Ts... > &tup) noexcept
Definition AMReX_Tuple.H:179
BottomSolver
Definition AMReX_MLLinOp.H:30
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:111
void Warning(const std::string &msg)
Print out warning message to cerr.
Definition AMReX.cpp:231
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:225
int verbose
Definition AMReX_DistributionMapping.cpp:36
std::unique_ptr< Hypre > makeHypre(const BoxArray &grids, const DistributionMapping &dmap, const Geometry &geom, MPI_Comm comm_, Hypre::Interface interface, const iMultiFab *overset_mask)
Definition AMReX_Hypre.cpp:12
std::array< T, N > Array
Definition AMReX_Array.H:24
Definition AMReX_Array4.H:61
Definition AMReX_FabDataType.H:9
Definition AMReX_Array.H:34
Definition AMReX_TypeTraits.H:29
Definition AMReX_MLLinOp.H:35
LPInfo & setConsolidationRatio(int x) noexcept
Definition AMReX_MLLinOp.H:54
int con_strategy
Definition AMReX_MLLinOp.H:42
bool do_semicoarsening
Definition AMReX_MLLinOp.H:38
bool has_metric_term
Definition AMReX_MLLinOp.H:43
LPInfo & setConsolidationGridSize(int x) noexcept
Definition AMReX_MLLinOp.H:53
int max_semicoarsening_level
Definition AMReX_MLLinOp.H:45
bool hasHiddenDimension() const noexcept
Definition AMReX_MLLinOp.H:62
int con_ratio
Definition AMReX_MLLinOp.H:41
bool do_consolidation
Definition AMReX_MLLinOp.H:37
int con_grid_size
Definition AMReX_MLLinOp.H:40
LPInfo & setSemicoarsening(bool x) noexcept
Definition AMReX_MLLinOp.H:51
LPInfo & setHiddenDirection(int n) noexcept
Definition AMReX_MLLinOp.H:60
LPInfo & setConsolidation(bool x) noexcept
Definition AMReX_MLLinOp.H:50
LPInfo & setSemicoarseningDirection(int n) noexcept
Definition AMReX_MLLinOp.H:59
LPInfo & setMaxSemicoarseningLevel(int n) noexcept
Definition AMReX_MLLinOp.H:58
bool do_agglomeration
Definition AMReX_MLLinOp.H:36
LPInfo & setMetricTerm(bool x) noexcept
Definition AMReX_MLLinOp.H:56
static constexpr int getDefaultConsolidationGridSize()
Definition AMReX_MLLinOp.H:74
LPInfo & setConsolidationStrategy(int x) noexcept
Definition AMReX_MLLinOp.H:55
int max_coarsening_level
Definition AMReX_MLLinOp.H:44
int agg_grid_size
Definition AMReX_MLLinOp.H:39
static constexpr int getDefaultAgglomerationGridSize()
Definition AMReX_MLLinOp.H:66
int hidden_direction
Definition AMReX_MLLinOp.H:47
LPInfo & setAgglomerationGridSize(int x) noexcept
Definition AMReX_MLLinOp.H:52
LPInfo & setMaxCoarseningLevel(int n) noexcept
Definition AMReX_MLLinOp.H:57
LPInfo & setAgglomeration(bool x) noexcept
Definition AMReX_MLLinOp.H:49
int semicoarsening_direction
Definition AMReX_MLLinOp.H:46
Definition AMReX_MLLinOp.H:84
StateMode
Definition AMReX_MLLinOp.H:86
BCMode
Definition AMReX_MLLinOp.H:85
Location
Definition AMReX_MLLinOp.H:87
FabArray memory allocation information.
Definition AMReX_FabArray.H:66
Definition AMReX_MLLinOp.H:615
void operator=(const CommContainer &)=delete
~CommContainer()
Definition AMReX_MLLinOp.H:622
MPI_Comm comm
Definition AMReX_MLLinOp.H:616
CommContainer(const CommContainer &)=delete
CommContainer(CommContainer &&)=delete
CommContainer(MPI_Comm m) noexcept
Definition AMReX_MLLinOp.H:617