TaggingInterface.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "DenseMatrix.h"
#include "MooseTypes.h"
#include "MultiMooseEnum.h"
#include "Assembly.h"
#include "MooseVariableFE.h"
#include "SystemBase.h"

#include "libmesh/dense_vector.h"
#include "metaphysicl/raw_type.h"

#include <vector>

// Forward declarations
class InputParameters;
class MooseObject;
class SubProblem;
class Assembly;

#ifdef MOOSE_KOKKOS_ENABLED
namespace Moose::Kokkos
{
class ResidualObject;
class LinearSystemContributionObject;
class FESystem;
}
#endif

template <typename T>
InputParameters validParams();

struct ADResidualsPacket
{
  const DenseVector<ADReal> & residuals;
  const std::vector<dof_id_type> & dof_indices;
  const Real scaling_factor;
};

class TaggingInterface
{
public:
  TaggingInterface(const MooseObject * moose_object);

#ifdef MOOSE_KOKKOS_ENABLED

  TaggingInterface(const TaggingInterface & object, const Moose::Kokkos::FunctorCopy & key);
#endif

  virtual ~TaggingInterface();

  static InputParameters validParams();

  class VectorTagsKey
  {
    friend class AttribVectorTags;
    friend class NonlinearEigenSystem;
    friend class LinearSystemContributionObject;
    template <typename>
    friend class MooseObjectTagWarehouse;
#ifdef MOOSE_KOKKOS_ENABLED
    friend class Moose::Kokkos::ResidualObject;
    friend class Moose::Kokkos::FESystem;
    friend class Moose::Kokkos::LinearSystemContributionObject;
#endif

    VectorTagsKey() {}
    VectorTagsKey(const VectorTagsKey &) {}
  };

  class MatrixTagsKey
  {
    friend class AttribMatrixTags;
    friend class NonlinearEigenSystem;
    friend class LinearSystemContributionObject;
    template <typename>
    friend class MooseObjectTagWarehouse;
#ifdef MOOSE_KOKKOS_ENABLED
    friend class Moose::Kokkos::ResidualObject;
    friend class Moose::Kokkos::FESystem;
    friend class Moose::Kokkos::LinearSystemContributionObject;
#endif

    MatrixTagsKey() {}
    MatrixTagsKey(const MatrixTagsKey &) {}
  };

  enum class ResidualTagType
  {
    NonReference,
    Reference
  };

  void useVectorTag(const TagName & tag_name, VectorTagsKey);

  void useMatrixTag(const TagName & tag_name, MatrixTagsKey);

  void useVectorTag(TagID tag_id, VectorTagsKey);

  void useMatrixTag(TagID tag_id, MatrixTagsKey);

  bool isVectorTagged() { return _vector_tags.size() > 0; }

  bool isMatrixTagged() { return _matrix_tags.size() > 0; }

  bool hasVectorTags() const { return !_vector_tags.empty(); }

  const std::set<TagID> & getVectorTags(VectorTagsKey) const { return _vector_tags; }

  const std::set<TagID> & getMatrixTags(MatrixTagsKey) const { return _matrix_tags; }

protected:
  void prepareVectorTag(Assembly & assembly, unsigned int ivar);

  void prepareVectorTag(Assembly & assembly, unsigned int ivar, ResidualTagType tag_type);

  void prepareVectorTagNeighbor(Assembly & assembly, unsigned int ivar);

  void prepareVectorTagLower(Assembly & assembly, unsigned int ivar);

  void prepareMatrixTag(Assembly & assembly, unsigned int ivar, unsigned int jvar);

  void prepareMatrixTag(Assembly & assembly,
                        unsigned int ivar,
                        unsigned int jvar,
                        DenseMatrix<Number> & k) const;

  void prepareMatrixTagNonlocal(Assembly & assembly, unsigned int ivar, unsigned int jvar);

  void prepareMatrixTagNeighbor(Assembly & assembly,
                                unsigned int ivar,
                                unsigned int jvar,
                                Moose::DGJacobianType type);

  void prepareMatrixTagNeighbor(Assembly & assembly,
                                unsigned int ivar,
                                unsigned int jvar,
                                Moose::DGJacobianType type,
                                DenseMatrix<Number> & k) const;

