ryujin 2.1.1 revision f751393c604e228afd3e8516e9e236ec609caae0
offline_data.template.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
10#include "offline_data.h"
11#include "scratch_data.h"
12#include "sparse_matrix_simd.template.h" /* instantiate read_in */
13
14#include <deal.II/base/graph_coloring.h>
15#include <deal.II/base/parallel.h>
16#include <deal.II/base/work_stream.h>
17#include <deal.II/dofs/dof_renumbering.h>
18#include <deal.II/dofs/dof_tools.h>
19#include <deal.II/fe/fe_values.h>
20#include <deal.II/grid/grid_tools.h>
21#include <deal.II/lac/dynamic_sparsity_pattern.h>
22#include <deal.II/lac/la_parallel_vector.h>
23#ifdef DEAL_II_WITH_TRILINOS
24#include <deal.II/lac/trilinos_sparse_matrix.h>
25#endif
26
27#ifdef FORCE_DEAL_II_SPARSE_MATRIX
28#undef DEAL_II_WITH_TRILINOS
29#endif
30
31namespace ryujin
32{
33 using namespace dealii;
34
35
36 template <int dim, typename Number>
38 const MPI_Comm &mpi_communicator,
39 const Discretization<dim> &discretization,
40 const std::string &subsection /*= "OfflineData"*/)
41 : ParameterAcceptor(subsection)
42 , discretization_(&discretization)
43 , mpi_communicator_(mpi_communicator)
44 {
45 }
46
47
48 template <int dim, typename Number>
50 {
51 /*
52 * First, we set up the locally_relevant index set, determine (globally
53 * indexed) affine constraints and create a (globally indexed) sparsity
54 * pattern:
55 */
56
57 auto &dof_handler = *dof_handler_;
58 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
59
60 IndexSet locally_relevant;
61 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
62
63 affine_constraints_.reinit(locally_relevant);
64 DoFTools::make_hanging_node_constraints(dof_handler, affine_constraints_);
65
66#ifndef DEAL_II_WITH_TRILINOS
67 AssertThrow(affine_constraints_.n_constraints() == 0,
68 ExcMessage("ryujin was built without Trilinos support - no "
69 "hanging node support available"));
70#endif
72 /*
73 * Enforce periodic boundary conditions. We assume that the mesh is in
74 * "normal configuration".
75 */
76
77 const auto &periodic_faces =
78 discretization_->triangulation().get_periodic_face_map();
79
80 for (const auto &[left, value] : periodic_faces) {
81 const auto &[right, orientation] = value;
82
83 typename DoFHandler<dim>::cell_iterator dof_cell_left(
84 &left.first->get_triangulation(),
85 left.first->level(),
86 left.first->index(),
87 &dof_handler);
88
89 typename DoFHandler<dim>::cell_iterator dof_cell_right(
90 &right.first->get_triangulation(),
91 right.first->level(),
92 right.first->index(),
93 &dof_handler);
94
95 if constexpr (dim != 1 && std::is_same<Number, double>::value) {
96 DoFTools::make_periodicity_constraints(
97 dof_cell_left->face(left.second),
98 dof_cell_right->face(right.second),
99 affine_constraints_,
100 ComponentMask(),
101#if DEAL_II_VERSION_GTE(9, 6, 0)
102 orientation);
103#else
104 /* orientation */ orientation[0],
105 /* flip */ orientation[1],
106 /* rotation */ orientation[2]);
107#endif
108 } else {
109 AssertThrow(false, dealii::ExcNotImplemented());
110 __builtin_trap();
111 }
112 }
113
114 affine_constraints_.close();
115
116 sparsity_pattern_.reinit(
117 dof_handler.n_dofs(), dof_handler.n_dofs(), locally_relevant);
118#ifdef DEAL_II_WITH_TRILINOS
119 DoFTools::make_sparsity_pattern(
120 dof_handler, sparsity_pattern_, affine_constraints_, false);
121#else
122 /*
123 * In case we use dealii::SparseMatrix<Number> for assembly we need a
124 * sparsity pattern that also includes the full locally relevant -
125 * locally relevant coupling block. This gets thrown out again later,
126 * but nevertheless we have to add it.
127 */
129 dof_handler, sparsity_pattern_, affine_constraints_, false);
130#endif
131
132 /*
133 * We have to complete the local stencil to have consistent size over
134 * all MPI ranks. Otherwise, MPI synchronization in our
135 * SparseMatrixSIMD class will fail.
136 */
137
138 SparsityTools::distribute_sparsity_pattern(
139 sparsity_pattern_, locally_owned, mpi_communicator_, locally_relevant);
140 }
141
142
143 template <int dim, typename Number>
144 void OfflineData<dim, Number>::setup(const unsigned int problem_dimension)
145 {
146#ifdef DEBUG_OUTPUT
147 std::cout << "OfflineData<dim, Number>::setup()" << std::endl;
148#endif
149
150 /*
151 * Initialize dof handler:
152 */
153
154 const auto &triangulation = discretization_->triangulation();
155 if (!dof_handler_)
156 dof_handler_ = std::make_unique<dealii::DoFHandler<dim>>(triangulation);
157 auto &dof_handler = *dof_handler_;
158
159 dof_handler.distribute_dofs(discretization_->finite_element());
160
161 n_locally_owned_ = dof_handler.locally_owned_dofs().n_elements();
162
163 /*
164 * Renumbering:
165 */
166
167 DoFRenumbering::Cuthill_McKee(dof_handler);
168
169 /*
170 * Reorder all (individual) export indices at the beginning of the
171 * locally_internal index range to achieve a better packing:
172 *
173 * Note: This function might miss export indices that come from
174 * eliminating hanging node and periodicity constraints (which we do
175 * not know at this point because they depend on the renumbering...).
176 */
178 dof_handler, mpi_communicator_, n_locally_owned_, 1);
179
180 /*
181 * Group degrees of freedom that have the same stencil size in groups
182 * of multiples of the VectorizedArray<Number>::size().
183 *
184 * In order to determine the stencil size we have to create a first,
185 * temporary sparsity pattern:
186 */
187 create_constraints_and_sparsity_pattern();
188 n_locally_internal_ = DoFRenumbering::internal_range(
189 dof_handler, sparsity_pattern_, VectorizedArray<Number>::size());
190
191 /*
192 * Reorder all (strides of) locally internal indices that contain
193 * export indices to the start of the index range. This reordering
194 * preserves the binning introduced by
195 * DoFRenumbering::internal_range().
196 *
197 * Note: This function might miss export indices that come from
198 * eliminating hanging node and periodicity constraints (which we do
199 * not know at this point because they depend on the renumbering...).
200 * We therefore have to update n_export_indices_ later again.
201 */
202 n_export_indices_ =
204 mpi_communicator_,
205 n_locally_internal_,
206 VectorizedArray<Number>::size());
207
208 /*
209 * A small lambda to check for stride-level consistency of the internal
210 * index range:
211 */
212 const auto consistent_stride_range [[maybe_unused]] = [&]() {
213 constexpr auto group_size = VectorizedArray<Number>::size();
214 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
215 const auto offset = n_locally_owned_ != 0 ? *locally_owned.begin() : 0;
216
217 unsigned int group_row_length = 0;
218 unsigned int i = 0;
219 for (; i < n_locally_internal_; ++i) {
220 if (i % group_size == 0) {
221 group_row_length = sparsity_pattern_.row_length(offset + i);
222 } else {
223 if (group_row_length != sparsity_pattern_.row_length(offset + i)) {
224 break;
225 }
226 }
227 }
228 return i / group_size * group_size;
229 };
230
231 /*
232 * Create final sparsity pattern:
233 */
234
235 create_constraints_and_sparsity_pattern();
236
237 const auto mpi_allreduce_logical_or = [&](const bool local_value) {
238 std::function<bool(const bool &, const bool &)> comparator =
239 [](const bool &left, const bool &right) -> bool {
240 return left || right;
241 };
242
243 return Utilities::MPI::all_reduce(
244 local_value, mpi_communicator_, comparator);
245 };
246
247 /*
248 * We have to ensure that the locally internal numbering range is still
249 * consistent, meaning that all strides have the same stencil size.
250 * This property might not hold any more after the elimination
251 * procedure of constrained degrees of freedom (periodicity, or hanging
252 * node constraints). Therefore, the following little dance:
253 */
254
255 if (mpi_allreduce_logical_or(affine_constraints_.n_constraints() > 0)) {
256
257#if DEAL_II_VERSION_GTE(9, 5, 0)
258 if (mpi_allreduce_logical_or( //
259 consistent_stride_range() != n_locally_internal_)) {
260 /*
261 * In this case we try to fix up the numbering by pushing affected
262 * strides to the end and slightly lowering the n_locally_internal_
263 * marker.
264 */
265 n_locally_internal_ = DoFRenumbering::inconsistent_strides_last(
266 dof_handler,
267 sparsity_pattern_,
268 n_locally_internal_,
269 VectorizedArray<Number>::size());
270 create_constraints_and_sparsity_pattern();
271 n_locally_internal_ = consistent_stride_range();
272 }
273#endif
274 }
275
276 /*
277 * Check that after all the dof manipulation and setup we still end up
278 * with indices in [0, locally_internal) that have uniform stencil size
279 * within a stride.
280 */
281 Assert(consistent_stride_range() == n_locally_internal_,
282 dealii::ExcInternalError());
283
284 /*
285 * Set up partitioner:
286 */
287
288 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
289 Assert(n_locally_owned_ == locally_owned.n_elements(),
290 dealii::ExcInternalError());
291
292 IndexSet locally_relevant;
293 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
294 /* Enlarge the locally relevant set to include all additional couplings: */
295 {
296 IndexSet additional_dofs(dof_handler.n_dofs());
297 for (auto &entry : sparsity_pattern_)
298 if (!locally_relevant.is_element(entry.column())) {
299 Assert(locally_owned.is_element(entry.row()), ExcInternalError());
300 additional_dofs.add_index(entry.column());
301 }
302 additional_dofs.compress();
303 locally_relevant.add_indices(additional_dofs);
304 locally_relevant.compress();
305 }
306
307 n_locally_relevant_ = locally_relevant.n_elements();
308
309 scalar_partitioner_ = std::make_shared<dealii::Utilities::MPI::Partitioner>(
310 locally_owned, locally_relevant, mpi_communicator_);
311
312 vector_partitioner_ =
313 create_vector_partitioner(scalar_partitioner_, problem_dimension);
314
315 /*
316 * After elminiating periodicity and hanging node constraints we need
317 * to update n_export_indices_ again. This happens because we need to
318 * call export_indices_first() with incomplete information (missing
319 * eliminated degrees of freedom).
320 */
321 if (affine_constraints_.n_constraints() > 0) {
322 /*
323 * Recalculate n_export_indices_:
324 */
325 n_export_indices_ = 0;
326 for (const auto &it : scalar_partitioner_->import_indices())
327 if (it.second <= n_locally_internal_)
328 n_export_indices_ = std::max(n_export_indices_, it.second);
329
330 constexpr auto simd_length = VectorizedArray<Number>::size();
331 n_export_indices_ =
332 (n_export_indices_ + simd_length - 1) / simd_length * simd_length;
333 }
334
335#ifdef DEBUG
336 /* Check that n_export_indices_ is valid: */
337 unsigned int control = 0;
338 for (const auto &it : scalar_partitioner_->import_indices())
339 if (it.second <= n_locally_internal_)
340 control = std::max(control, it.second);
341
342 Assert(control <= n_export_indices_, ExcInternalError());
343 Assert(n_export_indices_ <= n_locally_internal_, ExcInternalError());
344#endif
345
346 /*
347 * Set up SIMD sparsity pattern in local numbering. Nota bene: The
348 * SparsityPatternSIMD::reinit() function will translates the pattern
349 * from global deal.II (typical) dof indexing to local indices.
350 */
351
352 sparsity_pattern_simd_.reinit(
353 n_locally_internal_, sparsity_pattern_, scalar_partitioner_);
354
355 /*
356 * Next we can (re)initialize all local matrices:
357 */
358
359 lumped_mass_matrix_.reinit(scalar_partitioner_);
360 lumped_mass_matrix_inverse_.reinit(scalar_partitioner_);
361
362 mass_matrix_.reinit(sparsity_pattern_simd_);
363 betaij_matrix_.reinit(sparsity_pattern_simd_);
364 cij_matrix_.reinit(sparsity_pattern_simd_);
365 }
366
367
368 template <int dim, typename Number>
369 void OfflineData<dim, Number>::assemble()
370 {
371#ifdef DEBUG_OUTPUT
372 std::cout << "OfflineData<dim, Number>::assemble()" << std::endl;
373#endif
374
375 auto &dof_handler = *dof_handler_;
376
377 measure_of_omega_ = 0.;
378
379#ifdef DEAL_II_WITH_TRILINOS
380 /* Variant using TrilinosWrappers::SparseMatrix with global numbering */
381
382 AffineConstraints<double> affine_constraints_assembly;
383 affine_constraints_assembly.reinit(affine_constraints_.get_local_lines());
384 for (auto line : affine_constraints_.get_lines()) {
385 affine_constraints_assembly.add_line(line.index);
386 for (auto entry : line.entries)
387 affine_constraints_assembly.add_entry(
388 line.index, entry.first, entry.second);
389 affine_constraints_assembly.set_inhomogeneity(line.index,
390 line.inhomogeneity);
391 }
392 affine_constraints_assembly.close();
393
394 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
395 TrilinosWrappers::SparsityPattern trilinos_sparsity_pattern;
396 trilinos_sparsity_pattern.reinit(
397 locally_owned, sparsity_pattern_, mpi_communicator_);
398
399 TrilinosWrappers::SparseMatrix mass_matrix_tmp;
400 TrilinosWrappers::SparseMatrix betaij_matrix_tmp;
401 std::array<TrilinosWrappers::SparseMatrix, dim> cij_matrix_tmp;
402
403 mass_matrix_tmp.reinit(trilinos_sparsity_pattern);
404 betaij_matrix_tmp.reinit(trilinos_sparsity_pattern);
405 for (auto &matrix : cij_matrix_tmp)
406 matrix.reinit(trilinos_sparsity_pattern);
407
408#else
409 /* Variant using deal.II SparseMatrix with local numbering */
410
411 AffineConstraints<Number> affine_constraints_assembly;
412 affine_constraints_assembly.copy_from(affine_constraints_);
413 transform_to_local_range(*scalar_partitioner_, affine_constraints_assembly);
414
415 SparsityPattern sparsity_pattern_assembly;
416 {
417 DynamicSparsityPattern dsp(n_locally_relevant_, n_locally_relevant_);
418 for (const auto &entry : sparsity_pattern_) {
419 const auto i = scalar_partitioner_->global_to_local(entry.row());
420 const auto j = scalar_partitioner_->global_to_local(entry.column());
421 dsp.add(i, j);
422 }
423 sparsity_pattern_assembly.copy_from(dsp);
424 }
425
426 dealii::SparseMatrix<Number> mass_matrix_tmp;
427 dealii::SparseMatrix<Number> betaij_matrix_tmp;
428 std::array<dealii::SparseMatrix<Number>, dim> cij_matrix_tmp;
429
430 mass_matrix_tmp.reinit(sparsity_pattern_assembly);
431 betaij_matrix_tmp.reinit(sparsity_pattern_assembly);
432 for (auto &matrix : cij_matrix_tmp)
433 matrix.reinit(sparsity_pattern_assembly);
434#endif
435
436 const unsigned int dofs_per_cell =
437 discretization_->finite_element().dofs_per_cell;
438
439 const unsigned int n_q_points = discretization_->quadrature().size();
440
441 /*
442 * Now, assemble all matrices:
443 */
444
445 /* The local, per-cell assembly routine: */
446
447 const auto local_assemble_system =
448 [&](const auto &cell, auto &scratch, auto &copy) {
449 /* iterate over locally owned cells and the ghost layer */
450
451 auto &is_locally_owned = copy.is_locally_owned_;
452 auto &local_dof_indices = copy.local_dof_indices_;
453
454 auto &cell_mass_matrix = copy.cell_mass_matrix_;
455 auto &cell_betaij_matrix = copy.cell_betaij_matrix_;
456 auto &cell_cij_matrix = copy.cell_cij_matrix_;
457 auto &cell_measure = copy.cell_measure_;
458
459 auto &fe_values = scratch.fe_values_;
460
461#ifdef DEAL_II_WITH_TRILINOS
462 is_locally_owned = cell->is_locally_owned();
463#else
464 /*
465 * When using a local dealii::SparseMatrix<Number> we don not
466 * have a compress(VectorOperation::add) available. In this case
467 * we assemble contributions over all locally relevant (non
468 * artificial) cells.
469 */
470 is_locally_owned = !cell->is_artificial();
471#endif
472 if (!is_locally_owned)
473 return;
474
475 cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
476 cell_betaij_matrix.reinit(dofs_per_cell, dofs_per_cell);
477 for (auto &matrix : cell_cij_matrix)
478 matrix.reinit(dofs_per_cell, dofs_per_cell);
479
480 fe_values.reinit(cell);
481
482 local_dof_indices.resize(dofs_per_cell);
483 cell->get_dof_indices(local_dof_indices);
484
485 /* clear out copy data: */
486 cell_mass_matrix = 0.;
487 cell_betaij_matrix = 0.;
488 for (auto &matrix : cell_cij_matrix)
489 matrix = 0.;
490 cell_measure = 0.;
491
492 for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
493 const auto JxW = fe_values.JxW(q_point);
494
495 if (cell->is_locally_owned())
496 cell_measure += Number(JxW);
497
498 for (unsigned int j = 0; j < dofs_per_cell; ++j) {
499 const auto value_JxW = fe_values.shape_value(j, q_point) * JxW;
500 const auto grad_JxW = fe_values.shape_grad(j, q_point) * JxW;
501
502 for (unsigned int i = 0; i < dofs_per_cell; ++i) {
503
504 const auto value = fe_values.shape_value(i, q_point);
505 const auto grad = fe_values.shape_grad(i, q_point);
506
507 cell_mass_matrix(i, j) += Number(value * value_JxW);
508 cell_betaij_matrix(i, j) += Number(grad * grad_JxW);
509 for (unsigned int d = 0; d < dim; ++d)
510 cell_cij_matrix[d](i, j) += Number((value * grad_JxW)[d]);
511
512 } /* for i */
513 } /* for j */
514 } /* for q */
515 };
516
517 const auto copy_local_to_global = [&](const auto &copy) {
518 const auto &is_locally_owned = copy.is_locally_owned_;
519 auto local_dof_indices = copy.local_dof_indices_; /* make a copy */
520 const auto &cell_mass_matrix = copy.cell_mass_matrix_;
521 const auto &cell_cij_matrix = copy.cell_cij_matrix_;
522 const auto &cell_betaij_matrix = copy.cell_betaij_matrix_;
523 const auto &cell_measure = copy.cell_measure_;
524
525 if (!is_locally_owned)
526 return;
527
528#ifndef DEAL_II_WITH_TRILINOS
529 transform_to_local_range(*scalar_partitioner_, local_dof_indices);
530#endif
531
532 affine_constraints_assembly.distribute_local_to_global(
533 cell_mass_matrix, local_dof_indices, mass_matrix_tmp);
534
535 for (int k = 0; k < dim; ++k) {
536 affine_constraints_assembly.distribute_local_to_global(
537 cell_cij_matrix[k], local_dof_indices, cij_matrix_tmp[k]);
538 }
539
540 affine_constraints_assembly.distribute_local_to_global(
541 cell_betaij_matrix, local_dof_indices, betaij_matrix_tmp);
542
543 measure_of_omega_ += cell_measure;
544 };
545
546 WorkStream::run(dof_handler.begin_active(),
547 dof_handler.end(),
548 local_assemble_system,
549 copy_local_to_global,
550 AssemblyScratchData<dim>(*discretization_),
551#ifdef DEAL_II_WITH_TRILINOS
552 AssemblyCopyData<dim, double>());
553#else
554 AssemblyCopyData<dim, Number>());
555#endif
556
557 measure_of_omega_ =
558 Utilities::MPI::sum(measure_of_omega_, mpi_communicator_);
559
560#ifdef DEAL_II_WITH_TRILINOS
561 betaij_matrix_tmp.compress(VectorOperation::add);
562 mass_matrix_tmp.compress(VectorOperation::add);
563 for (auto &it : cij_matrix_tmp)
564 it.compress(VectorOperation::add);
565#endif
566
567 /*
568 * Create lumped mass matrix:
569 */
570
571 {
572#ifdef DEAL_II_WITH_TRILINOS
573 using scalar_type = dealii::LinearAlgebra::distributed::Vector<double>;
574 scalar_type one(scalar_partitioner_);
575 one = 1.;
576
577 scalar_type local_lumped_mass_matrix(scalar_partitioner_);
578 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
579 lumped_mass_matrix_.compress(VectorOperation::add);
580
581 for (unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
582 ++i) {
583 lumped_mass_matrix_.local_element(i) =
584 local_lumped_mass_matrix.local_element(i);
585 lumped_mass_matrix_inverse_.local_element(i) =
586 1. / lumped_mass_matrix_.local_element(i);
587 }
588 lumped_mass_matrix_.update_ghost_values();
589 lumped_mass_matrix_inverse_.update_ghost_values();
590
591#else
592
593 Vector<Number> one(mass_matrix_tmp.m());
594 one = 1.;
595
596 Vector<Number> local_lumped_mass_matrix(mass_matrix_tmp.m());
597 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
598
599 for (unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
600 ++i) {
601 lumped_mass_matrix_.local_element(i) = local_lumped_mass_matrix(i);
602 lumped_mass_matrix_inverse_.local_element(i) =
603 1. / lumped_mass_matrix_.local_element(i);
604 }
605 lumped_mass_matrix_.update_ghost_values();
606 lumped_mass_matrix_inverse_.update_ghost_values();
607#endif
608 }
609
610#ifdef DEAL_II_WITH_TRILINOS
611 betaij_matrix_.read_in(betaij_matrix_tmp, /*locally_indexed*/ false);
612 mass_matrix_.read_in(mass_matrix_tmp, /*locally_indexed*/ false);
613 cij_matrix_.read_in(cij_matrix_tmp, /*locally_indexed*/ false);
614#else
615 betaij_matrix_.read_in(betaij_matrix_tmp, /*locally_indexed*/ true);
616 mass_matrix_.read_in(mass_matrix_tmp, /*locally_indexed*/ true);
617 cij_matrix_.read_in(cij_matrix_tmp, /*locally_indexed*/ true);
618#endif
619 betaij_matrix_.update_ghost_rows();
620 mass_matrix_.update_ghost_rows();
621 cij_matrix_.update_ghost_rows();
622
623 /* Populate boundary map and collect coupling boundary pairs: */
624
625 boundary_map_ = construct_boundary_map(
626 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
627
628 coupling_boundary_pairs_ = collect_coupling_boundary_pairs(
629 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
630 }
631
632
633 template <int dim, typename Number>
634 void OfflineData<dim, Number>::create_multigrid_data()
635 {
636#ifdef DEBUG_OUTPUT
637 std::cout << "OfflineData<dim, Number>::create_multigrid_data()"
638 << std::endl;
639#endif
640
641 auto &dof_handler = *dof_handler_;
642
643 dof_handler.distribute_mg_dofs();
644
645 const auto n_levels = dof_handler.get_triangulation().n_global_levels();
646
647 AffineConstraints<float> level_constraints;
648 // TODO not yet thread-parallel and without periodicity
649
650 level_boundary_map_.resize(n_levels);
651 level_lumped_mass_matrix_.resize(n_levels);
652
653 for (unsigned int level = 0; level < n_levels; ++level) {
654 /* Assemble lumped mass matrix vector: */
655
656 IndexSet relevant_dofs;
657 dealii::DoFTools::extract_locally_relevant_level_dofs(
658 dof_handler, level, relevant_dofs);
659 const auto partitioner = std::make_shared<Utilities::MPI::Partitioner>(
660 dof_handler.locally_owned_mg_dofs(level),
661 relevant_dofs,
662 lumped_mass_matrix_.get_mpi_communicator());
663 level_lumped_mass_matrix_[level].reinit(partitioner);
664 std::vector<types::global_dof_index> dof_indices(
665 dof_handler.get_fe().dofs_per_cell);
666 Vector<Number> mass_values(dof_handler.get_fe().dofs_per_cell);
667 FEValues<dim> fe_values(discretization_->mapping(),
668 discretization_->finite_element(),
669 discretization_->quadrature(),
670 update_values | update_JxW_values);
671 for (const auto &cell : dof_handler.cell_iterators_on_level(level))
672 // TODO for assembly with dealii::SparseMatrix and local
673 // numbering this probably has to read !cell->is_artificial()
674 if (cell->is_locally_owned_on_level()) {
675 fe_values.reinit(cell);
676 for (unsigned int i = 0; i < mass_values.size(); ++i) {
677 double sum = 0;
678 for (unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
679 sum += fe_values.shape_value(i, q) * fe_values.JxW(q);
680 mass_values(i) = sum;
681 }
682 cell->get_mg_dof_indices(dof_indices);
683 level_constraints.distribute_local_to_global(
684 mass_values, dof_indices, level_lumped_mass_matrix_[level]);
685 }
686 level_lumped_mass_matrix_[level].compress(VectorOperation::add);
687
688 /* Populate boundary map: */
689
690 level_boundary_map_[level] = construct_boundary_map(
691 dof_handler.begin_mg(level), dof_handler.end_mg(level), *partitioner);
692 }
693 }
694
695
696 template <int dim, typename Number>
697 template <typename ITERATOR1, typename ITERATOR2>
698 typename OfflineData<dim, Number>::boundary_map_type
700 const ITERATOR1 &begin,
701 const ITERATOR2 &end,
702 const Utilities::MPI::Partitioner &partitioner) const
703 {
704#ifdef DEBUG_OUTPUT
705 std::cout << "OfflineData<dim, Number>::construct_boundary_map()"
706 << std::endl;
707#endif
708
709 decltype(boundary_map_) preliminary_map;
710
711 std::vector<dealii::types::global_dof_index> local_dof_indices;
712
713 const dealii::QGauss<dim - 1> face_quadrature(3);
714 dealii::FEFaceValues<dim> fe_face_values(discretization_->mapping(),
715 discretization_->finite_element(),
716 face_quadrature,
717 dealii::update_normal_vectors |
718 dealii::update_values |
719 dealii::update_JxW_values);
720
721 const unsigned int dofs_per_cell =
722 discretization_->finite_element().dofs_per_cell;
723
724 const auto support_points =
725 discretization_->finite_element().get_unit_support_points();
726
727 for (auto cell = begin; cell != end; ++cell) {
728
729 if (!cell->is_locally_owned_on_level())
730 continue;
731
732 local_dof_indices.resize(dofs_per_cell);
733 cell->get_active_or_mg_dof_indices(local_dof_indices);
734
735 for (auto f : GeometryInfo<dim>::face_indices()) {
736 const auto face = cell->face(f);
737 const auto id = face->boundary_id();
738
739 if (!face->at_boundary())
740 continue;
741
742 /*
743 * Skip periodic boundary faces. For our algorithm these are
744 * interior degrees of freedom (if not simultaneously located at
745 * another boundary as well).
746 */
747 if (id == Boundary::periodic)
748 continue;
749
750 fe_face_values.reinit(cell, f);
751 const unsigned int n_face_q_points = face_quadrature.size();
752
753 for (unsigned int j = 0; j < dofs_per_cell; ++j) {
754
755 if (!discretization_->finite_element().has_support_on_face(j, f))
756 continue;
757
758 Number boundary_mass = 0.;
759 dealii::Tensor<1, dim, Number> normal;
760
761 for (unsigned int q = 0; q < n_face_q_points; ++q) {
762 const auto JxW = fe_face_values.JxW(q);
763 const auto phi_i = fe_face_values.shape_value(j, q);
764
765 boundary_mass += phi_i * JxW;
766 normal += phi_i * fe_face_values.normal_vector(q) * JxW;
767 }
768
769 const auto global_index = local_dof_indices[j];
770 const auto index = partitioner.global_to_local(global_index);
771
772 /* Skip nonlocal degrees of freedom: */
773 if (index >= n_locally_owned_)
774 continue;
775
776 /* Skip constrained degrees of freedom: */
777 const unsigned int row_length =
778 sparsity_pattern_simd_.row_length(index);
779 if (row_length == 1)
780 continue;
781
782 Point<dim> position =
783 discretization_->mapping().transform_unit_to_real_cell(
784 cell, support_points[j]);
785
786 /*
787 * Temporarily insert a (wrong) boundary mass value for the
788 * normal mass. We'll fix this later.
789 */
790 preliminary_map.insert(
791 {index, {normal, boundary_mass, boundary_mass, id, position}});
792 } /* j */
793 } /* f */
794 } /* cell */
795
796 /*
797 * Filter boundary map:
798 *
799 * At this point we have collected multiple cell contributions for each
800 * boundary degree of freedom. We now merge all entries that have the
801 * same boundary id and whose normals describe an acute angle of about
802 * 85 degrees or less.
803 *
804 * FIXME: is this robust in 3D?
805 */
806
807 decltype(boundary_map_) filtered_map;
808 std::set<dealii::types::global_dof_index> boundary_dofs;
809 for (auto entry : preliminary_map) {
810 bool inserted = false;
811 const auto range = filtered_map.equal_range(entry.first);
812 for (auto it = range.first; it != range.second; ++it) {
813 const auto &[new_normal,
814 new_normal_mass,
815 new_boundary_mass,
816 new_id,
817 new_point] = entry.second;
818 auto &[normal, normal_mass, boundary_mass, id, point] = it->second;
819
820 if (id != new_id)
821 continue;
822
823 Assert(point.distance(new_point) < 1.0e-14, dealii::ExcInternalError());
824
825 if (normal * new_normal / normal.norm() / new_normal.norm() > 0.08) {
826 /* Both normals describe an acute angle of 85 degrees or less. */
827 normal += new_normal;
828 boundary_mass += new_boundary_mass;
829 inserted = true;
830 }
831 }
832 if (!inserted)
833 filtered_map.insert(entry);
834 }
835
836 /* Normalize all normal vectors: */
837
838 for (auto &it : filtered_map) {
839 auto &[normal, normal_mass, boundary_mass, id, point] = it.second;
840 const auto new_normal_mass =
841 normal.norm() + std::numeric_limits<Number>::epsilon();
842 /* Replace boundary mass with new definition: */
843 normal_mass = new_normal_mass;
844 normal /= new_normal_mass;
845 }
846
847 return filtered_map;
848 }
849
850
851 template <int dim, typename Number>
852 template <typename ITERATOR1, typename ITERATOR2>
855 const ITERATOR1 &begin,
856 const ITERATOR2 &end,
857 const Utilities::MPI::Partitioner &partitioner) const
858 {
859#ifdef DEBUG_OUTPUT
860 std::cout << "OfflineData<dim, Number>::collect_coupling_boundary_pairs()"
861 << std::endl;
862#endif
863
864 /*
865 * First, collect *all* locally relevant degrees of freedom that are
866 * located on a (non periodic) boundary. We also collect constrained
867 * degrees of freedom for the time being (and filter afterwards).
868 */
869
870 std::set<unsigned int> locally_relevant_boundary_indices;
871
872 std::vector<dealii::types::global_dof_index> local_dof_indices;
873
874 const unsigned int dofs_per_cell =
875 discretization_->finite_element().dofs_per_cell;
876
877 for (auto cell = begin; cell != end; ++cell) {
878
879 /* Make sure to iterate over the entire locally relevant set: */
880 if (cell->is_artificial())
881 continue;
882
883 local_dof_indices.resize(dofs_per_cell);
884 cell->get_active_or_mg_dof_indices(local_dof_indices);
885
886 for (auto f : GeometryInfo<dim>::face_indices()) {
887 const auto face = cell->face(f);
888 const auto id = face->boundary_id();
889
890 if (!face->at_boundary())
891 continue;
892
893 /* Skip periodic boundary faces; see above. */
894 if (id == Boundary::periodic)
895 continue;
896
897 for (unsigned int j = 0; j < dofs_per_cell; ++j) {
898
899 if (!discretization_->finite_element().has_support_on_face(j, f))
900 continue;
901
902 const auto global_index = local_dof_indices[j];
903 const auto index = partitioner.global_to_local(global_index);
904
905 /* Skip irrelevant degrees of freedom: */
906 if (index >= n_locally_relevant_)
907 continue;
908
909 locally_relevant_boundary_indices.insert(index);
910 } /* j */
911 } /* f */
912 } /* cell */
913
914 /*
915 * Now, collect all coupling boundary pairs:
916 */
917
918 coupling_boundary_pairs_type result;
919
920 for (const auto i : locally_relevant_boundary_indices) {
921
922 /* Only record pairs with a left index that is locally owned: */
923 if (i >= n_locally_owned_)
924 continue;
925
926 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
927
928 /* Skip all constrained degrees of freedom: */
929 if (row_length == 1)
930 continue;
931
932 const unsigned int *js = sparsity_pattern_simd_.columns(i);
933 constexpr auto simd_length = VectorizedArray<Number>::size();
934 /* skip diagonal: */
935 for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
936 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
937 : js + col_idx);
938
939 if (locally_relevant_boundary_indices.count(j) != 0) {
940 result.push_back({i, col_idx, j});
941 }
942 }
943 }
944
945 return result;
946 }
947
948} /* namespace ryujin */
OfflineData(const MPI_Comm &mpi_communicator, const Discretization< dim > &discretization, const std::string &subsection="/OfflineData")
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)
std::shared_ptr< const dealii::Utilities::MPI::Partitioner > create_vector_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &scalar_partitioner, const unsigned int n_components)
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)