adaptationmanager.hh

#ifndef DUNE_FEM_ADAPTATIONMANAGER_HH
#define DUNE_FEM_ADAPTATIONMANAGER_HH
 
//- system includes
#include <type_traits>
#include <string>
 
//- local includes
#include <dune/common/timer.hh>
#include <dune/fem/space/common/dofmanager.hh>
#include <dune/fem/operator/common/objpointer.hh>
 
#include <dune/fem/misc/capabilities.hh>
#include <dune/fem/space/common/communicationmanager.hh>
#include <dune/fem/space/common/loadbalancer.hh>
#include <dune/fem/space/common/restrictprolonginterface.hh>
#include <dune/fem/storage/singletonlist.hh>
 
#include <dune/fem/space/common/adaptcallbackhandle.hh>
#include <dune/fem/io/parameter.hh>
#include <dune/fem/io/file/persistencemanager.hh>
#include <dune/fem/misc/mpimanager.hh>
 
#include <dune/fem/space/common/dataprojection/dataprojection.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
  class AdaptationManagerInterface : virtual public LoadBalancerInterface
  {
  public:
    AdaptationManagerInterface () : am_ (0) {}
 
    virtual ~AdaptationManagerInterface () {}
 
    virtual void adapt ()
    {
      //std::cout << "called AdaptationManagerInterface::adapt()" << std::endl;
      if(am_) am_->adapt();
      else
      {
        std::cerr << "WARNING: AdaptationManagerInterface::adapt: no adaptation manager assigned! \n";
      }
    }
 
    virtual bool adaptive () const
    {
      return (am_) ? (am_->adaptive()) : false;
    }
 
    virtual bool isCallBackAdaptation() const { return (am_) ? (am_->isCallBackAdaptation()) : false; }
 
    virtual const char * methodName() const
    {
      return (am_) ? (am_->methodName()) : "unknown method";
    }
 
    AdaptationManagerInterface & operator = (const AdaptationManagerInterface & am)
    {
      am_ = const_cast<AdaptationManagerInterface *> (&am);
      return (*this);
    }
 
    virtual bool loadBalance ()
    {
      return (am_) ? (am_->loadBalance()) : false;
    }
 
    virtual int balanceCounter () const
    {
      return (am_) ? (am_->balanceCounter()) : 0;
    }
 
    virtual double adaptationTime () const
    {
      return 0.0;
    }
 
  private:
    AdaptationManagerInterface* am_;
  };
 
  template <class GridType>
  class AdaptationMethod : virtual public AdaptationManagerInterface
  {
  public:
    enum AdaptationMethodType { none = 0, 
                                generic = 1, 
                                callback = 2 
                              };
 
    static const int defaultMethod = 1;
 
  public:
    AdaptationMethod ( const GridType &grid,
                       const ParameterReader &parameter = Parameter::container(),
                       const bool noOutput = false )
      : adaptationMethod_(generic)
    {
      const bool output = ( ! noOutput && Parameter :: verbose( Parameter ::parameterOutput ) && MPIManager::isMainThread() );
      int am = defaultMethod;
      const std::string methodNames [] = { "none", "generic", "callback" };
      am = parameter.getEnum("fem.adaptation.method", methodNames, am);
      init(am, output);
    }
  private:
    void init(int am,const bool output)
    {
      // chose adaptation method
      if(am == 2) adaptationMethod_ = callback;
      else if(am == 1) adaptationMethod_ = generic;
      else adaptationMethod_ = none;
 
      // for structured grid adaptation is disabled
      if( ! Capabilities::isLocallyAdaptive<GridType>::v )
      {
        adaptationMethod_ = none;
        if( output )
        {
          std::cerr << "WARNING: AdaptationMethod: adaptation disabled for structured grid!"<< std::endl;
        }
      }
 
      static const bool noHierarchy = std::is_same< typename GridType::LevelGridView,
                                                    typename GridType::LeafGridView > :: value;
      // for grids without hierarchy generic adaptation is disabled
      // this is for example the case for P4estGrid. In this case we
      // assume that Leaf and Level iterators (views) are the same
      if( noHierarchy && (adaptationMethod_ == generic) )
      {
        adaptationMethod_ = callback;
        if( output )
        {
          std::cerr << "WARNING: AdaptationMethod: hierarchy not implemented, switching from 'generic' to 'callback'!" << std::endl;
        }
      }
 
      if( output )
      {
        std::cout << "Created AdaptationMethod: adaptation method = " << methodName() << std::endl;
      }
    }
  public:
    virtual ~AdaptationMethod () {}
 
    virtual const char * methodName() const
    {
      switch (adaptationMethod_) {
        case generic: return "generic";
        case callback: return "callback";
        case none: return "no adaptation";
        default:  return "unknown method";
      }
    }
 
    virtual bool adaptive () const { return adaptationMethod_ != none; }
 
    virtual bool isCallBackAdaptation() const { return adaptationMethod_ == callback; }
 
  protected:
    AdaptationMethodType adaptationMethod_;
  };
 
  template <class GridType, class RestProlOperatorImp >
  class AdaptationManagerBase
  : public AdaptationMethod< GridType >,
    public ObjPointerStorage
  {
    typedef AdaptationMethod< GridType > BaseType;
    typedef typename BaseType :: AdaptationMethodType AdaptationMethodType;
 
    template <class AdaptManager, class GridImp, bool isGoodGrid>
    struct CallAdaptationMethod
    {
      template <class DofManagerImp, class RPOpImp>
      static void adapt(const AdaptManager& am, GridImp & grid,
                        DofManagerImp& dm , RPOpImp& rpop,
                        AdaptationMethodType adaptMethod)
      {
        // use generic adapt method
        if( adaptMethod == BaseType :: generic )
        {
          am.template genericAdapt<All_Partition> ();
          return ;
        }
 
        // use grid call back adapt method
        if( adaptMethod == BaseType :: callback )
        {
          // combine dof manager and restrict prolong operator
          typedef RestrictProlongWrapper< GridImp, DofManagerType, RPOpImp > RPType;
 
          // create new handle
          RPType restrictProlongHandle ( dm , rpop );
 
          // reserve memory
          restrictProlongHandle.initialize();
 
          // call grid adaptation
          grid.adapt( restrictProlongHandle );
 
          // do compress (if not already called)
          restrictProlongHandle.finalize();
          return ;
        }
      }
    };
 
    template <class AdaptManager, class GridImp>
    struct CallAdaptationMethod<AdaptManager,GridImp,false>
    {
      template <class DofManagerImp, class RPOpImp>
      static void adapt(const AdaptManager& am, GridImp & grid,
                        DofManagerImp& dm , RPOpImp& rpop,
                        AdaptationMethodType adaptMethod)
      {
        // use generic adapt method
        if(adaptMethod != BaseType :: none )
        {
          // use partition type All_Partition,
          // since we also need to iterate on ghosts
          // for wasChanged information gathering
          am.template genericAdapt<All_Partition> ();
          return ;
        }
      }
    };
 
    typedef AdaptationManagerBase<GridType,RestProlOperatorImp> ThisType;
 
    typedef DofManager< GridType > DofManagerType;
 
  public:
    typedef typename GridType :: Traits :: LocalIdSet LocalIdSet;
 
    AdaptationManagerBase ( GridType &grid, RestProlOperatorImp &rpOp, const ParameterReader &parameter = Parameter::container() )
    : BaseType( grid, parameter ),
      grid_( grid ),
      dm_( DofManagerType::instance( grid_ ) ),
      rpOp_( rpOp ),
      adaptTime_( 0.0 ),
      wasChanged_( false )
    {}
 
    virtual ~AdaptationManagerBase () {}
 
    RestProlOperatorImp & getRestProlOp ()
    {
      return rpOp_;
    }
 
    virtual void adapt ()
    {
      // only call in single thread mode
      if( ! Fem :: MPIManager :: singleThreadMode() )
      {
        assert( Fem :: MPIManager :: singleThreadMode() );
        DUNE_THROW(InvalidStateException,"AdaptationManagerBase::adapt: only call in single thread mode!");
      }
 
      // get stopwatch
      Dune::Timer timer;
 
      const bool supportsCallback = Capabilities :: supportsCallbackAdaptation< GridType > :: v;
      CallAdaptationMethod< ThisType, GridType, supportsCallback >
        :: adapt(*this,grid_,dm_,rpOp_,this->adaptationMethod_);
 
      // take time
      adaptTime_ = timer.elapsed();
    }
 
    virtual bool loadBalance ()
    {
      return false;
    }
 
    virtual int balanceCounter () const
    {
      return 0;
    }
 
    virtual double adaptationTime() const
    {
      return adaptTime_;
    }
 
  protected:
    static DofManagerType& getDofManager(const GridType& grid)
    {
      return DofManagerType :: instance( grid );
    }
 
  private:
    template <PartitionIteratorType pitype>
    void genericAdapt () const
    {
      // initialize restrict prolong operator (e.g. PetscRestrictProlong... )
      rpOp_.initialize();
 
      // call pre-adapt, returns true if at least
      // one element is marked for coarsening
      bool restr = grid_.preAdapt();
 
      // get macro grid view
      typedef typename GridType::LevelGridView MacroGridView;
      typedef typename MacroGridView :: template Codim<0>::
              template Partition<pitype> :: Iterator MacroIterator;
 
      // reset flag
      wasChanged_ = false ;
 
      if(restr)
      {
 
        // get macro grid view
        MacroGridView macroView = grid_.levelGridView( 0 );
 
        // make a hierarchical to insert all elements
        // that are father of elements that might be coarsened
 
        {
          // get macro iterator
          MacroIterator endit  = macroView.template end<0,pitype>  ();
          for(MacroIterator it = macroView.template begin<0,pitype>();
              it != endit; ++it )
          {
            hierarchicRestrict( *it , dm_.indexSetRestrictProlongNoResize() );
          }
        }
 
        // if at least one element was found for restriction
        if( wasChanged_ )
        {
          // now resize memory
          dm_.resizeForRestrict();
 
          // now project all data to fathers
          {
            // get macro iterator
            MacroIterator endit  = macroView.template end<0,pitype>  ();
            for(MacroIterator it = macroView.template begin<0,pitype>();
                it != endit; ++it )
            {
              hierarchicRestrict( *it , rpOp_ );
            }
          }
        }
      }
 
      // adapt grid due to preset markers
      // returns true if at least one element was refined
      const bool refined = grid_.adapt();
 
      // if coarsening or refinement was done
      // adjust sizes
      if( refined || restr )
      {
        // resizes the index sets (insert all new indices)
        // and resizes the memory
        dm_.resize();
      }
 
      // in case elements were created do prolongation
      if( refined )
      {
        // get macro grid view
        MacroGridView macroView = grid_.levelGridView( 0 );
 
        // make run through grid to project data
        MacroIterator endit  = macroView.template end<0,pitype>  ();
        for(MacroIterator it = macroView.template begin<0,pitype>();
            it != endit; ++it )
        {
          hierarchicProlong( *it , rpOp_ );
        }
      }
 
      // notifyGlobalChange make wasChanged equal on all cores
      if( dm_.notifyGlobalChange( wasChanged_ ) )
      {
        // compress index sets and data
        // this will increase the sequence counter
        dm_.compress();
      }
 
      // do cleanup
      grid_.postAdapt();
 
      // finalize restrict prolong operator (e.g. PetscRestrictProlong... )
      rpOp_.finalize();
    }
 
  private:
    template <class EntityType, class RestrictOperatorType  >
    bool hierarchicRestrict ( const EntityType& entity, RestrictOperatorType & restop ) const
    {
      if( ! entity.isLeaf() )
      {
        // true means we are going to restrict data
        bool doRestrict = true;
 
        // check partition type
        const bool isGhost = entity.partitionType() == GhostEntity ;
 
        // if the children have children then we have to go deeper
        const int childLevel = entity.level() + 1;
        typedef typename EntityType::HierarchicIterator HierarchicIterator;
 
        // check all children first
        {
          const HierarchicIterator endit = entity.hend( childLevel );
          for(HierarchicIterator it = entity.hbegin( childLevel ); it != endit; ++it)
          {
            doRestrict &= hierarchicRestrict( *it , restop );
          }
        }
 
        // if doRestrict is still true, restrict data
        if(doRestrict)
        {
          // we did at least one restriction
          wasChanged_ = true;
 
          // do not restrict the solution on ghosts, this will
          // fail, but we still need the wasChanged info, so simply
          // calling hierarchicRestrict on interior won't work either
          if( ! isGhost )
          {
            // true for first child, otherwise false
            bool initialize = true;
            const HierarchicIterator endit = entity.hend( childLevel );
            for(HierarchicIterator it = entity.hbegin( childLevel ); it != endit; ++it)
            {
              // restrict solution
              restop.restrictLocal( entity, *it , initialize);
              // reset initialize flag
              initialize = false;
            }
            restop.restrictFinalize(entity);
          }
        }
      }
 
      // if all children return mightBeCoarsened,
      // then doRestrict on father remains true
      return entity.mightVanish();
    }
 
    template <class EntityType, class ProlongOperatorType >
    void hierarchicProlong ( const EntityType &entity, ProlongOperatorType & prolop ) const
    {
      typedef typename EntityType::HierarchicIterator HierarchicIterator;
 
      // NOTE: initialize not working here
      // because we call hierarchically
 
      // first call on this element
      bool initialize = true;
 
      // check partition type
      const bool isGhost = entity.partitionType() == GhostEntity ;
 
      const int maxLevel = grid_.maxLevel();
      const HierarchicIterator endit = entity.hend( maxLevel );
      for( HierarchicIterator it = entity.hbegin( maxLevel ); it != endit; ++it )
      {
        // should only get here on non-leaf entities
        assert( !entity.isLeaf() );
 
        const EntityType & son = *it;
        if( son.isNew() )
        {
          // the grid was obviously changed if we get here
          wasChanged_ = true ;
 
          // do not prolong the solution on ghosts, this will
          // fail, but we still need the wasChanged info, so simply
          // calling hierarchicRestrict on interior won't work either
          if( ! isGhost )
          {
            EntityType vati = son.father();
            prolop.prolongLocal( vati , son , initialize );
            initialize = false;
          }
        }
      }
    }
 
  protected:
    GridType &grid_;
 
    DofManagerType &dm_;
 
    RestProlOperatorImp &rpOp_;
 
    double adaptTime_;
 
    mutable bool wasChanged_;
  };
 
  template <class KeyType, class ObjectType>
  struct AdaptationManagerReferenceFactory
  {
    static ObjectType* createObject(const KeyType& key)
    {
      return new ObjectType(0);
    }
    static void deleteObject(ObjectType* obj)
    {
      delete obj;
    }
  };
 
  template <class GridType, class RestProlOperatorImp>
  class AdaptationManager :
    public AdaptationManagerBase<GridType,RestProlOperatorImp> ,
    public LoadBalancer<GridType>,
    public AutoPersistentObject
  {
    // type of key
    typedef const GridType* KeyType;
    // object type
    typedef size_t ObjectType;
    // type of factory
    typedef AdaptationManagerReferenceFactory<KeyType, ObjectType>  FactoryType;
 
    // type of singleton list
    typedef SingletonList< KeyType, ObjectType, FactoryType > ProviderType;
 
    typedef AdaptationManagerBase<GridType,RestProlOperatorImp> BaseType;
    typedef LoadBalancer<GridType> Base2Type;
 
    using BaseType :: rpOp_;
 
    // reference counter to ensure only one instance per grid exists
    ObjectType& referenceCounter_;
 
    // do not copy
    AdaptationManager(const AdaptationManager&);
 
  public:
    AdaptationManager ( GridType &grid, RestProlOperatorImp &rpOp, int balanceCounter, const ParameterReader &parameter = Parameter::container() )
      : BaseType(grid,rpOp, parameter)
      , Base2Type( grid, rpOp )
      , referenceCounter_( ProviderType :: getObject( &grid ) )
      , balanceStep_( parameter.getValue< int >( "fem.loadbalancing.step", 1 ) )
      , balanceCounter_( balanceCounter )
    {
      if( ++referenceCounter_ > 1 )
        DUNE_THROW(InvalidStateException,"Only one instance of AdaptationManager allowed per grid instance");
      if( Parameter::verbose( Parameter::parameterOutput ) )
        std::cout << "Created LoadBalancer: balanceStep = " << balanceStep_ << std::endl;
    }
 
    AdaptationManager ( GridType &grid, RestProlOperatorImp &rpOp, const ParameterReader &parameter = Parameter::container() )
      : AdaptationManager( grid, rpOp, 0, parameter )
    {
    }
 
    ~AdaptationManager()
    {
      -- referenceCounter_;
      ProviderType :: removeObject( referenceCounter_ );
    }
 
    virtual bool loadBalance ()
    {
      // same as for the adapt method
      rpOp_.initialize () ;
 
      // call load balance
      const bool result = Base2Type :: loadBalance( );
 
      // finalize rp object (mostly RestrictProlongDefault for PetscDF)
      rpOp_.finalize () ;
      return result ;
    }
 
    virtual double loadBalanceTime() const
    {
      return Base2Type::loadBalanceTime();
    }
 
    virtual void adapt ()
    {
      // adapt grid
      BaseType :: adapt ();
 
      // if adaptation is enabled
      if( this->adaptive() && (balanceStep_ > 0) )
      {
        // if balance counter has readed balanceStep do load balance
        const bool callBalance = (++balanceCounter_ >= balanceStep_);
 
#ifndef NDEBUG
        // make sure load balance is called on every process
        int willCall = (callBalance) ? 1 : 0;
        const int iCall = willCall;
 
        // send info from rank 0 to all other
        Base2Type::grid_.comm().broadcast(&willCall, 1 , 0);
 
        assert( willCall == iCall );
#endif
 
        if( callBalance )
        {
          // balance work load and restore consistency in the data
          loadBalance();
          balanceCounter_ = 0;
        }
        else
        {
          // only restore consistency in the data
          Base2Type::communicate();
        }
      }
    }
 
    int balanceCounter () const { return balanceCounter_; }
 
    void backup() const
    {
      std::tuple<const int& > value( balanceCounter_ );
      PersistenceManager::backupValue("loadbalancer",value);
    }
 
    void restore()
    {
      std::tuple< int& > value( balanceCounter_ );
      PersistenceManager::restoreValue("loadbalancer",value);
    }
 
  private:
    // call loadBalance ervery balanceStep_ step
    const int balanceStep_ ;
    // count actual balance call
    int balanceCounter_;
  };
 
  namespace hpDG
  {
 
    // AdaptationManager
    // -----------------
 
    template< class DiscreteFunctionSpace, class DataProjection >
    class AdaptationManager
    : public Dune::Fem::AdaptationManagerInterface
    {
      using ThisType = AdaptationManager< DiscreteFunctionSpace, DataProjection >;
 
    public:
      using DiscreteFunctionSpaceType = DiscreteFunctionSpace;
      using DataProjectionType = DataProjection;
 
    private:
      using GridType = typename DiscreteFunctionSpaceType::GridType;
      using DofManagerType = DofManager< GridType >;
 
      class DataProjectionWrapper;
 
    public:
      explicit AdaptationManager ( DiscreteFunctionSpaceType &space, DataProjectionType &&dataProjection )
        : space_( space ),
          dataProjection_( std::forward< DataProjectionType >( dataProjection ) ),
          dofManager_( DofManagerType::instance( space.gridPart().grid() ) ),
          commList_( dataProjection_ ),
          time_( 0. )
      {}
 
      AdaptationManager ( const ThisType & ) = delete;
 
      ThisType &operator= ( const ThisType & ) = delete;
 
      bool adaptive () const { return true; }
 
      void adapt ()
      {
        // only call in single thread mode
        if( ! Fem :: MPIManager :: singleThreadMode() )
        {
          assert( Fem :: MPIManager :: singleThreadMode() );
          DUNE_THROW(InvalidStateException,"AdaptationManager::adapt: only call in single thread mode!");
        }
 
        Dune::Timer timer;
 
        DataProjectionWrapper wrapper( dataProjection_, dofManager() );
        space().adapt( wrapper );
 
        if( dofManager().notifyGlobalChange( static_cast< bool >( wrapper ) ) )
          dofManager().compress();
 
        commList_.exchange();
 
        time_ = timer.elapsed();
      }
 
      const char *methodName () const { return "discrete function space adaptation"; }
 
      double adaptationTime () const { return time_; }
 
      bool loadBalance () { return false; }
 
      int balanceCounter () const { return 0; }
 
      double loadBalanceTime () const { return 0.; }
 
      DataProjection& dataProjection() { return dataProjection_; }
    private:
      DiscreteFunctionSpaceType &space () { return space_.get(); }
 
      const DiscreteFunctionSpaceType &space () const { return space_.get(); }
 
      DofManagerType &dofManager () { return dofManager_.get(); }
 
      const DofManagerType &dofManager () const { return dofManager_.get(); }
 
      std::reference_wrapper< DiscreteFunctionSpaceType > space_;
      DataProjectionType dataProjection_;
      std::reference_wrapper< DofManagerType > dofManager_;
      mutable CommunicationManagerList commList_;
      double time_;
    };
 
    // AdaptationManager::DataProjectionWrapper
    // ----------------------------------------
 
    template< class DiscreteFunctionSpace, class DataProjection >
    class AdaptationManager< DiscreteFunctionSpace, DataProjection >::DataProjectionWrapper
    : public Dune::Fem::hpDG::DataProjection< DiscreteFunctionSpace, DataProjectionWrapper >
    {
      using BaseType = Dune::Fem::hpDG::DataProjection< DiscreteFunctionSpace, DataProjectionWrapper >;
 
    public:
      using BasisFunctionSetType = typename BaseType::BasisFunctionSetType;
      using EntityType = typename BaseType::EntityType;
 
      DataProjectionWrapper ( DataProjectionType &dataProjection, DofManagerType &dofManager )
        : dataProjection_( dataProjection ),
          dofManager_( dofManager ),
          modified_( false )
      {}
 
      void operator() ( const EntityType &entity,
                        const BasisFunctionSetType &prior,
                        const BasisFunctionSetType &present,
                        const std::vector< std::size_t > &origin,
                        const std::vector< std::size_t > &destination )
      {
        dofManager_.get().resizeMemory();
        dataProjection_.get()( entity, prior, present, origin, destination );
        modified_ = true;
      }
 
      template <class TemporaryStorage>
      void operator() ( TemporaryStorage& tmp )
      {
        dataProjection_.get()( tmp );
        modified_ = true;
      }
 
      explicit operator bool () const
      {
        return modified_;
      }
 
    private:
      std::reference_wrapper< DataProjectionType > dataProjection_;
      std::reference_wrapper< DofManagerType > dofManager_;
      bool modified_;
    };
 
  } // namespace hpDG
 
 
 
  struct GlobalRefine
  {
    template <class GridType>
    static void apply(GridType& grid, const int step)
    {
      typedef DofManager< GridType > DofManagerType;
      DofManagerType& dm = DofManagerType :: instance(grid);
      grid.globalRefine(step);
      grid.loadBalance();
      dm.resize();
      dm.compress();
    }
  };
  struct LocalRefine
  {
    template <class GridType>
    static void apply(GridType& grid)
    {
      typedef DofManager< GridType > DofManagerType;
      DofManagerType& dm = DofManagerType :: instance(grid);
      grid.preAdapt();
      grid.adapt();
      grid.postAdapt();
      grid.loadBalance();
      dm.resize();
      dm.compress();
    }
  };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_ADAPTMANAGER_HH