AMReX_MLNodeLaplacian.H

Go to the documentation of this file.

#ifndef AMREX_ML_NODE_LAPLACIAN_H_
#define AMREX_ML_NODE_LAPLACIAN_H_
#include <AMReX_Config.H>
 
#include <AMReX_iMultiFab.H>
#include <AMReX_MLNodeLinOp.H>
#include <AMReX_MultiFab.H>
#include <AMReX_Vector.H>
 
// Feature macro: when defined, MLNodeLaplacian::updateVelocity and
// MLNodeLaplacian::getFluxes correctly route through the anisotropic
// (harmonic-average) mknewu kernel when setMapped(true) was used.  Prior
// to this, both functions only used sigma component 0 (sigma_x) for all
// velocity/flux components, producing the wrong correction for
// non-identity mesh mappings.  EB and RAP coarsening are not covered by
// this feature (see comments in AMReX_MLNodeLaplacian_misc.cpp).
#define AMREX_MLNODELAP_HAS_MKNEWU_HA 1
 
namespace amrex {

 
// del dot (sigma grad phi) = rhs
// where phi and rhs are nodal, and sigma is cell-centered.

class MLNodeLaplacian
    : public MLNodeLinOp
{
public :
 
    friend class HypreMLNodeLap;

    MLNodeLaplacian () = default;
    MLNodeLaplacian (const Vector<Geometry>& a_geom,
                     const Vector<BoxArray>& a_grids,
                     const Vector<DistributionMapping>& a_dmap,
                     const LPInfo& a_info = LPInfo(),
                     const Vector<FabFactory<FArrayBox> const*>& a_factory = {},
                     Real  a_const_sigma = Real(0.0));
#ifdef AMREX_USE_EB
    MLNodeLaplacian (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,
                     Real a_const_sigma = Real(0.0));
#endif
    ~MLNodeLaplacian () override = default;
 
    MLNodeLaplacian (const MLNodeLaplacian&) = delete;
    MLNodeLaplacian (MLNodeLaplacian&&) = delete;
    MLNodeLaplacian& operator= (const MLNodeLaplacian&) = delete;
    MLNodeLaplacian& operator= (MLNodeLaplacian&&) = 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<FArrayBox> const*>& a_factory = {},
                 Real  a_const_sigma = Real(0.0));
 
#ifdef AMREX_USE_EB
    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,
                 Real a_const_sigma = Real(0.0));
