Merge pull request #4750 from knelli2/reduce_allocations

Further reduce number of allocations in InternalMortarDataImpl

docs/DevGuide/ImplementingVectors.md

@@ -282,17 +282,18 @@ of `std::move`.
 This function intentionally generates an error when assigning values from one
 vector to a differently sized, non-owning vector (made non-owning by use of
 `set_data_ref`). The assertion test which calls this function should search for
-the string "Must copy into same size". Three forms of the test are provided,
-which are switched between using a value from the enum `RefSizeErrorTestKind` in
-the first function argument:
+the string "Must copy/move/assign into same size". Three forms of the test are
+provided, which are switched between using a value from the enum
+`RefSizeErrorTestKind` in the first function argument:
 - `RefSizeErrorTestKind::Copy`: tests that the size error is appropriately
-  generated when copying to a non-owning vector of the wrong size.
+  generated when copying to a non-owning vector of the wrong size. This has
+  "copy" in the message.
 - `RefSizeErrorTestKind::ExpressionAssign`: tests that the size error is
   appropriately generated when assigning the result of a mathematical expression
-  to a non-owning vector of the wrong size.
+  to a non-owning vector of the wrong size. This has "assign" in the message.
 - `RefSizeErrorTestKind::Move`: tests that the size error is appropriately
   generated when a vector is `std::move`d into a non-owning vector of the wrong
-  size
+  size. This has "move" in the message.
 
 ## `TestHelpers::VectorImpl::test_functions_with_vector_arguments()`
 

src/DataStructures/VectorImpl.hpp

