Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/MultiMatTerms.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Functions for computing volume integrals of multi-material terms
9 : : using DG methods
10 : : \details This file contains functionality for computing volume integrals of
11 : : non-conservative and pressure relaxation terms that appear in the
12 : : multi-material hydrodynamic equations, using the discontinuous Galerkin
13 : : method for various orders of numerical representation.
14 : : */
15 : : // *****************************************************************************
16 : :
17 : : #include "QuinoaConfig.hpp"
18 : :
19 : : #include "MultiMatTerms.hpp"
20 : : #include "Vector.hpp"
21 : : #include "Quadrature.hpp"
22 : : #include "MultiMat/MultiMatIndexing.hpp"
23 : : #include "Reconstruction.hpp"
24 : : #include "Inciter/InputDeck/InputDeck.hpp"
25 : : #include "EoS/GetMatProp.hpp"
26 : :
27 : : namespace inciter {
28 : : extern ctr::InputDeck g_inputdeck;
29 : : }
30 : :
31 : : // Lapacke forward declarations
32 : : extern "C" {
33 : :
34 : : using lapack_int = long;
35 : :
36 : : #define LAPACK_ROW_MAJOR 101
37 : :
38 : : lapack_int LAPACKE_dsysv( int, char, lapack_int, lapack_int, double*,
39 : : lapack_int, lapack_int*, double*, lapack_int );
40 : :
41 : : }
42 : :
43 : : namespace tk {
44 : :
45 : : void
46 : 5430 : nonConservativeInt( std::size_t nmat,
47 : : const std::vector< inciter::EOS >& mat_blk,
48 : : const std::size_t ndof,
49 : : const std::size_t rdof,
50 : : const std::size_t nelem,
51 : : const std::vector< std::size_t >& inpoel,
52 : : const UnsMesh::Coords& coord,
53 : : const Fields& geoElem,
54 : : const Fields& U,
55 : : const Fields& P,
56 : : const std::vector< std::vector< tk::real > >& riemannDeriv,
57 : : const std::vector< std::size_t >& ndofel,
58 : : Fields& R,
59 : : int intsharp )
60 : : // *****************************************************************************
61 : : // Compute volume integrals for multi-material DG
62 : : //! \details This is called for multi-material DG, computing volume integrals of
63 : : //! terms in the volume fraction and energy equations, which do not exist in
64 : : //! the single-material flow formulation (for `CompFlow` DG). For further
65 : : //! details see Pelanti, M., & Shyue, K. M. (2019). A numerical model for
66 : : //! multiphase liquid–vapor–gas flows with interfaces and cavitation.
67 : : //! International Journal of Multiphase Flow, 113, 208-230.
68 : : //! \param[in] nmat Number of materials in this PDE system
69 : : //! \param[in] mat_blk EOS material block
70 : : //! \param[in] ndof Maximum number of degrees of freedom
71 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
72 : : //! \param[in] nelem Total number of elements
73 : : //! \param[in] inpoel Element-node connectivity
74 : : //! \param[in] coord Array of nodal coordinates
75 : : //! \param[in] geoElem Element geometry array
76 : : //! \param[in] U Solution vector at recent time step
77 : : //! \param[in] P Vector of primitive quantities at recent time step
78 : : //! \param[in] riemannDeriv Derivatives of partial-pressures and velocities
79 : : //! computed from the Riemann solver for use in the non-conservative terms
80 : : //! \param[in] ndofel Vector of local number of degrees of freedome
81 : : //! \param[in,out] R Right-hand side vector added to
82 : : //! \param[in] intsharp Interface reconstruction indicator
83 : : // *****************************************************************************
84 : : {
85 : : using inciter::volfracIdx;
86 : : using inciter::densityIdx;
87 : : using inciter::momentumIdx;
88 : : using inciter::energyIdx;
89 : : using inciter::velocityIdx;
90 : : using inciter::deformIdx;
91 : : using inciter::newSolidsAccFn;
92 : :
93 : : const auto& solidx =
94 : 5430 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
95 : :
96 : 5430 : const auto& cx = coord[0];
97 : 5430 : const auto& cy = coord[1];
98 : 5430 : const auto& cz = coord[2];
99 : :
100 : 5430 : auto ncomp = U.nprop()/rdof;
101 : 5430 : auto nprim = P.nprop()/rdof;
102 : :
103 : : // compute volume integrals
104 [ + + ]: 2353170 : for (std::size_t e=0; e<nelem; ++e)
105 : : {
106 [ + - ]: 2347740 : auto ng = tk::NGvol(ndofel[e]);
107 : :
108 : : // arrays for quadrature points
109 : 4695480 : std::array< std::vector< real >, 3 > coordgp;
110 : 4695480 : std::vector< real > wgp;
111 : :
112 [ + - ]: 2347740 : coordgp[0].resize( ng );
113 [ + - ]: 2347740 : coordgp[1].resize( ng );
114 [ + - ]: 2347740 : coordgp[2].resize( ng );
115 [ + - ]: 2347740 : wgp.resize( ng );
116 : :
117 [ + - ]: 2347740 : GaussQuadratureTet( ng, coordgp, wgp );
118 : :
119 : : // Extract the element coordinates
120 : : std::array< std::array< real, 3>, 4 > coordel {{
121 : 7043220 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
122 : 7043220 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
123 : 7043220 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
124 : 7043220 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
125 : 25825140 : }};
126 : :
127 : : auto jacInv =
128 : 2347740 : inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
129 : :
130 : : // Compute the derivatives of basis function for DG(P1)
131 : 4695480 : std::array< std::vector<tk::real>, 3 > dBdx;
132 [ + + ]: 2347740 : if (ndofel[e] > 1)
133 [ + - ]: 500280 : dBdx = eval_dBdx_p1( ndofel[e], jacInv );
134 : :
135 : : // Gaussian quadrature
136 [ + + ]: 6696600 : for (std::size_t igp=0; igp<ng; ++igp)
137 : : {
138 [ - + ]: 4348860 : if (ndofel[e] > 4)
139 [ - - ]: 0 : eval_dBdx_p2( igp, coordgp, jacInv, dBdx );
140 : :
141 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
142 : : // basis functions and solutions by default. This implies that P0P1 is
143 : : // unsupported in the p-adaptive DG (PDG).
144 : : std::size_t dof_el;
145 [ + + ]: 4348860 : if (rdof > ndof)
146 : : {
147 : 527160 : dof_el = rdof;
148 : : }
149 : : else
150 : : {
151 : 3821700 : dof_el = ndofel[e];
152 : : }
153 : :
154 : : // Compute the basis function
155 : : auto B =
156 [ + - ]: 8697720 : eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
157 : :
158 [ + - ]: 4348860 : auto wt = wgp[igp] * geoElem(e, 0);
159 : :
160 : : auto state = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim,
161 : : rdof, nmat, e, dof_el, inpoel, coord, geoElem,
162 [ + - ]: 8697720 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, U, P);
163 : :
164 : : // get bulk properties
165 : 4348860 : tk::real rhob(0.0);
166 [ + + ]: 14139480 : for (std::size_t k=0; k<nmat; ++k)
167 : 9790620 : rhob += state[densityIdx(nmat, k)];
168 : :
169 : : // get the velocity vector
170 : 4348860 : std::array< tk::real, 3 > vel{{ state[ncomp+velocityIdx(nmat, 0)],
171 : 4348860 : state[ncomp+velocityIdx(nmat, 1)],
172 : 8697720 : state[ncomp+velocityIdx(nmat, 2)] }};
173 : :
174 [ + - ]: 8697720 : std::vector< tk::real > ymat(nmat, 0.0);
175 : 4348860 : std::array< tk::real, 3 > dap{{0.0, 0.0, 0.0}};
176 [ + + ]: 14139480 : for (std::size_t k=0; k<nmat; ++k)
177 : : {
178 : 9790620 : ymat[k] = state[densityIdx(nmat, k)]/rhob;
179 : :
180 : 9790620 : std::size_t mark(3*k);
181 [ + + ]: 9790620 : if (solidx[k] > 0) mark = 3*nmat+ndof+3*(solidx[k]-1);
182 : :
183 [ + + ]: 39162480 : for (std::size_t idir=0; idir<3; ++idir)
184 : 29371860 : dap[idir] += riemannDeriv[mark+idir][e];
185 : : }
186 : :
187 : : // compute non-conservative terms
188 : : std::vector< std::vector< tk::real > > ncf
189 [ + - ][ + - ]: 13046580 : (ncomp, std::vector<tk::real>(ndof,0.0));
190 : :
191 [ + + ]: 17395440 : for (std::size_t idir=0; idir<3; ++idir)
192 [ + + ]: 48605760 : for(std::size_t idof=0; idof<ndof; ++idof)
193 : 35559180 : ncf[momentumIdx(nmat, idir)][idof] = 0.0;
194 : :
195 [ + + ]: 14139480 : for (std::size_t k=0; k<nmat; ++k)
196 : : {
197 : : // evaluate non-conservative term for energy equation
198 : 9790620 : std::size_t mark(3*k);
199 [ + + ]: 9790620 : if (solidx[k] > 0) mark = 3*nmat+ndof+3*(solidx[k]-1);
200 : :
201 [ + + ]: 34589640 : for(std::size_t idof=0; idof<ndof; ++idof)
202 : : {
203 : 24799020 : ncf[densityIdx(nmat, k)][idof] = 0.0;
204 : :
205 [ + + ]: 99196080 : for (std::size_t idir=0; idir<3; ++idir)
206 : 148794120 : ncf[energyIdx(nmat, k)][idof] -= vel[idir] * ( ymat[k]*dap[idir]
207 : 74397060 : - riemannDeriv[mark+idir][e] );
208 : : }
209 : :
210 : : // Evaluate non-conservative term for volume fraction equation:
211 : : // Here we make an assumption that the derivative of Riemann velocity
212 : : // times the basis function is constant. Therefore, when P0P1/DGP1/DGP2
213 : : // are used for constant velocity problems, the discretization is
214 : : // consistent. However, for a general problem with varying velocity,
215 : : // there will be errors since the said derivative is not constant.
216 : : // A discretization that solves this issue has not been implemented yet.
217 : : // Nevertheless, this does not affect high-order accuracy in
218 : : // single material regions for problems with sharp interfaces. Since
219 : : // volume fractions are nearly constant in such regions, using
220 : : // high-order for volume fractions does not show any benefits over
221 : : // THINC. Therefore, for such problems, we only use FV for the volume
222 : : // fractions, and these non-conservative high-order terms do not need
223 : : // to be computed.
224 : : // In summary, high-order discretization for non-conservative terms in
225 : : // volume fraction equations is avoided for sharp interface problems.
226 [ - + ][ - - ]: 9790620 : if (ndof <= 4 || intsharp == 1) {
227 [ + + ]: 34589640 : for(std::size_t idof=0; idof<ndof; ++idof)
228 : 49598040 : ncf[volfracIdx(nmat, k)][idof] = state[volfracIdx(nmat, k)]
229 : 34589640 : * riemannDeriv[3*nmat+idof][e];
230 [ - - ]: 0 : } else if (intsharp == 0) { // If DGP2 without THINC
231 : : // DGP2 is discretized differently than DGP1/FV to guarantee 3rd order
232 : : // convergence for the testcases with uniform and constant velocity.
233 : :
234 : : // P0 contributions for all equations
235 [ - - ]: 0 : for(std::size_t idof=0; idof<ndof; ++idof)
236 : 0 : ncf[volfracIdx(nmat, k)][idof] = state[volfracIdx(nmat, k)]
237 : 0 : * riemannDeriv[3*nmat][e] * B[idof];
238 : : // High order contributions
239 [ - - ]: 0 : for(std::size_t idof=1; idof<ndof; ++idof)
240 [ - - ]: 0 : for(std::size_t idir=0; idir<3; ++idir)
241 : 0 : ncf[volfracIdx(nmat, k)][idof] += state[volfracIdx(nmat, k)]
242 : 0 : * vel[idir] * dBdx[idir][idof];
243 : : }
244 : : }
245 : :
246 [ + - ]: 4348860 : updateRhsNonCons( ncomp, nmat, ndof, ndofel[e], wt, e, B, dBdx, ncf, R );
247 : : }
248 : : }
249 : 5430 : }
250 : :
251 : : void
252 : 4348860 : updateRhsNonCons(
253 : : ncomp_t ncomp,
254 : : const std::size_t nmat,
255 : : const std::size_t ndof,
256 : : [[maybe_unused]] const std::size_t ndof_el,
257 : : const tk::real wt,
258 : : const std::size_t e,
259 : : const std::vector<tk::real>& B,
260 : : [[maybe_unused]] const std::array< std::vector<tk::real>, 3 >& dBdx,
261 : : const std::vector< std::vector< tk::real > >& ncf,
262 : : Fields& R )
263 : : // *****************************************************************************
264 : : // Update the rhs by adding the non-conservative term integrals
265 : : //! \param[in] ncomp Number of scalar components in this PDE system
266 : : //! \param[in] nmat Number of materials
267 : : //! \param[in] ndof Maximum number of degrees of freedom
268 : : //! \param[in] ndof_el Number of degrees of freedom for local element
269 : : //! \param[in] wt Weight of gauss quadrature point
270 : : //! \param[in] e Element index
271 : : //! \param[in] B Basis function evaluated at local quadrature point
272 : : //! \param[in] dBdx Vector of basis function derivatives
273 : : //! \param[in] ncf Vector of non-conservative terms
274 : : //! \param[in,out] R Right-hand side vector computed
275 : : // *****************************************************************************
276 : : {
277 : : using inciter::volfracIdx;
278 : : using inciter::energyIdx;
279 : : using inciter::volfracDofIdx;
280 : : using inciter::energyDofIdx;
281 : :
282 : : //Assert( dBdx[0].size() == ndof_el,
283 : : // "Size mismatch for basis function derivatives" );
284 : : //Assert( dBdx[1].size() == ndof_el,
285 : : // "Size mismatch for basis function derivatives" );
286 : : //Assert( dBdx[2].size() == ndof_el,
287 : : // "Size mismatch for basis function derivatives" );
288 : : //Assert( ncf.size() == ncomp,
289 : : // "Size mismatch for non-conservative term" );
290 [ - + ][ - - ]: 4348860 : Assert( ncf.size() == ncomp, "Size mismatch for non-conservative term" );
[ - - ][ - - ]
291 : :
292 [ + + ]: 47660460 : for (ncomp_t c=0; c<ncomp; ++c)
293 : : {
294 : 43311600 : auto mark = c*ndof;
295 : 43311600 : R(e, mark) += wt * ncf[c][0];
296 : : }
297 : :
298 [ + + ]: 4348860 : if( ndof_el > 1)
299 : : {
300 : : // Update rhs with distributions from volume fraction and energy equations
301 [ + + ]: 7504200 : for (std::size_t k=0; k<nmat; ++k)
302 : : {
303 [ + + ]: 20011200 : for(std::size_t idof = 1; idof < ndof; idof++)
304 : : {
305 : 15008400 : R(e, volfracDofIdx(nmat,k,ndof,idof)) +=
306 : 15008400 : wt * ncf[volfracIdx(nmat,k)][idof];
307 : 15008400 : R(e, energyDofIdx(nmat,k,ndof,idof)) +=
308 : 15008400 : wt * ncf[energyIdx(nmat,k)][idof] * B[idof];
309 : : }
310 : : }
311 : : }
312 : 4348860 : }
313 : :
314 : : std::vector< tk::real >
315 : 975256 : nonConservativeIntFV(
316 : : std::size_t nmat,
317 : : const std::size_t rdof,
318 : : const std::size_t e,
319 : : const std::array< tk::real, 3 >& fn,
320 : : const Fields& U,
321 : : const Fields& P,
322 : : const std::vector< tk::real >& var_riemann )
323 : : // *****************************************************************************
324 : : // Compute integrals of non-conservative terms for multi-material FV
325 : : //! \details This is called for multi-material FV, computing integrals of
326 : : //! terms in the volume fraction and energy equations, which do not exist in
327 : : //! the single-material flow formulation (for `CompFlow`). For further
328 : : //! details see Pelanti, M., & Shyue, K. M. (2019). A numerical model for
329 : : //! multiphase liquid–vapor–gas flows with interfaces and cavitation.
330 : : //! International Journal of Multiphase Flow, 113, 208-230.
331 : : //! \param[in] nmat Number of materials in this PDE system
332 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
333 : : //! \param[in] e Element for which contribution is to be calculated
334 : : //! \param[in] fn Face normal
335 : : //! \param[in] U Solution vector at recent time step
336 : : //! \param[in] P Vector of primitive quantities at recent time step
337 : : //! \param[in] var_riemann Riemann-values of partial-pressures and velocities
338 : : //! computed from the Riemann solver for use in the non-conservative terms
339 : : // *****************************************************************************
340 : : {
341 : : using inciter::volfracIdx;
342 : : using inciter::densityIdx;
343 : : using inciter::momentumIdx;
344 : : using inciter::energyIdx;
345 : : using inciter::velocityIdx;
346 : : using inciter::volfracDofIdx;
347 : : using inciter::densityDofIdx;
348 : : using inciter::velocityDofIdx;
349 : :
350 : 975256 : auto ncomp = U.nprop()/rdof;
351 : :
352 : : // get bulk properties
353 : 975256 : tk::real rhob(0.0), p_face(0.0);
354 [ + + ]: 3162824 : for (std::size_t k=0; k<nmat; ++k)
355 : : {
356 [ + - ]: 2187568 : rhob += U(e, densityDofIdx(nmat,k,rdof,0));
357 : 2187568 : p_face += var_riemann[k];
358 : : }
359 : :
360 [ + - ]: 975256 : std::array< tk::real, 3 > vel{{ P(e, velocityDofIdx(nmat,0,rdof,0)),
361 [ + - ]: 975256 : P(e, velocityDofIdx(nmat,1,rdof,0)),
362 [ + - ]: 1950512 : P(e, velocityDofIdx(nmat,2,rdof,0)) }};
363 : :
364 : : // compute non-conservative terms
365 [ + - ]: 975256 : std::vector< tk::real > ncf(ncomp, 0.0);
366 [ + - ]: 1950512 : std::vector< tk::real > ymat(nmat, 0.0);
367 [ + + ]: 3162824 : for (std::size_t k=0; k<nmat; ++k)
368 : : {
369 [ + - ]: 2187568 : ymat[k] = U(e, densityDofIdx(nmat,k,rdof,0))/rhob;
370 : :
371 : : // evaluate non-conservative term for energy equation
372 [ + + ]: 8750272 : for (std::size_t idir=0; idir<3; ++idir)
373 : 13125408 : ncf[energyIdx(nmat, k)] -= vel[idir] * ( ymat[k]*p_face*fn[idir]
374 : 6562704 : - var_riemann[k]*fn[idir] );
375 : :
376 : : // evaluate non-conservative term for volume fraction equation
377 [ + - ]: 4375136 : ncf[volfracIdx(nmat, k)] = U(e, volfracDofIdx(nmat,k,rdof,0))
378 : 2187568 : * var_riemann[nmat];
379 : : }
380 : :
381 : 1950512 : return ncf;
382 : : }
383 : :
384 : : void
385 : 450 : pressureRelaxationInt( std::size_t nmat,
386 : : const std::vector< inciter::EOS >& mat_blk,
387 : : const std::size_t ndof,
388 : : const std::size_t rdof,
389 : : const std::size_t nelem,
390 : : const std::vector< std::size_t >& inpoel,
391 : : const UnsMesh::Coords& coord,
392 : : const Fields& geoElem,
393 : : const Fields& U,
394 : : const Fields& P,
395 : : const std::vector< std::size_t >& ndofel,
396 : : const tk::real ct,
397 : : Fields& R,
398 : : int intsharp )
399 : : // *****************************************************************************
400 : : // Compute volume integrals of pressure relaxation terms in multi-material DG
401 : : //! \details This is called for multi-material DG to compute volume integrals of
402 : : //! finite pressure relaxation terms in the volume fraction and energy
403 : : //! equations, which do not exist in the single-material flow formulation (for
404 : : //! `CompFlow` DG). For further details see Dobrev, V. A., Kolev, T. V.,
405 : : //! Rieben, R. N., & Tomov, V. Z. (2016). Multi‐material closure model for
406 : : //! high‐order finite element Lagrangian hydrodynamics. International Journal
407 : : //! for Numerical Methods in Fluids, 82(10), 689-706.
408 : : //! \param[in] nmat Number of materials in this PDE system
409 : : //! \param[in] mat_blk EOS material block
410 : : //! \param[in] ndof Maximum number of degrees of freedom
411 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
412 : : //! \param[in] nelem Total number of elements
413 : : //! \param[in] inpoel Element-node connectivity
414 : : //! \param[in] coord Array of nodal coordinates
415 : : //! \param[in] geoElem Element geometry array
416 : : //! \param[in] U Solution vector at recent time step
417 : : //! \param[in] P Vector of primitive quantities at recent time step
418 : : //! \param[in] ndofel Vector of local number of degrees of freedome
419 : : //! \param[in] ct Pressure relaxation time-scale for this system
420 : : //! \param[in,out] R Right-hand side vector added to
421 : : //! \param[in] intsharp Interface reconstruction indicator
422 : : // *****************************************************************************
423 : : {
424 : : using inciter::volfracIdx;
425 : : using inciter::densityIdx;
426 : : using inciter::momentumIdx;
427 : : using inciter::energyIdx;
428 : : using inciter::pressureIdx;
429 : : using inciter::velocityIdx;
430 : : using inciter::deformIdx;
431 : :
432 : : const auto& solidx =
433 : 450 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
434 : :
435 : 450 : auto ncomp = U.nprop()/rdof;
436 : 450 : auto nprim = P.nprop()/rdof;
437 : :
438 : : // compute volume integrals
439 [ + + ]: 341550 : for (std::size_t e=0; e<nelem; ++e)
440 : : {
441 [ + - ]: 341100 : auto dx = geoElem(e,4)/2.0;
442 [ + - ]: 341100 : auto ng = NGvol(ndofel[e]);
443 : :
444 : : // arrays for quadrature points
445 : 682200 : std::array< std::vector< real >, 3 > coordgp;
446 : 682200 : std::vector< real > wgp;
447 : :
448 [ + - ]: 341100 : coordgp[0].resize( ng );
449 [ + - ]: 341100 : coordgp[1].resize( ng );
450 [ + - ]: 341100 : coordgp[2].resize( ng );
451 [ + - ]: 341100 : wgp.resize( ng );
452 : :
453 [ + - ]: 341100 : GaussQuadratureTet( ng, coordgp, wgp );
454 : :
455 : : // Compute the derivatives of basis function for DG(P1)
456 : 682200 : std::array< std::vector<real>, 3 > dBdx;
457 : :
458 : : // Gaussian quadrature
459 [ + + ]: 1591800 : for (std::size_t igp=0; igp<ng; ++igp)
460 : : {
461 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
462 : : // basis functions and solutions by default. This implies that P0P1 is
463 : : // unsupported in the p-adaptive DG (PDG).
464 : : std::size_t dof_el;
465 [ + + ]: 1250700 : if (rdof > ndof)
466 : : {
467 : 113700 : dof_el = rdof;
468 : : }
469 : : else
470 : : {
471 : 1137000 : dof_el = ndofel[e];
472 : : }
473 : :
474 : : // Compute the basis function
475 : : auto B =
476 [ + - ]: 2501400 : eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
477 : :
478 [ + - ]: 1250700 : auto wt = wgp[igp] * geoElem(e, 0);
479 : :
480 : : auto state = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim,
481 : : rdof, nmat, e, dof_el, inpoel, coord, geoElem,
482 [ + - ]: 2501400 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, U, P);
483 : :
484 : : // get bulk properties
485 : 1250700 : real rhob(0.0);
486 [ + + ]: 3752100 : for (std::size_t k=0; k<nmat; ++k)
487 : 2501400 : rhob += state[densityIdx(nmat, k)];
488 : :
489 : : // get pressures and bulk modulii
490 : 1250700 : real pb(0.0), nume(0.0), deno(0.0), trelax(0.0);
491 [ + - ][ + - ]: 2501400 : std::vector< real > apmat(nmat, 0.0), kmat(nmat, 0.0);
492 [ + + ]: 3752100 : for (std::size_t k=0; k<nmat; ++k)
493 : : {
494 : 2501400 : real arhomat = state[densityIdx(nmat, k)];
495 : 2501400 : real alphamat = state[volfracIdx(nmat, k)];
496 : 2501400 : apmat[k] = state[ncomp+pressureIdx(nmat, k)];
497 : 2501400 : real amat = 0.0;
498 [ - + ]: 2501400 : if (solidx[k] > 0)
499 : : {
500 : : std::array< std::array< tk::real, 3 >, 3 > g;
501 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
502 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
503 : 0 : g[i][j] = state[deformIdx(nmat,solidx[k],i,j)];
504 [ - - ]: 0 : auto grot = tk::rotateTensor(g, {{1.0, 0.0, 0.0}});
505 [ - - ]: 0 : amat = mat_blk[k].compute< inciter::EOS::soundspeed >( arhomat,
506 : 0 : apmat[k], alphamat, k, grot);
507 [ - - ]: 0 : grot = tk::rotateTensor(g, {{0.0, 1.0, 0.0}});
508 [ - - ]: 0 : amat = std::max(amat, mat_blk[k].compute< inciter::EOS::soundspeed >(
509 : 0 : arhomat, apmat[k], alphamat, k, grot));
510 [ - - ]: 0 : grot = tk::rotateTensor(g, {{0.0, 0.0, 1.0}});
511 [ - - ]: 0 : amat = std::max(amat, mat_blk[k].compute< inciter::EOS::soundspeed >(
512 : 0 : arhomat, apmat[k], alphamat, k, grot));
513 : : }
514 : : else
515 : : {
516 [ + - ]: 5002800 : amat = mat_blk[k].compute< inciter::EOS::soundspeed >( arhomat,
517 : 2501400 : apmat[k], alphamat, k );
518 : : }
519 : 2501400 : kmat[k] = arhomat * amat * amat;
520 : 2501400 : pb += apmat[k];
521 : :
522 : : // relaxation parameters
523 : 2501400 : trelax = std::max(trelax, ct*dx/amat);
524 : 2501400 : nume += alphamat * apmat[k] / kmat[k];
525 : 2501400 : deno += alphamat * alphamat / kmat[k];
526 : : }
527 : 1250700 : auto p_relax = nume/deno;
528 : :
529 : : // compute pressure relaxation terms
530 [ + - ]: 2501400 : std::vector< real > s_prelax(ncomp, 0.0);
531 [ + + ]: 3752100 : for (std::size_t k=0; k<nmat; ++k)
532 : : {
533 : 2501400 : auto s_alpha = (apmat[k]-p_relax*state[volfracIdx(nmat, k)])
534 : 2501400 : * (state[volfracIdx(nmat, k)]/kmat[k]) / trelax;
535 : 2501400 : s_prelax[volfracIdx(nmat, k)] = s_alpha;
536 : 2501400 : s_prelax[energyIdx(nmat, k)] = - pb*s_alpha;
537 : : }
538 : :
539 [ + - ]: 1250700 : updateRhsPre( ncomp, ndof, dof_el, wt, e, B, s_prelax, R );
540 : : }
541 : : }
542 : 450 : }
543 : :
544 : : void
545 : 1250700 : updateRhsPre(
546 : : ncomp_t ncomp,
547 : : const std::size_t ndof,
548 : : [[maybe_unused]] const std::size_t ndof_el,
549 : : const tk::real wt,
550 : : const std::size_t e,
551 : : const std::vector< tk::real >& B,
552 : : std::vector< tk::real >& ncf,
553 : : Fields& R )
554 : : // *****************************************************************************
555 : : // Update the rhs by adding the pressure relaxation integrals
556 : : //! \param[in] ncomp Number of scalar components in this PDE system
557 : : //! \param[in] ndof Maximum number of degrees of freedom
558 : : //! \param[in] ndof_el Number of degrees of freedom for local element
559 : : //! \param[in] wt Weight of gauss quadrature point
560 : : //! \param[in] e Element index
561 : : //! \param[in] B Basis function evaluated at local quadrature point
562 : : //! \param[in] ncf Vector of non-conservative terms
563 : : //! \param[in,out] R Right-hand side vector computed
564 : : // *****************************************************************************
565 : : {
566 : : //Assert( dBdx[0].size() == ndof_el,
567 : : // "Size mismatch for basis function derivatives" );
568 : : //Assert( dBdx[1].size() == ndof_el,
569 : : // "Size mismatch for basis function derivatives" );
570 : : //Assert( dBdx[2].size() == ndof_el,
571 : : // "Size mismatch for basis function derivatives" );
572 : : //Assert( ncf.size() == ncomp,
573 : : // "Size mismatch for non-conservative term" );
574 [ - + ][ - - ]: 1250700 : Assert( ncf.size() == ncomp, "Size mismatch for non-conservative term" );
[ - - ][ - - ]
575 : :
576 [ + + ]: 12507000 : for (ncomp_t c=0; c<ncomp; ++c)
577 : : {
578 : 11256300 : auto mark = c*ndof;
579 [ + + ]: 53211600 : for(std::size_t idof = 0; idof < ndof; idof++)
580 : 41955300 : R(e, mark+idof) += wt * ncf[c] * B[idof];
581 : : }
582 : 1250700 : }
583 : :
584 : : void
585 : 1568 : pressureRelaxationIntFV(
586 : : std::size_t nmat,
587 : : const std::vector< inciter::EOS >& mat_blk,
588 : : const std::size_t rdof,
589 : : const std::size_t nelem,
590 : : const std::vector< std::size_t >& inpoel,
591 : : const UnsMesh::Coords& coord,
592 : : const Fields& geoElem,
593 : : const Fields& U,
594 : : const Fields& P,
595 : : const tk::real ct,
596 : : Fields& R )
597 : : // *****************************************************************************
598 : : // Compute volume integrals of pressure relaxation terms in multi-material FV
599 : : //! \details This is called for multi-material FV to compute volume integrals of
600 : : //! finite pressure relaxation terms in the volume fraction and energy
601 : : //! equations, which do not exist in the single-material flow formulation (for
602 : : //! `CompFlow`). For further details see Dobrev, V. A., Kolev, T. V.,
603 : : //! Rieben, R. N., & Tomov, V. Z. (2016). Multi‐material closure model for
604 : : //! high‐order finite element Lagrangian hydrodynamics. International Journal
605 : : //! for Numerical Methods in Fluids, 82(10), 689-706.
606 : : //! \param[in] nmat Number of materials in this PDE system
607 : : //! \param[in] mat_blk EOS material block
608 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
609 : : //! \param[in] nelem Total number of elements
610 : : //! \param[in] geoElem Element geometry array
611 : : //! \param[in] U Solution vector at recent time step
612 : : //! \param[in] P Vector of primitive quantities at recent time step
613 : : //! \param[in] ct Pressure relaxation time-scale for this system
614 : : //! \param[in,out] R Right-hand side vector added to
615 : : // *****************************************************************************
616 : : {
617 : : using inciter::volfracIdx;
618 : : using inciter::energyIdx;
619 : : using inciter::pressureIdx;
620 : : using inciter::velocityIdx;
621 : : using inciter::densityIdx;
622 : :
623 : 1568 : auto ncomp = U.nprop()/rdof;
624 : 1568 : auto nprim = P.nprop()/rdof;
625 : :
626 : : // compute volume integrals
627 [ + + ]: 146528 : for (std::size_t e=0; e<nelem; ++e)
628 : : {
629 [ + - ]: 144960 : auto dx = geoElem(e,4)/2.0;
630 : :
631 : : // Compute the basis function
632 [ + - ]: 289920 : std::vector< tk::real > B(rdof, 0.0);
633 : 144960 : B[0] = 1.0;
634 : :
635 : : auto state = evalPolynomialSol(mat_blk, 0, ncomp, nprim,
636 : : rdof, nmat, e, rdof, inpoel, coord, geoElem,
637 [ + - ]: 289920 : {{0.25, 0.25, 0.25}}, B, U, P);
638 : :
639 : : // get bulk properties
640 : 144960 : real rhob(0.0);
641 [ + + ]: 487440 : for (std::size_t k=0; k<nmat; ++k)
642 : 342480 : rhob += state[densityIdx(nmat, k)];
643 : :
644 : : // get pressures and bulk modulii
645 : 144960 : real pb(0.0), nume(0.0), deno(0.0), trelax(0.0);
646 [ + - ][ + - ]: 289920 : std::vector< real > apmat(nmat, 0.0), kmat(nmat, 0.0);
647 [ + + ]: 487440 : for (std::size_t k=0; k<nmat; ++k)
648 : : {
649 : 342480 : real arhomat = state[densityIdx(nmat, k)];
650 : 342480 : real alphamat = state[volfracIdx(nmat, k)];
651 : 342480 : apmat[k] = state[ncomp+pressureIdx(nmat, k)];
652 [ + - ]: 684960 : real amat = mat_blk[k].compute< inciter::EOS::soundspeed >( arhomat,
653 : 342480 : apmat[k], alphamat, k );
654 : 342480 : kmat[k] = arhomat * amat * amat;
655 : 342480 : pb += apmat[k];
656 : :
657 : : // relaxation parameters
658 : 342480 : trelax = std::max(trelax, ct*dx/amat);
659 : 342480 : nume += alphamat * apmat[k] / kmat[k];
660 : 342480 : deno += alphamat * alphamat / kmat[k];
661 : : }
662 : 144960 : auto p_relax = nume/deno;
663 : :
664 : : // compute pressure relaxation terms
665 [ + - ]: 289920 : std::vector< real > s_prelax(ncomp, 0.0);
666 [ + + ]: 487440 : for (std::size_t k=0; k<nmat; ++k)
667 : : {
668 : 342480 : auto s_alpha = (apmat[k]-p_relax*state[volfracIdx(nmat, k)])
669 : 342480 : * (state[volfracIdx(nmat, k)]/kmat[k]) / trelax;
670 : 342480 : s_prelax[volfracIdx(nmat, k)] = s_alpha;
671 : 342480 : s_prelax[energyIdx(nmat, k)] = - pb*s_alpha;
672 : : }
673 : :
674 [ + + ]: 1607280 : for (ncomp_t c=0; c<ncomp; ++c)
675 : : {
676 [ + - ][ + - ]: 1462320 : R(e, c) += geoElem(e,0) * s_prelax[c];
677 : : }
678 : : }
679 : 1568 : }
680 : :
681 : : std::vector< std::vector< tk::real > >
682 : 0 : solvevriem( std::size_t nelem,
683 : : const std::vector< std::vector< tk::real > >& vriem,
684 : : const std::vector< std::vector< tk::real > >& riemannLoc )
685 : : // *****************************************************************************
686 : : // Solve the reconstruct velocity used for volume fraction equation by
687 : : // Least square method
688 : : //! \param[in] nelem Numer of elements
689 : : //! \param[in,out] vriem Vector of the riemann velocity
690 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
691 : : //! is available
692 : : //! \return Vector of Riemann velocity polynomial solution
693 : : // *****************************************************************************
694 : : {
695 : : std::vector< std::vector< tk::real > >
696 [ - - ][ - - ]: 0 : vriempoly( nelem, std::vector<tk::real>(12,0.0) );
697 : :
698 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
699 : : {
700 : : // Use the normal method to construct the linear system A^T * A * x = u
701 : 0 : auto numgp = riemannLoc[e].size()/3;
702 : : std::vector< std::vector< tk::real > > A(numgp,
703 [ - - ][ - - ]: 0 : std::vector<tk::real>(4, 1.0));
704 : :
705 [ - - ]: 0 : for(std::size_t k = 0; k < numgp; k++)
706 : : {
707 : 0 : auto mark = k * 3;
708 : 0 : A[k][1] = riemannLoc[e][mark];
709 : 0 : A[k][2] = riemannLoc[e][mark+1];
710 : 0 : A[k][3] = riemannLoc[e][mark+2];
711 : : }
712 : :
713 [ - - ]: 0 : for(std::size_t idir = 0; idir < 3; idir++)
714 : : {
715 : : double AA_T[4*4], u[4];
716 : :
717 [ - - ]: 0 : for(std::size_t i = 0; i < 4; i++)
718 [ - - ]: 0 : for(std::size_t j = 0; j < 4; j++)
719 : : {
720 : 0 : auto id = 4 * i + j;
721 : 0 : AA_T[id] = 0;
722 [ - - ]: 0 : for(std::size_t k = 0; k < numgp; k++)
723 : 0 : AA_T[id] += A[k][i] * A[k][j];
724 : : }
725 : :
726 [ - - ]: 0 : std::vector<tk::real> vel(numgp, 1.0);
727 [ - - ]: 0 : for(std::size_t k = 0; k < numgp; k++)
728 : : {
729 : 0 : auto mark = k * 3 + idir;
730 : 0 : vel[k] = vriem[e][mark];
731 : : }
732 [ - - ]: 0 : for(std::size_t k = 0; k < 4; k++)
733 : : {
734 : 0 : u[k] = 0;
735 [ - - ]: 0 : for(std::size_t i = 0; i < numgp; i++)
736 : 0 : u[k] += A[i][k] * vel[i];
737 : : }
738 : :
739 : : lapack_int IPIV[4];
740 : : #ifndef NDEBUG
741 : : lapack_int info =
742 : : #endif
743 [ - - ]: 0 : LAPACKE_dsysv( LAPACK_ROW_MAJOR, 'U', 4, 1, AA_T, 4, IPIV, u, 1 );
744 [ - - ][ - - ]: 0 : Assert( info == 0, "Error in linear system solver" );
[ - - ][ - - ]
745 : :
746 : 0 : auto idirmark = idir * 4;
747 [ - - ]: 0 : for(std::size_t k = 0; k < 4; k++)
748 : 0 : vriempoly[e][idirmark+k] = u[k];
749 : : }
750 : : }
751 : 0 : return vriempoly;
752 : : }
753 : :
754 : 0 : void evaluRiemann( ncomp_t ncomp,
755 : : const int e_left,
756 : : const int e_right,
757 : : const std::size_t nmat,
758 : : const std::vector< tk::real >& fl,
759 : : const std::array< tk::real, 3 >& fn,
760 : : const std::array< tk::real, 3 >& gp,
761 : : const std::array< std::vector< tk::real >, 2 >& state,
762 : : std::vector< std::vector< tk::real > >& vriem,
763 : : std::vector< std::vector< tk::real > >& riemannLoc )
764 : : // *****************************************************************************
765 : : // Compute the riemann velocity at the interface
766 : : //! \param[in] ncomp Number of scalar components in this PDE system
767 : : //! \param[in] e_left Index for the left element
768 : : //! \param[in] e_right Index for the right element
769 : : //! \param[in] nmat Number of materials in this PDE system
770 : : //! \param[in] fn Face/Surface normal
771 : : //! \param[in] gp Gauss points coordinates
772 : : //! \param[in] fl Surface flux
773 : : //! \param[in] state Vector of state variables for left and right side
774 : : //! \param[in,out] vriem Vector of the riemann velocity
775 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
776 : : //! is available
777 : : // *****************************************************************************
778 : : {
779 : : using inciter::densityIdx;
780 : : using inciter::momentumIdx;
781 : :
782 : 0 : std::size_t el(0), er(0);
783 : 0 : el = static_cast< std::size_t >(e_left);
784 [ - - ]: 0 : if(e_right != -1)
785 : 0 : er = static_cast< std::size_t >(e_right);
786 : :
787 [ - - ]: 0 : riemannLoc[el].push_back( gp[0] );
788 [ - - ]: 0 : riemannLoc[el].push_back( gp[1] );
789 [ - - ]: 0 : riemannLoc[el].push_back( gp[2] );
790 : :
791 [ - - ]: 0 : if(e_right != -1)
792 : : {
793 [ - - ]: 0 : riemannLoc[er].push_back( gp[0] );
794 [ - - ]: 0 : riemannLoc[er].push_back( gp[1] );
795 [ - - ]: 0 : riemannLoc[er].push_back( gp[2] );
796 : : }
797 : :
798 : 0 : tk::real rhobl(0.0), rhobr(0.0);
799 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
800 : : {
801 : 0 : rhobl += state[0][densityIdx(nmat, k)];
802 : 0 : rhobr += state[1][densityIdx(nmat, k)];
803 : : }
804 : :
805 : 0 : auto ul = state[0][momentumIdx(nmat, 0)] / rhobl;
806 : 0 : auto vl = state[0][momentumIdx(nmat, 1)] / rhobl;
807 : 0 : auto wl = state[0][momentumIdx(nmat, 2)] / rhobl;
808 : :
809 : 0 : auto ur = state[1][momentumIdx(nmat, 0)] / rhobr;
810 : 0 : auto vr = state[1][momentumIdx(nmat, 1)] / rhobr;
811 : 0 : auto wr = state[1][momentumIdx(nmat, 2)] / rhobr;
812 : :
813 : : // Compute the normal velocities from left and right cells
814 : 0 : auto vnl = ul * fn[0] + vl * fn[1] + wl * fn[2];
815 : 0 : auto vnr = ur * fn[0] + vr * fn[1] + wr * fn[2];
816 : :
817 : : // The interface velocity is evaluated by adding the normal velocity which
818 : : // is taken from the Riemann solver and the tangential velocity which is
819 : : // evaluated as an average of the left and right cells
820 : 0 : auto urie = 0.5 * ((ul + ur) - fn[0] * (vnl + vnr)) + fl[ncomp+nmat] * fn[0];
821 : 0 : auto vrie = 0.5 * ((vl + vr) - fn[1] * (vnl + vnr)) + fl[ncomp+nmat] * fn[1];
822 : 0 : auto wrie = 0.5 * ((wl + wr) - fn[2] * (vnl + vnr)) + fl[ncomp+nmat] * fn[2];
823 : :
824 [ - - ]: 0 : vriem[el].push_back(urie);
825 [ - - ]: 0 : vriem[el].push_back(vrie);
826 [ - - ]: 0 : vriem[el].push_back(wrie);
827 : :
828 [ - - ]: 0 : if(e_right != -1)
829 : : {
830 [ - - ]: 0 : vriem[er].push_back(urie);
831 [ - - ]: 0 : vriem[er].push_back(vrie);
832 [ - - ]: 0 : vriem[er].push_back(wrie);
833 : : }
834 : 0 : }
835 : :
836 : : std::vector< std::array< tk::real, 3 > >
837 : 7374000 : fluxTerms(
838 : : std::size_t ncomp,
839 : : std::size_t nmat,
840 : : const std::vector< inciter::EOS >& /*mat_blk*/,
841 : : const std::vector< tk::real >& ugp )
842 : : // *****************************************************************************
843 : : // Compute the flux-function for the multimaterial PDEs
844 : : //! \param[in] ncomp Number of components in this PDE system
845 : : //! \param[in] nmat Number of materials in this PDE system
846 : : // //! \param[in] mat_blk EOS material block
847 : : //! \param[in] ugp State vector at the Gauss point at which flux is required
848 : : //! \return Flux vectors for all components in multi-material PDE system
849 : : // *****************************************************************************
850 : : {
851 : : using inciter::volfracIdx;
852 : : using inciter::densityIdx;
853 : : using inciter::momentumIdx;
854 : : using inciter::energyIdx;
855 : : using inciter::velocityIdx;
856 : : using inciter::pressureIdx;
857 : : using inciter::deformIdx;
858 : :
859 : : const auto& solidx =
860 : 7374000 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
861 : :
862 [ + - ]: 7374000 : std::vector< std::array< tk::real, 3 > > fl( ncomp, {{0, 0, 0}} );
863 : :
864 : 7374000 : tk::real rho(0.0);
865 [ + + ]: 22122000 : for (std::size_t k=0; k<nmat; ++k)
866 : 14748000 : rho += ugp[densityIdx(nmat, k)];
867 : :
868 : 7374000 : auto u = ugp[ncomp+velocityIdx(nmat,0)];
869 : 7374000 : auto v = ugp[ncomp+velocityIdx(nmat,1)];
870 : 7374000 : auto w = ugp[ncomp+velocityIdx(nmat,2)];
871 : :
872 [ + - ][ - + ]: 7374000 : if (inciter::haveSolid(nmat, solidx))
873 : : {
874 [ - - ]: 0 : std::vector< tk::real > al(nmat, 0.0);
875 : 0 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > g, asig;
876 : : std::array< std::array< tk::real, 3 >, 3 >
877 : 0 : sig {{ {{0, 0, 0}}, {{0, 0, 0}}, {{0, 0, 0}} }};
878 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
879 : : {
880 : 0 : al[k] = ugp[volfracIdx(nmat, k)];
881 : : // inv deformation gradient and Cauchy stress tensors
882 [ - - ][ - - ]: 0 : g.push_back(inciter::getDeformGrad(nmat, k, ugp));
883 [ - - ][ - - ]: 0 : asig.push_back(inciter::getCauchyStress(nmat, k, ncomp, ugp));
884 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) asig[k][i][i] -= ugp[ncomp+pressureIdx(nmat,k)];
885 : :
886 [ - - ]: 0 : for (size_t i=0; i<3; ++i)
887 [ - - ]: 0 : for (size_t j=0; j<3; ++j)
888 : 0 : sig[i][j] += asig[k][i][j];
889 : : }
890 : :
891 : : // conservative part of momentum flux
892 : 0 : fl[momentumIdx(nmat, 0)][0] = ugp[momentumIdx(nmat, 0)] * u - sig[0][0];
893 : 0 : fl[momentumIdx(nmat, 1)][0] = ugp[momentumIdx(nmat, 1)] * u - sig[0][1];
894 : 0 : fl[momentumIdx(nmat, 2)][0] = ugp[momentumIdx(nmat, 2)] * u - sig[0][2];
895 : :
896 : 0 : fl[momentumIdx(nmat, 0)][1] = ugp[momentumIdx(nmat, 0)] * v - sig[1][0];
897 : 0 : fl[momentumIdx(nmat, 1)][1] = ugp[momentumIdx(nmat, 1)] * v - sig[1][1];
898 : 0 : fl[momentumIdx(nmat, 2)][1] = ugp[momentumIdx(nmat, 2)] * v - sig[1][2];
899 : :
900 : 0 : fl[momentumIdx(nmat, 0)][2] = ugp[momentumIdx(nmat, 0)] * w - sig[2][0];
901 : 0 : fl[momentumIdx(nmat, 1)][2] = ugp[momentumIdx(nmat, 1)] * w - sig[2][1];
902 : 0 : fl[momentumIdx(nmat, 2)][2] = ugp[momentumIdx(nmat, 2)] * w - sig[2][2];
903 : :
904 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
905 : : {
906 : : // conservative part of volume-fraction flux
907 : 0 : fl[volfracIdx(nmat, k)][0] = 0.0;
908 : 0 : fl[volfracIdx(nmat, k)][1] = 0.0;
909 : 0 : fl[volfracIdx(nmat, k)][2] = 0.0;
910 : :
911 : : // conservative part of material continuity flux
912 : 0 : fl[densityIdx(nmat, k)][0] = u * ugp[densityIdx(nmat, k)];
913 : 0 : fl[densityIdx(nmat, k)][1] = v * ugp[densityIdx(nmat, k)];
914 : 0 : fl[densityIdx(nmat, k)][2] = w * ugp[densityIdx(nmat, k)];
915 : :
916 : : // conservative part of material total-energy flux
917 : 0 : fl[energyIdx(nmat, k)][0] = u * ugp[energyIdx(nmat, k)]
918 : 0 : - u * asig[k][0][0] - v * asig[k][1][0] - w * asig[k][2][0];
919 : 0 : fl[energyIdx(nmat, k)][1] = v * ugp[energyIdx(nmat, k)]
920 : 0 : - u * asig[k][0][1] - v * asig[k][1][1] - w * asig[k][2][1];
921 : 0 : fl[energyIdx(nmat, k)][2] = w * ugp[energyIdx(nmat, k)]
922 : 0 : - u * asig[k][0][2] - v * asig[k][1][2] - w * asig[k][2][2];
923 : :
924 : : // conservative part of material inverse deformation gradient
925 : : // g_ij: \partial (g_il u_l) / \partial (x_j)
926 [ - - ]: 0 : if (solidx[k] > 0)
927 : : {
928 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
929 : : {
930 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
931 : : {
932 : 0 : fl[deformIdx(nmat,solidx[k],i,j)][j] =
933 : 0 : u*g[k][i][0] + v*g[k][i][1] + w*g[k][i][2];
934 : : }
935 : : // other components are zero
936 : : }
937 : : }
938 : : }
939 : : }
940 : : else
941 : : {
942 : 7374000 : tk::real p(0.0);
943 [ + - ]: 14748000 : std::vector< tk::real > apk( nmat, 0.0 );
944 [ + + ]: 22122000 : for (std::size_t k=0; k<nmat; ++k)
945 : : {
946 : 14748000 : apk[k] = ugp[ncomp+pressureIdx(nmat,k)];
947 : 14748000 : p += apk[k];
948 : : }
949 : :
950 : : // conservative part of momentum flux
951 : 7374000 : fl[momentumIdx(nmat, 0)][0] = ugp[momentumIdx(nmat, 0)] * u + p;
952 : 7374000 : fl[momentumIdx(nmat, 1)][0] = ugp[momentumIdx(nmat, 1)] * u;
953 : 7374000 : fl[momentumIdx(nmat, 2)][0] = ugp[momentumIdx(nmat, 2)] * u;
954 : :
955 : 7374000 : fl[momentumIdx(nmat, 0)][1] = ugp[momentumIdx(nmat, 0)] * v;
956 : 7374000 : fl[momentumIdx(nmat, 1)][1] = ugp[momentumIdx(nmat, 1)] * v + p;
957 : 7374000 : fl[momentumIdx(nmat, 2)][1] = ugp[momentumIdx(nmat, 2)] * v;
958 : :
959 : 7374000 : fl[momentumIdx(nmat, 0)][2] = ugp[momentumIdx(nmat, 0)] * w;
960 : 7374000 : fl[momentumIdx(nmat, 1)][2] = ugp[momentumIdx(nmat, 1)] * w;
961 : 7374000 : fl[momentumIdx(nmat, 2)][2] = ugp[momentumIdx(nmat, 2)] * w + p;
962 : :
963 [ + + ]: 22122000 : for (std::size_t k=0; k<nmat; ++k)
964 : : {
965 : : // conservative part of volume-fraction flux
966 : 14748000 : fl[volfracIdx(nmat, k)][0] = 0.0;
967 : 14748000 : fl[volfracIdx(nmat, k)][1] = 0.0;
968 : 14748000 : fl[volfracIdx(nmat, k)][2] = 0.0;
969 : :
970 : : // conservative part of material continuity flux
971 : 14748000 : fl[densityIdx(nmat, k)][0] = u * ugp[densityIdx(nmat, k)];
972 : 14748000 : fl[densityIdx(nmat, k)][1] = v * ugp[densityIdx(nmat, k)];
973 : 14748000 : fl[densityIdx(nmat, k)][2] = w * ugp[densityIdx(nmat, k)];
974 : :
975 : : // conservative part of material total-energy flux
976 : 14748000 : auto hmat = ugp[energyIdx(nmat, k)] + apk[k];
977 : 14748000 : fl[energyIdx(nmat, k)][0] = u * hmat;
978 : 14748000 : fl[energyIdx(nmat, k)][1] = v * hmat;
979 : 14748000 : fl[energyIdx(nmat, k)][2] = w * hmat;
980 : : }
981 : : }
982 : :
983 : 7374000 : return fl;
984 : : }
985 : :
986 : : }// tk::
|
