hypre_solver.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
#ifndef OOMPH_HYPRE_SOLVER_HEADER
#define OOMPH_HYPRE_SOLVER_HEADER
 
// Headers required for hypre
#include "_hypre_utilities.h"
#include "HYPRE.h"
#include "HYPRE_parcsr_ls.h"
#include "HYPRE_krylov.h"
#include "HYPRE_IJ_mv.h"
#include "HYPRE_parcsr_mv.h"
 
// OOMPH-LIB includes
#include "linear_algebra_distribution.h"
#include "linear_solver.h"
#include "preconditioner.h"
 
// hypre's global error flag
extern int hypre__global_error;
 
 
namespace oomph
{
  //==================================================================
  //==================================================================
  namespace HypreHelpers
  {
    extern double AMG_strength;
 
    extern unsigned AMG_coarsening;
 
    extern double AMG_truncation; // 0.0;
 
    int check_HYPRE_error_flag(std::ostringstream& message);
 
    extern void create_HYPRE_Vector(const DoubleVector& oomph_vec,
                                    const LinearAlgebraDistribution* dist_pt,
                                    HYPRE_IJVector& hypre_ij_vector,
                                    HYPRE_ParVector& hypre_par_vector);
 
    void create_HYPRE_Vector(const LinearAlgebraDistribution* oomph_vec,
                             HYPRE_IJVector& hypre_ij_vector,
                             HYPRE_ParVector& hypre_par_vector);
 
    void create_HYPRE_Matrix(CRDoubleMatrix* oomph_matrix,
                             HYPRE_IJMatrix& hypre_ij_matrix,
                             HYPRE_ParCSRMatrix& hypre_par_matrix,
                             LinearAlgebraDistribution* dist_pt);
 
    void euclid_settings_helper(const bool& use_block_jacobi,
                                const bool& use_row_scaling,
                                const bool& use_ilut,
                                const int& level,
                                const double& drop_tol,
                                const int& print_level,
                                HYPRE_Solver& euclid_object);
 
  } // namespace HypreHelpers
 
 
  //=====================================================================
  //====================================================================
  class HypreInterface
  {
  public:
    HypreInterface()
    {
#ifdef PARANOID
#ifndef HYPRE_SEQUENTIAL
      // For the MPI version of Hypre, check MPI_Helpers::setup has been called
      if (MPI_Helpers::mpi_has_been_initialised() == false)
      {
        std::ostringstream error_message;
        error_message << "When using the MPI version of Hypre please first\n"
                      << "call function MPI_Helpers::setup()\n";
        throw OomphLibError(error_message.str(),
                            OOMPH_CURRENT_FUNCTION,
                            OOMPH_EXCEPTION_LOCATION);
      }
#endif
#endif
 
      // setup the distribution
      Hypre_distribution_pt = new LinearAlgebraDistribution();
 
      // These keep track of which solver and preconditioner
      // (if any) currently exist
      Existing_solver = None;
      Existing_preconditioner = None;
 
      // Do we want to output info and results of timings?
      Output_info = true;
 
      // General control paramaters
      Tolerance = 1e-10;
      Max_iter = 100;
      Hypre_method = GMRES;
      Internal_preconditioner = None;
 
      // Default AMG control parameters -- these seem OK;
      // see hypre documenation for details.
      AMG_using_simple_smoothing = true;
      if (MPI_Helpers::communicator_pt()->nproc() > 1)
      {
        // Jacobi in parallel
        AMG_simple_smoother = 0;
      }
      else
      {
        // Gauss Seidel in serial
        AMG_simple_smoother = 1;
      }
      AMG_complex_smoother = 6;
      AMG_smoother_iterations = 2;
      AMG_damping = 1.0;
      AMG_strength = HypreHelpers::AMG_strength; // 0.25;
      AMG_truncation = HypreHelpers::AMG_truncation; // 0.0;
      AMG_coarsening = HypreHelpers::AMG_coarsening; // 0
      AMG_max_levels = 100;
      AMG_max_row_sum = 1.0;
      AMG_print_level = 0;
 
      // Parameters for using Euclicd as an AMG smoother (defaults copied
      // from the normal defaults listed in the manual).
      AMGEuclidSmoother_use_block_jacobi = false;
      AMGEuclidSmoother_use_row_scaling = false;
      AMGEuclidSmoother_use_ilut = false;
      AMGEuclidSmoother_level = 1;
      AMGEuclidSmoother_drop_tol = 0; // No dropping
      AMGEuclidSmoother_print_level = 0;
 
      // Print level for CG, GMRES and BiCGStab
      Krylov_print_level = 0;
 
      // Set ParaSails default values
      ParaSails_symmetry = 2;
      ParaSails_nlevel = 1;
      ParaSails_thresh = 0.1;
      ParaSails_filter = 0.1;
 
      // Set Euclid default values
      Euclid_droptol = 0.0;
      Euclid_rowScale = false;
      Euclid_using_ILUT = false;
      Euclid_using_BJ = false;
      Euclid_level = 1;
      Euclid_print_level = 0;
 
      // Set to true to periodically check the hypre error flag
      // and output any messages
      Hypre_error_messages = false;
    }
 
