]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add tests
authorBuğrahan Temür <bugrahan.temuer@tum.de>
Wed, 15 Mar 2023 14:47:15 +0000 (15:47 +0100)
committerBuğrahan Temür <bugrahan.temuer@tum.de>
Wed, 15 Mar 2023 14:47:15 +0000 (15:47 +0100)
- inverse_mass_10: tests cellwise inverse mass with random dyadic coefficients
- inverse_mass_11: same as 10, but with fe_degree=-1

Both tests use a GMRES solver (coefficients not symmetric) with a reduced accuracy (1e-10) compared to other tests in this test suite, which use 1e-12. Reason is to make the tests more robust as the problem can become poorly conditioned due to the random coefficient.

tests/matrix_free/inverse_mass_10.cc [new file with mode: 0644]
tests/matrix_free/inverse_mass_10.output [new file with mode: 0644]
tests/matrix_free/inverse_mass_11.cc [new file with mode: 0644]
tests/matrix_free/inverse_mass_11.output [new file with mode: 0644]

diff --git a/tests/matrix_free/inverse_mass_10.cc b/tests/matrix_free/inverse_mass_10.cc
new file mode 100644 (file)
index 0000000..8bd9ab6
--- /dev/null
@@ -0,0 +1,259 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 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 with random dyadic coefficients, otherwise
+// the same as inverse_mass_02.cc
+
+#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_gmres.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)
+  {
+    FEEvaluation<dim, fe_degree, fe_degree + 1, dim, Number> fe_eval(data);
+    const unsigned int n_q_points = fe_eval.n_q_points;
+
+    dyadic_coefficients.resize(n_q_points);
+    inverse_dyadic_coefficients.resize(n_q_points);
+    create_and_invert_dyadic_coefficients();
+  };
+
+  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(EvaluationFlags::values);
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            const auto value_flux =
+              dyadic_coefficients[q] * fe_eval.get_value(q);
+            fe_eval.submit_value(value_flux, q);
+          }
+        fe_eval.integrate(EvaluationFlags::values);
+        fe_eval.distribute_local_to_global(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;
+    AlignedVector<VectorizedArray<Number>> inverse_JxW_values(n_q_points);
+    AlignedVector<Tensor<2, dim, VectorizedArray<Number>>>
+      inverse_coefficients_with_JxW(n_q_points);
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        fe_eval.reinit(cell);
+        mass_inv.fill_inverse_JxW_values(inverse_JxW_values);
+        for (unsigned int q = 0; q < inverse_JxW_values.size(); ++q)
+          inverse_coefficients_with_JxW[q] =
+            inverse_JxW_values[q] * inverse_dyadic_coefficients[q];
+        fe_eval.read_dof_values(src);
+        mass_inv.apply(inverse_coefficients_with_JxW,
+                       fe_eval.begin_dof_values(),
+                       fe_eval.begin_dof_values());
+        fe_eval.distribute_local_to_global(dst);
+      }
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    dst = 0;
+    data.cell_loop(
+      &MatrixFreeTest<dim, fe_degree, Number, VectorType>::local_mass_operator,
+      this,
+      dst,
+      src);
+  };
+
+  void
+  apply_inverse(VectorType &dst, const VectorType &src) const
+  {
+    dst = 0;
+    data.cell_loop(&MatrixFreeTest<dim, fe_degree, Number, VectorType>::
+                     local_inverse_mass_operator,
+                   this,
+                   dst,
+                   src);
+  };
+
+private:
+  void
+  create_and_invert_dyadic_coefficients()
+  {
+    for (unsigned int q = 0; q < dyadic_coefficients.size(); ++q)
+      {
+        for (unsigned int d1 = 0; d1 < dim; ++d1)
+          {
+            for (unsigned int d2 = 0; d2 < dim; ++d2)
+              {
+                dyadic_coefficients[q][d1][d2] =
+                  make_vectorized_array(random_value<Number>(1.0e-2, 1.0e-1));
+                // make diagonal dominant because inverse is needed
+                if (d1 == d2)
+                  dyadic_coefficients[q][d1][d2] *= 1.0e2;
+              }
+          }
+        inverse_dyadic_coefficients[q] = invert(dyadic_coefficients[q]);
+      }
+  }
+
+  AlignedVector<Tensor<2, dim, VectorizedArray<Number>>> dyadic_coefficients;
+  AlignedVector<Tensor<2, dim, VectorizedArray<Number>>>
+                                 inverse_dyadic_coefficients;
+  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(MappingQ1<dim>{}, 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-10);
+  // GMRES solver as dyadic coefficients are not necessarily symmetric
+  SolverGMRES<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_10.output b/tests/matrix_free/inverse_mass_10.output
new file mode 100644 (file)
index 0000000..295b3ce
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2]
+DEAL:2d::Norm of difference: 4.55e-15
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2]
+DEAL:2d::Norm of difference: 2.16e-11
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2]
+DEAL:2d::Norm of difference: 4.98e-11
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3]
+DEAL:3d::Norm of difference: 4.22e-11
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3]
+DEAL:3d::Norm of difference: 2.12e-11
+DEAL:3d::
diff --git a/tests/matrix_free/inverse_mass_11.cc b/tests/matrix_free/inverse_mass_11.cc
new file mode 100644 (file)
index 0000000..4407c2e
--- /dev/null
@@ -0,0 +1,260 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 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 with random dyadic coefficients. Uses
+// implicit template argument fe_degree=-1 in FEEvaluation, otherwise the same
+// as inverse_mass_10.cc
+
+#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_gmres.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)
+  {
+    FEEvaluation<dim, -1, 0, dim, Number> fe_eval(data);
+    const unsigned int                    n_q_points = fe_eval.n_q_points;
+
+    dyadic_coefficients.resize(n_q_points);
+    inverse_dyadic_coefficients.resize(n_q_points);
+    create_and_invert_dyadic_coefficients();
+  };
+
+  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, -1, 0, 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(EvaluationFlags::values);
+        for (unsigned int q = 0; q < n_q_points; ++q)
+          {
+            const auto value_flux =
+              dyadic_coefficients[q] * fe_eval.get_value(q);
+            fe_eval.submit_value(value_flux, q);
+          }
+        fe_eval.integrate(EvaluationFlags::values);
+        fe_eval.distribute_local_to_global(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, -1, 0, dim, Number> fe_eval(data);
+    MatrixFreeOperators::CellwiseInverseMassMatrix<dim, -1, dim, Number>
+                                           mass_inv(fe_eval);
+    const unsigned int                     n_q_points = fe_eval.n_q_points;
+    AlignedVector<VectorizedArray<Number>> inverse_JxW_values(n_q_points);
+    AlignedVector<Tensor<2, dim, VectorizedArray<Number>>>
+      inverse_coefficients_with_JxW(n_q_points);
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        fe_eval.reinit(cell);
+        mass_inv.fill_inverse_JxW_values(inverse_JxW_values);
+        for (unsigned int q = 0; q < inverse_JxW_values.size(); ++q)
+          inverse_coefficients_with_JxW[q] =
+            inverse_JxW_values[q] * inverse_dyadic_coefficients[q];
+        fe_eval.read_dof_values(src);
+        mass_inv.apply(inverse_coefficients_with_JxW,
+                       fe_eval.begin_dof_values(),
+                       fe_eval.begin_dof_values());
+        fe_eval.distribute_local_to_global(dst);
+      }
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    dst = 0;
+    data.cell_loop(
+      &MatrixFreeTest<dim, fe_degree, Number, VectorType>::local_mass_operator,
+      this,
+      dst,
+      src);
+  };
+
+  void
+  apply_inverse(VectorType &dst, const VectorType &src) const
+  {
+    dst = 0;
+    data.cell_loop(&MatrixFreeTest<dim, fe_degree, Number, VectorType>::
+                     local_inverse_mass_operator,
+                   this,
+                   dst,
+                   src);
+  };
+
+private:
+  void
+  create_and_invert_dyadic_coefficients()
+  {
+    for (unsigned int q = 0; q < dyadic_coefficients.size(); ++q)
+      {
+        for (unsigned int d1 = 0; d1 < dim; ++d1)
+          {
+            for (unsigned int d2 = 0; d2 < dim; ++d2)
+              {
+                dyadic_coefficients[q][d1][d2] =
+                  make_vectorized_array(random_value<Number>(1.0e-2, 1.0e-1));
+                // make diagonal dominant because inverse is needed
+                if (d1 == d2)
+                  dyadic_coefficients[q][d1][d2] *= 1.0e2;
+              }
+          }
+        inverse_dyadic_coefficients[q] = invert(dyadic_coefficients[q]);
+      }
+  }
+
+  AlignedVector<Tensor<2, dim, VectorizedArray<Number>>> dyadic_coefficients;
+  AlignedVector<Tensor<2, dim, VectorizedArray<Number>>>
+                                 inverse_dyadic_coefficients;
+  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(MappingQ1<dim>{}, 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-10);
+  // GMRES solver as dyadic coefficients are not necessarily symmetric
+  SolverGMRES<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_11.output b/tests/matrix_free/inverse_mass_11.output
new file mode 100644 (file)
index 0000000..295b3ce
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2]
+DEAL:2d::Norm of difference: 4.55e-15
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2]
+DEAL:2d::Norm of difference: 2.16e-11
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2]
+DEAL:2d::Norm of difference: 4.98e-11
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3]
+DEAL:3d::Norm of difference: 4.22e-11
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3]
+DEAL:3d::Norm of difference: 2.12e-11
+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.