From 240f257cd2447616b8cd086202b7f51fce6a74ee Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 2 Oct 2022 16:38:59 +0200 Subject: [PATCH] TensorProductMatrixSymmetricSumCollection: allow to disable compression --- include/deal.II/lac/tensor_product_matrix.h | 164 ++++++++++--- .../numerics/tensor_product_matrix_creator.h | 40 ++-- tests/lac/tensor_product_matrix.h | 133 +++++++++++ tests/lac/tensor_product_matrix_08.cc | 217 ++++++++++++++++++ ..._product_matrix_08.with_lapack=true.output | 31 +++ 5 files changed, 543 insertions(+), 42 deletions(-) create mode 100644 tests/lac/tensor_product_matrix.h create mode 100644 tests/lac/tensor_product_matrix_08.cc create mode 100644 tests/lac/tensor_product_matrix_08.with_lapack=true.output diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index a59f68dae0..27f27169ab 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -301,8 +301,6 @@ namespace internal std::bitset mask; - using ScalarNumber = typename Number::value_type; - for (unsigned int v = 0; v < width; ++v) { ScalarNumber a = 0.0; @@ -370,6 +368,28 @@ class TensorProductMatrixSymmetricSumCollection using MatrixPairType = std::pair, Table<2, Number>>; public: + /** + * Struct to configure TensorProductMatrixSymmetricSumCollection. + */ + struct AdditionalData + { + /** + * Constructor. + */ + AdditionalData(const bool compress_matrices = true); + + /** + * Try to compress internal matrices. Default: true. + */ + bool compress_matrices; + }; + + /** + * Consturctor. + */ + TensorProductMatrixSymmetricSumCollection( + const AdditionalData &additional_data = AdditionalData()); + /** * Allocate memory. The parameter @p specifies the maximum value * of the index used in invert() and apply_inverse(). @@ -382,10 +402,9 @@ public: * stiffness matrices @p Ks to the stored data, looking out for * compression possibilities. */ + template void - insert(const unsigned int index, - const std::array, dim> &Ms, - const std::array, dim> &Ks); + insert(const unsigned int index, const T &Ms, const T &Ks); /** * Finalize setup. This function computes, e.g., the @@ -421,6 +440,17 @@ public: storage_size() const; private: + /** + * Try to compress matrices. + */ + const bool compress_matrices; + + /** + * Container used to collect 1D matrices if no compression is + * requested. The memory is freed during finalize(). + */ + std::vector mass_and_derivative_matrices; + /** * Container used during setup to determine the unique 1D * matrices. The memory is freed during finalize(). @@ -1041,35 +1071,69 @@ TensorProductMatrixSymmetricSum::reinit( +template +TensorProductMatrixSymmetricSumCollection:: + AdditionalData::AdditionalData(const bool compress_matrices) + : compress_matrices(compress_matrices) +{} + + + +template +TensorProductMatrixSymmetricSumCollection:: + TensorProductMatrixSymmetricSumCollection( + const AdditionalData &additional_data) + : compress_matrices(additional_data.compress_matrices) +{} + + + template void TensorProductMatrixSymmetricSumCollection::reserve( const unsigned int size) { - indices.assign(size * dim, numbers::invalid_unsigned_int); + if (compress_matrices == false) + mass_and_derivative_matrices.resize(size * dim); + else + indices.assign(size * dim, numbers::invalid_unsigned_int); } template +template void TensorProductMatrixSymmetricSumCollection::insert( - const unsigned int index, - const std::array, dim> &Ms, - const std::array, dim> &Ks) + const unsigned int index, + const T & Ms_in, + const T & Ks_in) { + const auto Ms = + internal::TensorProductMatrixSymmetricSum::convert(Ms_in); + const auto Ks = + internal::TensorProductMatrixSymmetricSum::convert(Ks_in); + for (unsigned int d = 0; d < dim; ++d) { const MatrixPairType matrix(Ms[d], Ks[d]); - const auto ptr = cache.find(matrix); - - if (ptr != cache.end()) - indices[index * dim + d] = ptr->second; + if (compress_matrices == false) + { + mass_and_derivative_matrices[index * dim + d] = matrix; + } else { - indices[index * dim + d] = cache.size(); - cache[matrix] = cache.size(); + const auto ptr = cache.find(matrix); + + if (ptr != cache.end()) + indices[index * dim + d] = ptr->second; + else + { + const auto size = cache.size(); + indices[index * dim + d] = size; + cache[matrix] = size; + } } } } @@ -1080,9 +1144,6 @@ template void TensorProductMatrixSymmetricSumCollection::finalize() { - this->vector_ptr.resize(cache.size() + 1); - this->matrix_ptr.resize(cache.size() + 1); - const auto store = [&](const unsigned int index, const MatrixPairType &M_and_K) { std::array, 1> mass_matrix; @@ -1114,8 +1175,47 @@ TensorProductMatrixSymmetricSumCollection::finalize() } }; - if (cache.size() == indices.size()) + if (compress_matrices == false) + { + // case 1) no compression requested + + AssertDimension(cache.size(), 0); + AssertDimension(indices.size(), 0); + + this->vector_ptr.resize(mass_and_derivative_matrices.size() + 1); + this->matrix_ptr.resize(mass_and_derivative_matrices.size() + 1); + + for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i) + { + const auto &M = mass_and_derivative_matrices[i].first; + + this->vector_ptr[i + 1] = M.n_rows(); + this->matrix_ptr[i + 1] = M.n_rows() * M.n_cols(); + } + + for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i) + { + this->vector_ptr[i + 1] += this->vector_ptr[i]; + this->matrix_ptr[i + 1] += this->matrix_ptr[i]; + } + + this->mass_matrices.resize_fast(matrix_ptr.back()); + this->derivative_matrices.resize_fast(matrix_ptr.back()); + this->eigenvectors.resize_fast(matrix_ptr.back()); + this->eigenvalues.resize_fast(vector_ptr.back()); + + for (unsigned int i = 0; i < mass_and_derivative_matrices.size(); ++i) + store(i, mass_and_derivative_matrices[i]); + + mass_and_derivative_matrices.clear(); + } + else if (cache.size() == indices.size()) { + // case 2) compression requested but none possible + + this->vector_ptr.resize(cache.size() + 1); + this->matrix_ptr.resize(cache.size() + 1); + std::map inverted_cache; for (const auto &i : cache) @@ -1135,18 +1235,24 @@ TensorProductMatrixSymmetricSumCollection::finalize() this->matrix_ptr[i + 1] += this->matrix_ptr[i]; } - this->mass_matrices.resize(matrix_ptr.back()); - this->derivative_matrices.resize(matrix_ptr.back()); - this->eigenvectors.resize(matrix_ptr.back()); - this->eigenvalues.resize(vector_ptr.back()); + this->mass_matrices.resize_fast(matrix_ptr.back()); + this->derivative_matrices.resize_fast(matrix_ptr.back()); + this->eigenvectors.resize_fast(matrix_ptr.back()); + this->eigenvalues.resize_fast(vector_ptr.back()); for (unsigned int i = 0; i < indices.size(); ++i) store(i, inverted_cache[indices[i]]); indices.clear(); + cache.clear(); } else { + // case 2) compression requested but none possible + + this->vector_ptr.resize(cache.size() + 1); + this->matrix_ptr.resize(cache.size() + 1); + for (const auto &i : cache) { const auto &M = i.first.first; @@ -1161,16 +1267,16 @@ TensorProductMatrixSymmetricSumCollection::finalize() this->matrix_ptr[i + 1] += this->matrix_ptr[i]; } - this->mass_matrices.resize(matrix_ptr.back()); - this->derivative_matrices.resize(matrix_ptr.back()); - this->eigenvectors.resize(matrix_ptr.back()); - this->eigenvalues.resize(vector_ptr.back()); + this->mass_matrices.resize_fast(matrix_ptr.back()); + this->derivative_matrices.resize_fast(matrix_ptr.back()); + this->eigenvectors.resize_fast(matrix_ptr.back()); + this->eigenvalues.resize_fast(vector_ptr.back()); for (const auto &i : cache) store(i.second, i.first); - } - cache.clear(); + cache.clear(); + } } diff --git a/include/deal.II/numerics/tensor_product_matrix_creator.h b/include/deal.II/numerics/tensor_product_matrix_creator.h index 1242e64b80..184851b04f 100644 --- a/include/deal.II/numerics/tensor_product_matrix_creator.h +++ b/include/deal.II/numerics/tensor_product_matrix_creator.h @@ -99,6 +99,19 @@ namespace TensorProductMatrixCreator { namespace internal { + template + void + clear_row_and_column(const unsigned int n_dofs_1D_with_overlap, + const unsigned int n, + FullMatrix &matrix) + { + for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i) + { + matrix[i][n] = 0.0; + matrix[n][i] = 0.0; + } + } + template std::tuple, FullMatrix, bool> create_reference_mass_and_stiffness_matrices( @@ -150,7 +163,8 @@ namespace TensorProductMatrixCreator fe_values.JxW(q_index)); } - return {mass_matrix_reference, derivative_matrix_reference, false}; + return std::tuple, FullMatrix, bool>{ + mass_matrix_reference, derivative_matrix_reference, false}; } } // namespace internal @@ -191,14 +205,6 @@ namespace TensorProductMatrixCreator std::array, dim> Ms; std::array, dim> Ks; - const auto clear_row_and_column = [&](const unsigned int n, auto &matrix) { - for (unsigned int i = 0; i < n_dofs_1D_with_overlap; ++i) - { - matrix[i][n] = 0.0; - matrix[n][i] = 0.0; - } - }; - for (unsigned int d = 0; d < dim; ++d) { Ms[d].reinit(n_dofs_1D_with_overlap, n_dofs_1D_with_overlap); @@ -235,8 +241,12 @@ namespace TensorProductMatrixCreator { // left DBC const unsigned i0 = n_overlap - 1; - clear_row_and_column(i0, Ms[d]); - clear_row_and_column(i0, Ks[d]); + internal::clear_row_and_column(n_dofs_1D_with_overlap, + i0, + Ms[d]); + internal::clear_row_and_column(n_dofs_1D_with_overlap, + i0, + Ks[d]); } else if (boundary_ids[d][0] == LaplaceBoundaryType::neumann) { @@ -268,8 +278,12 @@ namespace TensorProductMatrixCreator { // right DBC const unsigned i0 = n_overlap + n_dofs_1D - 2; - clear_row_and_column(i0, Ms[d]); - clear_row_and_column(i0, Ks[d]); + internal::clear_row_and_column(n_dofs_1D_with_overlap, + i0, + Ms[d]); + internal::clear_row_and_column(n_dofs_1D_with_overlap, + i0, + Ks[d]); } else if (boundary_ids[d][1] == LaplaceBoundaryType::neumann) { diff --git a/tests/lac/tensor_product_matrix.h b/tests/lac/tensor_product_matrix.h new file mode 100644 index 0000000000..9e311eb63f --- /dev/null +++ b/tests/lac/tensor_product_matrix.h @@ -0,0 +1,133 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + +#pragma once + +#include +#include + +namespace dealii +{ + namespace GridTools + { + // Compute harminic extent of all non-artificial cells. + template + std::vector> + compute_harmonic_cell_extent(const Mapping & mapping, + const Triangulation & triangulation, + const Quadrature &quadrature) + { + std::vector> result( + triangulation.n_active_cells()); + + FE_Nothing fe_nothing; + FEFaceValues fe_face_values_0(mapping, + fe_nothing, + quadrature, + update_quadrature_points); + FEFaceValues fe_face_values_1(mapping, + fe_nothing, + quadrature, + update_quadrature_points); + + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->is_artificial() == false) + { + for (unsigned int d = 0; d < dim; ++d) + { + fe_face_values_0.reinit(cell, 2 * d + 0); + fe_face_values_1.reinit(cell, 2 * d + 1); + + double extent = 0.0; + + for (unsigned int q = 0; q < quadrature.size(); ++q) + extent += fe_face_values_0.quadrature_point(q).distance( + fe_face_values_1.quadrature_point(q)) * + quadrature.weight(q); + + result[cell->active_cell_index()][d] = extent; + } + } + + return result; + } + + // Compute harmonic extent of each locally owned cell including of each + // of its neighbors. If there is no neigbor, its extent is zero. + template + std::vector> + compute_harmonic_patch_extent(const Mapping & mapping, + const Triangulation & triangulation, + const Quadrature &quadrature) + { + // 1) compute extent of each non-artificial cell + const auto harmonic_cell_extents = + GridTools::compute_harmonic_cell_extent(mapping, + triangulation, + quadrature); + + // 2) accumulate for each face the normal extent for the + // neigboring cell(s); here we also consider periodicies + std::vector face_extent(triangulation.n_faces(), 0.0); + + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->is_artificial() == false) + for (unsigned int d = 0; d < dim; ++d) + { + const auto extent = + harmonic_cell_extents[cell->active_cell_index()][d]; + + const auto add_extent_to_faces = [&](const unsigned int face_no) { + face_extent[cell->face(face_no)->index()] += extent; + + if (cell->has_periodic_neighbor(face_no) && + (cell->periodic_neighbor(face_no)->is_artificial() == + false)) + face_extent[cell->periodic_neighbor(face_no) + ->face(cell->periodic_neighbor_face_no(face_no)) + ->index()] += extent; + }; + + add_extent_to_faces(2 * d + 0); // face 0 + add_extent_to_faces(2 * d + 1); // face 1 + } + + // 3) collect cell extent including those of the neighboring + // cells, which corresponds to the difference of extent of the + // current cell and the face extent + std::vector> result( + triangulation.n_active_cells()); + + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->is_locally_owned()) + for (unsigned int d = 0; d < dim; ++d) + { + const auto cell_extent = + harmonic_cell_extents[cell->active_cell_index()][d]; + + const auto index = cell->active_cell_index(); + + result[index][d][0] = + face_extent[cell->face(2 * d + 0)->index()] - cell_extent; + result[index][d][1] = cell_extent; + result[index][d][2] = + face_extent[cell->face(2 * d + 1)->index()] - cell_extent; + } + + return result; + } + + } // namespace GridTools +} // namespace dealii diff --git a/tests/lac/tensor_product_matrix_08.cc b/tests/lac/tensor_product_matrix_08.cc new file mode 100644 index 0000000000..508869c3b9 --- /dev/null +++ b/tests/lac/tensor_product_matrix_08.cc @@ -0,0 +1,217 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + + +// Test compression of TensorProductMatrixSymmetricSumCollection. + +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +#include "./tensor_product_matrix.h" + +template +void +do_test_mesh(const Mapping &mapping, const Triangulation &tria) +{ + using FDM = TensorProductMatrixSymmetricSumCollection; + + const unsigned int fe_degree = 3; + const unsigned int n_overlap = 1; + + FE_Q fe(fe_degree); + FE_Q<1> fe_1D(fe_degree); + QGauss<1> quadrature_1D(fe_degree + 1); + QGauss quadrature_face(fe_degree + 1); + + const auto harmonic_patch_extent = + GridTools::compute_harmonic_patch_extent(mapping, tria, quadrature_face); + + FDM collection_0(typename FDM::AdditionalData(true)); + FDM collection_1(typename FDM::AdditionalData(false)); + + collection_0.reserve(tria.n_active_cells()); + collection_1.reserve(tria.n_active_cells()); + + for (const auto &cell : tria.active_cell_iterators()) + { + const auto &patch_extent = + harmonic_patch_extent[cell->active_cell_index()]; + + const auto M_and_K = TensorProductMatrixCreator:: + create_laplace_tensor_product_matrix( + cell, + std::set{0}, + std::set{}, + fe_1D, + quadrature_1D, + patch_extent, + n_overlap); + + collection_0.insert(cell->active_cell_index(), + M_and_K.first, + M_and_K.second); + collection_1.insert(cell->active_cell_index(), + M_and_K.first, + M_and_K.second); + } + + collection_0.finalize(); + collection_1.finalize(); + + deallog << "Storage sizes: " << collection_0.storage_size() << " " + << collection_1.storage_size() << std::endl; + + Vector src(fe.n_dofs_per_cell()); + Vector dst(fe.n_dofs_per_cell()); + AlignedVector tmp; + FullMatrix matrix_0(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); + FullMatrix matrix_1(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); + + for (unsigned int cell = 0; cell < tria.n_active_cells(); ++cell) + { + for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) + { + for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) + src[j] = i == j; + + collection_0.apply_inverse(cell, + make_array_view(dst), + make_array_view(src), + tmp); + for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) + matrix_0[j][i] = dst[j]; + + collection_1.apply_inverse(cell, + make_array_view(dst), + make_array_view(src), + tmp); + for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) + matrix_1[j][i] = dst[j]; + } + + + FloatingPointComparator comp(1e-5, false); + + Assert((comp.compare(matrix_0, matrix_1) == + FloatingPointComparator::ComparisonResult::equal), + ExcInternalError()); + } + + deallog << "OK!" << std::endl; + deallog << std::endl; +} + + +template +void +do_test() +{ + const int dim = 2; + + MappingQ1 mapping; + + const std::vector> settings = { + std::tuple{0, false, false}, // hyper-cube + std::tuple{0, true, false}, // hyper-cube with PBC + std::tuple{1, false, false}, // hyper-rectangle + std::tuple{1, + true, + false}, // hyper-rectangle with PBC + std::tuple{0, false, true} // deformed hyper-cube + }; + + for (const auto setting : settings) + { + Triangulation tria; + + const auto do_pbc = std::get<1>(setting); + + switch (std::get<0>(setting)) + { + case 0: + GridGenerator::hyper_cube(tria, 0.0, 1.0, do_pbc); + break; + case 1: + GridGenerator::subdivided_hyper_rectangle( + tria, {1, 2}, {0, 0}, {1.0, 1.0}, do_pbc); + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + + if (do_pbc) + { + std::vector::cell_iterator>> + periodic_faces; + for (unsigned int d = 0; d < dim; ++d) + GridTools::collect_periodic_faces( + tria, 2 * d, 2 * d + 1, d, periodic_faces); + tria.add_periodicity(periodic_faces); + } + + tria.refine_global(3); + + MappingQCache mapping(3); + + mapping.initialize( + MappingQ1(), + tria, + [&](const typename Triangulation::cell_iterator &, + const Point &point) -> Point { + Point result; + + if (std::get<2>(setting) == false) // Cartesian mesh + return result; + + for (unsigned int d = 0; d < dim; ++d) + result[d] = std::pow(point[d], d + 0.9) * + std::pow(point[(d + 1) % dim], d + 1.1); + + return result; + }, + true); + + + do_test_mesh(mapping, tria); + } +} + + +int +main() +{ + initlog(); + + do_test(); + do_test(); + + return 0; +} diff --git a/tests/lac/tensor_product_matrix_08.with_lapack=true.output b/tests/lac/tensor_product_matrix_08.with_lapack=true.output new file mode 100644 index 0000000000..d7f63ad173 --- /dev/null +++ b/tests/lac/tensor_product_matrix_08.with_lapack=true.output @@ -0,0 +1,31 @@ + +DEAL::Storage sizes: 3 128 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 1 128 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 6 256 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 2 256 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 128 128 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 3 128 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 1 128 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 6 256 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 2 256 +DEAL::OK! +DEAL:: +DEAL::Storage sizes: 128 128 +DEAL::OK! +DEAL:: -- 2.39.5