LinearFVFluxKernel.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 "LinearFVFluxKernel.h"
#include "LinearFVBoundaryCondition.h"
 
InputParameters
LinearFVFluxKernel::validParams()
{
  InputParameters params = LinearFVKernel::validParams();
  params.addParam<bool>("force_boundary_execution",
                        false,
                        "Whether to force execution of this object on all external boundaries.");
  params.registerSystemAttributeName("LinearFVFluxKernel");
  return params;
}
 
LinearFVFluxKernel::LinearFVFluxKernel(const InputParameters & params)
  : LinearFVKernel(params),
    FaceArgProducerInterface(),
    _current_face_info(nullptr),
    _current_face_type(FaceInfo::VarFaceNeighbors::NEITHER),
    _cached_matrix_contribution(false),
    _cached_rhs_contribution(false),
    _force_boundary_execution(getParam<bool>("force_boundary_execution")),
    _dof_indices(2, 0),
    _matrix_contribution(2, 2),
    _rhs_contribution(2, 0.0)
{
}
 
void
LinearFVFluxKernel::addMatrixContribution()
{
  // If we are on an internal face, we populate the four entries in the system matrix
  // which touch the face
  if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
  {
    // The dof ids of the variable corresponding to the element and neighbor
    _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
    _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
 
    // Compute the entries which will go to the neighbor (offdiagonal) and element
    // (diagonal).
    const auto elem_matrix_contribution = computeElemMatrixContribution();
    const auto neighbor_matrix_contribution = computeNeighborMatrixContribution();
 
    // Populate matrix
    if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()))
    {
      _matrix_contribution(0, 0) = elem_matrix_contribution;
      _matrix_contribution(0, 1) = neighbor_matrix_contribution;
    }
 
    if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()))
    {
      _matrix_contribution(1, 0) = -elem_matrix_contribution;
      _matrix_contribution(1, 1) = -neighbor_matrix_contribution;
    }
 
    // We add the contributions to every tagged matrix
    for (auto & matrix : _matrices)
      (*matrix).add_matrix(_matrix_contribution, _dof_indices.get_values());
  }
  // We are at a block boundary where the variable is not defined on one of the adjacent cells.
  // We check if we have a boundary condition here
  else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
           _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
  {
    if (_current_face_info->boundaryIDs().size() > 1)
      mooseError("We currently don't support multiple boundary conditions for the same variable on "
                 "the same face. Current face center : " +
                 Moose::stringify(_current_face_info->faceCentroid()) +
                 " boundaries specified: " + Moose::stringify(_current_face_info->boundaryIDs()));
 
    LinearFVBoundaryCondition * bc_pointer =
        _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
 
    if (bc_pointer || _force_boundary_execution)
    {
      if (bc_pointer)
        bc_pointer->setupFaceData(_current_face_info, _current_face_type);
      const auto matrix_contribution = computeBoundaryMatrixContribution(*bc_pointer);
 
      // We allow internal (for the mesh) boundaries too, so we have to check on which side we
      // are on (assuming that this is a boundary for the variable)
      if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
      {
        const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
 
        // We add the contributions to every tagged matrix
        for (auto & matrix : _matrices)
          (*matrix).add(dof_id_elem, dof_id_elem, matrix_contribution);
      }
      else if (_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
      {
        const auto dof_id_neighbor =
            _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
 
        // We add the contributions to every tagged matrix
        for (auto & matrix : _matrices)
          (*matrix).add(dof_id_neighbor, dof_id_neighbor, matrix_contribution);
      }
    }
  }
}
 
void
LinearFVFluxKernel::addRightHandSideContribution()
{
  // If we are on an internal face, we populate the two entries in the right hand side
  // which touch the face
  if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
  {
    // The dof ids of the variable corresponding to the element and neighbor
    _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
    _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
 
    // Compute the entries which will go to the neighbor and element positions.
    const auto elem_rhs_contribution = computeElemRightHandSideContribution();
    const auto neighbor_rhs_contribution = computeNeighborRightHandSideContribution();
 
    // Populate right hand side
    if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()))
      _rhs_contribution(0) = elem_rhs_contribution;
    if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()))
      _rhs_contribution(1) = neighbor_rhs_contribution;
 
    // We add the contributions to every tagged vector
    for (auto & vector : _vectors)
      (*vector).add_vector(_rhs_contribution.get_values().data(), _dof_indices.get_values());
  }
  // We are at a block boundary where the variable is not defined on one of the adjacent cells.
  // We check if we have a boundary condition here
  else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
           _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
  {
    if (_current_face_info->boundaryIDs().size() > 1)
      mooseError("We currently don't support multiple boundary conditions for the same variable on "
                 "the same face. Current face center : " +
                 Moose::stringify(_current_face_info->faceCentroid()) +
                 " boundaries specified: " + Moose::stringify(_current_face_info->boundaryIDs()));
    LinearFVBoundaryCondition * bc_pointer =
        _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
 
    if (bc_pointer || _force_boundary_execution)
    {
      if (bc_pointer)
        bc_pointer->setupFaceData(_current_face_info, _current_face_type);
 
      const auto rhs_contribution = computeBoundaryRHSContribution(*bc_pointer);
 
      // We allow internal (for the mesh) boundaries too, so we have to check on which side we
      // are on (assuming that this is a boundary for the variable)
      if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
      {
        const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
        // We add the contributions to every tagged vector
        for (auto & vector : _vectors)
          (*vector).add(dof_id_elem, rhs_contribution);
      }
      else if (_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
      {
        const auto dof_id_neighbor =
            _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
        // We add the contributions to every tagged matrix
        for (auto & vector : _vectors)
          (*vector).add(dof_id_neighbor, rhs_contribution);
      }
    }
  }
}
 
bool
LinearFVFluxKernel::hasFaceSide(const FaceInfo & fi, bool fi_elem_side) const
{
  const auto ft = fi.faceType(std::make_pair(_var_num, _sys_num));
  if (fi_elem_side)
    return ft == FaceInfo::VarFaceNeighbors::ELEM || ft == FaceInfo::VarFaceNeighbors::BOTH;
  else
    return ft == FaceInfo::VarFaceNeighbors::NEIGHBOR || ft == FaceInfo::VarFaceNeighbors::BOTH;
}
 
Moose::FaceArg
LinearFVFluxKernel::singleSidedFaceArg(const FaceInfo * fi,
                                       const Moose::FV::LimiterType limiter_type,
                                       const bool correct_skewness) const
{
  mooseAssert(fi, "FaceInfo should not be null!");
  return makeFace(*fi, limiter_type, true, correct_skewness);
}
 
Real
LinearFVFluxKernel::computeBoundaryFlux(const LinearFVBoundaryCondition & bc)
{
  const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
                             ? _current_face_info->elemInfo()
                             : _current_face_info->neighborInfo();
  return computeBoundaryMatrixContribution(bc) * _var.getElemValue(*elem_info, determineState()) -
         computeBoundaryRHSContribution(bc);
}
 
void
LinearFVFluxKernel::setupFaceData(const FaceInfo * face_info)
{
  _cached_matrix_contribution = false;
  _cached_rhs_contribution = false;
  _current_face_info = face_info;
  _current_face_type = _current_face_info->faceType(std::make_pair(_var_num, _sys_num));
  _matrix_contribution.zero();
  _dof_indices.zero();
  _rhs_contribution.zero();
}