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