The Shallow Water Equations: Solution

Initial conditions:

void exahype2::training::swe::SWESolver::initialCondition(
  [[maybe_unused]] double* const NOALIAS Q, // Q[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const bool                                   gridIsConstructed
) {
  Q[Shortcuts::hu + 0] = 0.0; // v_x
  Q[Shortcuts::hu + 1] = 0.0; // v_y
  Q[Shortcuts::h] = x[0] < 0.0 ? 1.0 : 2.0; // h
  Q[Shortcuts::b] = 0.0;
}

Boundary conditions:

void exahype2::training::swe::SWESolver::boundaryConditions(
  [[maybe_unused]] const double* const NOALIAS   Qinside, // Qinside[4+0]
  [[maybe_unused]] double* const NOALIAS         Qoutside, // Qoutside[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const double                                 t,
  [[maybe_unused]] const int                                    normal
) {
  Qoutside[Shortcuts::h]      = 1.0; // h
  Qoutside[Shortcuts::hu + 0] = 0.0; // v_x
  Qoutside[Shortcuts::hu + 1] = 0.0; // v_y
  Qoutside[Shortcuts::b]      = 0.0; // b
}

Max eigenvalue:

double exahype2::training::swe::SWESolver::maxEigenvalue(
  [[maybe_unused]] const double* const NOALIAS   Q, // Q[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const double                                 t,
  [[maybe_unused]] const double                                 dt,
  [[maybe_unused]] const int                                    normal
) {
  constexpr double g = 9.81;
  const double u = Q[Shortcuts::hu + normal] / Q[Shortcuts::h];
  const double c = std::sqrt(g * Q[Shortcuts::h]);
  return std::max(std::abs(u + c), std::abs(u - c));
}

Flux:

void exahype2::training::swe::SWESolver::flux(
  [[maybe_unused]] const double* const NOALIAS   Q, // Q[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const double                                 t,
  [[maybe_unused]] const double                                 dt,
  [[maybe_unused]] const int                                    normal,
  [[maybe_unused]] double* const NOALIAS F // F[4]
) {
  double ih = 1.0 / Q[Shortcuts::h];
  constexpr double g = 9.81;
  F[Shortcuts::h]      = Q[Shortcuts::hu + normal];
  F[Shortcuts::hu + 0] = Q[Shortcuts::hu + normal] * Q[Shortcuts::hu + 0] * ih;
  F[Shortcuts::hu + 1] = Q[Shortcuts::hu + normal] * Q[Shortcuts::hu + 1] * ih;
  F[Shortcuts::b]      = 0.0;
  F[Shortcuts::hu + normal] += 0.5 * g * Q[Shortcuts::h] * Q[Shortcuts::h];
}

Nonconservative product:

void exahype2::training::swe::SWESolver::nonconservativeProduct(
  [[maybe_unused]] const double* const NOALIAS   Q, // Q[4+0]
  [[maybe_unused]] const double* const NOALIAS   deltaQ, // deltaQ[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const double                                 t,
  [[maybe_unused]] const double                                 dt,
  [[maybe_unused]] const int                                    normal,
  [[maybe_unused]] double* const NOALIAS BTimesDeltaQ // BTimesDeltaQ[4]
) {
  BTimesDeltaQ[Shortcuts::h]      = 0.0;
  BTimesDeltaQ[Shortcuts::hu + 0] = 0.0;
  BTimesDeltaQ[Shortcuts::hu + 1] = 0.0;
  BTimesDeltaQ[Shortcuts::b]      = 0.0;
}

Refinement criterion:

::exahype2::RefinementCommand exahype2::training::swe::SWESolver::refinementCriterion(
  [[maybe_unused]] const double* const NOALIAS   Q, // Q[4+0]
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& x,
  [[maybe_unused]] const tarch::la::Vector<DIMENSIONS, double>& h,
  [[maybe_unused]] const double                                 t
) {
  return ::exahype2::RefinementCommand::Keep;
}

Overall these should yield following end results: