FrictionLaws.cpph

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
 
template <class T, class Shortcuts>
static inline GPU_CALLABLE_METHOD double applications::exahype2::ShallowWater::FrictionLaws::maxEigenvalue(
  const T* const NOALIAS                       Q,
  const tarch::la::Vector<DIMENSIONS, double>& x,
  const tarch::la::Vector<DIMENSIONS, double>& h,
  const double                                 t,
  const double                                 dt,
  const int                                    normal
) {
  if (Q[Shortcuts::h] <= hThreshold) {
    return 0.0;
  }
 
  // Wave speed estimation based on the flux Jacobian
  const auto Vn{Q[Shortcuts::hu + normal] / Q[Shortcuts::h]};
  const auto c{std::sqrt(g * Q[Shortcuts::h])};
  const auto sFlux{std::fmax(std::fabs(Vn + c), std::fabs(Vn - c))};
 
  // Wave speed estimation based on the stiff source
  const auto momentum{std::sqrt(Q[Shortcuts::hu] * Q[Shortcuts::hu] + Q[Shortcuts::hv] * Q[Shortcuts::hv])};
  const auto sSource{2.0 * g * momentum * invXi / (Q[Shortcuts::h] * Q[Shortcuts::h])};
 
  return std::fmax(sFlux, h[normal] * sSource);
}
 
template <class T, class Shortcuts, int NumberOfUnknowns>
static inline GPU_CALLABLE_METHOD void applications::exahype2::ShallowWater::FrictionLaws::nonconservativeProduct(
  const T* const NOALIAS                       Q,
  const T* const NOALIAS                       deltaQ,
  const tarch::la::Vector<DIMENSIONS, double>& x,
  const tarch::la::Vector<DIMENSIONS, double>& h,
  const double                                 t,
  const double                                 dt,
  const int                                    normal,
  T* const NOALIAS                             BTimesDeltaQ
) {
  for (int n = 0; n < NumberOfUnknowns; n++) {
    BTimesDeltaQ[n] = 0.0;
  }
 
  if (Q[Shortcuts::h] <= hThreshold) {
    return;
  }
 
  // Bathymetry/topography-related effects
  BTimesDeltaQ[Shortcuts::hu + normal] = g * Q[Shortcuts::h] * (deltaQ[Shortcuts::h] + deltaQ[Shortcuts::z]);
 
  // Compute stiff turbulent drag on the interface
  const auto momentumSQR{Q[Shortcuts::hu] * Q[Shortcuts::hu] + Q[Shortcuts::hv] * Q[Shortcuts::hv]};
  const auto momentum{std::sqrt(momentumSQR)};
  if (tarch::la::greater(momentum / Q[Shortcuts::h], 0.0)) {
    const auto turbulentDrag{g * momentumSQR * invXi / (Q[Shortcuts::h] * Q[Shortcuts::h])};
    const auto directionVect{Q[Shortcuts::hu + normal] / momentum};
    BTimesDeltaQ[Shortcuts::hu + normal] += h[normal] * turbulentDrag * directionVect;
  }
}
 
template <class T, class Shortcuts, int NumberOfUnknowns, int NumberOfAuxiliaryVariables>
static inline GPU_CALLABLE_METHOD void applications::exahype2::ShallowWater::FrictionLaws::sourceTermADERDG(
  const T* const NOALIAS                       Q,
  const T* const NOALIAS                       deltaQ,
  const tarch::la::Vector<DIMENSIONS, double>& x,
  const tarch::la::Vector<DIMENSIONS, double>& h,
  const double                                 t,
  const double                                 dt,
  T* const NOALIAS                             S
) {
  for (int n = 0; n < NumberOfUnknowns; n++) {
    S[n] = 0.0;
  }
 
  if (Q[Shortcuts::h] <= hThreshold) {
    return;
  }
 
  // Compute local slope cF
  const auto dzx{deltaQ[Shortcuts::z]};
  const auto dzy{deltaQ[Shortcuts::z + NumberOfUnknowns + NumberOfAuxiliaryVariables]};
  const auto cF{1.0 / std::sqrt(1.0 + dzx * dzx + dzy * dzy)};
 
  // Apply friction only if there is some movement
  const auto momentumSQR{Q[Shortcuts::hu] * Q[Shortcuts::hu] + Q[Shortcuts::hv] * Q[Shortcuts::hv]};
  const auto momentum{std::sqrt(momentumSQR)};
  if (tarch::la::greater(momentum / Q[Shortcuts::h], 0.0)) {
    const auto coulombFriction{mu * cF * Q[Shortcuts::h]};
    const auto turbulentDrag{momentumSQR * invXi / (Q[Shortcuts::h] * Q[Shortcuts::h])};
    const auto frictionTerm{-g * (coulombFriction + turbulentDrag) / momentum};
    S[Shortcuts::hu] = Q[Shortcuts::hu] * frictionTerm;
    S[Shortcuts::hv] = Q[Shortcuts::hv] * frictionTerm;
  }
}
 
template <class T, class Shortcuts, int NumberOfUnknowns, int NumberOfAuxiliaryVariables>
static inline GPU_CALLABLE_METHOD void applications::exahype2::ShallowWater::FrictionLaws::sourceTermFV(
  const T* const NOALIAS                       Q,
  const T* const NOALIAS                       deltaQ,
  const tarch::la::Vector<DIMENSIONS, double>& x,
  const tarch::la::Vector<DIMENSIONS, double>& h,
  const double                                 t,
  const double                                 dt,
  T* const NOALIAS                             S
) {
  for (int n = 0; n < NumberOfUnknowns; n++) {
    S[n] = 0.0;
  }
 
  if (Q[Shortcuts::h] <= hThreshold) {
    return;
  }
 
  // Compute local slope cF
  const auto dzx{deltaQ[Shortcuts::z]};
  const auto dzy{deltaQ[Shortcuts::z + NumberOfUnknowns + NumberOfAuxiliaryVariables]};
  const auto cF{1.0 / std::sqrt(1.0 + dzx * dzx + dzy * dzy)};
 
  // Apply friction only if there is some movement
  const auto momentumSQR{Q[Shortcuts::hu] * Q[Shortcuts::hu] + Q[Shortcuts::hv] * Q[Shortcuts::hv]};
  const auto momentum{std::sqrt(momentumSQR)};
  if (tarch::la::greater(momentum / Q[Shortcuts::h], 0.0)) {
    const auto coulombFriction{mu * cF * Q[Shortcuts::h]};
    const auto frictionTerm{-g * coulombFriction / momentum};
    S[Shortcuts::hu] = Q[Shortcuts::hu] * frictionTerm;
    S[Shortcuts::hv] = Q[Shortcuts::hv] * frictionTerm;
  }
}