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/DataBox/PrefixHelpers.hpp"
       9             : #include "DataStructures/DataBox/Prefixes.hpp"
      10             : #include "DataStructures/Tensor/EagerMath/OuterProduct.hpp"
      11             : #include "DataStructures/Tensor/Tensor.hpp"
      12             : #include "DataStructures/Variables.hpp"
      13             : #include "Utilities/Gsl.hpp"
      14             : #include "Utilities/MakeWithValue.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : /// @{
      18             : /*!
      19             :  * \brief Contract a surface normal covector with the first index of a flux
      20             :  * tensor or variables
      21             :  *
      22             :  * \details
      23             :  * Returns \f$n_i F^i_{j\ldots}\f$, where the flux tensor \f$F\f$ must have an
      24             :  * upper spatial first index and may have arbitrary extra indices.
      25             :  */
      26             : template <size_t VolumeDim, typename Fr, typename Symm,
      27             :           typename... RemainingIndices,
      28             :           typename ResultTensor = Tensor<DataVector, tmpl::pop_front<Symm>,
      29             :                                          index_list<RemainingIndices...>>>
      30           1 : void normal_dot_flux(
      31             :     const gsl::not_null<ResultTensor*> normal_dot_flux,
      32             :     const tnsr::i<DataVector, VolumeDim, Fr>& normal,
      33             :     const Tensor<DataVector, Symm,
      34             :                  index_list<SpatialIndex<VolumeDim, UpLo::Up, Fr>,
      35             :                             RemainingIndices...>>& flux_tensor) {
      36             :   for (auto it = normal_dot_flux->begin(); it != normal_dot_flux->end(); it++) {
      37             :     const auto result_indices = normal_dot_flux->get_tensor_index(it);
      38             :     *it = get<0>(normal) * flux_tensor.get(prepend(result_indices, size_t{0}));
      39             :     for (size_t d = 1; d < VolumeDim; d++) {
      40             :       *it += normal.get(d) * flux_tensor.get(prepend(result_indices, d));
      41             :     }
      42             :   }
      43             : }
      44             : 
      45             : template <typename... ReturnTags, typename... FluxTags, size_t VolumeDim,
      46             :           typename Fr>
      47           1 : void normal_dot_flux(
      48             :     const gsl::not_null<Variables<tmpl::list<ReturnTags...>>*> result,
      49             :     const tnsr::i<DataVector, VolumeDim, Fr>& normal,
      50             :     const Variables<tmpl::list<FluxTags...>>& fluxes) {
      51             :   if (result->number_of_grid_points() != fluxes.number_of_grid_points()) {
      52             :     result->initialize(fluxes.number_of_grid_points());
      53             :   }
      54             :   EXPAND_PACK_LEFT_TO_RIGHT(normal_dot_flux(
      55             :       make_not_null(&get<ReturnTags>(*result)), normal, get<FluxTags>(fluxes)));
      56             : }
      57             : 
      58             : template <typename TagsList, size_t VolumeDim, typename Fr>
      59           1 : auto normal_dot_flux(
      60             :     const tnsr::i<DataVector, VolumeDim, Fr>& normal,
      61             :     const Variables<db::wrap_tags_in<::Tags::Flux, TagsList,
      62             :                                      tmpl::size_t<VolumeDim>, Fr>>& fluxes) {
      63             :   Variables<db::wrap_tags_in<::Tags::NormalDotFlux, TagsList>> result{};
      64             :   normal_dot_flux(make_not_null(&result), normal, fluxes);
      65             :   return result;
      66             : }
      67             : /// @}
      68             : 
      69             : /// @{
      70             : /*!
      71             :  * \brief Multiplies a surface normal covector with a tensor or variables
      72             :  *
      73             :  * \details
      74             :  * Returns the outer product $n_j v_{k\ldots}$, where $n_j$ is the `normal` and
      75             :  * $v_{k\ldots}$ is the `rhs`.
      76             :  *
      77             :  * Note that this quantity is a "normal dot flux" where the flux involves
      78             :  * a Kronecker delta. For example:
      79             :  *
      80             :  * \f{equation}
      81             :  * n_j v_{k\ldots} = n_i \delta^i_j v_{k\ldots} = n_i F^i_{jk\ldots}
      82             :  * \f}
      83             :  *
      84             :  * This makes this quantity useful for optimizations of DG formulations for
      85             :  * "sparse" (i.e., Kronecker delta) fluxes.
      86             :  */
      87             : template <size_t VolumeDim, typename Fr, typename Symm, typename Indices>
      88           1 : void normal_times_flux(
      89             :     const gsl::not_null<TensorMetafunctions::prepend_spatial_index<
      90             :         Tensor<DataVector, Symm, Indices>, VolumeDim, UpLo::Lo, Fr>*>
      91             :         normal_times_flux,
      92             :     const tnsr::i<DataVector, VolumeDim, Fr>& normal,
      93             :     const Tensor<DataVector, Symm, Indices>& rhs) {
      94             :   outer_product(normal_times_flux, normal, rhs);
      95             : }
      96             : 
      97             : template <typename... ReturnTags, typename... FluxTags, size_t VolumeDim,
      98             :           typename Fr>
      99           1 : void normal_times_flux(
     100             :     const gsl::not_null<Variables<tmpl::list<ReturnTags...>>*> result,
     101             :     const tnsr::i<DataVector, VolumeDim, Fr>& normal,
     102             :     const Variables<tmpl::list<FluxTags...>>& fluxes) {
     103             :   if (result->number_of_grid_points() != fluxes.number_of_grid_points()) {
     104             :     result->initialize(fluxes.number_of_grid_points());
     105             :   }
     106             :   EXPAND_PACK_LEFT_TO_RIGHT(normal_times_flux(
     107             :       make_not_null(&get<ReturnTags>(*result)), normal, get<FluxTags>(fluxes)));
     108             : }
     109             : /// @}