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 <optional>
       8             : #include <utility>
       9             : 
      10             : #include "DataStructures/DataBox/DataBox.hpp"
      11             : #include "DataStructures/DataVector.hpp"
      12             : #include "DataStructures/Tensor/Tensor.hpp"
      13             : #include "Domain/Tags.hpp"
      14             : #include "NumericalAlgorithms/LinearOperators/Filters/Filter.hpp"
      15             : #include "NumericalAlgorithms/LinearOperators/Filters/None.hpp"
      16             : #include "NumericalAlgorithms/LinearOperators/Filters/Tag.hpp"
      17             : #include "Parallel/AlgorithmExecution.hpp"
      18             : #include "Parallel/GlobalCache.hpp"
      19             : #include "Time/Tags/StepNumberWithinSlab.hpp"
      20             : #include "Utilities/Gsl.hpp"
      21             : 
      22             : /// \cond
      23             : template <size_t Dim>
      24             : class Mesh;
      25             : namespace tuples {
      26             : template <typename... Tags>
      27             : class TaggedTuple;
      28             : }  // namespace tuples
      29             : /// \endcond
      30             : 
      31           0 : namespace dg::Actions {
      32             : /*!
      33             :  * \ingroup DiscontinuousGalerkinGroup
      34             :  * \brief Applies the element's spectral volume filter to the evolved variables.
      35             :  *
      36             :  * Retrieves the per-element filter from `Filters::Tags::SpectralFilter` from
      37             :  * the DataBox and applies the volume filter.
      38             :  *
      39             :  * The filter is skipped when:
      40             :  * - the filter is `Filters::None` (explicit no-op), or
      41             :  * - both `apply_volume_filter_on_substep()` and
      42             :  *   `apply_volume_filter_on_this_step(step_number)` return `false`.
      43             :  *
      44             :  * The grid-to-inertial Jacobian and its inverse are only retrieved from the
      45             :  * DataBox when `Filters::Filter::need_jacobians()` returns `true`; otherwise
      46             :  * `std::nullopt` is passed for both arguments.
      47             :  *
      48             :  * \note Currently the check for whether to filter on every `N` steps is done
      49             :  * relative to the start of the current Slab. This means that for GTS,
      50             :  * independent of the value of `N` for every `N` steps, every step has a
      51             :  * filter applied since GTS has one step per slab.
      52             :  *
      53             :  * Uses:
      54             :  * - DataBox:
      55             :  *   - `Filters::Tags::SpectralFilter<Dim, TagList>`
      56             :  *   - `Tags::StepNumberWithinSlab`
      57             :  *   - `domain::Tags::Mesh<Dim>`
      58             :  *   - `domain::Tags::InverseJacobian<Dim, Frame::Grid, Frame::Inertial>`
      59             :  *     (only when `need_jacobians()` is `true`)
      60             :  *   - `domain::Tags::Jacobian<Dim, Frame::Grid, Frame::Inertial>`
      61             :  *     (only when `need_jacobians()` is `true`)
      62             :  * - DataBox changes:
      63             :  *   - Modifies: `Metavariables::system::variables_tag`
      64             :  * - System:
      65             :  *   - `volume_dim`
      66             :  *   - `variables_tag`
      67             :  */
      68           1 : struct SpectralFilter {
      69             :   template <typename DbTags, typename... InboxTags, typename ArrayIndex,
      70             :             typename ActionList, typename ParallelComponent,
      71             :             typename Metavariables>
      72           0 :   static Parallel::iterable_action_return_t apply(
      73             :       db::DataBox<DbTags>& box,
      74             :       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
      75             :       const Parallel::GlobalCache<Metavariables>& /*cache*/,
      76             :       const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
      77             :       const ParallelComponent* const /*meta*/) {
      78             :     constexpr size_t Dim = Metavariables::system::volume_dim;
      79             :     using variables_tag = typename Metavariables::system::variables_tag;
      80             :     using TagList = typename variables_tag::tags_list;
      81             : 
      82             :     const auto& filter =
      83             :         db::get<Filters::Tags::SpectralFilter<Dim, TagList>>(box);
      84             : 
      85             :     const Mesh<Dim>& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
      86             :     if (not filter.supports_mesh(mesh)) {
      87             :       ERROR("The filter '"
      88             :             << filter.name() << "' on element "
      89             :             << db::get<domain::Tags::Element<Dim>>(box).id() << " with mesh "
      90             :             << mesh
      91             :             << " does not support that mesh. This is a bug in initialization "
      92             :                "or if AMR projection during h-refinement incorrectly adjusted "
      93             :                "the filter.");
      94             :     }
      95             : 
      96             :     // None is explicitly a no-op; skip immediately.
      97             :     if (dynamic_cast<const Filters::None<Dim, TagList>*>(&filter) != nullptr) {
      98             :       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
      99             :     }
     100             : 
     101             :     const auto step_number =
     102             :         static_cast<size_t>(db::get<::Tags::StepNumberWithinSlab>(box));
     103             :     if (not(filter.apply_volume_filter_on_substep() or
     104             :             filter.apply_volume_filter_on_this_step(step_number))) {
     105             :       return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     106             :     }
     107             : 
     108             :     const size_t num_pts = mesh.number_of_grid_points();
     109             : 
     110             :     std::optional<
     111             :         InverseJacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>>
     112             :         inv_jac{std::nullopt};
     113             :     std::optional<Jacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>> jac{
     114             :         std::nullopt};
     115             :     if (filter.need_jacobians()) {
     116             :       if constexpr (db::tag_is_retrievable_v<
     117             :                         domain::Tags::InverseJacobian<Dim, Frame::Grid,
     118             :                                                       Frame::Inertial>,
     119             :                         db::DataBox<DbTags>>) {
     120             :         const auto& src = db::get<
     121             :             domain::Tags::InverseJacobian<Dim, Frame::Grid, Frame::Inertial>>(
     122             :             box);
     123             :         inv_jac =
     124             :             InverseJacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>{};
     125             :         for (size_t i = 0; i < src.size(); ++i) {
     126             :           make_const_view(make_not_null(&std::as_const((*inv_jac)[i])), src[i],
     127             :                           0, num_pts);
     128             :         }
     129             :       } else {
     130             :         ERROR(
     131             :             "SpectralFilter: Missing "
     132             :             "domain::Tags::InverseJacobian<Dim, Grid, Inertial> in the "
     133             :             "DataBox. This should be present for both time-dependent and "
     134             :             "time-independent evolutions.");
     135             :       }
     136             :       if constexpr (db::tag_is_retrievable_v<
     137             :                         domain::Tags::Jacobian<Dim, Frame::Grid,
     138             :                                                Frame::Inertial>,
     139             :                         db::DataBox<DbTags>>) {
     140             :         const auto& src =
     141             :             db::get<domain::Tags::Jacobian<Dim, Frame::Grid, Frame::Inertial>>(
     142             :                 box);
     143             :         jac = Jacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>{};
     144             :         for (size_t i = 0; i < src.size(); ++i) {
     145             :           make_const_view(make_not_null(&std::as_const((*jac)[i])), src[i], 0,
     146             :                           num_pts);
     147             :         }
     148             :       } else {
     149             :         ERROR(
     150             :             "SpectralFilter: Missing "
     151             :             "domain::Tags::Jacobian<Dim, Grid, Inertial> in the "
     152             :             "DataBox. This should be present for both time-dependent and "
     153             :             "time-independent evolutions.");
     154             :       }
     155             :     }
     156             : 
     157             :     db::mutate<variables_tag>(
     158             :         [&filter, &mesh, &inv_jac,
     159             :          &jac](const gsl::not_null<typename variables_tag::type*> vars) {
     160             :           filter.apply_in_volume(vars, mesh, inv_jac, jac);
     161             :         },
     162             :         make_not_null(&box));
     163             : 
     164             :     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
     165             :   }
     166             : };
     167             : }  // namespace dg::Actions