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)] = w_l * 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 : + (Sm-vnl[0]) * (u[0][densityIdx(nmat, k)]*Sm
243 : 0 : - asignnlStar[k][0][0]/(Sl-vnl[0]));
244 : 0 : rholStar += uStar[0][densityIdx(nmat, k)];
245 : :
246 [ - - ]: 0 : auto amatl = mat_blk[k].compute< EOS::shearspeed >(
247 : 0 : uStar[0][densityIdx(nmat, k)], uStar[0][volfracIdx(nmat, k)], k );
248 : :
249 : : // Right
250 [ - - ]: 0 : if (solidx[k] > 0)
251 : : {
252 [ - - ]: 0 : gnrStar.push_back(gnr[k]);
253 : 0 : gnrStar[k][0][0] *= w_r;
254 : 0 : gnrStar[k][1][0] *= w_r;
255 : 0 : gnrStar[k][2][0] *= w_r;
256 : : // rotate g back to original frame of reference
257 [ - - ][ - - ]: 0 : grStar.push_back(tk::unrotateTensor(gnrStar[k], fn));
258 : : }
259 : 0 : uStar[1][volfracIdx(nmat, k)] = w_r * u[1][volfracIdx(nmat, k)];
260 : 0 : uStar[1][densityIdx(nmat, k)] = w_r * u[1][densityIdx(nmat, k)];
261 : 0 : uStar[1][energyIdx(nmat, k)] = w_r * u[1][energyIdx(nmat, k)]
262 : 0 : + (Sm-vnr[0]) * (u[1][densityIdx(nmat, k)]*Sm
263 : 0 : - asignnrStar[k][0][0]/(Sr-vnr[0]));
264 : 0 : rhorStar += uStar[1][densityIdx(nmat, k)];
265 : :
266 [ - - ]: 0 : auto amatr = mat_blk[k].compute< EOS::shearspeed >(
267 : 0 : uStar[1][densityIdx(nmat, k)], uStar[1][volfracIdx(nmat, k)], k );
268 : :
269 : : // Mixture speed of sound
270 : 0 : acsl += uStar[0][densityIdx(nmat, k)] * amatl * amatl;
271 : 0 : acsr += uStar[1][densityIdx(nmat, k)] * amatr * amatr;
272 : : }
273 : 0 : acsl = std::sqrt(acsl/rholStar);
274 : 0 : acsr = std::sqrt(acsr/rhorStar);
275 : :
276 : : // Signal velocities
277 : 0 : auto Ssl = vnlStar[0]-acsl;
278 : 0 : auto Ssr = vnrStar[0]+acsr;
279 : :
280 : : // Middle-zone (StarStar) variables. Only relevant if solids are present.
281 : : // -------------------------------------------------------------------------
282 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
283 : 0 : asignnlStarStar, asignnrStarStar;
284 : : std::array< std::array< tk::real, 3 >, 3 >
285 : 0 : signnlStarStar{{{0,0,0},{0,0,0},{0,0,0}}},
286 : 0 : signnrStarStar{{{0,0,0},{0,0,0},{0,0,0}}};
287 : : std::array< std::array< tk::real, 3 >, 3 >
288 : 0 : siglStarStar{{{0,0,0},{0,0,0},{0,0,0}}},
289 : 0 : sigrStarStar{{{0,0,0},{0,0,0},{0,0,0}}};
290 [ - - ]: 0 : asignnlStarStar.resize(nmat);
291 [ - - ]: 0 : asignnrStarStar.resize(nmat);
292 : 0 : std::array< tk::real, 3 > TnlStarStar{{0, 0, 0}}, TnrStarStar{{0, 0, 0}};
293 : 0 : std::vector< std::array< tk::real, 3 > > aTnlStarStar, aTnrStarStar;
294 : : std::array< tk::real, 3 > vnlStarStar, vnrStarStar;
295 : : std::array< tk::real, 3 > vlStarStar, vrStarStar;
296 : 0 : tk::real rholStarStar(0.0), rhorStarStar(0.0);
297 : : std::array< std::array< tk::real, 3 >, 3 >
298 : 0 : gnlStarStar{{{0,0,0},{0,0,0},{0,0,0}}};
299 : : std::array< std::array< tk::real, 3 >, 3 >
300 : 0 : gnrStarStar{{{0,0,0},{0,0,0},{0,0,0}}};
301 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
302 : 0 : glStarStar, grStarStar;
303 [ - - ]: 0 : auto uStarStar = uStar;
304 [ - - ]: 0 : if (nsld > 0) {
305 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
306 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
307 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
308 : : {
309 : 0 : asignnlStarStar[k][i][j] = asignnlStar[k][i][j];
310 : 0 : asignnrStarStar[k][i][j] = asignnrStar[k][i][j];
311 : : }
312 : :
313 [ - - ]: 0 : if (solidx[k] > 0)
314 : : {
315 [ - - ]: 0 : for (std::size_t i=1; i<3; ++i)
316 : : {
317 : 0 : asignnlStarStar[k][i][0] =
318 : 0 : ( (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]*asignnr[k][i][0]
319 : 0 : - (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]*asignnl[k][i][0]
320 : 0 : + (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]
321 : 0 : * (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]
322 : 0 : * (vnl[i]-vnr[i]) ) /
323 : 0 : ( (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]
324 : 0 : - (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]);
325 : 0 : asignnrStarStar[k][i][0] =
326 : 0 : ( (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]*asignnr[k][i][0]
327 : 0 : - (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]*asignnl[k][i][0]
328 : 0 : + (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]
329 : 0 : * (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]
330 : 0 : * (vnl[i]-vnr[i]) ) /
331 : 0 : ( (Sm-Ssl)*uStar[0][densityIdx(nmat,k)]
332 : 0 : - (Sm-Ssr)*uStar[1][densityIdx(nmat,k)]);
333 : : }
334 : : // Symmetry
335 : 0 : asignnlStarStar[k][0][1] = asignnlStarStar[k][1][0];
336 : 0 : asignnlStarStar[k][0][2] = asignnlStarStar[k][2][0];
337 : 0 : asignnrStarStar[k][0][1] = asignnrStarStar[k][1][0];
338 : 0 : asignnrStarStar[k][0][2] = asignnrStarStar[k][2][0];
339 : : }
340 : :
341 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
342 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
343 : : {
344 : 0 : signnlStarStar[i][j] += asignnlStarStar[k][i][j];
345 : 0 : signnrStarStar[i][j] += asignnrStarStar[k][i][j];
346 : : }
347 : :
348 [ - - ]: 0 : auto asiglStarStar = tk::unrotateTensor(asignnlStarStar[k], fn);
349 [ - - ]: 0 : auto asigrStarStar = tk::unrotateTensor(asignnrStarStar[k], fn);
350 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
351 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
352 : : {
353 : 0 : siglStarStar[i][j] += asiglStarStar[i][j];
354 : 0 : sigrStarStar[i][j] += asigrStarStar[i][j];
355 : : }
356 [ - - ]: 0 : aTnlStarStar.push_back(tk::matvec(asiglStarStar, fn));
357 [ - - ]: 0 : aTnrStarStar.push_back(tk::matvec(asigrStarStar, fn));
358 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
359 : : {
360 : 0 : TnlStarStar[i] += aTnlStarStar[k][i];
361 : 0 : TnrStarStar[i] += aTnrStarStar[k][i];
362 : : }
363 : : }
364 : :
365 : : // u*_L
366 : 0 : vnlStarStar[0] = Sm;
367 : 0 : vnlStarStar[1] = vnlStar[1]
368 : 0 : + (signnlStarStar[1][0] - signnl[1][0]) / (rholStar*(Sm-Ssl));
369 : 0 : vnlStarStar[2] = vnlStar[2]
370 : 0 : + (signnlStarStar[2][0] - signnl[2][0]) / (rholStar*(Sm-Ssl));
371 : : // u*_R
372 : 0 : vnrStarStar[0] = Sm;
373 : 0 : vnrStarStar[1] = vnrStar[1]
374 : 0 : + (signnrStarStar[1][0] - signnr[1][0]) / (rhorStar*(Sm-Ssr));
375 : 0 : vnrStarStar[2] = vnrStar[2]
376 : 0 : + (signnrStarStar[2][0] - signnr[2][0]) / (rhorStar*(Sm-Ssr));
377 : :
378 : 0 : vlStarStar = tk::unrotateVector(vnlStarStar, fn);
379 : 0 : vrStarStar = tk::unrotateVector(vnrStarStar, fn);
380 : :
381 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
382 : : // Left
383 [ - - ]: 0 : if (solidx[k] > 0)
384 : : {
385 : 0 : gnlStarStar = gnlStar[k];
386 : 0 : gnlStar[k][0][0] += gnlStar[k][0][1]*(vnlStarStar[1]-vnl[1])/(Sm-Ssl)
387 : 0 : + gnlStar[k][0][2]*(vnlStarStar[2]-vnl[2])/(Sm-Ssl);
388 : 0 : gnlStar[k][1][0] += gnlStar[k][1][1]*(vnlStarStar[1]-vnl[1])/(Sm-Ssl)
389 : 0 : + gnlStar[k][1][2]*(vnlStarStar[2]-vnl[2])/(Sm-Ssl);
390 : 0 : gnlStar[k][2][0] += gnlStar[k][2][1]*(vnlStarStar[1]-vnl[1])/(Sm-Ssl)
391 : 0 : + gnlStar[k][2][2]*(vnlStarStar[2]-vnl[2])/(Sm-Ssl);
392 : : // rotate g back to original frame of reference
393 [ - - ][ - - ]: 0 : glStarStar.push_back(tk::unrotateTensor(gnlStarStar, fn));
394 : : }
395 : 0 : uStarStar[0][volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)];
396 : 0 : uStarStar[0][densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)];
397 : 0 : uStarStar[0][energyIdx(nmat, k)] = uStar[0][energyIdx(nmat, k)]
398 : 0 : + ( - asignnl[k][1][0]*vnl[1]
399 : 0 : - asignnl[k][2][0]*vnl[2]
400 : 0 : + asignnlStarStar[k][1][0]*vnlStarStar[1]
401 : 0 : + asignnlStarStar[k][2][0]*vnlStarStar[2]
402 : 0 : ) / (Sm-Ssl);
403 : 0 : rholStarStar += uStarStar[0][densityIdx(nmat, k)];
404 : :
405 : : // Right
406 [ - - ]: 0 : if (solidx[k] > 0)
407 : : {
408 : 0 : gnrStarStar = gnrStar[k];
409 : 0 : gnrStar[k][0][0] += gnrStar[k][0][1]*(vnrStarStar[1]-vnr[1])/(Sm-Ssr)
410 : 0 : + gnrStar[k][0][2]*(vnrStarStar[2]-vnr[2])/(Sm-Ssr);
411 : 0 : gnrStar[k][1][0] += gnrStar[k][1][1]*(vnrStarStar[1]-vnr[1])/(Sm-Ssr)
412 : 0 : + gnrStar[k][1][2]*(vnrStarStar[2]-vnr[2])/(Sm-Ssr);
413 : 0 : gnrStar[k][2][0] += gnrStar[k][2][1]*(vnrStarStar[1]-vnr[1])/(Sm-Ssr)
414 : 0 : + gnrStar[k][2][2]*(vnrStarStar[2]-vnr[2])/(Sm-Ssr);
415 : : // rotate g back to original frame of reference
416 [ - - ][ - - ]: 0 : grStarStar.push_back(tk::unrotateTensor(gnrStarStar, fn));
417 : : }
418 : 0 : uStarStar[1][volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)];
419 : 0 : uStarStar[1][densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)];
420 : 0 : uStarStar[1][energyIdx(nmat, k)] = uStar[1][energyIdx(nmat, k)]
421 : 0 : + ( - asignnr[k][1][0]*vnr[1]
422 : 0 : - asignnr[k][2][0]*vnr[2]
423 : 0 : + asignnrStarStar[k][1][0]*vnrStarStar[1]
424 : 0 : + asignnrStarStar[k][2][0]*vnrStarStar[2]
425 : 0 : ) / (Sm-Ssr);
426 : 0 : rhorStarStar += uStarStar[1][densityIdx(nmat, k)];
427 : : }
428 : : }
429 : :
430 : : // Numerical fluxes
431 : : // -------------------------------------------------------------------------
432 [ - - ]: 0 : if (Sl > 0.0) {
433 : :
434 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
435 : 0 : flx[momentumIdx(nmat, idir)] =
436 : 0 : u[0][momentumIdx(nmat, idir)] * vnl[0] - Tnl[idir];
437 : :
438 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
439 : 0 : flx[volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)] * vnl[0];
440 : 0 : flx[densityIdx(nmat, k)] = u[0][densityIdx(nmat, k)] * vnl[0];
441 : 0 : flx[energyIdx(nmat, k)] = u[0][energyIdx(nmat, k)] * vnl[0];
442 : 0 : flx[energyIdx(nmat, k)] -= ul * aTnl[k][0];
443 : 0 : flx[energyIdx(nmat, k)] -= vl * aTnl[k][1];
444 : 0 : flx[energyIdx(nmat, k)] -= wl * aTnl[k][2];
445 [ - - ]: 0 : if (solidx[k] > 0) {
446 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
447 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
448 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
449 : 0 : gl[k][i][0] * ul +
450 : 0 : gl[k][i][1] * vl +
451 : 0 : gl[k][i][2] * wl ) * fn[j];
452 : : }
453 : : }
454 : :
455 : : // Quantities for non-conservative terms
456 : : // Store Riemann-advected partial pressures
457 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
458 : 0 : flx.push_back(std::sqrt((aTnl[k][0]*aTnl[k][0]
459 : 0 : +aTnl[k][1]*aTnl[k][1]
460 [ - - ]: 0 : +aTnl[k][2]*aTnl[k][2])));
461 : : // Store Riemann velocity
462 [ - - ]: 0 : flx.push_back(vnl[0]);
463 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
464 [ - - ]: 0 : if (solidx[k] > 0) {
465 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
466 [ - - ]: 0 : flx.push_back(aTnl[k][i]);
467 : : }
468 : : }
469 : :
470 : : }
471 : :
472 [ - - ][ - - ]: 0 : else if (Sl <= 0.0 && Ssl > 0.0) {
473 : :
474 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
475 : 0 : flx[momentumIdx(nmat, idir)] =
476 : 0 : vlStar[idir] * rholStar * Sm - TnlStar[idir];
477 : :
478 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
479 : 0 : flx[volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)] * Sm;
480 : 0 : flx[densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)] * Sm;
481 : 0 : flx[energyIdx(nmat, k)] = uStar[0][energyIdx(nmat, k)] * Sm;
482 : 0 : flx[energyIdx(nmat, k)] -= vlStar[0] * aTnlStar[k][0];
483 : 0 : flx[energyIdx(nmat, k)] -= vlStar[1] * aTnlStar[k][1];
484 : 0 : flx[energyIdx(nmat, k)] -= vlStar[2] * aTnlStar[k][2];
485 [ - - ]: 0 : if (solidx[k] > 0) {
486 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
487 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
488 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
489 : 0 : glStar[k][i][0] * vlStar[0] +
490 : 0 : glStar[k][i][1] * vlStar[1] +
491 : 0 : glStar[k][i][2] * vlStar[2] ) * fn[j];
492 : : }
493 : : }
494 : :
495 : : // Quantities for non-conservative terms
496 : : // Store Riemann-advected partial pressures
497 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
498 : 0 : flx.push_back(std::sqrt(aTnlStar[k][0]*aTnlStar[k][0]
499 : 0 : +aTnlStar[k][1]*aTnlStar[k][1]
500 [ - - ]: 0 : +aTnlStar[k][2]*aTnlStar[k][2]));
501 : : // Store Riemann velocity
502 [ - - ]: 0 : flx.push_back(Sm);
503 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
504 [ - - ]: 0 : if (solidx[k] > 0) {
505 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
506 [ - - ]: 0 : flx.push_back(aTnlStar[k][i]);
507 : : }
508 : 0 : }
509 : :
510 : : }
511 : :
512 [ - - ][ - - ]: 0 : else if (Ssl <= 0.0 && Sm > 0.0) {
513 : :
514 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
515 : 0 : flx[momentumIdx(nmat, idir)] =
516 : 0 : vlStarStar[idir] * rholStarStar * Sm - TnlStarStar[idir];
517 : :
518 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
519 : 0 : flx[volfracIdx(nmat, k)] = uStarStar[0][volfracIdx(nmat, k)] * Sm;
520 : 0 : flx[densityIdx(nmat, k)] = uStarStar[0][densityIdx(nmat, k)] * Sm;
521 : 0 : flx[energyIdx(nmat, k)] = uStarStar[0][energyIdx(nmat, k)] * Sm;
522 : 0 : flx[energyIdx(nmat, k)] -= vlStarStar[0] * aTnlStarStar[k][0];
523 : 0 : flx[energyIdx(nmat, k)] -= vlStarStar[1] * aTnlStarStar[k][1];
524 : 0 : flx[energyIdx(nmat, k)] -= vlStarStar[2] * aTnlStarStar[k][2];
525 [ - - ]: 0 : if (solidx[k] > 0) {
526 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
527 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
528 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
529 : 0 : glStarStar[k][i][0] * vlStarStar[0] +
530 : 0 : glStarStar[k][i][1] * vlStarStar[1] +
531 : 0 : glStarStar[k][i][2] * vlStarStar[2] ) * fn[j];
532 : : }
533 : : }
534 : :
535 : : // Quantities for non-conservative terms
536 : : // Store Riemann-advected partial pressures
537 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
538 : 0 : flx.push_back(std::sqrt(aTnlStarStar[k][0]*aTnlStarStar[k][0]
539 : 0 : +aTnlStarStar[k][1]*aTnlStarStar[k][1]
540 [ - - ]: 0 : +aTnlStarStar[k][2]*aTnlStarStar[k][2]));
541 : : // Store Riemann velocity
542 [ - - ]: 0 : flx.push_back(Sm);
543 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
544 [ - - ]: 0 : if (solidx[k] > 0) {
545 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
546 [ - - ]: 0 : flx.push_back(aTnlStarStar[k][i]);
547 : : }
548 : 0 : }
549 : :
550 : : }
551 : :
552 [ - - ][ - - ]: 0 : else if (Sm <= 0.0 && Ssr > 0.0) {
553 : :
554 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
555 : 0 : flx[momentumIdx(nmat, idir)] =
556 : 0 : vrStarStar[idir] * rhorStarStar * Sm - TnrStarStar[idir];
557 : :
558 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
559 : 0 : flx[volfracIdx(nmat, k)] = uStarStar[1][volfracIdx(nmat, k)] * Sm;
560 : 0 : flx[densityIdx(nmat, k)] = uStarStar[1][densityIdx(nmat, k)] * Sm;
561 : 0 : flx[energyIdx(nmat, k)] = uStarStar[1][energyIdx(nmat, k)] * Sm;
562 : 0 : flx[energyIdx(nmat, k)] -= vrStarStar[0] * aTnrStarStar[k][0];
563 : 0 : flx[energyIdx(nmat, k)] -= vrStarStar[1] * aTnrStarStar[k][1];
564 : 0 : flx[energyIdx(nmat, k)] -= vrStarStar[2] * aTnrStarStar[k][2];
565 [ - - ]: 0 : if (solidx[k] > 0) {
566 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
567 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
568 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
569 : 0 : grStarStar[k][i][0] * vrStarStar[0] +
570 : 0 : grStarStar[k][i][1] * vrStarStar[1] +
571 : 0 : grStarStar[k][i][2] * vrStarStar[2] ) * fn[j];
572 : : }
573 : : }
574 : :
575 : : // Quantities for non-conservative terms
576 : : // Store Riemann-advected partial pressures
577 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
578 : 0 : flx.push_back(std::sqrt(aTnrStarStar[k][0]*aTnrStarStar[k][0]
579 : 0 : +aTnrStarStar[k][1]*aTnrStarStar[k][1]
580 [ - - ]: 0 : +aTnrStarStar[k][2]*aTnrStarStar[k][2]));
581 : : // Store Riemann velocity
582 [ - - ]: 0 : flx.push_back(Sm);
583 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
584 [ - - ]: 0 : if (solidx[k] > 0) {
585 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
586 [ - - ]: 0 : flx.push_back(aTnrStarStar[k][i]);
587 : : }
588 : 0 : }
589 : :
590 : : }
591 : :
592 [ - - ][ - - ]: 0 : else if (Ssr <= 0.0 && Sr > 0.0) {
593 : :
594 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
595 : 0 : flx[momentumIdx(nmat, idir)] =
596 : 0 : vrStar[idir] * rhorStar * Sm - TnrStar[idir];
597 : :
598 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
599 : 0 : flx[volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)] * Sm;
600 : 0 : flx[densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)] * Sm;
601 : 0 : flx[energyIdx(nmat, k)] = uStar[1][energyIdx(nmat, k)] * Sm;
602 : 0 : flx[energyIdx(nmat, k)] -= vrStar[0] * aTnrStar[k][0];
603 : 0 : flx[energyIdx(nmat, k)] -= vrStar[1] * aTnrStar[k][1];
604 : 0 : flx[energyIdx(nmat, k)] -= vrStar[2] * aTnrStar[k][2];
605 [ - - ]: 0 : if (solidx[k] > 0) {
606 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
607 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
608 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
609 : 0 : grStar[k][i][0] * vrStar[0] +
610 : 0 : grStar[k][i][1] * vrStar[1] +
611 : 0 : grStar[k][i][2] * vrStar[2] ) * fn[j];
612 : : }
613 : : }
614 : :
615 : : // Quantities for non-conservative terms
616 : : // Store Riemann-advected partial pressures
617 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
618 : 0 : flx.push_back(std::sqrt(aTnrStar[k][0]*aTnrStar[k][0]
619 : 0 : +aTnrStar[k][1]*aTnrStar[k][1]
620 [ - - ]: 0 : +aTnrStar[k][2]*aTnrStar[k][2]));
621 : : // Store Riemann velocity
622 [ - - ]: 0 : flx.push_back(Sm);
623 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
624 [ - - ]: 0 : if (solidx[k] > 0) {
625 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
626 [ - - ]: 0 : flx.push_back(aTnrStar[k][i]);
627 : : }
628 : 0 : }
629 : :
630 : : }
631 : :
632 : : else {
633 : :
634 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
635 : 0 : flx[momentumIdx(nmat, idir)] =
636 : 0 : u[1][momentumIdx(nmat, idir)] * vnr[0] - Tnr[idir];
637 : :
638 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
639 : 0 : flx[volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)] * vnr[0];
640 : 0 : flx[densityIdx(nmat, k)] = u[1][densityIdx(nmat, k)] * vnr[0];
641 : 0 : flx[energyIdx(nmat, k)] = u[1][energyIdx(nmat, k)] * vnr[0];
642 : 0 : flx[energyIdx(nmat, k)] -= ur * aTnr[k][0];
643 : 0 : flx[energyIdx(nmat, k)] -= vr * aTnr[k][1];
644 : 0 : flx[energyIdx(nmat, k)] -= wr * aTnr[k][2];
645 [ - - ]: 0 : if (solidx[k] > 0) {
646 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
647 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
648 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
649 : 0 : gr[k][i][0] * ur +
650 : 0 : gr[k][i][1] * vr +
651 : 0 : gr[k][i][2] * wr ) * fn[j];
652 : : }
653 : : }
654 : :
655 : : // Quantities for non-conservative terms
656 : : // Store Riemann-advected partial pressures
657 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
658 : 0 : flx.push_back(std::sqrt(aTnr[k][0]*aTnr[k][0]
659 : 0 : +aTnr[k][1]*aTnr[k][1]
660 [ - - ]: 0 : +aTnr[k][2]*aTnr[k][2]));
661 : : // Store Riemann velocity
662 [ - - ]: 0 : flx.push_back(vnr[0]);
663 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
664 [ - - ]: 0 : if (solidx[k] > 0) {
665 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
666 [ - - ]: 0 : flx.push_back(aTnr[k][i]);
667 : : }
668 : : }
669 : :
670 : : }
671 : :
672 [ - - ][ - - ]: 0 : Assert( flx.size() == (ncomp+nmat+1+3*nsld), "Size of "
[ - - ][ - - ]
673 : : "multi-material flux vector incorrect" );
674 : :
675 : 0 : return flx;
676 : : }
677 : :
678 : : ////! Flux type accessor
679 : : ////! \return Flux type
680 : : //static ctr::FluxType type() noexcept {
681 : : // return ctr::FluxType::HLLDMultiMat; }
682 : : };
683 : :
684 : : } // inciter::
685 : :
686 : : #endif // HLLDMultiMat_h
|
