Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       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 "INSFVTurbulentAdvection.h"
      11             : #include "NavierStokesMethods.h"
      12             : #include "NS.h"
      13             : 
      14             : registerMooseObject("NavierStokesApp", INSFVTurbulentAdvection);
      15             : 
      16             : InputParameters
      17         914 : INSFVTurbulentAdvection::validParams()
      18             : {
      19         914 :   auto params = INSFVAdvectionKernel::validParams();
      20         914 :   params.addClassDescription(
      21             :       "Advects an arbitrary turbulent quantity, the associated nonlinear 'variable'.");
      22         914 :   params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
      23        1828 :   params.addParam<std::vector<BoundaryName>>(
      24             :       "walls", {}, "Boundaries that correspond to solid walls.");
      25        1828 :   params.addParam<bool>("neglect_advection_derivatives",
      26        1828 :                         true,
      27             :                         "Whether to remove automatic differentiation derivative terms "
      28             :                         "for velocity in the advection term");
      29         914 :   return params;
      30           0 : }
      31             : 
      32         514 : INSFVTurbulentAdvection::INSFVTurbulentAdvection(const InputParameters & params)
      33             :   : INSFVAdvectionKernel(params),
      34        1028 :     _rho(getFunctor<ADReal>(NS::density)),
      35        1028 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
      36        1542 :     _neglect_advection_derivatives(getParam<bool>("neglect_advection_derivatives"))
      37             : {
      38         514 : }
      39             : 
      40             : void
      41        1858 : INSFVTurbulentAdvection::initialSetup()
      42             : {
      43        1858 :   INSFVAdvectionKernel::initialSetup();
      44        3716 :   NS::getWallBoundedElements(
      45        1858 :       _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
      46        1858 : }
      47             : 
      48             : ADReal
      49     1320064 : INSFVTurbulentAdvection::computeQpResidual()
      50             : {
      51     1320064 :   const auto v = _rc_vel_provider.getVelocity(
      52     1320064 :       _velocity_interp_method, *_face_info, determineState(), _tid, false);
      53     1320064 :   const auto var_face = _var(makeFace(*_face_info,
      54             :                                       limiterType(_advected_interp_method),
      55     1320064 :                                       MetaPhysicL::raw_value(v) * _normal > 0),
      56     2640128 :                              determineState());
      57     1320064 :   const auto rho_face = _rho(makeFace(*_face_info,
      58             :                                       limiterType(_advected_interp_method),
      59     1320064 :                                       MetaPhysicL::raw_value(v) * _normal > 0),
      60     2640128 :                              determineState());
      61     1320064 :   if (!_neglect_advection_derivatives)
      62       51712 :     return _normal * v * rho_face * var_face;
      63             :   else
      64     2536704 :     return _normal * MetaPhysicL::raw_value(v) * rho_face * var_face;
      65             : }
      66             : 
      67             : void
      68      444592 : INSFVTurbulentAdvection::computeResidual(const FaceInfo & fi)
      69             : {
      70      444592 :   if (skipForBoundary(fi))
      71       74528 :     return;
      72             : 
      73      370064 :   _face_info = &fi;
      74      370064 :   _normal = fi.normal();
      75      370064 :   _face_type = _face_info->faceType(std::make_pair(_var.number(), _var.sys().number()));
      76      740128 :   auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual());
      77             : 
      78      370064 :   const Elem * elem = fi.elemPtr();
      79      370064 :   const Elem * neighbor = fi.neighborPtr();
      80             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
      81             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
      82             : 
      83      370064 :   if ((_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      84      370064 :        _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
      85             :       (!bounded_elem))
      86             :   {
      87             :     // residual contribution of this kernel to the elem element
      88      341152 :     prepareVectorTag(_assembly, _var.number());
      89      341152 :     _local_re(0) = r;
      90      341152 :     accumulateTaggedLocalResidual();
      91             :   }
      92      370064 :   if ((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
      93      370064 :        _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
      94             :       (!bounded_neigh))
      95             :   {
      96             :     // residual contribution of this kernel to the neighbor element
      97      341152 :     prepareVectorTagNeighbor(_assembly, _var.number());
      98      341152 :     _local_re(0) = -r;
      99      341152 :     accumulateTaggedLocalResidual();
     100             :   }
     101             : }
     102             : 
     103             : void
     104     1121224 : INSFVTurbulentAdvection::computeJacobian(const FaceInfo & fi)
     105             : {
     106     1121224 :   if (skipForBoundary(fi))
     107      171224 :     return;
     108             : 
     109      950000 :   _face_info = &fi;
     110      950000 :   _normal = fi.normal();
     111      950000 :   _face_type = _face_info->faceType(std::make_pair(_var.number(), _var.sys().number()));
     112     1900000 :   const ADReal r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
     113             : 
     114      950000 :   const Elem * elem = fi.elemPtr();
     115      950000 :   const Elem * neighbor = fi.neighborPtr();
     116             :   const auto bounded_elem = _wall_bounded.find(elem) != _wall_bounded.end();
     117             :   const auto bounded_neigh = _wall_bounded.find(neighbor) != _wall_bounded.end();
     118             : 
     119      950000 :   if ((_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     120      950000 :        _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
     121             :       (!bounded_elem))
     122             :   {
     123             :     mooseAssert(_var.dofIndices().size() == 1, "We're currently built to use CONSTANT MONOMIALS");
     124             : 
     125      850724 :     addResidualsAndJacobian(
     126     1701448 :         _assembly, std::array<ADReal, 1>{{r}}, _var.dofIndices(), _var.scalingFactor());
     127             :   }
     128             : 
     129      950000 :   if ((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
     130      874904 :        _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
     131             :       (!bounded_neigh))
     132             :   {
     133             :     mooseAssert((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR) ==
     134             :                     (_var.dofIndices().size() == 0),
     135             :                 "If the variable is only defined on the neighbor hand side of the face, then that "
     136             :                 "means it should have no dof indices on the elem element. Conversely if "
     137             :                 "the variable is defined on both sides of the face, then it should have a non-zero "
     138             :                 "number of degrees of freedom on the elem element");
     139             : 
     140             :     // We switch the sign for the neighbor residual
     141      788144 :     ADReal neighbor_r = -r;
     142             : 
     143             :     mooseAssert(_var.dofIndicesNeighbor().size() == 1,
     144             :                 "We're currently built to use CONSTANT MONOMIALS");
     145             : 
     146      788144 :     addResidualsAndJacobian(_assembly,
     147     1576288 :                             std::array<ADReal, 1>{{neighbor_r}},
     148             :                             _var.dofIndicesNeighbor(),
     149      788144 :                             _var.scalingFactor());
     150             :   }
     151             : }