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/utilities/utilities.h>
 
#include <prismspf/config.h>
 
#include "prismspf/core/matrix_free_manager.h"
 
#include <memory>
#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 || dim == 1,
                                   ScalarValue,
                                   dealii::Tensor<int(Rank), dim, ScalarValue>>;
 
  template <TensorRank Rank>
  static Value<Rank>
  identity()
  {
    if constexpr (Rank == TensorRank::Scalar || dim == 1)
      {
        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 obj;
        }();
        return ident;
      }
  }
 
  template <TensorRank Rank>
  static Value<Rank>
  zero()
  {
    if constexpr (Rank == TensorRank::Scalar || dim == 1)
      {
        return ScalarValue(0.0);
      }
    else
      {
        return Value<Rank>();
      }
  }

  MFOperator(const PDEOperatorBase<dim, degree, number> &operator_owner,
             Operator                                    oper,
             const std::vector<FieldAttributes>         &_field_attributes,
             const SolutionIndexer<dim, number>         &_solution_indexer,
             const MatrixFreeManager<dim, number>       &_matrix_free_manager,
             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)
    , matrix_free_manager(&_matrix_free_manager)
    , dependency_map(std::move(_dependency_map))
    , sim_timer(&_sim_timer)
  {
    set_relative_level(0);
  }

  void
  initialize(const GroupSolutionHandler<dim, number> &dst_solution)
  {
    solve_block          = dst_solution.get_solve_block();
    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
      }
  }
 
  void
  set_relative_level(unsigned int _relative_level)
  {
    relative_level = _relative_level;
    AssertThrow(matrix_free_manager != nullptr, dealii::ExcNotInitialized());
    data = &(matrix_free_manager->get_shared_matrix_free(relative_level));
  }
 
  // 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:
  void
  compute_diagonal(BlockVector<number> &dst, const BlockVector<number> &src) const;
 
private:
  void
  compute_local_diagonal(const MatrixFree<dim, number>               &_data,
                         BlockVector<number>                         &diagonal,
                         const BlockVector<number>                   &dummy_src,
                         const std::pair<unsigned int, unsigned int> &cell_range) const;
 
  template <TensorRank Rank>
  void
  compute_local_field_diagonal(FieldContainer<dim, degree, number> &variable_list,
                               BlockVector<number>                 &diagonal,
                               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();

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

  const std::shared_ptr<dealii::DiagonalMatrix<BlockVector<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)

  void
  reinit_matrix_diagonal(const BlockVector<number> &shape);

  void
  eval_matrix_diagonal();

  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 = nullptr;

  const MatrixFreeManager<dim, number> *matrix_free_manager = nullptr;

  const MatrixFree<dim, number> *data = nullptr;

  unsigned int relative_level = 0;

  DependencyMap dependency_map;

  const SimulationTimer *sim_timer = nullptr;

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

  bool scale_by_diagonal = false;

  std::vector<unsigned int> field_to_block_index;

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

  std::shared_ptr<dealii::DiagonalMatrix<BlockVector<number>>> diagonal_entries =
    std::make_shared<dealii::DiagonalMatrix<BlockVector<number>>>();

  std::shared_ptr<dealii::DiagonalMatrix<BlockVector<number>>> inverse_diagonal_entries =
    std::make_shared<dealii::DiagonalMatrix<BlockVector<number>>>();
};
 
PRISMS_PF_END_NAMESPACE