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     4486474 : 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             :   using std::sqrt, std::min, std::max;
      55             : 
      56     4486474 :   const ADReal rhoAL = UL[THMVACE3D::RHOA];
      57     4486474 :   const ADReal rhouAL = UL[THMVACE3D::RHOUA];
      58     4486474 :   const ADReal rhovAL = UL[THMVACE3D::RHOVA];
      59     4486474 :   const ADReal rhowAL = UL[THMVACE3D::RHOWA];
      60     4486474 :   const ADReal rhoEAL = UL[THMVACE3D::RHOEA];
      61     4486474 :   const ADReal AL = UL[THMVACE3D::AREA];
      62             : 
      63     4486474 :   const ADReal rhoAR = UR[THMVACE3D::RHOA];
      64     4486474 :   const ADReal rhouAR = UR[THMVACE3D::RHOUA];
      65     4486474 :   const ADReal rhovAR = UR[THMVACE3D::RHOVA];
      66     4486474 :   const ADReal rhowAR = UR[THMVACE3D::RHOWA];
      67     4486474 :   const ADReal rhoEAR = UR[THMVACE3D::RHOEA];
      68     4486474 :   const ADReal AR = UR[THMVACE3D::AREA];
      69             : 
      70             :   // compute the primitive variables
      71             : 
      72             :   const ADReal rhoL = rhoAL / AL;
      73     4486474 :   const ADRealVectorValue uvecL(rhouAL / rhoAL, rhovAL / rhoAL, rhowAL / rhoAL);
      74     4486474 :   const ADReal unL = uvecL * nLR;
      75     4486474 :   const ADReal ut1L = uvecL * t1;
      76     4486474 :   const ADReal ut2L = uvecL * t2;
      77             :   const ADReal rhoEL = rhoEAL / AL;
      78     8972948 :   const ADReal vL = 1.0 / rhoL;
      79             :   const ADReal EL = rhoEAL / rhoAL;
      80     4486474 :   const ADReal eL = EL - 0.5 * uvecL * uvecL;
      81     4486474 :   const ADReal pL = _fp.p_from_v_e(vL, eL);
      82     4486474 :   const ADReal cL = _fp.c_from_v_e(vL, eL);
      83             : 
      84             :   const ADReal rhoR = rhoAR / AR;
      85     4486474 :   const ADRealVectorValue uvecR(rhouAR / rhoAR, rhovAR / rhoAR, rhowAR / rhoAR);
      86     4486474 :   const ADReal unR = uvecR * nLR;
      87     4486474 :   const ADReal ut1R = uvecR * t1;
      88     4486474 :   const ADReal ut2R = uvecR * t2;
      89             :   const ADReal rhoER = rhoEAR / AR;
      90     8972948 :   const ADReal vR = 1.0 / rhoR;
      91             :   const ADReal ER = rhoEAR / rhoAR;
      92     4486474 :   const ADReal eR = ER - 0.5 * uvecR * uvecR;
      93     4486474 :   const ADReal pR = _fp.p_from_v_e(vR, eR);
      94     4486474 :   const ADReal cR = _fp.c_from_v_e(vR, eR);
      95             : 
      96             :   // compute left and right wave speeds
      97             :   ADReal sL, sR;
      98     4486474 :   if (_wave_speed_formulation == WaveSpeedFormulation::EINFELDT)
      99             :   {
     100             :     // compute Roe-averaged variables
     101     4482514 :     const ADReal sqrt_rhoL = sqrt(rhoL);
     102     4482514 :     const ADReal sqrt_rhoR = sqrt(rhoR);
     103     4482514 :     const ADReal un_roe = (sqrt_rhoL * unL + sqrt_rhoR * unR) / (sqrt_rhoL + sqrt_rhoR);
     104     4482514 :     const ADReal ut1_roe = (sqrt_rhoL * ut1L + sqrt_rhoR * ut1R) / (sqrt_rhoL + sqrt_rhoR);
     105     4482514 :     const ADReal ut2_roe = (sqrt_rhoL * ut2L + sqrt_rhoR * ut2R) / (sqrt_rhoL + sqrt_rhoR);
     106     4482514 :     const ADReal HL = EL + pL / rhoL;
     107     4482514 :     const ADReal HR = ER + pR / rhoR;
     108     4482514 :     const ADReal H_roe = (sqrt_rhoL * HL + sqrt_rhoR * HR) / (sqrt_rhoL + sqrt_rhoR);
     109             :     const ADRealVectorValue uvec_roe(un_roe, ut1_roe, ut2_roe);
     110     8965028 :     const ADReal h_roe = H_roe - 0.5 * uvec_roe * uvec_roe;
     111     4482514 :     const ADReal rho_roe = sqrt(rhoL * rhoR);
     112     4482514 :     const ADReal v_roe = 1.0 / rho_roe;
     113     4482514 :     const ADReal e_roe = _fp.e_from_v_h(v_roe, h_roe);
     114     4482514 :     const ADReal c_roe = _fp.c_from_v_e(v_roe, e_roe);
     115             : 
     116     4482514 :     sL = min(unL - cL, un_roe - c_roe);
     117     4482514 :     sR = max(unR + cR, un_roe + c_roe);
     118             :   }
     119        3960 :   else if (_wave_speed_formulation == WaveSpeedFormulation::DAVIS)
     120             :   {
     121        3960 :     sL = min(unL - cL, unR - cR);
     122        3960 :     sR = max(unL + cL, unR + cR);
     123             :   }
     124             :   else
     125             :   {
     126             :     mooseAssert(false, "Invalid 'wave_speed_formulation'.");
     127             :   }
     128             : 
     129             :   // compute middle wave speed
     130     4486474 :   const ADReal sm = (rhoR * unR * (sR - unR) - rhoL * unL * (sL - unL) + pL - pR) /
     131     4486474 :                     (rhoR * (sR - unR) - rhoL * (sL - unL));
     132             : 
     133             :   // compute Omega_L, Omega_R
     134     8972948 :   const ADReal omegL = 1.0 / (sL - sm);
     135     8972948 :   const ADReal omegR = 1.0 / (sR - sm);
     136             : 
     137             :   // compute p^*
     138     4486474 :   const ADReal ps = rhoL * (sL - unL) * (sm - unL) + pL;
     139             : 
     140             :   // compute U_L^*, U_R^*
     141             : 
     142     4486474 :   const ADReal rhoLs = omegL * (sL - unL) * rhoL;
     143     4486474 :   const ADReal rhounLs = omegL * ((sL - unL) * rhoL * unL + ps - pL);
     144     4486474 :   const ADReal rhoELs = omegL * ((sL - unL) * rhoEL - pL * unL + ps * sm);
     145             : 
     146     4486474 :   const ADReal rhoRs = omegR * (sR - unR) * rhoR;
     147     4486474 :   const ADReal rhounRs = omegR * ((sR - unR) * rhoR * unR + ps - pR);
     148     4486474 :   const ADReal rhoERs = omegR * ((sR - unR) * rhoER - pR * unR + ps * sm);
     149             : 
     150     4486474 :   std::vector<ADReal> UL_1d(THMVACE1D::N_FLUX_INPUTS);
     151     4486474 :   UL_1d[THMVACE1D::RHOA] = UL[THMVACE3D::RHOA];
     152     4486474 :   UL_1d[THMVACE1D::RHOUA] = UL[THMVACE3D::RHOUA];
     153     4486474 :   UL_1d[THMVACE1D::RHOEA] = UL[THMVACE3D::RHOEA];
     154     4486474 :   UL_1d[THMVACE1D::AREA] = UL[THMVACE3D::AREA];
     155             : 
     156     4486474 :   std::vector<ADReal> UR_1d(THMVACE1D::N_FLUX_INPUTS);
     157     4486474 :   UR_1d[THMVACE1D::RHOA] = UR[THMVACE3D::RHOA];
     158     4486474 :   UR_1d[THMVACE1D::RHOUA] = UR[THMVACE3D::RHOUA];
     159     4486474 :   UR_1d[THMVACE1D::RHOEA] = UR[THMVACE3D::RHOEA];
     160     4486474 :   UR_1d[THMVACE1D::AREA] = UR[THMVACE3D::AREA];
     161             : 
     162     4486474 :   const ADReal A_flow = computeFlowArea(UL_1d, UR_1d);
     163             : 
     164             :   // compute the fluxes
     165     4486474 :   FL.resize(THMVACE3D::N_FLUX_OUTPUTS);
     166     4486474 :   if (sL > 0.0)
     167             :   {
     168        3921 :     FL[THMVACE3D::MASS] = unL * rhoL * A_flow;
     169        3921 :     FL[THMVACE3D::MOM_NORM] = (unL * rhoL * unL + pL) * A_flow;
     170        3921 :     FL[THMVACE3D::MOM_TAN1] = rhoL * unL * ut1L * A_flow;
     171        3921 :     FL[THMVACE3D::MOM_TAN2] = rhoL * unL * ut2L * A_flow;
     172        3921 :     FL[THMVACE3D::ENERGY] = unL * (rhoEL + pL) * A_flow;
     173             : 
     174        3921 :     _last_region_index = 0;
     175             :   }
     176     4482553 :   else if (sL <= 0.0 && sm > 0.0)
     177             :   {
     178     3610911 :     FL[THMVACE3D::MASS] = sm * rhoLs * A_flow;
     179     3610911 :     FL[THMVACE3D::MOM_NORM] = (sm * rhounLs + ps) * A_flow;
     180     3610911 :     FL[THMVACE3D::MOM_TAN1] = rhounLs * ut1L * A_flow;
     181     3610911 :     FL[THMVACE3D::MOM_TAN2] = rhounLs * ut2L * A_flow;
     182     3610911 :     FL[THMVACE3D::ENERGY] = sm * (rhoELs + ps) * A_flow;
     183             : 
     184     3610911 :     _last_region_index = 1;
     185             :   }
     186      871642 :   else if (sm <= 0.0 && sR >= 0.0)
     187             :   {
     188      867663 :     FL[THMVACE3D::MASS] = sm * rhoRs * A_flow;
     189      867663 :     FL[THMVACE3D::MOM_NORM] = (sm * rhounRs + ps) * A_flow;
     190      867663 :     FL[THMVACE3D::MOM_TAN1] = rhounRs * ut1R * A_flow;
     191      867663 :     FL[THMVACE3D::MOM_TAN2] = rhounRs * ut2R * A_flow;
     192      867663 :     FL[THMVACE3D::ENERGY] = sm * (rhoERs + ps) * A_flow;
     193             : 
     194      867663 :     _last_region_index = 2;
     195             :   }
     196        3979 :   else if (sR < 0.0)
     197             :   {
     198        3929 :     FL[THMVACE3D::MASS] = unR * rhoR * A_flow;
     199        3929 :     FL[THMVACE3D::MOM_NORM] = (unR * rhoR * unR + pR) * A_flow;
     200        3929 :     FL[THMVACE3D::MOM_TAN1] = rhoR * unR * ut1R * A_flow;
     201        3929 :     FL[THMVACE3D::MOM_TAN2] = rhoR * unR * ut2R * A_flow;
     202        3929 :     FL[THMVACE3D::ENERGY] = unR * (rhoER + pR) * A_flow;
     203             : 
     204        3929 :     _last_region_index = 3;
     205             :   }
     206             :   else
     207          50 :     std::fill(FL.begin(), FL.end(), getNaN());
     208             : 
     209     4486474 :   FR = FL;
     210             : 
     211             :   const ADReal A_wall_L = AL - A_flow;
     212     4486474 :   FL[THMVACE3D::MOM_NORM] += pL * A_wall_L;
     213             : 
     214             :   const ADReal A_wall_R = AR - A_flow;
     215     4486474 :   FR[THMVACE3D::MOM_NORM] += pR * A_wall_R;
     216     4486474 : }
     217             : 
     218             : ADReal
     219     4486474 : ADNumericalFlux3EqnHLLC::computeFlowArea(const std::vector<ADReal> & UL,
     220             :                                          const std::vector<ADReal> & UR) const
     221             : {
     222     4486474 :   return std::min(UL[THMVACE1D::AREA], UR[THMVACE1D::AREA]);
     223             : }