]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add matrix-free interpolation test for Raviart-Thomas case 15948/head
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Thu, 31 Aug 2023 06:08:38 +0000 (08:08 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Thu, 31 Aug 2023 08:10:13 +0000 (10:10 +0200)
tests/matrix_free/interpolate_rt.cc [new file with mode: 0644]
tests/matrix_free/interpolate_rt.output [new file with mode: 0644]

diff --git a/tests/matrix_free/interpolate_rt.cc b/tests/matrix_free/interpolate_rt.cc
new file mode 100644 (file)
index 0000000..2c0c17c
--- /dev/null
@@ -0,0 +1,403 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 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
+// operations in getting the function values, the function gradients, and the
+// function Laplacians on a cartesian mesh (hyper cube). This tests whether
+// cartesian meshes are treated correctly. The test case is without any
+// constraints
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.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/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+
+template <int dim>
+class CompareFunction : public Function<dim>
+{
+public:
+  CompareFunction()
+    : Function<dim>(dim)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component) const
+  {
+    double value = (1.2 - 0.1 * component) * p[0] + 0.4;
+    for (unsigned int d = 1; d < dim; ++d)
+      value -= (2.7 - 0.2 * component) * d * p[d];
+    return value;
+  }
+
+  virtual Tensor<1, dim>
+  gradient(const Point<dim> &p, const unsigned int component) const
+  {
+    Tensor<1, dim> grad;
+    grad[0] = 1.2 - 0.1 * component;
+    for (unsigned int d = 1; d < dim; ++d)
+      grad[d] = -(2.7 - 0.2 * component) * d;
+    return grad;
+  }
+};
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d = fe_degree + 1,
+          typename Number   = double>
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in)
+    : data(data_in){};
+
+  MatrixFreeTest(const MatrixFree<dim, Number> &data_in,
+                 const Mapping<dim>            &mapping)
+    : data(data_in){};
+
+  virtual ~MatrixFreeTest()
+  {}
+
+  // make function virtual to allow derived classes to define a different
+  // function
+  virtual void
+  cell(const MatrixFree<dim, Number> &data,
+       Vector<Number> &,
+       const Vector<Number>                        &src,
+       const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_eval(data);
+
+    CompareFunction<dim> function;
+
+    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 j = 0; j < data.n_active_entries_per_cell_batch(cell);
+             ++j)
+          for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+            {
+              ++cell_times;
+              Point<dim> p;
+              for (unsigned int d = 0; d < dim; ++d)
+                p[d] = fe_eval.quadrature_point(q)[d][j];
+              for (unsigned int d = 0; d < dim; ++d)
+                {
+                  cell_errors[0][d] +=
+                    std::abs(fe_eval.get_value(q)[d][j] - function.value(p, d));
+                  for (unsigned int e = 0; e < dim; ++e)
+                    cell_errors[1][d] +=
+                      std::abs(fe_eval.get_gradient(q)[d][e][j] -
+                               function.gradient(p, d)[e]);
+                }
+              double divergence = 0;
+              for (unsigned int d = 0; d < dim; ++d)
+                divergence += function.gradient(p, d)[d];
+              cell_errors[2][0] +=
+                std::abs(fe_eval.get_divergence(q)[j] - divergence);
+            }
+      }
+  }
+
+  virtual void
+  face(const MatrixFree<dim, Number> &data,
+       Vector<Number> &,
+       const Vector<Number>                        &src,
+       const std::pair<unsigned int, unsigned int> &face_range) const
+  {
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_evalm(data,
+                                                                          true);
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_evalp(
+      data, false);
+
+    CompareFunction<dim> function;
+
+    for (unsigned int face = face_range.first; face < face_range.second; ++face)
+      {
+        fe_evalm.reinit(face);
+        fe_evalm.read_dof_values(src);
+        fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+                          EvaluationFlags::hessians);
+        fe_evalp.reinit(face);
+        fe_evalp.read_dof_values(src);
+        fe_evalp.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+                          EvaluationFlags::hessians);
+
+        for (unsigned int j = 0; j < VectorizedArray<Number>::size(); ++j)
+          {
+            // skip empty components in VectorizedArray
+            if (data.get_face_info(face).cells_interior[j] ==
+                numbers::invalid_unsigned_int)
+              break;
+            for (unsigned int q = 0; q < fe_evalm.n_q_points; ++q)
+              {
+                ++facem_times;
+                ++facep_times;
+                Point<dim> p;
+
+                // interior face
+                for (unsigned int d = 0; d < dim; ++d)
+                  p[d] = fe_evalm.quadrature_point(q)[d][j];
+                for (unsigned int d = 0; d < dim; ++d)
+                  {
+                    facem_errors[0][d] += std::abs(fe_evalm.get_value(q)[d][j] -
+                                                   function.value(p, d));
+                    for (unsigned int e = 0; e < dim; ++e)
+                      {
+                        facem_errors[1][d] +=
+                          std::abs(fe_evalm.get_gradient(q)[d][e][j] -
+                                   function.gradient(p, d)[e]);
+                      }
+                  }
+                double divergence = 0;
+                for (unsigned int d = 0; d < dim; ++d)
+                  divergence += function.gradient(p, d)[d];
+                facem_errors[2][0] +=
+                  std::abs(fe_evalm.get_divergence(q)[j] - divergence);
+
+                // exterior face
+                for (unsigned int d = 0; d < dim; ++d)
+                  {
+                    facep_errors[0][d] += std::abs(fe_evalp.get_value(q)[d][j] -
+                                                   function.value(p, d));
+                    for (unsigned int e = 0; e < dim; ++e)
+                      facep_errors[1][d] +=
+                        std::abs(fe_evalp.get_gradient(q)[d][e][j] -
+                                 function.gradient(p, d)[e]);
+                  }
+                facep_errors[2][0] +=
+                  std::abs(fe_evalp.get_divergence(q)[j] - divergence);
+              }
+          }
+      }
+  }
+
+  virtual void
+  boundary(const MatrixFree<dim, Number> &data,
+           Vector<Number> &,
+           const Vector<Number>                        &src,
+           const std::pair<unsigned int, unsigned int> &face_range) const
+  {
+    FEFaceEvaluation<dim, fe_degree, n_q_points_1d, dim, Number> fe_evalm(data,
+                                                                          true);
+
+    CompareFunction<dim> function;
+
+    for (unsigned int face = face_range.first; face < face_range.second; ++face)
+      {
+        fe_evalm.reinit(face);
+        fe_evalm.read_dof_values(src);
+        fe_evalm.evaluate(EvaluationFlags::values | EvaluationFlags::gradients |
+                          EvaluationFlags::hessians);
+
+        for (unsigned int j = 0; j < VectorizedArray<Number>::size(); ++j)
+          {
+            // skip empty components in VectorizedArray
+            if (data.get_face_info(face).cells_interior[j] ==
+                numbers::invalid_unsigned_int)
+              break;
+            for (unsigned int q = 0; q < fe_evalm.n_q_points; ++q)
+              {
+                ++boundary_times;
+                Point<dim> p;
+                for (unsigned int d = 0; d < dim; ++d)
+                  p[d] = fe_evalm.quadrature_point(q)[d][j];
+                for (unsigned int d = 0; d < dim; ++d)
+                  {
+                    boundary_errors[0][d] += std::abs(
+                      fe_evalm.get_value(q)[d][j] - function.value(p, d));
+                    for (unsigned int e = 0; e < dim; ++e)
+                      boundary_errors[1][d] +=
+                        std::abs(fe_evalm.get_gradient(q)[d][e][j] -
+                                 function.gradient(p, d)[e]);
+                  }
+                double divergence = 0;
+                for (unsigned int d = 0; d < dim; ++d)
+                  divergence += function.gradient(p, d)[d];
+                boundary_errors[2][0] +=
+                  std::abs(fe_evalm.get_divergence(q)[j] - divergence);
+              }
+          }
+      }
+  }
+
+
+
+  static void
+  print_error(const std::string                     &text,
+              const dealii::ndarray<double, 3, dim> &array,
+              const unsigned long long               n_times)
+  {
+    deallog << "Error " << std::left << std::setw(6) << text << " values:     ";
+    for (unsigned int d = 0; d < dim; ++d)
+      deallog << array[0][d] / n_times << " ";
+    deallog << std::endl;
+    deallog << "Error " << std::left << std::setw(6) << text << " gradients:  ";
+    for (unsigned int d = 0; d < dim; ++d)
+      deallog << array[1][d] / n_times << " ";
+    deallog << std::endl;
+    deallog << "Error " << std::left << std::setw(6) << text << " divergence: ";
+    deallog << array[2][0] / n_times << " ";
+    deallog << std::endl;
+  }
+
+  void
+  test_functions(const Vector<Number> &src) const
+  {
+    for (unsigned int d = 0; d < dim; ++d)
+      for (unsigned int i = 0; i < 3; ++i)
+        {
+          cell_errors[i][d]     = 0;
+          facem_errors[i][d]    = 0;
+          facep_errors[i][d]    = 0;
+          boundary_errors[i][d] = 0;
+        }
+    cell_times = facem_times = facep_times = boundary_times = 0;
+
+    Vector<Number> dst_dummy;
+    data.loop(&MatrixFreeTest::cell,
+              &MatrixFreeTest::face,
+              &MatrixFreeTest::boundary,
+              this,
+              dst_dummy,
+              src);
+
+    print_error("cell", cell_errors, cell_times);
+    print_error("face-", facem_errors, facem_times);
+    print_error("face+", facep_errors, facep_times);
+    print_error("face b", boundary_errors, boundary_times);
+    deallog << std::endl;
+  };
+
+protected:
+  const MatrixFree<dim, Number>          &data;
+  mutable dealii::ndarray<double, 3, dim> cell_errors, facem_errors,
+    facep_errors, boundary_errors;
+  mutable unsigned long long cell_times, facem_times, facep_times,
+    boundary_times;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void
+do_test(const DoFHandler<dim>           &dof,
+        const AffineConstraints<double> &constraints)
+{
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  // use this for info on problem
+  // 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;
+
+  Vector<number> interpolated(dof.n_dofs());
+  VectorTools::interpolate(dof, CompareFunction<dim>(), interpolated);
+
+  constraints.distribute(interpolated);
+  MatrixFree<dim, number> mf_data;
+  {
+    const QGauss<1> quad(dof.get_fe().degree + 1);
+    typename MatrixFree<dim, number>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, number>::AdditionalData::none;
+    data.mapping_update_flags  = update_gradients | update_quadrature_points;
+    data.mapping_update_flags_boundary_faces =
+      update_gradients | update_quadrature_points;
+    data.mapping_update_flags_inner_faces =
+      update_gradients | update_quadrature_points;
+    mf_data.reinit(MappingQ1<dim>{}, dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim, fe_degree, fe_degree + 1, number> mf(mf_data);
+  mf.test_functions(interpolated);
+}
+
+
+
+template <int dim, int fe_degree>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, -1, 1);
+  // 1 refinement will make sure that the Piola transform is 1
+  tria.refine_global(1);
+
+  {
+    FE_RaviartThomasNodal<dim> fe(fe_degree - 1);
+    DoFHandler<dim>            dof(tria);
+    dof.distribute_dofs(fe);
+
+    AffineConstraints<double> constraints;
+    constraints.close();
+    if (fe_degree > 2)
+      do_test<dim, -1, double>(dof, constraints);
+    else
+      do_test<dim, fe_degree, double>(dof, constraints);
+  }
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  deallog << std::setprecision(5);
+  {
+    deallog.push("2d");
+    test<2, 1>();
+    test<2, 2>();
+    test<2, 3>();
+    test<2, 4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3, 1>();
+    test<3, 2>();
+    test<3, 3>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/interpolate_rt.output b/tests/matrix_free/interpolate_rt.output
new file mode 100644 (file)
index 0000000..05c86cb
--- /dev/null
@@ -0,0 +1,99 @@
+
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(0)
+DEAL:2d::Error cell   values:     0.77942 0.31754 
+DEAL:2d::Error cell   gradients:  2.7000 1.1000 
+DEAL:2d::Error cell   divergence: 2.2204e-16 
+DEAL:2d::Error face-  values:     1.0647 0.43377 
+DEAL:2d::Error face-  gradients:  2.7000 1.1000 
+DEAL:2d::Error face-  divergence: 2.2204e-16 
+DEAL:2d::Error face+  values:     1.0647 0.43377 
+DEAL:2d::Error face+  gradients:  2.7000 1.1000 
+DEAL:2d::Error face+  divergence: 2.2204e-16 
+DEAL:2d::Error face b values:     1.0647 0.43377 
+DEAL:2d::Error face b gradients:  2.7000 1.1000 
+DEAL:2d::Error face b divergence: 2.2204e-16 
+DEAL:2d::
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(1)
+DEAL:2d::Error cell   values:     7.8641e-17 1.7193e-16 
+DEAL:2d::Error cell   gradients:  3.6391e-16 3.5774e-16 
+DEAL:2d::Error cell   divergence: 3.3923e-16 
+DEAL:2d::Error face-  values:     3.2382e-17 8.7893e-17 
+DEAL:2d::Error face-  gradients:  2.2204e-16 4.6259e-16 
+DEAL:2d::Error face-  divergence: 4.4409e-16 
+DEAL:2d::Error face+  values:     3.2382e-17 8.7893e-17 
+DEAL:2d::Error face+  gradients:  4.6259e-16 3.1456e-16 
+DEAL:2d::Error face+  divergence: 4.6259e-16 
+DEAL:2d::Error face b values:     1.0640e-16 7.8641e-17 
+DEAL:2d::Error face b gradients:  5.8287e-16 5.3661e-16 
+DEAL:2d::Error face b divergence: 6.9389e-16 
+DEAL:2d::
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(2)
+DEAL:2d::Error cell   values:     2.6498e-16 1.9418e-16 
+DEAL:2d::Error cell   gradients:  1.9949e-15 1.3982e-15 
+DEAL:2d::Error cell   divergence: 1.9186e-15 
+DEAL:2d::Error face-  values:     1.8757e-16 1.8388e-16 
+DEAL:2d::Error face-  gradients:  3.2058e-15 1.3323e-15 
+DEAL:2d::Error face-  divergence: 3.3445e-15 
+DEAL:2d::Error face+  values:     9.9096e-17 1.2750e-16 
+DEAL:2d::Error face+  gradients:  2.2482e-15 2.3870e-15 
+DEAL:2d::Error face+  divergence: 3.4417e-15 
+DEAL:2d::Error face b values:     3.0358e-16 3.7643e-16 
+DEAL:2d::Error face b gradients:  4.0870e-15 2.9213e-15 
+DEAL:2d::Error face b divergence: 5.1972e-15 
+DEAL:2d::
+DEAL:2d::Testing FE_RaviartThomasNodal<2>(3)
+DEAL:2d::Error cell   values:     2.5709e-16 2.3020e-16 
+DEAL:2d::Error cell   gradients:  3.3573e-15 3.1575e-15 
+DEAL:2d::Error cell   divergence: 2.8955e-15 
+DEAL:2d::Error face-  values:     1.7764e-16 2.3315e-16 
+DEAL:2d::Error face-  gradients:  2.7867e-15 3.4972e-15 
+DEAL:2d::Error face-  divergence: 2.2982e-15 
+DEAL:2d::Error face+  values:     1.4502e-16 1.7902e-16 
+DEAL:2d::Error face+  gradients:  3.9191e-15 3.1086e-15 
+DEAL:2d::Error face+  divergence: 2.8533e-15 
+DEAL:2d::Error face b values:     4.0662e-16 2.8276e-16 
+DEAL:2d::Error face b gradients:  6.0452e-15 6.1173e-15 
+DEAL:2d::Error face b divergence: 5.0460e-15 
+DEAL:2d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(0)
+DEAL:3d::Error cell   values:     1.5588 1.4434 0.66395 
+DEAL:3d::Error cell   gradients:  8.1000 6.1000 3.3000 
+DEAL:3d::Error cell   divergence: 4.4409e-16 
+DEAL:3d::Error face-  values:     1.9392 1.7956 0.82597 
+DEAL:3d::Error face-  gradients:  8.1000 6.1000 3.3000 
+DEAL:3d::Error face-  divergence: 3.7007e-16 
+DEAL:3d::Error face+  values:     1.9392 1.7956 0.82597 
+DEAL:3d::Error face+  gradients:  8.1000 6.1000 3.3000 
+DEAL:3d::Error face+  divergence: 4.4409e-16 
+DEAL:3d::Error face b values:     1.9392 1.7956 0.82597 
+DEAL:3d::Error face b gradients:  8.1000 6.1000 3.3000 
+DEAL:3d::Error face b divergence: 4.0708e-16 
+DEAL:3d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(1)
+DEAL:3d::Error cell   values:     2.9323e-16 3.6905e-16 2.5109e-16 
+DEAL:3d::Error cell   gradients:  1.0465e-15 1.0043e-15 1.2002e-15 
+DEAL:3d::Error cell   divergence: 1.1801e-15 
+DEAL:3d::Error face-  values:     1.8478e-16 2.6355e-16 1.8799e-16 
+DEAL:3d::Error face-  gradients:  7.1342e-16 7.8744e-16 1.0352e-15 
+DEAL:3d::Error face-  divergence: 1.0938e-15 
+DEAL:3d::Error face+  values:     1.8478e-16 2.6355e-16 1.8799e-16 
+DEAL:3d::Error face+  gradients:  9.7864e-16 7.4015e-16 1.1164e-15 
+DEAL:3d::Error face+  divergence: 1.3158e-15 
+DEAL:3d::Error face b values:     2.9555e-16 3.3975e-16 2.6985e-16 
+DEAL:3d::Error face b gradients:  1.1472e-15 1.1030e-15 1.4006e-15 
+DEAL:3d::Error face b divergence: 1.6324e-15 
+DEAL:3d::
+DEAL:3d::Testing FE_RaviartThomasNodal<3>(2)
+DEAL:3d::Error cell   values:     8.0941e-16 8.5099e-16 8.2424e-16 
+DEAL:3d::Error cell   gradients:  5.5828e-15 5.0810e-15 4.6376e-15 
+DEAL:3d::Error cell   divergence: 5.7940e-15 
+DEAL:3d::Error face-  values:     8.3182e-16 7.5558e-16 7.5558e-16 
+DEAL:3d::Error face-  gradients:  6.7793e-15 5.4089e-15 4.6999e-15 
+DEAL:3d::Error face-  divergence: 7.8224e-15 
+DEAL:3d::Error face+  values:     5.1581e-16 4.8728e-16 6.1290e-16 
+DEAL:3d::Error face+  gradients:  5.9964e-15 5.3950e-15 5.5135e-15 
+DEAL:3d::Error face+  divergence: 9.3305e-15 
+DEAL:3d::Error face b values:     1.1222e-15 1.2358e-15 1.0999e-15 
+DEAL:3d::Error face b gradients:  8.5892e-15 7.4477e-15 6.4367e-15 
+DEAL:3d::Error face b divergence: 1.1618e-14 
+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.