From 902ca5abdf0ad45635daf7e8747c436139c754d9 Mon Sep 17 00:00:00 2001
From: Bruno Turcksin <bruno.turcksin@gmail.com>
Date: Thu, 8 Aug 2019 21:28:22 +0000
Subject: [PATCH] Switch the tests to atomic operations and add a test for
 graph coloring

---
 tests/cuda/matrix_free_matrix_vector_03b.cu   | 77 +++++++++++++++++++
 .../cuda/matrix_free_matrix_vector_03b.output | 16 ++++
 tests/cuda/matrix_vector_common.h             |  4 +-
 3 files changed, 96 insertions(+), 1 deletion(-)
 create mode 100644 tests/cuda/matrix_free_matrix_vector_03b.cu
 create mode 100644 tests/cuda/matrix_free_matrix_vector_03b.output

diff --git a/tests/cuda/matrix_free_matrix_vector_03b.cu b/tests/cuda/matrix_free_matrix_vector_03b.cu
new file mode 100644
index 0000000000..772bcba081
--- /dev/null
+++ b/tests/cuda/matrix_free_matrix_vector_03b.cu
@@ -0,0 +1,77 @@
+// ---------------------------------------------------------------------
+//
+// 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 implementation of matrix free
+// matrix-vector products by comparing with the result of deal.II sparse
+// matrix. The mesh uses a hypercube mesh with hanging nodes (created by
+// randomly refining some cells) and zero Dirichlet conditions. Same as test
+// matrix_free_matrix_vector_03b but use coloring instead of atomics
+
+#include <deal.II/base/function.h>
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_common.h"
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
+                                                    endc = tria.end();
+  if (dim < 3 || fe_degree < 2)
+    tria.refine_global(2);
+  else
+    tria.refine_global(1);
+  tria.begin(tria.n_levels() - 1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  cell = tria.begin_active();
+  for (unsigned int i = 0; i < 10 - 3 * dim; ++i)
+    {
+      cell                 = tria.begin_active();
+      unsigned int counter = 0;
+      for (; cell != endc; ++cell, ++counter)
+        if (counter % (7 - i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+
+  VectorTools::interpolate_boundary_values(dof,
+                                           0,
+                                           Functions::ZeroFunction<dim>(),
+                                           constraints);
+  constraints.close();
+
+  // Skip 2D tests with even fe_degree
+  if ((dim == 3) || ((fe_degree % 2) == 1))
+    do_test<dim,
+            fe_degree,
+            double,
+            LinearAlgebra::CUDAWrappers::Vector<double>,
+            fe_degree + 1>(dof, constraints, tria.n_active_cells(), true, true);
+}
diff --git a/tests/cuda/matrix_free_matrix_vector_03b.output b/tests/cuda/matrix_free_matrix_vector_03b.output
new file mode 100644
index 0000000000..ab11bfed67
--- /dev/null
+++ b/tests/cuda/matrix_free_matrix_vector_03b.output
@@ -0,0 +1,16 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference: 0.
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Norm of difference: 0.
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Norm of difference: 0.
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Norm of difference: 0.
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(3)
+DEAL:3d::Norm of difference: 0.
+DEAL:3d::
diff --git a/tests/cuda/matrix_vector_common.h b/tests/cuda/matrix_vector_common.h
index 80dbb8aec5..2c3586b2d9 100644
--- a/tests/cuda/matrix_vector_common.h
+++ b/tests/cuda/matrix_vector_common.h
@@ -78,7 +78,8 @@ void
 do_test(const DoFHandler<dim> &          dof,
         const AffineConstraints<double> &constraints,
         const unsigned int               n_locally_owned_cells,
-        const bool                       constant_coefficient = true)
+        const bool                       constant_coefficient = true,
+        const bool                       coloring             = false)
 {
   deallog << "Testing " << dof.get_fe().get_name() << std::endl;
 
@@ -89,6 +90,7 @@ do_test(const DoFHandler<dim> &          dof,
   additional_data.mapping_update_flags = update_values | update_gradients |
                                          update_JxW_values |
                                          update_quadrature_points;
+  additional_data.use_coloring = coloring;
   const QGauss<1> quad(n_q_points_1d);
   mf_data.reinit(mapping, dof, constraints, quad, additional_data);
 
-- 
2.39.5