From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 17 Nov 2020 18:53:38 +0000 (+0100)
Subject: Add Poisson/Helmholtz test: MatrixFree + simplex mesh
X-Git-Tag: v9.3.0-rc1~858^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c66cf7ac57d765bc6a7b6fc114cfe5cadbcef442;p=dealii.git

Add Poisson/Helmholtz test: MatrixFree + simplex mesh
---

diff --git a/tests/simplex/matrix_free_01.cc b/tests/simplex/matrix_free_01.cc
new file mode 100644
index 0000000000..90741f8683
--- /dev/null
+++ b/tests/simplex/matrix_free_01.cc
@@ -0,0 +1,294 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Solve Poisson problem and Helmholtz problem on a simplex mesh with
+// continuous elements and compare results between matrix-free and matrix-based
+// implementations.
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+
+template <int dim>
+class PoissonOperator
+{
+public:
+  using VectorType = LinearAlgebra::distributed::Vector<double>;
+
+  PoissonOperator(const MatrixFree<dim, double> &matrix_free,
+                  const bool                     do_helmholtz)
+    : matrix_free(matrix_free)
+    , do_helmholtz(do_helmholtz)
+  {}
+
+  void
+  initialize_dof_vector(VectorType &vec)
+  {
+    matrix_free.initialize_dof_vector(vec);
+  }
+
+  void
+  rhs(VectorType &vec) const
+  {
+    const int dummy = 0;
+
+    matrix_free.template cell_loop<VectorType, int>(
+      [&](const auto &, auto &dst, const auto &, const auto cells) {
+        FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free);
+        for (unsigned int cell = cells.first; cell < cells.second; ++cell)
+          {
+            phi.reinit(cell);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              phi.submit_value(1.0, q);
+
+            phi.integrate_scatter(true, false, dst);
+          }
+      },
+      vec,
+      dummy,
+      true);
+  }
+
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    matrix_free.template cell_loop<VectorType, VectorType>(
+      [&](const auto &, auto &dst, const auto &src, const auto cells) {
+        FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free);
+        for (unsigned int cell = cells.first; cell < cells.second; ++cell)
+          {
+            phi.reinit(cell);
+            phi.gather_evaluate(src, do_helmholtz, true);
+
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              {
+                if (do_helmholtz)
+                  phi.submit_value(phi.get_value(q), q);
+
+                phi.submit_gradient(phi.get_gradient(q), q);
+              }
+
+            phi.integrate_scatter(do_helmholtz, true, dst);
+          }
+      },
+      dst,
+      src,
+      true);
+  }
+
+private:
+  const MatrixFree<dim, double> &matrix_free;
+  const bool                     do_helmholtz;
+};
+
+template <int dim>
+void
+test(const unsigned int degree, const bool do_helmholtz)
+{
+  Triangulation<dim> tria;
+
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 16 : 8);
+
+  Simplex::FE_P<dim>   fe(degree);
+  Simplex::QGauss<dim> quad(degree + 1);
+  MappingFE<dim>       mapping(Simplex::FE_P<dim>(1));
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  VectorTools::interpolate_boundary_values(
+    mapping, dof_handler, 0, Functions::ZeroFunction<dim>(), constraints);
+  constraints.close();
+
+  const auto solve_and_postprocess =
+    [&](const auto &poisson_operator,
+        auto &      x,
+        auto &      b) -> std::pair<unsigned int, double> {
+    ReductionControl reduction_control;
+    SolverCG<typename std::remove_reference<decltype(x)>::type> solver(
+      reduction_control);
+    solver.solve(poisson_operator, x, b, PreconditionIdentity());
+
+    if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+      printf("Solved in %d iterations.\n", reduction_control.last_step());
+
+    constraints.distribute(x);
+
+#if 0
+  DataOut<dim> data_out;
+  data_out.attach_dof_handler(dof_handler);
+  x.update_ghost_values();
+  data_out.add_data_vector(dof_handler, x, "solution");
+  data_out.build_patches(mapping, 2);
+  data_out.write_vtu_with_pvtu_record("./", "result", 0, MPI_COMM_WORLD);
+#endif
+
+    return {reduction_control.last_step(), reduction_control.last_value()};
+  };
+
+  const auto mf_algo = [&]() {
+    typename MatrixFree<dim, double>::AdditionalData additional_data;
+    additional_data.mapping_update_flags = update_gradients | update_values;
+
+    MatrixFree<dim, double> matrix_free;
+    matrix_free.reinit(
+      mapping, dof_handler, constraints, quad, additional_data);
+
+    PoissonOperator<dim> poisson_operator(matrix_free, do_helmholtz);
+
+    LinearAlgebra::distributed::Vector<double> x, b;
+    poisson_operator.initialize_dof_vector(x);
+    poisson_operator.initialize_dof_vector(b);
+
+    poisson_operator.rhs(b);
+
+    return solve_and_postprocess(poisson_operator, x, b);
+  };
+
+  const auto mb_algo = [&]() {
+    Vector<double> x, b;
+
+    x.reinit(dof_handler.n_dofs());
+    b.reinit(dof_handler.n_dofs());
+
+    SparseMatrix<double> A;
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
+
+    SparsityPattern sparsity_pattern;
+    sparsity_pattern.copy_from(dsp);
+    A.reinit(sparsity_pattern);
+
+    const auto flags = update_values | update_gradients | update_JxW_values;
+
+    FEValues<dim> fe_values(mapping, fe, quad, flags);
+
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      {
+        if (cell->is_locally_owned() == false)
+          continue;
+
+        fe_values.reinit(cell);
+
+        const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+        cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+        cell_rhs.reinit(dofs_per_cell);
+
+        for (const auto q : fe_values.quadrature_point_indices())
+          {
+            for (const auto i : fe_values.dof_indices())
+              for (const auto j : fe_values.dof_indices())
+                cell_matrix(i, j) += (fe_values.shape_grad(i, q) *        //
+                                        fe_values.shape_grad(j, q) +      //
+                                      static_cast<double>(do_helmholtz) * //
+                                        fe_values.shape_value(i, q) *     //
+                                        fe_values.shape_value(j, q)) *    //
+                                     fe_values.JxW(q);                    //
+
+            for (const unsigned int i : fe_values.dof_indices())
+              cell_rhs(i) += (fe_values.shape_value(i, q) * //
+                              1. *                          //
+                              fe_values.JxW(q));            //
+          }
+
+        local_dof_indices.resize(cell->get_fe().dofs_per_cell);
+        cell->get_dof_indices(local_dof_indices);
+
+        constraints.distribute_local_to_global(
+          cell_matrix, cell_rhs, local_dof_indices, A, b);
+      }
+
+    return solve_and_postprocess(A, x, b);
+  };
+
+  const auto compare = [&](const auto result_mf, const auto result_mb) {
+    AssertDimension(result_mf.first, result_mb.first);
+    Assert(std::abs(result_mf.second - result_mb.second) < 1e-8,
+           ExcNotImplemented());
+
+    deallog << "dim=" << dim << " ";
+    deallog << "degree=" << degree << " ";
+    deallog << "Type=";
+
+    if (do_helmholtz)
+      deallog << "Helmholtz"
+              << " : ";
+    else
+      deallog << "Possion  "
+              << " : ";
+
+    deallog << "Convergence step " << result_mf.first << " value "
+            << result_mf.second << "." << std::endl;
+  };
+
+  compare(mf_algo(), mb_algo());
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  deallog.depth_file(1);
+
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+
+  test<2>(/*degree=*/1, /*do_helmholtz*/ false);
+  test<2>(/*degree=*/1, /*do_helmholtz*/ true);
+  test<2>(/*degree=*/2, /*do_helmholtz*/ false);
+  test<2>(/*degree=*/2, /*do_helmholtz*/ true);
+
+  test<3>(/*degree=*/1, /*do_helmholtz*/ false);
+  test<3>(/*degree=*/1, /*do_helmholtz*/ true);
+  test<3>(/*degree=*/2, /*do_helmholtz*/ false);
+  test<3>(/*degree=*/2, /*do_helmholtz*/ true);
+}
diff --git a/tests/simplex/matrix_free_01.with_simplex_support=on.output b/tests/simplex/matrix_free_01.with_simplex_support=on.output
new file mode 100644
index 0000000000..ebff9c041c
--- /dev/null
+++ b/tests/simplex/matrix_free_01.with_simplex_support=on.output
@@ -0,0 +1,9 @@
+
+DEAL::dim=2 degree=1 Type=Possion   : Convergence step 15 value 0.000430255.
+DEAL::dim=2 degree=1 Type=Helmholtz : Convergence step 14 value 0.000558792.
+DEAL::dim=2 degree=2 Type=Possion   : Convergence step 36 value 0.000319106.
+DEAL::dim=2 degree=2 Type=Helmholtz : Convergence step 36 value 0.000294440.
+DEAL::dim=3 degree=1 Type=Possion   : Convergence step 7 value 0.000281643.
+DEAL::dim=3 degree=1 Type=Helmholtz : Convergence step 7 value 0.000268188.
+DEAL::dim=3 degree=2 Type=Possion   : Convergence step 20 value 0.000146315.
+DEAL::dim=3 degree=2 Type=Helmholtz : Convergence step 20 value 0.000137731.