deal.II version GIT relicensing-1982-gbb401394a0 2024-10-13 08:40:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
mg_level_global_transfer.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
18
20
21#include <deal.II/fe/fe.h>
22
23#include <deal.II/grid/tria.h>
25
32#include <deal.II/lac/vector.h>
33
36#include <deal.II/multigrid/mg_transfer.templates.h>
38
39#include <algorithm>
40
42
43
44/* ------------------ MGLevelGlobalTransfer<VectorType> ----------------- */
45
46
47template <typename VectorType>
48template <int dim, int spacedim>
49void
51 const DoFHandler<dim, spacedim> &mg_dof)
52{
54 mg_constrained_dofs,
55 copy_indices,
56 copy_indices_global_mine,
57 copy_indices_level_mine);
58
59 // check if we can run a plain copy operation between the global DoFs and
60 // the finest level.
61 bool my_perform_plain_copy =
62 (copy_indices.back().size() == mg_dof.locally_owned_dofs().n_elements()) &&
63 (mg_dof.locally_owned_dofs().n_elements() ==
64 mg_dof
66 .n_elements());
67 if (my_perform_plain_copy)
68 {
69 AssertDimension(copy_indices_global_mine.back().size(), 0);
70 AssertDimension(copy_indices_level_mine.back().size(), 0);
71
72 // check whether there is a renumbering of degrees of freedom on
73 // either the finest level or the global dofs, which means that we
74 // cannot apply a plain copy
75 for (const auto &copy_index : copy_indices.back())
76 if (copy_index.first != copy_index.second)
77 {
78 my_perform_plain_copy = false;
79 break;
80 }
81 }
82
83 // now do a global reduction over all processors to see what operation
84 // they can agree upon
87 &mg_dof.get_triangulation()))
88 perform_plain_copy =
89 (Utilities::MPI::min(my_perform_plain_copy ? 1 : 0,
90 ptria->get_mpi_communicator()) == 1);
91 else
92 perform_plain_copy = my_perform_plain_copy;
93}
94
95
96
97template <typename VectorType>
98void
100{
101 sizes.resize(0);
102 copy_indices.clear();
103 copy_indices_global_mine.clear();
104 copy_indices_level_mine.clear();
105 component_to_block_map.resize(0);
106 mg_constrained_dofs = nullptr;
107 perform_plain_copy = false;
108}
109
110
111
112template <typename VectorType>
113void
115{
116 for (unsigned int level = 0; level < copy_indices.size(); ++level)
117 {
118 for (unsigned int i = 0; i < copy_indices[level].size(); ++i)
119 os << "copy_indices[" << level << "]\t" << copy_indices[level][i].first
120 << '\t' << copy_indices[level][i].second << std::endl;
121 }
122
123 for (unsigned int level = 0; level < copy_indices_level_mine.size(); ++level)
124 {
125 for (unsigned int i = 0; i < copy_indices_level_mine[level].size(); ++i)
126 os << "copy_ifrom [" << level << "]\t"
127 << copy_indices_level_mine[level][i].first << '\t'
128 << copy_indices_level_mine[level][i].second << std::endl;
129 }
130 for (unsigned int level = 0; level < copy_indices_global_mine.size(); ++level)
131 {
132 for (unsigned int i = 0; i < copy_indices_global_mine[level].size(); ++i)
133 os << "copy_ito [" << level << "]\t"
134 << copy_indices_global_mine[level][i].first << '\t'
135 << copy_indices_global_mine[level][i].second << std::endl;
136 }
137}
138
139
140
141template <typename VectorType>
142std::size_t
144{
145 std::size_t result = sizeof(*this);
147 result += MemoryConsumption::memory_consumption(copy_indices);
148 result += MemoryConsumption::memory_consumption(copy_indices_global_mine);
149 result += MemoryConsumption::memory_consumption(copy_indices_level_mine);
150
151 return result;
152}
153
154
155
156/* ------------------ MGLevelGlobalTransfer<VectorType> ----------------- */
157
158namespace
159{
160 template <int dim, int spacedim, typename Number>
161 void
162 fill_internal(
163 const DoFHandler<dim, spacedim> &mg_dof,
165 const MPI_Comm mpi_communicator,
166 const bool transfer_solution_vectors,
167 std::vector<Table<2, unsigned int>> &copy_indices,
168 std::vector<Table<2, unsigned int>> &copy_indices_global_mine,
169 std::vector<Table<2, unsigned int>> &copy_indices_level_mine,
170 LinearAlgebra::distributed::Vector<Number> &ghosted_global_vector,
172 &ghosted_level_vector)
173 {
174 // first go to the usual routine...
175 std::vector<
176 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
177 my_copy_indices;
178 std::vector<
179 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
180 my_copy_indices_global_mine;
181 std::vector<
182 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
183 my_copy_indices_level_mine;
184
186 mg_constrained_dofs,
187 my_copy_indices,
188 my_copy_indices_global_mine,
189 my_copy_indices_level_mine,
190 !transfer_solution_vectors);
191
192 // get all degrees of freedom that we need read access to in copy_to_mg
193 // and copy_from_mg, respectively. We fill an IndexSet once on each level
194 // (for the global_mine indices accessing remote level indices) and once
195 // globally (for the level_mine indices accessing remote global indices).
196
197 // the variables index_set and level_index_set are going to define the
198 // ghost indices of the respective vectors (due to construction, these are
199 // precisely the indices that we need)
200
201 IndexSet index_set(mg_dof.locally_owned_dofs().size());
202 std::vector<types::global_dof_index> accessed_indices;
203 ghosted_level_vector.resize(0,
205 1);
206 std::vector<IndexSet> level_index_set(
208 for (unsigned int l = 0; l < mg_dof.get_triangulation().n_global_levels();
209 ++l)
210 {
211 for (const auto &indices : my_copy_indices_level_mine[l])
212 accessed_indices.push_back(indices.first);
213 std::vector<types::global_dof_index> accessed_level_indices;
214 for (const auto &indices : my_copy_indices_global_mine[l])
215 accessed_level_indices.push_back(indices.second);
216 std::sort(accessed_level_indices.begin(), accessed_level_indices.end());
217 level_index_set[l].set_size(mg_dof.locally_owned_mg_dofs(l).size());
218 level_index_set[l].add_indices(accessed_level_indices.begin(),
219 accessed_level_indices.end());
220 level_index_set[l].compress();
221 ghosted_level_vector[l].reinit(mg_dof.locally_owned_mg_dofs(l),
222 level_index_set[l],
223 mpi_communicator);
224 }
225 std::sort(accessed_indices.begin(), accessed_indices.end());
226 index_set.add_indices(accessed_indices.begin(), accessed_indices.end());
227 index_set.compress();
228 ghosted_global_vector.reinit(mg_dof.locally_owned_dofs(),
229 index_set,
230 mpi_communicator);
231
232 // localize the copy indices for faster access. Since all access will be
233 // through the ghosted vector in 'data', we can use this (much faster)
234 // option
235 copy_indices.resize(mg_dof.get_triangulation().n_global_levels());
236 copy_indices_level_mine.resize(
238 copy_indices_global_mine.resize(
240 for (unsigned int level = 0;
241 level < mg_dof.get_triangulation().n_global_levels();
242 ++level)
243 {
244 const Utilities::MPI::Partitioner &global_partitioner =
245 *ghosted_global_vector.get_partitioner();
246 const Utilities::MPI::Partitioner &level_partitioner =
247 *ghosted_level_vector[level].get_partitioner();
248
249 auto translate_indices =
250 [&](const std::vector<
251 std::pair<types::global_dof_index, types::global_dof_index>>
252 &global_copy_indices,
253 Table<2, unsigned int> &local_copy_indices) {
254 local_copy_indices.reinit(2, global_copy_indices.size());
255 for (unsigned int i = 0; i < global_copy_indices.size(); ++i)
256 {
257 local_copy_indices(0, i) = global_partitioner.global_to_local(
258 global_copy_indices[i].first);
259 local_copy_indices(1, i) = level_partitioner.global_to_local(
260 global_copy_indices[i].second);
261 }
262 };
263
264 // owned-owned case
265 translate_indices(my_copy_indices[level], copy_indices[level]);
266
267 // remote-owned case
268 translate_indices(my_copy_indices_level_mine[level],
269 copy_indices_level_mine[level]);
270
271 // owned-remote case
272 translate_indices(my_copy_indices_global_mine[level],
273 copy_indices_global_mine[level]);
274 }
275 }
276} // namespace
277
278template <typename Number>
279template <int dim, int spacedim>
280void
282 fill_and_communicate_copy_indices(const DoFHandler<dim, spacedim> &mg_dof)
283{
284 const MPI_Comm mpi_communicator = mg_dof.get_mpi_communicator();
285
286 fill_internal(mg_dof,
287 mg_constrained_dofs,
288 mpi_communicator,
289 false,
290 this->copy_indices,
291 this->copy_indices_global_mine,
292 this->copy_indices_level_mine,
293 ghosted_global_vector,
294 ghosted_level_vector);
295
296 // in case we have hanging nodes which imply different ownership of
297 // global and level dof indices, we must also fill the solution indices
298 // which contain additional indices, otherwise they can re-use the
299 // indices of the standard transfer.
300 int have_refinement_edge_dofs = 0;
301 if (mg_constrained_dofs != nullptr)
302 for (unsigned int level = 0;
303 level < mg_dof.get_triangulation().n_global_levels();
304 ++level)
305 if (mg_constrained_dofs->get_refinement_edge_indices(level).n_elements() >
306 0)
307 {
308 have_refinement_edge_dofs = 1;
309 break;
310 }
311 if (Utilities::MPI::max(have_refinement_edge_dofs, mpi_communicator) == 1)
312 {
313 // note: variables not needed
314 std::vector<Table<2, unsigned int>> solution_copy_indices_global_mine;
316 solution_ghosted_level_vector;
317
318 fill_internal(mg_dof,
319 mg_constrained_dofs,
320 mpi_communicator,
321 true,
322 this->solution_copy_indices,
323 solution_copy_indices_global_mine,
324 this->solution_copy_indices_level_mine,
325 solution_ghosted_global_vector,
326 solution_ghosted_level_vector);
327 }
328 else
329 {
330 this->solution_copy_indices = this->copy_indices;
331 this->solution_copy_indices_level_mine = this->copy_indices_level_mine;
332 solution_ghosted_global_vector = ghosted_global_vector;
333 }
334
335 // Check if we can perform a cheaper "plain copy" (with or without
336 // renumbering) instead of having to translate individual entries
337 // using copy_indices*. This only works if a) we don't have to send
338 // or receive any DoFs and we have all locally owned DoFs in our
339 // copy_indices (so no adaptive refinement) and b) all processors
340 // agree on the choice (see below).
341 const bool my_perform_renumbered_plain_copy =
342 (this->copy_indices.back().n_cols() ==
343 mg_dof.locally_owned_dofs().n_elements()) &&
344 (this->copy_indices_global_mine.back().n_rows() == 0) &&
345 (this->copy_indices_level_mine.back().n_rows() == 0);
346
347 bool my_perform_plain_copy = false;
348 if (my_perform_renumbered_plain_copy)
349 {
350 my_perform_plain_copy = true;
351 // check whether there is a renumbering of degrees of freedom on
352 // either the finest level or the global dofs, which means that we
353 // cannot apply a plain copy
354 for (unsigned int i = 0; i < this->copy_indices.back().n_cols(); ++i)
355 if (this->copy_indices.back()(0, i) != this->copy_indices.back()(1, i))
356 {
357 my_perform_plain_copy = false;
358 break;
359 }
360 }
361
362 // now do a global reduction over all processors to see what operation
363 // they can agree upon
364 perform_plain_copy =
365 Utilities::MPI::min(static_cast<int>(my_perform_plain_copy),
366 mpi_communicator);
367 perform_renumbered_plain_copy =
368 Utilities::MPI::min(static_cast<int>(my_perform_renumbered_plain_copy),
369 mpi_communicator);
370
371 // if we do a plain copy, no need to hold additional ghosted vectors
372 if (perform_renumbered_plain_copy)
373 {
374 for (unsigned int i = 0; i < this->copy_indices.back().n_cols(); ++i)
375 AssertDimension(this->copy_indices.back()(0, i), i);
376
377 ghosted_global_vector.reinit(0);
378 ghosted_level_vector.resize(0, 0);
379 solution_ghosted_global_vector.reinit(0);
380 }
381}
382
383
384
385template <typename Number>
386void
388{
389 sizes.resize(0);
390 copy_indices.clear();
391 solution_copy_indices.clear();
392 copy_indices_global_mine.clear();
393 copy_indices_level_mine.clear();
394 solution_copy_indices_level_mine.clear();
395 component_to_block_map.resize(0);
396 mg_constrained_dofs = nullptr;
397 ghosted_global_vector.reinit(0);
398 solution_ghosted_global_vector.reinit(0);
399 ghosted_level_vector.resize(0, 0);
400 perform_plain_copy = false;
401 perform_renumbered_plain_copy = false;
402 initialize_dof_vector = nullptr;
403}
404
405
406
407template <typename Number>
408void
410 print_indices(std::ostream &os) const
411{
412 for (unsigned int level = 0; level < copy_indices.size(); ++level)
413 {
414 for (unsigned int i = 0; i < copy_indices[level].n_cols(); ++i)
415 os << "copy_indices[" << level << "]\t" << copy_indices[level](0, i)
416 << '\t' << copy_indices[level](1, i) << std::endl;
417 }
418
419 for (unsigned int level = 0; level < copy_indices_level_mine.size(); ++level)
420 {
421 for (unsigned int i = 0; i < copy_indices_level_mine[level].n_cols(); ++i)
422 os << "copy_ifrom [" << level << "]\t"
423 << copy_indices_level_mine[level](0, i) << '\t'
424 << copy_indices_level_mine[level](1, i) << std::endl;
425 }
426 for (unsigned int level = 0; level < copy_indices_global_mine.size(); ++level)
427 {
428 for (unsigned int i = 0; i < copy_indices_global_mine[level].n_cols();
429 ++i)
430 os << "copy_ito [" << level << "]\t"
431 << copy_indices_global_mine[level](0, i) << '\t'
432 << copy_indices_global_mine[level](1, i) << std::endl;
433 }
434}
435
436
437
438template <typename Number>
439std::size_t
442{
443 std::size_t result = sizeof(*this);
445 result += MemoryConsumption::memory_consumption(copy_indices);
446 result += MemoryConsumption::memory_consumption(copy_indices_global_mine);
447 result += MemoryConsumption::memory_consumption(copy_indices_level_mine);
448 result += ghosted_global_vector.memory_consumption();
449 for (unsigned int i = ghosted_level_vector.min_level();
450 i <= ghosted_level_vector.max_level();
451 ++i)
452 result += ghosted_level_vector[i].memory_consumption();
453
454 return result;
455}
456
457
458
459// explicit instantiation
460#include "mg_level_global_transfer.inst"
461
462// create an additional instantiation currently not supported by the automatic
463// template instantiation scheme
465
466
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
const Triangulation< dim, spacedim > & get_triangulation() const
const IndexSet & locally_owned_dofs() const
MPI_Comm get_mpi_communicator() const
size_type size() const
Definition index_set.h:1776
size_type n_elements() const
Definition index_set.h:1934
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
std::size_t memory_consumption() const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
void print_indices(std::ostream &os) const
std::size_t memory_consumption() const
void fill_and_communicate_copy_indices(const DoFHandler< dim, spacedim > &dof_handler)
virtual unsigned int n_global_levels() const
unsigned int global_to_local(const types::global_dof_index global_index) const
virtual void reinit(const IndexSet &locally_owned_indices, const IndexSet &ghost_indices, const MPI_Comm communicator) override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
Point< 2 > second
Definition grid_out.cc:4624
Point< 2 > first
Definition grid_out.cc:4623
unsigned int level
Definition grid_out.cc:4626
#define AssertDimension(dim1, dim2)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
void fill_copy_indices(const DoFHandler< dim, spacedim > &dof_handler, const MGConstrainedDoFs *mg_constrained_dofs, std::vector< std::vector< std::pair< types::global_dof_index, types::global_dof_index > > > &copy_indices, std::vector< std::vector< std::pair< types::global_dof_index, types::global_dof_index > > > &copy_indices_global_mine, std::vector< std::vector< std::pair< types::global_dof_index, types::global_dof_index > > > &copy_indices_level_mine, const bool skip_interface_dofs=true)