multigrid.hpp

Go to the documentation of this file.

#pragma once
 
#include <vector>
 
#include <Eigen/Core>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
 
#include <amg/common.hpp>
#include <amg/grid.hpp>
#include <amg/interpolator.hpp>
#include <amg/smoother.hpp>
 
namespace AMG {
 
template <class EleType>
class Multigrid {
 private:
  InterpolatorBase<EleType>* interpolator;
 
  AMG::SmootherBase<EleType>* smoother;
 
  Eigen::SimplicialLDLT<Eigen::SparseMatrix<EleType>> coarse_direct_solver;
 
  EleType tolerance;
 
  size_t compute_error_every_n_iters;
 
  size_t n_iters;
 
  size_t n_fine_nodes;
 
  size_t n_levels;
 
  size_t finest_grid_ix = 0;
 
  size_t coarsest_grid_ix;
 
  std::vector<size_t> level_to_n_dofs;
 
  std::vector<Eigen::SparseMatrix<EleType>> level_to_coefficient_matrix;
 
  std::vector<Eigen::Matrix<EleType, -1, 1>> level_to_rhs;
 
  std::vector<Eigen::Matrix<EleType, -1, 1>> level_to_soln;
 
  std::vector<Eigen::Matrix<EleType, -1, 1>> level_to_residual;
 
  size_t n_H_dofs_from_n_h_dofs(size_t h_dofs) {
    size_t H_dofs = (h_dofs + 1) / 2 - 1;
    return H_dofs;
  }
 
  bool display_error{false};
 
 public:
  ~Multigrid() {}
 
  Multigrid() = delete;
 
  Multigrid(AMG::InterpolatorBase<EleType>* interpolator_,
            AMG::SmootherBase<EleType>* smoother_,
            const Eigen::SparseMatrix<EleType>& A,
            const Eigen::Matrix<EleType, -1, 1>& b, size_t n_levels_,
            EleType tolerance_ = 1e-9, size_t compute_error_every_n_iters_ = 10,
            size_t n_iters_ = 100)
      : interpolator(interpolator_),
        smoother(smoother_),
        n_levels(n_levels_),
        tolerance(tolerance_),
        compute_error_every_n_iters(compute_error_every_n_iters_),
        n_iters(n_iters_) {
 
    std::string errormsg;
    if (compute_error_every_n_iters > n_iters) {
      errormsg =
          "`compute_error_every_n_iters` must be leq to `n_iters`, got " +
          std::to_string(compute_error_every_n_iters) + " and " +
          std::to_string(n_iters);
      throw(std::invalid_argument(errormsg));
    }
 
    if (A.rows() != b.rows()) {
      errormsg =
          "`A` and `b` must have the same number of degrees of freedom, got " +
          std::to_string(A.rows()) + " and " + std::to_string(b.rows());
      throw(std::invalid_argument(errormsg));
    }
 
    // Initialize linear system info
    level_to_coefficient_matrix.resize(n_levels);
    level_to_soln.resize(n_levels);
    level_to_rhs.resize(n_levels);
    level_to_residual.resize(n_levels);
    level_to_n_dofs.resize(n_levels);
 
    // Initialize index for coarsest grid
    coarsest_grid_ix = n_levels - 1;
 
    // Initialize the finest coefficients matrix
    size_t n_fine_dofs = A.rows();
    level_to_n_dofs[finest_grid_ix] = n_fine_dofs;
    level_to_coefficient_matrix[finest_grid_ix] = A;
 
    // Initialize the finest solutions vector
    Eigen::Matrix<EleType, -1, 1> u(n_fine_dofs);
    u.setZero();
    level_to_soln[finest_grid_ix] = u;
 
    // Initialize the finest right hand side
    level_to_rhs[finest_grid_ix] = b;
 
    // Initialize the finest residual
    level_to_residual[finest_grid_ix] = b - A * u;
 
    // Use the finest level grid info to construct the linear interpolators
    // needed to construct the remaining grid stuff
    // TODO: could do while coarsest is greater than max coarse ndofs
    size_t n_H_dofs;
    size_t n_h_dofs;
    for (size_t level = 1; level < n_levels; ++level) {
      // Make restriction/prolongation matrices
      n_h_dofs = level_to_n_dofs[level - 1];
      n_H_dofs = n_H_dofs_from_n_h_dofs(n_h_dofs);
      level_to_n_dofs[level] = n_H_dofs;
      interpolator->make_operators(n_h_dofs, n_H_dofs, level - 1);
 
      // Make coefficient matrix using interpolation operators
      auto R_h = interpolator->get_R(level - 1);
      auto A_h = level_to_coefficient_matrix[level - 1];  // finer matrix
      auto P_h = interpolator->get_P(level - 1);
      auto A_H = R_h * (A_h * P_h);  // coarser matrix
      level_to_coefficient_matrix[level] = A_H;
 
      // Make solution vector
      Eigen::Matrix<EleType, -1, 1> u_H(n_H_dofs);
      u_H.setZero();
      level_to_soln[level] = u_H;
 
      // Make right hand side
      Eigen::Matrix<EleType, -1, 1> rhs_H(n_H_dofs);
      rhs_H.setZero();
      level_to_rhs[level] = rhs_H;
 
      // Make residual vector
      level_to_residual[level] = rhs_H - A_H * u_H;
    }
 
    // Initialize the coarse solver
    coarse_direct_solver.analyzePattern(
        level_to_coefficient_matrix[coarsest_grid_ix]);
    coarse_direct_solver.factorize(
        level_to_coefficient_matrix[coarsest_grid_ix]);
  }
 
