SpECTRE Documentation Coverage Report
Current view: top level - Evolution/DgSubcell/Events - ObserveFields.hpp Hit Total Coverage
Commit: b5ffa4904430ccef0b226f73dcd25c74cb5188f6 Lines: 2 21 9.5 %
Date: 2021-07-28 22:05:18
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <functional>
       8             : #include <initializer_list>
       9             : #include <optional>
      10             : #include <pup.h>
      11             : #include <string>
      12             : #include <type_traits>
      13             : #include <unordered_set>
      14             : #include <utility>
      15             : #include <vector>
      16             : 
      17             : #include "DataStructures/DataBox/Prefixes.hpp"
      18             : #include "DataStructures/DataBox/TagName.hpp"
      19             : #include "DataStructures/DataVector.hpp"
      20             : #include "DataStructures/Tensor/TensorData.hpp"
      21             : #include "Domain/CoordinateMaps/CoordinateMap.hpp"
      22             : #include "Domain/CoordinateMaps/Tags.hpp"
      23             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      24             : #include "Domain/FunctionsOfTime/Tags.hpp"
      25             : #include "Domain/Structure/ElementId.hpp"
      26             : #include "Domain/Tags.hpp"
      27             : #include "Evolution/DgSubcell/ActiveGrid.hpp"
      28             : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
      29             : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
      30             : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
      31             : #include "Evolution/TypeTraits.hpp"
      32             : #include "IO/Observer/ArrayComponentId.hpp"
      33             : #include "IO/Observer/ObservationId.hpp"
      34             : #include "IO/Observer/ObserverComponent.hpp"  // IWYU pragma: keep
      35             : #include "IO/Observer/VolumeActions.hpp"      // IWYU pragma: keep
      36             : #include "Options/Auto.hpp"
      37             : #include "Options/Options.hpp"
      38             : #include "Parallel/ArrayIndex.hpp"
      39             : #include "Parallel/CharmPupable.hpp"
      40             : #include "Parallel/GlobalCache.hpp"
      41             : #include "Parallel/Invoke.hpp"
      42             : #include "Parallel/PupStlCpp17.hpp"
      43             : #include "ParallelAlgorithms/Events/ObserveFields.hpp"
      44             : #include "ParallelAlgorithms/EventsAndTriggers/Event.hpp"
      45             : #include "PointwiseFunctions/AnalyticSolutions/Tags.hpp"
      46             : #include "Utilities/Algorithm.hpp"
      47             : #include "Utilities/Literals.hpp"
      48             : #include "Utilities/MakeString.hpp"
      49             : #include "Utilities/Numeric.hpp"
      50             : #include "Utilities/StdHelpers.hpp"
      51             : #include "Utilities/TMPL.hpp"
      52             : #include "Utilities/TypeTraits/IsA.hpp"
      53             : 
      54             : /// \cond
      55             : template <size_t Dim>
      56             : class Mesh;
      57             : namespace Frame {
      58             : struct Inertial;
      59             : }  // namespace Frame
      60             : /// \endcond
      61             : 
      62           0 : namespace evolution::dg::subcell::Events {
      63             : /// \cond
      64             : template <size_t VolumeDim, typename ObservationValueTag, typename Tensors,
      65             :           typename AnalyticSolutionTensors = tmpl::list<>,
      66             :           typename NonSolutionTensors =
      67             :               tmpl::list_difference<Tensors, AnalyticSolutionTensors>>
      68             : class ObserveFields;
      69             : /// \cond
      70             : 
      71             : /*!
      72             :  * \ingroup DiscontinuousGalerkinGroup
      73             :  * \brief %Observe volume tensor fields.
      74             :  *
      75             :  * A class that writes volume quantities to an h5 file during the simulation.
      76             :  * The observed quantitites are:
      77             :  * - `InertialCoordinates`
      78             :  * - Tensors listed in `Tensors` template parameter
      79             :  * - `Error(*)` = errors in `AnalyticSolutionTensors` =
      80             :  *   \f$\text{value} - \text{analytic solution}\f$
      81             :  *
      82             :  * The user may specify an `interpolation_mesh` to which the
      83             :  * data is interpolated.
      84             :  */
      85             : template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
      86             :           typename... AnalyticSolutionTensors, typename... NonSolutionTensors>
      87             : class ObserveFields<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
      88             :                     tmpl::list<AnalyticSolutionTensors...>,
      89             :                     tmpl::list<NonSolutionTensors...>> : public Event {
      90             :  private:
      91             :   static_assert(
      92             :       std::is_same_v<
      93             :           tmpl::list_difference<tmpl::list<AnalyticSolutionTensors...>,
      94             :                                 tmpl::list<Tensors...>>,
      95             :           tmpl::list<>>,
      96             :       "All AnalyticSolutionTensors must be listed in Tensors.");
      97             : 
      98             :   using dg_observe_fields =
      99             :       ::dg::Events::ObserveFields<VolumeDim, ObservationValueTag,
     100             :                                   tmpl::list<Tensors...>,
     101             :                                   tmpl::list<AnalyticSolutionTensors...>,
     102             :                                   tmpl::list<NonSolutionTensors...>>;
     103             : 
     104             :  public:
     105             :   /// The name of the subfile inside the HDF5 file
     106             :   struct SubfileName {
     107             :     using type = std::string;
     108             :     static constexpr Options::String help = {
     109             :         "The name of the subfile inside the HDF5 file without an extension and "
     110             :         "without a preceding '/'."};
     111             :   };
     112             : 
     113             :   /// \cond
     114             :   explicit ObserveFields(CkMigrateMessage* /*unused*/) noexcept {}
     115             :   using PUP::able::register_constructor;
     116             :   WRAPPED_PUPable_decl_template(ObserveFields);  // NOLINT
     117             :   /// \endcond
     118             : 
     119           0 :   using VariablesToObserve = typename dg_observe_fields::VariablesToObserve;
     120           0 :   using InterpolateToMesh = typename dg_observe_fields::InterpolateToMesh;
     121             : 
     122             :   /// The floating point type/precision with which to write the data to disk.
     123             :   ///
     124             :   /// Must be specified once for all data or individually for each variable
     125             :   /// being observed.
     126           1 :   using FloatingPointTypes = typename dg_observe_fields::FloatingPointTypes;
     127             : 
     128             :   /// The floating point type/precision with which to write the coordinates to
     129             :   /// disk.
     130           1 :   using CoordinatesFloatingPointType =
     131             :       typename dg_observe_fields::CoordinatesFloatingPointType;
     132             : 
     133           0 :   using options =
     134             :       tmpl::list<SubfileName, CoordinatesFloatingPointType, FloatingPointTypes,
     135             :                  VariablesToObserve, InterpolateToMesh>;
     136             : 
     137           0 :   static constexpr Options::String help =
     138             :       "Observe volume tensor fields.\n"
     139             :       "\n"
     140             :       "Writes volume quantities:\n"
     141             :       " * InertialCoordinates\n"
     142             :       " * Tensors listed in Tensors template parameter\n"
     143             :       " * Error(*) = errors in AnalyticSolutionTensors\n"
     144             :       "            = value - analytic solution\n"
     145             :       "\n";
     146             : 
     147           0 :   ObserveFields() = default;
     148             : 
     149           0 :   ObserveFields(const std::string& subfile_name,
     150             :                 FloatingPointType coordinates_floating_point_type,
     151             :                 const std::vector<FloatingPointType>& floating_point_types,
     152             :                 const std::vector<std::string>& variables_to_observe,
     153             :                 std::optional<Mesh<VolumeDim>> interpolation_mesh = {},
     154             :                 const Options::Context& context = {});
     155             : 
     156           0 :   using coordinates_tag =
     157             :       ::domain::Tags::Coordinates<VolumeDim, Frame::Inertial>;
     158             : 
     159           0 :   using argument_tags =
     160             :       tmpl::list<ObservationValueTag, ::domain::Tags::Mesh<VolumeDim>,
     161             :                  subcell::Tags::Mesh<VolumeDim>, subcell::Tags::ActiveGrid,
     162             :                  ::domain::Tags::Coordinates<VolumeDim, Frame::Grid>,
     163             :                  subcell::Tags::Coordinates<VolumeDim, Frame::Grid>,
     164             :                  ::domain::CoordinateMaps::Tags::CoordinateMap<
     165             :                      VolumeDim, Frame::Grid, Frame::Inertial>,
     166             :                  ::domain::Tags::FunctionsOfTime, AnalyticSolutionTensors...,
     167             :                  NonSolutionTensors...>;
     168             : 
     169             :   template <typename Metavariables, typename ParallelComponent>
     170           0 :   void operator()(
     171             :       const double observation_value, const Mesh<VolumeDim>& dg_mesh,
     172             :       const Mesh<VolumeDim>& subcell_mesh,
     173             :       const subcell::ActiveGrid active_grid,
     174             :       const tnsr::I<DataVector, VolumeDim, Frame::Grid>& dg_grid_coordinates,
     175             :       const tnsr::I<DataVector, VolumeDim, Frame::Grid>&
     176             :           subcell_grid_coordinates,
     177             :       const ::domain::CoordinateMapBase<Frame::Grid, Frame::Inertial,
     178             :                                         VolumeDim>& grid_to_inertial_map,
     179             :       const std::unordered_map<
     180             :           std::string,
     181             :           std::unique_ptr<::domain::FunctionsOfTime::FunctionOfTime>>&
     182             :           functions_of_time,
     183             :       const typename AnalyticSolutionTensors::
     184             :           type&... analytic_solution_tensors,
     185             :       const typename NonSolutionTensors::type&... non_solution_tensors,
     186             :       Parallel::GlobalCache<Metavariables>& cache,
     187             :       const ElementId<VolumeDim>& array_index,
     188             :       const ParallelComponent* const component) const noexcept {
     189             :     std::optional<
     190             :         Variables<tmpl::list<::Tags::Analytic<AnalyticSolutionTensors>...>>>
     191             :         analytic_solution_variables{};
     192             :     const auto set_analytic_soln =
     193             :         [&analytic_solution_variables, &cache, &observation_value](
     194             :             const Mesh<VolumeDim>& mesh,
     195             :             const tnsr::I<DataVector, VolumeDim, Frame::Inertial>&
     196             :                 inertial_coords) noexcept {
     197             :           if constexpr (evolution::is_analytic_solution_v<
     198             :                             typename Metavariables::initial_data>) {
     199             :             Variables<tmpl::list<AnalyticSolutionTensors...>> soln_vars{
     200             :                 mesh.number_of_grid_points()};
     201             :             soln_vars.assign_subset(
     202             :                 Parallel::get<::Tags::AnalyticSolutionBase>(cache).variables(
     203             :                     inertial_coords, observation_value,
     204             :                     tmpl::list<AnalyticSolutionTensors...>{}));
     205             :             analytic_solution_variables = std::move(soln_vars);
     206             :           } else {
     207             :             (void)analytic_solution_variables;
     208             :             (void)cache;
     209             :             (void)observation_value;
     210             :           }
     211             :         };
     212             :     if (active_grid == subcell::ActiveGrid::Dg) {
     213             :       const auto dg_inertial_coords = grid_to_inertial_map(
     214             :           dg_grid_coordinates, observation_value, functions_of_time);
     215             :       set_analytic_soln(dg_mesh, dg_inertial_coords);
     216             :       ::dg::Events::ObserveFields<VolumeDim, ObservationValueTag,
     217             :                                   tmpl::list<Tensors...>,
     218             :                                   tmpl::list<AnalyticSolutionTensors...>,
     219             :                                   tmpl::list<NonSolutionTensors...>>::
     220             :           call_operator_impl(
     221             :               subfile_path_, variables_to_observe_, interpolation_mesh_,
     222             :               observation_value, dg_mesh, dg_inertial_coords,
     223             :               analytic_solution_tensors..., non_solution_tensors...,
     224             :               analytic_solution_variables, cache, array_index, component);
     225             :     } else {
     226             :       ASSERT(active_grid == subcell::ActiveGrid::Subcell,
     227             :              "Active grid must be either Dg or Subcell");
     228             :       const auto subcell_inertial_coords = grid_to_inertial_map(
     229             :           subcell_grid_coordinates, observation_value, functions_of_time);
     230             :       set_analytic_soln(subcell_mesh, subcell_inertial_coords);
     231             :       dg_observe_fields::call_operator_impl(
     232             :           subfile_path_, variables_to_observe_, interpolation_mesh_,
     233             :           observation_value, subcell_mesh, subcell_inertial_coords,
     234             :           analytic_solution_tensors..., non_solution_tensors...,
     235             :           analytic_solution_variables, cache, array_index, component);
     236             :     }
     237             :   }
     238             : 
     239             :   // This overload is called when the list of analytic-solution tensors is
     240             :   // empty, i.e. it is clear at compile-time that no analytic solutions are
     241             :   // available
     242             :   template <typename Metavariables, typename ParallelComponent>
     243           0 :   void operator()(
     244             :       const double observation_value, const Mesh<VolumeDim>& mesh,
     245             :       const tnsr::I<DataVector, VolumeDim, Frame::Inertial>&
     246             :           inertial_coordinates,
     247             :       const typename NonSolutionTensors::type&... non_solution_tensors,
     248             :       Parallel::GlobalCache<Metavariables>& cache,
     249             :       const ElementId<VolumeDim>& array_index,
     250             :       const ParallelComponent* const component) const noexcept {
     251             :     this->operator()(observation_value, mesh, inertial_coordinates,
     252             :                      non_solution_tensors..., std::nullopt, cache, array_index,
     253             :                      component);
     254             :   }
     255             : 
     256           0 :   using observation_registration_tags = tmpl::list<>;
     257             :   std::pair<observers::TypeOfObservation, observers::ObservationKey>
     258           0 :   get_observation_type_and_key_for_registration() const noexcept {
     259             :     return {observers::TypeOfObservation::Volume,
     260             :             observers::ObservationKey(subfile_path_ + ".vol")};
     261             :   }
     262             : 
     263             :   // NOLINTNEXTLINE(google-runtime-references)
     264           0 :   void pup(PUP::er& p) noexcept override {
     265             :     Event::pup(p);
     266             :     p | subfile_path_;
     267             :     p | variables_to_observe_;
     268             :     p | interpolation_mesh_;
     269             :   }
     270             : 
     271           0 :   bool needs_evolved_variables() const noexcept override { return true; }
     272             : 
     273             :  private:
     274           0 :   std::string subfile_path_;
     275           0 :   std::unordered_map<std::string, FloatingPointType> variables_to_observe_{};
     276           0 :   std::optional<Mesh<VolumeDim>> interpolation_mesh_{};
     277             : };
     278             : 
     279             : /// \cond
     280             : template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
     281             :           typename... AnalyticSolutionTensors, typename... NonSolutionTensors>
     282             : ObserveFields<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
     283             :               tmpl::list<AnalyticSolutionTensors...>,
     284             :               tmpl::list<NonSolutionTensors...>>::
     285             :     ObserveFields(const std::string& subfile_name,
     286             :                   const FloatingPointType coordinates_floating_point_type,
     287             :                   const std::vector<FloatingPointType>& floating_point_types,
     288             :                   const std::vector<std::string>& variables_to_observe,
     289             :                   std::optional<Mesh<VolumeDim>> interpolation_mesh,
     290             :                   const Options::Context& context)
     291             :     : subfile_path_("/" + subfile_name),
     292             :       variables_to_observe_([&context, &floating_point_types,
     293             :                              &variables_to_observe]() {
     294             :         if (floating_point_types.size() != 1 and
     295             :             floating_point_types.size() != variables_to_observe.size()) {
     296             :           PARSE_ERROR(context, "The number of floating point types specified ("
     297             :                                    << floating_point_types.size()
     298             :                                    << ") must be 1 or the number of variables "
     299             :                                       "specified for observing ("
     300             :                                    << variables_to_observe.size() << ")");
     301             :         }
     302             :         std::unordered_map<std::string, FloatingPointType> result{};
     303             :         for (size_t i = 0; i < variables_to_observe.size(); ++i) {
     304             :           result[variables_to_observe[i]] = floating_point_types.size() == 1
     305             :                                                 ? floating_point_types[0]
     306             :                                                 : floating_point_types[i];
     307             :           ASSERT(
     308             :               result.at(variables_to_observe[i]) == FloatingPointType::Float or
     309             :                   result.at(variables_to_observe[i]) ==
     310             :                       FloatingPointType::Double,
     311             :               "Floating point type for variable '"
     312             :                   << variables_to_observe[i]
     313             :                   << "' must be either Float or Double.");
     314             :         }
     315             :         return result;
     316             :       }()),
     317             :       interpolation_mesh_(interpolation_mesh) {
     318             :   using ::operator<<;
     319             :   const std::unordered_set<std::string> valid_tensors{
     320             :       db::tag_name<Tensors>()...};
     321             :   for (const auto& [name, floating_point_type] : variables_to_observe_) {
     322             :     (void)floating_point_type;
     323             :     if (valid_tensors.count(name) != 1) {
     324             :       PARSE_ERROR(
     325             :           context,
     326             :           name << " is not an available variable.  Available variables:\n"
     327             :                << (std::vector<std::string>{db::tag_name<Tensors>()...}));
     328             :     }
     329             :     if (alg::count(variables_to_observe, name) != 1) {
     330             :       PARSE_ERROR(context, name << " specified multiple times");
     331             :     }
     332             :   }
     333             :   variables_to_observe_[coordinates_tag::name()] =
     334             :       coordinates_floating_point_type;
     335             :   if (interpolation_mesh.has_value()) {
     336             :     PARSE_ERROR(
     337             :         context,
     338             :         "We don't yet support interpolating to a mesh with DG-subcell because "
     339             :         "we don't have an interpolator on the subcell grid. Trilinear "
     340             :         "interpolation could be added to support interpolation.");
     341             :   }
     342             : }
     343             : /// \endcond
     344             : 
     345             : /// \cond
     346             : template <size_t VolumeDim, typename ObservationValueTag, typename... Tensors,
     347             :           typename... AnalyticSolutionTensors, typename... NonSolutionTensors>
     348             : PUP::able::PUP_ID
     349             :     ObserveFields<VolumeDim, ObservationValueTag, tmpl::list<Tensors...>,
     350             :                   tmpl::list<AnalyticSolutionTensors...>,
     351             :                   tmpl::list<NonSolutionTensors...>>::my_PUP_ID = 0;  // NOLINT
     352             : /// \endcond
     353             : }  // namespace evolution::dg::subcell::Events

Generated by: LCOV version 1.14