solver.h

Go to the documentation of this file.

#ifndef HYBRIDADRSOLVER_SOLVER_H
#define HYBRIDADRSOLVER_SOLVER_H
 
#include "types.h"
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/timer.h>
 
#include <deal.II/distributed/tria.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/lac/affine_constraints.h>
 
#include <utility>
 
namespace HybridADRSolver {
using namespace dealii;

template <int dim> class ParallelSolverBase {
public:
    ParallelSolverBase(MPI_Comm comm, SolverParameters params,
                       bool construct_mg_hierarchy = false);
    virtual ~ParallelSolverBase() = default;

    virtual void run(unsigned int n_refinements) = 0;

    const TimingResults& get_timing_results() const { return timing_results; }

    virtual SolverType get_solver_type() const = 0;

    virtual std::string get_name() const = 0;
 
protected:
    virtual void setup_grid(unsigned int n_refinements);

    virtual void setup_dofs() = 0;

    virtual void assemble_system() = 0;

    virtual void solve() = 0;

    virtual void output_results(unsigned int cycle) const = 0;
 
    // MPI communication
    MPI_Comm mpi_communicator;
    unsigned int n_mpi_processes;
    unsigned int this_mpi_process;
 
    // Mesh and finite element
    parallel::distributed::Triangulation<dim> triangulation;
    std::unique_ptr<FiniteElement<dim>> fe;
    DoFHandler<dim> dof_handler;
    std::unique_ptr<Mapping<dim>> mapping;
 
    // Constraints (Dirichlet BC, hanging nodes)
    AffineConstraints<double> constraints;
 
    // DoF distribution info
    IndexSet locally_owned_dofs;
    IndexSet locally_relevant_dofs;
 
    // Parameters
    SolverParameters parameters;
 
    // Timing
    TimingResults timing_results;
 
    // Output
    ConditionalOStream pcout;
    TimerOutput computing_timer;
 
    // Whether multigrid hierarchy was constructed
    bool has_mg_hierarchy;
 
private:
    static typename parallel::distributed::Triangulation<dim>::Settings
    get_triangulation_settings(bool construct_mg_hierarchy) {
        if (construct_mg_hierarchy) {
            return parallel::distributed::Triangulation<
                dim>::construct_multigrid_hierarchy;
        }
        return parallel::distributed::Triangulation<dim>::default_setting;
    }
};
 
template <int dim>
ParallelSolverBase<dim>::ParallelSolverBase(MPI_Comm comm,
                                            SolverParameters params,
                                            bool construct_mg_hierarchy)
    : mpi_communicator(comm),
      n_mpi_processes(Utilities::MPI::n_mpi_processes(comm)),
      this_mpi_process(Utilities::MPI::this_mpi_process(comm)),
      triangulation(comm,
                    typename Triangulation<dim>::MeshSmoothing(
                        Triangulation<dim>::limit_level_difference_at_vertices |
                        Triangulation<dim>::smoothing_on_refinement |
                        Triangulation<dim>::smoothing_on_coarsening),
                    get_triangulation_settings(construct_mg_hierarchy)),
      dof_handler(triangulation), parameters(std::move(params)),
      pcout(std::cout, this_mpi_process == 0),
      computing_timer(comm, pcout, TimerOutput::never, TimerOutput::wall_times),
      has_mg_hierarchy(construct_mg_hierarchy) {
    // Set thread limit
    if (parameters.n_threads != numbers::invalid_unsigned_int)
        MultithreadInfo::set_thread_limit(parameters.n_threads);
}
 
template <int dim>
void ParallelSolverBase<dim>::setup_grid(const unsigned int n_refinements) {
    // Create unit hypercube domain [0,1]^dim
    // The 'true' parameter sets boundary IDs:
    // - 2D: left=0, right=1, bottom=2, top=3
    // - 3D: left=0, right=1, front=2, back=3, bottom=4, top=5
    GridGenerator::hyper_cube(triangulation, 0.0, 1.0, true);
 
    // Refine globally
    triangulation.refine_global(n_refinements);
 
    if (parameters.verbose) {
        pcout << "   Active cells: " << triangulation.n_global_active_cells()
              << std::endl;
        if (has_mg_hierarchy) {
            pcout << "   MG levels:    " << triangulation.n_global_levels()
                  << std::endl;
        }
    }
}
 
} // namespace HybridADRSolver
 
#endif // HYBRIDADRSOLVER_SOLVER_H