AMReX_MLLinOp.H

Go to the documentation of this file.

#ifndef AMREX_ML_LINOP_H_
#define AMREX_ML_LINOP_H_
#include <AMReX_Config.H>
 
#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
#include <AMReX_Hypre.H>
#include <AMReX_HypreNodeLap.H>
#endif
 
#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
#include <AMReX_PETSc.H>
#endif
 
#ifdef AMREX_USE_EB
#include <AMReX_EBMultiFabUtil.H>
#include <AMReX_MultiCutFab.H>
#endif
 
#include <AMReX_Any.H>
#include <AMReX_BndryRegister.H>
#include <AMReX_FabDataType.H>
#include <AMReX_MLMGBndry.H>
#include <AMReX_MLAlgMG.H>
#include <AMReX_MultiFab.H>
#include <AMReX_MultiFabUtil.H>
 
#include <algorithm>
#include <iterator>
#include <limits>
#include <string>
 
namespace amrex {

 
enum class BottomSolver : int {
    Default, smoother, bicgstab, cg, bicgcg, cgbicg, hypre, petsc, custom, algmg
};

struct LPInfo
{
    bool do_agglomeration = true;
    bool do_consolidation = true;
    bool do_semicoarsening = false;
    int agg_grid_size = -1;
    int con_grid_size = -1;
    int con_ratio = 2;
    int con_strategy = 3;
    bool has_metric_term = true;
    int max_coarsening_level = 30;
    int max_semicoarsening_level = 0;
    int semicoarsening_direction = -1;
    int hidden_direction = -1;
    bool deterministic = false;   

    LPInfo& setAgglomeration (bool x) noexcept { do_agglomeration = x; return *this; }
    LPInfo& setConsolidation (bool x) noexcept { do_consolidation = x; return *this; }
    LPInfo& setSemicoarsening (bool x) noexcept { do_semicoarsening = x; return *this; }
    LPInfo& setAgglomerationGridSize (int x) noexcept { agg_grid_size = x; return *this; }
    LPInfo& setConsolidationGridSize (int x) noexcept { con_grid_size = x; return *this; }
    LPInfo& setConsolidationRatio (int x) noexcept { con_ratio = x; return *this; }
    LPInfo& setConsolidationStrategy (int x) noexcept { con_strategy = x; return *this; }
    LPInfo& setMetricTerm (bool x) noexcept { has_metric_term = x; return *this; }
    LPInfo& setMaxCoarseningLevel (int n) noexcept { max_coarsening_level = n; return *this; }
    LPInfo& setMaxSemicoarseningLevel (int n) noexcept { max_semicoarsening_level = n; return *this; }
    LPInfo& setSemicoarseningDirection (int n) noexcept { semicoarsening_direction = n; return *this; }
    LPInfo& setHiddenDirection (int n) noexcept { hidden_direction = n; return *this; }
    LPInfo& setDeterministic (bool x) noexcept { deterministic = x; return *this; }

    [[nodiscard]] bool hasHiddenDimension () const noexcept {
        return hidden_direction >=0 && hidden_direction < AMREX_SPACEDIM;
    }
 
    static constexpr int getDefaultAgglomerationGridSize () {
#ifdef AMREX_USE_GPU
        return 32;
#else
        return AMREX_D_PICK(32, 16, 8);
#endif
    }
 
    static constexpr int getDefaultConsolidationGridSize () {
#ifdef AMREX_USE_GPU
        return 32;
#else
        return AMREX_D_PICK(32, 16, 8);
#endif
    }
};
 
struct LinOpEnumType
{
    enum struct BCMode { Homogeneous, Inhomogeneous };
    enum struct StateMode { Solution, Correction };
    enum struct Location { FaceCenter, FaceCentroid, CellCenter, CellCentroid };
};
 
template <typename T> class MLMGT;
template <typename T> class MLCGSolverT;
template <typename T> class MLPoissonT;
template <typename T> class MLABecLaplacianT;
template <typename T> class GMRESMLMGT;


class HypreMLNodeLap;
 
template <typename MF>
class MLLinOpT
{
public:
 
    template <typename T> friend class MLMGT;
    template <typename T> friend class MLCGSolverT;
    template <typename T> friend class MLPoissonT;
    template <typename T> friend class MLABecLaplacianT;
    template <typename T> friend class GMRESMLMGT;
    friend class HypreMLNodeLap;
 
    using MFType = MF;
    using FAB = typename FabDataType<MF>::fab_type;
    using RT  = typename FabDataType<MF>::value_type;
 
    using BCType = LinOpBCType;
    using BCMode    = LinOpEnumType::BCMode;
    using StateMode = LinOpEnumType::StateMode;
    using Location  = LinOpEnumType::Location;
 
    MLLinOpT () = default;
    virtual ~MLLinOpT () = default;
 
    MLLinOpT (const MLLinOpT<MF>&) = delete;
    MLLinOpT (MLLinOpT<MF>&&) = delete;
    MLLinOpT<MF>& operator= (const MLLinOpT<MF>&) = delete;
    MLLinOpT<MF>& operator= (MLLinOpT<MF>&&) = delete;

    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);
 
    [[nodiscard]] virtual std::string name () const { return std::string("Unspecified"); }

    void setDomainBC (const Array<BCType,AMREX_SPACEDIM>& lobc,
                      const Array<BCType,AMREX_SPACEDIM>& hibc) noexcept;

    void setDomainBC (const Vector<Array<BCType,AMREX_SPACEDIM> >& lobc,
                      const Vector<Array<BCType,AMREX_SPACEDIM> >& hibc);

    void setDomainBCLoc (const Array<Real,AMREX_SPACEDIM>& lo_bcloc,
                         const Array<Real,AMREX_SPACEDIM>& hi_bcloc) noexcept;

    [[nodiscard]] bool needsCoarseDataForBC () const noexcept { return m_needs_coarse_data_for_bc; }

