HLLCRiemannSolver.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
 
namespace applications::exahype2::ShallowWater {
  constexpr int    NEWTON_ITER     = 1;
  constexpr double NEWTON_TOL      = 1e-5;
  constexpr double CONVERGENCE_TOL = 1e-5;
 
  template <class T>
  static inline GPU_CALLABLE_METHOD double middleHeightEstimation(
    const T& hLeft,
    const T& hRight,
    const T& uLeft,
    const T& uRight
  ) {
    T hMin = std::fmin(hLeft, hRight);
    T hMax = std::fmax(hLeft, hRight);
 
    T uDiff   = uRight - uLeft;
    T hMiddle = 0;
 
    if (hLeft <= hThreshold or hRight <= hThreshold) {
      hMiddle = 0;
      return hMiddle;
    }
 
    T phiMin            = uDiff + 2.0 * (std::sqrt(g * hMin) - std::sqrt(g * hMax));
    T phiMax            = uDiff + (hMax - hMin) * (std::sqrt(0.5 * g * (hMax + hMin) / (hMax * hMin)));
    T hZero             = 0;
    T gLeft             = 0;
    T gRight            = 0;
    T phiZero           = 0;
    T phiZeroDerivative = 0;
    T slope             = 0;
 
    if (phiMax <= 0.0) { // Shock-shock Riemann structure
      hZero = hMax;
      for (int i = 0; i < NEWTON_ITER; ++i) {
        gLeft   = std::sqrt(0.5 * g * (1 / hZero + 1 / hLeft));
        gRight  = std::sqrt(0.5 * g * (1 / hZero + 1 / hRight));
        phiZero = uDiff + (hZero - hLeft) * gLeft + (hZero - hRight) * gRight;
        if (std::fabs(phiZero) < NEWTON_TOL) {
          break;
        }
        phiZeroDerivative = gLeft - g * (hZero - hLeft) / (4.0 * (hZero * hZero) * gLeft) + gRight
                            - g * (hZero - hRight) / (4.0 * (hZero * hZero) * gRight);
        slope = 2.0 * std::sqrt(hZero) * phiZeroDerivative;
        hZero = std::pow((std::sqrt(hZero) - phiZero / slope), 2);
      }
      hMiddle = hZero;
      return hMiddle;
    } else if (phiMin > 0.0) { // Rarefaction-rarefaction
      hMiddle = std::max(0.0, uLeft - uRight + 2.0 * (std::sqrt(g * hLeft) + std::sqrt(g * hRight)));
      hMiddle = 1.0 / (16.0 * g) * hMiddle * hMiddle;
      return hMiddle;
    } else { // One rarefaction-one shock
      hZero = hMin;
      for (int i = 0; i < NEWTON_ITER; ++i) {
        phiZero = uDiff + 2.0 * (std::sqrt(g * hZero) - std::sqrt(g * hMax))
                  + (hZero - hMin) * std::sqrt(0.5 * g * (1 / hZero + 1 / hMin));
        slope = (phiMax - phiZero) / (hMax - hMin);
        hZero = hZero - phiZero / slope;
      }
      hMiddle = hZero;
      return hMiddle;
    }
  }
} // namespace applications::exahype2::ShallowWater
 
