ryujin 2.1.1 revision 46bf70e400e423a8ffffe8300887eeb35b8dfb2c
local_index_handling.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2022 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
10#include <deal.II/base/partitioner.h>
11#include <deal.II/dofs/dof_handler.h>
12#include <deal.II/dofs/dof_renumbering.h>
13#include <deal.II/dofs/dof_tools.h>
14#include <deal.II/lac/affine_constraints.h>
15#include <deal.II/lac/dynamic_sparsity_pattern.h>
16#include <deal.II/lac/sparsity_tools.h>
17
18namespace ryujin
19{
28 template <typename Number>
30 const dealii::Utilities::MPI::Partitioner &partitioner,
31 dealii::AffineConstraints<Number> &affine_constraints)
32 {
33 affine_constraints.close();
34
35 dealii::AffineConstraints<Number> temporary;
36
37 for (auto line : affine_constraints.get_lines()) {
38 /* translate into local index ranges: */
39 line.index = partitioner.global_to_local(line.index);
40 std::transform(line.entries.begin(),
41 line.entries.end(),
42 line.entries.begin(),
43 [&](auto entry) {
44 return std::make_pair(
45 partitioner.global_to_local(entry.first),
46 entry.second);
47 });
48
49 temporary.add_line(line.index);
50 temporary.add_entries(line.index, line.entries);
51 temporary.set_inhomogeneity(line.index, line.inhomogeneity);
52 }
53
54 temporary.close();
55
56 affine_constraints = std::move(temporary);
57 }
58
59
68 template <typename VECTOR>
70 const dealii::Utilities::MPI::Partitioner &partitioner, VECTOR &vector)
71 {
72 std::transform(
73 vector.begin(), vector.end(), vector.begin(), [&](auto index) {
74 return partitioner.global_to_local(index);
75 });
76 }
77
78
85 namespace DoFRenumbering
86 {
91 using dealii::DoFRenumbering::Cuthill_McKee;
92
102 template <int dim>
103 unsigned int export_indices_first(dealii::DoFHandler<dim> &dof_handler,
104 const MPI_Comm &mpi_communicator,
105 const unsigned int n_locally_internal,
106 const std::size_t group_size)
107 {
108 using namespace dealii;
109
110 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
111 const auto n_locally_owned = locally_owned.n_elements();
112
113 /* The locally owned index range has to be contiguous */
114 Assert(locally_owned.is_contiguous() == true,
115 dealii::ExcMessage(
116 "Need a contiguous set of locally owned indices."));
117
118 /* Offset to translate from global to local index range */
119 const auto offset = n_locally_owned != 0 ? *locally_owned.begin() : 0;
120
121#if DEAL_II_VERSION_GTE(9, 6, 0)
122 const auto locally_relevant =
123 DoFTools::extract_locally_relevant_dofs(dof_handler);
124#else
125 IndexSet locally_relevant;
126 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
127#endif
128
129 /* Create a temporary MPI partitioner: */
130
131 Utilities::MPI::Partitioner partitioner(
132 locally_owned, locally_relevant, mpi_communicator);
133
134 IndexSet export_indices(n_locally_owned);
135 for (const auto &it : partitioner.import_indices()) {
136 export_indices.add_range(it.first, it.second);
137 }
138
139 std::vector<dealii::types::global_dof_index> new_order(n_locally_owned);
140
141 /*
142 * First pass: reorder all strides containing export indices and mark
143 * all other indices with numbers::invalid_dof_index:
144 */
145
146 unsigned int n_export_indices = 0;
147
148 Assert(n_locally_internal <= n_locally_owned, dealii::ExcInternalError());
149
150 for (unsigned int i = 0; i < n_locally_internal; i += group_size) {
151 bool export_index_present = false;
152 for (unsigned int j = 0; j < group_size; ++j) {
153 if (export_indices.is_element(i + j)) {
154 export_index_present = true;
155 break;
156 }
157 }
158
159 if (export_index_present) {
160 Assert(n_export_indices % group_size == 0,
161 dealii::ExcInternalError());
162 for (unsigned int j = 0; j < group_size; ++j) {
163 new_order[i + j] = offset + n_export_indices++;
164 }
165 } else {
166 for (unsigned int j = 0; j < group_size; ++j)
167 new_order[i + j] = dealii::numbers::invalid_dof_index;
168 }
169 }
170
171#if DEBUG
172 unsigned int n_other = 0;
173 for (unsigned int i = n_locally_internal; i < n_locally_owned; ++i)
174 if (export_indices.is_element(i))
175 n_other++;
176
177 Assert(n_other + n_export_indices >= export_indices.n_elements(),
178 dealii::ExcInternalError());
179#endif
180
181 unsigned int running_index = n_export_indices;
182
183 /*
184 * Second pass: append the rest:
185 */
186
187 for (unsigned int i = 0; i < n_locally_internal; i += group_size) {
188 if (new_order[i] == dealii::numbers::invalid_dof_index) {
189 for (unsigned int j = 0; j < group_size; ++j) {
190 Assert(new_order[i + j] == dealii::numbers::invalid_dof_index,
191 dealii::ExcInternalError());
192 new_order[i + j] = offset + running_index++;
193 }
194 }
195 }
196
197 Assert(running_index == n_locally_internal, dealii::ExcInternalError());
198
199 for (unsigned int i = n_locally_internal; i < n_locally_owned; i++) {
200 new_order[i] = offset + running_index++;
201 }
202
203 Assert(running_index == n_locally_owned, dealii::ExcInternalError());
204
205 dof_handler.renumber_dofs(new_order);
206
207 Assert(n_export_indices % group_size == 0, dealii::ExcInternalError());
208 Assert(n_export_indices <= n_locally_internal,
209 dealii::ExcInternalError());
210 return n_export_indices;
211 }
212
213
220 template <int dim>
221 unsigned int
222 inconsistent_strides_last(dealii::DoFHandler<dim> &dof_handler,
223 const dealii::DynamicSparsityPattern &sparsity,
224 const unsigned int n_locally_internal,
225 const std::size_t group_size)
226 {
227 using namespace dealii;
228
229 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
230 const auto n_locally_owned = locally_owned.n_elements();
231
232 /* The locally owned index range has to be contiguous */
233 Assert(locally_owned.is_contiguous() == true,
234 dealii::ExcMessage(
235 "Need a contiguous set of locally owned indices."));
236
237 /* Offset to translate from global to local index range */
238 const auto offset = n_locally_owned != 0 ? *locally_owned.begin() : 0;
239
240 std::vector<dealii::types::global_dof_index> new_order(n_locally_owned);
241
242 /*
243 * First pass: keep all strides with consistent row length at the
244 * beginning of the locally internal index range and mark all other
245 * indices with numbers::invalid_dof_index:
246 */
247
248 unsigned int n_consistent_range = 0;
249
250 Assert(n_locally_internal <= n_locally_owned, dealii::ExcInternalError());
251
252 for (unsigned int i = 0; i < n_locally_internal; i += group_size) {
253
254 bool stride_is_consistent = true;
255 const auto group_row_length = sparsity.row_length(offset + i);
256 for (unsigned int j = 0; j < group_size; ++j) {
257 if (group_row_length != sparsity.row_length(offset + i + j)) {
258 stride_is_consistent = false;
259 break;
260 }
261 }
262
263 if (stride_is_consistent) {
264 for (unsigned int j = 0; j < group_size; ++j) {
265 new_order[i + j] = offset + n_consistent_range++;
266 }
267 } else {
268 for (unsigned int j = 0; j < group_size; ++j)
269 new_order[i + j] = dealii::numbers::invalid_dof_index;
270 }
271 }
272
273 /*
274 * Second pass: append the rest:
275 */
276
277 unsigned int running_index = n_consistent_range;
278
279 for (unsigned int i = 0; i < n_locally_internal; i += group_size) {
280 if (new_order[i] == dealii::numbers::invalid_dof_index) {
281 for (unsigned int j = 0; j < group_size; ++j) {
282 Assert(new_order[i + j] == dealii::numbers::invalid_dof_index,
283 dealii::ExcInternalError());
284 new_order[i + j] = offset + running_index++;
285 }
286 }
287 }
288
289 Assert(running_index == n_locally_internal, dealii::ExcInternalError());
290
291 for (unsigned int i = n_locally_internal; i < n_locally_owned; i++) {
292 new_order[i] = offset + running_index++;
293 }
294
295 Assert(running_index == n_locally_owned, dealii::ExcInternalError());
296
297 dof_handler.renumber_dofs(new_order);
298
299 Assert(n_consistent_range % group_size == 0, dealii::ExcInternalError());
300 Assert(n_consistent_range <= n_locally_internal,
301 dealii::ExcInternalError());
302 return n_consistent_range;
303 }
304
305
320 template <int dim>
321 unsigned int internal_range(dealii::DoFHandler<dim> &dof_handler,
322 const dealii::DynamicSparsityPattern &sparsity,
323 const std::size_t group_size)
324 {
325 using namespace dealii;
326
327 const auto &locally_owned = dof_handler.locally_owned_dofs();
328 const auto n_locally_owned = locally_owned.n_elements();
329
330 /* The locally owned index range has to be contiguous */
331
332 Assert(locally_owned.is_contiguous() == true,
333 dealii::ExcMessage(
334 "Need a contiguous set of locally owned indices."));
335
336 /* Offset to translate from global to local index range */
337 const auto offset = n_locally_owned != 0 ? *locally_owned.begin() : 0;
338
339 using dof_type = dealii::types::global_dof_index;
340 std::vector<dof_type> new_order(n_locally_owned);
341 dof_type current_index = offset;
342
343 /*
344 * Sort degrees of freedom into a map grouped by stencil size. Write
345 * out dof indices into the new_order vector in groups of group_size
346 * and with same stencil size.
347 */
348
349 std::map<unsigned int, std::set<dof_type>> bins;
350
351 for (unsigned int i = 0; i < n_locally_owned; ++i) {
352 const dof_type index = i;
353 const unsigned int row_length = sparsity.row_length(offset + index);
354 bins[row_length].insert(index);
355
356 if (bins[row_length].size() == group_size) {
357 for (const auto &index : bins[row_length])
358 new_order[index] = current_index++;
359 bins.erase(row_length);
360 }
361 }
362
363 unsigned int n_locally_internal = current_index - offset;
364
365 /* Write out the rest. */
366
367 for (const auto &entries : bins) {
368 Assert(entries.second.size() > 0, ExcInternalError());
369 for (const auto &index : entries.second)
370 new_order[index] = current_index++;
371 }
372 Assert(current_index == offset + n_locally_owned, ExcInternalError());
373
374 dof_handler.renumber_dofs(new_order);
375
376 Assert(n_locally_internal % group_size == 0, ExcInternalError());
377 return n_locally_internal;
378 }
379 } // namespace DoFRenumbering
380
381
388 namespace DoFTools
389 {
391 using dealii::DoFTools::extract_locally_relevant_dofs;
392
394 using dealii::DoFTools::make_hanging_node_constraints;
395
397 using dealii::DoFTools::make_periodicity_constraints;
398
400 using dealii::DoFTools::make_sparsity_pattern;
401
402
410 template <int dim, typename Number, typename SPARSITY>
412 const dealii::DoFHandler<dim> &dof_handler,
413 SPARSITY &dsp,
414 const dealii::AffineConstraints<Number> &affine_constraints,
415 bool keep_constrained)
416 {
417 std::vector<dealii::types::global_dof_index> dof_indices;
418
419 for (auto cell : dof_handler.active_cell_iterators()) {
420 /* iterate over locally owned cells and the ghost layer */
421 if (cell->is_artificial())
422 continue;
423
424 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
425 dof_indices.resize(dofs_per_cell);
426 cell->get_dof_indices(dof_indices);
427
428 affine_constraints.add_entries_local_to_global(
429 dof_indices, dsp, keep_constrained);
430 }
431 }
432
433
442 template <int dim, typename Number, typename SPARSITY>
444 const dealii::DoFHandler<dim> &dof_handler,
445 SPARSITY &dsp,
446 const dealii::AffineConstraints<Number> &affine_constraints,
447 bool keep_constrained)
448 {
449 Assert(affine_constraints.n_constraints() == 0,
450 dealii::ExcMessage("I don't think constraints make sense for dG"));
451
452 std::vector<dealii::types::global_dof_index> dof_indices;
453 std::vector<dealii::types::global_dof_index> neighbor_dof_indices;
454
455 /*
456 * We collect all coupling dof indices on a face and store the result
457 * in a vector.
458 */
459 std::vector<dealii::types::global_dof_index> coupling_indices;
460 std::vector<dealii::types::global_dof_index> neighbor_coupling_indices;
461
462 /* we iterate over locally owned cells and the ghost layer */
463 for (auto cell : dof_handler.active_cell_iterators()) {
464 if (cell->is_artificial())
465 continue;
466
467 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
468 dof_indices.resize(dofs_per_cell);
469 cell->get_dof_indices(dof_indices);
470
471 affine_constraints.add_entries_local_to_global(
472 dof_indices, dsp, keep_constrained);
473
474 for (const auto f_index : cell->face_indices()) {
475 const auto &face = cell->face(f_index);
476
477 /* Skip faces without neighbors... */
478 const bool has_neighbor =
479 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
480 if (!has_neighbor)
481 continue;
482
483 /* Avoid artificial cells: */
484 const auto neighbor_cell =
485 cell->neighbor_or_periodic_neighbor(f_index);
486 if (neighbor_cell->is_artificial())
487 continue;
488
489 const unsigned int neighbor_dofs_per_cell =
490 neighbor_cell->get_fe().n_dofs_per_cell();
491 neighbor_dof_indices.resize(neighbor_dofs_per_cell);
492 neighbor_cell->get_dof_indices(neighbor_dof_indices);
493
494 const unsigned int f_index_neighbor =
495 cell->has_periodic_neighbor(f_index)
496 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
497 : cell->neighbor_of_neighbor(f_index);
498
499 /*
500 * Construct all couplings between current and neighbor cell with
501 * DoFs located at the boundary:
502 */
503
504 coupling_indices.resize(0);
505 for (unsigned int i = 0; i < dofs_per_cell; ++i)
506 if (cell->get_fe().has_support_on_face(i, f_index))
507 coupling_indices.push_back(dof_indices[i]);
508
509 neighbor_coupling_indices.resize(0);
510 for (unsigned int j = 0; j < neighbor_dofs_per_cell; ++j)
511 if (neighbor_cell->get_fe().has_support_on_face(j,
512 f_index_neighbor))
513 neighbor_coupling_indices.push_back(neighbor_dof_indices[j]);
514
515 affine_constraints.add_entries_local_to_global(
516 coupling_indices,
517 neighbor_coupling_indices,
518 dsp,
519 keep_constrained);
520 }
521 }
522 }
523
524
525 } // namespace DoFTools
526
527} // namespace ryujin
void make_extended_sparsity_pattern(const dealii::DoFHandler< dim > &dof_handler, SPARSITY &dsp, const dealii::AffineConstraints< Number > &affine_constraints, bool keep_constrained)
void transform_to_local_range(const dealii::Utilities::MPI::Partitioner &partitioner, dealii::AffineConstraints< Number > &affine_constraints)
unsigned int internal_range(dealii::DoFHandler< dim > &dof_handler, const dealii::DynamicSparsityPattern &sparsity, const std::size_t group_size)
unsigned int export_indices_first(dealii::DoFHandler< dim > &dof_handler, const MPI_Comm &mpi_communicator, const unsigned int n_locally_internal, const std::size_t group_size)
unsigned int inconsistent_strides_last(dealii::DoFHandler< dim > &dof_handler, const dealii::DynamicSparsityPattern &sparsity, const unsigned int n_locally_internal, const std::size_t group_size)
void make_extended_sparsity_pattern_dg(const dealii::DoFHandler< dim > &dof_handler, SPARSITY &dsp, const dealii::AffineConstraints< Number > &affine_constraints, bool keep_constrained)