    void setCoarseFineBC (const MF* crse, int crse_ratio,
                          LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;
 
    void setCoarseFineBC (const MF* crse, IntVect const& crse_ratio,
                          LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;
 
    template <typename AMF>
    requires (!std::same_as<MF,AMF>)
    void setCoarseFineBC (const AMF* crse, int crse_ratio,
                          LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;
 
    template <typename AMF>
    requires (!std::same_as<MF,AMF>)
    void setCoarseFineBC (const AMF* crse, IntVect const& crse_ratio,
                          LinOpBCType bc_type = LinOpBCType::Dirichlet) noexcept;
 

    virtual void setLevelBC (int /*amrlev*/, const MF* /*levelbcdata*/,
                             const MF* /*robinbc_a*/ = nullptr,
                             const MF* /*robinbc_b*/ = nullptr,
                             const MF* /*robinbc_f*/ = nullptr) = 0;
 
    template <MultiFabLike AMF>
    requires (!std::same_as<MF,AMF>)
    void setLevelBC (int amrlev, const AMF* levelbcdata,
                     const AMF* robinbc_a = nullptr,
                     const AMF* robinbc_b = nullptr,
                     const AMF* robinbc_f = nullptr);

    void setVerbose (int v) noexcept { verbose = v; }

    void setPrintIndentation (std::string s) { print_ident = std::move(s); }

    void setMaxOrder (int o) noexcept { maxorder = o; }
    [[nodiscard]] int getMaxOrder () const noexcept { return maxorder; }

    void setEnforceSingularSolvable (bool o) noexcept { enforceSingularSolvable = o; }
    [[nodiscard]] bool getEnforceSingularSolvable () const noexcept { return enforceSingularSolvable; }
 
    [[nodiscard]] virtual BottomSolver getDefaultBottomSolver () const { return BottomSolver::bicgstab; }

    [[nodiscard]] virtual bool supportCustomBottomSolver () const { return false; }

    virtual void customBottomSolve (MLMGT<MF>* mlmg, MF& x, const MF& b,
                                    RT eps_rel, RT eps_abs, int maxiter)
    {
        amrex::ignore_unused(mlmg, x, b, eps_rel, eps_abs, maxiter);
        amrex::Abort("customBottomSolve not implemented");
    }

    [[nodiscard]] virtual int getNComp () const { return 1; }
 
    [[nodiscard]] virtual int getNGrow (int /*a_lev*/ = 0, int /*mg_lev*/ = 0) const { return 0; }

    [[nodiscard]] virtual bool needsUpdate () const { return false; }
    virtual void update () {}

    virtual void restriction (int amrlev, int cmglev, MF& crse, MF& fine) const = 0;

    virtual void interpolation (int amrlev, int fmglev, MF& fine, const MF& crse) const = 0;

    virtual void interpAssign (int amrlev, int fmglev, MF& fine, MF& crse) const
    {
        amrex::ignore_unused(amrlev, fmglev, fine, crse);
        amrex::Abort("MLLinOpT::interpAssign: Must be implemented for FMG cycle");
    }

    virtual void interpolationAmr (int famrlev, MF& fine, const MF& crse,
                                   IntVect const& nghost) const
    {
        amrex::ignore_unused(famrlev, fine, crse, nghost);
        amrex::Abort("MLLinOpT::interpolationAmr: Must be implemented for composite solves across multiple AMR levels");
    }

    virtual void interpCoarseFineSolution (int famrlev, MF& fine_sol, const MF& crse_sol) const
    {
        amrex::ignore_unused(famrlev, fine_sol, crse_sol);
    }

    virtual void averageDownSolutionRHS (int camrlev, MF& crse_sol, MF& crse_rhs,
                                         const MF& fine_sol, const MF& fine_rhs)
    {
        amrex::ignore_unused(camrlev, crse_sol, crse_rhs, fine_sol, fine_rhs);
        amrex::Abort("MLLinOpT::averageDownSolutionRHS: Must be implemented for composite solves across multiple AMR levels");
    }

    virtual void apply (int amrlev, int mglev, MF& out, MF& in, BCMode bc_mode,
                        StateMode s_mode, const MLMGBndryT<MF>* bndry=nullptr) const = 0;

    virtual void smooth (int amrlev, int mglev, MF& sol, const MF& rhs,
                         bool skip_fillboundary, int niter) const = 0;

    virtual void normalize (int amrlev, int mglev, MF& mf) const {
        amrex::ignore_unused(amrlev, mglev, mf);
    }

    virtual void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b,
                                   const MF* crse_bcdata=nullptr) = 0;

    virtual void prepareForFluxes (int amrlev, const MF* crse_bcdata = nullptr) {
        amrex::ignore_unused(amrlev, crse_bcdata);
    }

    virtual void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b,
                                     BCMode bc_mode, const MF* crse_bcdata=nullptr) = 0;

    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
    {
        amrex::ignore_unused(crse_amrlev, res, crse_sol, crse_rhs, fine_res,
                             fine_sol, fine_rhs);
        amrex::Abort("MLLinOpT::reflux: Must be implemented for composite solves across multiple AMR levels");
    }

    virtual void compFlux (int amrlev, const Array<MF*,AMREX_SPACEDIM>& fluxes,
                           MF& sol, Location loc) const
    {
        amrex::ignore_unused(amrlev, fluxes, sol, loc);
        amrex::Abort("AMReX_MLLinOp::compFlux::How did we get here?");
    }

    virtual void compGrad (int amrlev, const Array<MF*,AMREX_SPACEDIM>& grad,
                           MF& sol, Location loc) const
    {
        amrex::ignore_unused(amrlev, grad, sol, loc);
        amrex::Abort("AMReX_MLLinOp::compGrad::How did we get here?");
    }

    virtual void applyMetricTerm (int amrlev, int mglev, MF& rhs) const {
        amrex::ignore_unused(amrlev, mglev, rhs);
    }

    virtual void unapplyMetricTerm (int amrlev, int mglev, MF& rhs) const {
        amrex::ignore_unused(amrlev, mglev, rhs);
    }

    virtual void unimposeNeumannBC (int amrlev, MF& rhs) const {
        amrex::ignore_unused(amrlev, rhs);
    }

    virtual void applyInhomogNeumannTerm (int amrlev, MF& rhs) const {
        amrex::ignore_unused(amrlev, rhs);
    }

    virtual void applyOverset (int amrlev, MF& rhs) const {
        amrex::ignore_unused(amrlev, rhs);
    }

    [[nodiscard]] virtual bool scaleRHS (int amrlev, MF* rhs) const {
        amrex::ignore_unused(amrlev, rhs);
        return false;
    }

    virtual Vector<RT> getSolvabilityOffset (int amrlev, int mglev,
                                               MF const& rhs) const {
        amrex::ignore_unused(amrlev, mglev, rhs);
        return {};
    }

    virtual void fixSolvabilityByOffset (int amrlev, int mglev, MF& rhs,
                                         Vector<RT> const& offset) const {
        amrex::ignore_unused(amrlev, mglev, rhs, offset);
    }

    virtual void prepareForSolve () = 0;

    virtual void preparePrecond () {}

    virtual void setDirichletNodesToZero (int amrlev, int mglev,
                                          MF& mf) const
    {
        amrex::ignore_unused(amrlev, mglev, mf);
        amrex::Warning("This function might need to be implemented for GMRES to work with this LinOp.");
    }

    [[nodiscard]] virtual bool isSingular (int amrlev) const = 0;
    [[nodiscard]] virtual bool isBottomSingular () const = 0;

    virtual RT xdoty (int amrlev, int mglev, const MF& x, const MF& y, bool local) const = 0;

    virtual RT dotProductPrecond (Vector<MF const*> const& x, Vector<MF const*> const& y) const;

    virtual RT norm2Precond (Vector<MF const*> const& x) const;

    virtual std::unique_ptr<MLLinOpT<MF>> makeNLinOp (int grid_size) const
    {
        amrex::ignore_unused(grid_size);
        amrex::Abort("MLLinOp::makeNLinOp: NSolve not supported");
        return nullptr;
    }

    virtual void getFluxes (const Vector<Array<MF*,AMREX_SPACEDIM> >& a_flux,
                            const Vector<MF*>& a_sol,
                            Location a_loc) const {
        amrex::ignore_unused(a_flux, a_sol, a_loc);
        amrex::Abort("MLLinOp::getFluxes: How did we get here?");
    }

    virtual void getFluxes (const Vector<MF*>& a_flux,
                            const Vector<MF*>& a_sol) const {
        amrex::ignore_unused(a_flux, a_sol);
        amrex::Abort("MLLinOp::getFluxes: How did we get here?");
    }
 
#ifdef AMREX_USE_EB
    virtual void getEBFluxes (const Vector<MF*>& a_flux,
                              const Vector<MF*>& a_sol) const {
        amrex::ignore_unused(a_flux, a_sol);
        amrex::Abort("MLLinOp::getEBFluxes: How did we get here?");
    }
#endif

    [[nodiscard]] virtual std::unique_ptr<MLAlgMG> makeAlgMG (int mglev) const {
        amrex::ignore_unused(mglev);
        amrex::Abort("MLLinOp::makeAlgMG: not supported by this operator");
        return {nullptr};
    }

    [[nodiscard]] virtual bool supportsAlgMG () const { return false; }

    [[nodiscard]] virtual bool supportsAnisotropicCoarsening () const { return false; }

    [[nodiscard]] virtual GpuArray<Real,AMREX_SPACEDIM>
    anisotropicCoarseningCellSize (Geometry const& geom) const { return geom.CellSizeArray(); }
 
#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
    [[nodiscard]] virtual std::unique_ptr<Hypre> makeHypre (Hypre::Interface hypre_interface) const {
        amrex::ignore_unused(hypre_interface);
        amrex::Abort("MLLinOp::makeHypre: How did we get here?");
        return {nullptr};
    }
    [[nodiscard]] virtual std::unique_ptr<HypreNodeLap> makeHypreNodeLap(
        int bottom_verbose,
        const std::string& options_namespace) const
    {
        amrex::ignore_unused(bottom_verbose, options_namespace);
        amrex::Abort("MLLinOp::makeHypreNodeLap: How did we get here?");
        return {nullptr};
    }
#endif
 
#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
    [[nodiscard]] virtual std::unique_ptr<PETScABecLap> makePETSc () const {
        amrex::Abort("MLLinOp::makePETSc: How did we get here?");
        return {nullptr};
    }
#endif

    [[nodiscard]] virtual bool supportNSolve () const { return false; }

    virtual void copyNSolveSolution (MF& dst, MF const& src) const {
        amrex::ignore_unused(dst, src);
    }

    virtual void postSolve (Vector<MF*> const& sol) const {
        amrex::ignore_unused(sol);
    }

    [[nodiscard]] virtual RT normInf (int amrlev, MF const& mf, bool local) const = 0;

    virtual void averageDownAndSync (Vector<MF>& sol) const = 0;
 
    virtual void avgDownResAmr (int clev, MF& cres, MF const& fres) const
    {
        amrex::ignore_unused(clev, cres, fres);
        amrex::Abort("MLLinOpT::avgDownResAmr: Must be implemented for composite solves across multiple AMR levels");
    }

    virtual void avgDownResMG (int clev, MF& cres, MF const& fres) const;

    virtual void beginPrecondBC () { m_precond_mode = true; }
    virtual void endPrecondBC () { m_precond_mode = false; }

    [[nodiscard]] bool isMFIterSafe (int amrlev, int mglev1, int mglev2) const;

    [[nodiscard]] int NAMRLevels () const noexcept { return m_num_amr_levels; }

    [[nodiscard]] int NMGLevels (int amrlev) const noexcept { return m_num_mg_levels[amrlev]; }

    [[nodiscard]] const Geometry& Geom (int amr_lev, int mglev=0) const noexcept { return m_geom[amr_lev][mglev]; }
 
    // BC
    Vector<Array<BCType, AMREX_SPACEDIM> > m_lobc;
    Vector<Array<BCType, AMREX_SPACEDIM> > m_hibc;
    // Need to save the original copy because we change the BC type to
    // Neumann for inhomogeneous Neumann and Robin.
    Vector<Array<BCType, AMREX_SPACEDIM> > m_lobc_orig;
    Vector<Array<BCType, AMREX_SPACEDIM> > m_hibc_orig;
 
protected:
 
    static constexpr int mg_coarsen_ratio = 2;
    int mg_box_min_width = 2;
    int mg_domain_min_width = 2;
    int mg_agg_no_split_direction = -1; // keep agglomerated grids whole in it
    bool mg_odd_coarsening = false; // coarsen by 3 or 5 where 2 does not fit
    bool mg_independent_coarsening = false; // coarsen directions on their own
 
    LPInfo info;
 
    int verbose = 0;
    std::string print_ident;
 
    int maxorder = 3;
 
    bool enforceSingularSolvable = true;
    bool m_mg_deferred = false; 
    bool m_mg_built = false;
 
    int m_num_amr_levels = 0;
    Vector<int> m_amr_ref_ratio;
 
    Vector<int> m_num_mg_levels;
    const MLLinOpT<MF>* m_parent = nullptr;
 
    IntVect m_ixtype;
 
    bool m_do_agglomeration = false;
    bool m_do_consolidation = false;
 
    bool m_do_semicoarsening = false;
    Vector<IntVect> mg_coarsen_ratio_vec;

    Vector<Vector<Geometry> >            m_geom;
    Vector<Vector<BoxArray> >            m_grids;
    Vector<Vector<DistributionMapping> > m_dmap;
    Vector<Vector<std::unique_ptr<FabFactory<FAB> > > > m_factory;
    Vector<int>                          m_domain_covered;
 
    MPI_Comm m_default_comm = MPI_COMM_NULL;
    MPI_Comm m_bottom_comm = MPI_COMM_NULL;
    struct CommContainer {
        MPI_Comm comm;
        CommContainer (MPI_Comm m) noexcept : comm(m) {}
        CommContainer (const CommContainer&) = delete;
        CommContainer (CommContainer&&) = delete;
        void operator= (const CommContainer&) = delete;
        void operator= (CommContainer&&) = delete;
        ~CommContainer () { // NOLINT(modernize-use-equals-default)
#ifdef BL_USE_MPI
            if (comm != MPI_COMM_NULL) { MPI_Comm_free(&comm); }
#endif
        }
    };
    std::unique_ptr<CommContainer> m_raii_comm;
 
    Array<Real, AMREX_SPACEDIM> m_domain_bloc_lo {{AMREX_D_DECL(0._rt,0._rt,0._rt)}};
    Array<Real, AMREX_SPACEDIM> m_domain_bloc_hi {{AMREX_D_DECL(0._rt,0._rt,0._rt)}};
 
    bool m_needs_coarse_data_for_bc = false;
    LinOpBCType m_coarse_fine_bc_type = LinOpBCType::Dirichlet;
    IntVect m_coarse_data_crse_ratio = IntVect(-1);
    RealVect m_coarse_bc_loc;
    const MF* m_coarse_data_for_bc = nullptr;
    MF m_coarse_data_for_bc_raii;
 
    bool m_precond_mode = false;

    [[nodiscard]] const Vector<int>& AMRRefRatio () const noexcept { return m_amr_ref_ratio; }

    [[nodiscard]] int AMRRefRatio (int amr_lev) const noexcept { return m_amr_ref_ratio[amr_lev]; }

    [[nodiscard]] IntVect AMRRefRatioVect (int amr_lev) const noexcept {
        IntVect rr(m_amr_ref_ratio[amr_lev]);
        if (info.hasHiddenDimension()) { rr[info.hidden_direction] = 1; }
        return rr;
    }
 
    [[nodiscard]] FabFactory<FAB> const* Factory (int amr_lev, int mglev=0) const noexcept {
        return m_factory[amr_lev][mglev].get();
    }
 
    [[nodiscard]] GpuArray<BCType,AMREX_SPACEDIM> LoBC (int icomp = 0) const noexcept {
        return GpuArray<BCType,AMREX_SPACEDIM>{{AMREX_D_DECL(m_lobc[icomp][0],
                                                             m_lobc[icomp][1],
                                                             m_lobc[icomp][2])}};
    }
    [[nodiscard]] GpuArray<BCType,AMREX_SPACEDIM> HiBC (int icomp = 0) const noexcept {
        return GpuArray<BCType,AMREX_SPACEDIM>{{AMREX_D_DECL(m_hibc[icomp][0],
                                                             m_hibc[icomp][1],
                                                             m_hibc[icomp][2])}};
    }
 
    [[nodiscard]] bool hasBC (BCType bct) const noexcept;
    [[nodiscard]] bool hasInhomogNeumannBC () const noexcept;
    [[nodiscard]] bool hasRobinBC () const noexcept;
 
    [[nodiscard]] virtual bool supportRobinBC () const noexcept { return false; }
    [[nodiscard]] virtual bool supportInhomogNeumannBC () const noexcept { return false; }
 
#ifdef BL_USE_MPI
    [[nodiscard]] bool isBottomActive () const noexcept { return m_bottom_comm != MPI_COMM_NULL; }
#else
    [[nodiscard]] bool isBottomActive () const noexcept { return true; }
#endif
    [[nodiscard]] MPI_Comm BottomCommunicator () const noexcept { return m_bottom_comm; }
    [[nodiscard]] MPI_Comm Communicator () const noexcept { return m_default_comm; }
 
    void setCoarseFineBCLocation (const RealVect& cloc) noexcept { m_coarse_bc_loc = cloc; }
 
    [[nodiscard]] bool doAgglomeration () const noexcept { return m_do_agglomeration; }
    [[nodiscard]] bool doConsolidation () const noexcept { return m_do_consolidation; }
    [[nodiscard]] bool doSemicoarsening () const noexcept { return m_do_semicoarsening; }
 
    [[nodiscard]] bool isCellCentered () const noexcept { return m_ixtype == 0; }
 
    [[nodiscard]] virtual IntVect getNGrowVectRestriction () const {
        return isCellCentered() ? IntVect(0) : IntVect(1);
    }
 
    virtual void make (Vector<Vector<MF> >& mf, IntVect const& ng) const;
 
    [[nodiscard]] virtual MF make (int amrlev, int mglev, IntVect const& ng) const;

    [[nodiscard]] virtual MF make (int amrlev, int mglev, IntVect const& ng,
                                   MFInfo const& mf_info) const;
 
    [[nodiscard]] virtual MF makeAlias (MF const& mf) const;

    [[nodiscard]] virtual MF makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const;

    [[nodiscard]] virtual MF makeCoarseMG (int amrlev, int mglev, IntVect const& ng,
                                           MFInfo const& mf_info) const;

    [[nodiscard]] virtual MF makeCoarseAmr (int famrlev, IntVect const& ng) const;

    [[nodiscard]] virtual MF makeCoarseAmr (int famrlev, IntVect const& ng,
                                            MFInfo const& mf_info) const;
 
    [[nodiscard]] virtual std::unique_ptr<FabFactory<FAB> > makeFactory (int /*amrlev*/, int /*mglev*/) const {
        return std::make_unique<DefaultFabFactory<FAB>>();
    }
 
    virtual void resizeMultiGrid (int new_size);

    void prepareMGHierarchy ();

    virtual void buildMGHierarchy () {
        BL_PROFILE("MLLinOp::buildMGHierarchy()");
        defineCoarseMGLevels();
    }
 
    [[nodiscard]] bool hasHiddenDimension () const noexcept { return info.hasHiddenDimension(); }
    [[nodiscard]] int hiddenDirection () const noexcept { return info.hidden_direction; }
    [[nodiscard]] Box compactify (Box const& b) const noexcept;
 
    template <typename T>
    [[nodiscard]] Array4<T> compactify (Array4<T> const& a) const noexcept
    {
        if (info.hidden_direction == 0) {
            return Array4<T>(a.dataPtr(), {a.begin[1],a.begin[2],0}, {a.end[1],a.end[2],1}, a.nComp());
        } else if (info.hidden_direction == 1) {
            return Array4<T>(a.dataPtr(), {a.begin[0],a.begin[2],0}, {a.end[0],a.end[2],1}, a.nComp());
        } else if (info.hidden_direction == 2) {
            return Array4<T>(a.dataPtr(), {a.begin[0],a.begin[1],0}, {a.end[0],a.end[1],1}, a.nComp());
        } else {
            return a;
        }
    }
 
    template <typename T>
    [[nodiscard]] T get_d0 (T const& dx, T const& dy, T const&) const noexcept
    {
        if (info.hidden_direction == 0) {
            return dy;
        } else {
            return dx;
        }
    }
 
    template <typename T>
    [[nodiscard]] T get_d1 (T const&, T const& dy, T const& dz) const noexcept
    {
        if (info.hidden_direction == 0 || info.hidden_direction == 1) {
            return dz;
        } else {
            return dy;
        }
    }
 
private:
 
    void defineGrids (const Vector<Geometry>& a_geom,
                      const Vector<BoxArray>& a_grids,
                      const Vector<DistributionMapping>& a_dmap,
                      const Vector<FabFactory<FAB> const*>& a_factory);
    void defineBC ();
    void defineCoarseMGLevels ();
    static void makeAgglomeratedDMap (const Vector<BoxArray>& ba, Vector<DistributionMapping>& dm);
    static void makeConsolidatedDMap (const Vector<BoxArray>& ba, Vector<DistributionMapping>& dm,
                                      int ratio, int strategy);
    [[nodiscard]] MPI_Comm makeSubCommunicator (const DistributionMapping& dm);
 
    virtual void checkPoint (std::string const& /*file_name*/) const {
        amrex::Abort("MLLinOp:checkPoint: not implemented");
    }
 
    Vector<std::unique_ptr<MF>> levelbc_raii;
    Vector<std::unique_ptr<MF>> robin_a_raii;
    Vector<std::unique_ptr<MF>> robin_b_raii;
    Vector<std::unique_ptr<MF>> robin_f_raii;
};
 
template <typename MF>
void
MLLinOpT<MF>::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,
                      [[maybe_unused]] bool eb_limit_coarsening)
{
    BL_PROFILE("MLLinOp::define()");
 
    info = a_info;
#ifdef AMREX_USE_GPU
    if (Gpu::notInLaunchRegion())
    {
        if (info.agg_grid_size <= 0) { info.agg_grid_size = AMREX_D_PICK(32, 16, 8); }
        if (info.con_grid_size <= 0) { info.con_grid_size = AMREX_D_PICK(32, 16, 8); }
    }
    else
#endif
    {
        if (info.agg_grid_size <= 0) { info.agg_grid_size = LPInfo::getDefaultAgglomerationGridSize(); }
        if (info.con_grid_size <= 0) { info.con_grid_size = LPInfo::getDefaultConsolidationGridSize(); }
    }
 
#ifdef AMREX_USE_EB
    if (!a_factory.empty() && eb_limit_coarsening) {
        const auto *f = dynamic_cast<EBFArrayBoxFactory const*>(a_factory[0]);
        if (f) {
            info.max_coarsening_level = std::min(info.max_coarsening_level,
                                                 f->maxCoarseningLevel());
        }
    }
#endif
    m_mg_deferred = supportsAnisotropicCoarsening();
    m_mg_built = false;
    defineGrids(a_geom, a_grids, a_dmap, a_factory);
    defineBC();
}
 
