AMReX_HypreMLABecLap.H

Go to the documentation of this file.

#ifndef AMREX_HYPRE_ML_ABECLAP_H_
#define AMREX_HYPRE_ML_ABECLAP_H_
#include <AMReX_Config.H>
 
#include <AMReX_CompositeSystemUtil.H>
#include <AMReX_Geometry.H>
#include <AMReX_HypreIJIface.H>
#include <AMReX_LO_BCTYPES.H>
#include <AMReX_MLMGBndry.H>
#include <AMReX_MultiFab.H>
#include <AMReX_MultiMask.H>
 
#include <HYPRE_sstruct_ls.h>
 
#include <limits>
#include <unordered_map>
#include <utility>

 
namespace amrex {
 
enum struct HypreSolverID {
    BoomerAMG, SSAMG
};
 
// single component only, cell centered only

class HypreMLABecLap
{
public:
    HypreMLABecLap (HypreMLABecLap const&) = delete;
    HypreMLABecLap (HypreMLABecLap &&) = delete;
    HypreMLABecLap& operator= (HypreMLABecLap const&) = delete;
    HypreMLABecLap& operator= (HypreMLABecLap &&) = delete;

    HypreMLABecLap (Vector<Geometry> a_geom,
                    Vector<BoxArray> a_grids,
                    Vector<DistributionMapping> a_dmap,
                    HypreSolverID a_hypre_solver_id,
                    std::string a_parmparse_prefix = "hypre_mlabeclap");
 
    ~HypreMLABecLap ();

    void setVerbose (int v) { m_verbose = v; }
    void setMaxIter (int v) { m_maxiter = v; }
    void setIsSingular (bool v) { m_is_singular = v; }

    void setReuseSolverSetup (bool v) { m_reuse_setup = v; }
    void setReuseSetupMaxIter (int n) noexcept { m_reuse_max_iter = n; }
    void resetSolverSetup () {
        if (m_ij_solver) { m_ij_solver->setNeedSetup(); }
        m_algmg_reset = true;
    }

    void setThrowException (bool t) noexcept { m_throw_exception = t; }
    [[nodiscard]] int getNumIters () const noexcept { return m_num_iters; }
    [[nodiscard]] Real getFinalResidualNorm () const noexcept { return m_final_res; }

    void setup (Real a_ascalar, Real a_bscalar,
                Vector<MultiFab const*> const& a_acoefs,
                Vector<Array<MultiFab const*,AMREX_SPACEDIM>> const& a_bcoefs,
                Array<LinOpBCType,AMREX_SPACEDIM> const& a_lobc,
                Array<LinOpBCType,AMREX_SPACEDIM> const& a_hibc,
                Vector<MultiFab const*> const& a_levelbcdata,
                std::pair<MultiFab const*, IntVect> const& a_coarse_bc = {nullptr, IntVect(0)});

    void solve (Vector<MultiFab*> const& a_sol, Vector<MultiFab const*> const& a_rhs,
                Real a_reltol, Real a_abstol);
 
    // update? updateDirichleBC? or updateCoeffs?
 
    // public for cuda

    void commBCoefs (int flev, Array<MultiFab const*,AMREX_SPACEDIM> const& a_bcoefs);
    void commBCoefs_local (int flev, Array<MultiFab const*,AMREX_SPACEDIM> const& a_bcoefs,
                           Vector<FabArrayBase::CopyComTag> const& tags);

    void buildIJMaps ();

    void buildIJSlots (Gpu::DeviceVector<HYPRE_BigInt> const& rows,
                       Gpu::DeviceVector<HYPRE_BigInt> const& cols);
    void updateIJValues (Gpu::DeviceVector<Real> const& vals);

    void solveIJ (Vector<MultiFab*> const& a_sol, Vector<MultiFab const*> const& a_rhs,
                  Real a_reltol, Real a_abstol);
 
private:
 
    void addNonStencilEntriesToGraph ();

    [[nodiscard]] Long cellId (int lev, IntVect const& iv) const;
 
    int m_verbose = 0;
    int m_maxiter = 200;
    bool m_is_singular = false;
 
    Vector<Geometry> m_geom;
    Vector<BoxArray> m_grids;
    Vector<DistributionMapping> m_dmap;
    std::string m_parmparse_prefix;
    int m_nlevels = 0;
    MPI_Comm m_comm = MPI_COMM_NULL;
 
    Vector<IntVect> m_ref_ratio;
 
    Real m_ascalar = std::numeric_limits<Real>::max();
    Real m_bscalar = std::numeric_limits<Real>::max();
    Array<LinOpBCType,AMREX_SPACEDIM> m_lobc{AMREX_D_DECL(LinOpBCType::bogus,
                                                          LinOpBCType::bogus,
                                                          LinOpBCType::bogus)};
    Array<LinOpBCType,AMREX_SPACEDIM> m_hibc{AMREX_D_DECL(LinOpBCType::bogus,
                                                          LinOpBCType::bogus,
                                                          LinOpBCType::bogus)};
 
