AMReX_MLCellABecLap.H

Go to the documentation of this file.

#ifndef AMREX_ML_CELL_ABECLAP_H_
#define AMREX_ML_CELL_ABECLAP_H_
#include <AMReX_Config.H>
 
#include <AMReX_MLCellLinOp.H>
#include <AMReX_MLCellABecLap_K.H>
 
namespace amrex {
 
template <typename MF>
class MLCellABecLapT  // NOLINT(cppcoreguidelines-virtual-class-destructor)
    : public MLCellLinOpT<MF>
{
public:
 
    using FAB = typename MF::fab_type;
    using RT  = typename MF::value_type;
 
    using Location  = typename MLLinOpT<MF>::Location;
 
    MLCellABecLapT () = default;
    ~MLCellABecLapT () override = default;
 
    MLCellABecLapT (const MLCellABecLapT<MF>&) = delete;
    MLCellABecLapT (MLCellABecLapT<MF>&&) = delete;
    MLCellABecLapT<MF>& operator= (const MLCellABecLapT<MF>&) = delete;
    MLCellABecLapT<MF>& operator= (MLCellABecLapT<MF>&&) = delete;
 
    void define (const Vector<Geometry>& a_geom,
                 const Vector<BoxArray>& a_grids,
                 const Vector<DistributionMapping>& a_dmap,
                 const LPInfo& a_info = LPInfo(),
                 const Vector<FabFactory<FAB> const*>& a_factory = {});
 
    void define (const Vector<Geometry>& a_geom,
                 const Vector<BoxArray>& a_grids,
                 const Vector<DistributionMapping>& a_dmap,
                 const Vector<iMultiFab const*>& a_overset_mask,
                 const LPInfo& a_info = LPInfo(),
                 const Vector<FabFactory<FAB> const*>& a_factory = {});
 
    [[nodiscard]] iMultiFab const* getOversetMask (int amrlev, int mglev) const {
        return m_overset_mask[amrlev][mglev].get();
    }
 
    [[nodiscard]] bool needsUpdate () const override {
        return MLCellLinOpT<MF>::needsUpdate();
    }
    void update () override;
 
    void prepareForSolve () override;
 
    void setDirichletNodesToZero (int amrlev, int mglev, MF& mf) const override;
 
    void getFluxes (const Vector<Array<MF*,AMREX_SPACEDIM> >& a_flux,
                            const Vector<MF*>& a_sol,
                            Location a_loc) const final;
    void getFluxes (const Vector<MF*>& a_flux,
                            const Vector<MF*>& a_sol) const final {
        amrex::ignore_unused(a_flux, a_sol);
        amrex::Abort("MLCellABecLap::getFluxes: How did we get here?");
    }
 
    virtual RT getAScalar () const = 0;
    virtual RT getBScalar () const = 0;
    virtual MF const* getACoeffs (int amrlev, int mglev) const = 0;
    virtual Array<MF const*,AMREX_SPACEDIM> getBCoeffs (int amrlev, int mglev) const = 0;
 
    void applyInhomogNeumannTerm (int amrlev, MF& rhs) const final;
 
    void addInhomogNeumannFlux (
        int amrlev, const Array<MF*,AMREX_SPACEDIM>& grad,
        MF const& sol, bool mult_bcoef) const final;
 
    void applyOverset (int amrlev, MF& rhs) const override;
 
    [[nodiscard]] virtual MF const* getEBBCoeffs (int /*amrlev*/, int /*mglev*/) const { return nullptr; }
 
    [[nodiscard]] std::unique_ptr<MLAlgMG> makeAlgMG (int mglev) const override;
    [[nodiscard]] bool supportsAlgMG () const override { return std::is_same_v<MF,MultiFab>; }
 
    [[nodiscard]] Array<MF,AMREX_SPACEDIM> makeDefaultBCoeffs (int mglev) const;
 
#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
    [[nodiscard]] std::unique_ptr<Hypre> makeHypre (Hypre::Interface hypre_interface) const override;
#endif
 
#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
    [[nodiscard]] std::unique_ptr<PETScABecLap> makePETSc () const override;
#endif
 
protected:
    Vector<Vector<std::unique_ptr<iMultiFab> > > m_overset_mask;
 
    LPInfo m_lpinfo_arg;
 
    [[nodiscard]] bool supportInhomogNeumannBC () const noexcept override { return true; }
};
 
template <typename MF>
void
MLCellABecLapT<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)
{
    MLCellLinOpT<MF>::define(a_geom, a_grids, a_dmap, a_info, a_factory);
 
    this->m_overset_mask.resize(this->m_num_amr_levels);
    for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
        this->m_overset_mask[amrlev].resize(this->m_num_mg_levels[amrlev]);
    }
}
 