template <typename MF>
void
MLLinOpT<MF>::defineGrids (const Vector<Geometry>& a_geom,
                           const Vector<BoxArray>& a_grids,
                           const Vector<DistributionMapping>& a_dmap,
                           const Vector<FabFactory<FAB> const*>& a_factory)
{
    BL_PROFILE("MLLinOp::defineGrids()");
 
#ifdef AMREX_USE_EB
    if ( ! a_factory.empty() ) {
        auto const* ebf = dynamic_cast<EBFArrayBoxFactory const*>(a_factory[0]);
        if (ebf && !(ebf->isAllRegular())) { // Has non-trivial EB
            mg_domain_min_width = 4;
        }
    }
#endif
 
    m_num_amr_levels = 0;
    for (int amrlev = 0; amrlev < std::ssize(a_geom); amrlev++) {
        if (!a_grids[amrlev].empty()) {
            m_num_amr_levels++;
        }
    }
 
    m_amr_ref_ratio.resize(m_num_amr_levels);
    m_num_mg_levels.resize(m_num_amr_levels);
 
    m_geom.clear();
    m_grids.clear();
    m_dmap.clear();
    m_factory.clear();
    m_domain_covered.clear();
    mg_coarsen_ratio_vec.clear();
 
    m_geom.resize(m_num_amr_levels);
    m_grids.resize(m_num_amr_levels);
    m_dmap.resize(m_num_amr_levels);
    m_factory.resize(m_num_amr_levels);
 
    m_default_comm = ParallelContext::CommunicatorSub();
 
    const RealBox& rb = a_geom[0].ProbDomain();
    const int coord = a_geom[0].Coord();
    const Array<int,AMREX_SPACEDIM>& is_per = a_geom[0].isPeriodic();
 
    IntVect mg_coarsen_ratio_v(mg_coarsen_ratio);
    if (hasHiddenDimension()) {
        AMREX_ASSERT_WITH_MESSAGE(AMREX_SPACEDIM == 3,
                                  "Hidden direction only supported for 3d");
        mg_coarsen_ratio_v[info.hidden_direction] = 1;
    }
 
    // fine amr levels
    for (int amrlev = m_num_amr_levels-1; amrlev > 0; --amrlev)
    {
        m_num_mg_levels[amrlev] = 1;
        m_geom[amrlev].push_back(a_geom[amrlev]);
        m_grids[amrlev].push_back(a_grids[amrlev]);
        m_dmap[amrlev].push_back(a_dmap[amrlev]);
        if (amrlev < std::ssize(a_factory)) {
            m_factory[amrlev].emplace_back(a_factory[amrlev]->clone());
        } else {
            m_factory[amrlev].push_back(std::make_unique<DefaultFabFactory<FAB>>());
        }
 
        IntVect rr = mg_coarsen_ratio_v;
        const Box& dom = a_geom[amrlev].Domain();
        for (int i = 0; i < 2; ++i)
        {
            if (!dom.coarsenable(rr)) { amrex::Abort("MLLinOp: Uncoarsenable domain"); }
 
            const Box& cdom = amrex::coarsen(dom,rr);
            if (cdom == a_geom[amrlev-1].Domain()) { break; }
 
            ++(m_num_mg_levels[amrlev]);
 
            m_geom[amrlev].emplace_back(cdom, rb, coord, is_per);
 
            m_grids[amrlev].push_back(a_grids[amrlev]);
            AMREX_ASSERT(m_grids[amrlev].back().coarsenable(rr));
            m_grids[amrlev].back().coarsen(rr);
 
            m_dmap[amrlev].push_back(a_dmap[amrlev]);
 
            rr *= mg_coarsen_ratio_v;
        }
 
#if (AMREX_SPACEDIM > 1)
        if (hasHiddenDimension()) {
            m_amr_ref_ratio[amrlev-1] = rr[(info.hidden_direction+1) % AMREX_SPACEDIM];
        } else
#endif
        {
            m_amr_ref_ratio[amrlev-1] = rr[0];
        }
    }
 
    // coarsest amr level
    m_num_mg_levels[0] = 1;
    m_geom[0].push_back(a_geom[0]);
    m_grids[0].push_back(a_grids[0]);
    m_dmap[0].push_back(a_dmap[0]);
    if (!a_factory.empty()) {
        m_factory[0].emplace_back(a_factory[0]->clone());
    } else {
        m_factory[0].push_back(std::make_unique<DefaultFabFactory<FAB>>());
    }
 
    m_domain_covered.resize(m_num_amr_levels, false);
    auto npts0 = m_grids[0][0].numPts();
    m_domain_covered[0] = (npts0 == compactify(m_geom[0][0].Domain()).numPts());
    for (int amrlev = 1; amrlev < m_num_amr_levels; ++amrlev)
    {
        if (!m_domain_covered[amrlev-1]) { break; }
        m_domain_covered[amrlev] = (m_grids[amrlev][0].numPts() ==
                                    compactify(m_geom[amrlev][0].Domain()).numPts());
    }
 
    if (m_mg_deferred) {
        // MG levels > 0 of AMR level 0 are built by prepareMGHierarchy.
        m_bottom_comm = m_default_comm;
        m_do_agglomeration = false;
        m_do_consolidation = false;
    } else {
        defineCoarseMGLevels();
    }
 
    for (int amrlev = 1; amrlev < m_num_amr_levels; ++amrlev)
    {
        for (int mglev = 1; mglev < m_num_mg_levels[amrlev]; ++mglev)
        {
            m_factory[amrlev].emplace_back(makeFactory(amrlev,mglev));
        }
    }
 
    for (int amrlev = 1; amrlev < m_num_amr_levels; ++amrlev)
    {
        AMREX_ASSERT_WITH_MESSAGE(m_grids[amrlev][0].coarsenable(AMRRefRatioVect(amrlev-1)),
                                  "MLLinOp: grids not coarsenable between AMR levels");
    }
}
 
