From: marcfehling Date: Fri, 5 Jul 2019 13:44:59 +0000 (+0200) Subject: New tests: Parallel distributed version of step-27. X-Git-Tag: v9.2.0-rc1~1345^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F8421%2Fhead;p=dealii.git New tests: Parallel distributed version of step-27. --- diff --git a/tests/mpi/petsc_step-27.cc b/tests/mpi/petsc_step-27.cc new file mode 100644 index 0000000000..c36531d968 --- /dev/null +++ b/tests/mpi/petsc_step-27.cc @@ -0,0 +1,561 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2006 - 2018 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. + * + * --------------------------------------------------------------------- + */ + + + +// parallelized version of step-27 with PETSc + + +#include +namespace LA +{ + using namespace dealii::LinearAlgebraPETSc; +} // namespace LA + +#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 +#include + +#include +#include +#include + +#include "../tests.h" + + +namespace Step27 +{ + using namespace dealii; + + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + ~LaplaceProblem(); + + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + create_coarse_grid(); + void + estimate_smoothness(Vector &smoothness_indicators); + void + postprocess(); + std::pair + predicate(const TableIndices &indices); + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + hp::QCollection face_quadrature_collection; + + hp::QCollection fourier_q_collection; + std::shared_ptr> fourier; + std::vector ln_k; + Table> fourier_coefficients; + + AffineConstraints constraints; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + LA::MPI::SparseMatrix system_matrix; + + LA::MPI::Vector solution; + LA::MPI::Vector system_rhs; + + const unsigned int max_degree; + + ConditionalOStream pcout; + }; + + + + template + class RightHandSide : public Function + { + public: + RightHandSide() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component) const override; + }; + + + template + double + RightHandSide::value(const Point &p, + const unsigned int /*component*/) const + { + double product = 1; + for (unsigned int d = 0; d < dim; ++d) + product *= (p[d] + 1); + return product; + } + + + + template + void + resize(Table &coeff, const unsigned int N) + { + TableIndices size; + for (unsigned int d = 0; d < dim; d++) + size[d] = N; + coeff.reinit(size); + } + + + + template + LaplaceProblem::LaplaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator) + , dof_handler(triangulation) + , max_degree(dim <= 2 ? 7 : 5) + , pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + { + for (unsigned int degree = 2; degree <= max_degree; ++degree) + { + fe_collection.push_back(FE_Q(degree)); + quadrature_collection.push_back(QGauss(degree + 1)); + face_quadrature_collection.push_back(QGauss(degree + 1)); + } + + const unsigned int N = max_degree; + + QGauss<1> base_quadrature(2); + QIterated quadrature(base_quadrature, N); + for (unsigned int i = 0; i < fe_collection.size(); i++) + fourier_q_collection.push_back(quadrature); + + fourier = std::make_shared>(N, + fe_collection, + fourier_q_collection); + + resize(fourier_coefficients, N); + } + + + + template + LaplaceProblem::~LaplaceProblem() + { + dof_handler.clear(); + } + + + + template + void + LaplaceProblem::setup_system() + { + dof_handler.distribute_dofs(fe_collection); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + 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(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + + + template + void + LaplaceProblem::assemble_system() + { + hp::FEValues hp_fe_values(fe_collection, + quadrature_collection, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + + const RightHandSide rhs_function; + + 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()) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_matrix = 0; + + cell_rhs.reinit(dofs_per_cell); + cell_rhs = 0; + + hp_fe_values.reinit(cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + std::vector rhs_values(fe_values.n_quadrature_points); + rhs_function.value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q_point = 0; + q_point < fe_values.n_quadrature_points; + ++q_point) + 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_point) * // grad phi_i(x_q) + fe_values.shape_grad(j, q_point) * // grad phi_j(x_q) + fe_values.JxW(q_point)); // dx + + cell_rhs(i) += + (fe_values.shape_value(i, q_point) * // phi_i(x_q) + rhs_values[q_point] * // f(x_q) + fe_values.JxW(q_point)); // dx + } + + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + + + template + void + LaplaceProblem::solve() + { + LA::MPI::Vector completely_distributed_solution(locally_owned_dofs, + mpi_communicator); + + SolverControl solver_control(system_rhs.size(), + 1e-8 * system_rhs.l2_norm()); + // ^~~~ + // Loss of precision with a factor of 1e-12 with Trilinos + LA::SolverCG cg(solver_control, mpi_communicator); + + LA::MPI::PreconditionAMG preconditioner; + LA::MPI::PreconditionAMG::AdditionalData data; + data.symmetric_operator = true; + preconditioner.initialize(system_matrix, data); + + check_solver_within_range(cg.solve(system_matrix, + completely_distributed_solution, + system_rhs, + preconditioner), + solver_control.last_step(), + 5, + 40); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + solution = completely_distributed_solution; + } + + + + template + void + LaplaceProblem::postprocess() + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator::estimate( + dof_handler, + face_quadrature_collection, + std::map *>(), + solution, + estimated_error_per_cell); + + Vector smoothness_indicators(triangulation.n_active_cells()); + estimate_smoothness(smoothness_indicators); + + parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( + triangulation, estimated_error_per_cell, 0.3, 0.03); + + hp::Refinement::p_adaptivity_from_threshold(dof_handler, + smoothness_indicators, + 0.5, + 0.); + hp::Refinement::choose_p_over_h(dof_handler); + + triangulation.execute_coarsening_and_refinement(); + } + + + + template <> + void + LaplaceProblem<2>::create_coarse_grid() + { + const unsigned int dim = 2; + + const std::vector> vertices = { + {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0}, // + {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5}, // + {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0}, // + {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5}, // + {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}}; + + const std::vector::vertices_per_cell>> + cell_vertices = {{{0, 1, 5, 6}}, + {{1, 2, 6, 7}}, + {{2, 3, 7, 8}}, + {{3, 4, 8, 9}}, + {{5, 6, 10, 11}}, + {{8, 9, 12, 13}}, + {{10, 11, 14, 15}}, + {{12, 13, 17, 18}}, + {{14, 15, 19, 20}}, + {{15, 16, 20, 21}}, + {{16, 17, 21, 22}}, + {{17, 18, 22, 23}}}; + + const unsigned int n_cells = cell_vertices.size(); + + std::vector> cells(n_cells, CellData()); + for (unsigned int i = 0; i < n_cells; ++i) + { + for (unsigned int j = 0; j < GeometryInfo::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + triangulation.create_triangulation(vertices, cells, SubCellData()); + triangulation.refine_global(3); + } + + + + template + void + LaplaceProblem::run() + { + for (unsigned int cycle = 0; cycle < 5; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + create_coarse_grid(); + + setup_system(); + + pcout << " Number of active cells : " + << triangulation.n_global_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << " Number of constraints : " + << Utilities::MPI::sum(constraints.n_constraints(), + mpi_communicator) + << std::endl; + + assemble_system(); + solve(); + postprocess(); + } + } + + + + template + std::pair + LaplaceProblem::predicate(const TableIndices &ind) + { + unsigned int v = 0; + for (unsigned int i = 0; i < dim; i++) + v += ind[i] * ind[i]; + if (v > 0 && v < max_degree * max_degree) + return std::make_pair(true, v); + else + return std::make_pair(false, v); + } + + + + template + void + LaplaceProblem::estimate_smoothness(Vector &smoothness_indicators) + { + Vector local_dof_values; + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + local_dof_values.reinit(cell->get_fe().dofs_per_cell); + cell->get_dof_values(solution, local_dof_values); + + fourier->calculate(local_dof_values, + cell->active_fe_index(), + fourier_coefficients); + + std::pair, std::vector> res = + FESeries::process_coefficients( + fourier_coefficients, + std::bind(&LaplaceProblem::predicate, + this, + std::placeholders::_1), + VectorTools::Linfty_norm); + + Assert(res.first.size() == res.second.size(), ExcInternalError()); + + if (ln_k.size() == 0) + { + ln_k.resize(res.first.size(), 0); + for (unsigned int f = 0; f < ln_k.size(); f++) + ln_k[f] = + std::log(2.0 * numbers::PI * std::sqrt(1. * res.first[f])); + } + + for (double &residual_element : res.second) + residual_element = std::log(residual_element); + + std::pair fit = + FESeries::linear_regression(ln_k, res.second); + + smoothness_indicators(cell->active_cell_index()) = + -fit.first - 1. * dim / 2; + } + } +} // namespace Step27 + + + +int +main(int argc, char *argv[]) +{ + try + { + using namespace dealii; + using namespace Step27; + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + LaplaceProblem<2> laplace_problem; + laplace_problem.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/mpi/petsc_step-27.with_p4est=true.with_petsc=true.mpirun=1.output b/tests/mpi/petsc_step-27.with_p4est=true.with_petsc=true.mpirun=1.output new file mode 100644 index 0000000000..b94c448be7 --- /dev/null +++ b/tests/mpi/petsc_step-27.with_p4est=true.with_petsc=true.mpirun=1.output @@ -0,0 +1,25 @@ +Cycle 0: + Number of active cells : 768 + Number of degrees of freedom: 3264 + Number of constraints : 384 + Solved in 8 iterations. +Cycle 1: + Number of active cells : 966 + Number of degrees of freedom: 5229 + Number of constraints : 928 + Solved in 9 iterations. +Cycle 2: + Number of active cells : 1146 + Number of degrees of freedom: 8527 + Number of constraints : 1952 + Solved in 10 iterations. +Cycle 3: + Number of active cells : 1359 + Number of degrees of freedom: 12486 + Number of constraints : 3120 + Solved in 17 iterations. +Cycle 4: + Number of active cells : 1656 + Number of degrees of freedom: 18434 + Number of constraints : 4819 + Solved in 29 iterations. diff --git a/tests/mpi/petsc_step-27.with_p4est=true.with_petsc=true.mpirun=2.output b/tests/mpi/petsc_step-27.with_p4est=true.with_petsc=true.mpirun=2.output new file mode 100644 index 0000000000..14e48d0906 --- /dev/null +++ b/tests/mpi/petsc_step-27.with_p4est=true.with_petsc=true.mpirun=2.output @@ -0,0 +1,25 @@ +Cycle 0: + Number of active cells : 768 + Number of degrees of freedom: 3264 + Number of constraints : 404 + Solved in 8 iterations. +Cycle 1: + Number of active cells : 966 + Number of degrees of freedom: 5245 + Number of constraints : 956 + Solved in 9 iterations. +Cycle 2: + Number of active cells : 1146 + Number of degrees of freedom: 8549 + Number of constraints : 2005 + Solved in 11 iterations. +Cycle 3: + Number of active cells : 1359 + Number of degrees of freedom: 12537 + Number of constraints : 3256 + Solved in 17 iterations. +Cycle 4: + Number of active cells : 1647 + Number of degrees of freedom: 18544 + Number of constraints : 4954 + Solved in 29 iterations. diff --git a/tests/mpi/trilinos_step-27.cc b/tests/mpi/trilinos_step-27.cc new file mode 100644 index 0000000000..cca0b94db8 --- /dev/null +++ b/tests/mpi/trilinos_step-27.cc @@ -0,0 +1,560 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2006 - 2018 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. + * + * --------------------------------------------------------------------- + */ + + + +// parallelized version of step-27 with Trilinos + + +#include +namespace LA +{ + using namespace dealii::LinearAlgebraTrilinos; +} // namespace LA + +#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 +#include + +#include +#include +#include + +#include "../tests.h" + + +namespace Step27 +{ + using namespace dealii; + + + template + class LaplaceProblem + { + public: + LaplaceProblem(); + ~LaplaceProblem(); + + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + create_coarse_grid(); + void + estimate_smoothness(Vector &smoothness_indicators); + void + postprocess(); + std::pair + predicate(const TableIndices &indices); + + MPI_Comm mpi_communicator; + + parallel::distributed::Triangulation triangulation; + + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + hp::QCollection face_quadrature_collection; + + hp::QCollection fourier_q_collection; + std::shared_ptr> fourier; + std::vector ln_k; + Table> fourier_coefficients; + + AffineConstraints constraints; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + LA::MPI::SparseMatrix system_matrix; + + LA::MPI::Vector solution; + LA::MPI::Vector system_rhs; + + const unsigned int max_degree; + + ConditionalOStream pcout; + }; + + + + template + class RightHandSide : public Function + { + public: + RightHandSide() + : Function() + {} + + virtual double + value(const Point &p, const unsigned int component) const override; + }; + + + template + double + RightHandSide::value(const Point &p, + const unsigned int /*component*/) const + { + double product = 1; + for (unsigned int d = 0; d < dim; ++d) + product *= (p[d] + 1); + return product; + } + + + + template + void + resize(Table &coeff, const unsigned int N) + { + TableIndices size; + for (unsigned int d = 0; d < dim; d++) + size[d] = N; + coeff.reinit(size); + } + + + + template + LaplaceProblem::LaplaceProblem() + : mpi_communicator(MPI_COMM_WORLD) + , triangulation(mpi_communicator) + , dof_handler(triangulation) + , max_degree(dim <= 2 ? 7 : 5) + , pcout(std::cout, + (Utilities::MPI::this_mpi_process(mpi_communicator) == 0)) + { + for (unsigned int degree = 2; degree <= max_degree; ++degree) + { + fe_collection.push_back(FE_Q(degree)); + quadrature_collection.push_back(QGauss(degree + 1)); + face_quadrature_collection.push_back(QGauss(degree + 1)); + } + + const unsigned int N = max_degree; + + QGauss<1> base_quadrature(2); + QIterated quadrature(base_quadrature, N); + for (unsigned int i = 0; i < fe_collection.size(); i++) + fourier_q_collection.push_back(quadrature); + + fourier = std::make_shared>(N, + fe_collection, + fourier_q_collection); + + resize(fourier_coefficients, N); + } + + + + template + LaplaceProblem::~LaplaceProblem() + { + dof_handler.clear(); + } + + + + template + void + LaplaceProblem::setup_system() + { + dof_handler.distribute_dofs(fe_collection); + + locally_owned_dofs = dof_handler.locally_owned_dofs(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + solution.reinit(locally_owned_dofs, + locally_relevant_dofs, + mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + 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(), + mpi_communicator, + locally_relevant_dofs); + + system_matrix.reinit(locally_owned_dofs, + locally_owned_dofs, + dsp, + mpi_communicator); + } + + + + template + void + LaplaceProblem::assemble_system() + { + hp::FEValues hp_fe_values(fe_collection, + quadrature_collection, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + + const RightHandSide rhs_function; + + 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()) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_matrix = 0; + + cell_rhs.reinit(dofs_per_cell); + cell_rhs = 0; + + hp_fe_values.reinit(cell); + + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + std::vector rhs_values(fe_values.n_quadrature_points); + rhs_function.value_list(fe_values.get_quadrature_points(), + rhs_values); + + for (unsigned int q_point = 0; + q_point < fe_values.n_quadrature_points; + ++q_point) + 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_point) * // grad phi_i(x_q) + fe_values.shape_grad(j, q_point) * // grad phi_j(x_q) + fe_values.JxW(q_point)); // dx + + cell_rhs(i) += + (fe_values.shape_value(i, q_point) * // phi_i(x_q) + rhs_values[q_point] * // f(x_q) + fe_values.JxW(q_point)); // dx + } + + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + } + + + + template + void + LaplaceProblem::solve() + { + LA::MPI::Vector completely_distributed_solution(locally_owned_dofs, + mpi_communicator); + + SolverControl solver_control(system_rhs.size(), + 1e-8 * system_rhs.l2_norm()); + // ^~~~ + // Loss of precision with a factor of 1e-12 with Trilinos + LA::SolverCG cg(solver_control); + + LA::MPI::PreconditionAMG preconditioner; + preconditioner.initialize(system_matrix); + + check_solver_within_range(cg.solve(system_matrix, + completely_distributed_solution, + system_rhs, + preconditioner), + solver_control.last_step(), + 10, + 80); + + pcout << " Solved in " << solver_control.last_step() << " iterations." + << std::endl; + + constraints.distribute(completely_distributed_solution); + + solution = completely_distributed_solution; + } + + + + template + void + LaplaceProblem::postprocess() + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator::estimate( + dof_handler, + face_quadrature_collection, + std::map *>(), + solution, + estimated_error_per_cell); + + + Vector smoothness_indicators(triangulation.n_active_cells()); + estimate_smoothness(smoothness_indicators); + + parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number( + triangulation, estimated_error_per_cell, 0.3, 0.03); + + hp::Refinement::p_adaptivity_from_threshold(dof_handler, + smoothness_indicators, + 0.5, + 0.); + hp::Refinement::choose_p_over_h(dof_handler); + + triangulation.execute_coarsening_and_refinement(); + } + + + + template <> + void + LaplaceProblem<2>::create_coarse_grid() + { + const unsigned int dim = 2; + + const std::vector> vertices = { + {-1.0, -1.0}, {-0.5, -1.0}, {+0.0, -1.0}, {+0.5, -1.0}, {+1.0, -1.0}, // + {-1.0, -0.5}, {-0.5, -0.5}, {+0.0, -0.5}, {+0.5, -0.5}, {+1.0, -0.5}, // + {-1.0, +0.0}, {-0.5, +0.0}, {+0.5, +0.0}, {+1.0, +0.0}, // + {-1.0, +0.5}, {-0.5, +0.5}, {+0.0, +0.5}, {+0.5, +0.5}, {+1.0, +0.5}, // + {-1.0, +1.0}, {-0.5, +1.0}, {+0.0, +1.0}, {+0.5, +1.0}, {+1.0, +1.0}}; + + const std::vector::vertices_per_cell>> + cell_vertices = {{{0, 1, 5, 6}}, + {{1, 2, 6, 7}}, + {{2, 3, 7, 8}}, + {{3, 4, 8, 9}}, + {{5, 6, 10, 11}}, + {{8, 9, 12, 13}}, + {{10, 11, 14, 15}}, + {{12, 13, 17, 18}}, + {{14, 15, 19, 20}}, + {{15, 16, 20, 21}}, + {{16, 17, 21, 22}}, + {{17, 18, 22, 23}}}; + + const unsigned int n_cells = cell_vertices.size(); + + std::vector> cells(n_cells, CellData()); + for (unsigned int i = 0; i < n_cells; ++i) + { + for (unsigned int j = 0; j < GeometryInfo::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } + + triangulation.create_triangulation(vertices, cells, SubCellData()); + triangulation.refine_global(3); + } + + + + template + void + LaplaceProblem::run() + { + for (unsigned int cycle = 0; cycle < 5; ++cycle) + { + pcout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + create_coarse_grid(); + + setup_system(); + + pcout << " Number of active cells : " + << triangulation.n_global_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl + << " Number of constraints : " + << Utilities::MPI::sum(constraints.n_constraints(), + mpi_communicator) + << std::endl; + + assemble_system(); + solve(); + postprocess(); + } + } + + + + template + std::pair + LaplaceProblem::predicate(const TableIndices &ind) + { + unsigned int v = 0; + for (unsigned int i = 0; i < dim; i++) + v += ind[i] * ind[i]; + if (v > 0 && v < max_degree * max_degree) + return std::make_pair(true, v); + else + return std::make_pair(false, v); + } + + + + template + void + LaplaceProblem::estimate_smoothness(Vector &smoothness_indicators) + { + Vector local_dof_values; + + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + { + local_dof_values.reinit(cell->get_fe().dofs_per_cell); + cell->get_dof_values(solution, local_dof_values); + + fourier->calculate(local_dof_values, + cell->active_fe_index(), + fourier_coefficients); + + std::pair, std::vector> res = + FESeries::process_coefficients( + fourier_coefficients, + std::bind(&LaplaceProblem::predicate, + this, + std::placeholders::_1), + VectorTools::Linfty_norm); + + Assert(res.first.size() == res.second.size(), ExcInternalError()); + + if (ln_k.size() == 0) + { + ln_k.resize(res.first.size(), 0); + for (unsigned int f = 0; f < ln_k.size(); f++) + ln_k[f] = + std::log(2.0 * numbers::PI * std::sqrt(1. * res.first[f])); + } + + for (double &residual_element : res.second) + residual_element = std::log(residual_element); + + std::pair fit = + FESeries::linear_regression(ln_k, res.second); + + smoothness_indicators(cell->active_cell_index()) = + -fit.first - 1. * dim / 2; + } + } +} // namespace Step27 + + + +int +main(int argc, char *argv[]) +{ + try + { + using namespace dealii; + using namespace Step27; + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + LaplaceProblem<2> laplace_problem; + laplace_problem.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/mpi/trilinos_step-27.with_p4est=true.with_trilinos=true.mpirun=1.output b/tests/mpi/trilinos_step-27.with_p4est=true.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..653b0af9f9 --- /dev/null +++ b/tests/mpi/trilinos_step-27.with_p4est=true.with_trilinos=true.mpirun=1.output @@ -0,0 +1,25 @@ +Cycle 0: + Number of active cells : 768 + Number of degrees of freedom: 3264 + Number of constraints : 384 + Solved in 16 iterations. +Cycle 1: + Number of active cells : 966 + Number of degrees of freedom: 5237 + Number of constraints : 932 + Solved in 25 iterations. +Cycle 2: + Number of active cells : 1146 + Number of degrees of freedom: 8538 + Number of constraints : 1958 + Solved in 34 iterations. +Cycle 3: + Number of active cells : 1359 + Number of degrees of freedom: 12513 + Number of constraints : 3133 + Solved in 45 iterations. +Cycle 4: + Number of active cells : 1656 + Number of degrees of freedom: 18527 + Number of constraints : 4846 + Solved in 75 iterations. diff --git a/tests/mpi/trilinos_step-27.with_p4est=true.with_trilinos=true.mpirun=2.output b/tests/mpi/trilinos_step-27.with_p4est=true.with_trilinos=true.mpirun=2.output new file mode 100644 index 0000000000..1f605c13d5 --- /dev/null +++ b/tests/mpi/trilinos_step-27.with_p4est=true.with_trilinos=true.mpirun=2.output @@ -0,0 +1,25 @@ +Cycle 0: + Number of active cells : 768 + Number of degrees of freedom: 3264 + Number of constraints : 404 + Solved in 16 iterations. +Cycle 1: + Number of active cells : 966 + Number of degrees of freedom: 5237 + Number of constraints : 952 + Solved in 25 iterations. +Cycle 2: + Number of active cells : 1146 + Number of degrees of freedom: 8538 + Number of constraints : 1999 + Solved in 34 iterations. +Cycle 3: + Number of active cells : 1359 + Number of degrees of freedom: 12513 + Number of constraints : 3237 + Solved in 46 iterations. +Cycle 4: + Number of active cells : 1656 + Number of degrees of freedom: 18535 + Number of constraints : 4975 + Solved in 70 iterations.