]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add Poisson test: MatrixFree + DG + simplex mesh 11235/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 28 Nov 2020 19:16:07 +0000 (20:16 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 28 Nov 2020 19:16:07 +0000 (20:16 +0100)
tests/simplex/matrix_free_03.cc [new file with mode: 0644]
tests/simplex/matrix_free_03.with_simplex_support=on.output [new file with mode: 0644]

diff --git a/tests/simplex/matrix_free_03.cc b/tests/simplex/matrix_free_03.cc
new file mode 100644 (file)
index 0000000..5eee662
--- /dev/null
@@ -0,0 +1,611 @@
+// ---------------------------------------------------------------------
+//
+// 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_dgq.h>
+#include <deal.II/fe/fe_interface_values.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/meshworker/copy_data.h>
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.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;
+
+const double PENALTY = 8;
+
+
+template <int dim>
+class SmoothSolution : public Function<dim>
+{
+public:
+  SmoothSolution()
+    : Function<dim>()
+  {}
+  virtual void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<double> &          values,
+             const unsigned int             component = 0) const override;
+};
+
+template <int dim>
+void
+SmoothSolution<dim>::value_list(const std::vector<Point<dim>> &points,
+                                std::vector<double> &          values,
+                                const unsigned int /*component*/) const
+{
+  for (unsigned int i = 0; i < values.size(); ++i)
+    values[i] = 0.0;
+}
+
+template <int dim>
+class SmoothRightHandSide : public Function<dim>
+{
+public:
+  SmoothRightHandSide()
+    : Function<dim>()
+  {}
+  virtual void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<double> &          values,
+             const unsigned int /*component*/ = 0) const override;
+};
+
+template <int dim>
+void
+SmoothRightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
+                                     std::vector<double> &          values,
+                                     const unsigned int /*component*/) const
+{
+  for (unsigned int i = 0; i < values.size(); ++i)
+    values[i] = 1.0;
+}
+
+
+template <int dim>
+class PoissonOperator
+{
+public:
+  using VectorType = LinearAlgebra::distributed::Vector<double>;
+  using number     = double;
+
+  PoissonOperator(const MatrixFree<dim, double> &matrix_free)
+    : matrix_free(matrix_free)
+  {}
+
+  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 &data, auto &dst, const auto &, const auto cells) {
+        FEEvaluation<dim, -1, 0, 1, double> phi(data);
+        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);
+  }
+
+  const int fe_degree = 5; /*TODO*/
+
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    matrix_free.template loop<VectorType, VectorType>(
+      [&](const auto &data, auto &dst, const auto &src, const auto cell_range) {
+        FEEvaluation<dim, -1, 0, 1, number> phi(data);
+
+        for (unsigned int cell = cell_range.first; cell < cell_range.second;
+             ++cell)
+          {
+            phi.reinit(cell);
+            phi.read_dof_values(src);
+            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);
+            phi.set_dof_values(dst);
+          }
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval(data, true);
+        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval_neighbor(data, false);
+
+        for (unsigned int face = face_range.first; face < face_range.second;
+             face++)
+          {
+            fe_eval.reinit(face);
+            fe_eval_neighbor.reinit(face);
+
+            fe_eval.gather_evaluate(src,
+                                    EvaluationFlags::values |
+                                      EvaluationFlags::gradients);
+            fe_eval_neighbor.gather_evaluate(src,
+                                             EvaluationFlags::values |
+                                               EvaluationFlags::gradients);
+            VectorizedArray<number> sigmaF = PENALTY;
+            //  (std::abs((fe_eval.get_normal_vector(0) *
+            //             fe_eval.inverse_jacobian(0))[dim - 1]) +
+            //   std::abs((fe_eval.get_normal_vector(0) *
+            //             fe_eval_neighbor.inverse_jacobian(0))[dim - 1])) *
+            //  (number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
+
+            for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+              {
+                VectorizedArray<number> average_value =
+                  (fe_eval.get_value(q) - fe_eval_neighbor.get_value(q)) * 0.5;
+                VectorizedArray<number> average_valgrad =
+                  fe_eval.get_normal_derivative(q) +
+                  fe_eval_neighbor.get_normal_derivative(q);
+                average_valgrad =
+                  average_value * 2. * sigmaF - average_valgrad * 0.5;
+                fe_eval.submit_normal_derivative(-average_value, q);
+                fe_eval_neighbor.submit_normal_derivative(-average_value, q);
+                fe_eval.submit_value(average_valgrad, q);
+                fe_eval_neighbor.submit_value(-average_valgrad, q);
+              }
+            fe_eval.integrate_scatter(EvaluationFlags::values |
+                                        EvaluationFlags::gradients,
+                                      dst);
+            fe_eval_neighbor.integrate_scatter(EvaluationFlags::values |
+                                                 EvaluationFlags::gradients,
+                                               dst);
+          }
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval(data, true);
+        for (unsigned int face = face_range.first; face < face_range.second;
+             face++)
+          {
+            fe_eval.reinit(face);
+            fe_eval.read_dof_values(src);
+            fe_eval.evaluate(EvaluationFlags::values |
+                             EvaluationFlags::gradients);
+            VectorizedArray<number> sigmaF = PENALTY;
+            //  std::abs((fe_eval.get_normal_vector(0) *
+            //            fe_eval.inverse_jacobian(0))[dim - 1]) *
+            //  number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
+
+            for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+              {
+                VectorizedArray<number> average_value = fe_eval.get_value(q);
+                VectorizedArray<number> average_valgrad =
+                  -fe_eval.get_normal_derivative(q);
+                average_valgrad += average_value * sigmaF;
+                fe_eval.submit_normal_derivative(-average_value, q);
+                fe_eval.submit_value(average_valgrad, q);
+              }
+
+            fe_eval.integrate_scatter(EvaluationFlags::values |
+                                        EvaluationFlags::gradients,
+                                      dst);
+          }
+      },
+      dst,
+      src);
+  }
+
+private:
+  const MatrixFree<dim, double> &matrix_free;
+};
+
+
+struct CopyDataFace
+{
+  FullMatrix<double>                   cell_matrix;
+  std::vector<types::global_dof_index> joint_dof_indices;
+  std::array<double, 2>                values;
+  std::array<unsigned int, 2>          cell_indices;
+};
+
+struct CopyData
+{
+  FullMatrix<double>                   cell_matrix;
+  Vector<double>                       cell_rhs;
+  std::vector<types::global_dof_index> local_dof_indices;
+  std::vector<CopyDataFace>            face_data;
+  double                               value;
+  unsigned int                         cell_index;
+  template <class Iterator>
+  void
+  reinit(const Iterator &cell, unsigned int dofs_per_cell)
+  {
+    cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+    cell_rhs.reinit(dofs_per_cell);
+    local_dof_indices.resize(dofs_per_cell);
+    cell->get_dof_indices(local_dof_indices);
+  }
+};
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+  Triangulation<dim> tria;
+
+#if true
+  unsigned int n_subdivisions = 0;
+
+  if (dim == 2 && degree == 1)
+    n_subdivisions = 16;
+
+  if (dim == 2 && degree == 2)
+    n_subdivisions = 8;
+
+  if (dim == 3 && degree == 1)
+    n_subdivisions = 8;
+
+  if (dim == 3 && degree == 2)
+    n_subdivisions = 4;
+
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, n_subdivisions);
+
+  Simplex::FE_DGP<dim>     fe(degree);
+  Simplex::QGauss<dim>     quadrature(degree + 1);
+  Simplex::QGauss<dim - 1> face_quadrature(degree + 1);
+  MappingFE<dim>           mapping(Simplex::FE_P<dim>(1));
+#else
+  GridGenerator::subdivided_hyper_cube(tria, dim == 2 ? 16 : 8);
+
+  FE_DGQ<dim>    fe(degree);
+  QGauss<dim>    quadrature(degree + 1);
+  MappingFE<dim> mapping(FE_Q<dim>(1));
+
+#endif
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  const auto solve_and_postprocess =
+    [&](const auto &poisson_operator,
+        auto &      x,
+        auto &      b) -> std::pair<unsigned int, double> {
+    ReductionControl reduction_control(1000, 1e-7, 1e-3);
+    SolverCG<typename std::remove_reference<decltype(x)>::type> solver(
+      reduction_control);
+
+    try
+      {
+        solver.solve(poisson_operator, x, b, PreconditionIdentity());
+      }
+    catch (const std::exception &e)
+      {
+        deallog << e.what() << std::endl;
+      }
+
+    if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+      printf("Solved in %d iterations.\n", reduction_control.last_step());
+
+    constraints.distribute(x);
+
+#if 1
+    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
+
+    Vector<double> difference(tria.n_active_cells());
+
+    deallog << "dim=" << dim << " ";
+    deallog << "degree=" << degree << " ";
+
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      x,
+                                      Functions::ZeroFunction<dim>(),
+                                      difference,
+                                      quadrature,
+                                      VectorTools::NormType::L2_norm);
+
+    deallog << VectorTools::compute_global_error(tria,
+                                                 difference,
+                                                 VectorTools::NormType::L2_norm)
+            << std::endl;
+
+    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;
+    additional_data.mapping_update_flags_inner_faces =
+      update_gradients | update_values;
+    additional_data.mapping_update_flags_boundary_faces =
+      update_gradients | update_values;
+    additional_data.tasks_parallel_scheme =
+      MatrixFree<dim, double>::AdditionalData::none;
+
+    MatrixFree<dim, double> matrix_free;
+    matrix_free.reinit(
+      mapping, dof_handler, constraints, quadrature, additional_data);
+
+    PoissonOperator<dim> poisson_operator(matrix_free);
+
+    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> solution, system_rhs;
+
+    solution.reinit(dof_handler.n_dofs());
+    system_rhs.reinit(dof_handler.n_dofs());
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    SparsityPattern        sparsity_pattern;
+    DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
+    sparsity_pattern.copy_from(dsp);
+
+    SparseMatrix<double> system_matrix;
+    system_matrix.reinit(sparsity_pattern);
+
+    const double diffusion_coefficient = 1.0;
+
+    const auto exact_solution = std::make_shared<SmoothSolution<dim>>();
+    const auto rhs_function   = std::make_shared<SmoothRightHandSide<dim>>();
+
+    const auto cell_worker =
+      [&](const auto &cell, auto &scratch_data, auto &copy_data) {
+        const FEValues<dim> &fe_v          = scratch_data.reinit(cell);
+        const unsigned int   dofs_per_cell = fe_v.dofs_per_cell;
+        copy_data.reinit(cell, dofs_per_cell);
+
+        const auto &       q_points    = scratch_data.get_quadrature_points();
+        const unsigned int n_q_points  = q_points.size();
+        const std::vector<double> &JxW = scratch_data.get_JxW_values();
+
+        std::vector<double> rhs(n_q_points);
+        rhs_function->value_list(q_points, rhs);
+
+        for (unsigned int point = 0; point < n_q_points; ++point)
+          for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
+            {
+              for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
+                copy_data.cell_matrix(i, j) +=
+                  diffusion_coefficient *     // nu
+                  fe_v.shape_grad(i, point) * // grad v_h
+                  fe_v.shape_grad(j, point) * // grad u_h
+                  JxW[point];                 // dx
+
+              copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h
+                                       rhs[point] *                 // f
+                                       JxW[point];                  // dx
+            }
+      };
+
+    const auto boundary_worker = [&](const auto &        cell,
+                                     const unsigned int &face_no,
+                                     auto &              scratch_data,
+                                     auto &              copy_data) {
+      const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
+
+      const auto &       q_points      = scratch_data.get_quadrature_points();
+      const unsigned int n_q_points    = q_points.size();
+      const unsigned int dofs_per_cell = fe_fv.dofs_per_cell;
+
+      const std::vector<double> &        JxW = scratch_data.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals =
+        scratch_data.get_normal_vectors();
+
+      std::vector<double> g(n_q_points);
+      exact_solution->value_list(q_points, g);
+
+      const double penalty = PENALTY;
+
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              copy_data.cell_matrix(i, j) +=
+                (-diffusion_coefficient *        // - nu
+                   fe_fv.shape_value(i, point) * // v_h
+                   (fe_fv.shape_grad(j, point) * // (grad u_h .
+                    normals[point])              //  n)
+
+                 - diffusion_coefficient *         // - nu
+                     (fe_fv.shape_grad(i, point) * // (grad v_h .
+                      normals[point]) *            //  n)
+                     fe_fv.shape_value(j, point)   // u_h
+
+                 + diffusion_coefficient * penalty * // + nu sigma
+                     fe_fv.shape_value(i, point) *   // v_h
+                     fe_fv.shape_value(j, point)     // u_h
+
+                 ) *
+                JxW[point]; // dx
+
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            copy_data.cell_rhs(i) +=
+              (-diffusion_coefficient *        // - nu
+                 (fe_fv.shape_grad(i, point) * // (grad v_h .
+                  normals[point]) *            //  n)
+                 g[point]                      // g
+
+
+               + diffusion_coefficient * penalty *        // + nu sigma
+                   fe_fv.shape_value(i, point) * g[point] // v_h g
+
+               ) *
+              JxW[point]; // dx
+        }
+    };
+
+    const auto face_worker = [&](const auto &        cell,
+                                 const unsigned int &f,
+                                 const unsigned int &sf,
+                                 const auto &        ncell,
+                                 const unsigned int &nf,
+                                 const unsigned int &nsf,
+                                 auto &              scratch_data,
+                                 auto &              copy_data) {
+      const FEInterfaceValues<dim> &fe_iv =
+        scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
+
+      const auto &       q_points   = fe_iv.get_quadrature_points();
+      const unsigned int n_q_points = q_points.size();
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace &     copy_data_face = copy_data.face_data.back();
+      const unsigned int n_dofs_face    = fe_iv.n_current_interface_dofs();
+      copy_data_face.joint_dof_indices  = fe_iv.get_interface_dof_indices();
+      copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face);
+
+      const std::vector<double> &        JxW     = fe_iv.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+      const double penalty = PENALTY;
+
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        {
+          for (unsigned int i = 0; i < n_dofs_face; ++i)
+            for (unsigned int j = 0; j < n_dofs_face; ++j)
+              copy_data_face.cell_matrix(i, j) +=
+                (-diffusion_coefficient *              // - nu
+                   fe_iv.jump(i, point) *              // [v_h]
+                   (fe_iv.average_gradient(j, point) * // ({grad u_h} .
+                    normals[point])                    //  n)
+
+                 - diffusion_coefficient *               // - nu
+                     (fe_iv.average_gradient(i, point) * // (grad v_h .
+                      normals[point]) *                  //  n)
+                     fe_iv.jump(j, point)                // [u_h]
+
+                 + diffusion_coefficient * penalty * // + nu sigma
+                     fe_iv.jump(i, point) *          // [v_h]
+                     fe_iv.jump(j, point)            // [u_h]
+
+                 ) *
+                JxW[point]; // dx
+        }
+    };
+
+    AffineConstraints<double> constraints;
+    constraints.close();
+    const auto copier = [&](const auto &c) {
+      constraints.distribute_local_to_global(c.cell_matrix,
+                                             c.cell_rhs,
+                                             c.local_dof_indices,
+                                             system_matrix,
+                                             system_rhs);
+
+      for (auto &cdf : c.face_data)
+        {
+          constraints.distribute_local_to_global(cdf.cell_matrix,
+                                                 cdf.joint_dof_indices,
+                                                 system_matrix);
+        }
+    };
+
+    UpdateFlags cell_flags = update_values | update_gradients |
+                             update_quadrature_points | update_JxW_values;
+    UpdateFlags face_flags = update_values | update_gradients |
+                             update_quadrature_points | update_normal_vectors |
+                             update_JxW_values;
+
+    MeshWorker::ScratchData<dim> scratch_data(
+      mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
+    CopyData cd;
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          cd,
+                          MeshWorker::assemble_own_cells |
+                            MeshWorker::assemble_boundary_faces |
+                            MeshWorker::assemble_own_interior_faces_once,
+                          boundary_worker,
+                          face_worker);
+
+    return solve_and_postprocess(system_matrix, solution, system_rhs);
+  };
+
+  mb_algo();
+  mf_algo();
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  deallog.depth_file(1);
+
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+
+  test<2>(/*degree=*/1);
+  test<2>(/*degree=*/2);
+  test<3>(/*degree=*/1);
+  test<3>(/*degree=*/2);
+}
diff --git a/tests/simplex/matrix_free_03.with_simplex_support=on.output b/tests/simplex/matrix_free_03.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..730828d
--- /dev/null
@@ -0,0 +1,9 @@
+
+DEAL::dim=2 degree=1 0.0412453
+DEAL::dim=2 degree=1 0.0412453
+DEAL::dim=2 degree=2 0.0412685
+DEAL::dim=2 degree=2 0.0412685
+DEAL::dim=3 degree=1 0.0258555
+DEAL::dim=3 degree=1 0.0258555
+DEAL::dim=3 degree=2 0.0250053
+DEAL::dim=3 degree=2 0.0250053

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.