template <typename MF>
void
MLCellABecLapT<MF>::define (const Vector<Geometry>& a_geom,
                            const Vector<BoxArray>& a_grids,
                            const Vector<DistributionMapping>& a_dmap,
                            const Vector<iMultiFab const*>& a_overset_mask,
                            const LPInfo& a_info,
                            const Vector<FabFactory<FAB> const*>& a_factory)
{
    BL_PROFILE("MLCellABecLap::define(overset)");
 
    AMREX_ALWAYS_ASSERT(!this->hasHiddenDimension());
 
    this->m_lpinfo_arg = a_info;
 
    auto namrlevs = static_cast<int>(a_geom.size());
    this->m_overset_mask.resize(namrlevs);
    for (int amrlev = 0; amrlev < namrlevs; ++amrlev)
    {
        this->m_overset_mask[amrlev].push_back(std::make_unique<iMultiFab>(a_grids[amrlev],
                                                                     a_dmap[amrlev], 1, 1));
        iMultiFab::Copy(*(this->m_overset_mask[amrlev][0]), *a_overset_mask[amrlev], 0, 0, 1, 0);
        if (amrlev > 1) {
            AMREX_ALWAYS_ASSERT(amrex::refine(a_geom[amrlev-1].Domain(),2)
                                == a_geom[amrlev].Domain());
        }
    }
 
    int amrlev = 0;
    Box dom = a_geom[0].Domain();
    for (int mglev = 1; mglev <= a_info.max_coarsening_level; ++mglev)
    {
        AMREX_ALWAYS_ASSERT(this->mg_coarsen_ratio == 2);
        iMultiFab const& fine = *(this->m_overset_mask[amrlev][mglev-1]);
        if (dom.coarsenable(2) && fine.boxArray().coarsenable(2)) {
            dom.coarsen(2);
            auto crse = std::make_unique<iMultiFab>(amrex::coarsen(fine.boxArray(),2),
                                                    fine.DistributionMap(), 1, 1);
            ReduceOps<ReduceOpSum> reduce_op;
            ReduceData<int> reduce_data(reduce_op);
            using ReduceTuple = typename decltype(reduce_data)::Type;
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
            for (MFIter mfi(*crse, TilingIfNotGPU()); mfi.isValid(); ++mfi)
            {
                const Box& bx = mfi.tilebox();
                Array4<int const> const& fmsk = fine.const_array(mfi);
                Array4<int> const& cmsk = crse->array(mfi);
                reduce_op.eval(bx, reduce_data,
                [=] AMREX_GPU_HOST_DEVICE (Box const& b) -> ReduceTuple
                {
                    return { coarsen_overset_mask(b, cmsk, fmsk) };
                });
            }
            ReduceTuple hv = reduce_data.value(reduce_op);
            if (amrex::get<0>(hv) == 0) {
                this->m_overset_mask[amrlev].push_back(std::move(crse));
            } else {
                break;
            }
        } else {
            break;
        }
    }
    int max_overset_mask_coarsening_level = this->m_overset_mask[amrlev].size()-1;
    ParallelAllReduce::Min(max_overset_mask_coarsening_level, ParallelContext::CommunicatorSub());
    this->m_overset_mask[amrlev].resize(max_overset_mask_coarsening_level+1);
 
    LPInfo linfo = a_info;
    linfo.max_coarsening_level = std::min(a_info.max_coarsening_level,
                                          max_overset_mask_coarsening_level);
    // The overset masks above are coarsened isotropically by 2, so the MG
    // levels must be too, or the masks and the grids end up in different
    // index spaces.
    linfo.do_semicoarsening = false;
 
    MLCellLinOpT<MF>::define(a_geom, a_grids, a_dmap, linfo, a_factory);
 
    amrlev = 0;
    for (int mglev = 1; mglev < this->m_num_mg_levels[amrlev]; ++mglev) {
        MF foo(this->m_grids[amrlev][mglev], this->m_dmap[amrlev][mglev], 1, 0, MFInfo().SetAlloc(false));
        if (! amrex::isMFIterSafe(*(this->m_overset_mask[amrlev][mglev]), foo)) {
            auto osm = std::make_unique<iMultiFab>(this->m_grids[amrlev][mglev],
                                                   this->m_dmap[amrlev][mglev], 1, 1);
            osm->ParallelCopy(*(this->m_overset_mask[amrlev][mglev]));
            std::swap(osm, this->m_overset_mask[amrlev][mglev]);
        }
    }
 
    for (amrlev = 1; amrlev < this->m_num_amr_levels; ++amrlev) {
        for (int mglev = 1; mglev < this->m_num_mg_levels[amrlev]; ++mglev) { // for ref_ratio 4
            this->m_overset_mask[amrlev].push_back(std::make_unique<iMultiFab>(this->m_grids[amrlev][mglev],
                                                                         this->m_dmap[amrlev][mglev],
                                                                         1, 1));
 
#ifdef AMREX_USE_GPU
            if (Gpu::inLaunchRegion() && this->m_overset_mask[amrlev][mglev]->isFusingCandidate()) {
                auto const& crsema = this->m_overset_mask[amrlev][mglev]->arrays();
                auto const& finema = this->m_overset_mask[amrlev][mglev-1]->const_arrays();
                ParallelFor(*(this->m_overset_mask[amrlev][mglev]),
                [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
                {
                    coarsen_overset_mask(i,j,k, crsema[box_no], finema[box_no]);
                });
                if (!Gpu::inNoSyncRegion()) {
                    Gpu::streamSynchronize();
                }
            } else
#endif
            {
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
                for (MFIter mfi(*(this->m_overset_mask[amrlev][mglev]), TilingIfNotGPU()); mfi.isValid(); ++mfi)
                {
                    const Box& bx = mfi.tilebox();
                    Array4<int> const& cmsk = this->m_overset_mask[amrlev][mglev]->array(mfi);
                    Array4<int const> const fmsk = this->m_overset_mask[amrlev][mglev-1]->const_array(mfi);
                    AMREX_HOST_DEVICE_PARALLEL_FOR_3D(bx, i, j, k,
                    {
                        coarsen_overset_mask(i,j,k, cmsk, fmsk);
                    });
                }
            }
        }
    }
 
    for (amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
        for (int mglev = 0; mglev < this->m_num_mg_levels[amrlev]; ++mglev) {
            this->m_overset_mask[amrlev][mglev]->setBndry(1);
            this->m_overset_mask[amrlev][mglev]->FillBoundary(this->m_geom[amrlev][mglev].periodicity());
        }
    }
}
 
template <typename MF>
void
MLCellABecLapT<MF>::update ()
{
    if (MLCellLinOpT<MF>::needsUpdate()) { MLCellLinOpT<MF>::update(); }
}
 
template <typename MF>
void
MLCellABecLapT<MF>::prepareForSolve ()
{
    MLCellLinOpT<MF>::prepareForSolve();
}
 
template <typename MF>
void
MLCellABecLapT<MF>::setDirichletNodesToZero (int amrlev, int mglev, MF& mf) const
{
    auto const* omask = this->getOversetMask(amrlev, mglev);
    if (omask) {
        const int ncomp = this->getNComp();
        auto const& mskma = omask->const_arrays();
        auto const& ma = mf.arrays();
        ParallelFor(mf, IntVect(0), ncomp,
        [=] AMREX_GPU_DEVICE (int bno, int i, int j, int k, int n)
        {
            if (mskma[bno](i,j,k) == 0) { ma[bno](i,j,k,n) = RT(0.0); }
        });
        if (!Gpu::inNoSyncRegion()) {
            Gpu::streamSynchronize();
        }
    }
}
 
template <typename MF>
void
MLCellABecLapT<MF>::getFluxes (const Vector<Array<MF*,AMREX_SPACEDIM> >& a_flux,
                               const Vector<MF*>& a_sol,
                               Location a_loc) const
{
    BL_PROFILE("MLMG::getFluxes()");
 
    const int ncomp = this->getNComp();
    const RT betainv = RT(1.0) / getBScalar();
    const int nlevs = this->NAMRLevels();
    // Without b coefficients, the metric factor is in the flux kernels, not in
    // b, so it must be removed before the boundary faces are set using b = 1.
    const bool has_bcoef = (getBCoeffs(0,0)[0] != nullptr);
    for (int alev = 0; alev < nlevs; ++alev) {
        this->compFlux(alev, a_flux[alev], *a_sol[alev], a_loc);
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
            if (betainv != RT(1.0)) {
                a_flux[alev][idim]->mult(betainv, 0, ncomp);
            }
            if (!has_bcoef) {
                this->unapplyMetricTerm(alev, 0, *a_flux[alev][idim]);
            }
        }
        this->addInhomogNeumannFlux(alev, a_flux[alev], *a_sol[alev], true);
        if (has_bcoef) {
            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
                this->unapplyMetricTerm(alev, 0, *a_flux[alev][idim]);
            }
        }
    }
}
 
