Merge pull request #6611 from knelli2/intrp_specific_blocks

Only interpolate for horizon finds from certain blocks

Author: nilsdeppe

src/Domain/Creators/Disk.cpp

@@ -60,6 +60,17 @@ Disk::Disk(typename InnerRadius::type inner_radius,
   if (boundary_condition_ != nullptr and is_periodic(boundary_condition_)) {
     PARSE_ERROR(context, "Cannot have periodic boundary conditions on a disk.");
   }
+
+  block_groups_["Shell0"];
+  block_names_.reserve(5);
+  const std::array<std::string, 4> wedge_directions{"UpperX", "UpperY",
+                                                    "LowerX", "LowerY"};
+  for (const std::string& name : wedge_directions) {
+    block_names_.emplace_back(name);
+    block_groups_["Shell0"].insert(name);
+  }
+  block_names_.emplace_back("CenterSquare");
+  block_groups_["CenterSquare"].insert("CenterSquare");
 }
 
 Domain<2> Disk::create_domain() const {
@@ -131,7 +142,8 @@ Domain<2> Disk::create_domain() const {
     boundary_conditions_all_blocks.emplace_back();
   }
 
-  return Domain<2>{std::move(coord_maps), corners};
+  return Domain<2>{std::move(coord_maps), corners,      {}, {},
+                   block_names_,          block_groups_};
 }
 
 std::vector<DirectionMap<

src/Domain/Creators/Disk.hpp

@@ -6,6 +6,9 @@
 #include <array>
 #include <cstddef>
 #include <memory>
+#include <string>
+#include <unordered_map>
+#include <unordered_set>
 #include <vector>
 
 #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
@@ -138,6 +141,13 @@ class Disk : public DomainCreator<2> {
 
   std::vector<std::array<size_t, 2>> initial_refinement_levels() const override;
 
+  std::vector<std::string> block_names() const override { return block_names_; }
+
+  std::unordered_map<std::string, std::unordered_set<std::string>>
+  block_groups() const override {
+    return block_groups_;
+  }
+
  private:
   typename InnerRadius::type inner_radius_{};
   typename OuterRadius::type outer_radius_{};
@@ -146,6 +156,9 @@ class Disk : public DomainCreator<2> {
   typename UseEquiangularMap::type use_equiangular_map_{false};
   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
       boundary_condition_;
+  std::vector<std::string> block_names_;
+  std::unordered_map<std::string, std::unordered_set<std::string>>
+      block_groups_;
 };
 }  // namespace creators
 }  // namespace domain

src/Domain/Creators/Rectilinear.cpp

