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         212 : LinearFVTurbulentAdvection::validParams()
      20             : {
      21         212 :   InputParameters params = LinearFVFluxKernel::validParams();
      22         212 :   params.addClassDescription("Represents the matrix and right hand side contributions of an "
      23             :                              "advection term for a turbulence variable.");
      24             : 
      25         424 :   params.addRequiredParam<UserObjectName>(
      26             :       "rhie_chow_user_object",
      27             :       "The rhie-chow user-object which is used to determine the face velocity.");
      28             : 
      29         212 :   params += Moose::FV::advectedInterpolationParameter();
      30             : 
      31         424 :   params.addParam<std::vector<BoundaryName>>(
      32             :       "walls", {}, "Boundaries that correspond to solid walls.");
      33             : 
      34         212 :   return params;
      35           0 : }
      36             : 
      37         106 : LinearFVTurbulentAdvection::LinearFVTurbulentAdvection(const InputParameters & params)
      38             :   : LinearFVFluxKernel(params),
      39         106 :     _mass_flux_provider(getUserObject<RhieChowMassFlux>("rhie_chow_user_object")),
      40         106 :     _advected_interp_coeffs(std::make_pair<Real, Real>(0, 0)),
      41         106 :     _mass_face_flux(0.0),
      42         318 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls"))
      43             : {
      44         106 :   Moose::FV::setInterpolationMethod(*this, _advected_interp_method, "advected_interp_method");
      45         106 : }
      46             : 
      47             : void
      48         530 : LinearFVTurbulentAdvection::initialSetup()
      49             : {
      50         530 :   LinearFVFluxKernel::initialSetup();
      51        1060 :   NS::getWallBoundedElements(
      52         530 :       _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
      53         530 : }
      54             : 
      55             : void
      56     2821392 : LinearFVTurbulentAdvection::addMatrixContribution()
      57             : {
      58             :   // Coumputing bounding map
      59     2821392 :   const Elem * elem = _current_face_info->elemPtr();
      60             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
      61     4936848 :   const Elem * neighbor = _current_face_info->neighborPtr();
      62             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
      63             : 
      64             :   // If we are on an internal face, we populate the four entries in the system matrix
      65             :   // which touch the face
      66     2821392 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
      67             :   {
      68             :     // The dof ids of the variable corresponding to the element and neighbor
      69     2083312 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
      70     2083312 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
      71             : 
      72             :     // Compute the entries which will go to the neighbor (offdiagonal) and element
      73             :     // (diagonal).
      74     2083312 :     const auto elem_matrix_contribution = computeElemMatrixContribution();
      75     2083312 :     const auto neighbor_matrix_contribution = computeNeighborMatrixContribution();
      76             : 
      77             :     // Populate matrix
      78     2083312 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()) && !(bounded_elem))
      79             :     {
      80     1870572 :       _matrix_contribution(0, 0) = elem_matrix_contribution;
      81     1870572 :       _matrix_contribution(0, 1) = neighbor_matrix_contribution;
      82             :     }
      83             : 
      84     2083312 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()) && !(bounded_neigh))
      85             :     {
      86     1872692 :       _matrix_contribution(1, 0) = -elem_matrix_contribution;
      87     1872692 :       _matrix_contribution(1, 1) = -neighbor_matrix_contribution;
      88             :     }
      89     2083312 :     (*_linear_system.matrix).add_matrix(_matrix_contribution, _dof_indices.get_values());
      90             :   }
      91             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
      92             :   // We check if we have a boundary condition here
      93      738080 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      94             :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
      95             :   {
      96             :     mooseAssert(_current_face_info->boundaryIDs().size() == 1,
      97             :                 "We should only have one boundary on every face.");
      98             : 
      99             :     LinearFVBoundaryCondition * bc_pointer =
     100      738080 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     101             : 
     102      738080 :     if (bc_pointer || _force_boundary_execution)
     103             :     {
     104             :       if (bc_pointer)
     105      413808 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     106      413808 :       const auto matrix_contribution = computeBoundaryMatrixContribution(*bc_pointer);
     107             : 
     108             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     109             :       // are on (assuming that this is a boundary for the variable)
     110      413808 :       if ((_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) && !(bounded_elem))
     111             :       {
     112      355000 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     113      355000 :         (*_linear_system.matrix).add(dof_id_elem, dof_id_elem, matrix_contribution);
     114             :       }
     115       58808 :       else if ((_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) && !(bounded_neigh))
     116             :       {
     117             :         const auto dof_id_neighbor =
     118           0 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     119           0 :         (*_linear_system.matrix).add(dof_id_neighbor, dof_id_neighbor, matrix_contribution);
     120             :       }
     121             :     }
     122             :   }
     123     2821392 : }
     124             : 
     125             : void
     126     2821392 : LinearFVTurbulentAdvection::addRightHandSideContribution()
     127             : {
     128             :   // Coumputing bounding map
     129     2821392 :   const Elem * elem = _current_face_info->elemPtr();
     130             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
     131     4936848 :   const Elem * neighbor = _current_face_info->neighborPtr();
     132             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
     133             : 
     134             :   // If we are on an internal face, we populate the two entries in the right hand side
     135             :   // which touch the face
     136     2821392 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
     137             :   {
     138             :     // The dof ids of the variable corresponding to the element and neighbor
     139     2083312 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     140     2083312 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     141             : 
     142             :     // Compute the entries which will go to the neighbor and element positions.
     143     2083312 :     const auto elem_rhs_contribution = computeElemRightHandSideContribution();
     144     2083312 :     const auto neighbor_rhs_contribution = computeNeighborRightHandSideContribution();
     145             : 
     146             :     // Populate right hand side
     147     2083312 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()) && !(bounded_elem))
     148     1870572 :       _rhs_contribution(0) = elem_rhs_contribution;
     149     2083312 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()) && !(bounded_neigh))
     150     1872692 :       _rhs_contribution(1) = neighbor_rhs_contribution;
     151             : 
     152     2083312 :     (*_linear_system.rhs)
     153     2083312 :         .add_vector(_rhs_contribution.get_values().data(), _dof_indices.get_values());
     154             :   }
     155             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
     156             :   // We check if we have a boundary condition here
     157      738080 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     158             :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     159             :   {
     160             :     mooseAssert(_current_face_info->boundaryIDs().size() == 1,
     161             :                 "We should only have one boundary on every face.");
     162             :     LinearFVBoundaryCondition * bc_pointer =
     163      738080 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     164             : 
     165      738080 :     if (bc_pointer || _force_boundary_execution)
     166             :     {
     167             :       if (bc_pointer)
     168      413808 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     169             : 
     170      413808 :       const auto rhs_contribution = computeBoundaryRHSContribution(*bc_pointer);
     171             : 
     172             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     173             :       // are on (assuming that this is a boundary for the variable)
     174      413808 :       if ((_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) && !(bounded_elem))
     175             :       {
     176      355000 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     177      355000 :         (*_linear_system.rhs).add(dof_id_elem, rhs_contribution);
     178             :       }
     179       58808 :       else if ((_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) && !(bounded_neigh))
     180             :       {
     181             :         const auto dof_id_neighbor =
     182           0 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     183           0 :         (*_linear_system.rhs).add(dof_id_neighbor, rhs_contribution);
     184             :       }
     185             :     }
     186             :   }
     187     2821392 : }
     188             : 
     189             : Real
     190     2083312 : LinearFVTurbulentAdvection::computeElemMatrixContribution()
     191             : {
     192     2083312 :   return _advected_interp_coeffs.first * _mass_face_flux * _current_face_area;
     193             : }
     194             : 
     195             : Real
     196     2083312 : LinearFVTurbulentAdvection::computeNeighborMatrixContribution()
     197             : {
     198     2083312 :   return _advected_interp_coeffs.second * _mass_face_flux * _current_face_area;
     199             : }
     200             : 
     201             : Real
     202     2083312 : LinearFVTurbulentAdvection::computeElemRightHandSideContribution()
     203             : {
     204     2083312 :   return 0.0;
     205             : }
     206             : 
     207             : Real
     208     2083312 : LinearFVTurbulentAdvection::computeNeighborRightHandSideContribution()
     209             : {
     210     2083312 :   return 0.0;
     211             : }
     212             : 
     213             : Real
     214      413808 : LinearFVTurbulentAdvection::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc)
     215             : {
     216             :   const auto * const adv_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     217             :   mooseAssert(adv_bc, "This should be a valid BC!");
     218             : 
     219      413808 :   const auto boundary_value_matrix_contrib = adv_bc->computeBoundaryValueMatrixContribution();
     220             : 
     221             :   // We support internal boundaries too so we have to make sure the normal points always outward
     222      413808 :   const auto factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0;
     223             : 
     224      413808 :   return boundary_value_matrix_contrib * factor * _mass_face_flux * _current_face_area;
     225             : }
     226             : 
     227             : Real
     228      413808 : LinearFVTurbulentAdvection::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     229             : {
     230             :   const auto * const adv_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     231             :   mooseAssert(adv_bc, "This should be a valid BC!");
     232             : 
     233             :   // We support internal boundaries too so we have to make sure the normal points always outward
     234      413808 :   const auto factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ? 1.0 : -1.0);
     235             : 
     236      413808 :   const auto boundary_value_rhs_contrib = adv_bc->computeBoundaryValueRHSContribution();
     237      413808 :   return -boundary_value_rhs_contrib * factor * _mass_face_flux * _current_face_area;
     238             : }
     239             : 
     240             : void
     241     2821392 : LinearFVTurbulentAdvection::setupFaceData(const FaceInfo * face_info)
     242             : {
     243     2821392 :   LinearFVFluxKernel::setupFaceData(face_info);
     244             : 
     245             :   // Caching the mass flux on the face which will be reused in the advection term's matrix and right
     246             :   // hand side contributions
     247     2821392 :   _mass_face_flux = _mass_flux_provider.getMassFlux(*face_info);
     248             : 
     249             :   // Caching the interpolation coefficients so they will be reused for the matrix and right hand
     250             :   // side terms
     251             :   _advected_interp_coeffs =
     252     2821392 :       interpCoeffs(_advected_interp_method, *_current_face_info, true, _mass_face_flux);
     253     2821392 : }