]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add failing test
authorTimo Heister <timo.heister@gmail.com>
Tue, 5 Oct 2021 02:53:14 +0000 (22:53 -0400)
committerTimo Heister <timo.heister@gmail.com>
Tue, 5 Oct 2021 02:53:14 +0000 (22:53 -0400)
tests/mappings/mapping_q_eulerian_14.cc [new file with mode: 0644]
tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=3.output [new file with mode: 0644]
tests/mappings/mapping_q_eulerian_14.with_p4est=true.mpirun=4.output [new file with mode: 0644]

diff --git a/tests/mappings/mapping_q_eulerian_14.cc b/tests/mappings/mapping_q_eulerian_14.cc
new file mode 100644 (file)
index 0000000..be74c42
--- /dev/null
@@ -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<Number, MemorySpace>::operator()(dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>::size_type) const [with Number = double; MemorySpace = dealii::MemorySpace::Host; dealii::LinearAlgebra::distributed::Vector<Number, MemorySpace>::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 <deal.II/base/function.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/thread_management.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/identity_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/operators.h>
+
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+
+const std::vector<std::string> 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 <int dim>
+class Displacement : public Function<dim>
+{
+public:
+  Displacement()
+    : Function<dim>(dim)
+  {}
+
+  double
+  value(const Point<dim> &p, const unsigned int component) const
+  {
+    return (0.1 + 2. * component) + (0.2 + 3. * component) * p[component];
+  }
+
+  template <typename NumberType>
+  Tensor<1, dim, VectorizedArray<NumberType>>
+  shift_value(const Point<dim, VectorizedArray<NumberType>> &p_vec) const
+  {
+    Tensor<1, dim, VectorizedArray<NumberType>> shift_vec;
+    Point<dim>                                  p;
+    for (unsigned int v = 0; v < VectorizedArray<NumberType>::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 <int dim,
+          int fe_degree            = 2,
+          int n_q_points           = fe_degree + 1,
+          typename NumberType      = double,
+          typename LevelNumberType = NumberType>
+void
+test(const unsigned int n_ref = 0)
+{
+  Displacement<dim>  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<dim> triangulation(
+    mpi_communicator,
+    Triangulation<dim>::limit_level_difference_at_vertices,
+    parallel::distributed::Triangulation<dim>::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<dim>       fe(fe_degree);
+  DoFHandler<dim> 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<dim>   fe_euler(FE_Q<dim>(euler_fe_degree), dim);
+  DoFHandler<dim> 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<double> constraints_euler;
+  constraints_euler.reinit(locally_relevant_dofs_euler);
+  DoFTools::make_hanging_node_constraints(dof_handler_euler, constraints_euler);
+  constraints_euler.close();
+
+  AffineConstraints<double> 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<NumberType> displacement;
+  displacement.reinit(locally_owned_dofs_euler,
+                      locally_relevant_dofs_euler,
+                      mpi_communicator);
+
+  VectorTools::project<dim,
+                       LinearAlgebra::distributed::Vector<NumberType>,
+                       dim>(dof_handler_euler,
+                            constraints_euler,
+                            QGauss<dim>(n_q_points),
+                            displacement_function,
+                            displacement);
+  displacement.update_ghost_values();
+
+  MGTransferMatrixFree<dim, LevelNumberType> 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<LinearAlgebra::distributed::Vector<LevelNumberType>>
+    displacement_level(min_level, max_level);
+  mg_transfer_euler.interpolate_to_mg(dof_handler_euler,
+                                      displacement_level,
+                                      displacement);
+
+  // First, check fine-level only:
+  {
+    MappingQEulerian<dim, LinearAlgebra::distributed::Vector<NumberType>>
+      euler_fine(euler_fe_degree, dof_handler_euler, displacement);
+
+
+
+    MatrixFree<dim, NumberType>                          matrix_free_euler;
+    typename MatrixFree<dim, NumberType>::AdditionalData data;
+    data.tasks_parallel_scheme =
+      MatrixFree<dim, NumberType>::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<dim, NumberType> matrix_free;
+    matrix_free.reinit(dof_handler, constraints, quadrature_formula, data);
+
+
+    // test fine-level mapping:
+    {
+      FEEvaluation<dim, fe_degree, n_q_points, 1, NumberType> fe_eval_euler(
+        matrix_free_euler);
+      FEEvaluation<dim, fe_degree, n_q_points, 1, NumberType> 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<NumberType> dist = v1.distance(v2);
+              for (unsigned int v = 0; v < VectorizedArray<NumberType>::size();
+                   ++v)
+                AssertThrow(dist[v] < 1e-8,
+                            ExcMessage("distance: " + std::to_string(dist[v])));
+            }
+        }
+    }
+  }
+
+  // now go through all GMG levels:
+  std::set<types::boundary_id> 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<dim, LevelNumberType>::AdditionalData
+        mg_additional_data;
+      mg_additional_data.tasks_parallel_scheme =
+        MatrixFree<dim, LevelNumberType>::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<double> 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<dim, LevelNumberType> mg_level_euler;
+
+      MappingQEulerian<dim, LinearAlgebra::distributed::Vector<LevelNumberType>>
+        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<dim, LevelNumberType> mg_level;
+      mg_level.reinit(dof_handler,
+                      level_constraints,
+                      quadrature_formula,
+                      mg_additional_data);
+
+      // go through all cells and quadrature points:
+      {
+        FEEvaluation<dim, fe_degree, n_q_points, 1, NumberType> fe_eval_euler(
+          mg_level_euler);
+        FEEvaluation<dim, fe_degree, n_q_points, 1, NumberType> 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<NumberType> dist = v1.distance(v2);
+                for (unsigned int v = 0;
+                     v < VectorizedArray<NumberType>::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 (file)
index 0000000..1177cf3
--- /dev/null
@@ -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 (file)
index 0000000..15fea44
--- /dev/null
@@ -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 (file)
index 0000000..e69de29

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.