Skip to content

Commit

Permalink
Add volume data output to Schwarz solver
Browse files Browse the repository at this point in the history
  • Loading branch information
nilsvu committed Dec 10, 2024
1 parent 8cb1d9e commit def8f84
Show file tree
Hide file tree
Showing 18 changed files with 259 additions and 18 deletions.
53 changes: 44 additions & 9 deletions src/ParallelAlgorithms/LinearSolver/Schwarz/ElementActions.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -206,17 +206,21 @@ struct InitializeElement : tt::ConformsTo<amr::protocols::Projector> {
subdomain_solver<FieldsTag, SubdomainOperator, SubdomainPreconditioners>;
using subdomain_solver_tag =
Tags::SubdomainSolver<std::unique_ptr<SubdomainSolver>, OptionsGroup>;
using volume_data_tag =
Tags::VolumeDataForOutput<SubdomainData, OptionsGroup>;

public: // Iterable action
using simple_tags_from_options = tmpl::list<subdomain_solver_tag>;

using const_global_cache_tags = tmpl::list<Tags::MaxOverlap<OptionsGroup>>;
using const_global_cache_tags =
tmpl::list<Tags::MaxOverlap<OptionsGroup>,
LinearSolver::Tags::OutputVolumeData<OptionsGroup>>;

using simple_tags =
tmpl::list<Tags::IntrudingExtents<Dim, OptionsGroup>,
Tags::Weight<OptionsGroup>,
domain::Tags::Faces<Dim, Tags::Weight<OptionsGroup>>,
SubdomainDataBufferTag<SubdomainData, OptionsGroup>>;
using simple_tags = tmpl::list<
Tags::IntrudingExtents<Dim, OptionsGroup>, Tags::Weight<OptionsGroup>,
domain::Tags::Faces<Dim, Tags::Weight<OptionsGroup>>,
SubdomainDataBufferTag<SubdomainData, OptionsGroup>,
LinearSolver::Tags::ObservationId<OptionsGroup>, volume_data_tag>;
using compute_tags = tmpl::list<>;
template <typename DbTagsList, typename... InboxTags, typename Metavariables,
typename ActionList, typename ParallelComponent>
Expand Down Expand Up @@ -244,7 +248,9 @@ struct InitializeElement : tt::ConformsTo<amr::protocols::Projector> {
const gsl::not_null<DirectionMap<Dim, Scalar<DataVector>>*>
intruding_overlap_weights,
const gsl::not_null<SubdomainData*> subdomain_data,
[[maybe_unused]] const gsl::not_null<std::unique_ptr<SubdomainSolver>*>
const gsl::not_null<size_t*> observation_id,
const gsl::not_null<SubdomainData*> /*volume_data_for_output*/,
const gsl::not_null<std::unique_ptr<SubdomainSolver>*>
subdomain_solver,
const Element<Dim>& element, const Mesh<Dim>& mesh,
const tnsr::I<DataVector, Dim, Frame::ElementLogical>& logical_coords,
Expand Down Expand Up @@ -299,14 +305,23 @@ struct InitializeElement : tt::ConformsTo<amr::protocols::Projector> {
// 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<tuples::TaggedTuple, AmrData...>) {
// h-refinement: copy from the parent
*observation_id =
get<LinearSolver::Tags::ObservationId<OptionsGroup>>(amr_data...);
*subdomain_solver = get<subdomain_solver_tag>(amr_data...)->get_clone();
} else if constexpr (tt::is_a_v<std::unordered_map, AmrData...>) {
// h-coarsening, copy from one of the children (doesn't matter which)
*observation_id = get<LinearSolver::Tags::ObservationId<OptionsGroup>>(
amr_data.begin()->second...);
*subdomain_solver =
get<subdomain_solver_tag>(amr_data.begin()->second...)->get_clone();
} else {
(void)observation_id;
}
(*subdomain_solver)->reset();
}
Expand Down Expand Up @@ -453,6 +468,15 @@ struct SolveSubdomain {
db::get<Tags::ObservePerCoreReductions<OptionsGroup>>(box));
}

