From 8d2e388b67fd9b48f073a24e0416b680526f4799 Mon Sep 17 00:00:00 2001 From: Jooheon Yoo Date: Fri, 2 Feb 2024 11:22:24 -0800 Subject: [PATCH] Mass Accretion and other diagnostics --- .../EvolveValenciaDivClean.hpp | 7 +- src/IO/H5/VolumeData.cpp | 32 +- src/ParallelAlgorithms/Events/CMakeLists.txt | 1 + src/ParallelAlgorithms/Events/Factory.hpp | 3 + .../Events/ObserveInterpolatedData.hpp | 476 +++++++++++++++ .../ObserveInterpolatedIntegralData.hpp | 548 ++++++++++++++++++ 6 files changed, 1062 insertions(+), 5 deletions(-) create mode 100644 src/ParallelAlgorithms/Events/ObserveInterpolatedData.hpp create mode 100644 src/ParallelAlgorithms/Events/ObserveInterpolatedIntegralData.hpp diff --git a/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp b/src/Evolution/Executables/GrMhd/ValenciaDivClean/EvolveValenciaDivClean.hpp index 50043314bab4..5d0de9bc3f32 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 9bd9ac44254a..6ff10661846d 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 1d672e17977d..3f1af1b78f3f 100644 --- a/src/ParallelAlgorithms/Events/CMakeLists.txt +++ b/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 diff --git a/src/ParallelAlgorithms/Events/Factory.hpp b/src/ParallelAlgorithms/Events/Factory.hpp index 16f4dc725af1..d7da5b86d85d 100644 --- a/src/ParallelAlgorithms/Events/Factory.hpp +++ b/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, ObserveFields, + 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 000000000000..a7df2c762d6f --- /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 000000000000..3ab8a5670d76 --- /dev/null +++ b/src/ParallelAlgorithms/Events/ObserveInterpolatedIntegralData.hpp @@ -0,0 +1,548 @@ +// 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/DataBox.hpp" +#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 "ParallelAlgorithms/Interpolation/Tags.hpp" +#include "PointwiseFunctions/AnalyticSolutions/Tags.hpp" +#include "PointwiseFunctions/Hydro/Tags.hpp" +#include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp" +#include "Utilities/Algorithm.hpp" +#include "Utilities/CallWithDynamicType.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/TaggedTuple.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>, + // MassAccretionRate Mdot + Parallel::ReductionDatum>, + // MagneticFlux Phi_BH + Parallel::ReductionDatum>, + // Total Energy Flux Edot + Parallel::ReductionDatum>, + // AngularMomentum Flux Ldot + 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) { + // this will be surface integral of a given tag for debugging purpose. + double surface_integral_value = 0.0; + double magnetic_flux = 0.0; + double mdot = 0.0; + double ldot = 0.0; + double edot = 0.0; + + const DataVector det_jacobian = + 1. / + get(get<::Events::Tags::ObserverDetInvJacobian>(box)); + const auto jac = + get<::Events::Tags::ObserverJacobian>(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); + + const double x0 = element_id.segment_id(0).endpoint(Side::Lower); + const double y0 = element_id.segment_id(1).endpoint(Side::Lower); + const double z0 = element_id.segment_id(2).endpoint(Side::Lower); + + // we only need to work on the elements where the intepolation surface + // lies in between + if ((lower_bound_bl < interp_val) and (upper_bound_bl >= interp_val)) { + // (D-1) dimensional mesh with original mesh sliced away in the + // interpolation dimension. + const auto new_mesh = mesh.slice_away(interp_dim); + + // Actually need the element logical value corresponding to + // interpolation target for next step + // This step of course assumes, we are already in the element in which + // the interpolation target lies. + // checked this formula on Oct 11 once again!!. + const double elm_interp_val = + (interp_val - ((upper_bound_bl + lower_bound_bl) / 2.)) / + ((upper_bound_bl - lower_bound_bl) / 2.); + + // when sliced away the target_points has values -1 (lowest val) for + // interp_dim. set it to appropriate element value as computed above. + 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); + + // we need to convert the above target points in element logical frame + // to inertial frame so that we can evaluate the metric quantities + auto ones = + make_with_value>( + target_points, 1.0); + auto target_points_inertial = + make_with_value>( + target_points, x0); + + get<1>(target_points_inertial) = y0; + get<2>(target_points_inertial) = z0; + for (size_t i = 0; i < VolumeDim; ++i) { + for (size_t j = 0; j < VolumeDim; ++j) { + target_points_inertial.get(i) += + jac.get(i, j)[0] * (target_points.get(j) + ones.get(j)); + } + } + using derived_classes = + tmpl::at; + auto* initial_data_ptr = + &Parallel::get(cache); + // all the metric related tags that we need for computing + // MassAccretionRate: shift, lapse, sqrt_det_spatial_metric + using metric_tags = + tmpl::list, + gr::Tags::SqrtDetSpatialMetric, + gr::Tags::Lapse>; + using metric_tuples = tuples::tagged_tuple_from_typelist; + auto metric_quantities = + call_with_dynamic_type( + initial_data_ptr, + [&target_points_inertial](const auto* const data_or_solution) { + return evolution::Initialization::initial_data( + *data_or_solution, target_points_inertial, 0.0, + metric_tags{}); + }); + + const auto rho = get(get>(box)); + const auto vr = get<0>( + get>( + box)); + const auto br = get<0>( + get>(box)); + const auto lf = get(get>(box)); + const auto record_tensor_component_impl = + [&interpolant, &surface_integral_value, &magnetic_flux, &mdot, &edot, + &ldot, &new_mesh, &det_jacobian, &rho, &vr, &br, &lf, + &metric_quantities](const auto& tensor) { + const auto new_tensor = interpolant.interpolate(tensor[0]); + const auto new_rho = interpolant.interpolate(rho); + const auto new_det_jacobian = interpolant.interpolate(det_jacobian); + const auto new_vr = interpolant.interpolate(vr); + const auto new_br = interpolant.interpolate(br); + const auto new_lf = interpolant.interpolate(lf); + + // this is where we perform surface intgral of the "new_tensor" + // basically whichever tensor integrated at the surface + const double surface_integral_value_contribution = + definite_integral(new_tensor * new_det_jacobian, new_mesh); + surface_integral_value += surface_integral_value_contribution; + + // this is where we compute Mdot + // here we directly compute all the metric quantities and + // interpolate the rest + // integral rho * u^r sqrt(-g) dphi dtheta + const auto shift = + get<0>(get>( + metric_quantities)); + const auto lapse = + get(get>(metric_quantities)); + const auto detsqrtsm = + get(get>( + metric_quantities)); + // we zero out potential outflow + DataVector integrand = + new_rho * new_lf * detsqrtsm * (lapse * new_vr - shift); + for (size_t i = 0; i < integrand.size(); i++) { + if (integrand[i] > 0) { + integrand[i] = 0; + } + } + const double mdot_contribution = + definite_integral(integrand * new_det_jacobian, new_mesh); + mdot += mdot_contribution; + + // magnetic flux + // 1/2 integral *F^(tr) sqrt(-g) dphi dtheta + // note B^r = alpha * (*F^(tr)) and sqrt(-g) = alpha * gamma + // so we only need gamma * B^r + DataVector magnetic_flux_integrand = detsqrtsm * abs(new_br); + const double magnetic_flux_contribution = + 0.5 * definite_integral( + magnetic_flux_integrand * new_det_jacobian, new_mesh); + magnetic_flux += magnetic_flux_contribution; + }; + 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)); + } + }; + 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, "Mdot", "MagneticFlux", + "Edot", "Ldot", "SurfaceIntegral"}, + ReductionData{observation_value.value, std::move(mdot), + std::move(magnetic_flux), std::move(edot), + std::move(ldot), 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