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        6197 : PINSFVMomentumDiffusion::validParams()
      20             : {
      21        6197 :   auto params = INSFVMomentumDiffusion::validParams();
      22        6197 :   params.addClassDescription(
      23             :       "Viscous diffusion term, div(mu eps grad(u_d / eps)), in the porous media "
      24             :       "incompressible Navier-Stokes momentum equation.");
      25        6197 :   params.addRequiredParam<MooseFunctorName>(NS::porosity, "Porosity auxiliary variable");
      26        6197 :   return params;
      27           0 : }
      28             : 
      29        3399 : PINSFVMomentumDiffusion::PINSFVMomentumDiffusion(const InputParameters & params)
      30        3399 :   : INSFVMomentumDiffusion(params), _eps(getFunctor<ADReal>(NS::porosity))
      31             : {
      32        3399 :   if (!dynamic_cast<PINSFVSuperficialVelocityVariable *>(&_var))
      33           0 :     mooseError("PINSFVMomentumDiffusion may only be used with a superficial velocity "
      34             :                "variable, of variable type PINSFVSuperficialVelocityVariable.");
      35        3399 : }
      36             : 
      37             : ADReal
      38    30597676 : PINSFVMomentumDiffusion::computeStrongResidual(const bool populate_a_coeffs)
      39             : {
      40             :   using namespace Moose::FV;
      41             : 
      42    30597676 :   const bool has_elem = (_face_type == FaceInfo::VarFaceNeighbors::ELEM ||
      43    30597676 :                          _face_type == FaceInfo::VarFaceNeighbors::BOTH);
      44    30597676 :   const bool has_neighbor = (_face_type == FaceInfo::VarFaceNeighbors::NEIGHBOR ||
      45    30597676 :                              _face_type == FaceInfo::VarFaceNeighbors::BOTH);
      46             : 
      47    30597676 :   const auto elem_face = elemArg();
      48    30597676 :   const auto neighbor_face = neighborArg();
      49    30597676 :   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    30597676 :   const auto mu_elem = has_elem ? _mu(elem_face, state) : _mu(neighbor_face, state);
      56    30597676 :   const auto eps_elem = has_elem ? _eps(elem_face, state) : _eps(neighbor_face, state);
      57             : 
      58             :   ADReal mu_neighbor;
      59             :   ADReal eps_neighbor;
      60    30597676 :   if (onBoundary(*_face_info))
      61     1033043 :     mu_face = _mu(singleSidedFaceArg(), state);
      62             :   else
      63             :   {
      64    29564633 :     mu_neighbor = has_elem ? _mu(neighbor_face, state) : _mu(elem_face, state);
      65    29564633 :     eps_neighbor = has_neighbor ? _eps(neighbor_face, state) : _eps(elem_face, state);
      66    29564633 :     interpolate(Moose::FV::InterpMethod::Average, mu_face, mu_elem, mu_neighbor, *_face_info, true);
      67             :   }
      68             : 
      69             :   // Compute face superficial velocity gradient
      70    61195352 :   auto dudn = _var.gradient(makeCDFace(*_face_info), state) * _face_info->normal();
      71             : 
      72    30597676 :   if (populate_a_coeffs)
      73             :   {
      74    28684176 :     if (has_elem)
      75             :     {
      76    28684076 :       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    28684076 :       _ae = dudn.derivatives()[dof_number];
      80    28684076 :       _ae *= -mu_face;
      81             :     }
      82    28684176 :     if (has_neighbor)
      83             :     {
      84    27737463 :       const auto dof_number = _face_info->neighbor().dof_number(_sys.number(), _var.number(), 0);
      85    27737463 :       _an = dudn.derivatives()[dof_number];
      86    27737463 :       _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    61195352 :           ? MetaPhysicL::raw_value(_eps.gradient(makeCDFace(*_face_info), state))
      97    30597676 :           : MetaPhysicL::raw_value(_eps.gradient(
      98    31632551 :                 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    30597676 :       mu_elem / eps_elem * (has_elem ? _var(elem_face, state) : _var(neighbor_face, state));
     106    30597676 :   if (onBoundary(*_face_info))
     107     1033043 :     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    29564633 :     const auto coeff_neighbor = mu_neighbor / eps_neighbor *
     112    29564633 :                                 (has_elem ? _var(neighbor_face, state) : _var(elem_face, state));
     113    29564633 :     interpolate(Moose::FV::InterpMethod::Average,
     114             :                 coeff_face,
     115             :                 coeff_one_side,
     116             :                 coeff_neighbor,
     117    29564633 :                 *_face_info,
     118             :                 true);
     119             :   }
     120             : 
     121    30597676 :   residual -= coeff_face * grad_eps_face * _normal;
     122             : 
     123    30597676 :   return -residual;
     124             : }