Extending SWE: 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] = 1.0;
  Q[Shortcuts::b] = (::tarch::la::norm2(x) < 0.5 ? 0.2 : 0.0);
}

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];
  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;
}

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]
) {
  constexpr double g = 9.81;
  BTimesDeltaQ[Shortcuts::h] = 0.0;
  switch (normal) {
  case 0:
    BTimesDeltaQ[Shortcuts::hu + 0] = g * Q[Shortcuts::h] * (deltaQ[Shortcuts::h] + deltaQ[Shortcuts::b]);
    BTimesDeltaQ[Shortcuts::hu + 1] = 0.0;
    break;
  case 1:
    BTimesDeltaQ[Shortcuts::hu + 0] = 0.0;
    BTimesDeltaQ[Shortcuts::hu + 1] = g * Q[Shortcuts::h] * (deltaQ[Shortcuts::h] + deltaQ[Shortcuts::b]);
    break;
  }
  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 (std::abs(::tarch::la::norm2(x) - 0.5) < 0.2 ?
    ::exahype2::RefinementCommand::Refine :
    ::exahype2::RefinementCommand::Keep);
}

The boundary conditions and the maximal eigenvalue remain identical to the solution of exercise 1 The Shallow Water Equations: Solution

Overall these should yield following end results: