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 "ADNumericalFlux3EqnHLLC.h"
      11             : #include "THMIndicesVACE.h"
      12             : #include "Numerics.h"
      13             : 
      14             : registerMooseObject("ThermalHydraulicsApp", ADNumericalFlux3EqnHLLC);
      15             : 
      16             : InputParameters
      17        8881 : ADNumericalFlux3EqnHLLC::validParams()
      18             : {
      19        8881 :   InputParameters params = ADNumericalFlux3EqnBase::validParams();
      20        8881 :   params += NaNInterface::validParams();
      21             : 
      22       17762 :   MooseEnum wave_speed_formulation("einfeldt davis", "einfeldt");
      23       17762 :   params.addParam<MooseEnum>(
      24             :       "wave_speed_formulation", wave_speed_formulation, "Method for computing wave speeds");
      25             : 
      26       17762 :   params.addRequiredParam<UserObjectName>("fluid_properties",
      27             :                                           "Name for fluid properties user object");
      28             : 
      29        8881 :   params.addClassDescription("Computes internal side flux for the 1-D, 1-phase, variable-area "
      30             :                              "Euler equations using the HLLC approximate Riemann solver.");
      31             : 
      32        8881 :   return params;
      33        8881 : }
      34             : 
      35        4850 : ADNumericalFlux3EqnHLLC::ADNumericalFlux3EqnHLLC(const InputParameters & parameters)
      36             :   : ADNumericalFlux3EqnBase(parameters),
      37             :     NaNInterface(this),
      38        4850 :     _fp(getUserObject<SinglePhaseFluidProperties>("fluid_properties")),
      39        4850 :     _wave_speed_formulation(
      40        9700 :         getParam<MooseEnum>("wave_speed_formulation").getEnum<WaveSpeedFormulation>())
      41             : {
      42        4850 : }
      43             : 
      44             : void
      45     4495185 : ADNumericalFlux3EqnHLLC::calcFlux(const std::vector<ADReal> & UL,
      46             :                                   const std::vector<ADReal> & UR,
      47             :                                   const RealVectorValue & nLR,
      48             :                                   const RealVectorValue & t1,
      49             :                                   const RealVectorValue & t2,
      50             :                                   std::vector<ADReal> & FL,
      51             :                                   std::vector<ADReal> & FR) const
      52             : {
      53             :   // extract the conserved variables and area
      54             : 
      55     4495185 :   const ADReal rhoAL = UL[THMVACE3D::RHOA];
      56     4495185 :   const ADReal rhouAL = UL[THMVACE3D::RHOUA];
      57     4495185 :   const ADReal rhovAL = UL[THMVACE3D::RHOVA];
      58     4495185 :   const ADReal rhowAL = UL[THMVACE3D::RHOWA];
      59     4495185 :   const ADReal rhoEAL = UL[THMVACE3D::RHOEA];
      60     4495185 :   const ADReal AL = UL[THMVACE3D::AREA];
      61             : 
      62     4495185 :   const ADReal rhoAR = UR[THMVACE3D::RHOA];
      63     4495185 :   const ADReal rhouAR = UR[THMVACE3D::RHOUA];
      64     4495185 :   const ADReal rhovAR = UR[THMVACE3D::RHOVA];
      65     4495185 :   const ADReal rhowAR = UR[THMVACE3D::RHOWA];
      66     4495185 :   const ADReal rhoEAR = UR[THMVACE3D::RHOEA];
      67     4495185 :   const ADReal AR = UR[THMVACE3D::AREA];
      68             : 
      69             :   // compute the primitive variables
      70             : 
      71             :   const ADReal rhoL = rhoAL / AL;
      72     4495185 :   const ADRealVectorValue uvecL(rhouAL / rhoAL, rhovAL / rhoAL, rhowAL / rhoAL);
      73     4495185 :   const ADReal unL = uvecL * nLR;
      74     4495185 :   const ADReal ut1L = uvecL * t1;
      75     4495185 :   const ADReal ut2L = uvecL * t2;
      76             :   const ADReal rhoEL = rhoEAL / AL;
      77     8990370 :   const ADReal vL = 1.0 / rhoL;
      78             :   const ADReal EL = rhoEAL / rhoAL;
      79     4495185 :   const ADReal eL = EL - 0.5 * uvecL * uvecL;
      80     4495185 :   const ADReal pL = _fp.p_from_v_e(vL, eL);
      81     4495185 :   const ADReal cL = _fp.c_from_v_e(vL, eL);
      82             : 
      83             :   const ADReal rhoR = rhoAR / AR;
      84     4495185 :   const ADRealVectorValue uvecR(rhouAR / rhoAR, rhovAR / rhoAR, rhowAR / rhoAR);
      85     4495185 :   const ADReal unR = uvecR * nLR;
      86     4495185 :   const ADReal ut1R = uvecR * t1;
      87     4495185 :   const ADReal ut2R = uvecR * t2;
      88             :   const ADReal rhoER = rhoEAR / AR;
      89     8990370 :   const ADReal vR = 1.0 / rhoR;
      90             :   const ADReal ER = rhoEAR / rhoAR;
      91     4495185 :   const ADReal eR = ER - 0.5 * uvecR * uvecR;
      92     4495185 :   const ADReal pR = _fp.p_from_v_e(vR, eR);
      93     4495185 :   const ADReal cR = _fp.c_from_v_e(vR, eR);
      94             : 
      95             :   // compute left and right wave speeds
      96             :   ADReal sL, sR;
      97     4495185 :   if (_wave_speed_formulation == WaveSpeedFormulation::EINFELDT)
      98             :   {
      99             :     // compute Roe-averaged variables
     100     4491225 :     const ADReal sqrt_rhoL = std::sqrt(rhoL);
     101     4491225 :     const ADReal sqrt_rhoR = std::sqrt(rhoR);
     102     4491225 :     const ADReal un_roe = (sqrt_rhoL * unL + sqrt_rhoR * unR) / (sqrt_rhoL + sqrt_rhoR);
     103     4491225 :     const ADReal ut1_roe = (sqrt_rhoL * ut1L + sqrt_rhoR * ut1R) / (sqrt_rhoL + sqrt_rhoR);
     104     4491225 :     const ADReal ut2_roe = (sqrt_rhoL * ut2L + sqrt_rhoR * ut2R) / (sqrt_rhoL + sqrt_rhoR);
     105     4491225 :     const ADReal HL = EL + pL / rhoL;
     106     4491225 :     const ADReal HR = ER + pR / rhoR;
     107     4491225 :     const ADReal H_roe = (sqrt_rhoL * HL + sqrt_rhoR * HR) / (sqrt_rhoL + sqrt_rhoR);
     108             :     const ADRealVectorValue uvec_roe(un_roe, ut1_roe, ut2_roe);
     109     8982450 :     const ADReal h_roe = H_roe - 0.5 * uvec_roe * uvec_roe;
     110     4491225 :     const ADReal rho_roe = std::sqrt(rhoL * rhoR);
     111     4491225 :     const ADReal v_roe = 1.0 / rho_roe;
     112     4491225 :     const ADReal e_roe = _fp.e_from_v_h(v_roe, h_roe);
     113     4491225 :     const ADReal c_roe = _fp.c_from_v_e(v_roe, e_roe);
     114             : 
     115     4491225 :     sL = std::min(unL - cL, un_roe - c_roe);
     116     4491225 :     sR = std::max(unR + cR, un_roe + c_roe);
     117             :   }
     118        3960 :   else if (_wave_speed_formulation == WaveSpeedFormulation::DAVIS)
     119             :   {
     120        3960 :     sL = std::min(unL - cL, unR - cR);
     121        3960 :     sR = std::max(unL + cL, unR + cR);
     122             :   }
     123             :   else
     124             :   {
     125             :     mooseAssert(false, "Invalid 'wave_speed_formulation'.");
     126             :   }
     127             : 
     128             :   // compute middle wave speed
     129     4495185 :   const ADReal sm = (rhoR * unR * (sR - unR) - rhoL * unL * (sL - unL) + pL - pR) /
     130     4495185 :                     (rhoR * (sR - unR) - rhoL * (sL - unL));
     131             : 
     132             :   // compute Omega_L, Omega_R
     133     8990370 :   const ADReal omegL = 1.0 / (sL - sm);
     134     8990370 :   const ADReal omegR = 1.0 / (sR - sm);
     135             : 
     136             :   // compute p^*
     137     4495185 :   const ADReal ps = rhoL * (sL - unL) * (sm - unL) + pL;
     138             : 
     139             :   // compute U_L^*, U_R^*
     140             : 
     141     4495185 :   const ADReal rhoLs = omegL * (sL - unL) * rhoL;
     142     4495185 :   const ADReal rhounLs = omegL * ((sL - unL) * rhoL * unL + ps - pL);
     143     4495185 :   const ADReal rhoELs = omegL * ((sL - unL) * rhoEL - pL * unL + ps * sm);
     144             : 
     145     4495185 :   const ADReal rhoRs = omegR * (sR - unR) * rhoR;
     146     4495185 :   const ADReal rhounRs = omegR * ((sR - unR) * rhoR * unR + ps - pR);
     147     4495185 :   const ADReal rhoERs = omegR * ((sR - unR) * rhoER - pR * unR + ps * sm);
     148             : 
     149     4495185 :   std::vector<ADReal> UL_1d(THMVACE1D::N_FLUX_INPUTS);
     150     4495185 :   UL_1d[THMVACE1D::RHOA] = UL[THMVACE3D::RHOA];
     151     4495185 :   UL_1d[THMVACE1D::RHOUA] = UL[THMVACE3D::RHOUA];
     152     4495185 :   UL_1d[THMVACE1D::RHOEA] = UL[THMVACE3D::RHOEA];
     153     4495185 :   UL_1d[THMVACE1D::AREA] = UL[THMVACE3D::AREA];
     154             : 
     155     4495185 :   std::vector<ADReal> UR_1d(THMVACE1D::N_FLUX_INPUTS);
     156     4495185 :   UR_1d[THMVACE1D::RHOA] = UR[THMVACE3D::RHOA];
     157     4495185 :   UR_1d[THMVACE1D::RHOUA] = UR[THMVACE3D::RHOUA];
     158     4495185 :   UR_1d[THMVACE1D::RHOEA] = UR[THMVACE3D::RHOEA];
     159     4495185 :   UR_1d[THMVACE1D::AREA] = UR[THMVACE3D::AREA];
     160             : 
     161     4495185 :   const ADReal A_flow = computeFlowArea(UL_1d, UR_1d);
     162             : 
     163             :   // compute the fluxes
     164     4495185 :   FL.resize(THMVACE3D::N_FLUX_OUTPUTS);
     165     4495185 :   if (sL > 0.0)
     166             :   {
     167        3921 :     FL[THMVACE3D::MASS] = unL * rhoL * A_flow;
     168        3921 :     FL[THMVACE3D::MOM_NORM] = (unL * rhoL * unL + pL) * A_flow;
     169        3921 :     FL[THMVACE3D::MOM_TAN1] = rhoL * unL * ut1L * A_flow;
     170        3921 :     FL[THMVACE3D::MOM_TAN2] = rhoL * unL * ut2L * A_flow;
     171        3921 :     FL[THMVACE3D::ENERGY] = unL * (rhoEL + pL) * A_flow;
     172             : 
     173        3921 :     _last_region_index = 0;
     174             :   }
     175     4491264 :   else if (sL <= 0.0 && sm > 0.0)
     176             :   {
     177     3617927 :     FL[THMVACE3D::MASS] = sm * rhoLs * A_flow;
     178     3617927 :     FL[THMVACE3D::MOM_NORM] = (sm * rhounLs + ps) * A_flow;
     179     3617927 :     FL[THMVACE3D::MOM_TAN1] = rhounLs * ut1L * A_flow;
     180     3617927 :     FL[THMVACE3D::MOM_TAN2] = rhounLs * ut2L * A_flow;
     181     3617927 :     FL[THMVACE3D::ENERGY] = sm * (rhoELs + ps) * A_flow;
     182             : 
     183     3617927 :     _last_region_index = 1;
     184             :   }
     185      873337 :   else if (sm <= 0.0 && sR >= 0.0)
     186             :   {
     187      869358 :     FL[THMVACE3D::MASS] = sm * rhoRs * A_flow;
     188      869358 :     FL[THMVACE3D::MOM_NORM] = (sm * rhounRs + ps) * A_flow;
     189      869358 :     FL[THMVACE3D::MOM_TAN1] = rhounRs * ut1R * A_flow;
     190      869358 :     FL[THMVACE3D::MOM_TAN2] = rhounRs * ut2R * A_flow;
     191      869358 :     FL[THMVACE3D::ENERGY] = sm * (rhoERs + ps) * A_flow;
     192             : 
     193      869358 :     _last_region_index = 2;
     194             :   }
     195        3979 :   else if (sR < 0.0)
     196             :   {
     197        3929 :     FL[THMVACE3D::MASS] = unR * rhoR * A_flow;
     198        3929 :     FL[THMVACE3D::MOM_NORM] = (unR * rhoR * unR + pR) * A_flow;
     199        3929 :     FL[THMVACE3D::MOM_TAN1] = rhoR * unR * ut1R * A_flow;
     200        3929 :     FL[THMVACE3D::MOM_TAN2] = rhoR * unR * ut2R * A_flow;
     201        3929 :     FL[THMVACE3D::ENERGY] = unR * (rhoER + pR) * A_flow;
     202             : 
     203        3929 :     _last_region_index = 3;
     204             :   }
     205             :   else
     206          50 :     std::fill(FL.begin(), FL.end(), getNaN());
     207             : 
     208     4495185 :   FR = FL;
     209             : 
     210             :   const ADReal A_wall_L = AL - A_flow;
     211     4495185 :   FL[THMVACE3D::MOM_NORM] += pL * A_wall_L;
     212             : 
     213             :   const ADReal A_wall_R = AR - A_flow;
     214     4495185 :   FR[THMVACE3D::MOM_NORM] += pR * A_wall_R;
     215     4495185 : }
     216             : 
     217             : ADReal
     218     4495185 : ADNumericalFlux3EqnHLLC::computeFlowArea(const std::vector<ADReal> & UL,
     219             :                                          const std::vector<ADReal> & UR) const
     220             : {
     221     4495185 :   return std::min(UL[THMVACE1D::AREA], UR[THMVACE1D::AREA]);
     222             : }