LoopBody.cpph

Go to the documentation of this file.

#include "exahype2/Solver.h"
#include "exahype2/fd/PatchUtils.h"
 
#include "tarch/compiler/CompilerSpecificSettings.h"
#include "tarch/Assertions.h"
#include "tarch/multicore/otter.h"
 
#if defined(GPUOffloadingOMP)
#pragma omp declare target
#endif
template <class QInEnumeratorType, class QOutEnumeratorType>
void exahype2::fd::internal::copySolution_LoopBody(
  const double* __restrict__                QIn,
  const QInEnumeratorType&                  QInEnumerator,
  int                                       patchIndex,
  const tarch::la::Vector<Dimensions,int>&  volumeIndex,
  int                                       unknown,
  double* __restrict__                      QOut,
  const QOutEnumeratorType&                 QOutEnumerator
) {
  QOut[ QOutEnumerator(patchIndex,volumeIndex,unknown) ] = QIn[ QInEnumerator(patchIndex,volumeIndex,unknown) ];
}
#if defined(GPUOffloadingOMP)
#pragma omp end declare target
#endif
 
#if defined(GPUOffloadingOMP)
#pragma omp declare target
#endif
template <class QOutEnumeratorType>
void exahype2::fd::internal::clearSolution_LoopBody(
  int                                       patchIndex,
  const tarch::la::Vector<Dimensions,int>&  volumeIndex,
  int                                       unknown,
  double* __restrict__                      QOut,
  const QOutEnumeratorType&                 QOutEnumerator
) {
  QOut[ QOutEnumerator(patchIndex,volumeIndex,unknown) ] = 0.0;
}
#if defined(GPUOffloadingOMP)
#pragma omp end declare target
#endif
 
template <class QInEnumeratorType, class QOutEnumeratorType>
void exahype2::fd::internal::addAlgebraicSourceTerm_LoopBody(
  const double* __restrict__                   QIn,
  const QInEnumeratorType&                     QInEnumerator,
  exahype2::fd::Source                     AlgebraicSource,
  const tarch::la::Vector<Dimensions,double>&  patchCentre,
  const tarch::la::Vector<Dimensions,double>&  patchSize,
  int                                          patchIndex,
  const tarch::la::Vector<Dimensions,int>&     volumeIndex,
  double                                       t,
  double                                       dt,
  double* __restrict__                         QOut,
  const QOutEnumeratorType&                    QOutEnumerator
) {
  double* QInGathered              = new double[QInEnumerator._unknowns+QInEnumerator._numberOfAuxiliaryVariables];
  double* QAlgebraicSourceGathered = new double[QInEnumerator._unknowns];
 
  // gather
  for (int unknown=0; unknown<QInEnumerator._unknowns+QInEnumerator._numberOfAuxiliaryVariables; unknown++) {
    QInGathered[unknown] = QIn[ QInEnumerator(patchIndex,volumeIndex,unknown) ];
  }
 
  AlgebraicSource(
      QInGathered,
      ::exahype2::fd::getGridCellCentre( patchCentre, patchSize, QInEnumerator._numberOfDoFsPerAxisInCell, volumeIndex ),
      ::exahype2::fd::getGridCellSize( patchSize, QInEnumerator._numberOfDoFsPerAxisInCell ),
      t,
      dt,
      QAlgebraicSourceGathered
  );
 
  for (int unknown=0; unknown<QOutEnumerator._unknowns; unknown++) {
    QOut[ QOutEnumerator(patchIndex,volumeIndex,unknown) ] += dt * QAlgebraicSourceGathered[unknown];
    assertion4(
      QOut[ QOutEnumerator(patchIndex, volumeIndex,unknown) ]==QOut[ QOutEnumerator(patchIndex, volumeIndex,unknown) ],
      dt, unknown, volumeIndex,
      QOutEnumerator(patchIndex,volumeIndex,unknown)
    );
  }
 
  delete[] QInGathered;
  delete[] QAlgebraicSourceGathered;
}
 
template <typename Solver, class QInEnumeratorType, class QOutEnumeratorType>
void exahype2::fd::internal::addAlgebraicSourceTerm_LoopBody(
  const double* __restrict__                   QIn,
  const QInEnumeratorType&                     QInEnumerator,
  const tarch::la::Vector<Dimensions,double>&  patchCentre,
  const tarch::la::Vector<Dimensions,double>&  patchSize,
  int                                          patchIndex,
  const tarch::la::Vector<Dimensions,int>&     volumeIndex,
  double                                       t,
  double                                       dt,
  double* __restrict__                         QOut,
  const QOutEnumeratorType&                    QOutEnumerator
) {
  double* QInGathered              = new double[QInEnumerator._unknowns+QInEnumerator._numberOfAuxiliaryVariables];
  double* QAlgebraicSourceGathered = new double[QInEnumerator._unknowns];
 
  // gather
  for (int unknown=0; unknown<QInEnumerator._unknowns+QInEnumerator._numberOfAuxiliaryVariables; unknown++) {
    QInGathered[unknown] = QIn[ QInEnumerator(patchIndex,volumeIndex,unknown) ];
  }
 
  Solver::sourceTerm(
      QInGathered,
      ::exahype2::fd::getGridCellCentre( patchCentre, patchSize, QInEnumerator._numberOfDoFsPerAxisInCell, volumeIndex ),
      ::exahype2::fd::getGridCellSize( patchSize, QInEnumerator._numberOfDoFsPerAxisInCell ),
      t,
      dt,
      QAlgebraicSourceGathered,
      ::exahype2::Solver::Offloadable::Yes
  );
 
  for (int unknown=0; unknown<QOutEnumerator._unknowns; unknown++) {
    QOut[ QOutEnumerator(patchIndex,volumeIndex,unknown) ] += dt * QAlgebraicSourceGathered[unknown];
    assertion4(
      QOut[ QOutEnumerator(patchIndex, volumeIndex,unknown) ]==QOut[ QOutEnumerator(patchIndex, volumeIndex,unknown) ],
      dt, unknown, volumeIndex,
      QOutEnumerator(patchIndex,volumeIndex,unknown)
    );
  }
 
  delete[] QInGathered;
  delete[] QAlgebraicSourceGathered;
}
 