template <typename MF>
void
MLLinOpT<MF>::defineCoarseMGLevels ()
{
    BL_PROFILE("MLLinOp::defineCoarseMGLevels()");
 
    AMREX_ASSERT(m_grids[0].size() == 1 && mg_coarsen_ratio_vec.empty());
 
    // Copies, because m_geom[0], m_grids[0] and m_dmap[0] grow below.
    Geometry const geom0 = m_geom[0][0];
    BoxArray const ba0 = m_grids[0][0];
    DistributionMapping const dm0 = m_dmap[0][0];
 
    const RealBox& rb = geom0.ProbDomain();
    const int coord = geom0.Coord();
    const Array<int,AMREX_SPACEDIM>& is_per = geom0.isPeriodic();
    auto const npts0 = ba0.numPts();
 
    IntVect mg_coarsen_ratio_v(mg_coarsen_ratio);
    IntVect mg_box_min_width_v(mg_box_min_width);
    IntVect mg_domain_min_width_v(mg_domain_min_width);
    if (hasHiddenDimension()) {
        mg_coarsen_ratio_v[info.hidden_direction] = 1;
        mg_box_min_width_v[info.hidden_direction] = 0;
        mg_domain_min_width_v[info.hidden_direction] = 0;
    }
 
    Box aggbox;
    bool aggable = false;
 
    if (ba0.size() > 1 && info.do_agglomeration)
    {
        if (m_domain_covered[0])
        {
            aggbox = geom0.Domain();
            if (hasHiddenDimension()) {
                aggbox.makeSlab(hiddenDirection(), ba0[0].smallEnd(hiddenDirection()));
            }
            aggable = true;
        }
        else
        {
            aggbox = ba0.minimalBox();
            aggable = (aggbox.numPts() == npts0);
        }
    }
 
    bool agged = false;
    bool coned = false;
    int agg_lev = 0, con_lev = 0;
 
    AMREX_ALWAYS_ASSERT( ! (info.do_semicoarsening && info.hasHiddenDimension()
                            && (info.semicoarsening_direction == -1 ||
                                info.semicoarsening_direction == info.hidden_direction))
                         && info.semicoarsening_direction >= -1
                         && info.semicoarsening_direction < AMREX_SPACEDIM );
 
    // Directions to coarsen for stretched cells: those shorter than 1.5
    // times the shortest. A fixed semicoarsening direction takes precedence.
    bool const aniso_coarsening = supportsAnisotropicCoarsening() && !hasHiddenDimension()
        && !(info.do_semicoarsening && info.semicoarsening_direction != -1);
    auto const dx0 = anisotropicCoarseningCellSize(geom0);
    auto base_ratio = [&] (IntVect const& accum) -> IntVect
    {
        IntVect rr = mg_coarsen_ratio_v;
        if (aniso_coarsening) {
            Real dxmin = std::numeric_limits<Real>::max();
            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                dxmin = std::min(dxmin, dx0[idim]*Real(accum[idim]));
            }
            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                if (dx0[idim]*Real(accum[idim]) >= Real(1.5)*dxmin) { rr[idim] = 1; }
            }
        }
        return rr;
    };
 
    if (info.do_agglomeration && aggable)
    {
        Box dbx = geom0.Domain();
        Box bbx = aggbox;
        Real const nbxs = static_cast<Real>(ba0.size());
        Long const threshold_npts = AMREX_D_TERM(Long(info.agg_grid_size),
                                                     *info.agg_grid_size,
                                                     *info.agg_grid_size);
        Vector<Box> domainboxes{dbx};
        Vector<Box> boundboxes{bbx};
        Vector<int> agg_flag{false};
        Vector<IntVect> accum_coarsen_ratio{IntVect(1)};
        int numsclevs = 0;
 
        for (int lev = 0; lev < info.max_coarsening_level; ++lev)
        {
            IntVect const rr_base = base_ratio(accum_coarsen_ratio.back());
            IntVect rr_level = rr_base;
            bool const do_semicoarsening_level = info.do_semicoarsening
                && numsclevs < info.max_semicoarsening_level;
            if (do_semicoarsening_level
                && info.semicoarsening_direction != -1)
            {
                rr_level[info.semicoarsening_direction] = 1;
            }
            IntVect is_coarsenable;
            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                IntVect rr_dir(1);
                rr_dir[idim] = rr_level[idim];
                // With odd ratios, the first ratio must also fit the original boxes.
                bool const fit_ba0 = mg_odd_coarsening && lev == 0;
                is_coarsenable[idim] = dbx.coarsenable(rr_dir, mg_domain_min_width_v)
                    && bbx.coarsenable(rr_dir, mg_box_min_width_v)
                    && (!fit_ba0 || ba0.coarsenable(rr_dir, mg_box_min_width_v));
                for (int r : {3, 5}) {
                    if (is_coarsenable[idim] || !mg_odd_coarsening || rr_level[idim] != 2) {
                        break;
                    }
                    rr_dir[idim] = r;
                    if (dbx.coarsenable(rr_dir, mg_domain_min_width_v)
                        && bbx.coarsenable(rr_dir, mg_box_min_width_v)
                        && (!fit_ba0 || ba0.coarsenable(rr_dir, mg_box_min_width_v)))
                    {
                        is_coarsenable[idim] = true;
                        rr_level[idim] = r;
                    }
                }
                if (!is_coarsenable[idim] && do_semicoarsening_level
                    && info.semicoarsening_direction == -1)
                {
                    is_coarsenable[idim] = true;
                    rr_level[idim] = 1;
                }
            }
            if (mg_independent_coarsening) {
                // A direction that cannot be coarsened gets ratio 1, and one whose
                // cells are 1.5 times longer than the shortest coarsenable waits.
                Real hmin = std::numeric_limits<Real>::max();
                for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                    if (is_coarsenable[idim] && rr_level[idim] > 1) {
                        hmin = std::min(hmin, dx0[idim]*Real(accum_coarsen_ratio.back()[idim]));
                    }
                }
                for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                    if (!is_coarsenable[idim] || dx0[idim]*Real(accum_coarsen_ratio.back()[idim]) >= Real(1.5)*hmin) {
                        is_coarsenable[idim] = true;
                        rr_level[idim] = 1;
                    }
                }
            }
            if (is_coarsenable != IntVect(1) || rr_level == IntVect(1)) {
                break;
            }
            if (do_semicoarsening_level && info.semicoarsening_direction == -1) {
                // make sure there is at most one direction that is not coarsened
                int n_ones = AMREX_D_TERM(  static_cast<int>(rr_level[0] == 1),
                                          + static_cast<int>(rr_level[1] == 1),
                                          + static_cast<int>(rr_level[2] == 1));
                if (n_ones > 1 && !aniso_coarsening) { break; }
            }
            bool semi_level = false; // odd ratios do not count
            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                semi_level = semi_level || (rr_level[idim] == 1 && rr_base[idim] != 1);
            }
            if (semi_level) {
                ++numsclevs;
            }
 
            accum_coarsen_ratio.push_back(accum_coarsen_ratio.back()*rr_level);
            domainboxes.push_back(dbx.coarsen(rr_level));
            boundboxes.push_back(bbx.coarsen(rr_level));
            bool to_agg = (bbx.d_numPts() / nbxs)
                < Real(0.999)*static_cast<Real>(threshold_npts);
            agg_flag.push_back(to_agg);
        }
 
        for (int lev = 1, nlevs = static_cast<int>(domainboxes.size()); lev < nlevs; ++lev) {
            if (!agged && !agg_flag[lev] &&
                ba0.coarsenable(accum_coarsen_ratio[lev], mg_box_min_width_v))
            {
                m_grids[0].push_back(amrex::coarsen(ba0, accum_coarsen_ratio[lev]));
                m_dmap[0].push_back(dm0);
            } else {
                IntVect cr = domainboxes[lev-1].length() / domainboxes[lev].length();
                if (!agged && lev > 1 && !m_grids[0].back().coarsenable(cr)) {
                    // The boxes of level lev-1 cannot be coarsened further,
                    // so agglomerate from that level instead.
                    m_grids[0].pop_back();
                    m_dmap[0].pop_back();
                    m_geom[0].pop_back();
                    --lev;
                    cr = domainboxes[lev-1].length() / domainboxes[lev].length();
                }
                if (!m_grids[0].back().coarsenable(cr)) {
                    break; // average_down would fail if fine boxarray is not coarsenable.
                }
                m_grids[0].emplace_back(boundboxes[lev]);
                Box const cell_box = amrex::enclosedCells(boundboxes[lev]);
                if (cell_box.numPts() > threshold_npts) {
                    IntVect max_grid_size(info.agg_grid_size);
                    if (info.do_semicoarsening && info.max_semicoarsening_level >= lev
                        && info.semicoarsening_direction != -1)
                    {
                        IntVect blen = cell_box.size();
                        AMREX_D_TERM(int mgs_0 = (max_grid_size[0]+blen[0]-1) / blen[0];,
                                     int mgs_1 = (max_grid_size[1]+blen[1]-1) / blen[1];,
                                     int mgs_2 = (max_grid_size[2]+blen[2]-1) / blen[2]);
                        max_grid_size[info.semicoarsening_direction]
                            *= AMREX_D_TERM(mgs_0, *mgs_1, *mgs_2);
                    }
                    if (mg_agg_no_split_direction >= 0 &&
                        mg_agg_no_split_direction < AMREX_SPACEDIM)
                    {
                        int const d = mg_agg_no_split_direction;
                        max_grid_size[d] = std::max(max_grid_size[d], cell_box.length(d));
                    }
                    // Chop in units of the next coarsening ratio, so that
                    // the boxes can be coarsened.
                    IntVect const rn = (lev+1 < nlevs)
                        ? domainboxes[lev].length() / domainboxes[lev+1].length() : IntVect(1);
                    BoxArray& cba = m_grids[0].back();
                    cba.coarsen(rn);
                    cba.maxSize((max_grid_size+rn-1)/rn);
                    cba.refine(rn);
                }
                m_dmap[0].push_back(DistributionMapping());
                if (!agged) {
                    agged = true;
                    agg_lev = lev;
                }
            }
            m_geom[0].emplace_back(domainboxes[lev],rb,coord,is_per);
        }
    }
    else
    {
        Long consolidation_threshold = 0;
        Real avg_npts = 0.0;
        if (info.do_consolidation) {
            avg_npts = static_cast<Real>(ba0.d_numPts()) / static_cast<Real>(ParallelContext::NProcsSub());
            consolidation_threshold = AMREX_D_TERM(Long(info.con_grid_size),
                                                       *info.con_grid_size,
                                                       *info.con_grid_size);
        }
 
        Box const& dom0 = geom0.Domain();
        IntVect rr_vec(1);
        int numsclevs = 0;
        for (int lev = 0; lev < info.max_coarsening_level; ++lev)
        {
            IntVect const rr_base = base_ratio(rr_vec);
            IntVect rr_level = rr_base;
            bool do_semicoarsening_level = info.do_semicoarsening
                && numsclevs < info.max_semicoarsening_level;
            if (do_semicoarsening_level
                && info.semicoarsening_direction != -1)
            {
                rr_level[info.semicoarsening_direction] = 1;
            }
            IntVect is_coarsenable;
            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                IntVect rr_dir(1);
                rr_dir[idim] = rr_vec[idim] * rr_level[idim];
                is_coarsenable[idim] = dom0.coarsenable(rr_dir, mg_domain_min_width_v)
                    && ba0.coarsenable(rr_dir, mg_box_min_width_v);
                for (int r : {3, 5}) {
                    if (is_coarsenable[idim] || !mg_odd_coarsening || rr_level[idim] != 2) {
                        break;
                    }
                    rr_dir[idim] = rr_vec[idim] * r;
                    if (dom0.coarsenable(rr_dir, mg_domain_min_width_v)
                        && ba0.coarsenable(rr_dir, mg_box_min_width_v))
                    {
                        is_coarsenable[idim] = true;
                        rr_level[idim] = r;
                    }
                }
                if (!is_coarsenable[idim] && do_semicoarsening_level
                    && info.semicoarsening_direction == -1)
                {
                    is_coarsenable[idim] = true;
                    rr_level[idim] = 1;
                }
            }
            if (mg_independent_coarsening) {
                // A direction that cannot be coarsened gets ratio 1, and one whose
                // cells are 1.5 times longer than the shortest coarsenable waits.
                Real hmin = std::numeric_limits<Real>::max();
                for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                    if (is_coarsenable[idim] && rr_level[idim] > 1) {
                        hmin = std::min(hmin, dx0[idim]*Real(rr_vec[idim]));
                    }
                }
                for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                    if (!is_coarsenable[idim] || dx0[idim]*Real(rr_vec[idim]) >= Real(1.5)*hmin) {
                        is_coarsenable[idim] = true;
                        rr_level[idim] = 1;
                    }
                }
            }
            if (is_coarsenable != IntVect(1) || rr_level == IntVect(1)) {
                break;
            }
            if (do_semicoarsening_level && info.semicoarsening_direction == -1) {
                // make sure there is at most one direction that is not coarsened
                int n_ones = AMREX_D_TERM(  static_cast<int>(rr_level[0] == 1),
                                          + static_cast<int>(rr_level[1] == 1),
                                          + static_cast<int>(rr_level[2] == 1));
                if (n_ones > 1 && !aniso_coarsening) { break; }
            }
            bool semi_level = false; // odd ratios do not count
            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                semi_level = semi_level || (rr_level[idim] == 1 && rr_base[idim] != 1);
            }
            if (semi_level) {
                ++numsclevs;
            }
            rr_vec *= rr_level;
 
            m_geom[0].emplace_back(amrex::coarsen(dom0, rr_vec), rb, coord, is_per);
            m_grids[0].push_back(amrex::coarsen(ba0, rr_vec));
 
            if (info.do_consolidation)
            {
                if (avg_npts/static_cast<Real>(AMREX_D_TERM(rr_vec[0], *rr_vec[1], *rr_vec[2]))
                    < Real(0.999)*static_cast<Real>(consolidation_threshold))
                {
                    coned = true;
                    con_lev = m_dmap[0].size();
                    m_dmap[0].push_back(DistributionMapping());
                }
                else
                {
                    m_dmap[0].push_back(m_dmap[0].back());
                }
            }
            else
            {
                m_dmap[0].push_back(dm0);
            }
        }
    }
 
    m_num_mg_levels[0] = m_grids[0].size();
 
    for (int mglev = 0; mglev < m_num_mg_levels[0] - 1; mglev++){
        const Box& fine_domain = m_geom[0][mglev].Domain();
        const Box& crse_domain = m_geom[0][mglev+1].Domain();
        mg_coarsen_ratio_vec.push_back(fine_domain.length()/crse_domain.length());
    }
 
    for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) {
        if (AMRRefRatio(amrlev) == 4 && mg_coarsen_ratio_vec.empty()) {
            mg_coarsen_ratio_vec.push_back(IntVect(2));
        }
    }
 
    if (agged)
    {
        makeAgglomeratedDMap(m_grids[0], m_dmap[0]);
    }
    else if (coned)
    {
        makeConsolidatedDMap(m_grids[0], m_dmap[0], info.con_ratio, info.con_strategy);
    }
 
    if (agged || coned)
    {
        m_bottom_comm = makeSubCommunicator(m_dmap[0].back());
    }
    else
    {
        m_bottom_comm = m_default_comm;
    }
 
    m_do_agglomeration = agged;
    m_do_consolidation = coned;
 
    if (verbose > 1) {
        if (agged) {
            Print() << "MLLinOp::defineGrids(): agglomerated AMR level 0 starting at MG level "
                    << agg_lev << " of " << m_num_mg_levels[0] << "\n";
        } else if (coned) {
            Print() << "MLLinOp::defineGrids(): consolidated AMR level 0 starting at MG level "
                    << con_lev << " of " << m_num_mg_levels[0]
                    << " (ratio = " << info.con_ratio << ")" << "\n";
        } else {
            Print() << "MLLinOp::defineGrids(): no agglomeration or consolidation of AMR level 0\n";
        }
    }
 
    for (int mglev = 1; mglev < m_num_mg_levels[0]; ++mglev)
    {
        m_factory[0].emplace_back(makeFactory(0,mglev));
    }
}
 
