AMReX_HypreMLNodeLap.H

Go to the documentation of this file.

#ifndef AMREX_HYPRE_ML_NODELAP_H_
#define AMREX_HYPRE_ML_NODELAP_H_
#include <AMReX_Config.H>
 
#include <AMReX_CompositeSystemUtil.H>
#include <AMReX_HypreIJIface.H>
#include <AMReX_MultiFab.H>
 
#include <memory>
#include <string>

 
namespace amrex {
 
class MLNodeLaplacian;

class HypreMLNodeLap
{
public:
    using BigInt = HYPRE_BigInt;

    HypreMLNodeLap (MLNodeLaplacian const& a_linop,
                    std::string a_prefix = "hypre_mlnodelap",
                    int a_verbose = 0);
 
    ~HypreMLNodeLap ();
 
    HypreMLNodeLap (HypreMLNodeLap const&) = delete;
    HypreMLNodeLap (HypreMLNodeLap &&) = delete;
    HypreMLNodeLap& operator= (HypreMLNodeLap const&) = delete;
    HypreMLNodeLap& operator= (HypreMLNodeLap &&) = delete;

    void solve (Vector<MultiFab*> const& a_sol, Vector<MultiFab const*> const& a_rhs,
                Real a_reltol, Real a_abstol, int a_maxiter);

    void setThrowException (bool t) noexcept { m_throw_exception = t; }

    [[nodiscard]] int getNumIters () const { return m_num_iters; }
    [[nodiscard]] Real getFinalResidualNorm () const { return m_final_res; }

    [[nodiscard]] BigInt numUnknowns () const { return m_nglobal; }

    void buildMatrix ();
 
private:
    using IdFab = BaseFab<Long>;
    using IdFabArray = FabArray<IdFab>;
 
    void numberUnknowns ();
    void loadRHS (Vector<MultiFab const*> const& a_rhs, Vector<MultiFab*> const& a_sol);
    void getSolution (Vector<MultiFab*> const& a_sol);
    void dumpSystem (int k, bool solution);
    void fail (std::string const& msg) const;

    [[nodiscard]] int neumannScale (int lev, IntVect const& iv) const;
    [[nodiscard]] bool onDirichletFace (int lev, IntVect const& iv) const;
 
    [[nodiscard]] FabFactory<FArrayBox> const* factory (int lev) const;
 
    MLNodeLaplacian const* m_linop;
    std::string m_prefix;
    int m_verbose = 0;
    int m_nlevels = 0;
    MPI_Comm m_comm = MPI_COMM_NULL;
    int m_nprocs = 1;
    int m_myproc = 0;
 
    Vector<Geometry> m_geom;
    Vector<BoxArray> m_nd_grids;  // nodal BoxArrays
    Vector<DistributionMapping> m_dmap;
    Vector<int> m_ref_ratio;      // m_ref_ratio[lev]: ratio between lev and lev+1
 
    // Host-accessible (pinned in GPU builds) per-level data
    Vector<std::unique_ptr<IdFabArray>> m_gid;         // global id or -1
    Vector<std::unique_ptr<iMultiFab>> m_owner;        // node ownership
    Vector<std::unique_ptr<iMultiFab>> m_dmsk;         // Dirichlet (incl. c/f on fine levels)
    Vector<std::unique_ptr<iMultiFab>> m_ndfine;       // nodal fine mask (numbering only)
    Vector<std::unique_ptr<IdFabArray>> m_cgid_on_f;   // level lev-1 gid on coarsened level lev nodes
    Vector<std::unique_ptr<iMultiFab>> m_cndfine_on_f; // level lev-1 nodal fine mask, same layout
    Vector<std::unique_ptr<MultiFab>> m_ndvol;         // EB: nodal fluid volume fraction (or null)
    Vector<std::unique_ptr<MultiFab>> m_scale;         // EB: unknown scaling 1/sqrt(ndvol) (or null)
    Vector<std::unique_ptr<MultiFab>> m_cscale_on_f;   // level lev-1 scaling on coarsened level lev nodes
    bool m_eb_scaling = true;     // <prefix>.eb_scaling
    bool m_has_dirichlet = false; // any non-periodic Dirichlet domain face
    bool m_throw_exception = false;

    struct DirichletCoupling { Long row; Real coef; int lev; int li; IntVect node; };
    Vector<DirichletCoupling> m_dbc;
 
    Vector<Long> m_rank_offset;   // first gid of each rank, plus total at the end
    BigInt m_ilower = 0;
    BigInt m_iupper = -1;
    BigInt m_nglobal = 0;
    Long m_pin_gid = -1;          // pinned unknown for singular problems
 
    enum class Backend { hypre, algmg };
    Backend m_backend = Backend::hypre;
 
    // Matrix Market dumps
    int m_dump_max = 0;
    std::string m_dump_prefix;
    Long m_dump_gather_max = 50000000;
 
    // Owned rows of the assembled matrix (AlgMG backend and dumps)
    composite::OwnedCSR m_csr;
    std::unique_ptr<SpMatrix<Real>> m_spA;
    std::unique_ptr<AlgMG<Real>> m_algmg;
 
    // Right-hand side and solution of this rank's rows
    Vector<Real> m_bvec;
    Vector<Real> m_xvec;
 
    std::unique_ptr<HypreIJIface> m_ij;
    HYPRE_IJMatrix m_A = nullptr;
    HYPRE_IJVector m_b = nullptr;
    HYPRE_IJVector m_x = nullptr;
 
    int m_num_iters = 0;
    Real m_final_res = 0.0;
    Real m_bnorm = 0.0;           // 2-norm of the assembled right-hand side
};
 
}
 
#endif