Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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 "LinearFVTurbulentAdvection.h"
      11             : #include "MooseLinearVariableFV.h"
      12             : #include "NSFVUtils.h"
      13             : #include "NavierStokesMethods.h"
      14             : #include "NS.h"
      15             : 
      16             : registerMooseObject("NavierStokesApp", LinearFVTurbulentAdvection);
      17             : 
      18             : InputParameters
      19         176 : LinearFVTurbulentAdvection::validParams()
      20             : {
      21         176 :   InputParameters params = LinearFVScalarAdvection::validParams();
      22         176 :   params.addClassDescription("Represents the matrix and right hand side contributions of an "
      23             :                              "advection term for a turbulence variable.");
      24             : 
      25         352 :   params.addParam<std::vector<BoundaryName>>(
      26             :       "walls", {}, "Boundaries that correspond to solid walls.");
      27             : 
      28         176 :   return params;
      29           0 : }
      30             : 
      31          88 : LinearFVTurbulentAdvection::LinearFVTurbulentAdvection(const InputParameters & params)
      32             :   : LinearFVScalarAdvection(params),
      33          88 :     _advected_interp_coeffs(std::make_pair<Real, Real>(0, 0)),
      34         176 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls"))
      35             : {
      36          88 :   Moose::FV::setInterpolationMethod(*this, _advected_interp_method, "advected_interp_method");
      37          88 : }
      38             : 
      39             : void
      40         440 : LinearFVTurbulentAdvection::initialSetup()
      41             : {
      42         440 :   LinearFVScalarAdvection::initialSetup();
      43         880 :   NS::getWallBoundedElements(
      44         440 :       _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
      45         440 : }
      46             : 
      47             : void
      48     2516936 : LinearFVTurbulentAdvection::addMatrixContribution()
      49             : {
      50             :   // Coumputing bounding map
      51     2516936 :   const Elem * elem = _current_face_info->elemPtr();
      52             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
      53     4424080 :   const Elem * neighbor = _current_face_info->neighborPtr();
      54             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
      55             : 
      56             :   // If we are on an internal face, we populate the four entries in the system matrix
      57             :   // which touch the face
      58     2516936 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
      59             :   {
      60             :     // The dof ids of the variable corresponding to the element and neighbor
      61     1875000 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
      62     1875000 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
      63             : 
      64             :     // Compute the entries which will go to the neighbor (offdiagonal) and element
      65             :     // (diagonal).
      66     1875000 :     const auto elem_matrix_contribution = computeElemMatrixContribution();
      67     1875000 :     const auto neighbor_matrix_contribution = computeNeighborMatrixContribution();
      68             : 
      69             :     // Populate matrix
      70     1875000 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()) && !(bounded_elem))
      71             :     {
      72     1682290 :       _matrix_contribution(0, 0) = elem_matrix_contribution;
      73     1682290 :       _matrix_contribution(0, 1) = neighbor_matrix_contribution;
      74             :     }
      75             : 
      76     1875000 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()) && !(bounded_neigh))
      77             :     {
      78     1684410 :       _matrix_contribution(1, 0) = -elem_matrix_contribution;
      79     1684410 :       _matrix_contribution(1, 1) = -neighbor_matrix_contribution;
      80             :     }
      81     1875000 :     (*_linear_system.matrix).add_matrix(_matrix_contribution, _dof_indices.get_values());
      82             :   }
      83             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
      84             :   // We check if we have a boundary condition here
      85      641936 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      86             :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
      87             :   {
      88             :     mooseAssert(_current_face_info->boundaryIDs().size() == 1,
      89             :                 "We should only have one boundary on every face.");
      90             : 
      91             :     LinearFVBoundaryCondition * bc_pointer =
      92      641936 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
      93             : 
      94      641936 :     if (bc_pointer || _force_boundary_execution)
      95             :     {
      96             :       if (bc_pointer)
      97      349712 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
      98      349712 :       const auto matrix_contribution = computeBoundaryMatrixContribution(*bc_pointer);
      99             : 
     100             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     101             :       // are on (assuming that this is a boundary for the variable)
     102      349712 :       if ((_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) && !(bounded_elem))
     103             :       {
     104      298916 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     105      298916 :         (*_linear_system.matrix).add(dof_id_elem, dof_id_elem, matrix_contribution);
     106             :       }
     107       50796 :       else if ((_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) && !(bounded_neigh))
     108             :       {
     109             :         const auto dof_id_neighbor =
     110           0 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     111           0 :         (*_linear_system.matrix).add(dof_id_neighbor, dof_id_neighbor, matrix_contribution);
     112             :       }
     113             :     }
     114             :   }
     115     2516936 : }
     116             : 
     117             : void
     118     2516936 : LinearFVTurbulentAdvection::addRightHandSideContribution()
     119             : {
     120             :   // Coumputing bounding map
     121     2516936 :   const Elem * elem = _current_face_info->elemPtr();
     122             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
     123     4424080 :   const Elem * neighbor = _current_face_info->neighborPtr();
     124             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
     125             : 
     126             :   // If we are on an internal face, we populate the two entries in the right hand side
     127             :   // which touch the face
     128     2516936 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
     129             :   {
     130             :     // The dof ids of the variable corresponding to the element and neighbor
     131     1875000 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     132     1875000 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     133             : 
     134             :     // Compute the entries which will go to the neighbor and element positions.
     135     1875000 :     const auto elem_rhs_contribution = computeElemRightHandSideContribution();
     136     1875000 :     const auto neighbor_rhs_contribution = computeNeighborRightHandSideContribution();
     137             : 
     138             :     // Populate right hand side
     139     1875000 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()) && !(bounded_elem))
     140     1682290 :       _rhs_contribution(0) = elem_rhs_contribution;
     141     1875000 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()) && !(bounded_neigh))
     142     1684410 :       _rhs_contribution(1) = neighbor_rhs_contribution;
     143             : 
     144     1875000 :     (*_linear_system.rhs)
     145     1875000 :         .add_vector(_rhs_contribution.get_values().data(), _dof_indices.get_values());
     146             :   }
     147             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
     148             :   // We check if we have a boundary condition here
     149      641936 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     150             :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     151             :   {
     152             :     mooseAssert(_current_face_info->boundaryIDs().size() == 1,
     153             :                 "We should only have one boundary on every face.");
     154             :     LinearFVBoundaryCondition * bc_pointer =
     155      641936 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     156             : 
     157      641936 :     if (bc_pointer || _force_boundary_execution)
     158             :     {
     159             :       if (bc_pointer)
     160      349712 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     161             : 
     162      349712 :       const auto rhs_contribution = computeBoundaryRHSContribution(*bc_pointer);
     163             : 
     164             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     165             :       // are on (assuming that this is a boundary for the variable)
     166      349712 :       if ((_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) && !(bounded_elem))
     167             :       {
     168      298916 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     169      298916 :         (*_linear_system.rhs).add(dof_id_elem, rhs_contribution);
     170             :       }
     171       50796 :       else if ((_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) && !(bounded_neigh))
     172             :       {
     173             :         const auto dof_id_neighbor =
     174           0 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     175           0 :         (*_linear_system.rhs).add(dof_id_neighbor, rhs_contribution);
     176             :       }
     177             :     }
     178             :   }
     179     2516936 : }