Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <optional>
       9             : #include <string>
      10             : #include <tuple>
      11             : 
      12             : #include "DataStructures/ApplyMatrices.hpp"
      13             : #include "DataStructures/DataBox/DataBox.hpp"
      14             : #include "DataStructures/DataVector.hpp"
      15             : #include "DataStructures/Matrix.hpp"
      16             : #include "Domain/Creators/BlockGroups.hpp"
      17             : #include "Domain/Tags.hpp"
      18             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      19             : #include "Parallel/AlgorithmExecution.hpp"
      20             : #include "Parallel/GlobalCache.hpp"
      21             : #include "Utilities/Algorithm.hpp"
      22             : #include "Utilities/ErrorHandling/Error.hpp"
      23             : #include "Utilities/Gsl.hpp"
      24             : #include "Utilities/StdHelpers.hpp"
      25             : #include "Utilities/TMPL.hpp"
      26             : #include "Utilities/TypeTraits.hpp"
      27             : 
      28             : /// \cond
      29             : namespace domain::Tags {
      30             : template <size_t VolumeDim>
      31             : struct Domain;
      32             : }  // namespace domain::Tags
      33             : namespace Filters::Tags {
      34             : template <typename FilterType>
      35             : struct Filter;
      36             : }  // namespace Filters::Tags
      37             : /// \endcond
      38             : 
      39             : namespace dg {
      40           0 : namespace Actions {
      41             : namespace Filter_detail {
      42             : template <bool SameSize, bool SameList>
      43             : struct FilterAllEvolvedVars {
      44             :   template <typename EvolvedVarsTagList, typename... FilterTags>
      45             :   using f = std::integral_constant<bool, false>;
      46             : };
      47             : 
      48             : template <>
      49             : struct FilterAllEvolvedVars<true, false> {
      50             :   template <typename EvolvedVarsTagList, typename... FilterTags>
      51             :   using f =
      52             :       std::integral_constant<bool, tmpl2::flat_all_v<tmpl::list_contains_v<
      53             :                                        EvolvedVarsTagList, FilterTags>...>>;
      54             : };
      55             : 
      56             : template <>
      57             : struct FilterAllEvolvedVars<true, true> {
      58             :   template <typename EvolvedVarsTagList, typename... FilterTags>
      59             :   using f = std::integral_constant<bool, true>;
      60             : };
      61             : }  // namespace Filter_detail
      62             : 
      63             : /// \cond
      64             : template <typename FilterType, typename TagsToFilterList>
      65             : struct Filter;
      66             : /// \endcond
      67             : 
      68             : /*!
      69             :  * \ingroup DiscontinuousGalerkinGroup
      70             :  * \brief Applies a filter to the specified tags.
      71             :  *
      72             :  * If different Filters are desired for different tags then multiple `Filter`
      73             :  * actions must be inserted into the action list with different `FilterType`.
      74             :  * Here is an example of an action list with two different exponential filters:
      75             :  *
      76             :  * \snippet LinearOperators/Test_Filtering.cpp action_list_example
      77             :  *
      78             :  * Uses:
      79             :  * - GlobalCache:
      80             :  *   - `Filter`
      81             :  * - DataBox:
      82             :  *   - `Tags::Mesh`
      83             :  * - DataBox changes:
      84             :  *   - Adds: nothing
      85             :  *   - Removes: nothing
      86             :  *   - Modifies:
      87             :  *     - `TagsToFilter`
      88             :  * - System:
      89             :  *   - `volume_dim`
      90             :  *   - `variables_tag`
      91             :  *
      92             :  */
      93             : template <typename FilterType, typename... TagsToFilter>
      94           1 : class Filter<FilterType, tmpl::list<TagsToFilter...>> {
      95             :  public:
      96           0 :   using const_global_cache_tags =
      97             :       tmpl::list<::Filters::Tags::Filter<FilterType>>;
      98             : 
      99             :   template <typename DbTags, typename... InboxTags, typename ArrayIndex,
     100             :             typename ActionList, typename ParallelComponent,
     101             :             typename Metavariables>
     102           0 :   static Parallel::iterable_action_return_t apply(
     103             :       db::DataBox<DbTags>& box,
     104             :       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
     105             :       const Parallel::GlobalCache<Metavariables>& cache,
     106             :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
     107             :       const ParallelComponent* const /*meta*/) {
     108             :     constexpr size_t volume_dim = Metavariables::system::volume_dim;
     109             :     using evolved_vars_tag = typename Metavariables::system::variables_tag;
     110             :     using evolved_vars_tags_list = typename evolved_vars_tag::tags_list;
     111             :     const FilterType& filter_helper =
     112             :         Parallel::get<::Filters::Tags::Filter<FilterType>>(cache);
     113             :     const size_t block_id =
     114             :         db::get<domain::Tags::Element<volume_dim>>(box).id().block_id();
     115             :     const auto& domain = Parallel::get<domain::Tags::Domain<volume_dim>>(cache);
     116             :     const auto& block_groups = domain.block_groups();
     117             :     const std::string& block_name = domain.blocks()[block_id].name();
     118             : 
     119             :     // Technically this whole next block could be done on a single line, but
     120             :     // then it would be very dense and hard to understand. This way is easier to
     121             :     // read and understand
     122             :     bool enable = filter_helper.enable();
     123             :     // Only do this check if filtering is enabled. A `nullopt` means all blocks
     124             :     // are allowed to do filtering
     125             :     if (enable and filter_helper.blocks_to_filter().has_value()) {
     126             :       // Enable filtering for this block if it's in any of the listed groups
     127             :       enable = alg::any_of(
     128             :           filter_helper.blocks_to_filter().value(),
     129             :           [&block_name, &block_groups](const std::string& block_to_filter) {
     130             :             return domain::block_is_in_group(block_name, block_to_filter,
     131             :                                              block_groups);
     132             :           });
     133             :     }
     134             : 
     135             :     if (not enable) {
     136             :       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     137             :     }
     138             : 
     139             :     const Mesh<volume_dim> mesh = db::get<domain::Tags::Mesh<volume_dim>>(box);
     140             :     const Matrix empty{};
     141             :     auto filter = make_array<volume_dim>(std::cref(empty));
     142             :     for (size_t d = 0; d < volume_dim; d++) {
     143             :       gsl::at(filter, d) =
     144             :           std::cref(filter_helper.filter_matrix(mesh.slice_through(d)));
     145             :     }
     146             : 
     147             :     // In the case that the tags we are filtering are all the evolved variables
     148             :     // we filter the entire Variables at once to be more efficient. This case is
     149             :     // the first branch of the `if-else`.
     150             :     if (Filter_detail::FilterAllEvolvedVars<
     151             :             sizeof...(TagsToFilter) ==
     152             :                 tmpl::size<evolved_vars_tags_list>::value,
     153             :             std::is_same_v<evolved_vars_tags_list,
     154             :                            tmpl::list<TagsToFilter...>>>::
     155             :             template f<evolved_vars_tags_list, TagsToFilter...>::value) {
     156             :       db::mutate<typename Metavariables::system::variables_tag>(
     157             :           [&filter](const gsl::not_null<
     158             :                         typename Metavariables::system::variables_tag::type*>
     159             :                         vars,
     160             :                     const auto& local_mesh) {
     161             :             *vars = apply_matrices(filter, *vars, local_mesh.extents());
     162             :           },
     163             :           make_not_null(&box), mesh);
     164             :     } else {
     165             :       db::mutate<TagsToFilter...>(
     166             :           [&filter](const gsl::not_null<
     167             :                         typename TagsToFilter::type*>... tensors_to_filter,
     168             :                     const auto& local_mesh) {
     169             :             DataVector temp(local_mesh.number_of_grid_points(), 0.0);
     170             :             const auto helper = [&local_mesh, &filter,
     171             :                                  &temp](const auto tensor) {
     172             :               for (auto& component : *tensor) {
     173             :                 temp = 0.0;
     174             :                 apply_matrices(make_not_null(&temp), filter, component,
     175             :                                local_mesh.extents());
     176             :                 component = temp;
     177             :               }
     178             :             };
     179             :             EXPAND_PACK_LEFT_TO_RIGHT(helper(tensors_to_filter));
     180             :           },
     181             :           make_not_null(&box), mesh);
     182             :     }
     183             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     184             :   }
     185             : };
     186             : }  // namespace Actions
     187             : }  // namespace dg