template <typename MF>
void
MLCellABecLapT<MF>::applyInhomogNeumannTerm (int amrlev, MF& rhs) const
{
    bool has_inhomog_neumann = this->hasInhomogNeumannBC();
    bool has_robin = this->hasRobinBC();
 
    // A preconditioner must be linear, so it gets no inhomogeneous terms.
    if (this->m_precond_mode) { return; }
    if (!has_inhomog_neumann && !has_robin) { return; }
 
    int ncomp = this->getNComp();
    const int mglev = 0;
 
    const auto problo = this->m_geom[amrlev][mglev].ProbLoArray();
    const auto probhi = this->m_geom[amrlev][mglev].ProbHiArray();
    amrex::ignore_unused(probhi);
    const RT dxi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(0));
    const RT dyi = static_cast<RT>((AMREX_SPACEDIM >= 2) ? this->m_geom[amrlev][mglev].InvCellSize(1) : Real(1.0));
    const RT dzi = static_cast<RT>((AMREX_SPACEDIM == 3) ? this->m_geom[amrlev][mglev].InvCellSize(2) : Real(1.0));
    const RT xlo = static_cast<RT>(problo[0]);
    const RT dx = static_cast<RT>(this->m_geom[amrlev][mglev].CellSize(0));
    const Box& domain = this->m_geom[amrlev][mglev].Domain();
 
    const RT beta = getBScalar();
    Array<MF const*, AMREX_SPACEDIM> const& bcoef = getBCoeffs(amrlev,mglev);
    FAB foo(Box(IntVect(0),IntVect(1)));
    bool has_bcoef = (bcoef[0] != nullptr);
 
    const auto& maskvals = this->m_maskvals[amrlev][mglev];
    const auto& bcondloc = *(this->m_bcondloc[amrlev][mglev]);
    const auto& bndry = *(this->m_bndry_sol[amrlev]);
 
    MFItInfo mfi_info;
    if (Gpu::notInLaunchRegion()) { mfi_info.SetDynamic(true); }
 
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(rhs, mfi_info); mfi.isValid(); ++mfi)
    {
        const Box& vbx = mfi.validbox();
        auto const& rhsfab = rhs.array(mfi);
 
        const auto & bdlv = bcondloc.bndryLocs(mfi);
        const auto & bdcv = bcondloc.bndryConds(mfi);
 
        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
        {
            auto const bfab = (has_bcoef)
                ? bcoef[idim]->const_array(mfi) : foo.const_array();
            const Orientation olo(idim,Orientation::low);
            const Orientation ohi(idim,Orientation::high);
            const Box blo = amrex::adjCellLo(vbx, idim);
            const Box bhi = amrex::adjCellHi(vbx, idim);
            const auto& mlo = maskvals[olo].array(mfi);
            const auto& mhi = maskvals[ohi].array(mfi);
            const auto& bvlo = bndry.bndryValues(olo).array(mfi);
            const auto& bvhi = bndry.bndryValues(ohi).array(mfi);
            bool outside_domain_lo = !(domain.contains(blo));
            bool outside_domain_hi = !(domain.contains(bhi));
            if ((!outside_domain_lo) && (!outside_domain_hi)) { continue; }
            for (int icomp = 0; icomp < ncomp; ++icomp) {
                const BoundCond bctlo = bdcv[icomp][olo];
                const BoundCond bcthi = bdcv[icomp][ohi];
                const RT bcllo = bdlv[icomp][olo];
                const RT bclhi = bdlv[icomp][ohi];
                if (this->m_lobc_orig[icomp][idim] == LinOpBCType::inhomogNeumann && outside_domain_lo)
                {
                    if (idim == 0) {
                        RT fac = beta*dxi;
                        if (this->m_has_metric_term && !has_bcoef) {
#if (AMREX_SPACEDIM == 1)
                            fac *= static_cast<RT>(problo[0]*problo[0]);
#elif (AMREX_SPACEDIM == 2)
                            fac *= static_cast<RT>(problo[0]);
#endif
                        }
                        AMREX_HOST_DEVICE_FOR_3D(blo, i, j, k,
                        {
                            mllinop_apply_innu_xlo(i,j,k, rhsfab, mlo, bfab,
                                                   bctlo, bcllo, bvlo,
                                                   fac, has_bcoef, icomp);
                        });
                    } else if (idim == 1) {
                        RT fac = beta*dyi;
                        if (this->m_has_metric_term && !has_bcoef) {
                            AMREX_HOST_DEVICE_FOR_3D(blo, i, j, k,
                            {
                                mllinop_apply_innu_ylo_m(i,j,k, rhsfab, mlo,
                                                         bctlo, bcllo, bvlo,
                                                         fac, xlo, dx, icomp);
                            });
                        }
                        else {
                            AMREX_HOST_DEVICE_FOR_3D(blo, i, j, k,
                            {
                                mllinop_apply_innu_ylo(i,j,k, rhsfab, mlo, bfab,
                                                       bctlo, bcllo, bvlo,
                                                       fac, has_bcoef, icomp);
                            });
                        }
                    } else {
                        RT fac = beta*dzi;
                        AMREX_HOST_DEVICE_FOR_3D(blo, i, j, k,
                        {
                            mllinop_apply_innu_zlo(i,j,k, rhsfab, mlo, bfab,
                                                   bctlo, bcllo, bvlo,
                                                   fac, has_bcoef, icomp);
                        });
                    }
                }
                if (this->m_hibc_orig[icomp][idim] == LinOpBCType::inhomogNeumann && outside_domain_hi)
                {
                    if (idim == 0) {
                        RT fac = beta*dxi;
                        if (this->m_has_metric_term && !has_bcoef) {
#if (AMREX_SPACEDIM == 1)
                            fac *= static_cast<RT>(probhi[0]*probhi[0]);
#elif (AMREX_SPACEDIM == 2)
                            fac *= static_cast<RT>(probhi[0]);
#endif
                        }
                        AMREX_HOST_DEVICE_FOR_3D(bhi, i, j, k,
                        {
                            mllinop_apply_innu_xhi(i,j,k, rhsfab, mhi, bfab,
                                                   bcthi, bclhi, bvhi,
                                                   fac, has_bcoef, icomp);
                        });
                    } else if (idim == 1) {
                        RT fac = beta*dyi;
                        if (this->m_has_metric_term && !has_bcoef) {
                            AMREX_HOST_DEVICE_FOR_3D(bhi, i, j, k,
                            {
                                mllinop_apply_innu_yhi_m(i,j,k, rhsfab, mhi,
                                                         bcthi, bclhi, bvhi,
                                                         fac, xlo, dx, icomp);
                            });
                        } else {
                            AMREX_HOST_DEVICE_FOR_3D(bhi, i, j, k,
                            {
                                mllinop_apply_innu_yhi(i,j,k, rhsfab, mhi, bfab,
                                                       bcthi, bclhi, bvhi,
                                                       fac, has_bcoef, icomp);
                            });
                        }
                    } else {
                        RT fac = beta*dzi;
                        AMREX_HOST_DEVICE_FOR_3D(bhi, i, j, k,
                        {
                            mllinop_apply_innu_zhi(i,j,k, rhsfab, mhi, bfab,
                                                   bcthi, bclhi, bvhi,
                                                   fac, has_bcoef, icomp);
                        });
                    }
                }
 
                if (has_robin) {
                    // For Robin BC, see comments in AMReX_MLABecLaplacian.cpp above
                    // function applyRobinBCTermsCoeffs.
                    auto const& rbc = (*this->m_robin_bcval[amrlev])[mfi].const_array(icomp*3);
                    if (this->m_lobc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_lo)
                    {
                        if (idim == 0) {
                            RT fac = beta*dxi*dxi;
                            AMREX_HOST_DEVICE_FOR_3D(blo, i, j, k,
                            {
                                RT A = rbc(i,j,k,2)
                                    / (rbc(i,j,k,1)*dxi + rbc(i,j,k,0)*RT(0.5));
                                rhsfab(i+1,j,k,icomp) += fac*bfab(i+1,j,k,icomp)*A;
                            });
                        } else if (idim == 1) {
                            RT fac = beta*dyi*dyi;
                            AMREX_HOST_DEVICE_FOR_3D(blo, i, j, k,
                            {
                                RT A = rbc(i,j,k,2)
                                    / (rbc(i,j,k,1)*dyi + rbc(i,j,k,0)*RT(0.5));
                                rhsfab(i,j+1,k,icomp) += fac*bfab(i,j+1,k,icomp)*A;
                            });
                        } else {
                            RT fac = beta*dzi*dzi;
                            AMREX_HOST_DEVICE_FOR_3D(blo, i, j, k,
                            {
                                RT A = rbc(i,j,k,2)
                                    / (rbc(i,j,k,1)*dzi + rbc(i,j,k,0)*RT(0.5));
                                rhsfab(i,j,k+1,icomp) += fac*bfab(i,j,k+1,icomp)*A;
                            });
                        }
                    }
                    if (this->m_hibc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_hi)
                    {
                        if (idim == 0) {
                            RT fac = beta*dxi*dxi;
                            AMREX_HOST_DEVICE_FOR_3D(bhi, i, j, k,
                            {
                                RT A = rbc(i,j,k,2)
                                    / (rbc(i,j,k,1)*dxi + rbc(i,j,k,0)*RT(0.5));
                                rhsfab(i-1,j,k,icomp) += fac*bfab(i,j,k,icomp)*A;
                            });
                        } else if (idim == 1) {
                            RT fac = beta*dyi*dyi;
                            AMREX_HOST_DEVICE_FOR_3D(bhi, i, j, k,
                            {
                                RT A = rbc(i,j,k,2)
                                    / (rbc(i,j,k,1)*dyi + rbc(i,j,k,0)*RT(0.5));
                                rhsfab(i,j-1,k,icomp) += fac*bfab(i,j,k,icomp)*A;
                            });
                        } else {
                            RT fac = beta*dzi*dzi;
                            AMREX_HOST_DEVICE_FOR_3D(bhi, i, j, k,
                            {
                                RT A = rbc(i,j,k,2)
                                    / (rbc(i,j,k,1)*dzi + rbc(i,j,k,0)*RT(0.5));
                                rhsfab(i,j,k-1,icomp) += fac*bfab(i,j,k,icomp)*A;
                            });
                        }
                    }
                }
            }
        }
 
    }
}
 
