Add test

diff --git a/tests/matrix_free/loop_boundary.cc b/tests/matrix_free/loop_boundary.cc
new file mode 100644 (file)
index 0000000..a06429f
--- /dev/null
+++ b/tests/matrix_free/loop_boundary.cc
@@ -0,0 +1,235 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the matrix-free loop with continuous
+// elements when both a cell and a boundary function are given. In an initial
+// implementation, we used to miss to exchange some ghost entries upon the
+// compress() functionality within the matrix-free loop
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+template <int dim, int fe_degree>
+void
+do_test(const unsigned int n_refine, const bool overlap_communication)
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria, 0.1, 0.96);
+  tria.refine_global(n_refine);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  MatrixFree<dim>                          matrix_free;
+  typename MatrixFree<dim>::AdditionalData add_data;
+  add_data.mapping_update_flags = update_values | update_gradients |
+                                  update_JxW_values | update_quadrature_points;
+  add_data.mapping_update_flags_boundary_faces =
+    update_values | update_JxW_values | update_quadrature_points;
+  add_data.overlap_communication_computation = overlap_communication;
+  matrix_free.reinit(dof_handler,
+                     constraints,
+                     QGauss<1>(fe.degree + 1),
+                     add_data);
+
+  LinearAlgebra::distributed::Vector<double> in, ref, test;
+  matrix_free.initialize_dof_vector(in);
+  for (unsigned int i = 0; i < in.get_partitioner()->local_size(); ++i)
+    in.local_element(i) = in.get_partitioner()->local_to_global(i);
+  matrix_free.initialize_dof_vector(ref);
+  matrix_free.initialize_dof_vector(test);
+
+  std::function<void(const MatrixFree<dim> &,
+                     LinearAlgebra::distributed::Vector<double> &,
+                     const LinearAlgebra::distributed::Vector<double> &,
+                     const std::pair<unsigned int, unsigned int> &)>
+    cell_func = [](const MatrixFree<dim> &                           data,
+                   LinearAlgebra::distributed::Vector<double> &      out,
+                   const LinearAlgebra::distributed::Vector<double> &in,
+                   const std::pair<unsigned int, unsigned int> &range) -> void {
+    FEEvaluation<dim, fe_degree> eval(data);
+
+    for (unsigned int cell = range.first; cell < range.second; ++cell)
+      {
+        eval.reinit(cell);
+        eval.gather_evaluate(in, false, true);
+        for (unsigned int q = 0; q < eval.n_q_points; ++q)
+          {
+            eval.submit_gradient(eval.get_gradient(q), q);
+            eval.submit_value(eval.quadrature_point(q).square(), q);
+          }
+        eval.integrate_scatter(true, true, out);
+      }
+  };
+
+  std::function<void(const MatrixFree<dim> &,
+                     LinearAlgebra::distributed::Vector<double> &,
+                     const LinearAlgebra::distributed::Vector<double> &,
+                     const std::pair<unsigned int, unsigned int> &)>
+    boundary_func =
+      [](const MatrixFree<dim> &                           data,
+         LinearAlgebra::distributed::Vector<double> &      out,
+         const LinearAlgebra::distributed::Vector<double> &in,
+         const std::pair<unsigned int, unsigned int> &     range) -> void {
+    FEFaceEvaluation<dim, fe_degree> eval(data, true);
+
+    for (unsigned int face = range.first; face < range.second; ++face)
+      {
+        eval.reinit(face);
+        eval.gather_evaluate(in, true, false);
+        for (unsigned int q = 0; q < eval.n_q_points; ++q)
+          {
+            eval.submit_value(eval.quadrature_point(q).square() -
+                                6. * eval.get_value(q),
+                              q);
+          }
+        eval.integrate_scatter(true, false, out);
+      }
+  };
+
+  // compute reference result
+  in.update_ghost_values();
+  cell_func(matrix_free,
+            ref,
+            in,
+            std::make_pair(0U, matrix_free.n_cell_batches()));
+  boundary_func(matrix_free,
+                ref,
+                in,
+                std::make_pair(matrix_free.n_inner_face_batches(),
+                               matrix_free.n_inner_face_batches() +
+                                 matrix_free.n_boundary_face_batches()));
+  ref.compress(VectorOperation::add);
+  in.zero_out_ghosts();
+
+  std::function<void(const MatrixFree<dim> &,
+                     LinearAlgebra::distributed::Vector<double> &,
+                     const LinearAlgebra::distributed::Vector<double> &,
+                     const std::pair<unsigned int, unsigned int> &)>
+    inner_face_func;
+
+
+  // compute result through loop with various update settings. Note that we do
+  // no compute on inner faces, so MatrixFree<dim>::DataAccessOnFaces::none
+  // should be enough.
+  matrix_free.loop(cell_func,
+                   inner_face_func,
+                   boundary_func,
+                   test,
+                   in,
+                   true,
+                   MatrixFree<dim>::DataAccessOnFaces::values,
+                   MatrixFree<dim>::DataAccessOnFaces::values);
+
+  deallog << "Number of dofs: " << dof_handler.n_dofs()
+          << (overlap_communication ? " with overlap" : " without overlap")
+          << std::endl;
+
+  test -= ref;
+  deallog << "Error loop 1: " << test.linfty_norm() << std::endl;
+  matrix_free.loop(cell_func,
+                   inner_face_func,
+                   boundary_func,
+                   test,
+                   in,
+                   true,
+                   MatrixFree<dim>::DataAccessOnFaces::values,
+                   MatrixFree<dim>::DataAccessOnFaces::values);
+  test -= ref;
+  deallog << "Error loop 2: " << test.linfty_norm() << std::endl;
+  matrix_free.loop(cell_func,
+                   inner_face_func,
+                   boundary_func,
+                   test,
+                   in,
+                   true,
+                   MatrixFree<dim>::DataAccessOnFaces::none,
+                   MatrixFree<dim>::DataAccessOnFaces::none);
+  test -= ref;
+  deallog << "Error loop 3: " << test.linfty_norm() << std::endl;
+  matrix_free.loop(cell_func,
+                   inner_face_func,
+                   boundary_func,
+                   test,
+                   in,
+                   true,
+                   MatrixFree<dim>::DataAccessOnFaces::values,
+                   MatrixFree<dim>::DataAccessOnFaces::none);
+  test -= ref;
+  deallog << "Error loop 4: " << test.linfty_norm() << std::endl;
+  matrix_free.loop(cell_func,
+                   inner_face_func,
+                   boundary_func,
+                   test,
+                   in,
+                   true,
+                   MatrixFree<dim>::DataAccessOnFaces::unspecified,
+                   MatrixFree<dim>::DataAccessOnFaces::unspecified);
+  test -= ref;
+  deallog << "Error loop 5: " << test.linfty_norm() << std::endl;
+
+  // compute again, now only cell loop
+  ref = 0;
+  in.update_ghost_values();
+  cell_func(matrix_free,
+            ref,
+            in,
+            std::make_pair(0U, matrix_free.n_cell_batches()));
+  ref.compress(VectorOperation::add);
+  in.zero_out_ghosts();
+
+  matrix_free.cell_loop(cell_func, test, in, true);
+  test -= ref;
+  deallog << "Error cell loop: " << test.linfty_norm() << std::endl;
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+  mpi_initlog();
+
+  do_test<2, 1>(2, false);
+  do_test<2, 1>(2, true);
+  do_test<2, 1>(3, false);
+  do_test<2, 1>(3, true);
+  do_test<2, 2>(2, false);
+  do_test<2, 2>(2, true);
+  do_test<2, 2>(3, false);
+  do_test<2, 2>(3, true);
+  do_test<3, 2>(2, false);
+  do_test<3, 2>(2, true);
+}

diff --git a/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=1.output b/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..ea28743
--- /dev/null
+++ b/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=1.output
@@ -0,0 +1,71 @@
+
+DEAL::Number of dofs: 25 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 25 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 81 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 81 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 81 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 81 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 289 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 289 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 729 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 729 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000

diff --git a/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=2.output b/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..ea28743
--- /dev/null
+++ b/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=2.output
@@ -0,0 +1,71 @@
+
+DEAL::Number of dofs: 25 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 25 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 81 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 81 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 81 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 81 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 289 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 289 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 729 without overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000
+DEAL::Number of dofs: 729 with overlap
+DEAL::Error loop 1: 0.00000
+DEAL::Error loop 2: 0.00000
+DEAL::Error loop 3: 0.00000
+DEAL::Error loop 4: 0.00000
+DEAL::Error loop 5: 0.00000
+DEAL::Error cell loop: 0.00000