parabolicfemscheme.hh

#ifndef __DUNE_ACFEM_PARABOLIC_FEMSCHEME_HH__
#define __DUNE_ACFEM_PARABOLIC_FEMSCHEME_HH__
 
// iostream includes
#include <iostream>
 
// include grid part
#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
 
// include discrete function space
#include <dune/fem/space/lagrange.hh>
 
// adaptation ...
#include <dune/fem/function/adaptivefunction.hh>
#include <dune/fem/space/common/adaptationmanager.hh>
 
// include discrete function
#include <dune/fem/function/blockvectorfunction.hh>
 
// Newton's iteration
#include <dune/fem/solver/newtoninverseoperator.hh>
 
// natural interpolation
#include <dune/fem/space/common/interpolate.hh>
 
// include norms
#include <dune/fem/misc/l2norm.hh>
#include <dune/fem/misc/h1norm.hh>
 
// include parameter handling
#include <dune/fem/io/parameter.hh>
 
// include vtk output and checkpointing
#include <dune/fem/io/file/dataoutput.hh>
#include <dune/fem/io/file/datawriter.hh>
 
// local includes
#include "../algorithms/femschemeinterface.hh"
#include "../algorithms/operatorselector.hh"
#include "../operators/ellipticoperator.hh"
#include "../operators/constraints/dirichletconstraints.hh"
#include "../operators/functionals/functionals.hh"
#include "../functions/functions.hh"
#include "../estimators/parabolicestimator.hh"
#include "../estimators/trueerrorestimator.hh"
#include "../algorithms/marking.hh"
#include "../common/timeview.hh"
#include "../common/dataoutput.hh"
 
// Select an appropriate solver, depending on ModelType and solver-family (ISTL, PESC ...)
#include "../common/solverselector.hh"
 
namespace Dune {
 
  namespace ACFem {
 
