Add ProjectSpectralFilter AMR projector

src/Evolution/DiscontinuousGalerkin/Initialization/CMakeLists.txt

@@ -6,6 +6,7 @@ spectre_target_headers(
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
   Mortars.hpp
+  ProjectSpectralFilters.hpp
   QuadratureTag.hpp
   SpectralFilters.hpp
   SpectralFilters.tpp

src/Evolution/DiscontinuousGalerkin/Initialization/ProjectSpectralFilters.hpp

@@ -0,0 +1,79 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+#include <iterator>
+#include <memory>
+#include <unordered_map>
+#include <utility>
+
+#include "DataStructures/TaggedTuple.hpp"
+#include "Domain/Structure/Element.hpp"
+#include "Domain/Structure/ElementId.hpp"
+#include "NumericalAlgorithms/LinearOperators/Filters/Filter.hpp"
+#include "NumericalAlgorithms/LinearOperators/Filters/Tag.hpp"
+#include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
+#include "Utilities/ErrorHandling/Error.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace evolution::dg::Initialization {
+
+/// \brief AMR projector for `Filters::Tags::SpectralFilter<Dim, TagList>`.
+///
+/// \details
+/// - For p-refinement: leaves the filter unchanged.
+/// - For h-refinement (splitting): clones the parent's filter for the child.
+/// - For h-coarsening (joining): clones the first child's filter for the
+///   parent (all siblings share the same block and mesh basis/quadrature, so
+///   they carry identical filters).
+template <size_t Dim, typename TagList>
+struct ProjectSpectralFilters : tt::ConformsTo<amr::protocols::Projector> {
+  using return_tags = tmpl::list<Filters::Tags::SpectralFilter<Dim, TagList>>;
+  using argument_tags = tmpl::list<>;
+
+  // p-refinement: leave the filter unchanged.
+  static void apply(
+      const gsl::not_null<
+          std::unique_ptr<Filters::Filter<Dim, TagList>>*> /*filter*/,
+      const std::pair<Mesh<Dim>, Element<Dim>>& /*old_mesh_and_element*/) {}
+
+  // Splitting: clone the parent's filter for the new child.
+  template <typename... ParentTags>
+  static void apply(
+      const gsl::not_null<std::unique_ptr<Filters::Filter<Dim, TagList>>*>
+          filter,
+      const tuples::TaggedTuple<ParentTags...>& parent_items) {
+    *filter =
+        tuples::get<Filters::Tags::SpectralFilter<Dim, TagList>>(parent_items)
+            ->get_clone();
+  }
+
+  // Joining: clone the first child's filter for the parent.
+  template <typename... ChildrenTags>
+  static void apply(
+      const gsl::not_null<std::unique_ptr<Filters::Filter<Dim, TagList>>*>
+          filter,
+      const std::unordered_map<ElementId<Dim>,
+                               tuples::TaggedTuple<ChildrenTags...>>&
+          children_items) {
+    const auto first_child_items = children_items.begin();
+    const auto& first_filter =
+        tuples::get<Filters::Tags::SpectralFilter<Dim, TagList>>(
+            first_child_items->second);
+    for (auto it = std::next(first_child_items); it != children_items.end();
+         ++it) {
+      const auto& other_filter =
+          tuples::get<Filters::Tags::SpectralFilter<Dim, TagList>>(it->second);
+      if (not first_filter->is_equal(*other_filter)) {
+        ERROR("Children do not agree on all items!");
+      }
+    }
+    *filter = first_filter->get_clone();
+  }
+};
+
+}  // namespace evolution::dg::Initialization

src/NumericalAlgorithms/LinearOperators/Filters/Filter.hpp

