]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add test 17159/head
authorNils Much <nils.much@tum.de>
Fri, 21 Jun 2024 10:55:58 +0000 (12:55 +0200)
committerNils Much <nils.much@tum.de>
Fri, 21 Jun 2024 10:55:58 +0000 (12:55 +0200)
tests/matrix_free/compute_diagonal_09.cc [new file with mode: 0644]
tests/matrix_free/compute_diagonal_09.output [new file with mode: 0644]

diff --git a/tests/matrix_free/compute_diagonal_09.cc b/tests/matrix_free/compute_diagonal_09.cc
new file mode 100644 (file)
index 0000000..60480e1
--- /dev/null
@@ -0,0 +1,287 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2022 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// Similar to compute_diagonal_08 but with different dof_no and quad_no.
+
+#include <deal.II/fe/fe_dgq.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test()
+{
+  const int fe_degree       = 3;
+  const int n_points        = fe_degree + 1;
+  const int n_components    = 1;
+  using Number              = double;
+  using VectorizedArrayType = VectorizedArray<Number>;
+  using VectorType          = Vector<Number>;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball(tria);
+  tria.refine_global(0);
+
+  const FE_DGQ<dim> fe_q(fe_degree);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe_q);
+
+  AffineConstraints<Number> constraints;
+  constraints.close();
+
+  const unsigned int                         dof_index       = 1;
+  const std::vector<const DoFHandler<dim> *> dof_handler_vec = {&dof_handler,
+                                                                &dof_handler};
+  const std::vector<const AffineConstraints<Number> *> constraints_vec = {
+    &constraints, &constraints};
+
+  typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+    additional_data;
+  additional_data.mapping_update_flags = update_values | update_gradients;
+  additional_data.mapping_update_flags_inner_faces =
+    update_values | update_gradients;
+  additional_data.mapping_update_flags_boundary_faces =
+    update_values | update_gradients;
+
+  MappingQ<dim> mapping(1);
+  QGauss<1>     quad(fe_degree + 1);
+
+  MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+  matrix_free.reinit(
+    mapping, dof_handler_vec, constraints_vec, quad, additional_data);
+
+  const auto cell_operation = [&](auto &phi) {
+    phi.evaluate(EvaluationFlags::gradients);
+    for (unsigned int q = 0; q < phi.n_q_points; ++q)
+      phi.submit_gradient(phi.get_gradient(q), q);
+    phi.integrate(EvaluationFlags::gradients);
+  };
+
+  const auto face_operation = [&](auto &phi_m, auto &phi_p) {
+    phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+    phi_p.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+
+    VectorizedArrayType sigmaF =
+      (std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) +
+       std::abs(
+         (phi_m.normal_vector(0) * phi_p.inverse_jacobian(0))[dim - 1])) *
+      (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
+
+    for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+      {
+        VectorizedArrayType average_value =
+          (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
+        VectorizedArrayType average_valgrad =
+          phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q);
+        average_valgrad = average_value * 2. * sigmaF - average_valgrad * 0.5;
+        phi_m.submit_normal_derivative(-average_value, q);
+        phi_p.submit_normal_derivative(-average_value, q);
+        phi_m.submit_value(average_valgrad, q);
+        phi_p.submit_value(-average_valgrad, q);
+      }
+    phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+    phi_p.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+  };
+
+  const auto boundary_operation = [&](auto &phi_m) {
+    phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+    VectorizedArrayType sigmaF =
+      std::abs((phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) *
+      Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
+
+    for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+      {
+        VectorizedArrayType average_value   = phi_m.get_value(q);
+        VectorizedArrayType average_valgrad = -phi_m.get_normal_derivative(q);
+        average_valgrad += average_value * sigmaF * 2.0;
+        phi_m.submit_normal_derivative(-average_value, q);
+        phi_m.submit_value(average_valgrad, q);
+      }
+
+    phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+  };
+
+  const auto vmult = [&](auto &dst, const auto &src) {
+    matrix_free.template loop<VectorType, VectorType>(
+      [&](
+        const auto &matrix_free, auto &dst, const auto &src, const auto range) {
+        FEEvaluation<dim,
+                     fe_degree,
+                     n_points,
+                     n_components,
+                     Number,
+                     VectorizedArrayType>
+          phi(matrix_free, dof_index);
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            phi.reinit(cell);
+            phi.read_dof_values(src);
+            cell_operation(phi);
+            phi.set_dof_values(dst);
+          }
+      },
+      [&](
+        const auto &matrix_free, auto &dst, const auto &src, const auto range) {
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_points,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi_m(matrix_free, true, dof_index);
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_points,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi_p(matrix_free, false, dof_index);
+
+        for (unsigned int face = range.first; face < range.second; ++face)
+          {
+            phi_m.reinit(face);
+            phi_p.reinit(face);
+            phi_m.read_dof_values(src);
+            phi_p.read_dof_values(src);
+            face_operation(phi_m, phi_p);
+            phi_m.distribute_local_to_global(dst);
+            phi_p.distribute_local_to_global(dst);
+          }
+      },
+      [&](const auto &matrix_free,
+          auto       &dst,
+          const auto &src,
+          const auto  face_range) {
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_points,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi_m(matrix_free, true, dof_index);
+
+        for (unsigned int face = face_range.first; face < face_range.second;
+             face++)
+          {
+            phi_m.reinit(face);
+            phi_m.read_dof_values(src);
+            boundary_operation(phi_m);
+            phi_m.distribute_local_to_global(dst);
+          }
+      },
+      dst,
+      src,
+      true);
+  };
+
+  // Compute matrix and diagonal (brute-force)
+  FullMatrix<Number> full_matrix(dof_handler.n_dofs(), dof_handler.n_dofs());
+  VectorType         diagonal(dof_handler.n_dofs());
+
+  for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
+    {
+      VectorType src(dof_handler.n_dofs());
+      VectorType dst(dof_handler.n_dofs());
+
+      src[i] = 1.0;
+
+      vmult(dst, src);
+
+      diagonal[i] = dst[i];
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        full_matrix[j][i] = dst[j];
+    }
+
+
+  // Compute matrix via MatrixFreeTools
+  FullMatrix<Number> full_matrix_mft(dof_handler.n_dofs(),
+                                     dof_handler.n_dofs());
+
+  MatrixFreeTools::compute_matrix<dim,
+                                  fe_degree,
+                                  n_points,
+                                  n_components,
+                                  Number,
+                                  VectorizedArrayType>(matrix_free,
+                                                       constraints,
+                                                       full_matrix_mft,
+                                                       cell_operation,
+                                                       face_operation,
+                                                       boundary_operation,
+                                                       dof_index);
+
+  for (unsigned int i = 0; i < full_matrix.m(); ++i)
+    for (unsigned int j = 0; j < full_matrix.n(); ++j)
+      if (std::abs(full_matrix[i][j] - full_matrix_mft[i][j]) > 1e-6)
+        {
+          full_matrix.print_formatted(
+            deallog.get_file_stream(), 3, true, 0, "0.000e+00");
+          deallog << std::endl;
+
+          full_matrix_mft.print_formatted(
+            deallog.get_file_stream(), 3, true, 0, "0.000e+00");
+          deallog << std::endl;
+        }
+
+  // Compute diagonal via MatrixFreeTools
+  VectorType diagonal_mft;
+  matrix_free.initialize_dof_vector(diagonal_mft);
+
+  MatrixFreeTools::compute_diagonal<dim,
+                                    fe_degree,
+                                    n_points,
+                                    n_components,
+                                    Number,
+                                    VectorizedArrayType>(matrix_free,
+                                                         diagonal_mft,
+                                                         cell_operation,
+                                                         face_operation,
+                                                         boundary_operation,
+                                                         dof_index);
+
+  for (unsigned int i = 0; i < diagonal.size(); ++i)
+    if (std::abs(diagonal[i] - diagonal_mft[i]) > 1e-6)
+      {
+        diagonal.print(deallog.get_file_stream());
+        deallog << std::endl;
+
+        diagonal_mft.print(deallog.get_file_stream());
+        deallog << std::endl;
+      }
+
+  deallog << "OK!" << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  test<2>();
+}
diff --git a/tests/matrix_free/compute_diagonal_09.output b/tests/matrix_free/compute_diagonal_09.output
new file mode 100644 (file)
index 0000000..5055485
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL:0::OK!

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.