AmrMultiGrid.h

Go to the documentation of this file.

//
// Class AmrMultiGrid
//   Main class of the AMR Poisson multigrid solver.
//   It implements the multigrid solver described in https://doi.org/10.1016/j.cpc.2019.106912
//
// Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef AMR_MULTI_GRID_H
#define AMR_MULTI_GRID_H

#include <vector>
#include <memory>

#include <AMReX.H>
#include <AMReX_MultiFab.H>
#include <AMReX_Vector.H>

#include "AmrMultiGridCore.h"

#include "Solvers/AmrPoissonSolver.h"
#include "Amr/AmrBoxLib.h"

#include "AmrMultiGridLevel.h"

#include <fstream>

#define AMR_MG_TIMER true
#define AMR_MG_WRITE false

class AmrMultiGrid : public AmrPoissonSolver< AmrBoxLib > {
    
public:
    typedef amr::matrix_t         matrix_t;
    typedef amr::vector_t         vector_t;
    typedef amr::multivector_t    mv_t;
    typedef amr::dmap_t           dmap_t;
    typedef amr::comm_t           comm_t;
    typedef amr::local_ordinal_t  lo_t;
    typedef amr::global_ordinal_t go_t;
    typedef amr::scalar_t         scalar_t;
    
    typedef AmrMultiGridLevel<matrix_t, vector_t> AmrMultiGridLevel_t;
    
    typedef AmrMultiGridLevel_t::AmrField_t     AmrField_t;
    typedef AmrMultiGridLevel_t::AmrGeometry_t  AmrGeometry_t;
    typedef AmrMultiGridLevel_t::AmrField_u     AmrField_u;
    typedef AmrMultiGridLevel_t::AmrField_s     AmrField_s;
    typedef AmrMultiGridLevel_t::AmrIntVect_t   AmrIntVect_t;
    typedef AmrMultiGridLevel_t::indices_t      indices_t;
    typedef AmrMultiGridLevel_t::coefficients_t coefficients_t;
    typedef AmrMultiGridLevel_t::umap_t         umap_t;
    typedef AmrMultiGridLevel_t::boundary_t     boundary_t;
    
    typedef BottomSolver<
        Teuchos::RCP<matrix_t>,
        Teuchos::RCP<mv_t>,
        AmrMultiGridLevel_t
    > bsolver_t;
    
    typedef BelosBottomSolver<AmrMultiGridLevel_t>      BelosSolver_t;
    typedef Amesos2BottomSolver<AmrMultiGridLevel_t>    Amesos2Solver_t;
    typedef MueLuBottomSolver<AmrMultiGridLevel_t>      MueLuSolver_t;
    
    typedef AmrPreconditioner<matrix_t, AmrMultiGridLevel_t> preconditioner_t;
    
    typedef Ifpack2Preconditioner<AmrMultiGridLevel_t> Ifpack2Preconditioner_t;
    typedef MueLuPreconditioner<AmrMultiGridLevel_t> MueLuPreconditioner_t;
    
    typedef amrex::BoxArray boxarray_t;
    typedef amrex::Box box_t;
    typedef amrex::BaseFab<int> basefab_t;
    typedef amrex::FArrayBox farraybox_t;
            
    typedef AmrSmoother::Smoother Smoother;
    
    typedef amr::Preconditioner Preconditioner;
    
    enum Interpolater {
        TRILINEAR = 0,
        LAGRANGE,
        PIECEWISE_CONST
    };
    
    enum BaseSolver {
        // all Belos
        BICGSTAB,
        MINRES,
        PCPG,
        CG,
        GMRES,
        STOCHASTIC_CG,
        RECYCLING_CG,
        RECYCLING_GMRES
        // all Amesos2
#ifdef HAVE_AMESOS2_KLU2
        , KLU2
#endif
#ifdef HAVE_AMESOS2_SUPERLU
        , SUPERLU
#endif
#ifdef HAVE_AMESOS2_UMFPACK  
        , UMFPACK
#endif
#ifdef HAVE_AMESOS2_PARDISO_MKL
        , PARDISO_MKL
#endif
#ifdef HAVE_AMESOS2_MUMPS
        , MUMPS
#endif
#ifdef HAVE_AMESOS2_LAPACK
        , LAPACK
#endif
        // all MueLu
        , SA
        // add others ...
    };
    
    enum Boundary {
        DIRICHLET = 0,
        OPEN,
        PERIODIC
    };
    
    enum Norm {
        L1,
        L2,
        LINF
    };
    
public:
    
    AmrMultiGrid(AmrBoxLib* itsAmrObject_p,
                 const std::string& bsolver,
                 const std::string& prec,
                 const bool& rebalance,
                 const std::string& reuse,
                 const std::string& bcx,
                 const std::string& bcy,
                 const std::string& bcz,
                 const std::string& smoother,
                 const std::size_t& nSweeps,
                 const std::string& interp,
                 const std::string& norm);
    
    void solve(AmrScalarFieldContainer_t &rho,
               AmrScalarFieldContainer_t &phi,
               AmrVectorFieldContainer_t &efield,
               unsigned short baseLevel,
               unsigned short finestLevel,
               bool prevAsGuess = true);
    
    void setNumberOfSweeps(const std::size_t& nSweeps);
    
    void setMaxNumberOfIterations(const std::size_t& maxiter);
    
    std::size_t getNumIters();
    
    scalar_t getLevelResidualNorm(lo_t level);
    
    void setVerbose(bool verbose);


    void setTolerance(const scalar_t& eps);
    
    double getXRangeMin(unsigned short level = 0);
    double getXRangeMax(unsigned short level = 0);
    double getYRangeMin(unsigned short level = 0);
    double getYRangeMax(unsigned short level = 0);
    double getZRangeMin(unsigned short level = 0);
    double getZRangeMax(unsigned short level = 0);
    
    Inform &print(Inform &os) const;
    
private:
    
    void initPhysicalBoundary_m(const Boundary* bc);
    
    void initLevels_m(const amrex::Vector<AmrField_u>& rho,
                      const amrex::Vector<AmrGeometry_t>& geom,
                      bool regrid);
    
    void clearMasks_m();
    
    void initGuess_m(bool reset);
    
    scalar_t iterate_m();
    
    bool isConverged_m(std::vector<scalar_t>& rhsNorms,
                       std::vector<scalar_t>& resNorms);
    
    void residual_m(const lo_t& level,
                    Teuchos::RCP<vector_t>& r,
                    const Teuchos::RCP<vector_t>& b,
                    const Teuchos::RCP<vector_t>& x);
    
    void relax_m(const lo_t& level);
    
    void residual_no_fine_m(const lo_t& level,
                            Teuchos::RCP<vector_t>& result,
                            const Teuchos::RCP<vector_t>& rhs,
                            const Teuchos::RCP<vector_t>& crs_rhs,
                            const Teuchos::RCP<vector_t>& b);
    
#if AMR_MG_WRITE

    void writeResidualNorm_m();
#endif
    
    scalar_t evalNorm_m(const Teuchos::RCP<const vector_t>& x);
    
    void initResidual_m(std::vector<scalar_t>& rhsNorms,
                        std::vector<scalar_t>& resNorms);
    
    void computeEfield_m(AmrVectorFieldContainer_t& efield);
    
    void setup_m(const amrex::Vector<AmrField_u>& rho,
                 const amrex::Vector<AmrField_u>& phi,
                 const bool& matrices = true);
    
    void buildSingleLevel_m(const amrex::Vector<AmrField_u>& rho,
                            const amrex::Vector<AmrField_u>& phi,
                            const bool& matrices = true);
    
    void buildMultiLevel_m(const amrex::Vector<AmrField_u>& rho,
                           const amrex::Vector<AmrField_u>& phi,
                           const bool& matrices = true);

    void open_m(const lo_t& level, const bool& matrices);
    
    void close_m(const lo_t& level, const bool& matrices);
    
    void buildNoFinePoissonMatrix_m(const lo_t& level,
                                    const go_t& gidx,
                                    const AmrIntVect_t& iv,
                                    const basefab_t& mfab,
                                    const scalar_t* invdx2);
    
    void buildCompositePoissonMatrix_m(const lo_t& level,
                                       const go_t& gidx,
                                       const AmrIntVect_t& iv,
                                       const basefab_t& mfab,
                                       const basefab_t& rfab,
                                       const scalar_t* invdx2);
    
    void buildRestrictionMatrix_m(const lo_t& level,
                                  const go_t& gidx,
                                  const AmrIntVect_t& iv,
                                  D_DECL(const go_t& ii,
                                         const go_t& jj,
                                         const go_t& kk),
                                  const basefab_t& rfab);
    
    void buildInterpolationMatrix_m(const lo_t& level,
                                    const go_t& gidx,
                                    const AmrIntVect_t& iv,
                                    const basefab_t& cfab);
    
