Merge pull request sxs-collaboration#6349 from PunJustice/LinearFluxes

src/Elliptic/DiscontinuousGalerkin/Actions/ApplyOperator.hpp

@@ -27,6 +27,7 @@
 #include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
 #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
 #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
+#include "Elliptic/Systems/GetFluxesComputer.hpp"
 #include "Elliptic/Systems/GetSourcesComputer.hpp"
 #include "Elliptic/Utilities/ApplyAt.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/HasReceivedFromAllMortars.hpp"
@@ -123,7 +124,7 @@ template <typename System, bool Linearized, typename TemporalIdTag,
           typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
           typename PrimalMortarFluxesTag,
           typename FluxesArgsTags =
-              typename System::fluxes_computer::argument_tags,
+              elliptic::get_fluxes_argument_tags<System, Linearized>,
           typename SourcesArgsTags =
               elliptic::get_sources_argument_tags<System, Linearized>>
 struct PrepareAndSendMortarData;
@@ -179,28 +180,28 @@ struct PrepareAndSendMortarData<
     const auto& boundary_conditions =
         db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box).at(
             element_id.block_id());
-    const auto apply_boundary_condition =
-        [&box, &boundary_conditions, &element_id](
-            const Direction<Dim>& direction, auto&&... fields_and_fluxes) {
-          ASSERT(
-              boundary_conditions.contains(direction),
-              "No boundary condition is available in block "
-                  << element_id.block_id() << " in direction " << direction
-                  << ". Make sure you are setting up boundary conditions when "
-                     "creating the domain.");
-          ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
-                     boundary_conditions.at(direction).get()) != nullptr,
-                 "The boundary condition in block "
-                     << element_id.block_id() << " in direction " << direction
-                     << " has an unexpected type. Make sure it derives off the "
-                        "'boundary_conditions_base' class set in the system.");
-          const auto& boundary_condition =
-              dynamic_cast<const BoundaryConditionsBase&>(
-                  *boundary_conditions.at(direction));
-          elliptic::apply_boundary_condition<Linearized>(
-              boundary_condition, box, direction,
-              std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
-        };
+    const auto apply_boundary_condition = [&box, &boundary_conditions,
+                                           &element_id](
+                                              const Direction<Dim>& direction,
+                                              auto&&... fields_and_fluxes) {
+      ASSERT(boundary_conditions.contains(direction),
+             "No boundary condition is available in block "
+                 << element_id.block_id() << " in direction " << direction
+                 << ". Make sure you are setting up boundary conditions when "
+                    "creating the domain.");
+      ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
+                 boundary_conditions.at(direction).get()) != nullptr,
+             "The boundary condition in block "
+                 << element_id.block_id() << " in direction " << direction
+                 << " has an unexpected type. Make sure it derives off the "
+                    "'boundary_conditions_base' class set in the system.");
+      const auto& boundary_condition =
+          dynamic_cast<const BoundaryConditionsBase&>(
+              *boundary_conditions.at(direction));
+      elliptic::apply_boundary_condition<Linearized>(
+          boundary_condition, box, direction,
+          std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
+    };
 
     // 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
@@ -260,8 +261,7 @@ struct PrepareAndSendMortarData<
         // Make a copy of the local boundary data on the mortar to send to the
         // neighbor
         auto remote_boundary_data_on_mortar =
-            get<all_mortar_data_tag>(box).at(mortar_id).local_data(
-                temporal_id);
+            get<all_mortar_data_tag>(box).at(mortar_id).local_data(temporal_id);
         // Reorient the data to the neighbor orientation if necessary
         if (not orientation.is_aligned()) {
           remote_boundary_data_on_mortar.orient_on_slice(
@@ -274,7 +274,7 @@ struct PrepareAndSendMortarData<
                 ::dg::MortarId<Dim>{direction_from_neighbor, element.id()},
                 std::move(remote_boundary_data_on_mortar)));
       }  // loop over neighbors in direction
-    }    // loop over directions
+    }  // loop over directions
 
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
@@ -289,7 +289,7 @@ template <typename System, bool Linearized, typename TemporalIdTag,
           typename OperatorAppliedToFieldsTag, typename PrimalMortarFieldsTag,
           typename PrimalMortarFluxesTag,
           typename FluxesArgsTags =
