diff --git a/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp b/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp index 75ee98f13d0bd..43fc62e910202 100644 --- 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,11 @@ struct EvolutionMetavars, ::Events::Tags::ObserverCoordinates, hydro::Tags::TransportVelocity>, - hydro::Tags::InversePlasmaBetaCompute>; + hydro::Tags::InversePlasmaBetaCompute, + hydro::Tags::MassFluxCompute, + gr::Tags::SqrtDetSpatialMetric, + gr::Tags::SpatialMetric, + gr::Tags::Lapse>; using non_tensor_compute_tags = tmpl::list< tmpl::conditional_t< use_dg_subcell, diff --git a/src/IO/H5/VolumeData.cpp b/src/IO/H5/VolumeData.cpp index 9bd9ac44254a9..6ff10661846dd 100644 --- a/src/IO/H5/VolumeData.cpp +++ b/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,7 +68,18 @@ void append_element_extents_and_connectivity( // out size without computing all the connectivities. const std::vector connectivity = [&extents, &total_points_so_far]() { std::vector local_connectivity; - for (const auto& cell : vis::detail::compute_cells(extents)) { + std::vector 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(bounding_indices)); @@ -78,6 +87,7 @@ void append_element_extents_and_connectivity( } 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& elements, + const std::vector& elements2, const std::optional>& serialized_domain, const std::optional>& 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 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 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 = diff --git a/src/ParallelAlgorithms/Events/CMakeLists.txt b/src/ParallelAlgorithms/Events/CMakeLists.txt index 5a6dfacb696cb..4089bee86a996 100644 --- a/src/ParallelAlgorithms/Events/CMakeLists.txt +++ b/src/ParallelAlgorithms/Events/CMakeLists.txt @@ -24,6 +24,8 @@ spectre_target_headers( ObserveDataBox.hpp 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 index 16f4dc725af1b..6b653acaabc84 100644 --- a/src/ParallelAlgorithms/Events/Factory.hpp +++ b/src/ParallelAlgorithms/Events/Factory.hpp @@ -9,6 +9,8 @@ #include "ParallelAlgorithms/Events/ErrorIfDataTooBig.hpp" #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/ChangeSlabSize/Event.hpp" @@ -21,6 +23,10 @@ using field_observations = tmpl::flatten, ObserveFields, + ObserveInterpolatedData, + ObserveInterpolatedIntegralData, ::Events::ObserveNorms>>; } // namespace dg::Events diff --git a/src/ParallelAlgorithms/Events/ObserveInterpolatedData.hpp b/src/ParallelAlgorithms/Events/ObserveInterpolatedData.hpp new file mode 100644 index 0000000000000..a7df2c762d6f8 --- /dev/null +++ b/src/ParallelAlgorithms/Events/ObserveInterpolatedData.hpp @@ -0,0 +1,476 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/ObservationBox.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataBox/TagName.hpp" +#include "DataStructures/DataBox/ValidateSelection.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/FloatingPointType.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "IO/H5/TensorData.hpp" +#include "IO/Observer/GetSectionObservationKey.hpp" +#include "IO/Observer/ObservationId.hpp" +#include "IO/Observer/ObserverComponent.hpp" +#include "IO/Observer/Tags.hpp" +#include "IO/Observer/VolumeActions.hpp" +#include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Options/Auto.hpp" +#include "Options/String.hpp" +#include "Parallel/ArrayComponentId.hpp" +#include "Parallel/ArrayIndex.hpp" +#include "Parallel/GlobalCache.hpp" +#include "Parallel/Invoke.hpp" +#include "Parallel/Local.hpp" +#include "Parallel/Reduction.hpp" +#include "ParallelAlgorithms/Events/Tags.hpp" +#include "ParallelAlgorithms/EventsAndTriggers/Event.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Tags.hpp" +#include "Utilities/Algorithm.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/Literals.hpp" +#include "Utilities/MakeString.hpp" +#include "Utilities/Numeric.hpp" +#include "Utilities/OptionalHelpers.hpp" +#include "Utilities/Serialization/CharmPupable.hpp" +#include "Utilities/Serialization/PupStlCpp17.hpp" +#include "Utilities/StdHelpers.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TypeTraits/IsA.hpp" + +#include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp" +#include "Utilities/ErrorHandling/CaptureForError.hpp" //FIXME: + +/// \cond +template +class Mesh; +namespace Frame { +struct Inertial; +} // namespace Frame +/// \endcond + +// FIXME: got clang-tidy "nested namespace can be concatenated" +namespace dg::Events { +/// \cond +template , + typename ArraySectionIdTag = void> +class ObserveInterpolatedData; +/// \endcond + +/*! + * \ingroup DiscontinuousGalerkinGroup + * \brief %Observe volume tensor fields interpolated to a new mesh slice. + * + * A class that writes volume quantities to an h5 file during the simulation. + * The observed quantitites are specified in the `VariablesToObserve` option. + * Any `Tensor` in the `db::DataBox` can be observed but must be listed in the + * `Tensors` template parameter. Any additional compute tags that hold a + * `Tensor` can also be added to the `Tensors` template parameter. Finally, + * `Variables` and other non-tensor compute tags can be listed in the + * `NonTensorComputeTags` to facilitate observing. Note that the + * `InertialCoordinates` are always observed. + * + * FIXME: add more details here. + * + * + * \note The `NonTensorComputeTags` are intended to be used for `Variables` + * compute tags like `Tags::DerivCompute` + * + * \par Array sections + * This event supports sections (see `Parallel::Section`). Set the + * `ArraySectionIdTag` template parameter to split up observations into subsets + * of elements. The `observers::Tags::ObservationKey` must be + * available in the DataBox. It identifies the section and is used as a suffix + * for the path in the output file. + */ +template +class ObserveInterpolatedData, + tmpl::list, + ArraySectionIdTag> : public Event { + public: + /// The name of the subfile inside the HDF5 file + struct SubfileName { + using type = std::string; + static constexpr Options::String help = { + "The name of the subfile inside the HDF5 file without an extension and " + "without a preceding '/'."}; + }; + + /// \cond + explicit ObserveInterpolatedData(CkMigrateMessage* /*unused*/) {} + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(ObserveInterpolatedData); // NOLINT + /// \endcond + + struct VariablesToObserve { + static constexpr Options::String help = "Subset of variables to observe"; + using type = std::vector; + static size_t lower_bound_on_size() { return 1; } + }; + + struct InterpDim { + static constexpr Options::String help = "Dimension to take the interpolate"; + using type = size_t; + }; + + struct InterpVal { + static constexpr Options::String help = + "Value at which to do interpolation"; + using type = double; + }; + + /// The floating point type/precision with which to write the data to disk. + /// + /// Must be specified once for all data or individually for each variable + /// being observed. + struct FloatingPointTypes { + static constexpr Options::String help = + "The floating point type/precision with which to write the data to " + "disk.\n\n" + "Must be specified once for all data or individually for each " + "variable being observed."; + using type = std::vector; + static size_t upper_bound_on_size() { return sizeof...(Tensors); } + static size_t lower_bound_on_size() { return 1; } + }; + + /// The floating point type/precision with which to write the coordinates to + /// disk. + struct CoordinatesFloatingPointType { + static constexpr Options::String help = + "The floating point type/precision with which to write the coordinates " + "to disk."; + using type = FloatingPointType; + }; + + using options = + tmpl::list; + + static constexpr Options::String help = + "Observe volume tensor fields.\n" + "\n" + "Writes volume quantities:\n" + " * InertialCoordinates\n" + " * Tensors listed in the 'VariablesToObserve' option\n"; + + ObserveInterpolatedData() = default; + + ObserveInterpolatedData( + const std::string& subfile_name, + FloatingPointType coordinates_floating_point_type, + const std::vector& floating_point_types, + const std::vector& variables_to_observe, + const size_t interp_dim, const double interp_val, + const Options::Context& context = {}); + + using compute_tags_for_observation_box = + tmpl::list; + + using return_tags = tmpl::list<>; + using argument_tags = tmpl::list<::Tags::ObservationBox, + ::Events::Tags::ObserverMesh>; + + template + void operator()(const ObservationBox& box, + const Mesh& mesh, + Parallel::GlobalCache& cache, + const ElementId& array_index, + const ParallelComponent* const component, + const ObservationValue& observation_value) const { + // Skip observation on elements that are not part of a section + const std::optional section_observation_key = + observers::get_section_observation_key(box); + if (not section_observation_key.has_value()) { + return; + } + call_operator_impl(subfile_path_ + *section_observation_key, + variables_to_observe_, interp_dim_, interp_val_, mesh, + box, cache, array_index, component, observation_value); + } + + // We factor out the work into a static member function so it can be shared + // with other field observing events, like the one that deals with DG-subcell + // where there are two grids. This is to avoid copy-pasting all of the code. + template + static void call_operator_impl( + const std::string& subfile_path, + const std::unordered_map& + variables_to_observe, + const size_t interp_dim, const double interp_val, + const Mesh& mesh, + const ObservationBox& box, + Parallel::GlobalCache& cache, + const ElementId& element_id, + const ParallelComponent* const /*meta*/, + const ObservationValue& observation_value) { + // For each element, we find the block logical coordinate + // of the upper and lower side in the interpo_dim. + // We only need to use the element for which the interp_val + // falls within the range. + + const double lower_bound_bl = + element_id.segment_id(interp_dim).endpoint(Side::Lower); + const double upper_bound_bl = + element_id.segment_id(interp_dim).endpoint(Side::Upper); + // if ((lower_bound_bl < interp_val) and (upper_bound_bl >= interp_val)) { + if ((lower_bound_bl <= interp_val) and (upper_bound_bl > interp_val)) { + // Create a 3 dimensional new_mesh corresponding to a single point in + // interp_dim times 2-dimensional slice of the original Mesh. + // This is used for the sake of creating a RegularGrid Interpolant. + std::array target_extents{}; + std::array bases{}; + std::array quadratures{}; + for (size_t d = 0; d < VolumeDim; ++d) { + gsl::at(bases, d) = mesh.basis(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(quadratures, d) = mesh.quadrature(d); + } + } + const auto new_mesh = + Mesh{{target_extents}, {bases}, {quadratures}}; + + // We need to convert the block logical coordinates + // to corresponding element logical coordinate for the given + // element for the next step. + const double elm_interp_val = + (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(DataVector{}); + // gsl::at(override_target_mesh_with_1d_coords, interp_dim) = + // DataVector{elm_interp_val}; + // const intrp::RegularGrid interpolant( + // mesh, mesh, override_target_mesh_with_1d_coords); + auto target_points = data_on_slice(logical_coordinates(mesh), + mesh.extents(), interp_dim, 0); + target_points.get(interp_dim) = elm_interp_val; + + const intrp::Irregular interpolant(mesh, target_points); + + // // Remove tensor types, only storing individual components. + std::vector components; + // // This is larger than we need if we are only observing some + // // tensors, but that's not a big deal and calculating the correct + // // size is nontrivial. + components.reserve(alg::accumulate( + std::initializer_list{std::decay_t::size()...}, + 0_st)); + + const auto record_tensor_component_impl = + [&components, &interpolant, &elm_interp_val, &mesh]( + const auto& tensor, const FloatingPointType floating_point_type, + const std::string& tag_name) { + for (size_t i = 0; i < tensor.size(); ++i) { + const auto tensor_component = interpolant.interpolate(tensor[i]); + // const auto wow = logical_coordinates(mesh); + // CAPTURE_FOR_ERROR(tensor_component); + // CAPTURE_FOR_ERROR(tensor[i]); + // CAPTURE_FOR_ERROR(elm_interp_val); + // CAPTURE_FOR_ERROR(wow); + + if (tag_name == "RestMassDensity" and min(tensor_component) < 0) { + ERROR("negative rho"); + } + if (floating_point_type == FloatingPointType::Float) { + components.emplace_back( + tag_name + tensor.component_suffix(i), + std::vector{tensor_component.begin(), + tensor_component.end()}); + } else { + components.emplace_back(tag_name + tensor.component_suffix(i), + tensor_component); + } + } + }; + const auto record_tensor_components = + [&box, &record_tensor_component_impl, + &variables_to_observe](const auto tensor_tag_v) { + using tensor_tag = tmpl::type_from; + const std::string tag_name = db::tag_name(); + if (const auto var_to_observe = variables_to_observe.find(tag_name); + var_to_observe != variables_to_observe.end()) { + + const auto& tensor = get(box); + if (not has_value(tensor)) { + // This will only print a warning the first time it's called on + // a node. + [[maybe_unused]] static bool t = + ObserveInterpolatedData::print_warning_about_optional< + tensor_tag>(); + return; + } + const auto floating_point_type = var_to_observe->second; + record_tensor_component_impl(value(tensor), floating_point_type, + tag_name); + } + }; + EXPAND_PACK_LEFT_TO_RIGHT( + record_tensor_components(tmpl::type_{})); + + // // Send data to volume observer + auto& local_observer = *Parallel::local_branch( + Parallel::get_parallel_component>( + cache)); + + Parallel::simple_action( + local_observer, + observers::ObservationId(observation_value.value, + subfile_path + ".vol"), + subfile_path, + Parallel::make_array_component_id(element_id), + ElementVolumeData{element_id, std::move(components), new_mesh}); + } + } + + using observation_registration_tags = tmpl::list<::Tags::DataBox>; + + template + std::optional< + std::pair> + get_observation_type_and_key_for_registration( + const db::DataBox& box) const { + const std::optional section_observation_key = + observers::get_section_observation_key(box); + if (not section_observation_key.has_value()) { + return std::nullopt; + } + if ((db::get>(box) + .id() + .segment_id(interp_dim_) + .endpoint(Side::Lower) < interp_val_) and + (db::get>(box) + .id() + .segment_id(interp_dim_) + .endpoint(Side::Upper)) >= interp_val_) { + return {{observers::TypeOfObservation::Volume, + observers::ObservationKey( + subfile_path_ + section_observation_key.value() + ".vol")}}; + } else { + return std::nullopt; + } + } + + using is_ready_argument_tags = tmpl::list<>; + + template + bool is_ready(Parallel::GlobalCache& /*cache*/, + const ArrayIndex& /*array_index*/, + const Component* const /*meta*/) const { + return true; + } + + bool needs_evolved_variables() const override { return true; } + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& p) override { + Event::pup(p); + p | subfile_path_; + p | variables_to_observe_; + p | interp_dim_; + p | interp_val_; + } + + private: + template + static bool print_warning_about_optional() { + Parallel::printf( + "Warning: ObserveInterpolatedData is trying to dump the tag %s " + "but it is stored as a std::optional and has not been " + "evaluated. This most commonly occurs when you are " + "trying to either observe an analytic solution or errors when " + "no analytic solution is available.\n", + db::tag_name()); + return false; + } + + std::string subfile_path_; + std::unordered_map variables_to_observe_{}; + size_t interp_dim_{}; + double interp_val_{}; +}; + +template +ObserveInterpolatedData, + tmpl::list, + ArraySectionIdTag>:: + ObserveInterpolatedData( + const std::string& subfile_name, + const FloatingPointType coordinates_floating_point_type, + const std::vector& floating_point_types, + const std::vector& variables_to_observe, + const size_t interp_dim, const double interp_val, + const Options::Context& context) + : subfile_path_("/" + subfile_name), + variables_to_observe_([&context, &floating_point_types, + &variables_to_observe]() { + if (floating_point_types.size() != 1 and + floating_point_types.size() != variables_to_observe.size()) { + PARSE_ERROR(context, "The number of floating point types specified (" + << floating_point_types.size() + << ") must be 1 or the number of variables " + "specified for observing (" + << variables_to_observe.size() << ")"); + } + std::unordered_map result{}; + for (size_t i = 0; i < variables_to_observe.size(); ++i) { + result[variables_to_observe[i]] = floating_point_types.size() == 1 + ? floating_point_types[0] + : floating_point_types[i]; + ASSERT( + result.at(variables_to_observe[i]) == FloatingPointType::Float or + result.at(variables_to_observe[i]) == + FloatingPointType::Double, + "Floating point type for variable '" + << variables_to_observe[i] + << "' must be either Float or Double."); + } + return result; + }()), + interp_dim_(interp_dim), + interp_val_(interp_val) { + ASSERT( + (... or (db::tag_name() == "InertialCoordinates")), + "There is no tag with name 'InertialCoordinates' specified " + "for the observer. Please make sure you specify a tag in the 'Tensors' " + "list that has the 'db::tag_name()' 'InertialCoordinates'."); + db::validate_selection>(variables_to_observe, context); + variables_to_observe_["InertialCoordinates"] = + coordinates_floating_point_type; +} + +/// \cond +template +PUP::able::PUP_ID ObserveInterpolatedData, + tmpl::list, + ArraySectionIdTag>::my_PUP_ID = + 0; // NOLINT +/// \endcond +} // namespace dg::Events diff --git a/src/ParallelAlgorithms/Events/ObserveInterpolatedIntegralData.hpp b/src/ParallelAlgorithms/Events/ObserveInterpolatedIntegralData.hpp new file mode 100644 index 0000000000000..02903355b8d10 --- /dev/null +++ b/src/ParallelAlgorithms/Events/ObserveInterpolatedIntegralData.hpp @@ -0,0 +1,465 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/ObservationBox.hpp" +#include "DataStructures/DataBox/Prefixes.hpp" +#include "DataStructures/DataBox/TagName.hpp" +#include "DataStructures/DataBox/ValidateSelection.hpp" +#include "DataStructures/DataVector.hpp" +#include "DataStructures/FloatingPointType.hpp" +#include "Domain/Structure/ElementId.hpp" +#include "Domain/Tags.hpp" +#include "IO/H5/TensorData.hpp" +#include "IO/Observer/GetSectionObservationKey.hpp" +#include "IO/Observer/ObservationId.hpp" +#include "IO/Observer/ObserverComponent.hpp" +#include "IO/Observer/ReductionActions.hpp" +#include "IO/Observer/Tags.hpp" +#include "IO/Observer/VolumeActions.hpp" +#include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp" +#include "NumericalAlgorithms/Interpolation/RegularGridInterpolant.hpp" +#include "NumericalAlgorithms/LinearOperators/DefiniteIntegral.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Options/Auto.hpp" +#include "Options/String.hpp" +#include "Parallel/ArrayComponentId.hpp" +#include "Parallel/ArrayIndex.hpp" +#include "Parallel/GlobalCache.hpp" +#include "Parallel/Invoke.hpp" +#include "Parallel/Local.hpp" +#include "Parallel/Reduction.hpp" +#include "ParallelAlgorithms/Events/Tags.hpp" +#include "ParallelAlgorithms/EventsAndTriggers/Event.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Tags.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "Utilities/Algorithm.hpp" +#include "Utilities/ErrorHandling/Assert.hpp" +#include "Utilities/Functional.hpp" +#include "Utilities/Literals.hpp" +#include "Utilities/MakeString.hpp" +#include "Utilities/Numeric.hpp" +#include "Utilities/OptionalHelpers.hpp" +#include "Utilities/Serialization/CharmPupable.hpp" +#include "Utilities/Serialization/PupStlCpp17.hpp" +#include "Utilities/StdHelpers.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TypeTraits/IsA.hpp" + +/// \cond +template +class Mesh; +namespace Frame { +struct Inertial; +} // namespace Frame +/// \endcond + +namespace dg::Events { +namespace detail { +using ObserveInterpolatedReductionData = Parallel::ReductionData< + // Observation value (Time) + Parallel::ReductionDatum>, + // TotalSurfaceArea + Parallel::ReductionDatum>, + // MagneticFlux + Parallel::ReductionDatum>, + // SurfaceIntegral of the given tag + Parallel::ReductionDatum>>; +} // namespace detail +/// \cond +template , + typename ArraySectionIdTag = void> +class ObserveInterpolatedIntegralData; +/// \endcond + +/*! + * \ingroup DiscontinuousGalerkinGroup + * \brief %Observe volume tensor fields interpolated to a new mesh slice and + * the corresponding surface integral. + * + * A class that writes volume quantities to an h5 file during the simulation. + * The observed quantitites are specified in the `VariablesToObserve` option. + * Any `Tensor` in the `db::DataBox` can be observed but must be listed in the + * `Tensors` template parameter. Any additional compute tags that hold a + * `Tensor` can also be added to the `Tensors` template parameter. Finally, + * `Variables` and other non-tensor compute tags can be listed in the + * `NonTensorComputeTags` to facilitate observing. Note that the + * `InertialCoordinates` are always observed. + * + * FIXME: add more details here. + * + * + * \note The `NonTensorComputeTags` are intended to be used for `Variables` + * compute tags like `Tags::DerivCompute` + * + * \par Array sections + * This event supports sections (see `Parallel::Section`). Set the + * `ArraySectionIdTag` template parameter to split up observations into subsets + * of elements. The `observers::Tags::ObservationKey` must be + * available in the DataBox. It identifies the section and is used as a suffix + * for the path in the output file. + */ +template +class ObserveInterpolatedIntegralData, + tmpl::list, + ArraySectionIdTag> : public Event { + public: + using ReductionData = Events::detail::ObserveInterpolatedReductionData; + + /// The name of the subfile inside the HDF5 file + struct SubfileName { + using type = std::string; + static constexpr Options::String help = { + "The name of the subfile inside the HDF5 file without an extension and " + "without a preceding '/'."}; + }; + + /// \cond + explicit ObserveInterpolatedIntegralData(CkMigrateMessage* /*unused*/) {} + using PUP::able::register_constructor; + WRAPPED_PUPable_decl_template(ObserveInterpolatedIntegralData); // NOLINT + /// \endcond + + struct VariablesToObserve { + static constexpr Options::String help = "Subset of variables to observe"; + using type = std::vector; + static size_t lower_bound_on_size() { return 1; } + }; + + struct InterpDim { + static constexpr Options::String help = "Dimension to take the interpolate"; + using type = size_t; + }; + + struct InterpVal { + static constexpr Options::String help = + "Value at which to do interpolation"; + using type = double; + }; + + /// The floating point type/precision with which to write the data to disk. + /// + /// Must be specified once for all data or individually for each variable + /// being observed. + struct FloatingPointTypes { + static constexpr Options::String help = + "The floating point type/precision with which to write the data to " + "disk.\n\n" + "Must be specified once for all data or individually for each " + "variable being observed."; + using type = std::vector; + static size_t upper_bound_on_size() { return sizeof...(Tensors); } + static size_t lower_bound_on_size() { return 1; } + }; + + /// The floating point type/precision with which to write the coordinates to + /// disk. + struct CoordinatesFloatingPointType { + static constexpr Options::String help = + "The floating point type/precision with which to write the coordinates " + "to disk."; + using type = FloatingPointType; + }; + + using options = + tmpl::list; + + static constexpr Options::String help = + "Observe volume tensor fields.\n" + "\n" + "Writes volume quantities:\n" + " * InertialCoordinates\n" + " * Tensors listed in the 'VariablesToObserve' option\n"; + + ObserveInterpolatedIntegralData() = default; + + ObserveInterpolatedIntegralData( + const std::string& subfile_name, + FloatingPointType coordinates_floating_point_type, + const std::vector& floating_point_types, + const std::vector& variables_to_observe, + const size_t interp_dim, const double interp_val, + const Options::Context& context = {}); + + using observed_reduction_data_tags = + observers::make_reduction_data_tags>; + + using compute_tags_for_observation_box = + tmpl::list; + + using return_tags = tmpl::list<>; + using argument_tags = tmpl::list<::Tags::ObservationBox, + ::Events::Tags::ObserverMesh>; + + template + void operator()(const ObservationBox& box, + const Mesh& mesh, + Parallel::GlobalCache& cache, + const ElementId& array_index, + const ParallelComponent* const component, + const ObservationValue& observation_value) const { + // Skip observation on elements that are not part of a section + const std::optional section_observation_key = + observers::get_section_observation_key(box); + if (not section_observation_key.has_value()) { + return; + } + call_operator_impl(subfile_path_ + *section_observation_key, + variables_to_observe_, interp_dim_, interp_val_, mesh, + box, cache, array_index, component, observation_value); + } + + // We factor out the work into a static member function so it can be shared + // with other field observing events, like the one that deals with DG-subcell + // where there are two grids. This is to avoid copy-pasting all of the code. + template + static void call_operator_impl( + const std::string& subfile_path, + const std::unordered_map& + variables_to_observe, + const size_t interp_dim, const double interp_val, + const Mesh& mesh, + const ObservationBox& box, + Parallel::GlobalCache& cache, + const ElementId& element_id, + const ParallelComponent* const /*meta*/, + const ObservationValue& observation_value) { + double surface_integral_value = 0.0; + double magnetic_flux = 0.0; + double local_volume = 0.0; + const DataVector det_jacobian = + 1. / + get(get<::Events::Tags::ObserverDetInvJacobian>(box)); + const double lower_bound_bl = + element_id.segment_id(interp_dim).endpoint(Side::Lower); + const double upper_bound_bl = + element_id.segment_id(interp_dim).endpoint(Side::Upper); + if ((lower_bound_bl < interp_val) and (upper_bound_bl >= interp_val)) { + + const auto new_mesh = mesh.slice_away(interp_dim); + // 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.); + + // auto override_target_mesh_with_1d_coords = + // make_array(DataVector{}); + // gsl::at(override_target_mesh_with_1d_coords, interp_dim) = + // DataVector{elm_interp_val}; + // const intrp::RegularGrid interpolant( + // mesh, mesh, override_target_mesh_with_1d_coords); + auto target_points = data_on_slice(logical_coordinates(mesh), + mesh.extents(), interp_dim, 0); + target_points.get(interp_dim) = elm_interp_val; + const intrp::Irregular interpolant(mesh, target_points); + + // this was used to check 4*pi*r^2 + // const auto temp_metric = + // get>(box); + // auto det1 = (get<1, 1>(temp_metric)) * (get<2, 2>(temp_metric)); + // auto det2 = (get<1, 2>(temp_metric)) * (get<2, 1>(temp_metric)); + // auto sqrt_det = sqrt(det1 - det2); + // auto met_rr = get<0, 0>(temp_metric); + + const auto sqrt_det_spatial_metric = + get(get>(box)); + const auto lapse = get(get>(box)); + const auto sqrt_det_metric = (lapse) * (sqrt_det_spatial_metric); + const auto mag_field_r = get<0>( + get>( + box)); + + const auto record_tensor_component_impl = + [&interpolant, &local_volume, &magnetic_flux, &surface_integral_value, + &new_mesh, &det_jacobian, &mag_field_r, + &sqrt_det_metric](const auto& tensor, const std::string& tag_name) { + const auto new_tensor = interpolant.interpolate(tensor[0]); + const auto new_mag_field_r = interpolant.interpolate(mag_field_r); + const auto new_det_jacobian = interpolant.interpolate(det_jacobian); + const auto new_sqrt_det_metric = + interpolant.interpolate(sqrt_det_metric); + // const auto new_sqrt_det = interpolant.interpolate(sqrt_det); + // this one, right now, is the MassFlux first component. + surface_integral_value += + definite_integral(new_tensor * new_det_jacobian, new_mesh); + // surface area of the surface + local_volume += definite_integral( + new_sqrt_det_metric * new_det_jacobian, new_mesh); + // magnetic flux + magnetic_flux += definite_integral( + abs(new_mag_field_r) * new_sqrt_det_metric * new_det_jacobian, + new_mesh); + }; + const auto record_tensor_components = + [&box, &record_tensor_component_impl, + &variables_to_observe](const auto tensor_tag_v) { + using tensor_tag = tmpl::type_from; + const std::string tag_name = db::tag_name(); + if (const auto var_to_observe = variables_to_observe.find(tag_name); + var_to_observe != variables_to_observe.end()) { + const auto& tensor = get(box); + if (not has_value(tensor)) { + // This will only print a warning the first time it's called + // on a node. + [[maybe_unused]] static bool t = + ObserveInterpolatedIntegralData:: + print_warning_about_optional(); + return; + } + const auto floating_point_type = var_to_observe->second; + record_tensor_component_impl(value(tensor), tag_name); + } + }; + EXPAND_PACK_LEFT_TO_RIGHT( + record_tensor_components(tmpl::type_{})); + } + // Send data to volume observer + auto& local_observer = *Parallel::local_branch( + Parallel::get_parallel_component>( + cache)); + Parallel::simple_action( + local_observer, + observers::ObservationId(observation_value.value, + subfile_path + ".dat"), + Parallel::make_array_component_id(element_id), + subfile_path, + std::vector{observation_value.name, "TotalVolume", + "MagneticFlux", "SurfaceIntegral"}, + ReductionData{observation_value.value, std::move(local_volume), + std::move(magnetic_flux), + std::move(surface_integral_value)}); + } + + using observation_registration_tags = tmpl::list<::Tags::DataBox>; + + template + std::optional< + std::pair> + get_observation_type_and_key_for_registration( + const db::DataBox& box) const { + const std::optional section_observation_key = + observers::get_section_observation_key(box); + if (not section_observation_key.has_value()) { + return std::nullopt; + } + return {{observers::TypeOfObservation::Reduction, + observers::ObservationKey( + subfile_path_ + section_observation_key.value() + ".dat")}}; + } + + using is_ready_argument_tags = tmpl::list<>; + + template + bool is_ready(Parallel::GlobalCache& /*cache*/, + const ArrayIndex& /*array_index*/, + const Component* const /*meta*/) const { + return true; + } + + bool needs_evolved_variables() const override { return true; } + + // NOLINTNEXTLINE(google-runtime-references) + void pup(PUP::er& p) override { + Event::pup(p); + p | subfile_path_; + p | variables_to_observe_; + p | interp_dim_; + p | interp_val_; + } + + private: + template + static bool print_warning_about_optional() { + Parallel::printf( + "Warning: ObserveInterpolatedData is trying to dump the tag %s " + "but it is stored as a std::optional and has not been " + "evaluated. This most commonly occurs when you are " + "trying to either observe an analytic solution or errors when " + "no analytic solution is available.\n", + db::tag_name()); + return false; + } + + std::string subfile_path_; + std::unordered_map variables_to_observe_{}; + size_t interp_dim_{}; + double interp_val_{}; +}; + +template +ObserveInterpolatedIntegralData, + tmpl::list, + ArraySectionIdTag>:: + ObserveInterpolatedIntegralData( + const std::string& subfile_name, + const FloatingPointType coordinates_floating_point_type, + const std::vector& floating_point_types, + const std::vector& variables_to_observe, + const size_t interp_dim, const double interp_val, + const Options::Context& context) + : subfile_path_("/" + subfile_name), + variables_to_observe_([&context, &floating_point_types, + &variables_to_observe]() { + if (floating_point_types.size() != 1 and + floating_point_types.size() != variables_to_observe.size()) { + PARSE_ERROR(context, "The number of floating point types specified (" + << floating_point_types.size() + << ") must be 1 or the number of variables " + "specified for observing (" + << variables_to_observe.size() << ")"); + } + std::unordered_map result{}; + for (size_t i = 0; i < variables_to_observe.size(); ++i) { + result[variables_to_observe[i]] = floating_point_types.size() == 1 + ? floating_point_types[0] + : floating_point_types[i]; + ASSERT( + result.at(variables_to_observe[i]) == FloatingPointType::Float or + result.at(variables_to_observe[i]) == + FloatingPointType::Double, + "Floating point type for variable '" + << variables_to_observe[i] + << "' must be either Float or Double."); + } + return result; + }()), + interp_dim_(interp_dim), + interp_val_(interp_val) { + ASSERT( + (... or (db::tag_name() == "InertialCoordinates")), + "There is no tag with name 'InertialCoordinates' specified " + "for the observer. Please make sure you specify a tag in the 'Tensors' " + "list that has the 'db::tag_name()' 'InertialCoordinates'."); + db::validate_selection>(variables_to_observe, context); + // variables_to_observe_["InertialCoordinates"] = + // coordinates_floating_point_type; +} + +/// \cond +template +PUP::able::PUP_ID ObserveInterpolatedIntegralData< + VolumeDim, tmpl::list, tmpl::list, + ArraySectionIdTag>::my_PUP_ID = 0; // NOLINT +/// \endcond +} // namespace dg::Events