From 0178accb8f6f32fddbf9cc2a1a1b7af0c4ed21e9 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 15 Dec 2021 00:43:29 +0100 Subject: [PATCH] MatrixFree::gather_evalute(): allow mixed numbers --- .../deal.II/matrix_free/evaluation_kernels.h | 30 +- include/deal.II/matrix_free/fe_evaluation.h | 50 ++-- .../evaluation_template_factory.inst.in | 7 +- .../faces_value_optimization_03.cc | 265 ++++++++++++++++++ .../faces_value_optimization_03.output | 5 + 5 files changed, 318 insertions(+), 39 deletions(-) create mode 100644 tests/matrix_free/faces_value_optimization_03.cc create mode 100644 tests/matrix_free/faces_value_optimization_03.output diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index c85113a867..5ec7296c94 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3163,9 +3163,9 @@ namespace internal const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, typename Processor::Number2_ * global_vector_ptr, - const std::vector> *sm_ptr, - const EvaluationData & fe_eval, - typename Processor::VectorizedArrayType_ * temp1) + const std::vector> *sm_ptr, + const EvaluationData & fe_eval, + typename Processor::VectorizedArrayType_ * temp1) { constexpr int dim = Processor::dim_; constexpr int fe_degree = Processor::fe_degree_; @@ -3529,7 +3529,7 @@ namespace internal dof_info .dof_indices_contiguous_sm[dof_access_index] [cell * n_lanes + v]; - vector_ptrs[v] = const_cast( + vector_ptrs[v] = const_cast( sm_ptr->operator[](temp.first).data() + temp.second + index_offset); } @@ -3554,7 +3554,7 @@ namespace internal dof_info .dof_indices_contiguous_sm[dof_access_index] [cells[v]]; - vector_ptrs[v] = const_cast( + vector_ptrs[v] = const_cast( sm_ptr->operator[](temp.first).data() + temp.second + index_offset); } @@ -3681,18 +3681,17 @@ namespace internal - template + template struct FEFaceEvaluationImplGatherEvaluateSelector { + using Number = typename VectorizedArrayType::value_type; + template static bool run(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number2 * src_ptr, - const std::vector> * sm_ptr, + const std::vector> * sm_ptr, FEEvaluationData &fe_eval) { Assert(fe_degree > -1, ExcInternalError()); @@ -3805,7 +3804,7 @@ namespace internal supports( const EvaluationFlags::EvaluationFlags evaluation_flag, const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * vector_ptr, + const Number2 * vector_ptr, MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage) { const unsigned int fe_degree = shape_info.data.front().fe_degree; @@ -3883,7 +3882,7 @@ namespace internal hermite_grad(T0 & temp_1, T0 & temp_2, const T1 &src_ptr_1, - const T2 &src_ptr_2, + const T1 &src_ptr_2, const T2 &grad_weight) { // case 3a) @@ -3903,12 +3902,11 @@ namespace internal - template + template struct FEFaceEvaluationImplIntegrateScatterSelector { + using Number = typename VectorizedArrayType::value_type; + template static bool run(const unsigned int n_components, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 4840493292..61f9fff40e 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -8059,24 +8059,26 @@ FEFaceEvaluationdata->data.front().n_q_points_1d) && internal::FEFaceEvaluationImplGatherEvaluateSelector< dim, - Number, - VectorizedArrayType>::supports(evaluation_flag, - *this->data, - internal::get_beginning( - input_vector), - this->dof_info->index_storage_variants - [this->dof_access_index][this->cell])) + typename VectorType::value_type, + VectorizedArrayType>:: + supports(evaluation_flag, + *this->data, + internal::get_beginning( + input_vector), + this->dof_info->index_storage_variants[this->dof_access_index] + [this->cell])) { if (fe_degree > -1) { internal::FEFaceEvaluationImplGatherEvaluateSelector< dim, - Number, + typename VectorType::value_type, VectorizedArrayType>::template run( n_components, evaluation_flag, - internal::get_beginning(input_vector), + internal::get_beginning( + input_vector), shared_vector_data, *this); } @@ -8084,10 +8086,11 @@ FEFaceEvaluation::evaluate(n_components, evaluation_flag, - internal::get_beginning( + internal::get_beginning< + typename VectorType::value_type>( input_vector), shared_vector_data, *this); @@ -8173,24 +8176,26 @@ FEFaceEvaluationdata->data.front().n_q_points_1d) && internal::FEFaceEvaluationImplGatherEvaluateSelector< dim, - Number, - VectorizedArrayType>::supports(integration_flag, - *this->data, - internal::get_beginning( - destination), - this->dof_info->index_storage_variants - [this->dof_access_index][this->cell])) + typename VectorType::value_type, + VectorizedArrayType>:: + supports(integration_flag, + *this->data, + internal::get_beginning( + destination), + this->dof_info->index_storage_variants[this->dof_access_index] + [this->cell])) { if (fe_degree > -1) { internal::FEFaceEvaluationImplIntegrateScatterSelector< dim, - Number, + typename VectorType::value_type, VectorizedArrayType>::template run( n_components, integration_flag, - internal::get_beginning(destination), + internal::get_beginning( + destination), shared_vector_data, *this); } @@ -8198,10 +8203,11 @@ FEFaceEvaluation::integrate(n_components, integration_flag, - internal::get_beginning( + internal::get_beginning< + typename VectorType::value_type>( destination), shared_vector_data, *this); diff --git a/source/matrix_free/evaluation_template_factory.inst.in b/source/matrix_free/evaluation_template_factory.inst.in index 6e99a7faff..03b1de0a40 100644 --- a/source/matrix_free/evaluation_template_factory.inst.in +++ b/source/matrix_free/evaluation_template_factory.inst.in @@ -25,7 +25,12 @@ for (deal_II_dimension : DIMENSIONS; template struct dealii::internal::FEFaceEvaluationGatherFactory< deal_II_dimension, - deal_II_scalar_vectorized::value_type, + double, + deal_II_scalar_vectorized>; + + template struct dealii::internal::FEFaceEvaluationGatherFactory< + deal_II_dimension, + float, deal_II_scalar_vectorized>; // inverse mass diff --git a/tests/matrix_free/faces_value_optimization_03.cc b/tests/matrix_free/faces_value_optimization_03.cc new file mode 100644 index 0000000000..d7697fa589 --- /dev/null +++ b/tests/matrix_free/faces_value_optimization_03.cc @@ -0,0 +1,265 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2021 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. +// +// --------------------------------------------------------------------- + + + +// like faces_value_optimization but mixed types + +#include + +#include + +#include + +#include + +#include +#include + +#include "../tests.h" + + +template +class MatrixFreeTest +{ +public: + MatrixFreeTest(const MatrixFree &data) + : data(data) + , error(2) + {} + + ~MatrixFreeTest() + { + deallog << "Error between variants: " << error[0] / error[1] << std::endl; + } + + void + check_error(Vector &src) const + { + data.loop(&MatrixFreeTest::local_apply, + &MatrixFreeTest::local_apply_face, + &MatrixFreeTest::local_apply_face, + this, + src, + src); + } + +private: + void + local_apply(const MatrixFree &, + Vector &, + const Vector &, + const std::pair &) const + {} + + void + local_apply_face( + const MatrixFree &data, + Vector &, + const Vector & src_, + const std::pair &face_range) const + { + FEFaceEvaluation ref(data, true); + FEFaceEvaluation check(data, + true); + + Vector src; + src = src_; + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + ref.reinit(face); + check.reinit(face); + + ref.read_dof_values(src); + ref.evaluate(EvaluationFlags::values); + check.gather_evaluate(src, EvaluationFlags::values); + + for (unsigned int q = 0; q < ref.n_q_points; ++q) + { + VectorizedArray diff = + ref.get_value(q) - check.get_value(q); + for (unsigned int v = 0; v < VectorizedArray::size(); ++v) + { + if (std::abs(diff[v]) > 1e-12) + { + deallog << "Error detected on interior face " << face + << ", v=" << v << ", q=" << q << "!" << std::endl; + deallog << "ref: "; + for (unsigned int i = 0; i < ref.dofs_per_cell; ++i) + deallog << ref.get_dof_value(i)[v] << " "; + deallog << std::endl; + deallog << "done: " << check.get_value(q)[v] + << " instead of " << ref.get_value(q)[v] + << std::endl; + + deallog + << data.get_face_info(face).cells_interior[v] << " " + << (int)data.get_face_info(face).interior_face_no << " " + << (int)data.get_face_info(face).face_orientation << " " + << (int)data.get_face_info(face).subface_index + << std::endl; + deallog << std::endl; + } + error[0] += std::abs(diff[v]); + error[1] += std::abs(ref.get_value(q)[v]); + } + } + } + + FEFaceEvaluation refr(data, + false); + FEFaceEvaluation checkr(data, + false); + + for (unsigned int face = face_range.first; + face < std::min(data.n_inner_face_batches(), face_range.second); + face++) + { + refr.reinit(face); + checkr.reinit(face); + + refr.read_dof_values(src); + refr.evaluate(EvaluationFlags::values); + checkr.gather_evaluate(src, EvaluationFlags::values); + + for (unsigned int q = 0; q < ref.n_q_points; ++q) + { + VectorizedArray diff = + (refr.get_value(q) - checkr.get_value(q)); + for (unsigned int v = 0; v < VectorizedArray::size(); ++v) + { + if (std::abs(diff[v]) > 1e-12) + { + deallog << "Error detected on exterior face " << face + << ", v=" << v << ", q=" << q << "!" << std::endl; + deallog << "ref: "; + for (unsigned int i = 0; i < ref.dofs_per_cell; ++i) + deallog << refr.get_dof_value(i)[v] << " "; + deallog << std::endl; + deallog << "done: " << checkr.get_value(q)[v] + << " instead of " << refr.get_value(q)[v] + << std::endl; + + deallog + << data.get_face_info(face).cells_exterior[v] << " " + << (int)data.get_face_info(face).exterior_face_no << " " + << (int)data.get_face_info(face).face_orientation << " " + << (int)data.get_face_info(face).subface_index + << std::endl; + deallog << std::endl; + } + error[0] += std::abs(diff[v]); + error[1] += std::abs(refr.get_value(q)[v]); + } + } + } + } + + const MatrixFree &data; + mutable std::vector error; +}; + + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_ball(tria); + tria.begin_active()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + typename Triangulation::active_cell_iterator cell, endc; + cell = tria.begin_active(); + endc = tria.end(); + for (; cell != endc; ++cell) + if (cell->center().norm() < 0.5) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.begin(tria.n_levels() - 1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.refine_global(1); + cell = tria.begin_active(); + for (unsigned int i = 0; i < 10 - 3 * dim; ++i) + { + cell = tria.begin_active(); + endc = tria.end(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + FE_DGQ fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + using number = double; + using number2 = float; + + MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.tasks_block_size = 3; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + mf_data.reinit(dof, constraints, quad, data); + } + + MatrixFreeTest mf(mf_data); + + Vector in(dof.n_dofs()); + + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + if (constraints.is_constrained(i)) + continue; + const double entry = Testing::rand() / (double)RAND_MAX; + in(i) = entry; + } + + mf.check_error(in); +} + + + +int +main() +{ + initlog(); + + deallog << std::setprecision(14); + + { + deallog.push("2d"); + test<2, 1>(); + test<2, 2>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + test<3, 2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/faces_value_optimization_03.output b/tests/matrix_free/faces_value_optimization_03.output new file mode 100644 index 0000000000..0bb296849a --- /dev/null +++ b/tests/matrix_free/faces_value_optimization_03.output @@ -0,0 +1,5 @@ + +DEAL:2d::Error between variants: 0 +DEAL:2d::Error between variants: 0 +DEAL:3d::Error between variants: 0 +DEAL:3d::Error between variants: 0 -- 2.39.5