    ~HypreInterface()
    {
      // call function to delete solver data
      hypre_clean_up_memory();
 
      // delete teh oomph-lib distribution
      delete Hypre_distribution_pt;
    }
 
    HypreInterface(const HypreInterface&) = delete;
 
    void operator=(const HypreInterface&) = delete;
 
    void enable_hypre_error_messages()
    {
      Hypre_error_messages = true;
    }
 
    void disable_hypre_error_messages()
    {
      Hypre_error_messages = false;
    }
 
    enum Hypre_method_types
    {
      CG,
      GMRES,
      BiCGStab,
      BoomerAMG,
      Euclid,
      ParaSails,
      None
    };
 
    unsigned existing_solver()
    {
      return Existing_solver;
    }
 
    unsigned existing_preconditioner()
    {
      return Existing_preconditioner;
    }
 
    //??ds comment, write access functions and move to protected?
    bool AMGEuclidSmoother_use_block_jacobi;
    bool AMGEuclidSmoother_use_row_scaling;
    bool AMGEuclidSmoother_use_ilut;
    unsigned AMGEuclidSmoother_level;
    double AMGEuclidSmoother_drop_tol;
    unsigned AMGEuclidSmoother_print_level;
 
  protected:
    void hypre_clean_up_memory();
 
    void hypre_matrix_setup(CRDoubleMatrix* matrix_pt);
 
    void hypre_solver_setup();
 
    void hypre_solve(const DoubleVector& rhs, DoubleVector& solution);
 
    bool Output_info;
 
    unsigned Max_iter;
 
    double Tolerance;
 
    unsigned Hypre_method;
 
    unsigned Internal_preconditioner;
 
    unsigned AMG_print_level;
 
    unsigned AMG_max_levels;
 
    double AMG_max_row_sum;
 
    bool AMG_using_simple_smoothing;
 
    unsigned AMG_simple_smoother;
 
    unsigned AMG_complex_smoother;
 
    unsigned AMG_smoother_iterations;
 
    double AMG_damping;
 
    double AMG_strength;
 
    double AMG_truncation;
 
    unsigned AMG_coarsening;
 
    int ParaSails_symmetry;
 
    int ParaSails_nlevel;
 
    double ParaSails_thresh;
 
    double ParaSails_filter;
 
    double Euclid_droptol;
 
    bool Euclid_rowScale;
 
    bool Euclid_using_ILUT;
 
    bool Euclid_using_BJ;
 
    int Euclid_level;
 
    unsigned Euclid_print_level;
 
    unsigned Krylov_print_level;
 
    bool Hypre_error_messages;
 
    bool Delete_input_data;
 
#ifdef OOMPH_HAS_MPI
    bool Using_distributed_rhs;
 
    bool Returning_distributed_solution;
#endif
 
  private:
    bool Delete_matrix;
 
    HYPRE_IJMatrix Matrix_ij;
 
    HYPRE_ParCSRMatrix Matrix_par;
 
    HYPRE_Solver Solver;
 
    HYPRE_Solver Preconditioner;
 
    unsigned Existing_solver;
 
    unsigned Existing_preconditioner;
 
    LinearAlgebraDistribution* Hypre_distribution_pt;
  };
 
 
 
  //=====================================================================
  //====================================================================
  class HypreSolver : public LinearSolver, public HypreInterface
  {
  public:
    HypreSolver()
    {
      // Hypre copies matrix data from oomph-lib's CRDoubleMatrix
      // and Distributed CRDoubleMatrix into its own data structures,
      // doubling the memory requirements for the matrix.
      // As far as the Hypre solvers are concerned the oomph-lib matrix
      // is no longer required and could be deleted to save memory.
      // However, this will not be the expected behaviour and therefore
      // needs to be requested explicitly by the user by changing this
      // flag from false (its default) to true.
      Delete_matrix = false;
 
      // Do we want to output results of timings?
      Doc_time = true;
    }
 
    ~HypreSolver() {}
 
    HypreSolver(const HypreSolver&) = delete;
 
    void operator=(const HypreSolver&) = delete;
 
    void disable_resolve()
    {
      Enable_resolve = false;
      clean_up_memory();
    }
 
    unsigned& max_iter()
    {
      return Max_iter;
    }
 
    double& tolerance()
    {
      return Tolerance;
    }
 
    unsigned& hypre_method()
    {
      return Hypre_method;
    }
 
    unsigned& internal_preconditioner()
    {
      return Internal_preconditioner;
    }
 