    template<class DiscreteFunction,
             class TimeProvider,
             class ImplicitModel,
             class ExplicitModel,
             class InitialValue,
             class RHSFunctional,
             template<class> class QuadratureTraits = DefaultQuadratureTraits>
    class ParabolicFemScheme
      : public TransientAdaptiveFemScheme
    {
     public:
      typedef DiscreteFunction DiscreteFunctionType;
 
      typedef typename DiscreteFunctionType::GridType GridType;
 
      typedef typename DiscreteFunctionType::GridPartType GridPartType;
 
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
 
      typedef TimeProvider TimeProviderType;
 
      using DiscreteModelType = DiscreteModel<ImplicitModel, DiscreteFunctionSpaceType>;
      using ImplicitTraits = ModelTraits<DiscreteModelType>;
 
      using ImplicitModelType = ImplicitModel;
      using ExplicitModelType = ExplicitModel;
 
      typedef InitialValue InitialValueType;
 
      typedef RHSFunctional RHSFunctionalType; 
 
      typedef typename DiscreteModelType::FunctionSpaceType FunctionSpaceType;
 
      // or: (must coincide)
      //typedef typenema DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
 
      // choose type of discrete function, Matrix implementation and solver implementation
      typedef SolverSelector<DiscreteFunctionType, DiscreteModelType> SolverSelectorType;
      typedef typename SolverSelectorType::LinearOperatorType LinearOperatorType;
      typedef typename SolverSelectorType::LinearInverseOperatorType LinearInverseOperatorType;
      typedef typename LinearInverseOperatorType::SolverParameterType LinearSolverParameterType;
 
      typedef
      Fem::NewtonInverseOperator<LinearOperatorType, LinearInverseOperatorType>
      NonLinearInverseOperatorType;
 
      typedef Dune::Fem::RestrictProlongDefault<DiscreteFunctionType> RestrictionProlongationType;
 
      typedef Dune::Fem::AdaptationManager<GridType, RestrictionProlongationType> AdaptationManagerType;
 
      typedef EmptyBlockConstraints<DiscreteFunctionSpaceType>  EmptyConstraintsType;
 
      using ConstraintsOperatorType = DirichletConstraints<DiscreteFunctionSpaceType, DiscreteModelType>;
 
      typedef DifferentiableEllipticOperator<
        LinearOperatorType, DiscreteModelType, const ConstraintsOperatorType&, QuadratureTraits>
      ImplicitOperatorType;
 
      typedef EllipticOperator<ExplicitModelType,
                               DiscreteFunctionType, DiscreteFunctionType,
                               EmptyConstraintsType, QuadratureTraits> ExplicitOperatorType;
 
      typedef
      ParabolicEulerEstimator<DiscreteFunctionType, TimeProviderType,
                              DiscreteModelType,
                              ExplicitModelType>
      EstimatorType;
 
      typedef MarkingStrategy<GridPartType> MarkingStrategyType;
 
      typedef
      TrueErrorEstimator<InitialValueType, Dune::Fem::L2Norm<GridPartType> >
      InitialEstimatorType;
 
      typedef MarkingStrategy<GridPartType> InitialMarkingStrategyType;
 
      typedef InitialValueType ExactSolutionFunctionType; // slight abuse
 
      typedef std::tuple<DiscreteFunctionType *> IOTupleType;
 
      typedef DataOutput<GridType, IOTupleType> DataOutputType;
 
      typedef Dune::Fem::CheckPointer<GridType>   CheckPointerType;
 
      ParabolicFemScheme(DiscreteFunctionType& solution,
                         const TimeProviderType& timeProvider,
                         const ImplicitModelType& implicitModel,
                         const ExplicitModelType& explicitModel,
                         const InitialValueType& initialValue,
                         const RHSFunctionalType& rhsFunctional,
                         const std::string name = "heat")
        : grid_(solution.space().gridPart().grid()),
          timeProvider_(timeProvider),
        implicitModel_(asExpression(discreteModel(implicitModel, solution.space()))),
        explicitModel_(explicitModel),
        name_(name),
        gridPart_(solution.space().gridPart()),
        discreteSpace_(solution.space()),
        initialValue_(initialValue),
        rhsFunctional_(rhsFunctional),
        solution_(solution),
        oldSolution_("old solution", discreteSpace_),
        residual_("residual", discreteSpace_),
        update_("update", discreteSpace_),
        estimator_(timeProvider_, implicitModel_, explicitModel_, oldSolution_),
        markingStrategy_(gridPart_, estimator_, name_ + ".adaptation.space"),
        initialEstimator_(initialValue_),
        initialMarkingStrategy_(gridPart_, initialEstimator_, name_ + ".adaptation" + ".initial"),
        // restriction/prolongation operator
        restrictProlong_(0),
        // adaptation manager
        adaptationManager_(0),
        // create Dirichlet contraints
        constraints_(discreteSpace_, implicitModel_),
        // the elliptic operator (implicit)
        implicitOperator_(implicitModel_, constraints_),
        explicitOperator_(explicitModel_),
        // create linear operator (domainSpace,rangeSpace)
        linearOperator_("assembled elliptic operator", discreteSpace_, discreteSpace_),
        solverParams_(ImplicitTraits::isAffineLinear ? name_ + ".solver." : name_ + ".linear."),
        // other stuff
        sequence_(-1),
        ioTuple_(&solution_),
        dataOutput_(0),
        checkPointer_(grid_, timeProvider_)
      {
        // Make the solution of the current timestep persistent for check-pointing.
        Dune::Fem::persistenceManager << solution_;
 
        oldSolution_.clear();
        solution_.clear();
        chooseDefaultSolverMethod<SolverSelectorType>(solverParams_);
      }
 
      ParabolicFemScheme(const ParabolicFemScheme& other)
        : grid_(other.grid_),
          timeProvider_(other.timeProvider_),
          implicitModel_(other.implicitModel_),
          explicitModel_(other.explicitModel_),
          name_(other.name_),
          gridPart_(other.gridPart_),
          discreteSpace_(other.discreteSpace_),
          initialValue_(other.initialValue_),
          rhsFunctional_(other.rhsFunctional_),
          solution_(other.solution_),
          oldSolution_("old solution", discreteSpace_),
          residual_("residual", discreteSpace_),
          update_("update", discreteSpace_),
          estimator_(timeProvider_, implicitModel_, explicitModel_, oldSolution_),
          markingStrategy_(gridPart_, estimator_, name_ + ".adaptation.space"),
          initialEstimator_(other.initialEstimator_),
          initialMarkingStrategy_(other.initialMarkingStrategy_),
          // restriction/prolongation operator
          restrictProlong_(0),
          // adaptation manager
          adaptationManager_(0),
          // create Dirichlet contraints
          constraints_(discreteSpace_, implicitModel_),
          // the elliptic operator (implicit)
          implicitOperator_(implicitModel_, constraints_),
          explicitOperator_(explicitModel_),
          // create linear operator (domainSpace,rangeSpace)
          linearOperator_("assembled elliptic operator", discreteSpace_, discreteSpace_),
          solverParams_(other.solverParams_),
          // other stuff
          sequence_(-1),
          ioTuple_(&solution_),
          dataOutput_(0),
          checkPointer_(grid_, timeProvider_)
      {
        // Make the solution of the current timestep persistent for check-pointing.
        Dune::Fem::persistenceManager << solution_;
 
        oldSolution_.clear();
        solution_.clear();
        chooseDefaultSolverMethod<SolverSelectorType>(solverParams_);
      }
 
      virtual ~ParabolicFemScheme() {
        if (dataOutput_) {
          delete dataOutput_;
        }
        if (adaptationManager_) {
          delete adaptationManager_;
        }
        if (restrictProlong_) {
          delete restrictProlong_;
        }
      }
 
      virtual void initialize()
      {
        // initialValue_ already has to be a grid-function
 
        // apply natural interpolation
        interpolate(initialValue_, oldSolution_);
 
        // copy it over
        solution_.assign(oldSolution_);
 
        // std::cout << solution_ << std::endl;
      }
 
      virtual std::string name() const
      {
        return name_;
      }
 
      virtual void restart(const std::string& name = "") {
        std::string name_ = name;
        if (name_ == "") {
          name_ = Dune::Fem::CheckPointerParameters().checkPointPrefix();
        }
        checkPointer_.restoreData(grid_, name_);
        next();
      }
 
      virtual void next()
      {
        oldSolution_.assign(solution_);
      }
 
      virtual void solve(bool forceMatrixAssembling = true)
      {
        // Need to copy the data again in case the mesh was adapted in
        // between because the AdaptationManager only works on a single function ATM.
        solution_.assign(oldSolution_);
 
        if (ImplicitTraits::template Exists<PDEModel::dirichlet>::value) {
          // set Dirichlet constraints, new time
          constraints_.updateValues(); // DO NOT CHANGE THIS LINE. DON'T!!!!
          constraints_.constrain(solution_);
        }
 
        residual_.clear();
        residual_ += rhsFunctional_;
 
        // Compute the contribution from the time derivative
        explicitOperator_(oldSolution_, update_);
        residual_ -= update_;
 
        // do no spoil the new Dirichlet-values
        constraints_.zeroConstrain(residual_);
 
        if (ImplicitTraits::isAffineLinear) {
          linearSolve(residual_, forceMatrixAssembling);    // sove Lu = res
        } else {
          nonLinearSolve(residual_); // solve Lu = res
        }
      }
 
     protected:
 
      virtual void nonLinearSolve(DiscreteFunctionType& rhs)
      {
        assert(!ImplicitTraits::isAffineLinear);
 
        Dune::Fem::NewtonParameter newtonParam(solverParams_, name_ + ".newton.");
        NonLinearInverseOperatorType newton(newtonParam);
        newton.bind(implicitOperator_);
 
        newton(rhs, solution_);
 
        converged_ = newton.converged();
        assert(newton.converged());
      }
 
      virtual void linearSolve(DiscreteFunctionType& rhs, bool forceMatrixAssembling)
      {
        assert(ImplicitTraits::isAffineLinear);
 
        if (sequence_ != discreteSpace_.sequence() || forceMatrixAssembling) {
          // assemble linear operator (i.e. setup matrix) The jacobian
          // incorporates the constraints defined by the constraint class.
          implicitOperator_.jacobian(solution_, linearOperator_);
 
          // update sequence number
          sequence_ = discreteSpace_.sequence();
        }
 
        // inverse operator using linear operator
        //
        // Note: ATM, we use the difference quotient for the time
        // derivative, i.e. 1/tau as factor in from of the
        // mass-matrix. Consequently, the solver-tolerance should be
        // scaled by the time-step size, at least when we use any
        // preconditioner.
        //double scaledEps_ = solverEps_ * timeProvider_.inverseDeltaT();
        LinearInverseOperatorType solver(solverParams_);
        solver.bind(linearOperator_);
 
        // Apply the affine linear operator to the start value. This
        // computes the initial residual. In the linear case, Lagrange space,
        // Dirichlet data g, this computes
        //
        // update_ = A u - g
        //
        // where it is allowed that A is affine-linear, without having
        // to apply a non-linear solver for this trivial non-linearity
        implicitOperator_(solution_, update_);
        rhs -= update_; // rhs = rhs - A u ...
 
        // solve system. Zero initial values for the update_ vector are
        //  ok, because the information about the previous solution already
        // is contained in residual_.
        update_.clear();
        solver(rhs, update_);
        converged_ = (static_cast<int>(solver.iterations()) < solverParams_.maxIterations());
 
        // subtract the update, difference should be the solution.
        solution_ += update_; // ... so u = u + invA(rhs - Au)
      }
 
     public:
 
      virtual bool mark(const double tolerance)
      {
        return markingStrategy_.mark(tolerance);
      }
 
      virtual double estimate()
      {
        return estimator_.estimate(solution_);
      }
 
      virtual double timeEstimate()
      {
        return estimator_.timeEstimate();
      }
 
      virtual double initialEstimate()
      {
        return initialEstimator_.estimate(solution_);
      }
 
      virtual bool initialMarking(const double tolerance)
      {
        return initialMarkingStrategy_.mark(tolerance);
      }
 
      virtual void adapt()
      {
        // there can only one adaptation manager per grid, and the
        // RestrictionProlongationType which determines which
        // functions are restricted/prolongated is built into the
        // AdaptationManagerType. In order to allow for override by
        // derived classes, allocation that adaptationManager_
        // dynamically (and thus not at all, if ThisType::adapt() is
        // never called.
 
        if (!adaptationManager_) {
          // restriction/prolongation operator
          restrictProlong_ = new RestrictionProlongationType(oldSolution_);
 
          // adaptation manager
          adaptationManager_ = new AdaptationManagerType(grid_, *restrictProlong_);
        }
 
        // adaptation and load balancing
        adaptationManager_->adapt();
        if (adaptationManager_->loadBalance()) {
          // TODO: re-initialize stuff as needed.
        }
      }
 
      virtual int output()
      {
        // write checkpoint in order to be able to restart after a crash
        checkPointer_.write(timeProvider_);
 
        if (!dataOutput_) {
          // NOTE: this should allocated dynamically, otherwise a
          // derived class has problems to define completely different
          // IO-schemes Also DataOutputType likes to already generate
          // some files during construction, so make sure this never happens
          dataOutput_ = new DataOutputType(grid_, ioTuple_, timeProvider_, DataOutputParameters());
        }
 
        if (!dataOutput_->willWrite(timeProvider_)) {
          return -1;
        }
 
        dataOutput_->write(timeProvider_);
 
        return dataOutput_->writeStep() - 1;
      }
 
      virtual bool converged() const
      {
        return converged_;
      }
 
      virtual double residual() const
      {
        assert(0);
        return -1;
      }
 
      virtual double error() const
      {
        // can also use H1-norm
        typedef Dune::Fem::L2Norm<GridPartType> NormType;
 
        // the DiscreteFunctionSpaceAdapter sets the order per default
        // to 111 == \infty. We set therefore the order just high
        // enough that we see the asymptotic error. If the polynomial
        // order is D, then the error should decay with (h^D). In
        // principle it should thus suffice to use a quadrature of
        // degree D.
        NormType norm(gridPart_, 2*discreteSpace_.order()+2);
        return norm.distance(initialValue_, solution_);
      }
 
      virtual size_t size() const
      {
        // NOTE: slaveDofs.size() is always one TOO LARGE.
        size_t numberOfDofs = discreteSpace_.blockMapper().size() - discreteSpace_.slaveDofs().size() + 1;
 
        numberOfDofs = grid_.comm().sum(numberOfDofs);
 
        return numberOfDofs;
      }
 
     protected:
      GridType&     grid_; // hierarchical grid
      const TimeProviderType& timeProvider_; // maybe change to TimeView later ...
      DiscreteModelType implicitModel_; // new time-step
      ExplicitModelType explicitModel_; // old time-step
      std::string name_;
 
      GridPartType& gridPart_; // grid part(view), here the leaf grid the discrete space is build with
 
      const DiscreteFunctionSpaceType& discreteSpace_; // discrete function space
      InitialValueType initialValue_;
      RHSFunctionalType rhsFunctional_;
      DiscreteFunctionType& solution_;   // the unknown
      DiscreteFunctionType oldSolution_;   // the solution from the previous timestep
      DiscreteFunctionType residual_;   // initial residual
      DiscreteFunctionType update_;     // distance to solution
 
      EstimatorType        estimator_;  // residual error estimator
      MarkingStrategyType markingStrategy_;
 
      InitialEstimatorType initialEstimator_;
      InitialMarkingStrategyType initialMarkingStrategy_;
 
      RestrictionProlongationType *restrictProlong_; // local restriction/prolongation object
      AdaptationManagerType  *adaptationManager_;    // adaptation manager handling adaptation
 
      ConstraintsOperatorType  constraints_; // dirichlet boundary constraints
 
      ImplicitOperatorType implicitOperator_; // implicit non-linear operator
      ExplicitOperatorType explicitOperator_; // explicit operator
 
      LinearOperatorType linearOperator_;  // the linear operator (i.e. jacobian of the operator)
 
      LinearSolverParameterType solverParams_; // the Solver Parameters from Dune::Fem
      bool converged_ = false; // check whether solver has converged
 
      EmptyConstraintsType emptyBlockConstraints_;
 
      mutable int sequence_;            // sequence number
 
      IOTupleType ioTuple_; // tuple with pointers
      DataOutputType *dataOutput_; // data output class
      CheckPointerType checkPointer_;
 
    };
 
