DTR_mf.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/timer.h>
 
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/precondition.h>
 
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q1.h>
 
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_in.h>
 
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_matrix.h>
 
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
 
// Include the matrix-free headers
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/operators.h>
#include <deal.II/matrix_free/fe_evaluation.h>
 
#include <filesystem>
#include <iostream>
#include <fstream>
 
#include "problem_data.hpp"
 
const std::string output_dir = "./output_mf/";
 
namespace DTR_mf
{
  using namespace dealii;
 
  // To be efficient matrix-free implementation require knowledge of loop lengths at compile time
  const unsigned int dim = 2;
  
  template <int dim, int fe_degree, typename number>
  class DTROperation
      : public MatrixFreeOperators::
            Base<dim, LinearAlgebra::distributed::Vector<number>>
  {
  public:
    using value_type = number;
 
    DTROperation();
 
    void clear() override;
 
    void evaluate_coefficients(const problem_data::DiffusionCoefficient<dim> &diffusion_function,
                               const problem_data::TransportCoefficient<dim> &transport_function,
                               const problem_data::ReactionCoefficient<dim> &reaction_function,
                               const problem_data::ForcingTerm<dim> &forcing_term_function);
 
    virtual void compute_diagonal() override;
 
    const Table<2, VectorizedArray<number>> &get_diffusion_coefficient() const
    {
      return diffusion_coefficient;
    }
 
    const Table<2, Tensor<1, dim, VectorizedArray<number>>> &get_transport_coefficient() const
    {
      return transport_coefficient;
    }
 
    const Table<2, VectorizedArray<number>> &get_reaction_coefficient() const
    {
      return reaction_coefficient;
    }
 
    const Table<2, VectorizedArray<number>> &get_forcing_term_coefficient() const
    {
      return forcing_term_coefficient;
    }
 
  private:
    virtual void apply_add(
        LinearAlgebra::distributed::Vector<number> &dst,
        const LinearAlgebra::distributed::Vector<number> &src) const override;
 
    void
    local_apply(const MatrixFree<dim, number> &data,
                LinearAlgebra::distributed::Vector<number> &dst,
                const LinearAlgebra::distributed::Vector<number> &src,
                const std::pair<unsigned int, unsigned int> &cell_range) const;
 
    void local_compute_diagonal(
        const MatrixFree<dim, number> &data,
        LinearAlgebra::distributed::Vector<number> &dst,
        const unsigned int &dummy,
        const std::pair<unsigned int, unsigned int> &cell_range) const;
 
    Table<2, VectorizedArray<number>> diffusion_coefficient;
    Table<2, Tensor<1, dim, VectorizedArray<number>>> transport_coefficient;
    Table<2, VectorizedArray<number>> reaction_coefficient;
    Table<2, VectorizedArray<number>> forcing_term_coefficient;
  };
 
  template <int dim, int degree_finite_element = 2>
  class DTRProblem
  {
  public:
    DTRProblem(bool verbose = true);
 
    DTRProblem(std::ofstream& dimension_time_file, bool verbose = true);
 
    void run(unsigned int n_initial_refinements = 3, unsigned int n_cycles = 9);
 
    double compute_error(const VectorTools::NormType &norm_type) const;
 
    unsigned int get_cells() const { return triangulation.n_global_active_cells(); }
 
    unsigned int get_dofs() const { return dof_handler.n_dofs(); }
 
    unsigned int get_fe_degree() const { return degree_finite_element; }
 
  private:
    void setup_system();
 
    void assemble_rhs();
 
    void solve();
 
    void output_results(const unsigned int cycle) const;
 
#ifdef DEAL_II_WITH_P4EST
    parallel::distributed::Triangulation<dim> triangulation;
#else
    Triangulation<dim> triangulation;
#endif
 
    FE_Q<dim> fe;
 
    DoFHandler<dim> dof_handler;
 
    MappingQ1<dim> mapping;
 
    AffineConstraints<double> constraints;
 
    using SystemMatrixType =
        DTROperation<dim, degree_finite_element, double>;
 
    SystemMatrixType system_matrix;
 
    MGConstrainedDoFs mg_constrained_dofs;
 
    using LevelMatrixType = DTROperation<dim, degree_finite_element, float>;
 
    MGLevelObject<LevelMatrixType> mg_matrices;
 
    LinearAlgebra::distributed::Vector<double> solution;
 
    LinearAlgebra::distributed::Vector<double> lifting;
 
    LinearAlgebra::distributed::Vector<double> system_rhs;
 
    double setup_time;
 
    ConditionalOStream pcout;
 
    ConditionalOStream time_details;
 
    problem_data::DirichletBC1<dim> dirichletBC1;
 
    problem_data::DirichletBC2<dim> dirichletBC2;
 
    problem_data::NeumannBC1<dim> neumannBC1;
 
    problem_data::NeumannBC2<dim> neumannBC2;
  };
}
 
#include "DTR_mf.cpp"