mg_solver.h

Go to the documentation of this file.

/*
 * Copyright (c) 2010-2015:  G-CSC, Goethe University Frankfurt
 * Author: Andreas Vogel
 * 
 * This file is part of UG4.
 * 
 * UG4 is free software: you can redistribute it and/or modify it under the
 * terms of the GNU Lesser General Public License version 3 (as published by the
 * Free Software Foundation) with the following additional attribution
 * requirements (according to LGPL/GPL v3 §7):
 * 
 * (1) The following notice must be displayed in the Appropriate Legal Notices
 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (2) The following notice must be displayed at a prominent place in the
 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
 * 
 * (3) The following bibliography is recommended for citation and must be
 * preserved in all covered files:
 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
 *   parallel geometric multigrid solver on hierarchically distributed grids.
 *   Computing and visualization in science 16, 4 (2013), 151-164"
 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
 *   flexible software system for simulating pde based models on high performance
 *   computers. Computing and visualization in science 16, 4 (2013), 165-179"
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU Lesser General Public License for more details.
 */
 
#ifndef __H__UG__LIB_DISC__MULTI_GRID_SOLVER__MG_SOLVER__
#define __H__UG__LIB_DISC__MULTI_GRID_SOLVER__MG_SOLVER__
 
// extern includes
#include <vector>
#include <iostream>
 
// other ug4 modules
#include "common/common.h"
#ifdef UG_PARALLEL
  #include "lib_algebra/parallelization/parallel_storage_type.h"
  #include "lib_grid/parallelization/distributed_grid.h"
#endif
#include "lib_algebra/operator/interface/linear_iterator.h"
#include "lib_algebra/operator/interface/operator_inverse.h"
#include "lib_algebra/operator/interface/operator.h"
#include "lib_algebra/operator/preconditioner/jacobi.h"
#include "lib_algebra/operator/linear_solver/lu.h"
#include "lib_disc/dof_manager/dof_distribution.h"
#include "lib_disc/operator/linear_operator/transfer_interface.h"
//only for debugging!!!
#include "lib_grid/algorithms/debug_util.h"
 
// library intern headers
#include "lib_disc/function_spaces/grid_function_util.h"
#include "lib_disc/operator/linear_operator/assembled_linear_operator.h"
 
#include "mg_stats.h"
 
