semieuler.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) 2002-2021 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef SEMIEULERINNER_HH
#define SEMIEULERINNER_HH
 
#include <cmath>
 
#include "fem/fixdune.hh"
#include "fem/variables.hh"
#include "fem/functional_aux.hh"
namespace Kaskade
{
  namespace SemiEulerDetail
  {
    template <class ParabolicEquation, class SemiEuler, bool hasInnerBoundaryCache> 
    struct SemiImplicitEulerInnerBoundaryPolicyImpl
    {
    };
 
    template <class ParabolicEquation, class SemiEuler>
    struct SemiImplicitEulerInnerBoundaryPolicyImpl<ParabolicEquation,SemiEuler,true>
    {
    private:
      using SemiEulerStep = SemiEuler;
      
    public:
      class InnerBoundaryCache
      {
        using AnsatzVars = typename SemiEulerStep::AnsatzVars;
        using TestVars = typename SemiEulerStep::TestVars;
        using AnsatzVariableSet = typename AnsatzVars::VariableSet;
        using Scalar = typename ParabolicEquation::Scalar;
        
      public:
        InnerBoundaryCache(SemiEulerStep const& f_,
                           AnsatzVariableSet const& vars_,
                           AnsatzVariableSet const& varsJ_,
                           AnsatzVariableSet const& du_,
                           int flags)
      : f(f_)
      , peibc(f.parabolicEquation(),vars_,flags)
      , peibcJ(f.parabolicEquation(),varsJ_,flags)
      , du(du_)  
      {}
 
        
        template <class FaceIterator>
        void moveTo(FaceIterator const& face)
        {
        peibc.moveTo(face);
        peibcJ.moveTo(face);
        }
        
        template <class Position, class Evaluators>
        void evaluateAt(Position const& x, Evaluators const& evaluators, Evaluators const& neighbourEvaluators)
        {
          peibc.evaluateAt(x, evaluators, neighbourEvaluators);
          peibcJ.evaluateAt(x, evaluators,neighbourEvaluators);
        }
        
        
        template <int row, int dim>
        Dune::FieldVector<Scalar,TestVars::template Components<row>::m> d1(VariationalArg<Scalar,dim> const& argT) const
        {
          return f.getTau() * peibc.template d1<row>(argT);
        }
        
        template <int row, int col, int dim>
        Dune::FieldMatrix<Scalar,TestVars::template Components<row>::m,AnsatzVars::template Components<row>::m>
        d2(VariationalArg<Scalar,dim> const& argT, VariationalArg<Scalar,dim> const& argA, bool centerCell) const
        {
          Dune::FieldMatrix<Scalar,TestVars::template components<row>,
                                   AnsatzVars::template components<col>> result(0);
          
          
          if (ParabolicEquation::template B2<row,col>::present)
            result = peibcJ.template b2<row,col>(argT, argA, centerCell);
          
          if (ParabolicEquation::template D2<row,col>::present)
            result -= f.getTau()*peibcJ.template d2<row,col>(argT, argA, centerCell);
          
          return result;
        }
        
 
      private:
        SemiEulerStep const&                           f;
        typename ParabolicEquation::InnerBoundaryCache peibc;
        typename ParabolicEquation::InnerBoundaryCache peibcJ;
        AnsatzVariableSet const&                       du;
      };
 
