ellipticfemscheme.hh

#ifndef __ELLIPTIC_FEMSCHEME_HH__
#define __ELLIPTIC_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>
 
// lagrange 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 <dune/fem/solver/parameter.hh>
 
// include output
#include <dune/fem/io/file/dataoutput.hh>
 
// local includes
#include "../common/matrixhelper.hh"
#include "../algorithms/femschemeinterface.hh"
#include "../algorithms/operatorselector.hh"
#include "../models/modeltraits.hh"
#include "../operators/ellipticoperator.hh"
#include "../operators/functionals/functionals.hh"
#include "../operators/constraints/dirichletconstraints.hh"
#include "../functions/functions.hh"
 
#include "../estimators/ellipticestimator.hh"
#include "../algorithms/marking.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 Model, class InitialGuess, class RHSFunctional>
    class EllipticFemScheme
      : public AdaptiveFemScheme
    {
     public:
      typedef DiscreteFunction DiscreteFunctionType;
 
      typedef RHSFunctional RHSFunctionalType; 
 
      typedef typename DiscreteFunctionType::GridType GridType;
 
      typedef typename DiscreteFunctionType::GridPartType GridPartType;
 
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
 
      using ModelType = Model;
 
      using DiscreteModelType = DiscreteModel<Model, DiscreteFunctionSpaceType>;
 
      using TraitsType = ModelTraits<DiscreteModelType>;
 
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;
 
      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 Fem::RestrictProlongDefault<DiscreteFunctionType> RestrictionProlongationType;
 
      typedef Fem::AdaptationManager<GridType, RestrictionProlongationType> AdaptationManagerType;
 
      typedef DirichletConstraints<DiscreteFunctionSpaceType, DiscreteModelType> ConstraintsOperatorType;
 
      typedef DifferentiableEllipticOperator<LinearOperatorType, DiscreteModelType, const ConstraintsOperatorType&> OperatorType;
 
      typedef EllipticEstimator<DiscreteFunctionSpaceType, DiscreteModelType> EstimatorType;
 
      typedef MarkingStrategy<GridPartType> MarkingStrategyType;
 
      typedef InitialGuess InitialGuessType;
 
      typedef std::tuple<DiscreteFunctionType *, InitialGuessType *> IOTupleType;
 
      typedef DataOutput<GridType, IOTupleType> DataOutputType;
 
     public:
      EllipticFemScheme(DiscreteFunctionType& solution,
                        const ModelType& model,
                        const InitialGuessType& initialGuess,
                        const RHSFunctionalType& rhsFunctional,
                        const std::string& name)
        : grid_(solution.space().gridPart().grid()),
          model_(asExpression(discreteModel(model, solution.space()))),
          name_(name),
          gridPart_(solution.space().gridPart()),
          discreteSpace_(solution.space()),
          solution_(solution),
          // estimator
          estimator_(model_, discreteSpace_),
          markingStrategy_(gridPart_, estimator_, name_ + ".adaptation"),
          // restriction/prolongation operator
          restrictProlong_(0),
          // adaptation manager
          adaptationManager_(0),
          // create Dirichlet contraints
          constraints_(discreteSpace_, model_),
          // the elliptic operator (implicit)
          operator_(model_, constraints_),
          // create linear operator (domainSpace,rangeSpace)
          linearOperator_("assembled elliptic operator", discreteSpace_, discreteSpace_),
          solverParams_(TraitsType::isAffineLinear ? name_ + ".solver." : name_ + ".linear."),
          sequence_(-1),
          // initial guess or "exact" solution
          initialGuess_(initialGuess),
          // cast back to implementation
          rhsFunctional_(rhsFunctional),
          // io tuple
          ioTuple_(&solution_, &initialGuess_),
          // DataOutputType
          dataOutput_(0)
      {
        chooseDefaultSolverMethod<SolverSelectorType>(solverParams_);
      }
 
      EllipticFemScheme(const EllipticFemScheme& other)
        : grid_(other.grid_),
          model_(other.model_),
          name_(other.name_),
          gridPart_(other.gridPart_),
          discreteSpace_(other.discreteSpace_),
          solution_(other.solution_),
          // estimator
          estimator_(model_, discreteSpace_),
          markingStrategy_(other.markingStrategy_),
          // restriction/prolongation operator
          restrictProlong_(0),
          // adaptation manager
          adaptationManager_(0),
          // create Dirichlet contraints
          constraints_(discreteSpace_, model_),
          // the elliptic operator (implicit)
          operator_(model_, constraints_),
          // create linear operator (domainSpace,rangeSpace)
          linearOperator_("assembled elliptic operator", discreteSpace_, discreteSpace_),
          solverParams_(other.solverParams_),
          sequence_(-1),
          // initial guess or "exact" solution
          initialGuess_(other.initialGuess_),
          // optional additional rhs-functional
          rhsFunctional_(other.rhsFunctional_),
          // io tuple
          ioTuple_(&solution_, &initialGuess_),
          // DataOutputType
          dataOutput_(0)
      {
        chooseDefaultSolverMethod<SolverSelectorType>(solverParams_);
      }
 
     public:
      virtual ~EllipticFemScheme() {
        if (dataOutput_) {
          delete dataOutput_;
        }
        if (adaptationManager_) {
          delete adaptationManager_;
        }
        if (restrictProlong_) {
          delete restrictProlong_;
        }
      }
 
      virtual std::string name () const
      {
        return name_;
      }
 
      virtual void initialize()
      {
        if constexpr (ExpressionTraits<InitialGuessType>::isZero) {
          solution_.clear();
        } else {
          // apply natural interpolation
          interpolate(initialGuess_, solution_);
        }
      }
 
      virtual void assemble() {
        operator_.jacobian(solution_, linearOperator_);
      }
 
      virtual void solve(bool forceMatrixAssembling = true)
      {
        // set Dirichlet constraints in solution
        constraints_.constrain(solution_);
 
        // right hand side, if present
        DiscreteFunctionType rhs("rhs", discreteSpace_);
 
        //The Expression Template eliminate functional
        // whose type evaluates to the ZeroFunctional
        rhs.clear();
        rhs += rhsFunctional_;
        constraints_.zeroConstrain(rhs);
 
        if (TraitsType::isAffineLinear) {
          linearSolve(rhs, forceMatrixAssembling);
        } else {
          nonLinearSolve(rhs);
        }
      }
 
     protected:
 
      virtual void nonLinearSolve(DiscreteFunctionType& rhs)
      {
        assert(!TraitsType::isAffineLinear);
 
        Dune::Fem::NewtonParameter newtonParam(solverParams_, name_ + ".newton.");
        NonLinearInverseOperatorType newton(newtonParam);
        newton.bind(operator_);
 
        newton(rhs, solution_);
 
        converged_ = newton.converged();
        assert(newton.converged());
      }
 
      virtual void linearSolve(DiscreteFunctionType& rhs, bool forceMatrixAssembling)
      {
        assert(TraitsType::isAffineLinear);
 
        // if sequence number is outdated update linear operator
        if (sequence_ != discreteSpace_.sequence() || forceMatrixAssembling) {
          // assemble linear operator (i.e. setup matrix) The jacobian
          // incorporates the constraints defined by the constraint class.
          assemble();
 
          // update sequence number
          sequence_ = discreteSpace_.sequence();
        }
 
        // inverse operator using linear operator
        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
        DiscreteFunctionType update("update", discreteSpace_);
        operator_(solution_, update);
        rhs -= update; // rhs = rhs - A u ...
 
        // Zero initial values for the update_ vector are ok, because
        // the information about the previous solution is already is
        // contained in residual_.
        update.clear();
 
        solver(rhs, update);
        converged_ = (static_cast<int>(solver.iterations()) < solverParams_.maxIterations());
        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 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(solution_);
 
          // adaptation manager
          adaptationManager_ = new AdaptationManagerType(grid_, *restrictProlong_);
        }
 
        // apply adaptation and load balancing
        adaptationManager_->adapt();
        if (adaptationManager_->loadBalance()) {
          // TODO: re-initialize stuff as needed.
        }
      }
 
      virtual int output()
      {
        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_, DataOutputParameters());
        }
 
        if (!dataOutput_->willWrite()) {
          return -1;
        }
 
        // write data
        dataOutput_->write();
 
        return dataOutput_->writeStep() - 1;
      }
 
      virtual bool converged() const
      {
        return converged_;
      }
 
      virtual double residual() const
      {
        // right hand side, if present
        DiscreteFunctionType rhs("rhs", discreteSpace_);
 
        rhs.clear();
        rhs += rhsFunctional_;
        constraints_.zeroConstrain(rhs);
 
        DiscreteFunctionType update("update", discreteSpace_);
        operator_(solution_, update);
        rhs -= update; // rhs = rhs - A u ...
        return std::sqrt(rhs.scalarProductDofs(rhs));
      }
 
      virtual double error() const
      {
        // can also use L2-norm
        typedef Fem::H1Norm< 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 it 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());
        return norm.distance(initialGuess_, solution_);
      }
 
      virtual size_t size() const
      {
        // NOTE: slaveDofs.size() is always one TOO LARGE.
        size_t numberOfDofs = discreteSpace_.size() - discreteSpace_.slaveDofs().size() + 1;
 
        numberOfDofs = grid_.comm().sum(numberOfDofs);
 
        return numberOfDofs;
      }
 
      // TODO: make this const, and assemble somewhere else
      virtual void printSystemMatrix(bool sparse = false, std::ostream& out = std::clog) {
        assemble();
        MatrixHelper<LinearOperatorType> mh(linearOperator_);
        if (sparse) {
          mh.printSparse(out);
        } else {
          mh.printPretty(out);
        }
      }
 
     protected:
      GridType& grid_; // hierarchical grid
      DiscreteModelType model_; // always a copy
 
      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
      DiscreteFunctionType& solution_;   // the unknown
 
      EstimatorType       estimator_;  // residual error estimator
      MarkingStrategyType markingStrategy_;
 
      RestrictionProlongationType *restrictProlong_ ; // local restriction/prolongation object
      AdaptationManagerType  *adaptationManager_ ;    // adaptation manager handling adaptation
 
      ConstraintsOperatorType constraints_; // dirichlet boundary constraints
 
      OperatorType operator_; // the (potentially non-linear) operator
 
      LinearOperatorType linearOperator_;  // the linear operator (i.e. jacobian of the operator)
 
      LinearSolverParameterType solverParams_; // the solver Parameters of Dune::Fem
      bool converged_ = false; // check whether solver has converged
      mutable int sequence_;            // sequence number
 
      InitialGuessType initialGuess_;
      RHSFunctionalType rhsFunctional_;
      IOTupleType ioTuple_; // tuple with pointers
      DataOutputType *dataOutput_; // data output class
    };
 
    template<
      class DiscreteFunction, class Model, class InitialGuess, class RHSFunctional,
      std::enable_if_t<(IsWrappableByConstLocalFunction<InitialGuess>::value
                        && IsLinearFunctional<RHSFunctional>::value
      ), int> = 0>
    auto ellipticFemScheme(DiscreteFunction& solution,
                           const Model&model,
                           const InitialGuess& initialGuess,
                           const RHSFunctional& rhsFunctional,
                           const std::string name = "acfem.schemes.elliptic")
    {
      typedef EllipticFemScheme<DiscreteFunction, Model, InitialGuess, RHSFunctional> ReturnType;
 
      return ReturnType(solution, model, initialGuess, rhsFunctional, name);
    }
 
    template<class DiscreteFunction, class Model, class InitialGuess,
             std::enable_if_t<IsWrappableByConstLocalFunction<InitialGuess>::value, int> = 0>
    auto
    ellipticFemScheme(DiscreteFunction& solution,
                      const Model& model,
                      const InitialGuess& initialGuess,
                      const std::string name = "acfem.schemes.elliptic")
    {
      typedef
        EllipticFemScheme<DiscreteFunction,
                          Model,
                          InitialGuess,
                          ZeroLinearFunctional<typename DiscreteFunction::DiscreteFunctionSpaceType> >
        ReturnType;
 
      return ReturnType(solution,
                        model,
                        initialGuess,
                        zeroFunctional(solution.space()),
                        name);
    }
 
    template<class DiscreteFunction, class Model, class RHSFunctional,
             std::enable_if_t<IsLinearFunctional<RHSFunctional>::value, int> = 0>
    auto ellipticFemScheme(DiscreteFunction& solution,
                           const Model& model,
                           const RHSFunctional& rhsFunctional,
                           const std::string name = "acfem.schemes.elliptic")
    {
      typedef
        EllipticFemScheme<DiscreteFunction,
                          Model,
                          ZeroGridFunction<typename DiscreteFunction::FunctionSpaceType,
                                           typename DiscreteFunction::GridPartType>,
                          RHSFunctional>
        ReturnType;
 
      return ReturnType(solution,
                        model,
                        GridFunction::zeros(solution.space()),
                        rhsFunctional,
                        name);
    }
 
    template<class DiscreteFunction, class Model>
    auto ellipticFemScheme(DiscreteFunction& solution,
                           const Model& model,
                           const std::string name = "acfem.schemes.elliptic")
    {
      typedef
        EllipticFemScheme<DiscreteFunction,
                          Model,
                          ZeroGridFunction<typename DiscreteFunction::FunctionSpaceType,
                                           typename DiscreteFunction::GridPartType>,
                          ZeroLinearFunctional<typename DiscreteFunction::DiscreteFunctionSpaceType> >
        ReturnType;
 
      return ReturnType(solution,
                        model,
                        GridFunction::zeros(solution.space()),
                        zeroFunctional(solution.space()),
                        name);
    }
 
 
 
  } // ACFem::
 
} // Dune::
 
#endif // __ELLIPTIC_FEMSCHEME_HH__