From d00303e26bfececbca8d0840fd94724e1c08dc5f Mon Sep 17 00:00:00 2001 From: Magdalena Schreter Date: Mon, 31 Jan 2022 16:22:46 +0100 Subject: [PATCH] add typename VectorType to compute_diagonal Co-authored-by: peterrum --- doc/news/changes/minor/20220202Schreter | 4 + include/deal.II/matrix_free/tools.h | 99 +++++++++----- tests/matrix_free/compute_diagonal_07.cc | 137 +++++++++++++++++++ tests/matrix_free/compute_diagonal_07.output | 13 ++ 4 files changed, 222 insertions(+), 31 deletions(-) create mode 100644 doc/news/changes/minor/20220202Schreter create mode 100644 tests/matrix_free/compute_diagonal_07.cc create mode 100644 tests/matrix_free/compute_diagonal_07.output diff --git a/doc/news/changes/minor/20220202Schreter b/doc/news/changes/minor/20220202Schreter new file mode 100644 index 0000000000..9cea9e636e --- /dev/null +++ b/doc/news/changes/minor/20220202Schreter @@ -0,0 +1,4 @@ +Improved: In MatrixFreeTools::compute_diagonal a template argument VectorType +has been introduced to be applicable for arbitrary vector types. +
+(Magdalena Schreter, Peter Munch, 2022/02/02) diff --git a/include/deal.II/matrix_free/tools.h b/include/deal.II/matrix_free/tools.h index b8e608b083..a48808e3e9 100644 --- a/include/deal.II/matrix_free/tools.h +++ b/include/deal.II/matrix_free/tools.h @@ -56,11 +56,12 @@ namespace MatrixFreeTools int n_q_points_1d, int n_components, typename Number, - typename VectorizedArrayType> + typename VectorizedArrayType, + typename VectorType> void compute_diagonal( const MatrixFree &matrix_free, - LinearAlgebra::distributed::Vector & diagonal_global, + VectorType & diagonal_global, const std::function + typename VectorizedArrayType, + typename VectorType> void compute_diagonal( const MatrixFree &matrix_free, - LinearAlgebra::distributed::Vector & diagonal_global, + VectorType & diagonal_global, void (CLASS::*cell_operation)(FEEvaluationphi.reinit(cell); // STEP 1: get relevant information from FEEvaluation - const unsigned int first_selected_component = - phi.get_first_selected_component(); const auto & dof_info = phi.get_dof_info(); const unsigned int n_fe_components = dof_info.start_components.back(); const unsigned int dofs_per_component = phi.dofs_per_component; const auto & matrix_free = phi.get_matrix_free(); + // if we have a block vector with components with the same DoFHandler, + // each component is described with same set of constraints and + // we consider the shift in components only during access of the global + // vector + 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); @@ -322,9 +329,6 @@ namespace MatrixFreeTools if (n_components == 1 || n_fe_components == 1) { - AssertDimension(n_components, - 1); // TODO: currently no block vector supported - unsigned int ind_local = 0; for (; index_indicators != next_index_indicators; ++index_indicators, ++ind_local) @@ -420,7 +424,6 @@ namespace MatrixFreeTools } } } - // STEP 2b: sort and make unique // sort vector @@ -674,7 +677,9 @@ namespace MatrixFreeTools // set size locally-relevant diagonal for (unsigned int v = 0; v < n_lanes_filled; ++v) diagonals_local_constrained[v].assign( - c_pools[v].row_lid_to_gid.size(), Number(0.0)); + c_pools[v].row_lid_to_gid.size() * + (n_fe_components == 1 ? n_components : 1), + Number(0.0)); } void @@ -692,7 +697,15 @@ namespace MatrixFreeTools void submit() { - const auto ith_column = phi.begin_dof_values(); + // if we have a block vector with components with the same DoFHandler, + // we need to figure out which component and which DoF within the + // comonent are we currently considering + const unsigned int n_fe_components = + phi.get_dof_info().start_components.back(); + const unsigned int comp = + n_fe_components == 1 ? i / phi.dofs_per_component : 0; + const unsigned int i_comp = + n_fe_components == 1 ? (i % phi.dofs_per_component) : i; // apply local constraint matrix from left and from right: // loop over all rows of transposed constrained matrix @@ -710,43 +723,66 @@ namespace MatrixFreeTools const auto scale_iterator = std::lower_bound(c_pool.col.begin() + c_pool.row[j], c_pool.col.begin() + c_pool.row[j + 1], - i); + i_comp); // explanation: j-th row of C_e^T is empty (see above) if (scale_iterator == c_pool.col.begin() + c_pool.row[j + 1]) continue; // explanation: C_e^T(j,i) is zero (see above) - if (*scale_iterator != i) + if (*scale_iterator != i_comp) continue; // apply constraint matrix from the left Number temp = 0.0; for (unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k) - temp += c_pool.val[k] * ith_column[c_pool.col[k]][v]; + temp += c_pool.val[k] * + phi.begin_dof_values()[comp * phi.dofs_per_component + + c_pool.col[k]][v]; // apply constraint matrix from the right - diagonals_local_constrained[v][j] += + diagonals_local_constrained + [v][j + comp * c_pools[v].row_lid_to_gid.size()] += temp * c_pool.val[std::distance(c_pool.col.begin(), scale_iterator)]; } } } - void - distribute_local_to_global( - LinearAlgebra::distributed::Vector &diagonal_global) + template + inline void + distribute_local_to_global(VectorType &diagonal_global) { // STEP 4: assembly results: add into global vector + const unsigned int n_fe_components = + phi.get_dof_info().start_components.back(); + const unsigned int first_selected_component = + phi.get_first_selected_component(); + for (unsigned int v = 0; v < phi.get_matrix_free().n_active_entries_per_cell_batch( phi.get_current_cell_index()); ++v) + // if we have a block vector with components with the same + // DoFHandler, we need to loop over all components manully and + // need to apply the correct shift for (unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j) - ::dealii::internal::vector_access_add( - diagonal_global, - c_pools[v].row_lid_to_gid[j], - diagonals_local_constrained[v][j]); + for (unsigned int comp = 0; + comp < (n_fe_components == 1 ? + static_cast(n_components) : + 1); + ++comp) + ::dealii::internal::vector_access_add( + *::dealii::internal::BlockVectorSelector< + VectorType, + IsBlockVector::value>:: + get_vector_component(diagonal_global, + n_fe_components == 1 ? + comp + first_selected_component : + 0), + c_pools[v].row_lid_to_gid[j], + diagonals_local_constrained + [v][j + comp * c_pools[v].row_lid_to_gid.size()]); } private: @@ -778,11 +814,12 @@ namespace MatrixFreeTools int n_q_points_1d, int n_components, typename Number, - typename VectorizedArrayType> + typename VectorizedArrayType, + typename VectorType> void compute_diagonal( const MatrixFree &matrix_free, - LinearAlgebra::distributed::Vector & diagonal_global, + VectorType & diagonal_global, const std::function; - int dummy = 0; matrix_free.template cell_loop( [&](const MatrixFree &matrix_free, - LinearAlgebra::distributed::Vector & diagonal_global, + VectorType & diagonal_global, const int &, const std::pair &range) mutable { FEEvaluation + typename VectorizedArrayType, + typename VectorType> void compute_diagonal( const MatrixFree &matrix_free, - LinearAlgebra::distributed::Vector & diagonal_global, + VectorType & diagonal_global, void (CLASS::*cell_operation)(FEEvaluation( + VectorizedArrayType, + VectorType>( matrix_free, diagonal_global, [&](auto &feeval) { (owning_class->*cell_operation)(feeval); }, diff --git a/tests/matrix_free/compute_diagonal_07.cc b/tests/matrix_free/compute_diagonal_07.cc new file mode 100644 index 0000000000..eb0c3f087f --- /dev/null +++ b/tests/matrix_free/compute_diagonal_07.cc @@ -0,0 +1,137 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Similar to compute_diagonal_02 but testing block vectorss. + +#include + +#include "compute_diagonal_util.h" + +using namespace dealii; + +template > +void +test() +{ + Triangulation tria; + GridGenerator::hyper_ball(tria); + tria.refine_global(0); + + const FE_Q fe_q(fe_degree); + const FESystem fe_system(fe_q, n_components); + + DoFHandler dof_handler_q(tria); + dof_handler_q.distribute_dofs(fe_q); + + DoFHandler dof_handler_system(tria); + dof_handler_system.distribute_dofs(fe_system); + + AffineConstraints constraints_q; + DoFTools::make_hanging_node_constraints(dof_handler_q, constraints_q); + DoFTools::make_zero_boundary_constraints(dof_handler_q, constraints_q); + constraints_q.close(); + + constraints_q.print(std::cout); + + AffineConstraints constraints_system; + DoFTools::make_hanging_node_constraints(dof_handler_system, + constraints_system); + DoFTools::make_zero_boundary_constraints(dof_handler_system, + constraints_system); + constraints_system.close(); + + constraints_system.print(std::cout); + + typename MatrixFree::AdditionalData + additional_data; + additional_data.mapping_update_flags = update_values | update_gradients; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData:: + TasksParallelScheme::none; + + MappingQ mapping(1); + QGauss<1> quad(fe_degree + 1); + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, + std::vector *>{&dof_handler_system, &dof_handler_q}, + std::vector *>{&constraints_system, + &constraints_q}, + quad, + additional_data); + + const auto kernel = [](FEEvaluation &phi) { + phi.evaluate(false, true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + phi.submit_symmetric_gradient(2.0 * phi.get_symmetric_gradient(q), q); + } + phi.integrate(false, true); + }; + + LinearAlgebra::distributed::Vector diagonal_global; + matrix_free.initialize_dof_vector(diagonal_global, 0); + + MatrixFreeTools::compute_diagonal(matrix_free, + diagonal_global, + kernel, + 0); + + diagonal_global.print(deallog.get_file_stream()); + + LinearAlgebra::distributed::BlockVector diagonal_global_block(dim); + + for (unsigned int comp = 0; comp < dim; ++comp) + matrix_free.initialize_dof_vector(diagonal_global_block.block(comp), 1); + + MatrixFreeTools::compute_diagonal(matrix_free, + diagonal_global_block, + kernel, + 1); + + diagonal_global_block.print(deallog.get_file_stream()); +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + test<2, 1>(); +} diff --git a/tests/matrix_free/compute_diagonal_07.output b/tests/matrix_free/compute_diagonal_07.output new file mode 100644 index 0000000000..ecf8996667 --- /dev/null +++ b/tests/matrix_free/compute_diagonal_07.output @@ -0,0 +1,13 @@ + +Process #0 +Local range: [0, 16), global size: 16 +Vector data: +0.000e+00 0.000e+00 0.000e+00 0.000e+00 4.621e+00 4.621e+00 4.621e+00 4.621e+00 0.000e+00 0.000e+00 4.621e+00 4.621e+00 4.621e+00 4.621e+00 0.000e+00 0.000e+00 +Process #0 +Local range: [0, 8), global size: 8 +Vector data: +0.000e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 +Process #0 +Local range: [0, 8), global size: 8 +Vector data: +0.000e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 4.621e+00 4.621e+00 0.000e+00 -- 2.39.5