template <typename MF>
void
MLLinOpT<MF>::prepareMGHierarchy ()
{
    if (m_mg_deferred && !m_mg_built) {
        ParallelContext::push(m_default_comm);
        buildMGHierarchy();
        ParallelContext::pop();
        m_mg_built = true;
    }
}
 
template <typename MF>
void
MLLinOpT<MF>::defineBC ()
{
    m_needs_coarse_data_for_bc = !m_domain_covered[0];
 
    levelbc_raii.resize(m_num_amr_levels);
    robin_a_raii.resize(m_num_amr_levels);
    robin_b_raii.resize(m_num_amr_levels);
    robin_f_raii.resize(m_num_amr_levels);
}
 
template <typename MF>
void
MLLinOpT<MF>::setDomainBC (const Array<BCType,AMREX_SPACEDIM>& a_lobc,
                           const Array<BCType,AMREX_SPACEDIM>& a_hibc) noexcept
{
    const int ncomp = getNComp();
    setDomainBC(Vector<Array<BCType,AMREX_SPACEDIM> >(ncomp,a_lobc),
                Vector<Array<BCType,AMREX_SPACEDIM> >(ncomp,a_hibc));
}
 
template <typename MF>
void
MLLinOpT<MF>::setDomainBC (const Vector<Array<BCType,AMREX_SPACEDIM> >& a_lobc,
                           const Vector<Array<BCType,AMREX_SPACEDIM> >& a_hibc)
{
    const int ncomp = getNComp();
    AMREX_ASSERT_WITH_MESSAGE(ncomp == a_lobc.size() && ncomp == a_hibc.size(),
                              "MLLinOp::setDomainBC: wrong size");
    m_lobc = a_lobc;
    m_hibc = a_hibc;
    m_lobc_orig = m_lobc;
    m_hibc_orig = m_hibc;
    for (int icomp = 0; icomp < ncomp; ++icomp) {
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
            if (m_geom[0][0].isPeriodic(idim)) {
                AMREX_ALWAYS_ASSERT(m_lobc[icomp][idim] == BCType::Periodic &&
                                    m_hibc[icomp][idim] == BCType::Periodic);
            } else {
                AMREX_ALWAYS_ASSERT(m_lobc[icomp][idim] != BCType::Periodic &&
                                    m_hibc[icomp][idim] != BCType::Periodic);
            }
 
            if (m_lobc[icomp][idim] == LinOpBCType::inhomogNeumann ||
                m_lobc[icomp][idim] == LinOpBCType::Robin)
            {
                m_lobc[icomp][idim] = LinOpBCType::Neumann;
            }
 
            if (m_hibc[icomp][idim] == LinOpBCType::inhomogNeumann ||
                m_hibc[icomp][idim] == LinOpBCType::Robin)
            {
                m_hibc[icomp][idim] = LinOpBCType::Neumann;
            }
        }
    }
 
    if (hasHiddenDimension()) {
        const int hd = hiddenDirection();
        for (int n = 0; n < ncomp; ++n) {
            m_lobc[n][hd] = LinOpBCType::Neumann;
            m_hibc[n][hd] = LinOpBCType::Neumann;
        }
    }
 
    if (hasInhomogNeumannBC() && !supportInhomogNeumannBC()) {
        amrex::Abort("Inhomogeneous Neumann BC not supported");
    }
    if (hasRobinBC() && !supportRobinBC()) {
        amrex::Abort("Robin BC not supported");
    }
}
 
