FVInterfaceKernel.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 "FVInterfaceKernel.h"
#include "FEProblemBase.h"
#include "SystemBase.h"
#include "MooseVariableFV.h"
#include "Assembly.h"
#include "MathFVUtils.h"
 
InputParameters
FVInterfaceKernel::validParams()
{
  InputParameters params = MooseObject::validParams();
  params += BoundaryRestrictableRequired::validParams();
  params += SetupInterface::validParams();
  params += FunctionInterface::validParams();
  params += DistributionInterface::validParams();
  params += UserObjectInterface::validParams();
  params += TransientInterface::validParams();
  params += PostprocessorInterface::validParams();
  params += VectorPostprocessorInterface::validParams();
  params += GeometricSearchInterface::validParams();
  params += MeshChangedInterface::validParams();
  params += TaggingInterface::validParams();
  params += NeighborCoupleableMooseVariableDependencyIntermediateInterface::validParams();
  params += TwoMaterialPropertyInterface::validParams();
  params += ADFunctorInterface::validParams();
 
  params.addRequiredParam<std::vector<SubdomainName>>(
      "subdomain1", "The subdomains on the 1st side of the boundary.");
  params.addRequiredParam<std::vector<SubdomainName>>(
      "subdomain2", "The subdomains on the 2nd side of the boundary.");
  params.addRequiredParam<NonlinearVariableName>(
      "variable1", "The name of the first variable that this interface kernel applies to");
  params.addParam<NonlinearVariableName>(
      "variable2",
      "The name of the second variable that this interface kernel applies to. If not supplied, "
      "variable1 will be used.");
  params.addParam<bool>("use_displaced_mesh",
                        false,
                        "Whether or not this object should use the "
                        "displaced mesh for computation.  Note that "
                        "in the case this is true but no "
                        "displacements are provided in the Mesh block "
                        "the undisplaced mesh will still be used.");
 
  params.addParam<unsigned short>("ghost_layers", 1, "The number of layers of elements to ghost.");
  params.addParam<bool>("use_point_neighbors",
                        false,
                        "Whether to use point neighbors, which introduces additional ghosting to "
                        "that used for simple face neighbors.");
  params.addParamNamesToGroup("ghost_layers use_point_neighbors", "Parallel ghosting");
 
  // FV Interface Kernels always need one layer of ghosting because the elements
  // on each side of the interface may be on different MPI ranks, but we still
  // need to access them as a pair to compute the numerical face flux.
  params.addRelationshipManager(
      "ElementSideNeighborLayers",
      Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
          Moose::RelationshipManagerType::COUPLING,
      [](const InputParameters & obj_params, InputParameters & rm_params)
      {
        rm_params.set<unsigned short>("layers") = obj_params.get<unsigned short>("ghost_layers");
        rm_params.set<bool>("use_point_neighbors") = obj_params.get<bool>("use_point_neighbors");
      });
 
  params.addParamNamesToGroup("use_displaced_mesh", "Advanced");
  params.addCoupledVar("displacements", "The displacements");
  params.declareControllable("enable");
  params.registerBase("FVInterfaceKernel");
  params.registerSystemAttributeName("FVInterfaceKernel");
  params.set<bool>("_residual_object") = true;
  return params;
}
 
FVInterfaceKernel::FVInterfaceKernel(const InputParameters & parameters)
  : MooseObject(parameters),
    BoundaryRestrictableRequired(this, false),
    SetupInterface(this),
    FunctionInterface(this),
    DistributionInterface(this),
    UserObjectInterface(this),
    TransientInterface(this),
    PostprocessorInterface(this),
    VectorPostprocessorInterface(this),
    GeometricSearchInterface(this),
    MeshChangedInterface(parameters),
    TaggingInterface(this),
    NeighborCoupleableMooseVariableDependencyIntermediateInterface(
        this, /*nodal=*/false, /*neighbor_nodal=*/false, /*is_fv=*/true),
    TwoMaterialPropertyInterface(this, Moose::EMPTY_BLOCK_IDS, boundaryIDs()),
    ADFunctorInterface(this),
    _tid(getParam<THREAD_ID>("_tid")),
    _subproblem(*getCheckedPointerParam<SubProblem *>("_subproblem")),
    _var1(_subproblem.getVariable(_tid, getParam<NonlinearVariableName>("variable1"))
              .sys()
              .getFVVariable<Real>(_tid, getParam<NonlinearVariableName>("variable1"))),
    _var2(_subproblem
              .getVariable(_tid,
                           isParamValid("variable2") ? getParam<NonlinearVariableName>("variable2")
                                                     : getParam<NonlinearVariableName>("variable1"))
              .sys()
              .getFVVariable<Real>(_tid,
                                   isParamValid("variable2")
                                       ? getParam<NonlinearVariableName>("variable2")
                                       : getParam<NonlinearVariableName>("variable1"))),
    _assembly(_subproblem.assembly(_tid, _var1.sys().number())),
    _mesh(_subproblem.mesh())
{
  if (getParam<bool>("use_displaced_mesh"))
    paramError("use_displaced_mesh", "FV interface kernels do not yet support displaced mesh");
 
  _subproblem.haveADObjects(true);
  addMooseVariableDependency(&_var1);
  addMooseVariableDependency(&_var2);
 
  for (const auto & sub_name : getParam<std::vector<SubdomainName>>("subdomain1"))
    _subdomain1.insert(_mesh.getSubdomainID(sub_name));
  for (const auto & sub_name : getParam<std::vector<SubdomainName>>("subdomain2"))
    _subdomain2.insert(_mesh.getSubdomainID(sub_name));
 
  if (!_var1.hasBlocks(_subdomain1))
    paramError("variable1",
               "variable1 does not exist on all the blocks specified by the subdomain1 parameter");
  if (!_var2.hasBlocks(_subdomain2))
  {
    const bool var2_provided = isParamValid("variable2");
    const std::string var_name = var2_provided ? "variable2" : "variable1";
    paramError(var_name,
               var_name,
               " does not exist on all the blocks specified by the subdomain2 parameter",
               var2_provided ? ""
                             : ".\nNote that you did not provide the variable2 parameter, "
                               "so variable1 was implicitly used on subdomain2");
  }
}
 
