]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add test 11212/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 20 Nov 2020 08:05:37 +0000 (09:05 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 20 Nov 2020 09:58:09 +0000 (10:58 +0100)
tests/matrix_free/matrix_vector_hp_no_template.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_hp_no_template.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_mf.h

diff --git a/tests/matrix_free/matrix_vector_hp_no_template.cc b/tests/matrix_free/matrix_vector_hp_no_template.cc
new file mode 100644 (file)
index 0000000..4fad45a
--- /dev/null
@@ -0,0 +1,227 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2020 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 for hp DoFHandler on a hyperball mesh with hanging nodes and finite
+// elements orders distributed randomly.
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/hp/dof_handler.h>
+#include <deal.II/hp/fe_values.h>
+
+#include "../tests.h"
+
+#include "matrix_vector_common.h"
+
+
+
+template <int dim, typename Number>
+class MatrixFreeTestHP
+{
+public:
+  MatrixFreeTestHP(const MatrixFree<dim, Number> &data_in)
+    : data(data_in){};
+
+  void
+  local_apply(const MatrixFree<dim, Number> &              data,
+              Vector<Number> &                             dst,
+              const Vector<Number> &                       src,
+              const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    for (unsigned int i = 1; i <= 9 - 2 * dim; ++i) // skip FE_Nothing
+      {
+        const auto cell_subrange =
+          data.create_cell_subrange_hp_by_index(cell_range, i);
+
+        if (cell_subrange.second <= cell_subrange.first)
+          continue;
+
+        helmholtz_operator_no_template(data, dst, src, cell_subrange, i, i);
+      }
+  }
+
+  void
+  vmult(Vector<Number> &dst, const Vector<Number> &src) const
+  {
+    dst = 0;
+    data.cell_loop(&MatrixFreeTestHP<dim, Number>::local_apply, this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  if (fe_degree > 1)
+    return;
+
+  typedef double               number;
+  const SphericalManifold<dim> manifold;
+  Triangulation<dim>           tria;
+  GridGenerator::hyper_ball(tria);
+  typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
+                                                    endc = tria.end();
+  for (; cell != endc; ++cell)
+    for (const unsigned int f : GeometryInfo<dim>::face_indices())
+      if (cell->at_boundary(f))
+        cell->face(f)->set_all_manifold_ids(0);
+  tria.set_manifold(0, manifold);
+  tria.refine_global(1);
+
+  // refine a few cells
+  for (unsigned int i = 0; i < 11 - 3 * dim; ++i)
+    {
+      typename Triangulation<dim>::active_cell_iterator cell =
+                                                          tria.begin_active(),
+                                                        endc = tria.end();
+      unsigned int counter                                   = 0;
+      for (; cell != endc; ++cell, ++counter)
+        if (counter % (7 - i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  const unsigned int max_degree = 9 - 2 * dim;
+
+  hp::FECollection<dim> fe_collection;
+  hp::QCollection<dim>  quadrature_collection;
+  hp::QCollection<1>    quadrature_collection_mf;
+
+  fe_collection.push_back(FE_Nothing<dim>());
+  quadrature_collection.push_back(QGauss<dim>(1));
+  quadrature_collection_mf.push_back(QGauss<1>(1));
+  for (unsigned int deg = 1; deg <= max_degree; ++deg)
+    {
+      fe_collection.push_back(FE_Q<dim>(QGaussLobatto<1>(deg + 1)));
+      quadrature_collection.push_back(QGauss<dim>(deg + 1));
+      quadrature_collection_mf.push_back(QGauss<1>(deg + 1));
+    }
+
+  hp::DoFHandler<dim> dof(tria);
+  // set the active FE index in a random order
+  {
+    typename hp::DoFHandler<dim>::active_cell_iterator cell =
+                                                         dof.begin_active(),
+                                                       endc = dof.end();
+    for (; cell != endc; ++cell)
+      {
+        const unsigned int fe_index = Testing::rand() % max_degree;
+        cell->set_active_fe_index(fe_index + 1);
+      }
+    // We cannot set random cells to FE_Nothing. We get the following error
+    // The violated condition was: dominating_fe.n_dofs_per_face(face) <=
+    // subface_fe.n_dofs_per_face(face)
+    cell = dof.begin_active();
+    cell->set_active_fe_index(0);
+  }
+
+  // setup DoFs
+  dof.distribute_dofs(fe_collection);
+  AffineConstraints<double> constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values(dof,
+                                           0,
+                                           Functions::ZeroFunction<dim>(),
+                                           constraints);
+  constraints.close();
+  DynamicSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+  DoFTools::make_sparsity_pattern(dof, csp, constraints, false);
+  SparsityPattern sparsity;
+  sparsity.copy_from(csp);
+  SparseMatrix<double> system_matrix(sparsity);
+
+  // std::cout << "Number of cells: " <<
+  // dof.get_triangulation().n_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;
+
+  // set up MatrixFree
+  MatrixFree<dim, number>                          mf_data;
+  typename MatrixFree<dim, number>::AdditionalData data;
+  data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+  mf_data.reinit(dof, constraints, quadrature_collection_mf, data);
+  MatrixFreeTestHP<dim, number> mf(mf_data);
+
+  // assemble sparse matrix with (\nabla v,
+  // \nabla u) + (v, 10 * u)
+  {
+    hp::FEValues<dim>                    hp_fe_values(fe_collection,
+                                   quadrature_collection,
+                                   update_values | update_gradients |
+                                     update_JxW_values);
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> local_dof_indices;
+
+    typename hp::DoFHandler<dim>::active_cell_iterator cell =
+                                                         dof.begin_active(),
+                                                       endc = dof.end();
+    for (; cell != endc; ++cell)
+      {
+        const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+
+        cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+        cell_matrix = 0;
+        hp_fe_values.reinit(cell);
+        const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values();
+
+        for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_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));
+            }
+        local_dof_indices.resize(dofs_per_cell);
+        cell->get_dof_indices(local_dof_indices);
+
+        constraints.distribute_local_to_global(cell_matrix,
+                                               local_dof_indices,
+                                               system_matrix);
+      }
+  }
+
+  // fill a right hand side vector with random
+  // numbers in unconstrained degrees of freedom
+  Vector<double> src(dof.n_dofs());
+  Vector<double> result_spmv(src), result_mf(src);
+
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i) == false)
+        src(i) = random_value<double>();
+    }
+
+  // now perform matrix-vector product and check
+  // its correctness
+  system_matrix.vmult(result_spmv, src);
+  mf.vmult(result_mf, src);
+
+  result_mf -= result_spmv;
+  const double diff_norm = result_mf.linfty_norm();
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
diff --git a/tests/matrix_free/matrix_vector_hp_no_template.output b/tests/matrix_free/matrix_vector_hp_no_template.output
new file mode 100644 (file)
index 0000000..7252835
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
index afc1ebf6bf07521c12f2ac6555d62e9a467c7094..8a77fccac98c73404cf3f90aee0ad7682e647865 100644 (file)
@@ -56,6 +56,38 @@ helmholtz_operator(const MatrixFree<dim, typename VectorType::value_type> &data,
 
 
 
+template <int dim, typename VectorType>
+void
+helmholtz_operator_no_template(
+  const MatrixFree<dim, typename VectorType::value_type> &data,
+  VectorType &                                            dst,
+  const VectorType &                                      src,
+  const std::pair<unsigned int, unsigned int> &           cell_range,
+  const unsigned int                                      active_fe_index,
+  const unsigned int                                      active_quad_index)
+{
+  typedef typename VectorType::value_type Number;
+  FEEvaluation<dim, -1, 0, 1, Number>     fe_eval(
+    data, 0, 0, 0, active_fe_index, active_quad_index);
+  const unsigned int n_q_points = fe_eval.n_q_points;
+
+  for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+    {
+      fe_eval.reinit(cell);
+      fe_eval.read_dof_values(src);
+      fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+      for (unsigned int q = 0; q < n_q_points; ++q)
+        {
+          fe_eval.submit_value(Number(10) * fe_eval.get_value(q), q);
+          fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+        }
+      fe_eval.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+      fe_eval.distribute_local_to_global(dst);
+    }
+}
+
+
+
 template <int dim,
           int fe_degree,
           typename Number,

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.