timeprovider.hh

#ifndef DUNE_FEM_TIMEPROVIDER_HH
#define DUNE_FEM_TIMEPROVIDER_HH
 
#include <cassert>
#include <limits>
#include <tuple>
 
#include <dune/fem/io/file/persistencemanager.hh>
#include <dune/fem/io/parameter.hh>
#include <dune/fem/space/common/dofmanager.hh>
 
#include <dune/fem/misc/mpimanager.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
    class TimeProviderBase : public AutoPersistentObject
    {
      typedef TimeProviderBase ThisType;
 
    public:
      inline TimeProviderBase ( const ParameterReader &parameter = Parameter::container() )
      : time_( parameter.getValue( "fem.timeprovider.starttime",
                                      static_cast<double>(0.0) ) ),
        timeStep_( 0 ),
        dt_( 0.0 ),
        invdt_( HUGE_VAL ),
        valid_( false ),
        dtEstimateValid_( false ),
        parameter_( parameter )
      {
        initTimeStepEstimate();
      }
 
      inline explicit TimeProviderBase ( const double startTime, const ParameterReader &parameter = Parameter::container() )
      : time_( startTime ),
        timeStep_( 0 ),
        dt_( 0.0 ),
        invdt_( HUGE_VAL ),
        valid_( false ),
        dtEstimateValid_( false ),
        parameter_( parameter )
      {
        initTimeStepEstimate();
      }
 
      inline virtual ~TimeProviderBase()
      {}
 
      void backup() const
      {
        std::tuple<const double&,const int&,const double&,const bool&,const double&>
          values(time_,timeStep_,dt_,valid_,dtEstimate_);
        PersistenceManager::backupValue("timeprovider",values);
      }
 
      void restore()
      {
        std::tuple<double&,int&,double&,bool&,double&>
          values(time_,timeStep_,dt_,valid_,dtEstimate_);
        PersistenceManager::restoreValue("timeprovider",values);
        dtEstimateValid_ = true;
        invdt_ = 1.0 / dt_;
      }
 
 
      TimeProviderBase ( const ThisType & ) = delete;
 
      ThisType &operator= ( const ThisType & ) = delete;
 
 
      inline double time () const
      {
        return time_;
      }
 
      inline int timeStep () const
      {
        assert( timeStepValid() );
        return timeStep_;
      }
 
      inline double deltaT () const
      {
        assert( timeStepValid() );
        return dt_;
      }
 
      inline double inverseDeltaT () const
      {
        assert( timeStepValid() );
        return invdt_;
      }
 
      inline double timeStepEstimate () const
      {
        return dtEstimate_;
      }
 
      inline void provideTimeStepEstimate ( const double dtEstimate )
      {
        dtEstimate_ = std::min( dtEstimate_, dtEstimate );
        dtEstimateValid_ = true;
      }
      inline void provideTimeStepUpperBound ( const double upperBound )
      {
        dtUpperBound_ = std::min( dtUpperBound_, upperBound );
        dtEstimateValid_ = true;
      }
 
      inline void invalidateTimeStep ()
      {
        valid_ = false;
      }
 
      inline bool timeStepValid () const
      {
        return valid_;
      }
 
    protected:
      double time_;
      int timeStep_;
      double dt_;
      double invdt_;
      bool valid_;
      bool dtEstimateValid_;
      double dtEstimate_;
      double dtUpperBound_;
      ParameterReader parameter_;
 
      void advance ()
      {
        if( timeStepValid() )
        {
          time_ += deltaT();
          ++timeStep_;
        }
      }
 
      void initTimeStepEstimate ()
      {
        dtEstimate_ = std::numeric_limits< double >::max();
        dtUpperBound_ = std::numeric_limits< double >::max();
        dtEstimateValid_ = false;
      }
    };
 
 
    template< class Communication = typename MPIManager::Communication >
    class FixedStepTimeProvider
    : public TimeProviderBase
    {
      typedef FixedStepTimeProvider< Communication > ThisType;
      typedef TimeProviderBase BaseType;
 
    public:
      typedef Communication CommunicationType;
 
      explicit FixedStepTimeProvider ( const double startTime, const double timeStepSize,
                                       const CommunicationType &comm,
                                       const ParameterReader &parameter = Parameter::container() )
        : BaseType( startTime, parameter ), comm_( comm )
      {
        dt_ = timeStepSize;
        initTimeStep();
      }
 
      explicit FixedStepTimeProvider ( const double startTime, const double timeStepSize,
                                       const ParameterReader &parameter = Parameter::container() )
        : BaseType( startTime, parameter ), comm_( MPIManager::comm() )
      {
        dt_ = timeStepSize;
        initTimeStep();
      }
 
      explicit FixedStepTimeProvider ( const ParameterReader &parameter = Parameter::container() )
        : BaseType( parameter.getValue< double>( "fem.timeprovider.starttime", 0.0 ), parameter ), comm_( MPIManager::comm() )
      {
        dt_ = parameter.getValidValue< double >("fem.timeprovider.fixedtimestep", [] ( double v ) { return v > 0.0;} );
        initTimeStep();
      }
 
      explicit FixedStepTimeProvider ( const CommunicationType &comm,
                                       const ParameterReader &parameter = Parameter::container() )
        : BaseType( parameter.getValue< double>( "fem.timeprovider.starttime", 0.0 ), parameter ), comm_( comm )
      {
        dt_ = parameter.getValidValue< double >("fem.timeprovider.fixedtimestep", [] ( double v ) { return v > 0.0;} );
        initTimeStep();
      }
      virtual ~FixedStepTimeProvider () {}
 
      FixedStepTimeProvider ( const ThisType & ) = delete;
      FixedStepTimeProvider ( ThisType && ) = delete;
      ThisType &operator= ( const ThisType & ) = delete;
      ThisType &operator= ( ThisType && ) = delete;
 
      void next ()
      {
        if( !timeStepValid() )
          DUNE_THROW( InvalidStateException, "Invalid Time Step in FixedStepTimeProvider" );
        advance();
        initTimeStep();
      }
 
    protected:
      using BaseType::advance;
      using BaseType::initTimeStepEstimate;
      using BaseType::time_;
      using BaseType::dt_;
      using BaseType::valid_;
 
      inline void initTimeStep ()
      {
        valid_ = true;
        initTimeStepEstimate();
      }
 
      const CommunicationType &comm_;
    };
 
 
    template< class Communication = typename MPIManager::Communication >
    class TimeProvider
    : public TimeProviderBase
    {
      typedef TimeProvider< Communication > ThisType;
      typedef TimeProviderBase BaseType;
 
    public:
      typedef Communication CommunicationType;
 
    protected:
 
      using BaseType::parameter_;
 
      inline double getCflFactor() const
      {
        return parameter_.getValidValue( "fem.timeprovider.factor", static_cast<double>(1.0),
            [] ( double val ) { return val > 0.0; } );
      }
 
      inline int getUpdateStep () const
      {
        return parameter_.getValidValue( "fem.timeprovider.updatestep", static_cast<int>(1),
            [] ( int step ) { return step > 0; } );
      }
 
    public:
      explicit TimeProvider ( const ParameterReader &parameter = Parameter::container() )
      : BaseType( parameter ),
        comm_( MPIManager::comm() ),
        cfl_( getCflFactor() ),
        updateStep_( getUpdateStep() ),
        counter_( updateStep_ )
      {}
 
      explicit TimeProvider ( const CommunicationType &comm, const ParameterReader &parameter = Parameter::container() )
      : BaseType( parameter ),
        comm_( comm ),
        cfl_( getCflFactor() ),
        updateStep_( getUpdateStep() ),
        counter_( updateStep_ )
      {}
 
      explicit TimeProvider ( const double startTime,
                              const CommunicationType &comm = MPIManager::comm() )
      : BaseType( startTime ),
        comm_( comm ),
        cfl_( getCflFactor() ),
        updateStep_( getUpdateStep() ),
        counter_( updateStep_ )
      {}
 
      TimeProvider ( const double startTime,
                     const double cfl,
                     const CommunicationType &comm = MPIManager::comm() )
      : BaseType( startTime ),
        comm_( comm ),
        cfl_( cfl ),
        updateStep_( 1 ),
        counter_( updateStep_ )
      {}
 
      virtual ~TimeProvider()
      {}
 
      TimeProvider ( const ThisType & ) = delete;
 
      ThisType &operator= ( const ThisType & ) = delete;
 
      inline void init ()
      {
        initTimeStep( dtEstimate_ );
      }
 
      inline void init ( const double timeStep )
      {
        initTimeStep( timeStep );
      }
 
      inline void next ()
      {
        assert( this->dtEstimateValid_ );
        advance();
        initTimeStep( dtEstimate_ );
      }
 
      inline void next ( const double timeStep )
      {
        advance();
        initTimeStep(timeStep);
      }
 
      inline double factor () const
      {
        return cfl_;
      }
 
    protected:
      using BaseType::advance;
      using BaseType::initTimeStepEstimate;
 
      void initTimeStep ( const double dtEstimate )
      {
        // increase counter
        ++counter_ ;
 
        if( counter_ >= updateStep_ )
        {
          // set timestep estimate
          dt_ = std::min(cfl_ * dtEstimate,dtUpperBound_);
          dt_ = comm_.min( dt_ );
          invdt_ = 1.0 / dt_;
          valid_ = (dt_ > 0.0);
          // reset counter
          counter_ = 0;
        }
 
        initTimeStepEstimate();
      }
 
    public:
      inline void restore ( const double time, const int timeStep )
      {
        time_ = time;
        timeStep_ = timeStep;
      }
 
      inline virtual void backup () const
      {
        BaseType::backup();
      }
 
      inline virtual void restore ()
      {
        BaseType::restore();
        const_cast< double & >( cfl_ ) = getCflFactor();
      }
 
    protected:
      using BaseType::dt_;
      using BaseType::invdt_;
      using BaseType::dtEstimate_;
      using BaseType::dtUpperBound_;
      using BaseType::valid_;
      using BaseType::timeStep_;
 
      const CommunicationType& comm_;
      const double cfl_;
      const int updateStep_;
      int counter_;
    };
 
 
 
    template< class Grid >
    class GridTimeProvider
    : public TimeProvider< typename Grid::Traits::Communication >
    {
      typedef GridTimeProvider< Grid > ThisType;
      typedef TimeProvider< typename Grid::Traits::Communication > BaseType;
 
      // type of DofManager for sequence number
      typedef DofManager < Grid > DofManagerType ;
 
    public:
      typedef typename Grid::Traits::Communication CommunicationType;
 
      explicit GridTimeProvider ( const Grid &grid )
      : BaseType( grid.comm() ),
        dm_( DofManagerType ::instance( grid ) ),
        sequence_( -1 )
      {}
 
      GridTimeProvider ( const double startTime,
                         const Grid &grid )
      : BaseType( startTime, grid.comm() ),
        dm_( DofManagerType ::instance( grid ) ),
        sequence_( -1 )
      {}
 
      GridTimeProvider ( const double startTime,
                         const double cfl,
                         const Grid &grid )
      : BaseType( startTime, cfl, grid.comm() ),
        dm_( DofManagerType ::instance( grid ) ),
        sequence_( -1 )
      {}
 
      virtual ~GridTimeProvider() {}
 
    protected:
      using BaseType :: counter_ ;
      using BaseType :: updateStep_ ;
 
      // this initTimeStep method also check the sequence number
      // in case the grid has changed due to adaptivity
      void initTimeStep ( const double dtEstimate )
      {
        const int currentSequence = dm_.sequence();
        // check sequence number
        if( sequence_ != currentSequence )
        {
          // if sequence number changed, update in any case
          counter_  = updateStep_ ;
          sequence_ = currentSequence ;
        }
 
        // call initTimeStep on base class
        BaseType :: initTimeStep( dtEstimate );
      }
 
      const DofManagerType& dm_;
      int sequence_ ;
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #ifndef DUNE_FEM_TIMEPROVIDER_HH