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 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
466 [ - - ]: 0 : if (solidx[k] > 0) {
467 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
468 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
469 [ - - ]: 0 : flx.push_back(gl[k][i][j]);
470 : : }
471 : : }
472 : :
473 : : }
474 : :
475 [ - - ][ - - ]: 0 : else if (Sl < 0.0 && 0.0 <= Ssl) {
476 : :
477 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
478 : 0 : flx[momentumIdx(nmat, idir)] =
479 : 0 : vlStar[idir] * rholStar * Sm - TnlStar[idir];
480 : :
481 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
482 : 0 : flx[volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)] * Sm;
483 : 0 : flx[densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)] * Sm;
484 : 0 : flx[energyIdx(nmat, k)] = uStar[0][energyIdx(nmat, k)] * Sm
485 : 0 : - vlStar[0] * aTnlStar[k][0]
486 : 0 : - vlStar[1] * aTnlStar[k][1]
487 : 0 : - vlStar[2] * aTnlStar[k][2];
488 [ - - ]: 0 : if (solidx[k] > 0) {
489 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
490 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
491 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
492 : 0 : glStar[k][i][0] * vlStar[0] +
493 : 0 : glStar[k][i][1] * vlStar[1] +
494 : 0 : glStar[k][i][2] * vlStar[2] ) * fn[j];
495 : : }
496 : : }
497 : :
498 : : // Quantities for non-conservative terms
499 : : // Store Riemann-advected partial pressures
500 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
501 : 0 : flx.push_back(std::sqrt(aTnlStar[k][0]*aTnlStar[k][0]
502 : 0 : +aTnlStar[k][1]*aTnlStar[k][1]
503 [ - - ]: 0 : +aTnlStar[k][2]*aTnlStar[k][2]));
504 : : // Store Riemann velocity
505 [ - - ]: 0 : flx.push_back(Sm);
506 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
507 [ - - ]: 0 : if (solidx[k] > 0) {
508 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
509 [ - - ]: 0 : flx.push_back(aTnlStar[k][i]);
510 : : }
511 : : }
512 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
513 [ - - ]: 0 : if (solidx[k] > 0) {
514 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
515 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
516 [ - - ]: 0 : flx.push_back(glStar[k][i][j]);
517 : : }
518 : 0 : }
519 : :
520 : : }
521 : :
522 [ - - ][ - - ]: 0 : else if (Ssl < 0.0 && 0.0 <= Sm && nsld > 0) {
[ - - ]
523 : :
524 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
525 : 0 : flx[momentumIdx(nmat, idir)] =
526 : 0 : vlStarStar[idir] * rholStarStar * Sm - TnlStarStar[idir];
527 : :
528 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
529 : 0 : flx[volfracIdx(nmat, k)] = uStarStar[0][volfracIdx(nmat, k)] * Sm;
530 : 0 : flx[densityIdx(nmat, k)] = uStarStar[0][densityIdx(nmat, k)] * Sm;
531 : 0 : flx[energyIdx(nmat, k)] = uStarStar[0][energyIdx(nmat, k)] * Sm
532 : 0 : - vlStarStar[0] * aTnlStarStar[k][0]
533 : 0 : - vlStarStar[1] * aTnlStarStar[k][1]
534 : 0 : - vlStarStar[2] * aTnlStarStar[k][2];
535 [ - - ]: 0 : if (solidx[k] > 0) {
536 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
537 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
538 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
539 : 0 : glStarStar[k][i][0] * vlStarStar[0] +
540 : 0 : glStarStar[k][i][1] * vlStarStar[1] +
541 : 0 : glStarStar[k][i][2] * vlStarStar[2] ) * fn[j];
542 : : }
543 : : }
544 : :
545 : : // Quantities for non-conservative terms
546 : : // Store Riemann-advected partial pressures
547 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
548 : 0 : flx.push_back(std::sqrt(aTnlStarStar[k][0]*aTnlStarStar[k][0]
549 : 0 : +aTnlStarStar[k][1]*aTnlStarStar[k][1]
550 [ - - ]: 0 : +aTnlStarStar[k][2]*aTnlStarStar[k][2]));
551 : : // Store Riemann velocity
552 [ - - ]: 0 : flx.push_back(Sm);
553 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
554 [ - - ]: 0 : if (solidx[k] > 0) {
555 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
556 [ - - ]: 0 : flx.push_back(aTnlStarStar[k][i]);
557 : : }
558 : : }
559 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
560 [ - - ]: 0 : if (solidx[k] > 0) {
561 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
562 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
563 [ - - ]: 0 : flx.push_back(glStarStar[k][i][j]);
564 : : }
565 : 0 : }
566 : :
567 : : }
568 : :
569 [ - - ][ - - ]: 0 : else if (Sm < 0.0 && 0.0 <= Ssr && nsld > 0) {
[ - - ]
570 : :
571 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
572 : 0 : flx[momentumIdx(nmat, idir)] =
573 : 0 : vrStarStar[idir] * rhorStarStar * Sm - TnrStarStar[idir];
574 : :
575 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
576 : 0 : flx[volfracIdx(nmat, k)] = uStarStar[1][volfracIdx(nmat, k)] * Sm;
577 : 0 : flx[densityIdx(nmat, k)] = uStarStar[1][densityIdx(nmat, k)] * Sm;
578 : 0 : flx[energyIdx(nmat, k)] = uStarStar[1][energyIdx(nmat, k)] * Sm
579 : 0 : - vrStarStar[0] * aTnrStarStar[k][0]
580 : 0 : - vrStarStar[1] * aTnrStarStar[k][1]
581 : 0 : - vrStarStar[2] * aTnrStarStar[k][2];
582 [ - - ]: 0 : if (solidx[k] > 0) {
583 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
584 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
585 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
586 : 0 : grStarStar[k][i][0] * vrStarStar[0] +
587 : 0 : grStarStar[k][i][1] * vrStarStar[1] +
588 : 0 : grStarStar[k][i][2] * vrStarStar[2] ) * fn[j];
589 : : }
590 : : }
591 : :
592 : : // Quantities for non-conservative terms
593 : : // Store Riemann-advected partial pressures
594 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
595 : 0 : flx.push_back(std::sqrt(aTnrStarStar[k][0]*aTnrStarStar[k][0]
596 : 0 : +aTnrStarStar[k][1]*aTnrStarStar[k][1]
597 [ - - ]: 0 : +aTnrStarStar[k][2]*aTnrStarStar[k][2]));
598 : : // Store Riemann velocity
599 [ - - ]: 0 : flx.push_back(Sm);
600 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
601 [ - - ]: 0 : if (solidx[k] > 0) {
602 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
603 [ - - ]: 0 : flx.push_back(aTnrStarStar[k][i]);
604 : : }
605 : : }
606 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
607 [ - - ]: 0 : if (solidx[k] > 0) {
608 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
609 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
610 [ - - ]: 0 : flx.push_back(grStarStar[k][i][j]);
611 : : }
612 : 0 : }
613 : :
614 : : }
615 : :
616 [ - - ][ - - ]: 0 : else if (Ssr < 0.0 && 0.0 <= Sr) {
617 : :
618 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
619 : 0 : flx[momentumIdx(nmat, idir)] =
620 : 0 : vrStar[idir] * rhorStar * Sm - TnrStar[idir];
621 : :
622 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
623 : 0 : flx[volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)] * Sm;
624 : 0 : flx[densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)] * Sm;
625 : 0 : flx[energyIdx(nmat, k)] = uStar[1][energyIdx(nmat, k)] * Sm
626 : 0 : - vrStar[0] * aTnrStar[k][0]
627 : 0 : - vrStar[1] * aTnrStar[k][1]
628 : 0 : - vrStar[2] * aTnrStar[k][2];
629 [ - - ]: 0 : if (solidx[k] > 0) {
630 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
631 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
632 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
633 : 0 : grStar[k][i][0] * vrStar[0] +
634 : 0 : grStar[k][i][1] * vrStar[1] +
635 : 0 : grStar[k][i][2] * vrStar[2] ) * fn[j];
636 : : }
637 : : }
638 : :
639 : : // Quantities for non-conservative terms
640 : : // Store Riemann-advected partial pressures
641 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
642 : 0 : flx.push_back(std::sqrt(aTnrStar[k][0]*aTnrStar[k][0]
643 : 0 : +aTnrStar[k][1]*aTnrStar[k][1]
644 [ - - ]: 0 : +aTnrStar[k][2]*aTnrStar[k][2]));
645 : : // Store Riemann velocity
646 [ - - ]: 0 : flx.push_back(Sm);
647 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
648 [ - - ]: 0 : if (solidx[k] > 0) {
649 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
650 [ - - ]: 0 : flx.push_back(aTnrStar[k][i]);
651 : : }
652 : : }
653 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
654 [ - - ]: 0 : if (solidx[k] > 0) {
655 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
656 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
657 [ - - ]: 0 : flx.push_back(grStar[k][i][j]);
658 : : }
659 : 0 : }
660 : :
661 : : }
662 : :
663 : : else {
664 : :
665 [ - - ]: 0 : for (std::size_t idir=0; idir<3; ++idir)
666 : 0 : flx[momentumIdx(nmat, idir)] =
667 : 0 : u[1][momentumIdx(nmat, idir)] * vnr[0] - Tnr[idir];
668 : :
669 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
670 : 0 : flx[volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)] * vnr[0];
671 : 0 : flx[densityIdx(nmat, k)] = u[1][densityIdx(nmat, k)] * vnr[0];
672 : 0 : flx[energyIdx(nmat, k)] = u[1][energyIdx(nmat, k)] * vnr[0]
673 : 0 : - ur * aTnr[k][0] - vr * aTnr[k][1] - wr * aTnr[k][2];
674 [ - - ]: 0 : if (solidx[k] > 0) {
675 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
676 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
677 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
678 : 0 : gr[k][i][0] * ur +
679 : 0 : gr[k][i][1] * vr +
680 : 0 : gr[k][i][2] * wr ) * fn[j];
681 : : }
682 : : }
683 : :
684 : : // Quantities for non-conservative terms
685 : : // Store Riemann-advected partial pressures
686 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
687 : 0 : flx.push_back(std::sqrt(aTnr[k][0]*aTnr[k][0]
688 : 0 : +aTnr[k][1]*aTnr[k][1]
689 [ - - ]: 0 : +aTnr[k][2]*aTnr[k][2]));
690 : : // Store Riemann velocity
691 [ - - ]: 0 : flx.push_back(vnr[0]);
692 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
693 [ - - ]: 0 : if (solidx[k] > 0) {
694 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
695 [ - - ]: 0 : flx.push_back(aTnr[k][i]);
696 : : }
697 : : }
698 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
699 [ - - ]: 0 : if (solidx[k] > 0) {
700 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
701 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
702 [ - - ]: 0 : flx.push_back(gr[k][i][j]);
703 : : }
704 : : }
705 : :
706 : : }
707 : :
708 [ - - ][ - - ]: 0 : Assert( flx.size() == (ncomp+nmat+1+3*nsld), "Size of "
[ - - ][ - - ]
709 : : "multi-material flux vector incorrect" );
710 : :
711 : 0 : return flx;
712 : : }
713 : :
714 : : ////! Flux type accessor
715 : : ////! \return Flux type
716 : : //static ctr::FluxType type() noexcept {
717 : : // return ctr::FluxType::HLLDMultiMat; }
718 : : };
719 : :
720 : : } // inciter::
721 : :
722 : : #endif // HLLDMultiMat_h
|
