Add some elliptic AMR projectors

.clang-tidy

@@ -78,6 +78,8 @@ Checks: '*,
 -readability-magic-numbers,
 # Same as llvm-qualified-auto above.,
 -readability-qualified-auto,
+# Access specifiers can be useful to structure code,
+-readability-redundant-access-specifiers,
 # We can have two of: methods are static when possible, static,
 # methods are not called through instances, and methods of,
 # calling, e.g., x.size(), are consistent across classes.  We,

src/Elliptic/Actions/InitializeAnalyticSolution.hpp

@@ -12,14 +12,16 @@
 #include "DataStructures/DataBox/PrefixHelpers.hpp"
 #include "DataStructures/DataBox/Prefixes.hpp"
 #include "DataStructures/Tensor/Tensor.hpp"
-#include "DataStructures/VariablesTag.hpp"
+#include "DataStructures/Variables.hpp"
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Tags.hpp"
-#include "Elliptic/Utilities/GetAnalyticData.hpp"
 #include "Parallel/AlgorithmExecution.hpp"
 #include "Parallel/GlobalCache.hpp"
+#include "Parallel/Tags/Metavariables.hpp"
+#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
+#include "Utilities/CallWithDynamicType.hpp"
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
 
@@ -49,37 +51,69 @@ namespace elliptic::Actions {
  * - Adds:
  *   - `::Tags::AnalyticSolutionsBase`
  */
-template <typename BackgroundTag, typename AnalyticSolutionFields,
+template <size_t Dim, typename BackgroundTag, typename AnalyticSolutionFields,
           typename AnalyticSolutionType>
-struct InitializeOptionalAnalyticSolution {
+struct InitializeOptionalAnalyticSolution
+    : tt::ConformsTo<::amr::protocols::Projector> {
  private:
   using analytic_fields_tag = ::Tags::AnalyticSolutions<AnalyticSolutionFields>;
 
- public:
+ public:  // Iterable action
   using simple_tags = tmpl::list<analytic_fields_tag>;
   using compute_tags = tmpl::list<>;
   using const_global_cache_tags = tmpl::list<BackgroundTag>;
 
   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
-            size_t Dim, typename ActionList, typename ParallelComponent>
+            typename ActionList, typename ParallelComponent>
   static Parallel::iterable_action_return_t apply(
       db::DataBox<DbTagsList>& box,
       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
       const Parallel::GlobalCache<Metavariables>& /*cache*/,
       const ElementId<Dim>& /*array_index*/, const ActionList /*meta*/,
       const ParallelComponent* const /*meta*/) {
+    db::mutate_apply<InitializeOptionalAnalyticSolution>(make_not_null(&box));
+    return {Parallel::AlgorithmExecution::Continue, std::nullopt};
+  }
+
+ public:  // DataBox mutator, amr::protocols::Projector
+  using return_tags = tmpl::list<analytic_fields_tag>;
+  using argument_tags =
+      tmpl::list<domain::Tags::Mesh<Dim>,
+                 domain::Tags::Coordinates<Dim, Frame::Inertial>, BackgroundTag,
+                 Parallel::Tags::Metavariables>;
+
+  template <typename Background, typename Metavariables, typename... AmrData>
+  static void apply(const gsl::not_null<typename analytic_fields_tag::type*>
+                        analytic_solution_fields,
+                    const Mesh<Dim>& mesh,
+                    const tnsr::I<DataVector, Dim> inertial_coords,
+                    const Background& background, const Metavariables& /*meta*/,
+                    const AmrData&... amr_data) {
+    if constexpr (sizeof...(AmrData) == 1) {
+      if constexpr (std::is_same_v<AmrData...,
+                                   std::pair<Mesh<Dim>, Element<Dim>>>) {
+        if (((mesh == amr_data.first) and ...)) {
+          // This element hasn't changed during AMR. Nothing to do.
+          return;
+        }
+      }
+    }
+
     const auto analytic_solution =
-        dynamic_cast<const AnalyticSolutionType*>(&db::get<BackgroundTag>(box));
+        dynamic_cast<const AnalyticSolutionType*>(&background);
     if (analytic_solution != nullptr) {
-      const auto& inertial_coords =
-          get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box);
-      auto analytic_fields =
-          elliptic::util::get_analytic_data<AnalyticSolutionFields>(
-              *analytic_solution, box, inertial_coords);
-      Initialization::mutate_assign<simple_tags>(make_not_null(&box),
-                                                 std::move(analytic_fields));
+      using factory_classes = typename std::decay_t<
+          Metavariables>::factory_creation::factory_classes;
+      *analytic_solution_fields = call_with_dynamic_type<
+          Variables<AnalyticSolutionFields>,
+          tmpl::at<factory_classes, AnalyticSolutionType>>(
+          analytic_solution, [&inertial_coords](const auto* const derived) {
+            return variables_from_tagged_tuple(
+                derived->variables(inertial_coords, AnalyticSolutionFields{}));
+          });
+    } else {
+      *analytic_solution_fields = std::nullopt;
     }
-    return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
 };
 /// @}

src/Elliptic/Actions/InitializeFixedSources.hpp

@@ -14,10 +14,12 @@
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Tags.hpp"
 #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
-#include "Elliptic/Utilities/GetAnalyticData.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp"
 #include "Parallel/AlgorithmExecution.hpp"
+#include "Parallel/Tags/Metavariables.hpp"
+#include "ParallelAlgorithms/Amr/Protocols/Projector.hpp"
 #include "ParallelAlgorithms/Initialization/MutateAssign.hpp"
+#include "Utilities/CallWithDynamicType.hpp"
 #include "Utilities/MakeWithValue.hpp"
 #include "Utilities/TMPL.hpp"
 #include "Utilities/TaggedTuple.hpp"
@@ -49,50 +51,64 @@ namespace elliptic::Actions {
  * - Adds:
  *   - `db::wrap_tags_in<::Tags::FixedSource, primal_fields>`
  */
-template <typename System, typename BackgroundTag>
-struct InitializeFixedSources {
+template <typename System, typename BackgroundTag,
+          size_t Dim = System::volume_dim>
+struct InitializeFixedSources : tt::ConformsTo<::amr::protocols::Projector> {
  private:
   using fixed_sources_tag = ::Tags::Variables<
       db::wrap_tags_in<::Tags::FixedSource, typename System::primal_fields>>;
 
- public:
+ public:  // Iterable action
   using const_global_cache_tags =
       tmpl::list<elliptic::dg::Tags::Massive, BackgroundTag>;
   using simple_tags = tmpl::list<fixed_sources_tag>;
   using compute_tags = tmpl::list<>;
 
   template <typename DbTagsList, typename... InboxTags, typename Metavariables,
-            size_t Dim, typename ActionList, typename ParallelComponent>
+            typename ActionList, typename ParallelComponent>
   static Parallel::iterable_action_return_t apply(
       db::DataBox<DbTagsList>& box,
       const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
       const Parallel::GlobalCache<Metavariables>& /*cache*/,
       const ElementId<Dim>& /*array_index*/, const ActionList /*meta*/,
       const ParallelComponent* const /*meta*/) {
-    const auto& inertial_coords =
-        get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box);
-    const auto& background = db::get<BackgroundTag>(box);
+    db::mutate_apply<InitializeFixedSources>(make_not_null(&box));
+    return {Parallel::AlgorithmExecution::Continue, std::nullopt};
+  }
+
+ public:  // DataBox mutator, amr::protocols::Projector
+  using return_tags = tmpl::list<fixed_sources_tag>;
+  using argument_tags = tmpl::list<
+      domain::Tags::Coordinates<Dim, Frame::Inertial>, BackgroundTag,
+      elliptic::dg::Tags::Massive, domain::Tags::Mesh<Dim>,
+      domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>,
+      Parallel::Tags::Metavariables>;
 
+  template <typename Background, typename Metavariables, typename... AmrData>
+  static void apply(
+      const gsl::not_null<typename fixed_sources_tag::type*> fixed_sources,
+      const tnsr::I<DataVector, Dim>& inertial_coords,
+      const Background& background, const bool massive, const Mesh<Dim>& mesh,
+      const Scalar<DataVector>& det_inv_jacobian, const Metavariables& /*meta*/,
+      const AmrData&... /*amr_data*/) {
     // Retrieve the fixed-sources of the elliptic system from the background,
     // which (along with the boundary conditions) define the problem we want to
     // solve.
-    auto fixed_sources = elliptic::util::get_analytic_data<
-        typename fixed_sources_tag::tags_list>(background, box,
-                                               inertial_coords);
+    using factory_classes =
+        typename std::decay_t<Metavariables>::factory_creation::factory_classes;
+    *fixed_sources =
+        call_with_dynamic_type<Variables<typename fixed_sources_tag::tags_list>,
+                               tmpl::at<factory_classes, Background>>(
+            &background, [&inertial_coords](const auto* const derived) {
+              return variables_from_tagged_tuple(derived->variables(
+                  inertial_coords, typename fixed_sources_tag::tags_list{}));
+            });
 
     // Apply DG mass matrix to the fixed sources if the DG operator is massive
-    if (db::get<elliptic::dg::Tags::Massive>(box)) {
-      const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
-      const auto& det_inv_jacobian = db::get<
-          domain::Tags::DetInvJacobian<Frame::ElementLogical, Frame::Inertial>>(
-          box);
-      fixed_sources /= get(det_inv_jacobian);
-      ::dg::apply_mass_matrix(make_not_null(&fixed_sources), mesh);
+    if (massive) {
+      *fixed_sources /= get(det_inv_jacobian);
+      ::dg::apply_mass_matrix(fixed_sources, mesh);
     }
-
-    ::Initialization::mutate_assign<simple_tags>(make_not_null(&box),
-                                                 std::move(fixed_sources));
-    return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
 };
 

src/Elliptic/CMakeLists.txt

@@ -23,6 +23,7 @@ target_link_libraries(
   Domain
   EventsAndTriggers
   LinearOperators
+  ParallelAmr
   Serialization
   Utilities
   )

src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp

@@ -37,6 +37,7 @@
 #include "Parallel/GlobalCache.hpp"
 #include "Parallel/InboxInserters.hpp"
 #include "Parallel/Invoke.hpp"
+#include "ParallelAlgorithms/Amr/Projectors/DefaultInitialize.hpp"
 #include "Utilities/ErrorHandling/Assert.hpp"
 #include "Utilities/GetOutput.hpp"
 #include "Utilities/Gsl.hpp"
@@ -207,13 +208,16 @@ struct PrepareAndSendMortarData<
     // Can't `db::get` the arguments for the boundary conditions within
     // `db::mutate`, so we extract the data to mutate and move it back in when
     // we're done.
+    // Possible optimization: also keep memory for the mortar data around in the
+    // DataBox. Currently the mortar data is extracted and erased by
+    // `apply_operator` anyway, so we can just create a new map here to avoid
+    // dealing with AMR resizing the mortars. When we keep the memory around,
+    // its size has to be adjusted when the mesh changes during AMR.
     typename PrimalFluxesTag::type primal_fluxes;
-    typename all_mortar_data_tag::type all_mortar_data;
-    db::mutate<PrimalFluxesTag, all_mortar_data_tag>(
-        [&primal_fluxes, &all_mortar_data](const auto local_primal_fluxes,
-                                           const auto local_all_mortar_data) {
+    typename all_mortar_data_tag::type all_mortar_data{};
+    db::mutate<PrimalFluxesTag>(
+        [&primal_fluxes](const auto local_primal_fluxes) {
           primal_fluxes = std::move(*local_primal_fluxes);
-          all_mortar_data = std::move(*local_all_mortar_data);
         },
         make_not_null(&box));
 
@@ -446,13 +450,27 @@ template <typename System, typename BackgroundTag = void>
 using initialize_operator = tmpl::list<
     detail::InitializeFacesMortarsAndBackground<System, BackgroundTag>>;
 
+/*!
+ * \brief AMR projectors for the tags added by `initialize_operator`
+ */
+template <typename System, typename BackgroundTag = void>
+using amr_projectors = tmpl::append<
+    tmpl::list<elliptic::dg::InitializeFacesAndMortars<
+        System::volume_dim, typename System::inv_metric_tag, BackgroundTag>>,
+    tmpl::conditional_t<
+        std::is_same_v<typename System::background_fields, tmpl::list<>>,
+        tmpl::list<>,
+        tmpl::list<elliptic::dg::InitializeBackground<
+            System::volume_dim, typename System::background_fields,
+            BackgroundTag>>>>;
+
 /*!
  * \brief Apply the DG operator to the `PrimalFieldsTag` and write the result to
  * the `OperatorAppliedToFieldsTag`
  *
- * Add this list to the action list of a parallel component to compute the
- * elliptic DG operator or its linearization. The operator involves a
- * communication between nearest-neighbor elements. See `elliptic::dg` for
+ * Add the `apply_actions` list to the action list of a parallel component to
+ * compute the elliptic DG operator or its linearization. The operator involves
+ * a communication between nearest-neighbor elements. See `elliptic::dg` for
  * details on the elliptic DG operator. Make sure to add
  * `elliptic::dg::Actions::initialize_operator` to the initialization phase of
  * your parallel component so the required DataBox tags are set up before
@@ -468,21 +486,37 @@ using initialize_operator = tmpl::list<
  * example when applying the nonlinear and linearized operator. They default to
  * the `PrimalFieldsTag` and the `PrimalFluxesTag`, meaning memory buffers
  * corresponding to these tags are set up in the DataBox.
+ *
+ * \par AMR
+ * Also add the `amr_projectors` to the list of AMR projectors to support AMR.
  */
 template <typename System, bool Linearized, typename TemporalIdTag,
           typename PrimalFieldsTag, typename PrimalFluxesTag,
           typename OperatorAppliedToFieldsTag,
           typename PrimalMortarFieldsTag = PrimalFieldsTag,
           typename PrimalMortarFluxesTag = PrimalFluxesTag>
-using apply_operator =
-    tmpl::list<detail::PrepareAndSendMortarData<
-                   System, Linearized, TemporalIdTag, PrimalFieldsTag,
-                   PrimalFluxesTag, OperatorAppliedToFieldsTag,
-                   PrimalMortarFieldsTag, PrimalMortarFluxesTag>,
-               detail::ReceiveMortarDataAndApplyOperator<
-                   System, Linearized, TemporalIdTag, PrimalFieldsTag,
-                   PrimalFluxesTag, OperatorAppliedToFieldsTag,
-                   PrimalMortarFieldsTag, PrimalMortarFluxesTag>>;
+struct DgOperator {
+ private:
+  static constexpr size_t Dim = System::volume_dim;
+
+ public:
+  using apply_actions =
+      tmpl::list<detail::PrepareAndSendMortarData<
+                     System, Linearized, TemporalIdTag, PrimalFieldsTag,
+                     PrimalFluxesTag, OperatorAppliedToFieldsTag,
+                     PrimalMortarFieldsTag, PrimalMortarFluxesTag>,
+                 detail::ReceiveMortarDataAndApplyOperator<
+                     System, Linearized, TemporalIdTag, PrimalFieldsTag,
+                     PrimalFluxesTag, OperatorAppliedToFieldsTag,
+                     PrimalMortarFieldsTag, PrimalMortarFluxesTag>>;
+  using amr_projectors = tmpl::list<::amr::projectors::DefaultInitialize<
+      PrimalFluxesTag, OperatorAppliedToFieldsTag,
+      ::Tags::Mortars<elliptic::dg::Tags::MortarData<
+                          typename TemporalIdTag::type,
+                          typename PrimalMortarFieldsTag::tags_list,
+                          typename PrimalMortarFluxesTag::tags_list>,
+                      Dim>>>;
+};
 
 /*!
  * \brief For linear systems, impose inhomogeneous boundary conditions as

src/Elliptic/DiscontinuousGalerkin/CMakeLists.txt

@@ -33,6 +33,7 @@ target_link_libraries(
   ErrorHandling
   FunctionsOfTime
   LinearOperators
+  ParallelAmr
   Spectral
   Utilities
   INTERFACE