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        1888 : INSFVTurbulentAdvection::validParams()
      18             : {
      19        1888 :   auto params = INSFVAdvectionKernel::validParams();
      20        1888 :   params.addClassDescription(
      21             :       "Advects an arbitrary turbulent quantity, the associated nonlinear 'variable'.");
      22        1888 :   params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
      23        3776 :   params.addParam<std::vector<BoundaryName>>(
      24             :       "walls", {}, "Boundaries that correspond to solid walls.");
      25        3776 :   params.addParam<bool>("neglect_advection_derivatives",
      26        3776 :                         true,
      27             :                         "Whether to remove automatic differentiation derivative terms "
      28             :                         "for velocity in the advection term");
      29        1888 :   return params;
      30           0 : }
      31             : 
      32        1082 : INSFVTurbulentAdvection::INSFVTurbulentAdvection(const InputParameters & params)
      33             :   : INSFVAdvectionKernel(params),
      34        2164 :     _rho(getFunctor<ADReal>(NS::density)),
      35        2164 :     _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
      36        3246 :     _neglect_advection_derivatives(getParam<bool>("neglect_advection_derivatives"))
      37             : {
      38        1082 : }
      39             : 
      40             : void
      41        3994 : INSFVTurbulentAdvection::initialSetup()
      42             : {
      43        3994 :   INSFVAdvectionKernel::initialSetup();
      44        7988 :   NS::getWallBoundedElements(
      45        3994 :       _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
      46        3994 : }
      47             : 
      48             : ADReal
      49     1816800 : INSFVTurbulentAdvection::computeQpResidual()
      50             : {
      51     1816800 :   const auto v = _rc_vel_provider.getVelocity(
      52     1816800 :       _velocity_interp_method, *_face_info, determineState(), _tid, false);
      53     1816800 :   const auto var_face = _var(makeFace(*_face_info,
      54             :                                       limiterType(_advected_interp_method),
      55     1816800 :                                       MetaPhysicL::raw_value(v) * _normal > 0),
      56     3633600 :                              determineState());
      57     1816800 :   const auto rho_face = _rho(makeFace(*_face_info,
      58             :                                       limiterType(_advected_interp_method),
      59     1816800 :                                       MetaPhysicL::raw_value(v) * _normal > 0),
      60     3633600 :                              determineState());
      61     1816800 :   if (!_neglect_advection_derivatives)
      62      154624 :     return _normal * v * rho_face * var_face;
      63             :   else
      64     3478976 :     return _normal * MetaPhysicL::raw_value(v) * rho_face * var_face;
      65             : }
      66             : 
      67             : void
      68      498880 : INSFVTurbulentAdvection::computeResidual(const FaceInfo & fi)
      69             : {
      70      498880 :   if (skipForBoundary(fi))
      71       85440 :     return;
      72             : 
      73      413440 :   _face_info = &fi;
      74      413440 :   _normal = fi.normal();
      75      413440 :   _face_type = _face_info->faceType(std::make_pair(_var.number(), _var.sys().number()));
      76      826880 :   auto r = MetaPhysicL::raw_value(fi.faceArea() * fi.faceCoord() * computeQpResidual());
      77             : 
      78      413440 :   const Elem * elem = fi.elemPtr();
      79      413440 :   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      413440 :   if ((_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      84      413440 :        _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
      85             :       (!bounded_elem))
      86             :   {
      87             :     // residual contribution of this kernel to the elem element
      88      386900 :     prepareVectorTag(_assembly, _var.number());
      89      386900 :     _local_re(0) = r;
      90      386900 :     accumulateTaggedLocalResidual();
      91             :   }
      92      413440 :   if ((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
      93      413440 :        _face_type == FaceInfo::VarFaceNeighbors::BOTH) &&
      94             :       (!bounded_neigh))
      95             :   {
      96             :     // residual contribution of this kernel to the neighbor element
      97      386900 :     prepareVectorTagNeighbor(_assembly, _var.number());
      98      386900 :     _local_re(0) = -r;
      99      386900 :     accumulateTaggedLocalResidual();
     100             :   }
     101             : }
     102             : 
     103             : void
     104     1656976 : INSFVTurbulentAdvection::computeJacobian(const FaceInfo & fi)
     105             : {
     106     1656976 :   if (skipForBoundary(fi))
     107      253616 :     return;
     108             : 
     109     1403360 :   _face_info = &fi;
     110     1403360 :   _normal = fi.normal();
     111     1403360 :   _face_type = _face_info->faceType(std::make_pair(_var.number(), _var.sys().number()));
     112     2806720 :   const ADReal r = fi.faceArea() * fi.faceCoord() * computeQpResidual();
     113             : 
     114     1403360 :   const Elem * elem = fi.elemPtr();
     115     1403360 :   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     1403360 :   if ((_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
     120     1403360 :        _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     1256236 :     addResidualsAndJacobian(
     126     2512472 :         _assembly, std::array<ADReal, 1>{{r}}, _var.dofIndices(), _var.scalingFactor());
     127             :   }
     128             : 
     129     1403360 :   if ((_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
     130     1291376 :        _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     1162916 :     ADReal neighbor_r = -r;
     142             : 
     143             :     mooseAssert(_var.dofIndicesNeighbor().size() == 1,
     144             :                 "We're currently built to use CONSTANT MONOMIALS");
     145             : 
     146     1162916 :     addResidualsAndJacobian(_assembly,
     147     2325832 :                             std::array<ADReal, 1>{{neighbor_r}},
     148             :                             _var.dofIndicesNeighbor(),
     149     1162916 :                             _var.scalingFactor());
     150             :   }
     151             : }