template <typename MF>
void
MLCellABecLapT<MF>::addInhomogNeumannFlux (
    int amrlev, const Array<MF*,AMREX_SPACEDIM>& grad, MF const& sol,
    bool mult_bcoef) const
{
    /*
     * if (mult_bcoef == true)
     *     grad is -bceof*grad phi
     * else
     *     grad is grad phi
     */
    RT fac = mult_bcoef ? RT(-1.0) : RT(1.0);
 
    bool has_inhomog_neumann = this->hasInhomogNeumannBC();
    bool has_robin = this->hasRobinBC();
 
    if (!has_inhomog_neumann && !has_robin) { return; }
 
    int ncomp = this->getNComp();
    const int mglev = 0;
 
    const auto dxinv = this->m_geom[amrlev][mglev].InvCellSize();
    const Box domain = this->m_geom[amrlev][mglev].growPeriodicDomain(1);
 
    Array<MF const*, AMREX_SPACEDIM> bcoef = {AMREX_D_DECL(nullptr,nullptr,nullptr)};
    if (mult_bcoef) {
        bcoef = getBCoeffs(amrlev,mglev);
    }
 
    const auto& bndry = *this->m_bndry_sol[amrlev];
 
    MFItInfo mfi_info;
    if (Gpu::notInLaunchRegion()) { mfi_info.SetDynamic(true); }
 
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
    for (MFIter mfi(sol, mfi_info); mfi.isValid(); ++mfi)
    {
        Box const& vbx = mfi.validbox();
        for (OrientationIter orit; orit.isValid(); ++orit) {
            const Orientation ori = orit();
            const int idim = ori.coordDir();
            const Box& ccb = amrex::adjCell(vbx, ori);
            const Dim3 os = IntVect::TheDimensionVector(idim).dim3();
            const RT dxi = static_cast<RT>(dxinv[idim]);
            if (! domain.contains(ccb)) {
                for (int icomp = 0; icomp < ncomp; ++icomp) {
                    auto const& phi = sol.const_array(mfi,icomp);
                    auto const& bv = bndry.bndryValues(ori).multiFab().const_array(mfi,icomp);
                    auto const& bc = bcoef[idim] ? bcoef[idim]->const_array(mfi,icomp)
                        : Array4<RT const>{};
                    auto const& f = grad[idim]->array(mfi,icomp);
                    if (ori.isLow()) {
                        if (this->m_lobc_orig[icomp][idim] ==
                            LinOpBCType::inhomogNeumann) {
                            AMREX_HOST_DEVICE_FOR_3D(ccb, i, j, k,
                            {
                                int ii = i+os.x;
                                int jj = j+os.y;
                                int kk = k+os.z;
                                RT b = bc ? bc(ii,jj,kk) : RT(1.0);
                                f(ii,jj,kk) = fac*b*bv(i,j,k);
                            });
                        } else if (this->m_lobc_orig[icomp][idim] ==
                                   LinOpBCType::Robin) {
                            auto const& rbc = (*this->m_robin_bcval[amrlev])[mfi].const_array(icomp*3);
                            AMREX_HOST_DEVICE_FOR_3D(ccb, i, j, k,
                            {
                                int ii = i+os.x;
                                int jj = j+os.y;
                                int kk = k+os.z;
                                RT tmp = RT(1.0) /
                                    (rbc(i,j,k,1)*dxi + rbc(i,j,k,0)*RT(0.5));
                                RT RA = rbc(i,j,k,2) * tmp;
                                RT RB = (rbc(i,j,k,1)*dxi - rbc(i,j,k,0)*RT(0.5)) * tmp;
                                RT b = bc ? bc(ii,jj,kk) : RT(1.0);
                                f(ii,jj,kk) = fac*b*dxi*((RT(1.0)-RB)*phi(ii,jj,kk)-RA);
                            });
                        }
                    } else {
                        if (this->m_hibc_orig[icomp][idim] ==
                            LinOpBCType::inhomogNeumann) {
                            AMREX_HOST_DEVICE_FOR_3D(ccb, i, j, k,
                            {
                                RT b = bc ? bc(i,j,k) : RT(1.0);
                                f(i,j,k) = fac*b*bv(i,j,k);
                            });
                        } else if (this->m_hibc_orig[icomp][idim] ==
                                   LinOpBCType::Robin) {
                            auto const& rbc = (*this->m_robin_bcval[amrlev])[mfi].const_array(icomp*3);
                            AMREX_HOST_DEVICE_FOR_3D(ccb, i, j, k,
                            {
                                RT tmp = RT(1.0) /
                                    (rbc(i,j,k,1)*dxi + rbc(i,j,k,0)*RT(0.5));
                                RT RA = rbc(i,j,k,2) * tmp;
                                RT RB = (rbc(i,j,k,1)*dxi - rbc(i,j,k,0)*RT(0.5)) * tmp;
                                RT b = bc ? bc(i,j,k) : RT(1.0);
                                f(i,j,k) = fac*b*dxi*(RA+(RB-RT(1.0))*
                                                      phi(i-os.x,j-os.y,k-os.z));
                            });
                        }
                    }
                }
            }
        }
    }
}
 
