From be671725f715dee9fbd2c5eda790af0897a1dcea Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 10 Nov 2017 15:33:31 +0100 Subject: [PATCH] Add test. --- tests/trilinos/direct_solver_3.cc | 313 ++++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=2.output | 8 + 2 files changed, 321 insertions(+) create mode 100644 tests/trilinos/direct_solver_3.cc create mode 100644 tests/trilinos/direct_solver_3.with_mpi=true.with_p4est=true.mpirun=2.output diff --git a/tests/trilinos/direct_solver_3.cc b/tests/trilinos/direct_solver_3.cc new file mode 100644 index 0000000000..dd8daa82d2 --- /dev/null +++ b/tests/trilinos/direct_solver_3.cc @@ -0,0 +1,313 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Trilinos direct solvers on a 2D Poisson equation for linear elements with +// deal.II's parallel vector + +#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 + + +using namespace dealii; + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void assemble_system (); + void solve (); + + parallel::distributed::Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; + + TrilinosWrappers::SparseMatrix system_matrix; + + LinearAlgebra::distributed::Vector solution; + LinearAlgebra::distributed::Vector system_rhs; + LinearAlgebra::distributed::Vector system_rhs_two; +}; + + +template +class RightHandSide : public Function +{ +public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + +template +class BoundaryValues : public Function +{ +public: + BoundaryValues () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + +template +double BoundaryValues::value (const Point &p, + const unsigned int /*component*/) const +{ + return -0.5/dim*p.norm_square(); +} + + + +template +Step4::Step4 () + : + triangulation(MPI_COMM_WORLD, + typename Triangulation::MeshSmoothing + (Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)), + fe (1), + dof_handler (triangulation) +{} + + +template +void Step4::make_grid () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (6); +} + + + +template +void Step4::setup_system () +{ + dof_handler.distribute_dofs (fe); + + constraints.clear(); + std::map boundary_values; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + BoundaryValues(), + constraints); + constraints.close(); + + IndexSet locally_owned_dofs = dof_handler.locally_owned_dofs(); + IndexSet locally_relevant_dofs; + + DoFTools::extract_locally_relevant_dofs(dof_handler, + locally_relevant_dofs); + + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false); + SparsityTools::distribute_sparsity_pattern(dsp, + dof_handler.n_locally_owned_dofs_per_processor(), + MPI_COMM_WORLD, + locally_relevant_dofs); + + system_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + dsp, + MPI_COMM_WORLD); + + solution.reinit (locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD); + + system_rhs.reinit (locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD); + + system_rhs_two.reinit (locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD); + +} + + +template +void Step4::assemble_system () +{ + QGauss quadrature_formula(fe.degree+1); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + Vector cell_rhs_two (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell!=endc; ++cell) + { + if (cell->is_locally_owned()) + { + fe_values.reinit (cell); + cell_matrix = 0; + cell_rhs = 0; + cell_rhs_two = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + local_dof_indices, + system_matrix); + + constraints.distribute_local_to_global(cell_rhs, + local_dof_indices, + system_rhs, cell_matrix); + constraints.distribute_local_to_global(cell_rhs_two, + local_dof_indices, + system_rhs_two, cell_matrix); + } + } + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + system_rhs_two.compress(VectorOperation::add); +} + + + +template +void Step4::solve () +{ + SolverControl solver_control(100, 1e-12); + // factorize matrix for direct solver + solution = 0; + + deallog.push("DirectKLU"); + TrilinosWrappers::SolverDirect::AdditionalData data; + data.solver_type = "Amesos_Klu"; + TrilinosWrappers::SolverDirect direct_solver(solver_control, data); + direct_solver.initialize (system_matrix); + + // do solve 1 + direct_solver.solve (solution, system_rhs); + deallog << "Vector norm: " << solution.l2_norm(); + + constraints.distribute(solution); + deallog << " and " << solution.l2_norm() << std::endl; + solution.update_ghost_values(); + + Vector cellwise_error(triangulation.n_active_cells()); + VectorTools::integrate_difference(dof_handler, solution, + BoundaryValues(), + cellwise_error, + QGauss(3), + VectorTools::L2_norm); + deallog << " L2 error: " << VectorTools::compute_global_error(triangulation, cellwise_error, + VectorTools::L2_norm) << std::endl; + + // do solve 2 without refactorizing + solution = 0; + direct_solver.solve (solution, system_rhs_two); + deallog << "Vector norm: " << solution.l2_norm(); + + constraints.distribute(solution); + deallog << " and " << solution.l2_norm() << std::endl; + + deallog.pop(); + +} + + + +template +void Step4::run() +{ + make_grid(); + setup_system(); + assemble_system(); + solve(); +} + + +int main (int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads()); + mpi_initlog(); + deallog << std::setprecision(10); + + Step4<2> test; + test.run(); +} diff --git a/tests/trilinos/direct_solver_3.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/trilinos/direct_solver_3.with_mpi=true.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..6a93b169e8 --- /dev/null +++ b/tests/trilinos/direct_solver_3.with_mpi=true.with_p4est=true.mpirun=2.output @@ -0,0 +1,8 @@ + +DEAL:DirectKLU::Starting value 0.000000000 +DEAL:DirectKLU::Convergence step 0 value 0.000000000 +DEAL:DirectKLU::Vector norm: 12.03422482 and 13.21739951 +DEAL:DirectKLU:: L2 error: 0.0001707045651 +DEAL:DirectKLU::Starting value 0.000000000 +DEAL:DirectKLU::Convergence step 0 value 0.000000000 +DEAL:DirectKLU::Vector norm: 11.65872572 and 12.87645014 -- 2.39.5