SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/LinearOperators - WeakDivergence.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 3 4 75.0 %
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             : #include <type_traits>
       8             : 
       9             : #include "DataStructures/Matrix.hpp"
      10             : #include "DataStructures/Tensor/Tensor.hpp"
      11             : #include "DataStructures/Transpose.hpp"
      12             : #include "DataStructures/Variables.hpp"
      13             : #include "NumericalAlgorithms/LinearOperators/Divergence.hpp"
      14             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      15             : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
      16             : #include "Utilities/Blas.hpp"
      17             : #include "Utilities/ContainerHelpers.hpp"
      18             : #include "Utilities/Gsl.hpp"
      19             : #include "Utilities/Literals.hpp"
      20             : #include "Utilities/TMPL.hpp"
      21             : 
      22             : /*!
      23             :  * \ingroup NumericalAlgorithmsGroup
      24             :  * \brief Compute the weak form divergence of fluxes
      25             :  *
      26             :  * In a discontinuous Galerkin scheme we integrate the equations against the
      27             :  * basis functions over the element. For the flux divergence term this gives:
      28             :  *
      29             :  * \f{align*}{
      30             :  *  \int_{\Omega}d^n x \phi_{\breve{\imath}}\partial_i F^i,
      31             :  * \f}
      32             :  *
      33             :  * where the basis functions are denoted by \f$\phi_{\breve{\imath}}\f$.
      34             :  *
      35             :  * Integrating by parts we get
      36             :  *
      37             :  * \f{align*}{
      38             :  *  \int_{\Omega}d^n x\, \phi_{\breve{\imath}}\partial_i F^i
      39             :  *  = -\int_{\Omega}d^n x\,  F^i \partial_i \phi_{\breve{\imath}}
      40             :  *   + \int_{\partial\Omega} d^{(n-1)}\Sigma\, n_i F^i \phi_{\breve{\imath}}
      41             :  * \f}
      42             :  *
      43             :  * Next we expand the flux \f$F^i\f$ in terms of the basis functions, yielding
      44             :  *
      45             :  * \f{align*}{
      46             :  *  - \int_{\Omega}d^n x\,F^i_{\breve{\jmath}} \phi_{\breve{\jmath}} \partial_i
      47             :  *    \phi_{\breve{\imath}}
      48             :  *  + \int_{\partial\Omega} d^{(n-1)}\Sigma\, n_i F^i_{\breve{\jmath}}
      49             :  *    \phi_{\breve{\jmath}} \phi_{\breve{\imath}}
      50             :  * \f}
      51             :  *
      52             :  * This function computes the volume term:
      53             :  *
      54             :  * \f{align*}{
      55             :  * \frac{1}{w_\breve{\imath}} \int_{\Omega}d^n x \,
      56             :  *   F^i_{\breve{\jmath}} \phi_{\breve{\jmath}} \partial_i \phi_{\breve{\imath}}
      57             :  * \f}
      58             :  *
      59             :  * \note When using Gauss-Lobatto points the numerical values of the
      60             :  * divergence are the same for the strong and weak divergence at the interior
      61             :  * points. When using Gauss points they are only the same at the central grid
      62             :  * point and only when an odd number of grid points is used.
      63             :  */
      64             : template <typename... ResultTags, typename... FluxTags, size_t Dim>
      65           1 : void weak_divergence(
      66             :     const gsl::not_null<Variables<tmpl::list<ResultTags...>>*>
      67             :         divergence_of_fluxes,
      68             :     const Variables<tmpl::list<FluxTags...>>& fluxes, const Mesh<Dim>& mesh,
      69             :     const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
      70             :                           Frame::Inertial>& det_jac_times_inverse_jacobian) {
      71             :   if (UNLIKELY(divergence_of_fluxes->number_of_grid_points() !=
      72             :                fluxes.number_of_grid_points())) {
      73             :     divergence_of_fluxes->initialize(fluxes.number_of_grid_points());
      74             :   }
      75             : 
      76             :   const auto apply_matrix_in_first_dim =
      77             :       [](double* result, const double* const input, const Matrix& matrix,
      78             :          const size_t size, const bool add_to_result) {
      79             :         dgemm_<true>(
      80             :             'N', 'N',
      81             :             matrix.rows(),            // rows of matrix and result
      82             :             size / matrix.columns(),  // columns of result and input
      83             :             matrix.columns(),         // columns of matrix and rows of input
      84             :             1.0,                      // overall multiplier
      85             :             matrix.data(),            // matrix
      86             :             matrix.spacing(),         // rows of matrix including padding
      87             :             input,                    // input
      88             :             matrix.columns(),         // rows of input
      89             :             add_to_result
      90             :                 ? 1.0
      91             :                 : 0.0,  // 1.0 means add to result, 0.0 means overwrite result
      92             :             result,     // result
      93             :             matrix.rows());  // rows of result
      94             :       };
      95             : 
      96             :   if constexpr (Dim == 1) {
      97             :     (void)det_jac_times_inverse_jacobian;  // is identically 1.0 in 1d
      98             : 
      99             :     apply_matrix_in_first_dim(divergence_of_fluxes->data(), fluxes.data(),
     100             :                               Spectral::weak_flux_differentiation_matrix(mesh),
     101             :                               fluxes.size(), false);
     102             :   } else {
     103             :     // Multiplies the flux by det_jac_time_inverse_jacobian.
     104             :     const auto transform_to_logical_frame =
     105             :         [&det_jac_times_inverse_jacobian, &fluxes](
     106             :             auto flux_tag_v, auto div_tag_v,
     107             :             const gsl::not_null<Variables<tmpl::list<ResultTags...>>*>
     108             :                 result_buffer,
     109             :             auto logical_index_of_jacobian) {
     110             :           using flux_tag = tmpl::type_from<decltype(flux_tag_v)>;
     111             :           using div_tag = tmpl::type_from<decltype(div_tag_v)>;
     112             :           using TensorIndexType =
     113             :               std::array<size_t, div_tag::type::num_tensor_indices>;
     114             :           using FluxIndexType =
     115             :               std::array<size_t, flux_tag::type::num_tensor_indices>;
     116             : 
     117             :           auto& result = get<div_tag>(*result_buffer);
     118             :           const auto& flux = get<flux_tag>(fluxes);
     119             : 
     120             :           for (size_t result_storage_index = 0;
     121             :                result_storage_index < result.size(); ++result_storage_index) {
     122             :             const TensorIndexType result_tensor_index =
     123             :                 result.get_tensor_index(result_storage_index);
     124             :             const FluxIndexType flux_x_tensor_index =
     125             :                 prepend(result_tensor_index, 0_st);
     126             :             const FluxIndexType flux_y_tensor_index =
     127             :                 prepend(result_tensor_index, 1_st);
     128             :             if constexpr (Dim == 2) {
     129             :               result[result_storage_index] =
     130             :                   get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
     131             :                       0>(det_jac_times_inverse_jacobian) *
     132             :                       flux.get(flux_x_tensor_index) +
     133             :                   get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
     134             :                       1>(det_jac_times_inverse_jacobian) *
     135             :                       flux.get(flux_y_tensor_index);
     136             :             } else {
     137             :               const FluxIndexType flux_z_tensor_index =
     138             :                   prepend(result_tensor_index, 2_st);
     139             :               result[result_storage_index] =
     140             :                   get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
     141             :                       0>(det_jac_times_inverse_jacobian) *
     142             :                       flux.get(flux_x_tensor_index) +
     143             :                   get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
     144             :                       1>(det_jac_times_inverse_jacobian) *
     145             :                       flux.get(flux_y_tensor_index) +
     146             :                   get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
     147             :                       2>(det_jac_times_inverse_jacobian) *
     148             :                       flux.get(flux_z_tensor_index);
     149             :             }
     150             :           }
     151             :         };
     152             : 
     153             :     if constexpr (Dim == 2) {
     154             :       Variables<tmpl::list<ResultTags...>> data_buffer{
     155             :           divergence_of_fluxes->number_of_grid_points()};
     156             : 
     157             :       const Matrix& eta_weak_div_matrix =
     158             :           Spectral::weak_flux_differentiation_matrix(mesh.slice_through(1));
     159             :       const Matrix& xi_weak_div_matrix =
     160             :           Spectral::weak_flux_differentiation_matrix(mesh.slice_through(0));
     161             : 
     162             :       // Compute the eta divergence term. Since that also needs a transpose,
     163             :       // copy into result, then transpose into `data_buffer`
     164             :       EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
     165             :           tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
     166             :           make_not_null(&data_buffer), std::integral_constant<size_t, 1>{}));
     167             :       double* div_ptr = divergence_of_fluxes->data();
     168             :       raw_transpose(make_not_null(div_ptr), data_buffer.data(),
     169             :                     xi_weak_div_matrix.rows(),
     170             :                     divergence_of_fluxes->size() / xi_weak_div_matrix.rows());
     171             :       apply_matrix_in_first_dim(data_buffer.data(),
     172             :                                 divergence_of_fluxes->data(),
     173             :                                 eta_weak_div_matrix, data_buffer.size(), false);
     174             : 
     175             :       const size_t chunk_size = Variables<tmpl::list<Tags::div<FluxTags>...>>::
     176             :                                     number_of_independent_components *
     177             :                                 eta_weak_div_matrix.rows();
     178             :       raw_transpose(make_not_null(div_ptr), data_buffer.data(), chunk_size,
     179             :                     data_buffer.size() / chunk_size);
     180             : 
     181             :       // Now compute xi divergence and *add* to eta divergence
     182             :       EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
     183             :           tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
     184             :           make_not_null(&data_buffer), std::integral_constant<size_t, 0>{}));
     185             :       apply_matrix_in_first_dim(divergence_of_fluxes->data(),
     186             :                                 data_buffer.data(), xi_weak_div_matrix,
     187             :                                 data_buffer.size(), true);
     188             :     } else if constexpr (Dim == 3) {
     189             :       Variables<tmpl::list<ResultTags...>> data_buffer0{
     190             :           divergence_of_fluxes->number_of_grid_points()};
     191             :       Variables<tmpl::list<ResultTags...>> data_buffer1{
     192             :           divergence_of_fluxes->number_of_grid_points()};
     193             :       constexpr size_t number_of_independent_components =
     194             :           decltype(data_buffer1)::number_of_independent_components;
     195             : 
     196             :       const Matrix& zeta_weak_div_matrix =
     197             :           Spectral::weak_flux_differentiation_matrix(mesh.slice_through(2));
     198             :       const Matrix& eta_weak_div_matrix =
     199             :           Spectral::weak_flux_differentiation_matrix(mesh.slice_through(1));
     200             :       const Matrix& xi_weak_div_matrix =
     201             :           Spectral::weak_flux_differentiation_matrix(mesh.slice_through(0));
     202             : 
     203             :       // Compute the zeta divergence term. Since that also needs a transpose,
     204             :       // copy into data_buffer0, then transpose into `data_buffer1`.
     205             :       EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
     206             :           tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
     207             :           make_not_null(&data_buffer0), std::integral_constant<size_t, 2>{}));
     208             :       size_t chunk_size =
     209             :           xi_weak_div_matrix.rows() * eta_weak_div_matrix.rows();
     210             :       double* result_ptr = data_buffer1.data();
     211             :       raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
     212             :                     data_buffer0.size() / chunk_size);
     213             :       apply_matrix_in_first_dim(data_buffer0.data(), data_buffer1.data(),
     214             :                                 zeta_weak_div_matrix, data_buffer0.size(),
     215             :                                 false);
     216             :       chunk_size =
     217             :           number_of_independent_components * zeta_weak_div_matrix.rows();
     218             :       result_ptr = divergence_of_fluxes->data();
     219             :       raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
     220             :                     data_buffer1.size() / chunk_size);
     221             : 
     222             :       // Compute the eta divergence term. Since that also needs a transpose,
     223             :       // copy into data_buffer0, then transpose into `data_buffer1`.
     224             :       EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
     225             :           tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
     226             :           make_not_null(&data_buffer0), std::integral_constant<size_t, 1>{}));
     227             :       chunk_size = xi_weak_div_matrix.rows();
     228             :       result_ptr = data_buffer1.data();
     229             :       raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
     230             :                     data_buffer1.size() / chunk_size);
     231             :       apply_matrix_in_first_dim(data_buffer0.data(), data_buffer1.data(),
     232             :                                 eta_weak_div_matrix, data_buffer0.size(),
     233             :                                 false);
     234             :       chunk_size = number_of_independent_components *
     235             :                    eta_weak_div_matrix.rows() * zeta_weak_div_matrix.rows();
     236             :       result_ptr = data_buffer1.data();
     237             :       raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
     238             :                     data_buffer0.size() / chunk_size);
     239             :       *divergence_of_fluxes += data_buffer1;
     240             : 
     241             :       // Now compute xi divergence and *add* to eta divergence
     242             :       EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
     243             :           tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
     244             :           make_not_null(&data_buffer0), std::integral_constant<size_t, 0>{}));
     245             :       apply_matrix_in_first_dim(divergence_of_fluxes->data(),
     246             :                                 data_buffer0.data(), xi_weak_div_matrix,
     247             :                                 data_buffer0.size(), true);
     248             :     } else {
     249             :       static_assert(Dim == 1 or Dim == 2 or Dim == 3,
     250             :                     "Weak divergence only implemented in 1d, 2d, and 3d.");
     251             :     }
     252             :   }
     253             : }
     254             : 
     255             : /// @{
     256             : /*!
     257             :  * \brief Compute the weak divergence of fluxes in logical coordinates
     258             :  *
     259             :  * Applies the transpose logical differentiation matrix to the fluxes in each
     260             :  * dimension and sums over dimensions.
     261             :  *
     262             :  * \see weak_divergence
     263             :  */
     264             : template <typename ResultTensor, typename FluxTensor, size_t Dim>
     265           1 : void logical_weak_divergence(gsl::not_null<ResultTensor*> div_flux,
     266             :                              const FluxTensor& flux, const Mesh<Dim>& mesh);
     267             : 
     268             : /*!
     269             :  * \param divergence_of_fluxes Result buffer. Will be resized if necessary.
     270             :  * \param fluxes The fluxes to take the divergence of. Their first index must be
     271             :  * an upper spatial index in element-logical coordinates.
     272             :  * \param mesh The element mesh.
     273             :  * \param add_to_buffer If `true`, add the divergence to the existing values in
     274             :  * `divergence_of_fluxes`. Otherwise, overwrite the existing values.
     275             :  */
     276             : template <typename... ResultTags, typename... FluxTags, size_t Dim>
     277           1 : void logical_weak_divergence(
     278             :     const gsl::not_null<Variables<tmpl::list<ResultTags...>>*>
     279             :         divergence_of_fluxes,
     280             :     const Variables<tmpl::list<FluxTags...>>& fluxes, const Mesh<Dim>& mesh,
     281             :     const bool add_to_buffer = false) {
     282             :   static_assert(
     283             :       (... and
     284             :        std::is_same_v<tmpl::front<typename FluxTags::type::index_list>,
     285             :                       SpatialIndex<Dim, UpLo::Up, Frame::ElementLogical>>),
     286             :       "The first index of each flux must be an upper spatial index in "
     287             :       "element-logical coordinates.");
     288             :   static_assert(
     289             :       (... and
     290             :        std::is_same_v<
     291             :            typename ResultTags::type,
     292             :            TensorMetafunctions::remove_first_index<typename FluxTags::type>>),
     293             :       "The result tensors must have the same type as the flux tensors with "
     294             :       "their first index removed.");
     295             :   if (not add_to_buffer) {
     296             :     divergence_of_fluxes->initialize(mesh.number_of_grid_points(), 0.);
     297             :   }
     298             :   EXPAND_PACK_LEFT_TO_RIGHT(logical_weak_divergence(
     299             :       make_not_null(&get<ResultTags>(*divergence_of_fluxes)),
     300             :       get<FluxTags>(fluxes), mesh));
     301             : }
     302             : /// @}

Generated by: LCOV version 1.14