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 "LinearFVTurbulentDiffusion.h"
      11             : #include "Assembly.h"
      12             : #include "SubProblem.h"
      13             : #include "LinearFVAdvectionDiffusionBC.h"
      14             : #include "NavierStokesMethods.h"
      15             : #include "NS.h"
      16             : 
      17             : registerMooseObject("NavierStokesApp", LinearFVTurbulentDiffusion);
      18             : 
      19             : InputParameters
      20         424 : LinearFVTurbulentDiffusion::validParams()
      21             : {
      22         424 :   InputParameters params = LinearFVDiffusion::validParams();
      23         424 :   params.addClassDescription(
      24             :       "Represents the matrix and right hand side contributions of a "
      25             :       "diffusion term for a turbulent variable in a partial differential equation.");
      26             : 
      27         848 :   params.addParam<MooseFunctorName>(
      28         848 :       "scaling_coeff", 1.0, "The scaling coefficient for the diffusion term.");
      29             : 
      30         848 :   params.addParam<std::vector<BoundaryName>>(
      31             :       "walls", {}, "Boundaries that correspond to solid walls.");
      32         424 :   return params;
      33           0 : }
      34             : 
      35         224 : LinearFVTurbulentDiffusion::LinearFVTurbulentDiffusion(const InputParameters & params)
      36             :   : LinearFVDiffusion(params),
      37         224 :     _scaling_coeff(getFunctor<Real>("scaling_coeff")),
      38         672 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls"))
      39             : {
      40         224 :   if (_use_nonorthogonal_correction)
      41           0 :     _var.computeCellGradients();
      42         224 : }
      43             : 
      44             : void
      45        1156 : LinearFVTurbulentDiffusion::initialSetup()
      46             : {
      47        1156 :   LinearFVDiffusion::initialSetup();
      48        2312 :   NS::getWallBoundedElements(
      49        1156 :       _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
      50        1156 : }
      51             : 
      52             : Real
      53     3523552 : LinearFVTurbulentDiffusion::computeElemMatrixContribution()
      54             : {
      55     3523552 :   const auto face_arg = makeCDFace(*_current_face_info);
      56     3523552 :   return computeFluxMatrixContribution() / _scaling_coeff(face_arg, determineState());
      57             : }
      58             : 
      59             : Real
      60     3523552 : LinearFVTurbulentDiffusion::computeNeighborMatrixContribution()
      61             : {
      62     3523552 :   const auto face_arg = makeCDFace(*_current_face_info);
      63     3523552 :   return -computeFluxMatrixContribution() / _scaling_coeff(face_arg, determineState());
      64             : }
      65             : 
      66             : Real
      67     3523552 : LinearFVTurbulentDiffusion::computeElemRightHandSideContribution()
      68             : {
      69     3523552 :   const auto face_arg = makeCDFace(*_current_face_info);
      70     3523552 :   return computeFluxRHSContribution() / _scaling_coeff(face_arg, determineState());
      71             : }
      72             : 
      73             : Real
      74     3523552 : LinearFVTurbulentDiffusion::computeNeighborRightHandSideContribution()
      75             : {
      76     3523552 :   const auto face_arg = makeCDFace(*_current_face_info);
      77     3523552 :   return -computeFluxRHSContribution() / _scaling_coeff(face_arg, determineState());
      78             : }
      79             : 
      80             : Real
      81      813824 : LinearFVTurbulentDiffusion::computeBoundaryMatrixContribution(const LinearFVBoundaryCondition & bc)
      82             : {
      83             :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
      84             :   mooseAssert(diff_bc, "This should be a valid BC!");
      85             : 
      86      813824 :   auto grad_contrib = diff_bc->computeBoundaryGradientMatrixContribution() * _current_face_area;
      87             :   // If the boundary condition does not include the diffusivity contribution then
      88             :   // add it here.
      89      813824 :   if (!diff_bc->includesMaterialPropertyMultiplier())
      90             :   {
      91      404352 :     const auto face_arg = singleSidedFaceArg(_current_face_info);
      92      404352 :     grad_contrib *=
      93      404352 :         _diffusion_coeff(face_arg, determineState()) / _scaling_coeff(face_arg, determineState());
      94             :   }
      95             : 
      96      813824 :   return grad_contrib;
      97             : }
      98             : 
      99             : Real
     100      813824 : LinearFVTurbulentDiffusion::computeBoundaryRHSContribution(const LinearFVBoundaryCondition & bc)
     101             : {
     102             :   const auto * const diff_bc = static_cast<const LinearFVAdvectionDiffusionBC *>(&bc);
     103             :   mooseAssert(diff_bc, "This should be a valid BC!");
     104             : 
     105      813824 :   const auto face_arg = singleSidedFaceArg(_current_face_info);
     106      813824 :   const auto state = determineState();
     107      813824 :   auto grad_contrib = diff_bc->computeBoundaryGradientRHSContribution() * _current_face_area;
     108             : 
     109             :   // If the boundary condition does not include the diffusivity contribution then
     110             :   // add it here.
     111      813824 :   if (!diff_bc->includesMaterialPropertyMultiplier())
     112      404352 :     grad_contrib *= _diffusion_coeff(face_arg, state) / _scaling_coeff(face_arg, state);
     113             : 
     114             :   // We add the nonorthogonal corrector for the face here. Potential idea: we could do
     115             :   // this in the boundary condition too. For now, however, we keep it like this.
     116      813824 :   if (_use_nonorthogonal_correction)
     117             :   {
     118             :     const auto correction_vector =
     119             :         _current_face_info->normal() -
     120           0 :         1 / (_current_face_info->normal() * _current_face_info->eCN()) * _current_face_info->eCN();
     121             : 
     122           0 :     grad_contrib += _diffusion_coeff(face_arg, state) *
     123           0 :                     _var.gradSln(*_current_face_info->elemInfo(), state) * correction_vector *
     124           0 :                     _current_face_area;
     125             :   }
     126             : 
     127      813824 :   return grad_contrib;
     128             : }
     129             : 
     130             : void
     131     4912160 : LinearFVTurbulentDiffusion::addMatrixContribution()
     132             : {
     133             :   // Coumputing bounding map
     134     4912160 :   const Elem * elem = _current_face_info->elemPtr();
     135             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
     136     8478336 :   const Elem * neighbor = _current_face_info->neighborPtr();
     137             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
     138             : 
     139             :   // If we are on an internal face, we populate the four entries in the system matrix
     140             :   // which touch the face
     141     4912160 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
     142             :   {
     143             :     // The dof ids of the variable corresponding to the element and neighbor
     144     3523552 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     145     3523552 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     146             : 
     147             :     // Compute the entries which will go to the neighbor (offdiagonal) and element
     148             :     // (diagonal).
     149     3523552 :     const auto elem_matrix_contribution = computeElemMatrixContribution();
     150     3523552 :     const auto neighbor_matrix_contribution = computeNeighborMatrixContribution();
     151             : 
     152             :     // Populate matrix
     153     3523552 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()) && !(bounded_elem))
     154             :     {
     155     3169344 :       _matrix_contribution(0, 0) = elem_matrix_contribution;
     156     3169344 :       _matrix_contribution(0, 1) = neighbor_matrix_contribution;
     157             :     }
     158             : 
     159     3523552 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()) && !(bounded_neigh))
     160             :     {
     161     3162080 :       _matrix_contribution(1, 0) = -elem_matrix_contribution;
     162     3162080 :       _matrix_contribution(1, 1) = -neighbor_matrix_contribution;
     163             :     }
     164     3523552 :     (*_linear_system.matrix).add_matrix(_matrix_contribution, _dof_indices.get_values());
     165             :   }
     166             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
     167             :   // We check if we have a boundary condition here
     168     1388608 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     169             :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     170             :   {
     171             :     mooseAssert(_current_face_info->boundaryIDs().size() == 1,
     172             :                 "We should only have one boundary on every face.");
     173             : 
     174             :     LinearFVBoundaryCondition * bc_pointer =
     175     1388608 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     176             : 
     177     1388608 :     if (bc_pointer || _force_boundary_execution)
     178             :     {
     179             :       if (bc_pointer)
     180      813824 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     181      813824 :       const auto matrix_contribution = computeBoundaryMatrixContribution(*bc_pointer);
     182             : 
     183             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     184             :       // are on (assuming that this is a boundary for the variable)
     185      813824 :       if ((_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) && !(bounded_elem))
     186             :       {
     187      701920 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     188      701920 :         (*_linear_system.matrix).add(dof_id_elem, dof_id_elem, matrix_contribution);
     189             :       }
     190      111904 :       else if ((_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) && !(bounded_neigh))
     191             :       {
     192             :         const auto dof_id_neighbor =
     193           0 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     194           0 :         (*_linear_system.matrix).add(dof_id_neighbor, dof_id_neighbor, matrix_contribution);
     195             :       }
     196             :     }
     197             :   }
     198     4912160 : }
     199             : 
     200             : void
     201     4912160 : LinearFVTurbulentDiffusion::addRightHandSideContribution()
     202             : {
     203             :   // Coumputing bounding map
     204     4912160 :   const Elem * elem = _current_face_info->elemPtr();
     205             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
     206     8478336 :   const Elem * neighbor = _current_face_info->neighborPtr();
     207             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
     208             : 
     209             :   // If we are on an internal face, we populate the two entries in the right hand side
     210             :   // which touch the face
     211     4912160 :   if (_current_face_type == FaceInfo::VarFaceNeighbors::BOTH)
     212             :   {
     213             :     // The dof ids of the variable corresponding to the element and neighbor
     214     3523552 :     _dof_indices(0) = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     215     3523552 :     _dof_indices(1) = _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     216             : 
     217             :     // Compute the entries which will go to the neighbor and element positions.
     218     3523552 :     const auto elem_rhs_contribution = computeElemRightHandSideContribution();
     219     3523552 :     const auto neighbor_rhs_contribution = computeNeighborRightHandSideContribution();
     220             : 
     221             :     // Populate right hand side
     222     3523552 :     if (hasBlocks(_current_face_info->elemInfo()->subdomain_id()))
     223     3523552 :       _rhs_contribution(0) = elem_rhs_contribution;
     224     3523552 :     if (hasBlocks(_current_face_info->neighborInfo()->subdomain_id()))
     225     3523552 :       _rhs_contribution(1) = neighbor_rhs_contribution;
     226             : 
     227     3523552 :     (*_linear_system.rhs)
     228     3523552 :         .add_vector(_rhs_contribution.get_values().data(), _dof_indices.get_values());
     229             :   }
     230             :   // We are at a block boundary where the variable is not defined on one of the adjacent cells.
     231             :   // We check if we have a boundary condition here
     232     1388608 :   else if (_current_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     233             :            _current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR)
     234             :   {
     235             :     mooseAssert(_current_face_info->boundaryIDs().size() == 1,
     236             :                 "We should only have one boundary on every face.");
     237             :     LinearFVBoundaryCondition * bc_pointer =
     238     1388608 :         _var.getBoundaryCondition(*_current_face_info->boundaryIDs().begin());
     239             : 
     240     1388608 :     if (bc_pointer || _force_boundary_execution)
     241             :     {
     242             :       if (bc_pointer)
     243      813824 :         bc_pointer->setupFaceData(_current_face_info, _current_face_type);
     244             : 
     245      813824 :       const auto rhs_contribution = computeBoundaryRHSContribution(*bc_pointer);
     246             : 
     247             :       // We allow internal (for the mesh) boundaries too, so we have to check on which side we
     248             :       // are on (assuming that this is a boundary for the variable)
     249      813824 :       if ((_current_face_type == FaceInfo::VarFaceNeighbors::ELEM) && !(bounded_elem))
     250             :       {
     251      701920 :         const auto dof_id_elem = _current_face_info->elemInfo()->dofIndices()[_sys_num][_var_num];
     252      701920 :         (*_linear_system.rhs).add(dof_id_elem, rhs_contribution);
     253             :       }
     254      111904 :       else if ((_current_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) && !(bounded_neigh))
     255             :       {
     256             :         const auto dof_id_neighbor =
     257           0 :             _current_face_info->neighborInfo()->dofIndices()[_sys_num][_var_num];
     258           0 :         (*_linear_system.rhs).add(dof_id_neighbor, rhs_contribution);
     259             :       }
     260             :     }
     261             :   }
     262     4912160 : }