]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add tests for MatrixFree + SHMEM 13059/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 10 Dec 2021 19:10:59 +0000 (20:10 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 10 Dec 2021 19:50:16 +0000 (20:50 +0100)
include/deal.II/matrix_free/evaluation_kernels.h
tests/matrix_free/ecl.h [new file with mode: 0644]
tests/matrix_free/ecl_01.cc
tests/matrix_free/ecl_01.with_p4est=true.mpirun=1.output [moved from tests/matrix_free/ecl_01.output with 100% similarity]
tests/matrix_free/ecl_03.cc [new file with mode: 0644]
tests/matrix_free/ecl_03.with_p4est=true.mpirun=5.output [new file with mode: 0644]
tests/matrix_free/ecl_04.cc [new file with mode: 0644]
tests/matrix_free/ecl_04.with_p4est=true.mpirun=5.output [new file with mode: 0644]

index dc84931d6343e682a11b4d8cc2d6447d99743677..81ff490f76840f0c9b2128870c1eba25fab69121 100644 (file)
@@ -3996,6 +3996,8 @@ namespace internal
                               proc.value(
                                 temp1[reorientate(v, i)][v],
                                 vector_ptrs[v][index_array_nodal[v][i]]);
+                            else if (integrate == false)
+                              temp1[i][v] = 0.0;
                         }
                   }
               }
