From: Peter Munch Date: Wed, 8 May 2024 22:26:04 +0000 (+0200) Subject: Generalize MatrixFreeTools::compute_diagonal() and MatrixFreeTools::compute_matrix() X-Git-Tag: v9.6.0-rc1~200^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=94027f473305df22f7a151fe4b5cab7e77fcfef6;p=dealii.git Generalize MatrixFreeTools::compute_diagonal() and MatrixFreeTools::compute_matrix() --- diff --git a/doc/news/changes/minor/20240516Munch b/doc/news/changes/minor/20240516Munch new file mode 100644 index 0000000000..f141495076 --- /dev/null +++ b/doc/news/changes/minor/20240516Munch @@ -0,0 +1,4 @@ +New: MatrixFreeTools::compute_diagonal() and MatrixFreeTools::compute_matrix() +can now handle face and boundary integrals needed for DG. +
+(Peter Munch, Magdalena Schreter-Fleischhacker, 2024/05/16) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 6dae706111..adfce0cb6d 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -2744,6 +2744,20 @@ namespace internal return vec[component]; } }; + + template + struct BlockVectorSelector, false> + { + using BaseVectorType = VectorType; + + static BaseVectorType * + get_vector_component(std::array &vec, + const unsigned int component) + { + AssertIndexRange(component, vec.size()); + return vec[component]; + } + }; } // namespace internal diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index 73f77159f6..92e3f83166 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -32,6 +32,31 @@ DEAL_II_NAMESPACE_OPEN */ namespace MatrixFreeTools { + namespace internal + { + template + class ComputeMatrixScratchData + { + public: + using FEEvalType = FEEvaluationData; + + std::vector dof_numbers; + std::vector quad_numbers; + std::vector first_selected_components; + std::vector batch_type; + static const bool is_face = is_face_; + + std::function>( + const std::pair &)> + op_create; + std::function> &, + const unsigned int)> + op_reinit; + std::function> &)> + op_compute; + }; + } // namespace internal + /** * Modify @p additional_data so that cells are categorized * according to their boundary IDs, making face integrals in the case of @@ -42,9 +67,11 @@ namespace MatrixFreeTools categorize_by_boundary_ids(const Triangulation &tria, AdditionalData &additional_data); + + /** * Compute the diagonal of a linear operator (@p diagonal_global), given - * @p matrix_free and the local cell integral operation @p local_vmult. The + * @p matrix_free and the local cell integral operation @p cell_operation. The * vector is initialized to the right size in the function. * * The parameters @p dof_no, @p quad_no, and @p first_selected_component are @@ -66,11 +93,14 @@ namespace MatrixFreeTools n_q_points_1d, n_components, Number, - VectorizedArrayType> &)> &local_vmult, - const unsigned int dof_no = 0, - const unsigned int quad_no = 0, + VectorizedArrayType> &)> + &cell_operation, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, const unsigned int first_selected_component = 0); + + /** * Same as above but with a class and a function pointer. */ @@ -98,9 +128,112 @@ namespace MatrixFreeTools const unsigned int first_selected_component = 0); + + /** + * Compute the diagonal of a linear operator (@p diagonal_global), given + * @p matrix_free and the local cell integral operation @p cell_operation, + * interior face integral operation @p face_operation, and boundary face + * integral operation @p boundary_operation. The + * vector is initialized to the right size in the function. + * + * The parameters @p dof_no, @p quad_no, and @p first_selected_component are + * passed to the constructor of the FEEvaluation that is internally set up. + */ + template + void + compute_diagonal( + const MatrixFree &matrix_free, + VectorType &diagonal_global, + const std::function &)> + &cell_operation, + const std::function &, + FEFaceEvaluation &)> + &face_operation, + const std::function &)> + &boundary_operation, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, + const unsigned int first_selected_component = 0); + + + + /** + * Same as above but with a class and a function pointer. + */ + template + void + compute_diagonal( + const MatrixFree &matrix_free, + VectorType &diagonal_global, + void (CLASS::*cell_operation)(FEEvaluation &) const, + void (CLASS::*face_operation)(FEFaceEvaluation &, + FEFaceEvaluation &) + const, + void (CLASS::*boundary_operation)(FEFaceEvaluation &) + const, + const CLASS *owning_class, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, + const unsigned int first_selected_component = 0); + + + /** * Compute the matrix representation of a linear operator (@p matrix), given - * @p matrix_free and the local cell integral operation @p local_vmult. + * @p matrix_free and the local cell integral operation @p cell_operation. * Constrained entries on the diagonal are set to one. * * The parameters @p dof_no, @p quad_no, and @p first_selected_component are @@ -115,19 +248,22 @@ namespace MatrixFreeTools typename MatrixType> void compute_matrix( - const MatrixFree &matrix_free, - const AffineConstraints &constraints, - MatrixType &matrix, + const MatrixFree &matrix_free, + const AffineConstraints &constraints, + MatrixType &matrix, const std::function &)> &local_vmult, - const unsigned int dof_no = 0, - const unsigned int quad_no = 0, + VectorizedArrayType> &)> + &cell_operation, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, const unsigned int first_selected_component = 0); + + /** * Same as above but with a class and a function pointer. */ @@ -156,6 +292,133 @@ namespace MatrixFreeTools const unsigned int first_selected_component = 0); + namespace internal + { + /** + * Compute the matrix representation of a linear operator (@p matrix), given + * @p matrix_free and the local cell, face and boundary integral operation. + */ + template + void + compute_matrix( + const MatrixFree &matrix_free, + const AffineConstraints &constraints, + const internal::ComputeMatrixScratchData + &cell_operation, + const internal::ComputeMatrixScratchData + &face_operation, + const internal::ComputeMatrixScratchData + &boundary_operation, + MatrixType &matrix); + } // namespace internal + + + + /** + * Compute the matrix representation of a linear operator (@p matrix), given + * @p matrix_free and the local cell integral operation @p cell_operation, + * interior face integral operation @p face_operation, and boundary face + * integal operation @p boundary_operation. + * + * The parameters @p dof_no, @p quad_no, and @p first_selected_component are + * passed to the constructor of the FEEvaluation that is internally set up. + */ + template + void + compute_matrix( + const MatrixFree &matrix_free, + const AffineConstraints &constraints, + MatrixType &matrix, + const std::function &)> + &cell_operation, + const std::function &, + FEFaceEvaluation &)> + &face_operation, + const std::function &)> + &boundary_operation, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, + const unsigned int first_selected_component = 0); + + + + /** + * Same as above but with a class and a function pointer. + */ + template + void + compute_matrix( + const MatrixFree &matrix_free, + const AffineConstraints &constraints, + MatrixType &matrix, + void (CLASS::*cell_operation)(FEEvaluation &) const, + void (CLASS::*face_operation)(FEFaceEvaluation &, + FEFaceEvaluation &) + const, + void (CLASS::*boundary_operation)(FEFaceEvaluation &) + const, + const CLASS *owning_class, + const unsigned int dof_no = 0, + const unsigned int quad_no = 0, + const unsigned int first_selected_component = 0); + + /** * A wrapper around MatrixFree to help users to deal with DoFHandler @@ -430,16 +693,16 @@ namespace MatrixFreeTools std::vector> inverse_lookup_origins; }; - template + template class ComputeDiagonalHelper { public: - static const unsigned int n_lanes = VectorizedArrayType::size(); + using Number = typename FEEvaluationType::number_type; + using VectorizedArrayType = typename FEEvaluationType::NumberType; + + static const unsigned int dim = FEEvaluationType::dimension; + static const unsigned int n_components = FEEvaluationType::n_components; + static const unsigned int n_lanes = VectorizedArrayType::size(); ComputeDiagonalHelper() : phi(nullptr) @@ -452,12 +715,7 @@ namespace MatrixFreeTools {} void - initialize(FEEvaluation &phi) + initialize(FEEvaluationType &phi) { // if we are in hp mode and the number of unknowns changed, we must // clear the map of entries @@ -486,8 +744,9 @@ namespace MatrixFreeTools const unsigned int first_selected_component = n_fe_components == 1 ? 0 : phi->get_first_selected_component(); - const unsigned int n_lanes_filled = - matrix_free.n_active_entries_per_cell_batch(cell); + this->n_lanes_filled = + is_face ? matrix_free.n_active_entries_per_face_batch(cell) : + matrix_free.n_active_entries_per_cell_batch(cell); // STEP 2: setup CSR storage of transposed locally-relevant // constraint matrix @@ -500,13 +759,19 @@ namespace MatrixFreeTools std::vector constraint_position; std::vector is_constrained_hn; + const std::array &cells = + this->phi->get_cell_ids(); + for (unsigned int v = 0; v < n_lanes_filled; ++v) { + Assert(cells[v] != numbers::invalid_unsigned_int, + ExcInternalError()); + const unsigned int *dof_indices; unsigned int index_indicators, next_index_indicators; const unsigned int start = - (cell * n_lanes + v) * n_fe_components + first_selected_component; + cells[v] * n_fe_components + first_selected_component; dof_indices = dof_info.dof_indices.data() + dof_info.row_starts[start].first; index_indicators = dof_info.row_starts[start].second; @@ -625,7 +890,7 @@ namespace MatrixFreeTools [phi->get_active_fe_index()][first_selected_component]) { const auto mask = - dof_info.hanging_node_constraint_masks[cell * n_lanes + v]; + dof_info.hanging_node_constraint_masks[cells[v]]; // cell has hanging nodes if (mask != dealii::internal::MatrixFreeFunctions:: @@ -787,6 +1052,42 @@ namespace MatrixFreeTools c_pools[v].row_lid_to_gid.size() * (n_fe_components == 1 ? n_components : 1), Number(0.0)); + + // check if fast path can be taken via FEEvaluation + bool use_fast_path = true; + + for (unsigned int v = 0; v < n_lanes_filled; ++v) + { + auto &c_pool = c_pools[v]; + + for (unsigned int i = 0; i < c_pool.row.size() - 1; ++i) + { + if ((c_pool.row[i + 1] - c_pool.row[i]) > 1) + { + use_fast_path = false; + break; + } + else if (((c_pool.row[i + 1] - c_pool.row[i]) == 1) && + (c_pool.val[c_pool.row[i]] != 1.0)) + { + use_fast_path = false; + break; + } + } + + if (use_fast_path == false) + break; + } + + if (use_fast_path) + { + temp_values.resize(phi->dofs_per_cell); + return; + } + else + { + temp_values.clear(); + } } void @@ -866,9 +1167,23 @@ namespace MatrixFreeTools dof_values[i] = Number(1); } + void + zero_basis_vector() + { + VectorizedArrayType *dof_values = phi->begin_dof_values(); + for (const unsigned int j : phi->dof_indices()) + dof_values[j] = VectorizedArrayType(); + } + void submit() { + if (!temp_values.empty()) + { + temp_values[i] = phi->begin_dof_values()[i]; + return; + } + // if we have a block vector with components with the same DoFHandler, // we need to figure out which component and which DoF within the // component are we currently considering @@ -881,10 +1196,7 @@ namespace MatrixFreeTools // apply local constraint matrix from left and from right: // loop over all rows of transposed constrained matrix - for (unsigned int v = 0; - v < phi->get_matrix_free().n_active_entries_per_cell_batch( - phi->get_current_cell_index()); - ++v) + for (unsigned int v = 0; v < n_lanes_filled; ++v) { const auto &c_pool = c_pools[v]; @@ -913,14 +1225,21 @@ namespace MatrixFreeTools distribute_local_to_global( std::array &diagonal_global) { + if (!temp_values.empty()) + { + for (unsigned int j = 0; j < temp_values.size(); ++j) + phi->begin_dof_values()[j] = temp_values[j]; + + phi->distribute_local_to_global(diagonal_global); + + return; + } + // STEP 4: assembly results: add into global vector const unsigned int n_fe_components = phi->get_dof_info().start_components.back(); - for (unsigned int v = 0; - v < phi->get_matrix_free().n_active_entries_per_cell_batch( - phi->get_current_cell_index()); - ++v) + for (unsigned int v = 0; v < n_lanes_filled; ++v) // if we have a block vector with components with the same // DoFHandler, we need to loop over all components manually and // need to apply the correct shift @@ -937,18 +1256,21 @@ namespace MatrixFreeTools [v][j + comp * c_pools[v].row_lid_to_gid.size()]); } + bool + use_fast_path() const + { + return !temp_values.empty(); + } + private: - FEEvaluation *phi; + FEEvaluationType *phi; unsigned int dofs_per_component; unsigned int i; + unsigned int n_lanes_filled; + std::array, n_lanes> c_pools; // local storage: buffer so that we access the global vector once @@ -962,9 +1284,52 @@ namespace MatrixFreeTools locally_relevant_constraints_hn_map; // scratch array + AlignedVector temp_values; AlignedVector values_dofs; }; + template + bool + is_fe_nothing( + const MatrixFree &matrix_free, + const std::pair &range, + const unsigned int dof_no, + const unsigned int quad_no, + const unsigned int first_selected_component, + const unsigned int fe_degree, + const unsigned int n_q_points_1d, + const bool is_interior_face = true) + { + const unsigned int static_n_q_points = + is_face ? Utilities::pow(n_q_points_1d, dim - 1) : + Utilities::pow(n_q_points_1d, dim); + + unsigned int active_fe_index = 0; + if (!is_face) + active_fe_index = matrix_free.get_cell_active_fe_index(range); + else if (is_interior_face) + active_fe_index = matrix_free.get_face_range_category(range).first; + else + active_fe_index = matrix_free.get_face_range_category(range).second; + + const auto init_data = dealii::internal:: + extract_initialization_data( + matrix_free, + dof_no, + first_selected_component, + quad_no, + fe_degree, + static_n_q_points, + active_fe_index, + numbers::invalid_unsigned_int /*active_quad_index*/, + numbers::invalid_unsigned_int /*face_type*/); + + return init_data.shape_info->dofs_per_component_on_cell == 0; + } + } // namespace internal template &)> &local_vmult, - const unsigned int dof_no, - const unsigned int quad_no, + VectorizedArrayType> &)> + &cell_operation, + const unsigned int dof_no, + const unsigned int quad_no, + const unsigned int first_selected_component) + { + compute_diagonal(matrix_free, + diagonal_global, + cell_operation, + {}, + {}, + dof_no, + quad_no, + first_selected_component); + } + + template + void + compute_diagonal( + const MatrixFree &matrix_free, + VectorType &diagonal_global, + void (CLASS::*cell_operation)(FEEvaluation &) const, + const CLASS *owning_class, + const unsigned int dof_no, + const unsigned int quad_no, + const unsigned int first_selected_component) + { + compute_diagonal( + matrix_free, + diagonal_global, + [&](auto &phi) { (owning_class->*cell_operation)(phi); }, + dof_no, + quad_no, + first_selected_component); + } + + template + void + compute_diagonal( + const MatrixFree &matrix_free, + VectorType &diagonal_global, + const std::function &)> + &cell_operation, + const std::function &, + FEFaceEvaluation &)> + &face_operation, + const std::function &)> + &boundary_operation, + const unsigned int dof_no, + const unsigned int quad_no, const unsigned int first_selected_component) { int dummy = 0; @@ -1023,19 +1485,46 @@ namespace MatrixFreeTools *diagonal_global_components[0], matrix_free, dof_info); } - using Helper = internal::ComputeDiagonalHelper; - - Threads::ThreadLocalStorage scratch_data; - matrix_free.template cell_loop( + VectorizedArrayType>, + false>; + + using HelperFace = + internal::ComputeDiagonalHelper, + true>; + + Threads::ThreadLocalStorage scratch_data; + Threads::ThreadLocalStorage scratch_data_m; + Threads::ThreadLocalStorage scratch_data_p; + + const auto cell_operation_wrapped = [&](const MatrixFree &matrix_free, VectorType &, const int &, - const std::pair &range) mutable { + const std::pair &range) { + if (!cell_operation) + return; + + // shortcut for FE_Nothing cells + if (internal::is_fe_nothing(matrix_free, + range, + dof_no, + quad_no, + first_selected_component, + fe_degree, + n_q_points_1d)) + return; + Helper &helper = scratch_data.get(); FEEvaluation &matrix_free, + VectorType &, + const int &, + const std::pair &range) { + if (!face_operation) + return; + + // shortcut for faces containing purely FE_Nothing + if (internal::is_fe_nothing(matrix_free, + range, + dof_no, + quad_no, + first_selected_component, + fe_degree, + n_q_points_1d, + true) || + internal::is_fe_nothing(matrix_free, + range, + dof_no, + quad_no, + first_selected_component, + fe_degree, + n_q_points_1d, + false)) + return; + + HelperFace &helper_m = scratch_data_m.get(); + HelperFace &helper_p = scratch_data_p.get(); + + FEFaceEvaluation + phi_m(matrix_free, + range, + true, + dof_no, + quad_no, + first_selected_component); + FEFaceEvaluation + phi_p(matrix_free, + range, + false, + dof_no, + quad_no, + first_selected_component); + + helper_m.initialize(phi_m); + helper_p.initialize(phi_p); + + for (unsigned int face = range.first; face < range.second; ++face) + { + helper_m.reinit(face); + helper_p.reinit(face); + + // make check only if both adjacent cells have DoFs + Assert(helper_m.use_fast_path() && helper_p.use_fast_path(), + ExcNotImplemented()); + + for (unsigned int i = 0; i < phi_m.dofs_per_cell; ++i) + { + helper_m.prepare_basis_vector(i); + helper_p.zero_basis_vector(); + face_operation(phi_m, phi_p); + helper_m.submit(); + } + + helper_m.distribute_local_to_global(diagonal_global_components); + + for (unsigned int i = 0; i < phi_p.dofs_per_cell; ++i) + { + helper_m.zero_basis_vector(); + helper_p.prepare_basis_vector(i); + face_operation(phi_m, phi_p); + helper_p.submit(); + } + + helper_p.distribute_local_to_global(diagonal_global_components); + } + }; + + const auto boundary_operation_wrapped = + [&](const MatrixFree &matrix_free, + VectorType &, + const int &, + const std::pair &range) { + if (!boundary_operation) + return; + + // shortcut for faces containing purely FE_Nothing + if (internal::is_fe_nothing(matrix_free, + range, + dof_no, + quad_no, + first_selected_component, + fe_degree, + n_q_points_1d, + true)) + return; + + HelperFace &helper = scratch_data_m.get(); + + FEFaceEvaluation + phi(matrix_free, + range, + true, + dof_no, + quad_no, + first_selected_component); + helper.initialize(phi); + + for (unsigned int face = range.first; face < range.second; ++face) + { + helper.reinit(face); + + for (unsigned int i = 0; i < phi.dofs_per_cell; ++i) + { + helper.prepare_basis_vector(i); + boundary_operation(phi); + helper.submit(); + } + + helper.distribute_local_to_global(diagonal_global_components); + } + }; + + if (face_operation || boundary_operation) + matrix_free.template loop(cell_operation_wrapped, + face_operation_wrapped, + boundary_operation_wrapped, + diagonal_global, + dummy, + false); + else + matrix_free.template cell_loop(cell_operation_wrapped, + diagonal_global, + dummy, + false); } template &) const, + void (CLASS::*face_operation)(FEFaceEvaluation &, + FEFaceEvaluation &) + const, + void (CLASS::*boundary_operation)(FEFaceEvaluation &) + const, const CLASS *owning_class, const unsigned int dof_no, const unsigned int quad_no, @@ -1098,7 +1756,11 @@ namespace MatrixFreeTools VectorType>( matrix_free, diagonal_global, - [&](auto &feeval) { (owning_class->*cell_operation)(feeval); }, + [&](auto &phi) { (owning_class->*cell_operation)(phi); }, + [&](auto &phi_m, auto &phi_p) { + (owning_class->*face_operation)(phi_m, phi_p); + }, + [&](auto &phi) { (owning_class->*boundary_operation)(phi); }, dof_no, quad_no, first_selected_component); @@ -1170,99 +1832,435 @@ namespace MatrixFreeTools n_q_points_1d, n_components, Number, - VectorizedArrayType> &)> &local_vmult, - const unsigned int dof_no, - const unsigned int quad_no, + VectorizedArrayType> &)> + &cell_operation, + const unsigned int dof_no, + const unsigned int quad_no, const unsigned int first_selected_component) { - std::unique_ptr> - constraints_for_matrix; - const AffineConstraints &constraints = - internal::create_new_affine_constraints_if_needed(matrix, - constraints_in, - constraints_for_matrix); - - matrix_free.template cell_loop( - [&](const auto &, auto &dst, const auto &, const auto range) { - FEEvaluation - integrator( - matrix_free, range, dof_no, quad_no, first_selected_component); - - const unsigned int dofs_per_cell = integrator.dofs_per_cell; + compute_matrix(matrix_free, + constraints_in, + matrix, + cell_operation, + {}, + {}, + dof_no, + quad_no, + first_selected_component); + } - std::vector dof_indices(dofs_per_cell); - std::vector dof_indices_mf(dofs_per_cell); + template + void + compute_matrix( + const MatrixFree &matrix_free, + const AffineConstraints &constraints, + MatrixType &matrix, + void (CLASS::*cell_operation)(FEEvaluation &) const, + const CLASS *owning_class, + const unsigned int dof_no, + const unsigned int quad_no, + const unsigned int first_selected_component) + { + compute_matrix( + matrix_free, + constraints, + matrix, + [&](auto &phi) { (owning_class->*cell_operation)(phi); }, + dof_no, + quad_no, + first_selected_component); + } - std::array, - VectorizedArrayType::size()> - matrices; + template + void + compute_matrix( + const MatrixFree &matrix_free, + const AffineConstraints &constraints_in, + MatrixType &matrix, + const std::function &)> + &cell_operation, + const std::function &, + FEFaceEvaluation &)> + &face_operation, + const std::function &)> + &boundary_operation, + const unsigned int dof_no, + const unsigned int quad_no, + const unsigned int first_selected_component) + { + using FEEvalType = FEEvaluation; + + using FEFaceEvalType = FEFaceEvaluation; + + internal::ComputeMatrixScratchData + data_cell; + + data_cell.dof_numbers = {dof_no}; + data_cell.quad_numbers = {dof_no}; + data_cell.first_selected_components = {first_selected_component}; + data_cell.batch_type = {0}; + + data_cell.op_create = + [&](const std::pair &range) { + std::vector< + std::unique_ptr>> + phi; + + if (!internal::is_fe_nothing(matrix_free, + range, + dof_no, + quad_no, + first_selected_component, + fe_degree, + n_q_points_1d)) + phi.emplace_back(std::make_unique( + matrix_free, range, dof_no, quad_no, first_selected_component)); + + return phi; + }; - std::fill_n(matrices.begin(), - VectorizedArrayType::size(), - FullMatrix(dofs_per_cell, - dofs_per_cell)); + data_cell.op_reinit = [](auto &phi, const unsigned batch) { + static_cast(*phi[0]).reinit(batch); + }; - const auto lexicographic_numbering = - matrix_free - .get_shape_info(dof_no, - quad_no, - first_selected_component, - integrator.get_active_fe_index(), - integrator.get_active_quadrature_index()) - .lexicographic_numbering; + if (cell_operation) + data_cell.op_compute = [&](auto &phi) { + cell_operation(static_cast(*phi[0])); + }; - for (auto cell = range.first; cell < range.second; ++cell) + internal::ComputeMatrixScratchData + data_face; + + data_face.dof_numbers = {dof_no, dof_no}; + data_face.quad_numbers = {dof_no, quad_no}; + data_face.first_selected_components = {first_selected_component, + first_selected_component}; + data_face.batch_type = {1, 2}; + + data_face.op_create = + [&](const std::pair &range) { + std::vector< + std::unique_ptr>> + phi; + + if (!internal::is_fe_nothing(matrix_free, + range, + dof_no, + quad_no, + first_selected_component, + fe_degree, + n_q_points_1d, + true) && + !internal::is_fe_nothing(matrix_free, + range, + dof_no, + quad_no, + first_selected_component, + fe_degree, + n_q_points_1d, + false)) { - integrator.reinit(cell); + phi.emplace_back( + std::make_unique(matrix_free, + range, + true, + dof_no, + quad_no, + first_selected_component)); + phi.emplace_back( + std::make_unique(matrix_free, + range, + false, + dof_no, + quad_no, + first_selected_component)); + } - const unsigned int n_filled_lanes = - matrix_free.n_active_entries_per_cell_batch(cell); + return phi; + }; - for (unsigned int v = 0; v < n_filled_lanes; ++v) - matrices[v] = 0.0; + data_face.op_reinit = [](auto &phi, const unsigned batch) { + static_cast(*phi[0]).reinit(batch); + static_cast(*phi[1]).reinit(batch); + }; - for (unsigned int j = 0; j < dofs_per_cell; ++j) - { - for (unsigned int i = 0; i < dofs_per_cell; ++i) - integrator.begin_dof_values()[i] = - static_cast(i == j); + if (face_operation) + data_face.op_compute = [&](auto &phi) { + face_operation(static_cast(*phi[0]), + static_cast(*phi[1])); + }; + + internal::ComputeMatrixScratchData + data_boundary; + + data_boundary.dof_numbers = {dof_no}; + data_boundary.quad_numbers = {dof_no}; + data_boundary.first_selected_components = {first_selected_component}; + data_boundary.batch_type = {1}; + + data_boundary + .op_create = [&](const std::pair &range) { + std::vector< + std::unique_ptr>> + phi; + + if (!internal::is_fe_nothing(matrix_free, + range, + dof_no, + quad_no, + first_selected_component, + fe_degree, + n_q_points_1d, + true)) + phi.emplace_back(std::make_unique( + matrix_free, range, true, dof_no, quad_no, first_selected_component)); + + return phi; + }; + + data_boundary.op_reinit = [](auto &phi, const unsigned batch) { + static_cast(*phi[0]).reinit(batch); + }; + + if (boundary_operation) + data_boundary.op_compute = [&](auto &phi) { + boundary_operation(static_cast(*phi[0])); + }; + + internal::compute_matrix( + matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix); + } + + namespace internal + { + template + void + compute_matrix( + const MatrixFree &matrix_free, + const AffineConstraints &constraints_in, + const internal::ComputeMatrixScratchData + &data_cell, + const internal::ComputeMatrixScratchData + &data_face, + const internal::ComputeMatrixScratchData + &data_boundary, + MatrixType &matrix) + { + std::unique_ptr> + constraints_for_matrix; + const AffineConstraints &constraints = + internal::create_new_affine_constraints_if_needed( + matrix, constraints_in, constraints_for_matrix); + + const auto batch_operation = + [&matrix_free, &constraints, &matrix]( + auto &data, const std::pair &range) { + if (!data.op_compute) + return; // nothing to do + + auto phi = data.op_create(range); + + const unsigned int n_blocks = phi.size(); + + if (n_blocks == 0) + return; + + Table<1, unsigned int> dofs_per_cell(n_blocks); + + Table<1, std::vector> dof_indices(n_blocks); + Table<2, std::vector> dof_indices_mf( + n_blocks, VectorizedArrayType::size()); + Table<1, std::vector> lexicographic_numbering(n_blocks); + Table<2, + std::array, + VectorizedArrayType::size()>> + matrices(n_blocks, n_blocks); + + for (unsigned int b = 0; b < n_blocks; ++b) + { + const auto &shape_info = matrix_free.get_shape_info( + data.dof_numbers[b], + data.quad_numbers[b], + data.first_selected_components[b], + phi[b]->get_active_fe_index(), + phi[b]->get_active_quadrature_index()); + + dofs_per_cell[b] = + shape_info.dofs_per_component_on_cell * shape_info.n_components; + + dof_indices[b].resize(dofs_per_cell[b]); + + for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v) + dof_indices_mf[b][v].resize(dofs_per_cell[b]); + + lexicographic_numbering[b] = shape_info.lexicographic_numbering; + } + + for (unsigned int bj = 0; bj < n_blocks; ++bj) + for (unsigned int bi = 0; bi < n_blocks; ++bi) + std::fill_n(matrices[bi][bj].begin(), + VectorizedArrayType::size(), + FullMatrix( + dofs_per_cell[bi], dofs_per_cell[bj])); + + for (auto batch = range.first; batch < range.second; ++batch) + { + data.op_reinit(phi, batch); + + const unsigned int n_filled_lanes = + data.is_face ? + matrix_free.n_active_entries_per_face_batch(batch) : + matrix_free.n_active_entries_per_cell_batch(batch); + + for (unsigned int v = 0; v < n_filled_lanes; ++v) + for (unsigned int b = 0; b < n_blocks; ++b) + { + unsigned int const cell_index = + (data.batch_type[b] == 0) ? + (batch * VectorizedArrayType::size() + v) : + ((data.batch_type[b] == 1) ? + matrix_free.get_face_info(batch).cells_interior[v] : + matrix_free.get_face_info(batch).cells_exterior[v]); + + const auto cell_iterator = matrix_free.get_cell_iterator( + cell_index / VectorizedArrayType::size(), + cell_index % VectorizedArrayType::size(), + data.dof_numbers[b]); + + if (matrix_free.get_mg_level() != + numbers::invalid_unsigned_int) + cell_iterator->get_mg_dof_indices(dof_indices[b]); + else + cell_iterator->get_dof_indices(dof_indices[b]); + + for (unsigned int j = 0; j < dof_indices[b].size(); ++j) + dof_indices_mf[b][v][j] = + dof_indices[b][lexicographic_numbering[b][j]]; + } - local_vmult(integrator); + for (unsigned int bj = 0; bj < n_blocks; ++bj) + { + for (unsigned int j = 0; j < dofs_per_cell[bj]; ++j) + { + for (unsigned int bi = 0; bi < n_blocks; ++bi) + for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i) + phi[bi]->begin_dof_values()[i] = + (bj == bi) ? static_cast(i == j) : 0.0; + + data.op_compute(phi); + + for (unsigned int bi = 0; bi < n_blocks; ++bi) + for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i) + for (unsigned int v = 0; v < n_filled_lanes; ++v) + matrices[bi][bj][v](i, j) = + phi[bi]->begin_dof_values()[i][v]; + } - for (unsigned int i = 0; i < dofs_per_cell; ++i) for (unsigned int v = 0; v < n_filled_lanes; ++v) - matrices[v](i, j) = integrator.begin_dof_values()[i][v]; - } + for (unsigned int bi = 0; bi < n_blocks; ++bi) + constraints.distribute_local_to_global( + matrices[bi][bj][v], + dof_indices_mf[bi][v], + dof_indices_mf[bj][v], + matrix); + } + } + }; - for (unsigned int v = 0; v < n_filled_lanes; ++v) - { - const auto cell_v = - matrix_free.get_cell_iterator(cell, v, dof_no); + const auto cell_operation_wrapped = + [&](const auto &, auto &, const auto &, const auto range) { + batch_operation(data_cell, range); + }; - if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int) - cell_v->get_mg_dof_indices(dof_indices); - else - cell_v->get_dof_indices(dof_indices); + const auto face_operation_wrapped = + [&](const auto &, auto &, const auto &, const auto range) { + batch_operation(data_face, range); + }; - for (unsigned int j = 0; j < dof_indices.size(); ++j) - dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]]; + const auto boundary_operation_wrapped = + [&](const auto &, auto &, const auto &, const auto range) { + batch_operation(data_boundary, range); + }; - constraints.distribute_local_to_global(matrices[v], - dof_indices_mf, - dst); - } - } - }, - matrix, - matrix); + if (data_face.op_compute || data_boundary.op_compute) + { + matrix_free.template loop( + cell_operation_wrapped, + face_operation_wrapped, + boundary_operation_wrapped, + matrix, + matrix); + } + else + matrix_free.template cell_loop( + cell_operation_wrapped, matrix, matrix); - matrix.compress(VectorOperation::add); - } + matrix.compress(VectorOperation::add); + } + } // namespace internal template &) const, + void (CLASS::*face_operation)(FEFaceEvaluation &, + FEFaceEvaluation &) + const, + void (CLASS::*boundary_operation)(FEFaceEvaluation &) + const, const CLASS *owning_class, const unsigned int dof_no, const unsigned int quad_no, @@ -1298,7 +2316,11 @@ namespace MatrixFreeTools matrix_free, constraints, matrix, - [&](auto &feeval) { (owning_class->*cell_operation)(feeval); }, + [&](auto &phi) { (owning_class->*cell_operation)(phi); }, + [&](auto &phi_m, auto &phi_p) { + (owning_class->*face_operation)(phi_m, phi_p); + }, + [&](auto &phi) { (owning_class->*boundary_operation)(phi); }, dof_no, quad_no, first_selected_component); diff --git a/tests/matrix_free/compute_diagonal_08.cc b/tests/matrix_free/compute_diagonal_08.cc new file mode 100644 index 0000000000..8f2eecb09e --- /dev/null +++ b/tests/matrix_free/compute_diagonal_08.cc @@ -0,0 +1,276 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2022 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// Compute diagonal and matrix for SIP with MatrixFreeTools. + +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test() +{ + const int fe_degree = 3; + const int n_points = fe_degree + 1; + const int n_components = 1; + using Number = double; + using VectorizedArrayType = VectorizedArray; + using VectorType = Vector; + + Triangulation tria; + GridGenerator::hyper_ball(tria); + tria.refine_global(0); + + const FE_DGQ fe_q(fe_degree); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe_q); + + AffineConstraints constraints; + constraints.close(); + + typename MatrixFree::AdditionalData + additional_data; + additional_data.mapping_update_flags = update_values | update_gradients; + additional_data.mapping_update_flags_inner_faces = + update_values | update_gradients; + additional_data.mapping_update_flags_boundary_faces = + update_values | update_gradients; + + MappingQ mapping(1); + QGauss<1> quad(fe_degree + 1); + + MatrixFree matrix_free; + matrix_free.reinit(mapping, dof_handler, constraints, quad, additional_data); + + const auto cell_operation = [&](auto &phi) { + phi.evaluate(EvaluationFlags::gradients); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(EvaluationFlags::gradients); + }; + + const auto face_operation = [&](auto &phi_m, auto &phi_p) { + phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + phi_p.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + + VectorizedArrayType sigmaF = + (std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs( + (phi_m.normal_vector(0) * phi_p.inverse_jacobian(0))[dim - 1])) * + (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = + (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5; + VectorizedArrayType average_valgrad = + phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q); + average_valgrad = average_value * 2. * sigmaF - average_valgrad * 0.5; + phi_m.submit_normal_derivative(-average_value, q); + phi_p.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + phi_p.submit_value(-average_valgrad, q); + } + phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + phi_p.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + }; + + const auto boundary_operation = [&](auto &phi_m) { + phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + VectorizedArrayType sigmaF = + std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) * + Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = phi_m.get_value(q); + VectorizedArrayType average_valgrad = -phi_m.get_normal_derivative(q); + average_valgrad += average_value * sigmaF * 2.0; + phi_m.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + } + + phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + }; + + const auto vmult = [&](auto &dst, const auto &src) { + matrix_free.template loop( + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto range) { + FEEvaluation + phi(matrix_free); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + cell_operation(phi); + phi.set_dof_values(dst); + } + }, + [&]( + const auto &matrix_free, auto &dst, const auto &src, const auto range) { + FEFaceEvaluation + phi_m(matrix_free, true); + FEFaceEvaluation + phi_p(matrix_free, false); + + for (unsigned int face = range.first; face < range.second; ++face) + { + phi_m.reinit(face); + phi_p.reinit(face); + phi_m.read_dof_values(src); + phi_p.read_dof_values(src); + face_operation(phi_m, phi_p); + phi_m.distribute_local_to_global(dst); + phi_p.distribute_local_to_global(dst); + } + }, + [&](const auto &matrix_free, + auto &dst, + const auto &src, + const auto face_range) { + FEFaceEvaluation + phi_m(matrix_free, true); + + for (unsigned int face = face_range.first; face < face_range.second; + face++) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + boundary_operation(phi_m); + phi_m.distribute_local_to_global(dst); + } + }, + dst, + src, + true); + }; + + // Compute matrix and diagonal (brute-force) + FullMatrix full_matrix(dof_handler.n_dofs(), dof_handler.n_dofs()); + VectorType diagonal(dof_handler.n_dofs()); + + for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) + { + VectorType src(dof_handler.n_dofs()); + VectorType dst(dof_handler.n_dofs()); + + src[i] = 1.0; + + vmult(dst, src); + + diagonal[i] = dst[i]; + for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j) + full_matrix[j][i] = dst[j]; + } + + + // Compute matrix via MatrixFreeTools + FullMatrix full_matrix_mft(dof_handler.n_dofs(), + dof_handler.n_dofs()); + + MatrixFreeTools::compute_matrix(matrix_free, + constraints, + full_matrix_mft, + cell_operation, + face_operation, + boundary_operation); + + for (unsigned int i = 0; i < full_matrix.m(); ++i) + for (unsigned int j = 0; j < full_matrix.n(); ++j) + if (std::abs(full_matrix[i][j] - full_matrix_mft[i][j]) > 1e-6) + { + full_matrix.print_formatted( + deallog.get_file_stream(), 3, true, 0, "0.000e+00"); + deallog << std::endl; + + full_matrix_mft.print_formatted( + deallog.get_file_stream(), 3, true, 0, "0.000e+00"); + deallog << std::endl; + } + + // Compute diagonal via MatrixFreeTools + VectorType diagonal_mft; + matrix_free.initialize_dof_vector(diagonal_mft); + + MatrixFreeTools::compute_diagonal(matrix_free, + diagonal_mft, + cell_operation, + face_operation, + boundary_operation); + + for (unsigned int i = 0; i < diagonal.size(); ++i) + if (std::abs(diagonal[i] - diagonal_mft[i]) > 1e-6) + { + diagonal.print(deallog.get_file_stream()); + deallog << std::endl; + + diagonal_mft.print(deallog.get_file_stream()); + deallog << std::endl; + } + + deallog << "OK!" << std::endl; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + test<2>(); +} diff --git a/tests/matrix_free/compute_diagonal_08.output b/tests/matrix_free/compute_diagonal_08.output new file mode 100644 index 0000000000..505548540f --- /dev/null +++ b/tests/matrix_free/compute_diagonal_08.output @@ -0,0 +1,2 @@ + +DEAL:0::OK!