template <typename MF>
bool
MLLinOpT<MF>::hasBC (BCType bct) const noexcept
{
    int ncomp = m_lobc_orig.size();
    for (int n = 0; n < ncomp; ++n) {
        for (int idim = 0; idim <AMREX_SPACEDIM; ++idim) {
            if (m_lobc_orig[n][idim] == bct || m_hibc_orig[n][idim] == bct) {
                return true;
            }
        }
    }
    return false;
}
 
template <typename MF>
bool
MLLinOpT<MF>::hasInhomogNeumannBC () const noexcept
{
    return hasBC(BCType::inhomogNeumann);
}
 
template <typename MF>
bool
MLLinOpT<MF>::hasRobinBC () const noexcept
{
    return hasBC(BCType::Robin);
}
 
template <typename MF>
Box
MLLinOpT<MF>::compactify (Box const& b) const noexcept
{
#if (AMREX_SPACEDIM == 3)
    if (info.hasHiddenDimension()) {
        const auto& lo = b.smallEnd();
        const auto& hi = b.bigEnd();
        if (info.hidden_direction == 0) {
            return Box(IntVect(lo[1],lo[2],0), IntVect(hi[1],hi[2],0), b.ixType());
        } else if (info.hidden_direction == 1) {
            return Box(IntVect(lo[0],lo[2],0), IntVect(hi[0],hi[2],0), b.ixType());
        } else {
            return Box(IntVect(lo[0],lo[1],0), IntVect(hi[0],hi[1],0), b.ixType());
        }
    } else
#endif
    {
        return b;
    }
}
 
