From: Maximilian Bergbauer Date: Fri, 4 Oct 2024 12:45:51 +0000 (+0200) Subject: Allow access of exterior dof values for face-centric FEFaceEvaluation X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=47ef636ac5206256b3135cf2f8f668c1ed346381;p=dealii.git Allow access of exterior dof values for face-centric FEFaceEvaluation --- diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index d263fdfeee..9109769da1 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -262,7 +262,8 @@ namespace internal template void compute_face_index_compression( - const std::vector> &faces); + const std::vector> &faces, + bool hold_all_faces_to_owned_cells); /** * This function computes the connectivity of the currently stored diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 942aa7df2e..2c4b2dcc34 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -243,7 +243,8 @@ namespace internal template void DoFInfo::compute_face_index_compression( - const std::vector> &faces) + const std::vector> &faces, + const bool hold_all_faces_to_owned_cells) { AssertDimension(length, vectorization_length); @@ -256,20 +257,33 @@ namespace internal n_vectorization_lanes_filled[dof_access_face_interior].resize( faces.size()); - // all interior faces come before the boundary faces - unsigned int n_exterior_faces = 0; - for (; n_exterior_faces < faces.size(); ++n_exterior_faces) - if (faces[n_exterior_faces].cells_exterior[0] == + // all inner faces come before the boundary faces + unsigned int n_inner_faces = 0; + for (; n_inner_faces < faces.size(); ++n_inner_faces) + if (faces[n_inner_faces].cells_exterior[0] == + numbers::invalid_unsigned_int) + break; + + // all boundary faces come after the inner faces and before the ghosted + // inner faces + unsigned int n_boundary_faces = 0; + for (; n_inner_faces + n_boundary_faces < faces.size(); + ++n_boundary_faces) + if (faces[n_inner_faces + n_boundary_faces].cells_exterior[0] != numbers::invalid_unsigned_int) break; + + const unsigned int size_exterior_faces = + hold_all_faces_to_owned_cells ? faces.size() : n_inner_faces; + index_storage_variants[dof_access_face_exterior].resize( - n_exterior_faces, IndexStorageVariants::full); + size_exterior_faces, IndexStorageVariants::full); dof_indices_contiguous[dof_access_face_exterior].resize( - n_exterior_faces * length, numbers::invalid_unsigned_int); + size_exterior_faces * length, numbers::invalid_unsigned_int); dof_indices_interleave_strides[dof_access_face_exterior].resize( faces.size() * length, numbers::invalid_unsigned_int); n_vectorization_lanes_filled[dof_access_face_exterior].resize( - n_exterior_faces); + size_exterior_faces); for (unsigned int face = 0; face < faces.size(); ++face) { @@ -362,7 +376,9 @@ namespace internal face_computation(dof_access_face_interior, faces[face].cells_interior); - if (face < n_exterior_faces) + if (face < n_inner_faces || + (hold_all_faces_to_owned_cells && + face >= (n_inner_faces + n_boundary_faces))) face_computation(dof_access_face_exterior, faces[face].cells_exterior); } diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h index 6a1ef52ffb..6fe9850926 100644 --- a/include/deal.II/matrix_free/face_setup_internal.h +++ b/include/deal.II/matrix_free/face_setup_internal.h @@ -852,7 +852,10 @@ namespace internal else if (face_is_owned[dcell->face(f) ->child(c) ->index()] == - FaceCategory::ghosted) + FaceCategory::ghosted || + face_is_owned[dcell->face(f) + ->index()] == + FaceCategory::ghosted) { inner_ghost_faces.push_back(create_face( neighbor_face_no, @@ -863,15 +866,11 @@ namespace internal is_mixed_mesh)); } else - Assert( - face_is_owned[dcell->face(f) - ->index()] == - FaceCategory:: - locally_active_done_elsewhere || - face_is_owned[dcell->face(f) - ->index()] == - FaceCategory::ghosted, - ExcInternalError()); + Assert(face_is_owned[dcell->face(f) + ->index()] == + FaceCategory:: + locally_active_done_elsewhere, + ExcInternalError()); } else { diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index f7c6d6ff61..ee9430614e 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -2026,7 +2026,8 @@ MatrixFree::initialize_indices( task_info.refinement_edge_face_partition_data[0])); for (auto &di : dof_info) - di.compute_face_index_compression(face_info.faces); + di.compute_face_index_compression( + face_info.faces, additional_data.hold_all_faces_to_owned_cells); // build the inverse map back from the faces array to // cell_and_face_to_plain_faces diff --git a/source/matrix_free/dof_info.cc b/source/matrix_free/dof_info.cc index 71230646b8..4804548a36 100644 --- a/source/matrix_free/dof_info.cc +++ b/source/matrix_free/dof_info.cc @@ -1589,19 +1589,24 @@ namespace internal template void DoFInfo::compute_face_index_compression<1>( - const std::vector> &); + const std::vector> &, + bool); template void DoFInfo::compute_face_index_compression<2>( - const std::vector> &); + const std::vector> &, + bool); template void DoFInfo::compute_face_index_compression<4>( - const std::vector> &); + const std::vector> &, + bool); template void DoFInfo::compute_face_index_compression<8>( - const std::vector> &); + const std::vector> &, + bool); template void DoFInfo::compute_face_index_compression<16>( - const std::vector> &); + const std::vector> &, + bool); template void DoFInfo::compute_vector_zero_access_pattern<1>( diff --git a/tests/matrix_free/matrix_vector_faces_24.cc b/tests/matrix_free/matrix_vector_faces_24.cc index 20bbe9c8a6..e96250eaf3 100644 --- a/tests/matrix_free/matrix_vector_faces_24.cc +++ b/tests/matrix_free/matrix_vector_faces_24.cc @@ -103,4 +103,8 @@ test() mf.vmult(out, in); deallog << "Norm of result: " << out.l2_norm() << std::endl; + + mf.manual_loop_vmult(out, in); + + deallog << "Norm of result: " << out.l2_norm() << std::endl; } diff --git a/tests/matrix_free/matrix_vector_faces_24.with_p4est=true.mpirun=2.output b/tests/matrix_free/matrix_vector_faces_24.with_p4est=true.mpirun=2.output index 71f9902dab..ae81bd819f 100644 --- a/tests/matrix_free/matrix_vector_faces_24.with_p4est=true.mpirun=2.output +++ b/tests/matrix_free/matrix_vector_faces_24.with_p4est=true.mpirun=2.output @@ -1,5 +1,7 @@ DEAL:2d::Testing FE_DGQ<2>(1) DEAL:2d::Norm of result: 4.18 +DEAL:2d::Norm of result: 4.18 DEAL:3d::Testing FE_DGQ<3>(1) DEAL:3d::Norm of result: 3.39 +DEAL:3d::Norm of result: 3.39 diff --git a/tests/matrix_free/matrix_vector_faces_common.h b/tests/matrix_free/matrix_vector_faces_common.h index 052becdabd..97fb5ad5de 100644 --- a/tests/matrix_free/matrix_vector_faces_common.h +++ b/tests/matrix_free/matrix_vector_faces_common.h @@ -87,6 +87,32 @@ public: src); } + void + manual_loop_vmult(VectorType &dst, const VectorType &src) const + { + src.update_ghost_values(); + dst = 0; + local_apply(data, dst, src, std::make_pair(0, data.n_cell_batches())); + local_apply_face(data, + dst, + src, + std::make_pair(0, data.n_inner_face_batches())); + local_apply_boundary_face(data, + dst, + src, + std::make_pair(data.n_inner_face_batches(), + data.n_inner_face_batches() + + data.n_boundary_face_batches())); + local_apply_face(data, + dst, + src, + std::make_pair(data.n_inner_face_batches() + + data.n_boundary_face_batches(), + data.n_inner_face_batches() + + data.n_boundary_face_batches() + + data.n_ghost_inner_face_batches())); + } + private: void local_apply(const MatrixFree &data,