Adapt current tests to be able to switch between constant and varying coefficients

diff --git a/tests/cuda/matrix_free_matrix_vector_01.cu b/tests/cuda/matrix_free_matrix_vector_01.cu
index efd928d1193b724a0087eac38c2eedbf024d0fee..ce08bb12f1965ee53fb78b3951a7eca62bb88aa2 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_01.cu
+++ b/tests/cuda/matrix_free_matrix_vector_01.cu
@@ -45,5 +45,5 @@ test()
           fe_degree,
           double,
           LinearAlgebra::CUDAWrappers::Vector<double>,
-          fe_degree + 1>(dof, constraints);
+          fe_degree + 1>(dof, constraints, tria.n_active_cells());
 }

diff --git a/tests/cuda/matrix_free_matrix_vector_02.cu b/tests/cuda/matrix_free_matrix_vector_02.cu
index c6388d802418b5040b7c4f6c79598e5ca1c21182..166b73cf87044a4f5e2bd48c97c08ccebdad51d9 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_02.cu
+++ b/tests/cuda/matrix_free_matrix_vector_02.cu
@@ -49,5 +49,5 @@ test()
           fe_degree,
           double,
           LinearAlgebra::CUDAWrappers::Vector<double>,
-          fe_degree + 1>(dof, constraints);
+          fe_degree + 1>(dof, constraints, tria.n_active_cells());
 }

diff --git a/tests/cuda/matrix_free_matrix_vector_03.cu b/tests/cuda/matrix_free_matrix_vector_03.cu
index cd556e16ee71bc5d3c579078d07ad8962d6000e6..147febd3f506f138b8588d1223d170d10f7212ad 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_03.cu
+++ b/tests/cuda/matrix_free_matrix_vector_03.cu
@@ -72,5 +72,5 @@ test()
             fe_degree,
             double,
             LinearAlgebra::CUDAWrappers::Vector<double>,
-            fe_degree + 1>(dof, constraints);
+            fe_degree + 1>(dof, constraints, tria.n_active_cells());
 }

diff --git a/tests/cuda/matrix_free_matrix_vector_04.cu b/tests/cuda/matrix_free_matrix_vector_04.cu
index 61dd8b491d7d87306feb43fb850c8e9710fdc27c..c4c1085670d3b8270bbe880966b2a6358bec67f4 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_04.cu
+++ b/tests/cuda/matrix_free_matrix_vector_04.cu
@@ -48,5 +48,5 @@ test()
           fe_degree,
           double,
           LinearAlgebra::CUDAWrappers::Vector<double>,
-          fe_degree + 1>(dof, constraints);
+          fe_degree + 1>(dof, constraints, tria.n_active_cells());
 }

diff --git a/tests/cuda/matrix_free_matrix_vector_05.cu b/tests/cuda/matrix_free_matrix_vector_05.cu
index 54ce17e03231aa209b6538764c5826865d3486c4..48efc829d1cd2566de413bb995a8b3872495852d 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_05.cu
+++ b/tests/cuda/matrix_free_matrix_vector_05.cu
@@ -82,5 +82,5 @@ test()
             fe_degree,
             double,
             LinearAlgebra::CUDAWrappers::Vector<double>,
-            fe_degree + 1>(dof, constraints);
+            fe_degree + 1>(dof, constraints, tria.n_active_cells());
 }

diff --git a/tests/cuda/matrix_free_matrix_vector_06.cu b/tests/cuda/matrix_free_matrix_vector_06.cu
index 55a4540e3dc2f10b93384c5a06991bbc7c184758..704255a86ead47b1fb5a9ebd16198fbc3d3fab1e 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_06.cu
+++ b/tests/cuda/matrix_free_matrix_vector_06.cu
@@ -78,5 +78,5 @@ test()
           fe_degree,
           double,
           LinearAlgebra::CUDAWrappers::Vector<double>,
-          fe_degree + 1>(dof, constraints);
+          fe_degree + 1>(dof, constraints, tria.n_active_cells());
 }

diff --git a/tests/cuda/matrix_free_matrix_vector_06a.cu b/tests/cuda/matrix_free_matrix_vector_06a.cu
index c9c5d1cb42071b1dacfb6884be83f4d3bb36aedb..7aa0b084fda8972879db11ae9c946432f3536154 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_06a.cu
+++ b/tests/cuda/matrix_free_matrix_vector_06a.cu
@@ -74,5 +74,5 @@ test()
           fe_degree,
           double,
           LinearAlgebra::distributed::Vector<double, MemorySpace::CUDA>,
-          fe_degree + 1>(dof, constraints);
+          fe_degree + 1>(dof, constraints, tria.n_active_cells());
 }