template <typename MF>
void
MLLinOpT<MF>::makeAgglomeratedDMap (const Vector<BoxArray>& ba,
                                    Vector<DistributionMapping>& dm)
{
    BL_PROFILE("MLLinOp::makeAgglomeratedDMap");
 
    BL_ASSERT(!dm[0].empty());
    for (int i = 1, N=static_cast<int>(ba.size()); i < N; ++i)
    {
        if (dm[i].empty())
        {
            const std::vector< std::vector<int> >& sfc = DistributionMapping::makeSFC(ba[i]);
 
            const int nprocs = ParallelContext::NProcsSub();
            AMREX_ASSERT(std::ssize(sfc) == nprocs);
 
            Vector<int> pmap(ba[i].size());
            for (int iproc = 0; iproc < nprocs; ++iproc) {
                int grank = ParallelContext::local_to_global_rank(iproc);
                for (int ibox : sfc[iproc]) {
                    pmap[ibox] = grank;
                }
            }
            dm[i].define(std::move(pmap));
        }
    }
}
 
template <typename MF>
void
MLLinOpT<MF>::makeConsolidatedDMap (const Vector<BoxArray>& ba,
                                    Vector<DistributionMapping>& dm,
                                    int ratio, int strategy)
{
    BL_PROFILE("MLLinOp::makeConsolidatedDMap()");
 
    int factor = 1;
    BL_ASSERT(!dm[0].empty());
    for (int i = 1, N=static_cast<int>(ba.size()); i < N; ++i)
    {
        if (dm[i].empty())
        {
            factor *= ratio;
 
            const int nprocs = ParallelContext::NProcsSub();
            const auto& pmap_fine = dm[i-1].ProcessorMap();
            Vector<int> pmap(pmap_fine.size());
            ParallelContext::global_to_local_rank(pmap.data(), pmap_fine.data(), static_cast<int>(pmap.size()));
            if (strategy == 1) {
                for (auto& x: pmap) {
                    x /= ratio;
                }
            } else if (strategy == 2) {
                int nprocs_con = static_cast<int>(std::ceil(static_cast<Real>(nprocs)
                                                            / static_cast<Real>(factor)));
                for (auto& x: pmap) {
                    auto d = std::div(x,nprocs_con);
                    x = d.rem;
                }
            } else if (strategy == 3) {
                if (factor == ratio) {
                    const std::vector< std::vector<int> >& sfc = DistributionMapping::makeSFC(ba[i]);
                    for (int iproc = 0; iproc < nprocs; ++iproc) {
                        for (int ibox : sfc[iproc]) {
                            pmap[ibox] = iproc;
                        }
                    }
                }
                for (auto& x: pmap) {
                    x /= ratio;
                }
            }
 
            if (ParallelContext::CommunicatorSub() == ParallelDescriptor::Communicator()) {
                dm[i].define(std::move(pmap));
            } else {
                Vector<int> pmap_g(pmap.size());
                ParallelContext::local_to_global_rank(pmap_g.data(), pmap.data(), static_cast<int>(pmap.size()));
                dm[i].define(std::move(pmap_g));
            }
        }
    }
}
 
template <typename MF>
MPI_Comm
MLLinOpT<MF>::makeSubCommunicator (const DistributionMapping& dm)
{
    BL_PROFILE("MLLinOp::makeSubCommunicator()");
 
#ifdef BL_USE_MPI
 
    Vector<int> newgrp_ranks = dm.ProcessorMap();
    std::ranges::sort(newgrp_ranks);
    auto last = std::unique(newgrp_ranks.begin(), newgrp_ranks.end());
    newgrp_ranks.erase(last, newgrp_ranks.end());
 
    MPI_Comm newcomm;
    MPI_Group defgrp, newgrp;
    MPI_Comm_group(m_default_comm, &defgrp);
    if (ParallelContext::CommunicatorSub() == ParallelDescriptor::Communicator()) {
        MPI_Group_incl(defgrp, static_cast<int>(newgrp_ranks.size()), newgrp_ranks.data(), &newgrp);
    } else {
        Vector<int> local_newgrp_ranks(newgrp_ranks.size());
        ParallelContext::global_to_local_rank(local_newgrp_ranks.data(),
                                              newgrp_ranks.data(), static_cast<int>(newgrp_ranks.size()));
        MPI_Group_incl(defgrp, static_cast<int>(local_newgrp_ranks.size()), local_newgrp_ranks.data(), &newgrp);
    }
 
    MPI_Comm_create(m_default_comm, newgrp, &newcomm);
 
    m_raii_comm = std::make_unique<CommContainer>(newcomm);
 
    MPI_Group_free(&defgrp);
    MPI_Group_free(&newgrp);
 
    return newcomm;
#else
    amrex::ignore_unused(dm);
    return m_default_comm;
#endif
}
 
template <typename MF>
void
MLLinOpT<MF>::setDomainBCLoc (const Array<Real,AMREX_SPACEDIM>& lo_bcloc,
                              const Array<Real,AMREX_SPACEDIM>& hi_bcloc) noexcept
{
    m_domain_bloc_lo = lo_bcloc;
    m_domain_bloc_hi = hi_bcloc;
}
 
template <typename MF>
void
MLLinOpT<MF>::setCoarseFineBC (const MF* crse, int crse_ratio,
                               LinOpBCType bc_type) noexcept
{
    setCoarseFineBC(crse, IntVect(crse_ratio), bc_type);
}
 
template <typename MF>
void
MLLinOpT<MF>::setCoarseFineBC (const MF* crse, IntVect const& crse_ratio,
                               LinOpBCType bc_type) noexcept
{
    m_coarse_data_for_bc = crse;
    m_coarse_data_crse_ratio = crse_ratio;
    m_coarse_fine_bc_type = bc_type;
}
 
template <typename MF>
template <typename AMF>
requires (!std::same_as<MF,AMF>)
void
MLLinOpT<MF>::setCoarseFineBC (const AMF* crse, int crse_ratio,
                               LinOpBCType bc_type) noexcept
{
    setCoarseFineBC(crse, IntVect(crse_ratio), bc_type);
}
 
