From 672f97d263f2279a67a0518383ecd3d22d31dbe4 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 21 Nov 2020 14:45:05 +0100 Subject: [PATCH] Add Poisson/Helmholtz test: MatrixFree + mixed mesh --- tests/simplex/matrix_free_02.cc | 476 ++++++++++++++++++ ...rix_free_02.with_simplex_support=on.output | 17 + 2 files changed, 493 insertions(+) create mode 100644 tests/simplex/matrix_free_02.cc create mode 100644 tests/simplex/matrix_free_02.with_simplex_support=on.output diff --git a/tests/simplex/matrix_free_02.cc b/tests/simplex/matrix_free_02.cc new file mode 100644 index 0000000000..84f58b9d01 --- /dev/null +++ b/tests/simplex/matrix_free_02.cc @@ -0,0 +1,476 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// Same as matrix_free_01 but testing mixed meshes (and also pure simplex and +// hypercube mesh as special case of mixed meshs). + +#include + +#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; + +namespace dealii +{ + namespace GridGenerator + { + template + void + subdivided_hyper_rectangle_with_simplices_mix( + Triangulation & tria, + const std::vector &repetitions, + const Point & p1, + const Point & p2, + const bool colorize = false) + { + AssertDimension(dim, spacedim); + + AssertThrow(colorize == false, ExcNotImplemented()); + + std::vector> vertices; + std::vector> cells; + + if (dim == 2) + { + // determine cell sizes + const Point dx((p2[0] - p1[0]) / repetitions[0], + (p2[1] - p1[1]) / repetitions[1]); + + // create vertices + for (unsigned int j = 0; j <= repetitions[1]; ++j) + for (unsigned int i = 0; i <= repetitions[0]; ++i) + vertices.push_back( + Point(p1[0] + dx[0] * i, p1[1] + dx[1] * j)); + + // create cells + for (unsigned int j = 0; j < repetitions[1]; ++j) + for (unsigned int i = 0; i < repetitions[0]; ++i) + { + // create reference QUAD cell + std::array quad{{ + (j + 0) * (repetitions[0] + 1) + i + 0, // + (j + 0) * (repetitions[0] + 1) + i + 1, // + (j + 1) * (repetitions[0] + 1) + i + 0, // + (j + 1) * (repetitions[0] + 1) + i + 1 // + }}; // + + if (j < repetitions[1] / 2 && i < repetitions[0] / 2) + { + CellData quad_; + quad_.vertices = {quad[0], quad[1], quad[2], quad[3]}; + cells.push_back(quad_); + + continue; + } + + // TRI cell 0 + { + CellData tri; + tri.vertices = {quad[0], quad[1], quad[2]}; + cells.push_back(tri); + } + + // TRI cell 1 + { + CellData tri; + tri.vertices = {quad[3], quad[2], quad[1]}; + cells.push_back(tri); + } + } + } + else + { + AssertThrow(colorize == false, ExcNotImplemented()); + } + + // actually create triangulation + tria.create_triangulation(vertices, cells, SubCellData()); + } + + + template + void + subdivided_hyper_cube_with_simplices_mix(Triangulation &tria, + const unsigned int repetitions, + const double p1 = 0.0, + const double p2 = 1.0, + const bool colorize = false) + { + if (dim == 2) + { + subdivided_hyper_rectangle_with_simplices_mix( + tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize); + } + else if (dim == 3) + { + subdivided_hyper_rectangle_with_simplices_mix( + tria, + {{repetitions, repetitions, repetitions}}, + {p1, p1, p1}, + {p2, p2, p2}, + colorize); + } + else + { + AssertThrow(false, ExcNotImplemented()) + } + } + } // namespace GridGenerator +} // 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 &data, auto &dst, const auto &, const auto cell_range) { + for (unsigned int i = 0; i < 2; ++i) + { + const auto cell_subrange = + data.create_cell_subrange_hp_by_index(cell_range, i); + + FEEvaluation phi(matrix_free, 0, 0, 0, i, i); + + for (unsigned int cell = cell_subrange.first; + cell < cell_subrange.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 &data, auto &dst, const auto &src, const auto cell_range) { + for (unsigned int i = 0; i < 2; ++i) + { + const auto cell_subrange = + data.create_cell_subrange_hp_by_index(cell_range, i); + + FEEvaluation phi(matrix_free, 0, 0, 0, i, i); + for (unsigned int cell = cell_subrange.first; + cell < cell_subrange.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 version, const unsigned int degree, const bool do_helmholtz) +{ + Triangulation tria; + + const unsigned int subdivisions = dim == 2 ? 25 : 8; + + if (version == 0) + GridGenerator::subdivided_hyper_cube_with_simplices(tria, subdivisions); + else if (version == 1) + GridGenerator::subdivided_hyper_cube(tria, subdivisions); + else if (version == 2) + GridGenerator::subdivided_hyper_cube_with_simplices_mix(tria, subdivisions); + + Simplex::FE_P fe1(degree); + FE_Q fe2(degree); + hp::FECollection fes(fe1, fe2); + + Simplex::QGauss quad1(degree + 1); + QGauss quad2(degree + 1); + hp::QCollection quads(quad1, quad2); + + MappingFE mapping1(Simplex::FE_P(1)); + MappingQ mapping2(1); + hp::MappingCollection mappings(mapping1, mapping2); + + DoFHandler dof_handler(tria); + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->reference_cell_type() == ReferenceCell::Type::Tri || + cell->reference_cell_type() == ReferenceCell::Type::Tet) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); + + dof_handler.distribute_dofs(fes); + + AffineConstraints constraints; + DoFTools::make_zero_boundary_constraints(dof_handler, 0, constraints); + constraints.close(); + + const auto solve_and_postprocess = + [&](const auto &poisson_operator, + auto & x, + auto & b) -> std::tuple { + ReductionControl reduction_control(1000, 1e-10, 1e-4); + 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(mappings, 2); + data_out.write_vtu_with_pvtu_record("./", "result", 0, MPI_COMM_WORLD); +#endif + + Vector difference(tria.n_active_cells()); + + VectorTools::integrate_difference(mappings, + dof_handler, + x, + Functions::ZeroFunction(), + difference, + quads, + VectorTools::NormType::L2_norm); + + return {reduction_control.last_step(), + reduction_control.last_value(), + x.linfty_norm(), + VectorTools::compute_global_error(tria, + difference, + VectorTools::NormType::L2_norm)}; + }; + + const auto mf_algo = [&]() { + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags = update_gradients | update_values; + + MatrixFree matrix_free; + matrix_free.reinit( + mappings, dof_handler, constraints, quads, 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; + + hp::FEValues hp_fe_values(mappings, fes, quads, 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; + + hp_fe_values.reinit(cell); + + auto &fe_values = hp_fe_values.get_present_fe_values(); + + const unsigned int dofs_per_cell = cell->get_fe().n_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().n_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(std::get<0>(result_mf), std::get<0>(result_mb)); + Assert(std::abs(std::get<1>(result_mf) - std::get<1>(result_mb)) < 1e-6, + ExcNotImplemented()); + Assert(std::abs(std::get<2>(result_mf) - std::get<2>(result_mb)) < 1e-6, + ExcNotImplemented()); + Assert(std::abs(std::get<3>(result_mf) - std::get<3>(result_mb)) < 1e-6, + ExcNotImplemented()); + + deallog << "mesh="; + if (version == 0) + deallog << "P"; + else if (version == 1) + deallog << "Q"; + else if (version == 2) + deallog << "M"; + deallog << " : "; + + deallog << "dim=" << dim << " "; + deallog << "degree=" << degree << " "; + deallog << "Type="; + + if (do_helmholtz) + deallog << "Helmholtz"; + else + deallog << "Possion "; + deallog << " : "; + + deallog << "Convergence step " << std::get<0>(result_mf) << " value " + << std::get<1>(result_mf) << " max " << std::get<2>(result_mf) + << " norm " << std::get<3>(result_mf) << "." << 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); + + for (unsigned int i = 0; i < 3; ++i) + test<2>(i, /*degree=*/1, /*do_helmholtz*/ false); + deallog << std::endl; + + for (unsigned int i = 0; i < 3; ++i) + test<2>(i, /*degree=*/1, /*do_helmholtz*/ true); + deallog << std::endl; + + for (unsigned int i = 0; i < 3; ++i) + test<2>(i, /*degree=*/2, /*do_helmholtz*/ false); + deallog << std::endl; + + for (unsigned int i = 0; i < 3; ++i) + test<2>(i, /*degree=*/2, /*do_helmholtz*/ true); + deallog << std::endl; +} diff --git a/tests/simplex/matrix_free_02.with_simplex_support=on.output b/tests/simplex/matrix_free_02.with_simplex_support=on.output new file mode 100644 index 0000000000..d13a782577 --- /dev/null +++ b/tests/simplex/matrix_free_02.with_simplex_support=on.output @@ -0,0 +1,17 @@ + +DEAL::mesh=P : dim=2 degree=1 Type=Possion : Convergence step 32 value 1.97679e-06 max 0.0733790 norm 0.0410935. +DEAL::mesh=Q : dim=2 degree=1 Type=Possion : Convergence step 23 value 2.31438e-06 max 0.0735643 norm 0.0412027. +DEAL::mesh=M : dim=2 degree=1 Type=Possion : Convergence step 37 value 2.47295e-06 max 0.0735235 norm 0.0411386. +DEAL:: +DEAL::mesh=P : dim=2 degree=1 Type=Helmholtz : Convergence step 32 value 1.84146e-06 max 0.0695486 norm 0.0391345. +DEAL::mesh=Q : dim=2 degree=1 Type=Helmholtz : Convergence step 23 value 2.11901e-06 max 0.0697145 norm 0.0392337. +DEAL::mesh=M : dim=2 degree=1 Type=Helmholtz : Convergence step 33 value 3.66201e-06 max 0.0696810 norm 0.0391757. +DEAL:: +DEAL::mesh=P : dim=2 degree=2 Type=Possion : Convergence step 76 value 1.83468e-06 max 0.0736712 norm 0.0412614. +DEAL::mesh=Q : dim=2 degree=2 Type=Possion : Convergence step 64 value 2.13749e-06 max 0.0736712 norm 0.0412615. +DEAL::mesh=M : dim=2 degree=2 Type=Possion : Convergence step 77 value 1.77910e-06 max 0.0736712 norm 0.0412614. +DEAL:: +DEAL::mesh=P : dim=2 degree=2 Type=Helmholtz : Convergence step 76 value 1.76851e-06 max 0.0698084 norm 0.0392874. +DEAL::mesh=Q : dim=2 degree=2 Type=Helmholtz : Convergence step 64 value 1.95880e-06 max 0.0698084 norm 0.0392875. +DEAL::mesh=M : dim=2 degree=2 Type=Helmholtz : Convergence step 76 value 2.21571e-06 max 0.0698084 norm 0.0392875. +DEAL:: -- 2.39.5