template <class T, class Shortcuts>
static inline GPU_CALLABLE_METHOD double applications::exahype2::ShallowWater::hllcRiemannSolver(
  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
) {
  // Correction of particle velocities according to Kurganov (modified version).
  // The vectors contain the horizontal velocities u and v of the cells.
  T horizontalVelocitiesLeft[2] = {
    QL[Shortcuts::hu] * QL[Shortcuts::h] * std::sqrt(2)
      / std::sqrt(std::pow(QL[Shortcuts::h], 4) + std::pow(std::fmax(QL[Shortcuts::h], hThreshold), 4)),
    QL[Shortcuts::hv] * QL[Shortcuts::h] * std::sqrt(2)
      / std::sqrt(std::pow(QL[Shortcuts::h], 4) + std::pow(std::fmax(QL[Shortcuts::h], hThreshold), 4))
  };
 
  T horizontalVelocitiesRight[2] = {
    QR[Shortcuts::hu] * QR[Shortcuts::h] * std::sqrt(2)
      / std::sqrt(std::pow(QR[Shortcuts::h], 4) + std::pow(std::fmax(QR[Shortcuts::h], hThreshold), 4)),
    QR[Shortcuts::hv] * QR[Shortcuts::h] * std::sqrt(2)
      / std::sqrt(std::pow(QR[Shortcuts::h], 4) + std::pow(std::fmax(QR[Shortcuts::h], hThreshold), 4))
  };
 
  // Hydrostatic reconstruction of the heights
  const T zMax   = std::fmax(QL[Shortcuts::z], QR[Shortcuts::z]);
  const T hLeft  = std::fmax(0.0, QL[Shortcuts::h] + QL[Shortcuts::z] - zMax);
  const T hRight = std::fmax(0.0, QR[Shortcuts::h] + QR[Shortcuts::z] - zMax);
 
  // Check for dry land after the hydrostatic reconstruction
  // as it might set both heights to zero or close to zero
  // even if they weren't before the reconstruction.
  const bool leftWet  = tarch::la::greater(hLeft, hThreshold);
  const bool rightWet = tarch::la::greater(hRight, hThreshold);
  const bool leftDry  = not leftWet;
  const bool rightDry = not rightWet;
 
  for (int i = 0; i < 3; ++i) {
    FL[i] = 0.0;
    FR[i] = 0.0;
  }
 
  // Dry land, skip the Riemann problem
  if (leftDry and rightDry) {
    return 0.0;
  }
 
  // Also reconstruct the momenta based on the reconstructed height.
  // Compute the conserved momentum variables according to the Kurganov modified speeds.
  const T huLeft  = hLeft * horizontalVelocitiesLeft[0];
  const T hvLeft  = hLeft * horizontalVelocitiesLeft[1];
  const T huRight = hRight * horizontalVelocitiesRight[0];
  const T hvRight = hRight * horizontalVelocitiesRight[1];
 
  // Wave speed estimation process
  // 1) calculate the first h* estimate
  // celerity a = squareRoot(g * h)
  T celerityLeft  = std::sqrt(g * hLeft);
  T celerityRight = std::sqrt(g * hRight);
  T hStarEstimate = 0.0;
 
  // The following is a case distinction which determines which wave speed estimate is used:
  // If the option UseTwoRarefactionsEstimate is enabled, the middle height estimation is implemented as described in
  // Toro2019 and Toro2024.
  // If not, the middle height estimation process is implemented as in the Clawpack Riemann solver, which should be
  // more exact (this version is still experimental).
  constexpr bool UseTwoRarefactionsEstimate = true;
  if constexpr (UseTwoRarefactionsEstimate) {
    // Estimate of the celerity based on the assumption of two outer rarefaction waves
    T celerityEstimate = 0.5 * (celerityLeft + celerityRight)
                         - 0.25 * (horizontalVelocitiesRight[normal] - horizontalVelocitiesLeft[normal]);
    // Estimate for the middle state height
    hStarEstimate = std::pow(celerityEstimate, 2) / g;
    // This is an experimental idea which could be used for better capturing developing vacuum states
    // if ((2 * (celerityLeft + celerityRight)) <= (horizontalVelocitiesRight[normal] -
    // horizontalVelocitiesLeft[normal])) { // Vacuum condition, check for safety
    //  hStarEstimate = 0.0;
    //}
  } else {
    hStarEstimate = middleHeightEstimation(
      hLeft,
      hRight,
      horizontalVelocitiesLeft[normal],
      horizontalVelocitiesRight[normal]
    );
  }
 
  // Additionally, as described in the book Computational Algorithms for Shallow Water Equations, the
  // wave speeds estimates will be modified depending on whether there is a wet/dry front in order to better capture
  // the speed of that wet/dry front, described in the book section 11.10.3 Dry-Bed Approximate Riemann Solvers.
 
  // Outer wave speeds
  T sLeft  = 0;
  T sRight = 0;
  T qLeft  = 0;
  T qRight = 0;
 
  // 2) determine qL and qR, case distinction whether a shock wave or a rarefaction wave is expected && 3) compute wave
  // speed estimates sLeft, sRight and sStar.
  if (hLeft > hThreshold) {
    if (hRight > hThreshold) {
      qLeft = (hStarEstimate > hLeft)
                ? std::sqrt(0.5 * ((hStarEstimate + hLeft) * hStarEstimate) / (std::pow(hLeft, 2)))
                : 1.0;
      sLeft = horizontalVelocitiesLeft[normal] - qLeft * celerityLeft;
    } else {
      // If the right cell is dry, the left cell should always contain a rarefaction and not a shock.
      sLeft = horizontalVelocitiesLeft[normal] - 1.0 * celerityLeft;
    }
  } else {
    // Left cell is dry, its respective speed is set to the wet dry front characteristic of the single rarefaction fan.
    sLeft = horizontalVelocitiesRight[normal] - 2 * celerityRight;
  }
 
  if (hRight > hThreshold) {
    if (hLeft > hThreshold) {
      qRight = (hStarEstimate > hRight)
                 ? std::sqrt(0.5 * ((hStarEstimate + hRight) * hStarEstimate) / (std::pow(hRight, 2)))
                 : 1.0;
      sRight = horizontalVelocitiesRight[normal] + qRight * celerityRight;
    } else {
      // If the left cell is dry, the right cell should always contain a rarefaction and not a shock.
      sRight = horizontalVelocitiesRight[normal] + 1.0 * celerityRight;
    }
  } else {
    // Right cell is dry, its respective speed is set to the wet dry front characteristic of the single rarefaction fan.
    sRight = horizontalVelocitiesLeft[normal] + 2 * celerityLeft;
  }
 
  // Middle speed of the shear wave
  T sStar = 0;
 
  // If one takes the formula for sStar and the above speed estimates, sStar coincides with either the left or the
  // right wave speed respectively in the case of a dry state Riemann problem.
  if (hLeft <= hThreshold) {
    sStar = sLeft;
  } else if (hRight <= hThreshold) {
    sStar = sRight;
  } else {
    sStar
      = (sLeft * hRight * (horizontalVelocitiesRight[normal] - sRight)
         - sRight * hLeft * (horizontalVelocitiesLeft[normal] - sLeft))
        / (hRight * (horizontalVelocitiesRight[normal] - sRight) - hLeft * (horizontalVelocitiesLeft[normal] - sLeft)); // TODO: Possible instability
  }
 
  // The middle wave speed should not be higher than the outer wave speeds, so the left and right wave speeds.
  // Thus, it should be limited.
  sStar = std::max(sLeft, sStar);
  sStar = std::min(sRight, sStar);
 
  // 4) Compute qStarLeft and qStarRight (only differ in the transverse velocity component).
  // Compute hStar based on wave speeds estimates sLeft and sRight.
  T qStarLeft[3];
  T qStarRight[3];
  T hStarLeft  = 0;
  T hStarRight = 0;
 
  // Compute the conserved variable states Q in the star region
  if (hLeft > hThreshold) {
    hStarLeft = hLeft * ((sLeft - horizontalVelocitiesLeft[normal]) / (sLeft - sStar));
    // TODO: This commented-out code should hold mathematically if the left cell is wet and the right cell is dry but
    // still experimental
    if (hRight <= hThreshold) {
      hStarLeft = (1 / 3) * hLeft;
    }
    qStarLeft[0] = hStarLeft;
    qStarLeft[1] = (normal == 0) ? (hStarLeft * sStar) : (hStarLeft * horizontalVelocitiesLeft[0]);
    qStarLeft[2] = (normal == 0) ? (hStarLeft * horizontalVelocitiesLeft[1]) : (hStarLeft * sStar);
  } else {
    qStarLeft[0] = 0;
    qStarLeft[1] = 0;
    qStarLeft[2] = 0;
  }
 
  if (hRight > hThreshold) {
    if (hLeft <= hThreshold) {
      hStarRight = (1 / 3) * hRight;
    }
    hStarRight    = hRight * ((sRight - horizontalVelocitiesRight[normal]) / (sRight - sStar));
    qStarRight[0] = hStarRight;
    qStarRight[1] = (normal == 0) ? (hStarRight * sStar) : (hStarRight * horizontalVelocitiesRight[0]);
    qStarRight[2] = (normal == 0) ? (hStarRight * horizontalVelocitiesRight[1]) : (hStarRight * sStar);
  } else {
    qStarRight[0] = 0;
    qStarRight[1] = 0;
    qStarRight[2] = 0;
  }
 
  // Compute the fluxes in the left and right cell
  double fluxLeft[3] = {
    normal == 0 ? (hLeft * horizontalVelocitiesLeft[0]) : (hLeft * horizontalVelocitiesLeft[1]),
    normal == 0
      ? (hLeft * horizontalVelocitiesLeft[0] * horizontalVelocitiesLeft[0] + 0.5 * g * std::pow(hLeft, 2))
      : (hLeft * horizontalVelocitiesLeft[1] * horizontalVelocitiesLeft[0]),
    normal == 0
      ? (hLeft * horizontalVelocitiesLeft[0] * horizontalVelocitiesLeft[1])
      : (hLeft * horizontalVelocitiesLeft[1] * horizontalVelocitiesLeft[1] + 0.5 * g * std::pow(hLeft, 2))
  };
 
  double fluxRight[3] = {
    normal == 0 ? (hRight * horizontalVelocitiesRight[0]) : (hRight * horizontalVelocitiesRight[1]),
    normal == 0
      ? (hRight * horizontalVelocitiesRight[0] * horizontalVelocitiesRight[0] + 0.5 * g * std::pow(hRight, 2))
      : (hRight * horizontalVelocitiesRight[1] * horizontalVelocitiesRight[0]),
    normal == 0
      ? (hRight * horizontalVelocitiesRight[0] * horizontalVelocitiesRight[1])
      : (hRight * horizontalVelocitiesRight[1] * horizontalVelocitiesRight[1] + 0.5 * g * std::pow(hRight, 2))
  };
 
  // 5) Now the HLLC flux can be calculated, case distinction based on sLeft and sRight.
  double hllcFlux[3];
  if (hLeft > hThreshold && hRight > hThreshold) {
    if (sLeft >= 0) {
      hllcFlux[0] = fluxLeft[0];
      hllcFlux[1] = fluxLeft[1];
      hllcFlux[2] = fluxLeft[2];
    } else if (sStar >= 0) {
      hllcFlux[0] = fluxLeft[0] + sLeft * (qStarLeft[0] - hLeft);
      hllcFlux[1] = fluxLeft[1] + sLeft * (qStarLeft[1] - (hLeft * horizontalVelocitiesLeft[0]));
      hllcFlux[2] = fluxLeft[2] + sLeft * (qStarLeft[2] - (hLeft * horizontalVelocitiesLeft[1]));
    } else if (sRight >= 0) {
      hllcFlux[0] = fluxRight[0] + sRight * (qStarRight[0] - hRight);
      hllcFlux[1] = fluxRight[1] + sRight * (qStarRight[1] - (hRight * horizontalVelocitiesRight[0]));
      hllcFlux[2] = fluxRight[2] + sRight * (qStarRight[2] - (hRight * horizontalVelocitiesRight[1]));
    } else { // if sRight <= 0
      hllcFlux[0] = fluxRight[0];
      hllcFlux[1] = fluxRight[1];
      hllcFlux[2] = fluxRight[2];
    }
  } else if (hLeft <= hThreshold) { // Q_{Left}^* is not defined
    if (sLeft >= 0) {               // sLeft == sStar
      hllcFlux[0] = fluxLeft[0];
      hllcFlux[1] = fluxLeft[1];
      hllcFlux[2] = fluxLeft[2];
    } else if (sRight >= 0) {
      hllcFlux[0] = fluxRight[0] + sRight * (qStarRight[0] - hRight);
      hllcFlux[1] = fluxRight[1] + sRight * (qStarRight[1] - (hRight * horizontalVelocitiesRight[0]));
      hllcFlux[2] = fluxRight[2] + sRight * (qStarRight[2] - (hRight * horizontalVelocitiesRight[1]));
    } else {
      hllcFlux[0] = fluxRight[0];
      hllcFlux[1] = fluxRight[1];
      hllcFlux[2] = fluxRight[2];
    }
  } else if (hRight <= hThreshold) { // Q_{Right}^* is not defined
    if (sLeft >= 0) {
      hllcFlux[0] = fluxLeft[0];
      hllcFlux[1] = fluxLeft[1];
      hllcFlux[2] = fluxLeft[2];
    } else if (sRight >= 0) { // sRight == sStar
      hllcFlux[0] = fluxLeft[0] + sLeft * (qStarLeft[0] - hLeft);
      hllcFlux[1] = fluxLeft[1] + sLeft * (qStarLeft[1] - (hLeft * horizontalVelocitiesLeft[0]));
      hllcFlux[2] = fluxLeft[2] + sLeft * (qStarLeft[2] - (hLeft * horizontalVelocitiesLeft[1]));
    } else {
      hllcFlux[0] = fluxRight[0];
      hllcFlux[1] = fluxRight[1];
      hllcFlux[2] = fluxRight[2];
    }
  }
 
  // Now include the source term already in the fluctuation computation.
  // Source term is based on the reconstructed water heights (for well-balancedness).
  // Always the square of the original value minus the square of the reconstructed value.
  double sourceLeft[3] = {
    0,
    normal == 0 ? g * 0.5 * ((QL[Shortcuts::h] * QL[Shortcuts::h]) - (hLeft * hLeft)) : 0.0,
    normal == 0 ? 0.0 : g * 0.5 * ((QL[Shortcuts::h] * QL[Shortcuts::h]) - (hLeft * hLeft))
  };
 
  double sourceRight[3] = {
    0,
    normal == 0 ? g * 0.5 * ((QR[Shortcuts::h] * QR[Shortcuts::h]) - (hRight * hRight)) : 0.0,
    normal == 0 ? 0.0 : g * 0.5 * ((QR[Shortcuts::h] * QR[Shortcuts::h]) - (hRight * hRight))
  };
 
  // Compute left- and right-going fluctuations based on the HLLC-flux and the source term
  for (int i = 0; i < 3; ++i) {
    FL[i] = hllcFlux[i] + sourceLeft[i];
    FR[i] = hllcFlux[i] + sourceRight[i];
  }
 
  return (std::max(std::abs(sLeft), std::abs(sRight)));
}