Line data    Source code

       1            0 : // Distributed under the MIT License.
       2              : // See LICENSE.txt for details.
       3              : 
       4              : #pragma once
       5              : 
       6              : #include <array>
       7              : #include <cmath>
       8              : #include <cstddef>
       9              : #include <utility>
      10              : 
      11              : #include "NumericalAlgorithms/FiniteDifference/Reconstruct.hpp"
      12              : #include "Utilities/ForceInline.hpp"
      13              : #include "Utilities/Gsl.hpp"
      14              : #include "Utilities/Math.hpp"
      15              : #include "Utilities/TMPL.hpp"
      16              : 
      17              : /// \cond
      18              : template <size_t Dim>
      19              : class Direction;
      20              : template <size_t Dim>
      21              : class Index;
      22              : /// \endcond
      23              : 
      24              : namespace fd::reconstruction {
      25              : namespace detail {
      26              : struct MonotonisedCentralReconstructor {
      27              :   SPECTRE_ALWAYS_INLINE static std::array<double, 2> pointwise(
      28              :       const double* const q, const int stride) {
      29              :     using std::abs;
      30              : 
      31              :     const double a = q[stride] - q[0];
      32              :     const double b = q[0] - q[-stride];
      33              : 
      34              :     if (sgn(a) != sgn(b)) {
      35              :       return {{q[0], q[0]}};
      36              :     }
      37              :     if (3.0 * abs(a) <= abs(b)) {
      38              :       return {{q[0] - a, q[stride]}};
      39              :     } else if (3.0 * abs(b) <= abs(a)) {
      40              :       return {{q[-stride], q[0] + b}};
      41              :     } else {
      42              :       const double slope = 0.5 * (q[stride] - q[-stride]);
      43              :       return {{q[0] - 0.5 * slope, q[0] + 0.5 * slope}};
      44              :     }
      45              :   }
      46              : 
      47              :   SPECTRE_ALWAYS_INLINE static constexpr size_t stencil_width() { return 3; }
      48              : };
      49              : }  // namespace detail
      50              : 
      51              : /*!
      52              :  * \ingroup FiniteDifferenceGroup
      53              :  * \brief Performs monotonised central-difference reconstruction on the `vars`
      54              :  * in each direction.
      55              :  *
      56              :  * On a 1d mesh with spacing \f$\Delta x\f$ we denote the solution at the
      57              :  * \f$j\f$th point by \f$u_j\f$. The monotonised central-difference
      58              :  * reconstruction \cite RezzollaBook first computes:
      59              :  *
      60              :  * \f{align}
      61              :  * \sigma_j\equiv \textrm{minmod}
      62              :  *           \left(\frac{u_{j+1} - u_{j-1}}{2\Delta x},
      63              :  *                 2\frac{u_j-u_{j-1}}{\Delta x},
      64              :  *                 2\frac{u_{j+1}-u_j}{\Delta x}\right)
      65              :  * \f}
      66              :  *
      67              :  * where
      68              :  *
      69              :  * \f{align}
      70              :  *  \mathrm{minmod}(a,b,c)=
      71              :  * \left\{
      72              :  *   \begin{array}{ll}
      73              :  *    \mathrm{sgn}(a)\min(\lvert a\rvert, \lvert b\rvert, \lvert c\rvert)
      74              :  *      & \mathrm{if} \; \mathrm{sgn}(a)=\mathrm{sgn}(b)=\mathrm{sgn}(c) \\
      75              :  *    0 & \mathrm{otherwise}
      76              :  *   \end{array}\right.
      77              :  * \f}
      78              :  *
      79              :  * The reconstructed solution \f$u_j(x)\f$ in the \f$j\f$th cell is given by
      80              :  *
      81              :  * \f{align}
      82              :  * u_j(x)=u_j + \sigma_j (x-x_j)
      83              :  * \f}
      84              :  *
      85              :  * where \f$x_j\f$ is the coordinate \f$x\f$ at the center of the cell.
      86              :  */
      87              : template <size_t Dim>
      88            1 : void monotonised_central(
      89              :     const gsl::not_null<std::array<gsl::span<double>, Dim>*>
      90              :         reconstructed_upper_side_of_face_vars,
      91              :     const gsl::not_null<std::array<gsl::span<double>, Dim>*>
      92              :         reconstructed_lower_side_of_face_vars,
      93              :     const gsl::span<const double>& volume_vars,
      94              :     const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
      95              :     const Index<Dim>& volume_extents, const size_t number_of_variables) {
      96              :   detail::reconstruct<detail::MonotonisedCentralReconstructor>(
      97              :       reconstructed_upper_side_of_face_vars,
      98              :       reconstructed_lower_side_of_face_vars, volume_vars, ghost_cell_vars,
      99              :       volume_extents, number_of_variables);
     100              : }
     101              : }  // namespace fd::reconstruction