    void buildCrseBoundaryMatrix_m(const lo_t& level,
                                   const go_t& gidx,
                                   const AmrIntVect_t& iv,
                                   const basefab_t& mfab,
                                   const basefab_t& cfab,
                                   const scalar_t* invdx2);
    
    void buildFineBoundaryMatrix_m(const lo_t& level,
                                   const go_t& gidx,
                                   const AmrIntVect_t& iv,
                                   const basefab_t& mfab,
                                   const basefab_t& rfab);
    
    void buildDensityVector_m(const lo_t& level,
                              const AmrField_t& rho);
    
    void buildPotentialVector_m(const lo_t& level,
                                const AmrField_t& phi);
    
    void buildGradientMatrix_m(const lo_t& level,
                               const go_t& gidx,
                               const AmrIntVect_t& iv,
                               const basefab_t& mfab,
                               const scalar_t* invdx);
    
    void amrex2trilinos_m(const lo_t& level,
                          const lo_t& comp,
                          const AmrField_t& mf,
                          Teuchos::RCP<vector_t>& mv);
    
    void trilinos2amrex_m(const lo_t& level,
                          const lo_t& comp,
                          AmrField_t& mf,
                          const Teuchos::RCP<vector_t>& mv);
    
    void map2vector_m(umap_t& map, indices_t& indices,
                      coefficients_t& values);
    
    void smooth_m(const lo_t& level,
                  Teuchos::RCP<vector_t>& e,
                  Teuchos::RCP<vector_t>& r);
    
    void restrict_m(const lo_t& level);
    
    void prolongate_m(const lo_t& level);
    
    void averageDown_m(const lo_t& level);
    
    void initInterpolater_m(const Interpolater& interp);
    
    void initCrseFineInterp_m(const Interpolater& interface);
    
    void initBaseSolver_m(const BaseSolver& solver,
                          const bool& rebalance,
                          const std::string& reuse);
    
    void initPrec_m(const Preconditioner& prec,
                    const std::string& reuse);
    
    Boundary convertToEnumBoundary_m(const std::string& bc);
    
    Interpolater convertToEnumInterpolater_m(const std::string& interp);
    
    BaseSolver convertToEnumBaseSolver_m(const std::string& bsolver);
    
    Preconditioner convertToEnumPreconditioner_m(const std::string& prec);
    
    Smoother convertToEnumSmoother_m(const std::string& smoother);
    
    Norm convertToEnumNorm_m(const std::string& norm);
    
    void writeSDDSHeader_m(std::ofstream& outfile);
    
    void writeSDDSData_m(const scalar_t& error);
    
#if AMR_MG_TIMER

    void initTimer_m();
#endif
    
private:
    Teuchos::RCP<comm_t> comm_mp;       
    
    std::unique_ptr<AmrInterpolater<AmrMultiGridLevel_t> > interp_mp;
    
    std::unique_ptr<AmrInterpolater<AmrMultiGridLevel_t> > interface_mp;
    
    std::size_t nIter_m;            
    std::size_t bIter_m;            
    std::size_t maxiter_m;          
    std::size_t nSweeps_m;          
    Smoother smootherType_m;        
    
    std::vector<std::unique_ptr<AmrMultiGridLevel_t > > mglevel_m;
    
    std::shared_ptr<bsolver_t> solver_mp;
    
    std::vector<std::shared_ptr<AmrSmoother> > smoother_m;
    
    std::shared_ptr<preconditioner_t> prec_mp;
    
    int lbase_m;            
    int lfine_m;            
    int nlevel_m;           
    
    boundary_t bc_m[AMREX_SPACEDIM];    
    int nBcPoints_m;                    
    
    Norm norm_m;            
    std::string snorm_m;    
    
    scalar_t eps_m;         
    
    bool verbose_m;                 
    std::string fname_m;            
    std::ios_base::openmode flag_m; 
    
#if AMR_MG_TIMER
    IpplTimings::TimerRef buildTimer_m;         
    IpplTimings::TimerRef restrictTimer_m;      
    IpplTimings::TimerRef smoothTimer_m;        
    IpplTimings::TimerRef interpTimer_m;        
    IpplTimings::TimerRef efieldTimer_m;        
    IpplTimings::TimerRef averageTimer_m;       
    IpplTimings::TimerRef bottomTimer_m;        
    IpplTimings::TimerRef dumpTimer_m;          
#endif
};


inline Inform &operator<<(Inform &os, const AmrMultiGrid &fs) {
    return fs.print(os);
}

#endif