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