FWaveRiemannSolver.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::fWaveRiemannSolver(
  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
) {
  // Set speeds to zero (will be determined later)
  T uLeft  = 0.0;
  T uRight = 0.0;
 
  // Not const here as they may get reassigned due to boundary conditions
  T hLeft   = QL[Shortcuts::h];
  T hRight  = QR[Shortcuts::h];
  T huLeft  = QL[Shortcuts::hu + normal];
  T huRight = QR[Shortcuts::hu + normal];
  T zLeft   = QL[Shortcuts::z];
  T zRight  = QR[Shortcuts::z];
 
  WetDryState wetDryState = WetDryState::WetWet; // Assume everything is wet
 
  // Determine the wet/dry state
  if (hLeft <= hThreshold && hRight <= hThreshold) { // Both cells are dry
    wetDryState = WetDryState::DryDry;
  } else if (hLeft <= hThreshold) { // Left cell dry, right cell wet
    uRight = huRight / hRight;
    // Set wall boundary conditions (inundation is not supported!)
    hLeft       = hRight;
    zLeft       = zRight;
    huLeft      = -huRight;
    uLeft       = -uRight;
    wetDryState = WetDryState::DryWetWall;
  } else if (hRight <= hThreshold) { // Left cell wet, right cell dry
    uLeft       = huLeft / hLeft;
    hRight      = hLeft;
    zRight      = zLeft;
    huRight     = -huLeft;
    uRight      = -uLeft;
    wetDryState = WetDryState::WetDryWall;
  } else {
    uLeft       = huLeft / hLeft;
    uRight      = huRight / hRight;
    wetDryState = WetDryState::WetWet;
  }
 
  // Zero updates in the beginning
  for (int i = 0; i < 3; i++) {
    FL[i] = 0.0;
    FR[i] = 0.0;
  }
 
  // Return in the case of dry cells
  if (wetDryState == WetDryState::DryDry) {
    return 0.0;
  }
 
  // Compute eigenvalues of the Jacobian matrices in states Q_{i-1} and Q_{i}
  const T hRoe = 0.5 * (hRight + hLeft);
  const T uRoe = (uLeft * std::sqrt(hLeft) + uRight * std::sqrt(hRight)) / (std::sqrt(hLeft) + std::sqrt(hRight));
  const T roeSpeeds[2]            = {uRoe - std::sqrt(g * hRoe), uRoe + std::sqrt(g * hRoe)};
  const T characteristicSpeeds[2] = {uLeft - std::sqrt(g * hLeft), uRight + std::sqrt(g * hRight)};
  const T einfeldtSpeeds[2]       = {
    std::fmin(characteristicSpeeds[0], roeSpeeds[0]),
    std::fmax(characteristicSpeeds[1], roeSpeeds[1])
  };
 
  // Set wave speeds (lambda Roe)
  const T waveSpeeds[2]{einfeldtSpeeds[0], einfeldtSpeeds[1]};
  const T deltaLambda = waveSpeeds[1] - waveSpeeds[0];
 
#if !defined(WITH_GPU)
  // Wave speed of the 1st wave family should be less than the speed of the 2nd wave family.
  assertion(waveSpeeds[0] < waveSpeeds[1]);
  assertion(std::fabs(deltaLambda) > 0); // Make sure we do not divide by 0
#endif
 
  const T inverseDeltaLambda = 1.0 / deltaLambda;
 
  // Compute the inverse matrix R^{-1} to obtain Eigencoefficients (alpha)
  //
  // 1                   1             ^-1
  //
  // \lambda_{1}         \lambda_{2}
  //
  // Inverse: https://mathworld.wolfram.com/MatrixInverse.html
  //
  //       1          \lambda_{2}       -1
  //  ____________
  //   deltaLambda   -\lambda_{1}        1
  const T Rinv[2][2]{
    {inverseDeltaLambda * waveSpeeds[1], -inverseDeltaLambda},
    {inverseDeltaLambda * -waveSpeeds[0], inverseDeltaLambda}
  };
 
  // Right hand side
  const T fluxLeft[2]{huLeft, huLeft * uLeft + 0.5 * g * hLeft * hLeft};
  const T fluxRight[2]{huRight, huRight * uRight + 0.5 * g * hRight * hRight};
 
  // Add source term (bathymetry) into the decomposition of the jump in the flux function
  const T deltaxPsi = -g * (zRight - zLeft) * hRoe;
  // Calculate flux difference f(Q_i) - f(Q_{i-1})
  const T deltaFlux[2]{fluxRight[0] - fluxLeft[0], fluxRight[1] - fluxLeft[1] - deltaxPsi};
 
  // Solve linear equations
  const T alpha[2]{
    Rinv[0][0] * deltaFlux[0] + Rinv[0][1] * deltaFlux[1],
    Rinv[1][0] * deltaFlux[0] + Rinv[1][1] * deltaFlux[1]
  };
 
  // Z
  const T fWaves[2][2]{{alpha[0], alpha[0] * waveSpeeds[0]}, {alpha[1], alpha[1] * waveSpeeds[1]}};
 
  // Compute the net-updates
  T hUpdateLeft   = 0.0;
  T hUpdateRight  = 0.0;
  T huUpdateLeft  = 0.0;
  T huUpdateRight = 0.0;
 
  constexpr T ZERO_TOL = std::numeric_limits<T>::epsilon();
  // 1st wave family
  if (waveSpeeds[0] < ZERO_TOL) { // Left going
    hUpdateLeft += fWaves[0][0];
    huUpdateLeft += fWaves[0][1];
  } else if (waveSpeeds[0] > ZERO_TOL) { // Right going
    hUpdateRight -= fWaves[0][0];
    huUpdateRight -= fWaves[0][1];
  } else { // Split waves
    hUpdateLeft += 0.5 * fWaves[0][0];
    huUpdateLeft += 0.5 * fWaves[0][1];
    hUpdateRight -= 0.5 * fWaves[0][0];
    huUpdateRight -= 0.5 * fWaves[0][1];
  }
 
  // 2nd wave family
  if (waveSpeeds[1] < ZERO_TOL) { // Left going
    hUpdateLeft += fWaves[1][0];
    huUpdateLeft += fWaves[1][1];
  } else if (waveSpeeds[1] > ZERO_TOL) { // Right going
    hUpdateRight -= fWaves[1][0];
    huUpdateRight -= fWaves[1][1];
  } else { // Split waves
    hUpdateLeft += 0.5 * fWaves[1][0];
    huUpdateLeft += 0.5 * fWaves[1][1];
    hUpdateRight -= 0.5 * fWaves[1][0];
    huUpdateRight -= 0.5 * fWaves[1][1];
  }
 
  // Zero ghost updates (wall boundary)
  if (wetDryState == WetDryState::WetDryWall) {
    hUpdateRight  = 0;
    huUpdateRight = 0;
  } else if (wetDryState == WetDryState::DryWetWall) {
    hUpdateLeft  = 0;
    huUpdateLeft = 0;
  }
 
  // Compute maximum wave speed (-> CFL-condition)
  const T maxWaveSpeed = std::fmax(std::fabs(waveSpeeds[0]), std::fabs(waveSpeeds[1]));
 
#if !defined(WITH_GPU)
  assertion(!std::isnan(hUpdateLeft));
  assertion(!std::isnan(hUpdateRight));
  assertion(!std::isnan(huUpdateLeft));
  assertion(!std::isnan(huUpdateRight));
  assertion(!std::isnan(maxWaveSpeed));
#endif
 
  // Set net-updates
  FL[Shortcuts::h] = hUpdateLeft;
  FR[Shortcuts::h] = hUpdateRight;
  FL[normal + 1]   = huUpdateLeft;
  FR[normal + 1]   = huUpdateRight;
 
  return maxWaveSpeed;
}