surf integral

jyoo1042 · Feb 12, 2024 · e640758 · e640758
1 parent 6b39f51
commit e640758
Show file tree

Hide file tree

Showing 5 changed files with 493 additions and 5 deletions.
diff --git a/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp b/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,9 @@ struct EvolutionMetavars<tmpl::list<InterpolationTargetTags...>,
           hydro::Tags::TransportVelocity<DataVector, volume_dim,
                                          Frame::Inertial>>,
       hydro::Tags::InversePlasmaBetaCompute<DataVector>,
-      hydro::Tags::MassFluxCompute<DataVector, 3, ::Frame::Inertial>>;
+      hydro::Tags::MassFluxCompute<DataVector, 3, ::Frame::Inertial>,
+      gr::Tags::SqrtDetSpatialMetric<DataVector>,
+      gr::Tags::SpatialMetric<DataVector, 3, ::Frame::Inertial>>;
   using non_tensor_compute_tags = tmpl::list<
       tmpl::conditional_t<
           use_dg_subcell,

diff --git a/src/ParallelAlgorithms/Events/CMakeLists.txt b/src/ParallelAlgorithms/Events/CMakeLists.txt
@@ -23,6 +23,7 @@ spectre_target_headers(
   ObserveAtExtremum.hpp
   ObserveFields.hpp
   ObserveInterpolatedData.hpp
+  ObserveInterpolatedIntegralData.hpp
   ObserveNorms.hpp
   ObserveTimeStep.hpp
   Tags.hpp

diff --git a/src/ParallelAlgorithms/Events/Factory.hpp b/src/ParallelAlgorithms/Events/Factory.hpp
@@ -10,6 +10,7 @@
 #include "ParallelAlgorithms/Events/ObserveAdaptiveSteppingDiagnostics.hpp"
 #include "ParallelAlgorithms/Events/ObserveFields.hpp"
 #include "ParallelAlgorithms/Events/ObserveInterpolatedData.hpp"
+#include "ParallelAlgorithms/Events/ObserveInterpolatedIntegralData.hpp"
 #include "ParallelAlgorithms/Events/ObserveNorms.hpp"
 #include "ParallelAlgorithms/Events/ObserveTimeStep.hpp"
 #include "Time/Actions/ChangeSlabSize.hpp"
@@ -24,6 +25,8 @@ using field_observations = tmpl::flatten<tmpl::list<
                   ArraySectionIdTag>,
     ObserveInterpolatedData<VolumeDim, Fields, NonTensorComputeTagsList,
                             ArraySectionIdTag>,
+    ObserveInterpolatedIntegralData<VolumeDim, Fields, NonTensorComputeTagsList,
+                                    ArraySectionIdTag>,
     ::Events::ObserveNorms<Fields, NonTensorComputeTagsList,
                            ArraySectionIdTag>>>;
 }  // namespace dg::Events

diff --git a/src/ParallelAlgorithms/Events/ObserveInterpolatedData.hpp b/src/ParallelAlgorithms/Events/ObserveInterpolatedData.hpp
@@ -233,6 +233,7 @@ class ObserveInterpolatedData<VolumeDim, tmpl::list<Tensors...>,
       for (size_t d = 0; d < VolumeDim; ++d) {
         if (d == interp_dim) {
           gsl::at(target_extents, d) = 1;
+          gsl::at(quadratures, d) = Spectral::Quadrature::Gauss;
         } else {
           gsl::at(target_extents, d) = mesh.extents(d);
           gsl::at(bases, d) = mesh.basis(d);
@@ -244,14 +245,15 @@ class ObserveInterpolatedData<VolumeDim, tmpl::list<Tensors...>,
 
       // Actually need the element_logical_interp_val for next step
       const double elm_interp_val =
-          (interp_val + ((upper_bound_bl - lower_bound_bl) / 2.)) /
-          ((upper_bound_bl + lower_bound_bl) / 2.);
+          (interp_val - ((upper_bound_bl + lower_bound_bl) / 2.)) /
+          ((upper_bound_bl - lower_bound_bl) / 2.);
 
       auto override_target_mesh_with_1d_coords =
           make_array<VolumeDim>(DataVector{});
-      gsl::at(override_target_mesh_with_1d_coords, interp_dim) = elm_interp_val;
+      gsl::at(override_target_mesh_with_1d_coords, interp_dim) =
+          DataVector{elm_interp_val};
       const intrp::RegularGrid<VolumeDim> interpolant(
-          mesh, new_mesh, override_target_mesh_with_1d_coords);
+          mesh, mesh, override_target_mesh_with_1d_coords);
 
       // // Remove tensor types, only storing individual components.
       std::vector<TensorComponent> components;
@@ -287,6 +289,7 @@ class ObserveInterpolatedData<VolumeDim, tmpl::list<Tensors...>,
             const std::string tag_name = db::tag_name<tensor_tag>();
             if (const auto var_to_observe = variables_to_observe.find(tag_name);
                 var_to_observe != variables_to_observe.end()) {
+
               const auto& tensor = get<tensor_tag>(box);
               if (not has_value(tensor)) {
                 // This will only print a warning the first time it's called on