Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Riemann/HLLDMultiMat.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2023 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief HLLD Riemann flux function for solids
9 : : \details This file implements the HLLD Riemann solver for solids.
10 : : Ref. Barton, Philip T. "An interface-capturing Godunov method for
11 : : the simulation of compressible solid-fluid problems." Journal of
12 : : Computational Physics 390 (2019): 25-50.
13 : : */
14 : : // *****************************************************************************
15 : : #ifndef HLLDMultiMat_h
16 : : #define HLLDMultiMat_h
17 : :
18 : : #include <vector>
19 : :
20 : : #include "Fields.hpp"
21 : : #include "FunctionPrototypes.hpp"
22 : : #include "Inciter/Options/Flux.hpp"
23 : : #include "EoS/EOS.hpp"
24 : : #include "MultiMat/MultiMatIndexing.hpp"
25 : : #include "MultiMat/MiscMultiMatFns.hpp"
26 : :
27 : : namespace inciter {
28 : :
29 : : //! HLLD approximate Riemann solver for solids
30 : : struct HLLDMultiMat {
31 : :
32 : : //! HLLD approximate Riemann solver flux function
33 : : //! \param[in] fn Face/Surface normal
34 : : //! \param[in] u Left and right unknown/state vector
35 : : //! \return Riemann solution according to Harten-Lax-van Leer-Contact
36 : : //! \note The function signature must follow tk::RiemannFluxFn
37 : : static tk::RiemannFluxFn::result_type
38 : 0 : flux( const std::vector< EOS >& mat_blk,
39 : : const std::array< tk::real, 3 >& fn,
40 : : const std::array< std::vector< tk::real >, 2 >& u,
41 : : const std::vector< std::array< tk::real, 3 > >& = {} )
42 : : {
43 : 0 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
44 : 0 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
45 : :
46 [ - - ]: 0 : auto nsld = numSolids(nmat, solidx);
47 : 0 : auto ncomp = u[0].size()-(3+nmat+nsld*6);
48 [ - - ][ - - ]: 0 : std::vector< tk::real > flx(ncomp, 0), fl(ncomp, 0), fr(ncomp, 0),
[ - - ]
49 [ - - ][ - - ]: 0 : ftl(ncomp, 0), ftr(ncomp, 0);
50 : :
51 : : // Primitive variables
52 : : // -------------------------------------------------------------------------
53 : 0 : tk::real rhol(0.0), rhor(0.0);
54 [ - - ]: 0 : for (size_t k=0; k<nmat; ++k) {
55 : 0 : rhol += u[0][densityIdx(nmat, k)];
56 : 0 : rhor += u[1][densityIdx(nmat, k)];
57 : : }
58 : :
59 : 0 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
60 : 0 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
61 : 0 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
62 : 0 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
63 : 0 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
64 : 0 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
65 : :
66 : : // Outer states
67 : : // -------------------------------------------------------------------------
68 : 0 : tk::real pl(0.0), pr(0.0);
69 : 0 : tk::real acl(0.0), acr(0.0);
70 [ - - ][ - - ]: 0 : std::vector< tk::real > apl(nmat, 0.0), apr(nmat, 0.0);
71 : 0 : std::array< tk::real, 3 > Tnl{{0, 0, 0}}, Tnr{{0, 0, 0}};
72 : 0 : std::vector< std::array< tk::real, 3 > > aTnl, aTnr;
73 : : std::array< std::array< tk::real, 3 >, 3 > asigl, asigr;
74 : : std::array< std::array< tk::real, 3 >, 3 >
75 : 0 : signnl{{{0,0,0},{0,0,0},{0,0,0}}};
76 : : std::array< std::array< tk::real, 3 >, 3 >
77 : 0 : signnr{{{0,0,0},{0,0,0},{0,0,0}}};
78 : 0 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gl, gr;
79 : 0 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gnl, gnr;
80 : 0 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > asignnl, asignnr;
81 : :
82 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
83 : : // Left state
84 : 0 : apl[k] = u[0][ncomp+pressureIdx(nmat, k)];
85 : 0 : pl += apl[k];
86 : :
87 : : // inv deformation gradient and Cauchy stress tensors
88 [ - - ][ - - ]: 0 : gl.push_back(getDeformGrad(nmat, k, u[0]));
89 [ - - ]: 0 : asigl = getCauchyStress(nmat, k, ncomp, u[0]);
90 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) asigl[i][i] -= apl[k];
91 : :
92 : : // normal stress (traction) vector
93 [ - - ]: 0 : aTnl.push_back(tk::matvec(asigl, fn));
94 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
95 : 0 : Tnl[i] += aTnl[k][i];
96 : :
97 : : // rotate stress vector
98 [ - - ][ - - ]: 0 : asignnl.push_back(tk::rotateTensor(asigl, fn));
99 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
100 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
101 : 0 : signnl[i][j] += asignnl[k][i][j];
102 : :
103 : : // rotate deformation gradient tensor for speed of sound in normal dir
104 [ - - ][ - - ]: 0 : gnl.push_back(tk::rotateTensor(gl[k], fn));
105 [ - - ]: 0 : auto amatl = mat_blk[k].compute< EOS::soundspeed >(
106 : 0 : u[0][densityIdx(nmat, k)], apl[k],
107 : 0 : u[0][volfracIdx(nmat, k)], k, gnl[k] );
108 : :
109 : : // Right state
110 : 0 : apr[k] = u[1][ncomp+pressureIdx(nmat, k)];
111 : 0 : pr += apr[k];
112 : :
113 : : // inv deformation gradient and Cauchy stress tensors
114 [ - - ][ - - ]: 0 : gr.push_back(getDeformGrad(nmat, k, u[1]));
115 [ - - ]: 0 : asigr = getCauchyStress(nmat, k, ncomp, u[1]);
116 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) asigr[i][i] -= apr[k];
117 : :
118 : : // normal stress (traction) vector
119 [ - - ]: 0 : aTnr.push_back(tk::matvec(asigr, fn));
120 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
121 : 0 : Tnr[i] += aTnr[k][i];
122 : :
123 : : // rotate stress vector
124 [ - - ][ - - ]: 0 : asignnr.push_back(tk::rotateTensor(asigr, fn));
125 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
126 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
127 : 0 : signnr[i][j] += asignnr[k][i][j];
128 : :
129 : : // rotate deformation gradient tensor for speed of sound in normal dir
130 [ - - ][ - - ]: 0 : gnr.push_back(tk::rotateTensor(gr[k], fn));
131 [ - - ]: 0 : auto amatr = mat_blk[k].compute< EOS::soundspeed >(
132 : 0 : u[1][densityIdx(nmat, k)], apr[k],
133 : 0 : u[1][volfracIdx(nmat, k)], k, gnr[k] );
134 : :
135 : : // Mixture speed of sound
136 : 0 : acl += u[0][densityIdx(nmat, k)] * amatl * amatl;
137 : 0 : acr += u[1][densityIdx(nmat, k)] * amatr * amatr;
138 : : }
139 : 0 : acl = std::sqrt(acl/rhol);
140 : 0 : acr = std::sqrt(acr/rhor);
141 : :
142 : : // Rotated velocities from advective velocities
143 : 0 : auto vnl = tk::rotateVector({ul, vl, wl}, fn);
144 : 0 : auto vnr = tk::rotateVector({ur, vr, wr}, fn);
145 : :
146 : : // Signal velocities
147 : 0 : auto Sl = std::min((vnl[0]-acl), (vnr[0]-acr));
148 : 0 : auto Sr = std::max((vnl[0]+acl), (vnr[0]+acr));
149 : 0 : auto Sm = ( rhor*vnr[0]*(Sr-vnr[0]) - rhol*vnl[0]*(Sl-vnl[0])
150 : 0 : - signnl[0][0] + signnr[0][0] )
151 : 0 : /( rhor*(Sr-vnr[0]) - rhol*(Sl-vnl[0]) );
152 : :
153 : : // Next zone inwards (star) variables
154 : : // -------------------------------------------------------------------------
155 : 0 : tk::real acsl(0.0), acsr(0.0);
156 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
157 : 0 : asignnlStar, asignnrStar;
158 : : std::array< std::array< tk::real, 3 >, 3 >
159 : 0 : signnlStar{{{0,0,0},{0,0,0},{0,0,0}}},
160 : 0 : signnrStar{{{0,0,0},{0,0,0},{0,0,0}}};
161 : : std::array< std::array< tk::real, 3 >, 3 >
162 : 0 : siglStar{{{0,0,0},{0,0,0},{0,0,0}}}, sigrStar{{{0,0,0},{0,0,0},{0,0,0}}};
163 [ - - ]: 0 : asignnlStar.resize(nmat);
164 [ - - ]: 0 : asignnrStar.resize(nmat);
165 : 0 : std::array< tk::real, 3 > TnlStar{{0, 0, 0}}, TnrStar{{0, 0, 0}};
166 : 0 : std::vector< std::array< tk::real, 3 > > aTnlStar, aTnrStar;
167 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
168 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
169 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
170 : : {
171 : 0 : asignnlStar[k][i][j] = asignnl[k][i][j];
172 : 0 : asignnrStar[k][i][j] = asignnr[k][i][j];
173 : : }
174 : :
175 [ - - ]: 0 : for (std::size_t i=0; i<1; ++i)
176 : : {
177 : 0 : asignnlStar[k][i][i] -=
178 : 0 : u[0][densityIdx(nmat,k)]*(Sl-vnl[0])*(Sm-vnl[0]);
179 : 0 : asignnrStar[k][i][i] -=
180 : 0 : u[1][densityIdx(nmat,k)]*(Sr-vnr[0])*(Sm-vnr[0]);
181 : : }
182 : :
183 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
184 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
185 : : {
186 : 0 : signnlStar[i][j] += asignnlStar[k][i][j];
187 : 0 : signnrStar[i][j] += asignnrStar[k][i][j];
188 : : }
189 : :
190 [ - - ]: 0 : auto asiglStar = tk::unrotateTensor(asignnlStar[k], fn);
191 [ - - ]: 0 : auto asigrStar = tk::unrotateTensor(asignnrStar[k], fn);
192 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
193 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
194 : : {
195 : 0 : siglStar[i][j] += asiglStar[i][j];
196 : 0 : sigrStar[i][j] += asigrStar[i][j];
197 : : }
198 [ - - ]: 0 : aTnlStar.push_back(tk::matvec(asiglStar, fn));
199 [ - - ]: 0 : aTnrStar.push_back(tk::matvec(asigrStar, fn));
200 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
201 : : {
202 : 0 : TnlStar[i] += aTnlStar[k][i];
203 : 0 : TnrStar[i] += aTnrStar[k][i];
204 : : }
205 : : }
206 : :
207 : 0 : auto w_l = (vnl[0]-Sl)/(Sm-Sl);
208 : 0 : auto w_r = (vnr[0]-Sr)/(Sm-Sr);
209 : :
210 : : std::array< tk::real, 3 > vnlStar, vnrStar;
211 : : // u*_L
212 : 0 : vnlStar[0] = Sm;
213 : 0 : vnlStar[1] = vnl[1];
214 : 0 : vnlStar[2] = vnl[2];
215 : : // u*_R
216 : 0 : vnrStar[0] = Sm;
217 : 0 : vnrStar[1] = vnr[1];
218 : 0 : vnrStar[2] = vnr[2];
219 : :
220 : 0 : auto vlStar = tk::unrotateVector(vnlStar, fn);
221 : 0 : auto vrStar = tk::unrotateVector(vnrStar, fn);
222 : :
223 [ - - ]: 0 : auto uStar = u;
224 : :
225 : 0 : tk::real rholStar(0.0), rhorStar(0.0);
226 : 0 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gnlStar, gnrStar;
227 : 0 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > glStar, grStar;
228 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
229 : : // Left
230 [ - - ]: 0 : if (solidx[k] > 0)
231 : : {
232 [ - - ]: 0 : gnlStar.push_back(gnl[k]);
233 : 0 : gnlStar[k][0][0] *= w_l;
234 : 0 : gnlStar[k][1][0] *= w_l;
235 : 0 : gnlStar[k][2][0] *= w_l;
236 : : // rotate g back to original frame of reference
237 [ - - ][ - - ]: 0 : glStar.push_back(tk::unrotateTensor(gnlStar[k], fn));
238 : : }
239 : 0 : uStar[0][volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)];
240 : 0 : uStar[0][densityIdx(nmat, k)] = w_l * u[0][densityIdx(nmat, k)];
241 : 0 : uStar[0][energyIdx(nmat, k)] = w_l * u[0][energyIdx(nmat, k)]
242 : 0 : + (asignnlStar[k][0][0]*vnlStar[0] - asignnl[k][0][0]*vnl[0]) / (Sm-Sl);
243 : 0 : rholStar += uStar[0][densityIdx(nmat, k)];
244 : :
245 [ - - ]: 0 : auto amatl = mat_blk[k].compute< EOS::shearspeed >(
246 : 0 : uStar[0][densityIdx(nmat, k)], uStar[0][volfracIdx(nmat, k)], k );
247 : :
248 : : // Right
249 [ - - ]: 0 : if (solidx[k] > 0)
250 : : {
251 [ - - ]: 0 : gnrStar.push_back(gnr[k]);
252 : 0 : gnrStar[k][0][0] *= w_r;
253 : 0 : gnrStar[k][1][0] *= w_r;
254 : 0 : gnrStar[k][2][0] *= w_r;
255 : : // rotate g back to original frame of reference
256 [ - - ][ - - ]: 0 : grStar.push_back(tk::unrotateTensor(gnrStar[k], fn));
257 : : }
258 : 0 : uStar[1][volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)];
259 : 0 : uStar[1][densityIdx(nmat, k)] = w_r * u[1][densityIdx(nmat, k)];
260 : 0 : uStar[1][energyIdx(nmat, k)] = w_r * u[1][energyIdx(nmat, k)]
261 : 0 : + (asignnrStar[k][0][0]*vnrStar[0] - asignnr[k][0][0]*vnr[0]) / (Sm-Sr);
262 : 0 : rhorStar += uStar[1][densityIdx(nmat, k)];
263 : :
264 [ - - ]: 0 : auto amatr = mat_blk[k].compute< EOS::shearspeed >(
265 : 0 : uStar[1][densityIdx(nmat, k)], uStar[1][volfracIdx(nmat, k)], k );
266 : :
267 : : // Mixture speed of sound
268 : 0 : acsl += uStar[0][densityIdx(nmat, k)] * amatl * amatl;
269 : 0 : acsr += uStar[1][densityIdx(nmat, k)] * amatr * amatr;
270 : : }
271 : 0 : acsl = std::sqrt(acsl/rholStar);
272 : 0 : acsr = std::sqrt(acsr/rhorStar);
273 : :
274 : : // Signal velocities
275 : 0 : auto Ssl = vnlStar[0]-acsl;
276 : 0 : auto Ssr = vnrStar[0]+acsr;
277 : :
278 : : // Middle-zone (StarStar) variables. Only relevant if solids are present.
279 : : // -------------------------------------------------------------------------
280 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
281 : 0 : asignnlStarStar, asignnrStarStar;
282 : : std::array< std::array< tk::real, 3 >, 3 >
283 : 0 : signnlStarStar{{{0,0,0},{0,0,0},{0,0,0}}},
284 : 0 : signnrStarStar{{{0,0,0},{0,0,0},{0,0,0}}};
285 : : std::array< std::array< tk::real, 3 >, 3 >
286 : 0 : siglStarStar{{{0,0,0},{0,0,0},{0,0,0}}},
287 : 0 : sigrStarStar{{{0,0,0},{0,0,0},{0,0,0}}};
288 [ - - ]: 0 : asignnlStarStar.resize(nmat);
289 [ - - ]: 0 : asignnrStarStar.resize(nmat);
290 : 0 : std::array< tk::real, 3 > TnlStarStar{{0, 0, 0}}, TnrStarStar{{0, 0, 0}};
291 : 0 : std::vector< std::array< tk::real, 3 > > aTnlStarStar, aTnrStarStar;
292 : : std::array< tk::real, 3 > vnlStarStar, vnrStarStar;
293 : : std::array< tk::real, 3 > vlStarStar, vrStarStar;
294 : 0 : tk::real rholStarStar(0.0), rhorStarStar(0.0);
295 : : std::array< std::array< tk::real, 3 >, 3 >
296 : 0 : gnlStarStar{{{0,0,0},{0,0,0},{0,0,0}}};
297 : : std::array< std::array< tk::real, 3 >, 3 >
298 : 0 : gnrStarStar{{{0,0,0},{0,0,0},{0,0,0}}};
299 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
300 : 0 : glStarStar, grStarStar;
301 [ - - ]: 0 : auto uStarStar = uStar;
302 [ - - ]: 0 : if (nsld > 0) {
303 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
304 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
305 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
306 : : {
307 : 0 : asignnlStarStar[k][i][j] = asignnlStar[k][i][j];
308 : 0 : asignnrStarStar[k][i][j] = asignnrStar[k][i][j];
309 : : }
310 : :
311 [ - - ]: 0 : if (solidx[k] > 0)
312 : : {
313 [ - - ]: 0 : for (std::size_t i=1; i<3; ++i)
314 : : {
315 : 0 : asignnlStarStar[k][i][0] =
316 : 0 : ( (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]*asignnr[k][i][0]
317 : 0 : - (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]*asignnl[k][i][0]
318 : 0 : + (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]
319 : 0 : * (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]
320 : 0 : * (vnl[i]-vnr[i]) ) /
321 : 0 : ( (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]
322 : 0 : - (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]);
323 : 0 : asignnrStarStar[k][i][0] =
324 : 0 : ( (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]*asignnr[k][i][0]
325 : 0 : - (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]*asignnl[k][i][0]
326 : 0 : + (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]
327 : 0 : * (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]
328 : 0 : * (vnl[i]-vnr[i]) ) /
329 : 0 : ( (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]
330 : 0 : - (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]);
331 : : }
332 : : // Symmetry
333 : 0 : asignnlStarStar[k][0][1] = asignnlStarStar[k][1][0];
334 : 0 : asignnlStarStar[k][0][2] = asignnlStarStar[k][2][0];
335 : 0 : asignnrStarStar[k][0][1] = asignnrStarStar[k][1][0];
336 : 0 : asignnrStarStar[k][0][2] = asignnrStarStar[k][2][0];
337 : : }
338 : :
339 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
340 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
341 : : {
342 : 0 : signnlStarStar[i][j] += asignnlStarStar[k][i][j];
343 : 0 : signnrStarStar[i][j] += asignnrStarStar[k][i][j];
344 : : }
345 : :
346 [ - - ]: 0 : auto asiglStarStar = tk::unrotateTensor(asignnlStarStar[k], fn);
347 [ - - ]: 0 : auto asigrStarStar = tk::unrotateTensor(asignnrStarStar[k], fn);
348 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
349 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
350 : : {
351 : 0 : siglStarStar[i][j] += asiglStarStar[i][j];
352 : 0 : sigrStarStar[i][j] += asigrStarStar[i][j];
353 : : }
354 [ - - ]: 0 : aTnlStarStar.push_back(tk::matvec(asiglStarStar, fn));
355 [ - - ]: 0 : aTnrStarStar.push_back(tk::matvec(asigrStarStar, fn));
356 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
357 : : {
358 : 0 : TnlStarStar[i] += aTnlStarStar[k][i];
359 : 0 : TnrStarStar[i] += aTnrStarStar[k][i];
360 : : }
361 : : }
362 : :
363 : : // u*_L
364 : 0 : vnlStarStar[0] = Sm;
365 : 0 : vnlStarStar[1] = vnlStar[1]
366 : 0 : + (signnlStarStar[1][0] - signnl[1][0]) / (rholStar*(Sm-Ssl));
367 : 0 : vnlStarStar[2] = vnlStar[2]
368 : 0 : + (signnlStarStar[2][0] - signnl[2][0]) / (rholStar*(Sm-Ssl));
369 : : // u*_R
370 : 0 : vnrStarStar[0] = Sm;
371 : 0 : vnrStarStar[1] = vnrStar[1]
372 : 0 : + (signnrStarStar[1][0] - signnr[1][0]) / (rhorStar*(Sm-Ssr));
373 : 0 : vnrStarStar[2] = vnrStar[2]
374 : 0 : + (signnrStarStar[2][0] - signnr[2][0]) / (rhorStar*(Sm-Ssr));
375 : :
376 : 0 : vlStarStar = tk::unrotateVector(vnlStarStar, fn);
377 : 0 : vrStarStar = tk::unrotateVector(vnrStarStar, fn);
378 : :
379 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
380 : : // Left
381 [ - - ]: 0 : if (solidx[k] > 0)
382 : : {
383 : 0 : gnlStarStar = gnlStar[k];
384 : 0 : gnlStarStar[0][0] += gnlStar[k][0][1]*(vnl[1]-vnlStarStar[1])/(Sm-Ssl)
385 : 0 : + gnlStar[k][0][2]*(vnl[2]-vnlStarStar[2])/(Sm-Ssl);
386 : 0 : gnlStarStar[1][0] += gnlStar[k][1][1]*(vnl[1]-vnlStarStar[1])/(Sm-Ssl)
387 : 0 : + gnlStar[k][1][2]*(vnl[2]-vnlStarStar[2])/(Sm-Ssl);
388 : 0 : gnlStarStar[2][0] += gnlStar[k][2][1]*(vnl[1]-vnlStarStar[1])/(Sm-Ssl)
389 : 0 : + gnlStar[k][2][2]*(vnl[2]-vnlStarStar[2])/(Sm-Ssl);
390 : : // rotate g back to original frame of reference
391 [ - - ][ - - ]: 0 : glStarStar.push_back(tk::unrotateTensor(gnlStarStar, fn));
392 : : }
393 : 0 : uStarStar[0][volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)];
394 : 0 : uStarStar[0][densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)];
395 : 0 : uStarStar[0][energyIdx(nmat, k)] = uStar[0][energyIdx(nmat, k)]
396 : 0 : + ( - asignnl[k][1][0]*vnl[1]
397 : 0 : - asignnl[k][2][0]*vnl[2]
398 : 0 : + asignnlStarStar[k][1][0]*vnlStarStar[1]
399 : 0 : + asignnlStarStar[k][2][0]*vnlStarStar[2]
400 : 0 : ) / (Sm-Ssl);
401 : 0 : rholStarStar += uStarStar[0][densityIdx(nmat, k)];
402 : :
403 : : // Right
404 [ - - ]: 0 : if (solidx[k] > 0)
405 : : {
406 : 0 : gnrStarStar = gnrStar[k];
407 : 0 : gnrStarStar[0][0] += gnrStar[k][0][1]*(vnr[1]-vnrStarStar[1])/(Sm-Ssr)
408 : 0 : + gnrStar[k][0][2]*(vnr[2]-vnrStarStar[2])/(Sm-Ssr);
409 : 0 : gnrStarStar[1][0] += gnrStar[k][1][1]*(vnr[1]-vnrStarStar[1])/(Sm-Ssr)
410 : 0 : + gnrStar[k][1][2]*(vnr[2]-vnrStarStar[2])/(Sm-Ssr);
411 : 0 : gnrStarStar[2][0] += gnrStar[k][2][1]*(vnr[1]-vnrStarStar[1])/(Sm-Ssr)
412 : 0 : + gnrStar[k][2][2]*(vnr[2]-vnrStarStar[2])/(Sm-Ssr);
413 : : // rotate g back to original frame of reference
414 [ - - ][ - - ]: 0 : grStarStar.push_back(tk::unrotateTensor(gnrStarStar, fn));
415 : : }
416 : 0 : uStarStar[1][volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)];
417 : 0 : uStarStar[1][densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)];
418 : 0 : uStarStar[1][energyIdx(nmat, k)] = uStar[1][energyIdx(nmat, k)]
419 : 0 : + ( - asignnr[k][1][0]*vnr[1]
420 : 0 : - asignnr[k][2][0]*vnr[2]
421 : 0 : + asignnrStarStar[k][1][0]*vnrStarStar[1]
422 : 0 : + asignnrStarStar[k][2][0]*vnrStarStar[2]
423 : 0 : ) / (Sm-Ssr);
424 : 0 : rhorStarStar += uStarStar[1][densityIdx(nmat, k)];
425 : : }
426 : : }
427 : :
428 : : // Numerical fluxes
429 : : // -------------------------------------------------------------------------
430 [ - - ]: 0 : if (0.0 <= Sl) {
431 : :
432 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
433 : 0 : flx[momentumIdx(nmat, idir)] =
434 : 0 : u[0][momentumIdx(nmat, idir)] * vnl[0] - Tnl[idir];
435 : :
436 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
437 : 0 : flx[volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)] * vnl[0];
438 : 0 : flx[densityIdx(nmat, k)] = u[0][densityIdx(nmat, k)] * vnl[0];
439 : 0 : flx[energyIdx(nmat, k)] = u[0][energyIdx(nmat, k)] * vnl[0]
440 : 0 : - ul * aTnl[k][0] - vl * aTnl[k][1] - wl * aTnl[k][2];
441 [ - - ]: 0 : if (solidx[k] > 0) {
442 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
443 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
444 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
445 : 0 : gl[k][i][0] * ul +
446 : 0 : gl[k][i][1] * vl +
447 : 0 : gl[k][i][2] * wl ) * fn[j];
448 : : }
449 : : }
450 : :
451 : : // Quantities for non-conservative terms
452 : : // Store Riemann-advected partial pressures
453 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
454 : 0 : flx.push_back(std::sqrt((aTnl[k][0]*aTnl[k][0]
455 : 0 : +aTnl[k][1]*aTnl[k][1]
456 [ - - ]: 0 : +aTnl[k][2]*aTnl[k][2])));
457 : : // Store Riemann velocity
458 [ - - ]: 0 : flx.push_back(vnl[0]);
459 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
460 [ - - ]: 0 : if (solidx[k] > 0) {
461 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
462 [ - - ]: 0 : flx.push_back(aTnl[k][i]);
463 : : }
464 : : }
465 : :
466 : : }
467 : :
468 [ - - ][ - - ]: 0 : else if (Sl < 0.0 && 0.0 <= Ssl) {
469 : :
470 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
471 : 0 : flx[momentumIdx(nmat, idir)] =
472 : 0 : vlStar[idir] * rholStar * Sm - TnlStar[idir];
473 : :
474 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
475 : 0 : flx[volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)] * Sm;
476 : 0 : flx[densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)] * Sm;
477 : 0 : flx[energyIdx(nmat, k)] = uStar[0][energyIdx(nmat, k)] * Sm
478 : 0 : - vlStar[0] * aTnlStar[k][0]
479 : 0 : - vlStar[1] * aTnlStar[k][1]
480 : 0 : - vlStar[2] * aTnlStar[k][2];
481 [ - - ]: 0 : if (solidx[k] > 0) {
482 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
483 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
484 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
485 : 0 : glStar[k][i][0] * vlStar[0] +
486 : 0 : glStar[k][i][1] * vlStar[1] +
487 : 0 : glStar[k][i][2] * vlStar[2] ) * fn[j];
488 : : }
489 : : }
490 : :
491 : : // Quantities for non-conservative terms
492 : : // Store Riemann-advected partial pressures
493 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
494 : 0 : flx.push_back(std::sqrt(aTnlStar[k][0]*aTnlStar[k][0]
495 : 0 : +aTnlStar[k][1]*aTnlStar[k][1]
496 [ - - ]: 0 : +aTnlStar[k][2]*aTnlStar[k][2]));
497 : : // Store Riemann velocity
498 [ - - ]: 0 : flx.push_back(Sm);
499 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
500 [ - - ]: 0 : if (solidx[k] > 0) {
501 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
502 [ - - ]: 0 : flx.push_back(aTnlStar[k][i]);
503 : : }
504 : 0 : }
505 : :
506 : : }
507 : :
508 [ - - ][ - - ]: 0 : else if (Ssl < 0.0 && 0.0 <= Sm && nsld > 0) {
[ - - ]
509 : :
510 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
511 : 0 : flx[momentumIdx(nmat, idir)] =
512 : 0 : vlStarStar[idir] * rholStarStar * Sm - TnlStarStar[idir];
513 : :
514 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
515 : 0 : flx[volfracIdx(nmat, k)] = uStarStar[0][volfracIdx(nmat, k)] * Sm;
516 : 0 : flx[densityIdx(nmat, k)] = uStarStar[0][densityIdx(nmat, k)] * Sm;
517 : 0 : flx[energyIdx(nmat, k)] = uStarStar[0][energyIdx(nmat, k)] * Sm
518 : 0 : - vlStarStar[0] * aTnlStarStar[k][0]
519 : 0 : - vlStarStar[1] * aTnlStarStar[k][1]
520 : 0 : - vlStarStar[2] * aTnlStarStar[k][2];
521 [ - - ]: 0 : if (solidx[k] > 0) {
522 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
523 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
524 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
525 : 0 : glStarStar[k][i][0] * vlStarStar[0] +
526 : 0 : glStarStar[k][i][1] * vlStarStar[1] +
527 : 0 : glStarStar[k][i][2] * vlStarStar[2] ) * fn[j];
528 : : }
529 : : }
530 : :
531 : : // Quantities for non-conservative terms
532 : : // Store Riemann-advected partial pressures
533 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
534 : 0 : flx.push_back(std::sqrt(aTnlStarStar[k][0]*aTnlStarStar[k][0]
535 : 0 : +aTnlStarStar[k][1]*aTnlStarStar[k][1]
536 [ - - ]: 0 : +aTnlStarStar[k][2]*aTnlStarStar[k][2]));
537 : : // Store Riemann velocity
538 [ - - ]: 0 : flx.push_back(Sm);
539 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
540 [ - - ]: 0 : if (solidx[k] > 0) {
541 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
542 [ - - ]: 0 : flx.push_back(aTnlStarStar[k][i]);
543 : : }
544 : 0 : }
545 : :
546 : : }
547 : :
548 [ - - ][ - - ]: 0 : else if (Sm < 0.0 && 0.0 <= Ssr && nsld > 0) {
[ - - ]
549 : :
550 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
551 : 0 : flx[momentumIdx(nmat, idir)] =
552 : 0 : vrStarStar[idir] * rhorStarStar * Sm - TnrStarStar[idir];
553 : :
554 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
555 : 0 : flx[volfracIdx(nmat, k)] = uStarStar[1][volfracIdx(nmat, k)] * Sm;
556 : 0 : flx[densityIdx(nmat, k)] = uStarStar[1][densityIdx(nmat, k)] * Sm;
557 : 0 : flx[energyIdx(nmat, k)] = uStarStar[1][energyIdx(nmat, k)] * Sm
558 : 0 : - vrStarStar[0] * aTnrStarStar[k][0]
559 : 0 : - vrStarStar[1] * aTnrStarStar[k][1]
560 : 0 : - vrStarStar[2] * aTnrStarStar[k][2];
561 [ - - ]: 0 : if (solidx[k] > 0) {
562 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
563 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
564 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
565 : 0 : grStarStar[k][i][0] * vrStarStar[0] +
566 : 0 : grStarStar[k][i][1] * vrStarStar[1] +
567 : 0 : grStarStar[k][i][2] * vrStarStar[2] ) * fn[j];
568 : : }
569 : : }
570 : :
571 : : // Quantities for non-conservative terms
572 : : // Store Riemann-advected partial pressures
573 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
574 : 0 : flx.push_back(std::sqrt(aTnrStarStar[k][0]*aTnrStarStar[k][0]
575 : 0 : +aTnrStarStar[k][1]*aTnrStarStar[k][1]
576 [ - - ]: 0 : +aTnrStarStar[k][2]*aTnrStarStar[k][2]));
577 : : // Store Riemann velocity
578 [ - - ]: 0 : flx.push_back(Sm);
579 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
580 [ - - ]: 0 : if (solidx[k] > 0) {
581 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
582 [ - - ]: 0 : flx.push_back(aTnrStarStar[k][i]);
583 : : }
584 : 0 : }
585 : :
586 : : }
587 : :
588 [ - - ][ - - ]: 0 : else if (Ssr < 0.0 && 0.0 <= Sr) {
589 : :
590 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
591 : 0 : flx[momentumIdx(nmat, idir)] =
592 : 0 : vrStar[idir] * rhorStar * Sm - TnrStar[idir];
593 : :
594 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
595 : 0 : flx[volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)] * Sm;
596 : 0 : flx[densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)] * Sm;
597 : 0 : flx[energyIdx(nmat, k)] = uStar[1][energyIdx(nmat, k)] * Sm
598 : 0 : - vrStar[0] * aTnrStar[k][0]
599 : 0 : - vrStar[1] * aTnrStar[k][1]
600 : 0 : - vrStar[2] * aTnrStar[k][2];
601 [ - - ]: 0 : if (solidx[k] > 0) {
602 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
603 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
604 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
605 : 0 : grStar[k][i][0] * vrStar[0] +
606 : 0 : grStar[k][i][1] * vrStar[1] +
607 : 0 : grStar[k][i][2] * vrStar[2] ) * fn[j];
608 : : }
609 : : }
610 : :
611 : : // Quantities for non-conservative terms
612 : : // Store Riemann-advected partial pressures
613 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
614 : 0 : flx.push_back(std::sqrt(aTnrStar[k][0]*aTnrStar[k][0]
615 : 0 : +aTnrStar[k][1]*aTnrStar[k][1]
616 [ - - ]: 0 : +aTnrStar[k][2]*aTnrStar[k][2]));
617 : : // Store Riemann velocity
618 [ - - ]: 0 : flx.push_back(Sm);
619 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
620 [ - - ]: 0 : if (solidx[k] > 0) {
621 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
622 [ - - ]: 0 : flx.push_back(aTnrStar[k][i]);
623 : : }
624 : 0 : }
625 : :
626 : : }
627 : :
628 : : else {
629 : :
630 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
631 : 0 : flx[momentumIdx(nmat, idir)] =
632 : 0 : u[1][momentumIdx(nmat, idir)] * vnr[0] - Tnr[idir];
633 : :
634 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
635 : 0 : flx[volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)] * vnr[0];
636 : 0 : flx[densityIdx(nmat, k)] = u[1][densityIdx(nmat, k)] * vnr[0];
637 : 0 : flx[energyIdx(nmat, k)] = u[1][energyIdx(nmat, k)] * vnr[0]
638 : 0 : - ur * aTnr[k][0] - vr * aTnr[k][1] - wr * aTnr[k][2];
639 [ - - ]: 0 : if (solidx[k] > 0) {
640 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
641 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
642 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
643 : 0 : gr[k][i][0] * ur +
644 : 0 : gr[k][i][1] * vr +
645 : 0 : gr[k][i][2] * wr ) * fn[j];
646 : : }
647 : : }
648 : :
649 : : // Quantities for non-conservative terms
650 : : // Store Riemann-advected partial pressures
651 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
652 : 0 : flx.push_back(std::sqrt(aTnr[k][0]*aTnr[k][0]
653 : 0 : +aTnr[k][1]*aTnr[k][1]
654 [ - - ]: 0 : +aTnr[k][2]*aTnr[k][2]));
655 : : // Store Riemann velocity
656 [ - - ]: 0 : flx.push_back(vnr[0]);
657 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
658 [ - - ]: 0 : if (solidx[k] > 0) {
659 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
660 [ - - ]: 0 : flx.push_back(aTnr[k][i]);
661 : : }
662 : : }
663 : :
664 : : }
665 : :
666 [ - - ][ - - ]: 0 : Assert( flx.size() == (ncomp+nmat+1+3*nsld), "Size of "
[ - - ][ - - ]
667 : : "multi-material flux vector incorrect" );
668 : :
669 : 0 : return flx;
670 : : }
671 : :
672 : : ////! Flux type accessor
673 : : ////! \return Flux type
674 : : //static ctr::FluxType type() noexcept {
675 : : // return ctr::FluxType::HLLDMultiMat; }
676 : : };
677 : :
678 : : } // inciter::
679 : :
680 : : #endif // HLLDMultiMat_h
|