diff --git a/tests/matrix_free/ecl.h b/tests/matrix_free/ecl.h
new file mode 100644 (file)
index 0000000..01efcb6
--- /dev/null
@@ -0,0 +1,320 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+template <int dim, typename Number>
+class CosineFunction : public Function<dim, Number>
+{
+public:
+  Number
+  value(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    (void)component;
+    return std::cos(2 * numbers::PI * p[0]);
+  }
+};
+
+template <int dim,
+          int fe_degree,
+          int n_points                 = fe_degree + 1,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test(const unsigned int n_refinements = 1,
+     const bool         print_vector  = true,
+     const MPI_Comm     comm          = MPI_COMM_SELF)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria, 0.0, 1.0, true);
+
+  if (false)
+    {
+      std::vector<dealii::GridTools::PeriodicFacePair<
+        typename dealii::Triangulation<dim>::cell_iterator>>
+        periodic_faces;
+
+      if (dim >= 1)
+        dealii::GridTools::collect_periodic_faces(
+          tria, 0, 1, 0, periodic_faces);
+
+      if (dim >= 2)
+        dealii::GridTools::collect_periodic_faces(
+          tria, 2, 3, 1, periodic_faces);
+
+      if (dim >= 3)
+        dealii::GridTools::collect_periodic_faces(
+          tria, 4, 5, 2, periodic_faces);
+
+      tria.add_periodicity(periodic_faces);
+    }
+
+  tria.refine_global(n_refinements);
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  QGauss<1>     quad(n_points);
+
+  AffineConstraints<Number> constraint;
+
+  using MF = MatrixFree<dim, Number, VectorizedArrayType>;
+
+  typename MF::AdditionalData additional_data;
+  additional_data.mapping_update_flags = update_values | update_gradients;
+  additional_data.mapping_update_flags_inner_faces =
+    update_values | update_gradients;
+  additional_data.mapping_update_flags_boundary_faces =
+    update_values | update_gradients;
+  additional_data.tasks_parallel_scheme =
+    MF::AdditionalData::TasksParallelScheme::none;
+  // additional_data.mapping_update_flags_faces_by_cells = update_values;
+  // additional_data.hold_all_faces_to_owned_cells       = true;
+  additional_data.communicator_sm = comm;
+
+  MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data);
+
+  MF matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
+
+  VectorType src, dst;
+
+  matrix_free.initialize_dof_vector(src);
+  matrix_free.initialize_dof_vector(dst);
+
+  FEEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType> phi(
+    matrix_free);
+  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+    phi_m(matrix_free, true);
+  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+    phi_p(matrix_free, false);
+
+
+  CosineFunction<dim, Number> function;
+  VectorTools::interpolate(dof_handler, function, src);
+
+  dst = 0.0;
+
+  /**
+   * Face-centric loop
+   */
+  matrix_free.template loop<VectorType, VectorType>(
+    [&](const auto &, auto &dst, const auto &src, const auto range) {
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        {
+          phi.reinit(cell);
+          phi.read_dof_values(src);
+          phi.evaluate(false, true, false);
+          for (unsigned int q = 0; q < phi.n_q_points; ++q)
+            phi.submit_gradient(phi.get_gradient(q), q);
+          phi.integrate(false, true);
+          phi.set_dof_values(dst);
+        }
+    },
+    [&](const auto &, auto &dst, const auto &src, const auto range) {
+      for (unsigned int face = range.first; face < range.second; ++face)
+        {
+          phi_m.reinit(face);
+          phi_p.reinit(face);
+
+          phi_m.read_dof_values(src);
+          phi_m.evaluate(true, true);
+          phi_p.read_dof_values(src);
+          phi_p.evaluate(true, true);
+          VectorizedArrayType sigmaF =
+            (std::abs((phi_m.get_normal_vector(0) *
+                       phi_m.inverse_jacobian(0))[dim - 1]) +
+             std::abs((phi_m.get_normal_vector(0) *
+                       phi_p.inverse_jacobian(0))[dim - 1])) *
+            (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
+
+          for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+            {
+              VectorizedArrayType average_value =
+                (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
+              VectorizedArrayType average_valgrad =
+                phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q);
+              average_valgrad =
+                average_value * 2. * sigmaF - average_valgrad * 0.5;
+              phi_m.submit_normal_derivative(-average_value, q);
+              phi_p.submit_normal_derivative(-average_value, q);
+              phi_m.submit_value(average_valgrad, q);
+              phi_p.submit_value(-average_valgrad, q);
+            }
+          phi_m.integrate(true, true);
+          phi_m.distribute_local_to_global(dst);
+          phi_p.integrate(true, true);
+          phi_p.distribute_local_to_global(dst);
+        }
+    },
+    [&](const auto &, auto &dst, const auto &src, const auto face_range) {
+      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(EvaluationFlags::values | EvaluationFlags::gradients);
+          VectorizedArrayType sigmaF =
+            std::abs((phi_m.get_normal_vector(0) *
+                      phi_m.inverse_jacobian(0))[dim - 1]) *
+            Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
+
+          for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+            {
+              VectorizedArrayType average_value = phi_m.get_value(q);
+              VectorizedArrayType average_valgrad =
+                -phi_m.get_normal_derivative(q);
+              average_valgrad += average_value * sigmaF * 2.0;
+              phi_m.submit_normal_derivative(-average_value, q);
+              phi_m.submit_value(average_valgrad, q);
+            }
+
+          phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+          phi_m.distribute_local_to_global(dst);
+        }
+    },
+    dst,
+    src,
+    false,
+    MF::DataAccessOnFaces::gradients,
+    MF::DataAccessOnFaces::gradients);
+
+  deallog << dst.l2_norm() << std::endl;
+
+  if (print_vector)
+    dst.print(deallog.get_file_stream());
+
+  dst = 0.0;
+
+  /**
+   * Element-centric loop
+   */
+  matrix_free.template loop_cell_centric<VectorType, VectorType>(
+    [&](const auto &, auto &dst, const auto &src, const auto range) {
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        {
+          phi.reinit(cell);
+          phi.read_dof_values(src);
+          phi.evaluate(false, true, false);
+          for (unsigned int q = 0; q < phi.n_q_points; ++q)
+            phi.submit_gradient(phi.get_gradient(q), q);
+          phi.integrate(false, true);
+
+          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+               ++face)
+            {
+              auto bids =
+                matrix_free.get_faces_by_cells_boundary_id(cell, face);
+
+              if (bids[0] != numbers::internal_face_boundary_id)
+                {
+                  phi_m.reinit(cell, face);
+                  phi_m.read_dof_values(src);
+                  phi_m.evaluate(EvaluationFlags::values |
+                                 EvaluationFlags::gradients);
+                  VectorizedArrayType sigmaF =
+                    std::abs((phi_m.get_normal_vector(0) *
+                              phi_m.inverse_jacobian(0))[dim - 1]) *
+                    Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
+
+                  for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+                    {
+                      VectorizedArrayType average_value = phi_m.get_value(q);
+                      VectorizedArrayType average_valgrad =
+                        -phi_m.get_normal_derivative(q);
+                      average_valgrad += average_value * sigmaF * 2.0;
+                      phi_m.submit_normal_derivative(-average_value, q);
+                      phi_m.submit_value(average_valgrad, q);
+                    }
+
+                  phi_m.integrate(EvaluationFlags::values |
+                                  EvaluationFlags::gradients);
+
+                  for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
+                    phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q];
+                }
+              else
+                {
+                  phi_m.reinit(cell, face);
+                  phi_p.reinit(cell, face);
+
+                  phi_m.read_dof_values(src);
+                  phi_m.evaluate(EvaluationFlags::values |
+                                 EvaluationFlags::gradients);
+                  phi_p.read_dof_values(src);
+                  phi_p.evaluate(EvaluationFlags::values |
+                                 EvaluationFlags::gradients);
+
+                  VectorizedArrayType sigmaF =
+                    (std::abs((phi_m.get_normal_vector(0) *
+                               phi_m.inverse_jacobian(0))[dim - 1]) +
+                     std::abs((phi_m.get_normal_vector(0) *
+                               phi_p.inverse_jacobian(0))[dim - 1])) *
+                    (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
+
+                  for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+                    {
+                      VectorizedArrayType average_value =
+                        (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
+                      VectorizedArrayType average_valgrad =
+                        phi_m.get_normal_derivative(q) +
+                        phi_p.get_normal_derivative(q);
+                      average_valgrad =
+                        average_value * 2. * sigmaF - average_valgrad * 0.5;
+                      phi_m.submit_normal_derivative(-average_value, q);
+                      phi_m.submit_value(average_valgrad, q);
+                    }
+                  phi_m.integrate(EvaluationFlags::values |
+                                  EvaluationFlags::gradients);
+
+                  for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
+                    phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q];
+                }
+            }
+
+          phi.set_dof_values(dst);
+        }
+    },
+    dst,
+    src,
+    false,
+    MF::DataAccessOnFaces::gradients);
+
+  deallog << dst.l2_norm() << std::endl;
+
+  if (print_vector)
+    dst.print(deallog.get_file_stream());
+}
index ad85f1d5a20cf0df9b63e62fd5e67c4b8f6482da..a335c637d5654c3e91949affeb967497206b4b66 100644 (file)
 
 #include <deal.II/numerics/vector_tools.h>
 
