From: Martin Kronbichler Date: Tue, 23 Jul 2024 09:52:57 +0000 (+0200) Subject: FETools: Remove tasking from compute_embedding_matrices X-Git-Tag: v9.6.0-rc1~16^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F17358%2Fhead;p=dealii.git FETools: Remove tasking from compute_embedding_matrices --- diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index e7636e5370..0340f47bea 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -21,7 +21,6 @@ #include #include #include -#include #include #include @@ -1642,91 +1641,47 @@ namespace FETools - namespace internal + template + void + compute_embedding_matrices( + const FiniteElement &fe, + std::vector>> &matrices, + const bool isotropic_only, + const double threshold) { - namespace FEToolsComputeEmbeddingMatricesHelper - { - template - void - compute_embedding_for_shape_function( - const unsigned int i, - const FiniteElement &fe, - const FEValues &coarse, - const Householder &H, - FullMatrix &this_matrix, - const double threshold) - { - const unsigned int n = fe.n_dofs_per_cell(); - const unsigned int nd = fe.n_components(); - const unsigned int nq = coarse.n_quadrature_points; - - Vector v_coarse(nq * nd); - Vector v_fine(n); - - // The right hand side of - // the least squares - // problem consists of the - // function values of the - // coarse grid function in - // each quadrature point. - if (fe.is_primitive()) - { - const unsigned int d = fe.system_to_component_index(i).first; - const double *phi_i = &coarse.shape_value(i, 0); - - for (unsigned int k = 0; k < nq; ++k) - v_coarse(k * nd + d) = phi_i[k]; - } + // loop over all possible refinement cases + unsigned int ref_case_start, ref_case_end; - else - for (unsigned int d = 0; d < nd; ++d) - for (unsigned int k = 0; k < nq; ++k) - v_coarse(k * nd + d) = coarse.shape_value_component(i, k, d); - - // solve the least squares - // problem. - const double result = H.least_squares(v_fine, v_coarse); - Assert(result <= threshold, FETools::ExcLeastSquaresError(result)); - // Avoid warnings in release mode - (void)result; - (void)threshold; - - // Copy into the result - // matrix. Since the matrix - // maps a coarse grid - // function to a fine grid - // function, the columns - // are fine grid. - for (unsigned int j = 0; j < n; ++j) - this_matrix(j, i) = v_fine(j); + if (fe.reference_cell() == ReferenceCells::Tetrahedron) + { + ref_case_start = + static_cast(IsotropicRefinementChoice::cut_tet_68); + ref_case_end = + static_cast(IsotropicRefinementChoice::cut_tet_49); } - - - - template - void - compute_embedding_matrices_for_refinement_case( - const FiniteElement &fe, - std::vector> &matrices, - const unsigned int ref_case, - const double threshold) + else + { + ref_case_start = isotropic_only ? + RefinementCase::isotropic_refinement : + RefinementCase::cut_x; + ref_case_end = RefinementCase::isotropic_refinement; + } + for (unsigned int ref_case = ref_case_start; ref_case <= ref_case_end; + ++ref_case) { const unsigned int n = fe.n_dofs_per_cell(); const ReferenceCell reference_cell = fe.reference_cell(); const unsigned int nc = reference_cell.n_children(RefinementCase(ref_case)); - AssertDimension(matrices.size(), nc); + AssertDimension(matrices[ref_case - 1].size(), nc); for (unsigned int i = 0; i < nc; ++i) { - Assert(matrices[i].n() == n, - ExcDimensionMismatch(matrices[i].n(), n)); - Assert(matrices[i].m() == n, - ExcDimensionMismatch(matrices[i].m(), n)); + AssertDimension(matrices[ref_case - 1][i].n(), n); + AssertDimension(matrices[ref_case - 1][i].m(), n); } - // Set up meshes, one with a single // reference cell and refine it once Triangulation tria; @@ -1755,16 +1710,12 @@ namespace FETools update_quadrature_points | update_JxW_values | update_values); - // We search for the polynomial on - // the small cell, being equal to - // the coarse polynomial in all - // quadrature points. + // We search for the polynomial on the small cell, being equal to the + // coarse polynomial in all quadrature points. - // First build the matrix for this - // least squares problem. This - // contains the values of the fine - // cell polynomials in the fine - // cell grid points. + // First build the matrix for this least squares problem. This + // contains the values of the fine cell polynomials in the fine cell + // grid points. // This matrix is the same for all // children. @@ -1779,16 +1730,14 @@ namespace FETools Householder H(A); - Threads::TaskGroup task_group; - for (const auto &fine_cell : tria.active_cell_iterators()) { fine.reinit(fine_cell); - // evaluate on the coarse cell (which - // is the first -- inactive -- cell on - // the lowest level of the - // triangulation we have created) + // evaluate on the coarse cell (which is the first -- inactive -- + // cell on the lowest level of the triangulation we have created) + + // TODO: Convert this to FEPointEvaluation const std::vector> &q_points_fine = fine.get_quadrature_points(); std::vector> q_points_coarse(q_points_fine.size()); @@ -1805,93 +1754,52 @@ namespace FETools coarse.reinit(tria.begin(0)); FullMatrix &this_matrix = - matrices[fine_cell->active_cell_index()]; - - // Compute this once for each - // coarse grid basis function. can - // spawn subtasks if n is - // sufficiently large so that there - // are more than about 5000 - // operations in the inner loop - // (which is basically const * n^2 - // operations). - if (n > 30) - { - for (unsigned int i = 0; i < n; ++i) - { - task_group += Threads::new_task( - &compute_embedding_for_shape_function, - i, - fe, - coarse, - H, - this_matrix, - threshold); - } - task_group.join_all(); - } - else + matrices[ref_case - 1][fine_cell->active_cell_index()]; + + // Compute this once for each coarse grid basis function. + + Vector v_coarse(nq * nd); + Vector v_fine(n); + + for (unsigned int i = 0; i < n; ++i) { - for (unsigned int i = 0; i < n; ++i) + // The right hand side of the least squares problem consists + // of the function values of the coarse grid function in each + // quadrature point. + if (fe.is_primitive()) { - compute_embedding_for_shape_function( - i, fe, coarse, H, this_matrix, threshold); + const unsigned int d = + fe.system_to_component_index(i).first; + const double *phi_i = &coarse.shape_value(i, 0); + + for (unsigned int k = 0; k < nq; ++k) + v_coarse(k * nd + d) = phi_i[k]; } + + else + for (unsigned int d = 0; d < nd; ++d) + for (unsigned int k = 0; k < nq; ++k) + v_coarse(k * nd + d) = + coarse.shape_value_component(i, k, d); + + // solve the least squares problem. + const double result = H.least_squares(v_fine, v_coarse); + Assert(result <= threshold, + FETools::ExcLeastSquaresError(result)); + (void)result; + (void)threshold; + + for (unsigned int j = 0; j < n; ++j) + this_matrix(j, i) = v_fine(j); } - // Remove small entries from - // the matrix + // Remove small entries from the matrix for (unsigned int i = 0; i < this_matrix.m(); ++i) for (unsigned int j = 0; j < this_matrix.n(); ++j) if (std::fabs(this_matrix(i, j)) < 1e-12) this_matrix(i, j) = 0.; } } - } // namespace FEToolsComputeEmbeddingMatricesHelper - } // namespace internal - - - - template - void - compute_embedding_matrices( - const FiniteElement &fe, - std::vector>> &matrices, - const bool isotropic_only, - const double threshold) - { - Threads::TaskGroup task_group; - - // loop over all possible refinement cases - unsigned int ref_case_start, ref_case_end; - - if (fe.reference_cell() == ReferenceCells::Tetrahedron) - { - ref_case_start = - static_cast(IsotropicRefinementChoice::cut_tet_68); - ref_case_end = - static_cast(IsotropicRefinementChoice::cut_tet_49); - } - else - { - ref_case_start = isotropic_only ? - RefinementCase::isotropic_refinement : - RefinementCase::cut_x; - ref_case_end = RefinementCase::isotropic_refinement; - } - for (unsigned int ref_case = ref_case_start; ref_case <= ref_case_end; - ++ref_case) - task_group += Threads::new_task( - &internal::FEToolsComputeEmbeddingMatricesHelper:: - compute_embedding_matrices_for_refinement_case, - fe, - matrices[ref_case - 1], - ref_case, - threshold); - - task_group.join_all(); } @@ -2178,21 +2086,33 @@ namespace FETools mass.gauss_jordan(); } + // finally loop over all possible refinement cases + unsigned int ref_case_start, ref_case_end; + if (reference_cell == ReferenceCells::Tetrahedron) + { + ref_case_start = + static_cast(IsotropicRefinementChoice::cut_tet_68); + ref_case_end = + static_cast(IsotropicRefinementChoice::cut_tet_49); + } + else + { + ref_case_start = + (isotropic_only ? RefinementCase::isotropic_refinement : + RefinementCase::cut_x); + ref_case_end = RefinementCase::isotropic_refinement; + } - const auto compute_one_case = - [&reference_cell, &mapping, &fe, &q_fine, n, nd, nq]( - const unsigned int ref_case, - const FullMatrix &inverse_mass_matrix, - std::vector> &matrices) { + for (unsigned int ref_case = ref_case_start; ref_case <= ref_case_end; + ++ref_case) + { const unsigned int nc = reference_cell.n_children(RefinementCase(ref_case)); for (unsigned int i = 0; i < nc; ++i) { - Assert(matrices[i].n() == n, - ExcDimensionMismatch(matrices[i].n(), n)); - Assert(matrices[i].m() == n, - ExcDimensionMismatch(matrices[i].m(), n)); + AssertDimension(matrices[ref_case - 1][i].n(), n); + AssertDimension(matrices[ref_case - 1][i].m(), n); } // create a respective refinement on the triangulation @@ -2223,7 +2143,8 @@ namespace FETools for (unsigned int cell_number = 0; cell_number < nc; ++cell_number) { - FullMatrix &this_matrix = matrices[cell_number]; + FullMatrix &this_matrix = + matrices[ref_case - 1][cell_number]; // Compute right hand side, which is a fine level basis // function tested with the coarse level functions. @@ -2275,7 +2196,7 @@ namespace FETools } // RHS ready. Solve system and enter row into matrix - inverse_mass_matrix.vmult(v_coarse, v_fine); + mass.vmult(v_coarse, v_fine); for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) this_matrix(i, j) = v_coarse(i); } @@ -2286,33 +2207,7 @@ namespace FETools if (std::fabs(this_matrix(i, j)) < 1e-12) this_matrix(i, j) = 0.; } - }; - - - // finally loop over all possible refinement cases - Threads::TaskGroup<> tasks; - unsigned int ref_case_start, ref_case_end; - if (reference_cell == ReferenceCells::Tetrahedron) - { - ref_case_start = - static_cast(IsotropicRefinementChoice::cut_tet_68); - ref_case_end = - static_cast(IsotropicRefinementChoice::cut_tet_49); - } - else - { - ref_case_start = - (isotropic_only ? RefinementCase::isotropic_refinement : - RefinementCase::cut_x); - ref_case_end = RefinementCase::isotropic_refinement; } - for (unsigned int ref_case = ref_case_start; ref_case <= ref_case_end; - ++ref_case) - tasks += Threads::new_task([&, ref_case]() { - compute_one_case(ref_case, mass, matrices[ref_case - 1]); - }); - - tasks.join_all(); }