From d79b9acb470d5fbfdb420b72f77c5631f527add8 Mon Sep 17 00:00:00 2001 From: Johannes Heinz <johannes.heinz@tuwien.ac.at> Date: Mon, 27 Nov 2023 09:55:52 +0100 Subject: [PATCH] Add FERemoteEvaluation --- doc/news/changes/major/20231127HeinzMunch | 8 + .../matrix_free/fe_remote_evaluation.h | 1037 +++++++++++++++++ tests/matrix_free/fe_remote_evaluation_01.cc | 424 +++++++ ...luation_01.with_p4est=true.mpirun=1.output | 2 + ...luation_01.with_p4est=true.mpirun=2.output | 5 + tests/matrix_free/fe_remote_evaluation_02.cc | 428 +++++++ ...luation_02.with_p4est=true.mpirun=1.output | 2 + ...luation_02.with_p4est=true.mpirun=2.output | 5 + 8 files changed, 1911 insertions(+) create mode 100644 doc/news/changes/major/20231127HeinzMunch create mode 100644 include/deal.II/matrix_free/fe_remote_evaluation.h create mode 100644 tests/matrix_free/fe_remote_evaluation_01.cc create mode 100644 tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=1.output create mode 100644 tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=2.output create mode 100644 tests/matrix_free/fe_remote_evaluation_02.cc create mode 100644 tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=1.output create mode 100644 tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=2.output diff --git a/doc/news/changes/major/20231127HeinzMunch b/doc/news/changes/major/20231127HeinzMunch new file mode 100644 index 0000000000..2fadf547d2 --- /dev/null +++ b/doc/news/changes/major/20231127HeinzMunch @@ -0,0 +1,8 @@ +New: Introduce FERemoteEvaluation. This new class provides similar +interfaces as FEEvaluation, FEFaceEvaluation, and FEPointEvaluation. +It can be used to access values and/or gradients in quadrature points +at remote triangulations. In particular the class is designed such that +it can also be used for Nitsche-type mortaring or L2-Projection with +CGAL in the future. +<br> +(Johannes Heinz, Peter Munch, 2023/11/27) diff --git a/include/deal.II/matrix_free/fe_remote_evaluation.h b/include/deal.II/matrix_free/fe_remote_evaluation.h new file mode 100644 index 0000000000..94ef418901 --- /dev/null +++ b/include/deal.II/matrix_free/fe_remote_evaluation.h @@ -0,0 +1,1037 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 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. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_fe_remote_evaluation_h +#define dealii_matrix_free_fe_remote_evaluation_h + +#include <deal.II/base/mpi_remote_point_evaluation.h> + +#include <deal.II/matrix_free/fe_evaluation.h> +#include <deal.II/matrix_free/fe_point_evaluation.h> +#include <deal.II/matrix_free/matrix_free.h> + +#include <deal.II/numerics/vector_tools.h> + +#include <variant> + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + /** + * A class that stores values and/or gradients, e.g. at quadrature points that + * are accessed inside matrix-free loops. + */ + template <int dim, int n_components, typename value_type_> + struct PrecomputedEvaluationData + { + using value_type = typename internal::FEPointEvaluation:: + EvaluatorTypeTraits<dim, n_components, value_type_>::value_type; + + using gradient_type = typename internal::FEPointEvaluation:: + EvaluatorTypeTraits<dim, n_components, value_type_>::gradient_type; + + /** + * Values at quadrature points. + */ + AlignedVector<value_type> values; + + /** + * Gradients at quadrature points. + */ + AlignedVector<gradient_type> gradients; + }; + + /** + * A class that stores a CRS like structure to access + * PrecomputedEvaluationData. In most cases a one-level CRS structure + * is enough. In this case @c ptrs has to be constructed and the shift can be + * obtained with `get_shift(index)`. The field @c ptrs_ptrs stays empty. + * It is only filled if a two level structure is needed. In this case + * `get_shift(cell_index, face_number)` return the correct shift. + */ + struct PrecomputedEvaluationDataView + { + /** + * Get a pointer to data at @p index. + */ + unsigned int + get_shift(const unsigned int index) const; + + /** + * Get a pointer to data at (@p cell_index, @p face_number). + */ + unsigned int + get_shift(const unsigned int cell_index, + const unsigned int face_number) const; + + /** + * Get the number of stored values. + */ + unsigned int + size() const; + + /** + * This parameter can be used if indices do not start with 0. + */ + unsigned int start = 0; + + /** + * Pointers to @c ptrs. + */ + std::vector<unsigned int> ptrs_ptrs; + + /** + * Pointers to data. + */ + std::vector<unsigned int> ptrs; + }; + + /** + * A class helps to access PecomputedEvaluationData in a thread-safe + * manner. Each thread has to create this wrapper class on its own to + * avoid race-conditions during reinit(). + */ + template <int dim, int n_components, typename value_type_> + class PrecomputedEvaluationDataAccessor + { + using value_type = + typename PrecomputedEvaluationData<dim, n_components, value_type_>:: + value_type; + using gradient_type = + typename PrecomputedEvaluationData<dim, n_components, value_type_>:: + gradient_type; + + public: + /** + * Constructor. + */ + PrecomputedEvaluationDataAccessor( + const PrecomputedEvaluationData<dim, n_components, value_type_> &data, + const PrecomputedEvaluationDataView &view); + + /** + * Get the value at quadrature point @p q. The entity on which the values + * are defined is set via `reinit()`. + * + * @param[in] q Quadrature point at which the value is queried. + */ + const value_type + get_value(const unsigned int q) const; + + /** + * Get the gradients at quadrature point @p q. The entity on which the + * gradients are defined is set via `reinit()`. + * + * @param[in] q Quadrature point at which the gradient is queried. + */ + const gradient_type + get_gradient(const unsigned int q) const; + + /** + * This function has to be called before `get_value()` and/or + * `get_gradient()`. + * + * @param[in] index Entity index at which quadrature points are accessed. + * This can be, e.g., a cell index, a cell batch index, or a face batch + * index. + */ + void + reinit(const unsigned int index); + + /** + * Set cell and face number at which quadrature points are accessed. + * + * @param[in] index_0 cell index. + * @param[in] index_1 cell-local face number. + */ + void + reinit(const unsigned int index_0, const unsigned int index_1); + + private: + /** + * PrecomputedEvaluationDataView provides information where values are + * located. + */ + const PrecomputedEvaluationDataView &view; + + /** + * PrecomputedEvaluationData stores the actual values. + */ + const PrecomputedEvaluationData<dim, n_components, value_type_> &data; + + /** + * Offset to data after last call of `reinit()`. + */ + unsigned int data_offset; + }; + +} // namespace internal + + + +/** + * Communication objects know about the communication pattern. + * In case of (matrix-free) cells batches or faces batches + * a RemotePointEvaluation object stores the location of the + * remote points. @c batch_id_n_entities relates these points + * to the corresponding quadrature points of entity batches. + * For this the field stores batch IDs and the number of entities + * in the batch. + */ +template <int dim> +struct FERemoteCommunicationObjectEntityBatches +{ + /** + * Object that is reinitialized with the remote points we want to access. + */ + std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>> rpe; + + /** + * A vector that stores the batch IDs and the number of active entries in the + * batch related to the the remote points. + */ + std::vector<std::pair<unsigned int, unsigned int>> batch_id_n_entities; + + /** + * Function that gives access to @c batch_id_n_entities. The function is only + * used for a unified access to the vector that stores additional information + * on the points that are processed via @c rpe. + */ + std::vector<std::pair<unsigned int, unsigned int>> + get_communication_object_pntrs() const; +}; + +/** + * Similar as @c FERemoteCommunicationObjectEntityBatches. + * To relate the points from @c RemotePointEvaluation to + * quadrature points on corresponding cells, cell iterators + * have to be stored. + */ +template <int dim> +struct FERemoteCommunicationObjectCells +{ + /** + * Object that is reinitialized with the remote points we want to access. + */ + std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>> rpe; + + /** + * A vector that stores the cell iterators related to the the remote points. + */ + std::vector<typename Triangulation<dim>::cell_iterator> cells; + + /** + * Function that gives access to @c cells. The function is only used + * for a unified access to the vector that stores additional information on + * the points that are processed via @c rpe. + */ + std::vector<typename Triangulation<dim>::cell_iterator> + get_communication_object_pntrs() const; +}; + +/** + * Similar as @c FERemoteCommunicationObjectCells. + * To relate the points from @c RemotePointEvaluation to + * quadrature points on corresponding faces, cell iterators + * and face numbers have to be stored. + */ +template <int dim> +struct FERemoteCommunicationObjectFaces +{ + /** + * Object that is reinitialized with the remote points we want to access. + */ + std::shared_ptr<Utilities::MPI::RemotePointEvaluation<dim>> rpe; + + /** + * A vector that stores the cell iterators and the face number related to the + * the remote points. + */ + std::vector< + std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>> + cell_face_nos; + + /** + * Function that gives access to @c cell_face_nos. The function is only used + * for a unified access to the vector that stores additional information on + * the points that are processed via @c rpe. + */ + std::vector< + std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>> + get_communication_object_pntrs() const; +}; + +/** + * A class to fill the fields in PrecomputedEvaluationData. + */ +template <int dim> +class FERemoteEvaluationCommunicator : public Subscriptor +{ +public: + /** + * This function stores given communication objects + * and constructs a @c PrecomputedEvaluationDataView + * object if remote points are related to matrix-free face batches. + */ + void + reinit_faces(const std::vector<FERemoteCommunicationObjectEntityBatches<dim>> + &comm_objects, + const std::pair<unsigned int, unsigned int> &face_batch_range, + const std::vector<Quadrature<dim>> &quadrature_vector); + + /** + * This function stores given communication objects + * and constructs a @c PrecomputedEvaluationDataView + * object if remote points are related to faces of given + * cells. + */ + template <typename Iterator> + void + reinit_faces( + const std::vector<FERemoteCommunicationObjectFaces<dim>> &comm_objects, + const IteratorRange<Iterator> &cell_iterator_range, + const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector); + + /** + * Fill the fields stored in PrecomputedEvaluationData. + */ + template <int n_components, + typename PrecomputedEvaluationDataType, + typename MeshType, + typename VectorType> + void + update_ghost_values( + PrecomputedEvaluationDataType &dst, + const MeshType &mesh, + const VectorType &src, + const EvaluationFlags::EvaluationFlags eval_flags, + const unsigned int first_selected_component, + const VectorTools::EvaluationFlags::EvaluationFlags vec_flags) const; + + /** + * Provide access to @c PrecomputedEvaluationDataView. + */ + const internal::PrecomputedEvaluationDataView & + get_view() const; + +private: + /** + * CRS like data structure that describes the data positions at given + * indices. + */ + internal::PrecomputedEvaluationDataView view; + + /** + * A variant for all possible communication objects. + */ + std::vector<std::variant<FERemoteCommunicationObjectEntityBatches<dim>, + FERemoteCommunicationObjectCells<dim>, + FERemoteCommunicationObjectFaces<dim>>> + communication_objects; + + /** + * Functions in this class only deals with copying data to @c + * PrecomputedEvaluationData. + */ + class CopyInstructions + { + public: + /** + * Copy data from @c src to @c dst. Overload for @c + * FERemoteCommunicationObjectCells. + */ + template <typename T1, typename T2> + static void + copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector<T1> &dst, + const std::vector<T2> &src, + const std::vector<typename Triangulation<dim>::cell_iterator> &cells); + + /** + * Copy data from @c src to @c dst. Overload for @c + * FERemoteCommunicationObjectFaces. + */ + template <typename T1, typename T2> + static void + copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector<T1> &dst, + const std::vector<T2> &src, + const std::vector<std::pair<typename Triangulation<dim>::cell_iterator, + unsigned int>> &cell_face_nos); + + /** + * Copy data from @c src to @c dst. Overload for @c + * FERemoteCommunicationObjectEntityBatches. + */ + template <typename T1, typename T2> + static void + copy_data(const internal::PrecomputedEvaluationDataView &view, + AlignedVector<T1> &dst, + const std::vector<T2> &src, + const std::vector<std::pair<unsigned int, unsigned int>> + &batch_id_n_entities); + + private: + /** + * Copy data to the correct position in a @c VectorizedArray. + */ + template <typename T1, std::size_t n_lanes> + static void + copy_data_entries(VectorizedArray<T1, n_lanes> &dst, + const unsigned int v, + const T1 &src); + + /** + * Similar as @c copy_data_entries() above. + */ + template <typename T1, int rank_, std::size_t n_lanes, int dim_> + static void + copy_data_entries(Tensor<rank_, dim_, VectorizedArray<T1, n_lanes>> &dst, + const unsigned int v, + const Tensor<rank_, dim_, T1> &src); + + /** + * Similar as @c copy_data_entries() above. + */ + template <typename T1, + int rank_, + std::size_t n_lanes, + int n_components_, + int dim_> + static void + copy_data_entries( + Tensor<rank_, + n_components_, + Tensor<rank_, dim_, VectorizedArray<T1, n_lanes>>> &dst, + const unsigned int v, + const Tensor<rank_, n_components_, Tensor<rank_, dim_, T1>> &src); + + /** + * Throw a runtime exception if @c copy_data_entries() has not been + * implemented for a given type. + */ + template <typename T1, typename T2> + static void + copy_data_entries(T1 &, const unsigned int, const T2 &); + }; +}; + + + +/** + * Class to access data in a matrix-free loop for non-matching discretizations. + * Interfaces are named with FEEvaluation in mind. + * The main difference is, that `gather_evaluate()` updates and caches + * all values at once. Therefore, it has to be called on one thread before a + * matrix-free loop. + * + * To access values and gradients in a thread safe way, @c get_data_accessor() + * has to be called on every thread. It provides the functions `get_value()` and + * `get_gradient()`. + */ +template <int dim, int n_components, typename value_type> +class FERemoteEvaluation +{ +public: + /** + * The constructor needs a corresponding FERemoteEvaluationCommunicator + * which has to be setup outside of this class. This design choice is + * motivated since the same FERemoteEvaluationCommunicator can be used + * for different MeshTypes and number of components. + * + * @param[in] comm FERemoteEvaluationCommunicator. + * @param[in] mesh Triangulation or DoFHandler. + * @param[in] evaluation_flags Specify treatment of values at points which are + * found on multiple cells. + * @param[in] first_selected_component Select first component of evaluation in + * DoFHandlers with multiple components. + */ + template <typename MeshType> + FERemoteEvaluation(const FERemoteEvaluationCommunicator<dim> &comm, + const MeshType &mesh, + const unsigned int first_selected_component = 0, + const VectorTools::EvaluationFlags::EvaluationFlags + evaluation_flags = VectorTools::EvaluationFlags::avg); + + /** + * Update the data which can be accessed via `get_value()` and + * `get_gradient()`. + * + * @param[in] src Solution vector used to update data. + * @param[in] flags Evaluation flags. Currently supported are + * EvaluationFlags::values and EvaluationFlags::gradients. + */ + template <typename VectorType> + void + gather_evaluate(const VectorType &src, + const EvaluationFlags::EvaluationFlags flags); + + /** + * @c FERemoteEvaluation does not provide the functions `get_value()` and + *`get_gradient()`. To access values and/or gradients call @c + *get_data_accessor() on every thread, e.g., `auto remote_evaluator = + *get_data_accessor();` The returned object can be used as follows. + *@code + * for (unsigned int entity = range.first; entity < range.second; ++entity) + * { + * remote_evaluator.reinit(entity); + * for(unsigned int q : quadrature_point_indices()) + * remote_evaluator.get_value(q) + * } + *@endcode + */ + internal::PrecomputedEvaluationDataAccessor<dim, n_components, value_type> + get_data_accessor() const; + +private: + /** + * Use Triangulation as MeshType. + */ + void + set_mesh(const Triangulation<dim> &tria); + + /** + * Use DoFHandler as MeshType. + */ + void + set_mesh(const DoFHandler<dim> &dof_handler); + + /** + * Precomputed values and/or gradients at remote locations. + */ + internal::PrecomputedEvaluationData<dim, n_components, value_type> data; + + /** + * Underlying communicator which handles update of the ghost values. + */ + SmartPointer<const FERemoteEvaluationCommunicator<dim>> comm; + + /** + * Pointer to MeshType if used with Triangulation. + */ + SmartPointer<const Triangulation<dim>> tria; + + /** + * Pointer to MeshType if used with DoFHandler. + */ + SmartPointer<const DoFHandler<dim>> dof_handler; + + /** + * First selected component. + */ + const unsigned int first_selected_component; + + /** + * Flags that indicate which ghost values are updated. + */ + const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags; +}; + + + +namespace internal +{ + unsigned int + PrecomputedEvaluationDataView::get_shift(const unsigned int index) const + { + Assert(ptrs_ptrs.size() == 0, ExcMessage("Two level CRS set up")); + + Assert(index != numbers::invalid_unsigned_int, + ExcMessage("Index has to be valid!")); + + Assert(start <= index, ExcInternalError()); + AssertIndexRange(index - start, ptrs.size()); + return ptrs[index - start]; + } + + unsigned int + PrecomputedEvaluationDataView::get_shift(const unsigned int cell_index, + const unsigned int face_number) const + { + Assert(ptrs_ptrs.size() > 0, ExcMessage("No two level CRS set up")); + + Assert(cell_index != numbers::invalid_unsigned_int, + ExcMessage("Cell index has to be valid!")); + Assert(face_number != numbers::invalid_unsigned_int, + ExcMessage("Face number has to be valid!")); + + Assert(start <= cell_index, ExcInternalError()); + + AssertIndexRange(cell_index - start, ptrs_ptrs.size()); + const unsigned int face_index = ptrs_ptrs[cell_index - start] + face_number; + AssertIndexRange(face_index, ptrs.size()); + return ptrs[face_index]; + } + + unsigned int + PrecomputedEvaluationDataView::size() const + { + Assert(ptrs.size() > 0, ExcInternalError()); + return ptrs.back(); + } + + template <int dim, int n_components, typename value_type_> + PrecomputedEvaluationDataAccessor<dim, n_components, value_type_>:: + PrecomputedEvaluationDataAccessor( + const PrecomputedEvaluationData<dim, n_components, value_type_> &data, + const PrecomputedEvaluationDataView &view) + : view(view) + , data(data) + , data_offset(numbers::invalid_unsigned_int) + {} + + template <int dim, int n_components, typename value_type_> + const typename PrecomputedEvaluationData<dim, + n_components, + value_type_>::value_type + PrecomputedEvaluationDataAccessor<dim, n_components, value_type_>::get_value( + const unsigned int q) const + { + Assert(data_offset != numbers::invalid_unsigned_int, + ExcMessage("reinit() not called.")); + AssertIndexRange(data_offset + q, data.values.size()); + return data.values[data_offset + q]; + } + + template <int dim, int n_components, typename value_type_> + const typename PrecomputedEvaluationData<dim, n_components, value_type_>:: + gradient_type + PrecomputedEvaluationDataAccessor<dim, n_components, value_type_>:: + get_gradient(const unsigned int q) const + { + Assert(data_offset != numbers::invalid_unsigned_int, + ExcMessage("reinit() not called.")); + AssertIndexRange(data_offset + q, data.gradients.size()); + return data.gradients[data_offset + q]; + } + + template <int dim, int n_components, typename value_type_> + void + PrecomputedEvaluationDataAccessor<dim, n_components, value_type_>::reinit( + const unsigned int index) + { + data_offset = view.get_shift(index); + } + + template <int dim, int n_components, typename value_type_> + void + PrecomputedEvaluationDataAccessor<dim, n_components, value_type_>::reinit( + const unsigned int index_0, + const unsigned int index_1) + { + data_offset = view.get_shift(index_0, index_1); + } + +} // namespace internal + + + +template <int dim> +std::vector<std::pair<unsigned int, unsigned int>> +FERemoteCommunicationObjectEntityBatches<dim>::get_communication_object_pntrs() + const +{ + return batch_id_n_entities; +} + +template <int dim> +std::vector<typename Triangulation<dim>::cell_iterator> +FERemoteCommunicationObjectCells<dim>::get_communication_object_pntrs() const +{ + return cells; +} + +template <int dim> +std::vector<std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>> +FERemoteCommunicationObjectFaces<dim>::get_communication_object_pntrs() const +{ + return cell_face_nos; +} + + + +template <int dim> +void +FERemoteEvaluationCommunicator<dim>::reinit_faces( + const std::vector<FERemoteCommunicationObjectEntityBatches<dim>> + &comm_objects, + const std::pair<unsigned int, unsigned int> &face_batch_range, + const std::vector<Quadrature<dim>> &quadrature_vector) +{ + // erase type by converting to the base object + communication_objects.clear(); + for (const auto &co : comm_objects) + communication_objects.push_back(co); + + // fetch points and update communication patterns + const unsigned int n_cells = quadrature_vector.size(); + AssertDimension(n_cells, face_batch_range.second - face_batch_range.first); + + // construct view: + view.start = face_batch_range.first; + + view.ptrs.resize(n_cells + 1); + + view.ptrs[0] = 0; + for (unsigned int face = 0; face < n_cells; ++face) + { + view.ptrs[face + 1] = view.ptrs[face] + quadrature_vector[face].size(); + } +} + +template <int dim> +template <typename Iterator> +void +FERemoteEvaluationCommunicator<dim>::reinit_faces( + const std::vector<FERemoteCommunicationObjectFaces<dim>> &comm_objects, + const IteratorRange<Iterator> &cell_iterator_range, + const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector) +{ + // erase type + communication_objects.clear(); + for (const auto &co : comm_objects) + communication_objects.push_back(co); + + const unsigned int n_cells = quadrature_vector.size(); + AssertDimension(n_cells, + std::distance(cell_iterator_range.begin(), + cell_iterator_range.end())); + + // construct view: + auto &cell_ptrs = view.ptrs_ptrs; + auto &face_ptrs = view.ptrs; + + view.start = 0; + cell_ptrs.resize(n_cells); + unsigned int n_faces = 0; + for (const auto &cell : cell_iterator_range) + { + cell_ptrs[cell->active_cell_index()] = n_faces; + n_faces += cell->n_faces(); + } + + face_ptrs.resize(n_faces + 1); + face_ptrs[0] = 0; + for (const auto &cell : cell_iterator_range) + { + for (const auto &f : cell->face_indices()) + { + const unsigned int face_index = + cell_ptrs[cell->active_cell_index()] + f; + + face_ptrs[face_index + 1] = + face_ptrs[face_index] + + quadrature_vector[cell->active_cell_index()][f].size(); + } + } +} + +template <int dim> +template <int n_components, + typename PrecomputedEvaluationDataType, + typename MeshType, + typename VectorType> +void +FERemoteEvaluationCommunicator<dim>::update_ghost_values( + PrecomputedEvaluationDataType &dst, + const MeshType &mesh, + const VectorType &src, + const EvaluationFlags::EvaluationFlags eval_flags, + const unsigned int first_selected_component, + const VectorTools::EvaluationFlags::EvaluationFlags vec_flags) const +{ + const bool has_ghost_elements = src.has_ghost_elements(); + + if (has_ghost_elements == false) + src.update_ghost_values(); + + + for (const auto &communication_object : communication_objects) + { + if (eval_flags & EvaluationFlags::values) + { + std::visit( + [&](const auto &obj) { + CopyInstructions::copy_data( + view, + dst.values, + VectorTools::point_values<n_components>( + *obj.rpe, mesh, src, vec_flags, first_selected_component), + obj.get_communication_object_pntrs()); + }, + communication_object); + } + + if (eval_flags & EvaluationFlags::gradients) + { + std::visit( + [&](const auto &obj) { + CopyInstructions::copy_data( + view, + dst.gradients, + VectorTools::point_gradients<n_components>( + *obj.rpe, mesh, src, vec_flags, first_selected_component), + obj.get_communication_object_pntrs()); + }, + communication_object); + } + + Assert(!(eval_flags & EvaluationFlags::hessians), ExcNotImplemented()); + } + + if (has_ghost_elements == false) + src.zero_out_ghost_values(); +} + +template <int dim> +const internal::PrecomputedEvaluationDataView & +FERemoteEvaluationCommunicator<dim>::get_view() const +{ + return view; +} + +template <int dim> +template <typename T1, typename T2> +void +FERemoteEvaluationCommunicator<dim>::CopyInstructions::copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector<T1> &dst, + const std::vector<T2> &src, + const std::vector<typename Triangulation<dim>::cell_iterator> &cells) +{ + dst.resize(view.size()); + + unsigned int c = 0; + for (const auto &cell : cells) + { + for (unsigned int j = view.get_shift(cell->active_cell_index()); + j < view.get_shift(cell->active_cell_index() + 1); + ++j, ++c) + { + AssertIndexRange(j, dst.size()); + AssertIndexRange(c, src.size()); + dst[j] = src[c]; + } + } +} + +template <int dim> +template <typename T1, typename T2> +void +FERemoteEvaluationCommunicator<dim>::CopyInstructions::copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector<T1> &dst, + const std::vector<T2> &src, + const std::vector<std::pair<typename Triangulation<dim>::cell_iterator, + unsigned int>> &cell_face_nos) +{ + dst.resize(view.size()); + + unsigned int c = 0; + for (const auto &[cell, f] : cell_face_nos) + { + for (unsigned int j = view.get_shift(cell->active_cell_index(), f); + j < view.get_shift(cell->active_cell_index(), f + 1); + ++j, ++c) + { + AssertIndexRange(j, dst.size()); + AssertIndexRange(c, src.size()); + + dst[j] = src[c]; + } + } +} + +template <int dim> +template <typename T1, typename T2> +void +FERemoteEvaluationCommunicator<dim>::CopyInstructions::copy_data( + const internal::PrecomputedEvaluationDataView &view, + AlignedVector<T1> &dst, + const std::vector<T2> &src, + const std::vector<std::pair<unsigned int, unsigned int>> &batch_id_n_entities) +{ + dst.resize(view.size()); + + unsigned int c = 0; + for (const auto &[batch_id, n_entries] : batch_id_n_entities) + { + for (unsigned int v = 0; v < n_entries; ++v) + for (unsigned int j = view.get_shift(batch_id); + j < view.get_shift(batch_id + 1); + ++j, ++c) + { + AssertIndexRange(j, dst.size()); + AssertIndexRange(c, src.size()); + + copy_data_entries(dst[j], v, src[c]); + } + } +} + +template <int dim> +template <typename T1, std::size_t n_lanes> +void +FERemoteEvaluationCommunicator<dim>::CopyInstructions::copy_data_entries( + VectorizedArray<T1, n_lanes> &dst, + const unsigned int v, + const T1 &src) +{ + AssertIndexRange(v, n_lanes); + + dst[v] = src; +} + +template <int dim> +template <typename T1, int rank_, std::size_t n_lanes, int dim_> +void +FERemoteEvaluationCommunicator<dim>::CopyInstructions::copy_data_entries( + Tensor<rank_, dim_, VectorizedArray<T1, n_lanes>> &dst, + const unsigned int v, + const Tensor<rank_, dim_, T1> &src) +{ + AssertIndexRange(v, n_lanes); + + if constexpr (rank_ == 1) + { + for (unsigned int i = 0; i < dim_; ++i) + dst[i][v] = src[i]; + } + else + { + for (unsigned int i = 0; i < rank_; ++i) + for (unsigned int j = 0; j < dim_; ++j) + dst[i][j][v] = src[i][j]; + } +} + +template <int dim> +template <typename T1, + int rank_, + std::size_t n_lanes, + int n_components_, + int dim_> +void +FERemoteEvaluationCommunicator<dim>::CopyInstructions::copy_data_entries( + Tensor<rank_, + n_components_, + Tensor<rank_, dim_, VectorizedArray<T1, n_lanes>>> &dst, + const unsigned int v, + const Tensor<rank_, n_components_, Tensor<rank_, dim_, T1>> &src) +{ + if constexpr (rank_ == 1) + { + for (unsigned int i = 0; i < n_components_; ++i) + copy_data(dst[i], v, src[i]); + } + else + { + for (unsigned int i = 0; i < rank_; ++i) + for (unsigned int j = 0; j < n_components_; ++j) + dst[i][j][v] = src[i][j]; + } +} + +template <int dim> +template <typename T1, typename T2> +void +FERemoteEvaluationCommunicator<dim>::CopyInstructions::copy_data_entries( + T1 &, + const unsigned int, + const T2 &) +{ + Assert(false, + ExcMessage( + "copy_data_entries() not implemented for given arguments.")); +} + + + +template <int dim, int n_components, typename value_type> +template <typename MeshType> +FERemoteEvaluation<dim, n_components, value_type>::FERemoteEvaluation( + const FERemoteEvaluationCommunicator<dim> &comm, + const MeshType &mesh, + const unsigned int first_selected_component, + const VectorTools::EvaluationFlags::EvaluationFlags evaluation_flags) + : comm(&comm) + , first_selected_component(first_selected_component) + , evaluation_flags(evaluation_flags) +{ + set_mesh(mesh); +} + +template <int dim, int n_components, typename value_type> +template <typename VectorType> +void +FERemoteEvaluation<dim, n_components, value_type>::gather_evaluate( + const VectorType &src, + const EvaluationFlags::EvaluationFlags flags) +{ + if (tria) + { + Assert(n_components == 1, ExcNotImplemented()); + comm->template update_ghost_values<n_components>(this->data, + *tria, + src, + flags, + first_selected_component, + evaluation_flags); + } + else if (dof_handler) + { + comm->template update_ghost_values<n_components>(this->data, + *dof_handler, + src, + flags, + first_selected_component, + evaluation_flags); + } + else + Assert(false, ExcNotImplemented()); +} + +template <int dim, int n_components, typename value_type> +internal::PrecomputedEvaluationDataAccessor<dim, n_components, value_type> +FERemoteEvaluation<dim, n_components, value_type>::get_data_accessor() const +{ + internal::PrecomputedEvaluationDataAccessor data_accessor(data, + comm->get_view()); + return data_accessor; +} + +template <int dim, int n_components, typename value_type> +void +FERemoteEvaluation<dim, n_components, value_type>::set_mesh( + const Triangulation<dim> &tria) +{ + this->tria = &tria; +} + +template <int dim, int n_components, typename value_type> +void +FERemoteEvaluation<dim, n_components, value_type>::set_mesh( + const DoFHandler<dim> &dof_handler) +{ + this->dof_handler = &dof_handler; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/tests/matrix_free/fe_remote_evaluation_01.cc b/tests/matrix_free/fe_remote_evaluation_01.cc new file mode 100644 index 0000000000..278d1a6a02 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_01.cc @@ -0,0 +1,424 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 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. +// +// --------------------------------------------------------------------- + +// This test checks the access of values on boundary faces using +// FERemoteEvaluation. + +#include <deal.II/base/mpi.h> + +#include <deal.II/fe/fe_dgq.h> +#include <deal.II/fe/fe_system.h> + +#include <deal.II/grid/grid_generator.h> + +#include <deal.II/lac/la_parallel_block_vector.h> + +#include <deal.II/matrix_free/fe_remote_evaluation.h> + +#include "../tests.h" + +using namespace dealii; + +using Number = double; + +using VectorType = LinearAlgebra::distributed::Vector<Number>; + +constexpr unsigned int dim = 2; + +template <int dim> +struct Coordinates : public Function<dim> +{ + Coordinates(unsigned int offset) + : Function<dim>(dim + offset, 0.0) + , offset(offset) + {} + + // Function that returns the coordinates of given point. + // In case an offset is provided, 0.0 is returned for values + // below the offset and the coordinates are returned for + // components >= offset. This is used to check the functionality + // for first selected component. + double + value(const Point<dim> &p, const unsigned int component) const final + { + if (component < offset) + return 0.0; + return p[component - offset]; + } + + unsigned int offset; +}; + +FERemoteEvaluationCommunicator<dim> +construct_comm_for_face_batches(const MatrixFree<dim, Number> &matrix_free) +{ + // Setup Communication objects for all boundary faces + FERemoteCommunicationObjectEntityBatches<dim> co; + + // Get range of boundary face indices + const auto face_batch_range = + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + + // Fill a quadrature vector to keep track of the quadrature sizes on each + // face + std::vector<Quadrature<dim>> quadrature_vector( + matrix_free.n_boundary_face_batches()); + + // Points that are searched by rpe. + std::vector<Point<dim>> points; + + FEFaceEvaluation<dim, -1, 0, 1, Number> phi(matrix_free, true, 0, 0, 0); + for (unsigned int bface = 0; + bface < face_batch_range.second - face_batch_range.first; + ++bface) + { + const unsigned int face = face_batch_range.first + bface; + + phi.reinit(face); + + const unsigned int n_faces = + matrix_free.n_active_entries_per_face_batch(face); + + co.batch_id_n_entities.push_back(std::make_pair(face, n_faces)); + + // Append the quadrature points to the points we need to search + // for. + for (unsigned int v = 0; v < n_faces; ++v) + { + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Point<dim> temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + + points.push_back(temp); + } + } + + // append quadrature of correct size + quadrature_vector[bface] = Quadrature<dim>(phi.n_q_points); + } + + // use rpe to search for stored points + auto rpe = + std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>(1.0e-9); + + rpe->reinit(points, + matrix_free.get_dof_handler().get_triangulation(), + *matrix_free.get_mapping_info().mapping); + Assert(rpe->all_points_found(), ExcMessage("Not all remote points found.")); + Assert(rpe->is_map_unique(), ExcMessage("The map should be unique.")); + + co.rpe = rpe; + + // Renit the communicator `FERemoteEvaluationCommunicator` + // with the communication objects. + FERemoteEvaluationCommunicator<dim> remote_communicator; + remote_communicator.reinit_faces({co}, face_batch_range, quadrature_vector); + + return remote_communicator; +} + +void +test_face_batches(const MatrixFree<dim, Number> &matrix_free, + unsigned int max_first_selected_comp) +{ + // 0) Get remote communicator. + auto remote_communicator = construct_comm_for_face_batches(matrix_free); + + for (unsigned int first_selected_comp = 0; + first_selected_comp < max_first_selected_comp; + ++first_selected_comp) + { + // 1) Allocate vectors and interpolate grid coordinates in DoFs. + VectorType src; + VectorType dst; + matrix_free.initialize_dof_vector(src, first_selected_comp); + matrix_free.initialize_dof_vector(dst, first_selected_comp); + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(first_selected_comp), + Coordinates<dim>(first_selected_comp), + src); + + // 2) Construct remote evaluator... + FERemoteEvaluation<dim, dim, VectorizedArray<Number>> phi_remote_eval( + remote_communicator, + matrix_free.get_dof_handler(first_selected_comp), + first_selected_comp); + + // ...and precompute remote values. + phi_remote_eval.gather_evaluate(src, EvaluationFlags::values); + + // 3) Access values at quadrature points + const auto boundary_function = [&](const auto &data, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation<dim, -1, 0, dim, Number> phi( + matrix_free, true, first_selected_comp, 0, first_selected_comp); + auto phi_remote = phi_remote_eval.get_data_accessor(); + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + phi_remote.reinit(face); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto error = + phi.quadrature_point(q) - phi_remote.get_value(q); + + // only consider active entries in global error + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + AssertThrow(std::abs(error[d][v]) < 1e-13, + ExcMessage("Error too large.")); + } + } + }; + + matrix_free.template loop<VectorType, VectorType>( + {}, {}, boundary_function, dst, src, true); + } +} + + +FERemoteEvaluationCommunicator<dim> +construct_comm_for_cell_face_nos(const MatrixFree<dim, Number> &matrix_free) +{ + // Setup Communication objects for all boundary faces + FERemoteCommunicationObjectFaces<dim> co; + + // Get range of boundary face indices + const auto face_batch_range = + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + + // Fill a quadrature vector to keep track of the quadrature sizes on each + // cell, face pair + std::vector<std::vector<Quadrature<dim - 1>>> quadrature_vector; + for (const auto &cell : matrix_free.get_dof_handler() + .get_triangulation() + .active_cell_iterators()) + quadrature_vector.emplace_back( + std::vector<Quadrature<dim - 1>>(cell->n_faces())); + + // Points that are searched by rpe. + std::vector<Point<dim>> points; + + FEFaceEvaluation<dim, -1, 0, 1, Number> phi(matrix_free, true, 0, 0, 0); + for (unsigned int bface = 0; + bface < face_batch_range.second - face_batch_range.first; + ++bface) + { + const unsigned int face = face_batch_range.first + bface; + phi.reinit(face); + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto [cell, f] = matrix_free.get_face_iterator(face, v, true); + co.cell_face_nos.push_back(std::make_pair(cell, f)); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Point<dim> temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + points.push_back(temp); + } + + // append quadrature of correct size + quadrature_vector[cell->active_cell_index()][f] = + Quadrature<dim - 1>(phi.n_q_points); + } + } + + // use rpe to search for stored points + auto rpe = + std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>(1.0e-9); + + rpe->reinit(points, + matrix_free.get_dof_handler().get_triangulation(), + *matrix_free.get_mapping_info().mapping); + Assert(rpe->all_points_found(), ExcMessage("Not all remote points found.")); + Assert(rpe->is_map_unique(), ExcMessage("The map should be unique.")); + + co.rpe = rpe; + + // Renit the communicator `FERemoteEvaluationCommunicator` + // with the communication objects. + FERemoteEvaluationCommunicator<dim> remote_communicator; + remote_communicator.reinit_faces( + {co}, + matrix_free.get_dof_handler().get_triangulation().active_cell_iterators(), + quadrature_vector); + + return remote_communicator; +} + +void +test_cell_face_nos(const MatrixFree<dim, Number> &matrix_free, + unsigned int max_first_selected_comp) +{ + // 0) Get remote communicator. + auto remote_communicator = construct_comm_for_cell_face_nos(matrix_free); + + for (unsigned int first_selected_comp = 0; + first_selected_comp < max_first_selected_comp; + ++first_selected_comp) + { + // 1) Allocate vectors and interpolate grid coordinates in DoFs. + VectorType src; + VectorType dst; + matrix_free.initialize_dof_vector(src, first_selected_comp); + matrix_free.initialize_dof_vector(dst, first_selected_comp); + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(first_selected_comp), + Coordinates<dim>(first_selected_comp), + src); + + // 2) Construct remote evaluator... + FERemoteEvaluation<dim, dim, Number> phi_remote_eval( + remote_communicator, + matrix_free.get_dof_handler(first_selected_comp), + first_selected_comp); + + // ...and precompute remote values. + phi_remote_eval.gather_evaluate(src, EvaluationFlags::values); + + // 3) Access values at quadrature points + const auto boundary_function = [&](const auto &data, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation<dim, -1, 0, dim, Number> phi( + matrix_free, true, first_selected_comp, 0, first_selected_comp); + auto phi_remote = phi_remote_eval.get_data_accessor(); + + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto [cell, f] = + matrix_free.get_face_iterator(face, v, true); + phi_remote.reinit(cell->active_cell_index(), f); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Tensor<1, dim, Number> temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + + const auto error = temp - phi_remote.get_value(q); + for (unsigned int d = 0; d < dim; ++d) + AssertThrow(std::abs(error[d]) < 1e-13, + ExcMessage("Error too large.")); + } + } + } + }; + + matrix_free.template loop<VectorType, VectorType>( + {}, {}, boundary_function, dst, src, true); + } +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + MPILogInitAll all; + + // 1) setup MatrixFree for different finite elements, i.e. + // FESystem(FE_DGQ, dim + first_selected_component). + + const unsigned int refinements = 2; + const unsigned int degree = 5; + + // Construct triangulation. + parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(refinements); + + // Create DoFHandlers with different number of components to + // check different first_selected components + std::vector<std::unique_ptr<DoFHandler<dim>>> dof_handlers; + unsigned int max_first_selected_comp = 3; + for (unsigned int i = 0; i < max_first_selected_comp; ++i) + { + dof_handlers.push_back(std::make_unique<DoFHandler<dim>>(tria)); + dof_handlers.back()->distribute_dofs( + FESystem<dim>(FE_DGQ<dim>(degree), dim + i)); + } + + // Constraints are not considered in this test + AffineConstraints<Number> constraints; + constraints.close(); + + // Setup MatrixFree + typename MatrixFree<dim, Number>::AdditionalData data; + data.mapping_update_flags = update_values; + // We are not considering inner faces in this test. Therefore, + // data.mapping_update_flags_inner_faces is not set. + data.mapping_update_flags_boundary_faces = + update_quadrature_points | update_values; + + + // Construct vector of raw pointers needed to init matrix free + std::vector<const DoFHandler<dim> *> dof_handlers_mf; + std::vector<const AffineConstraints<Number> *> constraints_mf; + for (auto &dh : dof_handlers) + { + dof_handlers_mf.push_back(dh.get()); + constraints_mf.push_back(&constraints); + } + + MatrixFree<dim, Number> matrix_free; + matrix_free.reinit(MappingQ1<dim>(), + dof_handlers_mf, + constraints_mf, + QGauss<dim>(degree + 1), + data); + + // 2) test to access remote values at points in face batches + test_face_batches(matrix_free, max_first_selected_comp); + + // 3) test to access remote values at points on faces of cells + test_cell_face_nos(matrix_free, max_first_selected_comp); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=1.output b/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..be8d055f86 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=1.output @@ -0,0 +1,2 @@ + +DEAL:0::OK diff --git a/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=2.output b/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..52103938b0 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_01.with_p4est=true.mpirun=2.output @@ -0,0 +1,5 @@ + +DEAL:0::OK + +DEAL:1::OK + diff --git a/tests/matrix_free/fe_remote_evaluation_02.cc b/tests/matrix_free/fe_remote_evaluation_02.cc new file mode 100644 index 0000000000..d4b5fdb14a --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_02.cc @@ -0,0 +1,428 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 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. +// +// --------------------------------------------------------------------- + +// This test checks the access of gradients on boundary faces using +// FERemoteEvaluation. + +#include <deal.II/base/mpi.h> + +#include <deal.II/fe/fe_dgq.h> +#include <deal.II/fe/fe_system.h> + +#include <deal.II/grid/grid_generator.h> + +#include <deal.II/lac/la_parallel_block_vector.h> + +#include <deal.II/matrix_free/fe_remote_evaluation.h> + +#include "../tests.h" + +using namespace dealii; + +using Number = double; + +using VectorType = LinearAlgebra::distributed::Vector<Number>; + +constexpr unsigned int dim = 2; + +template <int dim> +struct QuadraticFunction : public Function<dim> +{ + QuadraticFunction(unsigned int offset) + : Function<dim>(dim + offset, 0.0) + , offset(offset) + {} + + // Function that returns the coordinates of given point. + // In case an offset is provided, 0.0 is returned for values + // below the offset and the coordinates are returned for + // components >= offset. This is used to check the functionality + // for first selected component. + double + value(const Point<dim> &p, const unsigned int component) const final + { + if (component < offset) + return 0.0; + return p[component - offset] * p[component - offset]; + } + + unsigned int offset; +}; + +FERemoteEvaluationCommunicator<dim> +construct_comm_for_face_batches(const MatrixFree<dim, Number> &matrix_free) +{ + // Setup Communication objects for all boundary faces + FERemoteCommunicationObjectEntityBatches<dim> co; + + // Get range of boundary face indices + const auto face_batch_range = + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + + // Fill a quadrature vector to keep track of the quadrature sizes on each + // face + std::vector<Quadrature<dim>> quadrature_vector( + matrix_free.n_boundary_face_batches()); + + // Points that are searched by rpe. + std::vector<Point<dim>> points; + + FEFaceEvaluation<dim, -1, 0, 1, Number> phi(matrix_free, true, 0, 0, 0); + for (unsigned int bface = 0; + bface < face_batch_range.second - face_batch_range.first; + ++bface) + { + const unsigned int face = face_batch_range.first + bface; + + phi.reinit(face); + + const unsigned int n_faces = + matrix_free.n_active_entries_per_face_batch(face); + + co.batch_id_n_entities.push_back(std::make_pair(face, n_faces)); + + // Append the quadrature points to the points we need to search + // for. + for (unsigned int v = 0; v < n_faces; ++v) + { + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Point<dim> temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + + points.push_back(temp); + } + } + + // append quadrature of correct size + quadrature_vector[bface] = Quadrature<dim>(phi.n_q_points); + } + + // use rpe to search for stored points + auto rpe = + std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>(1.0e-9); + + rpe->reinit(points, + matrix_free.get_dof_handler().get_triangulation(), + *matrix_free.get_mapping_info().mapping); + Assert(rpe->all_points_found(), ExcMessage("Not all remote points found.")); + Assert(rpe->is_map_unique(), ExcMessage("The map should be unique.")); + + co.rpe = rpe; + + // Renit the communicator `FERemoteEvaluationCommunicator` + // with the communication objects. + FERemoteEvaluationCommunicator<dim> remote_communicator; + remote_communicator.reinit_faces({co}, face_batch_range, quadrature_vector); + + return remote_communicator; +} + +void +test_face_batches(const MatrixFree<dim, Number> &matrix_free, + unsigned int max_first_selected_comp) +{ + // 0) Get remote communicator. + auto remote_communicator = construct_comm_for_face_batches(matrix_free); + + for (unsigned int first_selected_comp = 0; + first_selected_comp < max_first_selected_comp; + ++first_selected_comp) + { + // 1) Allocate vectors and interpolate grid coordinates in DoFs. + VectorType src; + VectorType dst; + matrix_free.initialize_dof_vector(src, first_selected_comp); + matrix_free.initialize_dof_vector(dst, first_selected_comp); + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(first_selected_comp), + QuadraticFunction<dim>(first_selected_comp), + src); + + // 2) Construct remote evaluator... + FERemoteEvaluation<dim, dim, VectorizedArray<Number>> phi_remote_eval( + remote_communicator, + matrix_free.get_dof_handler(first_selected_comp), + first_selected_comp); + + // ...and precompute remote values. + phi_remote_eval.gather_evaluate(src, EvaluationFlags::gradients); + + // 3) Access values at quadrature points + const auto boundary_function = [&](const auto &data, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation<dim, -1, 0, dim, Number> phi( + matrix_free, true, first_selected_comp, 0, first_selected_comp); + auto phi_remote = phi_remote_eval.get_data_accessor(); + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + phi.gather_evaluate(src, EvaluationFlags::gradients); + phi_remote.reinit(face); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto error = + phi.get_gradient(q) - phi_remote.get_gradient(q); + + // only consider active entries in global error + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + AssertThrow(std::abs(error[i][d][v]) < 1e-13, + ExcMessage("Error too large.")); + } + } + }; + + matrix_free.template loop<VectorType, VectorType>( + {}, {}, boundary_function, dst, src, true); + } +} + + +FERemoteEvaluationCommunicator<dim> +construct_comm_for_cell_face_nos(const MatrixFree<dim, Number> &matrix_free) +{ + // Setup Communication objects for all boundary faces + FERemoteCommunicationObjectFaces<dim> co; + + // Get range of boundary face indices + const auto face_batch_range = + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches()); + + // Fill a quadrature vector to keep track of the quadrature sizes on each + // cell, face pair + std::vector<std::vector<Quadrature<dim - 1>>> quadrature_vector; + for (const auto &cell : matrix_free.get_dof_handler() + .get_triangulation() + .active_cell_iterators()) + quadrature_vector.emplace_back( + std::vector<Quadrature<dim - 1>>(cell->n_faces())); + + // Points that are searched by rpe. + std::vector<Point<dim>> points; + + FEFaceEvaluation<dim, -1, 0, 1, Number> phi(matrix_free, true, 0, 0, 0); + for (unsigned int bface = 0; + bface < face_batch_range.second - face_batch_range.first; + ++bface) + { + const unsigned int face = face_batch_range.first + bface; + phi.reinit(face); + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto [cell, f] = matrix_free.get_face_iterator(face, v, true); + co.cell_face_nos.push_back(std::make_pair(cell, f)); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto point = phi.quadrature_point(q); + Point<dim> temp; + for (unsigned int i = 0; i < dim; ++i) + temp[i] = point[i][v]; + points.push_back(temp); + } + + // append quadrature of correct size + quadrature_vector[cell->active_cell_index()][f] = + Quadrature<dim - 1>(phi.n_q_points); + } + } + + // use rpe to search for stored points + auto rpe = + std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>(1.0e-9); + + rpe->reinit(points, + matrix_free.get_dof_handler().get_triangulation(), + *matrix_free.get_mapping_info().mapping); + Assert(rpe->all_points_found(), ExcMessage("Not all remote points found.")); + Assert(rpe->is_map_unique(), ExcMessage("The map should be unique.")); + + co.rpe = rpe; + + // Renit the communicator `FERemoteEvaluationCommunicator` + // with the communication objects. + FERemoteEvaluationCommunicator<dim> remote_communicator; + remote_communicator.reinit_faces( + {co}, + matrix_free.get_dof_handler().get_triangulation().active_cell_iterators(), + quadrature_vector); + + return remote_communicator; +} + +void +test_cell_face_nos(const MatrixFree<dim, Number> &matrix_free, + unsigned int max_first_selected_comp) +{ + // 0) Get remote communicator. + auto remote_communicator = construct_comm_for_cell_face_nos(matrix_free); + + for (unsigned int first_selected_comp = 0; + first_selected_comp < max_first_selected_comp; + ++first_selected_comp) + { + // 1) Allocate vectors and interpolate grid coordinates in DoFs. + VectorType src; + VectorType dst; + matrix_free.initialize_dof_vector(src, first_selected_comp); + matrix_free.initialize_dof_vector(dst, first_selected_comp); + + VectorTools::interpolate(*matrix_free.get_mapping_info().mapping, + matrix_free.get_dof_handler(first_selected_comp), + QuadraticFunction<dim>(first_selected_comp), + src); + + // 2) Construct remote evaluator... + FERemoteEvaluation<dim, dim, Number> phi_remote_eval( + remote_communicator, + matrix_free.get_dof_handler(first_selected_comp), + first_selected_comp); + + // ...and precompute remote values. + phi_remote_eval.gather_evaluate(src, EvaluationFlags::gradients); + + // 3) Access values at quadrature points + const auto boundary_function = [&](const auto &data, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation<dim, -1, 0, dim, Number> phi( + matrix_free, true, first_selected_comp, 0, first_selected_comp); + auto phi_remote = phi_remote_eval.get_data_accessor(); + + + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + phi.gather_evaluate(src, EvaluationFlags::gradients); + + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_face_batch(face); + ++v) + { + const auto [cell, f] = + matrix_free.get_face_iterator(face, v, true); + phi_remote.reinit(cell->active_cell_index(), f); + + for (unsigned int q : phi.quadrature_point_indices()) + { + const auto grad = phi.get_gradient(q); + Tensor<2, dim, Number> temp; + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < dim; ++j) + temp[i][j] = grad[i][j][v]; + + const auto error = temp - phi_remote.get_gradient(q); + for (unsigned int i = 0; i < 2; ++i) + for (unsigned int j = 0; j < dim; ++j) + AssertThrow(std::abs(error[i][j]) < 1e-13, + ExcMessage("Error too large.")); + } + } + } + }; + + matrix_free.template loop<VectorType, VectorType>( + {}, {}, boundary_function, dst, src, true); + } +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + MPILogInitAll all; + + // 1) setup MatrixFree for different finite elements, i.e. + // FESystem(FE_DGQ, dim + first_selected_component). + + const unsigned int refinements = 2; + const unsigned int degree = 5; + + // Construct triangulation. + parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria); + tria.refine_global(refinements); + + // Create DoFHandlers with different number of components to + // check different first_selected components + std::vector<std::unique_ptr<DoFHandler<dim>>> dof_handlers; + unsigned int max_first_selected_comp = 3; + for (unsigned int i = 0; i < max_first_selected_comp; ++i) + { + dof_handlers.push_back(std::make_unique<DoFHandler<dim>>(tria)); + dof_handlers.back()->distribute_dofs( + FESystem<dim>(FE_DGQ<dim>(degree), dim + i)); + } + + // Constraints are not considered in this test + AffineConstraints<Number> constraints; + constraints.close(); + + // Setup MatrixFree + typename MatrixFree<dim, Number>::AdditionalData data; + data.mapping_update_flags = update_gradients; + // We are not considering inner faces in this test. Therefore, + // data.mapping_update_flags_inner_faces is not set. + data.mapping_update_flags_boundary_faces = + update_quadrature_points | update_gradients; + + + // Construct vector of raw pointers needed to init matrix free + std::vector<const DoFHandler<dim> *> dof_handlers_mf; + std::vector<const AffineConstraints<Number> *> constraints_mf; + for (auto &dh : dof_handlers) + { + dof_handlers_mf.push_back(dh.get()); + constraints_mf.push_back(&constraints); + } + + MatrixFree<dim, Number> matrix_free; + matrix_free.reinit(MappingQ1<dim>(), + dof_handlers_mf, + constraints_mf, + QGauss<dim>(degree + 1), + data); + + // 2) test to access remote values at points in face batches + test_face_batches(matrix_free, max_first_selected_comp); + + // 3) test to access remote values at points on faces of cells + test_cell_face_nos(matrix_free, max_first_selected_comp); + + deallog << "OK" << std::endl; + + return 0; +} diff --git a/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=1.output b/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..be8d055f86 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=1.output @@ -0,0 +1,2 @@ + +DEAL:0::OK diff --git a/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=2.output b/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..52103938b0 --- /dev/null +++ b/tests/matrix_free/fe_remote_evaluation_02.with_p4est=true.mpirun=2.output @@ -0,0 +1,5 @@ + +DEAL:0::OK + +DEAL:1::OK + -- 2.39.5