template <typename MF>
void
MLCellABecLapT<MF>::applyOverset (int amrlev, MF& rhs) const
{
    if (m_overset_mask[amrlev][0]) {
        const int ncomp = this->getNComp();
#ifdef AMREX_USE_GPU
        if (Gpu::inLaunchRegion() && m_overset_mask[amrlev][0]->isFusingCandidate()) {
            auto const& osma = m_overset_mask[amrlev][0]->const_arrays();
            auto const& rhsa = rhs.arrays();
            ParallelFor(*m_overset_mask[amrlev][0], IntVect(0), ncomp,
            [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) noexcept
            {
                if (osma[box_no](i,j,k) == 0) {
                    rhsa[box_no](i,j,k,n) = RT(0.0);
                }
            });
            if (!Gpu::inNoSyncRegion()) {
                Gpu::streamSynchronize();
            }
        } else
#endif
        {
#ifdef AMREX_USE_OMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
            for (MFIter mfi(*m_overset_mask[amrlev][0],TilingIfNotGPU()); mfi.isValid(); ++mfi)
            {
                const Box& bx = mfi.tilebox();
                auto const& rfab = rhs.array(mfi);
                auto const& osm = m_overset_mask[amrlev][0]->const_array(mfi);
                AMREX_HOST_DEVICE_PARALLEL_FOR_4D(bx, ncomp, i, j, k, n,
                {
                    if (osm(i,j,k) == 0) { rfab(i,j,k,n) = RT(0.0); }
                });
            }
        }
    }
}
 
