From 56f92cadfce50bdcda9dae1119f954e6e197b373 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 22 May 2021 14:08:50 +0200 Subject: [PATCH] Add test [WIP] --- .../multigrid_p_02.cc | 165 ++++++++++++++++++ ...id_p_02.mpirun=1.with_trilinos=true.output | 10 ++ 2 files changed, 175 insertions(+) create mode 100644 tests/multigrid-global-coarsening/multigrid_p_02.cc create mode 100644 tests/multigrid-global-coarsening/multigrid_p_02.mpirun=1.with_trilinos=true.output diff --git a/tests/multigrid-global-coarsening/multigrid_p_02.cc b/tests/multigrid-global-coarsening/multigrid_p_02.cc new file mode 100644 index 0000000000..3bdfca8875 --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_p_02.cc @@ -0,0 +1,165 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2021 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. +// +// --------------------------------------------------------------------- + + +/** + * Test p-multigrid for a uniformly and locally refined meshes both for + * simplex and hypercube mesh. + */ + +#include "multigrid_util.h" + +template +void +test(const unsigned int n_refinements, const unsigned int fe_degree_fine) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + Triangulation tria( + Triangulation::MeshSmoothing::limit_level_difference_at_vertices); + GridGenerator::subdivided_hyper_cube( + tria, 2 * Utilities::pow(2, n_refinements)); + + const auto level_degrees = + MGTransferGlobalCoarseningTools::create_polynomial_coarsening_sequence( + fe_degree_fine, + MGTransferGlobalCoarseningTools::PolynomialCoarseningSequenceType:: + bisect); + + const unsigned int min_level = 0; + const unsigned int max_level = level_degrees.size() - 1; + + MGLevelObject> dof_handlers(min_level, max_level, tria); + MGLevelObject> constraints(min_level, max_level); + MGLevelObject> transfers(min_level, + max_level); + MGLevelObject> operators(min_level, max_level); + + std::unique_ptr> mapping_; + + // set up levels + for (auto l = min_level; l <= max_level; ++l) + { + auto &dof_handler = dof_handlers[l]; + auto &constraint = constraints[l]; + auto &op = operators[l]; + + std::unique_ptr> fe; + std::unique_ptr> quad; + std::unique_ptr> mapping; + + fe = std::make_unique>(level_degrees[l]); + quad = std::make_unique>(level_degrees[l] + 1); + mapping = std::make_unique>(FE_Q(1)); + + if (l == max_level) + mapping_ = mapping->clone(); + + // set up dofhandler + dof_handler.distribute_dofs(*fe); + dof_handler.distribute_mg_dofs(); + + // set up constraints + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(dof_handler, + locally_relevant_dofs); + constraint.reinit(locally_relevant_dofs); + VectorTools::interpolate_boundary_values( + *mapping, dof_handler, 0, Functions::ZeroFunction(), constraint); + DoFTools::make_hanging_node_constraints(dof_handler, constraint); + constraint.close(); + + // set up operator + op.reinit(*mapping, dof_handler, *quad, constraint); + } + + // set up transfer operator + for (unsigned int l = min_level; l < max_level; ++l) + transfers[l + 1].reinit_polynomial_transfer(dof_handlers[l + 1], + dof_handlers[l], + constraints[l + 1], + constraints[l], + 0, + 0); + + MGTransferGlobalCoarsening transfer( + transfers, + [&](const auto l, auto &vec) { operators[l].initialize_dof_vector(vec); }); + + GMGParameters mg_data; // TODO + + VectorType dst, src; + operators[max_level].initialize_dof_vector(dst); + operators[max_level].initialize_dof_vector(src); + + operators[max_level].rhs(src); + + ReductionControl solver_control( + mg_data.maxiter, mg_data.abstol, mg_data.reltol, false, false); + + mg_solve(solver_control, + dst, + src, + mg_data, + dof_handlers[max_level], + operators[max_level], + operators, + transfer); + + constraints[max_level].distribute(dst); + + deallog << dim << " " << fe_degree_fine << " " << n_refinements << " " + << solver_control.last_step() << std::endl; + + return; + + static unsigned int counter = 0; + + MGLevelObject results(min_level, max_level); + + transfer.interpolate_to_mg(dof_handlers[max_level], results, dst); + + for (unsigned int l = min_level; l <= max_level; ++l) + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handlers[l]); + data_out.add_data_vector( + results[l], + "solution", + DataOut_DoFData, dim>::DataVectorType::type_dof_data); + data_out.build_patches(*mapping_, 2); + + std::ofstream output("test." + std::to_string(dim) + "." + + std::to_string(counter) + "." + std::to_string(l) + + ".vtk"); + data_out.write_vtk(output); + } + + counter++; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + deallog.precision(8); + + for (unsigned int n_refinements = 2; n_refinements <= 4; ++n_refinements) + for (unsigned int degree = 2; degree <= 4; ++degree) + test<2>(n_refinements, degree); +} diff --git a/tests/multigrid-global-coarsening/multigrid_p_02.mpirun=1.with_trilinos=true.output b/tests/multigrid-global-coarsening/multigrid_p_02.mpirun=1.with_trilinos=true.output new file mode 100644 index 0000000000..24dc959efe --- /dev/null +++ b/tests/multigrid-global-coarsening/multigrid_p_02.mpirun=1.with_trilinos=true.output @@ -0,0 +1,10 @@ + +DEAL:0::2 2 2 4 +DEAL:0::2 3 2 6 +DEAL:0::2 4 2 8 +DEAL:0::2 2 3 6 +DEAL:0::2 3 3 11 +DEAL:0::2 4 3 17 +DEAL:0::2 2 4 13 +DEAL:0::2 3 4 22 +DEAL:0::2 4 4 32 -- 2.39.5