mf_operator.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/base/vectorization.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/operators.h>
 
#include <prismspf/core/field_attributes.h>
#include <prismspf/core/field_container.h>
#include <prismspf/core/group_solution_handler.h>
#include <prismspf/core/pde_operator_base.h>
#include <prismspf/core/simulation_timer.h>
#include <prismspf/core/solution_indexer.h>
#include <prismspf/core/solve_block.h>
#include <prismspf/core/types.h>
 
#include <prismspf/solvers/solver_base.h>
 
#include <prismspf/config.h>
 
#include <vector>
 
#if DEAL_II_VERSION_MAJOR >= 9 && DEAL_II_VERSION_MINOR >= 7
#  include <deal.II/base/enable_observer_pointer.h>
#  define MATRIX_FREE_OPERATOR_BASE dealii::EnableObserverPointer
#else
#  include <deal.II/base/subscriptor.h>
#  define MATRIX_FREE_OPERATOR_BASE dealii::Subscriptor
#endif
 
PRISMS_PF_BEGIN_NAMESPACE

template <unsigned int dim, unsigned int degree, typename number>
class MFOperator : public MATRIX_FREE_OPERATOR_BASE
{
public:
  using ScalarValue = dealii::VectorizedArray<number>;
  using VectorValue = dealii::Tensor<1, dim, ScalarValue>;
 
  using Operator = void (PDEOperatorBase<dim, degree, number>::*)(
    FieldContainer<dim, degree, number> &, /* variable_list */
    const SimulationTimer &,               /* sim_timer */
    unsigned int                           /* solve_block_id */
  ) const;
 
  template <TensorRank Rank>
  using Value = std::conditional_t<Rank == TensorRank::Scalar,
                                   ScalarValue,
                                   dealii::Tensor<int(Rank), dim, ScalarValue>>;
 
  template <TensorRank Rank>
  Value<Rank>
  identity()
  {
    if constexpr (Rank == TensorRank::Scalar)
      {
        return ScalarValue(1.0);
      }
    else
      {
        static Value<Rank> ident = []()
        {
          Value<Rank> obj;
          for (int i = 0; i < Value<Rank>::n_independent_components; ++i)
            {
              obj[Value<Rank>::unrolled_to_component_indices(i)] = 1.0;
            }
        }();
        return ident;
      }
  }
 
  template <TensorRank Rank>
  Value<Rank>
  zero()
  {
    if constexpr (Rank == TensorRank::Scalar)
      {
        return ScalarValue(0.0);
      }
    else
      {
        return Value<Rank>();
      }
  }

  explicit MFOperator(const PDEOperatorBase<dim, degree, number> &operator_owner,
                      Operator                                    oper,
                      const std::vector<FieldAttributes>         &_field_attributes,
                      const SolutionIndexer<dim, number>         &_solution_indexer,
                      unsigned int                                _relative_level,
                      DependencyMap                               _dependency_map,
                      const SimulationTimer                      &_sim_timer)
    : MATRIX_FREE_OPERATOR_BASE()
    , pde_operator(&operator_owner)
    , pde_op(oper)
    , field_attributes(_field_attributes)
    , solution_indexer(&_solution_indexer)
    , relative_level(_relative_level)
    , dependency_map(std::move(_dependency_map))
    , sim_timer(&_sim_timer)
  {}

  void
  initialize(const GroupSolutionHandler<dim, number> &dst_solution)
  {
    solve_block          = dst_solution.get_solve_block();
    data                 = &(dst_solution.get_matrix_free(relative_level));
    field_to_block_index = dst_solution.get_global_to_block_index();
    for (unsigned int field_index : solve_block.field_indices)
      {
        dependency_map[field_index]; // creates entry if not already present
      }
  }
 
  // public:
  void
  compute_operator(BlockVector<number>       &dst,
                   const BlockVector<number> &src = BlockVector<number>()) const;
 
private:
  void
  compute_local_operator(const MatrixFree<dim, number>               &_data,
                         BlockVector<number>                         &dst,
                         const BlockVector<number>                   &src,
                         const std::pair<unsigned int, unsigned int> &cell_range) const;
 
  // public:
  //   /**
  //    * @brief Compute the diagonal of this operator.
  //    */
  //   void
  //   compute_diagonal();
  //
  // private:
  //   /**
  //    * @brief Local computation of the diagonal of the operator.
  //    */
  //   void
  //   compute_local_diagonal(const MatrixFree<dim, number> &_data,
  //                          BlockVector<number>                                 &dst,
  //                          const unsigned int                          &dummy,
  //                          const std::pair<unsigned int, unsigned int> &cell_range)
  //                          const;
  //
  //   template <TensorRank Rank>
  //   dealii::AlignedVector<Value<Rank>>
  //   compute_field_diagonal(FieldContainer<dim, degree, number> &variable_list,
  //                          DSTContainer<dim, degree, number>   &dst_fields,
  //                          unsigned int                         field_index) const;
 
public:
  void
  set_scaling_diagonal(bool                                               scale,
                       const std::vector<const SolutionVector<number> *> &diagonal)
  {
    scaling_diagonal  = diagonal;
    scale_by_diagonal = scale;
  }

  dealii::types::global_dof_index
  m() const;

  number
  el(const unsigned int &row, const unsigned int &col) const;

  void
  clear();

  // void
  // set_constrained_entries_to_one(SolutionVector<number> &dst) const;

  const MatrixFree<dim, number> *
  get_matrix_free() const;

  const std::shared_ptr<dealii::DiagonalMatrix<SolutionVector<number>>> &
  get_matrix_diagonal_inverse() const;

  void
  vmult(BlockVector<number> &dst, const BlockVector<number> &src) const;
 
  // NOLINTBEGIN(readability-identifier-naming)

  void
  Tvmult(BlockVector<number> &dst, const BlockVector<number> &src) const;
 
  // NOLINTEND(readability-identifier-naming)

  bool read_plain = false;
 
private:
  std::vector<FieldAttributes> field_attributes;

  SolveBlock solve_block;

  const PDEOperatorBase<dim, degree, number> *pde_operator;
  Operator pde_op;

  const SolutionIndexer<dim, number> *solution_indexer;
  unsigned int relative_level;

  DependencyMap dependency_map;

  const SimulationTimer *sim_timer;

  std::vector<const SolutionVector<number> *> scaling_diagonal;

  bool scale_by_diagonal = false;

  std::vector<unsigned int> field_to_block_index;

  const MatrixFree<dim, number> *data;

  std::vector<std::vector<unsigned int>> edge_constrained_indices;

  std::shared_ptr<dealii::DiagonalMatrix<SolutionVector<number>>> diagonal_entries;

  std::shared_ptr<dealii::DiagonalMatrix<SolutionVector<number>>>
    inverse_diagonal_entries;
};
 
PRISMS_PF_END_NAMESPACE