      template <class Face>
      bool considerFace(Face const& face) const
      {
        return static_cast<SemiEuler const*>(this)->parabolicEquation().considerFace(face);
      }
 
    };
  }
  
  
  template <class PE>
  class SemiImplicitEulerStep
  : public SemiEulerDetail::SemiImplicitEulerInnerBoundaryPolicyImpl<PE,SemiImplicitEulerStep<PE>,
                                                                     hasInnerBoundaryCache<PE>>   
  {
    using Self = SemiImplicitEulerStep<PE>;
    using Base = SemiEulerDetail::SemiImplicitEulerInnerBoundaryPolicyImpl<PE,Self,hasInnerBoundaryCache<PE>>;
 
  public:
    typedef PE ParabolicEquation;
 
    typedef typename ParabolicEquation::Scalar  Scalar;
    typedef typename ParabolicEquation::AnsatzVars AnsatzVars;
    typedef typename ParabolicEquation::TestVars TestVars;
    typedef typename ParabolicEquation::OriginVars OriginVars;
    typedef typename AnsatzVars::Grid Grid;
    using GridView = typename AnsatzVars::GridView;
 
    static ProblemType const type = WeakFormulation;
 
    typedef PE  EvolutionEquation;
 
    typedef typename EvolutionEquation::AnsatzVars::template CoefficientVectorRepresentation<0,EvolutionEquation::TestVars::noOfVariables>::type CoefficientVectors;
 
    
    
    class DomainCache
    {
    public:
      DomainCache(SemiImplicitEulerStep<ParabolicEquation> const& f_, // TODO: implement tauRatio
                  typename AnsatzVars::VariableSet const& vars_,
                  typename AnsatzVars::VariableSet const& varsJ_,
                  typename AnsatzVars::VariableSet const& duVars_,
                  int flags):
        f(f_), pedc(*f_.eq,vars_,flags), pedcJ(*f_.eq,varsJ_,flags), duVars(duVars_) 
      {} 
 
      void moveTo(typename Grid::template Codim<0>::Entity const& entity) {
        pedc.moveTo(entity);
        pedcJ.moveTo(entity);
      }
 
      template <class Position, class Evaluators>
      void evaluateAt(Position const& x, Evaluators const& evaluators) 
      {
        pedc.evaluateAt(x, evaluators);
        pedcJ.evaluateAt(x, evaluators);
        du = evaluateVariables(duVars,evaluators,valueMethod);
      }
 
      template<int row, int dim>
      Dune::FieldVector<Scalar, TestVars::template components<row>>
      d1 (VariationalArg<Scalar,dim> const& arg) const 
      {
        Dune::FieldVector<Scalar, TestVars::template components<row>> result( f.getTau() * pedc.template d1<row>(arg) );
 
        boost::fusion::for_each(typename AnsatzVars::Variables(),MultiplyDifference<row,dim>(pedc,pedcJ,du,arg,result));
        return result;
      }
 
      template<int row, int col, int dim>
      Dune::FieldMatrix<Scalar, TestVars::template components<row>, AnsatzVars::template components<col>>
      d2 (VariationalArg<Scalar,dim> const &arg1, VariationalArg<Scalar,dim> const &arg2) const
      {
        Dune::FieldMatrix<Scalar,TestVars::template components<row>,AnsatzVars::template components<col>> result(0);
 
 
        if (ParabolicEquation::template B2<row,col>::present)
          result = pedcJ.template b2<row,col>(arg1, arg2);
 
        if (ParabolicEquation::template D2<row,col>::present)
          result -= f.getTau()*pedcJ.template d2<row,col>(arg1, arg2);
        
        return result;
      }
 
    private:
      SemiImplicitEulerStep<ParabolicEquation> const& f;
      typename ParabolicEquation::DomainCache pedc;  // for evaluating rhs at u_k
      typename ParabolicEquation::DomainCache pedcJ; // for evaluating Jacobian at u_0
      typename AnsatzVars::VariableSet const& duVars;
      EvaluateVariables<typename AnsatzVars::VariableSet,ValueMethod> du;
 
 
      template <int row, int dim>
      class MultiplyDifference
      {
      public:
        MultiplyDifference(typename ParabolicEquation::DomainCache pedc_,
                           typename ParabolicEquation::DomainCache pedcJ_,
                           EvaluateVariables<typename AnsatzVars::VariableSet,ValueMethod> const& du_,
                           VariationalArg<Scalar,dim> const& arg_,
                           Dune::FieldVector<Scalar, TestVars::template Components<row>::m>& result_)
        : pedc(pedc_), pedcJ(pedcJ_), du(du_), arg(arg_), result(result_)
        {}
 
        template <class VarDescription>
        void operator()(VarDescription const& vd) const 
        {
          // compute <(B(u0)-B(uk))*du,arg>
          int const col = VarDescription::id;
 
          if (ParabolicEquation::template B2<row,col>::present && !ParabolicEquation::template B2<row,col>::constant) 
          {
            // @TODO: should work with vector-valued shape functions
            VariationalArg<Scalar,dim> duArg;
            duArg.value[0] = 1;
            result += (pedcJ.template b2<row,col>(arg,duArg)-pedc.template b2<row,col>(arg,duArg)) * boost::fusion::at_c<col>(du);
          }
        }
 
      private:
        typename ParabolicEquation::DomainCache const& pedc;
        typename ParabolicEquation::DomainCache const& pedcJ;
        EvaluateVariables<typename AnsatzVars::VariableSet,ValueMethod> const& du;
        VariationalArg<Scalar,dim> const& arg;
        Dune::FieldVector<Scalar, TestVars::template Components<row>::m>& result;
      };
 
    };
 
    // TODO: currently  a dependence of the boundary part of B on u is not implemented
    class BoundaryCache
    {
    public:
      typedef typename AnsatzVars::GridView::IntersectionIterator FaceIterator;
 
      BoundaryCache(SemiImplicitEulerStep<ParabolicEquation> const& f_,
                    typename AnsatzVars::VariableSet const& vars_,
                    typename AnsatzVars::VariableSet const& varsJ_,
                    typename AnsatzVars::VariableSet const& du_,
                    int flags)
      : f(f_), pedc(*f_.eq,vars_,flags), pedcJ(*f_.eq,varsJ_,flags), du(du_)  {}
 
      void moveTo(FaceIterator const& entity) 
      {
        pedc.moveTo(entity);
        pedcJ.moveTo(entity);
      }
 
      template <class Evaluators>
      void evaluateAt(Dune::FieldVector<typename Grid::ctype,Grid::dimension-1> const& x,
                      Evaluators const& evaluators) 
      {
        pedc.evaluateAt(x, evaluators);
        pedcJ.evaluateAt(x, evaluators);
      }
 
      template<int row, int dim>
      Dune::FieldVector<Scalar, TestVars::template components<row>> d1 (VariationalArg<Scalar,dim> const& arg) const
      {
        return f.getTau() * pedc.template d1<row>(arg);
      }
 
      template<int row, int col, int dim>
      Dune::FieldMatrix<Scalar, TestVars::template components<row>, AnsatzVars::template components<col>>
      d2 (VariationalArg<Scalar,dim> const &arg1, VariationalArg<Scalar,dim> const &arg2) const 
      {
        // to be corrected - include b2
        if (ParabolicEquation::template D2<row,col>::present)
          return - f.getTau()*pedcJ.template d2<row,col>(arg1, arg2);
 
        return 0.0; // never get here
      }
 
    private:
      SemiImplicitEulerStep<ParabolicEquation> const& f;
      typename ParabolicEquation::BoundaryCache pedc;
      typename ParabolicEquation::BoundaryCache pedcJ;
      typename AnsatzVars::VariableSet const& du;
    };
 
 
 
  public:
    template <int row, int col>
    struct D2: public ParabolicEquation::template D2<row,col>
    {
      static bool const present = (ParabolicEquation::template D2<row,col>::present ||
                                   ParabolicEquation::template B2<row,col>::present);
      static bool const symmetric = (!ParabolicEquation::template D2<row,col>::present ||
                                     ParabolicEquation::template D2<row,col>::symmetric) &&
                                    (!ParabolicEquation::template B2<row,col>::present ||
                                     ParabolicEquation::template B2<row,col>::symmetric);
      //     static bool const constant = (!ParabolicEquation::template D2<row,col>::present ||
      //                                   ParabolicEquation::template D2<row,col>::constant) &&
      //                                  (!ParabolicEquation::template B2<row,col>::present ||
      //                                   ParabolicEquation::template B2<row,col>::constant);
      static bool const lumped = ((ParabolicEquation::template D2<row,col>::lumped ||
          !ParabolicEquation::template D2<row,col>::present) &&
          (ParabolicEquation::template B2<row,col>::lumped||
              !ParabolicEquation::template B2<row,col>::present));
 
    };
 
    template <int row>
    struct D1: public ParabolicEquation::template D1<row>
    {
    };
 
    SemiImplicitEulerStep(ParabolicEquation *eq_, Scalar dt) 
    : eq(eq_)
    , tau(dt)
    {}
 
    template <class Cell>
    int integrationOrder(Cell const& cell, int shapeFunctionOrder, bool boundary) const
    {
      return eq->integrationOrder(cell,shapeFunctionOrder,boundary);
    }
 
    Self& setTau(Scalar dt)
    { 
      tau = dt; 
      assert(tau >= 0);
      return *this;
    }
 
    Scalar getTau() const { return tau; }
 
 
    ParabolicEquation&       parabolicEquation()       { return *eq; }
    ParabolicEquation const& parabolicEquation() const { return *eq; }
    private:
      ParabolicEquation *eq;
      Scalar tau;                        // the time step
      bool onlyJacobian = false;         // 
      bool onlyMassMatrix = false;       // 
  };
} // namespace Kaskade
#endif