Hypre IJ interface: Enable access to additional solvers and preconditioners available in Hypre

Src/Extern/HYPRE/AMReX_Hypre.H

@@ -51,7 +51,10 @@ public:
                              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;}
@@ -73,6 +76,8 @@ protected:
     int num_sweeps = 2;  // Sweeeps on each level
     Real strong_threshold = 0.25; // Hypre default is 0.25
 
+    std::string options_namespace{"hypre"};
+
     MultiFab acoefs;
     Array<MultiFab,AMREX_SPACEDIM> bcoefs;
     Real scalar_a, scalar_b;

Src/Extern/HYPRE/AMReX_HypreABecLap3.H

@@ -5,9 +5,7 @@
 
 #include <AMReX_iMultiFab.H>
 #include <AMReX_LayoutData.H>
-
-#include "HYPRE_parcsr_ls.h"
-#include "_hypre_parcsr_mv.h"
+#include <AMReX_HypreIJIface.H>
 
 #include <algorithm>
 
@@ -31,11 +29,12 @@ public:
 #endif
 
 private :
+    std::unique_ptr<HypreIJIface> hypre_ij;
 
+    // Non-owning references to hypre matrix, rhs, and solution data
     HYPRE_IJMatrix A = NULL;
     HYPRE_IJVector b = NULL;
     HYPRE_IJVector x = NULL;
-    HYPRE_Solver solver = NULL;
 
     LayoutData<HYPRE_Int> ncells_grid;
     LayoutData<Gpu::ManagedDeviceVector<HYPRE_Int> > cell_id_vec;

Src/Extern/HYPRE/AMReX_HypreABecLap3.cpp

@@ -20,20 +20,10 @@ namespace amrex {
 HypreABecLap3::HypreABecLap3 (const BoxArray& grids, const DistributionMapping& dmap,
                               const Geometry& geom_, MPI_Comm comm_)
     : Hypre(grids, dmap, geom_, comm_)
-{
-}
-
+{}
+
 HypreABecLap3::~HypreABecLap3 ()
-{
-    HYPRE_IJMatrixDestroy(A);
-    A = NULL;
-    HYPRE_IJVectorDestroy(b);
-    b = NULL;
-    HYPRE_IJVectorDestroy(x);
-    x = NULL;
-    HYPRE_BoomerAMGDestroy(solver);
-    solver = NULL;
-}
+{}
 
 void
 HypreABecLap3::solve (MultiFab& soln, const MultiFab& rhs, Real rel_tol, Real abs_tol,
@@ -43,7 +33,7 @@ HypreABecLap3::solve (MultiFab& soln, const MultiFab& rhs, Real rel_tol, Real ab
 
     BL_PROFILE("HypreABecLap3::solve()");
 
-    if (solver == NULL || m_bndry != &bndry || m_maxorder != max_bndry_order)
+    if (!hypre_ij || m_bndry != &bndry || m_maxorder != max_bndry_order)
     {
         m_bndry = &bndry;
         m_maxorder = max_bndry_order;
@@ -54,52 +44,16 @@ HypreABecLap3::solve (MultiFab& soln, const MultiFab& rhs, Real rel_tol, Real ab
     {
         m_factory = &(rhs.Factory());
     }
-    
+
     HYPRE_IJVectorInitialize(b);
     HYPRE_IJVectorInitialize(x);
     //
     loadVectors(soln, rhs);
     //
     HYPRE_IJVectorAssemble(x);
     HYPRE_IJVectorAssemble(b);
-
-    HYPRE_ParCSRMatrix par_A = NULL;
-    HYPRE_ParVector par_b = NULL;
-    HYPRE_ParVector par_x = NULL;
-    HYPRE_IJMatrixGetObject(A, (void**)  &par_A);
-    HYPRE_IJVectorGetObject(b, (void **) &par_b);
-    HYPRE_IJVectorGetObject(x, (void **) &par_x);
-
-    HYPRE_BoomerAMGSetMinIter(solver, 1);
-    HYPRE_BoomerAMGSetMaxIter(solver, max_iter);
-    HYPRE_BoomerAMGSetTol(solver, rel_tol);
-    if (abs_tol > 0.0)
-    {
-        Real bnorm = hypre_ParVectorInnerProd(par_b, par_b);
-        bnorm = std::sqrt(bnorm);
-
-        const BoxArray& grids = rhs.boxArray();
-        Real volume = grids.numPts();
-        Real rel_tol_new = (bnorm > 0.0) ? (abs_tol / bnorm * std::sqrt(volume)) : rel_tol;
-
-        if (rel_tol_new > rel_tol) {
-            HYPRE_BoomerAMGSetTol(solver, rel_tol_new);
-        }
-    }
-
-    HYPRE_BoomerAMGSolve(solver, par_A, par_b, par_x);
 
-    if (verbose >= 2)
-    {
-        HYPRE_Int num_iterations;
-        Real res;
-        HYPRE_BoomerAMGGetNumIterations(solver, &num_iterations);
-        HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, &res);
-
-        amrex::Print() <<"\n" <<  num_iterations
-                       << " Hypre IJ BoomerAMG Iterations, Relative Residual "
-                       << res << std::endl;
-    }
+    hypre_ij->solve(rel_tol, abs_tol, max_iter);
 
     getSolution(soln);
 }
@@ -280,18 +234,13 @@ HypreABecLap3::prepareSolver ()
     HYPRE_Int ilower = proc_begin;
     HYPRE_Int iupper = proc_end-1;
 
-    //
-    HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper, &A);
-    HYPRE_IJMatrixSetObjectType(A, HYPRE_PARCSR);
-    HYPRE_IJMatrixInitialize(A);
-    //
-    HYPRE_IJVectorCreate(comm, ilower, iupper, &b);
-    HYPRE_IJVectorSetObjectType(b, HYPRE_PARCSR);
-    //
-    HYPRE_IJVectorCreate(comm, ilower, iupper, &x);
-    HYPRE_IJVectorSetObjectType(x, HYPRE_PARCSR);
-
-    // A.SetValues() & A.assemble()
+    hypre_ij.reset(new HypreIJIface(comm, ilower, iupper, verbose));
+    hypre_ij->parse_inputs(options_namespace);
+
+    // Obtain non-owning references to the matrix, rhs, and solution data
+    A = hypre_ij->A();
+    b = hypre_ij->b();
+    x = hypre_ij->x();
 
     const Real* dx = geom.CellSize();
     const int bho = (m_maxorder > 2) ? 1 : 0;
@@ -402,24 +351,6 @@ HypreABecLap3::prepareSolver ()
         }
     }
     HYPRE_IJMatrixAssemble(A);
