Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/Riemann/Rusanov.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     Rusanov Riemann flux function
       9                 :            :   \details   This file implements the Rusanov Riemann solver, specific to ALECG.
      10                 :            : */
      11                 :            : // *****************************************************************************
      12                 :            : #ifndef Rusanov_h
      13                 :            : #define Rusanov_h
      14                 :            : 
      15                 :            : #include <vector>
      16                 :            : 
      17                 :            : #include "Fields.hpp"
      18                 :            : #include "FunctionPrototypes.hpp"
      19                 :            : #include "Inciter/Options/Flux.hpp"
      20                 :            : #include "EoS/EOS.hpp"
      21                 :            : 
      22                 :            : namespace inciter {
      23                 :            : 
      24                 :            : //! Rusanov approximate Riemann solver
      25                 :            : struct Rusanov {
      26                 :            : 
      27                 :            :   using real = tk::real;
      28                 :            : 
      29                 :            :   //! Rusanov approximate Riemann solver flux function
      30                 :            :   //! \param[in] mat_blk EOS material block
      31                 :            :   //! \param[in] nx X component of the surface normal
      32                 :            :   //! \param[in] ny Y component of the surface normal
      33                 :            :   //! \param[in] nz Z component of the surface normal
      34                 :            :   //! \param[in] mx X component of the weighted surface normal on chare
      35                 :            :   //!   boundary, weighted by the number of contributions to the edge
      36                 :            :   //! \param[in] my Y component of the weighted surface normal on chare
      37                 :            :   //!   boundary, weighted by the number of contributions to the edge
      38                 :            :   //! \param[in] mz Z component of the weighted surface normal on chare
      39                 :            :   //!   boundary, weighted by the number of contributions to the edge
      40                 :            :   //! \param[in] rL Left density
      41                 :            :   //! \param[in] ruL Left X momentum
      42                 :            :   //! \param[in] rvL Left Y momentum
      43                 :            :   //! \param[in] rwL Left Z momentum
      44                 :            :   //! \param[in] reL Left total specific energy
      45                 :            :   //! \param[in] rR Right density
      46                 :            :   //! \param[in] ruR Right X momentum
      47                 :            :   //! \param[in] rvR Right Y momentum
      48                 :            :   //! \param[in] rwR Right Z momentum
      49                 :            :   //! \param[in] reR Right total specific energy
      50                 :            :   //! \param[in] w1L Left X mesh velocity
      51                 :            :   //! \param[in] w2L Left Y mesh velocity
      52                 :            :   //! \param[in] w3L Left Z mesh velocity
      53                 :            :   //! \param[in] w1R Right X mesh velocity
      54                 :            :   //! \param[in] w2R Right Y mesh velocity
      55                 :            :   //! \param[in] w3R Right Z mesh velocity
      56                 :            :   //! \param[in] pL Left pressure
      57                 :            :   //! \param[in] pR Right pressure
      58                 :            :   //! \param[in,out] fr Riemann solution for density according to Rusanov
      59                 :            :   //! \param[in,out] fru Riemann solution for X momenutm according to Rusanov
      60                 :            :   //! \param[in,out] frv Riemann solution for Y momenutm according to Rusanov
      61                 :            :   //! \param[in,out] frw Riemann solution for Z momenutm according to Rusanov
      62                 :            :   //! \param[in,out] fre Riemann solution for specific total energy according
      63                 :            :   //!   to Rusanov
      64                 :            :   #pragma omp declare simd
      65                 :            :   static void
      66                 :   22970664 :   flux( const std::vector< EOS >& mat_blk,
      67                 :            :         real nx, real ny, real nz,
      68                 :            :         real mx, real my, real mz,
      69                 :            :         real rL, real ruL, real rvL, real rwL, real reL,
      70                 :            :         real rR, real ruR, real rvR, real rwR, real reR,
      71                 :            :         real w1L, real w2L, real w3L, real w1R, real w2R, real w3R,
      72                 :            :         real pL, real pR,
      73                 :            :         real& fr, real& fru, real& frv, real& frw, real& fre )
      74                 :            :   {
      75                 :   22970664 :     auto ul = ruL/rL - w1L;
      76                 :   22970664 :     auto vl = rvL/rL - w2L;
      77                 :   22970664 :     auto wl = rwL/rL - w3L;
      78                 :            : 
      79                 :   22970664 :     auto ur = ruR/rR - w1R;
      80                 :   22970664 :     auto vr = rvR/rR - w2R;
      81                 :   22970664 :     auto wr = rwR/rR - w3R;
      82                 :            : 
      83         [ +  - ]:   22970664 :     auto al = mat_blk[0].compute< EOS::soundspeed >( rL, pL );
      84         [ +  - ]:   22970664 :     auto ar = mat_blk[0].compute< EOS::soundspeed >( rR, pR );
      85                 :            : 
      86                 :            :     // dissipation
      87                 :            :     real len = tk::length( {mx,my,mz} );
      88                 :   22970664 :     real vml = ul*mx + vl*my + wl*mz;
      89         [ +  + ]:   22970664 :     real vmr = ur*mx + vr*my + wr*mz;
      90                 :   22970664 :     auto sl = std::abs(vml) + al*len;
      91         [ +  + ]:   22970664 :     auto sr = std::abs(vmr) + ar*len;
      92                 :   22970664 :     auto smax = std::max( sl, sr );
      93                 :            : 
      94                 :            :     // face-normal velocities
      95                 :   22970664 :     real vnl = ul*nx + vl*ny + wl*nz;
      96                 :   22970664 :     real vnr = ur*nx + vr*ny + wr*nz;
      97                 :            : 
      98                 :            :     // numerical fluxes
      99                 :   22970664 :     fr  = 0.5*(rL*vnl + rR*vnr - smax*(rR - rL));
     100                 :   22970664 :     fru = 0.5*(ruL*vnl + pL*nx + ruR*vnr + pR*nx - smax*(ruR - ruL));
     101                 :   22970664 :     frv = 0.5*(rvL*vnl + pL*ny + rvR*vnr + pR*ny - smax*(rvR - rvL));
     102                 :   22970664 :     frw = 0.5*(rwL*vnl + pL*nz + rwR*vnr + pR*nz - smax*(rwR - rwL));
     103                 :   22970664 :     fre = 0.5*(reL*vnl + reR*vnr
     104                 :   22970664 :                + pL*(ruL*nx + rvL*ny + rwL*nz)/rL
     105                 :   22970664 :                + pR*(ruR*nx + rvR*ny + rwR*nz)/rR
     106                 :   22970664 :                - smax*(reR - reL));
     107                 :   22970664 :   }
     108                 :            : 
     109                 :            :   //! Flux type accessor
     110                 :            :   //! \return Flux type
     111                 :            :   static ctr::FluxType type() noexcept { return ctr::FluxType::Rusanov; }
     112                 :            : };
     113                 :            : 
     114                 :            : } // inciter::
     115                 :            : 
     116                 :            : #endif // Rusanov_h