embedded_errorest.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-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef EMBEDDED_ERROREST_HH
#define EMBEDDED_ERROREST_HH
 
#include <type_traits>
#include "boost/multi_array.hpp"
 
#include "fem/errorest.hh"
#include "fem/fixfusion.hh"
#include "fem/iterate_grid.hh"
#include "utilities/power.hh"
#include "utilities/timing.hh"
 
namespace Kaskade
{
  namespace EmbeddedErrorestimatorDetail
  {
 
    using namespace boost::fusion;
 
    struct GetProjections
    {
      template <class Evaluator>
      typename Evaluator::Space::Mapper::ShapeFunctionSet::Matrix operator()(Evaluator const& evaluator) const
      {
        typename Evaluator::Space::Mapper::ShapeFunctionSet::Matrix p = evaluator.shapeFunctions().hierarchicProjection();
        evaluator.combiner().rightTransform(p);
        evaluator.combiner().leftPseudoInverse(p);
        return p;
      }
    };
 
 
    template <class Pairs>
    struct ApplyProjections
    {
      ApplyProjections(Pairs const& data_): data(data_) {}
 
      template <class Pair>
      void operator()(Pair const& var) const
      {
        int const sid = std::remove_reference<typename boost::fusion::result_of::value_at_c<Pair,0>::type>::type::spaceIndex;
 
        typedef typename boost::fusion::result_of::value_at_c<Pairs,sid>::type DataPair;
        typename boost::fusion::result_of::value_at_c<DataPair,0>::type evaluator = at_c<0>(at_c<sid>(data));
        typename boost::fusion::result_of::value_at_c<DataPair,1>::type p = at_c<1>(at_c<sid>(data));
 
        typedef typename std::remove_reference<typename boost::fusion::result_of::value_at_c<Pair,1>::type>::type Function;
        Function& f = at_c<1>(var);
        typename Function::StorageType x(p.N());
 
        for (int i=0; i<p.N(); ++i) {
          x[i] = 0;
          for (int j=0; j<p.M(); ++j)
            x[i].axpy(p[i][j][0][0],f.coefficients()[evaluator.globalIndices()[j]]);
        }
        for (int i=0; i<p.N(); ++i) {
          f.coefficients()[evaluator.globalIndices()[i]] = x[i];
        }
      }
 
    private:
      Pairs const& data;
    };
 
    template <class Sequence>
    ApplyProjections<Sequence> applyProjections(Sequence const& data)
    {
      return ApplyProjections<Sequence>(data);
    }
 
 
 
  }; // end of namespace EmbeddedErrorestimatorDetail
 
  template <class VariableSetDescription>
  void projectHierarchically(VariableSet<VariableSetDescription>& f)
  {
    using namespace boost::fusion;
    using namespace EmbeddedErrorestimatorDetail;
 
    // Step through all cells.
    typedef typename VariableSetDescription::Spaces    Spaces;
    typedef typename SpaceType<Spaces,0>::type::Scalar Scalar;
    typedef typename VariableSetDescription::Grid      Grid;
    typedef typename VariableSetDescription::GridView  GridView;
 
    GridView const& gridView = f.descriptions.gridView ;
 
    // Evaluators for shape functions of all FE spaces. Remember that
    // every thread has to use its own collection of evaluators!
    typedef ShapeFunctionCache<Grid,Scalar> SfCache;
    auto evaluators = getEvaluators(f.descriptions.spaces,static_cast<SfCache*>(nullptr));
    using Evaluators = decltype(evaluators);
    
 
    // A sequence of projection matrices.
    typedef typename boost::fusion::result_of::as_vector<
                  typename boost::fusion::result_of::transform<Evaluators,GetProjections>::type
             >::type Projections;
    Projections projections;
 
    // Step through all cells. On each cell, compute K^+ P K a, where a
    // is the coefficient vector, K the combiner (often the identity),
    // and P the projection matrix that projects to a subspace of lower
    // order shape functions.
    for (auto ci=gridView.template begin<0,Dune::All_Partition>();
         ci!=gridView.template end<0,Dune::All_Partition>(); ++ci) {
      // for all spaces involved, compute the shape functions and
      // their global indices, which are needed for evaluating the
      // functional's derivative.
      moveEvaluatorsToCell(evaluators,*ci);
 
      // Compute projections.
      projections = transform(evaluators,GetProjections());
 
      // Apply projections.
      typename VariableSetDescription::Variables vars;
      for_each2(vars,f.data,
                applyProjections(zip(evaluators,projections)));
    }
  }
 