    void amg_using_simple_smoothing()
    {
      AMG_using_simple_smoothing = true;
    }
 
    unsigned& amg_simple_smoother()
    {
      return AMG_simple_smoother;
    }
 
    void amg_using_complex_smoothing()
    {
      AMG_using_simple_smoothing = false;
    }
 
    unsigned& amg_complex_smoother()
    {
      return AMG_complex_smoother;
    }
 
    unsigned& amg_print_level()
    {
      return AMG_print_level;
    }
 
    unsigned& amg_smoother_iterations()
    {
      return AMG_smoother_iterations;
    }
 
    unsigned& amg_coarsening()
    {
      return AMG_coarsening;
    }
 
    unsigned& amg_max_levels()
    {
      return AMG_max_levels;
    }
 
    double& amg_damping()
    {
      return AMG_damping;
    }
 
    double& amg_strength()
    {
      return AMG_strength;
    }
 
    double& amg_max_row_sum()
    {
      return AMG_max_row_sum;
    }
 
    double& amg_truncation()
    {
      return AMG_truncation;
    }
 
    int& parasails_symmetry()
    {
      return ParaSails_symmetry;
    }
 
    int& parasails_nlevel()
    {
      return ParaSails_nlevel;
    }
 
    double& parasails_thresh()
    {
      return ParaSails_thresh;
    }
 
    double& parasails_filter()
    {
      return ParaSails_filter;
    }
 
    double& euclid_droptol()
    {
      return Euclid_droptol;
    }
 
    int& euclid_level()
    {
      return Euclid_level;
    }
 
    void enable_euclid_rowScale()
    {
      Euclid_rowScale = true;
    }
 
    void disable_euclid_rowScale()
    {
      Euclid_rowScale = false;
    }
 
    void enable_euclid_using_BJ()
    {
      Euclid_using_BJ = true;
    }
 
    void disable_euclid_using_BJ()
    {
      Euclid_using_BJ = false;
    }
 
    void euclid_using_ILUK()
    {
      Euclid_using_ILUT = false;
    }
 
    void euclid_using_ILUT()
    {
      Euclid_using_ILUT = true;
    }
 
    unsigned& euclid_print_level()
    {
      return Euclid_print_level;
    }
 
    unsigned& krylov_print_level()
    {
      return Krylov_print_level;
    }
 
    void enable_delete_matrix()
    {
      Delete_matrix = true;
    }
 
    void disable_delete_matrix()
    {
      Delete_matrix = false;
    }
 
 
    void solve(Problem* const& problem_pt, DoubleVector& solution);
 
    void solve(DoubleMatrixBase* const& matrix_pt,
               const DoubleVector& rhs,
               DoubleVector& solution);
 
    void resolve(const DoubleVector& rhs, DoubleVector& solution);
 
    void clean_up_memory();
 
  private:
    bool Delete_matrix;
  };
 
 
 
  //=====================================================================
  //====================================================================
  class HyprePreconditioner : public Preconditioner, public HypreInterface
  {
  public:
    HyprePreconditioner(const std::string& context_string = "")
    {
      Context_string = context_string;
 
      // Hypre copies matrix data from oomph-lib's CRDoubleMatrix
      // or DistributedCRDoubleMatrix into its own data structures,
      // doubling the memory requirements for the matrix.
      // As far as the Hypre solvers are concerned the oomph-lib matrix
      // is no longer required and could be deleted to save memory.
      // However, this will not be the expected behaviour and therefore
      // needs to be requested explicitly by the user by changing this
      // flag from false (its default) to true.
      Delete_matrix = false;
 
      // Do we want to output results of timings?
      Doc_time = true;
 
      // General control paramaters for using HYPRE as preconditioner
      Tolerance = 0.0;
      Max_iter = 1;
      Hypre_method = BoomerAMG;
      Internal_preconditioner = None;
 
      // Initialise private double that accumulates the preconditioner
      // solve time of thi instantiation of this class. Is reported
      // in destructor if Report_my_cumulative_preconditioner_solve_time
      // is set to true
      My_cumulative_preconditioner_solve_time = 0.0;
      Report_my_cumulative_preconditioner_solve_time = false;
    }
 
    ~HyprePreconditioner()
    {
      if (Report_my_cumulative_preconditioner_solve_time)
      {
        oomph_info << "~HyprePreconditioner() in context \" " << Context_string
                   << "\": My_cumulative_preconditioner_solve_time = "
                   << My_cumulative_preconditioner_solve_time << std::endl;
      }
    }
 
    HyprePreconditioner(const HyprePreconditioner&) = delete;
 
    void operator=(const HyprePreconditioner&) = delete;
 
    static double Cumulative_preconditioner_solve_time;
 
    static std::map<std::string, double> Context_based_cumulative_solve_time;
 
 
    static unsigned Cumulative_npreconditioner_solve;
 