-#include "../tests.h"
+#include "ecl.h"
 
 
 // Compare the evaluation of a SIP Laplace operator with face- and element-
 // centric loops on a structured grid.
 
-template <int dim, typename Number>
-class CosineFunction : public Function<dim, Number>
-{
-public:
-  Number
-  value(const Point<dim> &p, const unsigned int component = 0) const
-  {
-    (void)component;
-    return std::cos(2 * numbers::PI * p[0]);
-  }
-};
-
-template <int dim,
-          int fe_degree,
-          int n_points                 = fe_degree + 1,
-          typename Number              = double,
-          typename VectorizedArrayType = VectorizedArray<Number>>
-void
-test(const unsigned int n_refinements = 1)
-{
-  using VectorType = LinearAlgebra::distributed::Vector<Number>;
-
-  Triangulation<dim> tria;
-  GridGenerator::hyper_cube(tria, 0.0, 1.0, true);
-
-  if (false)
-    {
-      std::vector<dealii::GridTools::PeriodicFacePair<
-        typename dealii::Triangulation<dim>::cell_iterator>>
-        periodic_faces;
-
-      if (dim >= 1)
-        dealii::GridTools::collect_periodic_faces(
-          tria, 0, 1, 0, periodic_faces);
-
-      if (dim >= 2)
-        dealii::GridTools::collect_periodic_faces(
-          tria, 2, 3, 1, periodic_faces);
-
-      if (dim >= 3)
-        dealii::GridTools::collect_periodic_faces(
-          tria, 4, 5, 2, periodic_faces);
-
-      tria.add_periodicity(periodic_faces);
-    }
-
-  tria.refine_global(n_refinements);
-
-  FE_DGQ<dim>     fe(fe_degree);
-  DoFHandler<dim> dof_handler(tria);
-  dof_handler.distribute_dofs(fe);
-
-  MappingQ<dim> mapping(1);
-  QGauss<1>     quad(n_points);
-
-  AffineConstraints<Number> constraint;
-
-  using MF = MatrixFree<dim, Number, VectorizedArrayType>;
-
-  typename MF::AdditionalData additional_data;
-  additional_data.mapping_update_flags = update_values | update_gradients;
-  additional_data.mapping_update_flags_inner_faces =
-    update_values | update_gradients;
-  additional_data.mapping_update_flags_boundary_faces =
-    update_values | update_gradients;
-  additional_data.tasks_parallel_scheme =
-    MF::AdditionalData::TasksParallelScheme::none;
-  // additional_data.mapping_update_flags_faces_by_cells = update_values;
-  // additional_data.hold_all_faces_to_owned_cells       = true;
-
-  MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data);
-
-  MF matrix_free;
-  matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
-
-  VectorType src, dst;
-
-  matrix_free.initialize_dof_vector(src);
-  matrix_free.initialize_dof_vector(dst);
-
-  FEEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType> phi(
-    matrix_free);
-  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
-    phi_m(matrix_free, true);
-  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
-    phi_p(matrix_free, false);
-
-
-  CosineFunction<dim, Number> function;
-  VectorTools::interpolate(dof_handler, function, src);
-
-  dst = 0.0;
-
-  /**
-   * Face-centric loop
-   */
-  matrix_free.template loop<VectorType, VectorType>(
-    [&](const auto &, auto &dst, const auto &src, const auto range) {
-      for (unsigned int cell = range.first; cell < range.second; ++cell)
-        {
-          phi.reinit(cell);
-          phi.read_dof_values(src);
-          phi.evaluate(false, true, false);
-          for (unsigned int q = 0; q < phi.n_q_points; ++q)
-            phi.submit_gradient(phi.get_gradient(q), q);
-          phi.integrate(false, true);
-          phi.set_dof_values(dst);
-        }
-    },
-    [&](const auto &, auto &dst, const auto &src, const auto range) {
-      for (unsigned int face = range.first; face < range.second; ++face)
-        {
-          phi_m.reinit(face);
-          phi_p.reinit(face);
-
-          phi_m.read_dof_values(src);
-          phi_m.evaluate(true, true);
-          phi_p.read_dof_values(src);
-          phi_p.evaluate(true, true);
-          VectorizedArrayType sigmaF =
-            (std::abs((phi_m.get_normal_vector(0) *
-                       phi_m.inverse_jacobian(0))[dim - 1]) +
-             std::abs((phi_m.get_normal_vector(0) *
-                       phi_p.inverse_jacobian(0))[dim - 1])) *
-            (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
-
-          for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
-            {
-              VectorizedArrayType average_value =
-                (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
-              VectorizedArrayType average_valgrad =
-                phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q);
-              average_valgrad =
-                average_value * 2. * sigmaF - average_valgrad * 0.5;
-              phi_m.submit_normal_derivative(-average_value, q);
-              phi_p.submit_normal_derivative(-average_value, q);
-              phi_m.submit_value(average_valgrad, q);
-              phi_p.submit_value(-average_valgrad, q);
-            }
-          phi_m.integrate(true, true);
-          phi_m.distribute_local_to_global(dst);
-          phi_p.integrate(true, true);
-          phi_p.distribute_local_to_global(dst);
-        }
-    },
-    [&](const auto &, auto &dst, const auto &src, const auto face_range) {
-      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(EvaluationFlags::values | EvaluationFlags::gradients);
-          VectorizedArrayType sigmaF =
-            std::abs((phi_m.get_normal_vector(0) *
-                      phi_m.inverse_jacobian(0))[dim - 1]) *
-            Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
-
-          for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
-            {
-              VectorizedArrayType average_value = phi_m.get_value(q);
-              VectorizedArrayType average_valgrad =
-                -phi_m.get_normal_derivative(q);
-              average_valgrad += average_value * sigmaF * 2.0;
-              phi_m.submit_normal_derivative(-average_value, q);
-              phi_m.submit_value(average_valgrad, q);
-            }
-
-          phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
-          phi_m.distribute_local_to_global(dst);
-        }
-    },
-    dst,
-    src,
-    false,
-    MF::DataAccessOnFaces::gradients,
-    MF::DataAccessOnFaces::gradients);
-
-  deallog << dst.l2_norm() << std::endl;
-  dst.print(deallog.get_file_stream());
-
-  dst = 0.0;
-
-  /**
-   * Element-centric loop
-   */
-  matrix_free.template loop_cell_centric<VectorType, VectorType>(
-    [&](const auto &, auto &dst, const auto &src, const auto range) {
-      for (unsigned int cell = range.first; cell < range.second; ++cell)
-        {
-          phi.reinit(cell);
-          phi.read_dof_values(src);
-          phi.evaluate(false, true, false);
-          for (unsigned int q = 0; q < phi.n_q_points; ++q)
-            phi.submit_gradient(phi.get_gradient(q), q);
-          phi.integrate(false, true);
-
-          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
-               ++face)
-            {
-              auto bids =
-                matrix_free.get_faces_by_cells_boundary_id(cell, face);
-
-              if (bids[0] != numbers::internal_face_boundary_id)
-                {
-                  phi_m.reinit(cell, face);
-                  phi_m.read_dof_values(src);
-                  phi_m.evaluate(EvaluationFlags::values |
-                                 EvaluationFlags::gradients);
-                  VectorizedArrayType sigmaF =
-                    std::abs((phi_m.get_normal_vector(0) *
-                              phi_m.inverse_jacobian(0))[dim - 1]) *
-                    Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
-
-                  for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
-                    {
-                      VectorizedArrayType average_value = phi_m.get_value(q);
-                      VectorizedArrayType average_valgrad =
-                        -phi_m.get_normal_derivative(q);
-                      average_valgrad += average_value * sigmaF * 2.0;
-                      phi_m.submit_normal_derivative(-average_value, q);
-                      phi_m.submit_value(average_valgrad, q);
-                    }
-
-                  phi_m.integrate(EvaluationFlags::values |
-                                  EvaluationFlags::gradients);
-
-                  for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
-                    phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q];
-                }
-              else
-                {
-                  phi_m.reinit(cell, face);
-                  phi_p.reinit(cell, face);
-
-                  phi_m.read_dof_values(src);
-                  phi_m.evaluate(EvaluationFlags::values |
-                                 EvaluationFlags::gradients);
-                  phi_p.read_dof_values(src);
-                  phi_p.evaluate(EvaluationFlags::values |
-                                 EvaluationFlags::gradients);
-
-                  VectorizedArrayType sigmaF =
-                    (std::abs((phi_m.get_normal_vector(0) *
-                               phi_m.inverse_jacobian(0))[dim - 1]) +
-                     std::abs((phi_m.get_normal_vector(0) *
-                               phi_p.inverse_jacobian(0))[dim - 1])) *
-                    (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
-
-                  for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
-                    {
-                      VectorizedArrayType average_value =
-                        (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
-                      VectorizedArrayType average_valgrad =
-                        phi_m.get_normal_derivative(q) +
-                        phi_p.get_normal_derivative(q);
-                      average_valgrad =
-                        average_value * 2. * sigmaF - average_valgrad * 0.5;
-                      phi_m.submit_normal_derivative(-average_value, q);
-                      phi_m.submit_value(average_valgrad, q);
-                    }
-                  phi_m.integrate(EvaluationFlags::values |
-                                  EvaluationFlags::gradients);
-
-                  for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
-                    phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q];
-                }
-            }
-
-          phi.set_dof_values(dst);
-        }
-    },
-    dst,
-    src,
-    false,
-    MF::DataAccessOnFaces::gradients);
-
-  deallog << dst.l2_norm() << std::endl;
-  dst.print(deallog.get_file_stream());
-}
-
 int
