SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/DiscontinuousGalerkin - ApplyMassMatrix.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 5 6 83.3 %
Date: 2024-03-29 00:33:31
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/Variables.hpp"
       9             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      10             : #include "Utilities/ErrorHandling/Assert.hpp"
      11             : #include "Utilities/Gsl.hpp"
      12             : 
      13           1 : namespace dg {
      14             : 
      15             : namespace detail {
      16             : template <size_t Dim>
      17             : void apply_mass_matrix_impl(gsl::not_null<double*> data, const Mesh<Dim>& mesh);
      18             : template <size_t Dim>
      19             : void apply_inverse_mass_matrix_impl(gsl::not_null<double*> data,
      20             :                                     const Mesh<Dim>& mesh);
      21             : }  // namespace detail
      22             : 
      23             : /*!
      24             :  * \brief Apply the DG mass matrix to the data, in the diagonal mass-matrix
      25             :  * approximation ("mass-lumping")
      26             :  *
      27             :  * The DG mass matrix is:
      28             :  *
      29             :  * \f{equation}
      30             :  * M_{pq} = \int_{\Omega_k} \psi_p(\xi) \psi_q(\xi) \mathrm{d}V
      31             :  * \f}
      32             :  *
      33             :  * where \f$\psi_p(\xi)\f$ are the basis functions on the element
      34             :  * \f$\Omega_k\f$. In the diagonal mass-matrix approximation ("mass-lumping") we
      35             :  * evaluate the integral directly on the collocation points, i.e. with a Gauss
      36             :  * or Gauss-Lobatto quadrature determined by the element mesh. Then it reduces
      37             :  * to:
      38             :  *
      39             :  * \f{equation}
      40             :  * M_{pq} \approx \delta_{pq} \prod_{i=1}^d w_{p_i}
      41             :  * \f}
      42             :  *
      43             :  * where \f$d\f$ is the spatial dimension and \f$w_{p_i}\f$ are the quadrature
      44             :  * weights in dimension \f$i\f$. To apply the mass matrix in coordinates
      45             :  * different than logical, or to account for a curved background metric, the
      46             :  * data can be pre-multiplied with the Jacobian determinant and/or a metric
      47             :  * determinant.
      48             :  *
      49             :  * \note The mass-lumping is exact on Legendre-Gauss meshes, but omits a
      50             :  * correction term on Legendre-Gauss-Lobatto meshes.
      51             :  */
      52             : /// @{
      53             : template <size_t Dim>
      54           1 : void apply_mass_matrix(const gsl::not_null<DataVector*> data,
      55             :                        const Mesh<Dim>& mesh) {
      56             :   ASSERT(data->size() == mesh.number_of_grid_points(),
      57             :          "The DataVector has size " << data->size() << ", but expected size "
      58             :                                     << mesh.number_of_grid_points()
      59             :                                     << " on the given mesh.");
      60             :   detail::apply_mass_matrix_impl(data->data(), mesh);
      61             : }
      62             : 
      63             : template <size_t Dim, typename TagsList>
      64           1 : void apply_mass_matrix(const gsl::not_null<Variables<TagsList>*> data,
      65             :                        const Mesh<Dim>& mesh) {
      66             :   const size_t num_points = data->number_of_grid_points();
      67             :   ASSERT(num_points == mesh.number_of_grid_points(),
      68             :          "The Variables data has "
      69             :              << num_points << " grid points, but expected "
      70             :              << mesh.number_of_grid_points() << " on the given mesh.");
      71             :   constexpr size_t num_comps =
      72             :       Variables<TagsList>::number_of_independent_components;
      73             :   for (size_t i = 0; i < num_comps; ++i) {
      74             :     detail::apply_mass_matrix_impl(data->data() + i * num_points, mesh);
      75             :   }
      76             : }
      77             : /// @}
      78             : 
      79             : /*!
      80             :  * \brief Apply the inverse DG mass matrix to the data
      81             :  *
      82             :  * \see dg::apply_mass_matrix
      83             :  */
      84             : /// @{
      85             : template <size_t Dim>
      86           1 : void apply_inverse_mass_matrix(const gsl::not_null<DataVector*> data,
      87             :                                const Mesh<Dim>& mesh) {
      88             :   ASSERT(data->size() == mesh.number_of_grid_points(),
      89             :          "The DataVector has size " << data->size() << ", but expected size "
      90             :                                     << mesh.number_of_grid_points()
      91             :                                     << " on the given mesh.");
      92             :   detail::apply_inverse_mass_matrix_impl(data->data(), mesh);
      93             : }
      94             : 
      95             : template <size_t Dim, typename TagsList>
      96           1 : void apply_inverse_mass_matrix(const gsl::not_null<Variables<TagsList>*> data,
      97             :                                const Mesh<Dim>& mesh) {
      98             :   const size_t num_points = data->number_of_grid_points();
      99             :   ASSERT(num_points == mesh.number_of_grid_points(),
     100             :          "The Variables data has "
     101             :              << num_points << " grid points, but expected "
     102             :              << mesh.number_of_grid_points() << " on the given mesh.");
     103             :   constexpr size_t num_comps =
     104             :       Variables<TagsList>::number_of_independent_components;
     105             :   for (size_t i = 0; i < num_comps; ++i) {
     106             :     detail::apply_inverse_mass_matrix_impl(data->data() + i * num_points, mesh);
     107             :   }
     108             : }
     109             : /// @}
     110             : 
     111             : }  // namespace dg

Generated by: LCOV version 1.14