From 5b431a74177484b27d4b0c88c1897a2bb656301c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 8 Dec 2019 14:59:54 +0100 Subject: [PATCH] Access neighbors during ECL --- include/deal.II/matrix_free/fe_evaluation.h | 139 ++++++++++---- include/deal.II/matrix_free/matrix_free.h | 19 ++ .../matrix_free/matrix_free.templates.h | 2 +- ...read_cell_data.cc => read_cell_data_01.cc} | 0 ...l_data.output => read_cell_data_01.output} | 0 tests/matrix_free/read_cell_data_02.cc | 180 ++++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=2.output | 7 + 7 files changed, 312 insertions(+), 35 deletions(-) rename tests/matrix_free/{read_cell_data.cc => read_cell_data_01.cc} (100%) rename tests/matrix_free/{read_cell_data.output => read_cell_data_01.output} (100%) create mode 100644 tests/matrix_free/read_cell_data_02.cc create mode 100644 tests/matrix_free/read_cell_data_02.with_mpi=true.with_p4est=true.mpirun=2.output diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index a654917d67..e8bb47f928 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -611,6 +611,13 @@ public: const AlignedVector> &array) const; + /** + * Return the id of the cells this FEEvaluation or FEFaceEvaluation is + * associated with. + */ + std::array + get_cell_ids() const; + //@} /** @@ -3555,6 +3562,82 @@ FEEvaluationBase:: return jacobian[q_index]; } +template +std::array +FEEvaluationBase:: + get_cell_ids() const +{ + const unsigned int v_len = VectorizedArrayType::n_array_elements; + std::array cells; + + // initialize array + for (unsigned int i = 0; i < v_len; ++i) + cells[i] = numbers::invalid_unsigned_int; + + if ((is_face == false) || + (is_face && + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face)) + { + // cell or interior face face (element-centric loop) + for (unsigned int i = 0; i < v_len; ++i) + cells[i] = cell * v_len + i; + } + else if (is_face && + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && + this->is_interior_face == false) + { + // exterior face (element-centric loop): for this case, we need to + // look into the FaceInfo field that collects information from both + // sides of a face once for the global mesh, and pick the face id that + // is not the local one (cell_this). + for (unsigned int i = 0; i < v_len; i++) + { + // compute actual (non vectorized) cell ID + const unsigned int cell_this = this->cell * v_len + i; + // compute face ID + unsigned int face_index = + this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell, + this->face_no, + i); + + if (face_index == numbers::invalid_unsigned_int) + continue; // invalid face ID: no neighbor on boundary + + // get cell ID on both sides of face + auto cell_m = this->matrix_info->get_face_info(face_index / v_len) + .cells_interior[face_index % v_len]; + auto cell_p = this->matrix_info->get_face_info(face_index / v_len) + .cells_exterior[face_index % v_len]; + + // compare the IDs with the given cell ID + if (cell_m == cell_this) + cells[i] = cell_p; // neighbor has the other ID + else if (cell_p == cell_this) + cells[i] = cell_m; + } + } + else if (is_face) + { + // face-centric faces + const unsigned int *cells_ = + is_interior_face ? + &this->matrix_info->get_face_info(cell).cells_interior[0] : + &this->matrix_info->get_face_info(cell).cells_exterior[0]; + for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i) + if (cells_[i] != numbers::invalid_unsigned_int) + cells[i] = cells_[i]; + } + + return cells; +} + template :: Assert(matrix_info != nullptr, ExcNotImplemented()); AssertDimension(array.size(), matrix_info->get_task_info().cell_partition_data.back()); - if (is_face) - { - VectorizedArrayType out = make_vectorized_array(Number(1.)); - const unsigned int *cells = - is_interior_face ? - &this->matrix_info->get_face_info(cell).cells_interior[0] : - &this->matrix_info->get_face_info(cell).cells_exterior[0]; - for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i) - if (cells[i] != numbers::invalid_unsigned_int) - out[i] = array[cells[i] / VectorizedArrayType::n_array_elements] - [cells[i] % VectorizedArrayType::n_array_elements]; - return out; - } - else - return array[cell]; + + // 1) collect ids of cell + const auto cells = this->get_cell_ids(); + + // 2) actually gather values + VectorizedArrayType out = make_vectorized_array(Number(1.)); + for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i) + if (cells[i] != numbers::invalid_unsigned_int) + out[i] = array[cells[i] / VectorizedArrayType::n_array_elements] + [cells[i] % VectorizedArrayType::n_array_elements]; + return out; } @@ -3603,21 +3682,17 @@ FEEvaluationBase:: Assert(matrix_info != nullptr, ExcNotImplemented()); AssertDimension(array.size(), matrix_info->get_task_info().cell_partition_data.back()); - if (is_face) - { - std::array out; - const unsigned int * cells = - is_interior_face ? - &this->matrix_info->get_face_info(cell).cells_interior[0] : - &this->matrix_info->get_face_info(cell).cells_exterior[0]; - for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i) - if (cells[i] != numbers::invalid_unsigned_int) - out[i] = array[cells[i] / VectorizedArrayType::n_array_elements] - [cells[i] % VectorizedArrayType::n_array_elements]; - return out; - } - else - return array[cell]; + + // 1) collect ids of cell + const auto cells = this->get_cell_ids(); + + // 2) actually gather values + std::array out; + for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i) + if (cells[i] != numbers::invalid_unsigned_int) + out[i] = array[cells[i] / VectorizedArrayType::n_array_elements] + [cells[i] % VectorizedArrayType::n_array_elements]; + return out; } @@ -7313,10 +7388,6 @@ FEFaceEvaluationmapped_geometry == nullptr, ExcMessage("FEEvaluation was initialized without a matrix-free object." " Integer indexing is not possible")); - Assert(this->is_interior_face == true, - ExcMessage( - "Cell-based FEFaceEvaluation::reinit only possible for the " - "interior face with second argument to constructor as true")); if (this->mapped_geometry != nullptr) return; Assert(this->matrix_info != nullptr, ExcNotInitialized()); diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index df97932f76..4907ffe549 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -1778,6 +1778,15 @@ public: VectorizedArrayType::n_array_elements> & get_face_info(const unsigned int face_batch_number) const; + + /** + * Return the table that translates a triple of the macro cell number, + * the index of a face within a cell and the index within the cell batch of + * vectorization into the index within the faces array. + */ + const Table<3, unsigned int> & + get_cell_and_face_to_plain_faces() const; + /** * Obtains a scratch data object for internal use. Make sure to release it * afterwards by passing the pointer you obtain from this object to the @@ -2465,6 +2474,16 @@ MatrixFree::get_face_info( +template +inline const Table<3, unsigned int> & +MatrixFree::get_cell_and_face_to_plain_faces() + const +{ + return face_info.cell_and_face_to_plain_faces; +} + + + template inline const Quadrature & MatrixFree::get_quadrature( diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 8db758b14e..b52d71be14 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1484,7 +1484,7 @@ MatrixFree::initialize_indices( true); face_info.cell_and_face_boundary_id.fill(numbers::invalid_boundary_id); - for (unsigned int f = 0; f < task_info.boundary_partition_data.back(); + for (unsigned int f = 0; f < task_info.ghost_face_partition_data.back(); ++f) for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements && face_info.faces[f].cells_interior[v] != diff --git a/tests/matrix_free/read_cell_data.cc b/tests/matrix_free/read_cell_data_01.cc similarity index 100% rename from tests/matrix_free/read_cell_data.cc rename to tests/matrix_free/read_cell_data_01.cc diff --git a/tests/matrix_free/read_cell_data.output b/tests/matrix_free/read_cell_data_01.output similarity index 100% rename from tests/matrix_free/read_cell_data.output rename to tests/matrix_free/read_cell_data_01.output diff --git a/tests/matrix_free/read_cell_data_02.cc b/tests/matrix_free/read_cell_data_02.cc new file mode 100644 index 0000000000..d7aad75e24 --- /dev/null +++ b/tests/matrix_free/read_cell_data_02.cc @@ -0,0 +1,180 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// tests matrix-free read_cell_data for each cell (incl. inner and outer faces, +// which is needed for element-centric-loops) + + +#include +#include + +#include + +#include + +#include +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template > +class Test +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + + void + run(unsigned int fe_degree) + { + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingQ mapping(1); + + QGauss<1> quad(fe_degree + 1); + + AffineConstraints constraint; + + typename MatrixFree::AdditionalData + additional_data; + additional_data.mapping_update_flags = update_values; + additional_data.mapping_update_flags_inner_faces = update_values; + additional_data.mapping_update_flags_boundary_faces = update_values; + additional_data.mapping_update_flags_faces_by_cells = update_values; + additional_data.hold_all_faces_to_owned_cells = true; + + MatrixFree matrix_free; + matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); + + VectorType src, dst; + + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + + this->setup_vector(matrix_free); + + matrix_free.cell_loop(&Test::cell_wise_operation, this, dst, src); + + deallog << std::endl; + } + +private: + void + setup_vector(const MatrixFree &data) + { + unsigned int n_cells = data.n_cell_batches() + data.n_ghost_cell_batches(); + cell_ids.resize(n_cells); + + for (unsigned int cell = 0; cell < n_cells; cell++) + for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell); + lane++) + cell_ids[cell][lane] = data.get_cell_iterator(cell, lane)->id(); + } + + void + cell_wise_operation(const MatrixFree &data, + VectorType &, + const VectorType &, + const std::pair &pair) const + { + FEEvaluation fe_eval(data); + FEFaceEvaluation fe_eval_m(data, true); + FEFaceEvaluation fe_eval_p(data, false); + + for (auto cell = pair.first; cell < pair.second; cell++) + { + for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell); + lane++) + { + fe_eval.reinit(cell); + const auto cell_data = fe_eval.read_cell_data(cell_ids); + Assert(cell_data[lane] == data.get_cell_iterator(cell, lane)->id(), + ExcInternalError()); + // deallog << "c " << cell_data[lane] << std::endl; + + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; + ++face) + { + fe_eval_m.reinit(cell, face); + fe_eval_p.reinit(cell, face); + + const auto cell_data_m = fe_eval_m.read_cell_data(cell_ids); + Assert(cell_data_m[lane] == + data.get_cell_iterator(cell, lane)->id(), + ExcInternalError()); + // deallog << "m" << cell_data_m[lane] << std::endl; + + const auto bids = + data.get_faces_by_cells_boundary_id(cell, face); + + if (bids[lane] != numbers::invalid_unsigned_int) + continue; + + const auto cell_data_p = fe_eval_p.read_cell_data(cell_ids); + + Assert( + cell_data_p[lane] == + data.get_cell_iterator(cell, lane)->neighbor(face)->id(), + ExcInternalError()); + // deallog << "p " << cell_data_p[lane] << std::endl; + } + } + } + } + + AlignedVector> + cell_ids; +}; + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + MPILogInitAll all; + + { + deallog.push("2d"); + Test<2> runner; + runner.run(1); + deallog.pop(); + } + { + deallog.push("3d"); + Test<3> runner; + runner.run(1); + deallog.pop(); + } +} diff --git a/tests/matrix_free/read_cell_data_02.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/matrix_free/read_cell_data_02.with_mpi=true.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..a924a83442 --- /dev/null +++ b/tests/matrix_free/read_cell_data_02.with_mpi=true.with_p4est=true.mpirun=2.output @@ -0,0 +1,7 @@ + +DEAL:0:2d:: +DEAL:0:3d:: + +DEAL:1:2d:: +DEAL:1:3d:: + -- 2.39.5