From: Denis Davydov Date: Fri, 21 Oct 2016 13:29:26 +0000 (+0200) Subject: add ArpackSolver::set_shift() X-Git-Tag: v8.5.0-rc1~544^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7c8b7214f720f0959b4ded99b3483c478f4c20e7;p=dealii.git add ArpackSolver::set_shift() --- diff --git a/doc/news/changes.h b/doc/news/changes.h index 59663ea22a..7edcfd2daa 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -392,6 +392,12 @@ inconvenience this causes.
    +
  1. New: Add ArpackSolver::set_shift() to set the shift value in spectral + transformation. +
    + (Denis Davydov, 2016/10/25) +
    2. +
  2. New: PreconditionChebyshev now offers a PreconditionChebyshev::step() and PreconditionChebyshev::Tstep() methods for usage in relaxation smoothers.
    diff --git a/include/deal.II/lac/arpack_solver.h b/include/deal.II/lac/arpack_solver.h index b3507eb7c3..e43e13000f 100644 --- a/include/deal.II/lac/arpack_solver.h +++ b/include/deal.II/lac/arpack_solver.h @@ -161,6 +161,14 @@ public: template void set_initial_vector(const VectorType &vec); + /** + * Set real @p sigmar and complex @p sigmai parts of the shift for + * shift-and-invert spectral transformation. + * + * If this function is not called, the shift is assumed to be zero. + */ + void set_shift(const double sigmar, const double sigmai = 0); + /** * Solve the generalized eigensprectrum problem $A x=\lambda B x$ by calling * the dsaupd and dseupd or @@ -238,6 +246,17 @@ protected: bool initial_vector_provided; std::vector resid; + /** + * Real part of the shift + */ + double sigmar; + + /** + * Imaginary part of the shift + */ + double sigmai; + + private: /** @@ -333,9 +352,24 @@ ArpackSolver::ArpackSolver (SolverControl &control, : solver_control (control), additional_data (data), - initial_vector_provided(false) + initial_vector_provided(false), + sigmar(0.0), + sigmai(0.0) {} + + + +inline +void +ArpackSolver::set_shift(const double r, const double i) +{ + sigmar = r; + sigmai = i; +} + + + template inline void ArpackSolver:: @@ -601,9 +635,6 @@ void ArpackSolver::solve (const MatrixType1 &/*system_matrix*/, int ldz = n; - double sigmar = 0.0; // real part of the shift - double sigmai = 0.0; // imaginary part of the shift - std::vector eigenvalues_real (nev+1, 0.); std::vector eigenvalues_im (nev+1, 0.); diff --git a/tests/arpack/step-36_ar_shift.cc b/tests/arpack/step-36_ar_shift.cc new file mode 100644 index 0000000000..7a235af75f --- /dev/null +++ b/tests/arpack/step-36_ar_shift.cc @@ -0,0 +1,375 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2016 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. + * + * --------------------------------------------------------------------- + + * + * Same as step-36_ar but with a non-zero shift. + */ + +#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 + +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; + + ConstraintMatrix 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 + 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", + 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); + } + } + + + + template + std::pair EigenvalueProblem::solve () + { + SolverControl solver_control (dof_handler.n_dofs(), 1e-10); + + const double shift = 4.0; + const auto op_H = linear_operator >(stiffness_matrix); + const auto op_M = linear_operator >(mass_matrix); + const auto op_shift = op_H - shift * op_M; + + SolverControl solver_control_lin (dof_handler.n_dofs(), 1e-12, false, false); + SolverCG > cg(solver_control_lin); + const auto op_shift_invert = inverse_operator(op_shift, cg, PreconditionIdentity ()); + + 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.set_shift(shift); + eigensolver.solve (stiffness_matrix, + mass_matrix, + op_shift_invert, + 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 + 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", + 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 (); + + 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; + } +} + +int main (int argc, char **argv) +{ + + try + { + using namespace dealii; + using namespace Step36; + + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + + 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_shift.with_cxx11=on.output b/tests/arpack/step-36_ar_shift.with_cxx11=on.output new file mode 100644 index 0000000000..684685d6fb --- /dev/null +++ b/tests/arpack/step-36_ar_shift.with_cxx11=on.output @@ -0,0 +1,6 @@ + +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)