diff --git a/tests/cuda/matrix_free_matrix_vector_09.cu b/tests/cuda/matrix_free_matrix_vector_09.cu
index 1224ed4b6854889ced1fa1a3e9e570a16b47e3da..24f36ed4f3ccb5057fd5faed89d4109dd3f9f5ed 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_09.cu
+++ b/tests/cuda/matrix_free_matrix_vector_09.cu
@@ -76,5 +76,5 @@ test()
           fe_degree,
           double,
           LinearAlgebra::CUDAWrappers::Vector<double>,
-          fe_degree + 1>(dof, constraints);
+          fe_degree + 1>(dof, constraints, tria.n_active_cells());
 }

diff --git a/tests/cuda/matrix_free_matrix_vector_10.cu b/tests/cuda/matrix_free_matrix_vector_10.cu
index fd26148e4e34f010ac2e72c1dc5700e66ce892ac..63a887d5c0e313ab5c61b2dc654f9592d6919021 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_10.cu
+++ b/tests/cuda/matrix_free_matrix_vector_10.cu
@@ -115,11 +115,13 @@ test()
   mf_data.reinit(
     mapping, dof, constraints, quad, MPI_COMM_WORLD, additional_data);
 
+  const unsigned int coef_size =
+    tria.n_locally_owned_active_cells() * std::pow(fe_degree + 1, dim);
   MatrixFreeTest<dim,
                  fe_degree,
                  Number,
                  LinearAlgebra::distributed::Vector<Number, MemorySpace::CUDA>>
