hypreinterface_sequential.hh

/* *************************************************** *
 *    Hypre solver interface
 * *************************************************** */
#ifndef HYPREINTERFACE_SEQUENTIAL_HH
#define HYPREINTERFACE_SEQUENTIAL_HH
 
#if HAVE_HYPRE
 
#include <vector>
#include <string>
#include <dune/common/parallel/communication.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/grid/common/indexidset.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/solvers.hh>
 
#include <HYPRE.h>
#include <HYPRE_IJ_mv.h>
#include <HYPRE_parcsr_ls.h>
#include <HYPRE_krylov.h>
 
#include "globalDOFNumbering.hh"
 
class SequentialHypreParameters {
    public:
        //TODO set to number of input EV
        int dim = 1; //For blocked aggregation, gives size of block
        /* ! Type of coarsening strategy
           0  CLJP-coarsening (a parallel coarsening algorithm using independent sets.
           1  classical Ruge-Stueben coarsening on each processor, no boundary treatment (not recommended!)
           3  classical Ruge-Stueben coarsening on each processor, followed by a third pass, which adds coarse points on the boundaries
           6  Falgout coarsening (uses 1 first, followed by CLJP using the interior coarse points generated by 1 as its first independent set)
           7  CLJP-coarsening (using a fixed random vector, for debugging purposes only)
           8   PMIS-coarsening (a parallel coarsening algorithm using independent sets, generating lower complexities than CLJP, might also lead to slower convergence)
           9  PMIS-coarsening (using a fixed random vector, for debugging purposes only)
           10 HMIS-coarsening (uses one pass Ruge-Stueben on each processor independently, followed by PMIS using the interior C-points generated as its first independent set)
           11 one-pass Ruge-Stueben coarsening on each processor, no boundary treatment (not recommended!)
           21 CGC coarsening by M. Griebel, B. Metsch and A. Schweitzer
           22 CGC-E coarsening by M. Griebel, B. Metsch and A.Schweitzer
         * */
        int coarsentype = 0;
 
        int interptype = 0;
        //Maximal number of elements per row for interpolation
        int pmaxelmts = 0;
        //number of levels of aggressive coarsening
        int aggnumlevels = 2;
        /* ! defines the smoother
           0 Jacobi
           1 Gauss-Seidel, sequential (very slow!)
           2 Gauss-Seidel, interior points in parallel, boundary sequential (slow!)
           3 hybrid Gauss-Seidel or SOR, forward solve
           4 hybrid Gauss-Seidel or SOR, backward solve
           5 hybrid chaotic Gauss-Seidel (works only with OpenMP)
           6 hybrid symmetric Gauss-Seidel or SSOR
           9 Gaussian elimination (only on coarsest level)
           */
        int relaxtype = 6;
        /* !
           defines which order the points are relaxed
           0 lexiographic
           1 CF relaxation
           */
        int relaxorder = 1;
        double strongthreshold = 0.9;
        int printlevel = 0;
        int maxlevel = 5;
        int coarsesolver = 6;
        int ncoarserelax = 1;
};
 
 
/* *************************************************** *
 * Solver class
 * *************************************************** */
template <typename GV, typename Matrix, typename U>
class SequentialHypreSolver : public Dune::PDELab::LinearResultStorage {
    public:
        // Number of rows and columns in each block of the matrix
        static constexpr int mBlock = Matrix::block_type::rows;
        static constexpr int nBlock = Matrix::block_type::cols;
 
