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 "PINSFVMomentumDiffusion.h"
      11             : #include "PINSFVSuperficialVelocityVariable.h"
      12             : #include "NS.h"
      13             : #include "INSFVRhieChowInterpolator.h"
      14             : #include "SystemBase.h"
      15             : 
      16             : registerMooseObject("NavierStokesApp", PINSFVMomentumDiffusion);
      17             : 
      18             : InputParameters
      19        9858 : PINSFVMomentumDiffusion::validParams()
      20             : {
      21        9858 :   auto params = INSFVMomentumDiffusion::validParams();
      22        9858 :   params.addClassDescription(
      23             :       "Viscous diffusion term, div(mu eps grad(u_d / eps)), in the porous media "
      24             :       "incompressible Navier-Stokes momentum equation.");
      25        9858 :   params.addRequiredParam<MooseFunctorName>(NS::porosity, "Porosity auxiliary variable");
      26        9858 :   return params;
      27           0 : }
      28             : 
      29        5575 : PINSFVMomentumDiffusion::PINSFVMomentumDiffusion(const InputParameters & params)
      30        5575 :   : INSFVMomentumDiffusion(params), _eps(getFunctor<ADReal>(NS::porosity))
      31             : {
      32        5575 :   if (!dynamic_cast<PINSFVSuperficialVelocityVariable *>(&_var))
      33           0 :     mooseError("PINSFVMomentumDiffusion may only be used with a superficial velocity "
      34             :                "variable, of variable type PINSFVSuperficialVelocityVariable.");
      35        5575 : }
      36             : 
      37             : ADReal
      38    39953472 : PINSFVMomentumDiffusion::computeStrongResidual(const bool populate_a_coeffs)
      39             : {
      40             :   using namespace Moose::FV;
      41             : 
      42    39953472 :   const bool has_elem = (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      43    39953472 :                          _face_type == FaceInfo::VarFaceNeighbors::BOTH);
      44    39953472 :   const bool has_neighbor = (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
      45    39953472 :                              _face_type == FaceInfo::VarFaceNeighbors::BOTH);
      46             : 
      47    39953472 :   const auto elem_face = elemArg();
      48    39953472 :   const auto neighbor_face = neighborArg();
      49    39953472 :   const auto state = determineState();
      50             : 
      51             :   // Compute the diffusion driven by the velocity gradient
      52             :   // Interpolate viscosity divided by porosity on the face
      53             :   ADReal mu_face;
      54             : 
      55    39953472 :   const auto mu_elem = has_elem ? _mu(elem_face, state) : _mu(neighbor_face, state);
      56    39953472 :   const auto eps_elem = has_elem ? _eps(elem_face, state) : _eps(neighbor_face, state);
      57             : 
      58             :   ADReal mu_neighbor;
      59             :   ADReal eps_neighbor;
      60    39953472 :   if (onBoundary(*_face_info))
      61     1456266 :     mu_face = _mu(singleSidedFaceArg(), state);
      62             :   else
      63             :   {
      64    38497206 :     mu_neighbor = has_elem ? _mu(neighbor_face, state) : _mu(elem_face, state);
      65    38497206 :     eps_neighbor = has_neighbor ? _eps(neighbor_face, state) : _eps(elem_face, state);
      66    38497206 :     interpolate(Moose::FV::InterpMethod::Average, mu_face, mu_elem, mu_neighbor, *_face_info, true);
      67             :   }
      68             : 
      69             :   // Compute face superficial velocity gradient
      70    79906944 :   auto dudn = _var.gradient(makeCDFace(*_face_info), state) * _face_info->normal();
      71             : 
      72    39953472 :   if (populate_a_coeffs)
      73             :   {
      74    37658052 :     if (has_elem)
      75             :     {
      76    37657932 :       const auto dof_number = _face_info->elem().dof_number(_sys.number(), _var.number(), 0);
      77             :       // A gradient is a linear combination of degrees of freedom so it's safe to straight-up index
      78             :       // into the derivatives vector at the dof we care about
      79    37657932 :       _ae = dudn.derivatives()[dof_number];
      80    37657932 :       _ae *= -mu_face;
      81             :     }
      82    37658052 :     if (has_neighbor)
      83             :     {
      84    36304940 :       const auto dof_number = _face_info->neighbor().dof_number(_sys.number(), _var.number(), 0);
      85    36304940 :       _an = dudn.derivatives()[dof_number];
      86    36304940 :       _an *= mu_face;
      87             :     }
      88             :   }
      89             : 
      90             :   // First term of residual
      91             :   ADReal residual = mu_face * dudn;
      92             : 
      93             :   // Get the face porosity gradient separately
      94             :   const auto & grad_eps_face =
      95             :       (has_elem && has_neighbor)
      96    79906944 :           ? MetaPhysicL::raw_value(_eps.gradient(makeCDFace(*_face_info), state))
      97    39953472 :           : MetaPhysicL::raw_value(_eps.gradient(
      98    41411846 :                 makeElemArg(has_elem ? &_face_info->elem() : _face_info->neighborPtr()), state));
      99             : 
     100             :   // Interpolate to get the face value
     101             :   ADReal coeff_face;
     102             :   // At this point, we already computed mu_elem/eps_elem by knowing which element owns the
     103             :   // face, so it is enough to switch between the variable evaluation here
     104             :   const auto coeff_one_side =
     105    39953472 :       mu_elem / eps_elem * (has_elem ? _var(elem_face, state) : _var(neighbor_face, state));
     106    39953472 :   if (onBoundary(*_face_info))
     107     1456266 :     coeff_face = coeff_one_side;
     108             :   else
     109             :   {
     110             :     mooseAssert(has_elem, "We should be defined on the element side if we're not on a boundary");
     111    38497206 :     const auto coeff_neighbor = mu_neighbor / eps_neighbor *
     112    38497206 :                                 (has_elem ? _var(neighbor_face, state) : _var(elem_face, state));
     113    38497206 :     interpolate(Moose::FV::InterpMethod::Average,
     114             :                 coeff_face,
     115             :                 coeff_one_side,
     116             :                 coeff_neighbor,
     117    38497206 :                 *_face_info,
     118             :                 true);
     119             :   }
     120             : 
     121    39953472 :   residual -= coeff_face * grad_eps_face * _normal;
     122             : 
     123    39953472 :   return -residual;
     124             : }