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