        /* Constructor */
        SequentialHypreSolver (GV& gv_, U& u_, Matrix& a_, U& b_,
                const double tolerance_) :
            gv(gv_),
            u(u_), b(b_), a(a_),
            tolerance(tolerance_),
            maxiter(200),
            ParCSRPCG_printlevel(0) {
                SequentialHypreParameters hypre_param;
 
                // Set up BoomerAMG preconditioner
                HYPRE_BoomerAMGCreate(&prec);
                HYPRE_BoomerAMGSetTol(prec,0);
                HYPRE_BoomerAMGSetMaxIter(prec,1);
 
                HYPRE_BoomerAMGSetNumFunctions(prec, hypre_param.dim);
 
                HYPRE_BoomerAMGSetInterpType(prec, hypre_param.interptype);
                HYPRE_BoomerAMGSetPMaxElmts(prec, hypre_param.pmaxelmts);
                HYPRE_BoomerAMGSetCoarsenType(prec, hypre_param.coarsentype);
                HYPRE_BoomerAMGSetAggNumLevels(prec, hypre_param.aggnumlevels);
                HYPRE_BoomerAMGSetRelaxType(prec, hypre_param.relaxtype);
                HYPRE_BoomerAMGSetRelaxOrder(prec, hypre_param.relaxorder);
                HYPRE_BoomerAMGSetMaxLevels (prec, hypre_param.maxlevel);
                HYPRE_BoomerAMGSetCycleNumSweeps(prec, hypre_param.ncoarserelax,3);
                HYPRE_BoomerAMGSetCycleRelaxType(prec, hypre_param.coarsesolver,3);
                HYPRE_BoomerAMGSetPrintLevel(prec,hypre_param.printlevel);
                HYPRE_BoomerAMGSetStrongThreshold(prec,hypre_param.strongthreshold);
 
                HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD,&solver);
                HYPRE_ParCSRPCGSetTol(solver,tolerance);
                HYPRE_ParCSRPCGSetMaxIter(solver,maxiter);
                HYPRE_ParCSRPCGSetPrintLevel(solver, ParCSRPCG_printlevel);
                HYPRE_PCGSetTwoNorm(solver,1);
 
                HYPRE_ParCSRPCGSetPrecond(solver,
                        HYPRE_BoomerAMGSolve,
                        HYPRE_BoomerAMGSetup,
                        prec);
 
                // Split the communicator into sequential parts and use the original rank for ordering
                MPI_Comm_split(gv.comm(), gv.comm().rank(), gv.comm().rank(), &row_comm);
 
                // Assemble hypre data structures
                AssembleHypreMatrix(a,a_hypre);
                AssembleHypreVector(b,b_hypre);
                AssembleHypreVector(u,u_hypre);
                // Extract ParCSR data
                HYPRE_IJMatrixGetObject(a_hypre,(void **) &a_parcsr);
                HYPRE_IJVectorGetObject(b_hypre,(void **) &b_parcsr);
                HYPRE_IJVectorGetObject(u_hypre,(void **) &u_parcsr);
            }
 
        /* Destructor */
        ~SequentialHypreSolver() {
            HYPRE_BoomerAMGDestroy(prec);
            HYPRE_ParCSRPCGDestroy(solver);
            HYPRE_IJVectorDestroy(b_hypre);
            HYPRE_IJMatrixDestroy(a_hypre);
            MPI_Comm_free(&row_comm);
        }
 
  // Solve for a given RHS
  void solve(U& u) {
      Dune::Timer watch;
 
      //Setup and solve
      watch.reset();
      HYPRE_ParCSRPCGSetup(solver, a_parcsr, b_parcsr, u_parcsr);
      double timeBoomerAMGSetup = watch.elapsed();
      timeBoomerAMGSetup = gv.comm().max(timeBoomerAMGSetup);
      watch.reset();
 
      HYPRE_ParCSRPCGSolve(solver, a_parcsr, b_parcsr,u_parcsr);
      double timeBoomerAMGSolve = watch.elapsed();
      timeBoomerAMGSolve = gv.comm().max(timeBoomerAMGSolve);
      MPI_Barrier(gv.comm());
 
      //Write back to dune vector
      ReadHypreVector(u_hypre,u);
      HYPRE_IJVectorDestroy(u_hypre);
      HYPRE_Int iterations;
      double resreduction;
      HYPRE_ParCSRPCGGetNumIterations(solver,&iterations);
      HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(solver,&resreduction);
 
      //Store solver statistics
      res.converged  = (resreduction <= tolerance);
      if (gv.comm().rank()==0)
          std::cout << " Converged to " << resreduction << " in " << iterations << " iterations." << std::endl;
      res.iterations = iterations;
      res.reduction = resreduction;
      res.elapsed    = timeBoomerAMGSolve;
      res.conv_rate  = pow(resreduction,1./iterations);
  }
 
  /* Read solver tolerance */
  double Tolerance() const { return tolerance; }
 
 private:
  void AssembleHypreMatrix(const Matrix& a, HYPRE_IJMatrix& hypre_matrix) {
      HYPRE_Int ilower, iupper;
      HYPRE_Int jlower, jupper;
      ilower = 0;
      iupper = nBlock*a.N();
      jlower = ilower;
      jupper = iupper;
 
      HYPRE_IJMatrixCreate(row_comm,
              ilower, iupper,
              jlower, jupper,
              &hypre_matrix);
      HYPRE_IJMatrixSetObjectType(hypre_matrix, HYPRE_PARCSR);
      HYPRE_IJMatrixInitialize(hypre_matrix);
      using Dune::PDELab::Backend::native;
      int nrows=1;
      HYPRE_Int ncol=1;
 
      for (auto row_iter = native(a).begin(); row_iter != native(a).end(); row_iter++) {
              for (auto col_iter = row_iter->begin(); col_iter != row_iter->end(); col_iter++) {
                  for(int block_i = 0; block_i < mBlock; block_i++){   //iterator over block values
                      for(int block_j = 0; block_j < mBlock; block_j++){   //iterator over block values
                          HYPRE_Int row_idx = mBlock*row_iter.index()+block_i;
                          HYPRE_Int col_idx = nBlock*col_iter.index()+block_j;
                          double value = (*col_iter)[block_i][block_j];
                          HYPRE_IJMatrixSetValues(hypre_matrix,nrows,&ncol,&row_idx,&col_idx,&value);
                      }
                  }
              }
      }
      HYPRE_IJMatrixAssemble(hypre_matrix);
  }
  /* *************************************************** *
   * Assemble hypre vector given vector u and return
   * *************************************************** */
