From: Svenja Schoeder <schoeder@lnm.mw.tum.de>
Date: Tue, 21 Aug 2018 12:24:37 +0000 (+0200)
Subject: bugfix for braces, added test for vectorization mask
X-Git-Tag: v9.1.0-rc1~734^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=719bbf3f2c6d8f8a80d464308d194244183a62cd;p=dealii.git

bugfix for braces, added test for vectorization mask
---

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 46c90653bc..7b509d085e 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -3986,12 +3986,12 @@ FEEvaluationBase<dim, n_components_, Number, is_face>::read_write_operation(
 
   // Case 3: standard operation with one index per degree of freedom -> go on
   // here
-
   constexpr unsigned int n_vectorization =
     VectorizedArray<Number>::n_array_elements;
   Assert(mask.count() == n_vectorization,
          ExcNotImplemented("Masking currently not implemented for "
                            "non-contiguous DoF storage"));
+
   const unsigned int dofs_per_component =
     this->data->dofs_per_component_on_cell;
   if (dof_info->index_storage_variants
@@ -4572,24 +4572,29 @@ FEEvaluationBase<dim, n_components_, Number, is_face>::
               contiguous)
           {
             if (n_components == 1 || n_fe_components == 1)
-              for (unsigned int v = 0; v < vectorization_populated; ++v)
-                if (mask[v] == true)
-                  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
-                       ++i)
-                    operation.process_dof(dof_indices[v] + i,
-                                          *src[comp],
-                                          values_dofs[comp][i][v]);
-                else
-                  for (unsigned int v = 0; v < vectorization_populated; ++v)
-                    if (mask[v] == true)
-                      for (unsigned int i = 0;
-                           i < data->dofs_per_component_on_cell;
-                           ++i)
-                        operation.process_dof(
-                          dof_indices[v] + i +
-                            comp * data->dofs_per_component_on_cell,
-                          *src[0],
-                          values_dofs[comp][i][v]);
+              {
+                for (unsigned int v = 0; v < vectorization_populated; ++v)
+                  if (mask[v] == true)
+                    for (unsigned int i = 0;
+                         i < data->dofs_per_component_on_cell;
+                         ++i)
+                      operation.process_dof(dof_indices[v] + i,
+                                            *src[comp],
+                                            values_dofs[comp][i][v]);
+              }
+            else
+              {
+                for (unsigned int v = 0; v < vectorization_populated; ++v)
+                  if (mask[v] == true)
+                    for (unsigned int i = 0;
+                         i < data->dofs_per_component_on_cell;
+                         ++i)
+                      operation.process_dof(
+                        dof_indices[v] + i +
+                          comp * data->dofs_per_component_on_cell,
+                        *src[0],
+                        values_dofs[comp][i][v]);
+              }
           }
         else
           {
@@ -4601,24 +4606,29 @@ FEEvaluationBase<dim, n_components_, Number, is_face>::
                                VectorizedArray<Number>::n_array_elements + 1);
             if (n_components == 1 || n_fe_components == 1)
               for (unsigned int v = 0; v < vectorization_populated; ++v)
-                if (mask[v] == true)
-                  for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
-                       ++i)
-                    operation.process_dof(dof_indices[v] + i * offsets[v],
-                                          *src[comp],
-                                          values_dofs[comp][i][v]);
-                else
-                  for (unsigned int v = 0; v < vectorization_populated; ++v)
-                    if (mask[v] == true)
-                      for (unsigned int i = 0;
-                           i < data->dofs_per_component_on_cell;
-                           ++i)
-                        operation.process_dof(
-                          dof_indices[v] +
-                            (i + comp * data->dofs_per_component_on_cell) *
-                              offsets[v],
-                          *src[0],
-                          values_dofs[comp][i][v]);
+                {
+                  if (mask[v] == true)
+                    for (unsigned int i = 0;
+                         i < data->dofs_per_component_on_cell;
+                         ++i)
+                      operation.process_dof(dof_indices[v] + i * offsets[v],
+                                            *src[comp],
+                                            values_dofs[comp][i][v]);
+                }
+            else
+              {
+                for (unsigned int v = 0; v < vectorization_populated; ++v)
+                  if (mask[v] == true)
+                    for (unsigned int i = 0;
+                         i < data->dofs_per_component_on_cell;
+                         ++i)
+                      operation.process_dof(
+                        dof_indices[v] +
+                          (i + comp * data->dofs_per_component_on_cell) *
+                            offsets[v],
+                        *src[0],
+                        values_dofs[comp][i][v]);
+              }
           }
       }
 }