  //---------------------------------------------------------------------
  // Embedded error estimation.
  namespace EmbeddedErrorestDetail 
  {
 
    // A policy class for the GroupedSummationCollector. This declares each cell
    // a group of its own.
    class CellByCell
    {
    public:
      template <class IndexSet>
      CellByCell(IndexSet const& is): nGroups(is.size(0)) {}
      
      size_t operator[](size_t idx) const { assert(idx<nGroups); return idx; }
      
      size_t nGroups;
    };
 
    // --------------------------------------------------------------------------------------------
 
    template <class F>
    class MaximumCollector
    {
    public:
      MaximumCollector(F const& f_, size_t n_)
      : f(f_)
      , n(n_)
      {}
 
      template <class Cell>
      int integrationOrder(Cell const&, int shapeFunctionOrder) const
      {
        return shapeFunctionOrder;
      }
 
      void join(MaximumCollector& c)
      {
        if (eps.empty())                          // no own data
          eps = std::move(c.eps);                 // -> just move here
        else if (!c.eps.empty())                  // both data exist
        {                                         // -> merge both
          assert(eps.size()==c.eps.size());       //    if coming from partitioned
          for (size_t i=0; i<eps.size(); ++i)     //    multithreading, they are complementary
            eps[i] = std::max(eps[i],c.eps[i]);   //    but for generality we extend taking the
        }                                         //    maximum nevertheless
      }                                           // own data but no other data -> nothing to do
 
      template <class Cell, class Index, class LocalPosition, class Values>
      void operator()(Cell const& cell, Index idx, LocalPosition const& /* xi */,
                      double /* weight */, Values const& u)
      {
        if (eps.empty())
          eps.resize(n);
        eps[idx] = std::max(eps[idx],double(f(u)));
      }
 
      std::vector<double> const& epsilon() const
      {
        return eps;
      }
 
    private:
      F f;
      size_t n;                       // number of cells
      std::vector<double> eps;    // cell maximum
    };
  } // End of namespace EmbeddedErrorestDetail
 
  template <class VariableSet, class GridManager>
  std::array<size_t,3>
  uniformEmbeddedErrorEstimation(VariableSet const& u,
                                 Dune::FieldVector<double,VariableSet::Descriptions::noOfVariables> const& tol,
                                 GridManager& gridManager,
                                 double rhoMin=0.1, double rhoMax=0.9)
  {
    ScopedTimingSection timing("uniform embedded error estimation");
 
    auto const& varDesc = u.descriptions;
 
    timing.timer().start("project hierarchically");
    auto du = u;                                  // Get the embedded error estimate by projecting
    projectHierarchically(du);                    // to lower order and taking the difference.
    du -= u;                                      // (Note that the sign doesn't matter, we use norms.)
    timing.timer().stop("project hierarchically");
 
    timing.timer().start("quantify error");
    auto const n = varDesc.indexSet.size(0);      // number of cells
    auto f = [&tol](auto const& us)               // This maps u to max |u_i|/tol_i.
    {
      int i = 0;
      double epsilon = 0;
      boost::fusion::for_each(us,[&](auto const& u)
                                 { epsilon = std::max(epsilon,u.two_norm() / tol[i++]); });
      return epsilon;
    };
 
    EmbeddedErrorestDetail::MaximumCollector max(f,n);  // Find the weighted maximum of all errors
    gridIterate(du,max);                                // over the cells.
    timing.timer().stop("quantify error");
 
    timing.timer().start("mark");
    std::array<size_t,3> counts{0,0,0};
    for (auto const& cell: Dune::elements(varDesc.gridView))
    {
      double e = max.epsilon()[varDesc.indexSet.index(cell)];
      gridManager.mark(e>rhoMax? 1: e<rhoMin? -1: 0, cell);
      if (e > 1)
        ++counts[0];
      if (e > rhoMax)
        ++counts[1];
      else if (e < rhoMin)
        ++counts[2];
    }
    timing.timer().stop("mark");
 
    return counts;
  }
 
  // ----------------------------------------------------------------------------------------------
 
  template <class VariableSetDescription, class Scaling=IdentityScaling>
  class EmbeddedErrorEstimator
  {
  public:
    typedef GridManager<typename VariableSetDescription::Grid> GManager;
    typedef EmbeddedErrorEstimator<VariableSetDescription,Scaling> Self;
    
