From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 28 Mar 2017 10:11:49 +0000 (+0200)
Subject: Fix bug in matrix-free evaluation kernel for degree 9 in 2D
X-Git-Tag: v8.5.0-rc1~1^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9af342b61e0a2de8067c7b44a6901194f907f19a;p=dealii.git

Fix bug in matrix-free evaluation kernel for degree 9 in 2D
---

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 15c0ac4b08..62695fea38 100644
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -157,14 +157,14 @@ namespace internal
         temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
                                  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
       }
-    else if (temp_size > 100)
+    else if (temp_size < 100)
       {
-        temp1 = scratch_data;
+        temp1 = &temp_data[0];
         temp2 = temp1 + temp_size;
       }
     else
       {
-        temp1 = &temp_data[0];
+        temp1 = scratch_data;
         temp2 = temp1 + temp_size;
       }
 
@@ -376,14 +376,14 @@ namespace internal
         temp2 = temp1 + std::max(Utilities::fixed_power<dim>(shape_info.fe_degree+1),
                                  Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
       }
-    else if (temp_size > 100)
+    else if (temp_size < 100)
       {
-        temp1 = scratch_data;
+        temp1 = &temp_data[0];
         temp2 = temp1 + temp_size;
       }
     else
       {
-        temp1 = &temp_data[0];
+        temp1 = scratch_data;
         temp2 = temp1 + temp_size;
       }
 
diff --git a/tests/matrix_free/matrix_vector_large_degree.cc b/tests/matrix_free/matrix_vector_large_degree_01.cc
similarity index 71%
rename from tests/matrix_free/matrix_vector_large_degree.cc
rename to tests/matrix_free/matrix_vector_large_degree_01.cc
index b327043175..d3b603a774 100644
--- a/tests/matrix_free/matrix_vector_large_degree.cc
+++ b/tests/matrix_free/matrix_vector_large_degree_01.cc
@@ -54,15 +54,6 @@ void test ()
 
   deallog << "Testing " << dof.get_fe().get_name() << std::endl;
 
-  MatrixFree<dim,double> mf_data;
-  {
-    const QGauss<1> quad (fe_degree+2);
-    typename MatrixFree<dim,double>::AdditionalData data;
-    data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::none;
-    mf_data.reinit (dof, constraints, quad, data);
-  }
-
-  MatrixFreeTest<dim,fe_degree,double,Vector<double>,fe_degree+2> mf (mf_data);
   Vector<double> in (dof.n_dofs()), out (dof.n_dofs());
 
   for (unsigned int i=0; i<dof.n_dofs(); ++i)
@@ -72,8 +63,32 @@ void test ()
       in(i) = 1.;
     }
 
-  mf.vmult (out, in);
-  deallog << "Result norm: " << out.l2_norm() << std::endl;
+  {
+    MatrixFree<dim,double> mf_data;
+    {
+      const QGauss<1> quad (fe_degree+1);
+      typename MatrixFree<dim,double>::AdditionalData data;
+      data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::none;
+      mf_data.reinit (dof, constraints, quad, data);
+    }
+
+    MatrixFreeTest<dim,fe_degree,double,Vector<double>,fe_degree+1> mf (mf_data);
+    mf.vmult (out, in);
+    deallog << "Result norm: " << out.l2_norm() << std::endl;
+  }
+  {
+    MatrixFree<dim,double> mf_data;
+    {
+      const QGauss<1> quad (fe_degree+2);
+      typename MatrixFree<dim,double>::AdditionalData data;
+      data.tasks_parallel_scheme = MatrixFree<dim,double>::AdditionalData::none;
+      mf_data.reinit (dof, constraints, quad, data);
+    }
+
+    MatrixFreeTest<dim,fe_degree,double,Vector<double>,fe_degree+2> mf (mf_data);
+    mf.vmult (out, in);
+    deallog << "Result norm: " << out.l2_norm() << std::endl;
+  }
 }
 
 
diff --git a/tests/matrix_free/matrix_vector_large_degree.output b/tests/matrix_free/matrix_vector_large_degree_01.output
similarity index 65%
rename from tests/matrix_free/matrix_vector_large_degree.output
rename to tests/matrix_free/matrix_vector_large_degree_01.output
index dc5e78dc50..24cc4aeca1 100644
--- a/tests/matrix_free/matrix_vector_large_degree.output
+++ b/tests/matrix_free/matrix_vector_large_degree_01.output
@@ -1,11 +1,16 @@
 
 DEAL:2d::Testing FE_Q<2>(5)
 DEAL:2d::Result norm: 1.14444
