RusanovRiemannSolver.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, int NumberOfUnknowns, int NumberOfAuxiliaryVariables>
static inline GPU_CALLABLE_METHOD double applications::exahype2::landslide::rusanovRiemannSolverFV(
  const T* const NOALIAS                       QL,
  const T* const NOALIAS                       QR,
  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                             FL,
  T* const NOALIAS                             FR
) {
  ShallowWater::flux<T, Shortcuts, NumberOfUnknowns>(QL, x, h, t, dt, normal, FL);
  ShallowWater::flux<T, Shortcuts, NumberOfUnknowns>(QR, x, h, t, dt, normal, FR);
 
  landslide::flux<T, Shortcuts>(QL, x, h, t, dt, normal, FL);
  landslide::flux<T, Shortcuts>(QR, x, h, t, dt, normal, FR);
 
  const auto smax{std::fmax(
    std::fmax(
      ShallowWater::maxEigenvalue<T, Shortcuts>(QL, x, h, t, dt, normal),
      ShallowWater::maxEigenvalue<T, Shortcuts>(QR, x, h, t, dt, normal)
    ),
    std::fmax(
      landslide::maxEigenvalue<T, Shortcuts>(QL, x, h, t, dt, normal),
      landslide::maxEigenvalue<T, Shortcuts>(QR, x, h, t, dt, normal)
    )
  )};
  const auto N{NumberOfAuxiliaryVariables + NumberOfUnknowns};
 
  double Qavg[N]   = {};
  double DeltaQ[N] = {};
  for (int i = 0; i < NumberOfUnknowns; ++i) {
    Qavg[i]   = 0.5 * (QL[i] + QR[i]);
    DeltaQ[i] = QR[i] - QL[i];
    FL[i]     = 0.5 * (FL[i] + FR[i] - smax * DeltaQ[i]);
    FR[i]     = FL[i];
  }
 
  for (int i = NumberOfUnknowns; i < N; ++i) {
    Qavg[i]   = 0.5 * (QL[i] + QR[i]);
    DeltaQ[i] = QR[i] - QL[i];
  }
 
  double ncp[NumberOfUnknowns] = {};
  ShallowWater::nonconservativeProduct<T, Shortcuts, NumberOfUnknowns>(Qavg, DeltaQ, x, h, t, dt, normal, ncp);
  landslide::nonconservativeProduct<T, Shortcuts>(Qavg, DeltaQ, x, h, t, dt, normal, ncp);
  for (int i = 0; i < NumberOfUnknowns; ++i) {
    FL[i] += 0.5 * ncp[i];
    FR[i] -= 0.5 * ncp[i];
  }
 
  auto sL{0.0};
  if (QL[Shortcuts::lsh] > hThreshold) {
    const auto hL{QL[Shortcuts::lsh] > hThreshold ? QL[Shortcuts::lsh] : hThreshold};
    const auto momentumL{
      std::sqrt(QL[Shortcuts::lshu] * QL[Shortcuts::lshu] + QL[Shortcuts::lshv] * QL[Shortcuts::lshv])
    };
    sL = 2.0 * g * momentumL * invXi / (hL * hL);
  }
 
  auto sR{0.0};
  if (QL[Shortcuts::lsh] > hThreshold) {
    const auto hR{QR[Shortcuts::lsh] > hThreshold ? QR[Shortcuts::lsh] : hThreshold};
    const auto momentumR{
      std::sqrt(QR[Shortcuts::lshu] * QR[Shortcuts::lshu] + QR[Shortcuts::lshv] * QR[Shortcuts::lshv])
    };
    sR = 2.0 * g * momentumR * invXi / (hR * hR);
  }
 
  return std::fmax(smax, h[normal] * std::fmax(sL, sR));
}