@@ -180,7 +180,8 @@ Domain<Dim> Rectilinear<Dim>::create_domain() const {
       {std::move(block_corners)},
       std::move(identifications),
       {},
-      block_names_};
+      block_names(),
+      block_groups()};
 
   if (not time_dependence_->is_none()) {
     domain.inject_time_dependent_map_for_block(

src/Domain/Creators/Rectilinear.hpp

@@ -7,6 +7,7 @@
 #include <cstddef>
 #include <memory>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
@@ -247,6 +248,11 @@ class Rectilinear : public DomainCreator<Dim> {
 
   std::vector<std::string> block_names() const override { return block_names_; }
 
+  std::unordered_map<std::string, std::unordered_set<std::string>>
+  block_groups() const override {
+    return {{name(), {name()}}};
+  }
+
   // Transforms from option-created boundary conditions to the type used in the
   // constructor
   template <typename BoundaryConditionsBase>

src/Evolution/Executables/GrMhd/GhValenciaDivClean/EvolveGhValenciaDivClean.hpp

@@ -24,10 +24,10 @@ template <bool UseControlSystems, typename... InterpolationTargetTags>
 struct EvolutionMetavars
     : public GhValenciaDivCleanTemplateBase<
           EvolutionMetavars<UseControlSystems, InterpolationTargetTags...>,
-          true, UseControlSystems> {
+          true, UseControlSystems, false> {
   using base = GhValenciaDivCleanTemplateBase<
       EvolutionMetavars<UseControlSystems, InterpolationTargetTags...>, true,
-      UseControlSystems>;
+      UseControlSystems, false>;
   using const_global_cache_tags = typename base::const_global_cache_tags;
   using observed_reduction_data_tags =
       typename base::observed_reduction_data_tags;

src/Evolution/Executables/GrMhd/GhValenciaDivClean/EvolveGhValenciaDivCleanWithHorizon.hpp

@@ -43,14 +43,14 @@ template <bool UseControlSystems, typename... InterpolationTargetTags>
 struct EvolutionMetavars
     : public GhValenciaDivCleanTemplateBase<
           EvolutionMetavars<UseControlSystems, InterpolationTargetTags...>,
-          false, false> {
+          false, false, true> {
   static_assert(not UseControlSystems,
                 "GhValenciaWithHorizon doesn't support control systems yet.");
   static constexpr bool use_dg_subcell = false;
 
   using defaults = GhValenciaDivCleanDefaults<use_dg_subcell>;
   using base = GhValenciaDivCleanTemplateBase<EvolutionMetavars, use_dg_subcell,
-                                              UseControlSystems>;
+                                              UseControlSystems, true>;
   static constexpr size_t volume_dim = defaults::volume_dim;
   using domain_frame = typename defaults::domain_frame;
   static constexpr bool use_damped_harmonic_rollon =

src/Evolution/Executables/GrMhd/GhValenciaDivClean/GhValenciaDivCleanBase.hpp

@@ -337,15 +337,15 @@ struct GhValenciaDivCleanDefaults {
 };
 
 template <typename EvolutionMetavarsDerived, bool UseDgSubcell,
-          bool UseControlSystems>
+          bool UseControlSystems, bool WithHorizon>
 struct GhValenciaDivCleanTemplateBase;
 
-template <bool UseDgSubcell, bool UseControlSystems,
+template <bool UseDgSubcell, bool UseControlSystems, bool WithHorizon,
           template <bool, typename...> class EvolutionMetavarsDerived,
           typename... InterpolationTargetTags>
 struct GhValenciaDivCleanTemplateBase<
     EvolutionMetavarsDerived<UseControlSystems, InterpolationTargetTags...>,
-    UseDgSubcell, UseControlSystems>
+    UseDgSubcell, UseControlSystems, WithHorizon>
     : public virtual GhValenciaDivCleanDefaults<UseDgSubcell> {
   using derived_metavars =
       EvolutionMetavarsDerived<UseControlSystems, InterpolationTargetTags...>;
@@ -690,9 +690,11 @@ struct GhValenciaDivCleanTemplateBase<
       gh::ConstraintDamping::Tags::DampingFunctionGamma2<volume_dim,
                                                          Frame::Grid>>>;
 
-  using dg_registration_list =
-      tmpl::list<intrp::Actions::RegisterElementWithInterpolator,
-                 observers::Actions::RegisterEventsWithObservers>;
+  using dg_registration_list = tmpl::flatten<tmpl::list<
+      tmpl::conditional_t<WithHorizon,
+                          intrp::Actions::RegisterElementWithInterpolator,
+                          tmpl::list<>>,
+      observers::Actions::RegisterEventsWithObservers>>;
 
   static constexpr auto default_phase_order = std::array<Parallel::Phase, 8>{
       {Parallel::Phase::Initialization,
@@ -928,7 +930,8 @@ struct GhValenciaDivCleanTemplateBase<
       observers::ObserverWriter<derived_metavars>,
       importers::ElementDataReader<derived_metavars>,
       control_system::control_components<derived_metavars, control_systems>,
-      intrp::Interpolator<derived_metavars>,
+      tmpl::conditional_t<WithHorizon, intrp::Interpolator<derived_metavars>,
+                          tmpl::list<>>,
       intrp::InterpolationTarget<derived_metavars, InterpolationTargetTags>...,
       dg_element_array_component>>;
 };

src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp

@@ -118,19 +118,12 @@
 #include "ParallelAlgorithms/EventsAndTriggers/EventsAndTriggers.hpp"
 #include "ParallelAlgorithms/EventsAndTriggers/LogicalTriggers.hpp"
 #include "ParallelAlgorithms/EventsAndTriggers/Trigger.hpp"
-#include "ParallelAlgorithms/Interpolation/Actions/CleanUpInterpolator.hpp"
 #include "ParallelAlgorithms/Interpolation/Actions/ElementInitInterpPoints.hpp"
 #include "ParallelAlgorithms/Interpolation/Actions/InitializeInterpolationTarget.hpp"
-#include "ParallelAlgorithms/Interpolation/Actions/InterpolationTargetReceiveVars.hpp"
-#include "ParallelAlgorithms/Interpolation/Actions/InterpolatorReceivePoints.hpp"
-#include "ParallelAlgorithms/Interpolation/Actions/InterpolatorReceiveVolumeData.hpp"
-#include "ParallelAlgorithms/Interpolation/Actions/InterpolatorRegisterElement.hpp"
-#include "ParallelAlgorithms/Interpolation/Actions/TryToInterpolate.hpp"
 #include "ParallelAlgorithms/Interpolation/Callbacks/ObserveTimeSeriesOnSurface.hpp"
 #include "ParallelAlgorithms/Interpolation/Events/Interpolate.hpp"
 #include "ParallelAlgorithms/Interpolation/Events/InterpolateWithoutInterpComponent.hpp"
 #include "ParallelAlgorithms/Interpolation/InterpolationTarget.hpp"
-#include "ParallelAlgorithms/Interpolation/Interpolator.hpp"
 #include "ParallelAlgorithms/Interpolation/Protocols/InterpolationTargetTag.hpp"
 #include "ParallelAlgorithms/Interpolation/Tags.hpp"
 #include "ParallelAlgorithms/Interpolation/Targets/KerrHorizon.hpp"
@@ -515,8 +508,7 @@ struct EvolutionMetavars<tmpl::list<InterpolationTargetTags...>,
                           dg_step_actions>;
 
   using dg_registration_list =
-      tmpl::list<intrp::Actions::RegisterElementWithInterpolator,
-                 observers::Actions::RegisterEventsWithObservers>;
+      tmpl::list<observers::Actions::RegisterEventsWithObservers>;
 
   using initialization_actions = tmpl::flatten<tmpl::list<
       Initialization::Actions::InitializeItems<
@@ -613,7 +605,6 @@ struct EvolutionMetavars<tmpl::list<InterpolationTargetTags...>,
   using component_list = tmpl::list<
       observers::Observer<EvolutionMetavars>,
       observers::ObserverWriter<EvolutionMetavars>,
-      intrp::Interpolator<EvolutionMetavars>,
       intrp::InterpolationTarget<EvolutionMetavars, InterpolationTargetTags>...,
       dg_element_array_component>;
 

src/ParallelAlgorithms/ApparentHorizonFinder/InterpolationTarget.cpp

@@ -5,8 +5,10 @@
 
 #include <algorithm>
 #include <cstddef>
+#include <unordered_set>
 
 #include "IO/Logging/Verbosity.hpp"
+#include "Options/Context.hpp"
 #include "Utilities/GenerateInstantiations.hpp"
 
 namespace Frame {
@@ -17,11 +19,13 @@ namespace intrp::OptionHolders {
 template <typename Frame>
 ApparentHorizon<Frame>::ApparentHorizon(
     ylm::Strahlkorper<Frame> initial_guess_in, ::FastFlow fast_flow_in,
-    ::Verbosity verbosity_in, const size_t max_interpolation_retries_in)
+    ::Verbosity verbosity_in, const size_t max_interpolation_retries_in,
+    std::optional<std::vector<std::string>> blocks_for_interpolation_in)
     : initial_guess(std::move(initial_guess_in)),
       fast_flow(std::move(fast_flow_in)),  // NOLINT
       verbosity(std::move(verbosity_in)),  // NOLINT
-      max_interpolation_retries(max_interpolation_retries_in) {}
+      max_interpolation_retries(max_interpolation_retries_in),
+      blocks_for_interpolation(std::move(blocks_for_interpolation_in)) {}
 // clang-tidy std::move of trivially copyable type.
 
 template <typename Frame>
@@ -30,14 +34,16 @@ void ApparentHorizon<Frame>::pup(PUP::er& p) {
   p | fast_flow;
   p | verbosity;
   p | max_interpolation_retries;
+  p | blocks_for_interpolation;
 }
 
 template <typename Frame>
 bool operator==(const ApparentHorizon<Frame>& lhs,
                 const ApparentHorizon<Frame>& rhs) {
   return lhs.initial_guess == rhs.initial_guess and
          lhs.fast_flow == rhs.fast_flow and lhs.verbosity == rhs.verbosity and
-         lhs.max_interpolation_retries == rhs.max_interpolation_retries;
+         lhs.max_interpolation_retries == rhs.max_interpolation_retries and
+         lhs.blocks_for_interpolation == rhs.blocks_for_interpolation;
 }
 
 template <typename Frame>