From 6351f255d16c3cf465752e2ce943726df6bf4ea4 Mon Sep 17 00:00:00 2001 From: Sebastian Proell Date: Mon, 27 Jun 2022 16:38:15 +0200 Subject: [PATCH] Fix some VectorizedArrayTypes for non-default vectorization --- include/deal.II/matrix_free/fe_evaluation.h | 18 +- .../matrix_vector_stokes_no_vectorization.cc | 338 ++++++++++++++++++ ...trix_vector_stokes_no_vectorization.output | 16 + 3 files changed, 362 insertions(+), 10 deletions(-) create mode 100644 tests/matrix_free/matrix_vector_stokes_no_vectorization.cc create mode 100644 tests/matrix_free/matrix_vector_stokes_no_vectorization.output diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 01362337cb..ed5d569267 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -5747,11 +5747,10 @@ FEEvaluationAccess::get_value( this->cell_type == internal::MatrixFreeFunctions::cartesian) { // Cartesian cell - const Tensor<2, dim, dealii::VectorizedArray> jac = - this->jacobian[1]; - const VectorizedArrayType inv_det = + const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1]; + const VectorizedArrayType inv_det = (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] : - this->jacobian[0][0][0] * this->jacobian[0][1][1] * + this->jacobian[0][0][0] * this->jacobian[0][1][1] * this->jacobian[0][2][2]; // J * u * det(J^-1) @@ -5762,7 +5761,7 @@ FEEvaluationAccess::get_value( else { // Affine or general cell - const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = (this->cell_type > internal::MatrixFreeFunctions::affine) ? this->jacobian[q_point] : this->jacobian[0]; @@ -5838,7 +5837,7 @@ FEEvaluationAccess:: else if (this->cell_type <= internal::MatrixFreeFunctions::affine) { // Affine cell - const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = this->jacobian[0]; const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; @@ -6081,8 +6080,7 @@ FEEvaluationAccess:: if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { - const Tensor<2, dim, dealii::VectorizedArray> jac = - this->jacobian[1]; + const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1]; const VectorizedArrayType weight = this->quadrature_weights[q_point]; for (unsigned int comp = 0; comp < n_components; ++comp) @@ -6092,7 +6090,7 @@ FEEvaluationAccess:: else { // Affine or general cell - const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = (this->cell_type > internal::MatrixFreeFunctions::affine) ? this->jacobian[q_point] : this->jacobian[0]; @@ -6173,7 +6171,7 @@ FEEvaluationAccess:: else if (this->cell_type <= internal::MatrixFreeFunctions::affine) { // Affine cell - const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = this->jacobian[0]; const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; diff --git a/tests/matrix_free/matrix_vector_stokes_no_vectorization.cc b/tests/matrix_free/matrix_vector_stokes_no_vectorization.cc new file mode 100644 index 0000000000..2223cb2784 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_no_vectorization.cc @@ -0,0 +1,338 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2021 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. +// +// --------------------------------------------------------------------- + + +// same as matrix_vector_stokes but explicitly test VectorizedArray + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include "../tests.h" + +#include "create_mesh.h" + + + +template +class MatrixFreeTest +{ +public: + using CellIterator = typename DoFHandler::active_cell_iterator; + using Number = double; + using VectorizedArrayType = VectorizedArray; + + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in){}; + + void + local_apply(const MatrixFree &data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + using vector_t = VectorizedArrayType; + FEEvaluation + velocity(data, 0); + FEEvaluation + pressure(data, 1); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + velocity.reinit(cell); + velocity.read_dof_values(src, 0); + velocity.evaluate(EvaluationFlags::gradients); + pressure.reinit(cell); + pressure.read_dof_values(src, dim); + pressure.evaluate(EvaluationFlags::values); + + for (unsigned int q = 0; q < velocity.n_q_points; ++q) + { + Tensor<2, dim, vector_t> grad_u = velocity.get_gradient(q); + vector_t pres = pressure.get_value(q); + vector_t div = -trace(grad_u); + pressure.submit_value(div, q); + + // subtract p * I + for (unsigned int d = 0; d < dim; ++d) + grad_u[d][d] -= pres; + + velocity.submit_gradient(grad_u, q); + } + + velocity.integrate(EvaluationFlags::gradients); + velocity.distribute_local_to_global(dst, 0); + pressure.integrate(EvaluationFlags::values); + pressure.distribute_local_to_global(dst, dim); + } + } + + + void + vmult(VectorType &dst, const VectorType &src) const + { + AssertDimension(dst.size(), dim + 1); + for (unsigned int d = 0; d < dim + 1; ++d) + dst[d] = 0; + data.cell_loop(&MatrixFreeTest::local_apply, + this, + dst, + src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + Triangulation triangulation; + create_mesh(triangulation); + if (fe_degree == 1) + triangulation.refine_global(4 - dim); + else + triangulation.refine_global(3 - dim); + + FE_Q fe_u(fe_degree + 1); + FE_Q fe_p(fe_degree); + FESystem fe(fe_u, dim, fe_p, 1); + DoFHandler dof_handler_u(triangulation); + DoFHandler dof_handler_p(triangulation); + DoFHandler dof_handler(triangulation); + + using MatrixFreeType = MatrixFree>; + MatrixFreeType mf_data; + + AffineConstraints constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + std::vector> vec1, vec2; + + dof_handler.distribute_dofs(fe); + dof_handler_u.distribute_dofs(fe_u); + dof_handler_p.distribute_dofs(fe_p); + DoFRenumbering::component_wise(dof_handler); + + constraints.close(); + + const std::vector dofs_per_block = + DoFTools::count_dofs_per_fe_component(dof_handler); + + // std::cout << " Number of active cells: " + // << triangulation.n_active_cells() + // << std::endl + // << " Number of degrees of freedom: " + // << dof_handler.n_dofs() + // << " (" << n_u << '+' << n_p << ')' + // << std::endl; + + { + BlockDynamicSparsityPattern csp(dim + 1, dim + 1); + + for (unsigned int d = 0; d < dim + 1; ++d) + for (unsigned int e = 0; e < dim + 1; ++e) + csp.block(d, e).reinit(dofs_per_block[d], dofs_per_block[e]); + + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern(dof_handler, csp, constraints, false); + sparsity_pattern.copy_from(csp); + } + + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dim + 1); + for (unsigned int i = 0; i < dim + 1; ++i) + solution.block(i).reinit(dofs_per_block[i]); + solution.collect_sizes(); + + system_rhs.reinit(solution); + + vec1.resize(dim + 1); + vec2.resize(dim + 1); + vec1[0].reinit(dofs_per_block[0]); + vec2[0].reinit(vec1[0]); + for (unsigned int i = 1; i < dim; ++i) + { + vec1[i].reinit(vec1[0]); + vec2[i].reinit(vec1[0]); + } + vec1[dim].reinit(dofs_per_block[dim]); + vec2[dim].reinit(vec1[dim]); + + // this is from step-22 + { + QGauss quadrature_formula(fe_degree + 2); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_JxW_values | + update_gradients); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + std::vector> phi_grad_u(dofs_per_cell); + std::vector div_phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + fe_values.reinit(cell); + local_matrix = 0; + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + phi_grad_u[k] = fe_values[velocities].gradient(k, q); + div_phi_u[k] = fe_values[velocities].divergence(k, q); + phi_p[k] = fe_values[pressure].value(k, q); + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j <= i; ++j) + { + local_matrix(i, j) += + (scalar_product(phi_grad_u[i], phi_grad_u[j]) - + div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) * + fe_values.JxW(q); + } + } + } + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + local_matrix(i, j) = local_matrix(j, i); + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(local_matrix, + local_dof_indices, + system_matrix); + } + } + + // first system_rhs with random numbers + for (unsigned int i = 0; i < dim + 1; ++i) + for (unsigned int j = 0; j < system_rhs.block(i).size(); ++j) + { + const double val = -1. + 2. * random_value(); + system_rhs.block(i)(j) = val; + vec1[i](j) = val; + } + + // setup matrix-free structure + { + std::vector *> dofs; + dofs.push_back(&dof_handler_u); + dofs.push_back(&dof_handler_p); + AffineConstraints dummy_constraints; + dummy_constraints.close(); + std::vector *> constraints; + constraints.push_back(&dummy_constraints); + constraints.push_back(&dummy_constraints); + QGauss<1> quad(fe_degree + 2); + mf_data.reinit(MappingQ1{}, + dofs, + constraints, + quad, + typename MatrixFreeType::AdditionalData( + MatrixFreeType::AdditionalData::none)); + } + + system_matrix.vmult(solution, system_rhs); + + using VectorType = std::vector>; + MatrixFreeTest mf(mf_data); + mf.vmult(vec2, vec1); + + // Verification + double error = 0.; + for (unsigned int i = 0; i < dim + 1; ++i) + for (unsigned int j = 0; j < system_rhs.block(i).size(); ++j) + error += std::fabs(solution.block(i)(j) - vec2[i](j)); + double relative = solution.block(0).l1_norm(); + deallog << " Verification fe degree " << fe_degree << ": " + << error / relative << std::endl + << std::endl; +} + + + +int +main() +{ + initlog(); + deallog << std::setprecision(3); + + { + deallog << std::endl << "Test with doubles" << std::endl << std::endl; + deallog.push("2d"); + test<2, 1>(); + test<2, 2>(); + test<2, 3>(); + test<2, 4>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + test<3, 2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_stokes_no_vectorization.output b/tests/matrix_free/matrix_vector_stokes_no_vectorization.output new file mode 100644 index 0000000000..fca0544652 --- /dev/null +++ b/tests/matrix_free/matrix_vector_stokes_no_vectorization.output @@ -0,0 +1,16 @@ + +DEAL:: +DEAL::Test with doubles +DEAL:: +DEAL:2d:: Verification fe degree 1: 0 +DEAL:2d:: +DEAL:2d:: Verification fe degree 2: 0 +DEAL:2d:: +DEAL:2d:: Verification fe degree 3: 0 +DEAL:2d:: +DEAL:2d:: Verification fe degree 4: 0 +DEAL:2d:: +DEAL:3d:: Verification fe degree 1: 0 +DEAL:3d:: +DEAL:3d:: Verification fe degree 2: 0 +DEAL:3d:: -- 2.39.5