// Record subdomain solution for volume data output
if (db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
db::mutate<Tags::VolumeDataForOutput<SubdomainData, OptionsGroup>>(
[&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 *=
Expand Down Expand Up @@ -537,11 +561,22 @@ struct ReceiveOverlapSolution {
pretty_type::name<OptionsGroup>(), 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<overlap_solution_inbox_tag>(inboxes)
.extract(iteration_id)
.mapped());

// Record overlap solutions for volume data output
if (db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
db::mutate<Tags::VolumeDataForOutput<SubdomainData, OptionsGroup>>(
[&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<fields_tag>(
[&received_overlap_solutions](
const auto fields, const Index<Dim>& full_extents,
Expand Down
179 changes: 179 additions & 0 deletions src/ParallelAlgorithms/LinearSolver/Schwarz/ObserveVolumeData.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,179 @@
// Distributed under the MIT License.
// See LICENSE.txt for details.

#pragma once

#include <cstddef>
#include <optional>
#include <string>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>

#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 <typename OptionsGroup, typename ArraySectionIdTag>
struct RegisterWithVolumeObserver {
template <typename ParallelComponent, typename DbTagsList,
typename ArrayIndex>
static std::pair<observers::TypeOfObservation, observers::ObservationKey>
register_info(const db::DataBox<DbTagsList>& box,
const ArrayIndex& /*array_index*/) {
const std::optional<std::string> section_observation_key =
observers::get_section_observation_key<ArraySectionIdTag>(box);
const std::string subfile_path = "/" + pretty_type::name<OptionsGroup>() +
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 <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator,
typename ArraySectionIdTag>
struct ObserveVolumeData {
private:
using fields_tag = FieldsTag;
using residual_tag =
db::add_tag_prefix<LinearSolver::Tags::Residual, fields_tag>;
static constexpr size_t Dim = SubdomainOperator::volume_dim;
using SubdomainData =
ElementCenteredSubdomainData<Dim, typename residual_tag::tags_list>;
using volume_data_tag =
Tags::VolumeDataForOutput<SubdomainData, OptionsGroup>;
using VolumeDataVars = typename volume_data_tag::type::ElementData;

public:
template <typename DbTagsList, typename... InboxTags, typename Metavariables,
size_t Dim, typename ActionList, typename ParallelComponent>
static Parallel::iterable_action_return_t apply(
db::DataBox<DbTagsList>& box,
const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
Parallel::GlobalCache<Metavariables>& cache,
const ElementId<Dim>& element_id, const ActionList /*meta*/,
const ParallelComponent* const /*meta*/) {
if (not db::get<LinearSolver::Tags::OutputVolumeData<OptionsGroup>>(box)) {
return {Parallel::AlgorithmExecution::Continue, std::nullopt};
}
const auto& volume_data = db::get<volume_data_tag>(box);
const auto& observation_id =
db::get<LinearSolver::Tags::ObservationId<OptionsGroup>>(box);
const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
const auto& inertial_coords =
db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box);
// Collect tensor components to observe
std::vector<TensorComponent> 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<decltype(tensor_tag_v)>;
using TensorType = std::decay_t<decltype(tensor)>;
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<tensor_tag>() + suffix + tensor.component_suffix(i);
if constexpr (std::is_same_v<ValueType, std::complex<double>>) {
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<Dim, Frame::Inertial>{},
inertial_coords);
const auto& all_intruding_extents =
db::get<Tags::IntrudingExtents<Dim, OptionsGroup>>(box);
const VolumeDataVars zero_vars{mesh.number_of_grid_points(), 0.};
tmpl::for_each<typename VolumeDataVars::tags_list>(
[&volume_data, &record_tensor_components, &mesh, &all_intruding_extents,
&zero_vars, &element_id](auto tag_v) {
using tag = tmpl::type_from<decltype(tag_v)>;
record_tensor_components(tag{}, get<tag>(volume_data.element_data),
"_Center");
for (const auto direction : Direction<Dim>::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<tag>(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<tag>(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<observers::Observer<Metavariables>>(
cache));
const std::optional<std::string> section_observation_key =
observers::get_section_observation_key<ArraySectionIdTag>(box);
const std::string subfile_path = "/" + pretty_type::name<OptionsGroup>() +
section_observation_key.value_or("");
Parallel::simple_action<observers::Actions::ContributeVolumeData>(
local_observer, observers::ObservationId(observation_id, subfile_path),
subfile_path,
Parallel::make_array_component_id<ParallelComponent>(element_id),
ElementVolumeData{element_id, std::move(components), mesh});

// Increment observation ID
db::mutate<LinearSolver::Tags::ObservationId<OptionsGroup>>(
[](const auto local_observation_id) { ++(*local_observation_id); },
make_not_null(&box));
return {Parallel::AlgorithmExecution::Continue, std::nullopt};
}
};

} // namespace LinearSolver::Schwarz::detail
15 changes: 10 additions & 5 deletions src/ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -174,11 +175,13 @@ struct Schwarz {

using amr_projectors = initialize_element;

using register_element =
tmpl::list<async_solvers::RegisterElement<FieldsTag, OptionsGroup,
SourceTag, ArraySectionIdTag>,
detail::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
ArraySectionIdTag>>;
using register_element = tmpl::list<
async_solvers::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
ArraySectionIdTag>,
detail::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
ArraySectionIdTag>,
observers::Actions::RegisterWithObservers<
detail::RegisterWithVolumeObserver<OptionsGroup, ArraySectionIdTag>>>;

template <typename ApplyOperatorActions, typename Label = OptionsGroup>
using solve = tmpl::list<
Expand All @@ -189,6 +192,8 @@ struct Schwarz {
ArraySectionIdTag>,
detail::ReceiveOverlapSolution<FieldsTag, OptionsGroup,
SubdomainOperator>,
detail::ObserveVolumeData<FieldsTag, OptionsGroup, SubdomainOperator,
ArraySectionIdTag>,
ApplyOperatorActions,
async_solvers::CompleteStep<FieldsTag, OptionsGroup, SourceTag, Label,
ArraySectionIdTag>>;
Expand Down
8 changes: 8 additions & 0 deletions src/ParallelAlgorithms/LinearSolver/Schwarz/Tags.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -185,6 +185,14 @@ struct SummedIntrudingOverlapWeights : db::SimpleTag {
using type = Scalar<DataVector>;
};

/// Buffer for recording volume data
template <typename SubdomainDataType, typename OptionsGroup>
struct VolumeDataForOutput : db::SimpleTag {
using type = SubdomainDataType;
static std::string name() {
return "VolumeDataForOutput(" + pretty_type::name<OptionsGroup>() + ")";
}
};
} // namespace Tags
} // namespace LinearSolver::Schwarz

Expand Down
1 change: 1 addition & 0 deletions support/Pipelines/Bbh/InitialData.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -225,6 +225,7 @@ LinearSolver:
BoundaryConditions: Auto
SkipResets: True
ObservePerCoreReductions: False
OutputVolumeData: False

RadiallyCompressedCoordinates:
InnerRadius: *outer_shell_inner_radius
Expand Down
1 change: 1 addition & 0 deletions tests/InputFiles/Elasticity/BentBeam.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@ LinearSolver:
ExplicitInverse:
WriteMatrixToFile: None
ObservePerCoreReductions: False
OutputVolumeData: False

EventsAndTriggersAtIterations:
- Trigger: HasConverged
Expand Down
1 change: 1 addition & 0 deletions tests/InputFiles/Elasticity/HalfSpaceMirror.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,7 @@ LinearSolver:
WriteMatrixToFile: None
BoundaryConditions: Auto
ObservePerCoreReductions: False
OutputVolumeData: False

EventsAndTriggersAtIterations:
- Trigger: HasConverged
Expand Down
1 change: 1 addition & 0 deletions tests/InputFiles/Elasticity/SingleCoatingMirror.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,7 @@ LinearSolver:
WriteMatrixToFile: None
BoundaryConditions: Auto
ObservePerCoreReductions: False
OutputVolumeData: False

EventsAndTriggersAtIterations:
- Trigger: HasConverged
Expand Down
1 change: 1 addition & 0 deletions tests/InputFiles/Poisson/Lorentzian.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,7 @@ LinearSolver:
ExplicitInverse:
WriteMatrixToFile: None
ObservePerCoreReductions: False
OutputVolumeData: False

RadiallyCompressedCoordinates:
InnerRadius: *outer_shell_inner_radius
Expand Down
1 change: 1 addition & 0 deletions tests/InputFiles/Poisson/ProductOfSinusoids1D.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -152,6 +152,7 @@ LinearSolver:
ExplicitInverse:
WriteMatrixToFile: "SubdomainMatrix"
ObservePerCoreReductions: False
OutputVolumeData: False

RadiallyCompressedCoordinates: None

Expand Down
Loading

0 comments on commit def8f84

Please sign in to comment.