Marker.inc.hpp

#include <cmath>

namespace AMDiS {

template <class Traits>
std::unique_ptr<Marker<Traits> > Marker<Traits>::
createMarker(std::string const& name, std::string const& component,
             Estimates const& est, std::shared_ptr<Grid> const& grid)
{
  int strategy = 0;
  Parameters::get(name + "->strategy", strategy);

  switch (strategy) {
  case 0:  // no refinement/coarsening
    break;
  case 1:
    msg("Creating global refinement (GR) marker");
    return std::make_unique<GRMarker<Traits> >(name, component, est, grid);
    break;
  case 2:
    msg("Creating maximum strategy (MS) marker");
    return std::make_unique<MSMarker<Traits> >(name, component, est, grid);
    break;
  case 3:
    msg("Creating equidistribution strategy (ES) marker");
    return std::make_unique<ESMarker<Traits> >(name, component, est, grid);
    break;
  case 4:
    msg("Creating guaranteed error reduction strategy (GERS) marker");
    return std::make_unique<GERSMarker<Traits> >(name, component, est, grid);
    break;
  default:
    error_exit("invalid strategy");
  }

  return {};
}


template <class Traits>
void Marker<Traits>::initMarking(AdaptInfo& adaptInfo)
{
  elMarkRefine_ = 0;
  elMarkCoarsen_ = 0;
  estSum_ = std::pow(adaptInfo.getEstSum(component_), p_);
  estMax_ = adaptInfo.getEstMax(component_);
  refineAllowed_ = adaptInfo.isRefinementAllowed(component_);
  coarsenAllowed_ = adaptInfo.isCoarseningAllowed(component_);
}


template <class Traits>
void Marker<Traits>::finishMarking(AdaptInfo& adaptInfo)
{
  msg("{} elements marked for refinement", elMarkRefine_);
  msg("{} elements marked for coarsening", elMarkCoarsen_);
}


template <class Traits>
void Marker<Traits>::markElement(AdaptInfo& adaptInfo, const Element& elem)
{
  const auto& index = grid_->leafIndexSet().index(elem);
  double lError = est_[index];

  if (lError > markRLimit_ && refineAllowed_
      && (maxRefineLevel_ == -1 || elem.level() < maxRefineLevel_)) {
    this->mark(elem, 1);
  } else  if (lError <= markCLimit_ && coarsenAllowed_
              && (minRefineLevel_ == -1 || elem.level() > minRefineLevel_)) {
    this->mark(elem, -1);
  }
}


template <class Traits>
Flag Marker<Traits>::markGrid(AdaptInfo& adaptInfo)
{
  if (!grid_)
    error_exit("No grid!");

  initMarking(adaptInfo);

  if (!coarsenAllowed_ && !refineAllowed_)
    return 0;

  for (const auto& elem : Dune::elements(grid_->leafGridView()))
    markElement(adaptInfo, elem);

  finishMarking(adaptInfo);

  Flag markFlag;
  if (elMarkRefine_)
    markFlag = 1;
  if (elMarkCoarsen_)
    markFlag |= 2;

  return markFlag;
}


template <class Traits>
void MSMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);

  double msGammaP = std::pow(msGamma_, this->p_);
  double msGammaCP = std::pow(msGammaC_, this->p_);

  this->markRLimit_ = msGammaP * adaptInfo.getEstMax(this->component_);
  this->markCLimit_ = msGammaCP * adaptInfo.getEstMax(this->component_);

  msg("start max_est: {}, mark_limits: {}, {}",
    adaptInfo.getEstMax(this->component_), this->markRLimit_ , this->markCLimit_);
}


template <class Traits>
void ESMarker<Traits>::initMarking(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);

  double esThetaP = std::pow(esTheta_, this->p_);
  double esThetaCP = std::pow(esThetaC_, this->p_);
  double epsP = std::pow(adaptInfo.getSpaceTolerance(this->component_), this->p_);

  int nLeaves = (this->grid_->leafGridView()).size(0);
/*#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
  Parallel::mpi::globalAdd(nLeaves);
#endif*/

  this->markRLimit_ = esThetaP * epsP / nLeaves;
  this->markCLimit_ = esThetaCP * epsP / nLeaves;

  msg("start mark_limits: {}, {}; nt = {}", this->markRLimit_, this->markCLimit_, nLeaves);
}


template <class Traits>
Flag GERSMarker<Traits>::markGrid(AdaptInfo& adaptInfo)
{
  Super::initMarking(adaptInfo);

  if (!this->coarsenAllowed_ && !this->refineAllowed_)
    return 0;

  gersSum_ = 0.0;

  double LTheta = std::pow(1.0 - gersThetaStar_, this->p_);
  double epsP = std::pow(adaptInfo.getSpaceTolerance(this->component_), this->p_);

  if (this->estSum_ < oldErrSum_) {
    double improv = this->estSum_ / oldErrSum_;
    double wanted = 0.8 * epsP / this->estSum_;
    double redfac = std::min((1.0 - wanted) / (1.0 - improv), 1.0);
    redfac = std::max(redfac, 0.0);

    if (redfac < 1.0) {
      LTheta *= redfac;
      msg("GERS: use extrapolated theta_star = {}", std::pow(LTheta, 1.0 / this->p_));
    }
  }

  oldErrSum_ = this->estSum_;
  double GERSGamma = 1.0;

  if (this->refineAllowed_) {
    if (LTheta > 0) {
      do {
        gersSum_ = 0.0;
        GERSGamma -= gersNu_;
        this->markRLimit_ = GERSGamma * this->estMax_;

        for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
          markElementForRefinement(adaptInfo, elem);

      } while((GERSGamma > 0) && (gersSum_ < LTheta * this->estSum_));
    }

    msg("GERS refinement with gamma = {}", GERSGamma);
  }

  if (this->coarsenAllowed_) {
    GERSGamma = 0.3;
    LTheta = gersThetaC_ * epsP;

    do {
      gersSum_ = 0.0;
      GERSGamma -= gersNu_;
      this->markCLimit_ = GERSGamma * this->estMax_;

      for (const auto& elem : Dune::elements(this->grid_->leafGridView()))
        markElementForCoarsening(adaptInfo, elem);

      msg("coarse loop: gamma = {}, sum = {}, limit = {}", GERSGamma, gersSum_, LTheta);
    } while(gersSum_ > LTheta);

    msg("GERS coarsening with gamma = {}", GERSGamma);
  }

  Super::finishMarking(adaptInfo);

  Flag markFlag;
  if (this->elMarkRefine_)
    markFlag = 1;
  if (this->elMarkCoarsen_)
    markFlag |= 2;

  return markFlag;
}


template <class Traits>
void GERSMarker<Traits>::markElementForRefinement(AdaptInfo& adaptInfo, const Element& elem)
{
  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];

  if (lError > this->markRLimit_) {
    gersSum_ += lError;
    this->mark(elem, 1);
  }
}


template <class Traits>
void GERSMarker<Traits>::markElementForCoarsening(AdaptInfo& adaptInfo, const Element& elem)
{
  double lError = this->est_[this->grid_->leafIndexSet().index(elem)];

  if (this->grid_->getMark(elem) <= 0) {
    if (lError <= this->markCLimit_) {
      gersSum_ += lError;
      this->mark(elem, -1);
    } else {
      this->mark(elem, 0);
    }
  }
}


} // end namespace AMDiS