  void prepareMatrixTagLower(Assembly & assembly,
                             unsigned int ivar,
                             unsigned int jvar,
                             Moose::ConstraintJacobianType type);

  void accumulateTaggedLocalResidual();

  void assignTaggedLocalResidual();

  void accumulateTaggedLocalMatrix();

  void accumulateTaggedLocalMatrix(Assembly & assembly,
                                   unsigned int ivar,
                                   unsigned int jvar,
                                   const DenseMatrix<Number> & k);

  void accumulateTaggedLocalMatrix(Assembly & assembly,
                                   unsigned int ivar,
                                   unsigned int jvar,
                                   Moose::DGJacobianType type,
                                   const DenseMatrix<Number> & k);

  void accumulateTaggedNonlocalMatrix();

  void assignTaggedLocalMatrix();

  template <typename Residuals, typename Indices>
  void addResiduals(Assembly & assembly,
                    const Residuals & residuals,
                    const Indices & dof_indices,
                    Real scaling_factor);

  template <typename Residuals, typename Indices>
  void addResiduals(Assembly & assembly,
                    const Residuals & residuals,
                    const Indices & dof_indices,
                    const std::vector<Real> & scaling_factors);

  template <typename T, typename Indices>
  void addResiduals(Assembly & assembly,
                    const DenseVector<T> & residuals,
                    const Indices & dof_indices,
                    Real scaling_factor);

  template <typename Residuals, typename Indices>
  void addResidualsAndJacobian(Assembly & assembly,
                               const Residuals & residuals,
                               const Indices & dof_indices,
                               Real scaling_factor);

  template <typename Residuals, typename Indices>
  void addJacobian(Assembly & assembly,
                   const Residuals & residuals,
                   const Indices & dof_indices,
                   Real scaling_factor);

  template <typename Residuals, typename Indices>
  void addJacobian(Assembly & assembly,
                   const Residuals & residuals,
                   const Indices & dof_indices,
                   const std::vector<Real> & scaling_factors);

  void addResiduals(Assembly & assembly, const ADResidualsPacket & packet);

  void addResidualsAndJacobian(Assembly & assembly, const ADResidualsPacket & packet);

  void addJacobian(Assembly & assembly, const ADResidualsPacket & packet);

  template <typename Residuals, typename Indices>
  void addResidualsWithoutConstraints(Assembly & assembly,
                                      const Residuals & residuals,
                                      const Indices & dof_indices,
                                      Real scaling_factor);

  template <typename Residuals, typename Indices>
  void addResidualsAndJacobianWithoutConstraints(Assembly & assembly,
                                                 const Residuals & residuals,
                                                 const Indices & dof_indices,
                                                 Real scaling_factor);

  template <typename Residuals, typename Indices>
  void addJacobianWithoutConstraints(Assembly & assembly,
                                     const Residuals & residuals,
                                     const Indices & dof_indices,
                                     Real scaling_factor);

  void addJacobianElement(Assembly & assembly,
                          Real value,
                          dof_id_type row_index,
                          dof_id_type column_index,
                          Real scaling_factor);

  void addJacobian(Assembly & assembly,
                   const DenseMatrix<Real> & local_k,
                   const std::vector<dof_id_type> & row_indices,
                   const std::vector<dof_id_type> & column_indices,
                   Real scaling_factor);

  template <typename T>
  void setResidual(SystemBase & sys, const T & residual, MooseVariableFE<T> & var);

  void setResidual(SystemBase & sys, Real residual, dof_id_type dof_index);

  template <typename SetResidualFunctor>
  void setResidual(SystemBase & sys, SetResidualFunctor set_residual_functor);

  SubProblem & _subproblem;

  DenseVector<Number> _local_re;

  DenseMatrix<Number> _local_ke;

  DenseMatrix<Number> _nonlocal_ke;

private:
  void prepareVectorTagInternal(Assembly & assembly,
                                unsigned int ivar,
                                const std::set<TagID> & vector_tags,
                                const std::set<TagID> & absolute_value_vector_tags);

#ifndef NDEBUG

  void checkForNans() const;
#endif