template <typename MF>
Array<MF,AMREX_SPACEDIM>
MLCellABecLapT<MF>::makeDefaultBCoeffs (int mglev) const
{
    const int amrlev = 0;
    Array<MF,AMREX_SPACEDIM> beta;
    for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
        beta[idim].define(amrex::convert(this->m_grids[amrlev][mglev],
                                         IntVect::TheDimensionVector(idim)),
                          this->m_dmap[amrlev][mglev], 1, 0,
                          MFInfo().SetArena(The_Async_Arena()),
                          *(this->m_factory[amrlev][mglev]));
        beta[idim].setVal(RT(1.0));
        this->applyMetricTerm(amrlev, mglev, beta[idim]);
    }
    return beta;
}
 
template <typename MF>
std::unique_ptr<MLAlgMG>
MLCellABecLapT<MF>::makeAlgMG (int mglev) const
{
    if constexpr (!std::is_same<MF,MultiFab>()) {
        amrex::ignore_unused(mglev);
        amrex::Abort("MLCellABecLap: AlgMG only supports MultiFab");
        return nullptr;
    } else {
        const int amrlev = 0;
        const BoxArray& ba = this->m_grids[amrlev][mglev];
        const DistributionMapping& dm = this->m_dmap[amrlev][mglev];
        const Geometry& geom = this->m_geom[amrlev][mglev];
        const auto& factory = *(this->m_factory[amrlev][mglev]);
 
        auto const& bcondloc = *(this->m_bcondloc[amrlev][mglev]);
        LayoutData<GpuArray<int,2*AMREX_SPACEDIM>> bctype(ba, dm);
        LayoutData<GpuArray<Real,2*AMREX_SPACEDIM>> bcl(ba, dm);
        for (MFIter mfi(bctype); mfi.isValid(); ++mfi) {
            for (OrientationIter oit; oit; oit++) {
                int const cdir = oit();
                bctype[mfi][cdir] = bcondloc.bndryConds(mfi, 0)[cdir];
                bcl[mfi][cdir] = bcondloc.bndryLocs(mfi, 0)[cdir];
            }
        }
 
        auto ac = getACoeffs(amrlev, mglev);
        auto bc = getBCoeffs(amrlev, mglev);
        Array<MF,AMREX_SPACEDIM> beta;
        if (bc[0] == nullptr) {
            beta = makeDefaultBCoeffs(mglev);
            bc = amrex::GetArrOfConstPtrs(beta);
        }
        return std::make_unique<MLAlgMG>(mglev, ba, dm, geom, factory,
                                         getOversetMask(amrlev, mglev),
                                         getAScalar(), getBScalar(), ac, bc,
                                         getEBBCoeffs(amrlev, mglev), bctype, bcl,
                                         this->getMaxOrder());
    }
}
 
