From 56090ff690503220c1cc617eca1a210725681681 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 14 Mar 2017 10:24:42 +0100 Subject: [PATCH] Use simpler code path for quadratic operator evaluation path. --- .../deal.II/matrix_free/evaluation_kernels.h | 91 ++++++++++++++++++- include/deal.II/matrix_free/shape_info.h | 12 +++ .../matrix_free/shape_info.templates.h | 27 ++++++ 3 files changed, 126 insertions(+), 4 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 15c0ac4b08..893bd2ac44 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -198,7 +198,7 @@ namespace internal AssertDimension(count_q, Utilities::fixed_power(shape_info.fe_degree+1)); } - // These avoid compiler errors; they are only used in sensible context but + // These avoid compiler warnings; they are only used in sensible context but // compilers typically cannot detect when we access something like // gradients_quad[2] only for dim==3. const unsigned int d1 = dim>1?1:0; @@ -207,6 +207,46 @@ namespace internal const unsigned int d4 = dim>2?4:0; const unsigned int d5 = dim>2?5:0; + // check if we can go through the spectral evaluation option which is + // faster than the standard one + if (fe_degree+1 == n_q_points_1d && + (type == MatrixFreeFunctions::tensor_symmetric || + type == MatrixFreeFunctions::tensor_general) && + evaluate_lapl == false) + { + Eval eval_grad(shape_info.shape_values, + variant == evaluate_evenodd ? shape_info.shape_grad_spectral_eo : + shape_info.shape_grad_spectral, + shape_info.shape_hessians, + shape_info.fe_degree, + shape_info.n_q_points_1d); + for (unsigned int c=0; c(values_dofs[c], values_quad[c]); + else if (dim == 2) + { + eval.template values<0,true,false>(values_dofs[c], gradients_quad[c][0]); + eval.template values<1,true,false>(gradients_quad[c][0], values_quad[c]); + } + else if (dim == 3) + { + eval.template values<0,true,false>(values_dofs[c], values_quad[c]); + eval.template values<1,true,false>(values_quad[c], gradients_quad[c][0]); + eval.template values<2,true,false>(gradients_quad[c][0], values_quad[c]); + } + if (evaluate_grad == true) + { + eval_grad.template gradients<0,true,false>(values_quad[c], gradients_quad[c][0]); + if (dim >= 2) + eval_grad.template gradients<1,true,false>(values_quad[c], gradients_quad[c][d1]); + if (dim >= 3) + eval_grad.template gradients<2,true,false>(values_quad[c], gradients_quad[c][d2]); + } + } + return; + } + switch (dim) { case 1: @@ -399,12 +439,54 @@ namespace internal c*Utilities::fixed_power(shape_info.fe_degree+1); } - // These avoid compiler errors; they are only used in sensible context but + // These avoid compiler warnings; they are only used in sensible context but // compilers typically cannot detect when we access something like // gradients_quad[2] only for dim==3. const unsigned int d1 = dim>1?1:0; const unsigned int d2 = dim>2?2:0; + // check if we can go through the spectral evaluation option which is + // faster than the standard one + if (fe_degree+1 == n_q_points_1d && + (type == MatrixFreeFunctions::tensor_symmetric || + type == MatrixFreeFunctions::tensor_general)) + { + Eval eval_grad(shape_info.shape_values, + variant == evaluate_evenodd ? shape_info.shape_grad_spectral_eo : + shape_info.shape_grad_spectral, + shape_info.shape_hessians, + shape_info.fe_degree, + shape_info.n_q_points_1d); + for (unsigned int c=0; c(gradients_quad[c][0], values_quad[c]); + else + eval_grad.template gradients<0,false,false>(gradients_quad[c][0], values_quad[c]); + if (dim >= 2) + eval_grad.template gradients<1,false,true>(gradients_quad[c][d1], values_quad[c]); + if (dim >= 3) + eval_grad.template gradients<2,false,true>(gradients_quad[c][d2], values_quad[c]); + } + if (dim == 1) + eval.template values<0,false,false>(values_quad[c], values_dofs[c]); + else if (dim == 2) + { + eval.template values<0,false,false>(values_quad[c], gradients_quad[c][0]); + eval.template values<1,false,false>(gradients_quad[c][0], values_dofs[c]); + } + else if (dim == 3) + { + eval.template values<0,false,false>(values_quad[c], gradients_quad[c][0]); + eval.template values<1,false,false>(gradients_quad[c][0], values_quad[c]); + eval.template values<2,false,false>(values_quad[c], values_dofs[c]); + } + } + return; + } + switch (dim) { case 1: @@ -527,8 +609,9 @@ namespace internal } } - // This a specialization for Gauss-Lobatto elements where the 'values' - // operation is identity, which allows us to write shorter code. + // This a specialization for "spectral" elements like Gauss-Lobatto elements + // where the 'values' operation is identity, which allows us to write + // shorter code. template struct FEEvaluationImpl diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index cc4962f206..2cda877895 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -144,6 +144,18 @@ namespace internal */ AlignedVector > shape_hes_evenodd; + /** + * Stores the shape gradients of the spectral element space + * FE_DGQ<1>(Quadrature<1>) for faster evaluation in even-odd format. + */ + AlignedVector > shape_grad_spectral_eo; + + /** + * Stores the shape gradients of the spectral element space + * FE_DGQ<1>(Quadrature<1>) for faster evaluation in standard format. + */ + AlignedVector > shape_grad_spectral; + /** * Stores the indices from cell DoFs to face DoFs. The rows go through * the 2*dim faces, and the columns the DoFs on the faces. diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 88b302494a..eb6a2e96b6 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -208,6 +209,16 @@ namespace internal this->face_gradient[1][i] = fe->shape_grad(my_i,q_point)[0]; } + // get spectral evaluation points + if (fe_degree+1 == n_q_points_1d) + { + FE_DGQArbitraryNodes<1> fe(quad.get_points()); + shape_grad_spectral.resize(n_dofs_1d*n_dofs_1d); + for (unsigned int i=0; i