@@ -392,7 +392,7 @@ VectorImpl<T, VectorType, StaticSize>::operator=(
       reset_pointer_vector(size());
       rhs.clear();
     } else {
-      ASSERT(rhs.size() == size(), "Must copy into same size, not "
+      ASSERT(rhs.size() == size(), "Must move into same size, not "
                                        << rhs.size() << " into " << size());
       if (LIKELY(data() != rhs.data())) {
         std::memcpy(data(), rhs.data(), size() * sizeof(value_type));
@@ -432,7 +432,7 @@ VectorImpl<T, VectorType, StaticSize>::operator=(
     owned_data_ = heap_alloc_if_necessary((*expression).size());
     reset_pointer_vector((*expression).size());
   } else if (not owning_) {
-    ASSERT((*expression).size() == size(), "Must copy into same size, not "
+    ASSERT((*expression).size() == size(), "Must assign into same size, not "
                                                << (*expression).size()
                                                << " into " << size());
   }

src/Evolution/DiscontinuousGalerkin/Actions/ComputeTimeDerivative.hpp

@@ -76,6 +76,10 @@ template <typename T>
 struct get_dg_package_temporary_tags {
   using type = typename T::dg_package_data_temporary_tags;
 };
+template <typename T>
+struct get_dg_package_field_tags {
+  using type = typename T::dg_package_field_tags;
+};
 template <typename System, typename T>
 struct get_primitive_tags_for_face {
   using type = typename get_primitive_vars<
@@ -420,6 +424,12 @@ ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers, LocalTimeStepping>::
           tmpl::list<dg_package_data_projected_tags,
                      detail::inverse_spatial_metric_tag<EvolutionSystem>>,
           detail::OneOverNormalVectorMagnitude, detail::NormalVector<Dim>>>>;
+  // To avoid additional allocations in internal_mortar_data, we provide a
+  // buffer used to compute the packaged data before it has to be projected to
+  // the mortar. We get all mortar tags for similar reasons as described above
+  using all_mortar_tags = tmpl::remove_duplicates<tmpl::flatten<
+      tmpl::transform<derived_boundary_corrections,
+                      detail::get_dg_package_field_tags<tmpl::_1>>>>;
 
   // We also don't use the number of volume mesh grid points. We instead use the
   // max number of grid points from each face. That way, our allocation will be
@@ -458,6 +468,7 @@ ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers, LocalTimeStepping>::
       ::Tags::div, db::wrap_tags_in<::Tags::Flux, flux_variables,
                                     tmpl::size_t<Dim>, Frame::Inertial>>>;
   using VarsFaceTemporaries = Variables<all_face_temporary_tags>;
+  using DgPackagedDataVarsOnFace = Variables<all_mortar_tags>;
   const size_t number_of_grid_points = mesh.number_of_grid_points();
   auto buffer = cpp20::make_unique_for_overwrite<double[]>(
       (VarsTemporaries::number_of_independent_components +
@@ -467,7 +478,8 @@ ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers, LocalTimeStepping>::
           number_of_grid_points +
       // Different number of grid points. See explanation above where
       // num_face_temporary_grid_points is defined
-      VarsFaceTemporaries::number_of_independent_components *
+      (VarsFaceTemporaries::number_of_independent_components +
+       DgPackagedDataVarsOnFace::number_of_independent_components) *
           num_face_temporary_grid_points);
   VarsTemporaries temporaries{
       &buffer[0], VarsTemporaries::number_of_independent_components *
@@ -500,6 +512,18 @@ ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers, LocalTimeStepping>::
       // num_face_temporary_grid_points is defined
       VarsFaceTemporaries::number_of_independent_components *
           num_face_temporary_grid_points);
+  gsl::span<double> packaged_data_buffer = gsl::make_span<double>(
+      &buffer[(VarsTemporaries::number_of_independent_components +
+               VarsFluxes::number_of_independent_components +
+               VarsPartialDerivatives::number_of_independent_components +
+               VarsDivFluxes::number_of_independent_components) *
+                  number_of_grid_points +
+              VarsFaceTemporaries::number_of_independent_components *
+                  num_face_temporary_grid_points],
+      // Different number of grid points. See explanation above where
+      // num_face_temporary_grid_points is defined
+      DgPackagedDataVarsOnFace::number_of_independent_components *
+          num_face_temporary_grid_points);
 
   const Scalar<DataVector>* det_inverse_jacobian = nullptr;
   if constexpr (tmpl::size<flux_variables>::value != 0) {
@@ -549,7 +573,7 @@ ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers, LocalTimeStepping>::
       "final.");
   tmpl::for_each<derived_boundary_corrections>(
       [&boundary_correction, &box, &partial_derivs, &primitive_vars,
-       &temporaries, &volume_fluxes,
+       &temporaries, &volume_fluxes, &packaged_data_buffer,
        &face_temporaries](auto derived_correction_v) {
         using DerivedCorrection =
             tmpl::type_from<decltype(derived_correction_v)>;
@@ -562,6 +586,7 @@ ComputeTimeDerivative<Dim, EvolutionSystem, DgStepChoosers, LocalTimeStepping>::
           //  - evolution::dg::Tags::MortarData<Dim>
           detail::internal_mortar_data<EvolutionSystem, Dim>(
               make_not_null(&box), make_not_null(&face_temporaries),
+              make_not_null(&packaged_data_buffer),
               dynamic_cast<const DerivedCorrection&>(boundary_correction),
               db::get<variables_tag>(box), volume_fluxes, temporaries,
               primitive_vars,

src/Evolution/DiscontinuousGalerkin/Actions/InternalMortarDataImpl.hpp

@@ -52,6 +52,7 @@ void internal_mortar_data_impl(
         boost::hash<std::pair<Direction<Dim>, ElementId<Dim>>>>*>
         mortar_data_ptr,
     const gsl::not_null<gsl::span<double>*> face_temporaries,
+    const gsl::not_null<gsl::span<double>*> packaged_data_buffer,
     const BoundaryCorrection& boundary_correction,
     const Variables<typename System::variables_tag::tags_list>&
         volume_evolved_vars,
@@ -97,8 +98,6 @@ void internal_mortar_data_impl(
   FieldsOnFace fields_on_face{};
   std::optional<tnsr::I<DataVector, Dim>> face_mesh_velocity{};
   for (const auto& [direction, neighbors_in_direction] : element.neighbors()) {
-    (void)neighbors_in_direction;  // unused variable
-
     const Mesh<Dim - 1> face_mesh =
         volume_mesh.slice_away(direction.dimension());
 
@@ -230,57 +229,108 @@ void internal_mortar_data_impl(
            "have been. Direction: "
                << direction);
 
-    // We point the Variables inside the DataVector because the DataVector
-    // will be moved later on so we want it owning its data
-    DataVector packaged_data_buffer{
+    const size_t total_face_size =
         face_mesh.number_of_grid_points() *
-        Variables<mortar_tags_list>::number_of_independent_components};
-    Variables<mortar_tags_list> packaged_data{packaged_data_buffer.data(),
-                                              packaged_data_buffer.size()};
+        Variables<mortar_tags_list>::number_of_independent_components;
+    Variables<mortar_tags_list> packaged_data{};
+
+    // This is the case where we only have one neighbor in this direction, so we
+    // may or may not have to do any projection. If we don't have to do
+    // projection, then we can use the local_mortar_data itself to calculate the
+    // dg_package_data. However, if we need to project, then we hae to use the
+    // packaged_data_buffer that was passed in.
+    if (neighbors_in_direction.size() == 1) {
+      const auto& neighbor = *neighbors_in_direction.begin();
+      const auto mortar_id = std::pair{direction, neighbor};
+      const auto& mortar_mesh = mortar_meshes.at(mortar_id);
+      const auto& mortar_size = mortar_sizes.at(mortar_id);
+
+      // Have to use packaged_data_buffer
+      if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
+        // The face mesh will be assigned below along with ensuring the size of
+        // the mortar data is correct
+        packaged_data.set_data_ref(packaged_data_buffer->data(),
+                                   total_face_size);
+      } else {
+        // Can use the local_mortar_data
+        auto& local_mortar_data_opt =
+            mortar_data_ptr->at(mortar_id).local_mortar_data();
+        // If this isn't the first time, set the face mesh
+        if (LIKELY(local_mortar_data_opt.has_value())) {
+          local_mortar_data_opt->first = face_mesh;
+        } else {
+          // Otherwise we need to initialize the pair. If we don't do this, then
+          // the DataVector will be non-owning which we don't want
+          local_mortar_data_opt =
+              std::optional{std::pair{face_mesh, DataVector{}}};
+        }
+
+        // Always set the time step ID
+        mortar_data_ptr->at(mortar_id).time_step_id() = temporal_id;
+
+        DataVector& local_mortar_data = local_mortar_data_opt->second;
+
+        // Do a destructive resize to account for potential p-refinement
+        local_mortar_data.destructive_resize(total_face_size);
+
+        packaged_data.set_data_ref(local_mortar_data.data(),
+                                   local_mortar_data.size());
+      }
+    } else {
+      // In this case, we have multiple neighbors in this direction so all will
+      // need to project their data which means we use the
+      // packaged_data_buffer to calculate the dg_package_data
+      packaged_data.set_data_ref(packaged_data_buffer->data(), total_face_size);
+    }
 
-    const double max_abs_char_speed_on_face = detail::dg_package_data<System>(
+    detail::dg_package_data<System>(
         make_not_null(&packaged_data), boundary_correction, fields_on_face,
         get<evolution::dg::Tags::NormalCovector<Dim>>(
             *normal_covector_and_magnitude_ptr->at(direction)),
         face_mesh_velocity, dg_package_data_projected_tags{},
         package_data_volume_args...);
-    (void)max_abs_char_speed_on_face;
 
     // Perform step 3
+    // This will only do something if
+    //  a) we have multiple neighbors in this direction
+    // or
+    //  b) the one (and only) neighbor in this direction needed projection
     for (const auto& neighbor : neighbors_in_direction) {
       const auto mortar_id = std::make_pair(direction, neighbor);
       const auto& mortar_mesh = mortar_meshes.at(mortar_id);
       const auto& mortar_size = mortar_sizes.at(mortar_id);
 
-      // Project the data from the face to the mortar if necessary.
-      // We can't move the data or modify it in-place when projecting, because
-      // in that case the face may touch two mortars so we need to keep the
-      // data around.
-      DataVector boundary_data_on_mortar{};
       if (Spectral::needs_projection(face_mesh, mortar_mesh, mortar_size)) {
-        boundary_data_on_mortar = DataVector{
+        auto& local_mortar_data_opt =
+            mortar_data_ptr->at(mortar_id).local_mortar_data();
+
+        // If this isn't the first time, set the face mesh
+        if (LIKELY(local_mortar_data_opt.has_value())) {
+          local_mortar_data_opt->first = face_mesh;
+        } else {
+          // If we don't do this, then the DataVector will be non-owning which
+          // we don't want
+          local_mortar_data_opt =
+              std::optional{std::pair{face_mesh, DataVector{}}};
+        }
+
+        // Set the time id since above we only set it for cases that didn't need
+        // projection
+        mortar_data_ptr->at(mortar_id).time_step_id() = temporal_id;
+
+        DataVector& local_mortar_data = local_mortar_data_opt->second;
+
+        // Do a destructive resize to account for potential p-refinement
+        local_mortar_data.destructive_resize(
             mortar_mesh.number_of_grid_points() *
-            Variables<mortar_tags_list>::number_of_independent_components};
-        Variables<mortar_tags_list> projected_packaged_data{
-            boundary_data_on_mortar.data(), boundary_data_on_mortar.size()};
-        // Since projected_packaged_data is non-owning, if we tried to
-        // move-assign it with the result of ::dg::project_to_mortar, it would
-        // then become owning and the buffer it previously pointed to wouldn't
-        // be filled (and we want it to be filled). So instead force a
-        // copy-assign by having an intermediary, data_to_copy.
-        const auto data_to_copy = ::dg::project_to_mortar(
-            packaged_data, face_mesh, mortar_mesh, mortar_size);
-        projected_packaged_data = data_to_copy;
+            Variables<mortar_tags_list>::number_of_independent_components);
 
-      } else {
-        ASSERT(neighbors_in_direction.size() == 1,
-               "Expected only 1 neighbor in the local direction but got "
-                   << neighbors_in_direction.size() << "instead");
-        boundary_data_on_mortar = std::move(packaged_data_buffer);
+        Variables<mortar_tags_list> projected_packaged_data{
+            local_mortar_data.data(), local_mortar_data.size()};
+        ::dg::project_to_mortar(make_not_null(&projected_packaged_data),
+                                packaged_data, face_mesh, mortar_mesh,
+                                mortar_size);
       }
-
-      mortar_data_ptr->at(mortar_id).insert_local_mortar_data(
-          temporal_id, face_mesh, std::move(boundary_data_on_mortar));
     }
   }
 }
@@ -290,6 +340,7 @@ template <typename System, size_t Dim, typename BoundaryCorrection,
 void internal_mortar_data(
     const gsl::not_null<db::DataBox<DbTagsList>*> box,
     const gsl::not_null<gsl::span<double>*> face_temporaries,
+    const gsl::not_null<gsl::span<double>*> packaged_data_buffer,
     const BoundaryCorrection& boundary_correction,
     const Variables<typename System::variables_tag::tags_list>
         evolved_variables,
@@ -305,7 +356,7 @@ void internal_mortar_data(
   db::mutate<evolution::dg::Tags::NormalCovectorAndMagnitude<Dim>,
              evolution::dg::Tags::MortarData<Dim>>(
       box,
-      [&boundary_correction, &face_temporaries,
+      [&boundary_correction, &face_temporaries, &packaged_data_buffer,
        &element = db::get<domain::Tags::Element<Dim>>(*box), &evolved_variables,
        &logical_to_inertial_inverse_jacobian =
            db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
@@ -322,11 +373,11 @@ void internal_mortar_data(
           const auto mortar_data_ptr, const auto&... package_data_volume_args) {
         detail::internal_mortar_data_impl<System>(
             normal_covector_and_magnitude_ptr, mortar_data_ptr,
-            face_temporaries, boundary_correction, evolved_variables,
-            volume_fluxes, temporaries, primitive_vars, element, mesh,
-            mortar_meshes, mortar_sizes, time_step_id, moving_mesh_map,
-            mesh_velocity, logical_to_inertial_inverse_jacobian,
-            package_data_volume_args...);
+            face_temporaries, packaged_data_buffer, boundary_correction,
+            evolved_variables, volume_fluxes, temporaries, primitive_vars,
+            element, mesh, mortar_meshes, mortar_sizes, time_step_id,
+            moving_mesh_map, mesh_velocity,
+            logical_to_inertial_inverse_jacobian, package_data_volume_args...);
       },
       db::get<PackageDataVolumeTags>(*box)...);
 }

src/Evolution/DiscontinuousGalerkin/MortarData.hpp

@@ -145,6 +145,8 @@ class MortarData {
 
   const TimeStepId& time_step_id() const { return time_step_id_; }
 
+  TimeStepId& time_step_id() { return time_step_id_; }
+
   auto local_mortar_data() const
       -> const std::optional<std::pair<Mesh<Dim - 1>, DataVector>>& {
     return local_mortar_data_;
@@ -155,6 +157,16 @@ class MortarData {
     return neighbor_mortar_data_;
   }
 
+  auto local_mortar_data()
+      -> std::optional<std::pair<Mesh<Dim - 1>, DataVector>>& {
+    return local_mortar_data_;
+  }
+
+  auto neighbor_mortar_data()
+      -> std::optional<std::pair<Mesh<Dim - 1>, DataVector>>& {
+    return neighbor_mortar_data_;
+  }
+
   // NOLINTNEXTLINE(google-runtime-references)
   void pup(PUP::er& p);