15#include <deal.II/base/graph_coloring.h>
16#include <deal.II/base/parallel.h>
17#include <deal.II/base/work_stream.h>
18#include <deal.II/dofs/dof_renumbering.h>
19#include <deal.II/dofs/dof_tools.h>
20#include <deal.II/fe/fe_values.h>
21#include <deal.II/grid/grid_tools.h>
22#include <deal.II/lac/dynamic_sparsity_pattern.h>
23#include <deal.II/lac/la_parallel_vector.h>
24#ifdef DEAL_II_WITH_TRILINOS
25#include <deal.II/lac/trilinos_sparse_matrix.h>
28#ifdef FORCE_DEAL_II_SPARSE_MATRIX
29#undef DEAL_II_WITH_TRILINOS
34 using namespace dealii;
37 template <
int dim,
typename Number>
41 const std::string &subsection )
42 : ParameterAcceptor(subsection)
43 , mpi_ensemble_(mpi_ensemble)
44 , discretization_(&discretization)
46 incidence_relaxation_even_ = 0.5;
47 add_parameter(
"incidence matrix relaxation even degree",
48 incidence_relaxation_even_,
49 "Scaling exponent for incidence matrix used for "
50 "discontinuous finite elements with even degree. The default "
51 "value 0.5 scales the jump penalization with (h_i+h_j)^0.5.");
53 incidence_relaxation_odd_ = 0.0;
54 add_parameter(
"incidence matrix relaxation odd degree",
55 incidence_relaxation_odd_,
56 "Scaling exponent for incidence matrix used for "
57 "discontinuous finite elements with even degree. The default "
58 "value of 0.0 sets the jump penalization to a constant 1.");
62 template <
int dim,
typename Number>
71 auto &dof_handler = *dof_handler_;
72 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
74 IndexSet locally_relevant;
75 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
77 affine_constraints_.reinit(locally_relevant);
78 DoFTools::make_hanging_node_constraints(dof_handler, affine_constraints_);
80#ifndef DEAL_II_WITH_TRILINOS
81 AssertThrow(affine_constraints_.n_constraints() == 0,
82 ExcMessage(
"ryujin was built without Trilinos support - no "
83 "hanging node support available"));
91 const auto &periodic_faces =
92 discretization_->triangulation().get_periodic_face_map();
94 for (
const auto &[left, value] : periodic_faces) {
95 const auto &[right, orientation] = value;
97 typename DoFHandler<dim>::cell_iterator dof_cell_left(
98 &left.first->get_triangulation(),
103 typename DoFHandler<dim>::cell_iterator dof_cell_right(
104 &right.first->get_triangulation(),
105 right.first->level(),
106 right.first->index(),
109 if constexpr (std::is_same_v<Number, double>) {
110 DoFTools::make_periodicity_constraints(
111 dof_cell_left->face(left.second),
112 dof_cell_right->face(right.second),
115#
if DEAL_II_VERSION_GTE(9, 6, 0)
123 AssertThrow(
false, dealii::ExcNotImplemented());
128 affine_constraints_.close();
133 const std::vector<IndexSet> &locally_owned_dofs =
134 Utilities::MPI::all_gather(mpi_ensemble_.ensemble_communicator(),
135 dof_handler.locally_owned_dofs());
136 const IndexSet locally_active =
137 dealii::DoFTools::extract_locally_active_dofs(dof_handler);
138 Assert(affine_constraints_.is_consistent_in_parallel(
141 mpi_ensemble_.ensemble_communicator(),
147 sparsity_pattern_.reinit(
148 dof_handler.n_dofs(), dof_handler.n_dofs(), locally_relevant);
150 if (discretization_->have_discontinuous_ansatz()) {
155 dof_handler, sparsity_pattern_, affine_constraints_,
false);
160#ifdef DEAL_II_WITH_TRILINOS
161 DoFTools::make_sparsity_pattern(
162 dof_handler, sparsity_pattern_, affine_constraints_,
false);
171 dof_handler, sparsity_pattern_, affine_constraints_,
false);
181 SparsityTools::distribute_sparsity_pattern(
184 mpi_ensemble_.ensemble_communicator(),
189 template <
int dim,
typename Number>
190 void OfflineData<dim, Number>::setup(
const unsigned int problem_dimension,
191 const unsigned int n_precomputed_values)
194 std::cout <<
"OfflineData<dim, Number>::setup()" << std::endl;
201 const auto &triangulation = discretization_->triangulation();
203 dof_handler_ = std::make_unique<dealii::DoFHandler<dim>>(triangulation);
204 auto &dof_handler = *dof_handler_;
206 dof_handler.distribute_dofs(discretization_->finite_element());
208 n_locally_owned_ = dof_handler.locally_owned_dofs().n_elements();
214 DoFRenumbering::Cuthill_McKee(dof_handler);
225 mpi_ensemble_.ensemble_communicator(),
236 create_constraints_and_sparsity_pattern();
238 dof_handler, sparsity_pattern_, VectorizedArray<Number>::size());
253 mpi_ensemble_.ensemble_communicator(),
255 VectorizedArray<Number>::size());
261 const auto consistent_stride_range [[maybe_unused]] = [&]() {
262 constexpr auto group_size = VectorizedArray<Number>::size();
263 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
264 const auto offset = n_locally_owned_ != 0 ? *locally_owned.begin() : 0;
266 unsigned int group_row_length = 0;
268 for (; i < n_locally_internal_; ++i) {
269 if (i % group_size == 0) {
270 group_row_length = sparsity_pattern_.row_length(offset + i);
272 if (group_row_length != sparsity_pattern_.row_length(offset + i)) {
277 return i / group_size * group_size;
283 const auto mpi_allreduce_logical_or = [&](
const bool local_value) {
284 std::function<bool(
const bool &,
const bool &)> comparator =
285 [](
const bool &left,
const bool &right) ->
bool {
286 return left || right;
288 return Utilities::MPI::all_reduce(
289 local_value, mpi_ensemble_.ensemble_communicator(), comparator);
296 create_constraints_and_sparsity_pattern();
306#if DEAL_II_VERSION_GTE(9, 5, 0)
307 if (mpi_allreduce_logical_or(affine_constraints_.n_constraints() > 0)) {
308 if (mpi_allreduce_logical_or(
309 consistent_stride_range() != n_locally_internal_)) {
319 VectorizedArray<Number>::size());
320 create_constraints_and_sparsity_pattern();
321 n_locally_internal_ = consistent_stride_range();
331 Assert(consistent_stride_range() == n_locally_internal_,
332 dealii::ExcInternalError());
338 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
339 Assert(n_locally_owned_ == locally_owned.n_elements(),
340 dealii::ExcInternalError());
342 IndexSet locally_relevant;
343 DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant);
346 IndexSet additional_dofs(dof_handler.n_dofs());
347 for (
auto &entry : sparsity_pattern_)
348 if (!locally_relevant.is_element(entry.column())) {
349 Assert(locally_owned.is_element(entry.row()), ExcInternalError());
350 additional_dofs.add_index(entry.column());
352 additional_dofs.compress();
353 locally_relevant.add_indices(additional_dofs);
354 locally_relevant.compress();
357 n_locally_relevant_ = locally_relevant.n_elements();
359 scalar_partitioner_ = std::make_shared<dealii::Utilities::MPI::Partitioner>(
360 locally_owned, locally_relevant, mpi_ensemble_.ensemble_communicator());
363 scalar_partitioner_, problem_dimension);
366 scalar_partitioner_, n_precomputed_values);
374 if (mpi_allreduce_logical_or(affine_constraints_.n_constraints() > 0)) {
378 n_export_indices_ = 0;
379 for (
const auto &it : scalar_partitioner_->import_indices())
380 if (it.second <= n_locally_internal_)
381 n_export_indices_ = std::max(n_export_indices_, it.second);
383 constexpr auto simd_length = VectorizedArray<Number>::size();
385 (n_export_indices_ + simd_length - 1) / simd_length * simd_length;
390 unsigned int control = 0;
391 for (
const auto &it : scalar_partitioner_->import_indices())
392 if (it.second <= n_locally_internal_)
393 control = std::max(control, it.second);
395 Assert(control <= n_export_indices_, ExcInternalError());
396 Assert(n_export_indices_ <= n_locally_internal_, ExcInternalError());
405 sparsity_pattern_simd_.reinit(
406 n_locally_internal_, sparsity_pattern_, scalar_partitioner_);
410 template <
int dim,
typename Number>
411 void OfflineData<dim, Number>::create_matrices()
414 std::cout <<
"OfflineData<dim, Number>::create_matrices()" << std::endl;
421 mass_matrix_.reinit(sparsity_pattern_simd_);
422 if (discretization_->have_discontinuous_ansatz())
423 mass_matrix_inverse_.reinit(sparsity_pattern_simd_);
425 lumped_mass_matrix_.reinit(scalar_partitioner_);
426 lumped_mass_matrix_inverse_.reinit(scalar_partitioner_);
428 cij_matrix_.reinit(sparsity_pattern_simd_);
429 if (discretization_->have_discontinuous_ansatz())
430 incidence_matrix_.reinit(sparsity_pattern_simd_);
436 auto &dof_handler = *dof_handler_;
438 measure_of_omega_ = 0.;
440#ifdef DEAL_II_WITH_TRILINOS
443 AffineConstraints<double> affine_constraints_assembly;
445 affine_constraints_assembly.reinit(affine_constraints_.get_local_lines());
446 for (
auto line : affine_constraints_.get_lines()) {
447 affine_constraints_assembly.add_line(line.index);
448 for (
auto entry : line.entries)
449 affine_constraints_assembly.add_entry(
450 line.index, entry.first, entry.second);
451 affine_constraints_assembly.set_inhomogeneity(line.index,
454 affine_constraints_assembly.close();
456 const IndexSet &locally_owned = dof_handler.locally_owned_dofs();
457 TrilinosWrappers::SparsityPattern trilinos_sparsity_pattern;
458 trilinos_sparsity_pattern.reinit(locally_owned,
460 mpi_ensemble_.ensemble_communicator());
462 TrilinosWrappers::SparseMatrix mass_matrix_tmp;
463 TrilinosWrappers::SparseMatrix mass_matrix_inverse_tmp;
464 if (discretization_->have_discontinuous_ansatz())
465 mass_matrix_inverse_tmp.reinit(trilinos_sparsity_pattern);
467 std::array<TrilinosWrappers::SparseMatrix, dim> cij_matrix_tmp;
469 mass_matrix_tmp.reinit(trilinos_sparsity_pattern);
470 for (
auto &matrix : cij_matrix_tmp)
471 matrix.reinit(trilinos_sparsity_pattern);
476 AffineConstraints<Number> affine_constraints_assembly;
477 affine_constraints_assembly.copy_from(affine_constraints_);
480 SparsityPattern sparsity_pattern_assembly;
482 DynamicSparsityPattern dsp(n_locally_relevant_, n_locally_relevant_);
483 for (
const auto &entry : sparsity_pattern_) {
484 const auto i = scalar_partitioner_->global_to_local(entry.row());
485 const auto j = scalar_partitioner_->global_to_local(entry.column());
488 sparsity_pattern_assembly.copy_from(dsp);
491 dealii::SparseMatrix<Number> mass_matrix_tmp;
492 dealii::SparseMatrix<Number> mass_matrix_inverse_tmp;
493 if (discretization_->have_discontinuous_ansatz())
494 mass_matrix_inverse_tmp.reinit(sparsity_pattern_assembly);
496 std::array<dealii::SparseMatrix<Number>, dim> cij_matrix_tmp;
498 mass_matrix_tmp.reinit(sparsity_pattern_assembly);
499 for (
auto &matrix : cij_matrix_tmp)
500 matrix.reinit(sparsity_pattern_assembly);
508 const auto local_assemble_system = [&](
const auto &cell,
513 auto &is_locally_owned = copy.is_locally_owned_;
514 auto &local_dof_indices = copy.local_dof_indices_;
515 auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
517 auto &cell_mass_matrix = copy.cell_mass_matrix_;
518 auto &cell_mass_matrix_inverse = copy.cell_mass_matrix_inverse_;
519 auto &cell_cij_matrix = copy.cell_cij_matrix_;
520 auto &interface_cij_matrix = copy.interface_cij_matrix_;
521 auto &cell_measure = copy.cell_measure_;
523 auto &hp_fe_values = scratch.hp_fe_values_;
524 auto &hp_fe_face_values = scratch.hp_fe_face_values_;
525 auto &hp_fe_neighbor_face_values = scratch.hp_fe_neighbor_face_values_;
527#ifdef DEAL_II_WITH_TRILINOS
528 is_locally_owned = cell->is_locally_owned();
536 is_locally_owned = !cell->is_artificial();
538 if (!is_locally_owned)
541 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
543 cell_mass_matrix.reinit(dofs_per_cell, dofs_per_cell);
544 for (
auto &matrix : cell_cij_matrix)
545 matrix.reinit(dofs_per_cell, dofs_per_cell);
546 if (discretization_->have_discontinuous_ansatz()) {
547 cell_mass_matrix_inverse.reinit(dofs_per_cell, dofs_per_cell);
548 for (
auto &it : interface_cij_matrix)
549 for (
auto &matrix : it)
550 matrix.reinit(dofs_per_cell, dofs_per_cell);
553 hp_fe_values.reinit(cell);
554 const auto &fe_values = hp_fe_values.get_present_fe_values();
556 local_dof_indices.resize(dofs_per_cell);
557 cell->get_dof_indices(local_dof_indices);
561 cell_mass_matrix = 0.;
562 for (
auto &matrix : cell_cij_matrix)
564 if (discretization_->have_discontinuous_ansatz()) {
565 cell_mass_matrix_inverse = 0.;
566 for (
auto &it : interface_cij_matrix)
567 for (
auto &matrix : it)
572 for (
unsigned int q : fe_values.quadrature_point_indices()) {
573 const auto JxW = fe_values.JxW(q);
575 if (cell->is_locally_owned())
576 cell_measure += Number(JxW);
578 for (
unsigned int j : fe_values.dof_indices()) {
579 const auto value_JxW = fe_values.shape_value(j, q) * JxW;
580 const auto grad_JxW = fe_values.shape_grad(j, q) * JxW;
582 for (
unsigned int i : fe_values.dof_indices()) {
583 const auto value = fe_values.shape_value(i, q);
585 cell_mass_matrix(i, j) += Number(value * value_JxW);
586 for (
unsigned int d = 0; d < dim; ++d)
587 cell_cij_matrix[d](i, j) += Number((value * grad_JxW)[d]);
597 if (!discretization_->have_discontinuous_ansatz())
600 for (
const auto f_index : cell->face_indices()) {
601 const auto &face = cell->face(f_index);
604 const bool has_neighbor =
605 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
609 neighbor_local_dof_indices[f_index].resize(0);
614 const auto neighbor_cell = cell->neighbor_or_periodic_neighbor(f_index);
615 if (neighbor_cell->is_artificial()) {
618 neighbor_local_dof_indices[f_index].resize(0);
622 hp_fe_face_values.reinit(cell, f_index);
623 const auto &fe_face_values = hp_fe_face_values.get_present_fe_values();
627 for (
unsigned int q : fe_face_values.quadrature_point_indices()) {
628 const auto JxW = fe_face_values.JxW(q);
629 const auto &normal = fe_face_values.get_normal_vectors()[q];
631 for (
unsigned int j : fe_face_values.dof_indices()) {
632 const auto value_JxW = fe_face_values.shape_value(j, q) * JxW;
634 for (
unsigned int i : fe_face_values.dof_indices()) {
635 const auto value = fe_face_values.shape_value(i, q);
637 for (
unsigned int d = 0; d < dim; ++d)
638 cell_cij_matrix[d](i, j) -=
639 Number(0.5 * normal[d] * value * value_JxW);
646 const unsigned int f_index_neighbor =
647 cell->has_periodic_neighbor(f_index)
648 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
649 : cell->neighbor_of_neighbor(f_index);
651 neighbor_local_dof_indices[f_index].resize(dofs_per_cell);
652 neighbor_cell->get_dof_indices(neighbor_local_dof_indices[f_index]);
654 hp_fe_neighbor_face_values.reinit(neighbor_cell, f_index_neighbor);
655 const auto &fe_neighbor_face_values =
656 hp_fe_neighbor_face_values.get_present_fe_values();
658 for (
unsigned int q : fe_face_values.quadrature_point_indices()) {
659 const auto JxW = fe_face_values.JxW(q);
660 const auto &normal = fe_face_values.get_normal_vectors()[q];
663 for (
unsigned int j : fe_face_values.dof_indices()) {
664 const auto value_JxW =
665 fe_neighbor_face_values.shape_value(j, q) * JxW;
667 for (
unsigned int i : fe_face_values.dof_indices()) {
668 const auto value = fe_face_values.shape_value(i, q);
670 for (
unsigned int d = 0; d < dim; ++d)
671 interface_cij_matrix[f_index][d](i, j) +=
672 Number(0.5 * normal[d] * value * value_JxW);
682 if (discretization_->have_discontinuous_ansatz()) {
684 cell_mass_matrix_inverse.invert(cell_mass_matrix);
688 const auto copy_local_to_global = [&](
const auto ©) {
689 const auto &is_locally_owned = copy.is_locally_owned_;
690#ifdef DEAL_II_WITH_TRILINOS
691 const auto &local_dof_indices = copy.local_dof_indices_;
692 const auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
699 auto local_dof_indices = copy.local_dof_indices_;
700 auto neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
702 const auto &cell_mass_matrix = copy.cell_mass_matrix_;
703 const auto &cell_mass_matrix_inverse = copy.cell_mass_matrix_inverse_;
704 const auto &cell_cij_matrix = copy.cell_cij_matrix_;
705 const auto &interface_cij_matrix = copy.interface_cij_matrix_;
706 const auto &cell_measure = copy.cell_measure_;
708 if (!is_locally_owned)
711#ifndef DEAL_II_WITH_TRILINOS
713 for (
auto &indices : neighbor_local_dof_indices)
717 affine_constraints_assembly.distribute_local_to_global(
718 cell_mass_matrix, local_dof_indices, mass_matrix_tmp);
720 for (
int k = 0; k < dim; ++k) {
721 affine_constraints_assembly.distribute_local_to_global(
722 cell_cij_matrix[k], local_dof_indices, cij_matrix_tmp[k]);
724 for (
unsigned int f_index = 0; f_index < copy.n_faces; ++f_index) {
725 if (neighbor_local_dof_indices[f_index].size() != 0) {
726 affine_constraints_assembly.distribute_local_to_global(
727 interface_cij_matrix[f_index][k],
729 neighbor_local_dof_indices[f_index],
735 if (discretization_->have_discontinuous_ansatz())
736 affine_constraints_assembly.distribute_local_to_global(
737 cell_mass_matrix_inverse,
739 mass_matrix_inverse_tmp);
741 measure_of_omega_ += cell_measure;
744 WorkStream::run(dof_handler.begin_active(),
746 local_assemble_system,
747 copy_local_to_global,
748 AssemblyScratchData<dim>(*discretization_),
749#ifdef DEAL_II_WITH_TRILINOS
750 AssemblyCopyData<dim, double>());
752 AssemblyCopyData<dim, Number>());
755#ifdef DEAL_II_WITH_TRILINOS
757 for (
auto &it : cij_matrix_tmp)
760 mass_matrix_.read_in(mass_matrix_tmp,
false);
761 if (discretization_->have_discontinuous_ansatz())
762 mass_matrix_inverse_.read_in(mass_matrix_inverse_tmp,
false);
763 cij_matrix_.read_in(cij_matrix_tmp,
false);
765 mass_matrix_.read_in(mass_matrix_tmp,
true);
766 if (discretization_->have_discontinuous_ansatz())
767 mass_matrix_inverse_.read_in(mass_matrix_inverse_tmp,
true);
768 cij_matrix_.read_in(cij_matrix_tmp,
true);
771 mass_matrix_.update_ghost_rows();
772 if (discretization_->have_discontinuous_ansatz())
773 mass_matrix_inverse_.update_ghost_rows();
774 cij_matrix_.update_ghost_rows();
776 measure_of_omega_ = Utilities::MPI::sum(
777 measure_of_omega_, mpi_ensemble_.ensemble_communicator());
784#ifdef DEAL_II_WITH_TRILINOS
787 affine_constraints_assembly.set_zero(one);
789 ScalarVector local_lumped_mass_matrix(scalar_partitioner_);
790 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
793 for (
unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
795 lumped_mass_matrix_.local_element(i) =
796 local_lumped_mass_matrix.local_element(i);
797 lumped_mass_matrix_inverse_.local_element(i) =
798 1. / lumped_mass_matrix_.local_element(i);
800 lumped_mass_matrix_.update_ghost_values();
801 lumped_mass_matrix_inverse_.update_ghost_values();
805 dealii::Vector<Number> one(mass_matrix_tmp.m());
807 affine_constraints_assembly.set_zero(one);
809 dealii::Vector<Number> local_lumped_mass_matrix(mass_matrix_tmp.m());
810 mass_matrix_tmp.vmult(local_lumped_mass_matrix, one);
812 for (
unsigned int i = 0; i < scalar_partitioner_->locally_owned_size();
814 lumped_mass_matrix_.local_element(i) = local_lumped_mass_matrix(i);
815 lumped_mass_matrix_inverse_.local_element(i) =
816 1. / lumped_mass_matrix_.local_element(i);
818 lumped_mass_matrix_.update_ghost_values();
819 lumped_mass_matrix_inverse_.update_ghost_values();
827 if (discretization_->have_discontinuous_ansatz()) {
828#ifdef DEAL_II_WITH_TRILINOS
829 TrilinosWrappers::SparseMatrix incidence_matrix_tmp;
830 incidence_matrix_tmp.reinit(trilinos_sparsity_pattern);
832 dealii::SparseMatrix<Number> incidence_matrix_tmp;
833 incidence_matrix_tmp.reinit(sparsity_pattern_assembly);
837 const auto local_assemble_system = [&](
const auto &cell,
842 auto &is_locally_owned = copy.is_locally_owned_;
843 auto &local_dof_indices = copy.local_dof_indices_;
844 auto &neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
845 auto &interface_incidence_matrix = copy.interface_incidence_matrix_;
846 auto &hp_fe_face_values_nodal = scratch.hp_fe_face_values_nodal_;
847 auto &hp_fe_neighbor_face_values_nodal =
848 scratch.hp_fe_neighbor_face_values_nodal_;
850#ifdef DEAL_II_WITH_TRILINOS
851 is_locally_owned = cell->is_locally_owned();
853 is_locally_owned = !cell->is_artificial();
855 if (!is_locally_owned)
858 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
860 for (
auto &matrix : interface_incidence_matrix)
861 matrix.reinit(dofs_per_cell, dofs_per_cell);
863 local_dof_indices.resize(dofs_per_cell);
864 cell->get_dof_indices(local_dof_indices);
867 for (
auto &matrix : interface_incidence_matrix)
870 for (
const auto f_index : cell->face_indices()) {
871 const auto &face = cell->face(f_index);
874 const bool has_neighbor =
875 !face->at_boundary() || cell->has_periodic_neighbor(f_index);
879 neighbor_local_dof_indices[f_index].resize(0);
884 const auto neighbor_cell =
885 cell->neighbor_or_periodic_neighbor(f_index);
886 if (neighbor_cell->is_artificial()) {
889 neighbor_local_dof_indices[f_index].resize(0);
893 const unsigned int f_index_neighbor =
894 cell->has_periodic_neighbor(f_index)
895 ? cell->periodic_neighbor_of_periodic_neighbor(f_index)
896 : cell->neighbor_of_neighbor(f_index);
898 neighbor_local_dof_indices[f_index].resize(dofs_per_cell);
899 neighbor_cell->get_dof_indices(neighbor_local_dof_indices[f_index]);
901 hp_fe_face_values_nodal.reinit(cell, f_index);
902 const auto &fe_face_values_nodal =
903 hp_fe_face_values_nodal.get_present_fe_values();
904 hp_fe_neighbor_face_values_nodal.reinit(neighbor_cell,
906 const auto &fe_neighbor_face_values_nodal =
907 hp_fe_neighbor_face_values_nodal.get_present_fe_values();
911 for (
unsigned int q :
912 fe_face_values_nodal.quadrature_point_indices()) {
914 for (
unsigned int j : fe_neighbor_face_values_nodal.dof_indices()) {
915 const auto v_j = fe_neighbor_face_values_nodal.shape_value(j, q);
916 for (
unsigned int i : fe_face_values_nodal.dof_indices()) {
917 const auto v_i = fe_face_values_nodal.shape_value(i, q);
918 constexpr auto eps = std::numeric_limits<Number>::epsilon();
919 if (std::abs(v_i * v_j) > 100. * eps) {
920 const auto &ansatz = discretization_->ansatz();
922 const auto glob_i = local_dof_indices[i];
923 const auto glob_j = neighbor_local_dof_indices[f_index][j];
924 const auto m_i = lumped_mass_matrix_[glob_i];
925 const auto m_j = lumped_mass_matrix_[glob_j];
927 Number(0.5) * (m_i + m_j) / measure_of_omega_;
938 r_ij =
std::pow(hd_ij, incidence_relaxation_even_ / dim);
945 r_ij =
std::pow(hd_ij, incidence_relaxation_odd_ / dim);
948 interface_incidence_matrix[f_index](i, j) += r_ij;
956 const auto copy_local_to_global = [&](
const auto ©) {
957 const auto &is_locally_owned = copy.is_locally_owned_;
958#ifdef DEAL_II_WITH_TRILINOS
959 const auto &local_dof_indices = copy.local_dof_indices_;
960 const auto &neighbor_local_dof_indices =
961 copy.neighbor_local_dof_indices_;
968 auto local_dof_indices = copy.local_dof_indices_;
969 auto neighbor_local_dof_indices = copy.neighbor_local_dof_indices_;
971 const auto &interface_incidence_matrix =
972 copy.interface_incidence_matrix_;
974 if (!is_locally_owned)
977#ifndef DEAL_II_WITH_TRILINOS
979 for (
auto &indices : neighbor_local_dof_indices)
983 for (
unsigned int f_index = 0; f_index < copy.n_faces; ++f_index) {
984 if (neighbor_local_dof_indices[f_index].size() != 0) {
985 affine_constraints_assembly.distribute_local_to_global(
986 interface_incidence_matrix[f_index],
988 neighbor_local_dof_indices[f_index],
989 incidence_matrix_tmp);
994 WorkStream::run(dof_handler.begin_active(),
996 local_assemble_system,
997 copy_local_to_global,
998 AssemblyScratchData<dim>(*discretization_),
999#ifdef DEAL_II_WITH_TRILINOS
1000 AssemblyCopyData<dim, double>());
1002 AssemblyCopyData<dim, Number>());
1005#ifdef DEAL_II_WITH_TRILINOS
1006 incidence_matrix_.read_in(incidence_matrix_tmp,
false);
1008 incidence_matrix_.read_in(incidence_matrix_tmp,
true);
1010 incidence_matrix_.update_ghost_rows();
1017 boundary_map_ = construct_boundary_map(
1018 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
1020 coupling_boundary_pairs_ = collect_coupling_boundary_pairs(
1021 dof_handler.begin_active(), dof_handler.end(), *scalar_partitioner_);
1024#ifdef DEBUG_SYMMETRY_CHECK
1029 double total_mass = 0.;
1030 for (
unsigned int i = 0; i < n_locally_owned_; ++i)
1031 total_mass += lumped_mass_matrix_.local_element(i);
1033 Utilities::MPI::sum(total_mass, mpi_ensemble_.ensemble_communicator());
1035 Assert(std::abs(measure_of_omega_ - total_mass) <
1036 1.e-12 * measure_of_omega_,
1038 "Total mass differs from the measure of the domain."));
1044 for (
unsigned int i = 0; i < n_locally_owned_; ++i) {
1046 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1047 if (row_length == 1)
1051 mass_matrix_.get_entry(i, 0) - lumped_mass_matrix_.local_element(i);
1054 constexpr auto simd_length = VectorizedArray<Number>::size();
1055 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1056 for (
unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1057 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1059 Assert(j < n_locally_relevant_, dealii::ExcInternalError());
1061 const auto m_ij = mass_matrix_.get_entry(i, col_idx);
1062 if (discretization_->have_discontinuous_ansatz()) {
1065 Assert(std::abs(m_ij) > -1.e-12, dealii::ExcInternalError());
1067 Assert(std::abs(m_ij) > 1.e-12, dealii::ExcInternalError());
1071 const auto m_ji = mass_matrix_.get_transposed_entry(i, col_idx);
1072 if (std::abs(m_ij - m_ji) >= 1.e-12) {
1074 std::stringstream ss;
1075 ss <<
"m_ij matrix is not symmetric: " << m_ij <<
" <-> " << m_ji;
1076 Assert(
false, dealii::ExcMessage(ss.str()));
1080 Assert(std::abs(sum) < 1.e-12, dealii::ExcInternalError());
1087 for (
unsigned int i = 0; i < n_locally_owned_; ++i) {
1089 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1090 if (row_length == 1)
1093 auto sum = cij_matrix_.get_tensor(i, 0);
1096 constexpr auto simd_length = VectorizedArray<Number>::size();
1097 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1098 for (
unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1099 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1101 Assert(j < n_locally_relevant_, dealii::ExcInternalError());
1103 const auto c_ij = cij_matrix_.get_tensor(i, col_idx);
1104 Assert(c_ij.norm() > 1.e-12, dealii::ExcInternalError());
1107 const auto c_ji = cij_matrix_.get_transposed_tensor(i, col_idx);
1108 if ((c_ij + c_ji).norm() >= 1.e-12) {
1112 CouplingDescription coupling{i, col_idx, j};
1113 const auto it = std::find(coupling_boundary_pairs_.begin(),
1114 coupling_boundary_pairs_.end(),
1116 if (it == coupling_boundary_pairs_.end()) {
1117 std::stringstream ss;
1118 ss <<
"c_ij matrix is not anti-symmetric: " << c_ij <<
" <-> "
1120 Assert(
false, dealii::ExcMessage(ss.str()));
1125 Assert(sum.norm() < 1.e-12, dealii::ExcInternalError());
1131 template <
int dim,
typename Number>
1132 void OfflineData<dim, Number>::create_multigrid_data()
1135 std::cout <<
"OfflineData<dim, Number>::create_multigrid_data()"
1139 auto &dof_handler = *dof_handler_;
1141 Assert(!dof_handler.has_hp_capabilities(), dealii::ExcInternalError());
1143 dof_handler.distribute_mg_dofs();
1145 const auto n_levels = dof_handler.get_triangulation().n_global_levels();
1147 AffineConstraints<float> level_constraints;
1150 level_boundary_map_.resize(n_levels);
1151 level_lumped_mass_matrix_.resize(n_levels);
1153 for (
unsigned int level = 0; level < n_levels; ++level) {
1156 IndexSet relevant_dofs;
1157 dealii::DoFTools::extract_locally_relevant_level_dofs(
1158 dof_handler, level, relevant_dofs);
1159 const auto partitioner = std::make_shared<Utilities::MPI::Partitioner>(
1160 dof_handler.locally_owned_mg_dofs(level),
1162 mpi_ensemble_.ensemble_communicator());
1163 level_lumped_mass_matrix_[level].reinit(partitioner);
1164 std::vector<types::global_dof_index> dof_indices(
1165 dof_handler.get_fe().dofs_per_cell);
1166 dealii::Vector<Number> mass_values(dof_handler.get_fe().dofs_per_cell);
1167 dealii::hp::FEValues<dim> hp_fe_values(discretization_->mapping(),
1168 discretization_->finite_element(),
1169 discretization_->quadrature(),
1170 update_values | update_JxW_values);
1171 for (
const auto &cell : dof_handler.cell_iterators_on_level(level))
1174 if (cell->is_locally_owned_on_level()) {
1175 hp_fe_values.reinit(cell);
1176 const auto &fe_values = hp_fe_values.get_present_fe_values();
1177 for (
unsigned int i = 0; i < mass_values.size(); ++i) {
1179 for (
unsigned int q = 0; q < fe_values.n_quadrature_points; ++q)
1180 sum += fe_values.shape_value(i, q) * fe_values.JxW(q);
1181 mass_values(i) = sum;
1183 cell->get_mg_dof_indices(dof_indices);
1184 level_constraints.distribute_local_to_global(
1185 mass_values, dof_indices, level_lumped_mass_matrix_[level]);
1191 level_boundary_map_[level] = construct_boundary_map(
1192 dof_handler.begin_mg(level), dof_handler.end_mg(level), *partitioner);
1197 template <
int dim,
typename Number>
1198 template <
typename ITERATOR1,
typename ITERATOR2>
1200 const ITERATOR1 &begin,
1201 const ITERATOR2 &end,
1202 const Utilities::MPI::Partitioner &partitioner)
const -> BoundaryMap
1205 std::cout <<
"OfflineData<dim, Number>::construct_boundary_map()"
1213 using BoundaryData = std::tuple<dealii::Tensor<1, dim, Number> ,
1216 dealii::types::boundary_id ,
1217 dealii::Point<dim>> ;
1218 std::multimap<unsigned int, BoundaryData> preliminary_map;
1220 std::vector<dealii::types::global_dof_index> local_dof_indices;
1222 dealii::hp::FEFaceValues<dim> hp_fe_face_values(
1223 discretization_->mapping(),
1224 discretization_->finite_element(),
1225 discretization_->face_quadrature(),
1226 dealii::update_normal_vectors | dealii::update_values |
1227 dealii::update_JxW_values);
1229 for (
auto cell = begin; cell != end; ++cell) {
1237 if ((cell->is_active() && cell->is_artificial()) ||
1238 (!cell->is_active() && cell->is_artificial_on_level()))
1241 const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
1242 const auto &support_points = cell->get_fe().get_unit_support_points();
1244 local_dof_indices.resize(dofs_per_cell);
1245 cell->get_active_or_mg_dof_indices(local_dof_indices);
1247 for (
auto f : cell->face_indices()) {
1248 const auto face = cell->face(f);
1249 const auto id = face->boundary_id();
1253 if (!face->at_boundary())
1264 hp_fe_face_values.reinit(cell, f);
1265 const auto &fe_face_values = hp_fe_face_values.get_present_fe_values();
1266 const auto &mapping =
1267 hp_fe_face_values.get_mapping_collection()[cell->active_fe_index()];
1269 for (
unsigned int j : fe_face_values.dof_indices()) {
1270 if (!cell->get_fe().has_support_on_face(j, f))
1273 Number boundary_mass = 0.;
1274 dealii::Tensor<1, dim, Number> normal;
1276 for (
unsigned int q : fe_face_values.quadrature_point_indices()) {
1277 const auto JxW = fe_face_values.JxW(q);
1278 const auto phi_i = fe_face_values.shape_value(j, q);
1280 boundary_mass += phi_i * JxW;
1281 normal += phi_i * fe_face_values.normal_vector(q) * JxW;
1284 const auto global_index = local_dof_indices[j];
1285 const auto index = partitioner.global_to_local(global_index);
1288 if (index >= n_locally_owned_)
1292 const unsigned int row_length =
1293 sparsity_pattern_simd_.row_length(index);
1294 if (row_length == 1)
1297 Point<dim> position =
1298 mapping.transform_unit_to_real_cell(cell, support_points[j]);
1304 preliminary_map.insert(
1305 {index, {normal, boundary_mass, boundary_mass, id, position}});
1321 std::multimap<unsigned int, BoundaryData> filtered_map;
1322 std::set<dealii::types::global_dof_index> boundary_dofs;
1323 for (
auto entry : preliminary_map) {
1324 bool inserted =
false;
1325 const auto range = filtered_map.equal_range(entry.first);
1326 for (
auto it = range.first; it != range.second; ++it) {
1331 new_point] = entry.second;
1332 auto &[normal, normal_mass, boundary_mass, id, point] = it->second;
1337 Assert(point.distance(new_point) < 1.0e-14, dealii::ExcInternalError());
1339 if (normal * new_normal / normal.norm() / new_normal.norm() > 0.50) {
1344 normal += new_normal;
1345 boundary_mass += new_boundary_mass;
1349 }
else if constexpr (dim == 2) {
1365 filtered_map.insert(entry);
1372 BoundaryMap boundary_map;
1374 std::begin(filtered_map),
1375 std::end(filtered_map),
1376 std::back_inserter(boundary_map),
1378 auto index = it.first;
1379 const auto &[normal, normal_mass, boundary_mass, id, point] =
1382 const auto new_normal_mass =
1383 normal.norm() + std::numeric_limits<Number>::epsilon();
1384 const auto new_normal = normal / new_normal_mass;
1386 return {index, new_normal, new_normal_mass, boundary_mass, id, point};
1389 return boundary_map;
1393 template <
int dim,
typename Number>
1394 template <
typename ITERATOR1,
typename ITERATOR2>
1396 const ITERATOR1 &begin,
1397 const ITERATOR2 &end,
1398 const Utilities::MPI::Partitioner &partitioner)
const
1399 -> CouplingBoundaryPairs
1402 std::cout <<
"OfflineData<dim, Number>::collect_coupling_boundary_pairs()"
1412 std::set<unsigned int> locally_relevant_boundary_indices;
1414 std::vector<dealii::types::global_dof_index> local_dof_indices;
1416 for (
auto cell = begin; cell != end; ++cell) {
1419 if (cell->is_artificial())
1422 const auto &finite_element = cell->get_fe();
1423 const unsigned int dofs_per_cell = finite_element.dofs_per_cell;
1424 local_dof_indices.resize(dofs_per_cell);
1425 cell->get_active_or_mg_dof_indices(local_dof_indices);
1427 for (
auto f : cell->face_indices()) {
1428 const auto face = cell->face(f);
1429 const auto id = face->boundary_id();
1431 if (!face->at_boundary())
1438 for (
unsigned int j = 0; j < dofs_per_cell; ++j) {
1440 if (!cell->get_fe().has_support_on_face(j, f))
1443 const auto global_index = local_dof_indices[j];
1444 const auto index = partitioner.global_to_local(global_index);
1447 if (index >= n_locally_relevant_)
1450 locally_relevant_boundary_indices.insert(index);
1459 CouplingBoundaryPairs result;
1461 for (
const auto i : locally_relevant_boundary_indices) {
1464 if (i >= n_locally_owned_)
1467 const unsigned int row_length = sparsity_pattern_simd_.row_length(i);
1470 if (row_length == 1)
1473 const unsigned int *js = sparsity_pattern_simd_.columns(i);
1474 constexpr auto simd_length = VectorizedArray<Number>::size();
1476 for (
unsigned int col_idx = 1; col_idx < row_length; ++col_idx) {
1477 const auto j = *(i < n_locally_internal_ ? js + col_idx * simd_length
1480 if (locally_relevant_boundary_indices.count(j) != 0) {
1481 result.push_back({i, col_idx, j});
std::tuple< unsigned int, dealii::Tensor< 1, dim, Number >, Number, Number, dealii::types::boundary_id, dealii::Point< dim > > BoundaryDescription
OfflineData(const MPIEnsemble &mpi_ensemble, 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)
void make_extended_sparsity_pattern_dg(const dealii::DoFHandler< dim > &dof_handler, SPARSITY &dsp, const dealii::AffineConstraints< Number > &affine_constraints, bool keep_constrained)
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)
T pow(const T x, const T b)
dealii::LinearAlgebra::distributed::Vector< Number > ScalarVector
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)