From a99365b7aa464b85d26aa94c5a97a0d086fd6738 Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bu=C4=9Frahan=20Tem=C3=BCr?= Date: Mon, 13 Mar 2023 18:09:28 +0100 Subject: [PATCH] Implement InverseCellwiseMassMatrix with tensor var. coeff. based on the `Flexible` implementation with its own new interface (vector of rank-2 tensors), see #14843 - Unify the flexible implementation as well as fill inverse_JxW_values() for all `fe_degree`s. - Add boolean template parameter `dyadic_coefficients` to the flexible implementation. - Flexible implementation no longer takes the shape values, but fe_eval. - Matrix-vector product for dyadic coefficients implemented in vmult. --- .../deal.II/matrix_free/evaluation_kernels.h | 161 +++++++++++++----- .../matrix_free/evaluation_template_factory.h | 12 +- .../evaluation_template_factory.templates.h | 29 ++-- include/deal.II/matrix_free/operators.h | 90 ++++++++-- 4 files changed, 212 insertions(+), 80 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index da52686722..f999a7eaa6 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -5705,24 +5705,42 @@ namespace internal { template static bool - run(const unsigned int n_desired_components, - const AlignedVector &inverse_shape, - const AlignedVector &inverse_coefficients, - const Number * in_array, - Number * out_array, - std::enable_if_t * = nullptr) + run(const unsigned int n_desired_components, + const FEEvaluationData &fe_eval, + const AlignedVector & inverse_coefficients, + const bool dyadic_coefficients, + const Number * in_array, + Number * out_array) { - constexpr unsigned int dofs_per_component = - Utilities::pow(fe_degree + 1, dim); + const unsigned int given_degree = + (fe_degree > -1) ? fe_degree : + fe_eval.get_shape_info().data.front().fe_degree; + + const unsigned int dofs_per_component = + Utilities::pow(given_degree + 1, dim); + Assert(inverse_coefficients.size() > 0 && inverse_coefficients.size() % dofs_per_component == 0, ExcMessage( "Expected diagonal to be a multiple of scalar dof per cells")); - if (inverse_coefficients.size() != dofs_per_component) - AssertDimension(n_desired_components * dofs_per_component, - inverse_coefficients.size()); + + if (!dyadic_coefficients) + { + if (inverse_coefficients.size() != dofs_per_component) + AssertDimension(n_desired_components * dofs_per_component, + inverse_coefficients.size()) + } + else + { + AssertDimension(n_desired_components * n_desired_components * + dofs_per_component, + inverse_coefficients.size()); + } Assert(dim >= 1 || dim <= 3, ExcNotImplemented()); + Assert(fe_eval.get_shape_info().element_type <= + MatrixFreeFunctions::tensor_symmetric_no_collocation, + ExcNotImplemented()); EvaluatorTensorProduct evaluator(AlignedVector(), AlignedVector(), - inverse_shape); + fe_eval.get_shape_info().data.front().inverse_shape_values_eo, + given_degree + 1, + given_degree + 1); + + const Number *in = in_array; + Number * out = out_array; + + const Number *inv_coefficient = inverse_coefficients.data(); const unsigned int shift_coefficient = inverse_coefficients.size() > dofs_per_component ? dofs_per_component : 0; - const Number *inv_coefficient = inverse_coefficients.data(); - for (unsigned int d = 0; d < n_desired_components; ++d) - { - const Number *in = in_array + d * dofs_per_component; - Number * out = out_array + d * dofs_per_component; - // Need to select 'apply' method with hessian slot because values - // assume symmetries that do not exist in the inverse shapes - evaluator.template hessians<0, true, false>(in, out); - if (dim > 1) - evaluator.template hessians<1, true, false>(out, out); - if (dim > 2) - evaluator.template hessians<2, true, false>(out, out); - for (unsigned int q = 0; q < dofs_per_component; ++q) - out[q] *= inv_coefficient[q]; + const auto n_comp_outer = dyadic_coefficients ? 1 : n_desired_components; + const auto n_comp_inner = dyadic_coefficients ? n_desired_components : 1; - if (dim > 2) - evaluator.template hessians<2, false, false>(out, out); - if (dim > 1) - evaluator.template hessians<1, false, false>(out, out); - evaluator.template hessians<0, false, false>(out, out); + for (unsigned int d = 0; d < n_comp_outer; ++d) + { + for (unsigned int di = 0; di < n_comp_inner; ++di) + { + const Number *in_ = in + di * dofs_per_component; + Number * out_ = out + di * dofs_per_component; + evaluator.template hessians<0, true, false>(in_, out_); + if (dim > 1) + evaluator.template hessians<1, true, false>(out_, out_); + if (dim > 2) + evaluator.template hessians<2, true, false>(out_, out_); + } + if (dyadic_coefficients) + { + const auto n_coeff_components = + n_desired_components * n_desired_components; + if (n_desired_components == dim) + { + for (unsigned int q = 0; q < dofs_per_component; ++q) + vmult(&inv_coefficient[q * n_coeff_components], + &in[q], + &out[q], + dofs_per_component); + } + else + { + for (unsigned int q = 0; q < dofs_per_component; ++q) + vmult<-1>(&inv_coefficient[q * n_coeff_components], + &in[q], + &out[q], + dofs_per_component, + n_desired_components); + } + } + else + for (unsigned int q = 0; q < dofs_per_component; ++q) + out[q] *= inv_coefficient[q]; + + for (unsigned int di = 0; di < n_comp_inner; ++di) + { + Number *out_ = out + di * dofs_per_component; + if (dim > 2) + evaluator.template hessians<2, false, false>(out_, out_); + if (dim > 1) + evaluator.template hessians<1, false, false>(out_, out_); + evaluator.template hessians<0, false, false>(out_, out_); + } + in += dofs_per_component; + out += dofs_per_component; inv_coefficient += shift_coefficient; } + return false; } - /** - * Version for degree = -1 - */ - template - static bool - run(const unsigned int, - const AlignedVector &, - const AlignedVector &, - const Number *, - Number *, - std::enable_if_t * = nullptr) + template + static void + vmult(const Number * inverse_coefficients, + const Number * src, + Number * dst, + const unsigned int dofs_per_component, + const unsigned int n_given_components = 0) { - static_assert(fe_degree == -1, "Only usable for degree -1"); - Assert(false, ExcNotImplemented()); - return false; + const unsigned int n_desired_components = + (n_components > -1) ? n_components : n_given_components; + + std::array tmp; + Assert(n_desired_components <= dim + 2, + ExcMessage( + "Number of components larger than dim+2 not supported.")); + + for (unsigned int d = 0; d < n_desired_components; ++d) + tmp[d] = src[d * dofs_per_component]; + + for (unsigned int d1 = 0; d1 < n_desired_components; ++d1) + { + const Number *inv_coeff_row = + &inverse_coefficients[d1 * n_desired_components]; + Number sum = inv_coeff_row[0] * tmp[0]; + for (unsigned int d2 = 1; d2 < n_desired_components; ++d2) + sum += inv_coeff_row[d2] * tmp[d2]; + dst[d1 * dofs_per_component] = sum; + } } }; diff --git a/include/deal.II/matrix_free/evaluation_template_factory.h b/include/deal.II/matrix_free/evaluation_template_factory.h index 64aaae2b4a..616386cb8b 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.h @@ -105,12 +105,12 @@ namespace internal Number * out_array); static void - apply(const unsigned int n_components, - const unsigned int fe_degree, - const AlignedVector &inverse_shape, - const AlignedVector &inverse_coefficients, - const Number * in_array, - Number * out_array); + apply(const unsigned int n_components, + const FEEvaluationData &fe_eval, + const AlignedVector & inverse_coefficients, + const bool dyadic_coefficients, + const Number * in_array, + Number * out_array); static void transform_from_q_points_to_basis( diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h index f02a7383e4..95b93c6f29 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h @@ -112,22 +112,25 @@ namespace internal template void CellwiseInverseMassFactory::apply( - const unsigned int n_components, - const unsigned int fe_degree, - const AlignedVector &inverse_shape, - const AlignedVector &inverse_coefficients, - const Number * in_array, - Number * out_array) + const unsigned int n_components, + const FEEvaluationData &fe_eval, + const AlignedVector & inverse_coefficients, + const bool dyadic_coefficients, + const Number * in_array, + Number * out_array) { + const unsigned int fe_degree = fe_eval.get_shape_info().data[0].fe_degree; instantiation_helper_run< 1, - CellwiseInverseMassMatrixImplFlexible>(fe_degree, - fe_degree + 1, - n_components, - inverse_shape, - inverse_coefficients, - in_array, - out_array); + CellwiseInverseMassMatrixImplFlexible>( + fe_degree, + fe_degree + 1, + n_components, + fe_eval, + inverse_coefficients, + dyadic_coefficients, + in_array, + out_array); } diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index f6bfcfe0d9..480aafb288 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -648,7 +648,8 @@ namespace MatrixFreeOperators apply(const AlignedVector &inverse_coefficient, const unsigned int n_actual_components, const VectorizedArrayType * in_array, - VectorizedArrayType * out_array) const; + VectorizedArrayType * out_array, + const bool dyadic_coefficients = false) const; /** * Applies the inverse @ref GlossMassMatrix "mass matrix" operation on an input array, using the @@ -665,6 +666,24 @@ namespace MatrixFreeOperators apply(const VectorizedArrayType *in_array, VectorizedArrayType * out_array) const; + /** + * This operation applies the inverse @ref GlossMassMatrix "mass matrix" + * operation on an input array with local dyadic-valued coefficients. + * The second-rank tensor at each quadrature point defines a linear operator + * on a vector holding the dof components. It is assumed that the passed + * input and output arrays are of correct size, namely + * FEEvaluation::dofs_per_cell long. + * `inverse_dyadic_coefficients` must be dofs_per_component long, and every + * element must be a second-rank tensor of dimension `n_components`. All + * entries should also contain the inverse JxW values. The `in_array` and + * `out_array` arguments may point to the same memory position. + */ + void + apply(const AlignedVector> + & inverse_dyadic_coefficients, + const VectorizedArrayType *in_array, + VectorizedArrayType * out_array) const; + /** * This operation performs a projection from the data given in quadrature * points to the actual basis underlying this object. This projection can @@ -1033,8 +1052,12 @@ namespace MatrixFreeOperators fill_inverse_JxW_values( AlignedVector &inverse_jxw) const { - constexpr unsigned int dofs_per_component_on_cell = - Utilities::pow(fe_degree + 1, dim); + const unsigned int dofs_per_component_on_cell = + (fe_degree > -1) ? + Utilities::pow(fe_degree + 1, dim) : + Utilities::pow(fe_eval.get_shape_info().data.front().fe_degree + 1, + dim - 1); + Assert(inverse_jxw.size() > 0 && inverse_jxw.size() % dofs_per_component_on_cell == 0, ExcMessage( @@ -1089,29 +1112,64 @@ namespace MatrixFreeOperators apply(const AlignedVector &inverse_coefficients, const unsigned int n_actual_components, const VectorizedArrayType * in_array, - VectorizedArrayType * out_array) const + VectorizedArrayType * out_array, + const bool dyadic_coefficients) const { - const unsigned int given_degree = - fe_eval.get_shape_info().data[0].fe_degree; if (fe_degree > -1) - internal::CellwiseInverseMassMatrixImplFlexible:: - template run( - n_actual_components, - fe_eval.get_shape_info().data.front().inverse_shape_values_eo, - inverse_coefficients, - in_array, - out_array); + internal::CellwiseInverseMassMatrixImplFlexible< + dim, + VectorizedArrayType>::template run(n_actual_components, + fe_eval, + inverse_coefficients, + dyadic_coefficients, + in_array, + out_array); else internal::CellwiseInverseMassFactory::apply( n_actual_components, - given_degree, - fe_eval.get_shape_info().data.front().inverse_shape_values_eo, + fe_eval, inverse_coefficients, + dyadic_coefficients, in_array, out_array); } + template + inline void + CellwiseInverseMassMatrix:: + apply(const AlignedVector> + & inverse_dyadic_coefficients, + const VectorizedArrayType *in_array, + VectorizedArrayType * out_array) const + { + const unsigned int dofs_per_component = inverse_dyadic_coefficients.size(); + constexpr unsigned int n_tensor_components = n_components * n_components; + + AlignedVector inverse_coefficients( + dofs_per_component * n_tensor_components); + + // Flatten the inverse dyadic coefficients into `inverse_coefficients` + { + auto begin = inverse_coefficients.begin(); + for (unsigned int q = 0; q < dofs_per_component; ++q) + { + const auto end = std::next(begin, n_tensor_components); + inverse_dyadic_coefficients[q].unroll(begin, end); + begin = end; + } + } + + apply(inverse_coefficients, n_components, in_array, out_array, true); + } + template