explicit.hh

#ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_EXPLICIT_HH
#define DUNE_FEM_SOLVER_RUNGEKUTTA_EXPLICIT_HH
 
//- system includes
#include <iostream>
#include <cmath>
#include <vector>
#include <cassert>
 
//- Dune includes
#include <dune/fem/io/parameter.hh>
#include <dune/fem/operator/common/spaceoperatorif.hh>
#include <dune/fem/solver/odesolverinterface.hh>
#include <dune/fem/solver/timeprovider.hh>
#include <dune/fem/solver/rungekutta/butchertable.hh>
 
namespace DuneODE
{
 
  using namespace Dune;
  using namespace Fem;
  using namespace std;
 
  template<class DestinationImp>
  class ExplicitRungeKuttaSolver
  : public OdeSolverInterface<DestinationImp>
  {
  public:
    typedef DestinationImp DestinationType;
    typedef SpaceOperatorInterface<DestinationImp> OperatorType;
    typedef typename DestinationType :: DiscreteFunctionSpaceType SpaceType;
 
    typedef typename OdeSolverInterface<DestinationImp> :: MonitorType MonitorType ;
 
    using OdeSolverInterface<DestinationImp> :: solve ;
  protected:
    SimpleButcherTable< double > defaultButcherTables( const int order ) const
    {
      switch( order )
      {
        case 1: return explicitEulerButcherTable();
        case 2: return tvd2ButcherTable();
        case 3: return tvd3ButcherTable();
        case 4: return rk4ButcherTable();
        case 5:
        case 6: return expl6ButcherTable();
 
        default:
               std::cerr<< "Warning: ExplicitRungeKutta of order "<< order << " not implemented, using order 6!" << std::endl;
               return expl6ButcherTable();
      }
    }
 
  public:
    ExplicitRungeKuttaSolver(OperatorType& op,
                             TimeProviderBase& tp,
                             const SimpleButcherTable< double >& butcherTable,
                             bool verbose )
      : A_( butcherTable.A() ),
        b_( butcherTable.b() ),
        c_( butcherTable.c() ),
        Upd(),
        op_(op),
        tp_(tp),
        ord_( butcherTable.order() ),
        stages_( butcherTable.stages() ),
        initialized_(false)
    {
      assert(ord_>0);
 
      // create update memory
      for (int i=0; i<stages_; ++i)
      {
        Upd.emplace_back( new DestinationType("URK",op_.space()) );
      }
      Upd.emplace_back(new DestinationType("Ustep",op_.space()) );
 
    }
 
    ExplicitRungeKuttaSolver(OperatorType& op,
                             TimeProviderBase& tp,
                             const int pord,
                             bool verbose )
      : ExplicitRungeKuttaSolver( op, tp, defaultButcherTables( pord ), verbose )
    {
    }
 
    ExplicitRungeKuttaSolver(OperatorType& op,
                             TimeProviderBase& tp,
                             const int pord,
                             const Dune::Fem::ParameterReader &parameter = Dune::Fem::Parameter::container() )
      : ExplicitRungeKuttaSolver( op, tp, pord, true )
    {}
 
    void initialize(const DestinationType& U0)
    {
      if( ! initialized_ )
      {
        // Compute Steps
        op_(U0, *(Upd[0]));
        initialized_ = true;
 
        // provide operators time step estimate
        tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
      }
    }
 
    void solve(DestinationType& U0, MonitorType& monitor )
    {
      // no information here
      monitor.reset();
 
      // initialize
      if( ! initialized_ )
      {
        DUNE_THROW(InvalidStateException,"ExplicitRungeKuttaSolver wasn't initialized before first call!");
      }
 
      // get cfl * timeStepEstimate
      const double dt = tp_.deltaT();
      // get time
      const double t = tp_.time();
 
      // set new time
      op_.setTime( t );
 
      // Compute Steps
      op_(U0, *(Upd[0]));
 
      // provide operators time step estimate
      tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
 
      for (int i=1; i<stages_; ++i)
      {
        (Upd[ord_])->assign(U0);
        for (int j=0; j<i ; ++j)
        {
          (Upd[ord_])->axpy((A_[i][j]*dt), *(Upd[j]));
        }
 
        // set new time
        op_.setTime( t + c_[i]*dt );
 
        // apply operator
        op_( *(Upd[ord_]), *(Upd[i]) );
 
        // provide operators time step estimate
        tp_.provideTimeStepEstimate( op_.timeStepEstimate() );
      }
 
      // Perform Update
      for (int j=0; j<stages_; ++j)
      {
        U0.axpy((b_[j]*dt), *(Upd[j]));
      }
    }
 
    void description(std::ostream& out) const
    {
      out << "ExplRungeKutta, steps: " << ord_
          //<< ", cfl: " << this->tp_.factor()
          << "\\\\" <<std::endl;
    }
 
  protected:
    // Butcher table A,b,c
    Dune::DynamicMatrix< double > A_;
    Dune::DynamicVector< double > b_;
    Dune::DynamicVector< double > c_;
 
    // stages of Runge-Kutta solver
    std::vector< std::unique_ptr< DestinationType > > Upd;
 
    // operator to solve for
    OperatorType& op_;
    // time provider
    TimeProviderBase& tp_;
 
    // order of RK solver
    const int ord_;
    // number of stages
    const int stages_;
    // init flag
    bool initialized_;
  };
 
} // namespace DuneODE
 
#endif // #ifndef DUNE_FEM_SOLVER_RUNGEKUTTA_EXPLICIT_HH