Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Riemann/HLLCMultiMat.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 HLLC Riemann flux function for solids
9 : : \details This file implements the HLLC Riemann solver for solids.
10 : : Ref. Ndanou, S., Favrie, N., & Gavrilyuk, S. (2015). Multi-solid
11 : : and multi-fluid diffuse interface model: Applications to dynamic
12 : : fracture and fragmentation. Journal of Computational Physics, 295,
13 : : 523-555.
14 : : */
15 : : // *****************************************************************************
16 : : #ifndef HLLCMultiMat_h
17 : : #define HLLCMultiMat_h
18 : :
19 : : #include <vector>
20 : :
21 : : #include "Fields.hpp"
22 : : #include "FunctionPrototypes.hpp"
23 : : #include "Inciter/Options/Flux.hpp"
24 : : #include "EoS/EOS.hpp"
25 : : #include "MultiMat/MultiMatIndexing.hpp"
26 : : #include "MultiMat/MiscMultiMatFns.hpp"
27 : :
28 : : namespace inciter {
29 : :
30 : : //! HLLC approximate Riemann solver for solids
31 : : struct HLLCMultiMat {
32 : :
33 : : //! HLLC approximate Riemann solver flux function
34 : : //! \param[in] fn Face/Surface normal
35 : : //! \param[in] u Left and right unknown/state vector
36 : : //! \return Riemann solution according to Harten-Lax-van Leer-Contact
37 : : //! \note The function signature must follow tk::RiemannFluxFn
38 : : static tk::RiemannFluxFn::result_type
39 : 868500 : flux( const std::vector< EOS >& mat_blk,
40 : : const std::array< tk::real, 3 >& fn,
41 : : const std::array< std::vector< tk::real >, 2 >& u,
42 : : const std::vector< std::array< tk::real, 3 > >& = {} )
43 : : {
44 : 868500 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
45 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
46 : :
47 : 868500 : auto nsld = numSolids(nmat, solidx);
48 : 868500 : auto ncomp = u[0].size()-(3+nmat+nsld*6);
49 : 868500 : std::vector< tk::real > flx(ncomp, 0);
50 : :
51 : : // Primitive variables
52 : : // -------------------------------------------------------------------------
53 : : tk::real rhol(0.0), rhor(0.0);
54 [ + + ]: 2605500 : for (size_t k=0; k<nmat; ++k) {
55 : 1737000 : rhol += u[0][densityIdx(nmat, k)];
56 : 1737000 : rhor += u[1][densityIdx(nmat, k)];
57 : : }
58 : :
59 [ + - ]: 868500 : auto ul = u[0][ncomp+velocityIdx(nmat, 0)];
60 : 868500 : auto vl = u[0][ncomp+velocityIdx(nmat, 1)];
61 : 868500 : auto wl = u[0][ncomp+velocityIdx(nmat, 2)];
62 : 868500 : auto ur = u[1][ncomp+velocityIdx(nmat, 0)];
63 : 868500 : auto vr = u[1][ncomp+velocityIdx(nmat, 1)];
64 : 868500 : auto wr = u[1][ncomp+velocityIdx(nmat, 2)];
65 : :
66 : : // Outer states
67 : : // -------------------------------------------------------------------------
68 : : tk::real pl(0.0), pr(0.0);
69 : : tk::real acl(0.0), acr(0.0);
70 [ + - ][ + - ]: 868500 : std::vector< tk::real > apl(nmat, 0.0), apr(nmat, 0.0);
[ - - ][ - - ]
71 : 868500 : std::array< tk::real, 3 > Tnl{{0, 0, 0}}, Tnr{{0, 0, 0}};
72 : : 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 : 868500 : signnl{{{0,0,0},{0,0,0},{0,0,0}}};
76 : : std::array< std::array< tk::real, 3 >, 3 >
77 : 868500 : signnr{{{0,0,0},{0,0,0},{0,0,0}}};
78 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gl, gr;
79 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gnl, gnr;
80 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > > asignnl, asignnr;
81 : :
82 [ + + ]: 2605500 : for (std::size_t k=0; k<nmat; ++k) {
83 : : // Left state
84 [ + - ]: 1737000 : apl[k] = u[0][ncomp+pressureIdx(nmat, k)];
85 : : pl += apl[k];
86 : :
87 : : // inv deformation gradient and Cauchy stress tensors
88 [ + - ]: 1737000 : gl.push_back(getDeformGrad(nmat, k, u[0]));
89 [ + - ]: 1737000 : asigl = getCauchyStress(nmat, k, ncomp, u[0]);
90 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i) asigl[i][i] -= apl[k];
91 : :
92 : : // normal stress (traction) vector
93 : 1737000 : aTnl.push_back(tk::matvec(asigl, fn));
94 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i)
95 : 5211000 : Tnl[i] += aTnl[k][i];
96 : :
97 : : // rotate stress vector
98 [ + - ]: 3474000 : asignnl.push_back(tk::rotateTensor(asigl, fn));
99 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i)
100 [ + + ]: 20844000 : for (std::size_t j=0; j<3; ++j)
101 : 15633000 : signnl[i][j] += asignnl[k][i][j];
102 : :
103 : : // rotate deformation gradient tensor for speed of sound in normal dir
104 [ + - ]: 1737000 : gnl.push_back(tk::rotateTensor(gl[k], fn));
105 [ + - ]: 1737000 : auto amatl = mat_blk[k].compute< EOS::soundspeed >(
106 [ + - ]: 1737000 : u[0][densityIdx(nmat, k)], apl[k],
107 [ + - ]: 1737000 : u[0][volfracIdx(nmat, k)], k, gnl[k] );
108 : :
109 : : // Right state
110 [ + - ]: 1737000 : apr[k] = u[1][ncomp+pressureIdx(nmat, k)];
111 : : pr += apr[k];
112 : :
113 : : // inv deformation gradient and Cauchy stress tensors
114 [ + - ]: 1737000 : gr.push_back(getDeformGrad(nmat, k, u[1]));
115 [ + - ]: 1737000 : asigr = getCauchyStress(nmat, k, ncomp, u[1]);
116 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i) asigr[i][i] -= apr[k];
117 : :
118 : : // normal stress (traction) vector
119 : 1737000 : aTnr.push_back(tk::matvec(asigr, fn));
120 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i)
121 : 5211000 : Tnr[i] += aTnr[k][i];
122 : :
123 : : // rotate stress vector
124 [ + - ]: 3474000 : asignnr.push_back(tk::rotateTensor(asigr, fn));
125 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i)
126 [ + + ]: 20844000 : for (std::size_t j=0; j<3; ++j)
127 : 15633000 : signnr[i][j] += asignnr[k][i][j];
128 : :
129 : : // rotate deformation gradient tensor for speed of sound in normal dir
130 [ + - ][ - - ]: 1737000 : gnr.push_back(tk::rotateTensor(gr[k], fn));
131 [ + - ]: 1737000 : auto amatr = mat_blk[k].compute< EOS::soundspeed >(
132 [ + - ]: 1737000 : u[1][densityIdx(nmat, k)], apr[k],
133 [ + - ]: 1737000 : u[1][volfracIdx(nmat, k)], k, gnr[k] );
134 : :
135 : : // Mixture speed of sound
136 : 1737000 : acl += u[0][densityIdx(nmat, k)] * amatl * amatl;
137 : 1737000 : acr += u[1][densityIdx(nmat, k)] * amatr * amatr;
138 : : }
139 : 868500 : acl = std::sqrt(acl/rhol);
140 : 868500 : acr = std::sqrt(acr/rhor);
141 : :
142 : : // Rotated velocities from advective velocities
143 : 868500 : auto vnl = tk::rotateVector({ul, vl, wl}, fn);
144 : 868500 : auto vnr = tk::rotateVector({ur, vr, wr}, fn);
145 : :
146 : : // Signal velocities
147 [ + + ]: 868500 : auto Sl = std::min((vnl[0]-acl), (vnr[0]-acr));
148 [ + + ]: 868500 : auto Sr = std::max((vnl[0]+acl), (vnr[0]+acr));
149 : 868500 : auto Sm = ( rhor*vnr[0]*(Sr-vnr[0]) - rhol*vnl[0]*(Sl-vnl[0])
150 : 868500 : - signnl[0][0] + signnr[0][0] )
151 [ + - ]: 868500 : /( rhor*(Sr-vnr[0]) - rhol*(Sl-vnl[0]) );
152 : :
153 : : // Middle-zone (star) variables
154 : : // -------------------------------------------------------------------------
155 : : // the stress in the star zone is theoretically equivalent when derived
156 : : // from the left or right state. However, there might be differences
157 : : // numerically due to truncation etc. Hence two separate evaluations
158 : : // are used.
159 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
160 : : asignnlStar, asignnrStar;
161 : : std::array< std::array< tk::real, 3 >, 3 >
162 : 868500 : signnlStar{{{0,0,0},{0,0,0},{0,0,0}}},
163 : 868500 : signnrStar{{{0,0,0},{0,0,0},{0,0,0}}};
164 : : std::array< std::array< tk::real, 3 >, 3 >
165 : 868500 : siglStar{{{0,0,0},{0,0,0},{0,0,0}}}, sigrStar{{{0,0,0},{0,0,0},{0,0,0}}};
166 [ + - ]: 868500 : asignnlStar.resize(nmat);
167 [ + - ]: 868500 : asignnrStar.resize(nmat);
168 : 868500 : std::array< tk::real, 3 > TnlStar{{0, 0, 0}}, TnrStar{{0, 0, 0}};
169 : : std::vector< std::array< tk::real, 3 > > aTnlStar, aTnrStar;
170 [ + + ]: 2605500 : for (std::size_t k=0; k<nmat; ++k) {
171 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i)
172 [ + + ]: 20844000 : for (std::size_t j=0; j<3; ++j)
173 : : {
174 : 15633000 : asignnlStar[k][i][j] = asignnl[k][i][j];
175 : 15633000 : asignnrStar[k][i][j] = asignnr[k][i][j];
176 : : }
177 : :
178 [ + + ]: 3474000 : for (std::size_t i=0; i<1; ++i)
179 : : {
180 : 1737000 : asignnlStar[k][i][i] -=
181 : 1737000 : u[0][densityIdx(nmat,k)]*(Sl-vnl[0])*(Sm-vnl[0]);
182 : 1737000 : asignnrStar[k][i][i] -=
183 : 1737000 : u[1][densityIdx(nmat,k)]*(Sr-vnr[0])*(Sm-vnr[0]);
184 : : }
185 : :
186 [ + + ]: 5211000 : for (std::size_t i=1; i<3; ++i)
187 : : {
188 : 3474000 : asignnlStar[k][i][0] =
189 : 3474000 : ( (Sm-Sl)*u[0][densityIdx(nmat,k)]*asignnr[k][i][0]
190 : 3474000 : - (Sm-Sr)*u[1][densityIdx(nmat,k)]*asignnl[k][i][0]
191 : 3474000 : + (Sm-Sl)*u[0][densityIdx(nmat,k)]
192 : 3474000 : * (Sm-Sr)*u[1][densityIdx(nmat,k)]
193 : 3474000 : * (vnl[i]-vnr[i]) ) /
194 : : ( (Sm-Sl)*u[0][densityIdx(nmat,k)]
195 : 3474000 : - (Sm-Sr)*u[1][densityIdx(nmat,k)]);
196 : 3474000 : asignnrStar[k][i][0] =
197 : 3474000 : ( (Sm-Sl)*u[0][densityIdx(nmat,k)]*asignnr[k][i][0]
198 : 3474000 : - (Sm-Sr)*u[1][densityIdx(nmat,k)]*asignnl[k][i][0]
199 : 3474000 : + (Sm-Sl)*u[0][densityIdx(nmat,k)]
200 : 3474000 : * (Sm-Sr)*u[1][densityIdx(nmat,k)]
201 : 3474000 : * (vnl[i]-vnr[i]) ) /
202 : : ( (Sm-Sl)*u[0][densityIdx(nmat,k)]
203 : 3474000 : - (Sm-Sr)*u[1][densityIdx(nmat,k)]);
204 : : }
205 : : // Symmetry
206 : 1737000 : asignnlStar[k][0][1] = asignnlStar[k][1][0];
207 : 1737000 : asignnlStar[k][0][2] = asignnlStar[k][2][0];
208 : 1737000 : asignnrStar[k][0][1] = asignnrStar[k][1][0];
209 : 1737000 : asignnrStar[k][0][2] = asignnrStar[k][2][0];
210 : :
211 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i)
212 [ + + ]: 20844000 : for (std::size_t j=0; j<3; ++j)
213 : : {
214 : 15633000 : signnlStar[i][j] += asignnlStar[k][i][j];
215 : 15633000 : signnrStar[i][j] += asignnrStar[k][i][j];
216 : : }
217 : :
218 [ + - ]: 1737000 : auto asiglStar = tk::unrotateTensor(asignnlStar[k], fn);
219 [ + - ]: 1737000 : auto asigrStar = tk::unrotateTensor(asignnrStar[k], fn);
220 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i)
221 [ + + ]: 20844000 : for (std::size_t j=0; j<3; ++j)
222 : : {
223 : 15633000 : siglStar[i][j] += asiglStar[i][j];
224 : 15633000 : sigrStar[i][j] += asigrStar[i][j];
225 : : }
226 : 1737000 : aTnlStar.push_back(tk::matvec(asiglStar, fn));
227 [ - - ]: 1737000 : aTnrStar.push_back(tk::matvec(asigrStar, fn));
228 [ + + ]: 6948000 : for (std::size_t i=0; i<3; ++i)
229 : : {
230 : 5211000 : TnlStar[i] += aTnlStar[k][i];
231 : 5211000 : TnrStar[i] += aTnrStar[k][i];
232 : : }
233 : : }
234 : :
235 : 868500 : auto w_l = (vnl[0]-Sl)/(Sm-Sl);
236 : 868500 : auto w_r = (vnr[0]-Sr)/(Sm-Sr);
237 : :
238 : : std::array< tk::real, 3 > vnlStar, vnrStar;
239 : : // u*_L
240 : 868500 : vnlStar[0] = Sm;
241 : 868500 : vnlStar[1] = vnl[1]
242 : 868500 : + (signnlStar[1][0] - signnl[1][0]) / (w_l*rhol*(vnl[0]-Sl));
243 : 868500 : vnlStar[2] = vnl[2]
244 : 868500 : + (signnlStar[2][0] - signnl[2][0]) / (w_l*rhol*(vnl[0]-Sl));
245 : : // u*_R
246 : 868500 : vnrStar[0] = Sm;
247 : 868500 : vnrStar[1] = vnr[1]
248 : 868500 : + (signnrStar[1][0] - signnr[1][0]) / (w_r*rhor*(vnr[0]-Sr));
249 : 868500 : vnrStar[2] = vnr[2]
250 : 868500 : + (signnrStar[2][0] - signnr[2][0]) / (w_r*rhor*(vnr[0]-Sr));
251 : :
252 : 868500 : auto vlStar = tk::unrotateVector(vnlStar, fn);
253 : 868500 : auto vrStar = tk::unrotateVector(vnrStar, fn);
254 : :
255 [ + - ]: 868500 : auto uStar = u;
256 : :
257 : : tk::real rholStar(0.0), rhorStar(0.0);
258 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gnlStar, gnrStar;
259 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > > glStar, grStar;
260 [ + + ]: 2605500 : for (std::size_t k=0; k<nmat; ++k) {
261 : : // Left
262 [ - + ]: 1737000 : if (solidx[k] > 0)
263 : : {
264 [ - - ]: 0 : gnlStar.push_back(gnl[k]);
265 [ - - ]: 0 : gnlStar[k][0][0] = w_l * gnl[k][0][0]
266 : 0 : + gnl[k][0][1]*(vnl[1]-vnlStar[1])/(Sm-Sl)
267 : 0 : + gnl[k][0][2]*(vnl[2]-vnlStar[2])/(Sm-Sl);
268 : 0 : gnlStar[k][1][0] = w_l * gnl[k][1][0]
269 : 0 : + gnl[k][1][1]*(vnl[1]-vnlStar[1])/(Sm-Sl)
270 : 0 : + gnl[k][1][2]*(vnl[2]-vnlStar[2])/(Sm-Sl);
271 : 0 : gnlStar[k][2][0] = w_l * gnl[k][2][0]
272 : 0 : + gnl[k][2][1]*(vnl[1]-vnlStar[1])/(Sm-Sl)
273 : 0 : + gnl[k][2][2]*(vnl[2]-vnlStar[2])/(Sm-Sl);
274 : : // rotate g back to original frame of reference
275 [ - - ]: 0 : glStar.push_back(tk::unrotateTensor(gnlStar[k], fn));
276 : : }
277 [ - + ]: 1737000 : uStar[0][volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)];
278 [ - + ]: 1737000 : uStar[0][densityIdx(nmat, k)] = w_l * u[0][densityIdx(nmat, k)];
279 : 1737000 : uStar[0][energyIdx(nmat, k)] = w_l * u[0][energyIdx(nmat, k)]
280 : 1737000 : + ( - asignnl[k][0][0]*vnl[0]
281 : 1737000 : - asignnl[k][1][0]*vnl[1]
282 : 1737000 : - asignnl[k][2][0]*vnl[2]
283 : 1737000 : + asignnlStar[k][0][0]*vnlStar[0]
284 : 1737000 : + asignnlStar[k][1][0]*vnlStar[1]
285 : 1737000 : + asignnlStar[k][2][0]*vnlStar[2]
286 : 1737000 : ) / (Sm-Sl);
287 : 1737000 : rholStar += uStar[0][densityIdx(nmat, k)];
288 : :
289 : : // Right
290 [ - + ]: 1737000 : if (solidx[k] > 0)
291 : : {
292 [ - - ]: 0 : gnrStar.push_back(gnr[k]);
293 [ - - ]: 0 : gnrStar[k][0][0] = w_r * gnr[k][0][0]
294 : 0 : + gnr[k][0][1]*(vnr[1]-vnrStar[1])/(Sm-Sr)
295 : 0 : + gnr[k][0][2]*(vnr[2]-vnrStar[2])/(Sm-Sr);
296 : 0 : gnrStar[k][1][0] = w_r * gnr[k][1][0]
297 : 0 : + gnr[k][1][1]*(vnr[1]-vnrStar[1])/(Sm-Sr)
298 : 0 : + gnr[k][1][2]*(vnr[2]-vnrStar[2])/(Sm-Sr);
299 : 0 : gnrStar[k][2][0] = w_r * gnr[k][2][0]
300 : 0 : + gnr[k][2][1]*(vnr[1]-vnrStar[1])/(Sm-Sr)
301 : 0 : + gnr[k][2][2]*(vnr[2]-vnrStar[2])/(Sm-Sr);
302 : : // rotate g back to original frame of reference
303 [ - - ]: 0 : grStar.push_back(tk::unrotateTensor(gnrStar[k], fn));
304 : : }
305 : 1737000 : uStar[1][volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)];
306 : 1737000 : uStar[1][densityIdx(nmat, k)] = w_r * u[1][densityIdx(nmat, k)];
307 : 1737000 : uStar[1][energyIdx(nmat, k)] = w_r * u[1][energyIdx(nmat, k)]
308 : 1737000 : + ( - asignnr[k][0][0]*vnr[0]
309 : 1737000 : - asignnr[k][1][0]*vnr[1]
310 : 1737000 : - asignnr[k][2][0]*vnr[2]
311 : 1737000 : + asignnrStar[k][0][0]*vnrStar[0]
312 : 1737000 : + asignnrStar[k][1][0]*vnrStar[1]
313 : 1737000 : + asignnrStar[k][2][0]*vnrStar[2]
314 : 1737000 : ) / (Sm-Sr);
315 : 1737000 : rhorStar += uStar[1][densityIdx(nmat, k)];
316 : : }
317 : :
318 : : // Numerical fluxes
319 : : // -------------------------------------------------------------------------
320 [ + + ]: 868500 : if (0.0 <= Sl) {
321 : :
322 [ + + ]: 375248 : for (std::size_t idir=0; idir<3; ++idir)
323 : 281436 : flx[momentumIdx(nmat, idir)] =
324 : 281436 : u[0][momentumIdx(nmat, idir)] * vnl[0] - Tnl[idir];
325 : :
326 [ + + ]: 281436 : for (std::size_t k=0; k<nmat; ++k) {
327 [ - + ]: 187624 : flx[volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)] * vnl[0];
328 [ - + ]: 187624 : flx[densityIdx(nmat, k)] = u[0][densityIdx(nmat, k)] * vnl[0];
329 : 187624 : flx[energyIdx(nmat, k)] = u[0][energyIdx(nmat, k)] * vnl[0]
330 : 187624 : - ul * aTnl[k][0] - vl * aTnl[k][1] - wl * aTnl[k][2];
331 [ - + ]: 187624 : if (solidx[k] > 0) {
332 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
333 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
334 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
335 : 0 : gl[k][i][0] * ul +
336 : 0 : gl[k][i][1] * vl +
337 : 0 : gl[k][i][2] * wl ) * fn[j];
338 : : }
339 : : }
340 : :
341 : : // Quantities for non-conservative terms
342 : : // Store Riemann-advected partial pressures
343 [ + + ]: 281436 : for (std::size_t k=0; k<nmat; ++k)
344 [ + - ]: 187624 : flx.push_back(std::sqrt((aTnl[k][0]*aTnl[k][0]
345 : 187624 : +aTnl[k][1]*aTnl[k][1]
346 [ + - ]: 187624 : +aTnl[k][2]*aTnl[k][2])));
347 : : // Store Riemann velocity
348 [ + - ]: 93812 : flx.push_back(vnl[0]);
349 [ + + ]: 281436 : for (std::size_t k=0; k<nmat; ++k) {
350 [ - + ]: 187624 : if (solidx[k] > 0) {
351 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
352 [ - - ]: 0 : flx.push_back(aTnl[k][i]);
353 : : }
354 : : }
355 : :
356 : : }
357 : :
358 [ + - ][ + + ]: 774688 : else if (Sl < 0.0 && 0.0 <= Sm) {
359 : :
360 [ + + ]: 1685780 : for (std::size_t idir=0; idir<3; ++idir)
361 : 1264335 : flx[momentumIdx(nmat, idir)] =
362 : 1264335 : vlStar[idir] * rholStar * Sm - TnlStar[idir];
363 : :
364 [ + + ]: 1264335 : for (std::size_t k=0; k<nmat; ++k) {
365 [ - + ]: 842890 : flx[volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)] * Sm;
366 [ - + ]: 842890 : flx[densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)] * Sm;
367 : 842890 : flx[energyIdx(nmat, k)] = uStar[0][energyIdx(nmat, k)] * Sm
368 : 842890 : - vlStar[0] * aTnlStar[k][0]
369 : 842890 : - vlStar[1] * aTnlStar[k][1]
370 : 842890 : - vlStar[2] * aTnlStar[k][2];
371 [ - + ]: 842890 : if (solidx[k] > 0) {
372 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
373 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
374 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
375 : 0 : glStar[k][i][0] * vlStar[0] +
376 : 0 : glStar[k][i][1] * vlStar[1] +
377 : 0 : glStar[k][i][2] * vlStar[2] ) * fn[j];
378 : : }
379 : : }
380 : :
381 : : // Quantities for non-conservative terms
382 : : // Store Riemann-advected partial pressures
383 [ + + ]: 1264335 : for (std::size_t k=0; k<nmat; ++k)
384 [ + - ]: 842890 : flx.push_back(std::sqrt(aTnlStar[k][0]*aTnlStar[k][0]
385 : 842890 : +aTnlStar[k][1]*aTnlStar[k][1]
386 [ + - ]: 842890 : +aTnlStar[k][2]*aTnlStar[k][2]));
387 : : // Store Riemann velocity
388 [ + - ]: 421445 : flx.push_back(Sm);
389 [ + + ]: 1264335 : for (std::size_t k=0; k<nmat; ++k) {
390 [ - + ]: 842890 : if (solidx[k] > 0) {
391 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
392 [ - - ]: 0 : flx.push_back(aTnlStar[k][i]);
393 : : }
394 : : }
395 : :
396 : : }
397 : :
398 [ - + ][ + + ]: 353243 : else if (Sm < 0.0 && 0.0 <= Sr) {
399 : :
400 [ + + ]: 883044 : for (std::size_t idir=0; idir<3; ++idir)
401 : 662283 : flx[momentumIdx(nmat, idir)] =
402 : 662283 : vrStar[idir] * rhorStar * Sm - TnrStar[idir];
403 : :
404 [ + + ]: 662283 : for (std::size_t k=0; k<nmat; ++k) {
405 [ - + ]: 441522 : flx[volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)] * Sm;
406 [ - + ]: 441522 : flx[densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)] * Sm;
407 : 441522 : flx[energyIdx(nmat, k)] = uStar[1][energyIdx(nmat, k)] * Sm
408 : 441522 : - vrStar[0] * aTnrStar[k][0]
409 : 441522 : - vrStar[1] * aTnrStar[k][1]
410 : 441522 : - vrStar[2] * aTnrStar[k][2];
411 [ - + ]: 441522 : if (solidx[k] > 0) {
412 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
413 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
414 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
415 : 0 : grStar[k][i][0] * vrStar[0] +
416 : 0 : grStar[k][i][1] * vrStar[1] +
417 : 0 : grStar[k][i][2] * vrStar[2] ) * fn[j];
418 : : }
419 : : }
420 : :
421 : : // Quantities for non-conservative terms
422 : : // Store Riemann-advected partial pressures
423 [ + + ]: 662283 : for (std::size_t k=0; k<nmat; ++k)
424 [ + - ]: 441522 : flx.push_back(std::sqrt(aTnrStar[k][0]*aTnrStar[k][0]
425 : 441522 : +aTnrStar[k][1]*aTnrStar[k][1]
426 [ + - ]: 441522 : +aTnrStar[k][2]*aTnrStar[k][2]));
427 : : // Store Riemann velocity
428 [ + - ]: 220761 : flx.push_back(Sm);
429 [ + + ]: 662283 : for (std::size_t k=0; k<nmat; ++k) {
430 [ - + ]: 441522 : if (solidx[k] > 0) {
431 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
432 [ - - ]: 0 : flx.push_back(aTnrStar[k][i]);
433 : : }
434 : : }
435 : :
436 : : }
437 : :
438 : : else {
439 : :
440 [ + + ]: 529928 : for (std::size_t idir=0; idir<3; ++idir)
441 : 397446 : flx[momentumIdx(nmat, idir)] =
442 : 397446 : u[1][momentumIdx(nmat, idir)] * vnr[0] - Tnr[idir];
443 : :
444 [ + + ]: 397446 : for (std::size_t k=0; k<nmat; ++k) {
445 [ - + ]: 264964 : flx[volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)] * vnr[0];
446 [ - + ]: 264964 : flx[densityIdx(nmat, k)] = u[1][densityIdx(nmat, k)] * vnr[0];
447 : 264964 : flx[energyIdx(nmat, k)] = u[1][energyIdx(nmat, k)] * vnr[0]
448 : 264964 : - ur * aTnr[k][0] - vr * aTnr[k][1] - wr * aTnr[k][2];
449 [ - + ]: 264964 : if (solidx[k] > 0) {
450 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
451 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
452 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
453 : 0 : gr[k][i][0] * ur +
454 : 0 : gr[k][i][1] * vr +
455 : 0 : gr[k][i][2] * wr ) * fn[j];
456 : : }
457 : : }
458 : :
459 : : // Quantities for non-conservative terms
460 : : // Store Riemann-advected partial pressures
461 [ + + ]: 397446 : for (std::size_t k=0; k<nmat; ++k)
462 [ + - ][ - - ]: 264964 : flx.push_back(std::sqrt(aTnr[k][0]*aTnr[k][0]
463 : 264964 : +aTnr[k][1]*aTnr[k][1]
464 [ + - ]: 264964 : +aTnr[k][2]*aTnr[k][2]));
465 : : // Store Riemann velocity
466 [ + - ]: 132482 : flx.push_back(vnr[0]);
467 [ + + ]: 397446 : for (std::size_t k=0; k<nmat; ++k) {
468 [ - + ]: 264964 : if (solidx[k] > 0) {
469 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
470 [ - - ]: 0 : flx.push_back(aTnr[k][i]);
471 : : }
472 : : }
473 : :
474 : : }
475 : :
476 : : Assert( flx.size() == (ncomp+nmat+1+3*nsld), "Size of "
477 : : "multi-material flux vector incorrect" );
478 : :
479 : 868500 : return flx;
480 : : }
481 : :
482 : : ////! Flux type accessor
483 : : ////! \return Flux type
484 : : //static ctr::FluxType type() noexcept {
485 : : // return ctr::FluxType::HLLCMultiMat; }
486 : : };
487 : :
488 : : } // inciter::
489 : :
490 : : #endif // HLLCMultiMat_h
|
