From 91c6744699d86ab2479e856718c5953d5a2bfa65 Mon Sep 17 00:00:00 2001 From: Nils Vu Date: Tue, 10 Dec 2024 15:12:01 -0800 Subject: [PATCH] Add volume data output to Schwarz solver --- .../LinearSolver/Schwarz/CMakeLists.txt | 1 + .../LinearSolver/Schwarz/ElementActions.hpp | 53 +++++- .../Schwarz/ObserveVolumeData.hpp | 179 ++++++++++++++++++ .../LinearSolver/Schwarz/Schwarz.hpp | 15 +- .../LinearSolver/Schwarz/Tags.hpp | 8 + support/Pipelines/Bbh/InitialData.yaml | 1 + tests/InputFiles/Elasticity/BentBeam.yaml | 1 + .../Elasticity/HalfSpaceMirror.yaml | 1 + .../Elasticity/SingleCoatingMirror.yaml | 1 + tests/InputFiles/Poisson/Lorentzian.yaml | 1 + .../Poisson/ProductOfSinusoids1D.yaml | 1 + .../Poisson/ProductOfSinusoids2D.yaml | 9 +- .../Poisson/ProductOfSinusoids3D.yaml | 1 + .../Punctures/MultiplePunctures.yaml | 1 + tests/InputFiles/Xcts/BinaryBlackHole.yaml | 1 + tests/InputFiles/Xcts/HeadOnBns.yaml | 1 + tests/InputFiles/Xcts/KerrSchild.yaml | 1 + tests/InputFiles/Xcts/TovStar.yaml | 1 + .../Schwarz/Test_SchwarzAlgorithm.yaml | 1 + 19 files changed, 260 insertions(+), 18 deletions(-) create mode 100644 src/ParallelAlgorithms/LinearSolver/Schwarz/ObserveVolumeData.hpp diff --git a/src/ParallelAlgorithms/LinearSolver/Schwarz/CMakeLists.txt b/src/ParallelAlgorithms/LinearSolver/Schwarz/CMakeLists.txt index d5c08087f3a5..b496c8fd2db4 100644 --- a/src/ParallelAlgorithms/LinearSolver/Schwarz/CMakeLists.txt +++ b/src/ParallelAlgorithms/LinearSolver/Schwarz/CMakeLists.txt @@ -19,6 +19,7 @@ spectre_target_headers( ComputeTags.hpp ElementActions.hpp ElementCenteredSubdomainData.hpp + ObserveVolumeData.hpp OverlapHelpers.hpp Schwarz.hpp SubdomainOperator.hpp diff --git a/src/ParallelAlgorithms/LinearSolver/Schwarz/ElementActions.hpp b/src/ParallelAlgorithms/LinearSolver/Schwarz/ElementActions.hpp index 8d13bd0c15a3..991cb6eef875 100644 --- a/src/ParallelAlgorithms/LinearSolver/Schwarz/ElementActions.hpp +++ b/src/ParallelAlgorithms/LinearSolver/Schwarz/ElementActions.hpp @@ -206,17 +206,21 @@ struct InitializeElement : tt::ConformsTo { subdomain_solver; using subdomain_solver_tag = Tags::SubdomainSolver, OptionsGroup>; + using volume_data_tag = + Tags::VolumeDataForOutput; public: // Iterable action using simple_tags_from_options = tmpl::list; - using const_global_cache_tags = tmpl::list>; + using const_global_cache_tags = + tmpl::list, + LinearSolver::Tags::OutputVolumeData>; - using simple_tags = - tmpl::list, - Tags::Weight, - domain::Tags::Faces>, - SubdomainDataBufferTag>; + using simple_tags = tmpl::list< + Tags::IntrudingExtents, Tags::Weight, + domain::Tags::Faces>, + SubdomainDataBufferTag, + LinearSolver::Tags::ObservationId, volume_data_tag>; using compute_tags = tmpl::list<>; template @@ -244,7 +248,9 @@ struct InitializeElement : tt::ConformsTo { const gsl::not_null>*> intruding_overlap_weights, const gsl::not_null subdomain_data, - [[maybe_unused]] const gsl::not_null*> + const gsl::not_null observation_id, + const gsl::not_null /*volume_data_for_output*/, + const gsl::not_null*> subdomain_solver, const Element& element, const Mesh& mesh, const tnsr::I& logical_coords, @@ -299,14 +305,23 @@ struct InitializeElement : tt::ConformsTo { // Subdomain solver // The subdomain solver initially gets created from options on each element. // Then we have to copy it around during AMR. - if constexpr (sizeof...(AmrData) == 1) { + if constexpr (sizeof...(AmrData) == 0) { + *observation_id = 0; + (void)subdomain_solver; + } else { if constexpr (tt::is_a_v) { // h-refinement: copy from the parent + *observation_id = + get>(amr_data...); *subdomain_solver = get(amr_data...)->get_clone(); } else if constexpr (tt::is_a_v) { // h-coarsening, copy from one of the children (doesn't matter which) + *observation_id = get>( + amr_data.begin()->second...); *subdomain_solver = get(amr_data.begin()->second...)->get_clone(); + } else { + (void)observation_id; } (*subdomain_solver)->reset(); } @@ -453,6 +468,15 @@ struct SolveSubdomain { db::get>(box)); } + // Record subdomain solution for volume data output + if (db::get>(box)) { + db::mutate>( + [&subdomain_solution](const auto volume_data) { + volume_data->element_data = subdomain_solution.element_data; + }, + make_not_null(&box)); + } + // Apply weighting if (LIKELY(max_overlap > 0)) { subdomain_solution.element_data *= @@ -537,11 +561,22 @@ struct ReceiveOverlapSolution { pretty_type::name(), iteration_id); } - // Add solutions on overlaps to this element's solution in a weighted sum + // Extract overlap solutions from inbox const auto received_overlap_solutions = std::move(tuples::get(inboxes) .extract(iteration_id) .mapped()); + + // Record overlap solutions for volume data output + if (db::get>(box)) { + db::mutate>( + [&received_overlap_solutions](const auto volume_data) { + volume_data->overlap_data = received_overlap_solutions; + }, + make_not_null(&box)); + } + + // Add overlap solutions to this element's solution in a weighted sum db::mutate( [&received_overlap_solutions]( const auto fields, const Index& full_extents, diff --git a/src/ParallelAlgorithms/LinearSolver/Schwarz/ObserveVolumeData.hpp b/src/ParallelAlgorithms/LinearSolver/Schwarz/ObserveVolumeData.hpp new file mode 100644 index 000000000000..f7445f14dd8d --- /dev/null +++ b/src/ParallelAlgorithms/LinearSolver/Schwarz/ObserveVolumeData.hpp @@ -0,0 +1,179 @@ +// Distributed under the MIT License. +// See LICENSE.txt for details. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "DataStructures/DataBox/DataBox.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/TypeOfObservation.hpp" +#include "IO/Observer/VolumeActions.hpp" +#include "NumericalAlgorithms/Spectral/Mesh.hpp" +#include "Parallel/AlgorithmExecution.hpp" +#include "Parallel/ArrayComponentId.hpp" +#include "Parallel/GlobalCache.hpp" +#include "Parallel/Invoke.hpp" +#include "Parallel/Local.hpp" +#include "ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp" +#include "Utilities/Algorithm.hpp" +#include "Utilities/Gsl.hpp" +#include "Utilities/MakeString.hpp" +#include "Utilities/PrettyType.hpp" +#include "Utilities/TMPL.hpp" +#include "Utilities/TaggedTuple.hpp" + +namespace LinearSolver::Schwarz::detail { + +template +struct RegisterWithVolumeObserver { + template + static std::pair + register_info(const db::DataBox& box, + const ArrayIndex& /*array_index*/) { + const std::optional section_observation_key = + observers::get_section_observation_key(box); + const std::string subfile_path = "/" + pretty_type::name() + + section_observation_key.value_or(""); + return {observers::TypeOfObservation::Volume, + observers::ObservationKey(subfile_path)}; + } +}; + +// Contribute the volume data recorded in the other actions to the observer at +// the end of a step. +template +struct ObserveVolumeData { + private: + using fields_tag = FieldsTag; + using residual_tag = + db::add_tag_prefix; + static constexpr size_t Dim = SubdomainOperator::volume_dim; + using SubdomainData = + ElementCenteredSubdomainData; + using volume_data_tag = + Tags::VolumeDataForOutput; + using VolumeDataVars = typename volume_data_tag::type::ElementData; + + public: + template + static Parallel::iterable_action_return_t apply( + db::DataBox& box, + const tuples::TaggedTuple& /*inboxes*/, + Parallel::GlobalCache& cache, + const ElementId& element_id, const ActionList /*meta*/, + const ParallelComponent* const /*meta*/) { + if (not db::get>(box)) { + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } + const auto& volume_data = db::get(box); + const auto& observation_id = + db::get>(box); + const auto& mesh = db::get>(box); + const auto& inertial_coords = + db::get>(box); + // Collect tensor components to observe + std::vector components{}; + const auto record_tensor_components = [&components]( + const auto tensor_tag_v, + const auto& tensor, + const std::string& suffix = "") { + using tensor_tag = std::decay_t; + using TensorType = std::decay_t; + using VectorType = typename TensorType::type; + using ValueType = typename VectorType::value_type; + for (size_t i = 0; i < tensor.size(); ++i) { + const std::string component_name = + db::tag_name() + suffix + tensor.component_suffix(i); + if constexpr (std::is_same_v>) { + components.emplace_back("Re(" + component_name + ")", + real(tensor[i])); + components.emplace_back("Im(" + component_name + ")", + imag(tensor[i])); + } else { + components.emplace_back(component_name, tensor[i]); + } + } + }; + record_tensor_components(domain::Tags::Coordinates{}, + inertial_coords); + const auto& all_intruding_extents = + db::get>(box); + const VolumeDataVars zero_vars{mesh.number_of_grid_points(), 0.}; + tmpl::for_each( + [&volume_data, &record_tensor_components, &mesh, &all_intruding_extents, + &zero_vars, &element_id](auto tag_v) { + using tag = tmpl::type_from; + record_tensor_components(tag{}, get(volume_data.element_data), + "_Center"); + for (const auto direction : Direction::all_directions()) { + const auto direction_predicate = + [&direction](const auto& overlap_data) { + return overlap_data.first.direction() == direction; + }; + const auto num_overlaps = + alg::count_if(volume_data.overlap_data, direction_predicate); + if (num_overlaps == 0) { + // No overlap data for this direction (e.g. external boundary), + // record zero + record_tensor_components(tag{}, get(zero_vars), + "_Overlap" + get_output(direction)); + } else if (num_overlaps == 1) { + // Overlap data from a single neighbor, record it + const auto& overlap_data = + alg::find_if(volume_data.overlap_data, direction_predicate) + ->second; + const auto& intruding_extents = + gsl::at(all_intruding_extents, direction.dimension()); + record_tensor_components( + tag{}, + get(extended_overlap_data(overlap_data, mesh.extents(), + intruding_extents, direction)), + "_Overlap" + get_output(direction)); + } else { + ERROR("Multiple neighbors (" + << num_overlaps << ") overlap with element " << element_id + << " in direction " << direction + << ", but we can record only one for volume data output."); + } + } + }); + + // Contribute tensor components to observer + auto& local_observer = *Parallel::local_branch( + Parallel::get_parallel_component>( + cache)); + const std::optional section_observation_key = + observers::get_section_observation_key(box); + const std::string subfile_path = "/" + pretty_type::name() + + section_observation_key.value_or(""); + Parallel::simple_action( + local_observer, observers::ObservationId(observation_id, subfile_path), + subfile_path, + Parallel::make_array_component_id(element_id), + ElementVolumeData{element_id, std::move(components), mesh}); + + // Increment observation ID + db::mutate>( + [](const auto local_observation_id) { ++(*local_observation_id); }, + make_not_null(&box)); + return {Parallel::AlgorithmExecution::Continue, std::nullopt}; + } +}; + +} // namespace LinearSolver::Schwarz::detail diff --git a/src/ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp b/src/ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp index 5c0dedd3708d..c7691f6a6aa0 100644 --- a/src/ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp +++ b/src/ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp @@ -6,6 +6,7 @@ #include "IO/Observer/Helpers.hpp" #include "ParallelAlgorithms/LinearSolver/AsynchronousSolvers/ElementActions.hpp" #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementActions.hpp" +#include "ParallelAlgorithms/LinearSolver/Schwarz/ObserveVolumeData.hpp" #include "Utilities/ProtocolHelpers.hpp" #include "Utilities/TMPL.hpp" @@ -174,11 +175,13 @@ struct Schwarz { using amr_projectors = initialize_element; - using register_element = - tmpl::list, - detail::RegisterElement>; + using register_element = tmpl::list< + async_solvers::RegisterElement, + detail::RegisterElement, + observers::Actions::RegisterWithObservers< + detail::RegisterWithVolumeObserver>>; template using solve = tmpl::list< @@ -189,6 +192,8 @@ struct Schwarz { ArraySectionIdTag>, detail::ReceiveOverlapSolution, + detail::ObserveVolumeData, ApplyOperatorActions, async_solvers::CompleteStep>; diff --git a/src/ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp b/src/ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp index c9382702f180..c02b1d52b972 100644 --- a/src/ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp +++ b/src/ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp @@ -185,6 +185,14 @@ struct SummedIntrudingOverlapWeights : db::SimpleTag { using type = Scalar; }; +/// Buffer for recording volume data +template +struct VolumeDataForOutput : db::SimpleTag { + using type = SubdomainDataType; + static std::string name() { + return "VolumeDataForOutput(" + pretty_type::name() + ")"; + } +}; } // namespace Tags } // namespace LinearSolver::Schwarz diff --git a/support/Pipelines/Bbh/InitialData.yaml b/support/Pipelines/Bbh/InitialData.yaml index 2adfdf724d27..aa9ad3a282e6 100644 --- a/support/Pipelines/Bbh/InitialData.yaml +++ b/support/Pipelines/Bbh/InitialData.yaml @@ -225,6 +225,7 @@ LinearSolver: BoundaryConditions: Auto SkipResets: True ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: InnerRadius: *outer_shell_inner_radius diff --git a/tests/InputFiles/Elasticity/BentBeam.yaml b/tests/InputFiles/Elasticity/BentBeam.yaml index f3739997c2b3..1e7065b002b9 100644 --- a/tests/InputFiles/Elasticity/BentBeam.yaml +++ b/tests/InputFiles/Elasticity/BentBeam.yaml @@ -90,6 +90,7 @@ LinearSolver: ExplicitInverse: WriteMatrixToFile: None ObservePerCoreReductions: False + OutputVolumeData: False EventsAndTriggersAtIterations: - Trigger: HasConverged diff --git a/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml b/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml index 1642d6e9e876..f7d3d7739bd6 100644 --- a/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml +++ b/tests/InputFiles/Elasticity/HalfSpaceMirror.yaml @@ -116,6 +116,7 @@ LinearSolver: WriteMatrixToFile: None BoundaryConditions: Auto ObservePerCoreReductions: False + OutputVolumeData: False EventsAndTriggersAtIterations: - Trigger: HasConverged diff --git a/tests/InputFiles/Elasticity/SingleCoatingMirror.yaml b/tests/InputFiles/Elasticity/SingleCoatingMirror.yaml index cdaa115754bc..55e8fc960109 100644 --- a/tests/InputFiles/Elasticity/SingleCoatingMirror.yaml +++ b/tests/InputFiles/Elasticity/SingleCoatingMirror.yaml @@ -131,6 +131,7 @@ LinearSolver: WriteMatrixToFile: None BoundaryConditions: Auto ObservePerCoreReductions: False + OutputVolumeData: False EventsAndTriggersAtIterations: - Trigger: HasConverged diff --git a/tests/InputFiles/Poisson/Lorentzian.yaml b/tests/InputFiles/Poisson/Lorentzian.yaml index a001d829376c..eb538cd47fce 100644 --- a/tests/InputFiles/Poisson/Lorentzian.yaml +++ b/tests/InputFiles/Poisson/Lorentzian.yaml @@ -108,6 +108,7 @@ LinearSolver: ExplicitInverse: WriteMatrixToFile: None ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: InnerRadius: *outer_shell_inner_radius diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml index 35ecae72fabf..b921e2c34d0a 100644 --- a/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml +++ b/tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml @@ -152,6 +152,7 @@ LinearSolver: ExplicitInverse: WriteMatrixToFile: "SubdomainMatrix" ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: None diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml index ed297d7661d5..8432633cb277 100644 --- a/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml +++ b/tests/InputFiles/Poisson/ProductOfSinusoids2D.yaml @@ -44,11 +44,11 @@ RandomizeInitialGuess: None DomainCreator: Rectangle: - LowerBound: [-1.570796326794896, 0.] - UpperBound: [3.141592653589793, 1.570796326794896] + LowerBound: [0., 0.] + UpperBound: [3.141592653589793, 3.141592653589793] Distribution: [Linear, Linear] - InitialRefinement: [1, 0] - InitialGridPoints: [4, 3] + InitialRefinement: [4, 4] + InitialGridPoints: [8, 8] TimeDependence: None BoundaryConditions: - AnalyticSolution: @@ -93,6 +93,7 @@ LinearSolver: ExplicitInverse: WriteMatrixToFile: None ObservePerCoreReductions: False + OutputVolumeData: True RadiallyCompressedCoordinates: None diff --git a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml b/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml index 29c94979d1d5..77ac0536097d 100644 --- a/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml +++ b/tests/InputFiles/Poisson/ProductOfSinusoids3D.yaml @@ -97,6 +97,7 @@ LinearSolver: ExplicitInverse: WriteMatrixToFile: None ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: None diff --git a/tests/InputFiles/Punctures/MultiplePunctures.yaml b/tests/InputFiles/Punctures/MultiplePunctures.yaml index 760bcd6b43eb..6de03b8eabe7 100644 --- a/tests/InputFiles/Punctures/MultiplePunctures.yaml +++ b/tests/InputFiles/Punctures/MultiplePunctures.yaml @@ -155,6 +155,7 @@ LinearSolver: BoundaryConditions: Auto SkipResets: True ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: None diff --git a/tests/InputFiles/Xcts/BinaryBlackHole.yaml b/tests/InputFiles/Xcts/BinaryBlackHole.yaml index 9fc167c771ae..b3aeeaffc303 100644 --- a/tests/InputFiles/Xcts/BinaryBlackHole.yaml +++ b/tests/InputFiles/Xcts/BinaryBlackHole.yaml @@ -228,6 +228,7 @@ LinearSolver: BoundaryConditions: Auto SkipResets: True ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: InnerRadius: *outer_shell_inner_radius diff --git a/tests/InputFiles/Xcts/HeadOnBns.yaml b/tests/InputFiles/Xcts/HeadOnBns.yaml index 7c8737888917..f84fc798d198 100644 --- a/tests/InputFiles/Xcts/HeadOnBns.yaml +++ b/tests/InputFiles/Xcts/HeadOnBns.yaml @@ -167,6 +167,7 @@ LinearSolver: BoundaryConditions: Auto SkipResets: True ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: InnerRadius: *outer_shell_inner_radius diff --git a/tests/InputFiles/Xcts/KerrSchild.yaml b/tests/InputFiles/Xcts/KerrSchild.yaml index d17b26cac044..3421b4d8aaa6 100644 --- a/tests/InputFiles/Xcts/KerrSchild.yaml +++ b/tests/InputFiles/Xcts/KerrSchild.yaml @@ -198,6 +198,7 @@ LinearSolver: BoundaryConditions: Auto SkipResets: True ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: None diff --git a/tests/InputFiles/Xcts/TovStar.yaml b/tests/InputFiles/Xcts/TovStar.yaml index 32316f79f4e6..bebb66b2a620 100644 --- a/tests/InputFiles/Xcts/TovStar.yaml +++ b/tests/InputFiles/Xcts/TovStar.yaml @@ -117,6 +117,7 @@ LinearSolver: BoundaryConditions: Auto SkipResets: True ObservePerCoreReductions: False + OutputVolumeData: False RadiallyCompressedCoordinates: None diff --git a/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Test_SchwarzAlgorithm.yaml b/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Test_SchwarzAlgorithm.yaml index 97df12da7bb7..8f7f3d49dcee 100644 --- a/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Test_SchwarzAlgorithm.yaml +++ b/tests/Unit/ParallelAlgorithms/LinearSolver/Schwarz/Test_SchwarzAlgorithm.yaml @@ -73,6 +73,7 @@ SchwarzSmoother: ExplicitInverse: WriteMatrixToFile: None ObservePerCoreReductions: False + OutputVolumeData: True ConvergenceReason: NumIterations