+DEAL:2d::Result norm: 1.14444
 DEAL:2d::Testing FE_Q<2>(15)
 DEAL:2d::Result norm: 0.399128
+DEAL:2d::Result norm: 0.399128
 DEAL:2d::Testing FE_Q<2>(35)
 DEAL:2d::Result norm: 0.173861
+DEAL:2d::Result norm: 0.173861
 DEAL:3d::Testing FE_Q<3>(10)
 DEAL:3d::Result norm: 0.143708
+DEAL:3d::Result norm: 0.143708
 DEAL:3d::Testing FE_Q<3>(20)
 DEAL:3d::Result norm: 0.0523289
+DEAL:3d::Result norm: 0.0523289
diff --git a/tests/matrix_free/matrix_vector_large_degree_02.cc b/tests/matrix_free/matrix_vector_large_degree_02.cc
new file mode 100644
index 0000000000..e88d90e79e
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_large_degree_02.cc
@@ -0,0 +1,173 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests matrix-vector product for larger polynomial degrees in 2D, otherwise
+// similar to matrix_vector_07.cc
+
+#include "../tests.h"
+
+#include "matrix_vector_mf.h"
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/function.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  typedef double number;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(1);
+  typename Triangulation<dim>::active_cell_iterator
+  cell = tria.begin_active (),
+  endc = tria.end();
+  cell = tria.begin_active ();
+  for (; cell!=endc; ++cell)
+    if (cell->is_locally_owned())
+      if (cell->center().norm()<0.2)
+        cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global(1);
+  tria.begin(tria.n_levels()-1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+
+  IndexSet owned_set = dof.locally_owned_dofs();
+  IndexSet relevant_set;
+  DoFTools::extract_locally_relevant_dofs (dof, relevant_set);
+
+  ConstraintMatrix constraints (relevant_set);
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction<dim>(),
+                                            constraints);
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  //std::cout << "Number of cells: " << tria.n_global_active_cells() << std::endl;
+  //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+  //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
+
+  MatrixFree<dim,number> mf_data;
+  {
+    const QGauss<1> quad (fe_degree+1);
+    typename MatrixFree<dim,number>::AdditionalData data;
+    data.tasks_parallel_scheme =
+      MatrixFree<dim,number>::AdditionalData::none;
+    mf_data.reinit (dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim,fe_degree,number,Vector<number> > mf (mf_data);
+  Vector<number> in(dof.n_dofs()), out(dof.n_dofs()), ref(dof.n_dofs());
+
+  for (unsigned int i=0; i<in.size(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      in(i) = (double)Testing::rand()/RAND_MAX;
+    }
+
+  mf.vmult (out, in);
+
+
+  // assemble trilinos sparse matrix with
+  // (\nabla v, \nabla u) + (v, 10 * u) for
+  // reference
+  SparsityPattern sparsity;
+  SparseMatrix<number> sparse_matrix;
+  {
+    DynamicSparsityPattern csp (dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
+    sparsity.copy_from(csp);
+    sparse_matrix.reinit (sparsity);
+  }
+  {
+    QGauss<dim>  quadrature_formula(fe_degree+1);
+
+    FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_JxW_values);
+
+    const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof.begin_active(),
+    endc = dof.end();
+    for (; cell!=endc; ++cell)
+      {
+        cell_matrix = 0;
+        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));
+            }
+
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global (cell_matrix,
+                                                local_dof_indices,
+                                                sparse_matrix);
+      }
+  }
+
+  sparse_matrix.vmult (ref, in);
+  out -= ref;
+  const double diff_norm = out.linfty_norm();
+
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+int main (int argc, char **argv)
+{
+  initlog();
+  deallog.threshold_double(1e-10);
+  test<2,9>();
+}
diff --git a/tests/matrix_free/matrix_vector_large_degree_02.output b/tests/matrix_free/matrix_vector_large_degree_02.output
new file mode 100644
index 0000000000..17a908ebde
--- /dev/null
+++ b/tests/matrix_free/matrix_vector_large_degree_02.output
@@ -0,0 +1,4 @@
+
+DEAL::Testing FE_Q<2>(9)
+DEAL::Norm of difference: 0
+DEAL::