-
-    // Create solver
-    HYPRE_BoomerAMGCreate(&solver);
-
-    if (old_default) HYPRE_BoomerAMGSetOldDefault(solver); // Falgout coarsening with modified classical interpolation
-//    HYPRE_BoomerAMGSetCoarsenType(solver, 6);
-//    HYPRE_BoomerAMGSetCycleType(solver, 1);
-    HYPRE_BoomerAMGSetRelaxType(solver, relax_type);
-    HYPRE_BoomerAMGSetRelaxOrder(solver, relax_order);
-    HYPRE_BoomerAMGSetNumSweeps(solver, num_sweeps);
-    HYPRE_BoomerAMGSetStrongThreshold(solver, strong_threshold);
-
-    int logging = (verbose >= 2) ? 1 : 0;
-    HYPRE_BoomerAMGSetLogging(solver, logging);
-
-    HYPRE_ParCSRMatrix par_A = NULL;
-    HYPRE_IJMatrixGetObject(A, (void**)  &par_A);
-    HYPRE_BoomerAMGSetup(solver, par_A, NULL, NULL);
 }
 
 void

Src/Extern/HYPRE/AMReX_HypreIJIface.H

@@ -0,0 +1,132 @@
+#ifndef AMREX_HYPREIJIFACE_H
+#define AMREX_HYPREIJIFACE_H
+
+#include <memory>
+
+#include <AMReX_MultiFab.H>
+
+#include <HYPRE.h>
+#include <HYPRE_parcsr_ls.h>
+#include <HYPRE_parcsr_mv.h>
+
+namespace amrex {
+
+class HypreIJIface
+{
+public:
+    using HypreIntType = HYPRE_Int;
+    using HypreRealType = HYPRE_Real;
+
+    HypreIJIface(MPI_Comm comm, const HypreIntType ilower, const HypreIntType iupper,
+                 int verbose);
+
+    ~HypreIJIface();
+
+    void parse_inputs(const std::string& prefix = "hypre");
+
+    void solve(const HypreRealType rel_tol, const HypreRealType abs_tol, const HypreIntType max_iter);
+
+    //! IJ matrix instance
+    HYPRE_IJMatrix A() { return m_mat; }
+
+    //! Right hand side IJ vector instance
+    HYPRE_IJVector b() { return m_rhs; }
+
+    //! Solution IJ vector instance
+    HYPRE_IJVector x() { return m_sln; }
+
+    //! Number of iterations taken by the solver to reach the desired tolerance
+    HypreIntType getNumIters() { return m_num_iterations; }
+
+    //! Final residual norm after a linear solve
+    HypreRealType getFinalResidualNorm() { return m_final_res_norm; }
+
+private:
+    void init_preconditioner(const std::string& prefix, const std::string& name);
+    void init_solver(const std::string& prefix, const std::string& name);
+
+    // Preconditioners
+    void boomeramg_precond_configure(const std::string& prefix);
+    void euclid_precond_configure(const std::string& prefix);
+
+    // Solvers
+    void boomeramg_solver_configure(const std::string& prefix);
+    void gmres_solver_configure(const std::string& prefix);
+    void cogmres_solver_configure(const std::string& prefix);
+    void lgmres_solver_configure(const std::string& prefix);
+    void flex_gmres_solver_configure(const std::string& prefix);
+    void bicgstab_solver_configure(const std::string& prefix);
+    void pcg_solver_configure(const std::string& prefix);
+    void hybrid_solver_configure(const std::string& prefix);
+
+
+    MPI_Comm m_comm{MPI_COMM_NULL};
+
+    HYPRE_IJMatrix m_mat{nullptr};
+    HYPRE_IJVector m_rhs{nullptr};
+    HYPRE_IJVector m_sln{nullptr};
+
+    HYPRE_ParCSRMatrix m_parA{nullptr};
+    HYPRE_ParVector m_parRhs{nullptr};
+    HYPRE_ParVector m_parSln{nullptr};
+
+    HYPRE_Solver m_solver{nullptr};
+    HYPRE_Solver m_precond{nullptr};
+
+    HypreIntType (*m_solverDestroyPtr)(HYPRE_Solver){nullptr};
+    HypreIntType (*m_solverSetupPtr)(
+        HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector){nullptr};
+    HypreIntType (*m_solverSolvePtr)(
+        HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector){nullptr};
+    HypreIntType (*m_solverPrecondPtr)(
+        HYPRE_Solver,
+        HYPRE_PtrToParSolverFcn,
+        HYPRE_PtrToParSolverFcn,
+        HYPRE_Solver){nullptr};
+
+    HypreIntType (*m_precondDestroyPtr)(HYPRE_Solver){nullptr};
+    HypreIntType (*m_precondSetupPtr)(
+        HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector){nullptr};
+    HypreIntType (*m_precondSolvePtr)(
+        HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector){nullptr};
+
+    HypreIntType (*m_solverSetTolPtr)(HYPRE_Solver, double){nullptr};
+    HypreIntType (*m_solverSetAbsTolPtr)(HYPRE_Solver, double){nullptr};
+    HypreIntType (*m_solverSetMaxIterPtr)(HYPRE_Solver, HypreIntType){nullptr};
+    HypreIntType (*m_solverNumItersPtr)(HYPRE_Solver, HypreIntType*){nullptr};
+    HypreIntType (*m_solverFinalResidualNormPtr)(HYPRE_Solver, double*){nullptr};
+
+    HypreIntType m_ilower{0};
+    HypreIntType m_iupper{0};
+
+    HypreRealType m_final_res_norm;
+    HypreIntType m_num_iterations;
+
+    std::string m_solver_name{"BoomerAMG"};
+    std::string m_preconditioner_name{"none"};
+    std::string m_file_prefix{"IJ"};
+
+    //! Verbosity of the HYPRE solvers
+    int m_verbose{0};
+
+    unsigned int m_write_counter{0};
+
+    //! Flag indicating whether a preconditioner has been set
+    bool m_has_preconditioner{false};
+
+    //! Flag indicating whether the solver/preconditioner has been setup
+    bool m_need_setup{true};
+
+    //! Flag indicating whether user has requested recomputation of preconditioner
+    bool m_recompute_preconditioner{true};
+
+    //! Flag indicating whether to dump matrix files
+    bool m_write_files{false};
+
+    //! Flag indicating whether the files are overwritten on subsequent writes
+    bool m_overwrite_files{true};
+};
+
+} // namespace amrex
+
+#endif /* AMREX_HYPREIJIFACE_H */