-                                                                mf(mf_data);
+                                                                mf(mf_data, coef_size);
   LinearAlgebra::distributed::Vector<Number, MemorySpace::CUDA> in_dev(
     owned_set, MPI_COMM_WORLD);
   LinearAlgebra::distributed::Vector<Number, MemorySpace::CUDA> out_dev(

diff --git a/tests/cuda/matrix_free_matrix_vector_19.cu b/tests/cuda/matrix_free_matrix_vector_19.cu
index 18c2f1b86be0f0b94d3077f4ad9b2570fefe4d6d..c60eb4058890ec191bf8fa7a06856728dbf8f82c 100644 (file)
--- a/tests/cuda/matrix_free_matrix_vector_19.cu
+++ b/tests/cuda/matrix_free_matrix_vector_19.cu
@@ -111,11 +111,13 @@ test()
   mf_data.reinit(
     mapping, dof, constraints, quad, MPI_COMM_WORLD, additional_data);
 
+  const unsigned int coef_size =
+    tria.n_locally_owned_active_cells() * std::pow(fe_degree + 1, dim);
   MatrixFreeTest<dim,
                  fe_degree,
                  Number,
                  LinearAlgebra::distributed::Vector<Number, MemorySpace::CUDA>>
-                                                                mf(mf_data);
+                                                                mf(mf_data, coef_size);
   LinearAlgebra::distributed::Vector<Number, MemorySpace::CUDA> in_dev(
     owned_set, MPI_COMM_WORLD);
   LinearAlgebra::distributed::Vector<Number, MemorySpace::CUDA> out_dev(

diff --git a/tests/cuda/matrix_vector_common.h b/tests/cuda/matrix_vector_common.h
index 662086f5950b995ec437027d33b7097a53c57505..80dbb8aec586cb0c73f8d72e277e80b3806e9da1 100644 (file)
--- a/tests/cuda/matrix_vector_common.h
+++ b/tests/cuda/matrix_vector_common.h
@@ -19,6 +19,7 @@
 // general, with and without hanging nodes). It also tests the multithreading
 
 #include <deal.II/base/logstream.h>
+#include <deal.II/base/point.h>
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/dofs/dof_handler.h>
@@ -52,6 +53,20 @@ template <int dim, int fe_degree>
 void
 test();
 
+template <int dim>
+double
+ConstantCoefficient(const Point<dim> & /*p*/)
+{
+  return 10.;
+}
+
+
+template <int dim>
+double
+VaryingCoefficient(const Point<dim> &p)
+{
+  return 10. / (0.05 + 2. * p.square());
+}
 
 
 template <int dim,
@@ -61,7 +76,9 @@ template <int dim,
           int n_q_points_1d>
 void
 do_test(const DoFHandler<dim> &          dof,
-        const AffineConstraints<double> &constraints)
+        const AffineConstraints<double> &constraints,
+        const unsigned int               n_locally_owned_cells,
+        const bool                       constant_coefficient = true)
 {
   deallog << "Testing " << dof.get_fe().get_name() << std::endl;
 
@@ -76,7 +93,10 @@ do_test(const DoFHandler<dim> &          dof,
   mf_data.reinit(mapping, dof, constraints, quad, additional_data);
 
   const unsigned int n_dofs = dof.n_dofs();
-  MatrixFreeTest<dim, fe_degree, Number, VectorType, n_q_points_1d> mf(mf_data);
+  MatrixFreeTest<dim, fe_degree, Number, VectorType, n_q_points_1d> mf(
+    mf_data,
+    n_locally_owned_cells * std::pow(n_q_points_1d, dim),
+    constant_coefficient);
   Vector<Number>                         in_host(n_dofs), out_host(n_dofs);
   LinearAlgebra::ReadWriteVector<Number> in(n_dofs), out(n_dofs);
   VectorType                             in_device(n_dofs);
@@ -110,7 +130,7 @@ do_test(const DoFHandler<dim> &          dof,
                             dof.get_fe(),
                             quadrature_formula,
                             update_values | update_gradients |
-                              update_JxW_values);
+                              update_quadrature_points | update_JxW_values);
 
     const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
     const unsigned int n_q_points    = quadrature_formula.size();
@@ -126,15 +146,22 @@ do_test(const DoFHandler<dim> &          dof,
         fe_values.reinit(cell);
 
         for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
-          for (unsigned int i = 0; i < dofs_per_cell; ++i)
-            {
-              for (unsigned int j = 0; j < dofs_per_cell; ++j)
-                cell_matrix(i, j) += ((fe_values.shape_grad(i, q_point) *
-                                         fe_values.shape_grad(j, q_point) +
-                                       10. * fe_values.shape_value(i, q_point) *
-                                         fe_values.shape_value(j, q_point)) *
-                                      fe_values.JxW(q_point));
-            }
+          {
+            const auto coef =
+              constant_coefficient ?
+                ConstantCoefficient<dim>(fe_values.quadrature_point(q_point)) :
+                VaryingCoefficient<dim>(fe_values.quadrature_point(q_point));
+            for (unsigned int i = 0; i < dofs_per_cell; ++i)
+              {
+                for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                  cell_matrix(i, j) +=
+                    ((fe_values.shape_grad(i, q_point) *
+                        fe_values.shape_grad(j, q_point) +
+                      coef * fe_values.shape_value(i, q_point) *
+                        fe_values.shape_value(j, q_point)) *
+                     fe_values.JxW(q_point));
+              }
+          }
 
         cell->get_dof_indices(local_dof_indices);
         constraints.distribute_local_to_global(cell_matrix,
@@ -169,7 +196,7 @@ main()
   {
     deallog.push("2d");
     test<2, 1>();
-    //  test<2, 2>();
+    // test<2, 2>();
     test<2, 3>();
     deallog.pop();
     deallog.push("3d");

diff --git a/tests/cuda/matrix_vector_mf.h b/tests/cuda/matrix_vector_mf.h
index 4deba79eaee869533199ca85d2c9ae583d7cf4bd..4d62d78cec3b63944fa89001fa0846203e7b7a84 100644 (file)
--- a/tests/cuda/matrix_vector_mf.h
+++ b/tests/cuda/matrix_vector_mf.h
@@ -26,14 +26,37 @@
 
 #include "../tests.h"
 
+__device__ unsigned int
+compute_thread_id()
+{
+  const unsigned int thread_num_in_block =
+    threadIdx.x + blockDim.x * threadIdx.y +
+    blockDim.x * blockDim.y * threadIdx.z;
+  const unsigned int n_threads_per_block = blockDim.x * blockDim.y * blockDim.z;
+  const unsigned int blocks_num_in_grid =
+    blockIdx.x + gridDim.x * blockIdx.y + gridDim.x * gridDim.y * blockIdx.z;
+
+  return thread_num_in_block + blocks_num_in_grid * n_threads_per_block;
+}
+
+
+
 template <int dim, int fe_degree, typename Number, int n_q_points_1d>
 class HelmholtzOperatorQuad
 {
 public:
+  __device__
+  HelmholtzOperatorQuad(Number coef)
+    : coef(coef)
+  {}
+
   __device__ void operator()(
     CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number>
       *                fe_eval,
     const unsigned int q_point) const;
+
+private:
+  Number coef;
 };
 
 
@@ -44,7 +67,7 @@ __device__ void HelmholtzOperatorQuad<dim, fe_degree, Number, n_q_points_1d>::
   CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number> *fe_eval,
   const unsigned int                                                    q) const
 {
-  fe_eval->submit_value(Number(10) * fe_eval->get_value(q), q);
+  fe_eval->submit_value(coef * fe_eval->get_value(q), q);
   fe_eval->submit_gradient(fe_eval->get_gradient(q), q);
 }
 
@@ -54,6 +77,10 @@ template <int dim, int fe_degree, typename Number, int n_q_points_1d>
 class HelmholtzOperator
 {
 public:
+  HelmholtzOperator(Number *coefficient)
+    : coef(coefficient)
+  {}
+
   __device__ void
   operator()(
     const unsigned int                                          cell,
@@ -67,6 +94,8 @@ public:
     dealii::Utilities::pow(fe_degree + 1, dim);
   static const unsigned int n_q_points =
     dealii::Utilities::pow(n_q_points_1d, dim);
+
+  Number *coef;
 };
 
 
@@ -80,18 +109,68 @@ operator()(const unsigned int                                          cell,
            const Number *                         src,
            Number *                               dst) const
 {
+  const unsigned int pos = CUDAWrappers::local_q_point_id<dim, Number>(
+    cell, gpu_data, n_dofs_1d, n_q_points);
+
   CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number> fe_eval(
     cell, gpu_data, shared_data);
   fe_eval.read_dof_values(src);
   fe_eval.evaluate(true, true);
   fe_eval.apply_quad_point_operations(
-    HelmholtzOperatorQuad<dim, fe_degree, Number, n_q_points_1d>());
+    HelmholtzOperatorQuad<dim, fe_degree, Number, n_q_points_1d>(coef[pos]));
   fe_eval.integrate(true, true);
   fe_eval.distribute_local_to_global(dst);
 }
 
 
 
+template <int dim, int fe_degree, typename Number, int n_q_points_1d>
+class VaryingCoefficientFunctor
+{
+public:
+  VaryingCoefficientFunctor(Number *coefficient)
+    : coef(coefficient)
+  {}
+
+  __device__ void
+  operator()(
+    const unsigned int                                          cell,
+    const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data);
+
+  static const unsigned int n_dofs_1d = fe_degree + 1;
+  static const unsigned int n_local_dofs =
+    dealii::Utilities::pow(fe_degree + 1, dim);
+  static const unsigned int n_q_points =
+    dealii::Utilities::pow(n_q_points_1d, dim);
+
+private:
+  Number *coef;
+};
+
+
+
+template <int dim, int fe_degree, typename Number, int n_q_points_1d>
+__device__ void
+VaryingCoefficientFunctor<dim, fe_degree, Number, n_q_points_1d>::
+operator()(const unsigned int                                          cell,
+           const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data)
+{
+  const unsigned int pos = CUDAWrappers::local_q_point_id<dim, Number>(
+    cell, gpu_data, n_dofs_1d, n_q_points);
+  const auto q_point =
+    CUDAWrappers::get_quadrature_point<dim, Number>(cell, gpu_data, n_dofs_1d);
+
+  Number p_square = 0.;
+  for (unsigned int i = 0; i < dim; ++i)
+    {
+      Number coord = q_point[i];
+      p_square += coord * coord;
+    }
+  coef[pos] = 10. / (0.05 + 2. * p_square);
+}
+
+
+
 template <int dim,
           int fe_degree,
           typename Number,
@@ -100,16 +179,42 @@ template <int dim,
 class MatrixFreeTest
 {
 public:
-  MatrixFreeTest(const CUDAWrappers::MatrixFree<dim, Number> &data_in)
-    : data(data_in){};
+  MatrixFreeTest(const CUDAWrappers::MatrixFree<dim, Number> &data_in,
+                 const unsigned int                           size,
+                 const bool constant_coeff = true);
 
   void
   vmult(VectorType &dst, const VectorType &src);
 
 private:
   const CUDAWrappers::MatrixFree<dim, Number> &data;
+  LinearAlgebra::CUDAWrappers::Vector<Number>  coef;
 };
 
+template <int dim,
+          int fe_degree,
+          typename Number,
+          typename VectorType,
+          int n_q_points_1d>
+MatrixFreeTest<dim, fe_degree, Number, VectorType, n_q_points_1d>::
+  MatrixFreeTest(const CUDAWrappers::MatrixFree<dim, Number> &data_in,
+                 const unsigned int                           size,
+                 const bool                                   constant_coeff)
+  : data(data_in)
+{
+  coef.reinit(size);
+  if (constant_coeff)
+    {
+      coef.add(10.);
+    }
+  else
+    {
+      VaryingCoefficientFunctor<dim, fe_degree, Number, n_q_points_1d> functor(
+        coef.get_values());
+      data.evaluate_coefficients(functor);
+    }
+}
+
 
 
 template <int dim,
@@ -123,7 +228,8 @@ MatrixFreeTest<dim, fe_degree, Number, VectorType, n_q_points_1d>::vmult(
   const VectorType &src)
 {
   dst = static_cast<Number>(0.);
-  HelmholtzOperator<dim, fe_degree, Number, n_q_points_1d> helmholtz_operator;
+  HelmholtzOperator<dim, fe_degree, Number, n_q_points_1d> helmholtz_operator(
+    coef.get_values());
   data.cell_loop(helmholtz_operator, src, dst);
   data.copy_constrained_values(src, dst);
 }