Mass Accretion and other diagnostics

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

@@ -158,6 +158,7 @@
 #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/RotatingStar.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/RelativisticEuler/TovStar.hpp"
 #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
+#include "PointwiseFunctions/GeneralRelativity/DetAndInverseSpatialMetric.hpp"
 #include "PointwiseFunctions/GeneralRelativity/Surfaces/Tags.hpp"
 #include "PointwiseFunctions/Hydro/EquationsOfState/Factory.hpp"
 #include "PointwiseFunctions/Hydro/EquationsOfState/RegisterDerivedWithCharm.hpp"
@@ -310,7 +311,11 @@ struct EvolutionMetavars<tmpl::list<InterpolationTargetTags...>,
           ::Events::Tags::ObserverCoordinates<volume_dim, Frame::Inertial>,
           hydro::Tags::TransportVelocity<DataVector, volume_dim,
                                          Frame::Inertial>>,
-      hydro::Tags::InversePlasmaBetaCompute<DataVector>>;
+      hydro::Tags::InversePlasmaBetaCompute<DataVector>,
+      hydro::Tags::MassFluxCompute<DataVector, 3, ::Frame::Inertial>,
+      gr::Tags::SqrtDetSpatialMetric<DataVector>,
+      gr::Tags::SpatialMetric<DataVector, 3, ::Frame::Inertial>,
+      gr::Tags::Lapse<DataVector>>;
   using non_tensor_compute_tags = tmpl::list<
       tmpl::conditional_t<
           use_dg_subcell,

src/IO/H5/VolumeData.cpp

@@ -54,8 +54,6 @@ void append_element_extents_and_connectivity(
     const ElementVolumeData& element) {
   // Process the element extents
   const auto& extents = element.extents;
-  ASSERT(alg::none_of(extents, [](const size_t extent) { return extent == 1; }),
-         "We cannot generate connectivity for any single grid point elements.");
   if (extents.size() != dim) {
     ERROR("Trying to write data of dimensionality"
           << extents.size() << "but the VolumeData file has dimensionality"
@@ -70,14 +68,26 @@ void append_element_extents_and_connectivity(
   // out size without computing all the connectivities.
   const std::vector<int> connectivity = [&extents, &total_points_so_far]() {
     std::vector<int> local_connectivity;
-    for (const auto& cell : vis::detail::compute_cells(extents)) {
+    std::vector<size_t> extents_for_compute_cells =
+        extents;  // FIXME: changed up to 124.
+    // basically what this is doing is it removes all occurence of 1_st extent
+    // store the iterator pointing to the new end of the container after removal
+    //"removed" item technically present but shifted to the end by std::shift
+    // and they are actually removed at line 121
+    const auto new_end_it = std::remove(extents_for_compute_cells.begin(),
+                                        extents_for_compute_cells.end(), 1_st);
+    extents_for_compute_cells.erase(new_end_it,
+                                    extents_for_compute_cells.end());
+    for (const auto& cell :
+         vis::detail::compute_cells(extents_for_compute_cells)) {
       for (const auto& bounding_indices : cell.bounding_indices) {
         local_connectivity.emplace_back(*total_points_so_far +
                                         static_cast<int>(bounding_indices));
       }
     }
     return local_connectivity;
   }();
+
   *total_points_so_far += element_num_points;
   total_connectivity->insert(total_connectivity->end(), connectivity.begin(),
                              connectivity.end());
@@ -186,7 +196,7 @@ VolumeData::VolumeData(const bool subfile_exists, detail::OpenGroup&& group,
 // an `observation_group` in a `VolumeData` file.
 void VolumeData::write_volume_data(
     const size_t observation_id, const double observation_value,
-    const std::vector<ElementVolumeData>& elements,
+    const std::vector<ElementVolumeData>& elements2,
     const std::optional<std::vector<char>>& serialized_domain,
     const std::optional<std::vector<char>>& serialized_functions_of_time) {
   const std::string path = "ObservationId" + std::to_string(observation_id);
@@ -209,6 +219,19 @@ void VolumeData::write_volume_data(
                << component.name << "'.");
     return component.name;
   };
+
+  std::vector<ElementVolumeData> elements = elements2;  // FIXME:284-293
+  // can I pop off the 1 grid point parts of the elements?
+  for (auto& element : elements) {
+    for (unsigned i = element.extents.size(); i-- > 0;) {
+      if (element.extents[i] == 1_st) {
+        element.extents.erase(element.extents.begin() + i);
+        element.basis.erase(element.basis.begin() + i);
+        element.quadrature.erase(element.quadrature.begin() + i);
+      }
+    }
+  }
+
   const std::vector<std::string> component_names(
       boost::make_transform_iterator(elements.front().tensor_components.begin(),
                                      get_component_name),
@@ -220,6 +243,7 @@ void VolumeData::write_volume_data(
   // should have the same dimensionality)
   if (not contains_attribute(volume_data_group_.id(), "", "dimension")) {
     h5::write_to_attribute(volume_data_group_.id(), "dimension",
+                           //  1);
                            elements.front().extents.size());
   }
   const auto dim =

src/ParallelAlgorithms/Events/CMakeLists.txt

@@ -25,6 +25,7 @@ spectre_target_headers(
   ObserveDataBox.hpp
   ObserveAtExtremum.hpp
   ObserveFields.hpp
+  ObserveInterpolatedIntegralData.hpp
   ObserveNorms.hpp
   ObserveTimeStep.hpp
   ObserveTimeStepVolume.hpp

src/ParallelAlgorithms/Events/Factory.hpp

@@ -9,6 +9,7 @@
 #include "ParallelAlgorithms/Events/ErrorIfDataTooBig.hpp"
 #include "ParallelAlgorithms/Events/ObserveAdaptiveSteppingDiagnostics.hpp"
 #include "ParallelAlgorithms/Events/ObserveFields.hpp"
+#include "ParallelAlgorithms/Events/ObserveInterpolatedIntegralData.hpp"
 #include "ParallelAlgorithms/Events/ObserveNorms.hpp"
 #include "ParallelAlgorithms/Events/ObserveTimeStep.hpp"
 #include "Time/ChangeSlabSize/Event.hpp"
@@ -21,6 +22,8 @@ using field_observations = tmpl::flatten<tmpl::list<
     ::Events::ErrorIfDataTooBig<VolumeDim, Fields, NonTensorComputeTagsList>,
     ObserveFields<VolumeDim, Fields, NonTensorComputeTagsList,
                   ArraySectionIdTag>,
+    ObserveInterpolatedIntegralData<VolumeDim, Fields, NonTensorComputeTagsList,
+                                    ArraySectionIdTag>,
     ::Events::ObserveNorms<Fields, NonTensorComputeTagsList,
                            ArraySectionIdTag>>>;
 }  // namespace dg::Events