From: Peter Munch Date: Thu, 3 Dec 2020 10:21:52 +0000 (+0100) Subject: MatrixFree: precompute subranges X-Git-Tag: v9.3.0-rc1~390^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9e9d7b4500489811576d0e992a2a77ac5db8646f;p=dealii.git MatrixFree: precompute subranges --- diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 6e407a2053..1c220afe3f 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -1472,32 +1472,7 @@ public: const unsigned int dof_handler_index = 0) const; /** - * In the hp-adaptive case, a subrange of internal faces as computed during - * loop() might contain internal faces with elements of different active - * FE indices. Use this function to compute what the subrange for a given pair - * of active FE indices is. - */ - std::pair - create_inner_face_subrange_hp_by_index( - const std::pair &range, - const unsigned int fe_index_interior, - const unsigned int fe_index_exterior, - const unsigned int dof_handler_index = 0) const; - - /** - * In the hp-adaptive case, a subrange of boundary faces as computed during - * loop() might contain boundary faces with elements of different active - * FE indices. Use this function to compute what the subrange for a given - * active FE indices is. - */ - std::pair - create_boundary_face_subrange_hp_by_index( - const std::pair &range, - const unsigned int fe_index, - const unsigned int dof_handler_index = 0) const; - - /** - * In the hp-adaptive case, return number of active_fe_indices. + * In the hp adaptive case, return number of active_fe_indices. */ unsigned int n_active_fe_indices() const; @@ -4642,49 +4617,51 @@ namespace internal } } - // Runs the assembler on interior faces. If no function is given, nothing - // is done virtual void - face(const std::pair &face_range) override + cell(const unsigned int range_index) override { - if (face_function != nullptr && face_range.second > face_range.first) - for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i) - for (unsigned int j = 0; j < matrix_free.n_active_fe_indices(); ++j) - { - const auto face_subrange = - matrix_free.create_inner_face_subrange_hp_by_index(face_range, - i, - j); - - if (face_subrange.second <= face_subrange.first) - continue; - (container.* - face_function)(matrix_free, this->dst, this->src, face_subrange); - } + process_range(cell_function, + matrix_free.get_task_info().cell_partition_data_hp_ptr, + matrix_free.get_task_info().cell_partition_data_hp, + range_index); } - // Runs the assembler on boundary faces. If no function is given, nothing - // is done virtual void - boundary(const std::pair &face_range) override + face(const unsigned int range_index) override { - if (boundary_function != nullptr && face_range.second > face_range.first) - for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i) - { - const auto face_subrange = - matrix_free.create_boundary_face_subrange_hp_by_index(face_range, - i); + process_range(face_function, + matrix_free.get_task_info().face_partition_data_hp_ptr, + matrix_free.get_task_info().face_partition_data_hp, + range_index); + } - if (face_subrange.second <= face_subrange.first) - continue; + virtual void + boundary(const unsigned int range_index) override + { + process_range(boundary_function, + matrix_free.get_task_info().boundary_partition_data_hp_ptr, + matrix_free.get_task_info().boundary_partition_data_hp, + range_index); + } - (container.*boundary_function)(matrix_free, - this->dst, - this->src, - face_subrange); - } + private: + void + process_range(const function_type & fu, + const std::vector &ptr, + const std::vector &data, + const unsigned int range_index) + { + if (fu == nullptr) + return; + + for (unsigned int i = ptr[range_index]; i < ptr[range_index + 1]; ++i) + (container.*fu)(matrix_free, + this->dst, + this->src, + std::make_pair(data[2 * i], data[2 * i + 1])); } + public: // Starts the communication for the update ghost values operation. We // cannot call this update if ghost and destination are the same because // that would introduce spurious entries in the destination (there is also diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index df3f382117..3a4efd3b8b 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -174,92 +174,6 @@ namespace -template -std::pair -MatrixFree:: - create_inner_face_subrange_hp_by_index( - const std::pair &range, - const unsigned int fe_index_interior, - const unsigned int fe_index_exterior, - const unsigned int dof_handler_index) const -{ - if (dof_info[dof_handler_index].max_fe_index == 0 || - dof_handlers[dof_handler_index]->get_fe_collection().size() == 1) - return range; - - AssertIndexRange(fe_index_interior, dof_info[dof_handler_index].max_fe_index); - AssertIndexRange(fe_index_exterior, dof_info[dof_handler_index].max_fe_index); - const std::vector &fe_indices = - dof_info[dof_handler_index].cell_active_fe_index; - if (fe_indices.empty() == true) - return range; - else - { - std::pair return_range; - return_range.first = - std::lower_bound(face_info.faces.begin() + range.first, - face_info.faces.begin() + range.second, - std::pair{ - fe_index_interior, fe_index_exterior}, - FaceRangeCompartor(fe_indices, false)) - - face_info.faces.begin(); - return_range.second = - std::lower_bound(face_info.faces.begin() + return_range.first, - face_info.faces.begin() + range.second, - std::pair{ - fe_index_interior, fe_index_exterior}, - FaceRangeCompartor(fe_indices, true)) - - face_info.faces.begin(); - Assert(return_range.first >= range.first && - return_range.second <= range.second, - ExcInternalError()); - return return_range; - } -} - - - -template -std::pair -MatrixFree:: - create_boundary_face_subrange_hp_by_index( - const std::pair &range, - const unsigned int fe_index, - const unsigned int dof_handler_index) const -{ - if (dof_info[dof_handler_index].max_fe_index == 0 || - dof_handlers[dof_handler_index]->get_fe_collection().size() == 1) - return range; - - AssertIndexRange(fe_index, dof_info[dof_handler_index].max_fe_index); - const std::vector &fe_indices = - dof_info[dof_handler_index].cell_active_fe_index; - if (fe_indices.empty() == true) - return range; - else - { - std::pair return_range; - return_range.first = - std::lower_bound(face_info.faces.begin() + range.first, - face_info.faces.begin() + range.second, - fe_index, - FaceRangeCompartor(fe_indices, false)) - - face_info.faces.begin(); - return_range.second = - std::lower_bound(face_info.faces.begin() + return_range.first, - face_info.faces.begin() + range.second, - fe_index, - FaceRangeCompartor(fe_indices, true)) - - face_info.faces.begin(); - Assert(return_range.first >= range.first && - return_range.second <= range.second, - ExcInternalError()); - return return_range; - } -} - - - template void MatrixFree::renumber_dofs( @@ -572,6 +486,212 @@ MatrixFree::internal_reinit( } } + + // subdivide cell, face and boundary face partitioner data, s.t., all + // ranges have the same active_fe_indices + if (task_info.scheme != internal::MatrixFreeFunctions::TaskInfo:: + TasksParallelScheme::partition_color) + { + // ... for cells + { + auto &ptr = task_info.cell_partition_data_hp_ptr; + auto &data = task_info.cell_partition_data_hp; + + ptr = {0}; + data = {}; + data.reserve(2 * task_info.cell_partition_data.size()); + + for (unsigned int i = 0; i < task_info.cell_partition_data.size() - 1; + ++i) + { + if (task_info.cell_partition_data[i + 1] > + task_info.cell_partition_data[i]) + { + std::pair range{ + task_info.cell_partition_data[i], + task_info.cell_partition_data[i + 1]}; + + if (range.second > range.first) + for (unsigned int i = 0; i < this->n_active_fe_indices(); ++i) + { + const auto cell_subrange = + this->create_cell_subrange_hp_by_index(range, i); + + if (cell_subrange.second <= cell_subrange.first) + continue; + + data.push_back(cell_subrange.first); + data.push_back(cell_subrange.second); + } + } + + ptr.push_back(data.size() / 2); + } + } + + // ... for inner faces + if (task_info.face_partition_data.empty() == false) + { + auto &ptr = task_info.face_partition_data_hp_ptr; + auto &data = task_info.face_partition_data_hp; + + ptr = {0}; + data = {}; + data.reserve(2 * task_info.face_partition_data.size()); + + const auto create_inner_face_subrange_hp_by_index = + [&](const std::pair &range, + const unsigned int fe_index_interior, + const unsigned int fe_index_exterior, + const unsigned int dof_handler_index = + 0) -> std::pair { + if (dof_info[dof_handler_index].max_fe_index == 0 || + dof_handlers[dof_handler_index]->get_fe_collection().size() == + 1) + return range; + + AssertIndexRange(fe_index_interior, + dof_info[dof_handler_index].max_fe_index); + AssertIndexRange(fe_index_exterior, + dof_info[dof_handler_index].max_fe_index); + const std::vector &fe_indices = + dof_info[dof_handler_index].cell_active_fe_index; + if (fe_indices.empty() == true) + return range; + else + { + std::pair return_range; + return_range.first = + std::lower_bound(face_info.faces.begin() + range.first, + face_info.faces.begin() + range.second, + std::pair{ + fe_index_interior, fe_index_exterior}, + FaceRangeCompartor(fe_indices, false)) - + face_info.faces.begin(); + return_range.second = + std::lower_bound(face_info.faces.begin() + return_range.first, + face_info.faces.begin() + range.second, + std::pair{ + fe_index_interior, fe_index_exterior}, + FaceRangeCompartor(fe_indices, true)) - + face_info.faces.begin(); + Assert(return_range.first >= range.first && + return_range.second <= range.second, + ExcInternalError()); + return return_range; + } + }; + + for (unsigned int i = 0; i < task_info.face_partition_data.size() - 1; + ++i) + { + if (task_info.face_partition_data[i + 1] > + task_info.face_partition_data[i]) + { + std::pair range{ + task_info.face_partition_data[i], + task_info.face_partition_data[i + 1]}; + + if (range.second > range.first) + for (unsigned int i = 0; i < this->n_active_fe_indices(); + ++i) + for (unsigned int j = 0; j < this->n_active_fe_indices(); + ++j) + { + const auto subrange = + create_inner_face_subrange_hp_by_index(range, i, j); + + if (subrange.second <= subrange.first) + continue; + + data.push_back(subrange.first); + data.push_back(subrange.second); + } + } + + ptr.push_back(data.size() / 2); + } + } + + // ... for boundary faces + if (task_info.boundary_partition_data.empty() == false) + { + auto &ptr = task_info.boundary_partition_data_hp_ptr; + auto &data = task_info.boundary_partition_data_hp; + + ptr = {0}; + data = {}; + data.reserve(2 * task_info.boundary_partition_data.size()); + + const auto create_boundary_face_subrange_hp_by_index = + [&](const std::pair &range, + const unsigned int fe_index, + const unsigned int dof_handler_index = + 0) -> std::pair { + if (dof_info[dof_handler_index].max_fe_index == 0 || + dof_handlers[dof_handler_index]->get_fe_collection().size() == + 1) + return range; + + AssertIndexRange(fe_index, + dof_info[dof_handler_index].max_fe_index); + const std::vector &fe_indices = + dof_info[dof_handler_index].cell_active_fe_index; + if (fe_indices.empty() == true) + return range; + else + { + std::pair return_range; + return_range.first = + std::lower_bound(face_info.faces.begin() + range.first, + face_info.faces.begin() + range.second, + fe_index, + FaceRangeCompartor(fe_indices, false)) - + face_info.faces.begin(); + return_range.second = + std::lower_bound(face_info.faces.begin() + return_range.first, + face_info.faces.begin() + range.second, + fe_index, + FaceRangeCompartor(fe_indices, true)) - + face_info.faces.begin(); + Assert(return_range.first >= range.first && + return_range.second <= range.second, + ExcInternalError()); + return return_range; + } + }; + + for (unsigned int i = 0; + i < task_info.boundary_partition_data.size() - 1; + ++i) + { + if (task_info.boundary_partition_data[i + 1] > + task_info.boundary_partition_data[i]) + { + std::pair range{ + task_info.boundary_partition_data[i], + task_info.boundary_partition_data[i + 1]}; + + if (range.second > range.first) + for (unsigned int i = 0; i < this->n_active_fe_indices(); + ++i) + { + const auto cell_subrange = + create_boundary_face_subrange_hp_by_index(range, i); + + if (cell_subrange.second <= cell_subrange.first) + continue; + + data.push_back(cell_subrange.first); + data.push_back(cell_subrange.second); + } + } + + ptr.push_back(data.size() / 2); + } + } + } + // Evaluates transformations from unit to real cell, Jacobian determinants, // quadrature points in real space, based on the ordering of the cells // determined in @p extract_local_to_global_indices. diff --git a/include/deal.II/matrix_free/task_info.h b/include/deal.II/matrix_free/task_info.h index d985c3a570..71899a2f33 100644 --- a/include/deal.II/matrix_free/task_info.h +++ b/include/deal.II/matrix_free/task_info.h @@ -78,15 +78,20 @@ namespace internal virtual void cell(const std::pair &cell_range) = 0; + /// Runs the cell work specified by MatrixFree::loop or + /// MatrixFree::cell_loop + virtual void + cell(const unsigned int range_index) = 0; + /// Runs the body of the work on interior faces specified by /// MatrixFree::loop virtual void - face(const std::pair &face_range) = 0; + face(const unsigned int range_index) = 0; /// Runs the body of the work on boundary faces specified by /// MatrixFree::loop virtual void - boundary(const std::pair &face_range) = 0; + boundary(const unsigned int range_index) = 0; }; @@ -468,6 +473,19 @@ namespace internal */ std::vector cell_partition_data; + /** + * Like cell_partition_data but with precomputed subranges for each + * active fe index. The start and end point of a partition is given + * by cell_partition_data_hp_ptr. + */ + std::vector cell_partition_data_hp; + + /** + * Pointers within cell_partition_data_hp, indicating the start and end + * of a partition. + */ + std::vector cell_partition_data_hp_ptr; + /** * This is a linear storage of all partitions of inner faces, building a * range of indices of the form face_partition_data[idx] to @@ -477,6 +495,19 @@ namespace internal */ std::vector face_partition_data; + /** + * Like face_partition_data but with precomputed subranges for each + * active fe index pair. The start and end point of a partition is given + * by face_partition_data_hp_ptr. + */ + std::vector face_partition_data_hp; + + /** + * Pointers within face_partition_data_hp, indicating the start and end + * of a partition. + */ + std::vector face_partition_data_hp_ptr; + /** * This is a linear storage of all partitions of boundary faces, * building a range of indices of the form boundary_partition_data[idx] @@ -486,6 +517,19 @@ namespace internal */ std::vector boundary_partition_data; + /** + * Like boundary_partition_data but with precomputed subranges for each + * active fe index. The start and end point of a partition is given + * by boundary_partition_data_hp_ptr. + */ + std::vector boundary_partition_data_hp; + + /** + * Pointers within boundary_partition_data_hp, indicating the start and + * end of a partition. + */ + std::vector boundary_partition_data_hp_ptr; + /** * This is a linear storage of all partitions of interior faces on * boundaries to other processors that are not locally used, building a diff --git a/source/matrix_free/task_info.cc b/source/matrix_free/task_info.cc index 9ab41f175b..881a444ff1 100644 --- a/source/matrix_free/task_info.cc +++ b/source/matrix_free/task_info.cc @@ -77,19 +77,12 @@ namespace internal MFWorkerInterface *used_worker = worker != nullptr ? worker : *worker_pointer; Assert(used_worker != nullptr, ExcInternalError()); - used_worker->cell( - std::make_pair(task_info.cell_partition_data[partition], - task_info.cell_partition_data[partition + 1])); + used_worker->cell(partition); if (task_info.face_partition_data.empty() == false) { - used_worker->face( - std::make_pair(task_info.face_partition_data[partition], - task_info.face_partition_data[partition + 1])); - - used_worker->boundary(std::make_pair( - task_info.boundary_partition_data[partition], - task_info.boundary_partition_data[partition + 1])); + used_worker->face(partition); + used_worker->boundary(partition); } } @@ -607,20 +600,16 @@ namespace internal AssertIndexRange(i + 1, cell_partition_data.size()); if (cell_partition_data[i + 1] > cell_partition_data[i]) { - funct.cell(std::make_pair(cell_partition_data[i], - cell_partition_data[i + 1])); + funct.cell(i); } if (face_partition_data.empty() == false) { if (face_partition_data[i + 1] > face_partition_data[i]) - funct.face(std::make_pair(face_partition_data[i], - face_partition_data[i + 1])); + funct.face(i); if (boundary_partition_data[i + 1] > boundary_partition_data[i]) - funct.boundary( - std::make_pair(boundary_partition_data[i], - boundary_partition_data[i + 1])); + funct.boundary(i); } funct.cell_loop_post_range(i); }