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 "NavierStokesMethods.h"
      13             : #include "NS.h"
      14             : 
      15             : registerMooseObject("NavierStokesApp", LinearFVTurbulentAdvection);
      16             : 
      17             : InputParameters
      18         302 : LinearFVTurbulentAdvection::validParams()
      19             : {
      20         302 :   InputParameters params = LinearFVFluxKernel::validParams();
      21         302 :   params.addClassDescription("Represents the matrix and right hand side contributions of an "
      22             :                              "advection term for a turbulence variable.");
      23             : 
      24         604 :   params.addRequiredParam<UserObjectName>(
      25             :       "rhie_chow_user_object",
      26             :       "The rhie-chow user-object which is used to determine the face velocity.");
      27             : 
      28         604 :   params.addRequiredParam<InterpolationMethodName>(
      29             :       "advected_interp_method_name",
      30             :       "Name of the FVInterpolationMethod to use for the advected quantity.");
      31             : 
      32         604 :   params.addParam<std::vector<BoundaryName>>(
      33             :       "walls", {}, "Boundaries that correspond to solid walls.");
      34             : 
      35         302 :   return params;
      36           0 : }
      37             : 
      38         166 : LinearFVTurbulentAdvection::LinearFVTurbulentAdvection(const InputParameters & params)
      39             :   : LinearFVFluxKernel(params),
      40             :     FVInterpolationMethodInterface(this),
      41         166 :     _mass_flux_provider(getUserObject<RhieChowMassFlux>("rhie_chow_user_object")),
      42         332 :     _adv_interp_method(getFVAdvectedInterpolationMethod(
      43             :         getParam<InterpolationMethodName>("advected_interp_method_name"))),
      44         166 :     _mass_face_flux(0.0),
      45         498 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls"))
      46             : {
      47         166 :   if (_adv_interp_method.needsGradients())
      48          24 :     _var.computeCellGradients(_adv_interp_method.gradientLimiter());
      49         166 : }
      50             : 
      51             : void
      52         866 : LinearFVTurbulentAdvection::initialSetup()
      53             : {
      54         866 :   LinearFVFluxKernel::initialSetup();
      55        1732 :   NS::getWallBoundedElements(
      56         866 :       _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
      57         866 : }
      58             : 
      59             : void
      60     3356984 : LinearFVTurbulentAdvection::addMatrixContribution()
      61             : {
      62             :   // Coumputing bounding map
      63     3356984 :   const Elem * elem = _current_face_info->elemPtr();
      64             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
      65     5756496 :   const Elem * neighbor = _current_face_info->neighborPtr();
      66             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
      67             : 
      68             :   // If we are on an internal face, we populate the four entries in the system matrix
      69             :   // which touch the face
      70     3356984 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
      71             :   {
      72             :     // The dof ids of the variable corresponding to the element and neighbor
      73     2378184 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
      74     2378184 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
      75             : 
      76             :     // Compute the entries which will go to the neighbor (offdiagonal) and element
      77             :     // (diagonal).
      78     2378184 :     const auto elem_matrix_contribution = computeElemMatrixContribution();
      79     2378184 :     const auto neighbor_matrix_contribution = computeNeighborMatrixContribution();
      80             : 
      81             :     // Populate matrix
      82     2378184 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()) && !(bounded_elem))
      83             :     {
      84     2141810 :       _matrix_contribution(0, 0) = elem_matrix_contribution;
      85     2141810 :       _matrix_contribution(0, 1) = neighbor_matrix_contribution;
      86             :     }
      87             : 
      88     2378184 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()) && !(bounded_neigh))
      89             :     {
      90     2138178 :       _matrix_contribution(1, 0) = -elem_matrix_contribution;
      91     2138178 :       _matrix_contribution(1, 1) = -neighbor_matrix_contribution;
      92             :     }
      93     2378184 :     (*_linear_system.matrix).add_matrix(_matrix_contribution, _dof_indices.get_values());
      94             :   }
      95             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
      96             :   // We check if we have a boundary condition here
      97      978800 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      98             :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
      99             :   {
     100             :     mooseAssert(_current_face_info->boundaryIDs().size() == 1,
     101             :                 "We should only have one boundary on every face.");
     102             : 
     103             :     LinearFVBoundaryCondition * bc_pointer =
     104      978800 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     105             : 
     106      978800 :     if (bc_pointer || _force_boundary_execution)
     107             :     {
     108             :       if (bc_pointer)
     109      596576 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     110      596576 :       const auto matrix_contribution = computeBoundaryMatrixContribution(*bc_pointer);
     111             : 
     112             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     113             :       // are on (assuming that this is a boundary for the variable)
     114      596576 :       if ((_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) && !(bounded_elem))
     115             :       {
     116      516916 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     117      516916 :         (*_linear_system.matrix).add(dof_id_elem, dof_id_elem, matrix_contribution);
     118             :       }
     119       79660 :       else if ((_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) && !(bounded_neigh))
     120             :       {
     121             :         const auto dof_id_neighbor =
     122           0 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     123           0 :         (*_linear_system.matrix).add(dof_id_neighbor, dof_id_neighbor, matrix_contribution);
     124             :       }
     125             :     }
     126             :   }
     127     3356984 : }
     128             : 
     129             : void
     130     3356984 : LinearFVTurbulentAdvection::addRightHandSideContribution()
     131             : {
     132             :   // Coumputing bounding map
     133     3356984 :   const Elem * elem = _current_face_info->elemPtr();
     134             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
     135     5756496 :   const Elem * neighbor = _current_face_info->neighborPtr();
     136             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
     137             : 
     138             :   // If we are on an internal face, we populate the two entries in the right hand side
     139             :   // which touch the face
     140     3356984 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
     141             :   {
     142             :     // The dof ids of the variable corresponding to the element and neighbor
     143     2378184 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     144     2378184 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     145             : 
     146             :     // Compute the entries which will go to the neighbor and element positions.
     147     2378184 :     const auto elem_rhs_contribution = computeElemRightHandSideContribution();
     148     2378184 :     const auto neighbor_rhs_contribution = computeNeighborRightHandSideContribution();
     149             : 
     150             :     // Populate right hand side
     151     2378184 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()) && !(bounded_elem))
     152     2141810 :       _rhs_contribution(0) = elem_rhs_contribution;
     153     2378184 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()) && !(bounded_neigh))
     154     2138178 :       _rhs_contribution(1) = neighbor_rhs_contribution;
     155             : 
     156     2378184 :     (*_linear_system.rhs)
     157     2378184 :         .add_vector(_rhs_contribution.get_values().data(), _dof_indices.get_values());
     158             :   }
     159             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
     160             :   // We check if we have a boundary condition here
     161      978800 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     162             :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     163             :   {
     164             :     mooseAssert(_current_face_info->boundaryIDs().size() == 1,
     165             :                 "We should only have one boundary on every face.");
     166             :     LinearFVBoundaryCondition * bc_pointer =
     167      978800 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     168             : 
     169      978800 :     if (bc_pointer || _force_boundary_execution)
     170             :     {
     171             :       if (bc_pointer)
     172      596576 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     173             : 
     174      596576 :       const auto rhs_contribution = computeBoundaryRHSContribution(*bc_pointer);
     175             : 
     176             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     177             :       // are on (assuming that this is a boundary for the variable)
     178      596576 :       if ((_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) && !(bounded_elem))
     179             :       {
     180      516916 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     181      516916 :         (*_linear_system.rhs).add(dof_id_elem, rhs_contribution);
     182             :       }
     183       79660 :       else if ((_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) && !(bounded_neigh))
     184             :       {
     185             :         const auto dof_id_neighbor =
     186           0 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     187           0 :         (*_linear_system.rhs).add(dof_id_neighbor, rhs_contribution);
     188             :       }
     189             :     }
     190             :   }
     191     3356984 : }
     192             : 
     193             : Real
     194     2378184 : LinearFVTurbulentAdvection::computeElemMatrixContribution()
     195             : {
     196             :   const auto & coeffs = _adv_interp_result.weights_matrix;
     197     2378184 :   return coeffs.first * _mass_face_flux * _current_face_area;
     198             : }
     199             : 
     200             : Real
     201     2378184 : LinearFVTurbulentAdvection::computeNeighborMatrixContribution()
     202             : {
     203             :   const auto & coeffs = _adv_interp_result.weights_matrix;
     204     2378184 :   return coeffs.second * _mass_face_flux * _current_face_area;
     205             : }
     206             : 
     207             : Real
     208     2378184 : LinearFVTurbulentAdvection::computeElemRightHandSideContribution()
     209             : {
     210     2378184 :   return _adv_interp_result.rhs_face_value * _mass_face_flux * _current_face_area;
     211             : }
     212             : 
     213             : Real
     214     2378184 : LinearFVTurbulentAdvection::computeNeighborRightHandSideContribution()
     215             : {
     216     2378184 :   return -_adv_interp_result.rhs_face_value * _mass_face_flux * _current_face_area;
     217             : }
     218             : 
     219             : Real
     220      596576 : LinearFVTurbulentAdvection::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc)
     221             : {
     222             :   const auto * const adv_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     223             :   mooseAssert(adv_bc, "This should be a valid BC!");
     224             : 
     225      596576 :   const auto boundary_value_matrix_contrib = adv_bc->computeBoundaryValueMatrixContribution();
     226             : 
     227             :   // We support internal boundaries too so we have to make sure the normal points always outward
     228      596576 :   const auto factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) ? 1.0 : -1.0;
     229             : 
     230      596576 :   return boundary_value_matrix_contrib * factor * _mass_face_flux * _current_face_area;
     231             : }
     232             : 
     233             : Real
     234      596576 : LinearFVTurbulentAdvection::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     235             : {
     236             :   const auto * const adv_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     237             :   mooseAssert(adv_bc, "This should be a valid BC!");
     238             : 
     239             :   // We support internal boundaries too so we have to make sure the normal points always outward
     240      596576 :   const auto factor = (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ? 1.0 : -1.0);
     241             : 
     242      596576 :   const auto boundary_value_rhs_contrib = adv_bc->computeBoundaryValueRHSContribution();
     243      596576 :   return -boundary_value_rhs_contrib * factor * _mass_face_flux * _current_face_area;
     244             : }
     245             : 
     246             : void
     247     3356984 : LinearFVTurbulentAdvection::setupFaceData(const FaceInfo * face_info)
     248             : {
     249     3356984 :   LinearFVFluxKernel::setupFaceData(face_info);
     250             : 
     251             :   // Caching the mass flux on the face which will be reused in the advection term's matrix and right
     252             :   // hand side contributions
     253     3356984 :   _mass_face_flux = _mass_flux_provider.getMassFlux(*face_info);
     254             : 
     255             :   // Only internal faces need advected interpolation results; boundary contributions are handled
     256             :   // through the linear FV boundary conditions.
     257     3356984 :   if (_current_face_type != FaceInfo::VarFaceNeighbors::BOTH)
     258      978800 :     return;
     259             : 
     260     2378184 :   const auto state = determineState();
     261     2378184 :   const auto & elem_info = *_current_face_info->elemInfo();
     262             :   const auto & neighbor_info = *_current_face_info->neighborInfo();
     263             : 
     264     2378184 :   const Real elem_value = _var.getElemValue(elem_info, state);
     265     2378184 :   const Real neighbor_value = _var.getElemValue(neighbor_info, state);
     266     2378184 :   if (_adv_interp_method.needsGradients())
     267             :   {
     268      269984 :     const auto limiter_type = _adv_interp_method.gradientLimiter();
     269      269984 :     _elem_grad_storage = _var.gradSln(elem_info, state, limiter_type);
     270      269984 :     _neighbor_grad_storage = _var.gradSln(neighbor_info, state, limiter_type);
     271             :   }
     272             : 
     273     2378184 :   _adv_interp_result = _adv_interp_method.advectedInterpolate(*_current_face_info,
     274             :                                                               elem_value,
     275             :                                                               neighbor_value,
     276     2378184 :                                                               &_elem_grad_storage,
     277     2378184 :                                                               &_neighbor_grad_storage,
     278             :                                                               _mass_face_flux);
     279             : }