From 85caea0a100c0259466e78768fa1013fbc724c16 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Sat, 22 Nov 2014 07:44:54 +0100 Subject: [PATCH] added a step-36 test to arpack category --- tests/CMakeLists.txt | 2 +- tests/arpack/step-36_ar.cc | 374 +++++++++++++++++++++++++++++++++ tests/arpack/step-36_ar.output | 5 + 3 files changed, 380 insertions(+), 1 deletion(-) create mode 100644 tests/arpack/step-36_ar.cc create mode 100644 tests/arpack/step-36_ar.output diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 35bfe8b82d..a59fcf8108 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -62,7 +62,7 @@ SET(_categories a-framework algorithms all-headers aniso base bits build_tests codim_one deal.II distributed_grids fe gla grid hp integrators lac lapack manifold matrix_free mesh_converter metis mpi multigrid mumps opencascade petsc serialization - slepc trilinos umfpack + slepc trilinos umfpack arpack ) # diff --git a/tests/arpack/step-36_ar.cc b/tests/arpack/step-36_ar.cc new file mode 100644 index 0000000000..db792017de --- /dev/null +++ b/tests/arpack/step-36_ar.cc @@ -0,0 +1,374 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2009 - 2014 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. + * + * --------------------------------------------------------------------- + + * + * Authors: Toby D. Young, Polish Academy of Sciences, + * Wolfgang Bangerth, Texas A&M University + */ + +#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 + +namespace Step36 +{ + using namespace dealii; + + + template + class EigenvalueProblem + { + public: + EigenvalueProblem (const std::string &prm_file); + void run (); + + private: + void make_grid_and_dofs (); + void assemble_system (); + std::pair solve (); + void output_results () const; + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix stiffness_matrix, mass_matrix; + std::vector > eigenfunctions; + std::vector> eigenvalues; + + //ParameterHandler parameters; + + ConstraintMatrix constraints; + }; + + + + template + EigenvalueProblem::EigenvalueProblem (const std::string &prm_file) + : + fe (1), + dof_handler (triangulation) + { +// parameters.declare_entry ("Global mesh refinement steps", "5", +// Patterns::Integer (0, 20), +// "The number of times the 1-cell coarse mesh should " +// "be refined globally for our computations."); +// parameters.declare_entry ("Number of eigenvalues/eigenfunctions", "5", +// Patterns::Integer (0, 100), +// "The number of eigenvalues/eigenfunctions " +// "to be computed."); +// parameters.declare_entry ("Potential", "0", +// Patterns::Anything(), +// "A functional description of the potential."); + +// parameters.read_input (prm_file); + } + + + + template + void EigenvalueProblem::make_grid_and_dofs () + { + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (5/*parameters.get_integer ("Global mesh refinement steps")*/); + dof_handler.distribute_dofs (fe); + + DoFTools::make_zero_boundary_constraints (dof_handler, constraints); + constraints.close (); + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); + constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + stiffness_matrix.reinit (sparsity_pattern); + mass_matrix.reinit (sparsity_pattern); + + eigenfunctions + .resize (5/*parameters.get_integer ("Number of eigenvalues/eigenfunctions")*/); + for (unsigned int i=0; i + void EigenvalueProblem::assemble_system () + { + QGauss quadrature_formula(2); + + 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_stiffness_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix cell_mass_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + FunctionParser potential; + potential.initialize (FunctionParser::default_variable_names (), + "0"/*parameters.get ("Potential")*/, + typename FunctionParser::ConstMap()); + + std::vector potential_values (n_q_points); + + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + cell_stiffness_matrix = 0; + cell_mass_matrix = 0; + + potential.value_list (fe_values.get_quadrature_points(), + potential_values); + + 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 (VectorOperation::add); + mass_matrix.compress (VectorOperation::add); + + + double min_spurious_eigenvalue = std::numeric_limits::max(), + max_spurious_eigenvalue = -std::numeric_limits::max(); + + for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i) + if (constraints.is_constrained(i)) + { + const double ev = stiffness_matrix(i,i)/mass_matrix(i,i); + min_spurious_eigenvalue = std::min (min_spurious_eigenvalue, ev); + max_spurious_eigenvalue = std::max (max_spurious_eigenvalue, ev); + } + +// std::cout << " Spurious eigenvalues are all in the interval " +// << "[" << min_spurious_eigenvalue << "," << max_spurious_eigenvalue << "]" +// << std::endl; + + } + + + + template + std::pair EigenvalueProblem::solve () + { + SolverControl solver_control (dof_handler.n_dofs(), 1e-9); + SparseDirectUMFPACK inverse; + inverse.initialize (stiffness_matrix); + const unsigned int num_arnoldi_vectors = 2*eigenvalues.size() + 2; + ArpackSolver::AdditionalData additional_data(num_arnoldi_vectors, + ArpackSolver::largest_real_part); + ArpackSolver eigensolver (solver_control, additional_data); + eigensolver.solve (stiffness_matrix, + mass_matrix, + inverse, + eigenvalues, + eigenfunctions, + eigenvalues.size()); + + for (unsigned int i=0; i + void EigenvalueProblem::output_results () const + { + DataOut data_out; + + data_out.attach_dof_handler (dof_handler); + + for (unsigned int i=0; i projected_potential (dof_handler.n_dofs()); + { + FunctionParser potential; + potential.initialize (FunctionParser::default_variable_names (), + "0"/*parameters.get ("Potential")*/, + typename FunctionParser::ConstMap()); + VectorTools::interpolate (dof_handler, potential, projected_potential); + } + data_out.add_data_vector (projected_potential, "interpolated_potential"); + + data_out.build_patches (); + + std::ofstream output ("eigenvectors.vtk"); + data_out.write_vtk (output); + } + + + + template + void EigenvalueProblem::run () + { + make_grid_and_dofs (); + +// std::cout << " Number of active cells: " +// << triangulation.n_active_cells () +// << std::endl +// << " Number of degrees of freedom: " +// << dof_handler.n_dofs () +// << std::endl; + + assemble_system (); + + const std::pair res = solve (); +// std::cout << " Solver converged in " << res.first +// << " iterations. Residual is " << res.second << std::endl; + +// output_results (); +// +// std::cout << std::endl; +// for (unsigned int i=0; i problem (""); + problem.run (); + } + } + + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + +// std::cout << std::endl +// << " Job done." +// << std::endl; + + return 0; +} diff --git a/tests/arpack/step-36_ar.output b/tests/arpack/step-36_ar.output new file mode 100644 index 0000000000..6c321a9872 --- /dev/null +++ b/tests/arpack/step-36_ar.output @@ -0,0 +1,5 @@ +DEAL:: Eigenvalue 0 : (4.93877,0.00000) +DEAL:: Eigenvalue 1 : (19.8027,0.00000) +DEAL:: Eigenvalue 2 : (12.3707,0.00000) +DEAL:: Eigenvalue 3 : (12.3707,0.00000) +DEAL:: Eigenvalue 4 : (24.8370,0.00000) -- 2.39.5