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/adaptmanager.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>
 
// blah blah blah if only Dune were as hip as it ought to be
 
// local includes
#include "../algorithms/femschemeinterface.hh"
#include "../models/modelinterface.hh"
#include "../functions/gridfunctionwrapper.hh"
#include "../operators/ellipticoperator.hh"
#include "../operators/functionals/l2innerproductfunctional.hh"
#include "../operators/functionals/l2boundaryfunctional.hh"
#include "../operators/functionals/functionalexpression.hh"
#include "../operators/constraints/dirichletconstraints.hh"
#include "../operators/constraints/emptyblockconstraints.hh"
#include "../estimators/parabolicestimator.hh"
#include "../estimators/trueerrorestimator.hh"
#include "../algorithms/marking.hh"
#include "../models/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,
             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;
 
      typedef ImplicitModel ImplicitModelImplementationType;
      typedef ModelInterface<ImplicitModelImplementationType> ImplicitModelType;
      typedef typename ImplicitModelType::OperatorPartsType ImplicitOperatorPartsType;
      typedef ModelConstituents<ImplicitModelType> ImplicitModelConstituents;
      typedef ExplicitModel ExplicitModelImplementationType;
      typedef ModelInterface<ExplicitModelImplementationType> ExplicitModelType;
      typedef typename ExplicitModelType::OperatorPartsType ExplicitOperatorPartsType;
      typedef ModelConstituents<ExplicitModelType> ExplicitModelConstituents;
 
      typedef
      decltype(std::declval<ImplicitModelType>() + std::declval<ExplicitModelType>())
        SumModelType;
      typedef ModelConstituents<SumModelType> SumModelConstituents;
 
      typedef InitialValue InitialValueType;
      typedef
      typename GridFunctionConverter<InitialValueType, GridPartType>::WrappedGridFunctionType
      InitialValueFunctionType;
 
      typedef typename ImplicitModelType::FunctionSpaceType FunctionSpaceType;
 
      // or: (must coincide)
      //typedef typenema DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
 
      // choose type of discrete function, Matrix implementation and solver implementation
      typedef SolverSelector<DiscreteFunctionType, ImplicitOperatorPartsType> SolverSelectorType;
      typedef typename SolverSelectorType::LinearOperatorType LinearOperatorType;
      typedef typename SolverSelectorType::LinearInverseOperatorType LinearInverseOperatorType;
 
      typedef
      Fem::NewtonInverseOperator<LinearOperatorType, LinearInverseOperatorType>
      NonLinearInverseOperatorType;
 
      typedef Dune::Fem::RestrictProlongDefault<DiscreteFunctionType> RestrictionProlongationType;
 
      typedef Dune::Fem::AdaptationManager<GridType, RestrictionProlongationType> AdaptationManagerType;
 
      typedef typename ImplicitModelType::DirichletIndicatorType DirichletIndicatorType;
 
      typedef typename ImplicitModelType::DirichletBoundaryFunctionType DirichletBoundaryFunctionType;
 
      typedef typename ImplicitModelType::DirichletWeightFunctionType DirichletWeightFunctionType;
 
      typedef
      decltype(std::declval<DirichletBoundaryFunctionType>()
               /
               std::declval<typename DirichletWeightFunctionType::GridFunctionType>())
        EffectiveDirichletFunctionType;
 
      typedef DirichletConstraints<LinearOperatorType, EffectiveDirichletFunctionType, DirichletIndicatorType> ConstraintsOperatorType;
 
      typedef DifferentiableEmptyBlockConstraints<LinearOperatorType>  EmptyConstraintsType;
 
      typedef DifferentiableEllipticOperator<LinearOperatorType, ImplicitOperatorPartsType, ConstraintsOperatorType, QuadratureTraits> ImplicitOperatorType;
 
      typedef EllipticOperator<ExplicitOperatorPartsType,
                               DiscreteFunctionType, DiscreteFunctionType,
                               EmptyConstraintsType, QuadratureTraits> ExplicitOperatorType;
 
      // Isn't ther a more simpler way for doing this?
      typedef typename SumModelType::BulkForcesFunctionType BulkForcesFunctionType;
      typedef
      L2InnerProductFunctional<DiscreteFunctionSpaceType, BulkForcesFunctionType>
      BulkForcesFunctional;
      typedef typename SumModelType::NeumannBoundaryFunctionType BoundaryFluxFunctionType;
      typedef
      L2BoundaryFunctional<DiscreteFunctionSpaceType,
                           BoundaryFluxFunctionType,
                           typename ImplicitModelType::NeumannIndicatorType>
      BoundaryFluxFunctional;
 
      typedef
      ParabolicEulerEstimator<DiscreteFunctionType, TimeProviderType,
                              ImplicitModelImplementationType,
                              ExplicitModelImplementationType>
      EstimatorType;
 
      typedef MarkingStrategy<EstimatorType> MarkingStrategyType;
 
      typedef
      TrueErrorEstimator<InitialValueFunctionType, Dune::Fem::L2Norm<GridPartType> >
      InitialEstimatorType;
 
      typedef MarkingStrategy<InitialEstimatorType> InitialMarkingStrategyType;
 
      typedef InitialValueFunctionType ExactSolutionFunctionType; // slight abuse
 
      typedef std::tuple<DiscreteFunctionType *, const ExactSolutionFunctionType *> 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 std::string name = "heat")
        : grid_(solution.space().gridPart().grid()),
          timeProvider_(timeProvider),
          implicitModel_(implicitModel),
        explicitModel_(explicitModel),
        name_(name),
        gridPart_(solution.space().gridPart()),
        discreteSpace_(solution.space()),
        initialValue_(wrapToGridFunction("Initial Value", initialValue, gridPart_, discreteSpace_.order()+1)),
        solution_(solution),
        oldSolution_("old solution", discreteSpace_),
        residual_("residual", discreteSpace_),
        update_("update", discreteSpace_),
        estimator_(timeProvider_, implicitModel_, explicitModel_, oldSolution_),
        markingStrategy_(estimator_, name_ + ".adaptation.space"),
        initialEstimator_(initialValue_()),
        initialMarkingStrategy_(initialEstimator_, name_ + ".adaptation" + ".initial"),
        // restriction/prolongation operator
        restrictProlong_(0),
        // adaptation manager
        adaptationManager_(0),
        // create Dirichlet contraints
        boundaryValues_(implicitModel_.dirichletBoundaryFunction(gridPart_)),
        boundaryWeight_(implicitModel_.dirichletWeightFunction(gridPart_)),
        dirichletConstraints_(discreteSpace_,
                              boundaryValues_ / boundaryWeight_.function(),
                              implicitModel_.dirichletIndicator()),
        // the elliptic operator (implicit)
        implicitOperator_(implicitModel_.operatorParts(), dirichletConstraints_),
        explicitOperator_(explicitModel_.operatorParts()),
        // create linear operator (domainSpace,rangeSpace)
        linearOperator_("assembled elliptic operator", discreteSpace_, discreteSpace_),
        solverEps_(Dune::Fem::Parameter::getValue<double>(name_ + ".solvereps", 1e-8)),
        // the right-hand-side
        sumModel_(implicitModel_ + explicitModel_),
        bulkForces_(sumModel_.bulkForcesFunction(gridPart_)),
        boundaryFlux_(sumModel_.neumannBoundaryFunction(gridPart_)),
        // other stuff
        sequence_(-1),
        ioTuple_(&solution_, &initialValue_()),
        dataOutput_(0),
        checkPointer_(grid_, timeProvider_)
      {
        // Make the solution of the current timestep persistent for check-pointing.
        Dune::Fem::persistenceManager << solution_;
 
        oldSolution_.clear();
        solution_.clear();
      }
 
      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 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 (ImplicitModelConstituents::hasDirichletBoundary) {
          // set Dirichlet constraints, new time
          dirichletConstraints_.rebuildValues(); // force update
          dirichletConstraints_.constrain(solution_);
        }
 
        // If we have a RHS _function_ then also take that into
        // account. Note that it would be comparatively cheap to
        // omit the enclosing if-claus, because the default
        // configuration results into an efficient no-op. There is
        // also no point to optimize further into the case
        // Forces-but-no-Neumann or Neumann-but-no-forces for the
        // same reason.
        auto rhsFunctional =
          BulkForcesFunctional(discreteSpace_, bulkForces_)
          +
          BoundaryFluxFunctional(discreteSpace_, boundaryFlux_, implicitModel_.neumannIndicator())
          +
          sumModel_.forcesFunctional(discreteSpace_);
 
        if (!isZero(rhsFunctional)) { // avoid copying thousands of zero values
          rhsFunctional.coefficients(residual_);
        } else {
          residual_.clear();
        }
 
        // Compute the contribution from the time derivative
        explicitOperator_(oldSolution_, update_);
        residual_ -= update_;
 
        // do no spoil the new Dirichlet-values
        dirichletConstraints_.zeroConstrain(residual_);
 
        if (ImplicitOperatorPartsType::isLinear) {
          linearSolve(residual_, forceMatrixAssembling);    // sove Lu = res
        } else {
          nonLinearSolve(residual_); // solve Lu = res
        }
      }
 
     protected:
 
      virtual void nonLinearSolve(DiscreteFunctionType& rhs)
      {
        assert(!ImplicitOperatorPartsType::isLinear);
 
        NonLinearInverseOperatorType newton(implicitOperator_);
 
        newton(rhs, solution_);
 
        assert(newton.converged());
      }
 
      virtual void linearSolve(DiscreteFunctionType& rhs, bool forceMatrixAssembling)
      {
        assert(ImplicitOperatorPartsType::isLinear);
 
        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(linearOperator_, solverEps_, solverEps_);
 
        // 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_);
 
        // 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 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 ...
      const ImplicitModelType& implicitModel_; // new time-step
      const ExplicitModelType& explicitModel_; // old time-step
      const std::string name_;
 
      GridPartType& gridPart_; // grid part(view), here the leaf grid the discrete space is build with
 
      const DiscreteFunctionSpaceType& discreteSpace_; // discrete function space
      ExpressionStorage<InitialValueFunctionType> initialValue_;  // this is already the GRID function
      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
 
      DirichletBoundaryFunctionType boundaryValues_;
      DirichletWeightFunctionType boundaryWeight_;
      ConstraintsOperatorType  dirichletConstraints_; // dirichlet boundary constraints
 
      ImplicitOperatorType implicitOperator_; // the affine-linear operator.
      ExplicitOperatorType explicitOperator_;
 
      LinearOperatorType linearOperator_;  // the linear operator (i.e. jacobian of the operator)
 
      const double solverEps_ ; // eps for linear solver
 
      DifferentiableEmptyBlockConstraints<LinearOperatorType> emptyBlockConstraints_;
 
      SumModelType sumModel_;
      BulkForcesFunctionType bulkForces_;
      BoundaryFluxFunctionType boundaryFlux_;
 
      mutable int sequence_;            // sequence number
 
      IOTupleType ioTuple_; // tuple with pointers
      DataOutputType *dataOutput_; // data output class
      CheckPointerType checkPointer_;
 
    };
 
 
  } // ACFem
 
} // Dune
 
#endif // __PARABOLIC_FEMSCHEME_HH__