@@ -87,6 +87,11 @@ class Filter : public PUP::able {
   virtual bool apply_volume_filter_on_substep() const = 0;
 
   /// Whether the volume filter should be applied at the given `step_number`.
+  ///
+  /// \note Currently the check for whether to filter on every `N` steps is done
+  /// relative to the start of the current Slab. This means that for GTS,
+  /// independent of the value of `N` for every `N` steps, every step has a
+  /// filter applied since GTS has one step per slab.
   virtual bool apply_volume_filter_on_this_step(size_t step_number) const = 0;
 
   /// Whether the boundary filter should be applied inside every Runge-Kutta
@@ -95,6 +100,11 @@ class Filter : public PUP::able {
 
   /// Whether the boundary filter should be applied at the given
   /// `step_number`.
+  ///
+  /// \note Currently the check for whether to filter on every `N` steps is done
+  /// relative to the start of the current Slab. This means that for GTS,
+  /// independent of the value of `N` for every `N` steps, every step has a
+  /// filter applied since GTS has one step per slab.
   virtual bool apply_boundary_filter_on_this_step(size_t step_number) const = 0;
 
   /// \brief Whether the filter needs the grid-to-inertial Jacobian and its
@@ -108,6 +118,10 @@ class Filter : public PUP::able {
   /// \brief Returns `true` if this filter can filter the `mesh`.
   virtual bool supports_mesh(const Mesh<Dim>& mesh) const = 0;
 
+  /// \brief A human-readable name for the concrete filter type, used in
+  /// diagnostics.
+  virtual std::string name() const = 0;
+
   /// \brief Returns `true` if `other` is the same concrete filter type and is
   /// equivalent to `*this`.
   virtual bool is_equal(const Filter& other) const = 0;

src/NumericalAlgorithms/LinearOperators/Filters/Hypercube.hpp

@@ -115,6 +115,11 @@ class Hypercube : public Filter<Dim, TagList> {
 
   /// \brief Apply the volume filter once every `N` steps. `None`
   /// (`std::nullopt`) disables the every-N-steps trigger.
+  ///
+  /// \note Currently the check for whether to filter on every `N` steps is done
+  /// relative to the start of the current Slab. This means that for GTS,
+  /// independent of the value of `N` for every `N` steps, every step has a
+  /// filter applied since GTS has one step per slab.
   struct VolumeFilterEveryNSteps {
     using type = Options::Auto<size_t, Options::AutoLabel::None>;
     static constexpr Options::String help = {
@@ -123,6 +128,11 @@ class Hypercube : public Filter<Dim, TagList> {
 
   /// \brief Apply the boundary correction filter once every `N` steps. `None`
   /// (`std::nullopt`) disables the every-N-steps trigger.
+  ///
+  /// \note Currently the check for whether to filter on every `N` steps is done
+  /// relative to the start of the current Slab. This means that for GTS,
+  /// independent of the value of `N` for every `N` steps, every step has a
+  /// filter applied since GTS has one step per slab.
   struct BoundaryCorrectionFilterEveryNSteps {
     using type = Options::Auto<size_t, Options::AutoLabel::None>;
     static constexpr Options::String help = {
@@ -166,6 +176,8 @@ class Hypercube : public Filter<Dim, TagList> {
 
   bool supports_mesh(const Mesh<Dim>& mesh) const override;
 
+  std::string name() const override { return "Hypercube"; }
+
   const std::optional<std::vector<size_t>>& blocks_to_filter() const override;
 
   void set_blocks_to_filter(

src/NumericalAlgorithms/LinearOperators/Filters/None.hpp

@@ -85,6 +85,8 @@ class None : public Filter<Dim, TagList> {
 
   bool supports_mesh(const Mesh<Dim>& /*mesh*/) const override { return true; }
 
+  std::string name() const override { return "None"; }
+
   const std::optional<std::vector<size_t>>& blocks_to_filter() const override;
 
   void set_blocks_to_filter(

src/NumericalAlgorithms/LinearOperators/Filters/SphericalShell.hpp

@@ -158,6 +158,11 @@ class SphericalShell : public Filter<3, TagList> {
 
   /// \brief Apply the volume filter once every `N` steps. `None`
   /// (`std::nullopt`) disables the every-N-steps trigger.
+  ///
+  /// \note Currently the check for whether to filter on every `N` steps is done
+  /// relative to the start of the current Slab. This means that for GTS,
+  /// independent of the value of `N` for every `N` steps, every step has a
+  /// filter applied since GTS has one step per slab.
   struct VolumeFilterEveryNSteps {
     using type = Options::Auto<size_t, Options::AutoLabel::None>;
     static constexpr Options::String help = {
@@ -166,6 +171,11 @@ class SphericalShell : public Filter<3, TagList> {
 
   /// \brief Apply the boundary correction filter once every `N` steps. `None`
   /// (`std::nullopt`) disables the every-N-steps trigger.
+  ///
+  /// \note Currently the check for whether to filter on every `N` steps is done
+  /// relative to the start of the current Slab. This means that for GTS,
+  /// independent of the value of `N` for every `N` steps, every step has a
+  /// filter applied since GTS has one step per slab.
   struct BoundaryCorrectionFilterEveryNSteps {
     using type = Options::Auto<size_t, Options::AutoLabel::None>;
     static constexpr Options::String help = {
@@ -213,6 +223,8 @@ class SphericalShell : public Filter<3, TagList> {
 
   bool supports_mesh(const Mesh<3>& mesh) const override;
 
+  std::string name() const override { return "SphericalShell"; }
+
   const std::optional<std::vector<size_t>>& blocks_to_filter() const override;
 
   void set_blocks_to_filter(

src/ParallelAlgorithms/Actions/CMakeLists.txt

@@ -21,6 +21,7 @@ spectre_target_headers(
   MutateApply.hpp
   RandomizeVariables.hpp
   SetData.hpp
+  SpectralFilter.hpp
   TerminatePhase.hpp
   UpdateMessageQueue.hpp
   )
@@ -34,9 +35,11 @@ target_link_libraries(
   DomainStructure
   ErrorHandling
   FunctionsOfTime
+  LinearOperators
   Parallel
   Serialization
   Spectral
+  Time
   Utilities
   )
 

src/ParallelAlgorithms/Actions/SpectralFilter.hpp

@@ -0,0 +1,167 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include <cstddef>
+#include <optional>
+#include <utility>
+
+#include "DataStructures/DataBox/DataBox.hpp"
+#include "DataStructures/DataVector.hpp"
+#include "DataStructures/Tensor/Tensor.hpp"
+#include "Domain/Tags.hpp"
+#include "NumericalAlgorithms/LinearOperators/Filters/Filter.hpp"
+#include "NumericalAlgorithms/LinearOperators/Filters/None.hpp"
+#include "NumericalAlgorithms/LinearOperators/Filters/Tag.hpp"
+#include "Parallel/AlgorithmExecution.hpp"
+#include "Parallel/GlobalCache.hpp"
+#include "Time/Tags/StepNumberWithinSlab.hpp"
+#include "Utilities/Gsl.hpp"
+
+/// \cond
+template <size_t Dim>
+class Mesh;
+namespace tuples {
+template <typename... Tags>
+class TaggedTuple;
+}  // namespace tuples
+/// \endcond
+
+namespace dg::Actions {
+/*!
+ * \ingroup DiscontinuousGalerkinGroup
+ * \brief Applies the element's spectral volume filter to the evolved variables.
+ *
+ * Retrieves the per-element filter from `Filters::Tags::SpectralFilter` from
+ * the DataBox and applies the volume filter.
+ *
+ * The filter is skipped when:
+ * - the filter is `Filters::None` (explicit no-op), or
+ * - both `apply_volume_filter_on_substep()` and
+ *   `apply_volume_filter_on_this_step(step_number)` return `false`.
+ *
+ * The grid-to-inertial Jacobian and its inverse are only retrieved from the
+ * DataBox when `Filters::Filter::need_jacobians()` returns `true`; otherwise
+ * `std::nullopt` is passed for both arguments.
+ *
+ * \note Currently the check for whether to filter on every `N` steps is done
+ * relative to the start of the current Slab. This means that for GTS,
+ * independent of the value of `N` for every `N` steps, every step has a
+ * filter applied since GTS has one step per slab.
+ *
+ * Uses:
+ * - DataBox:
+ *   - `Filters::Tags::SpectralFilter<Dim, TagList>`
+ *   - `Tags::StepNumberWithinSlab`
+ *   - `domain::Tags::Mesh<Dim>`
+ *   - `domain::Tags::InverseJacobian<Dim, Frame::Grid, Frame::Inertial>`
+ *     (only when `need_jacobians()` is `true`)
+ *   - `domain::Tags::Jacobian<Dim, Frame::Grid, Frame::Inertial>`
+ *     (only when `need_jacobians()` is `true`)
+ * - DataBox changes:
+ *   - Modifies: `Metavariables::system::variables_tag`
+ * - System:
+ *   - `volume_dim`
+ *   - `variables_tag`
+ */
+struct SpectralFilter {
+  template <typename DbTags, typename... InboxTags, typename ArrayIndex,
+            typename ActionList, typename ParallelComponent,
+            typename Metavariables>
+  static Parallel::iterable_action_return_t apply(
+      db::DataBox<DbTags>& box,
+      const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
+      const Parallel::GlobalCache<Metavariables>& /*cache*/,
+      const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
+      const ParallelComponent* const /*meta*/) {
+    constexpr size_t Dim = Metavariables::system::volume_dim;
+    using variables_tag = typename Metavariables::system::variables_tag;
+    using TagList = typename variables_tag::tags_list;
+
+    const auto& filter =
+        db::get<Filters::Tags::SpectralFilter<Dim, TagList>>(box);
+
+    const Mesh<Dim>& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
+    if (not filter.supports_mesh(mesh)) {
+      ERROR("The filter '"
+            << filter.name() << "' on element "
+            << db::get<domain::Tags::Element<Dim>>(box).id() << " with mesh "
+            << mesh
+            << " does not support that mesh. This is a bug in initialization "
+               "or if AMR projection during h-refinement incorrectly adjusted "
+               "the filter.");
+    }
+
+    // None is explicitly a no-op; skip immediately.
+    if (dynamic_cast<const Filters::None<Dim, TagList>*>(&filter) != nullptr) {
+      return {Parallel::AlgorithmExecution::Continue, std::nullopt};
+    }
+
+    const auto step_number =
+        static_cast<size_t>(db::get<::Tags::StepNumberWithinSlab>(box));
+    if (not(filter.apply_volume_filter_on_substep() or
+            filter.apply_volume_filter_on_this_step(step_number))) {
+      return {Parallel::AlgorithmExecution::Continue, std::nullopt};
+    }
+
+    const size_t num_pts = mesh.number_of_grid_points();
+
+    std::optional<
+        InverseJacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>>
+        inv_jac{std::nullopt};
+    std::optional<Jacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>> jac{
+        std::nullopt};
+    if (filter.need_jacobians()) {
+      if constexpr (db::tag_is_retrievable_v<
+                        domain::Tags::InverseJacobian<Dim, Frame::Grid,
+                                                      Frame::Inertial>,
+                        db::DataBox<DbTags>>) {
+        const auto& src = db::get<
+            domain::Tags::InverseJacobian<Dim, Frame::Grid, Frame::Inertial>>(
+            box);
+        inv_jac =
+            InverseJacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>{};
+        for (size_t i = 0; i < src.size(); ++i) {
+          make_const_view(make_not_null(&std::as_const((*inv_jac)[i])), src[i],
+                          0, num_pts);
+        }
+      } else {
+        ERROR(
+            "SpectralFilter: Missing "
+            "domain::Tags::InverseJacobian<Dim, Grid, Inertial> in the "
+            "DataBox. This should be present for both time-dependent and "
+            "time-independent evolutions.");
+      }
+      if constexpr (db::tag_is_retrievable_v<
+                        domain::Tags::Jacobian<Dim, Frame::Grid,
+                                               Frame::Inertial>,
+                        db::DataBox<DbTags>>) {
+        const auto& src =
+            db::get<domain::Tags::Jacobian<Dim, Frame::Grid, Frame::Inertial>>(
+                box);
+        jac = Jacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>{};
+        for (size_t i = 0; i < src.size(); ++i) {
+          make_const_view(make_not_null(&std::as_const((*jac)[i])), src[i], 0,
+                          num_pts);
+        }
+      } else {
+        ERROR(
+            "SpectralFilter: Missing "
+            "domain::Tags::Jacobian<Dim, Grid, Inertial> in the "
+            "DataBox. This should be present for both time-dependent and "
+            "time-independent evolutions.");
+      }
+    }
+
+    db::mutate<variables_tag>(
+        [&filter, &mesh, &inv_jac,
+         &jac](const gsl::not_null<typename variables_tag::type*> vars) {
+          filter.apply_in_volume(vars, mesh, inv_jac, jac);
+        },
+        make_not_null(&box));
+
+    return {Parallel::AlgorithmExecution::Continue, std::nullopt};
+  }
+};
+}  // namespace dg::Actions

tests/Unit/Evolution/DiscontinuousGalerkin/CMakeLists.txt

@@ -12,6 +12,7 @@ set(LIBRARY_SOURCES
   Actions/Test_ComputeTimeDerivative.cpp
   Actions/Test_NormalCovectorAndMagnitude.cpp
   Initialization/Test_Mortars.cpp
+  Initialization/Test_ProjectSpectralFilters.cpp
   Initialization/Test_QuadratureTag.cpp
   Initialization/Test_SpectralFilters.cpp
   Test_AtomicInboxBoundaryData.cpp