From: Benjamin Brands Date: Sun, 15 Dec 2019 21:01:38 +0000 (+0100) Subject: add test derived from step-4 to check interplay of hp::DoFHandler and parallel::distr... X-Git-Tag: v9.2.0-rc1~778^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F9036%2Fhead;p=dealii.git add test derived from step-4 to check interplay of hp::DoFHandler and parallel::distributed::Triangulation for a locally refined mesh --- diff --git a/tests/mpi/hp_step-4.cc b/tests/mpi/hp_step-4.cc new file mode 100644 index 0000000000..595ca96a56 --- /dev/null +++ b/tests/mpi/hp_step-4.cc @@ -0,0 +1,374 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + +// A lightly adapted version of the step-4 tutorial program. Compared +// to step-4, this test uses hp::DoFHandler and +// parallel::distributed::Triangulation. In the upper half of the domain +// FENothing is used to save effort, as these elements would only be used later. +// The test checks whether the hanging-node constraints are setup correctly. + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +namespace Step4 +{ + using namespace dealii; + + template + class Step4 + { + public: + Step4(); + void + run(); + + private: + void + make_grid(); + void + setup_system(); + void + assemble_system(); + void + solve(); + void + output_results() const; + + MPI_Comm communicator; + + parallel::distributed::Triangulation triangulation; + hp::FECollection fe; + hp::DoFHandler dof_handler; + hp::QCollection q_collection; + + IndexSet locally_owned_dofs; + + IndexSet locally_relevant_dofs; + + TrilinosWrappers::SparseMatrix system_matrix; + + TrilinosWrappers::MPI::Vector solution, rhs; + + AffineConstraints constraints; + + ConditionalOStream pcout; + }; + + + + template + class RightHandSide : public Function + { + public: + RightHandSide() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const override; + }; + + + + template + double + RightHandSide::value(const Point &p, + const unsigned int /*component*/) const + { + double return_value = 0.0; + for (unsigned int i = 0; i < dim; ++i) + return_value += 4.0 * std::pow(p(i), 4.0); + + return return_value; + } + + + + template + Step4::Step4() + : communicator(MPI_COMM_WORLD) + , triangulation(communicator, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)) + , dof_handler(triangulation) + , pcout(std::cout, (Utilities::MPI::this_mpi_process(communicator) == 0)) + { + fe.push_back(FE_Q(2)); + fe.push_back(FE_Nothing()); + + q_collection.push_back(QGauss(fe[0].degree + 1)); + } + + + + template + void + Step4::make_grid() + { + GridGenerator::hyper_cube(triangulation, -1, 1); + + if (dim == 3) + triangulation.refine_global(3); + else + triangulation.refine_global(4); + + for (auto &cell : triangulation.active_cell_iterators()) + if (cell->is_locally_owned()) + if ((cell->center()[dim - 1] < 0.) && (cell->center()[dim - 1] > -0.2)) + cell->set_refine_flag(); + + triangulation.execute_coarsening_and_refinement(); + } + + + + template + void + Step4::setup_system() + { + for (auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + Point center = cell->center(); + + if (center[dim - 1] >= 0) + cell->set_active_fe_index(1); + + else + cell->set_active_fe_index(0); + } + dof_handler.distribute_dofs(fe); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + locally_relevant_dofs.clear(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + locally_relevant_dofs.compress(); + + constraints.clear(); + + constraints.reinit(locally_relevant_dofs); + + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + DoFTools::make_zero_boundary_constraints(dof_handler, constraints); + +#ifdef DEBUG + // We did not think about hp constraints on ghost cells yet. + // Thus, we are content with verifying their consistency for now. + IndexSet locally_active_dofs; + DoFTools::extract_locally_active_dofs(dof_handler, locally_active_dofs); + AssertThrow(constraints.is_consistent_in_parallel( + dof_handler.locally_owned_dofs_per_processor(), + locally_active_dofs, + communicator, + /*verbose=*/true), + ExcMessage( + "AffineConstraints object contains inconsistencies!")); +#endif + + constraints.close(); + + DynamicSparsityPattern dsp(locally_relevant_dofs); + + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + + SparsityTools::distribute_sparsity_pattern( + dsp, + dof_handler.compute_n_locally_owned_dofs_per_processor(), + communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + communicator); + + solution.reinit(locally_relevant_dofs, communicator); + + rhs.reinit(locally_owned_dofs, communicator); + } + + + + template + void + Step4::assemble_system() + { + const RightHandSide right_hand_side; + + system_matrix = 0; + rhs = 0; + + hp::FEValues hp_fe_v(dof_handler.get_fe_collection(), + q_collection, + update_JxW_values | update_gradients | + update_values | update_quadrature_points); + + std::vector local_dof_indices; + + FullMatrix cell_matrix; + Vector cell_rhs; + + for (auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->is_locally_owned() && cell->active_fe_index() == 0) + { + const unsigned int active_fe_index = cell->active_fe_index(); + + if (active_fe_index == 0) + { + hp_fe_v.reinit(cell); + + const FEValues &fe_values = + hp_fe_v.get_present_fe_values(); + + cell_matrix.reinit(fe[active_fe_index].dofs_per_cell, + fe[active_fe_index].dofs_per_cell); + + cell_rhs.reinit(fe[active_fe_index].dofs_per_cell); + + local_dof_indices.resize(fe[active_fe_index].dofs_per_cell); + + cell->get_dof_indices(local_dof_indices); + + const unsigned int dofs_per_cell = + fe[active_fe_index].dofs_per_cell; + + const unsigned int n_q_points = + q_collection[active_fe_index].size(); + + for (unsigned int q_index = 0; q_index < n_q_points; ++q_index) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += (fe_values.shape_grad(i, q_index) * + fe_values.shape_grad(j, q_index) * + fe_values.JxW(q_index)); + + const auto x_q = fe_values.quadrature_point(q_index); + cell_rhs(i) += + (fe_values.shape_value(i, q_index) * + right_hand_side.value(x_q) * fe_values.JxW(q_index)); + } + + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, rhs); + } + } + } + system_matrix.compress(VectorOperation::add); + + rhs.compress(VectorOperation::add); + } + + + + template + void + Step4::solve() + { + SolverControl solver_control(system_matrix.m(), 1e-12 * rhs.l2_norm()); + + TrilinosWrappers::SolverDirect::AdditionalData add_data(false, + "Amesos_Mumps"); + + TrilinosWrappers::SolverDirect direct_solver(solver_control, add_data); + + TrilinosWrappers::MPI::Vector completely_distributed_solution( + locally_owned_dofs, communicator); + + direct_solver.solve(system_matrix, completely_distributed_solution, rhs); + + constraints.distribute(completely_distributed_solution); + + solution = completely_distributed_solution; + } + + + + template + void + Step4::run() + { + pcout << "Solving problem in " << dim << " space dimensions." << std::endl; + + make_grid(); + + setup_system(); + + pcout << " Number of active cells: " + << triangulation.n_global_active_cells() << std::endl; + + pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + assemble_system(); + solve(); + } +} // namespace Step4 + + + +int +main(int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + Step4::Step4<2> laplace_problem_2d; + + laplace_problem_2d.run(); + + Step4::Step4<3> laplace_problem_3d; + + laplace_problem_3d.run(); + + return 0; +} diff --git a/tests/mpi/hp_step-4.with_trilinos=true.mpirun=1.output b/tests/mpi/hp_step-4.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..7940d91d53 --- /dev/null +++ b/tests/mpi/hp_step-4.with_trilinos=true.mpirun=1.output @@ -0,0 +1,6 @@ +Solving problem in 2 space dimensions. + Number of active cells: 352 + Number of degrees of freedom: 997 +Solving problem in 3 space dimensions. + Number of active cells: 960 + Number of degrees of freedom: 7387 diff --git a/tests/mpi/hp_step-4.with_trilinos=true.mpirun=2.output b/tests/mpi/hp_step-4.with_trilinos=true.mpirun=2.output new file mode 100644 index 0000000000..7940d91d53 --- /dev/null +++ b/tests/mpi/hp_step-4.with_trilinos=true.mpirun=2.output @@ -0,0 +1,6 @@ +Solving problem in 2 space dimensions. + Number of active cells: 352 + Number of degrees of freedom: 997 +Solving problem in 3 space dimensions. + Number of active cells: 960 + Number of degrees of freedom: 7387 diff --git a/tests/mpi/hp_step-4.with_trilinos=true.mpirun=4.output b/tests/mpi/hp_step-4.with_trilinos=true.mpirun=4.output new file mode 100644 index 0000000000..7940d91d53 --- /dev/null +++ b/tests/mpi/hp_step-4.with_trilinos=true.mpirun=4.output @@ -0,0 +1,6 @@ +Solving problem in 2 space dimensions. + Number of active cells: 352 + Number of degrees of freedom: 997 +Solving problem in 3 space dimensions. + Number of active cells: 960 + Number of degrees of freedom: 7387