sdc_euler.hh

Go to the documentation of this file.

#ifndef SDC_EULER_HH
#define SDC_EULER_HH
 
#include <algorithm>
#include <functional>
#include <vector>
#include "timestepping/sdc.hh"
#include "fem/fixdune.hh"
 
namespace Kaskade
{
  //FIXME:Adapt and move the norms to Kaskade::LinAlg::Norm
  template<class Vector>
  typename Vector::field_type oneNorm(std::vector<Vector> const& values)
  {
    return std::accumulate(values.begin(),values.end(),0.0,[](double init, Vector const& v)
    {
        return init + v.one_norm();
      });
//    auto norm = 0.0;
 
//    for (auto i = 0u; i < values.size(); ++i)
//      //one_norm: public member function of Dune::DenseVector
//      norm += values[i].one_norm();
 
//    return norm;
  }
 
  template<class Vector>
  typename Vector::field_type maxNorm(std::vector<Vector> const& values)
  {
    typename Vector::field_type max = values[0].infinity_norm();
 
    for (auto i = 1u; i < values.size(); ++i)
      //infinity_norm: public member function of Dune::DenseVector
      max = std::max(max, values[i].infinity_norm());
 
    return max;
  }
 
 
  //FIXME: make the sdcExEulerIterationStep generic, instead of std::vector<Vector> replace with generic class Vectors.
  //single SDC iteration step, which returns
  template <class Vector>
  typename Vector::field_type sdcExEulerIterationStep(Kaskade::SDCTimeGrid const& grid, std::vector<Vector> & yi, std::vector<Vector> & dyi,
                                                      std::function<Vector(typename Vector::field_type, Vector const&)> rhsFunc,
                                                      std::function<typename Vector::field_type(std::vector<Vector> const&)> normFunc,
                                                      bool verbose = false)
  {
    //extracting time points from the time grid
    auto const& pts = grid.points();
    //extracting values of the integration matrix corresponding to the time grid
    auto const& integ = grid.integrationMatrix();
    //number of subintervals
    int const n = pts.size() - 1;
 
    //including starting point for yi and dyi
    assert(yi.size() == n+1);
    assert(dyi.size() == n+1);
 
    //initialize the correction at starting point to zero, where dyi[0] is a vector in general
    //and every component of dyi is set to 0.0 by assignment below.
    dyi[0] = 0.0;
 
    //declare Vector total to be used inside the for loop
    Vector total(0.0);
    typename Vector::field_type norm = 0.0;
 
    //The Sdc step for explicit Euler
    //perform n Euler steps
    boost::timer::cpu_timer timer;
    for (int i = 1; i <= n; ++i)
      {
        //Initialize all entries of Vector total to 0.
        for (size_t t = 0; t < total.size(); ++t)
          total[t] = 0.0;
        //Computation of the Lagrange interpolation
        for (int k = 0; k < n; ++k)
          total += integ[i-1][k] * rhsFunc(pts[k], yi[k]);
        //approximate correction computation.
        dyi[i] = dyi[i-1] + (pts[i] - pts[i-1]) *(rhsFunc(pts[i-1],yi[i-1]+dyi[i-1]) - rhsFunc(pts[i-1],yi[i-1])) + total - yi[i] + yi[i-1];
 
      }
 
    //compute the next iteration of yi.
    for (int i = 0; i <= n; ++i)
      {
        yi[i] += dyi[i];
        if(verbose) std::cout << "yi[" << i << "] = " << yi[i] << "\n";
      }
 
    norm = normFunc(dyi);
    return norm;
  }
 
  template <class Vector, class TimeGrid=LobattoTimeGrid>
  class EulerSDC
  {
  public:
    EulerSDC(double t0, double t1, size_t nIntervals, bool verbose_=false)
      : tbegin(t0), tend(t1), verbose(verbose_), grid(nIntervals,tbegin,tend)
    {}
 
    EulerSDC(double t0, double t1, size_t nIntervals, size_t maxSDCIter, bool verbose_=false)
      : tbegin(t0), tend(t1), maxSDCiterations(maxSDCIter), verbose(verbose_), grid(nIntervals,tbegin,tend)
    {}
 
    EulerSDC(double t0, double t1, size_t nIntervals, size_t maxSDCIter, double tol, bool verbose_=false)
      : tbegin(t0), tend(t1), maxSDCiterations(maxSDCIter), tolerance(tol), verbose(verbose_), grid(nIntervals, tbegin, tend)
    {}
 
 
    void setMaxSDCIterations(size_t maxIter)
    {
      maxSDCiterations = maxIter;
    }
 
    void setTolerance(double tau)
    {
      tolerance = tau;
    }
 
 
    //integrate function to carry out SDC iterations based on explicit Euler method based on the number of iterations(maxSDCiterations) given by user.
    Vector integrate( Vector const& initialValue, std::function<Vector (double,Vector const&)> rightHandSide,
                      std::function<typename Vector::field_type(std::vector<Vector const&>)> normFunc)
    {
      //initialize the vector yi
      std::vector<Vector> yi(grid.points().size(),initialValue), dyi(grid.points().size(),Vector(0.0));
 
      //print the entries of the vector yi
      if(verbose)
        {
          std::cout << "print entries to the vector yi" << std::endl;
          for (int i = 0; i < yi.size(); ++i)
            std::cout << "[ " << yi[i] << " ]" << std::endl;
        }
 
      for (int i = 0; i < maxSDCiterations; ++i)
        {
          if(verbose) std::cout << "Iteration : " << i << "\n";
          auto stepNorm = sdcExEulerIterationStep(grid, yi, dyi, rightHandSide, normFunc, verbose);
          if(verbose) std::cout << "norm = " << stepNorm << std::endl;
        }
 
      return yi.back();
    }
 
    //integrate function to carry out SDC iterations based on explicit Euler method based on the tolerance given by the user.
    //In case the tolerance is not reached within the default value of maxSDCiterations, the last value of yi at maxSDCiterations is returned.
    Vector integrateTOL( Vector const& initialValue, std::function<Vector(double,Vector const&)> rightHandSide,
                         std::function<typename Vector::field_type(std::vector<Vector> const&)> normFunc)
    {
      //initialize the vector yi with initial value and the correction vector dyi to 0 at all grid points.
      std::vector<Vector> yi(grid.points().size(),initialValue), dyi(grid.points().size(),Vector(0.0));
      //initialize loop variable to keep track of when loop exceeds maxSDCiterations
      int loop = 1;
 
      //print the entries of the vector yi
      if(verbose)
        {
          std::cout << "print entries to the vector yi" << std::endl;
          for (int i = 0; i < yi.size(); ++i)
            std::cout << "[ " << yi[i] << " ]" << std::endl;
        }
 
      //while loop to carry out sdc iterations till the given tolerance is reached failing which it loops till maxSDCiterations is reached.
      auto stepNorm = sdcExEulerIterationStep(grid, yi, dyi, rightHandSide, normFunc, verbose);
      while (stepNorm > tolerance  && loop < maxSDCiterations)
        {
          //std::cout << "step norm = " << stepNorm << "\n";
          if(verbose) std::cout << "Iteration: " << loop << "\n";
          stepNorm = sdcExEulerIterationStep(grid, yi, dyi, rightHandSide, normFunc, verbose);
          if(verbose) std::cout << "norm = " << stepNorm << std::endl;
          ++loop;
        }
      return yi.back();
    }
 
  private:
    double tbegin = 0, tend = 1;
    size_t maxSDCiterations = 100;
    double tolerance = 0.1;
    bool verbose = false;
    TimeGrid grid;
  };
 
}
 
#endif // SDC_EULER_HH