Skip to content

Commit b15abf3

Browse files
nilsvunilsvu
authored
Merge pull request #3096 from nilsleiffischer/multigrid_hierarchy
Add multigrid hierarchy functions
2 parents 03f16ef + c69b22a commit b15abf3

11 files changed

Lines changed: 432 additions & 2 deletions

File tree

‎src/Domain/Structure/SegmentId.cpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,3 +65,10 @@ size_t hash_value(const SegmentId& segment_id) noexcept {
6565
return hash;
6666
}
6767
// LCOV_EXCL_STOP
68+
69+
// NOLINTNEXTLINE(cert-dcl58-cpp)
70+
namespace std {
71+
size_t hash<SegmentId>::operator()(const SegmentId& segment_id) const noexcept {
72+
return hash_value(segment_id);
73+
}
74+
} // namespace std

‎src/Domain/Structure/SegmentId.hpp‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -229,6 +229,13 @@ inline bool SegmentId::overlaps(const SegmentId& other) const noexcept {
229229
// hash_value is called by boost::hash and related functions.
230230
size_t hash_value(const SegmentId& s) noexcept;
231231

232+
namespace std {
233+
template <>
234+
struct hash<SegmentId> {
235+
size_t operator()(const SegmentId& segment_id) const noexcept;
236+
};
237+
} // namespace std
238+
232239
inline bool operator==(const SegmentId& lhs, const SegmentId& rhs) noexcept {
233240
return (lhs.refinement_level() == rhs.refinement_level() and
234241
lhs.index() == rhs.index());

‎src/ParallelAlgorithms/LinearSolver/CMakeLists.txt‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,6 @@
11
# Distributed under the MIT License.
22
# See LICENSE.txt for details.
33

4-
add_subdirectory(Schwarz)
5-
64
set(LIBRARY ParallelLinearSolver)
75

86
add_spectre_library(${LIBRARY} INTERFACE)
@@ -33,4 +31,6 @@ add_subdirectory(Actions)
3331
add_subdirectory(AsynchronousSolvers)
3432
add_subdirectory(ConjugateGradient)
3533
add_subdirectory(Gmres)
34+
add_subdirectory(Multigrid)
3635
add_subdirectory(Richardson)
36+
add_subdirectory(Schwarz)

‎src/ParallelAlgorithms/LinearSolver/Multigrid/CMakeLists.txt‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
# Distributed under the MIT License.
2+
# See LICENSE.txt for details.
3+
4+
set(LIBRARY ParallelMultigrid)
5+
6+
add_spectre_library(${LIBRARY})
7+
8+
spectre_target_sources(
9+
${LIBRARY}
10+
PRIVATE
11+
Hierarchy.cpp
12+
)
13+
14+
spectre_target_headers(
15+
${LIBRARY}
16+
INCLUDE_DIRECTORY ${CMAKE_SOURCE_DIR}/src
17+
HEADERS
18+
Hierarchy.hpp
19+
Multigrid.hpp
20+
)
21+
22+
target_link_libraries(
23+
${LIBRARY}
24+
PUBLIC
25+
DomainStructure
26+
PRIVATE
27+
ErrorHandling
28+
Utilities
29+
)

‎src/ParallelAlgorithms/LinearSolver/Multigrid/Hierarchy.cpp‎

Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
1+
// Distributed under the MIT License.
2+
// See LICENSE.txt for details.
3+
4+
#include "ParallelAlgorithms/LinearSolver/Multigrid/Hierarchy.hpp"
5+
6+
#include <array>
7+
#include <cstddef>
8+
#include <unordered_set>
9+
#include <vector>
10+
11+
#include "Domain/Structure/ElementId.hpp"
12+
#include "Domain/Structure/SegmentId.hpp"
13+
#include "Domain/Structure/Side.hpp"
14+
#include "Utilities/ErrorHandling/Assert.hpp"
15+
#include "Utilities/GenerateInstantiations.hpp"
16+
17+
namespace LinearSolver::multigrid {
18+
19+
template <size_t Dim>
20+
std::vector<std::array<size_t, Dim>> coarsen(
21+
std::vector<std::array<size_t, Dim>> initial_refinement_levels) noexcept {
22+
for (auto& block_refinement : initial_refinement_levels) {
23+
for (size_t d = 0; d < Dim; ++d) {
24+
auto& refinement_level = gsl::at(block_refinement, d);
25+
if (refinement_level > 0) {
26+
--refinement_level;
27+
}
28+
}
29+
}
30+
return initial_refinement_levels;
31+
}
32+
33+
template <size_t Dim>
34+
ElementId<Dim> parent_id(const ElementId<Dim>& child_id) noexcept {
35+
std::array<SegmentId, Dim> parent_segment_ids = child_id.segment_ids();
36+
for (size_t d = 0; d < Dim; ++d) {
37+
auto& segment_id = gsl::at(parent_segment_ids, d);
38+
if (segment_id.refinement_level() > 0) {
39+
segment_id = segment_id.id_of_parent();
40+
}
41+
}
42+
return {child_id.block_id(), std::move(parent_segment_ids),
43+
child_id.grid_index() + 1};
44+
}
45+
46+
namespace {
47+
template <size_t Dim>
48+
void assert_finest_grid(
49+
const std::array<SegmentId, Dim>& parent_segment_ids,
50+
const std::array<size_t, Dim>& children_refinement_levels) noexcept {
51+
for (size_t d = 0; d < Dim; ++d) {
52+
ASSERT(gsl::at(children_refinement_levels, d) ==
53+
gsl::at(parent_segment_ids, d).refinement_level(),
54+
"On the finest grid, expected the children refinement levels "
55+
<< children_refinement_levels
56+
<< " to equal the refinement levels of the parent segment IDs "
57+
<< parent_segment_ids << " in dimension " << d << ".");
58+
}
59+
}
60+
61+
std::unordered_set<SegmentId> child_segment_ids_impl(
62+
const SegmentId& parent_segment_id,
63+
const size_t children_refinement_level) noexcept {
64+
ASSERT(parent_segment_id.refinement_level() == children_refinement_level or
65+
(children_refinement_level > 0 and
66+
parent_segment_id.refinement_level() ==
67+
children_refinement_level - 1),
68+
"The parent refinement level must be exactly 1 smaller than the "
69+
"children refinement level, or the same. Parent refinement level: "
70+
<< parent_segment_id.refinement_level()
71+
<< ", children refinement level: " << children_refinement_level);
72+
if (parent_segment_id.refinement_level() < children_refinement_level) {
73+
return {parent_segment_id.id_of_child(Side::Lower),
74+
parent_segment_id.id_of_child(Side::Upper)};
75+
} else {
76+
return {parent_segment_id};
77+
}
78+
}
79+
} // namespace
80+
81+
template <>
82+
std::unordered_set<ElementId<1>> child_ids<1>(
83+
const ElementId<1>& parent_id,
84+
const std::array<size_t, 1>& children_refinement_levels) noexcept {
85+
if (parent_id.grid_index() == 0) {
86+
#ifdef SPECTRE_DEBUG
87+
assert_finest_grid(parent_id.segment_ids(), children_refinement_levels);
88+
#endif // SPECTRE_DEBUG
89+
return {};
90+
}
91+
const std::unordered_set<SegmentId> child_segment_ids =
92+
child_segment_ids_impl(parent_id.segment_ids()[0],
93+
children_refinement_levels[0]);
94+
std::unordered_set<ElementId<1>> child_ids{};
95+
for (const auto& child_segment_id : child_segment_ids) {
96+
child_ids.emplace(parent_id.block_id(),
97+
std::array<SegmentId, 1>{child_segment_id},
98+
parent_id.grid_index() - 1);
99+
}
100+
return child_ids;
101+
}
102+
103+
template <>
104+
std::unordered_set<ElementId<2>> child_ids<2>(
105+
const ElementId<2>& parent_id,
106+
const std::array<size_t, 2>& children_refinement_levels) noexcept {
107+
if (parent_id.grid_index() == 0) {
108+
#ifdef SPECTRE_DEBUG
109+
assert_finest_grid(parent_id.segment_ids(), children_refinement_levels);
110+
#endif // SPECTRE_DEBUG
111+
return {};
112+
}
113+
std::array<std::unordered_set<SegmentId>, 2> child_segment_ids{};
114+
for (size_t d = 0; d < 2; ++d) {
115+
gsl::at(child_segment_ids, d) =
116+
child_segment_ids_impl(gsl::at(parent_id.segment_ids(), d),
117+
gsl::at(children_refinement_levels, d));
118+
}
119+
std::unordered_set<ElementId<2>> child_ids{};
120+
for (const auto& child_segment_id_x : child_segment_ids[0]) {
121+
for (const auto& child_segment_id_y : child_segment_ids[1]) {
122+
child_ids.emplace(
123+
parent_id.block_id(),
124+
std::array<SegmentId, 2>{{child_segment_id_x, child_segment_id_y}},
125+
parent_id.grid_index() - 1);
126+
}
127+
}
128+
return child_ids;
129+
}
130+
131+
template <>
132+
std::unordered_set<ElementId<3>> child_ids<3>(
133+
const ElementId<3>& parent_id,
134+
const std::array<size_t, 3>& children_refinement_levels) noexcept {
135+
if (parent_id.grid_index() == 0) {
136+
#ifdef SPECTRE_DEBUG
137+
assert_finest_grid(parent_id.segment_ids(), children_refinement_levels);
138+
#endif // SPECTRE_DEBUG
139+
return {};
140+
}
141+
std::array<std::unordered_set<SegmentId>, 3> child_segment_ids{};
142+
for (size_t d = 0; d < 3; ++d) {
143+
gsl::at(child_segment_ids, d) =
144+
child_segment_ids_impl(gsl::at(parent_id.segment_ids(), d),
145+
gsl::at(children_refinement_levels, d));
146+
}
147+
std::unordered_set<ElementId<3>> child_ids{};
148+
for (const auto& child_segment_id_x : child_segment_ids[0]) {
149+
for (const auto& child_segment_id_y : child_segment_ids[1]) {
150+
for (const auto& child_segment_id_z : child_segment_ids[2]) {
151+
child_ids.emplace(
152+
parent_id.block_id(),
153+
std::array<SegmentId, 3>{
154+
{child_segment_id_x, child_segment_id_y, child_segment_id_z}},
155+
parent_id.grid_index() - 1);
156+
}
157+
}
158+
}
159+
return child_ids;
160+
}
161+
162+
#define DIM(data) BOOST_PP_TUPLE_ELEM(0, data)
163+
#define INSTANTIATE(r, data) \
164+
template std::vector<std::array<size_t, DIM(data)>> coarsen( \
165+
std::vector<std::array<size_t, DIM(data)>> \
166+
initial_refinement_levels) noexcept; \
167+
template ElementId<DIM(data)> parent_id( \
168+
const ElementId<DIM(data)>& child_id) noexcept;
169+
170+
GENERATE_INSTANTIATIONS(INSTANTIATE, (1, 2, 3))
171+
172+
#undef DIM
173+
#undef INSTANTIATE
174+
175+
} // namespace LinearSolver::multigrid

‎src/ParallelAlgorithms/LinearSolver/Multigrid/Hierarchy.hpp‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
// Distributed under the MIT License.
2+
// See LICENSE.txt for details.
3+
4+
#pragma once
5+
6+
#include <array>
7+
#include <cstddef>
8+
#include <unordered_set>
9+
#include <vector>
10+
11+
/// \cond
12+
template <size_t Dim>
13+
struct ElementId;
14+
/// \endcond
15+
16+
namespace LinearSolver::multigrid {
17+
18+
/*!
19+
* \brief Coarsen the initial refinement levels of all blocks in the domain
20+
*
21+
* Simply decrement the refinement level uniformly over the entire domain.
22+
* Doesn't do anything for blocks that are already fully coarsened, so if the
23+
* return value equals the input argument the entire domain is fully coarsened.
24+
* Decrementing the refinement level means combining two elements into one,
25+
* thereby halving the number of elements per dimension.
26+
*
27+
* \tparam Dim The spatial dimension of the domain
28+
* \param initial_refinement_levels The refinement level in each block of the
29+
* domain and in every dimension.
30+
* \return std::vector<std::array<size_t, Dim>> The coarsened refinement levels
31+
* by decrementing every entry in `initial_refinement_levels` unless it is
32+
* already zero.
33+
*/
34+
template <size_t Dim>
35+
std::vector<std::array<size_t, Dim>> coarsen(
36+
std::vector<std::array<size_t, Dim>> initial_refinement_levels) noexcept;
37+
38+
/*!
39+
* \brief The element covering the `child_id` on the coarser grid.
40+
*
41+
* \tparam Dim The spatial dimension of the domain
42+
* \param child_id The ID of an element on the finer grid
43+
* \return ElementId<Dim> The ID of the element on the coarser grid that
44+
* covers the `child_id`. This parent element covers at most two child elements
45+
* per dimension.
46+
*/
47+
template <size_t Dim>
48+
ElementId<Dim> parent_id(const ElementId<Dim>& child_id) noexcept;
49+
50+
/*!
51+
* \brief The elements covering the `parent_id` on the finer grid.
52+
*
53+
* \tparam Dim The spatial dimension of the domain
54+
* \param parent_id The ID of an element on the coarser grid
55+
* \param children_refinement_levels The refinement level of the finer grid in
56+
* this block
57+
* \return std::unordered_set<ElementId<Dim>> The IDs of the elements on the
58+
* finer grid that cover the `parent_id`. Returns an empty set if the
59+
* `parent_id` is already on the finest grid. Else, returns at least one child
60+
* (if the grids have the same refinement levels) and at most
61+
* \f$2^\mathrm{Dim}\f$ children (if the grid is finer in every dimension).
62+
*/
63+
template <size_t Dim>
64+
std::unordered_set<ElementId<Dim>> child_ids(
65+
const ElementId<Dim>& parent_id,
66+
const std::array<size_t, Dim>& children_refinement_levels) noexcept;
67+
68+
} // namespace LinearSolver::multigrid

‎src/ParallelAlgorithms/LinearSolver/Multigrid/Multigrid.hpp‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
// Distributed under the MIT License.
2+
// See LICENSE.txt for details.
3+
4+
#pragma once
5+
6+
/// \ingroup LinearSolverGroup
7+
/// Items related to the multigrid linear solver
8+
namespace LinearSolver::multigrid {}

‎tests/Unit/Domain/Structure/Test_SegmentId.cpp‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#include <cstddef>
77
#include <limits>
88
#include <string>
9+
#include <unordered_set>
910

1011
#include "Domain/Structure/SegmentId.hpp"
1112
#include "Domain/Structure/Side.hpp"
@@ -72,6 +73,13 @@ SPECTRE_TEST_CASE("Unit.Domain.Structure.SegmentId", "[Domain][Unit]") {
7273
// Test output operator:
7374
SegmentId level_3_index_2(3, 2);
7475
CHECK(get_output(level_3_index_2) == "L3I2");
76+
77+
{
78+
INFO("Hash");
79+
CHECK(std::unordered_set<SegmentId>{SegmentId{2, 3}, SegmentId{2, 3},
80+
SegmentId{3, 2}} ==
81+
std::unordered_set<SegmentId>{SegmentId{3, 2}, SegmentId{2, 3}});
82+
}
7583
}
7684

7785
// [[OutputRegex, index = 8, refinement_level = 3]]

‎tests/Unit/ParallelAlgorithms/LinearSolver/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -79,5 +79,6 @@ add_subdirectory(Actions)
7979
add_subdirectory(AsynchronousSolvers)
8080
add_subdirectory(ConjugateGradient)
8181
add_subdirectory(Gmres)
82+
add_subdirectory(Multigrid)
8283
add_subdirectory(Richardson)
8384
add_subdirectory(Schwarz)

‎tests/Unit/ParallelAlgorithms/LinearSolver/Multigrid/CMakeLists.txt‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
# Distributed under the MIT License.
2+
# See LICENSE.txt for details.
3+
4+
set(LIBRARY "Test_ParallelMultigrid")
5+
6+
set(LIBRARY_SOURCES
7+
Test_Hierarchy.cpp
8+
)
9+
10+
add_test_library(
11+
${LIBRARY}
12+
"ParallelAlgorithms/LinearSolver/Multigrid"
13+
"${LIBRARY_SOURCES}"
14+
""
15+
)
16+
17+
target_link_libraries(
18+
${LIBRARY}
19+
PRIVATE
20+
DomainStructure
21+
ParallelMultigrid
22+
Utilities
23+
)

0 commit comments

Comments
 (0)