Line data Source code
1 1 : // Distributed under the MIT License. 2 : // See LICENSE.txt for details. 3 : 4 : /// \file 5 : /// Defines function lift_flux. 6 : 7 : #pragma once 8 : 9 : #include <cstddef> 10 : #include <utility> 11 : 12 : #include "DataStructures/DataBox/PrefixHelpers.hpp" 13 : #include "DataStructures/DataVector.hpp" 14 : #include "DataStructures/Tensor/Tensor.hpp" 15 : #include "DataStructures/Variables.hpp" 16 : #include "Utilities/Gsl.hpp" 17 : #include "Utilities/TMPL.hpp" 18 : 19 : namespace dg { 20 : /// @{ 21 : /// \ingroup DiscontinuousGalerkinGroup 22 : /// \brief Lifts the flux contribution from an interface to the volume. 23 : /// 24 : /// The lifting operation takes the (d-1)-dimensional flux term at the 25 : /// interface and computes the corresponding d-dimensional term in the 26 : /// volume. SpECTRE implements an efficient DG method in which each 27 : /// interface grid point contributes only to that same grid point of the 28 : /// volume. 29 : /// 30 : /// \details 31 : /// SpECTRE implements a DG method with a diagonalized mass matrix (also 32 : /// known as a mass-lumping scheme). This choice gives a large 33 : /// reduction in the computational cost of the lifting operation, however, 34 : /// the scheme is slightly less accurate, especially when the grid is 35 : /// deformed by non-trivial Jacobians. For more details on the 36 : /// diagonalization of the mass matrix and its implications, 37 : /// \cite Teukolsky2015ega, especially Section 3. 38 : /// 39 : /// \note The result is still provided only on the boundary grid. The 40 : /// values away from the boundary are zero and are not stored. 41 : template <typename... BoundaryCorrectionTags> 42 1 : void lift_flux( 43 : const gsl::not_null<Variables<tmpl::list<BoundaryCorrectionTags...>>*> 44 : boundary_correction_terms, 45 : const size_t extent_perpendicular_to_boundary, 46 : const Scalar<DataVector>& magnitude_of_face_normal) { 47 : // For an Nth degree basis (i.e., one with N+1 basis functions), the LGL 48 : // weights are: 49 : // w_i = 2 / ((N + 1) * N * (P_{N}(xi_i))^2) 50 : // and so at the end points (xi = +/- 1) we get: 51 : // w_0 = 2 / ((N + 1) * N) 52 : // w_N = 2 / ((N + 1) * N) 53 : // The negative sign comes from bringing the boundary correction over to the 54 : // RHS of the equal sign (e.g. `du/dt=Source - div Flux - boundary corr`), 55 : // while the magnitude of the normal vector above accounts for the ratios of 56 : // spatial metrics and Jacobians. 57 : *boundary_correction_terms *= 58 : -0.5 * 59 : static_cast<double>((extent_perpendicular_to_boundary * 60 : (extent_perpendicular_to_boundary - 1))) * 61 : get(magnitude_of_face_normal); 62 : } 63 : 64 : template <typename... FluxTags> 65 1 : auto lift_flux(Variables<tmpl::list<FluxTags...>> flux, 66 : const size_t extent_perpendicular_to_boundary, 67 : const Scalar<DataVector>& magnitude_of_face_normal) 68 : -> Variables<tmpl::list<db::remove_tag_prefix<FluxTags>...>> { 69 : Variables<tmpl::list<db::remove_tag_prefix<FluxTags>...>> lifted_data( 70 : std::move(flux)); 71 : lift_flux(make_not_null(&lifted_data), extent_perpendicular_to_boundary, 72 : magnitude_of_face_normal); 73 : return lifted_data; 74 : } 75 : /// @} 76 : } // namespace dg