From: Peter Munch Date: Tue, 20 Oct 2020 11:15:16 +0000 (+0200) Subject: Fix DoFInfo::assign_ghosts for hp in parallel X-Git-Tag: v9.3.0-rc1~991^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=96598756ce9aa098acdc881e68a689621ee90766;p=dealii.git Fix DoFInfo::assign_ghosts for hp in parallel --- diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index f8d115f6d8..c94c42819c 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -473,7 +473,7 @@ namespace internal (cell_active_fe_index.size() == 0 || dofs_per_cell.size() == 1) ? 0 : - cell_active_fe_index[i]; + cell_active_fe_index[boundary_cells[i]]; AssertIndexRange(fe_index, dofs_per_cell.size()); const unsigned int *row_end = data_ptr + dofs_per_cell[fe_index]; diff --git a/tests/matrix_free/matrix_vector_hp_parallel_01.cc b/tests/matrix_free/matrix_vector_hp_parallel_01.cc new file mode 100644 index 0000000000..53f81eb8bd --- /dev/null +++ b/tests/matrix_free/matrix_vector_hp_parallel_01.cc @@ -0,0 +1,293 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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 matrix_vector_hp.cc but in parallel. + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_mf.h" + + + +template +class MatrixFreeTestHP +{ +public: + MatrixFreeTestHP(const MatrixFree &data_in) + : data(data_in){}; + + void + local_apply(const MatrixFree & data, + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const std::pair &cell_range) const + { + // ask MatrixFree for cell_range for different orders + std::pair subrange_deg = + data.create_cell_subrange_hp(cell_range, 1); + if (subrange_deg.second > subrange_deg.first) + helmholtz_operator, 2>( + data, dst, src, subrange_deg); + subrange_deg = data.create_cell_subrange_hp(cell_range, 2); + if (subrange_deg.second > subrange_deg.first) + helmholtz_operator, 3>( + data, dst, src, subrange_deg); + subrange_deg = data.create_cell_subrange_hp(cell_range, 3); + if (subrange_deg.second > subrange_deg.first) + helmholtz_operator, 4>( + data, dst, src, subrange_deg); + subrange_deg = data.create_cell_subrange_hp(cell_range, 4); + if (subrange_deg.second > subrange_deg.first) + helmholtz_operator, 5>( + data, dst, src, subrange_deg); + subrange_deg = data.create_cell_subrange_hp(cell_range, 5); + if (subrange_deg.second > subrange_deg.first) + helmholtz_operator, 6>( + data, dst, src, subrange_deg); + subrange_deg = data.create_cell_subrange_hp(cell_range, 6); + if (subrange_deg.second > subrange_deg.first) + helmholtz_operator, 7>( + data, dst, src, subrange_deg); + subrange_deg = data.create_cell_subrange_hp(cell_range, 7); + if (subrange_deg.second > subrange_deg.first) + helmholtz_operator, 8>( + data, dst, src, subrange_deg); + } + + void + vmult(LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src) const + { + dst = 0; + data.cell_loop(&MatrixFreeTestHP::local_apply, this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + if (fe_degree > 1) + return; + + typedef double number; + const SphericalManifold manifold; + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_ball(tria); + typename Triangulation::active_cell_iterator cell = tria.begin_active(), + endc = tria.end(); + for (; cell != endc; ++cell) + for (const unsigned int f : GeometryInfo::face_indices()) + if (cell->at_boundary(f)) + cell->face(f)->set_all_manifold_ids(0); + tria.set_manifold(0, manifold); + tria.refine_global(1); + + // refine a few cells + for (unsigned int i = 0; i < 11 - 3 * dim; ++i) + { + typename Triangulation::active_cell_iterator cell = + tria.begin_active(), + endc = tria.end(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (cell->is_locally_owned() && (counter % (7 - i) == 0)) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + const unsigned int max_degree = 9 - 2 * dim; + + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + hp::QCollection<1> quadrature_collection_mf; + + for (unsigned int deg = 1; deg <= max_degree; ++deg) + { + fe_collection.push_back(FE_Q(QGaussLobatto<1>(deg + 1))); + quadrature_collection.push_back(QGauss(deg + 1)); + quadrature_collection_mf.push_back(QGauss<1>(deg + 1)); + } + + hp::DoFHandler dof(tria); + // set the active FE index in a random order + { + typename hp::DoFHandler::active_cell_iterator cell = + dof.begin_active(), + endc = dof.end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned() == false) + continue; + + const unsigned int fe_index = Testing::rand() % max_degree; + cell->set_active_fe_index(fe_index); + } + } + + // setup DoFs + dof.distribute_dofs(fe_collection); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values(dof, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + TrilinosWrappers::SparsityPattern dsp(dof.locally_owned_dofs(), + MPI_COMM_WORLD); + DoFTools::make_sparsity_pattern(dof, dsp, constraints, false); + dsp.compress(); + + TrilinosWrappers::SparseMatrix system_matrix; + system_matrix.reinit(dsp); + + // std::cout << "Number of cells: " << + // dof.get_triangulation().n_active_cells() << std::endl; std::cout << "Number + // of degrees of freedom: " << dof.n_dofs() << std::endl; std::cout << "Number + // of constraints: " << constraints.n_constraints() << std::endl; + + // set up MatrixFree + MatrixFree mf_data; + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + mf_data.reinit(dof, constraints, quadrature_collection_mf, data); + MatrixFreeTestHP mf(mf_data); + + // assemble sparse matrix with (\nabla v, + // \nabla u) + (v, 10 * u) + { + hp::FEValues hp_fe_values(fe_collection, + quadrature_collection, + update_values | update_gradients | + update_JxW_values); + FullMatrix cell_matrix; + std::vector local_dof_indices; + + typename hp::DoFHandler::active_cell_iterator cell = + dof.begin_active(), + endc = dof.end(); + for (; cell != endc; ++cell) + { + if (cell->is_locally_owned() == false) + continue; + + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_matrix = 0; + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_points; + ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += ((fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) + + 10. * fe_values.shape_value(i, q_point) * + fe_values.shape_value(j, q_point)) * + fe_values.JxW(q_point)); + } + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + constraints.distribute_local_to_global(cell_matrix, + local_dof_indices, + system_matrix); + } + } + system_matrix.compress(VectorOperation::values::add); + + // fill a right hand side vector with random + // numbers in unconstrained degrees of freedom + LinearAlgebra::distributed::Vector src; + mf_data.initialize_dof_vector(src); + + + LinearAlgebra::distributed::Vector result_spmv(src), result_mf(src); + + for (unsigned int i = 0; i < dof.n_locally_owned_dofs(); ++i) + { + if (constraints.is_constrained( + src.get_partitioner()->local_to_global(i)) == false) + src.local_element(i) = random_value(); + } + + // now perform matrix-vector product and check + // its correctness + system_matrix.vmult(result_spmv, src); + mf.vmult(result_mf, src); + + result_mf -= result_spmv; + const double diff_norm = result_mf.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + mpi_initlog(); + + { + deallog.push("2d"); + test<2, 1>(); + test<2, 2>(); + deallog.pop(); + deallog.push("3d"); + test<3, 1>(); + test<3, 2>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=1.with_trilinos=true.output b/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=1.with_trilinos=true.output new file mode 100644 index 0000000000..7252835322 --- /dev/null +++ b/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=1.with_trilinos=true.output @@ -0,0 +1,5 @@ + +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=5.with_trilinos=true.output b/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=5.with_trilinos=true.output new file mode 100644 index 0000000000..7252835322 --- /dev/null +++ b/tests/matrix_free/matrix_vector_hp_parallel_01.with_mpi=true.with_p4est=true.mpirun=5.with_trilinos=true.output @@ -0,0 +1,5 @@ + +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Norm of difference: 0 +DEAL:3d::