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 : 6705 : 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 : 6705 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
97 : :
98 : 6705 : const auto& cx = coord[0];
99 : 6705 : const auto& cy = coord[1];
100 : 6705 : const auto& cz = coord[2];
101 : :
102 : 6705 : auto ncomp = U.nprop()/rdof;
103 : 6705 : auto nprim = P.nprop()/rdof;
104 : :
105 : : // compute volume integrals
106 [ + + ]: 2504385 : for (std::size_t e=0; e<nelem; ++e)
107 : : {
108 [ + - ]: 2497680 : auto ng = tk::NGvol(ndofel[e]);
109 : :
110 : : // arrays for quadrature points
111 : 4995360 : std::array< std::vector< real >, 3 > coordgp;
112 : 4995360 : std::vector< real > wgp;
113 : :
114 [ + - ]: 2497680 : coordgp[0].resize( ng );
115 [ + - ]: 2497680 : coordgp[1].resize( ng );
116 [ + - ]: 2497680 : coordgp[2].resize( ng );
117 [ + - ]: 2497680 : wgp.resize( ng );
118 : :
119 [ + - ]: 2497680 : GaussQuadratureTet( ng, coordgp, wgp );
120 : :
121 : : // Extract the element coordinates
122 : : std::array< std::array< real, 3>, 4 > coordel {{
123 : 7493040 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
124 : 7493040 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
125 : 7493040 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
126 : 7493040 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
127 : 27474480 : }};
128 : :
129 : : auto jacInv =
130 : 2497680 : 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 [ + + ]: 2497680 : if (rdof > ndof)
137 : : {
138 : 341100 : dof_el = rdof;
139 : : }
140 : : else
141 : : {
142 : 2156580 : 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 [ + + ][ - + ]: 2497680 : if(dof_el == 1 && pref)
148 : 0 : dof_el = 4;
149 : :
150 : : // Compute the derivatives of basis function for second order terms
151 : 4995360 : std::array< std::vector<tk::real>, 3 > dBdx;
152 [ + + ]: 2497680 : if (ndofel[e] > 1)
153 [ + - ]: 500280 : dBdx = eval_dBdx_p1( ndofel[e], jacInv );
154 : :
155 : : // Gaussian quadrature
156 [ + + ]: 6996480 : for (std::size_t igp=0; igp<ng; ++igp)
157 : : {
158 [ - + ]: 4498800 : if (ndofel[e] > 4)
159 [ - - ]: 0 : eval_dBdx_p2( igp, coordgp, jacInv, dBdx );
160 : :
161 : : // Compute the basis function
162 : : auto B =
163 [ + - ]: 8997600 : eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
164 : :
165 [ + - ]: 4498800 : 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 [ + - ]: 8997600 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, U, P);
170 : :
171 : : // get bulk properties
172 : 4498800 : tk::real rhob(0.0);
173 [ + + ]: 14589300 : for (std::size_t k=0; k<nmat; ++k)
174 : 10090500 : rhob += state[densityIdx(nmat, k)];
175 : :
176 : : // get the velocity vector
177 : 4498800 : std::array< tk::real, 3 > vel{{ state[ncomp+velocityIdx(nmat, 0)],
178 : 4498800 : state[ncomp+velocityIdx(nmat, 1)],
179 : 8997600 : state[ncomp+velocityIdx(nmat, 2)] }};
180 : :
181 [ + - ]: 8997600 : std::vector< tk::real > ymat(nmat, 0.0);
182 : 4498800 : std::array< tk::real, 3 > dap{{0.0, 0.0, 0.0}};
183 [ + + ]: 14589300 : for (std::size_t k=0; k<nmat; ++k)
184 : : {
185 : 10090500 : ymat[k] = state[densityIdx(nmat, k)]/rhob;
186 : :
187 : 10090500 : std::size_t mark(3*k);
188 [ - + ]: 10090500 : if (solidx[k] > 0) mark = 3*nmat+ndof+3*(solidx[k]-1);
189 : :
190 [ + + ]: 40362000 : for (std::size_t idir=0; idir<3; ++idir)
191 : 30271500 : dap[idir] += riemannDeriv[mark+idir][e];
192 : : }
193 : :
194 : : // compute non-conservative terms
195 : : std::vector< std::vector< tk::real > > ncf
196 [ + - ][ + - ]: 13496400 : (ncomp, std::vector<tk::real>(ndofel[e],0.0));
197 : :
198 [ + + ]: 17995200 : for (std::size_t idir=0; idir<3; ++idir)
199 [ + + ]: 49505400 : for(std::size_t idof=0; idof<ndofel[e]; ++idof)
200 : 36009000 : ncf[momentumIdx(nmat, idir)][idof] = 0.0;
201 : :
202 [ + + ]: 14589300 : for (std::size_t k=0; k<nmat; ++k)
203 : : {
204 : : // evaluate non-conservative term for energy equation
205 : 10090500 : std::size_t mark(3*k);
206 [ - + ]: 10090500 : if (solidx[k] > 0) mark = 3*nmat+ndof+3*(solidx[k]-1);
207 : :
208 [ + + ]: 35189400 : for(std::size_t idof=0; idof<ndofel[e]; ++idof)
209 : : {
210 : 25098900 : ncf[densityIdx(nmat, k)][idof] = 0.0;
211 : :
212 [ + + ]: 100395600 : for (std::size_t idir=0; idir<3; ++idir)
213 : 150593400 : ncf[energyIdx(nmat, k)][idof] -= vel[idir] * ( ymat[k]*dap[idir]
214 : 75296700 : - riemannDeriv[mark+idir][e] );
215 : : }
216 : :
217 : : // evaluate non-conservative terms for g equation
218 [ - + ]: 10090500 : 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 [ - + ][ - - ]: 10090500 : if (ndof <= 4 || intsharp == 1) {
248 [ + + ]: 35189400 : for(std::size_t idof=0; idof<ndofel[e]; ++idof)
249 : 50197800 : ncf[volfracIdx(nmat, k)][idof] = state[volfracIdx(nmat, k)]
250 : 35189400 : * 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 [ + - ]: 4498800 : updateRhsNonCons( ncomp, nmat, ndof, ndofel[e], wt, e, B, dBdx, ncf, R );
268 : : }
269 : : }
270 : 6705 : }
271 : :
272 : : void
273 : 4498800 : 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 [ - + ][ - - ]: 4498800 : Assert( ncf.size() == ncomp, "Size mismatch for non-conservative term" );
[ - - ][ - - ]
314 : :
315 [ + + ]: 48266700 : for (ncomp_t c=0; c<ncomp; ++c)
316 : : {
317 : 43767900 : auto mark = c*ndof;
318 : 43767900 : R(e, mark) += wt * ncf[c][0];
319 : : }
320 : :
321 [ + + ]: 4498800 : if( ndof_el > 1)
322 : : {
323 : : // Update rhs with distributions from volume fraction and energy equations
324 [ + + ]: 7504200 : for (std::size_t k=0; k<nmat; ++k)
325 : : {
326 [ + + ]: 20011200 : for(std::size_t idof = 1; idof < ndof_el; idof++)
327 : : {
328 : 15008400 : R(e, volfracDofIdx(nmat,k,ndof,idof)) +=
329 : 15008400 : wt * ncf[volfracIdx(nmat,k)][idof];
330 : 15008400 : R(e, energyDofIdx(nmat,k,ndof,idof)) +=
331 : 15008400 : wt * ncf[energyIdx(nmat,k)][idof] * B[idof];
332 : : // Note: High order non-conservative g terms not implemented yet!
333 : : }
334 : : }
335 : : }
336 : 4498800 : }
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 : :
388 : : // compute non-conservative terms
389 [ + - ]: 1091856 : std::vector< tk::real > ncf(ncomp, 0.0);
390 [ + - ]: 2183712 : std::vector< tk::real > ymat(nmat, 0.0);
391 [ + + ]: 3512624 : for (std::size_t k=0; k<nmat; ++k)
392 : : {
393 [ + - ]: 2420768 : ymat[k] = U(e, densityDofIdx(nmat,k,rdof,0))/rhob;
394 : :
395 : : // evaluate non-conservative term for energy equation
396 [ + + ]: 9683072 : for (std::size_t idir=0; idir<3; ++idir)
397 : 14524608 : ncf[energyIdx(nmat, k)] -= vel[idir] * ( ymat[k]*p_face*fn[idir]
398 : 7262304 : - var_riemann[k]*fn[idir] );
399 : :
400 : : // evaluate non-conservative term for volume fraction equation
401 [ + - ]: 4841536 : ncf[volfracIdx(nmat, k)] = U(e, volfracDofIdx(nmat,k,rdof,0))
402 : 2420768 : * var_riemann[nmat];
403 : : }
404 : :
405 : 2183712 : return ncf;
406 : : }
407 : :
408 : : void
409 : 375 : pressureRelaxationInt( const bool pref,
410 : : std::size_t nmat,
411 : : const std::vector< inciter::EOS >& mat_blk,
412 : : const std::size_t ndof,
413 : : const std::size_t rdof,
414 : : const std::size_t nelem,
415 : : const std::vector< std::size_t >& inpoel,
416 : : const UnsMesh::Coords& coord,
417 : : const Fields& geoElem,
418 : : const Fields& U,
419 : : const Fields& P,
420 : : const std::vector< std::size_t >& ndofel,
421 : : const tk::real ct,
422 : : Fields& R,
423 : : int intsharp )
424 : : // *****************************************************************************
425 : : // Compute volume integrals of pressure relaxation terms in multi-material DG
426 : : //! \details This is called for multi-material DG to compute volume integrals of
427 : : //! finite pressure relaxation terms in the volume fraction and energy
428 : : //! equations, which do not exist in the single-material flow formulation (for
429 : : //! `CompFlow` DG). For further details see Dobrev, V. A., Kolev, T. V.,
430 : : //! Rieben, R. N., & Tomov, V. Z. (2016). Multi‐material closure model for
431 : : //! high‐order finite element Lagrangian hydrodynamics. International Journal
432 : : //! for Numerical Methods in Fluids, 82(10), 689-706.
433 : : //! \param[in] pref Indicator for p-adaptive algorithm
434 : : //! \param[in] nmat Number of materials in this PDE system
435 : : //! \param[in] mat_blk EOS material block
436 : : //! \param[in] ndof Maximum number of degrees of freedom
437 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
438 : : //! \param[in] nelem Total number of elements
439 : : //! \param[in] inpoel Element-node connectivity
440 : : //! \param[in] coord Array of nodal coordinates
441 : : //! \param[in] geoElem Element geometry array
442 : : //! \param[in] U Solution vector at recent time step
443 : : //! \param[in] P Vector of primitive quantities at recent time step
444 : : //! \param[in] ndofel Vector of local number of degrees of freedome
445 : : //! \param[in] ct Pressure relaxation time-scale for this system
446 : : //! \param[in,out] R Right-hand side vector added to
447 : : //! \param[in] intsharp Interface reconstruction indicator
448 : : // *****************************************************************************
449 : : {
450 : : using inciter::volfracIdx;
451 : : using inciter::densityIdx;
452 : : using inciter::momentumIdx;
453 : : using inciter::energyIdx;
454 : : using inciter::pressureIdx;
455 : : using inciter::velocityIdx;
456 : : using inciter::deformIdx;
457 : :
458 : : const auto& solidx =
459 : 375 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
460 : :
461 : 375 : auto ncomp = U.nprop()/rdof;
462 : 375 : auto nprim = P.nprop()/rdof;
463 : :
464 : : // compute volume integrals
465 [ + + ]: 227775 : for (std::size_t e=0; e<nelem; ++e)
466 : : {
467 [ + - ]: 227400 : auto dx = geoElem(e,4)/2.0;
468 [ + - ]: 227400 : auto ng = NGvol(ndofel[e]);
469 : :
470 : : // arrays for quadrature points
471 : 454800 : std::array< std::vector< real >, 3 > coordgp;
472 : 454800 : std::vector< real > wgp;
473 : :
474 [ + - ]: 227400 : coordgp[0].resize( ng );
475 [ + - ]: 227400 : coordgp[1].resize( ng );
476 [ + - ]: 227400 : coordgp[2].resize( ng );
477 [ + - ]: 227400 : wgp.resize( ng );
478 : :
479 [ + - ]: 227400 : GaussQuadratureTet( ng, coordgp, wgp );
480 : :
481 : : // Compute the derivatives of basis function for DG(P1)
482 : 454800 : std::array< std::vector<real>, 3 > dBdx;
483 : :
484 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
485 : : // basis functions and solutions by default. This implies that P0P1 is
486 : : // unsupported in the p-adaptive DG (PDG).
487 : : std::size_t dof_el;
488 [ - + ]: 227400 : if (rdof > ndof)
489 : : {
490 : 0 : dof_el = rdof;
491 : : }
492 : : else
493 : : {
494 : 227400 : dof_el = ndofel[e];
495 : : }
496 : :
497 : : // For multi-material p-adaptive simulation, when dofel = 1, p0p1 is applied
498 : : // and ndof for solution evaluation should be 4
499 [ - + ][ - - ]: 227400 : if(dof_el == 1 && pref)
500 : 0 : dof_el = 4;
501 : :
502 : : // Gaussian quadrature
503 [ + + ]: 1364400 : for (std::size_t igp=0; igp<ng; ++igp)
504 : : {
505 : : // Compute the basis function
506 : : auto B =
507 [ + - ]: 2274000 : eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
508 : :
509 [ + - ]: 1137000 : auto wt = wgp[igp] * geoElem(e, 0);
510 : :
511 : : auto state = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim,
512 : : rdof, nmat, e, dof_el, inpoel, coord, geoElem,
513 [ + - ]: 2274000 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, U, P);
514 : :
515 : : // get bulk properties
516 : 1137000 : real rhob(0.0);
517 [ + + ]: 3411000 : for (std::size_t k=0; k<nmat; ++k)
518 : 2274000 : rhob += state[densityIdx(nmat, k)];
519 : :
520 : : // get pressures and bulk modulii
521 : 1137000 : real pb(0.0), nume(0.0), deno(0.0), trelax(0.0);
522 [ + - ][ + - ]: 2274000 : std::vector< real > apmat(nmat, 0.0), kmat(nmat, 0.0);
523 [ + - ]: 2274000 : std::vector< int > do_relax(nmat, 1);
524 : 1137000 : bool is_relax(false);
525 [ + + ]: 3411000 : for (std::size_t k=0; k<nmat; ++k)
526 : : {
527 : 2274000 : real arhomat = state[densityIdx(nmat, k)];
528 : 2274000 : real alphamat = state[volfracIdx(nmat, k)];
529 : 2274000 : apmat[k] = state[ncomp+pressureIdx(nmat, k)];
530 : 2274000 : real amat = 0.0;
531 [ + - ][ + + ]: 2274000 : if (solidx[k] == 0 && alphamat >= inciter::volfracPRelaxLim()) {
[ + + ]
532 [ + - ]: 2337396 : amat = mat_blk[k].compute< inciter::EOS::soundspeed >( arhomat,
533 : 1168698 : apmat[k], alphamat, k );
534 : 1168698 : kmat[k] = arhomat * amat * amat;
535 : 1168698 : pb += apmat[k];
536 : :
537 : : // relaxation parameters
538 : 1168698 : trelax = std::max(trelax, ct*dx/amat);
539 : 1168698 : nume += alphamat * apmat[k] / kmat[k];
540 : 1168698 : deno += alphamat * alphamat / kmat[k];
541 : :
542 : 1168698 : is_relax = true;
543 : : }
544 : 1105302 : else do_relax[k] = 0;
545 : : }
546 : 1137000 : real p_relax(0.0);
547 [ + - ]: 1137000 : if (is_relax) p_relax = nume/deno;
548 : :
549 : : // compute pressure relaxation terms
550 [ + - ]: 2274000 : std::vector< real > s_prelax(ncomp, 0.0);
551 [ + + ]: 3411000 : for (std::size_t k=0; k<nmat; ++k)
552 : : {
553 : : // only perform prelax on existing quantities
554 [ + + ]: 2274000 : if (do_relax[k] == 1) {
555 : 1168698 : auto s_alpha = (apmat[k]-p_relax*state[volfracIdx(nmat, k)])
556 : 1168698 : * (state[volfracIdx(nmat, k)]/kmat[k]) / trelax;
557 : 1168698 : s_prelax[volfracIdx(nmat, k)] = s_alpha;
558 : 1168698 : s_prelax[energyIdx(nmat, k)] = - pb*s_alpha;
559 : : }
560 : : }
561 : :
562 [ + - ]: 1137000 : updateRhsPre( ncomp, ndof, ndofel[e], wt, e, B, s_prelax, R );
563 : : }
564 : : }
565 : 375 : }
566 : :
567 : : void
568 : 1137000 : updateRhsPre(
569 : : ncomp_t ncomp,
570 : : const std::size_t ndof,
571 : : [[maybe_unused]] const std::size_t ndof_el,
572 : : const tk::real wt,
573 : : const std::size_t e,
574 : : const std::vector< tk::real >& B,
575 : : std::vector< tk::real >& ncf,
576 : : Fields& R )
577 : : // *****************************************************************************
578 : : // Update the rhs by adding the pressure relaxation integrals
579 : : //! \param[in] ncomp Number of scalar components in this PDE system
580 : : //! \param[in] ndof Maximum number of degrees of freedom
581 : : //! \param[in] ndof_el Number of degrees of freedom for local element
582 : : //! \param[in] wt Weight of gauss quadrature point
583 : : //! \param[in] e Element index
584 : : //! \param[in] B Basis function evaluated at local quadrature point
585 : : //! \param[in] ncf Vector of non-conservative terms
586 : : //! \param[in,out] R Right-hand side vector computed
587 : : // *****************************************************************************
588 : : {
589 : : //Assert( dBdx[0].size() == ndof_el,
590 : : // "Size mismatch for basis function derivatives" );
591 : : //Assert( dBdx[1].size() == ndof_el,
592 : : // "Size mismatch for basis function derivatives" );
593 : : //Assert( dBdx[2].size() == ndof_el,
594 : : // "Size mismatch for basis function derivatives" );
595 : : //Assert( ncf.size() == ncomp,
596 : : // "Size mismatch for non-conservative term" );
597 [ - + ][ - - ]: 1137000 : Assert( ncf.size() == ncomp, "Size mismatch for non-conservative term" );
[ - - ][ - - ]
598 : :
599 [ + + ]: 11370000 : for (ncomp_t c=0; c<ncomp; ++c)
600 : : {
601 : 10233000 : auto mark = c*ndof;
602 [ + + ]: 51165000 : for(std::size_t idof = 0; idof < ndof; idof++)
603 : 40932000 : R(e, mark+idof) += wt * ncf[c] * B[idof];
604 : : }
605 : 1137000 : }
606 : :
607 : : void
608 : 1218 : pressureRelaxationIntFV(
609 : : std::size_t nmat,
610 : : const std::vector< inciter::EOS >& mat_blk,
611 : : const std::size_t rdof,
612 : : const std::size_t nelem,
613 : : const std::vector< std::size_t >& inpoel,
614 : : const UnsMesh::Coords& coord,
615 : : const Fields& geoElem,
616 : : const Fields& U,
617 : : const Fields& P,
618 : : const tk::real ct,
619 : : Fields& R )
620 : : // *****************************************************************************
621 : : // Compute volume integrals of pressure relaxation terms in multi-material FV
622 : : //! \details This is called for multi-material FV to compute volume integrals of
623 : : //! finite pressure relaxation terms in the volume fraction and energy
624 : : //! equations, which do not exist in the single-material flow formulation (for
625 : : //! `CompFlow`). For further details see Dobrev, V. A., Kolev, T. V.,
626 : : //! Rieben, R. N., & Tomov, V. Z. (2016). Multi‐material closure model for
627 : : //! high‐order finite element Lagrangian hydrodynamics. International Journal
628 : : //! for Numerical Methods in Fluids, 82(10), 689-706.
629 : : //! \param[in] nmat Number of materials in this PDE system
630 : : //! \param[in] mat_blk EOS material block
631 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
632 : : //! \param[in] nelem Total number of elements
633 : : //! \param[in] geoElem Element geometry array
634 : : //! \param[in] U Solution vector at recent time step
635 : : //! \param[in] P Vector of primitive quantities at recent time step
636 : : //! \param[in] ct Pressure relaxation time-scale for this system
637 : : //! \param[in,out] R Right-hand side vector added to
638 : : // *****************************************************************************
639 : : {
640 : : using inciter::volfracIdx;
641 : : using inciter::energyIdx;
642 : : using inciter::pressureIdx;
643 : : using inciter::velocityIdx;
644 : : using inciter::densityIdx;
645 : :
646 : 1218 : auto ncomp = U.nprop()/rdof;
647 : 1218 : auto nprim = P.nprop()/rdof;
648 : :
649 : : // compute volume integrals
650 [ + + ]: 129578 : for (std::size_t e=0; e<nelem; ++e)
651 : : {
652 [ + - ]: 128360 : auto dx = geoElem(e,4)/2.0;
653 : :
654 : : // Compute the basis function
655 [ + - ]: 256720 : std::vector< tk::real > B(rdof, 0.0);
656 : 128360 : B[0] = 1.0;
657 : :
658 : : auto state = evalPolynomialSol(mat_blk, 0, ncomp, nprim,
659 : : rdof, nmat, e, rdof, inpoel, coord, geoElem,
660 [ + - ]: 256720 : {{0.25, 0.25, 0.25}}, B, U, P);
661 : :
662 : : // get bulk properties
663 : 128360 : real rhob(0.0);
664 [ + + ]: 437640 : for (std::size_t k=0; k<nmat; ++k)
665 : 309280 : rhob += state[densityIdx(nmat, k)];
666 : :
667 : : // get pressures and bulk modulii
668 : 128360 : real pb(0.0), nume(0.0), deno(0.0), trelax(0.0);
669 [ + - ][ + - ]: 256720 : std::vector< real > apmat(nmat, 0.0), kmat(nmat, 0.0);
670 [ + - ]: 256720 : std::vector< int > do_relax(nmat, 1);
671 : 128360 : bool is_relax(false);
672 [ + + ]: 437640 : for (std::size_t k=0; k<nmat; ++k)
673 : : {
674 : 309280 : real arhomat = state[densityIdx(nmat, k)];
675 : 309280 : real alphamat = state[volfracIdx(nmat, k)];
676 [ + + ]: 309280 : if (alphamat >= inciter::volfracPRelaxLim()) {
677 : 128874 : apmat[k] = state[ncomp+pressureIdx(nmat, k)];
678 [ + - ]: 257748 : real amat = mat_blk[k].compute< inciter::EOS::soundspeed >( arhomat,
679 : 128874 : apmat[k], alphamat, k );
680 : 128874 : kmat[k] = arhomat * amat * amat;
681 : 128874 : pb += apmat[k];
682 : :
683 : : // relaxation parameters
684 : 128874 : trelax = std::max(trelax, ct*dx/amat);
685 : 128874 : nume += alphamat * apmat[k] / kmat[k];
686 : 128874 : deno += alphamat * alphamat / kmat[k];
687 : :
688 : 128874 : is_relax = true;
689 : : }
690 : 180406 : else do_relax[k] = 0;
691 : : }
692 : 128360 : real p_relax(0.0);
693 [ + - ]: 128360 : if (is_relax) p_relax = nume/deno;
694 : :
695 : : // compute pressure relaxation terms
696 [ + - ]: 256720 : std::vector< real > s_prelax(ncomp, 0.0);
697 : :
698 [ + + ]: 437640 : for (std::size_t k=0; k<nmat; ++k)
699 : : {
700 : : // only perform prelax on existing quantities
701 [ + + ]: 309280 : if (do_relax[k] == 1) {
702 : 128874 : auto s_alpha = (apmat[k]-p_relax*state[volfracIdx(nmat, k)])
703 : 128874 : * (state[volfracIdx(nmat, k)]/kmat[k]) / trelax;
704 : 128874 : s_prelax[volfracIdx(nmat, k)] = s_alpha;
705 : 128874 : s_prelax[energyIdx(nmat, k)] = - pb*s_alpha;
706 : : }
707 : : }
708 : :
709 [ + + ]: 1441280 : for (ncomp_t c=0; c<ncomp; ++c)
710 : : {
711 [ + - ][ + - ]: 1312920 : R(e, c) += geoElem(e,0) * s_prelax[c];
712 : : }
713 : : }
714 : 1218 : }
715 : :
716 : : std::vector< std::vector< tk::real > >
717 : 0 : solvevriem( std::size_t nelem,
718 : : const std::vector< std::vector< tk::real > >& vriem,
719 : : const std::vector< std::vector< tk::real > >& riemannLoc )
720 : : // *****************************************************************************
721 : : // Solve the reconstruct velocity used for volume fraction equation by
722 : : // Least square method
723 : : //! \param[in] nelem Numer of elements
724 : : //! \param[in,out] vriem Vector of the riemann velocity
725 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
726 : : //! is available
727 : : //! \return Vector of Riemann velocity polynomial solution
728 : : // *****************************************************************************
729 : : {
730 : : std::vector< std::vector< tk::real > >
731 [ - - ][ - - ]: 0 : vriempoly( nelem, std::vector<tk::real>(12,0.0) );
732 : :
733 [ - - ]: 0 : for (std::size_t e=0; e<nelem; ++e)
734 : : {
735 : : // Use the normal method to construct the linear system A^T * A * x = u
736 : 0 : auto numgp = riemannLoc[e].size()/3;
737 : : std::vector< std::vector< tk::real > > A(numgp,
738 [ - - ][ - - ]: 0 : std::vector<tk::real>(4, 1.0));
739 : :
740 [ - - ]: 0 : for(std::size_t k = 0; k < numgp; k++)
741 : : {
742 : 0 : auto mark = k * 3;
743 : 0 : A[k][1] = riemannLoc[e][mark];
744 : 0 : A[k][2] = riemannLoc[e][mark+1];
745 : 0 : A[k][3] = riemannLoc[e][mark+2];
746 : : }
747 : :
748 [ - - ]: 0 : for(std::size_t idir = 0; idir < 3; idir++)
749 : : {
750 : : double AA_T[4*4], u[4];
751 : :
752 [ - - ]: 0 : for(std::size_t i = 0; i < 4; i++)
753 [ - - ]: 0 : for(std::size_t j = 0; j < 4; j++)
754 : : {
755 : 0 : auto id = 4 * i + j;
756 : 0 : AA_T[id] = 0;
757 [ - - ]: 0 : for(std::size_t k = 0; k < numgp; k++)
758 : 0 : AA_T[id] += A[k][i] * A[k][j];
759 : : }
760 : :
761 [ - - ]: 0 : std::vector<tk::real> vel(numgp, 1.0);
762 [ - - ]: 0 : for(std::size_t k = 0; k < numgp; k++)
763 : : {
764 : 0 : auto mark = k * 3 + idir;
765 : 0 : vel[k] = vriem[e][mark];
766 : : }
767 [ - - ]: 0 : for(std::size_t k = 0; k < 4; k++)
768 : : {
769 : 0 : u[k] = 0;
770 [ - - ]: 0 : for(std::size_t i = 0; i < numgp; i++)
771 : 0 : u[k] += A[i][k] * vel[i];
772 : : }
773 : :
774 : : lapack_int IPIV[4];
775 : : #ifndef NDEBUG
776 : : lapack_int info =
777 : : #endif
778 [ - - ]: 0 : LAPACKE_dsysv( LAPACK_ROW_MAJOR, 'U', 4, 1, AA_T, 4, IPIV, u, 1 );
779 [ - - ][ - - ]: 0 : Assert( info == 0, "Error in linear system solver" );
[ - - ][ - - ]
780 : :
781 : 0 : auto idirmark = idir * 4;
782 [ - - ]: 0 : for(std::size_t k = 0; k < 4; k++)
783 : 0 : vriempoly[e][idirmark+k] = u[k];
784 : : }
785 : : }
786 : 0 : return vriempoly;
787 : : }
788 : :
789 : 0 : void evaluRiemann( ncomp_t ncomp,
790 : : const int e_left,
791 : : const int e_right,
792 : : const std::size_t nmat,
793 : : const std::vector< tk::real >& fl,
794 : : const std::array< tk::real, 3 >& fn,
795 : : const std::array< tk::real, 3 >& gp,
796 : : const std::array< std::vector< tk::real >, 2 >& state,
797 : : std::vector< std::vector< tk::real > >& vriem,
798 : : std::vector< std::vector< tk::real > >& riemannLoc )
799 : : // *****************************************************************************
800 : : // Compute the riemann velocity at the interface
801 : : //! \param[in] ncomp Number of scalar components in this PDE system
802 : : //! \param[in] e_left Index for the left element
803 : : //! \param[in] e_right Index for the right element
804 : : //! \param[in] nmat Number of materials in this PDE system
805 : : //! \param[in] fn Face/Surface normal
806 : : //! \param[in] gp Gauss points coordinates
807 : : //! \param[in] fl Surface flux
808 : : //! \param[in] state Vector of state variables for left and right side
809 : : //! \param[in,out] vriem Vector of the riemann velocity
810 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
811 : : //! is available
812 : : // *****************************************************************************
813 : : {
814 : : using inciter::densityIdx;
815 : : using inciter::momentumIdx;
816 : :
817 : 0 : std::size_t el(0), er(0);
818 : 0 : el = static_cast< std::size_t >(e_left);
819 [ - - ]: 0 : if(e_right != -1)
820 : 0 : er = static_cast< std::size_t >(e_right);
821 : :
822 [ - - ]: 0 : riemannLoc[el].push_back( gp[0] );
823 [ - - ]: 0 : riemannLoc[el].push_back( gp[1] );
824 [ - - ]: 0 : riemannLoc[el].push_back( gp[2] );
825 : :
826 [ - - ]: 0 : if(e_right != -1)
827 : : {
828 [ - - ]: 0 : riemannLoc[er].push_back( gp[0] );
829 [ - - ]: 0 : riemannLoc[er].push_back( gp[1] );
830 [ - - ]: 0 : riemannLoc[er].push_back( gp[2] );
831 : : }
832 : :
833 : 0 : tk::real rhobl(0.0), rhobr(0.0);
834 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
835 : : {
836 : 0 : rhobl += state[0][densityIdx(nmat, k)];
837 : 0 : rhobr += state[1][densityIdx(nmat, k)];
838 : : }
839 : :
840 : 0 : auto ul = state[0][momentumIdx(nmat, 0)] / rhobl;
841 : 0 : auto vl = state[0][momentumIdx(nmat, 1)] / rhobl;
842 : 0 : auto wl = state[0][momentumIdx(nmat, 2)] / rhobl;
843 : :
844 : 0 : auto ur = state[1][momentumIdx(nmat, 0)] / rhobr;
845 : 0 : auto vr = state[1][momentumIdx(nmat, 1)] / rhobr;
846 : 0 : auto wr = state[1][momentumIdx(nmat, 2)] / rhobr;
847 : :
848 : : // Compute the normal velocities from left and right cells
849 : 0 : auto vnl = ul * fn[0] + vl * fn[1] + wl * fn[2];
850 : 0 : auto vnr = ur * fn[0] + vr * fn[1] + wr * fn[2];
851 : :
852 : : // The interface velocity is evaluated by adding the normal velocity which
853 : : // is taken from the Riemann solver and the tangential velocity which is
854 : : // evaluated as an average of the left and right cells
855 : 0 : auto urie = 0.5 * ((ul + ur) - fn[0] * (vnl + vnr)) + fl[ncomp+nmat] * fn[0];
856 : 0 : auto vrie = 0.5 * ((vl + vr) - fn[1] * (vnl + vnr)) + fl[ncomp+nmat] * fn[1];
857 : 0 : auto wrie = 0.5 * ((wl + wr) - fn[2] * (vnl + vnr)) + fl[ncomp+nmat] * fn[2];
858 : :
859 [ - - ]: 0 : vriem[el].push_back(urie);
860 [ - - ]: 0 : vriem[el].push_back(vrie);
861 [ - - ]: 0 : vriem[el].push_back(wrie);
862 : :
863 [ - - ]: 0 : if(e_right != -1)
864 : : {
865 [ - - ]: 0 : vriem[er].push_back(urie);
866 [ - - ]: 0 : vriem[er].push_back(vrie);
867 [ - - ]: 0 : vriem[er].push_back(wrie);
868 : : }
869 : 0 : }
870 : :
871 : : std::vector< std::array< tk::real, 3 > >
872 : 7374000 : fluxTerms(
873 : : std::size_t ncomp,
874 : : std::size_t nmat,
875 : : const std::vector< inciter::EOS >& /*mat_blk*/,
876 : : const std::vector< tk::real >& ugp )
877 : : // *****************************************************************************
878 : : // Compute the flux-function for the multimaterial PDEs
879 : : //! \param[in] ncomp Number of components in this PDE system
880 : : //! \param[in] nmat Number of materials in this PDE system
881 : : // //! \param[in] mat_blk EOS material block
882 : : //! \param[in] ugp State vector at the Gauss point at which flux is required
883 : : //! \return Flux vectors for all components in multi-material PDE system
884 : : // *****************************************************************************
885 : : {
886 : : using inciter::volfracIdx;
887 : : using inciter::densityIdx;
888 : : using inciter::momentumIdx;
889 : : using inciter::energyIdx;
890 : : using inciter::velocityIdx;
891 : : using inciter::pressureIdx;
892 : : using inciter::deformIdx;
893 : :
894 : : const auto& solidx =
895 : 7374000 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
896 : :
897 [ + - ]: 7374000 : std::vector< std::array< tk::real, 3 > > fl( ncomp, {{0, 0, 0}} );
898 : :
899 : 7374000 : tk::real rho(0.0);
900 [ + + ]: 22122000 : for (std::size_t k=0; k<nmat; ++k)
901 : 14748000 : rho += ugp[densityIdx(nmat, k)];
902 : :
903 : 7374000 : auto u = ugp[ncomp+velocityIdx(nmat,0)];
904 : 7374000 : auto v = ugp[ncomp+velocityIdx(nmat,1)];
905 : 7374000 : auto w = ugp[ncomp+velocityIdx(nmat,2)];
906 : :
907 [ + - ][ - + ]: 7374000 : if (inciter::haveSolid(nmat, solidx))
908 : : {
909 [ - - ]: 0 : std::vector< tk::real > al(nmat, 0.0);
910 : 0 : std::vector< std::array< std::array< tk::real, 3 >, 3 > > g, asig;
911 : : std::array< std::array< tk::real, 3 >, 3 >
912 : 0 : sig {{ {{0, 0, 0}}, {{0, 0, 0}}, {{0, 0, 0}} }};
913 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
914 : : {
915 : 0 : al[k] = ugp[volfracIdx(nmat, k)];
916 : : // inv deformation gradient and Cauchy stress tensors
917 [ - - ][ - - ]: 0 : g.push_back(inciter::getDeformGrad(nmat, k, ugp));
918 [ - - ][ - - ]: 0 : asig.push_back(inciter::getCauchyStress(nmat, k, ncomp, ugp));
919 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) asig[k][i][i] -= ugp[ncomp+pressureIdx(nmat,k)];
920 : :
921 [ - - ]: 0 : for (size_t i=0; i<3; ++i)
922 [ - - ]: 0 : for (size_t j=0; j<3; ++j)
923 : 0 : sig[i][j] += asig[k][i][j];
924 : : }
925 : :
926 : : // conservative part of momentum flux
927 : 0 : fl[momentumIdx(nmat, 0)][0] = ugp[momentumIdx(nmat, 0)] * u - sig[0][0];
928 : 0 : fl[momentumIdx(nmat, 1)][0] = ugp[momentumIdx(nmat, 1)] * u - sig[0][1];
929 : 0 : fl[momentumIdx(nmat, 2)][0] = ugp[momentumIdx(nmat, 2)] * u - sig[0][2];
930 : :
931 : 0 : fl[momentumIdx(nmat, 0)][1] = ugp[momentumIdx(nmat, 0)] * v - sig[1][0];
932 : 0 : fl[momentumIdx(nmat, 1)][1] = ugp[momentumIdx(nmat, 1)] * v - sig[1][1];
933 : 0 : fl[momentumIdx(nmat, 2)][1] = ugp[momentumIdx(nmat, 2)] * v - sig[1][2];
934 : :
935 : 0 : fl[momentumIdx(nmat, 0)][2] = ugp[momentumIdx(nmat, 0)] * w - sig[2][0];
936 : 0 : fl[momentumIdx(nmat, 1)][2] = ugp[momentumIdx(nmat, 1)] * w - sig[2][1];
937 : 0 : fl[momentumIdx(nmat, 2)][2] = ugp[momentumIdx(nmat, 2)] * w - sig[2][2];
938 : :
939 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
940 : : {
941 : : // conservative part of volume-fraction flux
942 : 0 : fl[volfracIdx(nmat, k)][0] = 0.0;
943 : 0 : fl[volfracIdx(nmat, k)][1] = 0.0;
944 : 0 : fl[volfracIdx(nmat, k)][2] = 0.0;
945 : :
946 : : // conservative part of material continuity flux
947 : 0 : fl[densityIdx(nmat, k)][0] = u * ugp[densityIdx(nmat, k)];
948 : 0 : fl[densityIdx(nmat, k)][1] = v * ugp[densityIdx(nmat, k)];
949 : 0 : fl[densityIdx(nmat, k)][2] = w * ugp[densityIdx(nmat, k)];
950 : :
951 : : // conservative part of material total-energy flux
952 : 0 : fl[energyIdx(nmat, k)][0] = u * ugp[energyIdx(nmat, k)]
953 : 0 : - u * asig[k][0][0] - v * asig[k][1][0] - w * asig[k][2][0];
954 : 0 : fl[energyIdx(nmat, k)][1] = v * ugp[energyIdx(nmat, k)]
955 : 0 : - u * asig[k][0][1] - v * asig[k][1][1] - w * asig[k][2][1];
956 : 0 : fl[energyIdx(nmat, k)][2] = w * ugp[energyIdx(nmat, k)]
957 : 0 : - u * asig[k][0][2] - v * asig[k][1][2] - w * asig[k][2][2];
958 : :
959 : : // conservative part of material inverse deformation gradient
960 : : // g_ij: \partial (g_il u_l) / \partial (x_j)
961 [ - - ]: 0 : if (solidx[k] > 0)
962 : : {
963 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
964 : : {
965 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
966 : : {
967 : 0 : fl[deformIdx(nmat,solidx[k],i,j)][j] =
968 : 0 : u*g[k][i][0] + v*g[k][i][1] + w*g[k][i][2];
969 : : }
970 : : // other components are zero
971 : : }
972 : : }
973 : : }
974 : : }
975 : : else
976 : : {
977 : 7374000 : tk::real p(0.0);
978 [ + - ]: 14748000 : std::vector< tk::real > apk( nmat, 0.0 );
979 [ + + ]: 22122000 : for (std::size_t k=0; k<nmat; ++k)
980 : : {
981 : 14748000 : apk[k] = ugp[ncomp+pressureIdx(nmat,k)];
982 : 14748000 : p += apk[k];
983 : : }
984 : :
985 : : // conservative part of momentum flux
986 : 7374000 : fl[momentumIdx(nmat, 0)][0] = ugp[momentumIdx(nmat, 0)] * u + p;
987 : 7374000 : fl[momentumIdx(nmat, 1)][0] = ugp[momentumIdx(nmat, 1)] * u;
988 : 7374000 : fl[momentumIdx(nmat, 2)][0] = ugp[momentumIdx(nmat, 2)] * u;
989 : :
990 : 7374000 : fl[momentumIdx(nmat, 0)][1] = ugp[momentumIdx(nmat, 0)] * v;
991 : 7374000 : fl[momentumIdx(nmat, 1)][1] = ugp[momentumIdx(nmat, 1)] * v + p;
992 : 7374000 : fl[momentumIdx(nmat, 2)][1] = ugp[momentumIdx(nmat, 2)] * v;
993 : :
994 : 7374000 : fl[momentumIdx(nmat, 0)][2] = ugp[momentumIdx(nmat, 0)] * w;
995 : 7374000 : fl[momentumIdx(nmat, 1)][2] = ugp[momentumIdx(nmat, 1)] * w;
996 : 7374000 : fl[momentumIdx(nmat, 2)][2] = ugp[momentumIdx(nmat, 2)] * w + p;
997 : :
998 [ + + ]: 22122000 : for (std::size_t k=0; k<nmat; ++k)
999 : : {
1000 : : // conservative part of volume-fraction flux
1001 : 14748000 : fl[volfracIdx(nmat, k)][0] = 0.0;
1002 : 14748000 : fl[volfracIdx(nmat, k)][1] = 0.0;
1003 : 14748000 : fl[volfracIdx(nmat, k)][2] = 0.0;
1004 : :
1005 : : // conservative part of material continuity flux
1006 : 14748000 : fl[densityIdx(nmat, k)][0] = u * ugp[densityIdx(nmat, k)];
1007 : 14748000 : fl[densityIdx(nmat, k)][1] = v * ugp[densityIdx(nmat, k)];
1008 : 14748000 : fl[densityIdx(nmat, k)][2] = w * ugp[densityIdx(nmat, k)];
1009 : :
1010 : : // conservative part of material total-energy flux
1011 : 14748000 : auto hmat = ugp[energyIdx(nmat, k)] + apk[k];
1012 : 14748000 : fl[energyIdx(nmat, k)][0] = u * hmat;
1013 : 14748000 : fl[energyIdx(nmat, k)][1] = v * hmat;
1014 : 14748000 : fl[energyIdx(nmat, k)][2] = w * hmat;
1015 : : }
1016 : : }
1017 : :
1018 : 7374000 : return fl;
1019 : : }
1020 : :
1021 : : }// tk::
|
