Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "LinearFVFluxKernel.h"
      11             : #include "LinearFVBoundaryCondition.h"
      12             : 
      13             : InputParameters
      14       11589 : LinearFVFluxKernel::validParams()
      15             : {
      16       11589 :   InputParameters params = LinearFVKernel::validParams();
      17       34767 :   params.addParam<bool>("force_boundary_execution",
      18       23178 :                         false,
      19             :                         "Whether to force execution of this object on all external boundaries.");
      20       11589 :   params.registerSystemAttributeName("LinearFVFluxKernel");
      21       11589 :   return params;
      22           0 : }
      23             : 
      24        1161 : LinearFVFluxKernel::LinearFVFluxKernel(const InputParameters & params)
      25             :   : LinearFVKernel(params),
      26             :     FaceArgProducerInterface(),
      27        1161 :     _current_face_info(nullptr),
      28        1161 :     _current_face_type(FaceInfo::VarFaceNeighbors::NEITHER),
      29        1161 :     _cached_matrix_contribution(false),
      30        1161 :     _cached_rhs_contribution(false),
      31        1161 :     _force_boundary_execution(getParam<bool>("force_boundary_execution")),
      32        1161 :     _dof_indices(2, 0),
      33        1161 :     _matrix_contribution(2, 2),
      34        2322 :     _rhs_contribution(2, 0.0)
      35             : {
      36        1161 : }
      37             : 
      38             : void
      39    57669530 : LinearFVFluxKernel::addMatrixContribution()
      40             : {
      41             :   // If we are on an internal face, we populate the four entries in the system matrix
      42             :   // which touch the face
      43    57669530 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
      44             :   {
      45             :     // The dof ids of the variable corresponding to the element and neighbor
      46    56307129 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
      47    56307129 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
      48             : 
      49             :     // Compute the entries which will go to the neighbor (offdiagonal) and element
      50             :     // (diagonal).
      51    56307129 :     const auto elem_matrix_contribution = computeElemMatrixContribution();
      52    56307129 :     const auto neighbor_matrix_contribution = computeNeighborMatrixContribution();
      53             : 
      54             :     // Populate matrix
      55    56307129 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()))
      56             :     {
      57    56307126 :       _matrix_contribution(0, 0) = elem_matrix_contribution;
      58    56307126 :       _matrix_contribution(0, 1) = neighbor_matrix_contribution;
      59             :     }
      60             : 
      61    56307129 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()))
      62             :     {
      63    56307119 :       _matrix_contribution(1, 0) = -elem_matrix_contribution;
      64    56307119 :       _matrix_contribution(1, 1) = -neighbor_matrix_contribution;
      65             :     }
      66             : 
      67             :     // We add the contributions to every tagged matrix
      68   112614258 :     for (auto & matrix : _matrices)
      69    56307129 :       (*matrix).add_matrix(_matrix_contribution, _dof_indices.get_values());
      70             :   }
      71             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
      72             :   // We check if we have a boundary condition here
      73     1362401 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      74         844 :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
      75             :   {
      76     1362401 :     if (_current_face_info->boundaryIDs().size() > 1)
      77           0 :       mooseError("We currently don't support multiple boundary conditions for the same variable on "
      78           0 :                  "the same face. Current face center : " +
      79           0 :                  Moose::stringify(_current_face_info->faceCentroid()) +
      80           0 :                  " boundaries specified: " + Moose::stringify(_current_face_info->boundaryIDs()));
      81             : 
      82             :     LinearFVBoundaryCondition * bc_pointer =
      83     1362401 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
      84             : 
      85     1362401 :     if (bc_pointer || _force_boundary_execution)
      86             :     {
      87     1356561 :       if (bc_pointer)
      88     1356561 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
      89     1356561 :       const auto matrix_contribution = computeBoundaryMatrixContribution(*bc_pointer);
      90             : 
      91             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
      92             :       // are on (assuming that this is a boundary for the variable)
      93     1356561 :       if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
      94             :       {
      95     1355717 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
      96             : 
      97             :         // We add the contributions to every tagged matrix
      98     2711434 :         for (auto & matrix : _matrices)
      99     1355717 :           (*matrix).add(dof_id_elem, dof_id_elem, matrix_contribution);
     100             :       }
     101         844 :       else if (_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     102             :       {
     103             :         const auto dof_id_neighbor =
     104         844 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     105             : 
     106             :         // We add the contributions to every tagged matrix
     107        1688 :         for (auto & matrix : _matrices)
     108         844 :           (*matrix).add(dof_id_neighbor, dof_id_neighbor, matrix_contribution);
     109             :       }
     110             :     }
     111             :   }
     112    57669530 : }
     113             : 
     114             : void
     115    57669530 : LinearFVFluxKernel::addRightHandSideContribution()
     116             : {
     117             :   // If we are on an internal face, we populate the two entries in the right hand side
     118             :   // which touch the face
     119    57669530 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
     120             :   {
     121             :     // The dof ids of the variable corresponding to the element and neighbor
     122    56307129 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     123    56307129 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     124             : 
     125             :     // Compute the entries which will go to the neighbor and element positions.
     126    56307129 :     const auto elem_rhs_contribution = computeElemRightHandSideContribution();
     127    56307129 :     const auto neighbor_rhs_contribution = computeNeighborRightHandSideContribution();
     128             : 
     129             :     // Populate right hand side
     130    56307129 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()))
     131    56307126 :       _rhs_contribution(0) = elem_rhs_contribution;
     132    56307129 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()))
     133    56307119 :       _rhs_contribution(1) = neighbor_rhs_contribution;
     134             : 
     135             :     // We add the contributions to every tagged vector
     136   112614258 :     for (auto & vector : _vectors)
     137    56307129 :       (*vector).add_vector(_rhs_contribution.get_values().data(), _dof_indices.get_values());
     138             :   }
     139             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
     140             :   // We check if we have a boundary condition here
     141     1362401 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     142         844 :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     143             :   {
     144     1362401 :     if (_current_face_info->boundaryIDs().size() > 1)
     145           0 :       mooseError("We currently don't support multiple boundary conditions for the same variable on "
     146           0 :                  "the same face. Current face center : " +
     147           0 :                  Moose::stringify(_current_face_info->faceCentroid()) +
     148           0 :                  " boundaries specified: " + Moose::stringify(_current_face_info->boundaryIDs()));
     149             :     LinearFVBoundaryCondition * bc_pointer =
     150     1362401 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     151             : 
     152     1362401 :     if (bc_pointer || _force_boundary_execution)
     153             :     {
     154     1356561 :       if (bc_pointer)
     155     1356561 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     156             : 
     157     1356561 :       const auto rhs_contribution = computeBoundaryRHSContribution(*bc_pointer);
     158             : 
     159             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     160             :       // are on (assuming that this is a boundary for the variable)
     161     1356561 :       if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
     162             :       {
     163     1355717 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     164             :         // We add the contributions to every tagged vector
     165     2711434 :         for (auto & vector : _vectors)
     166     1355717 :           (*vector).add(dof_id_elem, rhs_contribution);
     167             :       }
     168         844 :       else if (_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     169             :       {
     170             :         const auto dof_id_neighbor =
     171         844 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     172             :         // We add the contributions to every tagged matrix
     173        1688 :         for (auto & vector : _vectors)
     174         844 :           (*vector).add(dof_id_neighbor, rhs_contribution);
     175             :       }
     176             :     }
     177             :   }
     178    57669530 : }
     179             : 
     180             : bool
     181     1254866 : LinearFVFluxKernel::hasFaceSide(const FaceInfo & fi, bool fi_elem_side) const
     182             : {
     183     1254866 :   const auto ft = fi.faceType(std::make_pair(_var_num, _sys_num));
     184     1254866 :   if (fi_elem_side)
     185      627433 :     return ft == FaceInfo::VarFaceNeighbors::ELEM || ft == FaceInfo::VarFaceNeighbors::BOTH;
     186             :   else
     187      627433 :     return ft == FaceInfo::VarFaceNeighbors::NEIGHBOR || ft == FaceInfo::VarFaceNeighbors::BOTH;
     188             : }
     189             : 
     190             : Moose::FaceArg
     191       83742 : LinearFVFluxKernel::singleSidedFaceArg(const FaceInfo * fi,
     192             :                                        const Moose::FV::LimiterType limiter_type,
     193             :                                        const bool correct_skewness) const
     194             : {
     195             :   mooseAssert(fi, "FaceInfo should not be null!");
     196       83742 :   return makeFace(*fi, limiter_type, true, correct_skewness);
     197             : }
     198             : 
     199             : Real
     200          56 : LinearFVFluxKernel::computeBoundaryFlux(const LinearFVBoundaryCondition & bc)
     201             : {
     202          56 :   const auto elem_info = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
     203          56 :                              ? _current_face_info->elemInfo()
     204           0 :                              : _current_face_info->neighborInfo();
     205         112 :   return computeBoundaryMatrixContribution(bc) * _var.getElemValue(*elem_info, determineState()) -
     206         112 :          computeBoundaryRHSContribution(bc);
     207             : }
     208             : 
     209             : void
     210    57669586 : LinearFVFluxKernel::setupFaceData(const FaceInfo * face_info)
     211             : {
     212    57669586 :   _cached_matrix_contribution = false;
     213    57669586 :   _cached_rhs_contribution = false;
     214    57669586 :   _current_face_info = face_info;
     215    57669586 :   _current_face_type = _current_face_info->faceType(std::make_pair(_var_num, _sys_num));
     216    57669586 :   _matrix_contribution.zero();
     217    57669586 :   _dof_indices.zero();
     218    57669586 :   _rhs_contribution.zero();
     219    57669586 : }