ADNodalBC.C

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
 
#include "ADNodalBC.h"
#include "NodalBCBase.h"
#include "ADDirichletBCBase.h"
 
// MOOSE includes
#include "Assembly.h"
#include "SubProblem.h"
#include "SystemBase.h"
#include "MooseVariableFE.h"
#include "MooseVariableScalar.h"
#include "FEProblemBase.h"
 
#include "libmesh/quadrature.h"
 
template <typename T, typename Base>
InputParameters
ADNodalBCTempl<T, Base>::validParams()
{
  return Base::validParams();
}
 
template <>
InputParameters
ADNodalBCTempl<RealVectorValue, NodalBCBase>::validParams()
{
  InputParameters params = NodalBCBase::validParams();
  // The below parameters are useful for vector Nodal BCs
  params.addParam<bool>("set_x_comp", true, "Whether to set the x-component of the variable");
  params.addParam<bool>("set_y_comp", true, "Whether to set the y-component of the variable");
  params.addParam<bool>("set_z_comp", true, "Whether to set the z-component of the variable");
  return params;
}
 
template <>
InputParameters
ADNodalBCTempl<RealVectorValue, ADDirichletBCBase>::validParams()
{
  InputParameters params = ADDirichletBCBase::validParams();
  // The below parameters are useful for vector Nodal BCs
  params.addParam<bool>("set_x_comp", true, "Whether to set the x-component of the variable");
  params.addParam<bool>("set_y_comp", true, "Whether to set the y-component of the variable");
  params.addParam<bool>("set_z_comp", true, "Whether to set the z-component of the variable");
  return params;
}
 
template <typename T, typename Base>
ADNodalBCTempl<T, Base>::ADNodalBCTempl(const InputParameters & parameters)
  : Base(parameters),
    MooseVariableInterface<T>(this,
                              true,
                              "variable",
                              Moose::VarKindType::VAR_SOLVER,
                              std::is_same_v<T, Real> ? Moose::VarFieldType::VAR_FIELD_STANDARD
                              : std::is_same_v<T, RealVectorValue>
                                  ? Moose::VarFieldType::VAR_FIELD_VECTOR
                                  : Moose::VarFieldType::VAR_FIELD_ARRAY),
    ADFunctorInterface(this),
    _var(*this->mooseVariable()),
    _current_node(_var.node()),
    _u(_var.adNodalValue()),
    _undisplaced_assembly(_fe_problem.assembly(_tid, _sys.number()))
{
  if constexpr (std::is_same_v<T, RealVectorValue>)
  {
    _set_components.resize(Moose::dim);
    _set_components[0] = this->template getParam<bool>("set_x_comp");
    _set_components[1] = this->template getParam<bool>("set_y_comp");
    _set_components[2] = this->template getParam<bool>("set_z_comp");
  }
  else
    _set_components.resize(_var.count(), true);
 
  _subproblem.haveADObjects(true);
 
  addMooseVariableDependency(this->mooseVariable());
}
 
namespace
{
template <typename T>
const T &
conversionHelper(const T & value, const unsigned int)
{
  return value;
}
 
template <typename T>
const T &
conversionHelper(const VectorValue<T> & value, const unsigned int i)
{
  return value(i);
}
 
template <typename T>
const T &
conversionHelper(const Eigen::Matrix<T, Eigen::Dynamic, 1> & value, const unsigned int i)
{
  return value(i);
}
}
 
template <typename T, typename Base>
void
ADNodalBCTempl<T, Base>::addResidual(const T & residual,
                                     const std::vector<dof_id_type> & dof_indices)
{
  mooseAssert(dof_indices.size() <= _set_components.size(),
              "The number of dof indices must be less than the number of settable components");
 
  for (const auto i : index_range(dof_indices))
    if (_set_components[i])
      setResidual(_sys, conversionHelper(residual, i), dof_indices[i]);
}
 
template <typename T, typename Base>
template <typename ADResidual>
void
ADNodalBCTempl<T, Base>::addJacobian(const ADResidual & residual,
                                     const std::vector<dof_id_type> & dof_indices)
{
  mooseAssert(dof_indices.size() <= _set_components.size(),
              "The number of dof indices must be less than the number of settable components");
  if (!std::is_same_v<T, RealVectorValue>)
    mooseAssert(dof_indices.size() == _var.count(),
                "The number of dof indices should match the variable count");
 
  for (const auto i : index_range(dof_indices))
    if (_set_components[i])
      // If we store into the displaced assembly for nodal bc objects the data never actually makes
      // it into the global Jacobian. A nodal bc replaces its own row, so the row must not be
      // distributed by constraints (e.g. a hanging node) onto the rows of other dofs
      this->addJacobianWithoutConstraints(_undisplaced_assembly,
                                          Moose::Span(&conversionHelper(residual, i), 1),
                                          Moose::Span(&dof_indices[i], 1),
                                          /*scaling_factor=*/1);
}
 
template <typename T, typename Base>
void
ADNodalBCTempl<T, Base>::computeResidual()
{
  const std::vector<dof_id_type> & dof_indices = _var.dofIndices();
  if (dof_indices.empty())
    return;
 
  const auto residual = MetaPhysicL::raw_value(computeQpResidual());
 
  addResidual(residual, dof_indices);
}
 
template <typename T, typename Base>
void
ADNodalBCTempl<T, Base>::computeJacobian()
{
  const std::vector<dof_id_type> & dof_indices = _var.dofIndices();
  if (dof_indices.empty())
    return;
 
  const auto residual = computeQpResidual();
 
  addJacobian(residual, dof_indices);
}
 
template <typename T, typename Base>
void
ADNodalBCTempl<T, Base>::computeResidualAndJacobian()
{
  const std::vector<dof_id_type> & dof_indices = _var.dofIndices();
  if (dof_indices.empty())
    return;
 
  const auto residual = computeQpResidual();
 
  addResidual(MetaPhysicL::raw_value(residual), dof_indices);
  addJacobian(residual, dof_indices);
}
 
template <typename T, typename Base>
void
ADNodalBCTempl<T, Base>::computeOffDiagJacobian(const unsigned int jvar_num)
{
  // Only need to do this once because AD does all the derivatives at once
  if (jvar_num == _var.number())
    computeJacobian();
}
 
template <typename T, typename Base>
void
ADNodalBCTempl<T, Base>::computeOffDiagJacobianScalar(unsigned int)
{
  // scalar coupling will have been included in the all-at-once handling in computeOffDiagJacobian
}
 
template class ADNodalBCTempl<Real, NodalBCBase>;
template class ADNodalBCTempl<RealVectorValue, NodalBCBase>;
template class ADNodalBCTempl<RealEigenVector, NodalBCBase>;
template class ADNodalBCTempl<Real, ADDirichletBCBase>;
template class ADNodalBCTempl<RealVectorValue, ADDirichletBCBase>;