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 DataType, typename Symm,
      27              :           typename... RemainingIndices,
      28              :           typename ResultTensor = Tensor<DataType, 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<DataType, 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 DataType, typename Symm,
      88              :           typename Indices>
      89            1 : void normal_times_flux(
      90              :     const gsl::not_null<TensorMetafunctions::prepend_spatial_index<
      91              :         Tensor<DataType, Symm, Indices>, VolumeDim, UpLo::Lo, Fr>*>
      92              :         normal_times_flux,
      93              :     const tnsr::i<DataVector, VolumeDim, Fr>& normal,
      94              :     const Tensor<DataType, Symm, Indices>& rhs) {
      95              :   outer_product(normal_times_flux, normal, rhs);
      96              : }
      97              : 
      98              : template <typename... ReturnTags, typename... FluxTags, size_t VolumeDim,
      99              :           typename Fr>
     100            1 : void normal_times_flux(
     101              :     const gsl::not_null<Variables<tmpl::list<ReturnTags...>>*> result,
     102              :     const tnsr::i<DataVector, VolumeDim, Fr>& normal,
     103              :     const Variables<tmpl::list<FluxTags...>>& fluxes) {
     104              :   if (result->number_of_grid_points() != fluxes.number_of_grid_points()) {
     105              :     result->initialize(fluxes.number_of_grid_points());
     106              :   }
     107              :   EXPAND_PACK_LEFT_TO_RIGHT(normal_times_flux(
     108              :       make_not_null(&get<ReturnTags>(*result)), normal, get<FluxTags>(fluxes)));
     109              : }
     110              : /// @}