sdc.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2012-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef SDC_HH
#define SDC_HH
 
#include <boost/timer/timer.hpp>
 
#include "dune/grid/config.h"
#include "dune/istl/operators.hh"
#include "dune/istl/solvers.hh"
 
#include "dune/common/dynvector.hh"
 
#include "linalg/dynamicMatrix.hh"
#include "utilities/timing.hh"
 
#include <cassert>
 
 
namespace Kaskade {
 
  class SDCTimeGrid
  {
  public:
    typedef Dune::DynamicVector<double> RealVector;
    
    typedef DynamicMatrix<double> RealMatrix;
    
    virtual RealVector const& points() const = 0;
    
    double point(int k) const
    {
      return points()[k];
    }
    
    virtual RealMatrix const& integrationMatrix() const = 0;
    
    virtual RealMatrix const& differentiationMatrix() const = 0;
    
    virtual void refine(RealMatrix& p) = 0;
    
    virtual RealMatrix interpolate(RealVector const& x) const = 0;
  };
  
  // ---------------------------------------------------------------------------------------------------
 
  class SDCGridImplementationBase: public SDCTimeGrid
  {
  public:
   
    virtual RealVector const& points() const 
    {
      return pts;
    }
    
    virtual RealMatrix const& integrationMatrix() const
    {
      return integ;
    }
    
    virtual RealMatrix const& differentiationMatrix() const
    {
      return diff;
    }
    
    
  protected:
    RealVector pts;
    RealMatrix integ;
    RealMatrix diff;
  };
  
  // ---------------------------------------------------------------------------------------------------
 
  void eulerIntegrationMatrix(SDCTimeGrid const& grid, SDCTimeGrid::RealMatrix& Shat);
  
  SDCTimeGrid::RealMatrix eulerIntegrationMatrix(SDCTimeGrid const& grid);
 
  void luIntegrationMatrix(SDCTimeGrid const& grid, SDCTimeGrid::RealMatrix& Shat);
 
  SDCTimeGrid::RealMatrix luIntegrationMatrix(SDCTimeGrid const& grid);
 
  SDCTimeGrid::RealMatrix singleEulerIntegrationMatrix(SDCTimeGrid const& grid);
 
  SDCTimeGrid::RealMatrix SDIRKIntegrationMatrix(SDCTimeGrid const& grid);
 
 
  
//   /**
//    * \ingroup sdc
//    * \brief Types of weight functions for optimized integration matrices.
//    */
//   enum IntegrationMatrixOptimizationWeight
//   {
//     /// optimization for a contraction factor uniform on the real spectrum
//     OW_FLAT,
//     /// optimization for a contraction factor uniform on the real negative spectrum, but for SDIRK
//     OW_FLAT_SDIRK
//   };
//
//   /**
//    * \ingroup sdc
//    * \brief Triangular approximate differentiation and integration matrix \f$ \hat D, \hat S \f$ corresponding to specialized multi-step integrators
//    *
//    * \param grid the collocation time grid
//    * \param k    the sweep number
//    * \param w    the type of weight function for the optimization
//    *
//    * This returns one of a set of precomputed, optimized integration matrices, if available for the given grid.
//    * If no precomputed matrix is available, a lookup exception is thrown.
//    */
//   std::pair<SDCTimeGrid::RealMatrix,SDCTimeGrid::RealMatrix>
//   optimizedMatrices(SDCTimeGrid const& grid, int k, IntegrationMatrixOptimizationWeight w=OW_FLAT);
//
//   /**
//    * \ingroup sdc
//    * \brief Triangular approximate differentiation and integration matrix \f$ \hat D, \hat S \f$ corresponding
//    *        to specialized multi-step integrators
//    *
//    * \param grid     the time grid
//    * \param k        the sweep number
//    * \param Df       the differentiation matrix returned in case no precomputed optimized matrix is available
//    * \param Sf       the integration matrix returned in case no precomputed optimized matrix is available
//    * \param w        the type of weight function for the optimization
//    *
//    * This returns one of a set of precomputed, optimized integration matrices, if available for the given grid.
//    * If no precomputed matrix is available, the provided fallback matrix is returned.
//    * This is a convenience overload.
//    */
//   std::pair<SDCTimeGrid::RealMatrix,SDCTimeGrid::RealMatrix>
//   optimizedIntegrationMatrix(SDCTimeGrid const& grid, int k, SDCTimeGrid::RealMatrix const& Df, SDCTimeGrid::RealMatrix const& Sf,
//                              IntegrationMatrixOptimizationWeight w=OW_FLAT);
 
  
  // ---------------------------------------------------------------------------------------------------
  
