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 "PCNSFVHLLC.h"
      11             : #include "NS.h"
      12             : #include "SinglePhaseFluidProperties.h"
      13             : 
      14             : using MetaPhysicL::raw_value;
      15             : 
      16             : InputParameters
      17         663 : PCNSFVHLLC::validParams()
      18             : {
      19         663 :   InputParameters params = FVFluxKernel::validParams();
      20         663 :   params.addRequiredParam<UserObjectName>(NS::fluid, "Fluid properties userobject");
      21         663 :   params.set<unsigned short>("ghost_layers") = 2;
      22         663 :   return params;
      23           0 : }
      24             : 
      25         336 : PCNSFVHLLC::PCNSFVHLLC(const InputParameters & params)
      26             :   : FVFluxKernel(params),
      27         336 :     _fluid(UserObjectInterface::getUserObject<SinglePhaseFluidProperties>(NS::fluid)),
      28         336 :     _specific_internal_energy_elem(getADMaterialProperty<Real>(NS::specific_internal_energy)),
      29         336 :     _specific_internal_energy_neighbor(
      30             :         getNeighborADMaterialProperty<Real>(NS::specific_internal_energy)),
      31         336 :     _rho_et_elem(getADMaterialProperty<Real>(NS::total_energy_density)),
      32         336 :     _rho_et_neighbor(getNeighborADMaterialProperty<Real>(NS::total_energy_density)),
      33         336 :     _vel_elem(getADMaterialProperty<RealVectorValue>(NS::velocity)),
      34         336 :     _vel_neighbor(getNeighborADMaterialProperty<RealVectorValue>(NS::velocity)),
      35         336 :     _speed_elem(getADMaterialProperty<Real>(NS::speed)),
      36         336 :     _speed_neighbor(getNeighborADMaterialProperty<Real>(NS::speed)),
      37         336 :     _rho_elem(getADMaterialProperty<Real>(NS::density)),
      38         336 :     _rho_neighbor(getNeighborADMaterialProperty<Real>(NS::density)),
      39         336 :     _pressure_elem(getADMaterialProperty<Real>(NS::pressure)),
      40         336 :     _pressure_neighbor(getNeighborADMaterialProperty<Real>(NS::pressure)),
      41         672 :     _eps_elem(getMaterialProperty<Real>(NS::porosity)),
      42         336 :     _eps_neighbor(getNeighborMaterialProperty<Real>(NS::porosity))
      43             : {
      44         336 : }
      45             : 
      46             : std::array<ADReal, 3>
      47      244125 : PCNSFVHLLC::waveSpeed(const ADReal & rho_elem,
      48             :                       const ADRealVectorValue & vel_elem,
      49             :                       const ADReal & e_elem,
      50             :                       const Real eps_elem,
      51             :                       const ADReal & rho_neighbor,
      52             :                       const ADRealVectorValue & vel_neighbor,
      53             :                       const ADReal & e_neighbor,
      54             :                       const Real eps_neighbor,
      55             :                       const SinglePhaseFluidProperties & fluid,
      56             :                       const ADRealVectorValue & normal)
      57             : {
      58             :   const auto & rho1 = rho_elem;
      59      244125 :   const auto u1 = vel_elem.norm();
      60      244125 :   const auto q1 = normal * vel_elem;
      61      244125 :   const auto v1 = 1.0 / rho1;
      62             :   const auto & e1 = e_elem;
      63      244125 :   const auto et1 = e1 + 0.5 * u1 * u1;
      64      244125 :   const auto p1 = fluid.p_from_v_e(v1, e1);
      65      244125 :   const auto ht1 = et1 + p1 / rho1;
      66      244125 :   const auto c1 = fluid.c_from_v_e(v1, e1);
      67      244125 :   const auto eps1 = eps_elem;
      68             : 
      69             :   const auto & rho2 = rho_neighbor;
      70      244125 :   const auto u2 = vel_neighbor.norm();
      71      244125 :   const auto q2 = normal * vel_neighbor;
      72      244125 :   const auto v2 = 1.0 / rho2;
      73             :   const auto & e2 = e_neighbor;
      74      244125 :   const auto et2 = e2 + 0.5 * u2 * u2;
      75      244125 :   const auto p2 = fluid.p_from_v_e(v2, e2);
      76      244125 :   const auto ht2 = et2 + p2 / rho2;
      77      244125 :   const auto c2 = fluid.c_from_v_e(v2, e2);
      78      244125 :   const auto eps2 = eps_neighbor;
      79             : 
      80             :   // compute Roe-averaged variables
      81      244125 :   const auto sqrt_rho1 = std::sqrt(rho1);
      82      244125 :   const auto sqrt_rho2 = std::sqrt(rho2);
      83      244125 :   const auto u_roe = (sqrt_rho1 * u1 + sqrt_rho2 * u2) / (sqrt_rho1 + sqrt_rho2);
      84      244125 :   const auto q_roe = (sqrt_rho1 * q1 + sqrt_rho2 * q2) / (sqrt_rho1 + sqrt_rho2);
      85      244125 :   const auto e_roe = (sqrt_rho1 * e1 + sqrt_rho2 * e2) / (sqrt_rho1 + sqrt_rho2);
      86      244125 :   const auto rho_roe = std::sqrt(rho1 * rho2);
      87      244125 :   const auto v_roe = 1.0 / rho_roe;
      88      244125 :   const auto c_roe = fluid.c_from_v_e(v_roe, e_roe);
      89             : 
      90             :   // compute wave speeds
      91             :   // I may want to change the estimate of these wave speeds!
      92      244125 :   auto SL = std::min(q1 - c1, q_roe - c_roe);
      93      244125 :   auto SR = std::max(q2 + c2, q_roe + c_roe);
      94      244125 :   auto SM = (eps2 * rho2 * q2 * (SR - q2) - eps1 * rho1 * q1 * (SL - q1) + eps1 * p1 - eps2 * p2) /
      95      244125 :             (eps2 * rho2 * (SR - q2) - eps1 * rho1 * (SL - q1));
      96             : 
      97      244125 :   return {{std::move(SL), std::move(SM), std::move(SR)}};
      98      244125 : }
      99             : 
     100             : ADReal
     101      241635 : PCNSFVHLLC::computeQpResidual()
     102             : {
     103      241635 :   _normal_speed_elem = _normal * _vel_elem[_qp];
     104      241635 :   _normal_speed_neighbor = _normal * _vel_neighbor[_qp];
     105      241635 :   auto wave_speeds = waveSpeed(_rho_elem[_qp],
     106      241635 :                                _vel_elem[_qp],
     107      241635 :                                _specific_internal_energy_elem[_qp],
     108      241635 :                                _eps_elem[_qp],
     109      241635 :                                _rho_neighbor[_qp],
     110      241635 :                                _vel_neighbor[_qp],
     111      241635 :                                _specific_internal_energy_neighbor[_qp],
     112      241635 :                                _eps_neighbor[_qp],
     113             :                                _fluid,
     114      241635 :                                _normal);
     115      241635 :   _SL = std::move(wave_speeds[0]);
     116      241635 :   _SM = std::move(wave_speeds[1]);
     117      241635 :   _SR = std::move(wave_speeds[2]);
     118      241635 :   if (_SL >= 0)
     119           0 :     return fluxElem();
     120      241635 :   else if (_SR <= 0)
     121           0 :     return fluxNeighbor();
     122             :   else
     123             :   {
     124      241635 :     if (_SM >= 0)
     125             :     {
     126      241635 :       ADReal f = _eps_elem[_qp] * _rho_elem[_qp] * (_SL - _normal_speed_elem) / (_SL - _SM);
     127      724905 :       return fluxElem() + _SL * (f * hllcElem() - conservedVariableElem());
     128             :     }
     129             :     else
     130             :     {
     131             :       ADReal f =
     132           0 :           _eps_neighbor[_qp] * _rho_neighbor[_qp] * (_SR - _normal_speed_neighbor) / (_SR - _SM);
     133           0 :       return fluxNeighbor() + _SR * (f * hllcNeighbor() - conservedVariableNeighbor());
     134             :     }
     135             :   }
     136             :   mooseError("Should never get here");
     137             : }