-main()
+main(int argc, char **argv)
 {
-  initlog();
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
   test<2, 2, 3, double, VectorizedArray<double>>();
 }
diff --git a/tests/matrix_free/ecl_03.cc b/tests/matrix_free/ecl_03.cc
new file mode 100644 (file)
index 0000000..62ce5fb
--- /dev/null
@@ -0,0 +1,65 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "ecl.h"
+
+
+// Like ecl_01.cc but running with shared-memory MPI enabled.
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+
+  // test without shared memory
+  {
+    test<2, 2, 3, double, VectorizedArray<double>>(6, false, MPI_COMM_SELF);
+  }
+
+  // test with shared memory
+  {
+    const auto rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+    MPI_Comm subcommunicator;
+    MPI_Comm_split_type(MPI_COMM_WORLD,
+                        MPI_COMM_TYPE_SHARED,
+                        rank,
+                        MPI_INFO_NULL,
+                        &subcommunicator);
+
+    test<2, 2, 3, double, VectorizedArray<double>>(6, false, subcommunicator);
+
+    MPI_Comm_free(&subcommunicator);
+  }
+}
diff --git a/tests/matrix_free/ecl_03.with_p4est=true.mpirun=5.output b/tests/matrix_free/ecl_03.with_p4est=true.mpirun=5.output
new file mode 100644 (file)
index 0000000..62b728c
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::210.386
+DEAL::210.386
+DEAL::210.386
+DEAL::210.386
diff --git a/tests/matrix_free/ecl_04.cc b/tests/matrix_free/ecl_04.cc
new file mode 100644 (file)
index 0000000..e746eca
--- /dev/null
@@ -0,0 +1,205 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+// Test FEFaceEvaluation::gather_evaluate() for values, ECL, shared-memory MPI,
+// and different geometries (i.e., different orientations).
+
+template <int dim>
+class ExactSolution : public Function<dim>
+{
+public:
+  ExactSolution()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int /*component*/ = 0) const
+  {
+    double result = p[0];
+
+    for (unsigned i = 1; i < dim; ++i)
+      result *= p[i];
+
+    return result;
+  }
+};
+
+
+template <int dim,
+          int fe_degree,
+          int n_points                 = fe_degree + 1,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test(const unsigned int geometry, const MPI_Comm comm = MPI_COMM_SELF)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+  if (geometry == 0)
+    {
+      GridGenerator::hyper_cube(tria);
+
+      if (dim == 2)
+        tria.refine_global(4);
+      else
+        tria.refine_global(3);
+    }
+  else if (geometry == 1)
+    {
+      GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1.0);
+
+      if (dim == 2)
+        tria.refine_global(4);
+      else
+        tria.refine_global(3);
+    }
+  else if (geometry == 2)
+    {
+      GridGenerator::hyper_ball(tria, Point<dim>(), 1.0);
+
+      if (dim == 2)
+        tria.refine_global(4);
+      else
+        tria.refine_global(3);
+    }
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  QGauss<1>     quad(n_points);
+
+  AffineConstraints<Number> constraint;
+
+  using MF = MatrixFree<dim, Number, VectorizedArrayType>;
+
+  typename MF::AdditionalData additional_data;
+  additional_data.mapping_update_flags                = update_values;
+  additional_data.mapping_update_flags_inner_faces    = update_values;
+  additional_data.mapping_update_flags_boundary_faces = update_values;
+  additional_data.mapping_update_flags_faces_by_cells = update_values;
+  additional_data.hold_all_faces_to_owned_cells       = true;
+  additional_data.communicator_sm                     = comm;
+
+  MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data);
+
+  MF matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
+
+  VectorType src, dst;
+
+  matrix_free.initialize_dof_vector(src);
+  matrix_free.initialize_dof_vector(dst);
+
+  VectorTools::interpolate(dof_handler, ExactSolution<dim>(), src);
+
+  dst = 0.0;
+
+  /**
+   * Element-centric loop
+   */
+  matrix_free.template loop_cell_centric<VectorType, VectorType>(
+    [&](const auto &, auto &dst, const auto &src, const auto range) {
+      FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+        phi_m(matrix_free, true);
+      FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+        phi_p(matrix_free, false);
+
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        {
+          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+               face++)
+            {
+              auto bids =
+                matrix_free.get_faces_by_cells_boundary_id(cell, face);
+
+              if (bids[0] != numbers::internal_face_boundary_id)
+                continue;
+
+              phi_m.reinit(cell, face);
+              phi_p.reinit(cell, face);
+
+              phi_m.gather_evaluate(src, EvaluationFlags::values);
+              phi_p.gather_evaluate(src, EvaluationFlags::values);
+
+              for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+                {
+                  const auto u_minus = phi_m.get_value(q);
+                  const auto u_plus  = phi_p.get_value(q);
+
+                  for (unsigned int v = 0; v < VectorizedArray<double>::size();
+                       ++v)
+                    {
+                      Assert(std::abs(u_minus[v] - u_plus[v]) < 1e-10,
+                             ExcMessage("Entries do not match!"));
+                    }
+                }
+            }
+        }
+    },
+    dst,
+    src,
+    false,
+    MF::DataAccessOnFaces::values);
+
+  deallog << "OK" << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+
+  const auto rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  MPI_Comm subcommunicator;
+  MPI_Comm_split_type(MPI_COMM_WORLD,
+                      MPI_COMM_TYPE_SHARED,
+                      rank,
+                      MPI_INFO_NULL,
+                      &subcommunicator);
+
+  for (unsigned int i = 0; i < 3; ++i)
+    test<2, 3, 4, double>(i, subcommunicator);
+
+  for (unsigned int i = 0; i < 3; ++i)
+    test<3, 3, 4, double>(i, subcommunicator);
+
+  MPI_Comm_free(&subcommunicator);
+}
diff --git a/tests/matrix_free/ecl_04.with_p4est=true.mpirun=5.output b/tests/matrix_free/ecl_04.with_p4est=true.mpirun=5.output
new file mode 100644 (file)
index 0000000..0aa61ff
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK

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.