  class LobattoTimeGrid: public SDCTimeGrid
  {
  public:
    LobattoTimeGrid(int n, double a, double b);
    
    virtual RealVector const& points() const 
    {
      return pts;
    }
    
    virtual RealMatrix const& integrationMatrix() const
    {
      return integ;
    }
    
    virtual RealMatrix const& differentiationMatrix() const
    {
      return diff;
    }
    
    virtual void refine(RealMatrix& p);
    
    virtual RealMatrix interpolate(RealVector const& x) const;
 
  private:
    RealVector pts;
    RealMatrix integ;
    RealMatrix diff;
  };
  
  class RadauTimeGrid: public SDCGridImplementationBase
  {
  public:
    RadauTimeGrid(int n, double a, double b);
    
    virtual void refine(RealMatrix& p);
    
    virtual RealMatrix interpolate(RealVector const& x) const;
 
  private:
    void computeMatrices();
  };
  
  
  class GaussTimeGrid: public SDCGridImplementationBase
  {
  public:
    GaussTimeGrid(int n, double a, double b);
    
    virtual void refine(RealMatrix& p);
    
    virtual RealMatrix interpolate(RealVector const& x) const;
 
  private:
    void computeMatrices();
    double b; // store the end point of the interval
  };
  
  
  // -----------------------------------------------------------------------------------------------
 
  namespace SdcDetail
  {
    template <class Matrix, class Vectors>
    std::vector<typename Vectors::value_type> computeMdu(Matrix const& M, Vectors const& u)
    {
      int const n = u.size()-1;      // number of subintervals
 
      using Vector = typename Vectors::value_type;
 
      Vector tmp = u[0];
      std::vector<Vector> Mdu(n,tmp);
 
      for (int i=0; i<n; ++i)
      {
        // compute M (u_{i}-u_{i+1})  TODO: loop fusion
        tmp = u[i];
        tmp.axpy(-1.0,u[i+1]);
        M.mv(tmp,Mdu[i]);
      }
 
      return Mdu;
    }
 
 
    // Construct J = M - s*(A+f_u), where f_u may have a sparsity pattern that is a subset
    // of the sparsity pattern of J, M, and A (which need to be equal.
    // f_u may be nullptr, in which case f_u is completely ignored.
    template <class Matrix, class ReactionDerivative>
    void buildEulerStepMatrix(Matrix const& M, double s, Matrix const& A,
                              ReactionDerivative const& f_u, Matrix& J)
    {
      for (size_t row=0; row<J.N(); ++row)
      {
        auto colJ = J[row].begin();
        auto end = J[row].end();
        auto colM = M[row].begin();
        auto colA = A[row].begin();
 
        if constexpr (std::is_same_v<ReactionDerivative,std::nullptr_t>)  // no reaction derivative
        {                                                                 // to be considered
          while (colJ != end)
          {
            *colJ = *colM - s * *colA;
            ++colJ; ++colM; ++colA;
          }
        }
        else
        {
          auto colF = f_u[row].begin();
          auto endF = f_u[row].end();
 
          while (colJ != end)
          {
            *colJ = *colM - s * *colA;
 
            if (colF != endF && colJ.index() == colF.index()) // f_u can have subset of sparsity pattern
            {
              // guarantee M - Shat*fu is nonnegative -- reduce by at most 50%
              *colJ -= std::min(0.5 * *colM, s * *colF);
              ++colF;
            }
            ++colJ; ++colM; ++colA;
          }
        }
 
        assert(colM==M[row].end());
        assert(colA==A[row].end());
      }
    }
 
