From a36c5bc5fe032abff1ded245d3bfb2bd4af1bbb8 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 31 Aug 2020 16:27:05 +0200 Subject: [PATCH] Enabel BlockVector for FEEvaluation::gather_evalute() --- doc/news/changes/minor/20200831Munch | 4 + .../deal.II/matrix_free/evaluation_kernels.h | 108 ++++++- include/deal.II/matrix_free/fe_evaluation.h | 304 ++++++++++++------ tests/matrix_free/matrix_vector_faces_29.cc | 148 +++++++++ ...h_mpi=true.with_p4est=true.mpirun=2.output | 17 + 5 files changed, 484 insertions(+), 97 deletions(-) create mode 100644 doc/news/changes/minor/20200831Munch create mode 100644 tests/matrix_free/matrix_vector_faces_29.cc create mode 100644 tests/matrix_free/matrix_vector_faces_29.with_mpi=true.with_p4est=true.mpirun=2.output diff --git a/doc/news/changes/minor/20200831Munch b/doc/news/changes/minor/20200831Munch new file mode 100644 index 0000000000..bb36682850 --- /dev/null +++ b/doc/news/changes/minor/20200831Munch @@ -0,0 +1,4 @@ +New: The methods FEEvaluation/FEEFacevaluation::gather_evaluate/integrate_scatter() +can now also accept block vectors. +
+(Peter Munch, Magdalena Schreter, Martin Kronbichler, 2020/08/31) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 4451027de4..120f8fe705 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -26,6 +26,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -2041,10 +2042,44 @@ namespace internal data, temp1, values_array, integrate_gradients, face_no); } - template + template ::value || + !std::is_same().begin()), + Number2 *>::value, + VectorType>::type * = nullptr> + static bool + gather_evaluate(const VectorType &, + const MatrixFreeFunctions::ShapeInfo &, + const MatrixFreeFunctions::DoFInfo &, + VectorizedArrayType *, + VectorizedArrayType *, + VectorizedArrayType *, + const bool, + const bool, + const unsigned int, + const unsigned int, + const std::array, + const std::array, + const unsigned int, + const MatrixFreeFunctions::DoFInfo::DoFAccessIndex, + const std::array, + const Table<2, unsigned int> &) + { + return false; + } + + template ::value && + std::is_same().begin()), + Number2 *>::value, + VectorType>::type * = nullptr> static bool gather_evaluate( - const Number2 * src_ptr, + const VectorType & src, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_quad, @@ -2061,9 +2096,13 @@ namespace internal const std::array face_orientations, const Table<2, unsigned int> & orientation_map) { + static_assert((std::is_same().begin()), + Number2 *>::value), + "Types do not match."); + return process_and_io( // false /*=evaluate*/, - src_ptr, + src.begin(), data, dof_info, values_quad, @@ -2165,10 +2204,63 @@ namespace internal }); } - template + template ::value || + !std::is_same().begin()), + Number2 *>::value, + VectorType>::type * = nullptr> + static bool + integrate_scatter( + VectorType &, + const MatrixFreeFunctions::ShapeInfo &data, + const MatrixFreeFunctions::DoFInfo &, + VectorizedArrayType *values_array, + VectorizedArrayType *values_quad, + VectorizedArrayType *gradients_quad, + VectorizedArrayType *scratch_data, + const bool integrate_values, + const bool integrate_gradients, + const unsigned int, + const unsigned int, + const std::array, + const std::array face_nos, + const unsigned int subface_index, + const MatrixFreeFunctions::DoFInfo::DoFAccessIndex, + const std::array face_orientations, + const Table<2, unsigned int> & orientation_map) + { + AssertDimension(face_nos.size(), 1); + AssertDimension(face_orientations.size(), 1); + + // for block vectors simply integrate + integrate(data, + values_array, + values_quad, + gradients_quad, + scratch_data, + integrate_values, + integrate_gradients, + face_nos[0], + subface_index, + face_orientations[0], + orientation_map); + + // default vector access + return false; + } + + template ::value && + std::is_same().begin()), + Number2 *>::value, + VectorType>::type * = nullptr> static bool integrate_scatter( - Number2 * dst_ptr, + VectorType & dst, const MatrixFreeFunctions::ShapeInfo &data, const MatrixFreeFunctions::DoFInfo & dof_info, VectorizedArrayType * values_array, @@ -2186,9 +2278,13 @@ namespace internal const std::array face_orientations, const Table<2, unsigned int> & orientation_map) { + static_assert((std::is_same().begin()), + Number2 *>::value), + "Types do not match."); + return process_and_io( // true /*=integrate*/, - dst_ptr, + dst.begin(), data, dof_info, values_quad, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 560d7954ef..b91ba87259 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -7315,6 +7315,192 @@ FEEvaluation< gather_evaluate(input_vector, flag); } + +namespace internal +{ + /** + * Implementation for vectors that have the begin() methods. + */ + template ::value && + std::is_same().begin()), + Number *>::value, + VectorType>::type * = nullptr> + bool + try_gather_evaluate_inplace( + T phi, + const VectorType & input_vector, + const unsigned int cell, + const unsigned int active_fe_index, + const unsigned int first_selected_component, + const internal::MatrixFreeFunctions::DoFInfo *dof_info, + const EvaluationFlags::EvaluationFlags evaluation_flag) + { + // If the index storage is interleaved and contiguous and the vector storage + // has the correct alignment, we can directly pass the pointer into the + // vector to the evaluate() call, without reading the vector entries into a + // separate data field. This saves some operations. + if (std::is_same::value && + dof_info->index_storage_variants + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous && + reinterpret_cast( + input_vector.begin() + + dof_info->dof_indices_contiguous + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [cell * VectorizedArrayType::size()]) % + sizeof(VectorizedArrayType) == + 0) + { + const VectorizedArrayType *vec_values = + reinterpret_cast( + input_vector.begin() + + dof_info->dof_indices_contiguous + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [cell * VectorizedArrayType::size()] + + dof_info->component_dof_indices_offset[active_fe_index] + [first_selected_component] * + VectorizedArrayType::size()); + + phi->evaluate(vec_values, evaluation_flag); + + return true; + } + + return false; + } + + /** + * Implementation for all other vectors like block vectors. + */ + template ::value || + !std::is_same().begin()), + Number *>::value, + VectorType>::type * = nullptr> + bool + try_gather_evaluate_inplace(T, + const VectorType &, + const unsigned int, + const unsigned int, + const unsigned int, + const internal::MatrixFreeFunctions::DoFInfo *, + const EvaluationFlags::EvaluationFlags) + { + return false; + } + + /** + * Implementation for vectors that have the begin() methods. + */ + template ::value && + std::is_same().begin()), + Number *>::value, + VectorType>::type * = nullptr> + bool + try_integrate_scatter_inplace( + VectorType & destination, + const unsigned int cell, + const unsigned int n_components, + const unsigned int active_fe_index, + const unsigned int first_selected_component, + const internal::MatrixFreeFunctions::DoFInfo *dof_info, + VectorizedArrayType * values_quad, + VectorizedArrayType * gradients_quad, + VectorizedArrayType * scratch_data, + const internal::MatrixFreeFunctions::ShapeInfo *data, + const EvaluationFlags::EvaluationFlags integration_flag) + { + // If the index storage is interleaved and contiguous and the vector storage + // has the correct alignment, we can directly pass the pointer into the + // vector to the integrate() call, without writing temporary results into a + // separate data field that will later be added into the vector. This saves + // some operations. + if (std::is_same::value && + dof_info->index_storage_variants + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved_contiguous && + reinterpret_cast( + destination.begin() + + dof_info->dof_indices_contiguous + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [cell * VectorizedArrayType::size()]) % + sizeof(VectorizedArrayType) == + 0) + { + VectorizedArrayType *vec_values = + reinterpret_cast( + destination.begin() + + dof_info->dof_indices_contiguous + [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] + [cell * VectorizedArrayType::size()] + + dof_info->component_dof_indices_offset[active_fe_index] + [first_selected_component] * + VectorizedArrayType::size()); + SelectEvaluator:: + integrate(n_components, + integration_flag, + *data, + vec_values, + values_quad, + gradients_quad, + scratch_data, + true); + return true; + } + + return false; + } + + /** + * Implementation for all other vectors like block vectors. + */ + template ::value || + !std::is_same().begin()), + Number *>::value, + VectorType>::type * = nullptr> + bool + try_integrate_scatter_inplace( + VectorType &, + const unsigned int, + const unsigned int, + const unsigned int, + const unsigned int, + const internal::MatrixFreeFunctions::DoFInfo *, + const VectorizedArrayType *, + const VectorizedArrayType *, + const VectorizedArrayType *, + const internal::MatrixFreeFunctions::ShapeInfo *, + const EvaluationFlags::EvaluationFlags) + { + return false; + } +} // namespace internal + template ::value && - this->dof_info->index_storage_variants - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell] == internal::MatrixFreeFunctions::DoFInfo:: - IndexStorageVariants::interleaved_contiguous && - reinterpret_cast( - input_vector.begin() + - this->dof_info->dof_indices_contiguous - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell * VectorizedArrayType::size()]) % - sizeof(VectorizedArrayType) == - 0) - { - const VectorizedArrayType *vec_values = - reinterpret_cast( - input_vector.begin() + - this->dof_info->dof_indices_contiguous - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell * VectorizedArrayType::size()] + - this->dof_info - ->component_dof_indices_offset[this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::size()); - - evaluate(vec_values, evaluation_flag); - } - else + if (internal::try_gather_evaluate_inplace( + this, + input_vector, + this->cell, + this->active_fe_index, + this->first_selected_component, + this->dof_info, + evaluation_flag) == false) { this->read_dof_values(input_vector); evaluate(this->begin_dof_values(), evaluation_flag); @@ -7540,44 +7703,22 @@ FEEvaluation::value && - this->dof_info->index_storage_variants - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell] == internal::MatrixFreeFunctions::DoFInfo:: - IndexStorageVariants::interleaved_contiguous && - reinterpret_cast( - destination.begin() + - this->dof_info->dof_indices_contiguous - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell * VectorizedArrayType::size()]) % - sizeof(VectorizedArrayType) == - 0) - { - VectorizedArrayType *vec_values = reinterpret_cast( - destination.begin() + - this->dof_info->dof_indices_contiguous - [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [this->cell * VectorizedArrayType::size()] + - this->dof_info - ->component_dof_indices_offset[this->active_fe_index] - [this->first_selected_component] * - VectorizedArrayType::size()); - SelectEvaluator:: - integrate(n_components, - integration_flag, - *this->data, - vec_values, - this->values_quad, - this->gradients_quad, - this->scratch_data, - true); - } - else + if (internal::try_integrate_scatter_inplace( + destination, + this->cell, + n_components, + this->active_fe_index, + this->first_selected_component, + this->dof_info, + this->values_quad, + this->gradients_quad, + this->scratch_data, + this->data, + integration_flag) == false) { integrate(integration_flag, this->begin_dof_values()); this->distribute_local_to_global(destination); @@ -8065,15 +8206,6 @@ FEFaceEvaluation::value && - (std::is_same().begin()), - double *>::value || - std::is_same().begin()), - float *>::value), - "This function requires a vector type with begin() function " - "evaluating to a pointer to basic number (float,double). " - "Use read_dof_values() followed by evaluate() instead."); - Assert( (evaluation_flag & ~(EvaluationFlags::values | EvaluationFlags::gradients)) == 0, @@ -8095,7 +8227,7 @@ FEFaceEvaluation:: - gather_evaluate(input_vector.begin(), + gather_evaluate(input_vector, *this->data, *this->dof_info, this->begin_values(), @@ -8121,7 +8253,7 @@ FEFaceEvaluation:: - gather_evaluate(input_vector.begin(), + gather_evaluate(input_vector, *this->data, *this->dof_info, this->begin_values(), @@ -8202,16 +8334,6 @@ FEFaceEvaluation::value && - (std::is_same().begin()), - double *>::value || - std::is_same().begin()), - float *>::value), - "This function requires a vector type with begin() function " - "evaluating to a pointer to basic number (float,double). " - "Use integrate() followed by distribute_local_to_global() " - "instead."); - Assert((this->dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell && this->is_interior_face == false) == false, @@ -8223,7 +8345,7 @@ FEFaceEvaluation:: - integrate_scatter(destination.begin(), + integrate_scatter(destination, *this->data, *this->dof_info, this->begin_dof_values(), diff --git a/tests/matrix_free/matrix_vector_faces_29.cc b/tests/matrix_free/matrix_vector_faces_29.cc new file mode 100644 index 0000000000..e1706bbd0a --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_29.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2020 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. +// +// --------------------------------------------------------------------- + + + +// similar to matrix_vector_faces_13 (renumbering of degrees of freedom to +// better vectorized access, gather_evaluate and integrate_scatter) but with +// block vectors + +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_faces_common.h" + + + +template +void +test() +{ + // raise element degree by one to test quadratic and cubic shape functions + // rather than linears and quadratics according to the main function in + // matrix_vector_faces_common.h + + const unsigned int fe_degree = fe_degree_ + 1; + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria, -1, 1); + + for (typename Triangulation::cell_iterator cell = tria.begin(); + cell != tria.end(); + ++cell) + for (const unsigned int f : GeometryInfo::face_indices()) + if (cell->at_boundary(f)) + cell->face(f)->set_all_boundary_ids(f); + std::vector< + GridTools::PeriodicFacePair::cell_iterator>> + periodic_faces; + for (unsigned int d = 0; d < dim; ++d) + GridTools::collect_periodic_faces( + tria, 2 * d, 2 * d + 1, d, periodic_faces); + tria.add_periodicity(periodic_faces); + + tria.refine_global(3 - dim); + + FE_DGQHermite fe(fe_degree); + DoFHandler dof(tria); + AffineConstraints constraints; + constraints.close(); + + for (unsigned int test = 0; test < 2; ++test) + { + tria.refine_global(1); + + dof.distribute_dofs(fe); + + deallog << "Testing " << dof.get_fe().get_name(); + deallog << " on " << dof.n_dofs() << " DoFs"; + deallog << std::endl; + + MappingQGeneric mapping(dof.get_fe().degree + 1); + + LinearAlgebra::distributed::BlockVector in(1); + LinearAlgebra::distributed::BlockVector out(1); + LinearAlgebra::distributed::BlockVector out_dist(1); + + MatrixFree mf_data; + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + data.initialize_mapping = false; + + mf_data.reinit(mapping, dof, constraints, quad, data); + + std::vector renumbering; + mf_data.renumber_dofs(renumbering); + dof.renumber_dofs(renumbering); + + data.initialize_mapping = true; + mf_data.reinit(mapping, dof, constraints, quad, data); + + mf_data.initialize_dof_vector(in.block(0)); + mf_data.initialize_dof_vector(out.block(0)); + mf_data.initialize_dof_vector(out_dist.block(0)); + + // Set random seed for reproducibility + Testing::srand(42); + for (unsigned int i = 0; i < in.block(0).local_size(); ++i) + { + const double entry = Testing::rand() / (double)RAND_MAX; + in.block(0).local_element(i) = entry; + } + + MatrixFreeTest> + mf(mf_data); + mf.vmult(out, in); + + MatrixFreeVariant> + mf2(mf_data); + mf2.vmult(out_dist, in); + + out_dist -= out; + + double diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << " "; + + // test again, now doing matrix-vector product twice + mf2.vmult(out_dist, in); + mf2.vmult(out_dist, in); + out_dist -= out; + diff_norm = out_dist.linfty_norm() / out.linfty_norm(); + deallog << diff_norm << std::endl; + } +} diff --git a/tests/matrix_free/matrix_vector_faces_29.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/matrix_free/matrix_vector_faces_29.with_mpi=true.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..bf136d0816 --- /dev/null +++ b/tests/matrix_free/matrix_vector_faces_29.with_mpi=true.with_p4est=true.mpirun=2.output @@ -0,0 +1,17 @@ + +DEAL:2d::Testing FE_DGQHermite<2>(2) on 144 DoFs +DEAL:2d::Norm of difference: 0 0 +DEAL:2d::Testing FE_DGQHermite<2>(2) on 576 DoFs +DEAL:2d::Norm of difference: 0 0 +DEAL:2d::Testing FE_DGQHermite<2>(3) on 256 DoFs +DEAL:2d::Norm of difference: 0 0 +DEAL:2d::Testing FE_DGQHermite<2>(3) on 1024 DoFs +DEAL:2d::Norm of difference: 0 0 +DEAL:3d::Testing FE_DGQHermite<3>(2) on 216 DoFs +DEAL:3d::Norm of difference: 0 0 +DEAL:3d::Testing FE_DGQHermite<3>(2) on 1728 DoFs +DEAL:3d::Norm of difference: 0 0 +DEAL:3d::Testing FE_DGQHermite<3>(3) on 512 DoFs +DEAL:3d::Norm of difference: 0 0 +DEAL:3d::Testing FE_DGQHermite<3>(3) on 4096 DoFs +DEAL:3d::Norm of difference: 0 0 -- 2.39.5