diff --git a/tests/matrix_free/advect_1d_vectorization_mask.cc b/tests/matrix_free/advect_1d_vectorization_mask.cc
new file mode 100644
index 0000000000..4706713b2f
--- /dev/null
+++ b/tests/matrix_free/advect_1d_vectorization_mask.cc
@@ -0,0 +1,353 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// similar to  advect_1d but testing the mask for distribute_local_to_global
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.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 <deal.II/matrix_free/operators.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d   = fe_degree + 1,
+          typename number     = double,
+          typename VectorType = Vector<number>,
+          int n_components    = 1>
+class MatrixFreeAdvectionBasic
+{
+public:
+  MatrixFreeAdvectionBasic(const MatrixFree<dim, number> &data,
+                           const bool         zero_within_loop       = true,
+                           const unsigned int start_vector_component = 0)
+    : data(data)
+    , zero_within_loop(zero_within_loop)
+    , start_vector_component(start_vector_component)
+  {
+    for (unsigned int d = 0; d < dim; ++d)
+      advection[d] = 0.4 + 0.12 * d;
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    if (!zero_within_loop)
+      dst = 0;
+    data.loop(&MatrixFreeAdvectionBasic::local_apply,
+              &MatrixFreeAdvectionBasic::local_apply_face,
+              &MatrixFreeAdvectionBasic::local_apply_boundary_face,
+              this,
+              dst,
+              src,
+              zero_within_loop,
+              MatrixFree<dim, number>::DataAccessOnFaces::values,
+              MatrixFree<dim, number>::DataAccessOnFaces::values);
+
+    FEEvaluation<dim, fe_degree, fe_degree + 1, n_components, number> phi(data);
+
+    const unsigned int dofs_per_cell = phi.dofs_per_cell;
+
+    AlignedVector<VectorizedArray<number>> coefficients(phi.dofs_per_cell);
+    MatrixFreeOperators::
+      CellwiseInverseMassMatrix<dim, fe_degree, n_components, number>
+        inverse(phi);
+
+    for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
+      {
+        phi.reinit(cell);
+        phi.read_dof_values(dst);
+
+        inverse.fill_inverse_JxW_values(coefficients);
+        inverse.apply(coefficients,
+                      n_components,
+                      phi.begin_dof_values(),
+                      phi.begin_dof_values());
+
+        phi.set_dof_values(dst);
+      }
+  }
+
+private:
+  void
+  local_apply(const MatrixFree<dim, number> &              data,
+              VectorType &                                 dst,
+              const VectorType &                           src,
+              const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, number> phi(
+      data, 0, 0, start_vector_component);
+
+    const unsigned int n_vect = VectorizedArray<number>::n_array_elements;
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        phi.reinit(cell);
+        phi.read_dof_values(src);
+        phi.evaluate(true, false);
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          phi.submit_gradient(advection * phi.get_value(q), q);
+        phi.integrate(false, true);
+        for (unsigned int v = 0; v < n_vect; ++v)
+          {
+            std::bitset<n_vect> mask;
+            mask[v] = true;
+            phi.distribute_local_to_global(dst, 0, mask);
+          }
+      }
+  }
+
+  void
+  local_apply_face(
+    const MatrixFree<dim, number> &              data,
+    VectorType &                                 dst,
+    const VectorType &                           src,
+    const std::pair<unsigned int, unsigned int> &face_range) const
+  {
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, number> phi_m(
+      data, true, 0, 0, start_vector_component);
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, number> phi_p(
+      data, false, 0, 0, start_vector_component);
+    typedef typename FEFaceEvaluation<dim,
+                                      fe_degree,
+                                      n_q_points_1d,
+                                      n_components,
+                                      number>::value_type value_type;
+
+    const unsigned int n_vect = VectorizedArray<number>::n_array_elements;
+
+    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(true, false);
+        phi_p.reinit(face);
+        phi_p.read_dof_values(src);
+        phi_p.evaluate(true, false);
+
+        for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+          {
+            value_type u_minus = phi_m.get_value(q),
+                       u_plus  = phi_p.get_value(q);
+            const VectorizedArray<number> normal_times_advection =
+              advection * phi_m.get_normal_vector(q);
+            const value_type flux_times_normal =
+              make_vectorized_array<number>(0.5) *
+              ((u_minus + u_plus) * normal_times_advection +
+               std::abs(normal_times_advection) * (u_minus - u_plus));
+            phi_m.submit_value(-flux_times_normal, q);
+            phi_p.submit_value(flux_times_normal, q);
+          }
+
+        phi_m.integrate(true, false);
+        for (unsigned int v = 0; v < n_vect; ++v)
+          {
+            std::bitset<n_vect> mask;
+            mask[v] = true;
+            phi_m.distribute_local_to_global(dst, 0, mask);
+          }
+        phi_p.integrate(true, false);
+        for (unsigned int v = 0; v < n_vect; ++v)
+          {
+            std::bitset<n_vect> mask;
+            mask[v] = true;
+            phi_p.distribute_local_to_global(dst, 0, mask);
+          }
+      }
+  }
+
+  void
+  local_apply_boundary_face(
+    const MatrixFree<dim, number> &              data,
+    VectorType &                                 dst,
+    const VectorType &                           src,
+    const std::pair<unsigned int, unsigned int> &face_range) const
+  {
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components, number>
+                                                          fe_eval(data, true, 0, 0, start_vector_component);
+    typedef typename FEFaceEvaluation<dim,
+                                      fe_degree,
+                                      n_q_points_1d,
+                                      n_components,
+                                      number>::value_type value_type;
+
+    const unsigned int n_vect = VectorizedArray<number>::n_array_elements;
+
+    for (unsigned int face = face_range.first; face < face_range.second; face++)
+      {
+        fe_eval.reinit(face);
+        fe_eval.read_dof_values(src);
+        fe_eval.evaluate(true, false);
+
+        for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+          {
+            value_type                    u_minus = fe_eval.get_value(q);
+            value_type                    u_plus  = -u_minus;
+            const VectorizedArray<number> normal_times_advection =
+              advection * fe_eval.get_normal_vector(q);
+            const value_type flux_times_normal =
+              make_vectorized_array<number>(0.5) *
+              ((u_minus + u_plus) * normal_times_advection +
+               std::abs(normal_times_advection) * (u_minus - u_plus));
+            fe_eval.submit_value(-flux_times_normal, q);
+          }
+
+        fe_eval.integrate(true, false);
+        for (unsigned int v = 0; v < n_vect; ++v)
+          {
+            std::bitset<n_vect> mask =
+              std::bitset<VectorizedArray<number>::n_array_elements>();
+            mask[v] = true;
+            fe_eval.distribute_local_to_global(dst, 0, mask);
+          }
+      }
+  }
+
+  const MatrixFree<dim, number> &         data;
+  const bool                              zero_within_loop;
+  const unsigned int                      start_vector_component;
+  Tensor<1, dim, VectorizedArray<number>> advection;
+};
+
+
+
+template <int dim>
+class AnalyticFunction : public Function<dim>
+{
+public:
+  static_assert(dim == 1, "Only 1D implemented");
+  AnalyticFunction()
+    : Function<dim>(1)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int) const override
+  {
+    return std::sin(3 * numbers::PI * p[0] / 0.8);
+  }
+};
+
+
+
+template <int dim>
+class AnalyticDerivative : public Function<dim>
+{
+public:
+  static_assert(dim == 1, "Only 1D implemented");
+  AnalyticDerivative()
+    : Function<dim>(1)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int) const override
+  {
+    Tensor<1, dim> advection;
+    for (unsigned int d = 0; d < dim; ++d)
+      advection[d] = 0.4 + 0.12 * d;
+
+    return -std::cos(3 * numbers::PI * p[0] / 0.8) * advection[0] * 3 *
+           numbers::PI / 0.8;
+  }
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test(const unsigned int n_refine)
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, 0, 0.8);
+  tria.refine_global(n_refine);
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  if (n_refine == 3)
+    {
+      deallog << "Testing " << dof.get_fe().get_name();
+      deallog << std::endl;
+    }
+
+  LinearAlgebra::distributed::Vector<double> in, out;
+
+  const QGauss<1>                                  quad(fe_degree + 1);
+  typename MatrixFree<dim, double>::AdditionalData data;
+  data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none;
+  data.mapping_update_flags_inner_faces =
+    (update_gradients | update_JxW_values);
+  data.mapping_update_flags_boundary_faces =
+    (update_gradients | update_JxW_values);
+
+  MatrixFree<dim, double> mf_data;
+  mf_data.reinit(dof, constraints, quad, data);
+
+  mf_data.initialize_dof_vector(in);
+  mf_data.initialize_dof_vector(out);
+
+  VectorTools::interpolate(dof, AnalyticFunction<dim>(), in);
+
+  MatrixFreeAdvectionBasic<dim,
+                           fe_degree,
+                           fe_degree + 1,
+                           double,
+                           LinearAlgebra::distributed::Vector<double>>
+    mf2(mf_data);
+  mf2.vmult(out, in);
+
+  VectorTools::interpolate(dof, AnalyticDerivative<dim>(), in);
+  out -= in;
+
+  double diff_norm = out.linfty_norm();
+  deallog << "Norm of difference:          " << diff_norm << " " << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  for (unsigned int r = 3; r < 9; ++r)
+    test<1, 2>(r);
+
+  for (unsigned int r = 3; r < 9; ++r)
+    test<1, 4>(r);
+
+  for (unsigned int r = 3; r < 9; ++r)
+    test<1, 5>(r);
+}
diff --git a/tests/matrix_free/advect_1d_vectorization_mask.output b/tests/matrix_free/advect_1d_vectorization_mask.output
new file mode 100644
index 0000000000..b213520047
--- /dev/null
+++ b/tests/matrix_free/advect_1d_vectorization_mask.output
@@ -0,0 +1,22 @@
+
+DEAL::Testing FE_DGQ<1>(2)
+DEAL::Norm of difference:          0.531539 
+DEAL::Norm of difference:          0.135407 
+DEAL::Norm of difference:          0.0340112 
+DEAL::Norm of difference:          0.00851279 
+DEAL::Norm of difference:          0.00212882 
+DEAL::Norm of difference:          0.000532245 
+DEAL::Testing FE_DGQ<1>(4)
+DEAL::Norm of difference:          0.00529591 
+DEAL::Norm of difference:          0.000336014 
+DEAL::Norm of difference:          2.10801e-05 
+DEAL::Norm of difference:          1.31874e-06 
+DEAL::Norm of difference:          8.24413e-08 
+DEAL::Norm of difference:          5.15331e-09 
+DEAL::Testing FE_DGQ<1>(5)
+DEAL::Norm of difference:          0.000346887 
+DEAL::Norm of difference:          1.09982e-05 
+DEAL::Norm of difference:          3.44938e-07 
+DEAL::Norm of difference:          1.07881e-08 
+DEAL::Norm of difference:          3.35085e-10 
+DEAL::Norm of difference:          9.72689e-12