AMReX_Hypre.H

Go to the documentation of this file.

#ifndef AMREX_HYPRE_H_
#define AMREX_HYPRE_H_
#include <AMReX_Config.H>
 
#include <AMReX_Geometry.H>
#include <AMReX_MultiFab.H>
#include <AMReX_BndryData.H>
 
#include "HYPRE.h"
#include "_hypre_utilities.h"
#include "HYPRE_parcsr_ls.h"
 
#include <memory>
 
namespace amrex
{
 
namespace HypreDefaults {
#ifdef HYPRE_USING_GPU
    inline constexpr bool gpu = true;
#else
    inline constexpr bool gpu = false;
#endif
    inline constexpr bool old_default = !gpu;
    inline constexpr int coarsen_type = gpu ? 8 : 6;
    inline constexpr int interp_type = gpu ? 6 : 0;
    inline constexpr int pmax_elmts = 4;
    inline constexpr int keep_transpose = gpu ? 1 : 0;
    inline constexpr int relax_type = gpu ? 18 : 6;
    inline constexpr int relax_order = gpu ? 0 : 1;
 
    inline void setGpuOptions (HYPRE_Solver solver)
    {
        HYPRE_BoomerAMGSetCoarsenType(solver, coarsen_type);
        HYPRE_BoomerAMGSetInterpType(solver, interp_type);
        HYPRE_BoomerAMGSetPMaxElmts(solver, pmax_elmts);
        HYPRE_BoomerAMGSetKeepTranspose(solver, keep_transpose);
    }
}
 
class Hypre
{
public:
 
    enum class Interface : int { structed, semi_structed, ij };
 
    Hypre (const BoxArray& grids,
           const DistributionMapping& dmap,
           const Geometry& geom,
           MPI_Comm comm_);
 
    virtual ~Hypre ();
    Hypre (Hypre &&) noexcept = default;
    Hypre& operator= (Hypre &&) noexcept = default;
 
    Hypre (Hypre const&) = delete;
    Hypre& operator= (Hypre const&) = delete;
 
    void setScalars (Real sa, Real sb);
    void setACoeffs (const MultiFab& alpha);
    void setBCoeffs (const Array<const MultiFab*,BL_SPACEDIM>& beta);
    void setVerbose (int _verbose);
    void setIsMatrixSingular(bool flag) { is_matrix_singular = flag; }
    virtual void solve (MultiFab& soln, const MultiFab& rhs, Real rel_tol, Real abs_tol,
                        int max_iter, const BndryData& bndry, int max_bndry_order) = 0;
 
    static HYPRE_Int ispow2 (HYPRE_Int i)
    {
        if (i == 1) {
            return 1;
        } else if ((i <= 0) || (i & 1)) {
            return 0;
        } else {
            return ispow2(i/2);
        }
    }
 
    static Array<HYPRE_Int,AMREX_SPACEDIM> loV (const Box& b) {
        const auto& v = b.loVect();
        return {AMREX_D_DECL(static_cast<HYPRE_Int>(v[0]),
                             static_cast<HYPRE_Int>(v[1]),
                             static_cast<HYPRE_Int>(v[2]))};
    }
 
    static Array<HYPRE_Int,AMREX_SPACEDIM> hiV (const Box& b) {
        const auto& v = b.hiVect();
        return {AMREX_D_DECL(static_cast<HYPRE_Int>(v[0]),
                             static_cast<HYPRE_Int>(v[1]),
                             static_cast<HYPRE_Int>(v[2]))};
    }
 
 
    void setHypreOptionsNamespace (const std::string& ns) noexcept
    { options_namespace = ns; }
    void setHypreOldDefault (bool l) noexcept {old_default = l;}
    void setHypreRelaxType (int n) noexcept {relax_type = n;}
    void setHypreRelaxOrder (int n) noexcept {relax_order = n;}
    void setHypreNumSweeps (int n) noexcept {num_sweeps = n;}
    void setHypreStrongThreshold (Real t) noexcept {strong_threshold = t;}
 
    static constexpr HYPRE_Int regular_stencil_size = 2*AMREX_SPACEDIM + 1;
    static constexpr HYPRE_Int eb_stencil_size = AMREX_D_TERM(3, *3, *3);
 
protected:
 
    MPI_Comm comm = MPI_COMM_NULL;
    Geometry geom;
 
    int verbose = 0;
    bool old_default = HypreDefaults::old_default;
    int relax_type = HypreDefaults::relax_type;
    int relax_order = HypreDefaults::relax_order;
    int num_sweeps = 2;  // Sweeeps on each level
    Real strong_threshold = Real(0.25); // HYPRE default is 0.25
 
    std::string options_namespace{"hypre"};
 
    MultiFab acoefs;
    Array<MultiFab,AMREX_SPACEDIM> bcoefs;
    Real scalar_a, scalar_b;
 
    FabFactory<FArrayBox> const* m_factory = nullptr;
    BndryData const* m_bndry = nullptr;
    int m_maxorder = -1;
 
    bool is_matrix_singular { false };
};
 
[[nodiscard]] std::unique_ptr<Hypre>
makeHypre (const BoxArray& grids, const DistributionMapping& dmap,
           const Geometry& geom, MPI_Comm comm_, Hypre::Interface interface,
           const iMultiFab* overset_mask = nullptr);
 
}
 
#endif