    // Checks if all diagonal values (that lead to implicit behavior)
    // are the same, such that the same system has to be solved for all stages.
    bool isSinglyDiagonal(SDCTimeGrid::RealMatrix const& Shat, double relError=1e-6);
 
 
  }
 
  // -----------------------------------------------------------------------------------------------
 
  template <class Matrix, class Vectors, class ReactionDerivatives, class Solver>
  typename Matrix::field_type
  sdcIterationStep(SDCTimeGrid const& grid, SDCTimeGrid::RealMatrix const& Shat, Solver const& solve,
                   Matrix const& M, Matrix const& A,
                   Vectors const& r, ReactionDerivatives const& f_u, Vectors const& u, Vectors& du)
  {
    using namespace SdcDetail;
    auto Mdu = computeMdu(M,u);
    return sdcIterationStep2(grid,Shat,solve,M,A,r,f_u,Mdu,du);
  } 
 
  template <class Matrix, class Vectors, class Solver>
  typename Matrix::field_type
  sdcIMEXStep(SDCTimeGrid const& grid, SDCTimeGrid::RealMatrix const& Shat, Solver const& solve,
              Matrix const& M, Matrix const& A,
              Vectors const& r, Vectors const& u, Vectors& du)
  {
    using namespace SdcDetail;
    assert(isSinglyDiagonal(Shat));
 
    auto Mdu = computeMdu(M,u);
    auto solver = [&](auto const& /* J */, auto& x, auto const& b)
    {
      solve(x,b);
    };
    return sdcIterationStep2(grid,Shat,solver,M,A,r,nullptr,Mdu,du);
  }
  
