From c66cf7ac57d765bc6a7b6fc114cfe5cadbcef442 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 17 Nov 2020 19:53:38 +0100 Subject: [PATCH] Add Poisson/Helmholtz test: MatrixFree + simplex mesh --- tests/simplex/matrix_free_01.cc | 294 ++++++++++++++++++ ...rix_free_01.with_simplex_support=on.output | 9 + 2 files changed, 303 insertions(+) create mode 100644 tests/simplex/matrix_free_01.cc create mode 100644 tests/simplex/matrix_free_01.with_simplex_support=on.output 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 + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + + +template +class PoissonOperator +{ +public: + using VectorType = LinearAlgebra::distributed::Vector; + + PoissonOperator(const MatrixFree &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( + [&](const auto &, auto &dst, const auto &, const auto cells) { + FEEvaluation 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( + [&](const auto &, auto &dst, const auto &src, const auto cells) { + FEEvaluation 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 &matrix_free; + const bool do_helmholtz; +}; + +template +void +test(const unsigned int degree, const bool do_helmholtz) +{ + Triangulation tria; + + GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 16 : 8); + + Simplex::FE_P fe(degree); + Simplex::QGauss quad(degree + 1); + MappingFE mapping(Simplex::FE_P(1)); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + AffineConstraints constraints; + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 0, Functions::ZeroFunction(), constraints); + constraints.close(); + + const auto solve_and_postprocess = + [&](const auto &poisson_operator, + auto & x, + auto & b) -> std::pair { + ReductionControl reduction_control; + SolverCG::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 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::AdditionalData additional_data; + additional_data.mapping_update_flags = update_gradients | update_values; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, quad, additional_data); + + PoissonOperator poisson_operator(matrix_free, do_helmholtz); + + LinearAlgebra::distributed::Vector 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 x, b; + + x.reinit(dof_handler.n_dofs()); + b.reinit(dof_handler.n_dofs()); + + SparseMatrix 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 fe_values(mapping, fe, quad, flags); + + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector 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(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. -- 2.39.5