template <typename Solver>
void exahype2::fd::internal::addAlgebraicSourceTerm_LoopBody(
  const double* __restrict__                   QIn,
  const exahype2::enumerator::AoSLexicographicEnumerator& QInEnumerator,
  const tarch::la::Vector<Dimensions,double>&  patchCentre,
  const tarch::la::Vector<Dimensions,double>&  patchSize,
  int                                          patchIndex,
  const tarch::la::Vector<Dimensions,int>&     volumeIndex,
  double                                       t,
  double                                       dt,
  double* __restrict__                         QOut,
  const exahype2::enumerator::AoSLexicographicEnumerator& QOutEnumerator
) {
  double* QAlgebraicSourceGathered = new double[QInEnumerator._unknowns];
 
  Solver::sourceTerm(
      & (QIn[ QInEnumerator(patchIndex,volumeIndex,0) ]),
      ::exahype2::fd::getGridCellCentre( patchCentre, patchSize, QInEnumerator._numberOfDoFsPerAxisInCell, volumeIndex ),
      ::exahype2::fd::getGridCellSize( patchSize, QInEnumerator._numberOfDoFsPerAxisInCell ),
      t,
      dt,
      QAlgebraicSourceGathered,
      ::exahype2::Solver::Offloadable::Yes
  );
 
  for (int unknown=0; unknown<QOutEnumerator._unknowns; unknown++) {
    QOut[ QOutEnumerator(patchIndex,volumeIndex,unknown) ] += dt * QAlgebraicSourceGathered[unknown];
    assertion4(
      QOut[ QOutEnumerator(patchIndex, volumeIndex,unknown) ]==QOut[ QOutEnumerator(patchIndex, volumeIndex,unknown) ],
      dt, unknown, volumeIndex,
      QOutEnumerator(patchIndex,volumeIndex,unknown)
    );
  }
 
  delete[] QAlgebraicSourceGathered;
}
 
template <class QOutEnumeratorType>
double exahype2::fd::internal::reduceMaxEigenvalue_LoopBody(
    const double* __restrict__                   QOut,
    const QOutEnumeratorType&                    QOutEnumerator,
    exahype2::fd::MaxEigenvalue                  maxEigenvalue,
    const tarch::la::Vector<Dimensions,double>&  patchCentre,
    const tarch::la::Vector<Dimensions,double>&  patchSize,
    int                                          patchIndex,
    const tarch::la::Vector<Dimensions,int>&     volumeIndex,
    double                                       t,
    double                                       dt
) {
    double* QGathered = new double[QOutEnumerator._unknowns+QOutEnumerator._numberOfAuxiliaryVariables];
 
    // gather
    #if defined(SharedOMP)
    #pragma omp simd
    #endif
    for (int unknown=0; unknown<QOutEnumerator._unknowns+QOutEnumerator._numberOfAuxiliaryVariables; unknown++) {
      QGathered[unknown] = QOut[ QOutEnumerator(patchIndex,volumeIndex,unknown) ];
    }
 
    double result = 0.0;
    for (int normal=0; normal<Dimensions; normal++) {
      result = std::max(
        result,
        maxEigenvalue(
          QGathered,
          ::exahype2::fd::getGridCellCentre( patchCentre, patchSize, QOutEnumerator._numberOfDoFsPerAxisInCell, volumeIndex ),
          ::exahype2::fd::getGridCellSize( patchSize, QOutEnumerator._numberOfDoFsPerAxisInCell ),
          t,
          dt,
          normal,
          ::exahype2::Solver::Offloadable::Yes
        )
      );
    }
 
    delete[] QGathered;
 
    return result;
}
 
template <typename Solver>
double exahype2::fd::internal::reduceMaxEigenvalue_LoopBody(
    const double* __restrict__                   QOut,
    const exahype2::enumerator::AoSLexicographicEnumerator& QOutEnumerator,
    const tarch::la::Vector<Dimensions,double>&  patchCentre,
    const tarch::la::Vector<Dimensions,double>&  patchSize,
    int                                          patchIndex,
    const tarch::la::Vector<Dimensions,int>&     volumeIndex,
    double                                       t,
    double                                       dt
) {
    double result = 0.0;
    for (int normal=0; normal<Dimensions; normal++) {
      result = std::max(
        result,
        Solver::maxEigenvalue(
          &( QOut[ QOutEnumerator(patchIndex,volumeIndex,0) ] ),
          ::exahype2::fd::getGridCellCentre( patchCentre, patchSize, QOutEnumerator._numberOfDoFsPerAxisInCell, volumeIndex ),
          ::exahype2::fd::getGridCellSize( patchSize, QOutEnumerator._numberOfDoFsPerAxisInCell ),
          t,
          dt,
          normal,
          ::exahype2::Solver::Offloadable::Yes
        )
      );
    }
    return result;
}