  template <class Matrix, class Vectors, class ReactionDerivatives, class Solver>
  typename Matrix::field_type
  sdcIterationStep2(SDCTimeGrid const& grid, SDCTimeGrid::RealMatrix const& Shat, Solver const& solve,
                    Matrix const& M, Matrix const& A,
                    Vectors const& rUi, ReactionDerivatives const& f_u, Vectors const& Mdu, Vectors& du)
  {
    using namespace SdcDetail;
    auto& timer = Timings::instance();
    ScopedTimingSection("SDC sweep",timer);
 
    auto const& pts = grid.points();
    int const n = pts.size()-1;      // number of subintervals
 
    constexpr bool reactionExplicit = std::is_same_v<ReactionDerivatives,std::nullptr_t>;
    bool const matrixIsFixed = reactionExplicit && isSinglyDiagonal(Shat);
 
    assert(Mdu.size()>=n); 
    assert(du.size()>=n+1);  // including start point
    
    size_t const dofs = Mdu[0].size();
 
    // compute exact integration matrix
    auto const& S = grid.integrationMatrix();
 
 
    // initialize correction at starting point to zero
    du[0] = 0.0; 
 
    typedef typename Vectors::value_type Vector;
    Vector rhs(dofs), tmp(dofs);  // declare here to prevent frequent reallocation
    
    // perform n "Euler" steps
    typename Matrix::field_type norm = 0;
    Matrix J = M;
 
    for (int i=1; i<=n; i++)
    {
      timer.start("building matrix");
      // matrix J = M - Shat_i-1,i*(A+f_u)
      if constexpr (reactionExplicit)
      {
        if (i==1 || !matrixIsFixed)                         // only rebuild if it might change
          buildEulerStepMatrix(M,Shat[i-1][i],A,nullptr,J);
      }
      else
        buildEulerStepMatrix(M,Shat[i-1][i],A,f_u[i],J);
      timer.stop("building matrix");
 
      // build right-hand side for linear system
      timer.start("building rhs");
      // M * ( u_i^{k} - u_{i+1}^k + du_i)
      rhs = Mdu[i-1];
      M.umv(du[i-1],rhs);
 
      // add sum_j S_ij r_j to right hand side
      for (int j=0; j<=n; ++j)
        rhs.axpy(S[i-1][j],rUi[j]);
         
      // add sum_j Shat_ij r'(u_j) du_j with r' = A + f_u, but ignore
      // f_u if that is not provided
      tmp = 0;
      for (int j=1; j<i; ++j) // we can start at 1 since du[0] is zero
      {
        if constexpr ( ! reactionExplicit)
          f_u[j].usmv(Shat[i-1][j],du[j],rhs);
 
        tmp.axpy(Shat[i-1][j],du[j]);
      }
      A.umv(tmp,rhs);
      timer.stop("building rhs");
 
      // solve linear system
      timer.start("solve");
      du[i] = du[i-1]; // previous increment is probably a good starting value
      solve(J,du[i],rhs);
      timer.stop("solve");
      
      // evaluate norm of correction
      norm += (pts[i]-pts[i-1]) * (du[i]*rhs);
    }   // end i - loop
    
    return std::sqrt(norm/(pts[n]-pts[0]));
  } 
  
// THIS HAS BEEN A TESTBED FOR PLAYING WITH SDC VARIANTS. WHILE INTERESTING, IN ITS CURRENT FORM
// IT DOES NOT BELONG INTO THE KASKADE CORE LIBRARY. [MW 2024]
//
//   template <class Matrix, class Vectors, class Reaction, class Solver>
//   typename Matrix::field_type
//   sdcIterationStep3(SDCTimeGrid const& grid, SDCTimeGrid::RealMatrix const& Shat, Solver const& solve,
//                     Matrix const& M, Matrix const& Stiff,
//                     Vectors const& rUi, Reaction const& fAt, Vectors const& Mdu, Vectors& du,
//                     int nReactionSweeps)
//   {
//     auto const& pts = grid.points();
//     int const n = pts.size()-1;      // number of subintervals
//
//     assert(Mdu.size()>=n);
//     assert(du.size()>=n+1);  // including start point
//
//     size_t const dofs = Mdu[0].size();
//
//     // compute exact integration matrix
//     auto const& S = grid.integrationMatrix();
// // auto const& S = Shat;
//
//     // initialize correction at starting point to zero
//     du[0] = 0.0;
//
//     typedef typename Vectors::value_type Vector;
//     Vector rhs(dofs), tmp(dofs);  // declare here to prevent frequent reallocation
//
//
//     // perform n basic steps
//     typename Matrix::field_type norm = 0;
//     Matrix J = M;
//     for (int i=1; i<=n; i++)
//     {
//       // Each basic step consists of a splitting method, separating a linearly implicit Euler step
//       // from the pointwise nonlinearity of the right hand side. First we perform the linearly implicit
//       // Euler step.
//
//       //  right-hand side for linear system
//
//       // M * ( u_{i-1}^{k} - u_{i}^k + du_{i-1})
//       rhs = Mdu[i-1];
//       M.umv(du[i-1],rhs);
//
//       // add sum_j S_ij r_j to right hand side
//       for (int j=0; j<=n; ++j)
//         rhs.axpy(S[i-1][j],rUi[j]);
//
//       // add sum_j Shat_ij r'(u_j) du_j with r' = A + f_u, i.e. A (sum_j Shat_ij du_j) + sum_j Shat_ij f_uj du_j
//       // addition of f_u is postponed to matrix loop
//       tmp = 0;
//       for (int j=0; j<i; ++j) // TODO: start at 1 instead of 0? du[0] is zero anyway...
//         tmp.axpy(Shat[i-1][j],du[j]);
//       Stiff.umv(tmp,rhs);
//
//       // matrix J = M - Shat_i-1,i*(A+f_u)
//       auto const f = fAt( pts[i] );
//       for (size_t row=0; row<J.N(); ++row)
//       {
//         auto colJ = J[row].begin();
//         auto end = J[row].end();
//         auto colM = M[row].begin();
//         auto colA = Stiff[row].begin();
//
//         while (colJ != end)
//         {
//           auto fu = *colM * (f(row,0.0,1) + f(colJ.index(),0.0,1)) / 2;
// //           *colJ = *colM - Shat[i-1][i] * (*colA + fu); // exact Newton, but for large step sizes the Jacobi matrix becomes singular...
// //           *colJ = *colM - Shat[i-1][i] * (*colA + std::min(0.0,fu));  // this modification guarantees an invertible Jacobian
//           *colJ = *colM - Shat[i-1][i] * (*colA);         // explicit reaction
//           if (colJ.index()==row)
//             *colJ -= Shat[i-1][i] * fu; // only diagonal of reaction mass matrix
//
// //           // add f_u contribution to right hand side (no-op for Euler SDC as Shat[i-1][j]=0 for j<i)
// //           for (int j=0; j<i; ++j) // TODO: start at 1 instead of 0? du[0] is zero anyway...
// //           {
// //             auto f = fAt(pts[j]);
// //             rhs[row] += *colM*Shat[i-1][j]*(f(row,0.0,1)+f(colJ.index(),0.0,1))/2 * du[j][colJ.index()];
// //           }
//
//           ++colJ; ++colM; ++colA;
//         }
//       }
//
//       // solve linear system
//       du[i] = du[i-1]; // previous increment is probably a good starting value
//       solve(J,du[i],rhs);
//
// std::cerr << "du=" << du[i] << "\n";
// // std::cerr << "|du| = " << du[i].two_norm2() << "  du[i]=" << du[i][0] << "  du[i-1]=" << du[i-1][0] << "\n";
//       // Second part is the remaining nonlinearity, which we address by a subsequent run of a couple of Euler steps
//       double const tau = pts[i]-pts[i-1];
//       double snorm = 0;
//       for (size_t row=0; row<J.N(); ++row)
//       {
//         auto fend = fAt(pts[i]);
//         double dfdu = std::min(0.0,fend(row,0.0,1)) * du[i][row];
//
//         double s = 0;
//         int l = nReactionSweeps;
//         for (int j=0; j<l; ++j)
//         {
//           double theta = 0.0; // where in the subinterval to evaluate (theta=0 explicit Euler, theta=0.5 lin. impl. midpoint, theta=1 lin. impl Euler)
//           auto f = fAt(pts[i-1]+(j+theta)*tau/l);
//           double dut = ((l-j-theta)*du[i-1][row]+(j+theta)*du[i][row])/l;
//
// // if (row==0)
// // std::cerr << "j=" << j << " du+s=" << dut+s << "   f(u+du+s)=" <<  f(row,dut+s,0) << "  f(u)=" <<  f(row,0.0,0) << "  dfdu=" <<  dfdu  <<  "  sum rhs=" << (f(row,dut+s,0) - f(row,0.0,0) - dfdu);
// //           s += tau/l * (f(row,dut+s,0) - f(row,0.0,0) - dfdu) / (1-1.0*tau/l*std::min(0.0,f(row,dut+s,1)));
// if (row==0)
// std::cerr << "j=" << j << " du+s=" << dut+s << "   f(u+du+s)=" <<  f(row,dut+s,0) << "  f(u)=" <<  f(row,0.0,0) << "  dfdu=" <<  f(row,dut+s,1)  <<  "  sum rhs=" << (f(row,dut+s,0) - f(row,0.0,0));
//           s += tau/l * (f(row,dut+s,0) - f(row,0.0,0) - f(row,0.0,1)*(dut+s)) / (1-theta*tau/l*std::min(0.0,f(row,dut+s,1)));
// if (row==0) std::cerr << "  -> s=" << s << "\n";
//         }
//         du[i][row] += s;
//         snorm += s*s;
//       }
// // std::cerr << "|s| = " << snorm << "\n";
//
//       // evaluate norm of correction
//       norm += (pts[i]-pts[i-1]) * (du[i]*rhs);
//     }   // end i - loop
//
//     return std::sqrt(norm/(pts[n]-pts[0]));
//   }
  