    EmbeddedErrorEstimator(GManager& gridManager_, VariableSetDescription const& varDesc_, Scaling const& scaling_=Scaling())
    : gridManager(gridManager_), varDesc(varDesc_), tol(varDesc.noOfVariables,std::make_pair(0.0,0.01)), scaling(scaling_) {}
 
 
    bool estimate(typename VariableSetDescription::VariableSet const& err,
                  typename VariableSetDescription::VariableSet const& sol) const
    {
      using namespace boost::fusion;
      using namespace EmbeddedErrorestDetail;
      
      assert(size(err.data)==size(sol.data));
 
      ErrorestDetail::GroupedSummationCollector<CellByCell> sum(CellByCell(varDesc.indexSet));
      scaledTwoNormSquared(join(typename VariableSetDescription::Variables(),typename VariableSetDescription::Variables()),
                           join(err.data,sol.data),varDesc.spaces,scaling,sum);
 
      int const s = varDesc.noOfVariables;
      int const n = sum.sums.shape()[0];
      
      assert(sum.sums.shape()[1]==2*s);
      
      // Compute scaled L2 norms of solution and errors.
      std::vector<double> norm2(s), error2(s);
      for (int idx=0; idx<n; ++idx)
        for (int j=0; j<s; ++j) {
          error2[j] += sum.sums[idx][j];
          norm2[j]  += sum.sums[idx][j+s];
        }
        
        
      if (verbosity>0)
      {
        std::cout << "   ||solution||^2 = " << std::scientific << std::setprecision(3); std::copy(norm2.begin(),norm2.end(),std::ostream_iterator<double>(std::cout," "));
        std::cout << std::endl;
        std::cout << "     ||errors||^2 = "; std::copy(error2.begin(),error2.end(),std::ostream_iterator<double>(std::cout," "));
        std::cout << std::resetiosflags( ::std::ios::scientific ) << std::setprecision(6) << std::endl;
      };
        
      // Compute error limit and normalized errors
      std::vector<double> errorLimit(s);
      boost::multi_array<double,2> normalizedErrors(boost::extents[n][s]);
        
      for (int j=0; j<s; ++j) 
      {
        errorLimit[j] = square(tol[j].first) + square(tol[j].second)*norm2[j];
        for (int i=0; i<n; ++i) 
          normalizedErrors[i][j] = sum.sums[i][j] / errorLimit[j];
      }
        
      // Get refinement thresholds
      FixedFractionCriterion marker(0.1);
      auto threshold = marker.threshold(normalizedErrors);
        
      // Check for refinement necessity
      if (threshold.empty()) // all variables accurate enough
        return true;
      
      // Refine the mesh
      markAndRefine(gridManager,varDesc.gridView,normalizedErrors,threshold);
      return false;
    }
    
    Self& setTolerances(std::vector<std::pair<double,double> > const& tol_)
    {
      tol = tol_;
      return *this;
    }
    
    Self& setScaling(Scaling const& scaling_)
    {
      scaling = scaling_;
      return *this;
    }
    
    Self& setCriterion(RefinementCriterion const& criterion_)
    {
      criterion = &criterion_;
      return *this;
    }
    
    Self& setVerbosity(int verbosity_)
    {
      verbosity = verbosity_;
      return *this;
    }
    
    
  private:
    GManager&                              gridManager;
    VariableSetDescription const&          varDesc;
    std::vector<std::pair<double,double> > tol;
    Scaling                                scaling;
    RefinementCriterion const*             criterion;
    int                                    verbosity;
  };
 
  template <class VariableSetDescription, class Scaling>
  bool embeddedErrorEstimator(VariableSetDescription const& varDesc,
                              typename VariableSetDescription::VariableSet const& err,
                              typename VariableSetDescription::VariableSet const& sol,
                              Scaling const& scaling,
                              std::vector<std::pair<double,double> > const& tol,
                              GridManager<typename VariableSetDescription::Grid>& gridManager,
                              int verbosity=1)
  {
    EmbeddedErrorEstimator<VariableSetDescription,Scaling> estimator(gridManager,varDesc,scaling);
    FixedFractionCriterion marker(0.1);
    estimator.setTolerances(tol).setCriterion(marker).setVerbosity(verbosity);
    return estimator.estimate(err,sol);
  }
} /* end of namespace Kaskade */
 
#endif