Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2003 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16
19
21
22#include <deal.II/fe/fe.h>
23
24#include <deal.II/grid/tria.h>
26
33#include <deal.II/lac/vector.h>
34
37#include <deal.II/multigrid/mg_transfer.templates.h>
39
40#include <algorithm>
41
43
44
45/* ------------------ MGLevelGlobalTransfer<VectorType> ----------------- */
46
47
48template <typename VectorType>
49template <int dim, int spacedim>
50void
52 const DoFHandler<dim, spacedim> &mg_dof)
53{
55 mg_constrained_dofs,
56 copy_indices,
57 copy_indices_global_mine,
58 copy_indices_level_mine);
59
60 // check if we can run a plain copy operation between the global DoFs and
61 // the finest level.
62 bool my_perform_plain_copy =
63 (copy_indices.back().size() == mg_dof.locally_owned_dofs().n_elements()) &&
64 (mg_dof.locally_owned_dofs().n_elements() ==
65 mg_dof
67 .n_elements());
68 if (my_perform_plain_copy)
69 {
70 AssertDimension(copy_indices_global_mine.back().size(), 0);
71 AssertDimension(copy_indices_level_mine.back().size(), 0);
72
73 // check whether there is a renumbering of degrees of freedom on
74 // either the finest level or the global dofs, which means that we
75 // cannot apply a plain copy
76 for (const auto &copy_index : copy_indices.back())
77 if (copy_index.first != copy_index.second)
78 {
79 my_perform_plain_copy = false;
80 break;
81 }
82 }
83
84 // now do a global reduction over all processors to see what operation
85 // they can agree upon
88 &mg_dof.get_triangulation()))
89 perform_plain_copy = (Utilities::MPI::min(my_perform_plain_copy ? 1 : 0,
90 ptria->get_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 MGConstrainedDoFs,
167 mg_constrained_dofs,
168 const MPI_Comm mpi_communicator,
169 const bool transfer_solution_vectors,
170 std::vector<Table<2, unsigned int>> & copy_indices,
171 std::vector<Table<2, unsigned int>> & copy_indices_global_mine,
172 std::vector<Table<2, unsigned int>> & copy_indices_level_mine,
173 LinearAlgebra::distributed::Vector<Number> &ghosted_global_vector,
175 &ghosted_level_vector)
176 {
177 // first go to the usual routine...
178 std::vector<
179 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
180 my_copy_indices;
181 std::vector<
182 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
183 my_copy_indices_global_mine;
184 std::vector<
185 std::vector<std::pair<types::global_dof_index, types::global_dof_index>>>
186 my_copy_indices_level_mine;
187
189 mg_constrained_dofs,
190 my_copy_indices,
191 my_copy_indices_global_mine,
192 my_copy_indices_level_mine,
193 !transfer_solution_vectors);
194
195 // get all degrees of freedom that we need read access to in copy_to_mg
196 // and copy_from_mg, respectively. We fill an IndexSet once on each level
197 // (for the global_mine indices accessing remote level indices) and once
198 // globally (for the level_mine indices accessing remote global indices).
199
200 // the variables index_set and level_index_set are going to define the
201 // ghost indices of the respective vectors (due to construction, these are
202 // precisely the indices that we need)
203
204 IndexSet index_set(mg_dof.locally_owned_dofs().size());
205 std::vector<types::global_dof_index> accessed_indices;
206 ghosted_level_vector.resize(0,
208 1);
209 std::vector<IndexSet> level_index_set(
211 for (unsigned int l = 0; l < mg_dof.get_triangulation().n_global_levels();
212 ++l)
213 {
214 for (const auto &indices : my_copy_indices_level_mine[l])
215 accessed_indices.push_back(indices.first);
216 std::vector<types::global_dof_index> accessed_level_indices;
217 for (const auto &indices : my_copy_indices_global_mine[l])
218 accessed_level_indices.push_back(indices.second);
219 std::sort(accessed_level_indices.begin(), accessed_level_indices.end());
220 level_index_set[l].set_size(mg_dof.locally_owned_mg_dofs(l).size());
221 level_index_set[l].add_indices(accessed_level_indices.begin(),
222 accessed_level_indices.end());
223 level_index_set[l].compress();
224 ghosted_level_vector[l].reinit(mg_dof.locally_owned_mg_dofs(l),
225 level_index_set[l],
226 mpi_communicator);
227 }
228 std::sort(accessed_indices.begin(), accessed_indices.end());
229 index_set.add_indices(accessed_indices.begin(), accessed_indices.end());
230 index_set.compress();
231 ghosted_global_vector.reinit(mg_dof.locally_owned_dofs(),
232 index_set,
233 mpi_communicator);
234
235 // localize the copy indices for faster access. Since all access will be
236 // through the ghosted vector in 'data', we can use this (much faster)
237 // option
238 copy_indices.resize(mg_dof.get_triangulation().n_global_levels());
239 copy_indices_level_mine.resize(
241 copy_indices_global_mine.resize(
243 for (unsigned int level = 0;
244 level < mg_dof.get_triangulation().n_global_levels();
245 ++level)
246 {
247 const Utilities::MPI::Partitioner &global_partitioner =
248 *ghosted_global_vector.get_partitioner();
249 const Utilities::MPI::Partitioner &level_partitioner =
250 *ghosted_level_vector[level].get_partitioner();
251
252 auto translate_indices =
253 [&](const std::vector<
254 std::pair<types::global_dof_index, types::global_dof_index>>
255 & global_copy_indices,
256 Table<2, unsigned int> &local_copy_indices) {
257 local_copy_indices.reinit(2, global_copy_indices.size());
258 for (unsigned int i = 0; i < global_copy_indices.size(); ++i)
259 {
260 local_copy_indices(0, i) = global_partitioner.global_to_local(
261 global_copy_indices[i].first);
262 local_copy_indices(1, i) = level_partitioner.global_to_local(
263 global_copy_indices[i].second);
264 }
265 };
266
267 // owned-owned case
268 translate_indices(my_copy_indices[level], copy_indices[level]);
269
270 // remote-owned case
271 translate_indices(my_copy_indices_level_mine[level],
272 copy_indices_level_mine[level]);
273
274 // owned-remote case
275 translate_indices(my_copy_indices_global_mine[level],
276 copy_indices_global_mine[level]);
277 }
278 }
279} // namespace
280
281template <typename Number>
282template <int dim, int spacedim>
283void
285 fill_and_communicate_copy_indices(const DoFHandler<dim, spacedim> &mg_dof)
286{
287 const MPI_Comm mpi_communicator = mg_dof.get_communicator();
288
289 fill_internal(mg_dof,
290 mg_constrained_dofs,
291 mpi_communicator,
292 false,
293 this->copy_indices,
294 this->copy_indices_global_mine,
295 this->copy_indices_level_mine,
296 ghosted_global_vector,
297 ghosted_level_vector);
298
299 // in case we have hanging nodes which imply different ownership of
300 // global and level dof indices, we must also fill the solution indices
301 // which contain additional indices, otherwise they can re-use the
302 // indices of the standard transfer.
303 int have_refinement_edge_dofs = 0;
304 if (mg_constrained_dofs != nullptr)
305 for (unsigned int level = 0;
306 level < mg_dof.get_triangulation().n_global_levels();
307 ++level)
308 if (mg_constrained_dofs->get_refinement_edge_indices(level).n_elements() >
309 0)
310 {
311 have_refinement_edge_dofs = 1;
312 break;
313 }
314 if (Utilities::MPI::max(have_refinement_edge_dofs, mpi_communicator) == 1)
315 fill_internal(mg_dof,
316 mg_constrained_dofs,
317 mpi_communicator,
318 true,
319 this->solution_copy_indices,
320 this->solution_copy_indices_global_mine,
321 this->solution_copy_indices_level_mine,
322 solution_ghosted_global_vector,
323 solution_ghosted_level_vector);
324 else
325 {
326 this->solution_copy_indices = this->copy_indices;
327 this->solution_copy_indices_global_mine = this->copy_indices_global_mine;
328 this->solution_copy_indices_level_mine = this->copy_indices_level_mine;
329 solution_ghosted_global_vector = ghosted_global_vector;
330 solution_ghosted_level_vector = ghosted_level_vector;
331 }
332
333 // Check if we can perform a cheaper "plain copy" (with or without
334 // renumbering) instead of having to translate individual entries
335 // using copy_indices*. This only works if a) we don't have to send
336 // or receive any DoFs and we have all locally owned DoFs in our
337 // copy_indices (so no adaptive refinement) and b) all processors
338 // agree on the choice (see below).
339 const bool my_perform_renumbered_plain_copy =
340 (this->copy_indices.back().n_cols() ==
341 mg_dof.locally_owned_dofs().n_elements()) &&
342 (this->copy_indices_global_mine.back().n_rows() == 0) &&
343 (this->copy_indices_level_mine.back().n_rows() == 0);
344
345 bool my_perform_plain_copy = false;
346 if (my_perform_renumbered_plain_copy)
347 {
348 my_perform_plain_copy = true;
349 // check whether there is a renumbering of degrees of freedom on
350 // either the finest level or the global dofs, which means that we
351 // cannot apply a plain copy
352 for (unsigned int i = 0; i < this->copy_indices.back().n_cols(); ++i)
353 if (this->copy_indices.back()(0, i) != this->copy_indices.back()(1, i))
354 {
355 my_perform_plain_copy = false;
356 break;
357 }
358 }
359
360 // now do a global reduction over all processors to see what operation
361 // they can agree upon
362 perform_plain_copy =
363 Utilities::MPI::min(static_cast<int>(my_perform_plain_copy),
364 mpi_communicator);
365 perform_renumbered_plain_copy =
366 Utilities::MPI::min(static_cast<int>(my_perform_renumbered_plain_copy),
367 mpi_communicator);
368
369 // if we do a plain copy, no need to hold additional ghosted vectors
370 if (perform_renumbered_plain_copy)
371 {
372 for (unsigned int i = 0; i < this->copy_indices.back().n_cols(); ++i)
373 AssertDimension(this->copy_indices.back()(0, i), i);
374
375 ghosted_global_vector.reinit(0);
376 ghosted_level_vector.resize(0, 0);
377 solution_ghosted_global_vector.reinit(0);
378 solution_ghosted_level_vector.resize(0, 0);
379 }
380}
381
382
383
384template <typename Number>
385void
387{
388 sizes.resize(0);
389 copy_indices.clear();
390 solution_copy_indices.clear();
391 copy_indices_global_mine.clear();
392 solution_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 solution_ghosted_level_vector.resize(0, 0);
401 perform_plain_copy = false;
402 perform_renumbered_plain_copy = false;
403 initialize_dof_vector = nullptr;
404}
405
406
407
408template <typename Number>
409void
411 print_indices(std::ostream &os) const
412{
413 for (unsigned int level = 0; level < copy_indices.size(); ++level)
414 {
415 for (unsigned int i = 0; i < copy_indices[level].n_cols(); ++i)
416 os << "copy_indices[" << level << "]\t" << copy_indices[level](0, i)
417 << '\t' << copy_indices[level](1, i) << std::endl;
418 }
419
420 for (unsigned int level = 0; level < copy_indices_level_mine.size(); ++level)
421 {
422 for (unsigned int i = 0; i < copy_indices_level_mine[level].n_cols(); ++i)
423 os << "copy_ifrom [" << level << "]\t"
424 << copy_indices_level_mine[level](0, i) << '\t'
425 << copy_indices_level_mine[level](1, i) << std::endl;
426 }
427 for (unsigned int level = 0; level < copy_indices_global_mine.size(); ++level)
428 {
429 for (unsigned int i = 0; i < copy_indices_global_mine[level].n_cols();
430 ++i)
431 os << "copy_ito [" << level << "]\t"
432 << copy_indices_global_mine[level](0, i) << '\t'
433 << copy_indices_global_mine[level](1, i) << std::endl;
434 }
435}
436
437
438
439template <typename Number>
440std::size_t
443{
444 std::size_t result = sizeof(*this);
446 result += MemoryConsumption::memory_consumption(copy_indices);
447 result += MemoryConsumption::memory_consumption(copy_indices_global_mine);
448 result += MemoryConsumption::memory_consumption(copy_indices_level_mine);
449 result += ghosted_global_vector.memory_consumption();
450 for (unsigned int i = ghosted_level_vector.min_level();
451 i <= ghosted_level_vector.max_level();
452 ++i)
453 result += ghosted_level_vector[i].memory_consumption();
454
455 return result;
456}
457
458
459
460// explicit instantiation
461#include "mg_level_global_transfer.inst"
462
463// create an additional instantiation currently not supported by the automatic
464// template instantiation scheme
466
467
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_communicator() const
size_type size() const
Definition index_set.h:1661
size_type n_elements() const
Definition index_set.h:1816
virtual std::size_t memory_consumption() const override
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() 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 &vector_space_vector_index_set, const IndexSet &read_write_vector_index_set, const MPI_Comm communicator) override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 2 > second
Definition grid_out.cc:4616
Point< 2 > first
Definition grid_out.cc:4615
unsigned int level
Definition grid_out.cc:4618
#define AssertDimension(dim1, dim2)
std::enable_if_t< std::is_fundamental< T >::value, 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)