From 719bbf3f2c6d8f8a80d464308d194244183a62cd Mon Sep 17 00:00:00 2001 From: Svenja Schoeder Date: Tue, 21 Aug 2018 14:24:37 +0200 Subject: [PATCH] bugfix for braces, added test for vectorization mask --- include/deal.II/matrix_free/fe_evaluation.h | 84 +++-- .../advect_1d_vectorization_mask.cc | 353 ++++++++++++++++++ .../advect_1d_vectorization_mask.output | 22 ++ 3 files changed, 422 insertions(+), 37 deletions(-) create mode 100644 tests/matrix_free/advect_1d_vectorization_mask.cc create mode 100644 tests/matrix_free/advect_1d_vectorization_mask.output diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 46c90653bc..7b509d085e 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -3986,12 +3986,12 @@ FEEvaluationBase::read_write_operation( // Case 3: standard operation with one index per degree of freedom -> go on // here - constexpr unsigned int n_vectorization = VectorizedArray::n_array_elements; Assert(mask.count() == n_vectorization, ExcNotImplemented("Masking currently not implemented for " "non-contiguous DoF storage")); + const unsigned int dofs_per_component = this->data->dofs_per_component_on_cell; if (dof_info->index_storage_variants @@ -4572,24 +4572,29 @@ FEEvaluationBase:: contiguous) { if (n_components == 1 || n_fe_components == 1) - for (unsigned int v = 0; v < vectorization_populated; ++v) - if (mask[v] == true) - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; - ++i) - operation.process_dof(dof_indices[v] + i, - *src[comp], - values_dofs[comp][i][v]); - else - for (unsigned int v = 0; v < vectorization_populated; ++v) - if (mask[v] == true) - for (unsigned int i = 0; - i < data->dofs_per_component_on_cell; - ++i) - operation.process_dof( - dof_indices[v] + i + - comp * data->dofs_per_component_on_cell, - *src[0], - values_dofs[comp][i][v]); + { + for (unsigned int v = 0; v < vectorization_populated; ++v) + if (mask[v] == true) + for (unsigned int i = 0; + i < data->dofs_per_component_on_cell; + ++i) + operation.process_dof(dof_indices[v] + i, + *src[comp], + values_dofs[comp][i][v]); + } + else + { + for (unsigned int v = 0; v < vectorization_populated; ++v) + if (mask[v] == true) + for (unsigned int i = 0; + i < data->dofs_per_component_on_cell; + ++i) + operation.process_dof( + dof_indices[v] + i + + comp * data->dofs_per_component_on_cell, + *src[0], + values_dofs[comp][i][v]); + } } else { @@ -4601,24 +4606,29 @@ FEEvaluationBase:: VectorizedArray::n_array_elements + 1); if (n_components == 1 || n_fe_components == 1) for (unsigned int v = 0; v < vectorization_populated; ++v) - if (mask[v] == true) - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; - ++i) - operation.process_dof(dof_indices[v] + i * offsets[v], - *src[comp], - values_dofs[comp][i][v]); - else - for (unsigned int v = 0; v < vectorization_populated; ++v) - if (mask[v] == true) - for (unsigned int i = 0; - i < data->dofs_per_component_on_cell; - ++i) - operation.process_dof( - dof_indices[v] + - (i + comp * data->dofs_per_component_on_cell) * - offsets[v], - *src[0], - values_dofs[comp][i][v]); + { + if (mask[v] == true) + for (unsigned int i = 0; + i < data->dofs_per_component_on_cell; + ++i) + operation.process_dof(dof_indices[v] + i * offsets[v], + *src[comp], + values_dofs[comp][i][v]); + } + else + { + for (unsigned int v = 0; v < vectorization_populated; ++v) + if (mask[v] == true) + for (unsigned int i = 0; + i < data->dofs_per_component_on_cell; + ++i) + operation.process_dof( + dof_indices[v] + + (i + comp * data->dofs_per_component_on_cell) * + offsets[v], + *src[0], + values_dofs[comp][i][v]); + } } } } diff --git a/tests/matrix_free/advect_1d_vectorization_mask.cc b/tests/matrix_free/advect_1d_vectorization_mask.cc new file mode 100644 index 0000000000..4706713b2f --- /dev/null +++ b/tests/matrix_free/advect_1d_vectorization_mask.cc @@ -0,0 +1,353 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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 advect_1d but testing the mask for distribute_local_to_global + +#include +#include + +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +std::ofstream logfile("output"); + + +template , + int n_components = 1> +class MatrixFreeAdvectionBasic +{ +public: + MatrixFreeAdvectionBasic(const MatrixFree &data, + const bool zero_within_loop = true, + const unsigned int start_vector_component = 0) + : data(data) + , zero_within_loop(zero_within_loop) + , start_vector_component(start_vector_component) + { + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + if (!zero_within_loop) + dst = 0; + data.loop(&MatrixFreeAdvectionBasic::local_apply, + &MatrixFreeAdvectionBasic::local_apply_face, + &MatrixFreeAdvectionBasic::local_apply_boundary_face, + this, + dst, + src, + zero_within_loop, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + + FEEvaluation phi(data); + + const unsigned int dofs_per_cell = phi.dofs_per_cell; + + AlignedVector> coefficients(phi.dofs_per_cell); + MatrixFreeOperators:: + CellwiseInverseMassMatrix + inverse(phi); + + for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell) + { + phi.reinit(cell); + phi.read_dof_values(dst); + + inverse.fill_inverse_JxW_values(coefficients); + inverse.apply(coefficients, + n_components, + phi.begin_dof_values(), + phi.begin_dof_values()); + + phi.set_dof_values(dst); + } + } + +private: + void + local_apply(const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &cell_range) const + { + FEEvaluation phi( + data, 0, 0, start_vector_component); + + const unsigned int n_vect = VectorizedArray::n_array_elements; + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(advection * phi.get_value(q), q); + phi.integrate(false, true); + for (unsigned int v = 0; v < n_vect; ++v) + { + std::bitset mask; + mask[v] = true; + phi.distribute_local_to_global(dst, 0, mask); + } + } + } + + void + local_apply_face( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &face_range) const + { + FEFaceEvaluation phi_m( + data, true, 0, 0, start_vector_component); + FEFaceEvaluation phi_p( + data, false, 0, 0, start_vector_component); + typedef typename FEFaceEvaluation::value_type value_type; + + const unsigned int n_vect = VectorizedArray::n_array_elements; + + for (unsigned int face = face_range.first; face < face_range.second; face++) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + phi_m.evaluate(true, false); + phi_p.reinit(face); + phi_p.read_dof_values(src); + phi_p.evaluate(true, false); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + value_type u_minus = phi_m.get_value(q), + u_plus = phi_p.get_value(q); + const VectorizedArray normal_times_advection = + advection * phi_m.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + phi_m.submit_value(-flux_times_normal, q); + phi_p.submit_value(flux_times_normal, q); + } + + phi_m.integrate(true, false); + for (unsigned int v = 0; v < n_vect; ++v) + { + std::bitset mask; + mask[v] = true; + phi_m.distribute_local_to_global(dst, 0, mask); + } + phi_p.integrate(true, false); + for (unsigned int v = 0; v < n_vect; ++v) + { + std::bitset mask; + mask[v] = true; + phi_p.distribute_local_to_global(dst, 0, mask); + } + } + } + + void + local_apply_boundary_face( + const MatrixFree & data, + VectorType & dst, + const VectorType & src, + const std::pair &face_range) const + { + FEFaceEvaluation + fe_eval(data, true, 0, 0, start_vector_component); + typedef typename FEFaceEvaluation::value_type value_type; + + const unsigned int n_vect = VectorizedArray::n_array_elements; + + for (unsigned int face = face_range.first; face < face_range.second; face++) + { + fe_eval.reinit(face); + fe_eval.read_dof_values(src); + fe_eval.evaluate(true, false); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + value_type u_minus = fe_eval.get_value(q); + value_type u_plus = -u_minus; + const VectorizedArray normal_times_advection = + advection * fe_eval.get_normal_vector(q); + const value_type flux_times_normal = + make_vectorized_array(0.5) * + ((u_minus + u_plus) * normal_times_advection + + std::abs(normal_times_advection) * (u_minus - u_plus)); + fe_eval.submit_value(-flux_times_normal, q); + } + + fe_eval.integrate(true, false); + for (unsigned int v = 0; v < n_vect; ++v) + { + std::bitset mask = + std::bitset::n_array_elements>(); + mask[v] = true; + fe_eval.distribute_local_to_global(dst, 0, mask); + } + } + } + + const MatrixFree & data; + const bool zero_within_loop; + const unsigned int start_vector_component; + Tensor<1, dim, VectorizedArray> advection; +}; + + + +template +class AnalyticFunction : public Function +{ +public: + static_assert(dim == 1, "Only 1D implemented"); + AnalyticFunction() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int) const override + { + return std::sin(3 * numbers::PI * p[0] / 0.8); + } +}; + + + +template +class AnalyticDerivative : public Function +{ +public: + static_assert(dim == 1, "Only 1D implemented"); + AnalyticDerivative() + : Function(1) + {} + + virtual double + value(const Point &p, const unsigned int) const override + { + Tensor<1, dim> advection; + for (unsigned int d = 0; d < dim; ++d) + advection[d] = 0.4 + 0.12 * d; + + return -std::cos(3 * numbers::PI * p[0] / 0.8) * advection[0] * 3 * + numbers::PI / 0.8; + } +}; + + + +template +void +test(const unsigned int n_refine) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria, 0, 0.8); + tria.refine_global(n_refine); + + FE_DGQ fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + constraints.close(); + + if (n_refine == 3) + { + deallog << "Testing " << dof.get_fe().get_name(); + deallog << std::endl; + } + + LinearAlgebra::distributed::Vector in, out; + + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + + MatrixFree mf_data; + mf_data.reinit(dof, constraints, quad, data); + + mf_data.initialize_dof_vector(in); + mf_data.initialize_dof_vector(out); + + VectorTools::interpolate(dof, AnalyticFunction(), in); + + MatrixFreeAdvectionBasic> + mf2(mf_data); + mf2.vmult(out, in); + + VectorTools::interpolate(dof, AnalyticDerivative(), in); + out -= in; + + double diff_norm = out.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << " " << std::endl; +} + + + +int +main() +{ + initlog(); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 2>(r); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 4>(r); + + for (unsigned int r = 3; r < 9; ++r) + test<1, 5>(r); +} diff --git a/tests/matrix_free/advect_1d_vectorization_mask.output b/tests/matrix_free/advect_1d_vectorization_mask.output new file mode 100644 index 0000000000..b213520047 --- /dev/null +++ b/tests/matrix_free/advect_1d_vectorization_mask.output @@ -0,0 +1,22 @@ + +DEAL::Testing FE_DGQ<1>(2) +DEAL::Norm of difference: 0.531539 +DEAL::Norm of difference: 0.135407 +DEAL::Norm of difference: 0.0340112 +DEAL::Norm of difference: 0.00851279 +DEAL::Norm of difference: 0.00212882 +DEAL::Norm of difference: 0.000532245 +DEAL::Testing FE_DGQ<1>(4) +DEAL::Norm of difference: 0.00529591 +DEAL::Norm of difference: 0.000336014 +DEAL::Norm of difference: 2.10801e-05 +DEAL::Norm of difference: 1.31874e-06 +DEAL::Norm of difference: 8.24413e-08 +DEAL::Norm of difference: 5.15331e-09 +DEAL::Testing FE_DGQ<1>(5) +DEAL::Norm of difference: 0.000346887 +DEAL::Norm of difference: 1.09982e-05 +DEAL::Norm of difference: 3.44938e-07 +DEAL::Norm of difference: 1.07881e-08 +DEAL::Norm of difference: 3.35085e-10 +DEAL::Norm of difference: 9.72689e-12 -- 2.39.5