From: Martin Kronbichler Date: Sun, 5 May 2024 06:03:17 +0000 (+0200) Subject: MatrixFree: Switch ShapeInfo from VectorizedArray to Number type X-Git-Tag: v9.6.0-rc1~277^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16985%2Fhead;p=dealii.git MatrixFree: Switch ShapeInfo from VectorizedArray to Number type --- diff --git a/examples/step-76/step-76.cc b/examples/step-76/step-76.cc index f01c40479d..90f67c66e0 100644 --- a/examples/step-76/step-76.cc +++ b/examples/step-76/step-76.cc @@ -663,10 +663,11 @@ namespace Euler_DG dim, n_points_1d, n_points_1d, - VectorizedArrayType> - eval(AlignedVector(), + VectorizedArrayType, + Number> + eval({}, data.get_shape_info().data[0].shape_gradients_collocation_eo, - AlignedVector()); + {}); AlignedVector buffer(phi.static_n_q_points * phi.n_components); diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index a440434ae6..69b77a776e 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -405,7 +405,7 @@ public: * we need to store the evaluations of the 1d polynomials at * the 1d quadrature points. That is what this variable is for. */ - internal::MatrixFreeFunctions::ShapeInfo shape_info; + internal::MatrixFreeFunctions::ShapeInfo shape_info; /** * In case the quadrature rule given represents a tensor product diff --git a/include/deal.II/matrix_free/constraint_info.h b/include/deal.II/matrix_free/constraint_info.h index 6348f96eae..3dde53aabd 100644 --- a/include/deal.II/matrix_free/constraint_info.h +++ b/include/deal.II/matrix_free/constraint_info.h @@ -187,7 +187,7 @@ namespace internal std::vector constraint_pool_data; std::vector constraint_pool_row_index; - std::vector> shape_infos; + std::vector> shape_infos; std::vector hanging_node_constraint_masks; std::vector active_fe_indices; diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 25419d5ded..3ec800f16f 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -1285,7 +1285,7 @@ namespace internal // might have non-symmetric quadrature formula, so use the more // conservative 'evaluate_general' scheme rather than 'even_odd' as the // Hessians are not used very often - const MatrixFreeFunctions::UnivariateShapeData &data = + const MatrixFreeFunctions::UnivariateShapeData &data = fe_eval.get_shape_info().data[0]; AssertDimension(data.shape_gradients_collocation.size(), data.n_q_points_1d * data.n_q_points_1d); @@ -1356,7 +1356,7 @@ namespace internal using Number2 = typename FEEvaluationData::shape_info_number_type; - const MatrixFreeFunctions::UnivariateShapeData &data = + const MatrixFreeFunctions::UnivariateShapeData &data = fe_eval.get_shape_info().data[0]; AssertDimension(data.shape_gradients_collocation.size(), data.n_q_points_1d * data.n_q_points_1d); diff --git a/include/deal.II/matrix_free/evaluation_kernels_face.h b/include/deal.II/matrix_free/evaluation_kernels_face.h index 974d273bbc..f6468be33c 100644 --- a/include/deal.II/matrix_free/evaluation_kernels_face.h +++ b/include/deal.II/matrix_free/evaluation_kernels_face.h @@ -2429,10 +2429,10 @@ namespace internal const auto face_no = fe_eval.get_face_no(v); grad_weight[v] = - shape_data.shape_data_on_face[0][fe_degree + - (integrate ? - (2 - (face_no % 2)) : - (1 + (face_no % 2)))][0]; + shape_data + .shape_data_on_face[0][fe_degree + (integrate ? + (2 - (face_no % 2)) : + (1 + (face_no % 2)))]; index_array_hermite[v] = &fe_eval.get_shape_info().face_to_cell_index_hermite(face_no, diff --git a/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h index 15959f9991..bbd49888b1 100644 --- a/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h +++ b/include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h @@ -88,15 +88,18 @@ namespace internal Number> { private: - template + template static void - interpolate(const unsigned int offset, - const unsigned int outer_stride, - const unsigned int given_degree, - const Number mask_weight, - const Number mask_write, - const Number *DEAL_II_RESTRICT weights, - Number *DEAL_II_RESTRICT values) + interpolate(const unsigned int offset, + const unsigned int outer_stride, + const unsigned int given_degree, + const Number mask_weight, + const Number mask_write, + const Number2 *DEAL_II_RESTRICT weights, + Number *DEAL_II_RESTRICT values) { static constexpr unsigned int max_n_points_1D = 40; @@ -178,19 +181,19 @@ namespace internal } public: - template + template static void run_internal( - const unsigned int n_components, - const MatrixFreeFunctions::ShapeInfo &shape_info, + const unsigned int n_components, + const MatrixFreeFunctions::ShapeInfo &shape_info, const std::array &constraint_mask, - Number *values) + Number::size()> &constraint_mask, + Number *values) { const unsigned int given_degree = fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree; - const Number *DEAL_II_RESTRICT weights = + const Number2 *DEAL_II_RESTRICT weights = shape_info.data.front().subface_interpolation_matrices[0].data(); const unsigned int points = given_degree + 1; @@ -1524,14 +1527,14 @@ namespace internal } public: - template + template static void run_internal( - const unsigned int n_desired_components, - const MatrixFreeFunctions::ShapeInfo &shape_info, + const unsigned int n_desired_components, + const MatrixFreeFunctions::ShapeInfo &shape_info, const std::array &constraint_mask, - Number *values) + Number::size()> &constraint_mask, + Number *values) { const unsigned int given_degree = fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree; @@ -1703,14 +1706,14 @@ namespace internal struct FEEvaluationImplHangingNodes { public: - template + template static bool - run(const unsigned int n_desired_components, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const bool transpose, + run(const unsigned int n_desired_components, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const bool transpose, const std::array &c_mask, - Number *values) + Number::size()> &c_mask, + Number *values) { using RunnerType = FEEvaluationImplHangingNodesRunner(), diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index 7f8f7c3603..b50575ecea 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -153,13 +153,13 @@ namespace internal struct FEEvaluationHangingNodesFactory { static void - apply(const unsigned int n_components, - const unsigned int fe_degree, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const bool transpose, + apply(const unsigned int n_components, + const unsigned int fe_degree, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const bool transpose, const std::array &c_mask, - VectorizedArrayType *values); + VectorizedArrayType::size()> &c_mask, + VectorizedArrayType *values); }; } // end of namespace internal diff --git a/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h index bcfcb992ba..ed41d8fab4 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h @@ -31,13 +31,13 @@ namespace internal template void FEEvaluationHangingNodesFactory::apply( - const unsigned int n_components, - const unsigned int fe_degree, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const bool transpose, + const unsigned int n_components, + const unsigned int fe_degree, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const bool transpose, const std::array &c_mask, - VectorizedArrayType *values) + VectorizedArrayType::size()> &c_mask, + VectorizedArrayType *values) { instantiation_helper_degree_run< 1, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 58fc2bae7f..4db72c6a37 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2439,11 +2439,10 @@ inline FEEvaluationBasedata == nullptr, ExcInternalError()); - this->data = - new internal::MatrixFreeFunctions::ShapeInfo( - Quadrature<(is_face ? dim - 1 : dim)>(quadrature), - fe, - fe.component_to_base_index(first_selected_component).first); + this->data = new internal::MatrixFreeFunctions::ShapeInfo( + Quadrature<(is_face ? dim - 1 : dim)>(quadrature), + fe, + fe.component_to_base_index(first_selected_component).first); this->set_data_pointers(scratch_data_array, n_components_); } @@ -2475,8 +2474,7 @@ inline FEEvaluationBasedata = - new internal::MatrixFreeFunctions::ShapeInfo( - *other.data); + new internal::MatrixFreeFunctions::ShapeInfo(*other.data); // Create deep copy of mapped geometry for use in parallel this->mapped_geometry = @@ -2543,8 +2541,7 @@ operator=(const FEEvaluationBasedata = - new internal::MatrixFreeFunctions::ShapeInfo( - *other.data); + new internal::MatrixFreeFunctions::ShapeInfo(*other.data); scratch_data_array = new AlignedVector(); // Create deep copy of mapped geometry for use in parallel diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h index f9dcc3b8d7..9e955ee19b 100644 --- a/include/deal.II/matrix_free/fe_evaluation_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -112,7 +112,8 @@ namespace internal template class FEEvaluationData { - using ShapeInfoType = internal::MatrixFreeFunctions::ShapeInfo; + using ShapeInfoType = internal::MatrixFreeFunctions::ShapeInfo< + typename internal::VectorizedArrayTrait::value_type>; using MappingInfoStorageType = internal::MatrixFreeFunctions:: MappingInfoStorage<(is_face ? dim - 1 : dim), dim, Number>; using DoFInfo = internal::MatrixFreeFunctions::DoFInfo; @@ -120,10 +121,10 @@ class FEEvaluationData public: static constexpr unsigned int dimension = dim; - using NumberType = Number; - using shape_info_number_type = Number; + using NumberType = Number; using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; + using shape_info_number_type = ScalarNumber; static constexpr unsigned int n_lanes = internal::VectorizedArrayTrait::width(); diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 2d11f50e2f..44d38e9c6e 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -1151,13 +1151,13 @@ namespace internal typename VectorizedDouble> void mapping_q_compute_range( - const unsigned int begin_cell, - const unsigned int end_cell, - const std::vector &cell_type, - const std::vector &process_cell, - const UpdateFlags update_flags_cells, - const AlignedVector &plain_quadrature_points, - const ShapeInfo &shape_info, + const unsigned int begin_cell, + const unsigned int end_cell, + const std::vector &cell_type, + const std::vector &process_cell, + const UpdateFlags update_flags_cells, + const AlignedVector &plain_quadrature_points, + const ShapeInfo &shape_info, MappingInfoStorage &my_data) { constexpr unsigned int n_lanes = VectorizedArrayType::size(); @@ -2083,12 +2083,12 @@ namespace internal const unsigned int begin_face, const unsigned int end_face, const std::vector> - &faces, - const std::vector &face_type, - const std::vector &process_face, - const UpdateFlags update_flags_faces, - const AlignedVector &plain_quadrature_points, - const ShapeInfo &shape_info, + &faces, + const std::vector &face_type, + const std::vector &process_face, + const UpdateFlags update_flags_faces, + const AlignedVector &plain_quadrature_points, + const ShapeInfo &shape_info, MappingInfoStorage &my_data) { constexpr unsigned int n_lanes = VectorizedArrayType::size(); @@ -2696,7 +2696,7 @@ namespace internal // functions or the quadrature points; shape info is merely a vehicle to // return us the right interpolation matrices from the cell support // points to the cell and face quadrature points. - std::vector> shape_infos(cell_data.size()); + std::vector> shape_infos(cell_data.size()); { FE_DGQ fe_geometry(mapping_degree); for (unsigned int my_q = 0; my_q < cell_data.size(); ++my_q) diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index df1deeb17d..78410ba48d 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -2147,7 +2147,7 @@ public: /** * Return the unit cell information for given hp-index. */ - const internal::MatrixFreeFunctions::ShapeInfo & + const internal::MatrixFreeFunctions::ShapeInfo & get_shape_info(const unsigned int dof_handler_index_component = 0, const unsigned int quad_index = 0, const unsigned int fe_base_element = 0, @@ -2292,8 +2292,7 @@ private: /** * Contains shape value information on the unit cell. */ - Table<4, internal::MatrixFreeFunctions::ShapeInfo> - shape_info; + Table<4, internal::MatrixFreeFunctions::ShapeInfo> shape_info; /** * Describes how the cells are gone through. With the cell level (first @@ -2856,7 +2855,7 @@ MatrixFree::get_ghost_set( template -inline const internal::MatrixFreeFunctions::ShapeInfo & +inline const internal::MatrixFreeFunctions::ShapeInfo & MatrixFree::get_shape_info( const unsigned int dof_handler_index, const unsigned int index_quad, @@ -3920,7 +3919,17 @@ namespace internal zero_vector_region(const unsigned int range_index, VectorType &vec) const { if (range_index == numbers::invalid_unsigned_int || range_index == 0) - vec = typename VectorType::value_type(); + { + if constexpr (std::is_same_v< + ArrayView, + VectorType>) + { + for (unsigned int i = 0; i < vec.size(); ++i) + vec[i] = typename VectorType::value_type(); + } + else + vec = typename VectorType::value_type(); + } } diff --git a/tests/matrix_free/fe_evaluation_shift.cc b/tests/matrix_free/fe_evaluation_shift.cc index abaddb0fdc..ae41d86066 100644 --- a/tests/matrix_free/fe_evaluation_shift.cc +++ b/tests/matrix_free/fe_evaluation_shift.cc @@ -74,19 +74,18 @@ namespace dealii }; - const internal::MatrixFreeFunctions::ShapeInfo - *shape_info_base; + const internal::MatrixFreeFunctions::ShapeInfo *shape_info_base; void reinit(unsigned int cell_batch_index, - const internal::MatrixFreeFunctions::ShapeInfo - *shape_info = nullptr); + const internal::MatrixFreeFunctions::ShapeInfo *shape_info = + nullptr); void reinit(unsigned int cell_batch_index, unsigned int face_number, - const internal::MatrixFreeFunctions::ShapeInfo - *shape_info = nullptr); + const internal::MatrixFreeFunctions::ShapeInfo *shape_info = + nullptr); }; template :: reinit(const unsigned int cell_batch_index, - const internal::MatrixFreeFunctions::ShapeInfo - *shape_info) + const internal::MatrixFreeFunctions::ShapeInfo *shape_info) { Assert(this->mapped_geometry == nullptr, ExcMessage( @@ -163,8 +161,7 @@ namespace dealii VectorizedArrayType>:: reinit(const unsigned int cell_batch_index, const unsigned int face_number, - const internal::MatrixFreeFunctions::ShapeInfo - *shape_info) + const internal::MatrixFreeFunctions::ShapeInfo *shape_info) { Assert( this->quad_no < @@ -344,8 +341,7 @@ test(const dealii::FE_Poly &fe) matrix_free.initialize_dof_vector(dst); - std::array, - 2 * dim> + std::array, 2 * dim> shape_info_shift; const QGauss<1> quadrature_1D(n_points); @@ -355,11 +351,11 @@ test(const dealii::FE_Poly &fe) const Quadrature<1> quadrature_1D_shift_plus = shift_1d_quadrature(quadrature_1D, +1.); - dealii::internal::MatrixFreeFunctions::ShapeInfo - shape_info_base(quadrature_1D, fe); - dealii::internal::MatrixFreeFunctions::ShapeInfo + dealii::internal::MatrixFreeFunctions::ShapeInfo shape_info_base( + quadrature_1D, fe); + dealii::internal::MatrixFreeFunctions::ShapeInfo shape_info_shift_plus(quadrature_1D_shift_plus, fe); - dealii::internal::MatrixFreeFunctions::ShapeInfo + dealii::internal::MatrixFreeFunctions::ShapeInfo shape_info_shift_minus(quadrature_1D_shift_minus, fe); for (unsigned int f = 0; f < 2 * dim; ++f) diff --git a/tests/matrix_free/hanging_node_kernels_01.cc b/tests/matrix_free/hanging_node_kernels_01.cc index 49a6520476..fc66737187 100644 --- a/tests/matrix_free/hanging_node_kernels_01.cc +++ b/tests/matrix_free/hanging_node_kernels_01.cc @@ -202,10 +202,10 @@ namespace dealii const auto w = transpose ? constraint_weights[i * (fe_degree + 1) + - interp_idx][v] : + interp_idx] : constraint_weights[interp_idx * (fe_degree + 1) + - i][v]; + i]; t += w * values_temp[real_idx][v]; } } @@ -222,11 +222,10 @@ namespace dealii transpose ? constraint_weights[(fe_degree - i) * (fe_degree + 1) + - fe_degree - interp_idx] - [v] : + fe_degree - interp_idx] : constraint_weights[(fe_degree - interp_idx) * (fe_degree + 1) + - fe_degree - i][v]; + fe_degree - i]; t += w * values_temp[real_idx][v]; } } @@ -375,10 +374,10 @@ namespace dealii const auto w = transpose ? constraint_weights[i * (fe_degree + 1) + - interp_idx][v] : + interp_idx] : constraint_weights[interp_idx * (fe_degree + 1) + - i][v]; + i]; t += w * values_temp[real_idx][v]; } } @@ -397,12 +396,11 @@ namespace dealii transpose ? constraint_weights[(fe_degree - i) * (fe_degree + 1) + - fe_degree - interp_idx] - [v] : + fe_degree - interp_idx] : constraint_weights[(fe_degree - interp_idx) * (fe_degree + 1) + - fe_degree - i][v]; + fe_degree - i]; t += w * values_temp[real_idx][v]; } } diff --git a/tests/performance/timing_navier_stokes.cc b/tests/performance/timing_navier_stokes.cc index bf322a51f9..849779097b 100644 --- a/tests/performance/timing_navier_stokes.cc +++ b/tests/performance/timing_navier_stokes.cc @@ -674,17 +674,19 @@ namespace NavierStokes_DG dim, n_points_1d, n_points_1d, - VectorizedArrayType> - eval(AlignedVector(), + VectorizedArrayType, + Number> + eval({}, data.get_shape_info().data[0].shape_gradients_collocation_eo, - AlignedVector()); + {}); internal::EvaluatorTensorProduct< internal::EvaluatorVariant::evaluate_evenodd, dim - 1, n_q_points_1d, n_q_points_1d, - VectorizedArrayType> + VectorizedArrayType, + Number> eval_face({}, data.get_shape_info().data[0].shape_gradients_collocation_eo, {}); @@ -789,7 +791,7 @@ namespace NavierStokes_DG phi_m.reinit(cell, face); - const AlignedVector &shape_data = + const AlignedVector &shape_data = data.get_shape_info().data.front().quadrature_data_on_face[face % 2]; const std::array n_blocks{