From 6dc1d32c4d16cbfd1ec5643907b43e8a9055a298 Mon Sep 17 00:00:00 2001 From: Roland Date: Mon, 14 Oct 2019 15:18:55 +0200 Subject: [PATCH] Simple test added, including test output --- tests/arpack/step-36_ar_with_iterations.cc | 373 ++++++++++++++++++ ...r_with_iterations.with_umfpack=true.output | 8 + 2 files changed, 381 insertions(+) create mode 100644 tests/arpack/step-36_ar_with_iterations.cc create mode 100644 tests/arpack/step-36_ar_with_iterations.with_umfpack=true.output diff --git a/tests/arpack/step-36_ar_with_iterations.cc b/tests/arpack/step-36_ar_with_iterations.cc new file mode 100644 index 0000000000..bf6cb2658b --- /dev/null +++ b/tests/arpack/step-36_ar_with_iterations.cc @@ -0,0 +1,373 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2009 - 2018 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.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * + * Authors: Toby D. Young, Polish Academy of Sciences, + * Wolfgang Bangerth, Texas A&M University + * + * This file tests the ARPACK interface for a symmetric operator taken from + step-36. + * + * We test that the computed vectors are eigenvectors and mass-orthonormal, i.e. + * a) (A*x_i-\lambda*B*x_i).L2() == 0 + * b) x_j*B*x_i = \delta_{i,j} + * + */ + +#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 "../tests.h" + +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; + + AffineConstraints constraints; + }; + + + + template + EigenvalueProblem::EigenvalueProblem(const std::string &prm_file) + : fe(1) + , dof_handler(triangulation) + {} + + + + template + void + EigenvalueProblem::make_grid_and_dofs() + { + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global(5); + 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(8); + for (unsigned int i = 0; i < eigenfunctions.size(); ++i) + eigenfunctions[i].reinit(dof_handler.n_dofs()); + + eigenvalues.resize(eigenfunctions.size()); + } + + + + template + 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); + + Functions::ZeroFunction potential; + + 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_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + cell_stiffness_matrix(i, j) += + (fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) + + potential_values[q_point] * + fe_values.shape_value(i, q_point) * + fe_values.shape_value(j, q_point)) * + fe_values.JxW(q_point); + + cell_mass_matrix(i, j) += (fe_values.shape_value(i, q_point) * + fe_values.shape_value(j, q_point)) * + fe_values.JxW(q_point); + } + + cell->get_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); + } + } + + + + template + std::pair + EigenvalueProblem::solve() + { + SolverControl solver_control(dof_handler.n_dofs(), 1e-10); + 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_magnitude, true); + ArpackSolver eigensolver(solver_control, additional_data); + eigensolver.solve(stiffness_matrix, + mass_matrix, + inverse, + eigenvalues, + eigenfunctions, + eigenvalues.size()); + + // 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} + { + 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++) + if (std::abs(eigenfunctions[j] * Bx - (i == j)) > 1e-8) + deallog << "Eigenvectors " + Utilities::int_to_string(i) + + " and " + Utilities::int_to_string(j) + + " are not orthonormal!" + " failing norm is " + + Utilities::to_string( + std::abs(eigenfunctions[j] * Bx - (i == j))) + << std::endl; + + stiffness_matrix.vmult(Ax, eigenfunctions[i]); + Ax.add(-1.0 * std::real(eigenvalues[i]), Bx); + if (Ax.l2_norm() > 1e-8) + deallog << "Returned vector " + Utilities::int_to_string(i) + + " is not an eigenvector!" + " L2 norm of the residual is " + + Utilities::to_string(Ax.l2_norm()) + << std::endl; + } + } + for (unsigned int i = 0; i < eigenfunctions.size(); ++i) + eigenfunctions[i] /= eigenfunctions[i].linfty_norm(); + + return std::make_pair(solver_control.last_step(), + solver_control.last_value()); + } + + + + template + void + EigenvalueProblem::output_results() const + { + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + + for (unsigned int i = 0; i < eigenfunctions.size(); ++i) + data_out.add_data_vector(eigenfunctions[i], + std::string("eigenfunction_") + + Utilities::int_to_string(i)); + + Vector projected_potential(dof_handler.n_dofs()); + { + FunctionParser potential; + potential.initialize(FunctionParser::default_variable_names(), + "0", + 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); + } + + + bool + my_compare(std::complex a, std::complex b) + { + return a.real() < b.real(); + } + + template + void + EigenvalueProblem::run() + { + make_grid_and_dofs(); + + assemble_system(); + + const std::pair res = solve(); + + deallog << res.first << " iterations used" << std::endl; + + std::sort(eigenvalues.begin(), eigenvalues.end(), my_compare); + + for (unsigned int i = 0; i < 5 && i < eigenvalues.size(); ++i) + deallog << " Eigenvalue " << i << " : " << eigenvalues[i] + << std::endl; + } +} // namespace Step36 + +int +main(int argc, char **argv) +{ + try + { + using namespace dealii; + using namespace Step36; + + initlog(); + + + EigenvalueProblem<2> 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; + } + + return 0; +} diff --git a/tests/arpack/step-36_ar_with_iterations.with_umfpack=true.output b/tests/arpack/step-36_ar_with_iterations.with_umfpack=true.output new file mode 100644 index 0000000000..89431ddcc6 --- /dev/null +++ b/tests/arpack/step-36_ar_with_iterations.with_umfpack=true.output @@ -0,0 +1,8 @@ + +DEAL::Convergence step 8 value 0.00000 +DEAL::8 iterations used +DEAL:: Eigenvalue 0 : (4.93877,0.00000) +DEAL:: Eigenvalue 1 : (12.3707,0.00000) +DEAL:: Eigenvalue 2 : (12.3707,0.00000) +DEAL:: Eigenvalue 3 : (19.8027,0.00000) +DEAL:: Eigenvalue 4 : (24.8370,0.00000) -- 2.39.5