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 : 171500 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
42 : :
43 [ + - ]: 171500 : auto nsld = numSolids(nmat, solidx);
44 : 171500 : auto ncomp = u[0].size()-(3+nmat+nsld*6);
45 [ + - ][ + - ]: 343000 : std::vector< tk::real > flx(ncomp, 0), fl(ncomp, 0), fr(ncomp, 0);
[ + - ]
46 : :
47 : : // Primitive quantities
48 : 171500 : tk::real rhol(0.0), rhor(0.0);
49 : 171500 : tk::real amatl(0.0), amatr(0.0), ac_l(0.0), ac_r(0.0);
50 [ + - ][ + - ]: 343000 : std::vector< tk::real > al_l(nmat, 0.0), al_r(nmat, 0.0),
51 [ + - ][ + - ]: 343000 : pml(nmat, 0.0), pmr(nmat, 0.0);
52 : 343000 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > g_l, g_r,
53 : 343000 : asig_l, asig_r;
54 : 343000 : std::vector< std::array< tk::real, 3 > > asign_l, asign_r;
55 : 343000 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gn_l, gn_r;
56 : 171500 : std::array< tk::real, 3 > sign_l {{0, 0, 0}}, sign_r {{0, 0, 0}};
57 [ + + ]: 514500 : for (std::size_t k=0; k<nmat; ++k)
58 : : {
59 : : // Left state
60 : : // -----------------------------------------------------------------------
61 : 343000 : al_l[k] = u[0][volfracIdx(nmat, k)];
62 : 343000 : pml[k] = u[0][ncomp+pressureIdx(nmat, k)];
63 : 343000 : rhol += u[0][densityIdx(nmat, k)];
64 : :
65 : : // inv deformation gradient and Cauchy stress tensors
66 [ + - ][ + - ]: 343000 : g_l.push_back(getDeformGrad(nmat, k, u[0]));
67 [ + - ][ + - ]: 343000 : asig_l.push_back(getCauchyStress(nmat, k, ncomp, u[0]));
68 [ + + ]: 1372000 : for (std::size_t i=0; i<3; ++i) asig_l[k][i][i] -= pml[k];
69 : :
70 : : // normal stress (traction) vector
71 [ + - ]: 343000 : asign_l.push_back(tk::matvec(asig_l[k], fn));
72 [ + + ]: 1372000 : for (std::size_t i=0; i<3; ++i)
73 : 1029000 : sign_l[i] += asign_l[k][i];
74 : :
75 : : // rotate deformation gradient tensor for speed of sound in normal dir
76 [ + - ][ + - ]: 343000 : gn_l.push_back(tk::rotateTensor(g_l[k], fn));
77 [ + - ]: 686000 : amatl = mat_blk[k].compute< EOS::soundspeed >(
78 : 343000 : u[0][densityIdx(nmat, k)], pml[k], al_l[k], k, gn_l[k] );
79 : :
80 : : // Right state
81 : : // -----------------------------------------------------------------------
82 : 343000 : al_r[k] = u[1][volfracIdx(nmat, k)];
83 : 343000 : pmr[k] = u[1][ncomp+pressureIdx(nmat, k)];
84 : 343000 : rhor += u[1][densityIdx(nmat, k)];
85 : :
86 : : // inv deformation gradient and Cauchy stress tensors
87 [ + - ][ + - ]: 343000 : g_r.push_back(getDeformGrad(nmat, k, u[1]));
88 [ + - ][ + - ]: 343000 : asig_r.push_back(getCauchyStress(nmat, k, ncomp, u[1]));
89 [ + + ]: 1372000 : for (std::size_t i=0; i<3; ++i) asig_r[k][i][i] -= pmr[k];
90 : :
91 : : // normal stress (traction) vector
92 [ + - ]: 343000 : asign_r.push_back(tk::matvec(asig_r[k], fn));
93 [ + + ]: 1372000 : for (std::size_t i=0; i<3; ++i)
94 : 1029000 : sign_r[i] += asign_r[k][i];
95 : :
96 : : // rotate deformation gradient tensor for speed of sound in normal dir
97 [ + - ][ + - ]: 343000 : gn_r.push_back(tk::rotateTensor(g_r[k], fn));
98 [ + - ]: 686000 : amatr = mat_blk[k].compute< EOS::soundspeed >(
99 : 343000 : u[1][densityIdx(nmat, k)], pmr[k], al_r[k], k, gn_r[k] );
100 : :
101 : : // Mixture speed of sound
102 : : // -----------------------------------------------------------------------
103 : 343000 : ac_l += u[0][densityIdx(nmat, k)] * amatl * amatl;
104 : 343000 : ac_r += u[1][densityIdx(nmat, k)] * amatr * amatr;
105 : : }
106 : :
107 : 171500 : ac_l = std::sqrt(ac_l/rhol);
108 : 171500 : ac_r = std::sqrt(ac_r/rhor);
109 : :
110 : : // Independently limited velocities for advection
111 : 171500 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
112 : 171500 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
113 : 171500 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
114 : 171500 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
115 : 171500 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
116 : 171500 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
117 : :
118 : : // Face-normal velocities from advective velocities
119 : 171500 : auto vnl = ul*fn[0] + vl*fn[1] + wl*fn[2];
120 : 171500 : auto vnr = ur*fn[0] + vr*fn[1] + wr*fn[2];
121 : :
122 : : // Conservative flux functions
123 : : // -------------------------------------------------------------------------
124 [ + + ]: 514500 : for (std::size_t k=0; k<nmat; ++k)
125 : : {
126 : : // Left fluxes
127 : 343000 : fl[volfracIdx(nmat, k)] = vnl * al_l[k];
128 : 343000 : fl[densityIdx(nmat, k)] = vnl * u[0][densityIdx(nmat, k)];
129 : 343000 : fl[energyIdx(nmat, k)] = vnl * u[0][energyIdx(nmat, k)];
130 [ + + ]: 1372000 : for (std::size_t i=0; i<3; ++i) {
131 : 2058000 : fl[energyIdx(nmat, k)] -= u[0][ncomp+velocityIdx(nmat,i)] *
132 : 1029000 : asign_l[k][i];
133 : : }
134 : :
135 : : // inv deformation gradient tensor
136 [ - + ]: 343000 : if (solidx[k] > 0) {
137 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
138 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
139 : 0 : fl[deformIdx(nmat,solidx[k],i,j)] = (
140 : 0 : g_l[k][i][0] * ul +
141 : 0 : g_l[k][i][1] * vl +
142 : 0 : g_l[k][i][2] * wl ) * fn[j];
143 : : }
144 : :
145 : : // Right fluxes
146 : 343000 : fr[volfracIdx(nmat, k)] = vnr * al_r[k];
147 : 343000 : fr[densityIdx(nmat, k)] = vnr * u[1][densityIdx(nmat, k)];
148 : 343000 : fr[energyIdx(nmat, k)] = vnr * u[1][energyIdx(nmat, k)];
149 [ + + ]: 1372000 : for (std::size_t i=0; i<3; ++i) {
150 : 2058000 : fr[energyIdx(nmat, k)] -= u[1][ncomp+velocityIdx(nmat,i)] *
151 : 1029000 : asign_r[k][i];
152 : : }
153 : :
154 : : // inv deformation gradient tensor
155 [ - + ]: 343000 : if (solidx[k] > 0) {
156 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
157 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
158 : 0 : fr[deformIdx(nmat,solidx[k],i,j)] = (
159 : 0 : g_r[k][i][0] * ur +
160 : 0 : g_r[k][i][1] * vr +
161 : 0 : g_r[k][i][2] * wr ) * fn[j];
162 : : }
163 : : }
164 : :
165 : : // bulk momentum
166 [ + + ]: 686000 : for (std::size_t idir=0; idir<3; ++idir)
167 : : {
168 : 1029000 : fl[momentumIdx(nmat, idir)] = vnl*u[0][momentumIdx(nmat, idir)]
169 : 514500 : - sign_l[idir];
170 : 1029000 : fr[momentumIdx(nmat, idir)] = vnr*u[1][momentumIdx(nmat, idir)]
171 : 514500 : - sign_r[idir];
172 : : }
173 : :
174 : : // Numerical fluxes
175 : : // -------------------------------------------------------------------------
176 : :
177 : : // Signal velocities
178 : 171500 : auto Sl = std::min((vnl-ac_l), (vnr-ac_r));
179 : 171500 : auto Sr = std::max((vnl+ac_l), (vnr+ac_r));
180 : : //// Signal velocities by Einfeldt (HLLE)
181 : : //auto Sl = std::min( 0.0, std::min((vnl-ac_l), 0.5*((vnl+vnr)-(ac_l+ac_r))) );
182 : : //auto Sr = std::max( 0.0, std::max((vnr+ac_r), 0.5*((vnl+vnr)+(ac_l+ac_r))) );
183 : :
184 : : // Numerical flux functions and wave-speeds
185 : 171500 : auto c_plus(0.0), c_minus(0.0), p_plus(0.0), p_minus(0.0);
186 [ - + ]: 171500 : if (Sl >= 0.0)
187 : : {
188 [ - - ]: 0 : flx = fl;
189 : 0 : c_plus = vnl;
190 : 0 : p_plus = 1.0;
191 : : }
192 [ - + ]: 171500 : else if (Sr <= 0.0)
193 : : {
194 [ - - ]: 0 : flx = fr;
195 : 0 : c_minus = vnr;
196 : 0 : p_minus = 1.0;
197 : : }
198 : : else
199 : : {
200 [ + + ]: 1715000 : for (std::size_t k=0; k<flx.size(); ++k)
201 : 1543500 : flx[k] = (Sr*fl[k] - Sl*fr[k] + Sl*Sr*(u[1][k]-u[0][k])) / (Sr-Sl);
202 : 171500 : c_plus = (Sr*vnl - Sr*Sl) / (Sr-Sl);
203 : 171500 : c_minus = (Sr*Sl - Sl*vnr) / (Sr-Sl);
204 : 171500 : p_plus = Sr / (Sr-Sl);
205 : 171500 : p_minus = -Sl / (Sr-Sl);
206 : : }
207 : :
208 : : // Quantities for non-conservative terms
209 : : // -------------------------------------------------------------------------
210 : :
211 : 171500 : auto vriem = c_plus+c_minus;
212 : :
213 : : // Store Riemann-advected partial pressures
214 [ + + ]: 514500 : for (std::size_t k=0; k<nmat; ++k)
215 [ + - ]: 343000 : flx.push_back(p_plus*pml[k] + p_minus*pmr[k]);
216 : :
217 : : // Store Riemann velocity
218 [ + - ]: 171500 : flx.push_back( vriem );
219 : :
220 : : // Store Riemann asign_ij (3*nsld)
221 [ + + ]: 514500 : for (std::size_t k=0; k<nmat; ++k) {
222 [ - + ]: 343000 : if (solidx[k] > 0) {
223 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
224 [ - - ]: 0 : flx.push_back(p_plus*asign_l[k][i] + p_minus*asign_r[k][i]);
225 : : }
226 : : }
227 : :
228 [ - + ][ - - ]: 171500 : Assert( flx.size() == (ncomp+nmat+1+3*nsld), "Size of "
[ - - ][ - - ]
229 : : "multi-material flux vector incorrect" );
230 : :
231 : 343000 : return flx;
232 : : }
233 : :
234 : : //! Flux type accessor
235 : : //! \return Flux type
236 : : static ctr::FluxType type() noexcept { return ctr::FluxType::HLL; }
237 : :
238 : : };
239 : :
240 : : } // inciter::
241 : :
242 : : #endif // HLL_h
|
