eigenvalues.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_gmres.h>
 
namespace MeltPoolDG
{
  template <typename number, typename VectorType, typename Operator>
  std::vector<std::complex<number>>
  estimate_eigenvalues_gmres(const Operator   &op,
                             const VectorType &b,
                             const unsigned    max_eigenvalues = 100);
 
  // --- template and inline functions --- //
 
  template <typename number, typename VectorType, typename Operator>
  std::vector<std::complex<number>>
  estimate_eigenvalues_gmres(const Operator    &op,
                             const VectorType  &b,
                             const unsigned int max_eigenvalues)
  {
    VectorType x(b);
    x.update_ghost_values();
 
    const unsigned iterations = std::min(b.size(), max_eigenvalues);
 
    // Set up the SolverControl to perform the specified number of iterations by
    // requiring an unreachable tolerance.
    dealii::SolverControl control_eigen(iterations, 0.0);
 
    // Make sure that GMRES never restarts by setting its basis size
    // accordingly.
    typename dealii::SolverGMRES<VectorType>::AdditionalData data(iterations + 2);
    data.orthogonalization_strategy =
      dealii::LinearAlgebra::OrthogonalizationStrategy::modified_gram_schmidt;
 
    dealii::SolverGMRES<VectorType>   eigen_solver(control_eigen, data);
    std::vector<std::complex<number>> eigen_values;
 
    eigen_solver.connect_eigenvalues_slot(
      [&eigen_values](const std::vector<std::complex<number>> &evs) { eigen_values = evs; }, true);
 
    try
      {
        // We do not precondition the system because GMRES estimates the
        // eigenvalues of the preconditioned system.
        eigen_solver.solve(op, x, b, dealii::PreconditionIdentity());
      }
    catch (const dealii::SolverControl::NoConvergence &exc)
      {
        // We don't expect convergence and this is fine. However, we can check
        // that we indeed got as many eigenvectors as iterations.
        AssertDimension(eigen_values.size(), iterations);
      }
 
    return eigen_values;
  }
} // namespace MeltPoolDG