AMReX_MLEBNodeFDLaplacian.H

Go to the documentation of this file.

#ifndef AMREX_MLEBNODEFDLAPLACIAN_H_
#define AMREX_MLEBNODEFDLAPLACIAN_H_
#include <AMReX_Config.H>
 
#include <AMReX_Array.H>
#ifdef AMREX_USE_EB
#include <AMReX_EBFabFactory.H>
#endif
#include <AMReX_LayoutData.H>
#include <AMReX_MLNodeLinOp.H>
 
#include <limits>
 
namespace amrex {
 
// Although the class has EB in the name, it works for non-EB build too.
//
// del dot (sigma grad phi) = rhs, for non-RZ
// where phi and rhs are nodal multifab, and sigma is a tensor constant
// with only diagonal components.  The EB is assumed to be Dirichlet.
//
// del dot (sigma grad phi) - alpha/r^2 phi = rhs, for RZ where alpha is a
// scalar constant that is zero by default.
//
// New feature: sigma can also be a single-component cell-centered multifab.
//
// The EB information needed on the coarse MG levels (covered nodes and open
// EB positions on edges) is built here from the level-0 factory, so multigrid
// coarsening is not limited by the EB index space.  Coarsening stops at the
// last level that still has a handful of open nodes.
 
class MLEBNodeFDLaplacian
    : public MLNodeLinOp
{
public:
 
    MLEBNodeFDLaplacian () = default;
 
#ifdef AMREX_USE_EB
    MLEBNodeFDLaplacian (const Vector<Geometry>& a_geom,
                         const Vector<BoxArray>& a_grids,
                         const Vector<DistributionMapping>& a_dmap,
                         const LPInfo& a_info,
                         const Vector<EBFArrayBoxFactory const*>& a_factory);
#endif
 
    MLEBNodeFDLaplacian (const Vector<Geometry>& a_geom,
                         const Vector<BoxArray>& a_grids,
                         const Vector<DistributionMapping>& a_dmap,
                         const LPInfo& a_info);
 
    ~MLEBNodeFDLaplacian () override = default;
 
    MLEBNodeFDLaplacian (const MLEBNodeFDLaplacian&) = delete;
    MLEBNodeFDLaplacian (MLEBNodeFDLaplacian&&) = delete;
    MLEBNodeFDLaplacian& operator= (const MLEBNodeFDLaplacian&) = delete;
    MLEBNodeFDLaplacian& operator= (MLEBNodeFDLaplacian&&) = delete;
 
    void setSigma (Array<Real,AMREX_SPACEDIM> const& a_sigma) noexcept;
 
    void setSigma (int amrlev, MultiFab const& a_sigma);
 
    void setRZ (bool flag);
 
    void setAlpha (Real a_alpha);
 
#ifdef AMREX_USE_EB
 
    void setEBDirichlet (Real a_phi_eb);
    template <typename F>
    requires (IsCallableR<Real,F,AMREX_D_DECL(Real,Real,Real)>::value)
    void setEBDirichlet (F const& f);
 
    void define (const Vector<Geometry>& a_geom,
                 const Vector<BoxArray>& a_grids,
                 const Vector<DistributionMapping>& a_dmap,
                 const LPInfo& a_info,
                 const Vector<EBFArrayBoxFactory const*>& a_factory);
 
    [[nodiscard]] bool scaleRHS (int amrlev, MultiFab* rhs) const final;
 
#endif
 
    void define (const Vector<Geometry>& a_geom,
                 const Vector<BoxArray>& a_grids,
                 const Vector<DistributionMapping>& a_dmap,
                 const LPInfo& a_info);
 
    [[nodiscard]] std::string name () const override { return std::string("MLEBNodeFDLaplacian"); }
 
    void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const final;
    void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const final;
 
    [[nodiscard]] bool needsUpdate () const override {
        return (m_needs_update || MLNodeLinOp::needsUpdate());
    }
    void update () override;
 
    void prepareForSolve () final;
    void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const final;
    void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs) const final;
    void normalize (int amrlev, int mglev, MultiFab& mf) const final;
 
    void fixUpResidualMask (int amrlev, iMultiFab& resmsk) final;
 
    [[nodiscard]] bool isSingular (int) const final { return false; }
    [[nodiscard]] bool isBottomSingular () const final { return false; }
 
    void compGrad (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>& grad,
                   MultiFab& sol, Location /*loc*/) const override;
 
    void compGrad_doit (int amrlev, const Array<MultiFab*,AMREX_SPACEDIM>& grad,
                        MultiFab& sol) const;
 
#if (AMREX_SPACEDIM > 1)
#if defined(AMREX_USE_HYPRE)
    void fillIJMatrix (MFIter const& mfi,
                       Array4<HypreNodeLap::AtomicInt const> const& gid,
                       Array4<int const> const& lid,
                       HypreNodeLap::Int* ncols,
                       HypreNodeLap::Int* cols,
                       Real* mat) const override;
#endif
 