void AssembleHypreVector(const U& u, HYPRE_IJVector& u_hypre) {
 
    HYPRE_Int jlower = 0;
    HYPRE_Int jupper = nBlock*a.N();
    HYPRE_IJVectorCreate(row_comm, jlower, jupper, &u_hypre);
    HYPRE_IJVectorSetObjectType(u_hypre, HYPRE_PARCSR);
    HYPRE_IJVectorInitialize(u_hypre);
    using Dune::PDELab::Backend::native;
    int nvalues=1;
    for (unsigned int i = 0; i < u.N(); i++) {
            for(int block_i = 0; block_i < mBlock; block_i++){   //iterator over block values
                HYPRE_Int indices = nBlock*i+block_i;
                double values = native(u)[i][block_i];
                HYPRE_IJVectorSetValues(u_hypre, nvalues, &indices, &values);
            }
    }
    HYPRE_IJVectorAssemble(u_hypre);
}
 
/* *************************************************** *
 * Extract from hypre vector u and return
 * *************************************************** */
void ReadHypreVector(const HYPRE_IJVector& u_hypre, U& u) {
    int nvalues=1;
    double values;
    using Dune::PDELab::Backend::native;
    for (unsigned int i = 0; i < u.N(); i++) {
            for(int block_i = 0; block_i < mBlock; block_i++){   //iterator over block values
                HYPRE_Int indices = nBlock*i + block_i;
                HYPRE_IJVectorGetValues(u_hypre, nvalues, &indices, &values);
                native(u)[i][block_i] = values;
            }
        }
}
/* *************************************************** *
 * Class data
 * *************************************************** */
 
const GV gv;    // gridview of trial grid function space
U& u;           // Point where the Jacobian is evaluated
U b;            // RHS
Matrix& a;
 
double tolerance; // ParCSRPCG Solver tolerance
HYPRE_Int maxiter; // ParCSRPCG maxiter
int ParCSRPCG_printlevel;
 
MPI_Comm row_comm;
 
HYPRE_IJMatrix a_hypre;  // Hypre matrix storage
HYPRE_IJVector b_hypre;
HYPRE_IJVector u_hypre;
 
HYPRE_Solver prec;
HYPRE_Solver solver;
 
HYPRE_ParCSRMatrix a_parcsr;
HYPRE_ParVector b_parcsr;
HYPRE_ParVector u_parcsr;
 
};
 
#endif
 
#endif