#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
template <typename MF>
std::unique_ptr<Hypre>
MLCellABecLapT<MF>::makeHypre (Hypre::Interface hypre_interface) const
{
    if constexpr (!std::is_same<MF,MultiFab>()) {
        amrex::Abort("MLCellABecLap Hypre interface only supports MultiFab");
    } else {
        const BoxArray& ba = this->m_grids[0].back();
        const DistributionMapping& dm = this->m_dmap[0].back();
        const Geometry& geom = this->m_geom[0].back();
        const auto& factory = *(this->m_factory[0].back());
        MPI_Comm comm = this->BottomCommunicator();
 
        const int mglev = this->NMGLevels(0)-1;
 
        auto om = getOversetMask(0, mglev);
 
        auto hypre_solver = amrex::makeHypre(ba, dm, geom, comm, hypre_interface, om);
 
        hypre_solver->setScalars(getAScalar(), getBScalar());
 
        auto ac = getACoeffs(0, mglev);
        if (ac)
        {
            hypre_solver->setACoeffs(*ac);
        }
        else
        {
            MultiFab alpha(ba,dm,1,0,MFInfo().SetArena(The_Async_Arena()),factory);
            alpha.setVal(0.0);
            hypre_solver->setACoeffs(alpha);
        }
 
        auto bc = getBCoeffs(0, mglev);
        if (bc[0])
        {
            hypre_solver->setBCoeffs(bc);
        }
        else
        {
            auto beta = makeDefaultBCoeffs(mglev);
            hypre_solver->setBCoeffs(amrex::GetArrOfConstPtrs(beta));
        }
        hypre_solver->setIsMatrixSingular(this->isBottomSingular());
 
        return hypre_solver;
    }
    return nullptr;
}
#endif
 
