compute_eigenvalues.h

/*  ______________________________________________________________________
 *
 *  ExaDG - High-Order Discontinuous Galerkin for the Exa-Scale
 *
 *  Copyright (C) 2021 by the ExaDG authors
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 *  ______________________________________________________________________
 */
 
#ifndef EXADG_SOLVERS_AND_PRECONDITIONERS_UTILITIES_COMPUTE_EIGENVALUES_H_
#define EXADG_SOLVERS_AND_PRECONDITIONERS_UTILITIES_COMPUTE_EIGENVALUES_H_
 
// deal.II
#include <deal.II/numerics/vector_tools_mean_value.h>
 
namespace ExaDG
{
// manually compute eigenvalues for the coarsest level for proper setup of the
// Chebyshev iteration
template<typename Operator, typename VectorType>
std::pair<double, double>
compute_eigenvalues(Operator const &   op,
                    VectorType const & inverse_diagonal,
                    bool const         operator_is_singular,
                    unsigned int const eig_n_iter = 10000)
{
  VectorType solution, rhs;
  solution.reinit(inverse_diagonal);
  rhs.reinit(inverse_diagonal, true);
  // NB: initialize rand in order to obtain "reproducible" results!
  srand(1);
  for(unsigned int i = 0; i < rhs.locally_owned_size(); ++i)
    rhs.local_element(i) = (double)rand() / RAND_MAX;
  if(operator_is_singular)
    dealii::VectorTools::subtract_mean_value(rhs);
 
  dealii::SolverControl               control(eig_n_iter, rhs.l2_norm() * 1e-5);
  dealii::internal::EigenvalueTracker eigenvalue_tracker;
 
  dealii::SolverCG<VectorType> solver(control);
 
  solver.connect_eigenvalues_slot(std::bind(&dealii::internal::EigenvalueTracker::slot,
                                            &eigenvalue_tracker,
                                            std::placeholders::_1));
 
  JacobiPreconditioner<Operator> preconditioner(op, true /* initialize */);
 
  try
  {
    solver.solve(op, solution, rhs, preconditioner);
  }
  catch(dealii::SolverControl::NoConvergence &)
  {
  }
 
  std::pair<double, double> eigenvalues;
  if(eigenvalue_tracker.values.empty())
  {
    eigenvalues.first = eigenvalues.second = 1.;
  }
  else
  {
    eigenvalues.first  = eigenvalue_tracker.values.front();
    eigenvalues.second = eigenvalue_tracker.values.back();
  }
 
  return eigenvalues;
}
 
template<typename Number>
struct EigenvalueTracker
{
public:
  void
  slot(const std::vector<Number> & eigenvalues)
  {
    values = eigenvalues;
  }
 
  std::vector<Number> values;
};
 
} // namespace ExaDG
 
#endif /* EXADG_SOLVERS_AND_PRECONDITIONERS_UTILITIES_COMPUTE_EIGENVALUES_H_ */