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 "HLLCUserObject.h"
      11             : #include "NS.h"
      12             : #include "SinglePhaseFluidProperties.h"
      13             : 
      14             : using MetaPhysicL::raw_value;
      15             : 
      16             : registerMooseObject("NavierStokesApp", HLLCUserObject);
      17             : 
      18             : InputParameters
      19          82 : HLLCUserObject::validParams()
      20             : {
      21          82 :   InputParameters params = InternalSideUserObject::validParams();
      22          82 :   params.addRequiredParam<UserObjectName>(NS::fluid, "Fluid properties userobject");
      23          82 :   params.addParam<MaterialPropertyName>(
      24             :       NS::porosity,
      25             :       "An optional parameter for specifying a porosity material property. If not specified "
      26             :       "free-flow conditions will be assumed.");
      27          82 :   params.addClassDescription(
      28             :       "Computes free-flow wave speeds on internal sides, useful in HLLC contexts");
      29          82 :   return params;
      30           0 : }
      31             : 
      32          44 : HLLCUserObject::HLLCUserObject(const InputParameters & parameters)
      33             :   : InternalSideUserObject(parameters),
      34          88 :     _fluid(UserObjectInterface::getUserObject<SinglePhaseFluidProperties>(NS::fluid)),
      35          44 :     _face_info(_mesh.faceInfo()),
      36          44 :     _vel_elem(getADMaterialProperty<RealVectorValue>(NS::velocity)),
      37          44 :     _vel_neighbor(getNeighborADMaterialProperty<RealVectorValue>(NS::velocity)),
      38          44 :     _speed_elem(getADMaterialProperty<Real>(NS::speed)),
      39          44 :     _speed_neighbor(getNeighborADMaterialProperty<Real>(NS::speed)),
      40          44 :     _pressure_elem(getADMaterialProperty<Real>(NS::pressure)),
      41          44 :     _pressure_neighbor(getNeighborADMaterialProperty<Real>(NS::pressure)),
      42          44 :     _rho_elem(getADMaterialProperty<Real>(NS::density)),
      43          44 :     _rho_neighbor(getNeighborADMaterialProperty<Real>(NS::density)),
      44          44 :     _specific_internal_energy_elem(getADMaterialProperty<Real>(NS::specific_internal_energy)),
      45          44 :     _specific_internal_energy_neighbor(
      46             :         getNeighborADMaterialProperty<Real>(NS::specific_internal_energy)),
      47          44 :     _eps_elem(isParamValid(NS::porosity) ? &getMaterialProperty<Real>(NS::porosity) : nullptr),
      48          44 :     _eps_neighbor(isParamValid(NS::porosity) ? &getNeighborMaterialProperty<Real>(NS::porosity)
      49          44 :                                              : nullptr)
      50             : {
      51          44 : }
      52             : 
      53             : void
      54          44 : HLLCUserObject::initialSetup()
      55             : {
      56             :   // This means that we have no finite volume variables, so we need to build the
      57             :   // finite volume information on the mesh so that we can still use this object
      58          44 :   if (!_face_info.size())
      59          38 :     _mesh.setupFiniteVolumeMeshData();
      60             : 
      61         172 :   for (unsigned int j = 0; j < _face_info.size(); ++j)
      62             :   {
      63         128 :     const Elem * elem = &_face_info[j]->elem();
      64             :     unsigned int side = _face_info[j]->elemSideID();
      65             :     side_type elem_side(elem, side);
      66         128 :     _side_to_face_info[elem_side] = j;
      67             :   }
      68          44 : }
      69             : 
      70             : void
      71          18 : HLLCUserObject::execute()
      72             : {
      73             :   using std::sqrt, std::min, std::max;
      74             : 
      75          18 :   const FaceInfo & fi = faceInfoHelper(_current_elem, _current_side);
      76             : 
      77          18 :   const ADReal & rho1 = _rho_elem[_qp];
      78          18 :   const ADReal & u1 = _speed_elem[_qp];
      79          18 :   const ADReal q1 = fi.normal() * _vel_elem[_qp];
      80          18 :   const ADReal v1 = 1.0 / rho1;
      81          18 :   const ADReal & e1 = _specific_internal_energy_elem[_qp];
      82          18 :   const ADReal E1 = e1 + 0.5 * u1 * u1;
      83          18 :   const ADReal & p1 = _pressure_elem[_qp];
      84          18 :   const ADReal H1 = E1 + p1 / rho1;
      85          18 :   const ADReal c1 = _fluid.c_from_v_e(v1, e1);
      86          18 :   const Real eps1 = _eps_elem ? (*_eps_elem)[_qp] : 1.;
      87             : 
      88          18 :   const ADReal & rho2 = _rho_neighbor[_qp];
      89          18 :   const ADReal & u2 = _speed_neighbor[_qp];
      90          18 :   const ADReal q2 = fi.normal() * _vel_neighbor[_qp];
      91          18 :   const ADReal v2 = 1.0 / rho2;
      92          18 :   const ADReal & e2 = _specific_internal_energy_neighbor[_qp];
      93          18 :   const ADReal E2 = e2 + 0.5 * u2 * u2;
      94          18 :   const ADReal & p2 = _pressure_neighbor[_qp];
      95          18 :   const ADReal H2 = E2 + p2 / rho2;
      96          18 :   const ADReal c2 = _fluid.c_from_v_e(v2, e2);
      97          18 :   const Real eps2 = _eps_neighbor ? (*_eps_neighbor)[_qp] : 1.;
      98             : 
      99             :   // compute Roe-averaged variables
     100          18 :   const ADReal sqrt_rho1 = sqrt(rho1);
     101          18 :   const ADReal sqrt_rho2 = sqrt(rho2);
     102          18 :   const ADReal u_roe = (sqrt_rho1 * u1 + sqrt_rho2 * u2) / (sqrt_rho1 + sqrt_rho2);
     103          18 :   const ADReal q_roe = (sqrt_rho1 * q1 + sqrt_rho2 * q2) / (sqrt_rho1 + sqrt_rho2);
     104          18 :   const ADReal H_roe = (sqrt_rho1 * H1 + sqrt_rho2 * H2) / (sqrt_rho1 + sqrt_rho2);
     105          36 :   const ADReal h_roe = H_roe - 0.5 * u_roe * u_roe;
     106          18 :   const ADReal rho_roe = sqrt(rho1 * rho2);
     107          18 :   const ADReal v_roe = 1.0 / rho_roe;
     108          18 :   const ADReal e_roe = _fluid.e_from_v_h(v_roe, h_roe);
     109          18 :   const ADReal c_roe = _fluid.c_from_v_e(v_roe, e_roe);
     110             : 
     111             :   // compute wave speeds
     112          18 :   const ADReal SL = min(q1 - c1, q_roe - c_roe);
     113          18 :   const ADReal SR = max(q2 + c2, q_roe + c_roe);
     114             :   const ADReal SM =
     115          18 :       (eps2 * rho2 * q2 * (SR - q2) - eps1 * rho1 * q1 * (SL - q1) + eps1 * p1 - eps2 * p2) /
     116          18 :       (eps2 * rho2 * (SR - q2) - eps1 * rho1 * (SL - q1));
     117             : 
     118             :   // store these results in _wave_speed
     119          18 :   side_type elem_side(_current_elem, _current_side);
     120          18 :   _wave_speed[elem_side] = {SL, SM, SR};
     121          36 : }
     122             : 
     123             : bool
     124          28 : HLLCUserObject::hasData(const Elem * const elem, const unsigned int side) const
     125             : {
     126             :   side_type elem_side(elem, side);
     127          28 :   return (_wave_speed.find(elem_side) != _wave_speed.end());
     128             : }
     129             : 
     130             : void
     131           6 : HLLCUserObject::threadJoin(const UserObject & y)
     132             : {
     133             :   const auto & pps = static_cast<const HLLCUserObject &>(y);
     134           6 :   for (auto & ws : pps._wave_speed)
     135             :   {
     136           0 :     const auto & it = _wave_speed.find(ws.first);
     137           0 :     if (it == _wave_speed.end())
     138           0 :       _wave_speed[ws.first] = ws.second;
     139             :     else
     140           0 :       for (unsigned int j = 0; j < 3; ++j)
     141           0 :         _wave_speed[ws.first][j] = ws.second[j];
     142             :   }
     143           6 : }
     144             : 
     145             : std::vector<ADReal>
     146          18 : HLLCUserObject::waveSpeed(const Elem * elem, unsigned int side) const
     147             : {
     148             :   side_type elem_side(elem, side);
     149             :   auto it = _wave_speed.find(elem_side);
     150          18 :   if (it == _wave_speed.end())
     151           0 :     mooseError("what the heck (version 2)?");
     152          18 :   return it->second;
     153             : }
     154             : 
     155             : const FaceInfo &
     156          18 : HLLCUserObject::faceInfoHelper(const Elem * elem, unsigned int side) const
     157             : {
     158             :   side_type elem_side(elem, side);
     159             :   auto it = _side_to_face_info.find(elem_side);
     160          18 :   if (it == _side_to_face_info.end())
     161           0 :     mooseError("what the heck?");
     162          18 :   return *_face_info[it->second];
     163             : }