    template<class DiscreteFunction,
             class TimeProvider,
             class ImplicitModel,
             class ExplicitModel,
             class InitialValue,
             class RHSFunctional,
             template<class> class QuadratureTraits>
    std::enable_if_t<IsLinearFunctional<RHSFunctional>::value,
      ParabolicFemScheme<DiscreteFunction,
                         TimeProvider,
                         ImplicitModel,
                         ExplicitModel,
                         InitialValue,
                         RHSFunctional,
                         QuadratureTraits> >
    parabolicFemScheme(DiscreteFunction& solution,
                       const TimeProvider& timeProvider,
                       const ImplicitModel& implicitModel,
                       const ExplicitModel& explicitModel,
                       const InitialValue& initialValue,
                       const RHSFunctional& rhsFunctional,
                       const QuadratureTraits<typename DiscreteFunction::DiscreteFunctionSpaceType::GridPartType>& quadTraits,
                       const std::string& name = "acfem.schemes.parabolic")
    {
      typedef ParabolicFemScheme<DiscreteFunction,
                                 TimeProvider,
                                 ImplicitModel,
                                 ExplicitModel,
                                 InitialValue,
                                 RHSFunctional,
                                 QuadratureTraits>
        ReturnType;
 
      return ReturnType(solution,
                        timeProvider,
                        implicitModel,
                        explicitModel,
                        initialValue,
                        rhsFunctional,
                        name);
    }
 
