Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Riemann/HLL.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 Harten-Lax-vanLeer's (HLL) Riemann flux function
9 : : \details This file implements Harten-Lax-vanLeer's (HLL) Riemann solver.
10 : : */
11 : : // *****************************************************************************
12 : : #ifndef HLL_h
13 : : #define HLL_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 : : #include "MultiMat/MultiMatIndexing.hpp"
22 : :
23 : : namespace inciter {
24 : :
25 : : //! HLL approximate Riemann solver
26 : : struct HLL {
27 : :
28 : : //! HLL approximate Riemann solver flux function
29 : : //! \param[in] fn Face/Surface normal
30 : : //! \param[in] u Left and right unknown/state vector
31 : : //! \return Riemann flux solution according to HLL, appended by Riemann
32 : : //! velocities and volume-fractions.
33 : : //! \note The function signature must follow tk::RiemannFluxFn
34 : : static tk::RiemannFluxFn::result_type
35 : 171500 : flux( const std::vector< EOS >& mat_blk,
36 : : const std::array< tk::real, 3 >& fn,
37 : : const std::array< std::vector< tk::real >, 2 >& u,
38 : : const std::vector< std::array< tk::real, 3 > >& )
39 : : {
40 : 171500 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
41 : :
42 : 171500 : auto ncomp = u[0].size()-(3+nmat);
43 [ + - ][ + - ]: 343000 : std::vector< tk::real > flx(ncomp, 0), fl(ncomp, 0), fr(ncomp, 0);
[ + - ]
44 : :
45 : : // Primitive quantities
46 : 171500 : tk::real rhol(0.0), rhor(0.0);
47 : 171500 : tk::real pl(0.0), pr(0.0), amatl(0.0), amatr(0.0), ac_l(0.0), ac_r(0.0);
48 [ + - ][ + - ]: 343000 : std::vector< tk::real > al_l(nmat, 0.0), al_r(nmat, 0.0),
49 [ + - ][ + - ]: 343000 : hml(nmat, 0.0), hmr(nmat, 0.0),
50 [ + - ][ + - ]: 343000 : pml(nmat, 0.0), pmr(nmat, 0.0);
51 [ + + ]: 514500 : for (std::size_t k=0; k<nmat; ++k)
52 : : {
53 : 343000 : al_l[k] = u[0][volfracIdx(nmat, k)];
54 : 343000 : pml[k] = u[0][ncomp+pressureIdx(nmat, k)];
55 : 343000 : pl += pml[k];
56 : 343000 : hml[k] = u[0][energyIdx(nmat, k)] + pml[k];
57 [ + - ]: 686000 : amatl = mat_blk[k].compute< EOS::soundspeed >(
58 : 343000 : u[0][densityIdx(nmat, k)], pml[k], al_l[k], k );
59 : 343000 : rhol += u[0][densityIdx(nmat, k)];
60 : :
61 : 343000 : al_r[k] = u[1][volfracIdx(nmat, k)];
62 : 343000 : pmr[k] = u[1][ncomp+pressureIdx(nmat, k)];
63 : 343000 : pr += pmr[k];
64 : 343000 : hmr[k] = u[1][energyIdx(nmat, k)] + pmr[k];
65 [ + - ]: 686000 : amatr = mat_blk[k].compute< EOS::soundspeed >(
66 : 343000 : u[1][densityIdx(nmat, k)], pmr[k], al_r[k], k );
67 : 343000 : rhor += u[1][densityIdx(nmat, k)];
68 : :
69 : : // Mixture speed of sound
70 : 343000 : ac_l += u[0][densityIdx(nmat, k)] * amatl * amatl;
71 : 343000 : ac_r += u[1][densityIdx(nmat, k)] * amatr * amatr;
72 : : }
73 : :
74 : 171500 : ac_l = std::sqrt(ac_l/rhol);
75 : 171500 : ac_r = std::sqrt(ac_r/rhor);
76 : :
77 : : // Independently limited velocities for advection
78 : 171500 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
79 : 171500 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
80 : 171500 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
81 : 171500 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
82 : 171500 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
83 : 171500 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
84 : :
85 : : // Face-normal velocities from advective velocities
86 : 171500 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
87 : 171500 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
88 : :
89 : : // Flux functions
90 [ + + ]: 514500 : for (std::size_t k=0; k<nmat; ++k)
91 : : {
92 : 343000 : fl[volfracIdx(nmat, k)] = vnl * al_l[k];
93 : 343000 : fl[densityIdx(nmat, k)] = vnl * u[0][densityIdx(nmat, k)];
94 : 343000 : fl[energyIdx(nmat, k)] = vnl * hml[k];
95 : :
96 : 343000 : fr[volfracIdx(nmat, k)] = vnr * al_r[k];
97 : 343000 : fr[densityIdx(nmat, k)] = vnr * u[1][densityIdx(nmat, k)];
98 : 343000 : fr[energyIdx(nmat, k)] = vnr * hmr[k];
99 : : }
100 : :
101 [ + + ]: 686000 : for (std::size_t idir=0; idir<3; ++idir)
102 : : {
103 : 1029000 : fl[momentumIdx(nmat, idir)] = vnl * u[0][momentumIdx(nmat, idir)]
104 : 514500 : + pl*fn[idir];
105 : :
106 : 1029000 : fr[momentumIdx(nmat, idir)] = vnr * u[1][momentumIdx(nmat, idir)]
107 : 514500 : + pr*fn[idir];
108 : : }
109 : :
110 : : // Signal velocities
111 : 171500 : auto Sl = std::min((vnl-ac_l), (vnr-ac_r));
112 : 171500 : auto Sr = std::max((vnl+ac_l), (vnr+ac_r));
113 : :
114 : : // Numerical flux functions and wave-speeds
115 : 171500 : auto c_plus(0.0), c_minus(0.0), p_plus(0.0), p_minus(0.0);
116 [ - + ]: 171500 : if (Sl >= 0.0)
117 : : {
118 [ - - ]: 0 : for (std::size_t k=0; k<flx.size(); ++k)
119 : 0 : flx[k] = fl[k];
120 : 0 : c_plus = vnl;
121 : 0 : p_plus = vnl;
122 : : }
123 [ - + ]: 171500 : else if (Sr <= 0.0)
124 : : {
125 [ - - ]: 0 : for (std::size_t k=0; k<flx.size(); ++k)
126 : 0 : flx[k] = fr[k];
127 : 0 : c_minus = vnr;
128 : 0 : p_minus = vnr;
129 : : }
130 : : else
131 : : {
132 [ + + ]: 1715000 : for (std::size_t k=0; k<flx.size(); ++k)
133 : 1543500 : flx[k] = (Sr*fl[k] - Sl*fr[k] + Sl*Sr*(u[1][k]-u[0][k])) / (Sr-Sl);
134 : 171500 : c_plus = (Sr*vnl - Sr*Sl) / (Sr-Sl);
135 : 171500 : c_minus = (Sr*Sl - Sl*vnr) / (Sr-Sl);
136 : 171500 : p_plus = Sr * vnl / (Sr-Sl);
137 : 171500 : p_minus = -Sl * vnr / (Sr-Sl);
138 : : }
139 : :
140 : 171500 : auto vriem = c_plus+c_minus;
141 : :
142 : 171500 : p_plus = p_plus/( vriem + std::copysign(1.0e-16,vriem) );
143 : 171500 : p_minus = p_minus/( vriem + std::copysign(1.0e-16,vriem) );
144 : :
145 : : // Store Riemann-advected partial pressures
146 [ + + ]: 514500 : for (std::size_t k=0; k<nmat; ++k)
147 [ + - ]: 343000 : flx.push_back(p_plus*pml[k] + p_minus*pmr[k]);
148 : :
149 : : // Store Riemann velocity
150 [ + - ]: 171500 : flx.push_back( vriem );
151 : :
152 [ - + ][ - - ]: 171500 : Assert( flx.size() == (3*nmat+3+nmat+1), "Size of multi-material flux "
[ - - ][ - - ]
153 : : "vector incorrect" );
154 : :
155 : 343000 : return flx;
156 : : }
157 : :
158 : : //! Flux type accessor
159 : : //! \return Flux type
160 : : static ctr::FluxType type() noexcept { return ctr::FluxType::HLL; }
161 : :
162 : : };
163 : :
164 : : } // inciter::
165 : :
166 : : #endif // HLL_h
|