bool
FVInterfaceKernel::checkFaceIntegrity(const FaceInfo & fi) const
{
  if (!fi.neighborPtr())
    mooseError("FVInterfaceKernel '",
               name(),
               "' is being triggered on an external face with centroid ",
               fi.faceCentroid(),
               ". An interface kernel requires elements on both sides of the face.");
 
  const auto elem_subdomain = fi.elem().subdomain_id();
  const auto neighbor_subdomain = fi.neighbor().subdomain_id();
 
  const bool elem_on_side1 =
      _subdomain1.count(elem_subdomain) && _subdomain2.count(neighbor_subdomain);
  const bool elem_on_side2 =
      _subdomain2.count(elem_subdomain) && _subdomain1.count(neighbor_subdomain);
 
  if (!elem_on_side1 && !elem_on_side2)
    mooseError("FVInterfaceKernel '",
               name(),
               "' is being triggered on a face with centroid ",
               fi.faceCentroid(),
               " connecting subdomains ",
               elem_subdomain,
               " and ",
               neighbor_subdomain,
               ", but this face does not connect the subdomains specified by 'subdomain1' and "
               "'subdomain2'.");
 
  if (elem_on_side1 && elem_on_side2)
    mooseError("FVInterfaceKernel '",
               name(),
               "' is being triggered on a face with centroid ",
               fi.faceCentroid(),
               " connecting subdomains ",
               elem_subdomain,
               " and ",
               neighbor_subdomain,
               ", but both orientations satisfy the user-defined 'subdomain1' and 'subdomain2' "
               "restrictions, making the interface orientation ambiguous.");
 
  return elem_on_side1;
}
 
void
FVInterfaceKernel::setupData(const FaceInfo & fi)
{
  _face_info = &fi;
  _face_info_elem_on_side1 = checkFaceIntegrity(fi);
 
#ifndef NDEBUG
  const auto ft1 = fi.faceType(std::make_pair(_var1.number(), _var1.sys().number()));
  const auto ft2 = fi.faceType(std::make_pair(_var2.number(), _var2.sys().number()));
  constexpr auto ft_both = FaceInfo::VarFaceNeighbors::BOTH;
  constexpr auto ft_elem = FaceInfo::VarFaceNeighbors::ELEM;
  constexpr auto ft_neigh = FaceInfo::VarFaceNeighbors::NEIGHBOR;
  mooseAssert((_face_info_elem_on_side1 && (ft1 == ft_elem || ft1 == ft_both) &&
               (ft2 == ft_neigh || ft2 == ft_both)) ||
                  (!_face_info_elem_on_side1 && (ft1 == ft_neigh || ft1 == ft_both) &&
                   (ft2 == ft_elem || ft2 == ft_both)),
              "Face type was not recognized. Check that the specified boundaries are interfaces.");
#endif
}
 
void
FVInterfaceKernel::addResidual(const Real resid, const unsigned int var_num, const bool neighbor)
{
  neighbor ? prepareVectorTagNeighbor(_assembly, var_num) : prepareVectorTag(_assembly, var_num);
  _local_re(0) = resid;
  accumulateTaggedLocalResidual();
}
 
void
FVInterfaceKernel::addResidualToVariable1(const Real residual)
{
  addResidual(residual, _var1.number(), !_face_info_elem_on_side1);
}
 
void
FVInterfaceKernel::addResidualToVariable2(const Real residual)
{
  addResidual(residual, _var2.number(), _face_info_elem_on_side1);
}
 
