From 979d00057703ba0ed8d4f0c55edc6ca148cc13cb Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 10 Jun 2022 12:18:00 +0200 Subject: [PATCH] MatrixFree: Element-centric loops with different cell geometry compressions --- include/deal.II/matrix_free/fe_evaluation.h | 5 +- include/deal.II/matrix_free/mapping_info.h | 21 ++-- .../matrix_free/mapping_info.templates.h | 100 +++++++++++++----- tests/matrix_free/compare_faces_by_cells.cc | 8 +- 4 files changed, 93 insertions(+), 41 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 7987649224..84c26effb4 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -8189,8 +8189,9 @@ FEFaceEvaluationmatrix_free != nullptr, ExcNotInitialized()); - this->cell_type = this->matrix_free->get_mapping_info().cell_type[cell_index]; - this->cell = cell_index; + this->cell_type = this->matrix_free->get_mapping_info() + .faces_by_cells_type[cell_index][face_number]; + this->cell = cell_index; this->subface_index = GeometryInfo::max_children_per_cell; this->dof_access_index = internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h index 7087e5fda3..bd46a1ce5f 100644 --- a/include/deal.II/matrix_free/mapping_info.h +++ b/include/deal.II/matrix_free/mapping_info.h @@ -67,7 +67,7 @@ namespace internal initialize( const dealii::Triangulation & tria, const std::vector> &cells, - const FaceInfo & faces, + const FaceInfo & face_info, const std::vector &active_fe_index, const std::shared_ptr> &mapping, const std::vector> & quad, @@ -86,7 +86,7 @@ namespace internal update_mapping( const dealii::Triangulation & tria, const std::vector> &cells, - const FaceInfo & faces, + const FaceInfo & face_info, const std::vector &active_fe_index, const std::shared_ptr> &mapping); @@ -157,6 +157,15 @@ namespace internal */ std::vector face_type; + /** + * Store whether the face geometry for the face_data_by_cells data + * structure is represented as Cartesian (cell type 0), with constant + * transform data (Jacobians) (cell type 1), or a general type (cell + * type 3). Note that both the interior and exterior agree on the type + * of the data structure, using the more general of the two. + */ + std::vector> faces_by_cells_type; + /** * The data cache for the cells. */ @@ -209,15 +218,14 @@ namespace internal * given as a tuple of the level and index within the level as used in * the main initialization of the class * - * @param faces The description of the connectivity from faces to cells - * as filled in the MatrixFree class + * @param face_info The description of the connectivity from faces to + * cells as filled in the MatrixFree class */ void compute_mapping_q( const dealii::Triangulation & tria, const std::vector> &cells, - const std::vector> - &faces); + const FaceInfo & face_info); /** * Computes the information in the given cells, called within @@ -251,6 +259,7 @@ namespace internal initialize_faces_by_cells( const dealii::Triangulation & tria, const std::vector> &cells, + const FaceInfo & face_info, const dealii::hp::MappingCollection & mapping); }; diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 5aa4fb3409..e9178bb5e2 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -181,7 +181,7 @@ namespace internal // use the fast method. if (active_fe_index.empty() && !cells.empty() && mapping->size() == 1 && dynamic_cast *>(&mapping->operator[](0))) - compute_mapping_q(tria, cells, face_info.faces); + compute_mapping_q(tria, cells, face_info); else { // Could call these functions in parallel, but not useful because @@ -189,7 +189,7 @@ namespace internal initialize_cells(tria, cells, active_fe_index, *mapping); initialize_faces( tria, cells, face_info.faces, active_fe_index, *mapping); - initialize_faces_by_cells(tria, cells, *mapping); + initialize_faces_by_cells(tria, cells, face_info, *mapping); } } @@ -219,7 +219,7 @@ namespace internal if (active_fe_index.empty() && !cells.empty() && mapping->size() == 1 && dynamic_cast *>(&mapping->operator[](0))) - compute_mapping_q(tria, cells, face_info.faces); + compute_mapping_q(tria, cells, face_info); else { // Could call these functions in parallel, but not useful because @@ -227,7 +227,7 @@ namespace internal initialize_cells(tria, cells, active_fe_index, *mapping); initialize_faces( tria, cells, face_info.faces, active_fe_index, *mapping); - initialize_faces_by_cells(tria, cells, *mapping); + initialize_faces_by_cells(tria, cells, face_info, *mapping); } } @@ -2530,7 +2530,7 @@ namespace internal MappingInfo::compute_mapping_q( const dealii::Triangulation & tria, const std::vector> &cell_array, - const std::vector> &faces) + const FaceInfo & face_info) { // step 1: extract quadrature point data with the data appropriate for // MappingQ @@ -2714,6 +2714,8 @@ namespace internal std::size_t(2U))); } + const std::vector> + &faces = face_info.faces; if (faces.empty()) return; @@ -2882,7 +2884,10 @@ namespace internal // transitioned to extracting the information from cell quadrature // points but we need to figure out the correct indices of neighbors // within the list of arrays still - initialize_faces_by_cells(tria, cell_array, *this->mapping_collection); + initialize_faces_by_cells(tria, + cell_array, + face_info, + *this->mapping_collection); } @@ -2892,6 +2897,7 @@ namespace internal MappingInfo::initialize_faces_by_cells( const dealii::Triangulation & tria, const std::vector> &cells, + const FaceInfo & face_info, const dealii::hp::MappingCollection & mapping_in) { if (update_flags_faces_by_cells == update_default) @@ -2901,14 +2907,49 @@ namespace internal AssertDimension(mapping_in.size(), 1); const auto &mapping = mapping_in[0]; - const unsigned int n_quads = face_data_by_cells.size(); - const unsigned int n_lanes = VectorizedArrayType::size(); - UpdateFlags update_flags = + const unsigned int n_quads = face_data_by_cells.size(); + constexpr unsigned int n_lanes = VectorizedArrayType::size(); + UpdateFlags update_flags = (update_flags_faces_by_cells & update_quadrature_points ? update_quadrature_points : update_default) | update_normal_vectors | update_JxW_values | update_jacobians; + const auto compute_neighbor_index = + [&face_info](const unsigned int cell, + const unsigned int face, + const unsigned int lane) { + constexpr unsigned int n_lanes = VectorizedArrayType::size(); + const unsigned int cell_this = cell * n_lanes + lane; + unsigned int face_index = + face_info.cell_and_face_to_plain_faces(cell, face, lane); + if (face_index == numbers::invalid_unsigned_int) + return numbers::invalid_unsigned_int; + + unsigned int cell_neighbor = face_info.faces[face_index / n_lanes] + .cells_interior[face_index % n_lanes]; + if (cell_neighbor == cell_this) + cell_neighbor = face_info.faces[face_index / n_lanes] + .cells_exterior[face_index % n_lanes]; + return cell_neighbor; + }; + + faces_by_cells_type.resize(cell_type.size()); + for (unsigned int cell = 0; cell < cell_type.size(); ++cell) + for (const unsigned int face : GeometryInfo::face_indices()) + { + faces_by_cells_type[cell][face] = cell_type[cell]; + for (unsigned int v = 0; v < n_lanes; ++v) + { + const unsigned int cell_neighbor = + compute_neighbor_index(cell, face, v); + if (cell_neighbor != numbers::invalid_unsigned_int) + faces_by_cells_type[cell][face] = + std::max(faces_by_cells_type[cell][face], + cell_type[cell_neighbor / n_lanes]); + } + } + for (unsigned int my_q = 0; my_q < n_quads; ++my_q) { // since we already know the cell type, we can pre-allocate the right @@ -2924,7 +2965,7 @@ namespace internal for (unsigned int i = 0; i < cell_type.size(); ++i) for (const unsigned int face : GeometryInfo::face_indices()) { - if (cell_type[i] <= affine) + if (faces_by_cells_type[i][face] <= affine) { face_data_by_cells[my_q].data_index_offsets [i * GeometryInfo::faces_per_cell + face] = @@ -3010,6 +3051,8 @@ namespace internal .data_index_offsets[cell * GeometryInfo::faces_per_cell + face]; + const GeometryType my_cell_type = faces_by_cells_type[cell][face]; + for (unsigned int v = 0; v < n_lanes; ++v) { typename dealii::Triangulation::cell_iterator cell_it( @@ -3018,18 +3061,15 @@ namespace internal cells[cell * n_lanes + v].second); fe_val.reinit(cell_it, face); - const bool is_local = - (cell_it->is_active() ? - cell_it->is_locally_owned() : - cell_it->is_locally_owned_on_level()) && - (!cell_it->at_boundary(face) || - (cell_it->at_boundary(face) && - cell_it->has_periodic_neighbor(face))); + const unsigned int cell_neighbor = + compute_neighbor_index(cell, face, v); - if (is_local) + if (cell_neighbor != numbers::invalid_unsigned_int) { - auto cell_it_neigh = - cell_it->neighbor_or_periodic_neighbor(face); + typename dealii::Triangulation::cell_iterator + cell_it_neigh(&tria, + cells[cell_neighbor].first, + cells[cell_neighbor].second); fe_val_neigh.reinit(cell_it_neigh, cell_it->at_boundary(face) ? cell_it->periodic_neighbor_face_no( @@ -3038,7 +3078,7 @@ namespace internal } // copy data for affine data type - if (cell_type[cell] <= affine) + if (my_cell_type <= affine) { if (update_flags & update_JxW_values) face_data_by_cells[my_q].JxW_values[offset][v] = @@ -3057,7 +3097,8 @@ namespace internal inv_jac[d][ee]; } } - if (is_local && (update_flags & update_jacobians)) + if (cell_neighbor != numbers::invalid_unsigned_int && + (update_flags & update_jacobians)) for (unsigned int q = 0; q < fe_val.n_quadrature_points; ++q) { @@ -3109,7 +3150,8 @@ namespace internal inv_jac[d][ee]; } } - if (is_local && (update_flags & update_jacobians)) + if (cell_neighbor != numbers::invalid_unsigned_int && + (update_flags & update_jacobians)) for (unsigned int q = 0; q < fe_val.n_quadrature_points; ++q) { @@ -3149,9 +3191,9 @@ namespace internal } if (update_flags & update_normal_vectors && update_flags & update_jacobians) - for (unsigned int q = 0; q < (cell_type[cell] <= affine ? - 1 : - fe_val.n_quadrature_points); + for (unsigned int q = 0; + q < + (my_cell_type <= affine ? 1 : fe_val.n_quadrature_points); ++q) face_data_by_cells[my_q] .normals_times_jacobians[0][offset + q] = @@ -3159,9 +3201,9 @@ namespace internal face_data_by_cells[my_q].jacobians[0][offset + q]; if (update_flags & update_normal_vectors && update_flags & update_jacobians) - for (unsigned int q = 0; q < (cell_type[cell] <= affine ? - 1 : - fe_val.n_quadrature_points); + for (unsigned int q = 0; + q < + (my_cell_type <= affine ? 1 : fe_val.n_quadrature_points); ++q) face_data_by_cells[my_q] .normals_times_jacobians[1][offset + q] = diff --git a/tests/matrix_free/compare_faces_by_cells.cc b/tests/matrix_free/compare_faces_by_cells.cc index a30a4d3397..8e23091004 100644 --- a/tests/matrix_free/compare_faces_by_cells.cc +++ b/tests/matrix_free/compare_faces_by_cells.cc @@ -15,7 +15,7 @@ -// compares the computation of the diagonal using a face integration +// compares the computation of the diagonal using the face integration // facilities with the alternative reinit(cell_index, face_number) approach #include @@ -82,8 +82,8 @@ public: compute_diagonal_by_face( LinearAlgebra::distributed::Vector &result) const { - int dummy; - result = 0; + int dummy = 0; + result = 0; data.loop(&LaplaceOperator::local_diagonal_dummy, &LaplaceOperator::local_diagonal_face, &LaplaceOperator::local_diagonal_boundary, @@ -96,7 +96,7 @@ public: compute_diagonal_by_cell( LinearAlgebra::distributed::Vector &result) const { - int dummy; + int dummy = 0; result.zero_out_ghost_values(); data.cell_loop(&LaplaceOperator::local_diagonal_by_cell, this, -- 2.39.5