  template <class Matrix, class Vectors, class ReactionDerivatives>
  double waveformRelaxationStep2(SDCTimeGrid const& grid, 
                                 Matrix const& M, Matrix const& A,
                                 Vectors const& rUi, ReactionDerivatives const& f_u, Vectors& du)
  {
    auto const& pts = grid.points();
    int const n = pts.size()-1;      // number of subintervals -- left interval boundary always included in points
    
 
    assert(rUi.size()>=n); 
    assert(du.size()>=n+1);  // including start point - usually 0
    
    size_t const dofs = rUi[0].size();
    
    // We need to solve the scalar ODE
    // Mii y_t = Aii y + Rii y + ri, y(0) = 0
    // We approximate y_t by the spectral differentiation matrix D as Dy and obtain, as system in the 
    // collocation points (excluding the initial value 0)
    // (MiiD - diag(Aii+Rii)) y = r. 
    // The system matrix is called S here.
    
    DynamicMatrix<double> S(n,n);  // n x n system matrix TODO: use Kaskade::DynamicMatrix?
    Dune::DynamicVector<double> r(n), y(n);
    auto const& D = grid.differentiationMatrix();
    
    std::cerr << "using D = \n" << D << "\n";
    
    // a function to extract diagonal entries 
    auto diagonal = [](Matrix const& M, size_t i) -> double
    {
      auto const& arow = M[i];
      auto first = arow.begin();
      auto last = arow.end();
      while (first!=last && first.index()<i)
        ++first;
      if (first==last)
      {
        std::cerr << "first==last\n"; 
        std::cerr.flush();
      }
      return *first;
    };
    
    double retval = 0;
    
    for (size_t i=0; i<dofs; ++i)
    {
      std::cerr.flush();
      auto mii = diagonal(M,i);
      auto aii = diagonal(A,i);
      
      for (int j=0; j<n; ++j)
      {
        for (int k=0; k<n; ++k)
          S[j][k] = mii*D[j+1][k+1]; // remember D starts with left interval point
        S[j][j] -= aii+diagonal(f_u[j+1],i);
        r[j] = rUi[j+1][i];
      }
      S.solve(y,r);
      for (int j=0; j<n; ++j)
        du[j+1][i] = 0.5*y[j]; // Jacobi damping (otherwise we get spatially high-frequent errors)
      retval += 0.25 * y.two_norm2();
    }
    
    return std::sqrt(retval);
  }
  
  template <class Vector>
  class Sdc 
  {
  public:
    Sdc(SDCTimeGrid const& ts_, Vector const& u0)
    : ts(ts_), us(ts.points().size(),u0)
    {
    }
    
    void setApproximateQuadratureMatrix(SDCTimeGrid::RealMatrix const& Shat_)
    {
      Shat = Shat_;
    }
    
    void setInitialValue(Vector const& u0, double decayRate=0.0)
    {
      Vector du = u0;
      du -= us[0];
      for (int k=0; k<us.size(); ++k)
        us[k].axpy(std::exp(-(ts.points()[k]-ts.points()[0])*decayRate),du);
    }
    
  protected:
    
  private:
    SDCTimeGrid             ts;
    std::vector<Vector>     us;
    SDCTimeGrid::RealMatrix Shat;
  };
 
}
 
 
#endif