#endif
 
    std::string name () const override { return std::string("MLNodeLaplacian"); }

    void setRZCorrection (bool rz) noexcept { m_is_rz = rz; }

    void setNormalizationThreshold (Real t) noexcept {
        if (m_normalization_threshold != t) {
            m_normalization_threshold = t;
            m_needs_update = true;
        }
    }

    void setSigma (int amrlev, const MultiFab& a_sigma);

    void compDivergence (const Vector<MultiFab*>& rhs, const Vector<MultiFab*>& vel);

    void compRHS (const Vector<MultiFab*>& rhs, const Vector<MultiFab*>& vel,
                  const Vector<const MultiFab*>& rhnd,
                  const Vector<MultiFab*>& rhcc);

    void updateVelocity (const Vector<MultiFab*>& vel, const Vector<MultiFab const*>& sol) const;

    void compSyncResidualCoarse (MultiFab& sync_resid, const MultiFab& phi,
                                 const MultiFab& vold, const MultiFab* rhcc,
                                 const BoxArray& fine_grids, const IntVect& ref_ratio);

    void compSyncResidualFine (MultiFab& sync_resid, const MultiFab& phi, const MultiFab& vold,
                               const MultiFab* rhcc);

    void setGaussSeidel (bool flag) noexcept { m_use_gauss_seidel = flag; }
    void setHarmonicAverage (bool flag) noexcept {
        if (m_use_harmonic_average != flag) {
            m_use_harmonic_average = flag;
            m_needs_update = true;
        }
    }

    void setMapped (bool flag) noexcept {
        if (m_use_mapped != flag) {
            m_use_mapped = flag;
            m_needs_update = true;
        }
    }

    void setCoarseningStrategy (CoarseningStrategy cs) noexcept {
        if (m_const_sigma == Real(0.0) && m_coarsening_strategy != cs) {
            m_coarsening_strategy = cs;
            m_masks_built = false;
            m_precond_weight_mask.clear();
            m_needs_update = true;
        }
    }
 
    BottomSolver getDefaultBottomSolver () const final {
        return (m_coarsening_strategy == CoarseningStrategy::RAP) ?
            BottomSolver::bicgcg : BottomSolver::bicgstab;
    }

    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;
    void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs,
                                         const MultiFab& fine_sol, const MultiFab& fine_rhs) final;

    void interpCoarseFineSolution (int famrlev, MultiFab& fine_sol,
                                   const MultiFab& crse_sol) const final;

    void reflux (int crse_amrlev,
                         MultiFab& res, const MultiFab& crse_sol, const MultiFab& crse_rhs,
                         MultiFab& fine_res, MultiFab& fine_sol, const MultiFab& fine_rhs) const final;

    void prepareForSolve () final;
    [[nodiscard]] bool needsUpdate () const final { return m_needs_update; }
    void update () 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;

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

    void getFluxes (const Vector<MultiFab*>& a_flux,
                            const Vector<MultiFab*>& a_sol) const final;
    void unimposeNeumannBC (int amrlev, MultiFab& rhs) const final;
    Vector<Real> getSolvabilityOffset (int amrlev, int mglev,
                                       MultiFab const& rhs) const override;
    void fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs,
                                 Vector<Real> const& offset) const override;

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

    void compGrad (int amrlev, MultiFab& grad, MultiFab& sol) const;

    void averageDownCoeffs ();
    void averageDownCoeffsToCoarseAmrLevel (int flev);
    void averageDownCoeffsSameAmrLevel (int amrlev);

    void restrictInteriorNodes (int camrlev, MultiFab& crhs, MultiFab& frhs) const;

    void FillBoundaryCoeff (MultiFab& sigma, const Geometry& geom);

    void buildStencil ();
 
#ifdef AMREX_USE_EB
    void buildIntegral ();
    void buildSurfaceIntegral ();
 
    // Tells the solver that EB boundaries have inflow specified by "eb_vel"
    void setEBInflowVelocity (int amrlev, const MultiFab& eb_vel);
 
#endif
 
#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 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. AlgGid is the global row id type: Long for AlgMG,
    // HypreNodeLap::AtomicInt for hypre.
    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;
#ifdef AMREX_USE_GPU
    template <typename AlgInt, typename AlgGid>
    void fillMatrix_gpu (int mglev, MFIter const& mfi,
                         Array4<AlgGid const> const& gid,
                         Array4<int const> const& lid,
                         AlgInt* ncols, AlgInt* cols, Real* mat) const;
#endif
    template <typename AlgInt, typename AlgGid>
    void fillMatrix_cpu (int mglev, MFIter const& mfi,
                         Array4<AlgGid const> const& gid,
                         Array4<int const> const& lid,
                         AlgInt* ncols, AlgInt* cols, Real* mat) const;
#endif
 
protected:
 
    void resizeMultiGrid (int new_size) final;
 
private:
 
    int m_is_rz = 0;
 
    Real m_const_sigma = Real(0.0);
    Vector<Vector<Array<std::unique_ptr<MultiFab>,AMREX_SPACEDIM> > > m_sigma;
    Vector<Vector<std::unique_ptr<MultiFab> > > m_stencil;
    Vector<std::unique_ptr<MultiFab> > m_nosigma_stencil;
    Vector<Vector<Real> > m_s0_norm0;
 
    Real m_normalization_threshold = Real(1.e-8);
 
#ifdef AMREX_USE_EB
    // they could be MultiCutFab
    Vector<std::unique_ptr<MultiFab> > m_integral;
    bool m_integral_built = false;
 
    Vector<std::unique_ptr<MultiFab> > m_surface_integral;
    bool m_surface_integral_built = false;
    bool m_build_surface_integral = false;
    Vector<std::unique_ptr<MultiFab> > m_eb_vel_dot_n;
#endif
 
    bool m_use_gauss_seidel     = true;
    bool m_use_harmonic_average = false;
    bool m_use_mapped           = false;
    bool m_needs_update         = true;
 
    void checkPoint (std::string const& file_name) const final;
};
 
}
 
#endif