matrix_free.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2020 by the hyper.deal authors
//
// This file is part of the hyper.deal library.
//
// The hyper.deal library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 3.0 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.MD at
// the top level directory of hyper.deal.
//
// ---------------------------------------------------------------------
 
#ifndef HYPERDEAL_NDIM_MATRIXFREE
#define HYPERDEAL_NDIM_MATRIXFREE
 
#include <hyper.deal/base/config.h>
 
#include <deal.II/lac/la_parallel_vector.h>
 
#include <deal.II/matrix_free/matrix_free.h>
 
#include <hyper.deal/base/memory_consumption.h>
#include <hyper.deal/base/timers.h>
#include <hyper.deal/matrix_free/dof_info.h>
#include <hyper.deal/matrix_free/face_info.h>
#include <hyper.deal/matrix_free/id.h>
#include <hyper.deal/matrix_free/shape_info.h>
#include <hyper.deal/matrix_free/vector_partitioner.h>
 
namespace hyperdeal
{
  template <int dim_x, int dim_v, typename Number, typename VectorizedArrayType>
  class MatrixFree
  {
  public:
    using ID = TensorID;
 
    using VectorizedArrayTypeX = VectorizedArrayType;
    using VectorizedArrayTypeV = dealii::VectorizedArray<Number, 1>;
 
    enum class DataAccessOnFaces
    {
      none,
      values,
      gradients,
      unspecified
    };
 
    struct AdditionalData
    {
      AdditionalData()
        : do_ghost_faces(true)
        , do_buffering(false)
        , use_ecl(true)
        , overlapping_level(0 /*no overlapping communication-computation*/)
      {}
 
      bool do_ghost_faces;
 
      bool do_buffering;
 
      bool use_ecl;
 
      unsigned int overlapping_level;
    };
 
 
    MatrixFree(const MPI_Comm comm,
               const MPI_Comm comm_sm,
               const dealii::MatrixFree<dim_x, Number, VectorizedArrayTypeX>
                 &matrix_free_x,
               const dealii::MatrixFree<dim_v, Number, VectorizedArrayTypeV>
                 &matrix_free_v);
 
    void
    reinit(const AdditionalData &ad = AdditionalData());
 
    template <typename OutVector, typename InVector>
    void
    cell_loop(
      const std::function<
        void(const MatrixFree &, OutVector &, const InVector &, const ID)>
        &             cell_operation,
      OutVector &     dst,
      const InVector &src) const;
 
    template <typename CLASS, typename OutVector, typename InVector>
    void
    cell_loop(void (CLASS::*cell_operation)(const MatrixFree &,
                                            OutVector &,
                                            const InVector &,
                                            const ID),
              CLASS *         owning_class,
              OutVector &     dst,
              const InVector &src) const;
 
    template <typename OutVector, typename InVector>
    void
    loop_cell_centric(
      const std::function<
        void(const MatrixFree &, OutVector &, const InVector &, const ID)>
        &                     cell_operation,
      OutVector &             dst,
      const InVector &        src,
      const DataAccessOnFaces src_vector_face_access =
        DataAccessOnFaces::unspecified,
      Timers *timers = nullptr) const;
 
    template <typename CLASS, typename OutVector, typename InVector>
    void
    loop_cell_centric(void (CLASS::*cell_operation)(const MatrixFree &,
                                                    OutVector &,
                                                    const InVector &,
                                                    const ID),
                      CLASS *                 owning_class,
                      OutVector &             dst,
                      const InVector &        src,
                      const DataAccessOnFaces src_vector_face_access =
                        DataAccessOnFaces::unspecified,
                      Timers *timers = nullptr) const;
 
    template <typename OutVector, typename InVector>
    void
    loop(const std::function<
           void(const MatrixFree &, OutVector &, const InVector &, const ID)>
           &cell_operation,
         const std::function<
           void(const MatrixFree &, OutVector &, const InVector &, const ID)>
           &face_operation,
         const std::function<
           void(const MatrixFree &, OutVector &, const InVector &, const ID)>
           &                     boundary_operation,
         OutVector &             dst,
         const InVector &        src,
         const DataAccessOnFaces dst_vector_face_access =
           DataAccessOnFaces::unspecified,
         const DataAccessOnFaces src_vector_face_access =
           DataAccessOnFaces::unspecified,
         Timers *timers = nullptr) const;
 
    template <typename CLASS, typename OutVector, typename InVector>
    void
    loop(void (CLASS::*cell_operation)(const MatrixFree &,
                                       OutVector &,
                                       const InVector &,
                                       const ID),
         void (CLASS::*face_operation)(const MatrixFree &,
                                       OutVector &,
                                       const InVector &,
                                       const ID),
         void (CLASS::*boundary_operation)(const MatrixFree &,
                                           OutVector &,
                                           const InVector &,
                                           const ID),
         CLASS *                 owning_class,
         OutVector &             dst,
         const InVector &        src,
         const DataAccessOnFaces dst_vector_face_access =
           DataAccessOnFaces::unspecified,
         const DataAccessOnFaces src_vector_face_access =
           DataAccessOnFaces::unspecified,
         Timers *timers = nullptr) const;
 
    const MPI_Comm &
    get_communicator() const;
 
    void
    initialize_dof_vector(
      dealii::LinearAlgebra::distributed::Vector<Number> &vec,
      const unsigned int                                  dof_handler_index = 0,
      const bool                                          do_ghosts = true,
      const bool zero_out_values = true) const;
 
    dealii::types::boundary_id
    get_boundary_id(const ID macro_face) const;
 
    dealii::types::boundary_id
    get_faces_by_cells_boundary_id(const TensorID &   macro_cell,
                                   const unsigned int face_number) const;
 
    bool
    is_ecl_supported() const;
 
    bool
    are_ghost_faces_supported() const;
 
    const dealii::MatrixFree<dim_x, Number, VectorizedArrayTypeX> &
    get_matrix_free_x() const;
 
    const dealii::MatrixFree<dim_v, Number, VectorizedArrayTypeV> &
    get_matrix_free_v() const;
 
    const internal::MatrixFreeFunctions::DoFInfo &
    get_dof_info() const;
 
    const internal::MatrixFreeFunctions::FaceInfo &
    get_face_info() const;
 
    const internal::MatrixFreeFunctions::ShapeInfo<Number> &
    get_shape_info() const;
 
    MemoryConsumption
    memory_consumption() const;
 
    const std::shared_ptr<
      const dealii::internal::MatrixFreeFunctions::VectorDataExchange::Base> &
    get_vector_partitioner() const;
 
 
  private:
    const MPI_Comm comm;
 
    const MPI_Comm comm_sm;
 
    const dealii::MatrixFree<dim_x, Number, VectorizedArrayTypeX>
      &matrix_free_x;
 
    const dealii::MatrixFree<dim_v, Number, VectorizedArrayTypeV>
      &matrix_free_v;
 
    bool do_buffering;
 
    bool do_ghost_faces;
 
    bool use_ecl;
 
    std::shared_ptr<
      const dealii::internal::MatrixFreeFunctions::VectorDataExchange::Base>
      partitioner;
 
    internal::MatrixFreeFunctions::DoFInfo dof_info;
 
    internal::MatrixFreeFunctions::FaceInfo face_info;
 
    internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info;
 
    std::array<std::vector<ID>, 3> partitions;
  };
 
} // namespace hyperdeal
 
#include <hyper.deal/matrix_free/matrix_free.templates.h>
 
#endif