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 <cstddef>
       8              : 
       9              : #include "NumericalAlgorithms/FiniteDifference/Reconstruct.hpp"
      10              : #include "Utilities/ForceInline.hpp"
      11              : #include "Utilities/Gsl.hpp"
      12              : 
      13              : /// \cond
      14              : template <size_t Dim, typename T>
      15              : class DirectionMap;
      16              : template <size_t Dim>
      17              : class Index;
      18              : /// \endcond
      19              : 
      20              : namespace fd::reconstruction {
      21              : namespace detail {
      22              : template <size_t Degree>
      23              : struct UnlimitedReconstructor {
      24              :   static_assert(Degree == 2 or Degree == 4 or Degree == 6 or Degree == 8);
      25              :   SPECTRE_ALWAYS_INLINE static std::array<double, 2> pointwise(
      26              :       const double* const q, const int stride) {
      27              :     if constexpr (Degree == 2) {
      28              :       // quadratic polynomial
      29              :       return {{0.375 * q[-stride] + 0.75 * q[0] - 0.125 * q[stride],
      30              :                -0.125 * q[-stride] + 0.75 * q[0] + 0.375 * q[stride]}};
      31              :     } else if constexpr (Degree == 4) {
      32              :       // quintic polynomial
      33              :       return {{-0.0390625 * q[-2 * stride] + 0.46875 * q[-stride] +
      34              :                0.703125 * q[0] - 0.15625 * q[stride] +
      35              :                0.0234375 * q[2 * stride],
      36              :                0.0234375 * q[-2 * stride] - 0.15625 * q[-stride] +
      37              :                0.703125 * q[0] + 0.46875 * q[stride] -
      38              :                0.0390625 * q[2 * stride]}};
      39              :     } else if constexpr (Degree == 6) {
      40              :       return {{-0.1708984375 * q[stride] + 0.041015625 * q[2 * stride] -
      41              :                    0.0048828125 * q[3 * stride] + 0.5126953125 * q[-stride] -
      42              :                    0.068359375 * q[-2 * stride] +
      43              :                    0.0068359375 * q[-3 * stride] + 0.68359375 * q[0],
      44              :                0.5126953125 * q[stride] - 0.068359375 * q[2 * stride] +
      45              :                    0.0068359375 * q[3 * stride] - 0.1708984375 * q[-stride] +
      46              :                    0.041015625 * q[-2 * stride] -
      47              :                    0.0048828125 * q[-3 * stride] + 0.68359375 * q[0]}};
      48              :     } else if constexpr (Degree == 8) {
      49              :       return {
      50              :           {-0.179443359375 * q[stride] + 0.0538330078125 * q[2 * stride] -
      51              :                0.010986328125 * q[3 * stride] +
      52              :                0.001068115234375 * q[4 * stride] + 0.538330078125 * q[-stride] -
      53              :                0.0897216796875 * q[-2 * stride] +
      54              :                0.015380859375 * q[-3 * stride] -
      55              :                0.001373291015625 * q[-4 * stride] + 0.67291259765625 * q[0],
      56              :            0.538330078125 * q[stride] - 0.0897216796875 * q[2 * stride] +
      57              :                0.015380859375 * q[3 * stride] -
      58              :                0.001373291015625 * q[4 * stride] - 0.179443359375 * q[-stride] +
      59              :                0.0538330078125 * q[-2 * stride] -
      60              :                0.010986328125 * q[-3 * stride] +
      61              :                0.001068115234375 * q[-4 * stride] + 0.67291259765625 * q[0]}};
      62              :     }
      63              :   }
      64              : 
      65              :   SPECTRE_ALWAYS_INLINE static constexpr size_t stencil_width() {
      66              :     return Degree + 1;
      67              :   }
      68              : };
      69              : }  // namespace detail
      70              : 
      71              : /*!
      72              :  * \ingroup FiniteDifferenceGroup
      73              :  * \brief Performs unlimited reconstruction on the `vars` in each direction.
      74              :  *
      75              :  * On a 1d mesh with spacing \f$\Delta x\f$ we denote the solution at the
      76              :  * \f$j\f$th point by \f$u_j\f$. The degree 2 reconstruction is
      77              :  *
      78              :  * \f{align}{
      79              :  * u_{j+1/2} &= -\frac{1}{8} u_{j-1} + \frac{3}{4} u_j + \frac{3}{8} u_{j+1}, \\
      80              :  * u_{j-1/2} &= \frac{3}{8} u_{j-1} + \frac{3}{4} u_j - \frac{1}{8} u_{j+1}.
      81              :  * \f}
      82              :  *
      83              :  * The degree 4 reconstruction is
      84              :  *
      85              :  * \f{align}{
      86              :  *  u_{j-1/2} &= - \frac{5}{128}u_{j-2} + \frac{15}{32}u_{j-1} +
      87              :  *            \frac{45}{64}u_{j} -\frac{5}{32}u_{j+1} + \frac{3}{128}u_{j+2}, \\
      88              :  *  u_{j+1/2} &= \frac{3}{128}u_{j-2} - \frac{5}{32}u_{j-1} + \frac{45}{64}u_{j}
      89              :  *            + \frac{15}{32}u_{j+1} - \frac{5}{128}u_{j+2}.
      90              :  * \f}
      91              :  *
      92              :  * Degree 6 and 8 reconstruction is also supported.
      93              :  */
      94              : template <size_t Degree, size_t Dim>
      95            1 : void unlimited(
      96              :     const gsl::not_null<std::array<gsl::span<double>, Dim>*>
      97              :     reconstructed_upper_side_of_face_vars,
      98              :     const gsl::not_null<std::array<gsl::span<double>, Dim>*>
      99              :     reconstructed_lower_side_of_face_vars,
     100              :     const gsl::span<const double>& volume_vars,
     101              :     const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
     102              :     const Index<Dim>& volume_extents, const size_t number_of_variables) {
     103              :   detail::reconstruct<detail::UnlimitedReconstructor<Degree>>(
     104              :       reconstructed_upper_side_of_face_vars,
     105              :       reconstructed_lower_side_of_face_vars, volume_vars, ghost_cell_vars,
     106              :       volume_extents, number_of_variables);
     107              : }
     108              : }  // namespace fd::reconstruction