    [[nodiscard]] bool supportsAnisotropicCoarsening () const override { return true; }
 
    [[nodiscard]] GpuArray<Real,AMREX_SPACEDIM>
    anisotropicCoarseningCellSize (Geometry const& geom) const override;
 
    [[nodiscard]] bool supportsAlgMG () const override { return true; }
 
    void fillAlgMatrix (int mglev, MFIter const& mfi,
                        Array4<Long const> const& gid,
                        Array4<int const> const& lid,
                        Long* ncols, Long* cols, Real* mat) const override;
 
    void fillRHS (int mglev, MFIter const& mfi,
                  Array4<int const> const& lid,
                  Real* rhs,
                  Array4<Real const> const& bfab) const override;
 
    // Public for nvcc.
    template <typename AlgInt, typename AlgGid>
    void fillMatrix_doit (int mglev, MFIter const& mfi,
                          Array4<AlgGid const> const& gid,
                          Array4<int const> const& lid,
                          AlgInt* ncols, AlgInt* cols, Real* mat) const;
#endif
 
    void postSolve (Vector<MultiFab*> const& sol) const override;
 
    [[nodiscard]] BottomSolver getDefaultBottomSolver () const override {
        return BottomSolver::custom;
    }
 
    [[nodiscard]] bool supportCustomBottomSolver () const override { return true; }
 
    void customBottomSolve (MLMGT<MultiFab>* mlmg, MultiFab& x, const MultiFab& b,
                            Real eps_rel, Real eps_abs, int maxiter) override;
 
protected:
    void buildMGHierarchy () override;
 
private:
    GpuArray<Real,AMREX_SPACEDIM> m_sigma{{AMREX_D_DECL(1_rt,1_rt,1_rt)}};
    Vector<std::unique_ptr<MultiFab>> m_sigma_mf;
    Vector<Vector<Array<MultiFab,AMREX_SPACEDIM>>> m_sigma_edge;
    bool m_has_sigma_mf = false;
    bool m_needs_update = true;
    Real m_s_phi_eb = std::numeric_limits<Real>::lowest();
    Vector<MultiFab> m_phi_eb;
    int m_rz = false;
    Real m_rz_alpha = 0._rt;
 
#ifdef AMREX_USE_EB
    Vector<Vector<MultiFab>> m_levset;
    Vector<Vector<Array<MultiFab,AMREX_SPACEDIM>>> m_eb_pos;
    Vector<Vector<LayoutData<int>>> m_has_eb;
#ifdef AMREX_USE_GPU
    Vector<Vector<Gpu::DeviceVector<int>>> m_has_eb_d;
#endif
    Vector<Vector<int>> m_eb_lost;
    Vector<Vector<MultiFab>> m_row_scale;
    Vector<Vector<MultiFab>> m_row_scale_crse;
#endif
 
public:
    // These are public only because they launch GPU kernels.
#ifdef AMREX_USE_EB
    void build_eb_data ();
    void limit_coarsening ();
#endif
    void update_sigma ();
};
 
#ifdef AMREX_USE_EB
 
template <typename F>
requires (IsCallableR<Real,F,AMREX_D_DECL(Real,Real,Real)>::value)
void MLEBNodeFDLaplacian::setEBDirichlet (F const& f)
{
    // Select m_phi_eb again, in case a previous prepareForSolve defaulted
    // to homogeneous Dirichlet or a scalar value was set earlier.
    m_s_phi_eb = std::numeric_limits<Real>::lowest();
    m_phi_eb.resize(m_num_amr_levels);
    for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) {
        auto const* factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][0].get());
        if (factory) {
            Geometry const& geom = m_geom[amrlev][0];
            auto const problo = geom.ProbLoArray();
            auto const cellsize = geom.CellSizeArray();
            if (m_phi_eb[amrlev].empty()) {
                m_phi_eb[amrlev].define(amrex::convert(m_grids[amrlev][0],IntVect(1)),
                                        m_dmap[amrlev][0], 1, 1);
                m_phi_eb[amrlev].setVal(0.0);
            }
            auto const& flags = factory->getMultiEBCellFlagFab();
            auto const& levset = factory->getLevelSet();
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
            for (MFIter mfi(m_phi_eb[amrlev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
            {
                const Box& ndbx = mfi.growntilebox();
                const auto& flag = flags[mfi];
                if (flag.getType() != FabType::regular) {
                    Array4<Real const> const lstarr = levset.const_array(mfi);
                    Array4<Real> const& phi = m_phi_eb[amrlev].array(mfi);
                    AMREX_HOST_DEVICE_FOR_3D(ndbx, i, j, k,
                    {
                        if (lstarr(i,j,k) >= Real(0.0)) {
                            phi(i,j,k) = f(AMREX_D_DECL(problo[0]+Real(i)*cellsize[0],
                                                        problo[1]+Real(j)*cellsize[1],
                                                        problo[2]+Real(k)*cellsize[2]));
                        }
                    });
                }
            }
        }
    }
}
 
#endif
 
}
 
#endif