#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
template <typename MF>
std::unique_ptr<PETScABecLap>
MLCellABecLapT<MF>::makePETSc () const
{
    if constexpr (!std::is_same<MF,MultiFab>()) {
        amrex::Abort("MLCellABecLap PETSc interface only supports MultiFab");
    } else {
        const BoxArray& ba = this->m_grids[0].back();
        const DistributionMapping& dm = this->m_dmap[0].back();
        const Geometry& geom = this->m_geom[0].back();
        const auto& factory = *(this->m_factory[0].back());
        MPI_Comm comm = this->BottomCommunicator();
 
        auto petsc_solver = makePetsc(ba, dm, geom, comm);
 
        petsc_solver->setScalars(getAScalar(), getBScalar());
 
        const int mglev = this->NMGLevels(0)-1;
        auto ac = getACoeffs(0, mglev);
        if (ac)
        {
            petsc_solver->setACoeffs(*ac);
        }
        else
        {
            MultiFab alpha(ba,dm,1,0,MFInfo().SetArena(The_Async_Arena()),factory);
            alpha.setVal(0.0);
            petsc_solver->setACoeffs(alpha);
        }
 
        auto bc = getBCoeffs(0, mglev);
        if (bc[0])
        {
            petsc_solver->setBCoeffs(bc);
        }
        else
        {
            auto beta = makeDefaultBCoeffs(mglev);
            petsc_solver->setBCoeffs(amrex::GetArrOfConstPtrs(beta));
        }
        return petsc_solver;
    }
    return nullptr;
}
#endif
 
extern template class MLCellABecLapT<MultiFab>;
 
using MLCellABecLap = MLCellABecLapT<MultiFab>;
 
}
 
#endif