  std::set<TagID> _vector_tags;

  std::set<TagID> _abs_vector_tags;

  std::set<TagID> _matrix_tags;

  std::set<TagID> _non_ref_vector_tags;

  std::set<TagID> _non_ref_abs_vector_tags;

  std::set<TagID> _ref_vector_tags;

  std::set<TagID> _ref_abs_vector_tags;

  const MooseObject & _moose_object;

  const InputParameters & _tag_params;

  std::vector<DenseVector<Number> *> _re_blocks;

  std::vector<DenseVector<Number> *> _absre_blocks;

  std::vector<DenseMatrix<Number> *> _ke_blocks;

  std::vector<Real> _absolute_residuals;

  friend class NonlinearSystemBase;
};

#define usingTaggingInterfaceMembers                                                               \
  using TaggingInterface::_subproblem;                                                             \
  using TaggingInterface::accumulateTaggedLocalResidual;                                           \
  using TaggingInterface::accumulateTaggedLocalMatrix;                                             \
  using TaggingInterface::prepareVectorTag;                                                        \
  using TaggingInterface::prepareMatrixTag;                                                        \
  using TaggingInterface::prepareVectorTagNeighbor;                                                \
  using TaggingInterface::_local_re;                                                               \
  using TaggingInterface::prepareVectorTagLower;                                                   \
  using TaggingInterface::prepareMatrixTagNeighbor;                                                \
  using TaggingInterface::prepareMatrixTagLower;                                                   \
  using TaggingInterface::_local_ke

template <typename Residuals, typename Indices>
void
TaggingInterface::addResiduals(Assembly & assembly,
                               const Residuals & residuals,
                               const Indices & dof_indices,
                               const Real scaling_factor)
{
  assembly.cacheResiduals(
      residuals, dof_indices, scaling_factor, Assembly::LocalDataKey{}, _vector_tags);
  if (!_abs_vector_tags.empty())
  {
    _absolute_residuals.resize(residuals.size());
    for (const auto i : index_range(residuals))
      _absolute_residuals[i] = std::abs(MetaPhysicL::raw_value(residuals[i]));

    assembly.cacheResiduals(_absolute_residuals,
                            dof_indices,
                            scaling_factor,
                            Assembly::LocalDataKey{},
                            _abs_vector_tags);
  }
}

template <typename Residuals, typename Indices>
void
TaggingInterface::addResiduals(Assembly & assembly,
                               const Residuals & residuals,
                               const Indices & dof_indices,
                               const std::vector<Real> & scaling_factors)
{
  const auto count = scaling_factors.size();
  mooseAssert(dof_indices.size() % count == 0,
              "The number of dof indices should be divided cleanly by the variable count");
  const auto nshapes = dof_indices.size() / count;

  for (const auto j : make_range(count))
    // The Residuals type may not offer operator[] (e.g. eigen vectors) but more commonly it
    // should offer data()
    addResiduals(assembly,
                 Moose::makeSpan(residuals, j * nshapes, nshapes),
                 Moose::makeSpan(dof_indices, j * nshapes, nshapes),
                 scaling_factors[j]);
}

template <typename T, typename Indices>
void
TaggingInterface::addResiduals(Assembly & assembly,
                               const DenseVector<T> & residuals,
                               const Indices & dof_indices,
                               const Real scaling_factor)
{
  addResiduals(assembly, residuals.get_values(), dof_indices, scaling_factor);
}

template <typename Residuals, typename Indices>
void
TaggingInterface::addResidualsWithoutConstraints(Assembly & assembly,
                                                 const Residuals & residuals,
                                                 const Indices & dof_indices,
                                                 const Real scaling_factor)
{
  assembly.cacheResidualsWithoutConstraints(
      residuals, dof_indices, scaling_factor, Assembly::LocalDataKey{}, _vector_tags);
  if (!_abs_vector_tags.empty())
  {
    _absolute_residuals.resize(residuals.size());
    for (const auto i : index_range(residuals))
      _absolute_residuals[i] = std::abs(MetaPhysicL::raw_value(residuals[i]));

    assembly.cacheResidualsWithoutConstraints(_absolute_residuals,
                                              dof_indices,
                                              scaling_factor,
                                              Assembly::LocalDataKey{},
                                              _abs_vector_tags);
  }
}

