Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/Riemann/AUSMCompFlow.hpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     Advection Upstream Splitting Method (AUSM+) Riemann flux function
       9                 :            :   \details   This file implements the Advection Upstream Splitting Method (AUSM)
      10                 :            :              Riemann solver, with the all-speed corrections.
      11                 :            :              Ref. Liou, M. S. (2006). A sequel to AUSM, Part II: AUSM+-up for
      12                 :            :              all speeds. Journal of computational physics, 214(1), 137-170.
      13                 :            : */
      14                 :            : // *****************************************************************************
      15                 :            : #ifndef AUSMCompFlow_h
      16                 :            : #define AUSMCompFlow_h
      17                 :            : 
      18                 :            : #include <vector>
      19                 :            : 
      20                 :            : #include "Fields.hpp"
      21                 :            : #include "FunctionPrototypes.hpp"
      22                 :            : #include "Inciter/Options/Flux.hpp"
      23                 :            : #include "SplitMachFns.hpp"
      24                 :            : #include "EoS/EOS.hpp"
      25                 :            : 
      26                 :            : namespace inciter {
      27                 :            : 
      28                 :            : //! AUSM+up approximate Riemann solver
      29                 :            : struct AUSMCompFlow {
      30                 :            : 
      31                 :            :   //! AUSM+up approximate Riemann solver flux function
      32                 :            :   //! \param[in] fn Face/Surface normal
      33                 :            :   //! \param[in] u Left and right unknown/state vector
      34                 :            :   //! \return Riemann flux solution according to AUSM+up, appended by Riemann
      35                 :            :   //!   velocities and volume-fractions.
      36                 :            :   //! \note The function signature must follow tk::RiemannFluxFn
      37                 :            :   static tk::RiemannFluxFn::result_type
      38                 :          0 :   flux( const std::vector< EOS >& mat_blk,
      39                 :            :         const std::array< tk::real, 3 >& fn,
      40                 :            :         const std::array< std::vector< tk::real >, 2 >& u,
      41                 :            :         const std::vector< std::array< tk::real, 3 > >& = {} )
      42                 :            :   {
      43                 :            :     // All-speed parameters
      44                 :            :     // These parameters control the amount of all-speed diffusion necessary for
      45                 :            :     // low-Mach flows. Setting k_u and k_p to zero does not add any all-speed
      46                 :            :     // diffusion, whereas setting k_u and k_p to 1 adds maximum recommended
      47                 :            :     // all-speed diffusion. See "Liou, M. S. (2006). A sequel to AUSM, Part II:
      48                 :            :     // AUSM+-up for all speeds. Journal of computational physics, 214(1),
      49                 :            :     // 137-170" for more mathematical explanation. k_u is the velocity diffusion
      50                 :            :     // term and k_p is the pressure diffusion term. These two terms reduce
      51                 :            :     // pressure-velocity decoupling (chequerboarding/odd-even oscillations).
      52                 :          0 :     auto k_u = g_inputdeck.get< tag::lowspeed_ku >();
      53                 :          0 :     auto k_p = g_inputdeck.get< tag::lowspeed_kp >();
      54                 :            : 
      55                 :          0 :     auto ncomp = u[0].size();
      56                 :          0 :     std::vector< tk::real > flx( ncomp, 0 );
      57                 :            : 
      58                 :            :     // Primitive variables
      59                 :          0 :     auto rhol = u[0][0];
      60                 :          0 :     auto ul = u[0][1]/rhol;
      61                 :          0 :     auto vl = u[0][2]/rhol;
      62                 :          0 :     auto wl = u[0][3]/rhol;
      63                 :          0 :     auto rhor = u[1][0];
      64                 :          0 :     auto ur = u[1][1]/rhor;
      65                 :          0 :     auto vr = u[1][2]/rhor;
      66                 :          0 :     auto wr = u[1][3]/rhor;
      67                 :            : 
      68                 :          0 :     tk::real pl(0.0), pr(0.0), amatl(0.0), amatr(0.0);
      69                 :            : 
      70         [ -  - ]:          0 :     pl = mat_blk[0].compute< EOS::pressure >(rhol, ul, vl, wl, u[0][4]);
      71         [ -  - ]:          0 :     amatl = mat_blk[0].compute< EOS::soundspeed >( rhol, pl );
      72                 :            : 
      73         [ -  - ]:          0 :     pr = mat_blk[0].compute< EOS::pressure >(rhor, ur, vr, wr, u[1][4]);
      74         [ -  - ]:          0 :     amatr = mat_blk[0].compute< EOS::soundspeed >( rhor, pr );
      75                 :            : 
      76                 :            :     // Average states for mixture speed of sound
      77                 :            :     auto ac12 = 0.5*(amatl+amatr);
      78                 :          0 :     auto rho12 = 0.5*(rhol+rhor);
      79                 :            : 
      80                 :            :     // Face-normal velocities from advective velocities
      81                 :          0 :     auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
      82                 :          0 :     auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
      83                 :            : 
      84                 :            :     // Mach numbers
      85                 :          0 :     auto ml = vnl/ac12;
      86                 :          0 :     auto mr = vnr/ac12;
      87                 :            : 
      88                 :            :     tk::real f_a(1.0);
      89                 :            : 
      90                 :            :     // Split Mach polynomials
      91                 :          0 :     auto msl = splitmach_ausm( ml, f_a );
      92                 :          0 :     auto msr = splitmach_ausm( mr, f_a );
      93                 :            : 
      94                 :            :     // Riemann Mach number
      95                 :          0 :     auto m0 = 1.0 - (0.5*(vnl*vnl + vnr*vnr)/(ac12*ac12));
      96         [ -  - ]:          0 :     auto mp = -k_p* std::max(m0, 0.0) * (pr-pl) / (f_a*rho12*ac12*ac12);
      97                 :          0 :     auto m12 = msl[0] + msr[1] + mp;
      98                 :          0 :     auto vriem = ac12 * m12;
      99                 :            : 
     100                 :            :     // Riemann pressure
     101                 :          0 :     auto pu = -k_u* msl[2] * msr[3] * f_a * rho12 * ac12 * (vnr-vnl);
     102                 :          0 :     auto p12 = msl[2]*pl + msr[3]*pr + pu;
     103                 :            : 
     104                 :            :     // Flux vector splitting
     105                 :          0 :     auto l_plus = 0.5 * (vriem + std::fabs(vriem));
     106                 :          0 :     auto l_minus = 0.5 * (vriem - std::fabs(vriem));
     107                 :            : 
     108                 :            :     // Conservative fluxes
     109                 :          0 :     flx[0] = l_plus*u[0][0] + l_minus*u[1][0];
     110                 :            : 
     111                 :          0 :     flx[1] = l_plus*u[0][1] + l_minus*u[1][1] + p12*fn[0];
     112                 :          0 :     flx[2] = l_plus*u[0][2] + l_minus*u[1][2] + p12*fn[1];
     113                 :          0 :     flx[3] = l_plus*u[0][3] + l_minus*u[1][3] + p12*fn[2];
     114                 :            : 
     115                 :          0 :     flx[4] = l_plus*(u[0][4] + pl) + l_minus*(u[1][4] + pr);
     116                 :            : 
     117                 :          0 :     return flx;
     118                 :            :   }
     119                 :            : 
     120                 :            :   //! Flux type accessor
     121                 :            :   //! \return Flux type
     122                 :            :   static ctr::FluxType type() noexcept { return ctr::FluxType::AUSM; }
     123                 :            : };
     124                 :            : 
     125                 :            : } // inciter::
     126                 :            : 
     127                 :            : #endif // AUSMCompFlow_h