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