preconditioner_trilinos_wrapper.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <deal.II/base/exceptions.h>
 
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
 
#include <meltpooldg/core/scratch_data.hpp>
 
#include <any>
#include <type_traits>
 
namespace MeltPoolDG
{
  template <typename OperatorType, typename VectorType>
  concept DealiipreconditionerWrapperOperatorType =
    requires(const OperatorType op, dealii::TrilinosWrappers::SparseMatrix &matrix) {
      op.compute_system_matrix_from_matrixfree(matrix);
    };
 
  template <int dim,
            typename number,
            typename VectorType,
            typename DealiiPreconditionerType,
            typename OperatorType>
  class DealiiPreconditionerWrapper
  {
  public:
    explicit DealiiPreconditionerWrapper(const OperatorType                  *operator_in,
                                         const ScratchData<dim, dim, number> &scratch_data_in,
                                         const unsigned                       dof_idx_in,
                                         const bool                           do_matrix_free = true)
      : eq_operator(operator_in)
      , scratch_data(scratch_data_in)
      , dof_idx(dof_idx_in)
      , do_matrix_free(do_matrix_free)
    {}
 
    void
    vmult(VectorType &dst, const VectorType &src) const
    {
      // deal.II trilinos preconditioner do not support vmult(BlockVectorType& BlockVectorType&).
      // This is done manually if required.
      // TODO: Can this be done more elegant?
      if constexpr (std::is_same_v<dealii::LinearAlgebra::distributed::BlockVector<
                                     typename VectorType::value_type>,
                                   VectorType>)
        for (unsigned int b = 0; b < dst.n_blocks(); ++b)
          preconditioner.vmult(dst.block(b), src.block(b));
      else
        preconditioner.vmult(dst, src);
    }
 
    void
    update(const std::any &external_preconditioner_matrix = std::any())
    {
      dealii::TrilinosWrappers::SparseMatrix *preconditioner_matrix_ptr;
      if (not external_preconditioner_matrix.has_value())
        {
          if constexpr (DealiipreconditionerWrapperOperatorType<OperatorType, VectorType>)
            {
              Assert(eq_operator != nullptr,
                     dealii::ExcMessage("The provided operator has no valid address (nullptr)."));
              eq_operator->compute_system_matrix_from_matrixfree(preconditioner_system_matrix);
              preconditioner_matrix_ptr = &preconditioner_system_matrix;
            }
          else
            {
              AssertThrow(false,
                          dealii::ExcMessage(
                            "Neither a rule to compute the preconditioner matrix nor an external"
                            " preconditioner matrix are provided."));
            }
        }
      else
        {
          Assert(external_preconditioner_matrix.has_value(),
                 dealii::ExcMessage(
                   "No external system matrix has been passed to the preconditioner!"));
          preconditioner_matrix_ptr =
            std::any_cast<dealii::TrilinosWrappers::SparseMatrix *>(external_preconditioner_matrix);
        }
      preconditioner.initialize(*preconditioner_matrix_ptr,
                                typename DealiiPreconditionerType::AdditionalData());
    }
 
    void
    reinit()
    {
      if (eq_operator != nullptr && do_matrix_free)
        {
          dsp.reinit(scratch_data.get_dof_handler(dof_idx).n_dofs(),
                     scratch_data.get_dof_handler(dof_idx).n_dofs(),
                     scratch_data.get_locally_relevant_dofs(dof_idx));
          if (scratch_data.is_FE_DGQ(dof_idx) or
              not scratch_data.get_cut_type(dof_idx) == CutUtil::CutPhaseType::not_cut)
            {
              dealii::DoFTools::make_flux_sparsity_pattern(scratch_data.get_dof_handler(dof_idx),
                                                           dsp,
                                                           scratch_data.get_constraint(dof_idx),
                                                           false);
            }
          else
            {
              dealii::DoFTools::make_sparsity_pattern(scratch_data.get_dof_handler(dof_idx),
                                                      dsp,
                                                      scratch_data.get_constraint(dof_idx));
            }
          dealii::SparsityTools::distribute_sparsity_pattern(
            dsp,
            scratch_data.get_locally_owned_dofs(dof_idx),
            scratch_data.get_mpi_comm(),
            scratch_data.get_locally_relevant_dofs(dof_idx));
          preconditioner_system_matrix.reinit(scratch_data.get_locally_owned_dofs(dof_idx),
                                              scratch_data.get_locally_owned_dofs(dof_idx),
                                              dsp,
                                              scratch_data.get_mpi_comm());
        }
    }
 
  private:
    DealiiPreconditionerType               preconditioner;
    dealii::TrilinosWrappers::SparseMatrix preconditioner_system_matrix;
    dealii::DynamicSparsityPattern         dsp;
 
    const OperatorType *eq_operator;
 
    const ScratchData<dim, dim, number> &scratch_data;
 
    const unsigned int dof_idx;
 
    const bool do_matrix_free;
  };
 
  template <int dim, typename VectorType, typename number>
  class IdentityPreconditioner
  {
  public:
    void
    vmult(VectorType &dst, const VectorType &src) const
    {
      dst = src;
    }
 
    void
    update(const std::any &) const
    {}
 
    void
    reinit() const
    {}
  };
} // namespace MeltPoolDG