namespace ug{
 
// AssembledMultiGridCycle
 
 
template <typename TDomain, typename TAlgebra>
class AssembledMultiGridCycle :
 public ILinearIterator<  typename TAlgebra::vector_type>
{
  public:
    typedef TDomain domain_type;
 
    typedef TAlgebra algebra_type;
 
    typedef typename algebra_type::vector_type vector_type;
 
    typedef typename algebra_type::matrix_type matrix_type;
 
    typedef GridFunction<TDomain, TAlgebra> GF;
 
    typedef MGStats<TDomain, TAlgebra> mg_stats_type;
 
  //  Setup
 
  public:
    AssembledMultiGridCycle();
 
    AssembledMultiGridCycle(SmartPtr<ApproximationSpace<TDomain> > approxSpace);
 
 
    void set_mg_stats(SmartPtr<mg_stats_type> mgstats)
      {m_mgstats = mgstats;}
 
    void set_approximation_space(SmartPtr<ApproximationSpace<TDomain> > approxSpace);
 
    void set_discretization(SmartPtr<IAssemble<TAlgebra> > spAss)
      {m_spAss = spAss;}
 
    void set_base_level(int baseLevel) {m_baseLev = baseLevel;}
 
    void set_surface_level(int topLevel) {m_surfaceLev = topLevel;}
 
    void set_base_solver(SmartPtr<ILinearOperatorInverse<vector_type> > baseSolver)
      {m_spBaseSolver = baseSolver;}
 
    void set_gathered_base_solver_if_ambiguous(bool bGathered) {m_bGatheredBaseIfAmbiguous = bGathered;}
 
    void set_emulate_full_refined_grid(bool bEmulate){
      if(bEmulate) m_GridLevelType = GridLevel::SURFACE;
      else m_GridLevelType = GridLevel::LEVEL;
    }
 
    void set_rap(bool bRAP) {m_bUseRAP = bRAP;}
 
    void set_smooth_on_surface_rim(bool bSmooth) {m_bSmoothOnSurfaceRim = bSmooth;}
 
    void set_cycle_type(int type) {m_cycleType = type;}
 
    void set_cycle_type(const std::string& type) {
      if(TrimString(type) == "V") {m_cycleType = _V_;}
      else if(TrimString(type) == "W") {m_cycleType = _W_;}
      else if(TrimString(type) == "F") {m_cycleType = _F_;}
      else {UG_THROW("GMG::set_cycle_type: option '"<<type<<"' not supported.");}
    }
 
    void set_comm_comp_overlap(bool bOverlap) {m_bCommCompOverlap = bOverlap;}
 
    void set_num_presmooth(int num) {m_numPreSmooth = num;}
 
    void set_num_postsmooth(int num) {m_numPostSmooth = num;}
 
    void set_smoother(SmartPtr<ILinearIterator<vector_type> > smoother)
      {set_presmoother(smoother); set_postsmoother(smoother);}
 
    void set_presmoother(SmartPtr<ILinearIterator<vector_type> > smoother)
      {m_spPreSmootherPrototype = smoother;}
 
    void set_postsmoother(SmartPtr<ILinearIterator<vector_type> > smoother)
      {m_spPostSmootherPrototype = smoother;}
 
    void set_transfer(SmartPtr<ITransferOperator<TDomain, TAlgebra> > P)
      {set_prolongation(P); set_restriction(P);}
 
    void set_prolongation(SmartPtr<ITransferOperator<TDomain, TAlgebra> > P)
      {m_spProlongationPrototype = P;}
 
    void set_restriction(SmartPtr<ITransferOperator<TDomain, TAlgebra> > P)
      {m_spRestrictionPrototype = P;}
 
    void set_projection(SmartPtr<ITransferOperator<TDomain, TAlgebra> > P)
      {m_spProjectionPrototype = P;}
 
    void clear_transfer_post_process()
      {m_vspProlongationPostProcess.clear(); m_vspRestrictionPostProcess.clear();}
 
    void add_prolongation_post_process(SmartPtr<ITransferPostProcess<TDomain, TAlgebra> > PP)
      {m_vspProlongationPostProcess.push_back(PP);}
 
    void add_restriction_post_process(SmartPtr<ITransferPostProcess<TDomain, TAlgebra> > PP)
      {m_vspRestrictionPostProcess.push_back(PP);}
 
  //  Linear Solver interface methods
 
    virtual const char* name() const {return "Geometric MultiGrid";}
 
    virtual bool supports_parallel() const
    {
      if(!m_spPreSmootherPrototype->supports_parallel())
        return false;
      if(!m_spPostSmootherPrototype->supports_parallel())
        return false;
      return true;
    }
 
    virtual std::string config_string() const;
 
    virtual bool init(SmartPtr<ILinearOperator<vector_type> > J, const vector_type& u);
 
    virtual bool init(SmartPtr<ILinearOperator<vector_type> > L);
 
 
    virtual void ignore_init_for_base_solver(bool ignore);
    virtual bool ignore_init_for_base_solver() const;
    void set_matrix_structure_is_const(bool b) {m_bMatrixStructureIsConst = b;}
 
    void force_reinit();
 
    virtual bool apply(vector_type& c, const vector_type& d);
 
    virtual bool apply_update_defect(vector_type& c, vector_type& d);
 
    SmartPtr<ILinearIterator<vector_type> > clone();
 
    ~AssembledMultiGridCycle();
 
  protected:
    void lmgc(int lev, int cycleType);
 
  //  The methods in this section rely on each other and should be called in sequence
    void presmooth_and_restriction(int lev);
 
    void prolongation_and_postsmooth(int lev);
 
    void base_solve(int lev);
  //  end of section
 
    void init_level_memory(int baseLev, int topLev);
 
    void init();
 
    void init_smoother();
 
    void assemble_level_operator();
    void init_rap_operator();
 
    void init_base_solver();
 
    void init_transfer();
 
    void init_projection();
 
    void assemble_rim_cpl(const vector_type* u);
    void init_rap_rim_cpl();
 
  protected:
    ConstSmartPtr<matrix_type> m_spSurfaceMat;
 
    const vector_type* m_pSurfaceSol;
 
    SmartPtr<IAssemble<TAlgebra> > m_spAss;
 
    SmartPtr<ApproximationSpace<TDomain> > m_spApproxSpace;
 
    int m_topLev;
    int m_surfaceLev;
 
    int m_baseLev;
 
    int m_cycleType;
    static const int _V_ = 1;
    static const int _W_ = 2;
    static const int _F_ = -1;
 
    int m_numPreSmooth;
 
    int m_numPostSmooth;
 
    int m_LocalFullRefLevel;
 
    GridLevel::ViewType m_GridLevelType;
 
    bool m_bUseRAP;
 
    bool m_bSmoothOnSurfaceRim;
 
    bool m_bCommCompOverlap;
 
    RevisionCounter m_ApproxSpaceRevision;
 
    SmartPtr<ILinearIterator<vector_type> > m_spPreSmootherPrototype;
 
    SmartPtr<ILinearIterator<vector_type> > m_spPostSmootherPrototype;
 
    SmartPtr<ITransferOperator<TDomain, TAlgebra> > m_spProjectionPrototype;
 
    SmartPtr<ITransferOperator<TDomain, TAlgebra> > m_spProlongationPrototype;
 
    SmartPtr<ITransferOperator<TDomain, TAlgebra> > m_spRestrictionPrototype;
 
    std::vector<SmartPtr<ITransferPostProcess<TDomain, TAlgebra> > > m_vspProlongationPostProcess;
    std::vector<SmartPtr<ITransferPostProcess<TDomain, TAlgebra> > > m_vspRestrictionPostProcess;
 
    SmartPtr<ILinearOperatorInverse<vector_type> > m_spBaseSolver;
 
    // Storage for each grid level
 
    struct LevelIndex{
      LevelIndex() : index(-1), indexLower(-1), level(-1), levelLower(-1) {}
      size_t index, indexLower;
      int level, levelLower;
    };
    std::vector<LevelIndex> m_vSurfToLevelMap;
 
    template <typename TElem>
    void init_index_mappings();
    void init_index_mappings();
 
    struct SurfLevelMap{
      SurfLevelMap() : levIndex(-1), surfIndex(-1) {}
      SurfLevelMap(size_t levIndex_, size_t surfIndex_)
        : levIndex(levIndex_), surfIndex(surfIndex_) {}
      size_t levIndex, surfIndex;
    };
 
    struct LevData
    {
      SmartPtr<MatrixOperator<matrix_type, vector_type> > A;
 
      SmartPtr<ILinearIterator<vector_type> > PreSmoother;
      SmartPtr<ILinearIterator<vector_type> > PostSmoother;
 
      SmartPtr<ITransferOperator<TDomain, TAlgebra> > Projection;
 
      SmartPtr<ITransferOperator<TDomain, TAlgebra> > Prolongation;
      SmartPtr<ITransferOperator<TDomain, TAlgebra> > Restriction;
 
      SmartPtr<GF> sc, sd, st, t;
 
      std::vector<size_t> vMapPatchToGlobal;
 
      std::vector<size_t> vShadowing;
 
      std::vector<size_t> vSurfShadowing;
 
      std::vector<SurfLevelMap> vSurfLevelMap;
 
      matrix_type RimCpl_Fine_Coarse;
 
      matrix_type RimCpl_Coarse_Fine;
      
      int n_pre_calls, n_post_calls, n_base_calls, n_restr_calls, n_prolong_calls;
    };
 
    std::vector<SmartPtr<LevData> > m_vLevData;
 
    bool m_bGatheredBaseIfAmbiguous;
 
    bool m_bGatheredBaseUsed;
 
    bool m_ignoreInitForBaseSolver;
 
    bool m_bMatrixStructureIsConst;
 
    SmartPtr<MatrixOperator<matrix_type, vector_type> > spGatheredBaseMat;
 
    SmartPtr<GF> spGatheredBaseCorr;
 
    bool gathered_base_master() const;
 
    GF* m_pC;
 
    template <typename TElem>
    void init_noghost_to_ghost_mapping( std::vector<size_t>& vNoGhostToGhostMap,
                      ConstSmartPtr<DoFDistribution> spNoGhostDD,
                      ConstSmartPtr<DoFDistribution> spGhostDD);
    void init_noghost_to_ghost_mapping(int lev);
 
    void copy_ghost_to_noghost(SmartPtr<GF> spVecTo,
                               ConstSmartPtr<GF> spVecFrom,
                               const std::vector<size_t>& vMapPatchToGlobal);
 
    void copy_noghost_to_ghost(SmartPtr<GF> spVecTo,
                   ConstSmartPtr<GF> spVecFrom,
                   const std::vector<size_t>& vMapPatchToGlobal);
 
    void copy_noghost_to_ghost(SmartPtr<matrix_type> spMatTo,
                               ConstSmartPtr<matrix_type> spMatFrom,
                               const std::vector<size_t>& vMapPatchToGlobal);
 
#ifdef UG_PARALLEL
    pcl::InterfaceCommunicator<IndexLayout> m_Com;
#endif
 
  public:
 
    void set_debug(SmartPtr<GridFunctionDebugWriter<TDomain, TAlgebra> > spDebugWriter)
    {
      m_spDebugWriter = spDebugWriter;
    }
 
  protected:
 
    inline void write_debug(ConstSmartPtr<GF> spGF, std::string name, int cycleNo = -1);
    void write_debug(const GF& rGF, std::string name, int cycleNo = -1);
 
 
    void write_debug(const matrix_type& mat, std::string name,
                     const GridLevel& glTo, const GridLevel& glFrom);
    void write_debug(const matrix_type& mat, std::string name,
                     const GF& rTo, const GF& rFrom);
  
    void enter_debug_writer_section(GridLevel& orig_gl, const char * sec_name, int lev, int cycleNo = -1, int callNo = -1);
    
    void leave_debug_writer_section(GridLevel& orig_gl);
 
    void log_debug_data(int lvl, int cycleNo, std::string name);
 
    void mg_stats_defect(GF& gf, int lvl, typename mg_stats_type::Stage stage);
 
    SmartPtr<GridFunctionDebugWriter<TDomain, TAlgebra> > m_spDebugWriter;
 
    int m_dbgIterCnt;
 
    SmartPtr<mg_stats_type> m_mgstats;
};
 
// Selections
 
void SelectNonShadowsAdjacentToShadowsOnLevel(BoolMarker& sel,
                       const SurfaceView& surfView,
                       int level);
 
} // end namespace ug
 
// include implementation
#include "mg_solver_impl.hpp"
 
#endif /* __H__UG__LIB_DISC__MULTI_GRID_SOLVER__MG_SOLVER__ */