From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 19 Jun 2022 07:50:04 +0000 (+0200)
Subject: Add test (MG+DG+hp)
X-Git-Tag: v9.5.0-rc1~1153^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=531222e76a1eeb32e5e4f8ac5b432f167450b78f;p=dealii.git

Add test (MG+DG+hp)
---

diff --git a/tests/matrix_free/poisson_dg_hp.cc b/tests/matrix_free/poisson_dg_hp.cc
new file mode 100644
index 0000000000..7907f02386
--- /dev/null
+++ b/tests/matrix_free/poisson_dg_hp.cc
@@ -0,0 +1,363 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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 problem on a with DG, MatrixFree and hp.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#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_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_wedge_p.h>
+#include <deal.II/fe/mapping_fe.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/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 "../tests.h"
+
+using namespace dealii;
+
+
+double
+get_penalty_parameter(const unsigned int i,
+                      const unsigned int j,
+                      const unsigned int degree)
+{
+  if (degree == 1)
+    {
+      if (i != j)
+        return 32.0;
+      if (i == 0)
+        return 32.0;
+      if (i == 1)
+        return 64.0;
+    }
+  else if (degree == 2)
+    {
+      if (i != j)
+        return 32.0;
+      if (i == 0)
+        return 32.0;
+      if (i == 1)
+        return 64.0;
+    }
+
+  Assert(false, ExcNotImplemented());
+
+  return 0.0;
+}
+
+
+
+template <int dim>
+class PoissonOperator
+{
+public:
+  using VectorType = LinearAlgebra::distributed::Vector<double>;
+  using number     = double;
+
+  using FECellIntegrator = FEEvaluation<dim, -1, 0, 1, number>;
+  using FEFaceIntegrator = FEFaceEvaluation<dim, -1, 0, 1, number>;
+
+  PoissonOperator(const MatrixFree<dim, double> &matrix_free,
+                  const unsigned int             degree)
+    : matrix_free(matrix_free)
+    , degree(degree)
+  {}
+
+  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 range) {
+        FECellIntegrator phi(matrix_free, range);
+
+        for (unsigned int cell = range.first; cell < range.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 loop<VectorType, VectorType>(
+      [&](const auto &data, auto &dst, const auto &src, const auto range) {
+        FECellIntegrator phi(matrix_free, range);
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            phi.reinit(cell);
+            phi.gather_evaluate(src, EvaluationFlags::gradients);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              phi.submit_gradient(phi.get_gradient(q), q);
+            phi.integrate_scatter(EvaluationFlags::gradients, dst);
+          }
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto range) {
+        FEFaceIntegrator fe_eval(data, range, true);
+        FEFaceIntegrator fe_eval_neighbor(data, range, false);
+
+        for (unsigned int face = range.first; 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 =
+              get_penalty_parameter(data.get_face_active_fe_index(range, true),
+                                    data.get_face_active_fe_index(range, false),
+                                    degree);
+
+            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 range) {
+        FEFaceIntegrator fe_eval(data, range, true);
+
+        for (unsigned int face = range.first; face < range.second; ++face)
+          {
+            fe_eval.reinit(face);
+            fe_eval.gather_evaluate(src,
+                                    EvaluationFlags::values |
+                                      EvaluationFlags::gradients);
+            VectorizedArray<number> sigmaF =
+              get_penalty_parameter(data.get_face_active_fe_index(range),
+                                    data.get_face_active_fe_index(range),
+                                    degree);
+
+            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,
+      true);
+  }
+
+private:
+  const MatrixFree<dim, double> &matrix_free;
+  const unsigned int             degree;
+};
+
+template <int dim>
+void
+test(const unsigned int degree)
+{
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+  unsigned int subdivisions = degree == 1 ? 16 : 8;
+
+  GridGenerator::subdivided_hyper_cube(tria, subdivisions);
+
+  FE_DGQ<dim>           fe1(degree);
+  FE_DGQ<dim>           fe2(degree + 1);
+  hp::FECollection<dim> fes(fe1, fe2);
+
+  QGauss<dim>   quad(degree + 2);
+  MappingQ<dim> mapping(1);
+
+  DoFHandler<dim> dof_handler(tria);
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    if (cell->is_locally_owned())
+      {
+        if (cell->center()[0] < 0.5)
+          cell->set_active_fe_index(0);
+        else
+          cell->set_active_fe_index(1);
+      }
+
+  dof_handler.distribute_dofs(fes);
+
+  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(2000, 1e-7, 1e-2);
+    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 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, degree + 1);
+    data_out.write_vtu_with_pvtu_record("./",
+                                        "result-" + std::to_string(dim) + "-" +
+                                          std::to_string(degree),
+                                        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,
+                                      quad,
+                                      VectorTools::NormType::L2_norm);
+
+    const double error =
+      VectorTools::compute_global_error(tria,
+                                        difference,
+                                        VectorTools::NormType::L2_norm);
+
+    if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+      printf("Error %f.\n", error);
+
+    if (error < 0.042)
+      deallog << "OK" << std::endl;
+    else
+      deallog << "FAIL" << 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, quad, additional_data);
+
+    PoissonOperator<dim> poisson_operator(matrix_free, degree);
+
+    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);
+  };
+
+  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);
+}
diff --git a/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=1.output b/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=1.output
new file mode 100644
index 0000000000..9bb3640efd
--- /dev/null
+++ b/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=1.output
@@ -0,0 +1,3 @@
+
+DEAL::dim=2 degree=1 OK
+DEAL::dim=2 degree=2 OK
diff --git a/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=2.output b/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=2.output
new file mode 100644
index 0000000000..9bb3640efd
--- /dev/null
+++ b/tests/matrix_free/poisson_dg_hp.with_p4est=true.mpirun=2.output
@@ -0,0 +1,3 @@
+
+DEAL::dim=2 degree=1 OK
+DEAL::dim=2 degree=2 OK