From 818d94dd426cb8e9662bf38334f359cdfe9ae0ab Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Mon, 4 Oct 2021 22:53:14 -0400 Subject: [PATCH] add failing test --- tests/mappings/mapping_q_eulerian_14.cc | 395 ++++++++++++++++++ ...ulerian_14.with_p4est=true.mpirun=1.output | 15 + ...ulerian_14.with_p4est=true.mpirun=3.output | 15 + ...ulerian_14.with_p4est=true.mpirun=4.output | 0 4 files changed, 425 insertions(+) create mode 100644 tests/mappings/mapping_q_eulerian_14.cc create mode 100644 tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=1.output create mode 100644 tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=3.output create mode 100644 tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=4.output diff --git a/tests/mappings/mapping_q_eulerian_14.cc b/tests/mappings/mapping_q_eulerian_14.cc new file mode 100644 index 0000000000..be74c4291b --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_14.cc @@ -0,0 +1,395 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 - 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. +// +// --------------------------------------------------------------------- + + +// test matrix-free Laplace and Mass operators with MappingQEulerian +// on a refined mesh that triggered the following bug with 4 MPI ranks: + +// clang-format off +/* +An error occurred in line <1638> of file <../include/deal.II/lac/la_parallel_vector.h> in function + Number dealii::LinearAlgebra::distributed::Vector::operator()(dealii::LinearAlgebra::distributed::Vector::size_type) const [with Number = double; MemorySpace = dealii::MemorySpace::Host; dealii::LinearAlgebra::distributed::Vector::size_type = unsigned int] +The violated condition was: + partitioner->in_local_range(global_index) || partitioner->ghost_indices().is_element(global_index) +*/ +// clang-format on + +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +const std::vector refined = { + "0_0:", "10_0:", "11_0:", "12_0:", "13_0:", "14_0:", "15_0:", + "1_0:", "21_0:", "23_0:", "26_0:", "27_0:", "28_0:", "29_0:", + "2_0:", "30_0:", "31_0:", "32_0:", "33_0:", "34_0:", "34_1:0", + "34_1:1", "34_1:2", "34_1:3", "35_0:", "35_1:0", "35_1:2", "35_1:3", + "36_0:", "37_0:", "38_0:", "38_1:2", "39_0:", "3_0:", "40_0:", + "40_1:0", "40_1:1", "40_1:2", "40_1:3", "41_0:", "41_1:0", "41_1:1", + "41_1:2", "41_1:3", "42_0:", "42_1:0", "42_1:1", "42_1:2", "42_1:3", + "43_0:", "43_1:0", "43_1:1", "43_1:2", "43_1:3", "44_0:", "44_1:0", + "44_1:1", "44_1:2", "44_1:3", "45_0:", "46_0:", "46_1:0", "46_1:1", + "46_1:2", "46_1:3", "47_0:", "47_1:0", "47_1:2", "48_0:", "49_0:", + "4_0:", "50_0:", "51_0:", "52_0:", "53_0:", "54_0:", "55_0:", + "56_0:", "57_0:", "58_0:", "59_0:", "60_0:", "61_0:", "62_0:", + "63_0:", "64_0:", "65_0:", "66_0:", "67_0:", "68_0:", "69_0:", + "6_0:", "70_0:", "71_0:", "72_0:", "73_0:", "74_0:", "75_0:", + "76_0:", "77_0:", "78_0:", "79_0:", "80_0:", "80_1:0", "80_1:1", + "80_1:2", "80_1:3", "81_0:", "81_1:0", "81_1:1", "81_1:2", "81_1:3", + "82_0:", "82_1:0", "82_1:1", "82_1:2", "82_1:3", "83_0:", "83_1:0", + "83_1:1", "83_1:2", "83_1:3", "84_0:", "84_1:0", "84_1:1", "84_1:2", + "84_1:3", "85_0:", "85_1:0", "86_0:", "86_1:0", "86_1:1", "86_1:2", + "86_1:3", "87_0:", "88_0:", "89_0:", "8_0:", "92_0:", "9_0:"}; + + +// Displacement to be applied +template +class Displacement : public Function +{ +public: + Displacement() + : Function(dim) + {} + + double + value(const Point &p, const unsigned int component) const + { + return (0.1 + 2. * component) + (0.2 + 3. * component) * p[component]; + } + + template + Tensor<1, dim, VectorizedArray> + shift_value(const Point> &p_vec) const + { + Tensor<1, dim, VectorizedArray> shift_vec; + Point p; + for (unsigned int v = 0; v < VectorizedArray::size(); ++v) + { + for (unsigned int d = 0; d < dim; ++d) + p[d] = p_vec[d][v]; + + for (unsigned int d = 0; d < dim; ++d) + shift_vec[d][v] = this->value(p, d); + } + + return shift_vec; + } +}; + + +template +void +test(const unsigned int n_ref = 0) +{ + Displacement displacement_function; + const unsigned int euler_fe_degree = 2; + + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + unsigned int myid = Utilities::MPI::this_mpi_process(mpi_communicator); + unsigned int numproc = Utilities::MPI::n_mpi_processes(mpi_communicator); + + parallel::distributed::Triangulation triangulation( + mpi_communicator, + Triangulation::limit_level_difference_at_vertices, + parallel::distributed::Triangulation::construct_multigrid_hierarchy); + + GridGenerator::subdivided_hyper_cube(triangulation, 10, -5, 5); + while (true) + { + deallog << "refine..." << std::endl; + bool changes = false; + + for (const auto &cell : triangulation.active_cell_iterators()) + + if (!cell->is_artificial() && + std::find(std::begin(refined), + std::end(refined), + cell->id().to_string()) != std::end(refined)) + { + cell->set_refine_flag(); + changes = true; + } + + if (Utilities::MPI::max(changes ? 1 : 0, MPI_COMM_WORLD) == 0) + break; + + triangulation.execute_coarsening_and_refinement(); + } + deallog << "Checksum: " << triangulation.get_checksum() << '\n' + << "n_active_cells: " << triangulation.n_active_cells() << '\n' + << "n_levels: " << triangulation.n_levels() << std::endl; + + FE_Q fe(fe_degree); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + dof_handler.distribute_mg_dofs(); + + // quadrature for MatrixFree, not related to the degree in the + // FE displacement field. + QGauss<1> quadrature_formula(n_q_points); + + + FESystem fe_euler(FE_Q(euler_fe_degree), dim); + DoFHandler dof_handler_euler(triangulation); + dof_handler_euler.distribute_dofs(fe_euler); + dof_handler_euler.distribute_mg_dofs(); + + // IndexSets and constraints + const IndexSet &locally_owned_dofs_euler = + dof_handler_euler.locally_owned_dofs(); + IndexSet locally_relevant_dofs_euler; + DoFTools::extract_locally_relevant_dofs(dof_handler_euler, + locally_relevant_dofs_euler); + + const IndexSet &locally_owned_dofs = dof_handler.locally_owned_dofs(); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + // constraints: + AffineConstraints constraints_euler; + constraints_euler.reinit(locally_relevant_dofs_euler); + DoFTools::make_hanging_node_constraints(dof_handler_euler, constraints_euler); + constraints_euler.close(); + + AffineConstraints constraints; + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints_euler.close(); + + // MG constraints: + MGConstrainedDoFs mg_constrained_dofs_euler; + mg_constrained_dofs_euler.initialize(dof_handler_euler); + + // Displacement vector + LinearAlgebra::distributed::Vector displacement; + displacement.reinit(locally_owned_dofs_euler, + locally_relevant_dofs_euler, + mpi_communicator); + + VectorTools::project, + dim>(dof_handler_euler, + constraints_euler, + QGauss(n_q_points), + displacement_function, + displacement); + displacement.update_ghost_values(); + + MGTransferMatrixFree mg_transfer_euler( + mg_constrained_dofs_euler); + mg_transfer_euler.build(dof_handler_euler); + + // now the core of the test: + const unsigned int max_level = + dof_handler.get_triangulation().n_global_levels() - 1; + const unsigned int min_level = 0; + MGLevelObject> + displacement_level(min_level, max_level); + mg_transfer_euler.interpolate_to_mg(dof_handler_euler, + displacement_level, + displacement); + + // First, check fine-level only: + { + MappingQEulerian> + euler_fine(euler_fe_degree, dof_handler_euler, displacement); + + + + MatrixFree matrix_free_euler; + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + data.mapping_update_flags = update_values | update_gradients | + update_JxW_values | update_quadrature_points; + matrix_free_euler.reinit( + euler_fine, dof_handler, constraints, quadrature_formula, data); + + MatrixFree matrix_free; + matrix_free.reinit(dof_handler, constraints, quadrature_formula, data); + + + // test fine-level mapping: + { + FEEvaluation fe_eval_euler( + matrix_free_euler); + FEEvaluation fe_eval( + matrix_free); + const unsigned int n_cells = matrix_free_euler.n_cell_batches(); + Assert(matrix_free_euler.n_cell_batches() == matrix_free.n_cell_batches(), + ExcInternalError()); + const unsigned int nqp = fe_eval.n_q_points; + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + fe_eval_euler.reinit(cell); + fe_eval.reinit(cell); + for (unsigned int q = 0; q < nqp; ++q) + { + const auto &v1 = fe_eval_euler.quadrature_point(q); + const auto &qp = fe_eval.quadrature_point(q); + const auto v2 = qp + displacement_function.shift_value(qp); + VectorizedArray dist = v1.distance(v2); + for (unsigned int v = 0; v < VectorizedArray::size(); + ++v) + AssertThrow(dist[v] < 1e-8, + ExcMessage("distance: " + std::to_string(dist[v]))); + } + } + } + } + + // now go through all GMG levels: + std::set dirichlet_boundary_ids; + dirichlet_boundary_ids.insert(0); + MGConstrainedDoFs mg_constrained_dofs; + mg_constrained_dofs.initialize(dof_handler); + mg_constrained_dofs.make_zero_boundary_constraints(dof_handler, + dirichlet_boundary_ids); + + for (unsigned int level = min_level; level <= max_level; ++level) + { + typename MatrixFree::AdditionalData + mg_additional_data; + mg_additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::partition_color; + mg_additional_data.mg_level = level; + mg_additional_data.mapping_update_flags = + update_values | update_gradients | update_JxW_values | + update_quadrature_points; + + AffineConstraints level_constraints; + IndexSet relevant_dofs; + DoFTools::extract_locally_relevant_level_dofs(dof_handler, + level, + relevant_dofs); + level_constraints.reinit(relevant_dofs); + level_constraints.add_lines( + mg_constrained_dofs.get_boundary_indices(level)); + level_constraints.close(); + + MatrixFree mg_level_euler; + + MappingQEulerian> + euler_level(euler_fe_degree, + dof_handler_euler, + displacement_level[level], + level); + + mg_level_euler.reinit(euler_level, + dof_handler, + level_constraints, + quadrature_formula, + mg_additional_data); + + MatrixFree mg_level; + mg_level.reinit(dof_handler, + level_constraints, + quadrature_formula, + mg_additional_data); + + // go through all cells and quadrature points: + { + FEEvaluation fe_eval_euler( + mg_level_euler); + FEEvaluation fe_eval( + mg_level); + const unsigned int n_cells = mg_level_euler.n_cell_batches(); + Assert(mg_level_euler.n_cell_batches() == mg_level.n_cell_batches(), + ExcInternalError()); + const unsigned int nqp = fe_eval.n_q_points; + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + fe_eval_euler.reinit(cell); + fe_eval.reinit(cell); + for (unsigned int q = 0; q < nqp; ++q) + { + const auto &v1 = fe_eval_euler.quadrature_point(q); + const auto &qp = fe_eval.quadrature_point(q); + const auto v2 = qp + displacement_function.shift_value(qp); + VectorizedArray dist = v1.distance(v2); + for (unsigned int v = 0; + v < VectorizedArray::size(); + ++v) + AssertThrow(dist[v] < 1e-8, + ExcMessage( + "Level " + std::to_string(level) + + " distance: " + std::to_string(dist[v]))); + } + } + } + } + + deallog << "Ok" << std::endl; +} + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll log; + + test<2>(); +} diff --git a/tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=1.output b/tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..1177cf30be --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=1.output @@ -0,0 +1,15 @@ + +DEAL:2d::refine... +DEAL:2d::refine... +DEAL:2d::refine... +DEAL:2d::Checksum: 4107161474 +n_active_cells: 520 +n_levels: 3 +DEAL:2d::Mass operator: +DEAL:2d::l2_norm: 0.00000 +DEAL:2d::l1_norm: 0.00000 +DEAL:2d::linfty_norm: 0.00000 +DEAL:2d::Laplace operator: +DEAL:2d::l2_norm: 0.00000 +DEAL:2d::l1_norm: 0.00000 +DEAL:2d::linfty_norm: 0.00000 diff --git a/tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=3.output b/tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..15fea44d2a --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=3.output @@ -0,0 +1,15 @@ + +DEAL:2d::refine... +DEAL:2d::refine... +DEAL:2d::refine... +DEAL:2d::Checksum: 4107161474 +n_active_cells: 307 +n_levels: 3 +DEAL:2d::Mass operator: +DEAL:2d::l2_norm: 0.00000 +DEAL:2d::l1_norm: 0.00000 +DEAL:2d::linfty_norm: 0.00000 +DEAL:2d::Laplace operator: +DEAL:2d::l2_norm: 0.00000 +DEAL:2d::l1_norm: 0.00000 +DEAL:2d::linfty_norm: 0.00000 diff --git a/tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=4.output b/tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=4.output new file mode 100644 index 0000000000..e69de29bb2 -- 2.39.5