Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/Riemann/SplitMachFns.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     Split-Mach functions for the AUSM+ Riemann solver
       9                 :            :   \details   This file implements the split-Mach functions used by the 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 SplitMachFns_h
      16                 :            : #define SplitMachFns_h
      17                 :            : 
      18                 :            : namespace inciter {
      19                 :            : 
      20                 :            : //! Split Mach polynomials for AUSM+ flux
      21                 :            : //! \param[in] mach Local Mach number
      22                 :            : //! \param[in] fa All-speed parameter. Default is 1
      23                 :            : //! \param[in] alpha High-order pressure polynomial coefficient. Default is 3/16
      24                 :            : //! \param[in] beta High-order Mach polynomial coefficient. Default is 1/8
      25                 :            : //! \return Values of the positive and negative split Mach and pressure
      26                 :            : //!   polynomials.
      27                 :            : //! \details This function returns a vector with positive and negative Mach
      28                 :            : //!   and pressure polynomials, as:
      29                 :            : //!   ms[0] = M_4(+),
      30                 :            : //!   ms[1] = M_4(-),
      31                 :            : //!   ms[2] = P_5(+), and
      32                 :            : //!   ms[3] = P_5(-).
      33                 :            : //!   For more details, ref. Liou, M. S. (2006). A sequel to AUSM, Part II:
      34                 :            : //!   AUSM+-up for all speeds. J. Comp. Phys., 214(1), 137-170.
      35                 :   18950082 : static std::array< tk::real, 4 > splitmach_ausm( tk::real mach,
      36                 :            :                                                  tk::real fa = 1.0,
      37                 :            :                                                  tk::real alpha = 3.0/16.0,
      38                 :            :                                                  tk::real beta = 1.0/8.0 )
      39                 :            : {
      40                 :            :   std::array< tk::real, 4 > ms;
      41                 :            : 
      42                 :            :   std::array< tk::real, 3 > msplus, msminus;
      43                 :            :   tk::real psplus, psminus;
      44                 :            : 
      45                 :   18950082 :   msplus[0] = 0.5*(mach + std::fabs(mach));
      46                 :   18950082 :   msminus[0]= 0.5*(mach - std::fabs(mach));
      47                 :            : 
      48                 :   18950082 :   msplus[1] = +0.25*(mach + 1.0)*(mach + 1.0);
      49                 :   18950082 :   msminus[1]= -0.25*(mach - 1.0)*(mach - 1.0);
      50                 :            : 
      51                 :   18950082 :   auto alph_fa = alpha * (-4.0 + 5.0*fa*fa);
      52                 :            : 
      53         [ +  + ]:   18950082 :   if (std::fabs(mach) >= 1.0)
      54                 :            :   {
      55                 :            :       msplus[2] = msplus[0];
      56                 :            :       msminus[2]= msminus[0];
      57                 :     530400 :       psplus    = msplus[0]/mach;
      58                 :     530400 :       psminus   = msminus[0]/mach;
      59                 :            :   }
      60                 :            :   else
      61                 :            :   {
      62                 :   18419682 :       msplus[2] = msplus[1]* (1.0 - 16.0*beta*msminus[1]);
      63                 :   18419682 :       msminus[2]= msminus[1]* (1.0 + 16.0*beta*msplus[1]);
      64                 :   18419682 :       psplus    = msplus[1]*
      65                 :   18419682 :                   ((+2.0 - mach) - (16.0 * alph_fa)*mach*msminus[1]);
      66                 :   18419682 :       psminus   = msminus[1]*
      67                 :   18419682 :                   ((-2.0 - mach) + (16.0 * alph_fa)*mach*msplus[1]);
      68                 :            :   }
      69                 :            : 
      70                 :   18950082 :   ms[0] = msplus[2];
      71                 :   18950082 :   ms[1] = msminus[2];
      72                 :   18950082 :   ms[2] = psplus;
      73                 :   18950082 :   ms[3] = psminus;
      74                 :            : 
      75                 :   18950082 :   return ms;
      76                 :            : }
      77                 :            : 
      78                 :            : } // inciter::
      79                 :            : 
      80                 :            : #endif // SplitMachFns_h