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