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 "Types.hpp"
18 : : #include "Fields.hpp"
19 : : #include "Tags.hpp"
20 : : #include "FunctionPrototypes.hpp"
21 : : #include "Inciter/Options/Flux.hpp"
22 : : #include "EoS/EoS.hpp"
23 : : #include "MultiMat/MultiMatIndexing.hpp"
24 : :
25 : : namespace inciter {
26 : :
27 : : //! HLL approximate Riemann solver
28 : : //! \details This class can be used polymorphically with inciter::RiemannSolver
29 : : struct HLL {
30 : :
31 : : //! HLL 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 HLL, 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::array< tk::real, 3 >& fn,
39 : : const std::array< std::vector< tk::real >, 2 >& u,
40 : : const std::vector< std::array< tk::real, 3 > >& )
41 : : {
42 : : const auto nmat =
43 : 0 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[0];
44 : :
45 : 0 : auto ncomp = u[0].size()-(3+nmat);
46 [ - - ][ - - ]: 0 : std::vector< tk::real > flx(ncomp, 0), fl(ncomp, 0), fr(ncomp, 0);
[ - - ][ - - ]
47 : :
48 : : // Primitive variables
49 : : tk::real rhol(0.0), rhor(0.0);
50 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
51 : : {
52 : 0 : rhol += u[0][densityIdx(nmat, k)];
53 : 0 : rhor += u[1][densityIdx(nmat, k)];
54 : : }
55 : :
56 : : tk::real pl(0.0), pr(0.0), amatl(0.0), amatr(0.0), ac_l(0.0), ac_r(0.0);
57 [ - - ][ - - ]: 0 : std::vector< tk::real > al_l(nmat, 0.0), al_r(nmat, 0.0),
[ - - ][ - - ]
58 [ - - ][ - - ]: 0 : hml(nmat, 0.0), hmr(nmat, 0.0),
[ - - ][ - - ]
59 [ - - ][ - - ]: 0 : pml(nmat, 0.0), pmr(nmat, 0.0);
[ - - ][ - - ]
60 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
61 : : {
62 [ - - ]: 0 : al_l[k] = u[0][volfracIdx(nmat, k)];
63 : 0 : pml[k] = u[0][ncomp+pressureIdx(nmat, k)];
64 [ - - ]: 0 : pl += pml[k];
65 : 0 : hml[k] = u[0][energyIdx(nmat, k)] + pml[k];
66 [ - - ]: 0 : amatl = eos_soundspeed< tag::multimat >( 0, u[0][densityIdx(nmat, k)],
67 : : pml[k], al_l[k], k );
68 : :
69 [ - - ]: 0 : al_r[k] = u[1][volfracIdx(nmat, k)];
70 : 0 : pmr[k] = u[1][ncomp+pressureIdx(nmat, k)];
71 : 0 : pr += pmr[k];
72 : 0 : hmr[k] = u[1][energyIdx(nmat, k)] + pmr[k];
73 [ - - ]: 0 : amatr = eos_soundspeed< tag::multimat >( 0, u[1][densityIdx(nmat, k)],
74 : : pmr[k], al_r[k], k );
75 : :
76 : : // Mixture speed of sound
77 : 0 : ac_l += u[0][densityIdx(nmat, k)] * amatl * amatl;
78 : 0 : ac_r += u[1][densityIdx(nmat, k)] * amatr * amatr;
79 : : }
80 : :
81 : 0 : ac_l = std::sqrt(ac_l/rhol);
82 : 0 : ac_r = std::sqrt(ac_r/rhor);
83 : :
84 : : // Independently limited velocities for advection
85 : 0 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
86 : 0 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
87 : 0 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
88 : 0 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
89 : 0 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
90 : 0 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
91 : :
92 : : // Face-normal velocities from advective velocities
93 : 0 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
94 : 0 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
95 : :
96 : : // Flux functions
97 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
98 : : {
99 : 0 : fl[volfracIdx(nmat, k)] = vnl * al_l[k];
100 : 0 : fl[densityIdx(nmat, k)] = vnl * u[0][densityIdx(nmat, k)];
101 : 0 : fl[energyIdx(nmat, k)] = vnl * hml[k];
102 : :
103 : 0 : fr[volfracIdx(nmat, k)] = vnr * al_r[k];
104 : 0 : fr[densityIdx(nmat, k)] = vnr * u[1][densityIdx(nmat, k)];
105 : 0 : fr[energyIdx(nmat, k)] = vnr * hmr[k];
106 : : }
107 : :
108 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
109 : : {
110 : 0 : fl[momentumIdx(nmat, idir)] = vnl * u[0][momentumIdx(nmat, idir)]
111 : 0 : + pl*fn[idir];
112 : :
113 : 0 : fr[momentumIdx(nmat, idir)] = vnr * u[1][momentumIdx(nmat, idir)]
114 : 0 : + pr*fn[idir];
115 : : }
116 : :
117 : : // Signal velocities
118 [ - - ]: 0 : auto Sl = std::min((vnl-ac_l), (vnr-ac_r));
119 [ - - ]: 0 : auto Sr = std::min((vnl+ac_l), (vnr+ac_r));
120 : :
121 : : // Numerical flux functions and wave-speeds
122 : : auto c_plus(0.0), c_minus(0.0), p_plus(0.0), p_minus(0.0);
123 [ - - ]: 0 : if (Sl >= 0.0)
124 : : {
125 [ - - ]: 0 : for (std::size_t k=0; k<flx.size(); ++k)
126 : 0 : flx[k] = fl[k];
127 : : c_plus = vnl;
128 : : p_plus = 1.0;
129 : : }
130 [ - - ]: 0 : else if (Sr <= 0.0)
131 : : {
132 [ - - ]: 0 : for (std::size_t k=0; k<flx.size(); ++k)
133 : 0 : flx[k] = fr[k];
134 : : c_minus = vnr;
135 : : p_minus = 1.0;
136 : : }
137 : : else
138 : : {
139 [ - - ]: 0 : for (std::size_t k=0; k<flx.size(); ++k)
140 : 0 : flx[k] = (Sr*fl[k] - Sl*fr[k] + Sl*Sr*(u[1][k]-u[0][k])) / (Sr-Sl);
141 : 0 : c_plus = (Sr*vnl - Sr*Sl) / (Sr-Sl);
142 : 0 : c_minus = (Sr*Sl - Sl*vnr) / (Sr-Sl);
143 : 0 : p_plus = Sr / (Sr-Sl);
144 : 0 : p_minus = -Sl / (Sr-Sl);
145 : : }
146 : :
147 : 0 : auto vriem = c_plus+c_minus;
148 : :
149 : : c_plus = c_plus/( std::fabs(vriem) + 1.0e-16 );
150 : : c_minus = c_minus/( std::fabs(vriem) + 1.0e-16 );
151 : :
152 : : // Store Riemann-advected partial pressures
153 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
154 [ - - ][ - - ]: 0 : flx.push_back(p_plus*pml[k] + p_minus*pmr[k]);
155 : :
156 : : // Store Riemann velocity
157 [ - - ]: 0 : flx.push_back( vriem );
158 : :
159 : : Assert( flx.size() == (3*nmat+3+nmat+1), "Size of multi-material flux "
160 : : "vector incorrect" );
161 : :
162 : 0 : return flx;
163 : : }
164 : :
165 : : //! Flux type accessor
166 : : //! \return Flux type
167 : : static ctr::FluxType type() noexcept { return ctr::FluxType::HLL; }
168 : :
169 : : };
170 : :
171 : : } // inciter::
172 : :
173 : : #endif // HLL_h
|
