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 : 868500 : 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 : 868500 : 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 : 868500 : tk::real pl(0.0), pr(0.0);
69 : 868500 : tk::real acl(0.0), acr(0.0);
70 [ + - ][ + - ]: 1737000 : 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 : 1737000 : 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 : 1737000 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gl, gr;
79 : 1737000 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gnl, gnr;
80 : 1737000 : 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 : 1737000 : 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 [ + - ][ + - ]: 1737000 : 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 [ + - ]: 3474000 : 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 : 1737000 : 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 [ + - ][ + - ]: 1737000 : 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 [ + - ]: 3474000 : 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 : 1737000 : 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 : 1737000 : 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 : 3474000 : ( (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 : 3474000 : ( (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 [ + - ]: 1737000 : auto uStar = u;
256 : :
257 : 868500 : tk::real rholStar(0.0), rhorStar(0.0);
258 : 1737000 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > gnlStar, gnrStar;
259 : 1737000 : 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 : 3474000 : 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 : 3474000 : 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 [ + + ]: 373748 : for (std::size_t idir=0; idir<3; ++idir)
323 : 280311 : flx[momentumIdx(nmat, idir)] =
324 : 280311 : u[0][momentumIdx(nmat, idir)] * vnl[0] - Tnl[idir];
325 : :
326 [ + + ]: 280311 : for (std::size_t k=0; k<nmat; ++k) {
327 : 186874 : flx[volfracIdx(nmat, k)] = u[0][volfracIdx(nmat, k)] * vnl[0];
328 : 186874 : flx[densityIdx(nmat, k)] = u[0][densityIdx(nmat, k)] * vnl[0];
329 : 373748 : flx[energyIdx(nmat, k)] = u[0][energyIdx(nmat, k)] * vnl[0]
330 : 186874 : - ul * aTnl[k][0] - vl * aTnl[k][1] - wl * aTnl[k][2];
331 [ - + ]: 186874 : 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 [ + + ]: 280311 : for (std::size_t k=0; k<nmat; ++k)
344 : 373748 : flx.push_back(std::sqrt((aTnl[k][0]*aTnl[k][0]
345 : 186874 : +aTnl[k][1]*aTnl[k][1]
346 [ + - ]: 186874 : +aTnl[k][2]*aTnl[k][2])));
347 : : // Store Riemann velocity
348 [ + - ]: 93437 : flx.push_back(vnl[0]);
349 [ + + ]: 280311 : for (std::size_t k=0; k<nmat; ++k) {
350 [ - + ]: 186874 : 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 [ + + ]: 280311 : for (std::size_t k=0; k<nmat; ++k) {
356 [ - + ]: 186874 : if (solidx[k] > 0) {
357 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
358 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
359 [ - - ]: 0 : flx.push_back(gl[k][i][j]);
360 : : }
361 : : }
362 : :
363 : : }
364 : :
365 [ + - ][ + + ]: 775063 : else if (Sl < 0.0 && 0.0 <= Sm) {
366 : :
367 [ + + ]: 1687280 : for (std::size_t idir=0; idir<3; ++idir)
368 : 1265460 : flx[momentumIdx(nmat, idir)] =
369 : 1265460 : vlStar[idir] * rholStar * Sm - TnlStar[idir];
370 : :
371 [ + + ]: 1265460 : for (std::size_t k=0; k<nmat; ++k) {
372 : 843640 : flx[volfracIdx(nmat, k)] = uStar[0][volfracIdx(nmat, k)] * Sm;
373 : 843640 : flx[densityIdx(nmat, k)] = uStar[0][densityIdx(nmat, k)] * Sm;
374 : 1687280 : flx[energyIdx(nmat, k)] = uStar[0][energyIdx(nmat, k)] * Sm
375 : 843640 : - vlStar[0] * aTnlStar[k][0]
376 : 843640 : - vlStar[1] * aTnlStar[k][1]
377 : 843640 : - vlStar[2] * aTnlStar[k][2];
378 [ - + ]: 843640 : if (solidx[k] > 0) {
379 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
380 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
381 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
382 : 0 : glStar[k][i][0] * vlStar[0] +
383 : 0 : glStar[k][i][1] * vlStar[1] +
384 : 0 : glStar[k][i][2] * vlStar[2] ) * fn[j];
385 : : }
386 : : }
387 : :
388 : : // Quantities for non-conservative terms
389 : : // Store Riemann-advected partial pressures
390 [ + + ]: 1265460 : for (std::size_t k=0; k<nmat; ++k)
391 : 1687280 : flx.push_back(std::sqrt(aTnlStar[k][0]*aTnlStar[k][0]
392 : 843640 : +aTnlStar[k][1]*aTnlStar[k][1]
393 [ + - ]: 843640 : +aTnlStar[k][2]*aTnlStar[k][2]));
394 : : // Store Riemann velocity
395 [ + - ]: 421820 : flx.push_back(Sm);
396 [ + + ]: 1265460 : for (std::size_t k=0; k<nmat; ++k) {
397 [ - + ]: 843640 : if (solidx[k] > 0) {
398 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
399 [ - - ]: 0 : flx.push_back(aTnlStar[k][i]);
400 : : }
401 : : }
402 [ + + ]: 1265460 : for (std::size_t k=0; k<nmat; ++k) {
403 [ - + ]: 843640 : if (solidx[k] > 0) {
404 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
405 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
406 [ - - ]: 0 : flx.push_back(glStar[k][i][j]);
407 : : }
408 : 421820 : }
409 : :
410 : : }
411 : :
412 [ + - ][ + + ]: 353243 : else if (Sm < 0.0 && 0.0 <= Sr) {
413 : :
414 [ + + ]: 881544 : for (std::size_t idir=0; idir<3; ++idir)
415 : 661158 : flx[momentumIdx(nmat, idir)] =
416 : 661158 : vrStar[idir] * rhorStar * Sm - TnrStar[idir];
417 : :
418 [ + + ]: 661158 : for (std::size_t k=0; k<nmat; ++k) {
419 : 440772 : flx[volfracIdx(nmat, k)] = uStar[1][volfracIdx(nmat, k)] * Sm;
420 : 440772 : flx[densityIdx(nmat, k)] = uStar[1][densityIdx(nmat, k)] * Sm;
421 : 881544 : flx[energyIdx(nmat, k)] = uStar[1][energyIdx(nmat, k)] * Sm
422 : 440772 : - vrStar[0] * aTnrStar[k][0]
423 : 440772 : - vrStar[1] * aTnrStar[k][1]
424 : 440772 : - vrStar[2] * aTnrStar[k][2];
425 [ - + ]: 440772 : if (solidx[k] > 0) {
426 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
427 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
428 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
429 : 0 : grStar[k][i][0] * vrStar[0] +
430 : 0 : grStar[k][i][1] * vrStar[1] +
431 : 0 : grStar[k][i][2] * vrStar[2] ) * fn[j];
432 : : }
433 : : }
434 : :
435 : : // Quantities for non-conservative terms
436 : : // Store Riemann-advected partial pressures
437 [ + + ]: 661158 : for (std::size_t k=0; k<nmat; ++k)
438 : 881544 : flx.push_back(std::sqrt(aTnrStar[k][0]*aTnrStar[k][0]
439 : 440772 : +aTnrStar[k][1]*aTnrStar[k][1]
440 [ + - ]: 440772 : +aTnrStar[k][2]*aTnrStar[k][2]));
441 : : // Store Riemann velocity
442 [ + - ]: 220386 : flx.push_back(Sm);
443 [ + + ]: 661158 : for (std::size_t k=0; k<nmat; ++k) {
444 [ - + ]: 440772 : if (solidx[k] > 0) {
445 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
446 [ - - ]: 0 : flx.push_back(aTnrStar[k][i]);
447 : : }
448 : : }
449 [ + + ]: 661158 : for (std::size_t k=0; k<nmat; ++k) {
450 [ - + ]: 440772 : if (solidx[k] > 0) {
451 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
452 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
453 [ - - ]: 0 : flx.push_back(grStar[k][i][j]);
454 : : }
455 : 220386 : }
456 : :
457 : : }
458 : :
459 : : else {
460 : :
461 [ + + ]: 531428 : for (std::size_t idir=0; idir<3; ++idir)
462 : 398571 : flx[momentumIdx(nmat, idir)] =
463 : 398571 : u[1][momentumIdx(nmat, idir)] * vnr[0] - Tnr[idir];
464 : :
465 [ + + ]: 398571 : for (std::size_t k=0; k<nmat; ++k) {
466 : 265714 : flx[volfracIdx(nmat, k)] = u[1][volfracIdx(nmat, k)] * vnr[0];
467 : 265714 : flx[densityIdx(nmat, k)] = u[1][densityIdx(nmat, k)] * vnr[0];
468 : 531428 : flx[energyIdx(nmat, k)] = u[1][energyIdx(nmat, k)] * vnr[0]
469 : 265714 : - ur * aTnr[k][0] - vr * aTnr[k][1] - wr * aTnr[k][2];
470 [ - + ]: 265714 : if (solidx[k] > 0) {
471 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
472 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
473 : 0 : flx[deformIdx(nmat,solidx[k],i,j)] = (
474 : 0 : gr[k][i][0] * ur +
475 : 0 : gr[k][i][1] * vr +
476 : 0 : gr[k][i][2] * wr ) * fn[j];
477 : : }
478 : : }
479 : :
480 : : // Quantities for non-conservative terms
481 : : // Store Riemann-advected partial pressures
482 [ + + ]: 398571 : for (std::size_t k=0; k<nmat; ++k)
483 : 531428 : flx.push_back(std::sqrt(aTnr[k][0]*aTnr[k][0]
484 : 265714 : +aTnr[k][1]*aTnr[k][1]
485 [ + - ]: 265714 : +aTnr[k][2]*aTnr[k][2]));
486 : : // Store Riemann velocity
487 [ + - ]: 132857 : flx.push_back(vnr[0]);
488 [ + + ]: 398571 : for (std::size_t k=0; k<nmat; ++k) {
489 [ - + ]: 265714 : if (solidx[k] > 0) {
490 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
491 [ - - ]: 0 : flx.push_back(aTnr[k][i]);
492 : : }
493 : : }
494 [ + + ]: 398571 : for (std::size_t k=0; k<nmat; ++k) {
495 [ - + ]: 265714 : if (solidx[k] > 0) {
496 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
497 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
498 [ - - ]: 0 : flx.push_back(gr[k][i][j]);
499 : : }
500 : : }
501 : :
502 : : }
503 : :
504 [ - + ][ - - ]: 868500 : Assert( flx.size() == (ncomp+nmat+1+3*nsld+9*nsld), "Size of "
[ - - ][ - - ]
505 : : "multi-material flux vector incorrect" );
506 : :
507 : 1737000 : return flx;
508 : : }
509 : :
510 : : ////! Flux type accessor
511 : : ////! \return Flux type
512 : : //static ctr::FluxType type() noexcept {
513 : : // return ctr::FluxType::HLLCMultiMat; }
514 : : };
515 : :
516 : : } // inciter::
517 : :
518 : : #endif // HLLCMultiMat_h
|
