From: Martin Kronbichler Date: Mon, 15 May 2023 14:32:26 +0000 (+0200) Subject: Make MatrixFree evaluation kernels more flexible w.r.t. number type X-Git-Tag: v9.5.0-rc1~175^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9da569fdc531088d1833811ab9812857c89403f9;p=dealii.git Make MatrixFree evaluation kernels more flexible w.r.t. number type --- diff --git a/examples/step-76/step-76.cc b/examples/step-76/step-76.cc index 79b888b576..c66a6ac44b 100644 --- a/examples/step-76/step-76.cc +++ b/examples/step-76/step-76.cc @@ -918,15 +918,13 @@ namespace Euler_DG internal::EvaluatorQuantity::hessian, dim, degree + 1, - n_points_1d, - VectorizedArrayType, - VectorizedArrayType>::do_backward(dim + 2, - data.get_shape_info() - .data[0] - .inverse_shape_values_eo, - false, - phi.begin_values(), - phi.begin_dof_values()); + n_points_1d>::do_backward(dim + 2, + data.get_shape_info() + .data[0] + .inverse_shape_values_eo, + false, + phi.begin_values(), + phi.begin_dof_values()); // Perform Runge-Kutta update and write results back to global // vectors: diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 7167750c2d..4d6fa24347 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -118,12 +118,15 @@ namespace internal { static const EvaluatorVariant variant = EvaluatorSelector 4)>::variant; + using Number2 = + typename FEEvaluationData::shape_info_number_type; using Eval = EvaluatorTensorProduct; + Number, + Number2>; static void evaluate(const unsigned int n_components, @@ -140,7 +143,7 @@ namespace internal static Eval create_evaluator_tensor_product( - const MatrixFreeFunctions::UnivariateShapeData + const MatrixFreeFunctions::UnivariateShapeData *univariate_shape_data) { if (variant == evaluate_evenodd) @@ -205,10 +208,10 @@ namespace internal const bool add_into_values_array = false); private: - template + template static EvalType create_evaluator_tensor_product( - const MatrixFreeFunctions::UnivariateShapeData &shape_data) + const MatrixFreeFunctions::UnivariateShapeData &shape_data) { return EvalType(shape_data.shape_values, shape_data.shape_gradients, @@ -251,7 +254,7 @@ namespace internal if (evaluation_flag == EvaluationFlags::nothing) return; - std::array *, 3> + std::array *, 3> univariate_shape_data; const auto &shape_data = fe_eval.get_shape_info().data; @@ -522,7 +525,7 @@ namespace internal FEEvaluationData & fe_eval, const bool add_into_values_array) { - std::array *, 3> + std::array *, 3> univariate_shape_data; const auto &shape_data = fe_eval.get_shape_info().data; @@ -843,8 +846,10 @@ namespace internal const auto &shape_data = fe_eval.get_shape_info().data; + using Number2 = + typename FEEvaluationData::shape_info_number_type; using Eval = - EvaluatorTensorProduct; + EvaluatorTensorProduct; if (evaluation_flag & EvaluationFlags::values) { @@ -913,8 +918,10 @@ namespace internal const auto &shape_data = fe_eval.get_shape_info().data; + using Number2 = + typename FEEvaluationData::shape_info_number_type; using Eval = - EvaluatorTensorProduct; + EvaluatorTensorProduct; if (integration_flag & EvaluationFlags::values) { @@ -1034,13 +1041,16 @@ namespace internal { (void)add_into_values_array; + using Number2 = + typename FEEvaluationData::shape_info_number_type; using EvalNormal = EvaluatorTensorProductAnisotropic; + normal_dir, + Number2>; using EvalTangent = EvaluatorTensorProductAnisotropic; + normal_dir, + Number2>; using Eval0 = typename std::conditional::type; using Eval1 = @@ -1056,9 +1067,8 @@ namespace internal using Eval2 = typename std::conditional::type; - const MatrixFreeFunctions::ShapeInfo &shape_info = - fe_eval.get_shape_info(); - Eval0 eval0 = create_evaluator_tensor_product( + const auto &shape_info = fe_eval.get_shape_info(); + Eval0 eval0 = create_evaluator_tensor_product( ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1])); Eval1 eval1 = create_evaluator_tensor_product( ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1])); @@ -1227,13 +1237,16 @@ namespace internal const bool add_into_values_array, std::integral_constant) { + using Number2 = + typename FEEvaluationData::shape_info_number_type; using EvalNormal = EvaluatorTensorProductAnisotropic; + normal_dir, + Number2>; using EvalTangent = EvaluatorTensorProductAnisotropic; + normal_dir, + Number2>; using Eval0 = typename std::conditional::type; using Eval1 = @@ -1249,9 +1263,8 @@ namespace internal using Eval2 = typename std::conditional::type; - const MatrixFreeFunctions::ShapeInfo &shape_info = - fe_eval.get_shape_info(); - Eval0 eval0 = create_evaluator_tensor_product( + const auto &shape_info = fe_eval.get_shape_info(); + Eval0 eval0 = create_evaluator_tensor_product( ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1])); Eval1 eval1 = create_evaluator_tensor_product( ((normal_dir == 1) ? shape_info.data[0] : shape_info.data[1])); @@ -1423,15 +1436,15 @@ namespace internal * derivative through collocation techniques. * * This class allows for dimension-independent application of the operation, - * implemented by template recursion. It has been tested up to 6D. + * implemented by template recursion. It has been tested up to 6d. */ template + typename = bool, + typename = bool> struct FEEvaluationImplBasisChange { static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2, @@ -1459,17 +1472,19 @@ namespace internal * is zero, the size of the second basis can alternatively be passed in as a * run time argument. */ + template #ifndef DEBUG DEAL_II_ALWAYS_INLINE #endif - static void - do_forward( - const unsigned int n_components, - const AlignedVector &transformation_matrix, - const Number * values_in, - Number * values_out, - const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int, - const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int) + static void + do_forward(const unsigned int n_components, + const AlignedVector &transformation_matrix, + const Number * values_in, + Number * values_out, + const unsigned int basis_size_1_variable = + numbers::invalid_unsigned_int, + const unsigned int basis_size_2_variable = + numbers::invalid_unsigned_int) { Assert( basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable, @@ -1494,8 +1509,8 @@ namespace internal Number, Number2> eval_val(transformation_matrix, - AlignedVector(), - AlignedVector(), + {}, + {}, basis_size_1_variable, basis_size_2_variable); const unsigned int np_1 = @@ -1518,23 +1533,19 @@ namespace internal values_out -= Utilities::fixed_power(np_2); if (next_dim < dim) for (unsigned int q = np_1; q != 0; --q) - FEEvaluationImplBasisChange< - variant, - quantity, - next_dim, - basis_size_1, - basis_size_2, - Number, - Number2>::do_forward(1, - transformation_matrix, - values_in + - (q - 1) * - Utilities::fixed_power(np_1), - values_out + - (q - 1) * - Utilities::fixed_power(np_2), - basis_size_1_variable, - basis_size_2_variable); + FEEvaluationImplBasisChange:: + do_forward(1, + transformation_matrix, + values_in + + (q - 1) * Utilities::fixed_power(np_1), + values_out + + (q - 1) * Utilities::fixed_power(np_2), + basis_size_1_variable, + basis_size_2_variable); // the recursion stops if dim==1 or if dim==2 and // basis_size_1==basis_size_2 (the latter is used because the @@ -1583,18 +1594,20 @@ namespace internal * is zero, the size of the second basis can alternatively be passed in as a * run time argument. */ + template #ifndef DEBUG DEAL_II_ALWAYS_INLINE #endif - static void - do_backward( - const unsigned int n_components, - const AlignedVector &transformation_matrix, - const bool add_into_result, - Number * values_in, - Number * values_out, - const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int, - const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int) + static void + do_backward(const unsigned int n_components, + const AlignedVector &transformation_matrix, + const bool add_into_result, + Number * values_in, + Number * values_out, + const unsigned int basis_size_1_variable = + numbers::invalid_unsigned_int, + const unsigned int basis_size_2_variable = + numbers::invalid_unsigned_int) { Assert( basis_size_1 != 0 || basis_size_1_variable <= basis_size_2_variable, @@ -1698,9 +1711,7 @@ namespace internal quantity, next_dim, basis_size_1, - basis_size_2, - Number, - Number2>:: + basis_size_2>:: do_backward(1, transformation_matrix, add_into_result, @@ -1736,6 +1747,7 @@ namespace internal * of the transformation are stored. It may alias with * the values_in array. */ + template static void do_mass(const unsigned int n_components, const AlignedVector &transformation_matrix, @@ -1761,16 +1773,14 @@ namespace internal EvaluatorQuantity::value, next_dim, basis_size_1, - basis_size_2, - Number, - Number2>::do_forward(n_components, - transformation_matrix, - values_in + - (q - 1) * - Utilities::pow(basis_size_1, dim - 1), - my_scratch + - (q - 1) * - Utilities::pow(basis_size_2, dim - 1)); + basis_size_2>::do_forward(n_components, + transformation_matrix, + values_in + + (q - 1) * + Utilities::pow(basis_size_1, dim - 1), + my_scratch + + (q - 1) * + Utilities::pow(basis_size_2, dim - 1)); EvaluatorTensorProduct::do_backward(n_components, - transformation_matrix, - false, - my_scratch + - q * Utilities::pow(basis_size_2, dim - 1), - values_out + - q * Utilities::pow(basis_size_1, dim - 1)); + FEEvaluationImplBasisChange:: + do_backward(n_components, + transformation_matrix, + false, + my_scratch + q * Utilities::pow(basis_size_2, dim - 1), + values_out + q * Utilities::pow(basis_size_1, dim - 1)); } }; @@ -1831,6 +1837,15 @@ namespace internal template struct FEEvaluationImplCollocation { + using Number2 = + typename FEEvaluationData::shape_info_number_type; + using Eval = EvaluatorTensorProduct; + static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, @@ -1838,7 +1853,7 @@ namespace internal FEEvaluationData & fe_eval); static void - do_evaluate(const MatrixFreeFunctions::UnivariateShapeData &shape, + do_evaluate(const MatrixFreeFunctions::UnivariateShapeData &shape, const EvaluationFlags::EvaluationFlags evaluation_flag, const Number * values_dofs, Number * gradients_quad, @@ -1852,7 +1867,7 @@ namespace internal const bool add_into_values_array); static void - do_integrate(const MatrixFreeFunctions::UnivariateShapeData &shape, + do_integrate(const MatrixFreeFunctions::UnivariateShapeData &shape, const EvaluationFlags::EvaluationFlags integration_flag, Number * values_dofs, Number * gradients_quad, @@ -1893,24 +1908,19 @@ namespace internal template inline void FEEvaluationImplCollocation::do_evaluate( - const MatrixFreeFunctions::UnivariateShapeData &shape, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number * values_dofs, - Number * gradients_quad, - Number * hessians_quad) + const MatrixFreeFunctions::UnivariateShapeData &shape, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + Number * gradients_quad, + Number * hessians_quad) { AssertDimension(shape.shape_gradients_collocation_eo.size(), (fe_degree + 2) / 2 * (fe_degree + 1)); constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); - EvaluatorTensorProduct - eval(AlignedVector(), - shape.shape_gradients_collocation_eo, - shape.shape_hessians_collocation_eo); + Eval eval({}, + shape.shape_gradients_collocation_eo, + shape.shape_hessians_collocation_eo); if ((evaluation_flag & (EvaluationFlags::gradients | EvaluationFlags::hessians)) != 0u) { @@ -1992,24 +2002,19 @@ namespace internal template inline void FEEvaluationImplCollocation::do_integrate( - const MatrixFreeFunctions::UnivariateShapeData &shape, - const EvaluationFlags::EvaluationFlags integration_flag, - Number * values_dofs, - Number * gradients_quad, - const Number * hessians_quad, + const MatrixFreeFunctions::UnivariateShapeData &shape, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + Number * gradients_quad, + const Number * hessians_quad, const bool add_into_values_array) { AssertDimension(shape.shape_gradients_collocation_eo.size(), (fe_degree + 2) / 2 * (fe_degree + 1)); - EvaluatorTensorProduct - eval(AlignedVector(), - shape.shape_gradients_collocation_eo, - shape.shape_hessians_collocation_eo); + Eval eval({}, + shape.shape_gradients_collocation_eo, + shape.shape_hessians_collocation_eo); constexpr std::size_t n_points = Utilities::pow(fe_degree + 1, dim); if ((integration_flag & EvaluationFlags::hessians) != 0u) @@ -2151,12 +2156,10 @@ namespace internal EvaluatorQuantity::value, dim, (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1), - n_q_points_1d, - Number, - Number>::do_forward(1, - shape_data.shape_values_eo, - values_dofs + c * n_dofs, - fe_eval.begin_values() + c * n_q_points); + n_q_points_1d>::do_forward(1, + shape_data.shape_values_eo, + values_dofs + c * n_dofs, + fe_eval.begin_values() + c * n_q_points); // apply derivatives in the collocation space if (evaluation_flag & @@ -2217,14 +2220,12 @@ namespace internal EvaluatorQuantity::value, dim, (fe_degree >= n_q_points_1d ? n_q_points_1d : fe_degree + 1), - n_q_points_1d, - Number, - Number>::do_backward(1, - shape_data.shape_values_eo, - add_into_values_array, - fe_eval.begin_values() + c * n_q_points, - values_dofs + - c * Utilities::pow(fe_degree + 1, dim)); + n_q_points_1d>::do_backward(1, + shape_data.shape_values_eo, + add_into_values_array, + fe_eval.begin_values() + c * n_q_points, + values_dofs + + c * Utilities::pow(fe_degree + 1, dim)); } } @@ -2464,18 +2465,22 @@ namespace internal // choice in terms of operation counts (third condition) and if we were // able to initialize the fields in shape_info.templates.h from the // polynomials (fourth condition). + using Number2 = + typename FEEvaluationData::shape_info_number_type; + using Eval = EvaluatorTensorProduct; + Number, + Number2>; static Eval create_evaluator_tensor_product( - const MatrixFreeFunctions::UnivariateShapeData &data, - const unsigned int subface_index, - const unsigned int direction) + const MatrixFreeFunctions::UnivariateShapeData &data, + const unsigned int subface_index, + const unsigned int direction) { if (symmetric_evaluate) return Eval(data.shape_values_eo, @@ -2503,15 +2508,15 @@ namespace internal static void evaluate_in_face( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::UnivariateShapeData &data, - Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const unsigned int subface_index) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const MatrixFreeFunctions::UnivariateShapeData &data, + Number * values_dofs, + Number * values_quad, + Number * gradients_quad, + Number * hessians_quad, + Number * scratch_data, + const unsigned int subface_index) { Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0); Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1); @@ -2570,10 +2575,9 @@ namespace internal dim - 1, n_q_points_1d, n_q_points_1d, - Number> - eval_grad(AlignedVector(), - data.shape_gradients_collocation_eo, - AlignedVector()); + Number, + Number2> + eval_grad({}, data.shape_gradients_collocation_eo, {}); eval_grad.template gradients<0, true, false>( values_quad, gradients_quad); eval_grad.template gradients<1, true, false>( @@ -2704,15 +2708,15 @@ namespace internal static void integrate_in_face( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - const MatrixFreeFunctions::UnivariateShapeData &data, - Number * values_dofs, - Number * values_quad, - Number * gradients_quad, - Number * hessians_quad, - Number * scratch_data, - const unsigned int subface_index) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + const MatrixFreeFunctions::UnivariateShapeData &data, + Number * values_dofs, + Number * values_quad, + Number * gradients_quad, + Number * hessians_quad, + Number * scratch_data, + const unsigned int subface_index) { Eval eval0 = create_evaluator_tensor_product(data, subface_index, 0); Eval eval1 = create_evaluator_tensor_product(data, subface_index, 1); @@ -2774,10 +2778,9 @@ namespace internal dim - 1, n_q_points_1d, n_q_points_1d, - Number> - eval_grad(AlignedVector(), - data.shape_gradients_collocation_eo, - AlignedVector()); + Number, + Number2> + eval_grad({}, data.shape_gradients_collocation_eo, {}); if ((integration_flag & EvaluationFlags::values) != 0u) eval_grad.template gradients<1, false, true>( gradients_quad + n_q_points, values_quad); @@ -2940,17 +2943,21 @@ namespace internal template struct FEFaceEvaluationImplRaviartThomas { + using Number2 = + typename FEEvaluationData::shape_info_number_type; using EvalGeneral = EvaluatorTensorProduct; + Number, + Number2>; + template static EvalType create_evaluator_tensor_product( - const MatrixFreeFunctions::UnivariateShapeData &data, - const unsigned int subface_index, - const unsigned int direction) + const MatrixFreeFunctions::UnivariateShapeData &data, + const unsigned int subface_index, + const unsigned int direction) { if (subface_index >= GeometryInfo::max_children_per_cell) return EvalType(data.shape_values, @@ -3047,14 +3054,16 @@ namespace internal (fe_degree == -1) ? 1 : fe_degree + 1, n_q_points_1d, Number, - normal_dir>; + normal_dir, + Number2>; using EvalTangent = EvaluatorTensorProductAnisotropic; + normal_dir, + Number2>; using TempEval0 = typename std:: conditional::type; @@ -3065,9 +3074,8 @@ namespace internal using Eval1 = typename std:: conditional::type; - const MatrixFreeFunctions::ShapeInfo &shape_info = - fe_eval.get_shape_info(); - Eval0 eval0 = create_evaluator_tensor_product( + const auto &shape_info = fe_eval.get_shape_info(); + Eval0 eval0 = create_evaluator_tensor_product( ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]), subface_index, 0); @@ -3251,14 +3259,16 @@ namespace internal (fe_degree == -1) ? 1 : fe_degree + 1, n_q_points_1d, Number, - normal_dir>; + normal_dir, + Number2>; using EvalTangent = EvaluatorTensorProductAnisotropic; + normal_dir, + Number2>; using TempEval0 = typename std:: conditional::type; @@ -3269,9 +3279,8 @@ namespace internal using Eval1 = typename std:: conditional::type; - const MatrixFreeFunctions::ShapeInfo &shape_info = - fe_eval.get_shape_info(); - Eval0 eval0 = create_evaluator_tensor_product( + const auto &shape_info = fe_eval.get_shape_info(); + Eval0 eval0 = create_evaluator_tensor_product( ((normal_dir == 0) ? shape_info.data[0] : shape_info.data[1]), subface_index, 0); @@ -3473,14 +3482,17 @@ namespace internal template struct FEFaceNormalEvaluationImpl { + using Number2 = + typename FEEvaluationData::shape_info_number_type; + template static void - interpolate(const unsigned int n_components, - const EvaluationFlags::EvaluationFlags flags, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * input, - Number * output, - const unsigned int face_no) + interpolate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags flags, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const Number * input, + Number * output, + const unsigned int face_no) { Assert(static_cast(fe_degree) == shape_info.data.front().fe_degree || @@ -3508,12 +3520,12 @@ namespace internal template static void interpolate_quadrature( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags flags, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number * input, - Number * output, - const unsigned int face_no) + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags flags, + const MatrixFreeFunctions::ShapeInfo &shape_info, + const Number * input, + Number * output, + const unsigned int face_no) { Assert(static_cast(fe_degree + 1) == shape_info.data.front().n_q_points_1d || @@ -3541,7 +3553,7 @@ namespace internal const EvaluationFlags::EvaluationFlags flag, const unsigned int face_no, const unsigned int n_points_1d, - const std::array, 2> &shape_data, + const std::array, 2> &shape_data, const unsigned int dofs_per_component_on_cell, const unsigned int dofs_per_component_on_face) { @@ -3551,7 +3563,8 @@ namespace internal dim, fe_degree + 1, 0, - Number> + Number, + Number2> evalf(shape_data[face_no % 2].begin(), nullptr, nullptr, @@ -3606,12 +3619,10 @@ namespace internal template static EvalType create_evaluator_tensor_product( - const MatrixFreeFunctions::UnivariateShapeData &data, - const unsigned int face_no) + const MatrixFreeFunctions::UnivariateShapeData &data, + const unsigned int face_no) { - return EvalType(data.shape_data_on_face[face_no % 2], - AlignedVector(), - AlignedVector()); + return EvalType(data.shape_data_on_face[face_no % 2], {}, {}); } template static void interpolate_generic_raviart_thomas( - const unsigned int n_components, - const Number * input, - Number * output, - const EvaluationFlags::EvaluationFlags flag, - const unsigned int face_no, - const MatrixFreeFunctions::ShapeInfo &shape_info) + const unsigned int n_components, + const Number * input, + Number * output, + const EvaluationFlags::EvaluationFlags flag, + const unsigned int face_no, + const MatrixFreeFunctions::ShapeInfo &shape_info) { if (dim == 1) { @@ -3688,10 +3699,10 @@ namespace internal int max_derivative> static inline void interpolate_generic_raviart_thomas_apply_face( - const MatrixFreeFunctions::ShapeInfo &shape_info, - const unsigned int face_no, - const Number * input, - Number * output) + const MatrixFreeFunctions::ShapeInfo &shape_info, + const unsigned int face_no, + const Number * input, + Number * output) { // These types are evaluators in either normal or tangential direction // depending on the face direction, with different normal directions for @@ -3703,13 +3714,15 @@ namespace internal (fe_degree == -1) ? 1 : fe_degree + 1, 0, Number, - 0>, + 0, + Number2>, EvaluatorTensorProductAnisotropic>::type; + 0, + Number2>>::type; using Evalf1 = typename std::conditional< face_direction == 1, EvaluatorTensorProductAnisotropic, + 1, + Number2>, EvaluatorTensorProductAnisotropic>::type; + 1, + Number2>>::type; using Evalf2 = typename std::conditional< face_direction == 2, EvaluatorTensorProductAnisotropic, + 2, + Number2>, EvaluatorTensorProductAnisotropic>::type; + 2, + Number2>>::type; Evalf0 evalf0 = create_evaluator_tensor_product((face_direction == 0) ? @@ -4055,6 +4072,8 @@ namespace internal { const auto &shape_info = fe_eval.get_shape_info(); const auto &shape_data = shape_info.data.front(); + using Number2 = + typename FEEvaluationData::shape_info_number_type; if (shape_info.element_type == MatrixFreeFunctions::tensor_none) { @@ -4069,7 +4088,7 @@ namespace internal const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; using Eval = - EvaluatorTensorProduct; + EvaluatorTensorProduct; if (evaluation_flag & EvaluationFlags::values) { @@ -4095,7 +4114,7 @@ namespace internal auto gradients_quad_ptr = fe_eval.begin_gradients(); auto values_dofs_actual_ptr = values_dofs; - std::array shape_gradients; + std::array shape_gradients; for (unsigned int d = 0; d < dim; ++d) shape_gradients[d] = &shape_data.shape_gradients_face( face_no, face_orientation, d, 0); @@ -4292,6 +4311,8 @@ namespace internal { const auto &shape_info = fe_eval.get_shape_info(); const auto &shape_data = shape_info.data.front(); + using Number2 = + typename FEEvaluationData::shape_info_number_type; if (shape_info.element_type == MatrixFreeFunctions::tensor_none) { @@ -4306,7 +4327,7 @@ namespace internal const std::size_t n_q_points = shape_info.n_q_points_faces[face_no]; using Eval = - EvaluatorTensorProduct; + EvaluatorTensorProduct; if (integration_flag & EvaluationFlags::values) { @@ -4332,7 +4353,7 @@ namespace internal auto gradients_quad_ptr = fe_eval.begin_gradients(); auto values_dofs_actual_ptr = values_dofs; - std::array shape_gradients; + std::array shape_gradients; for (unsigned int d = 0; d < dim; ++d) shape_gradients[d] = &shape_data.shape_gradients_face( face_no, face_orientation, d, 0); @@ -5207,12 +5228,12 @@ namespace internal return false; } + template static bool - supports( - const EvaluationFlags::EvaluationFlags evaluation_flag, - const MatrixFreeFunctions::ShapeInfo &shape_info, - const Number2 * vector_ptr, - MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage) + supports(const EvaluationFlags::EvaluationFlags evaluation_flag, + const MatrixFreeFunctions::ShapeInfo & shape_info, + const Number2 * vector_ptr, + MatrixFreeFunctions::DoFInfo::IndexStorageVariants storage) { const unsigned int fe_degree = shape_info.data.front().fe_degree; if (fe_degree < 1 || !shape_info.data.front().nodal_at_cell_boundaries || @@ -5518,6 +5539,9 @@ namespace internal template struct CellwiseInverseMassMatrixImplBasic { + using Number2 = + typename FEEvaluationData::shape_info_number_type; + template static bool run(const unsigned int n_components, @@ -5538,10 +5562,11 @@ namespace internal dim, fe_degree + 1, fe_degree + 1, - Number> + Number, + Number2> evaluator( - AlignedVector(), - AlignedVector(), + {}, + {}, fe_eval.get_shape_info().data.front().inverse_shape_values_eo); for (unsigned int d = 0; d < n_components; ++d) @@ -5588,12 +5613,12 @@ namespace internal Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); - EvaluatorTensorProduct evaluator( - fe_eval.get_shape_info().data.front().inverse_shape_values, - AlignedVector(), - AlignedVector(), - fe_eval.get_shape_info().data.front().fe_degree + 1, - fe_eval.get_shape_info().data.front().fe_degree + 1); + EvaluatorTensorProduct + evaluator(fe_eval.get_shape_info().data.front().inverse_shape_values, + {}, + {}, + fe_eval.get_shape_info().data.front().fe_degree + 1, + fe_eval.get_shape_info().data.front().fe_degree + 1); for (unsigned int d = 0; d < n_components; ++d) { @@ -5637,6 +5662,9 @@ namespace internal template struct CellwiseInverseMassMatrixImplFlexible { + using Number2 = + typename FEEvaluationData::shape_info_number_type; + template static bool run(const unsigned int n_desired_components, @@ -5680,9 +5708,10 @@ namespace internal dim, fe_degree + 1, fe_degree + 1, - Number> - evaluator(AlignedVector(), - AlignedVector(), + Number, + Number2> + evaluator({}, + {}, fe_eval.get_shape_info().data.front().inverse_shape_values_eo, given_degree + 1, given_degree + 1); @@ -5824,13 +5853,16 @@ namespace internal Utilities::pow(fe_degree + 1, dim) : fe_eval.get_shape_info().n_q_points; + using Number2 = + typename FEEvaluationData::shape_info_number_type; EvaluatorTensorProduct - evaluator(AlignedVector(), - AlignedVector(), + Number, + Number2> + evaluator({}, + {}, inverse_shape, fe_eval.get_shape_info().data.front().fe_degree + 1, fe_eval.get_shape_info().data.front().n_q_points_1d); diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h index aad1d02521..605ee88f4e 100644 --- a/include/deal.II/matrix_free/fe_evaluation_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -121,8 +121,8 @@ class FEEvaluationData public: static constexpr unsigned int dimension = dim; - using NumberType = Number; - + using NumberType = Number; + using shape_info_number_type = Number; using ScalarNumber = typename internal::VectorizedArrayTrait::value_type; @@ -1502,7 +1502,7 @@ FEEvaluationData::get_cell_type() const template -inline const internal::MatrixFreeFunctions::ShapeInfo & +inline const typename FEEvaluationData::ShapeInfoType & FEEvaluationData::get_shape_info() const { Assert(data != nullptr, ExcInternalError()); diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h index 2b6115c832..a46843b837 100644 --- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h +++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h @@ -137,14 +137,12 @@ namespace internal::EvaluatorQuantity::value, dim, degree_coarse + 1, - degree_fine + 1, - Number, - Number>::do_forward(1, - prolongation_matrix_1d, - evaluation_data_coarse, - evaluation_data_fine, - degree_coarse_ + 1, - degree_fine_ + 1); + degree_fine + 1>::do_forward(1, + prolongation_matrix_1d, + evaluation_data_coarse, + evaluation_data_fine, + degree_coarse_ + 1, + degree_fine_ + 1); } void @@ -157,14 +155,12 @@ namespace internal::EvaluatorQuantity::value, 1, 0, - 0, - Number, - Number>::do_forward(1, - prolongation_matrix, - evaluation_data_coarse, - evaluation_data_fine, - n_dofs_coarse, - n_dofs_fine); + 0>::do_forward(1, + prolongation_matrix, + evaluation_data_coarse, + evaluation_data_fine, + n_dofs_coarse, + n_dofs_fine); } private: @@ -203,15 +199,14 @@ namespace internal::EvaluatorQuantity::value, dim, degree_coarse == 0 ? -1 : (degree_coarse + 1), - degree_fine == 0 ? -1 : (degree_fine + 1), - Number, - Number>::do_backward(1, - prolongation_matrix_1d, - false, - evaluation_data_fine, - evaluation_data_coarse, - degree_coarse_ + 1, - degree_fine_ + 1); + degree_fine == 0 ? -1 : (degree_fine + 1)>:: + do_backward(1, + prolongation_matrix_1d, + false, + evaluation_data_fine, + evaluation_data_coarse, + degree_coarse_ + 1, + degree_fine_ + 1); } void @@ -224,15 +219,13 @@ namespace internal::EvaluatorQuantity::value, 1, 0, - 0, - Number, - Number>::do_backward(1, - prolongation_matrix, - false, - evaluation_data_fine, - evaluation_data_coarse, - n_dofs_coarse, - n_dofs_fine); + 0>::do_backward(1, + prolongation_matrix, + false, + evaluation_data_fine, + evaluation_data_coarse, + n_dofs_coarse, + n_dofs_fine); } private: diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc index 1859913180..8b59f142f8 100644 --- a/source/multigrid/mg_transfer_matrix_free.cc +++ b/source/multigrid/mg_transfer_matrix_free.cc @@ -496,17 +496,15 @@ MGTransferMatrixFree::do_prolongate_add( internal::EvaluatorQuantity::value, dim, degree + 1, - 2 * degree + 1, - VectorizedArray, - VectorizedArray>::do_forward(1, - prolongation_matrix_1d, - evaluation_data.begin() + - c * Utilities::fixed_power< - dim>(degree_size), - evaluation_data.begin() + - c * n_scalar_cell_dofs, - fe_degree + 1, - 2 * fe_degree + 1); + 2 * degree + 1>::do_forward(1, + prolongation_matrix_1d, + evaluation_data.begin() + + c * Utilities::fixed_power( + degree_size), + evaluation_data.begin() + + c * n_scalar_cell_dofs, + fe_degree + 1, + 2 * fe_degree + 1); internal::weight_fe_q_dofs_by_entity::do_prolongate_add( internal::EvaluatorQuantity::value, dim, degree + 1, - 2 * degree + 2, - VectorizedArray, - VectorizedArray>::do_forward(1, - prolongation_matrix_1d, - evaluation_data.begin() + - c * Utilities::fixed_power< - dim>(degree_size), - evaluation_data.begin() + - c * n_scalar_cell_dofs, - fe_degree + 1, - 2 * fe_degree + 2); + 2 * degree + 2>::do_forward(1, + prolongation_matrix_1d, + evaluation_data.begin() + + c * Utilities::fixed_power( + degree_size), + evaluation_data.begin() + + c * n_scalar_cell_dofs, + fe_degree + 1, + 2 * fe_degree + 2); } // write into dst vector @@ -606,19 +602,16 @@ MGTransferMatrixFree::do_restrict_add( internal::EvaluatorQuantity::value, dim, degree + 1, - 2 * degree + 1, - VectorizedArray, - VectorizedArray>::do_backward(1, - prolongation_matrix_1d, - false, - evaluation_data.begin() + - c * n_scalar_cell_dofs, - evaluation_data.begin() + - c * - Utilities::fixed_power< - dim>(degree_size), - fe_degree + 1, - 2 * fe_degree + 1); + 2 * degree + 1>::do_backward(1, + prolongation_matrix_1d, + false, + evaluation_data.begin() + + c * n_scalar_cell_dofs, + evaluation_data.begin() + + c * Utilities::fixed_power( + degree_size), + fe_degree + 1, + 2 * fe_degree + 1); } else { @@ -628,19 +621,16 @@ MGTransferMatrixFree::do_restrict_add( internal::EvaluatorQuantity::value, dim, degree + 1, - 2 * degree + 2, - VectorizedArray, - VectorizedArray>::do_backward(1, - prolongation_matrix_1d, - false, - evaluation_data.begin() + - c * n_scalar_cell_dofs, - evaluation_data.begin() + - c * - Utilities::fixed_power< - dim>(degree_size), - fe_degree + 1, - 2 * fe_degree + 2); + 2 * degree + 2>::do_backward(1, + prolongation_matrix_1d, + false, + evaluation_data.begin() + + c * n_scalar_cell_dofs, + evaluation_data.begin() + + c * Utilities::fixed_power( + degree_size), + fe_degree + 1, + 2 * fe_degree + 2); } // write into dst vector diff --git a/tests/performance/timing_navier_stokes.cc b/tests/performance/timing_navier_stokes.cc index 0b2e54544f..4480b17b2c 100644 --- a/tests/performance/timing_navier_stokes.cc +++ b/tests/performance/timing_navier_stokes.cc @@ -1014,15 +1014,13 @@ namespace NavierStokes_DG internal::EvaluatorQuantity::hessian, dim, degree + 1, - n_points_1d, - VectorizedArrayType, - VectorizedArrayType>::do_backward(dim + 2, - data.get_shape_info() - .data[0] - .inverse_shape_values_eo, - false, - phi.begin_values(), - phi.begin_dof_values()); + n_points_1d>::do_backward(dim + 2, + data.get_shape_info() + .data[0] + .inverse_shape_values_eo, + false, + phi.begin_values(), + phi.begin_dof_values()); if (ai == Number()) {