]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New test cases
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sun, 5 Jan 2020 17:21:09 +0000 (18:21 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sun, 5 Jan 2020 17:21:09 +0000 (18:21 +0100)
tests/matrix_free/inverse_mass_06.cc [new file with mode: 0644]
tests/matrix_free/inverse_mass_06.output [new file with mode: 0644]
tests/matrix_free/inverse_mass_07.cc [new file with mode: 0644]
tests/matrix_free/inverse_mass_07.output [new file with mode: 0644]

diff --git a/tests/matrix_free/inverse_mass_06.cc b/tests/matrix_free/inverse_mass_06.cc
new file mode 100644 (file)
index 0000000..0c16331
--- /dev/null
@@ -0,0 +1,228 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests CellwiseInverseMassMatrix::apply() with implicit definition from
+// FEEvaluationBase::JxW(), otherwise the same as inverse_mass_01
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/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 "../tests.h"
+
+
+
+template <int dim,
+          int fe_degree,
+          typename Number,
+          typename VectorType = Vector<Number>>
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+    : data(data_in){};
+
+  void
+  local_mass_operator(
+    const MatrixFree<dim, Number> &              data,
+    VectorType &                                 dst,
+    const VectorType &                           src,
+    const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, fe_degree, fe_degree + 1, 1, Number> fe_eval(data);
+    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(true, false);
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          fe_eval.submit_value(fe_eval.get_value(q), q);
+        fe_eval.integrate(true, false);
+        fe_eval.set_dof_values(dst);
+      }
+  }
+
+  void
+  local_inverse_mass_operator(
+    const MatrixFree<dim, Number> &              data,
+    VectorType &                                 dst,
+    const VectorType &                           src,
+    const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, fe_degree, fe_degree + 1, 1, Number> fe_eval(data);
+    MatrixFreeOperators::CellwiseInverseMassMatrix<dim, fe_degree, 1, Number>
+      mass_inv(fe_eval);
+
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        fe_eval.reinit(cell);
+        fe_eval.read_dof_values(src);
+        mass_inv.apply(fe_eval.begin_dof_values(), fe_eval.begin_dof_values());
+        fe_eval.set_dof_values(dst);
+      }
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    data.cell_loop(
+      &MatrixFreeTest<dim, fe_degree, Number, VectorType>::local_mass_operator,
+      this,
+      dst,
+      src);
+  };
+
+  void
+  apply_inverse(VectorType &dst, const VectorType &src) const
+  {
+    data.cell_loop(&MatrixFreeTest<dim, fe_degree, Number, VectorType>::
+                     local_inverse_mass_operator,
+                   this,
+                   dst,
+                   src);
+  };
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+do_test(const DoFHandler<dim> &dof)
+{
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  MappingQ<dim> mapping(4);
+
+  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::partition_color;
+    data.tasks_block_size = 3;
+    AffineConstraints<double> constraints;
+
+    mf_data.reinit(mapping, dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim, fe_degree, number> mf(mf_data);
+  Vector<number> in(dof.n_dofs()), inverse(dof.n_dofs()),
+    reference(dof.n_dofs());
+
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    {
+      const double entry = random_value<double>();
+      in(i)              = entry;
+    }
+
+  mf.apply_inverse(inverse, in);
+
+  SolverControl            control(1000, 1e-12);
+  SolverCG<Vector<number>> solver(control);
+  solver.solve(mf, reference, in, PreconditionIdentity());
+
+  inverse -= reference;
+  const double diff_norm = inverse.linfty_norm() / reference.linfty_norm();
+
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  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 (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+      if (cell->at_boundary(f))
+        cell->face(f)->set_all_manifold_ids(0);
+  tria.set_manifold(0, manifold);
+
+  if (dim < 3 || fe_degree < 2)
+    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 (; cell != endc; ++cell)
+    if (cell->center().norm() < 1e-8)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+
+  do_test<dim, fe_degree, double>(dof);
+
+  if (dim == 2)
+    {
+      deallog.push("float");
+      do_test<dim, fe_degree, float>(dof);
+      deallog.pop();
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  deallog << std::setprecision(3);
+  // do not output CG iteration numbers to log file because these are too
+  // sensitive
+  deallog.depth_file(2);
+
+  {
+    deallog.push("2d");
+    test<2, 1>();
+    test<2, 2>();
+    test<2, 4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1>();
+    test<3, 2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/inverse_mass_06.output b/tests/matrix_free/inverse_mass_06.output
new file mode 100644 (file)
index 0000000..19887e1
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL:2d::Testing FE_DGQ<2>(1)
+DEAL:2d::Norm of difference: 2.21e-13
+DEAL:2d::
+DEAL:2d::Testing FE_DGQ<2>(2)
+DEAL:2d::Norm of difference: 1.29e-13
+DEAL:2d::
+DEAL:2d::Testing FE_DGQ<2>(4)
+DEAL:2d::Norm of difference: 4.60e-14
+DEAL:2d::
+DEAL:3d::Testing FE_DGQ<3>(1)
+DEAL:3d::Norm of difference: 3.27e-14
+DEAL:3d::
+DEAL:3d::Testing FE_DGQ<3>(2)
+DEAL:3d::Norm of difference: 1.03e-15
+DEAL:3d::
diff --git a/tests/matrix_free/inverse_mass_07.cc b/tests/matrix_free/inverse_mass_07.cc
new file mode 100644 (file)
index 0000000..27243c0
--- /dev/null
@@ -0,0 +1,214 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests CellwiseInverseMassMatrix::apply() with implicit definition from
+// FEEvaluationBase::JxW() for vector DG elements, otherwise the same as
+// inverse_mass_02
+
+#include <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/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 "../tests.h"
+
+
+
+template <int dim,
+          int fe_degree,
+          typename Number,
+          typename VectorType = Vector<Number>>
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+    : data(data_in){};
+
+  void
+  local_mass_operator(
+    const MatrixFree<dim, Number> &              data,
+    VectorType &                                 dst,
+    const VectorType &                           src,
+    const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, fe_degree, fe_degree + 1, dim, Number> fe_eval(data);
+    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(true, false);
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          fe_eval.submit_value(fe_eval.get_value(q), q);
+        fe_eval.integrate(true, false);
+        fe_eval.set_dof_values(dst);
+      }
+  }
+
+  void
+  local_inverse_mass_operator(
+    const MatrixFree<dim, Number> &              data,
+    VectorType &                                 dst,
+    const VectorType &                           src,
+    const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, fe_degree, fe_degree + 1, dim, Number> fe_eval(data);
+    MatrixFreeOperators::CellwiseInverseMassMatrix<dim, fe_degree, dim, Number>
+                       mass_inv(fe_eval);
+    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);
+        mass_inv.apply(fe_eval.begin_dof_values(), fe_eval.begin_dof_values());
+        fe_eval.set_dof_values(dst);
+      }
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    data.cell_loop(
+      &MatrixFreeTest<dim, fe_degree, Number, VectorType>::local_mass_operator,
+      this,
+      dst,
+      src);
+  };
+
+  void
+  apply_inverse(VectorType &dst, const VectorType &src) const
+  {
+    data.cell_loop(&MatrixFreeTest<dim, fe_degree, Number, VectorType>::
+                     local_inverse_mass_operator,
+                   this,
+                   dst,
+                   src);
+  };
+
+private:
+  const MatrixFree<dim, Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+do_test(const DoFHandler<dim> &dof)
+{
+  deallog << "Testing " << dof.get_fe().get_name() << 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::partition_color;
+    data.tasks_block_size = 3;
+    AffineConstraints<double> constraints;
+
+    mf_data.reinit(dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim, fe_degree, number> mf(mf_data);
+  Vector<number> in(dof.n_dofs()), inverse(dof.n_dofs()),
+    reference(dof.n_dofs());
+
+  for (unsigned int i = 0; i < dof.n_dofs(); ++i)
+    {
+      const double entry = random_value<double>();
+      in(i)              = entry;
+    }
+
+  mf.apply_inverse(inverse, in);
+
+  SolverControl            control(1000, 1e-12);
+  SolverCG<Vector<number>> solver(control);
+  solver.solve(mf, reference, in, PreconditionIdentity());
+
+  inverse -= reference;
+  const double diff_norm = inverse.linfty_norm() / reference.linfty_norm();
+
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, -1, 1);
+  tria.refine_global(1);
+  if (dim < 3 || fe_degree < 2)
+    tria.refine_global(1);
+  tria.begin(tria.n_levels() - 1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  typename Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
+                                                    endc = tria.end();
+  for (; cell != endc; ++cell)
+    if (cell->center().norm() < 1e-8)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  const unsigned int degree = fe_degree;
+  FESystem<dim>      fe(FE_DGQ<dim>(degree), dim);
+  DoFHandler<dim>    dof(tria);
+  dof.distribute_dofs(fe);
+
+  do_test<dim, fe_degree, double>(dof);
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  deallog << std::setprecision(3);
+  deallog.depth_file(2);
+
+  {
+    deallog.push("2d");
+    test<2, 1>();
+    test<2, 2>();
+    test<2, 4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1>();
+    test<3, 2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/inverse_mass_07.output b/tests/matrix_free/inverse_mass_07.output
new file mode 100644 (file)
index 0000000..07daa53
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2]
+DEAL:2d::Norm of difference: 5.67e-16
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2]
+DEAL:2d::Norm of difference: 1.01e-14
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2]
+DEAL:2d::Norm of difference: 9.43e-16
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3]
+DEAL:3d::Norm of difference: 1.23e-15
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3]
+DEAL:3d::Norm of difference: 2.30e-15
+DEAL:3d::

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.