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 "NumericalAlgorithms/Spectral/Basis.hpp"
17 : #include "Utilities/Gsl.hpp"
18 : #include "Utilities/TMPL.hpp"
19 :
20 : namespace dg {
21 : /// @{
22 : /// \ingroup DiscontinuousGalerkinGroup
23 : /// \brief Lifts the flux contribution from an interface to the volume.
24 : ///
25 : /// The lifting operation takes the (d-1)-dimensional flux term at the
26 : /// interface and computes the corresponding d-dimensional term in the
27 : /// volume. SpECTRE implements an efficient DG method in which each
28 : /// interface grid point contributes only to that same grid point of the
29 : /// volume.
30 : ///
31 : /// \details
32 : /// SpECTRE implements a DG method with a diagonalized mass matrix (also
33 : /// known as a mass-lumping scheme). This choice gives a large
34 : /// reduction in the computational cost of the lifting operation, however,
35 : /// the scheme is slightly less accurate, especially when the grid is
36 : /// deformed by non-trivial Jacobians. For more details on the
37 : /// diagonalization of the mass matrix and its implications,
38 : /// \cite Teukolsky2015ega, especially Section 3.
39 : ///
40 : /// This function is implemented to handle Legendre, Chebyshev, and Zernike
41 : /// bases.
42 : ///
43 : /// \note The result is still provided only on the boundary grid. The
44 : /// values away from the boundary are zero and are not stored.
45 : template <typename... BoundaryCorrectionTags>
46 1 : void lift_flux(
47 : const gsl::not_null<Variables<tmpl::list<BoundaryCorrectionTags...>>*>
48 : boundary_correction_terms,
49 : const size_t extent_perpendicular_to_boundary,
50 : const Scalar<DataVector>& magnitude_of_face_normal,
51 : const Spectral::Basis basis = Spectral::Basis::Legendre) {
52 : if (basis == Spectral::Basis::Legendre) {
53 : // For an Nth degree basis (i.e., one with N+1 basis functions), the LGL
54 : // weights are:
55 : // w_i = 2 / ((N + 1) * N * (P_{N}(xi_i))^2)
56 : // and so at the end points (xi = +/- 1) we get:
57 : // w_0 = 2 / ((N + 1) * N)
58 : // w_N = 2 / ((N + 1) * N)
59 : // The negative sign comes from bringing the boundary correction over to the
60 : // RHS of the equal sign (e.g. `du/dt=Source - div Flux - boundary corr`),
61 : // while the magnitude of the normal vector above accounts for the ratios of
62 : // spatial metrics and Jacobians.
63 : *boundary_correction_terms *=
64 : -0.5 *
65 : static_cast<double>((extent_perpendicular_to_boundary *
66 : (extent_perpendicular_to_boundary - 1))) *
67 : get(magnitude_of_face_normal);
68 : } else if (basis == Spectral::Basis::Chebyshev) {
69 : // Chebyshev with GaussLobotto, with N grid points, has the weight at
70 : // the upper and lower sides
71 : // w_0 = 1 / ((N - 1)^2 - a)
72 : // w_N = 1 / ((N - 1)^2 - a)
73 : // where a = 0 if N is even, else it is 1
74 : // The negative sign follows from the Legendre discussion
75 : *boundary_correction_terms *=
76 : -static_cast<double>((extent_perpendicular_to_boundary - 1) *
77 : (extent_perpendicular_to_boundary - 1) -
78 : extent_perpendicular_to_boundary % 2) *
79 : get(magnitude_of_face_normal);
80 : } else if (basis == Spectral::Basis::ZernikeB1) {
81 : // ZernikeB1 with GaussRadauUpper, with N grid points, has the weight at
82 : // the upper side
83 : // w_N = 1 / (N * (2 * N - 1))
84 : // The negative sign follows from the Legendre discussion
85 : *boundary_correction_terms *=
86 : -1.0 *
87 : static_cast<double>((extent_perpendicular_to_boundary *
88 : (2 * extent_perpendicular_to_boundary - 1))) *
89 : get(magnitude_of_face_normal);
90 : } else if (basis == Spectral::Basis::ZernikeB2) {
91 : // ZernikeB2 with GaussRadauUpper, with N grid points, has the weight at
92 : // the upper side
93 : // w_N = 1 / (2 * N^2)
94 : // The negative sign follows from the Legendre discussion
95 : *boundary_correction_terms *=
96 : -1.0 *
97 : static_cast<double>((2 * extent_perpendicular_to_boundary *
98 : extent_perpendicular_to_boundary)) *
99 : get(magnitude_of_face_normal);
100 : } else if (basis == Spectral::Basis::ZernikeB3) {
101 : // ZernikeB3 with GaussRadauUpper, with N grid points, has the weight at
102 : // the upper side
103 : // w_N = 1 / (N * (2 * N + 1))
104 : // The negative sign follows from the Legendre discussion
105 : *boundary_correction_terms *=
106 : -1.0 *
107 : static_cast<double>((extent_perpendicular_to_boundary *
108 : (2 * extent_perpendicular_to_boundary + 1))) *
109 : get(magnitude_of_face_normal);
110 : } else {
111 : ERROR("Got unexpected basis " << basis);
112 : }
113 : }
114 :
115 : template <typename... FluxTags>
116 1 : auto lift_flux(Variables<tmpl::list<FluxTags...>> flux,
117 : const size_t extent_perpendicular_to_boundary,
118 : const Scalar<DataVector>& magnitude_of_face_normal,
119 : const Spectral::Basis basis = Spectral::Basis::Legendre)
120 : -> Variables<tmpl::list<db::remove_tag_prefix<FluxTags>...>> {
121 : Variables<tmpl::list<db::remove_tag_prefix<FluxTags>...>> lifted_data(
122 : std::move(flux));
123 : lift_flux(make_not_null(&lifted_data), extent_perpendicular_to_boundary,
124 : magnitude_of_face_normal, basis);
125 : return lifted_data;
126 : }
127 : /// @}
128 : } // namespace dg
|