template <typename Residuals, typename Indices>
void
TaggingInterface::addResidualsAndJacobian(Assembly & assembly,
                                          const Residuals & residuals,
                                          const Indices & dof_indices,
                                          Real scaling_factor)
{
  addResiduals(assembly, residuals, dof_indices, scaling_factor);
  addJacobian(assembly, residuals, dof_indices, scaling_factor);
}

template <typename Residuals, typename Indices>
void
TaggingInterface::addJacobian(Assembly & assembly,
                              const Residuals & residuals,
                              const Indices & dof_indices,
                              Real scaling_factor)
{
  assembly.cacheJacobian(
      residuals, dof_indices, scaling_factor, Assembly::LocalDataKey{}, _matrix_tags);
}

template <typename Residuals, typename Indices>
void
TaggingInterface::addJacobian(Assembly & assembly,
                              const Residuals & residuals,
                              const Indices & dof_indices,
                              const std::vector<Real> & scaling_factors)
{
  const auto count = scaling_factors.size();
  mooseAssert(dof_indices.size() % count == 0,
              "The number of dof indices should be divided cleanly by the variable count");
  const auto nshapes = dof_indices.size() / count;

  for (const auto j : make_range(count))
    // The Residuals type may not offer operator[] (e.g. eigen vectors) but more commonly it
    // should offer data()
    addJacobian(assembly,
                Moose::makeSpan(residuals, j * nshapes, nshapes),
                Moose::makeSpan(dof_indices, j * nshapes, nshapes),
                scaling_factors[j]);
}

template <typename Residuals, typename Indices>
void
TaggingInterface::addResidualsAndJacobianWithoutConstraints(Assembly & assembly,
                                                            const Residuals & residuals,
                                                            const Indices & dof_indices,
                                                            Real scaling_factor)
{
  addResidualsWithoutConstraints(assembly, residuals, dof_indices, scaling_factor);
  addJacobianWithoutConstraints(assembly, residuals, dof_indices, scaling_factor);
}

template <typename Residuals, typename Indices>
void
TaggingInterface::addJacobianWithoutConstraints(Assembly & assembly,
                                                const Residuals & residuals,
                                                const Indices & dof_indices,
                                                Real scaling_factor)
{
  assembly.cacheJacobianWithoutConstraints(
      residuals, dof_indices, scaling_factor, Assembly::LocalDataKey{}, _matrix_tags);
}

inline void
TaggingInterface::addJacobianElement(Assembly & assembly,
                                     const Real value,
                                     const dof_id_type row_index,
                                     const dof_id_type column_index,
                                     const Real scaling_factor)
{
  assembly.cacheJacobian(
      row_index, column_index, value * scaling_factor, Assembly::LocalDataKey{}, _matrix_tags);
}

inline void
TaggingInterface::addJacobian(Assembly & assembly,
                              const DenseMatrix<Real> & local_k,
                              const std::vector<dof_id_type> & row_indices,
                              const std::vector<dof_id_type> & column_indices,
                              const Real scaling_factor)
{
  assembly.cacheJacobianBlock(
      local_k, row_indices, column_indices, scaling_factor, Assembly::LocalDataKey{}, _matrix_tags);
}

template <typename T>
void
TaggingInterface::setResidual(SystemBase & sys, const T & residual, MooseVariableFE<T> & var)
{
  for (const auto tag_id : _vector_tags)
    if (sys.hasVector(tag_id))
      var.insertNodalValue(sys.getVector(tag_id), residual);
}

inline void
TaggingInterface::setResidual(SystemBase & sys, const Real residual, const dof_id_type dof_index)
{
  for (const auto tag_id : _vector_tags)
    if (sys.hasVector(tag_id))
      sys.getVector(tag_id).set(dof_index, residual);
}

template <typename SetResidualFunctor>
void
TaggingInterface::setResidual(SystemBase & sys, const SetResidualFunctor set_residual_functor)
{
  for (const auto tag_id : _vector_tags)
    if (sys.hasVector(tag_id))
      set_residual_functor(sys.getVector(tag_id));
}