Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\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\}}\)
mg_level_global_transfer.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2003 - 2021 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 
17 #include <deal.II/base/function.h>
18 #include <deal.II/base/logstream.h>
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 
48 template <typename VectorType>
49 template <int dim, int spacedim>
50 void
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
66  .locally_owned_mg_dofs(mg_dof.get_triangulation().n_global_levels() - 1)
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
87  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
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 
97 template <typename VectorType>
98 void
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 
112 template <typename VectorType>
113 void
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 
141 template <typename VectorType>
142 std::size_t
144 {
145  std::size_t result = sizeof(*this);
146  result += MemoryConsumption::memory_consumption(sizes);
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 
158 namespace
159 {
160  template <int dim, int spacedim, typename Number>
161  void
162  fill_internal(
163  const DoFHandler<dim, spacedim> &mg_dof,
164  SmartPointer<
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,
207  mg_dof.get_triangulation().n_global_levels() -
208  1);
209  std::vector<IndexSet> level_index_set(
210  mg_dof.get_triangulation().n_global_levels());
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(
240  mg_dof.get_triangulation().n_global_levels());
241  copy_indices_global_mine.resize(
242  mg_dof.get_triangulation().n_global_levels());
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 
281 template <typename Number>
282 template <int dim, int spacedim>
283 void
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 
384 template <typename Number>
385 void
387 {
388  sizes.resize(0);
389  copy_indices.clear();
390  copy_indices_global_mine.clear();
391  copy_indices_level_mine.clear();
392  component_to_block_map.resize(0);
393  mg_constrained_dofs = nullptr;
394  ghosted_global_vector.reinit(0);
395  ghosted_level_vector.resize(0, 0);
396  perform_plain_copy = false;
397  perform_renumbered_plain_copy = false;
398 }
399 
400 
401 
402 template <typename Number>
403 void
405  print_indices(std::ostream &os) const
406 {
407  for (unsigned int level = 0; level < copy_indices.size(); ++level)
408  {
409  for (unsigned int i = 0; i < copy_indices[level].n_cols(); ++i)
410  os << "copy_indices[" << level << "]\t" << copy_indices[level](0, i)
411  << '\t' << copy_indices[level](1, i) << std::endl;
412  }
413 
414  for (unsigned int level = 0; level < copy_indices_level_mine.size(); ++level)
415  {
416  for (unsigned int i = 0; i < copy_indices_level_mine[level].n_cols(); ++i)
417  os << "copy_ifrom [" << level << "]\t"
418  << copy_indices_level_mine[level](0, i) << '\t'
419  << copy_indices_level_mine[level](1, i) << std::endl;
420  }
421  for (unsigned int level = 0; level < copy_indices_global_mine.size(); ++level)
422  {
423  for (unsigned int i = 0; i < copy_indices_global_mine[level].n_cols();
424  ++i)
425  os << "copy_ito [" << level << "]\t"
426  << copy_indices_global_mine[level](0, i) << '\t'
427  << copy_indices_global_mine[level](1, i) << std::endl;
428  }
429 }
430 
431 
432 
433 template <typename Number>
434 std::size_t
437 {
438  std::size_t result = sizeof(*this);
439  result += MemoryConsumption::memory_consumption(sizes);
440  result += MemoryConsumption::memory_consumption(copy_indices);
441  result += MemoryConsumption::memory_consumption(copy_indices_global_mine);
442  result += MemoryConsumption::memory_consumption(copy_indices_level_mine);
443  result += ghosted_global_vector.memory_consumption();
444  for (unsigned int i = ghosted_level_vector.min_level();
445  i <= ghosted_level_vector.max_level();
446  ++i)
447  result += ghosted_level_vector[i].memory_consumption();
448 
449  return result;
450 }
451 
452 
453 
454 // explicit instantiation
455 #include "mg_level_global_transfer.inst"
456 
457 // create an additional instantiation currently not supported by the automatic
458 // template instantiation scheme
460 
461 
const Triangulation< dim, spacedim > & get_triangulation() const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
std::size_t memory_consumption() const
MPI_Comm get_communicator() const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Point< 2 > second
Definition: grid_out.cc:4587
unsigned int global_to_local(const types::global_dof_index global_index) const
void print_indices(std::ostream &os) const
size_type size() const
Definition: index_set.h:1639
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
unsigned int level
Definition: grid_out.cc:4589
const IndexSet & locally_owned_mg_dofs(const unsigned int level) const
Point< 2 > first
Definition: grid_out.cc:4586
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
T min(const T &t, const MPI_Comm &mpi_communicator)
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_partitioner() const
void fill_and_communicate_copy_indices(const DoFHandler< dim, spacedim > &dof_handler)
const IndexSet & locally_owned_dofs() const
virtual void reinit(const IndexSet &vector_space_vector_index_set, const IndexSet &read_write_vector_index_set, const MPI_Comm &communicator) override
Definition: partitioner.cc:134
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)
size_type n_elements() const
Definition: index_set.h:1837
T max(const T &t, const MPI_Comm &mpi_communicator)
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)