TimeStepping.h

Go to the documentation of this file.

// This file is part of the ExaHyPE2 project. For conditions of distribution and
// use, please see the copyright notice at www.peano-framework.org
#pragma once
 
 
#include "tarch/la/ScalarOperations.h"
#include "tarch/logging/Log.h"
 
#include "peano4/datamanagement/FaceEnumerator.h"
 
#include <utility>
 
 
namespace exahype2 {
  namespace internal {
    template <typename FaceLabel>
    double getMinTimeStampOfNeighbours(
      const peano4::datamanagement::FaceEnumerator< FaceLabel >& faceLabelEnumerator
    ) {
      double result = std::numeric_limits<double>::max();
 
      for (int d=0; d<Dimensions; d++) {
        result = std::min(result, faceLabelEnumerator(d).getNewTimeStamp(0) );
        result = std::min(result, faceLabelEnumerator(d+Dimensions).getNewTimeStamp(1) );
      }
 
      return result;
    }
  }
 
  template <typename CellLabel, typename FaceLabel>
  double getMinTimeStampOfNeighboursAhead(
    const CellLabel& cellLabel,
    const peano4::datamanagement::FaceEnumerator< FaceLabel >& faceLabelEnumerator
  ) {
    double result = std::numeric_limits<double>::max();
    bool   oneIsAhead  = false;
    bool   oneIsBehind = false;
 
    for (int d=0; d<Dimensions; d++) {
      double leftNeighbourTimeStamp = faceLabelEnumerator(d).getNewTimeStamp(0);
      if (
        not faceLabelEnumerator(d).getBoundary()
        and
        tarch::la::greater(leftNeighbourTimeStamp,cellLabel.getTimeStamp())
      ) {
        result     = std::min(result,leftNeighbourTimeStamp);
        oneIsAhead = true;
      }
      if (
        not faceLabelEnumerator(d).getBoundary()
        and
        tarch::la::smaller(leftNeighbourTimeStamp,cellLabel.getTimeStamp())
      ) {
        oneIsBehind = true;
      }
 
      double rightNeighbourTimeStamp = faceLabelEnumerator(d+Dimensions).getNewTimeStamp(1);
      if (
        not faceLabelEnumerator(d+Dimensions).getBoundary()
        and
        tarch::la::greater(rightNeighbourTimeStamp,cellLabel.getTimeStamp())
      ) {
        result     = std::min(result,rightNeighbourTimeStamp);
        oneIsAhead = true;
      }
      if (
        not faceLabelEnumerator(d+Dimensions).getBoundary()
        and
        tarch::la::smaller(rightNeighbourTimeStamp,cellLabel.getTimeStamp())
      ) {
        oneIsBehind = true;
      }
    }
 
    if (oneIsAhead and not oneIsBehind) {
      assertion(result>cellLabel.getTimeStamp());
      return result;
    }
    else {
      return cellLabel.getTimeStamp();
    }
  }
 
  template <typename CellLabel, typename FaceLabel>
  bool runTimeStepOnCell(
    const CellLabel&                                            cellLabel,
    const peano4::datamanagement::FaceEnumerator< FaceLabel >&  faceLabelEnumerator,
    double                                                      globalMinTimeStamp
  ) {
    double cellTimeStamp    = cellLabel.getTimeStamp();
    double cellTimeStepSize = cellLabel.getTimeStepSize();
 
    const double Tolerance = 0.01;
 
    bool result = tarch::la::greaterEquals( ::exahype2::internal::getMinTimeStampOfNeighbours(faceLabelEnumerator), cellTimeStamp, Tolerance*cellTimeStepSize )
               or tarch::la::equals(cellTimeStamp, globalMinTimeStamp);
 
    // @todo das ist net das Problem. Das Problem ist, dass ich u.U. keine Eigenwerte habe
    if (not result) {
//      std::cout << "skip " << cellLabel.toString() << " as neighbours have t_min=" << ::exahype2::internal::getMinTimeStampOfNeighbours(faceLabelEnumerator) << std::endl;
    }
    return result;
  }
 
 
  template <typename CellLabel, typename FaceLabel>
  double removeTimeStepAccumulationErrorsFromCell(
    const CellLabel&  cellLabel,
    const peano4::datamanagement::FaceEnumerator< FaceLabel >&  faceLabelEnumerator,
    double  timeStepSize
  ) {
    static tarch::logging::Log _log( "exahype2" );
 
    const double Tolerance = 0.01;
    double relativeTolerance = tarch::la::equals(timeStepSize,0.0) ? tarch::la::NUMERICAL_ZERO_DIFFERENCE : Tolerance * timeStepSize;
 
    double       maxOfMatchingTimeStamps = -1.0;
    const double currentFutureTimeStamp  = cellLabel.getTimeStamp() + timeStepSize;
    for (int d=0; d<Dimensions; d++) {
      if (
        not faceLabelEnumerator(d).getBoundary()
        and
        tarch::la::equals( faceLabelEnumerator(d).getNewTimeStamp(0), currentFutureTimeStamp, relativeTolerance)
      ) {
        maxOfMatchingTimeStamps = std::max(maxOfMatchingTimeStamps,faceLabelEnumerator(d).getNewTimeStamp(0));
      }
      if (
        not faceLabelEnumerator(d).getBoundary()
        and
        tarch::la::equals( faceLabelEnumerator(d+Dimensions).getNewTimeStamp(1), currentFutureTimeStamp, relativeTolerance)
      ) {
        maxOfMatchingTimeStamps = std::max(maxOfMatchingTimeStamps,faceLabelEnumerator(d+Dimensions).getNewTimeStamp(1));
      }
    }
 
    if ( maxOfMatchingTimeStamps>0.0 ) {
      double tweakedTimeStepSize = maxOfMatchingTimeStamps - cellLabel.getTimeStamp();
      assertion2( tarch::la::greaterEquals( tweakedTimeStepSize,timeStepSize ), tweakedTimeStepSize, timeStepSize );
      if (tweakedTimeStepSize>timeStepSize) {
        logDebug( "removeTimeStepAccummulationErrorsFromCell(...)", "adopt local time step from " << timeStepSize << " to " << tweakedTimeStepSize );
      }
      timeStepSize = tweakedTimeStepSize;
    }
 
    logDebug( "removeTimeStepAccummulationErrorsFromCell(...)", "dt=" << timeStepSize << ", max-match=" << maxOfMatchingTimeStamps );
    return timeStepSize;
  }
 
 
  std::pair<double,double> getInterpolationWeights( double oldTimeStampOnFace, double newTimeStampOnFace, double cellTimeStamp );
 
  double discretiseAndTruncateTimeStepSizes(
    double cellTimeStepSize,
    double maxGlobalTimeStepSize,
    int    discretisationStepsSize
  );
}