-              typename System::fluxes_computer::argument_tags,
+              elliptic::get_fluxes_argument_tags<System, Linearized>,
           typename SourcesArgsTags =
               elliptic::get_sources_argument_tags<System, Linearized>>
 struct ReceiveMortarDataAndApplyOperator;
@@ -359,9 +359,10 @@ struct ReceiveMortarDataAndApplyOperator<
     };
 
     // Apply DG operator
-    using fluxes_args_tags = typename System::fluxes_computer::argument_tags;
+    using fluxes_args_tags =
+        typename elliptic::get_fluxes_argument_tags<System, Linearized>;
     using fluxes_args_volume_tags =
-        typename System::fluxes_computer::volume_tags;
+        typename elliptic::get_fluxes_volume_tags<System, Linearized>;
     DirectionMap<Dim, std::tuple<decltype(db::get<FluxesArgsTags>(box))...>>
         fluxes_args_on_faces{};
     for (const auto& direction : Direction<Dim>::all_directions()) {
@@ -527,7 +528,7 @@ struct DgOperator {
  */
 template <typename System, typename FixedSourcesTag,
           typename FluxesArgsTags =
-              typename System::fluxes_computer::argument_tags,
+              elliptic::get_fluxes_argument_tags<System, false>,
           typename SourcesArgsTags =
               elliptic::get_sources_argument_tags<System, false>>
 struct ImposeInhomogeneousBoundaryConditionsOnSource;
@@ -563,28 +564,28 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
     const auto& boundary_conditions =
         db::get<domain::Tags::ExternalBoundaryConditions<Dim>>(box).at(
             element_id.block_id());
-    const auto apply_boundary_condition =
-        [&box, &boundary_conditions, &element_id](
-            const Direction<Dim>& direction, auto&&... fields_and_fluxes) {
-          ASSERT(
-              boundary_conditions.contains(direction),
-              "No boundary condition is available in block "
-                  << element_id.block_id() << " in direction " << direction
-                  << ". Make sure you are setting up boundary conditions when "
-                     "creating the domain.");
-          ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
-                     boundary_conditions.at(direction).get()) != nullptr,
-                 "The boundary condition in block "
-                     << element_id.block_id() << " in direction " << direction
-                     << " has an unexpected type. Make sure it derives off the "
-                        "'boundary_conditions_base' class set in the system.");
-          const auto& boundary_condition =
-              dynamic_cast<const BoundaryConditionsBase&>(
-                  *boundary_conditions.at(direction));
-          elliptic::apply_boundary_condition<false>(
-              boundary_condition, box, direction,
-              std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
-        };
+    const auto apply_boundary_condition = [&box, &boundary_conditions,
+                                           &element_id](
+                                              const Direction<Dim>& direction,
+                                              auto&&... fields_and_fluxes) {
+      ASSERT(boundary_conditions.contains(direction),
+             "No boundary condition is available in block "
+                 << element_id.block_id() << " in direction " << direction
+                 << ". Make sure you are setting up boundary conditions when "
+                    "creating the domain.");
+      ASSERT(dynamic_cast<const BoundaryConditionsBase*>(
+                 boundary_conditions.at(direction).get()) != nullptr,
+             "The boundary condition in block "
+                 << element_id.block_id() << " in direction " << direction
+                 << " has an unexpected type. Make sure it derives off the "
+                    "'boundary_conditions_base' class set in the system.");
+      const auto& boundary_condition =
+          dynamic_cast<const BoundaryConditionsBase&>(
+              *boundary_conditions.at(direction));
+      elliptic::apply_boundary_condition<false>(
+          boundary_condition, box, direction,
+          std::forward<decltype(fields_and_fluxes)>(fields_and_fluxes)...);
+    };
 
     // 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
@@ -596,9 +597,10 @@ struct ImposeInhomogeneousBoundaryConditionsOnSource<
         },
         make_not_null(&box));
 
-    using fluxes_args_tags = typename System::fluxes_computer::argument_tags;
+    using fluxes_args_tags =
+        typename elliptic::get_fluxes_argument_tags<System, false>;
     using fluxes_args_volume_tags =
