MonotisedCentral.hpp
1 // 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"
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 MonotisedCentralReconstructor {
28  const double* const q, const int stride) noexcept {
29  using std::abs;
30  using std::min;
31  // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
32  const double a = q[stride] - q[0];
33  // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
34  const double b = q[0] - q[-stride];
35  const double slope = 0.5 * (sgn(a) + sgn(b)) *
36  min(0.5 * abs(a + b), 2.0 * min(abs(a), abs(b)));
37  // NOLINTNEXTLINE(cppcoreguidelines-pro-bounds-pointer-arithmetic)
38  return {{q[0] - 0.5 * slope, q[0] + 0.5 * slope}};
39  }
40 
41  SPECTRE_ALWAYS_INLINE static constexpr size_t stencil_width() noexcept {
42  return 3;
43  }
44 };
45 } // namespace detail
46 
47 /*!
48  * \ingroup FiniteDifferenceGroup
49  * \brief Performs monotised central-difference reconstruction on the `vars` in
50  * each direction.
51  *
52  * On a 1d mesh with spacing \f$\Delta x\f$ we denote the solution at the
53  * \f$j\f$th point by \f$u_j\f$. The monotised central-difference reconstruction
54  * \cite RezzollaBook first computes:
55  *
56  * \f{align}
57  * \sigma_j\equiv \textrm{minmod}
58  * \left(\frac{u_{j+1} - u_{j-1}}{2\Delta x},
59  * 2\frac{u_j-u_{j-1}}{\Delta x},
60  * 2\frac{u_{j+1}-u_j}{\Delta x}\right)
61  * \f}
62  *
63  * where
64  *
65  * \f{align}
66  * \mathrm{minmod}(a,b,c)=
67  * \left\{
68  * \begin{array}{ll}
69  * \mathrm{sgn}(a)\min(\lvert a\rvert, \lvert b\rvert, \lvert c\rvert)
70  * & \mathrm{if} \; \mathrm{sgn}(a)=\mathrm{sgn}(b)=\mathrm{sgn}(c) \\
71  * 0 & \mathrm{otherwise}
72  * \end{array}\right.
73  * \f}
74  *
75  * The reconstructed solution \f$u_j(x)\f$ in the \f$j\f$th cell is given by
76  *
77  * \f{align}
78  * u_j(x)=u_j + \sigma_j (x-x_j)
79  * \f}
80  *
81  * where \f$x_j\f$ is the coordinate \f$x\f$ at the center of the cell.
82  */
83 template <size_t Dim>
86  reconstructed_upper_side_of_face_vars,
88  reconstructed_lower_side_of_face_vars,
89  const gsl::span<const double>& volume_vars,
90  const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
91  const Index<Dim>& volume_extents,
92  const size_t number_of_variables) noexcept {
93  detail::reconstruct<detail::MonotisedCentralReconstructor>(
94  reconstructed_upper_side_of_face_vars,
95  reconstructed_lower_side_of_face_vars, volume_vars, ghost_cell_vars,
96  volume_extents, number_of_variables);
97 }
98 } // namespace fd::reconstruction
utility
fd::reconstruction
Variable and flux vector splitting reconstruction schemes for finite difference methods.
Definition: AoWeno.hpp:25
fd::reconstruction::monotised_central
void monotised_central(const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_upper_side_of_face_vars, const gsl::not_null< std::array< gsl::span< double >, Dim > * > reconstructed_lower_side_of_face_vars, const gsl::span< const double > &volume_vars, const DirectionMap< Dim, gsl::span< const double >> &ghost_cell_vars, const Index< Dim > &volume_extents, const size_t number_of_variables) noexcept
Performs monotised central-difference reconstruction on the vars in each direction.
Definition: MonotisedCentral.hpp:84
Index
Definition: Index.hpp:31
cmath
sgn
constexpr T sgn(const T &val) noexcept
Compute the sign function of val defined as 1 if val > 0, 0 if val == 0, and -1 if val < 0.
Definition: Math.hpp:163
Direction
Definition: Direction.hpp:23
SPECTRE_ALWAYS_INLINE
#define SPECTRE_ALWAYS_INLINE
Definition: ForceInline.hpp:16
cstddef
array
std::min
T min(T... args)
gsl::span
Create a span/view on a range, which is cheap to copy (one pointer).
Definition: Gsl.hpp:292
DirectionMap
Definition: DirectionMap.hpp:15
Gsl.hpp
ForceInline.hpp
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecThirdOrderPiecewisePolynomial.hpp:13