util.hh

Go to the documentation of this file.

/*
 * util.hh
 *
 *  Created on: May 3, 2015
 *      Author: sunayana_ghosh
 */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2015-2015 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
//In this file we implement some of the utility functions necessary for
//the classes WorkModel, ErrorEstimate and the main integrate function of the class EulerSDC.
 
#ifndef KASKADE_TIMESTEPPING_EULERSDC_UTIL_HH_
#define KASKADE_TIMESTEPPING_EULERSDC_UTIL_HH_
 
#include "norm.hh"
 
#include <cmath>          //std::abs
#include <functional>     //std::function
#include <vector>
 
#include <boost/timer/timer.hpp>
 
#include "dune/common/dynmatrix.hh"
#include "dune/common/dynvector.hh"
 
#include "fem/fixdune.hh"
#include "timestepping/sdc.hh"
 
namespace Kaskade {
//=================================================================================================
//  Utility function to find norm of a difference of vector of two Vectors
//=================================================================================================
 
template<class Vector, class Norm>
typename Vector::value_type normVecDiff(std::vector<Vector> const& y0, std::vector<Vector> const& y1, Norm norm);
 
 
//=================================================================================================
//  Implementation for single sdc iteration step based on explicit Euler method.
//=================================================================================================
  template<class Vector>
  void sdcExplicitEulerIterationStep(Kaskade::SDCTimeGrid const& grid,
                                     std::function<Vector(typename Vector::value_type, Vector const&)> rhsFunc,
                                     std::vector<Vector> & yi, std::vector<Vector> & dyi,
                                     bool verbose = false)
  {
    //extract time points from the grid
    auto const& pts = grid.points();
    //extracting values of the integration matrix corresponding to the time grid.
    auto const& integ = grid.integrationMatrix();
    //number of intervals
    int const n = pts.size() - 1;
 
    //including starting point for yi and dyi
    assert(yi.size() == n+1);
    assert(dyi.size() == n+1);
 
    //initialize all elements of dyi to zero.
    for (int i = 0; i <= n; ++i)
    {
      dyi[i] = Vector(0.0);
    }
 
    //declare Vector total to be used inside the for loop
    Vector total(0.0);
 
    //the SDC step for explicit Euler
    //perform n Euler steps
    boost::timer::cpu_timer timer;
    for (auto i = 1; i <= n; ++i)
    {
      //Initialize all the entries of total to zero
      total = Vector(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 (auto i = 0; i <= n; ++i)
    {
      yi[i] += dyi[i];
      if (verbose)
        std::cout << "yi[" << i << "] = " << yi[i] << "\n";
    }
 
  }
 
  //=================================================================================================
   //  Implementation for single inexact sdc iteration step based on explicit Euler method.
   //=================================================================================================
 
 
   template<class Vector, class RealVector>
   void inexactSDCExplicitEulerIterationStep(Kaskade::SDCTimeGrid const& grid,
                                             std::function<Vector(typename Vector::value_type, Vector const&)> rhsFunc,
                                             RealVector const& toleranceVector, typename Vector::value_type rho,
                                             Kaskade::NormType norm_t,
                                             std::vector<Vector> & yi, std::vector<Vector> & dyi,
                                             bool verbose = false)
     {
       //extract time points from the grid
       auto const& tpts = grid.points();
       //extracting values of the integration matrix corresponding to the time grid.
       auto const& integ = grid.integrationMatrix();
       //number of intervals
       unsigned int const n = tpts.size() - 1;
 
       //including starting point for yi and dyi
       assert(yi.size() == n+1);
       assert(dyi.size() == n+1);
 
       //initialize all elements of dyi to zero.
       for (auto i = 0u; i <= n; ++i)
       {
         dyi[i] = Vector(0.0);
       }
 
       //declare Vector total to be used inside the for loop
       Vector total(0.0);
 
       //convert the toleranceVector which is a vector of doubles to a vector of Vectors.
       std::vector<Vector> errVector(n+1, Vector(0.0));
       //depending on the normType fill in err vector. This step is taken since the tolerance vector is a vector of scalars
       switch(norm_t)
       {
       case Kaskade::NormType::ONE_NORM:
       {
         auto sz = dyi[0].size();
         for (auto i = 1u; i <= n; ++i)
           errVector[i] = Vector(toleranceVector[i]/sz);
         break;
       }
       case Kaskade::NormType::MAX_NORM:
       {
         for (auto i = 1u; i <= n; ++i)
           errVector[i] = Vector(toleranceVector[i]);
         break;
        }
       }
 
       //the perturbed SDC step for explicit Euler
       //perform n Euler steps
       boost::timer::cpu_timer timer;
       for (auto i = 1u; i <= n; ++i)
       {
         //Initialize all the entries of total to zero
         total = Vector(0.0);
 
         //computation of the Lagrange interpolation
         for (auto k = 0u; k <= n; ++k)
         {
           total += integ[i-1][k] * (rhsFunc(tpts[k], yi[k]) + errVector[k]);
         }
         //approximate correction computation
         dyi[i] = dyi[i-1] + (tpts[i] - tpts[i-1]) *((rhsFunc(tpts[i-1],yi[i-1]+dyi[i-1]) - rhsFunc(tpts[i-1],yi[i-1]))- (1 - rho)*errVector[i-1]) + total - yi[i] + yi[i-1];
       }
 
       //Compute the next iteration of yi.
       for (auto i = 0u; i <= n; ++i)
       {
         yi[i] += dyi[i];
         if (verbose)
           std::cout << "yi[" << i << "] = " << yi[i] << "\n";
       }
     }
 
 
  //===================================================================================================
  //   Implementation of the abstract base class SDCUtil for all utility methods depending on the norms.
  //===================================================================================================
 
  template <class Vector, class Norm>
  class SDCUtil
  {
  public:
    using field_type = typename Vector::value_type;
    using RealVector = Dune::DynamicVector<double>;
    using RealMatrix = Dune::DynamicMatrix<double>;
 
    //virtual function
    field_type sdcContractionFactor(std::vector<Vector> const& yPrev,
                                    std::vector<Vector> const& yCurrent,
                                    std::vector<Vector> const& yNext,
                                    Norm& norm);
 
    //pure virtual functions
    virtual RealVector const& computeAlphaVec(RealVector const& timePoints,
                                              RealMatrix const& integrationMatrix,
                                              std::vector<Vector> const& yPrev,
                                              std::vector<Vector> const& yCurrent,
                                              std::vector<Vector> const& yNext) = 0;
 
 
    //virtual destructor
    virtual ~SDCUtil(){}
  };
 
  //===================================================================================================
  //   Implementation of the derived class SDCUtilOneNorm w.r.t. 1-norm.
  //===================================================================================================
  template<class Vector, class Norm>
  class SDCUtilOneNorm : public SDCUtil<Vector, Norm>
  {
  public:
    using field_type = typename Vector::value_type;
    using RealVector = Dune::DynamicVector<double>;
    using RealMatrix = Dune::DynamicMatrix<double>;
 
    virtual RealVector const& computeAlphaVec(RealVector const& timePts,
                                              RealMatrix const& integrationMatrix,
                                              std::vector<Vector> const& yPrev,
                                              std::vector<Vector> const& yCurrent,
                                              std::vector<Vector> const& yNext)
    {
      //number of collocation points
      auto nCollocationPts = timePts.size();
      alphaVec = RealVector(nCollocationPts-1, 0.0);
 
      //create one norm object
      Kaskade::OneNorm<Vector> on;
      //compute sdc contraction rate
      auto rho = Kaskade::SDCUtil<Vector, Norm>::sdcContractionFactor(yPrev, yCurrent, yNext, on);
      for (auto i = 0u; i < nCollocationPts-1; ++i)
      {
        field_type total = 0.0;
        for (auto k = 0u; k < nCollocationPts-1; ++k)
          total += std::abs(integrationMatrix[k][i]);
        if (i != nCollocationPts-2)
          alphaVec[i] = (total + (timePts[i+1] - timePts[i]) * (1 + rho));
        else
          alphaVec[i] = total;
      }
      return alphaVec;
    }
 
 
    virtual ~SDCUtilOneNorm(){}
 
  private:
    RealVector alphaVec;
    RealMatrix alphaMat;
  };
 
  //===================================================================================================
  //   Implementation of the derived class SDCUtilMaxNorm w.r.t. max-norm.
  //===================================================================================================
  template<class Vector, class Norm>
    class SDCUtilMaxNorm : public SDCUtil<Vector, Norm>
    {
    public:
      using field_type = typename Vector::value_type;
      using RealVector = Dune::DynamicVector<double>;
      using RealMatrix = Dune::DynamicMatrix<double>;
 
 
      //TODO: Doc me
      virtual RealVector const& computeAlphaVec(RealVector const& timePts,
                                                RealMatrix const& integrationMatrix,
                                                std::vector<Vector> const& yPrev,
                                                std::vector<Vector> const& yCurrent,
                                                std::vector<Vector> const& yNext)
      {
        //number of collocation points.
        auto nCollocationPts = timePts.size();
        gamma = RealVector(nCollocationPts, 0.0);
 
        //initialize the matrix alpha with size (nCollocationPts - 1) x nCollocationPts
        auto alphaMat = computeAlphaMat(timePts, integrationMatrix, yPrev, yCurrent, yNext, nCollocationPts);
        //Compute the vector Gamma.
        for (int i = 0; i < nCollocationPts -1; ++i)
        {
          for (int j = 0; j < nCollocationPts; ++j)
          {
            if (alphaMat[i][j] > gamma[j])
              gamma[j] = alphaMat[i][j];
          }
        }
        return gamma;
       }
 
 
 
      virtual ~SDCUtilMaxNorm() {};
    private:
      RealVector gamma;
      RealMatrix alphaMat;
 
      RealMatrix const& computeAlphaMat(RealVector const& timePts,
                                        RealMatrix const& integrationMatrix,
                                        std::vector<Vector> const& yPrev,
                                        std::vector<Vector> const& yCurrent,
                                        std::vector<Vector> const& yNext,
                                        int nCollocationPts);
 
 
    };
 
 
  //=================================================================================================
  //   FUNCTION DEFINITION BEGINS FROM HERE
  //=================================================================================================
 
  //=================================================================================================
  //  Utility function to find norm of a difference of vector of two Vectors
  //=================================================================================================
 
  template<class Vector, class Norm>
  typename Vector::value_type normVecDiff(std::vector<Vector> const& y0, std::vector<Vector> const& y1, Norm norm)
   {
     //Initialize the all the entries of dy to 0.0
     std::vector<Vector> dy(y0.size(), Vector(0.0));
     for (auto i = 0u; i < y0.size(); ++i)
     {
       dy[i] = y1[i] - y0[i];
     }
     return norm.value(dy);
   }
 
  //===================================================================================================
  //   Implementation of the abstract base class SDCUtil for all utility methods depending on the norms.
  //===================================================================================================
 
  template <class Vector, class Norm>
  typename Vector::value_type SDCUtil<Vector, Norm>::sdcContractionFactor(std::vector<Vector> const& yPrev,
                                                                          std::vector<Vector> const& yCurrent,
                                                                          std::vector<Vector> const& yNext,
                                                                          Norm& norm)
   {
    size_t n = yPrev.size();
    //Compute the differences: yPrev - yCurrent and yCurrent - yNext
    std::vector<Vector> y01(n, Vector(0.0));
    std::vector<Vector> y12(n, Vector(0.0));
    for (auto i = 0u; i < n; ++i)
    {
      y01[i] = yCurrent[i] - yPrev[i] ;
      y12[i] = yNext[i] - yCurrent[i];
    }
    auto val01 = norm.value(y01);
    auto val12 = norm.value(y12);
    auto beta = val12/(val01 + val12);
    auto val = val12/val01;
    if (val > 1 || val != val) //check if val > 1 or val is NaN.
    {
      if (beta != beta) //check if beta = NaN
        val = 0.70;
      else
        val = beta;
    }
    return val;
   }
 
 
  //===================================================================================================
  //   Implementation of the private method computeAlphaMat for the derived class SDCUtilMaxNorm
  //===================================================================================================
  using RealVector = Dune::DynamicVector<double>;
  using RealMatrix = Dune::DynamicMatrix<double>;
 
  template<class Vector, class Norm>
  RealMatrix const& SDCUtilMaxNorm<Vector, Norm>::computeAlphaMat(RealVector const& timePts,
                                                                  RealMatrix const& integrationMatrix,
                                                                  std::vector<Vector> const& yPrev,
                                                                  std::vector<Vector> const& yCurrent,
                                                                  std::vector<Vector> const& yNext,
                                                                  int nCollocationPts)
    {
    //create max norm object
    Kaskade::MaxNorm<Vector> mn;
    //compute sdc contraction rate
    auto rho = Kaskade::SDCUtil<Vector, Norm>::sdcContractionFactor(yPrev, yCurrent, yNext, mn);
    alphaMat = RealMatrix(nCollocationPts-1, nCollocationPts, 0.0);
    for (auto i = 0; i < nCollocationPts-1; ++i)
    {
      for (auto k = 0; k < nCollocationPts; ++k)
      {
        if (i == k)
          alphaMat[i][k] = std::abs(integrationMatrix[i][k]) + (timePts[i+1] - timePts[i]) * (1 + rho);
        else
          alphaMat[i][k] = std::abs(integrationMatrix[i][k]);
      }
    }
    return alphaMat;
    }
 
 
}//end namespace Kaskade
 
 
 
 
#endif /* KASKADE_TIMESTEPPING_EULERSDC_UTIL_HH_ */