From f54ddf3774749b71f3666fafbf541501c2eef6b0 Mon Sep 17 00:00:00 2001 From: Mike H Date: Tue, 28 Jun 2016 11:06:03 -0500 Subject: [PATCH] Added new feactures of trilinos direct solver to changes.h --- doc/news/changes.h | 6 + tests/trilinos/direct_solver_2.cc | 312 ++++++++++++++++---------- tests/trilinos/direct_solver_2.output | 16 +- 3 files changed, 211 insertions(+), 123 deletions(-) diff --git a/doc/news/changes.h b/doc/news/changes.h index cef4132d1a..7c4ec9d931 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -139,6 +139,12 @@ inconvenience this causes.

General

    +
  1. New: Added TrilinosWrappers::SolveDirect::Initialize and + TrilinosWrappers::SolverDirect::Solve to solve distributed linear systems + with multiple right hand sides without needing to refactorize the matrix + everytime. Also, added unit test for testing the new functionality. + (Michael Harmon, 2016/06/30) +
  2. New: Add new classes to expand a scalar finite element solution into the orthogonal bases FESeries::Fourier and FESeries::Legendre. Also diff --git a/tests/trilinos/direct_solver_2.cc b/tests/trilinos/direct_solver_2.cc index 55c0866662..8ae8c32863 100644 --- a/tests/trilinos/direct_solver_2.cc +++ b/tests/trilinos/direct_solver_2.cc @@ -13,30 +13,43 @@ // // --------------------------------------------------------------------- - - // tests Trilinos direct solvers on a 2D Poisson equation for linear elements #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 -#include -#include -#include +#include +#include + +#include #include #include +using namespace dealii; template class Step4 @@ -51,16 +64,18 @@ private: void assemble_system (); void solve (); - Triangulation triangulation; - FE_Q fe; - DoFHandler dof_handler; + parallel::distributed::Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; - ConstraintMatrix constraints; + ConstraintMatrix constraints; + SparsityPattern sparsity_pattern; - TrilinosWrappers::SparseMatrix system_matrix; + TrilinosWrappers::SparseMatrix system_matrix; - Vector solution; - Vector system_rhs; + TrilinosWrappers::MPI::Vector solution; + TrilinosWrappers::MPI::Vector system_rhs; + TrilinosWrappers::MPI::Vector system_rhs_two; }; @@ -75,6 +90,18 @@ public: }; +template +class RightHandSideTwo : public Function +{ +public: + RightHandSideTwo () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + + template class BoundaryValues : public Function @@ -95,12 +122,23 @@ double RightHandSide::value (const Point &p, { double return_value = 0; for (unsigned int i=0; i +double RightHandSideTwo::value (const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0; + for (unsigned int i=0; i double BoundaryValues::value (const Point &p, @@ -114,6 +152,10 @@ double BoundaryValues::value (const Point &p, template Step4::Step4 () : + triangulation(MPI_COMM_WORLD, + typename Triangulation::MeshSmoothing + (Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening)), fe (1), dof_handler (triangulation) {} @@ -141,12 +183,38 @@ void Step4::setup_system () constraints); constraints.close(); - CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false); - system_matrix.reinit (c_sparsity); + 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_relevant_dofs, + MPI_COMM_WORLD); + + system_rhs.reinit (locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD, + true); + + system_rhs_two.reinit (locally_owned_dofs, + locally_relevant_dofs, + MPI_COMM_WORLD, + true); - solution.reinit (dof_handler.n_dofs()); - system_rhs.reinit (dof_handler.n_dofs()); } @@ -166,6 +234,7 @@ void Step4::assemble_system () 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); @@ -175,29 +244,44 @@ void Step4::assemble_system () for (; cell!=endc; ++cell) { - fe_values.reinit (cell); - cell_matrix = 0; - cell_rhs = 0; - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - constraints.distribute_local_to_global(cell_matrix, cell_rhs, - local_dof_indices, - system_matrix, system_rhs); + 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, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + + constraints.distribute_local_to_global(cell_rhs_two, + local_dof_indices, + system_rhs_two); + } } system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + system_rhs_two.compress(VectorOperation::add); } @@ -207,79 +291,79 @@ void Step4::solve () { // Compute 'reference' solution with CG solver and SSOR preconditioner TrilinosWrappers::PreconditionSSOR preconditioner; - solution = 0; - SolverControl solver_control (1000, 1e-12); - SolverCG<> solver (solver_control); preconditioner.initialize(system_matrix); - solver.solve (system_matrix, solution, system_rhs, + TrilinosWrappers::MPI::Vector temp_solution(system_rhs); + temp_solution = 0; + SolverControl solver_control (1000, 1e-12); + SolverCG solver (solver_control); + + solver.solve (system_matrix, temp_solution, system_rhs, + preconditioner); + + constraints.distribute(temp_solution); + solution = temp_solution; + + TrilinosWrappers::MPI::Vector output(temp_solution); + + // do CG solve for new rhs + temp_solution = 0; + solution = 0; + solver.solve (system_matrix, temp_solution, system_rhs_two, preconditioner); - Vector output(solution); - { - deallog.push("DirectKLU"); - TrilinosWrappers::SolverDirect::AdditionalData data; - data.solver_type = "Amesos_Klu"; - SolverControl solver_control (1000, 1e-10); - TrilinosWrappers::SolverDirect solver(solver_control, data); - solver.initialize (system_matrix); - solver.solve (system_matrix, output, system_rhs); - output -= solution; - deallog << "Norm of error in direct solve: " << output.l2_norm() - << std::endl; - deallog.pop(); - } - - { - deallog.push("DirectLAPACK"); - TrilinosWrappers::SolverDirect::AdditionalData data; - data.solver_type = "Amesos_Lapack"; - SolverControl solver_control (1000, 1e-10); - TrilinosWrappers::SolverDirect solver(solver_control, data); - solver.initialize (system_matrix); - solver.solve (output, system_rhs); - output -= solution; - deallog << "Norm of error in direct solve: " << output.l2_norm() - << std::endl; - deallog.pop(); - } - - { - deallog.push("DirectDscpack"); - TrilinosWrappers::SolverDirect::AdditionalData data; - data.solver_type = "Amesos_Dscpack"; - SolverControl solver_control (1000, 1e-10); - TrilinosWrappers::SolverDirect solver(solver_control, data); - try - { - solver.initialize (system_matrix); - solver.solve (output, system_rhs); - } - catch (dealii::ExceptionBase &exc) - { - deallog << "Error: " << std::endl; - exc.print_info(deallog.get_file_stream()); - } - deallog.pop(); - } - - { - deallog.push("Dummy"); - TrilinosWrappers::SolverDirect::AdditionalData data; - data.solver_type = "DummySolver"; - SolverControl solver_control (1000, 1e-10); - TrilinosWrappers::SolverDirect solver(solver_control, data); - try - { - solver.initialize (system_matrix); - solver.solve (output, system_rhs); - } - catch (dealii::ExceptionBase &exc) - { - deallog << "Error: " << std::endl; - exc.print_info(deallog.get_file_stream()); - } - deallog.pop(); - } + constraints.distribute(temp_solution); + solution = temp_solution; + + TrilinosWrappers::MPI::Vector output_two(temp_solution); + + // factorize matrix for direct solver + temp_solution = 0; + 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 (temp_solution, system_rhs); + + constraints.distribute(temp_solution); + solution = temp_solution; + + // calculate l2 errors + output.add(-1.0,temp_solution); + + const double local_error = output.l2_norm(); + const double global_error = std::sqrt(Utilities::MPI::sum( + local_error * local_error, + MPI_COMM_WORLD)); + + deallog << "Norm of error in direct solve 1: " << global_error + << std::endl; + + // do solve 2 without refactorizing + temp_solution = 0; + solution = 0; + direct_solver.solve (temp_solution, system_rhs_two); + + constraints.distribute(temp_solution); + solution = temp_solution; + + // calculate l2 errors + output_two.add(-1.0,temp_solution); + + const double local_error_two = output_two.l2_norm(); + const double global_error_two = std::sqrt(Utilities::MPI::sum( + local_error_two * local_error_two, + MPI_COMM_WORLD)); + + deallog << "Norm of error in direct solve 2: " << global_error_two + << std::endl; + + deallog.pop(); + } diff --git a/tests/trilinos/direct_solver_2.output b/tests/trilinos/direct_solver_2.output index 9b38539243..cef7644c01 100644 --- a/tests/trilinos/direct_solver_2.output +++ b/tests/trilinos/direct_solver_2.output @@ -1,13 +1,11 @@ -DEAL:cg::Starting value 21.8299 +DEAL:cg::Starting value 21.8027 DEAL:cg::Convergence step 86 value 0 +DEAL:cg::Starting value 0.0940512 +DEAL:cg::Convergence step 77 value 0 DEAL:DirectKLU::Starting value 0 DEAL:DirectKLU::Convergence step 0 value 0 -DEAL:DirectKLU::Norm of error in direct solve: 0 -DEAL:DirectLAPACK::Starting value 0 -DEAL:DirectLAPACK::Convergence step 0 value 0 -DEAL:DirectLAPACK::Norm of error in direct solve: 0 -DEAL:DirectDscpack::Error: -You tried to select the solver type but this solver is not supported by Trilinos either because it does not exist, or because Trilinos was not configured for its use. -DEAL:Dummy::Error: -You tried to select the solver type but this solver is not supported by Trilinos either because it does not exist, or because Trilinos was not configured for its use. +DEAL:DirectKLU::Norm of error in direct solve 1: 0 +DEAL:DirectKLU::Starting value 0 +DEAL:DirectKLU::Convergence step 0 value 0 +DEAL:DirectKLU::Norm of error in direct solve 2: 0 -- 2.39.5