    template<class DiscreteFunction,
             class TimeProvider,
             class ImplicitModel,
             class ExplicitModel,
             class InitialValue,
             template<class> class QuadratureTraits>
    auto
    parabolicFemScheme(DiscreteFunction& solution,
                       const TimeProvider& timeProvider,
                       const ImplicitModel& implicitModel,
                       const ExplicitModel& explicitModel,
                       const InitialValue& initialValue,
                       const QuadratureTraits<typename DiscreteFunction::DiscreteFunctionSpaceType::GridPartType>& quadTraits,
                       const std::string& name = "acfem.schemes.parabolic")
    {
      typedef ParabolicFemScheme<DiscreteFunction,
                                 TimeProvider,
                                 ImplicitModel,
                                 ExplicitModel,
                                 InitialValue,
                                 ZeroLinearFunctional<typename DiscreteFunction::DiscreteFunctionSpaceType>,
                                 QuadratureTraits>
        ReturnType;
 
      return ReturnType(solution,
                        timeProvider,
                        implicitModel,
                        explicitModel,
                        initialValue,
                        zeroFunctional(solution.space()),
                        name);
    }
 
    template<class DiscreteFunction,
             class TimeProvider,
             class ImplicitModel,
             class ExplicitModel,
             class InitialValue,
             class RHSFunctional>
    std::enable_if_t<IsLinearFunctional<RHSFunctional>::value,
      ParabolicFemScheme<DiscreteFunction,
                         TimeProvider,
                         ImplicitModel,
                         ExplicitModel,
                         InitialValue,
                         RHSFunctional,
                         DefaultQuadratureTraits> >
    parabolicFemScheme(DiscreteFunction& solution,
                       const TimeProvider& timeProvider,
                       const ImplicitModel& implicitModel,
                       const ExplicitModel& explicitModel,
                       const InitialValue& initialValue,
                       const RHSFunctional& rhsFunctional,
                       const std::string& name = "acfem.schemes.parabolic")
    {
      typedef ParabolicFemScheme<DiscreteFunction,
                                 TimeProvider,
                                 ImplicitModel,
                                 ExplicitModel,
                                 InitialValue,
                                 RHSFunctional,
                                 DefaultQuadratureTraits>
        ReturnType;
 
      return ReturnType(solution,
                        timeProvider,
                        implicitModel,
                        explicitModel,
                        initialValue,
                        rhsFunctional,
                        name);
    }
 
    template<class DiscreteFunction,
             class TimeProvider,
             class ImplicitModel,
             class ExplicitModel,
             class InitialValue>
    auto
    parabolicFemScheme(DiscreteFunction& solution,
                       const TimeProvider& timeProvider,
                       const ImplicitModel& implicitModel,
                       const ExplicitModel& explicitModel,
                       const InitialValue& initialValue,
                       const std::string& name = "acfem.schemes.parabolic")
    {
      typedef ParabolicFemScheme<DiscreteFunction,
                                 TimeProvider,
                                 ImplicitModel,
                                 ExplicitModel,
                                 InitialValue,
                                 ZeroLinearFunctional<typename DiscreteFunction::DiscreteFunctionSpaceType>,
                                 DefaultQuadratureTraits>
        ReturnType;
 
      return ReturnType(solution,
                        timeProvider,
                        implicitModel,
                        explicitModel,
                        initialValue,
                        zeroFunctional(solution.space()),
                        name);
    }
 
 
 
  } // ACFem
 
} // Dune
 
#endif // __PARABOLIC_FEMSCHEME_HH__