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));
}