matrix_free_solver.h

Go to the documentation of this file.

 
#ifndef HYBRIDADRSOLVER_MATRIX_FREE_SOLVER_H
#define HYBRIDADRSOLVER_MATRIX_FREE_SOLVER_H
 
#include "adr_operator.h"
#include "core/problem_definition.h"
#include "core/solver.h"
#include "core/types.h"
 
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/numerics/vector_tools.h>
 
namespace HybridADRSolver {
 
using namespace dealii;

template <int dim, int fe_degree>
class MatrixFreeSolver : public ParallelSolverBase<dim> {
public:
    using Number = double;
    using LevelNumber = float; // Use float for level matrices (memory efficiency)
    using VectorType = LinearAlgebra::distributed::Vector<Number>;
    using LevelVectorType = LinearAlgebra::distributed::Vector<LevelNumber>;
 
    // Type aliases for multigrid components
    using SystemMatrixType = ADROperator<dim, fe_degree, Number>;
    using LevelMatrixType = ADROperator<dim, fe_degree, LevelNumber>;
    using SmootherType =
        PreconditionChebyshev<LevelMatrixType, LevelVectorType>;
    using SmootherPreconditionerType =
        DiagonalMatrix<LevelVectorType>;

    MatrixFreeSolver(const ProblemInterface<dim>& problem, MPI_Comm comm,
                     const SolverParameters& params = SolverParameters());

    void run(unsigned int n_refinements) override;

    SolverType get_solver_type() const override {
        return SolverType::MatrixFree;
    }

    std::string get_name() const override {
        return "Matrix-Free Solver (GMG + Hybrid MPI+Threading)";
    }

    const SystemMatrixType& get_system_operator() const {
        return system_operator;
    }

    const VectorType& get_solution() const { return solution; }

    double compute_memory_usage() const;

    double compute_l2_error() const;
 
protected:
    void setup_dofs() override;
    void assemble_system() override;
    void solve() override;
    void output_results(unsigned int cycle) const override;

    void setup_matrix_free();

    void setup_multigrid();

    void assemble_rhs();

    void solve_gmres_gmg();

    void solve_cg_gmg();

    void solve_gmres_jacobi();

    void solve_cg_jacobi();
 
private:
    const ProblemInterface<dim>& problem;
 
    // MatrixFree data structure for finest level
    std::shared_ptr<MatrixFree<dim, Number>> matrix_free_data;
 
    // System operator (matrix-free) for finest level
    SystemMatrixType system_operator;
 
    // Vectors
    VectorType solution;
    VectorType system_rhs;
 
    // Multigrid components
    MGConstrainedDoFs mg_constrained_dofs;
    MGLevelObject<LevelMatrixType> mg_matrices;
    MGLevelObject<std::shared_ptr<MatrixFree<dim, LevelNumber>>>
        mg_matrix_free_storage;
};
 
// Explicit instantiation declarations
extern template class MatrixFreeSolver<2, 1>;
extern template class MatrixFreeSolver<2, 2>;
extern template class MatrixFreeSolver<2, 3>;
extern template class MatrixFreeSolver<3, 1>;
extern template class MatrixFreeSolver<3, 2>;
extern template class MatrixFreeSolver<3, 3>;
 
} // namespace HybridADRSolver
 
#endif // HYBRIDADRSOLVER_MATRIX_FREE_SOLVER_H