-        typename System::fluxes_computer::volume_tags;
+        elliptic::get_fluxes_volume_tags<System, false>;
     DirectionMap<Dim, std::tuple<decltype(db::get<FluxesArgsTags>(box))...>>
         fluxes_args_on_faces{};
     for (const auto& direction : Direction<Dim>::all_directions()) {

src/Elliptic/DiscontinuousGalerkin/DgOperator.hpp

@@ -21,6 +21,7 @@
 #include "Domain/Structure/ElementId.hpp"
 #include "Domain/Structure/IndexToSliceAt.hpp"
 #include "Elliptic/Protocols/FirstOrderSystem.hpp"
+#include "Elliptic/Systems/GetFluxesComputer.hpp"
 #include "Elliptic/Systems/GetSourcesComputer.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/ApplyMassMatrix.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
@@ -222,7 +223,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
       tt::assert_conforms_to_v<System, elliptic::protocols::FirstOrderSystem>);
 
   static constexpr size_t Dim = System::volume_dim;
-  using FluxesComputer = typename System::fluxes_computer;
+  using FluxesComputer = elliptic::get_fluxes_computer<System, Linearized>;
   using SourcesComputer = elliptic::get_sources_computer<System, Linearized>;
 
   struct AllDirections {
@@ -511,7 +512,7 @@ struct DgOperatorImpl<System, Linearized, tmpl::list<PrimalFields...>,
         all_mortar_data->at(mortar_id).remote_insert(temporal_id,
                                                      std::move(boundary_data));
       }  // if (is_internal)
-    }    // loop directions
+    }  // loop directions
   }
 
   // --- This is essentially a break to communicate the mortar data ---

src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/InitializeSubdomain.hpp

@@ -37,6 +37,7 @@
 #include "Domain/Tags/NeighborMesh.hpp"
 #include "Elliptic/DiscontinuousGalerkin/Initialization.hpp"
 #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp"
+#include "Elliptic/Systems/GetFluxesComputer.hpp"
 #include "Elliptic/Utilities/ApplyAt.hpp"
 #include "Elliptic/Utilities/GetAnalyticData.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
@@ -160,10 +161,10 @@ struct InitializeSubdomain {
   // necessary for the DG operator, i.e. the background fields in the
   // System::fluxes_computer::argument_tags
   using fluxes_non_background_args =
-      tmpl::list_difference<typename System::fluxes_computer::argument_tags,
+      tmpl::list_difference<elliptic::get_fluxes_argument_tags<System, true>,
                             typename System::background_fields>;
   using background_fields_internal =
-      tmpl::list_difference<typename System::fluxes_computer::argument_tags,
+      tmpl::list_difference<elliptic::get_fluxes_argument_tags<System, true>,
                             fluxes_non_background_args>;
   // Slice all background fields to external boundaries for use in boundary
   // conditions
@@ -245,7 +246,7 @@ struct InitializeSubdomain {
         // Faces on the other side of the overlapped element's mortars
         initialize_remote_faces(make_not_null(&box), overlap_id);
       }  // neighbors in direction
-    }    // directions
+    }  // directions
     return {Parallel::AlgorithmExecution::Continue, std::nullopt};
   }
 
@@ -382,7 +383,7 @@ struct InitializeSubdomain {
               box, std::make_tuple(overlap_id, mortar_id));
         }
       }  // neighbors
-    }    // internal directions
+    }  // internal directions
   }
 };
 

src/Elliptic/DiscontinuousGalerkin/SubdomainOperator/SubdomainOperator.hpp

@@ -30,6 +30,7 @@
 #include "Elliptic/DiscontinuousGalerkin/DgOperator.hpp"
 #include "Elliptic/DiscontinuousGalerkin/SubdomainOperator/Tags.hpp"
 #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
+#include "Elliptic/Systems/GetFluxesComputer.hpp"
 #include "Elliptic/Systems/GetSourcesComputer.hpp"
 #include "Elliptic/Utilities/ApplyAt.hpp"
 #include "NumericalAlgorithms/DiscontinuousGalerkin/MortarHelpers.hpp"
