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