From: Matthias Maier Date: Fri, 11 Feb 2022 00:44:50 +0000 (-0600) Subject: add two performance tests X-Git-Tag: v9.4.0-rc1~508^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=124896a82daa8ceb4827f676aff7ca8eedb0c06d;p=dealii.git add two performance tests - Step 3 with SSOR preconditioner - Step 22 with increased workload --- diff --git a/tests/performance/timing_step_22.cc b/tests/performance/timing_step_22.cc new file mode 100644 index 0000000000..c5cccbecba --- /dev/null +++ b/tests/performance/timing_step_22.cc @@ -0,0 +1,675 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + +// +// Description: +// +// A performance benchmark based on step 3 that measures timings for system +// setup, assembly, solve and postprocessing for a Stokes problem. +// + +#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 "performance_test_driver.h" + +using namespace dealii; + +dealii::ConditionalOStream debug_output(std::cout, false); + +template +struct InnerPreconditioner; + +template <> +struct InnerPreconditioner<2> +{ + using type = SparseDirectUMFPACK; +}; + +template <> +struct InnerPreconditioner<3> +{ + using type = SparseILU; +}; + +template +class StokesProblem +{ +public: + StokesProblem(const unsigned int degree); + std::vector + run(); + +private: + void + setup_dofs(); + void + assemble_system(); + void + solve(); + void + output_results(const unsigned int refinement_cycle) const; + void + refine_mesh(); + + const unsigned int degree; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + AffineConstraints constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockSparsityPattern preconditioner_sparsity_pattern; + BlockSparseMatrix preconditioner_matrix; + + BlockVector solution; + BlockVector system_rhs; + + std::shared_ptr::type> A_preconditioner; +}; + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues() + : Function(dim + 1) + {} + + virtual double + value(const Point &p, const unsigned int component = 0) const override; + + virtual void + vector_value(const Point &p, Vector &value) const override; +}; + + +template +double +BoundaryValues::value(const Point & p, + const unsigned int component) const +{ + Assert(component < this->n_components, + ExcIndexRange(component, 0, this->n_components)); + + if (component == 0) + return (p[0] < 0 ? -1 : (p[0] > 0 ? 1 : 0)); + return 0; +} + + +template +void +BoundaryValues::vector_value(const Point &p, + Vector & values) const +{ + for (unsigned int c = 0; c < this->n_components; ++c) + values(c) = BoundaryValues::value(p, c); +} + + +template +class RightHandSide : public TensorFunction<1, dim> +{ +public: + RightHandSide() + : TensorFunction<1, dim>() + {} + + virtual Tensor<1, dim> + value(const Point &p) const override; + + virtual void + value_list(const std::vector> &p, + std::vector> & value) const override; +}; + + +template +Tensor<1, dim> +RightHandSide::value(const Point & /*p*/) const +{ + return Tensor<1, dim>(); +} + + +template +void +RightHandSide::value_list(const std::vector> &vp, + std::vector> & values) const +{ + for (unsigned int c = 0; c < vp.size(); ++c) + { + values[c] = RightHandSide::value(vp[c]); + } +} + + +template +class InverseMatrix : public Subscriptor +{ +public: + InverseMatrix(const MatrixType &m, const PreconditionerType &preconditioner); + + void + vmult(Vector &dst, const Vector &src) const; + +private: + const SmartPointer matrix; + const SmartPointer preconditioner; +}; + + +template +InverseMatrix::InverseMatrix( + const MatrixType & m, + const PreconditionerType &preconditioner) + : matrix(&m) + , preconditioner(&preconditioner) +{} + + +template +void +InverseMatrix::vmult( + Vector & dst, + const Vector &src) const +{ + SolverControl solver_control(src.size(), 1e-6 * src.l2_norm()); + SolverCG> cg(solver_control); + + dst = 0; + + cg.solve(*matrix, dst, src, *preconditioner); +} + + +template +class SchurComplement : public Subscriptor +{ +public: + SchurComplement( + const BlockSparseMatrix &system_matrix, + const InverseMatrix, PreconditionerType> &A_inverse); + + void + vmult(Vector &dst, const Vector &src) const; + +private: + const SmartPointer> system_matrix; + const SmartPointer< + const InverseMatrix, PreconditionerType>> + A_inverse; + + mutable Vector tmp1, tmp2; +}; + + +template +SchurComplement::SchurComplement( + const BlockSparseMatrix & system_matrix, + const InverseMatrix, PreconditionerType> &A_inverse) + : system_matrix(&system_matrix) + , A_inverse(&A_inverse) + , tmp1(system_matrix.block(0, 0).m()) + , tmp2(system_matrix.block(0, 0).m()) +{} + + +template +void +SchurComplement::vmult(Vector & dst, + const Vector &src) const +{ + system_matrix->block(0, 1).vmult(tmp1, src); + A_inverse->vmult(tmp2, tmp1); + system_matrix->block(1, 0).vmult(dst, tmp2); +} + + +template +StokesProblem::StokesProblem(const unsigned int degree) + : degree(degree) + , triangulation(Triangulation::maximum_smoothing) + , fe(FE_Q(degree + 1), dim, FE_Q(degree), 1) + , dof_handler(triangulation) +{} + + +template +void +StokesProblem::setup_dofs() +{ + A_preconditioner.reset(); + system_matrix.clear(); + preconditioner_matrix.clear(); + + dof_handler.distribute_dofs(fe); + DoFRenumbering::Cuthill_McKee(dof_handler); + + std::vector block_component(dim + 1, 0); + block_component[dim] = 1; + DoFRenumbering::component_wise(dof_handler, block_component); + + { + constraints.clear(); + + const FEValuesExtractors::Vector velocities(0); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values(dof_handler, + 1, + BoundaryValues(), + constraints, + fe.component_mask(velocities)); + } + + constraints.close(); + + const std::vector dofs_per_block = + DoFTools::count_dofs_per_fe_block(dof_handler, block_component); + const types::global_dof_index n_u = dofs_per_block[0]; + const types::global_dof_index n_p = dofs_per_block[1]; + + debug_output << " Number of active cells: " + << triangulation.n_active_cells() << std::endl + << " Number of degrees of freedom: " << dof_handler.n_dofs() + << " (" << n_u << '+' << n_p << ')' << std::endl; + + { + BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block); + + Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1); + for (unsigned int c = 0; c < dim + 1; ++c) + for (unsigned int d = 0; d < dim + 1; ++d) + if (!((c == dim) && (d == dim))) + coupling[c][d] = DoFTools::always; + else + coupling[c][d] = DoFTools::none; + + DoFTools::make_sparsity_pattern( + dof_handler, coupling, dsp, constraints, false); + + sparsity_pattern.copy_from(dsp); + } + + { + BlockDynamicSparsityPattern preconditioner_dsp(dofs_per_block, + dofs_per_block); + + Table<2, DoFTools::Coupling> preconditioner_coupling(dim + 1, dim + 1); + for (unsigned int c = 0; c < dim + 1; ++c) + for (unsigned int d = 0; d < dim + 1; ++d) + if (((c == dim) && (d == dim))) + preconditioner_coupling[c][d] = DoFTools::always; + else + preconditioner_coupling[c][d] = DoFTools::none; + + DoFTools::make_sparsity_pattern(dof_handler, + preconditioner_coupling, + preconditioner_dsp, + constraints, + false); + + preconditioner_sparsity_pattern.copy_from(preconditioner_dsp); + } + + system_matrix.reinit(sparsity_pattern); + preconditioner_matrix.reinit(preconditioner_sparsity_pattern); + + solution.reinit(dofs_per_block); + system_rhs.reinit(dofs_per_block); +} + + +template +void +StokesProblem::assemble_system() +{ + system_matrix = 0; + system_rhs = 0; + preconditioner_matrix = 0; + + QGauss quadrature_formula(degree + 2); + + FEValues fe_values(fe, + quadrature_formula, + update_values | update_quadrature_points | + update_JxW_values | update_gradients); + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix local_preconditioner_matrix(dofs_per_cell, dofs_per_cell); + Vector local_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + const RightHandSide right_hand_side; + std::vector> rhs_values(n_q_points, Tensor<1, dim>()); + + const FEValuesExtractors::Vector velocities(0); + const FEValuesExtractors::Scalar pressure(dim); + + std::vector> symgrad_phi_u(dofs_per_cell); + std::vector div_phi_u(dofs_per_cell); + std::vector> phi_u(dofs_per_cell); + std::vector phi_p(dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + local_matrix = 0; + local_preconditioner_matrix = 0; + local_rhs = 0; + + right_hand_side.value_list(fe_values.get_quadrature_points(), rhs_values); + + for (unsigned int q = 0; q < n_q_points; ++q) + { + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + symgrad_phi_u[k] = fe_values[velocities].symmetric_gradient(k, q); + div_phi_u[k] = fe_values[velocities].divergence(k, q); + phi_u[k] = fe_values[velocities].value(k, q); + phi_p[k] = fe_values[pressure].value(k, q); + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j <= i; ++j) + { + local_matrix(i, j) += + (2 * (symgrad_phi_u[i] * symgrad_phi_u[j]) // (1) + - div_phi_u[i] * phi_p[j] // (2) + - phi_p[i] * div_phi_u[j]) // (3) + * fe_values.JxW(q); // * dx + + local_preconditioner_matrix(i, j) += + (phi_p[i] * phi_p[j]) // (4) + * fe_values.JxW(q); // * dx + } + local_rhs(i) += phi_u[i] // phi_u_i(x_q) + * rhs_values[q] // * f(x_q) + * fe_values.JxW(q); // * dx + } + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + { + local_matrix(i, j) = local_matrix(j, i); + local_preconditioner_matrix(i, j) = + local_preconditioner_matrix(j, i); + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global( + local_matrix, local_rhs, local_dof_indices, system_matrix, system_rhs); + constraints.distribute_local_to_global(local_preconditioner_matrix, + local_dof_indices, + preconditioner_matrix); + } + + debug_output << " Computing preconditioner..." << std::endl << std::flush; + + A_preconditioner = + std::make_shared::type>(); + A_preconditioner->initialize( + system_matrix.block(0, 0), + typename InnerPreconditioner::type::AdditionalData()); +} + + +template +void +StokesProblem::solve() +{ + const InverseMatrix, + typename InnerPreconditioner::type> + A_inverse(system_matrix.block(0, 0), *A_preconditioner); + Vector tmp(solution.block(0).size()); + + { + Vector schur_rhs(solution.block(1).size()); + A_inverse.vmult(tmp, system_rhs.block(0)); + system_matrix.block(1, 0).vmult(schur_rhs, tmp); + schur_rhs -= system_rhs.block(1); + + SchurComplement::type> schur_complement( + system_matrix, A_inverse); + + SolverControl solver_control(solution.block(1).size(), + 1e-6 * schur_rhs.l2_norm()); + SolverCG> cg(solver_control); + + SparseILU preconditioner; + preconditioner.initialize(preconditioner_matrix.block(1, 1), + SparseILU::AdditionalData()); + + InverseMatrix, SparseILU> m_inverse( + preconditioner_matrix.block(1, 1), preconditioner); + + cg.solve(schur_complement, solution.block(1), schur_rhs, m_inverse); + + constraints.distribute(solution); + + debug_output << " " << solver_control.last_step() + << " outer CG Schur complement iterations for pressure" + << std::endl; + } + + { + system_matrix.block(0, 1).vmult(tmp, solution.block(1)); + tmp *= -1; + tmp += system_rhs.block(0); + + A_inverse.vmult(solution.block(0), tmp); + + constraints.distribute(solution); + } +} + + +template +void +StokesProblem::output_results(const unsigned int refinement_cycle) const +{ + std::vector solution_names(dim, "velocity"); + solution_names.emplace_back("pressure"); + + std::vector + data_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation.push_back( + DataComponentInterpretation::component_is_scalar); + + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, + solution_names, + DataOut::type_dof_data, + data_component_interpretation); + data_out.build_patches(); + + std::ofstream output("solution-" + + Utilities::int_to_string(refinement_cycle, 2) + ".vtk"); + data_out.write_vtk(output); +} + + +template +void +StokesProblem::refine_mesh() +{ + Vector estimated_error_per_cell(triangulation.n_active_cells()); + + const FEValuesExtractors::Scalar pressure(dim); + KellyErrorEstimator::estimate( + dof_handler, + QGauss(degree + 1), + std::map *>(), + solution, + estimated_error_per_cell, + fe.component_mask(pressure)); + + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + estimated_error_per_cell, + 0.3, + 0.0); + triangulation.execute_coarsening_and_refinement(); +} + + +template +std::vector +StokesProblem::run() +{ + std::map timer; + + { + std::vector subdivisions(dim, 1); + subdivisions[0] = 4; + + const Point bottom_left = (dim == 2 ? // + Point(-2, -1) : // 2d case + Point(-2, 0, -1)); // 3d case + + const Point top_right = (dim == 2 ? // + Point(2, 0) : // 2d case + Point(2, 1, 0)); // 3d case + + GridGenerator::subdivided_hyper_rectangle(triangulation, + subdivisions, + bottom_left, + top_right); + } + + for (const auto &cell : triangulation.active_cell_iterators()) + for (const auto &face : cell->face_iterators()) + if (face->center()[dim - 1] == 0) + face->set_all_boundary_ids(1); + + + switch (get_testing_envrionment()) + { + case TestingEnvironment::light: + triangulation.refine_global(5 - dim); + break; + case TestingEnvironment::medium: + /* fallthrough */ + case TestingEnvironment::heavy: + triangulation.refine_global(6 - dim); + break; + } + + for (unsigned int refinement_cycle = 0; refinement_cycle < 6; + ++refinement_cycle) + { + debug_output << "Refinement cycle " << refinement_cycle << std::endl; + + if (refinement_cycle > 0) + { + timer["refinement"].start(); + refine_mesh(); + timer["refinement"].stop(); + } + + timer["setup_system"].start(); + setup_dofs(); + timer["setup_system"].stop(); + + debug_output << " Assembling..." << std::endl << std::flush; + timer["assemble_system"].start(); + assemble_system(); + timer["assemble_system"].stop(); + + timer["solve"].start(); + solve(); + timer["solve"].stop(); + + timer["output_results"].start(); + output_results(refinement_cycle); + timer["output_results"].stop(); + + debug_output << std::endl; + } + return {timer["refinement"].cpu_time(), + timer["setup_system"].cpu_time(), + timer["assemble_system"].cpu_time(), + timer["solve"].cpu_time(), + timer["output_results"].cpu_time()}; +} + + +std::tuple> +describe_measurements() +{ + return {Metric::timing, + 4, + {"refinement", + "setup_system", + "assemble_system", + "solve", + "output_results"}}; +} + + +Measurement +perform_single_measurement() +{ + return StokesProblem<2>(1).run(); +} diff --git a/tests/performance/timing_step_22.threads=2.release.run_only b/tests/performance/timing_step_22.threads=2.release.run_only new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/performance/timing_step_3.cc b/tests/performance/timing_step_3.cc new file mode 100644 index 0000000000..4c6186873f --- /dev/null +++ b/tests/performance/timing_step_3.cc @@ -0,0 +1,263 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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. +// +// --------------------------------------------------------------------- + +// +// Description: +// +// A performance benchmark based on step 3 that measures timings for system +// setup, assembly, solve and postprocessing for a simple Poisson problem. +// + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "performance_test_driver.h" + +using namespace dealii; + +dealii::ConditionalOStream debug_output(std::cout, false); + +class Step3 +{ +public: + Step3(); + + std::vector + run(); + +private: + void + make_grid(); + void + setup_system(); + void + assemble_system(); + void + solve(); + void + output_results() const; + + Triangulation<2> triangulation; + FE_Q<2> fe; + DoFHandler<2> dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +Step3::Step3() + : fe(1) + , dof_handler(triangulation) +{} + + +void +Step3::make_grid() +{ + GridGenerator::hyper_cube(triangulation, -1, 1); + + switch (get_testing_envrionment()) + { + case TestingEnvironment::light: + triangulation.refine_global(9); + break; + case TestingEnvironment::medium: + /* fallthrough */ + case TestingEnvironment::heavy: + triangulation.refine_global(10); + break; + } + + debug_output << "Number of active cells: " << triangulation.n_active_cells() + << std::endl; +} + + +void +Step3::setup_system() +{ + dof_handler.distribute_dofs(fe); + debug_output << "Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); +} + + +void +Step3::assemble_system() +{ + QGauss<2> quadrature_formula(fe.degree + 1); + FEValues<2> fe_values(fe, + quadrature_formula, + update_values | update_gradients | update_JxW_values); + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + + cell_matrix = 0; + cell_rhs = 0; + + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + { + for (const unsigned int i : fe_values.dof_indices()) + for (const unsigned int j : fe_values.dof_indices()) + cell_matrix(i, j) += + (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q) + fe_values.shape_grad(j, q_index) * // grad phi_j(x_q) + fe_values.JxW(q_index)); // dx + + for (const unsigned int i : fe_values.dof_indices()) + cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q) + 1. * // f(x_q) + fe_values.JxW(q_index)); // dx + } + cell->get_dof_indices(local_dof_indices); + + for (const unsigned int i : fe_values.dof_indices()) + for (const unsigned int j : fe_values.dof_indices()) + system_matrix.add(local_dof_indices[i], + local_dof_indices[j], + cell_matrix(i, j)); + + for (const unsigned int i : fe_values.dof_indices()) + system_rhs(local_dof_indices[i]) += cell_rhs(i); + } + + + std::map boundary_values; + VectorTools::interpolate_boundary_values(dof_handler, + 0, + Functions::ZeroFunction<2>(), + boundary_values); + MatrixTools::apply_boundary_values(boundary_values, + system_matrix, + solution, + system_rhs); +} + + +void +Step3::solve() +{ + SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm()); + SolverCG> solver(solver_control); + PreconditionSSOR> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + solver.solve(system_matrix, solution, system_rhs, preconditioner); +} + + +void +Step3::output_results() const +{ + DataOut<2> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + data_out.build_patches(); +} + + +std::vector +Step3::run() +{ + std::map timer; + + timer["make_grid"].start(); + make_grid(); + timer["make_grid"].stop(); + + timer["setup_system"].start(); + setup_system(); + timer["setup_system"].stop(); + + timer["assemble_system"].start(); + assemble_system(); + timer["assemble_system"].stop(); + + timer["solve"].start(); + solve(); + timer["solve"].stop(); + + timer["output_results"].start(); + output_results(); + timer["output_results"].stop(); + + return {timer["make_grid"].cpu_time(), + timer["setup_system"].cpu_time(), + timer["assemble_system"].cpu_time(), + timer["solve"].cpu_time(), + timer["output_results"].cpu_time()}; +} + + +std::tuple> +describe_measurements() +{ + return {Metric::timing, + 4, + {"make_grid", + "setup_system", + "assemble_system", + "solve", + "output_results"}}; +} + + +Measurement +perform_single_measurement() +{ + return Step3().run(); +} diff --git a/tests/performance/timing_step_3.threads=2.release.run_only b/tests/performance/timing_step_3.threads=2.release.run_only new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/run_testsuite.cmake b/tests/run_testsuite.cmake index 52553a1e74..6652789304 100644 --- a/tests/run_testsuite.cmake +++ b/tests/run_testsuite.cmake @@ -271,7 +271,7 @@ ENDIF() # Pass all relevant variables down to configure: GET_CMAKE_PROPERTY(_variables VARIABLES) FOREACH(_var ${_variables}) - IF( _var MATCHES "^(ENABLE|TEST|DEAL_II|ALLOW|WITH|FORCE|COMPONENT)_" OR + IF( _var MATCHES "^(ENABLE|TEST|TESTING|DEAL_II|ALLOW|WITH|FORCE|COMPONENT)_" OR _var MATCHES "^(DOCUMENTATION|EXAMPLES)" OR _var MATCHES "^(ADOLC|ARBORX|ARPACK|BOOST|OPENCASCADE|MUPARSER|HDF5|KOKKOS|METIS|MPI)_" OR _var MATCHES "^(GINKGO|P4EST|PETSC|SCALAPACK|SLEPC|THREADS|TBB|TRILINOS)_" OR