@@ -82,12 +83,16 @@ struct make_neighbor_mortars_tag_impl {
  *   `elliptic::dg::subdomain_operator::InitializeSubdomain`.
  * - All `System::fluxes_computer::argument_tags` and
  *   `System::sources_computer::argument_tags` (or
+ *   `System::fluxes_computer_linearized::argument_tags` and
  *   `System::sources_computer_linearized::argument_tags` for nonlinear
  *   systems), except those listed in `ArgsTagsFromCenter`. The latter will be
  *   taken from the central element's DataBox, so they don't need to be made
  *   available on overlaps.
- * - The `System::fluxes_computer::argument_tags` on internal and external
- *   interfaces, except those listed in `System::fluxes_computer::volume_tags`.
+ * - The `System::fluxes_computer::argument_tags` (or
+ *   `System::fluxes_computer_linearized::argument_tags`) on internal and
+ *   external interfaces, except those listed in
+ *   `System::fluxes_computer::volume_tags` (or
+ *   `System::fluxes_computer_linearized::volume_tags`).
  *
  * Some of these tags may require communication between elements. For example,
  * nonlinear system fields are constant background fields for the linearized DG
@@ -174,18 +179,20 @@ struct SubdomainOperator
                       Dim>,
       ::Tags::Mortars<elliptic::dg::Tags::PenaltyFactor, Dim>,
       elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>;
-  using fluxes_args_tags = typename System::fluxes_computer::argument_tags;
+  using fluxes_args_tags =
+      elliptic::get_fluxes_argument_tags<System, linearized>;
   using sources_args_tags =
       elliptic::get_sources_argument_tags<System, linearized>;
 
   // We need the fluxes args also on interfaces (internal and external). The
   // volume tags are the subset that don't have to be taken from interfaces.
-  using fluxes_args_volume_tags = typename System::fluxes_computer::volume_tags;
+  using fluxes_args_volume_tags =
+      elliptic::get_fluxes_volume_tags<System, linearized>;
 
   // These tags can be taken directly from the central element's DataBox, even
   // when evaluating neighbors
   using args_tags_from_center = tmpl::remove_duplicates<tmpl::push_back<
-      typename System::fluxes_computer::const_global_cache_tags,
+      elliptic::get_fluxes_const_global_cache_tags<System, linearized>,
       elliptic::dg::Tags::Massive, elliptic::dg::Tags::Formulation>>;
 
   // Data on neighbors is stored in the central element's DataBox in
@@ -600,8 +607,8 @@ struct SubdomainOperator
                 .remote_insert(temporal_id, std::move(zero_mortar_data));
           }
         }  // loop over neighbor's mortars
-      }    // loop over neighbors
-    }      // loop over directions
+      }  // loop over neighbors
+    }  // loop over directions
 
     // 2. Apply the operator on all elements in the subdomain
     //
@@ -613,8 +620,7 @@ struct SubdomainOperator
               make_not_null(&central_mortar_data_), operand.element_data,
               central_primal_fluxes_, args...);
         },
-        box, temporal_id, fluxes_args_on_faces, sources_args,
-        data_is_zero);
+        box, temporal_id, fluxes_args_on_faces, sources_args, data_is_zero);
     // Apply on neighbors
     for (const auto& [direction, neighbors] : central_element.neighbors()) {
       const auto& orientation = neighbors.orientation();
@@ -673,7 +679,7 @@ struct SubdomainOperator
             extended_results_.at(overlap_id), neighbor_mesh.extents(),
             overlap_extent, direction_from_neighbor);
       }  // loop over neighbors
-    }    // loop over directions
+    }  // loop over directions
   }
 
   // NOLINTNEXTLINE(google-runtime-references)

src/Elliptic/Protocols/FirstOrderSystem.hpp

@@ -178,8 +178,8 @@ struct FirstOrderSystem {
 
     using fluxes_computer = typename ConformingType::fluxes_computer;
     using sources_computer = typename ConformingType::sources_computer;
-    using fluxes_argument_tags = typename fluxes_computer::argument_tags;
-    static_assert(tt::is_a_v<tmpl::list, fluxes_argument_tags>);
+    static_assert(
+        tt::is_a_v<tmpl::list, typename fluxes_computer::argument_tags>);
 
     using boundary_conditions_base =
         typename ConformingType::boundary_conditions_base;

src/Elliptic/Systems/CMakeLists.txt

@@ -9,6 +9,7 @@ spectre_target_headers(
   ${LIBRARY}
   INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
   HEADERS
+  GetFluxesComputer.hpp
   GetSourcesComputer.hpp
   )