From f77ab6e095321f47f691d44cc890dbff073ca984 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Sat, 12 Dec 2015 21:08:37 +0100 Subject: [PATCH] SLEPc unit test which uses PETSc preconditioner and solver inside SLEPc --- tests/slepc/step-36_parallel.cc | 332 ++++++++++++++++++ ...petsc=true.with_slepc=true.mpirun=3.output | 6 + 2 files changed, 338 insertions(+) create mode 100644 tests/slepc/step-36_parallel.cc create mode 100644 tests/slepc/step-36_parallel.with_mpi=true.with_petsc=true.with_slepc=true.mpirun=3.output diff --git a/tests/slepc/step-36_parallel.cc b/tests/slepc/step-36_parallel.cc new file mode 100644 index 0000000000..ab8ae7bed1 --- /dev/null +++ b/tests/slepc/step-36_parallel.cc @@ -0,0 +1,332 @@ +#include "../tests.h" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + + +#include +#include +#include +#include + +#include +#include + +// test parallel (MPI) version of Step-36 + +const unsigned int dim = 2;//run in 2d to save time + +const double eps = 1e-10; + +void test (std::string solver_name) +{ + const unsigned int global_mesh_refinement_steps = 5; + const unsigned int number_of_eigenvalues = 5; + + MPI_Comm mpi_communicator = MPI_COMM_WORLD; + const unsigned int n_mpi_processes = dealii::Utilities::MPI::n_mpi_processes(mpi_communicator); + const unsigned int this_mpi_process = dealii::Utilities::MPI::this_mpi_process(mpi_communicator); + + + dealii::Triangulation triangulation; + dealii::DoFHandler dof_handler(triangulation); + dealii::FE_Q fe(1); + dealii::ConstraintMatrix constraints; + dealii::IndexSet locally_owned_dofs; + dealii::IndexSet locally_relevant_dofs; + + std::vector eigenfunctions; + std::vector eigenvalues; + dealii::PETScWrappers::MPI::SparseMatrix stiffness_matrix, mass_matrix; + + dealii::GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (global_mesh_refinement_steps); + + // we do not use metis but rather partition by hand below. + //dealii::GridTools::partition_triangulation (n_mpi_processes, triangulation); + { + const double x0 = -1.0; + const double x1 = 1.0; + const double dL = (x1-x0) / n_mpi_processes; + + dealii::Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + { + const dealii::Point ¢er = cell->center(); + const double x = center[0]; + + const unsigned int id = std::floor ( (x-x0)/dL); + cell->set_subdomain_id (id); + } + } + + dof_handler.distribute_dofs (fe); + dealii::DoFRenumbering::subdomain_wise (dof_handler); + std::vector locally_owned_dofs_per_processor + = DoFTools::locally_owned_dofs_per_subdomain (dof_handler); + locally_owned_dofs = locally_owned_dofs_per_processor[this_mpi_process]; + locally_relevant_dofs.clear(); + dealii::DoFTools::extract_locally_relevant_dofs (dof_handler, + locally_relevant_dofs); + + constraints.clear(); + constraints.reinit (locally_relevant_dofs); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + dealii::VectorTools::interpolate_boundary_values (dof_handler, + 0, + dealii::ZeroFunction (), + constraints); + constraints.close (); + + dealii::CompressedSimpleSparsityPattern csp (locally_relevant_dofs); + // Fill in ignoring all cells that are not locally owned + dealii::DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, + /* keep constrained dofs */ true); + std::vector n_locally_owned_dofs(n_mpi_processes); + for (unsigned int i = 0; i < n_mpi_processes; i++) + n_locally_owned_dofs[i] = locally_owned_dofs_per_processor[i].n_elements(); + + dealii::SparsityTools::distribute_sparsity_pattern + (csp, + n_locally_owned_dofs, + mpi_communicator, + locally_relevant_dofs); + + // Initialise the stiffness and mass matrices + stiffness_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + mass_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + eigenfunctions.resize (5); + for (unsigned int i=0; i quadrature_formula(2); + dealii::FEValues fe_values (fe, quadrature_formula, + dealii::update_values | + dealii::update_gradients | + dealii::update_quadrature_points | + dealii::update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + dealii::FullMatrix cell_stiffness_matrix (dofs_per_cell, dofs_per_cell); + dealii::FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename dealii::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + if (cell->subdomain_id() == this_mpi_process) + { + fe_values.reinit (cell); + cell_stiffness_matrix = 0; + cell_mass_matrix = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + constraints + .distribute_local_to_global (cell_stiffness_matrix, + local_dof_indices, + stiffness_matrix); + constraints + .distribute_local_to_global (cell_mass_matrix, + local_dof_indices, + mass_matrix); + } + + stiffness_matrix.compress (dealii::VectorOperation::add); + mass_matrix.compress (dealii::VectorOperation::add); + + // test SLEPc by + { + PETScWrappers::PreconditionJacobi preconditioner(mpi_communicator); + dealii::SolverControl linear_solver_control (dof_handler.n_dofs(), 1e-12,/*log_history*/false,/*log_results*/false); + PETScWrappers::SolverCG linear_solver(linear_solver_control,mpi_communicator); + linear_solver.initialise(preconditioner); + + for (unsigned int i=0; i < eigenvalues.size(); i++) + eigenfunctions[i] = PetscScalar(); + + dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-9,/*log_history*/false,/*log_results*/false); + + dealii::SLEPcWrappers::SolverBase *eigensolver; + + // Get a handle on the wanted eigenspectrum solver + if (solver_name == "KrylovSchur") + { + eigensolver = new dealii::SLEPcWrappers::SolverKrylovSchur (solver_control, + mpi_communicator); + } + + else if (solver_name == "GeneralizedDavidson") + { + eigensolver = new dealii::SLEPcWrappers::SolverGeneralizedDavidson (solver_control, + mpi_communicator); + } + else if (solver_name == "JacobiDavidson") + { + eigensolver = new dealii::SLEPcWrappers::SolverJacobiDavidson (solver_control, + mpi_communicator); + } + else if (solver_name == "Lanczos") + { + eigensolver = new dealii::SLEPcWrappers::SolverLanczos (solver_control, + mpi_communicator); + } + else + { + AssertThrow (false, ExcMessage ("not supported eigensolver")); + + // Make compiler happy and not complaining about non + // uninitialized variables + eigensolver = new dealii::SLEPcWrappers::SolverKrylovSchur (solver_control, + mpi_communicator); + } + + SLEPcWrappers::TransformationShift spectral_transformation(mpi_communicator); + spectral_transformation.set_solver(linear_solver); + eigensolver->set_transformation(spectral_transformation); + //eigensolver->set_initial_vector(eigenfunctions[0]); + + eigensolver->set_which_eigenpairs (EPS_SMALLEST_REAL); + eigensolver->set_problem_type (EPS_GHEP); + + //zero constrained DoFs + for (unsigned int i=0; isolve (stiffness_matrix, + mass_matrix, + eigenvalues, + eigenfunctions, + eigenfunctions.size()); + + for (unsigned int i=0; i < eigenvalues.size(); i++) + dealii::deallog << eigenvalues[i] << std::endl; + + // make sure that we have eigenvectors and they are mass-orthonormal: + // a) (A*x_i-\lambda*B*x_i).L2() == 0 + // b) x_j*B*x_i=\delta_{ij} + { + const double precision = 1e-5; // ridiculously small... + PETScWrappers::MPI::Vector Ax(eigenfunctions[0]), Bx(eigenfunctions[0]); + for (unsigned int i=0; i < eigenfunctions.size(); ++i) + { + mass_matrix.vmult(Bx,eigenfunctions[i]); + + for (unsigned int j=0; j < eigenfunctions.size(); j++) + Assert( std::abs( eigenfunctions[j] * Bx - (i==j))< precision, + ExcMessage("Eigenvectors " + + Utilities::int_to_string(i) + + " and " + + Utilities::int_to_string(j) + + " are not orthonormal!")); + + stiffness_matrix.vmult(Ax,eigenfunctions[i]); + Ax.add(-1.0*eigenvalues[i],Bx); + Assert (Ax.l2_norm() < precision, + ExcMessage(std::to_string(Ax.l2_norm()))); + } + } + } + + + dof_handler.clear (); + dealii::deallog << "Ok"<