  void vcycle() {
    // At each level of the grid hiearchy, finest-to-coarsest:
    for (size_t level = 0; level < n_levels; ++level) {
      //  1. Apply a couple of smoothing iterations (pre-relaxation) to the
      // current solution ui=Si(Ai,fi,ui)
      smoother->smooth(level_to_coefficient_matrix[level], level_to_soln[level],
                       level_to_rhs[level]);
 
      //  2. Find residual ei=fi−Aiui and restrict it to the coarser level
      level_to_residual[level] =
          level_to_rhs[level] -
          level_to_coefficient_matrix[level] * level_to_soln[level];
 
      if (level + 1 != n_levels) {
        // initialize coarse solution to zero here refs [2-4].
        level_to_soln[level + 1].setZero();
 
        // f_{i+1} = R_{i} e_{i}
        level_to_rhs[level + 1] =
            interpolator->restriction(level_to_residual[level], level);
      }
    }
 
    // Solve the coarsest system directly: uL=A−1LfL
    level_to_soln[coarsest_grid_ix] =
        coarse_direct_solver.solve(level_to_rhs[coarsest_grid_ix]);
 
    // At each level of the grid hierarchy, coarsest-to-finest:
    for (int level = coarsest_grid_ix - 1; level >= 0; --level) {
      //  1. Update the current solution with the interpolated solution from the
      // coarser level: ui=ui+Piui+1
      level_to_soln[level] =
          level_to_soln[level] +
          interpolator->prolongation(level_to_soln[level + 1], level);
 
      //  2. Apply a couple of smoothing iterations (post-relaxation) to the
      // updated solution: ui=Si(Ai,fi,ui)
      smoother->smooth(level_to_coefficient_matrix[level], level_to_soln[level],
                       level_to_rhs[level]);
    }
 
    return;
  }
 
  const Eigen::Matrix<EleType, -1, 1>& solve() {
    size_t iter = 0;
    EleType error = 100;
    auto A = get_coefficient_matrix(finest_grid_ix);
    auto b = get_rhs(finest_grid_ix);
    while (iter < n_iters && error > tolerance) {
      vcycle();
      iter += 1;
      if ((iter % compute_error_every_n_iters) == 0) {
        auto u = get_soln(finest_grid_ix);
        error = rss(A, u, b);
 
        if (display_error) {
          std::cout << "Iter: " << iter << " | Error: " << error << std::endl;
        }
      }
    }
 
    if (error <= tolerance)
      std::cout << "AMG converged after " << iter << " iterations."
                << std::endl;
    else
      std::cout << "AMG did not converge after " << iter << " iterations."
                << std::endl;
 
    return level_to_soln[finest_grid_ix];
  }
 
  const Eigen::SparseMatrix<EleType>& get_coefficient_matrix(
      size_t level) const {
    return level_to_coefficient_matrix[level];
  }
 
  const Eigen::Matrix<EleType, -1, 1>& get_soln(size_t level) const {
    return level_to_soln[level];
  }
 
  const Eigen::Matrix<EleType, -1, 1>& get_rhs(size_t level) const {
    return level_to_rhs[level];
  }
 
  const size_t get_n_dofs(size_t level) const { return level_to_n_dofs[level]; }
 
  const EleType get_tolerance() const { return tolerance; }
 
  void display_error_on() {
    display_error = true;
    return;
  }
 
  void display_error_off() {
    display_error = true;
    return;
  }
};
 
}  // end namespace AMG