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