template <typename MF>
template <typename AMF>
requires (!std::same_as<MF,AMF>)
void
MLLinOpT<MF>::setCoarseFineBC (const AMF* crse, IntVect const& crse_ratio,
                               LinOpBCType bc_type) noexcept
{
    if (crse) {
        m_coarse_data_for_bc_raii = MF(crse->boxArray(), crse->DistributionMap(),
                                       crse->nComp(), crse->nGrowVect());
        m_coarse_data_for_bc_raii.LocalCopy(*crse, 0, 0, crse->nComp(),
                                            crse->nGrowVect());
        m_coarse_data_for_bc = &m_coarse_data_for_bc_raii;
    } else {
        m_coarse_data_for_bc = nullptr;
    }
    m_coarse_data_crse_ratio = crse_ratio;
    m_coarse_fine_bc_type = bc_type;
}
 
template <typename MF>
void
MLLinOpT<MF>::make (Vector<Vector<MF> >& mf, IntVect const& ng) const
{
    mf.clear();
    mf.resize(m_num_amr_levels);
    for (int alev = 0; alev < m_num_amr_levels; ++alev) {
        mf[alev].resize(m_num_mg_levels[alev]);
        for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) {
            mf[alev][mlev] = make(alev, mlev, ng);
        }
    }
}
 
template <typename MF>
MF
MLLinOpT<MF>::make (int amrlev, int mglev, IntVect const& ng) const
{
    return make(amrlev, mglev, ng, MFInfo());
}
 
template <typename MF>
MF
MLLinOpT<MF>::make (int amrlev, int mglev, IntVect const& ng, MFInfo const& mf_info) const
{
    if constexpr (IsMultiFabLike_v<MF>) {
        return MF(amrex::convert(m_grids[amrlev][mglev], m_ixtype),
                  m_dmap[amrlev][mglev], getNComp(), ng, mf_info,
                  *m_factory[amrlev][mglev]);
    } else {
        amrex::ignore_unused(amrlev, mglev, ng);
        amrex::Abort("MLLinOpT::make: how did we get here?");
        return {};
    }
}
 
template <typename MF>
MF
MLLinOpT<MF>::makeAlias (MF const& mf) const
{
    if constexpr (IsMultiFabLike_v<MF>) {
        return MF(mf, amrex::make_alias, 0, mf.nComp());
    } else {
        amrex::ignore_unused(mf);
        amrex::Abort("MLLinOpT::makeAlias: how did we get here?");
        return {};
    }
}
 
template <typename MF>
MF
MLLinOpT<MF>::makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const
{
    return makeCoarseMG(amrlev, mglev, ng, MFInfo());
}
 
template <typename MF>
MF
MLLinOpT<MF>::makeCoarseMG (int amrlev, int mglev, IntVect const& ng,
                            MFInfo const& mf_info) const
{
    if constexpr (IsMultiFabLike_v<MF>) {
        BoxArray cba = m_grids[amrlev][mglev];
        IntVect ratio = (amrlev > 0) ? IntVect(2) : mg_coarsen_ratio_vec[mglev];
        if (hasHiddenDimension()) { ratio[hiddenDirection()] = 1; }
        cba.coarsen(ratio);
        cba.convert(m_ixtype);
        return MF(cba, m_dmap[amrlev][mglev], getNComp(), ng, mf_info);
    } else {
        amrex::ignore_unused(amrlev, mglev, ng);
        amrex::Abort("MLLinOpT::makeCoarseMG: how did we get here?");
        return {};
    }
}
 
template <typename MF>
MF
MLLinOpT<MF>::makeCoarseAmr (int famrlev, IntVect const& ng) const
{
    return makeCoarseAmr(famrlev, ng, MFInfo());
}
 
template <typename MF>
MF
MLLinOpT<MF>::makeCoarseAmr (int famrlev, IntVect const& ng, MFInfo const& mf_info) const
{
    if constexpr (IsMultiFabLike_v<MF>) {
        BoxArray cba = m_grids[famrlev][0];
        IntVect ratio(AMRRefRatioVect(famrlev-1));
        cba.coarsen(ratio);
        cba.convert(m_ixtype);
        return MF(cba, m_dmap[famrlev][0], getNComp(), ng, mf_info);
    } else {
        amrex::ignore_unused(famrlev, ng);
        amrex::Abort("MLLinOpT::makeCoarseAmr: how did we get here?");
        return {};
    }
}
 
template <typename MF>
void
MLLinOpT<MF>::resizeMultiGrid (int new_size)
{
    if (new_size <= 0 || new_size >= m_num_mg_levels[0]) { return; }
 
    m_num_mg_levels[0] = new_size;
 
    m_geom[0].resize(new_size);
    m_grids[0].resize(new_size);
    m_dmap[0].resize(new_size);
    m_factory[0].resize(new_size);
 
    if (m_bottom_comm != m_default_comm) {
        m_bottom_comm = makeSubCommunicator(m_dmap[0].back());
    }
}
 
template <typename MF>
void
MLLinOpT<MF>::avgDownResMG (int clev, MF& cres, MF const& fres) const
{
    amrex::ignore_unused(clev, cres, fres);
    if constexpr (amrex::IsFabArray<MF>::value) {
        const int ncomp = this->getNComp();
#ifdef AMREX_USE_EB
        if (!fres.isAllRegular()) {
            if constexpr (std::is_same<MF,MultiFab>()) {
                amrex::EB_average_down(fres, cres, 0, ncomp,
                                       mg_coarsen_ratio_vec[clev-1]);
            } else {
                amrex::Abort("EB_average_down only works with MultiFab");
            }
        } else
#endif
        {
            amrex::average_down(fres, cres, 0, ncomp, mg_coarsen_ratio_vec[clev-1]);
        }
    } else {
        amrex::Abort("For non-FabArray, MLLinOpT<MF>::avgDownResMG should be overridden.");
    }
}
 
template <typename MF>
bool
MLLinOpT<MF>::isMFIterSafe (int amrlev, int mglev1, int mglev2) const
{
    return m_dmap[amrlev][mglev1] == m_dmap[amrlev][mglev2]
        && BoxArray::SameRefs(m_grids[amrlev][mglev1], m_grids[amrlev][mglev2]);
}
 
template <typename MF>
template <MultiFabLike AMF>
requires (!std::same_as<MF,AMF>)
void
MLLinOpT<MF>::setLevelBC (int amrlev, const AMF* levelbcdata,
                          const AMF* robinbc_a, const AMF* robinbc_b,
                          const AMF* robinbc_f)
{
    const int ncomp = this->getNComp();
    if (levelbcdata) {
        levelbc_raii[amrlev] = std::make_unique<MF>(levelbcdata->boxArray(),
                                                    levelbcdata->DistributionMap(),
                                                    ncomp, levelbcdata->nGrowVect());
        levelbc_raii[amrlev]->LocalCopy(*levelbcdata, 0, 0, ncomp,
                                        levelbcdata->nGrowVect());
    } else {
        levelbc_raii[amrlev].reset();
    }
 
    if (robinbc_a) {
        robin_a_raii[amrlev] = std::make_unique<MF>(robinbc_a->boxArray(),
                                                    robinbc_a->DistributionMap(),
                                                    ncomp, robinbc_a->nGrowVect());
        robin_a_raii[amrlev]->LocalCopy(*robinbc_a, 0, 0, ncomp,
                                        robinbc_a->nGrowVect());
    } else {
        robin_a_raii[amrlev].reset();
    }
 
    if (robinbc_b) {
        robin_b_raii[amrlev] = std::make_unique<MF>(robinbc_b->boxArray(),
                                                    robinbc_b->DistributionMap(),
                                                    ncomp, robinbc_b->nGrowVect());
        robin_b_raii[amrlev]->LocalCopy(*robinbc_b, 0, 0, ncomp,
                                        robinbc_b->nGrowVect());
    } else {
        robin_b_raii[amrlev].reset();
    }
 
    if (robinbc_f) {
        robin_f_raii[amrlev] = std::make_unique<MF>(robinbc_f->boxArray(),
                                                    robinbc_f->DistributionMap(),
                                                    ncomp, robinbc_f->nGrowVect());
        robin_f_raii[amrlev]->LocalCopy(*robinbc_f, 0, 0, ncomp,
                                        robinbc_f->nGrowVect());
    } else {
        robin_f_raii[amrlev].reset();
    }
 
    this->setLevelBC(amrlev, levelbc_raii[amrlev].get(), robin_a_raii[amrlev].get(),
                     robin_b_raii[amrlev].get(), robin_f_raii[amrlev].get());
}
 
template <typename MF>
auto
MLLinOpT<MF>::dotProductPrecond (Vector<MF const*> const& x, Vector<MF const*> const& y) const -> RT
{
    AMREX_ALWAYS_ASSERT(NAMRLevels() == 1);
    return xdoty(0,0,*x[0],*y[0],false);
}
 
template <typename MF>
auto
MLLinOpT<MF>::norm2Precond (Vector<MF const*> const& x) const -> RT
{
    AMREX_ALWAYS_ASSERT(NAMRLevels() == 1);
    auto r = xdoty(0,0,*x[0],*x[0],false);
    return std::sqrt(r);
}
 
extern template class MLLinOpT<MultiFab>;

using MLLinOp = MLLinOpT<MultiFab>;
 
}
 
#endif