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       11475 : LinearFVFluxKernel::validParams()
      15             : {
      16       11475 :   InputParameters params = LinearFVKernel::validParams();
      17       34425 :   params.addParam<bool>("force_boundary_execution",
      18       22950 :                         false,
      19             :                         "Whether to force execution of this object on all external boundaries.");
      20       11475 :   params.registerSystemAttributeName("LinearFVFluxKernel");
      21       11475 :   return params;
      22           0 : }
      23             : 
      24        1146 : LinearFVFluxKernel::LinearFVFluxKernel(const InputParameters & params)
      25             :   : LinearFVKernel(params),
      26             :     FaceArgProducerInterface(),
      27        1146 :     _current_face_info(nullptr),
      28        1146 :     _current_face_type(FaceInfo::VarFaceNeighbors::NEITHER),
      29        1146 :     _cached_matrix_contribution(false),
      30        1146 :     _cached_rhs_contribution(false),
      31        1146 :     _force_boundary_execution(getParam<bool>("force_boundary_execution")),
      32        1146 :     _dof_indices(2, 0),
      33        1146 :     _matrix_contribution(2, 2),
      34        2292 :     _rhs_contribution(2, 0.0)
      35             : {
      36        1146 : }
      37             : 
      38             : void
      39    58426149 : 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    58426149 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
      44             :   {
      45             :     // The dof ids of the variable corresponding to the element and neighbor
      46    57028585 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
      47    57028585 :     _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    57028585 :     const auto elem_matrix_contribution = computeElemMatrixContribution();
      52    57028585 :     const auto neighbor_matrix_contribution = computeNeighborMatrixContribution();
      53             : 
      54             :     // Populate matrix
      55    57028585 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()))
      56             :     {
      57    57028582 :       _matrix_contribution(0, 0) = elem_matrix_contribution;
      58    57028582 :       _matrix_contribution(0, 1) = neighbor_matrix_contribution;
      59             :     }
      60             : 
      61    57028585 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()))
      62             :     {
      63    57028568 :       _matrix_contribution(1, 0) = -elem_matrix_contribution;
      64    57028568 :       _matrix_contribution(1, 1) = -neighbor_matrix_contribution;
      65             :     }
      66             : 
      67             :     // We add the contributions to every tagged matrix
      68   114057170 :     for (auto & matrix : _matrices)
      69    57028585 :       (*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     1397564 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      74        1681 :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
      75             :   {
      76     1397564 :     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     1397564 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
      84             : 
      85     1397564 :     if (bc_pointer || _force_boundary_execution)
      86             :     {
      87     1391724 :       if (bc_pointer)
      88     1391724 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
      89     1391724 :       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     1391724 :       if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
      94             :       {
      95     1390043 :         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     2780086 :         for (auto & matrix : _matrices)
      99     1390043 :           (*matrix).add(dof_id_elem, dof_id_elem, matrix_contribution);
     100             :       }
     101        1681 :       else if (_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     102             :       {
     103             :         const auto dof_id_neighbor =
     104        1681 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     105             : 
     106             :         // We add the contributions to every tagged matrix
     107        3362 :         for (auto & matrix : _matrices)
     108        1681 :           (*matrix).add(dof_id_neighbor, dof_id_neighbor, matrix_contribution);
     109             :       }
     110             :     }
     111             :   }
     112    58426149 : }
     113             : 
     114             : void
     115    58426149 : 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    58426149 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
     120             :   {
     121             :     // The dof ids of the variable corresponding to the element and neighbor
     122    57028585 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     123    57028585 :     _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    57028585 :     const auto elem_rhs_contribution = computeElemRightHandSideContribution();
     127    57028585 :     const auto neighbor_rhs_contribution = computeNeighborRightHandSideContribution();
     128             : 
     129             :     // Populate right hand side
     130    57028585 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()))
     131    57028582 :       _rhs_contribution(0) = elem_rhs_contribution;
     132    57028585 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()))
     133    57028568 :       _rhs_contribution(1) = neighbor_rhs_contribution;
     134             : 
     135             :     // We add the contributions to every tagged vector
     136   114057170 :     for (auto & vector : _vectors)
     137    57028585 :       (*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     1397564 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     142        1681 :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     143             :   {
     144     1397564 :     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     1397564 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     151             : 
     152     1397564 :     if (bc_pointer || _force_boundary_execution)
     153             :     {
     154     1391724 :       if (bc_pointer)
     155     1391724 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     156             : 
     157     1391724 :       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     1391724 :       if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM)
     162             :       {
     163     1390043 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     164             :         // We add the contributions to every tagged vector
     165     2780086 :         for (auto & vector : _vectors)
     166     1390043 :           (*vector).add(dof_id_elem, rhs_contribution);
     167             :       }
     168        1681 :       else if (_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     169             :       {
     170             :         const auto dof_id_neighbor =
     171        1681 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     172             :         // We add the contributions to every tagged matrix
     173        3362 :         for (auto & vector : _vectors)
     174        1681 :           (*vector).add(dof_id_neighbor, rhs_contribution);
     175             :       }
     176             :     }
     177             :   }
     178    58426149 : }
     179             : 
     180             : bool
     181     1536160 : LinearFVFluxKernel::hasFaceSide(const FaceInfo & fi, bool fi_elem_side) const
     182             : {
     183     1536160 :   const auto ft = fi.faceType(std::make_pair(_var_num, _sys_num));
     184     1536160 :   if (fi_elem_side)
     185      768080 :     return ft == FaceInfo::VarFaceNeighbors::ELEM || ft == FaceInfo::VarFaceNeighbors::BOTH;
     186             :   else
     187      768080 :     return ft == FaceInfo::VarFaceNeighbors::NEIGHBOR || ft == FaceInfo::VarFaceNeighbors::BOTH;
     188             : }
     189             : 
     190             : Moose::FaceArg
     191      102271 : 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      102271 :   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    58426205 : LinearFVFluxKernel::setupFaceData(const FaceInfo * face_info)
     211             : {
     212    58426205 :   _cached_matrix_contribution = false;
     213    58426205 :   _cached_rhs_contribution = false;
     214    58426205 :   _current_face_info = face_info;
     215    58426205 :   _current_face_type = _current_face_info->faceType(std::make_pair(_var_num, _sys_num));
     216    58426205 :   _matrix_contribution.zero();
     217    58426205 :   _dof_indices.zero();
     218    58426205 :   _rhs_contribution.zero();
     219    58426205 : }