void
FVInterfaceKernel::addResidualAndJacobian(const ADReal & residual,
                                          const MooseVariableFV<Real> & variable,
                                          const bool side_is_one)
{
  const bool variable_is_on_elem = side_is_one == _face_info_elem_on_side1;
  const auto & dof_indices =
      variable_is_on_elem ? variable.dofIndices() : variable.dofIndicesNeighbor();
  mooseAssert(dof_indices.size() == 1, "We're currently built to use CONSTANT MONOMIALS");
  addResidualsAndJacobian(
      _assembly, std::array<ADReal, 1>{{residual}}, dof_indices, variable.scalingFactor());
}
 
void
FVInterfaceKernel::addResidualAndJacobianToVariable1(const ADReal & residual)
{
  addResidualAndJacobian(residual, _var1, true);
}
 
void
FVInterfaceKernel::addResidualAndJacobianToVariable2(const ADReal & residual)
{
  addResidualAndJacobian(residual, _var2, false);
}
 
void
FVInterfaceKernel::computeResidual(const FaceInfo & fi)
{
  setupData(fi);
 
  const auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual());
 
  // If the two variables belong to two different nonlinear systems, we only contribute to the one
  // which is being assembled right now
  mooseAssert(_var1.sys().number() == _subproblem.currentNlSysNum(),
              "The interface kernel should contribute to the system which variable1 belongs to!");
  addResidualToVariable1(r);
  if (_var1.sys().number() == _var2.sys().number())
    addResidualToVariable2(-r);
}
 
void
FVInterfaceKernel::computeResidualAndJacobian(const FaceInfo & fi)
{
  computeJacobian(fi);
}
 
void
FVInterfaceKernel::computeJacobian(const FaceInfo & fi)
{
  setupData(fi);
 
  const auto r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
 
  // If the two variables belong to two different nonlinear systems, we only contribute to the one
  // which is being assembled right now
  mooseAssert(_var1.sys().number() == _subproblem.currentNlSysNum(),
              "The interface kernel should contribute to the system which variable1 belongs to!");
  addResidualAndJacobianToVariable1(r);
  if (_var1.sys().number() == _var2.sys().number())
    addResidualAndJacobianToVariable2(-r);
}
 
Point
FVInterfaceKernel::normal() const
{
  Point oriented_normal = _face_info->normal();
  if (!_face_info_elem_on_side1)
    oriented_normal *= -1;
  return oriented_normal;
}
 
const Elem &
FVInterfaceKernel::elem1() const
{
  return _face_info_elem_on_side1 ? _face_info->elem() : _face_info->neighbor();
}
 
const Elem &
FVInterfaceKernel::elem2() const
{
  return _face_info_elem_on_side1 ? _face_info->neighbor() : _face_info->elem();
}
 
const Point &
FVInterfaceKernel::centroid1() const
{
  return _face_info_elem_on_side1 ? _face_info->elemCentroid() : _face_info->neighborCentroid();
}
 
const Point &
FVInterfaceKernel::centroid2() const
{
  return _face_info_elem_on_side1 ? _face_info->neighborCentroid() : _face_info->elemCentroid();
}
 
Moose::ElemArg
FVInterfaceKernel::elemArg1(const bool correct_skewness) const
{
  return {&elem1(), correct_skewness};
}
 
Moose::ElemArg
FVInterfaceKernel::elemArg2(const bool correct_skewness) const
{
  return {&elem2(), correct_skewness};
}
 
Moose::FaceArg
FVInterfaceKernel::faceArg1(const Moose::FV::LimiterType limiter_type,
                            const bool correct_skewness,
                            const Moose::StateArg * state_limiter) const
{
  return {_face_info,
          limiter_type,
          _face_info_elem_on_side1,
          correct_skewness,
          &elem1(),
          state_limiter};
}
 
Moose::FaceArg
FVInterfaceKernel::faceArg2(const Moose::FV::LimiterType limiter_type,
                            const bool correct_skewness,
                            const Moose::StateArg * state_limiter) const
{
  return {_face_info,
          limiter_type,
          !_face_info_elem_on_side1,
          correct_skewness,
          &elem2(),
          state_limiter};
}
 
ADReal
FVInterfaceKernel::interpolateValue(const Moose::FV::InterpMethod method,
                                    const ADReal & value1,
                                    const ADReal & value2) const
{
  ADReal result;
  Moose::FV::interpolate(method, result, value1, value2, *_face_info, _face_info_elem_on_side1);
  return result;
}
 
bool
FVInterfaceKernel::hasFaceSide(const FaceInfo &, bool) const
{
  // Our default interface kernel treats elem and neighbor sides equivalently so we will assume for
  // now that we will happily consume functor evaluations on either side of a face and any
  // interpolation between said evaluations
  return true;
}