    static std::map<std::string, unsigned>
      Context_based_cumulative_npreconditioner_solve;
 
    static std::map<std::string, unsigned> Context_based_nrow;
 
    static void report_cumulative_solve_times();
 
    static void reset_cumulative_solve_times();
 
    void enable_report_my_cumulative_preconditioner_solve_time()
    {
      Report_my_cumulative_preconditioner_solve_time = true;
    }
 
    void disable_report_my_cumulative_preconditioner_solve_time()
    {
      Report_my_cumulative_preconditioner_solve_time = false;
    }
 
    void enable_doc_time()
    {
      Doc_time = true;
    }
 
    void disable_doc_time()
    {
      Doc_time = false;
    }
 
    unsigned& hypre_method()
    {
      return Hypre_method;
    }
 
    unsigned& internal_preconditioner()
    {
      return Internal_preconditioner;
    }
 
    void use_BoomerAMG()
    {
      Hypre_method = BoomerAMG;
    }
 
    void set_amg_iterations(const unsigned& amg_iterations)
    {
      Max_iter = amg_iterations;
    }
 
    unsigned& amg_iterations()
    {
      return Max_iter;
    }
 
    void amg_using_simple_smoothing()
    {
      AMG_using_simple_smoothing = true;
    }
 
    unsigned& amg_simple_smoother()
    {
      return AMG_simple_smoother;
    }
 
    void amg_using_complex_smoothing()
    {
      AMG_using_simple_smoothing = false;
    }
 
    unsigned& amg_complex_smoother()
    {
      return AMG_complex_smoother;
    }
 
    bool& amg_using_simple_smoothing_flag()
    {
      return AMG_using_simple_smoothing;
    }
 
    unsigned& amg_print_level()
    {
      return AMG_print_level;
    }
 
    unsigned& amg_smoother_iterations()
    {
      return AMG_smoother_iterations;
    }
 
    unsigned& amg_coarsening()
    {
      return AMG_coarsening;
    }
 
    unsigned& amg_max_levels()
    {
      return AMG_max_levels;
    }
 
    double& amg_damping()
    {
      return AMG_damping;
    }
 
    double& amg_strength()
    {
      return AMG_strength;
    }
 
    double& amg_max_row_sum()
    {
      return AMG_max_row_sum;
    }
 
    double& amg_truncation()
    {
      return AMG_truncation;
    }
 
    void use_ParaSails()
    {
      Hypre_method = ParaSails;
    }
 
    int& parasails_symmetry()
    {
      return ParaSails_symmetry;
    }
 
    int& parasails_nlevel()
    {
      return ParaSails_nlevel;
    }
 
    double& parasails_thresh()
    {
      return ParaSails_thresh;
    }
 
    double& parasails_filter()
    {
      return ParaSails_filter;
    }
 
    void use_Euclid()
    {
      Hypre_method = Euclid;
    }
 
    double& euclid_droptol()
    {
      return Euclid_droptol;
    }
 
    int& euclid_level()
    {
      return Euclid_level;
    }
 
    void enable_euclid_rowScale()
    {
      Euclid_rowScale = true;
    }
 
    void disable_euclid_rowScale()
    {
      Euclid_rowScale = false;
    }
 
    void enable_euclid_using_BJ()
    {
      Euclid_using_BJ = true;
    }
 
    void disable_euclid_using_BJ()
    {
      Euclid_using_BJ = false;
    }
 
    void euclid_using_ILUK()
    {
      Euclid_using_ILUT = false;
    }
 
    void euclid_using_ILUT()
    {
      Euclid_using_ILUT = true;
    }
 
    unsigned& euclid_print_level()
    {
      return Euclid_print_level;
    }
 
    void enable_delete_matrix()
    {
      Delete_matrix = true;
    }
 
    void disable_delete_matrix()
    {
      Delete_matrix = false;
    }
 
    void setup();
 
    void preconditioner_solve(const DoubleVector& r, DoubleVector& z);
 
    void clean_up_memory();
 
  private:
    bool Delete_matrix;
 
    // Flag is true to output results of timings
    bool Doc_time;
 
    double My_cumulative_preconditioner_solve_time;
 
    bool Report_my_cumulative_preconditioner_solve_time;
 
    std::string Context_string;
  };
 
 
 
  //==================================================================
  //==================================================================
  namespace Hypre_default_settings
  {
    extern void set_defaults_for_navier_stokes_momentum_block(
      HyprePreconditioner* hypre_preconditioner_pt);
 
    extern void set_defaults_for_2D_poisson_problem(
      HyprePreconditioner* hypre_preconditioner_pt);
 
    extern void set_defaults_for_3D_poisson_problem(
      HyprePreconditioner* hypre_preconditioner_pt);
 
  } // namespace Hypre_default_settings
} // namespace oomph
 
#endif