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        4364 : ADNumericalFlux3EqnHLLC::validParams()
      18             : {
      19        4364 :   InputParameters params = ADNumericalFlux3EqnBase::validParams();
      20        4364 :   params += NaNInterface::validParams();
      21             : 
      22        8728 :   MooseEnum wave_speed_formulation("einfeldt davis", "einfeldt");
      23        8728 :   params.addParam<MooseEnum>(
      24             :       "wave_speed_formulation", wave_speed_formulation, "Method for computing wave speeds");
      25             : 
      26        8728 :   params.addRequiredParam<UserObjectName>("fluid_properties",
      27             :                                           "Name for fluid properties user object");
      28             : 
      29        4364 :   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        4364 :   return params;
      33        4364 : }
      34             : 
      35        2325 : ADNumericalFlux3EqnHLLC::ADNumericalFlux3EqnHLLC(const InputParameters & parameters)
      36             :   : ADNumericalFlux3EqnBase(parameters),
      37             :     NaNInterface(this),
      38        2325 :     _fp(getUserObject<SinglePhaseFluidProperties>("fluid_properties")),
      39        2325 :     _wave_speed_formulation(
      40        4650 :         getParam<MooseEnum>("wave_speed_formulation").getEnum<WaveSpeedFormulation>())
      41             : {
      42        2325 : }
      43             : 
      44             : void
      45     2556747 : 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     2556747 :   const ADReal rhoAL = UL[THMVACE3D::RHOA];
      57     2556747 :   const ADReal rhouAL = UL[THMVACE3D::RHOUA];
      58     2556747 :   const ADReal rhovAL = UL[THMVACE3D::RHOVA];
      59     2556747 :   const ADReal rhowAL = UL[THMVACE3D::RHOWA];
      60     2556747 :   const ADReal rhoEAL = UL[THMVACE3D::RHOEA];
      61     2556747 :   const ADReal AL = UL[THMVACE3D::AREA];
      62             : 
      63     2556747 :   const ADReal rhoAR = UR[THMVACE3D::RHOA];
      64     2556747 :   const ADReal rhouAR = UR[THMVACE3D::RHOUA];
      65     2556747 :   const ADReal rhovAR = UR[THMVACE3D::RHOVA];
      66     2556747 :   const ADReal rhowAR = UR[THMVACE3D::RHOWA];
      67     2556747 :   const ADReal rhoEAR = UR[THMVACE3D::RHOEA];
      68     2556747 :   const ADReal AR = UR[THMVACE3D::AREA];
      69             : 
      70     2556747 :   const auto n_passives = UL.size() - THMVACE3D::N_FLUX_INPUTS;
      71     2556747 :   std::vector<ADReal> passivesL(n_passives, 0.0), passivesR(n_passives, 0.0);
      72     2576341 :   for (const auto i : make_range(n_passives))
      73             :   {
      74       39188 :     passivesL[i] = UL[THMVACE3D::N_FLUX_INPUTS + i] / AL;
      75       19594 :     passivesR[i] = UR[THMVACE3D::N_FLUX_INPUTS + i] / AR;
      76             :   }
      77             : 
      78             :   // compute the primitive variables
      79             : 
      80             :   const ADReal rhoL = rhoAL / AL;
      81           0 :   const ADRealVectorValue uvecL(rhouAL / rhoAL, rhovAL / rhoAL, rhowAL / rhoAL);
      82     2556747 :   const ADReal unL = uvecL * nLR;
      83     2556747 :   const ADReal ut1L = uvecL * t1;
      84     2556747 :   const ADReal ut2L = uvecL * t2;
      85             :   const ADReal rhoEL = rhoEAL / AL;
      86     5113494 :   const ADReal vL = 1.0 / rhoL;
      87             :   const ADReal EL = rhoEAL / rhoAL;
      88     2556747 :   const ADReal eL = EL - 0.5 * uvecL * uvecL;
      89     2556747 :   const ADReal pL = _fp.p_from_v_e(vL, eL);
      90     2556747 :   const ADReal cL = _fp.c_from_v_e(vL, eL);
      91             : 
      92             :   const ADReal rhoR = rhoAR / AR;
      93     2556747 :   const ADRealVectorValue uvecR(rhouAR / rhoAR, rhovAR / rhoAR, rhowAR / rhoAR);
      94     2556747 :   const ADReal unR = uvecR * nLR;
      95     2556747 :   const ADReal ut1R = uvecR * t1;
      96     2556747 :   const ADReal ut2R = uvecR * t2;
      97             :   const ADReal rhoER = rhoEAR / AR;
      98     5113494 :   const ADReal vR = 1.0 / rhoR;
      99             :   const ADReal ER = rhoEAR / rhoAR;
     100     2556747 :   const ADReal eR = ER - 0.5 * uvecR * uvecR;
     101     2556747 :   const ADReal pR = _fp.p_from_v_e(vR, eR);
     102     2556747 :   const ADReal cR = _fp.c_from_v_e(vR, eR);
     103             : 
     104             :   // compute left and right wave speeds
     105             :   ADReal sL, sR;
     106     2556747 :   if (_wave_speed_formulation == WaveSpeedFormulation::EINFELDT)
     107             :   {
     108             :     // compute Roe-averaged variables
     109     2554767 :     const ADReal sqrt_rhoL = sqrt(rhoL);
     110     2554767 :     const ADReal sqrt_rhoR = sqrt(rhoR);
     111     2554767 :     const ADReal un_roe = (sqrt_rhoL * unL + sqrt_rhoR * unR) / (sqrt_rhoL + sqrt_rhoR);
     112     2554767 :     const ADReal ut1_roe = (sqrt_rhoL * ut1L + sqrt_rhoR * ut1R) / (sqrt_rhoL + sqrt_rhoR);
     113     2554767 :     const ADReal ut2_roe = (sqrt_rhoL * ut2L + sqrt_rhoR * ut2R) / (sqrt_rhoL + sqrt_rhoR);
     114     2554767 :     const ADReal HL = EL + pL / rhoL;
     115     2554767 :     const ADReal HR = ER + pR / rhoR;
     116     2554767 :     const ADReal H_roe = (sqrt_rhoL * HL + sqrt_rhoR * HR) / (sqrt_rhoL + sqrt_rhoR);
     117             :     const ADRealVectorValue uvec_roe(un_roe, ut1_roe, ut2_roe);
     118     5109534 :     const ADReal h_roe = H_roe - 0.5 * uvec_roe * uvec_roe;
     119     2554767 :     const ADReal rho_roe = sqrt(rhoL * rhoR);
     120     2554767 :     const ADReal v_roe = 1.0 / rho_roe;
     121     2554767 :     const ADReal e_roe = _fp.e_from_v_h(v_roe, h_roe);
     122     2554767 :     const ADReal c_roe = _fp.c_from_v_e(v_roe, e_roe);
     123             : 
     124     2554767 :     sL = min(unL - cL, un_roe - c_roe);
     125     2554767 :     sR = max(unR + cR, un_roe + c_roe);
     126             :   }
     127        1980 :   else if (_wave_speed_formulation == WaveSpeedFormulation::DAVIS)
     128             :   {
     129        1980 :     sL = min(unL - cL, unR - cR);
     130        1980 :     sR = max(unL + cL, unR + cR);
     131             :   }
     132             :   else
     133             :   {
     134             :     mooseAssert(false, "Invalid 'wave_speed_formulation'.");
     135             :   }
     136             : 
     137             :   // compute middle wave speed
     138     2556747 :   const ADReal sm = (rhoR * unR * (sR - unR) - rhoL * unL * (sL - unL) + pL - pR) /
     139     2556747 :                     (rhoR * (sR - unR) - rhoL * (sL - unL));
     140             : 
     141             :   // compute Omega_L, Omega_R
     142     5113494 :   const ADReal omegL = 1.0 / (sL - sm);
     143     5113494 :   const ADReal omegR = 1.0 / (sR - sm);
     144             : 
     145             :   // compute p^*
     146     2556747 :   const ADReal ps = rhoL * (sL - unL) * (sm - unL) + pL;
     147             : 
     148             :   // compute U_L^*, U_R^*
     149             : 
     150     2556747 :   const ADReal rhoLs = omegL * (sL - unL) * rhoL;
     151     2556747 :   const ADReal rhounLs = omegL * ((sL - unL) * rhoL * unL + ps - pL);
     152     2556747 :   const ADReal rhoELs = omegL * ((sL - unL) * rhoEL - pL * unL + ps * sm);
     153             : 
     154     2556747 :   const ADReal rhoRs = omegR * (sR - unR) * rhoR;
     155     2556747 :   const ADReal rhounRs = omegR * ((sR - unR) * rhoR * unR + ps - pR);
     156     2556747 :   const ADReal rhoERs = omegR * ((sR - unR) * rhoER - pR * unR + ps * sm);
     157             : 
     158     2556747 :   std::vector<ADReal> UL_1d(THMVACE1D::N_FLUX_INPUTS);
     159     2556747 :   UL_1d[THMVACE1D::RHOA] = UL[THMVACE3D::RHOA];
     160     2556747 :   UL_1d[THMVACE1D::RHOUA] = UL[THMVACE3D::RHOUA];
     161     2556747 :   UL_1d[THMVACE1D::RHOEA] = UL[THMVACE3D::RHOEA];
     162     2556747 :   UL_1d[THMVACE1D::AREA] = UL[THMVACE3D::AREA];
     163             : 
     164     2556747 :   std::vector<ADReal> UR_1d(THMVACE1D::N_FLUX_INPUTS);
     165     2556747 :   UR_1d[THMVACE1D::RHOA] = UR[THMVACE3D::RHOA];
     166     2556747 :   UR_1d[THMVACE1D::RHOUA] = UR[THMVACE3D::RHOUA];
     167     2556747 :   UR_1d[THMVACE1D::RHOEA] = UR[THMVACE3D::RHOEA];
     168     2556747 :   UR_1d[THMVACE1D::AREA] = UR[THMVACE3D::AREA];
     169             : 
     170     2556747 :   const ADReal A_flow = computeFlowArea(UL_1d, UR_1d);
     171             : 
     172             :   // compute the fluxes
     173     2556747 :   FL.resize(THMVACE3D::N_FLUX_OUTPUTS + n_passives);
     174     2556747 :   if (sL > 0.0)
     175             :   {
     176        1961 :     FL[THMVACE3D::MASS] = unL * rhoL * A_flow;
     177        1961 :     FL[THMVACE3D::MOM_NORM] = (unL * rhoL * unL + pL) * A_flow;
     178        1961 :     FL[THMVACE3D::MOM_TAN1] = rhoL * unL * ut1L * A_flow;
     179        1961 :     FL[THMVACE3D::MOM_TAN2] = rhoL * unL * ut2L * A_flow;
     180        1961 :     FL[THMVACE3D::ENERGY] = unL * (rhoEL + pL) * A_flow;
     181        1961 :     for (const auto i : make_range(n_passives))
     182           0 :       FL[THMVACE3D::N_FLUX_OUTPUTS + i] = passivesL[i] * unL * A_flow;
     183             : 
     184        1961 :     _last_region_index = 0;
     185             :   }
     186     2554786 :   else if (sL <= 0.0 && sm > 0.0)
     187             :   {
     188     2022683 :     FL[THMVACE3D::MASS] = sm * rhoLs * A_flow;
     189     2022683 :     FL[THMVACE3D::MOM_NORM] = (sm * rhounLs + ps) * A_flow;
     190     2022683 :     FL[THMVACE3D::MOM_TAN1] = rhounLs * ut1L * A_flow;
     191     2022683 :     FL[THMVACE3D::MOM_TAN2] = rhounLs * ut2L * A_flow;
     192     2022683 :     FL[THMVACE3D::ENERGY] = sm * (rhoELs + ps) * A_flow;
     193     2042083 :     for (const auto i : make_range(n_passives))
     194             :     {
     195       19400 :       const auto passiveLs = omegL * (sL - unL) * passivesL[i];
     196       19400 :       FL[THMVACE3D::N_FLUX_OUTPUTS + i] = passiveLs * sm * A_flow;
     197             :     }
     198             : 
     199     2022683 :     _last_region_index = 1;
     200             :   }
     201      532103 :   else if (sm <= 0.0 && sR >= 0.0)
     202             :   {
     203      530114 :     FL[THMVACE3D::MASS] = sm * rhoRs * A_flow;
     204      530114 :     FL[THMVACE3D::MOM_NORM] = (sm * rhounRs + ps) * A_flow;
     205      530114 :     FL[THMVACE3D::MOM_TAN1] = rhounRs * ut1R * A_flow;
     206      530114 :     FL[THMVACE3D::MOM_TAN2] = rhounRs * ut2R * A_flow;
     207      530114 :     FL[THMVACE3D::ENERGY] = sm * (rhoERs + ps) * A_flow;
     208      530308 :     for (const auto i : make_range(n_passives))
     209             :     {
     210         194 :       const auto passiveRs = omegR * (sR - unR) * passivesR[i];
     211         194 :       FL[THMVACE3D::N_FLUX_OUTPUTS + i] = passiveRs * sm * A_flow;
     212             :     }
     213             : 
     214      530114 :     _last_region_index = 2;
     215             :   }
     216        1989 :   else if (sR < 0.0)
     217             :   {
     218        1965 :     FL[THMVACE3D::MASS] = unR * rhoR * A_flow;
     219        1965 :     FL[THMVACE3D::MOM_NORM] = (unR * rhoR * unR + pR) * A_flow;
     220        1965 :     FL[THMVACE3D::MOM_TAN1] = rhoR * unR * ut1R * A_flow;
     221        1965 :     FL[THMVACE3D::MOM_TAN2] = rhoR * unR * ut2R * A_flow;
     222        1965 :     FL[THMVACE3D::ENERGY] = unR * (rhoER + pR) * A_flow;
     223        1965 :     for (const auto i : make_range(n_passives))
     224           0 :       FL[THMVACE3D::N_FLUX_OUTPUTS + i] = passivesR[i] * unR * A_flow;
     225             : 
     226        1965 :     _last_region_index = 3;
     227             :   }
     228             :   else
     229          24 :     std::fill(FL.begin(), FL.end(), getNaN());
     230             : 
     231     2556747 :   FR = FL;
     232             : 
     233             :   const ADReal A_wall_L = AL - A_flow;
     234     2556747 :   FL[THMVACE3D::MOM_NORM] += pL * A_wall_L;
     235             : 
     236             :   const ADReal A_wall_R = AR - A_flow;
     237     2556747 :   FR[THMVACE3D::MOM_NORM] += pR * A_wall_R;
     238     2556747 : }
     239             : 
     240             : ADReal
     241     2556747 : ADNumericalFlux3EqnHLLC::computeFlowArea(const std::vector<ADReal> & UL,
     242             :                                          const std::vector<ADReal> & UR) const
     243             : {
     244     2556747 :   return std::min(UL[THMVACE1D::AREA], UR[THMVACE1D::AREA]);
     245             : }