From eafe25f27007da1bdbdc7338559ed7a596912050 Mon Sep 17 00:00:00 2001 From: Mike H Date: Fri, 24 Jun 2016 14:47:41 -0500 Subject: [PATCH] Added a unit test for the refactored trilinos direct solver. --- tests/trilinos/direct_solver_2.cc | 334 ++++++++++++++++++++++++++ tests/trilinos/direct_solver_2.output | 13 + 2 files changed, 347 insertions(+) create mode 100644 tests/trilinos/direct_solver_2.cc create mode 100644 tests/trilinos/direct_solver_2.output diff --git a/tests/trilinos/direct_solver_2.cc b/tests/trilinos/direct_solver_2.cc new file mode 100644 index 0000000000..55c0866662 --- /dev/null +++ b/tests/trilinos/direct_solver_2.cc @@ -0,0 +1,334 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2015 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. +// +// --------------------------------------------------------------------- + + + +// 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 + + +template +class Step4 +{ +public: + Step4 (); + void run (); + +private: + void make_grid (); + void setup_system(); + void assemble_system (); + void solve (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + ConstraintMatrix constraints; + + TrilinosWrappers::SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +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 RightHandSide::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, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4 () + : + 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(); + + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false); + system_matrix.reinit (c_sparsity); + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); +} + + +template +void Step4::assemble_system () +{ + QGauss quadrature_formula(fe.degree+1); + + const RightHandSide right_hand_side; + + 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); + + 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) + { + 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); + } + system_matrix.compress(VectorOperation::add); +} + + + +template +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, + 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(); + } +} + + + +template +void Step4::run() +{ + make_grid(); + setup_system(); + assemble_system(); + solve(); +} + + +int main (int argc, char **argv) +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, testing_max_num_threads()); + + try + { + Step4<2> test; + test.run(); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/trilinos/direct_solver_2.output b/tests/trilinos/direct_solver_2.output new file mode 100644 index 0000000000..9b38539243 --- /dev/null +++ b/tests/trilinos/direct_solver_2.output @@ -0,0 +1,13 @@ + +DEAL:cg::Starting value 21.8299 +DEAL:cg::Convergence step 86 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. -- 2.39.5