    Vector<std::unique_ptr<MLMGBndry>> m_bndry;
    Vector<std::unique_ptr<BndryRegister>> m_bndry_rhs;
    Vector<iMultiFab> m_fine_masks;
    Vector<iMultiFab> m_crse_masks;
 
    // For coarse cells at coarse/fine interface. The vector is for AMR
    // levels.
    Vector<iMultiFab> m_c2f_offset_from; // offset for sparse coarse from-cells
    Vector<LayoutData<int>> m_c2f_total_from; // # of coarse from-cells w/ c2f entries
    Vector<iMultiFab> m_c2f_nentries; // # of non-stencil entries
    Vector<iMultiFab> m_c2f_offset_to; // offset for sparse to-cells, including fine (and coarse in 3d) cells
    Vector<LayoutData<int>> m_c2f_total_to; // total sum of non-stencil entries in a Box
 
    // B coefficients at coarse/fine interface
    Vector<Array<iMultiFab,AMREX_SPACEDIM>> m_offset_cf_bcoefs;
    Vector<Array<LayoutData<std::unique_ptr<Gpu::DeviceVector<Real>>>,AMREX_SPACEDIM>> m_cf_bcoefs;
 
#ifdef AMREX_USE_GPU
    template <class T> using HostVector = Gpu::PinnedVector<T>;
#else
    template <class T> using HostVector = Vector<T>;
#endif
 
    // For fine cells at coarse/fine interface. The non-stencil entries are
    // from fine to coarse. The outer vector is for AMR levels.
    Vector<HostVector<int>> m_f2c_bno; // local box number
    Vector<HostVector<IntVect>> m_f2c_cell; // fine cell
    Vector<Vector<HYPRE_Int>> m_f2c_nentries; // # of non-stencil entries
    Vector<HostVector<std::size_t>> m_f2c_offset; // offset into m_f2c_values
    Vector<HostVector<Real>> m_f2c_values; // values for non-stencil entries
 
    HYPRE_SStructGrid    m_ss_grid = nullptr;
    HYPRE_SStructStencil m_ss_stencil = nullptr;
    HYPRE_SStructGraph   m_ss_graph = nullptr;
    HYPRE_SStructSolver  m_ss_solver = nullptr;
    HYPRE_SStructSolver  m_ss_precond = nullptr;
    HYPRE_SStructMatrix  m_ss_A = nullptr;
    HYPRE_SStructVector  m_ss_x = nullptr;
    HYPRE_SStructVector  m_ss_b = nullptr;
 
    HYPRE_Solver m_solver = nullptr;

    std::unique_ptr<HypreIJIface> m_ij_solver;

    bool m_use_ij = false;
    bool m_reuse_setup = false;
    bool m_throw_exception = false;
    bool m_ij_vectors_assembled = false;
    bool m_ij_gather_failed = false; // hypre dropped an entry: no in-place updates
    bool m_algmg_reset = false;
    bool m_algmg_kept = false;
    int m_reuse_max_iter = 0;
    Vector<Vector<Long>> m_box_offset;  // first row id of each box, per level
    Long m_ij_lower = 0;
    Long m_ij_upper = -1;
    std::unordered_map<Long,Vector<Long>> m_ns_targets;
    Vector<std::unique_ptr<FabArray<BaseFab<Long>>>> m_cell_id;
    Vector<LayoutData<Gpu::DeviceVector<HYPRE_BigInt>>> m_c2f_targets;
    Vector<Gpu::DeviceVector<HYPRE_BigInt>> m_f2c_rows;
    Vector<Gpu::DeviceVector<HYPRE_BigInt>> m_f2c_cols;
    Gpu::DeviceVector<HYPRE_BigInt> m_ij_rowids; // the local rows
    struct EntryGather
    {
        Gpu::DeviceVector<Long> offset;
        Gpu::DeviceVector<Long> trip;
        [[nodiscard]] bool empty () const { return offset.empty(); }
        [[nodiscard]] Long numTriplets () const { return Long(trip.size()); }
        void clear () { offset.clear(); trip.clear(); }
    };
    EntryGather m_ij_gather;

    enum class Backend { hypre, algmg };
    Backend m_backend = Backend::hypre;
    int m_num_iters = 0;                     // last solve (IJ mode)
    Real m_final_res = 0.0;
    Vector<Long> m_rank_offset;              // first row of each rank, plus the total
    composite::OwnedCSR m_csr;               // owned rows (AlgMG backend, dumps)
    std::unique_ptr<SpMatrix<Real>> m_spA;
    std::unique_ptr<AlgMG<Real>> m_algmg;
    Gpu::DeviceVector<Long> m_alg_row_offset;
    Gpu::DeviceVector<Long> m_alg_col_index;
    EntryGather m_alg_gather;

    int m_dump_max = 0;
    std::string m_dump_prefix;
    Long m_dump_gather_max = 50000000;
    void dumpSystem (int k, Gpu::DeviceVector<Real> const& v, bool solution);
    void fail (std::string const& msg) const;
 
    HypreSolverID m_hypre_solver_id = HypreSolverID::BoomerAMG;
    HYPRE_Int m_hypre_object_type = HYPRE_PARCSR;
};
 
}
 
#endif