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 : : #ifdef HAS_MKL
20 : : #include <mkl_lapacke.h>
21 : : #else
22 : : #include <lapacke.h>
23 : : #endif
24 : :
25 : : #include "MultiMatTerms.hpp"
26 : : #include "Vector.hpp"
27 : : #include "Quadrature.hpp"
28 : : #include "EoS/EoS.hpp"
29 : : #include "MultiMat/MultiMatIndexing.hpp"
30 : : #include "Reconstruction.hpp"
31 : :
32 : : namespace tk {
33 : :
34 : : void
35 : 5250 : nonConservativeInt( [[maybe_unused]] ncomp_t system,
36 : : std::size_t nmat,
37 : : ncomp_t offset,
38 : : const std::size_t ndof,
39 : : const std::size_t rdof,
40 : : const std::size_t nelem,
41 : : const std::vector< std::size_t >& inpoel,
42 : : const UnsMesh::Coords& coord,
43 : : const Fields& geoElem,
44 : : const Fields& U,
45 : : const Fields& P,
46 : : const std::vector< std::vector< tk::real > >& riemannDeriv,
47 : : const std::vector< std::vector< tk::real > >& vriempoly,
48 : : const std::vector< std::size_t >& ndofel,
49 : : Fields& R,
50 : : int intsharp )
51 : : // *****************************************************************************
52 : : // Compute volume integrals for multi-material DG
53 : : //! \details This is called for multi-material DG, computing volume integrals of
54 : : //! terms in the volume fraction and energy equations, which do not exist in
55 : : //! the single-material flow formulation (for `CompFlow` DG). For further
56 : : //! details see Pelanti, M., & Shyue, K. M. (2019). A numerical model for
57 : : //! multiphase liquid–vapor–gas flows with interfaces and cavitation.
58 : : //! International Journal of Multiphase Flow, 113, 208-230.
59 : : //! \param[in] system Equation system index
60 : : //! \param[in] nmat Number of materials in this PDE system
61 : : //! \param[in] offset Offset this PDE system operates from
62 : : //! \param[in] ndof Maximum number of degrees of freedom
63 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
64 : : //! \param[in] nelem Total number of elements
65 : : //! \param[in] inpoel Element-node connectivity
66 : : //! \param[in] coord Array of nodal coordinates
67 : : //! \param[in] geoElem Element geometry array
68 : : //! \param[in] U Solution vector at recent time step
69 : : //! \param[in] P Vector of primitive quantities at recent time step
70 : : //! \param[in] riemannDeriv Derivatives of partial-pressures and velocities
71 : : //! computed from the Riemann solver for use in the non-conservative terms
72 : : //! \param[in] vriempoly Vector of Riemann velocity polynomial
73 : : //! \param[in] ndofel Vector of local number of degrees of freedome
74 : : //! \param[in,out] R Right-hand side vector added to
75 : : //! \param[in] intsharp Interface reconstruction indicator
76 : : // *****************************************************************************
77 : : {
78 : : using inciter::volfracIdx;
79 : : using inciter::densityIdx;
80 : : using inciter::momentumIdx;
81 : : using inciter::energyIdx;
82 : : using inciter::velocityIdx;
83 : :
84 : : const auto& cx = coord[0];
85 : : const auto& cy = coord[1];
86 : : const auto& cz = coord[2];
87 : :
88 : 5250 : auto ncomp = U.nprop()/rdof;
89 : 5250 : auto nprim = P.nprop()/rdof;
90 : :
91 : : // compute volume integrals
92 [ + + ]: 2235150 : for (std::size_t e=0; e<nelem; ++e)
93 : : {
94 : 2229900 : auto ng = tk::NGvol(ndofel[e]);
95 : :
96 : : // arrays for quadrature points
97 : : std::array< std::vector< real >, 3 > coordgp;
98 : : std::vector< real > wgp;
99 : :
100 [ + - ]: 2229900 : coordgp[0].resize( ng );
101 [ + - ]: 2229900 : coordgp[1].resize( ng );
102 [ + - ]: 2229900 : coordgp[2].resize( ng );
103 [ + - ]: 2229900 : wgp.resize( ng );
104 : :
105 [ + - ]: 2229900 : GaussQuadratureTet( ng, coordgp, wgp );
106 : :
107 : : // Extract the element coordinates
108 : : std::array< std::array< real, 3>, 4 > coordel {{
109 : 2229900 : {{ cx[ inpoel[4*e ] ], cy[ inpoel[4*e ] ], cz[ inpoel[4*e ] ] }},
110 : 2229900 : {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
111 : 2229900 : {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
112 : 2229900 : {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
113 : 2229900 : }};
114 : :
115 : : auto jacInv =
116 : 2229900 : inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );
117 : :
118 : : // Compute the derivatives of basis function for DG(P1)
119 : : std::array< std::vector<tk::real>, 3 > dBdx;
120 [ + + ]: 2229900 : if (ndofel[e] > 1)
121 [ + - ]: 909600 : dBdx = eval_dBdx_p1( ndofel[e], jacInv );
122 : :
123 : : // Gaussian quadrature
124 [ + + ]: 6279000 : for (std::size_t igp=0; igp<ng; ++igp)
125 : : {
126 [ - + ]: 4049100 : if (ndofel[e] > 4)
127 [ - - ]: 0 : eval_dBdx_p2( igp, coordgp, jacInv, dBdx );
128 : :
129 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
130 : : // basis functions and solutions by default. This implies that P0P1 is
131 : : // unsupported in the p-adaptive DG (PDG).
132 : : std::size_t dof_el;
133 [ + + ]: 4049100 : if (rdof > ndof)
134 : : {
135 : : dof_el = rdof;
136 : : }
137 : : else
138 : : {
139 : 3594300 : dof_el = ndofel[e];
140 : : }
141 : :
142 : : // Compute the basis function
143 : : auto B =
144 [ + - ]: 4049100 : eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
145 : :
146 [ + - ]: 4049100 : auto wt = wgp[igp] * geoElem(e, 0, 0);
147 : :
148 : : auto state = evalPolynomialSol(system, offset, intsharp, ncomp, nprim,
149 : : rdof, nmat, e, dof_el, inpoel, coord, geoElem,
150 [ + - ][ - - ]: 4049100 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, U, P);
151 : :
152 : : // get bulk properties
153 : : tk::real rhob(0.0);
154 [ + + ]: 13240200 : for (std::size_t k=0; k<nmat; ++k)
155 : 9191100 : rhob += state[densityIdx(nmat, k)];
156 : :
157 : : // get the velocity vector
158 [ + - ]: 4049100 : std::array< tk::real, 3 > vel{{ state[ncomp+velocityIdx(nmat, 0)],
159 : 4049100 : state[ncomp+velocityIdx(nmat, 1)],
160 [ + - ]: 4049100 : state[ncomp+velocityIdx(nmat, 2)] }};
161 : :
162 [ + - ][ - - ]: 4049100 : std::vector< tk::real > ymat(nmat, 0.0);
163 : 4049100 : std::array< tk::real, 3 > dap{{0.0, 0.0, 0.0}};
164 [ + + ]: 13240200 : for (std::size_t k=0; k<nmat; ++k)
165 : : {
166 : 9191100 : ymat[k] = state[densityIdx(nmat, k)]/rhob;
167 : :
168 [ + + ]: 36764400 : for (std::size_t idir=0; idir<3; ++idir)
169 : 27573300 : dap[idir] += riemannDeriv[3*k+idir][e];
170 : : }
171 : :
172 : : // Evaluate the velocity used for the multi-material term integration for
173 : : // volume fraction equation
174 [ + - ][ - - ]: 4049100 : std::vector< tk::real> vriem(3, 0.0);
175 [ + + ]: 4049100 : if(ndofel[e] > 1)
176 : : {
177 [ + - ]: 2274000 : auto gp = eval_gp( igp, coordel, coordgp );
178 [ + + ]: 9096000 : for(std::size_t idir = 0; idir < 3; idir++)
179 : : {
180 : 6822000 : auto mark = idir * 4;
181 : 6822000 : vriem[idir] = vriempoly[e][mark];
182 [ + + ]: 27288000 : for(std::size_t k = 1; k < 4; k++)
183 : 20466000 : vriem[idir] += vriempoly[e][mark+k] * gp[k-1];
184 : : }
185 : : }
186 : :
187 : : // compute non-conservative terms
188 : : std::vector< std::vector< tk::real > > ncf
189 [ + - ][ + - ]: 12147300 : (ncomp, std::vector<tk::real>(ndof,0.0));
[ - - ]
190 : :
191 [ + + ]: 16196400 : for (std::size_t idir=0; idir<3; ++idir)
192 [ + + ]: 44760600 : for(std::size_t idof=0; idof<ndof; ++idof)
193 : 32613300 : ncf[momentumIdx(nmat, idir)][idof] = 0.0;
194 : :
195 [ + + ]: 13240200 : for (std::size_t k=0; k<nmat; ++k)
196 : : {
197 : : // evaluate non-conservative term for energy equation
198 [ + + ]: 32026200 : for(std::size_t idof=0; idof<ndof; ++idof)
199 : : {
200 : 22835100 : ncf[densityIdx(nmat, k)][idof] = 0.0;
201 [ + + ]: 91340400 : for (std::size_t idir=0; idir<3; ++idir)
202 : 68505300 : ncf[energyIdx(nmat, k)][idof] -= vel[idir] * ( ymat[k]*dap[idir]
203 : 68505300 : - riemannDeriv[3*k+idir][e] );
204 : : }
205 : :
206 : : // evaluate non-conservative term for volume fraction equation
207 [ + + ]: 32026200 : for(std::size_t idof=0; idof<ndof; ++idof)
208 : 22835100 : ncf[volfracIdx(nmat, k)][idof] = state[volfracIdx(nmat, k)]
209 : 22835100 : * riemannDeriv[3*nmat][e] * B[idof];
210 : :
211 : : // evaluate the non-conservative term for volume fraction equation in
212 : : // high order discretization. The following code compuetes the
213 : : // non-conservative term:
214 : : // alpha * (d(u*B)/dx) = alpha * (u*dBdx + B*dudx)
215 [ + + ]: 9191100 : if (ndof > 1 && intsharp == 0)
216 [ + + ]: 9096000 : for(std::size_t idof=1; idof<ndof; ++idof)
217 [ + + ]: 27288000 : for(std::size_t idir=0; idir<3; ++idir)
218 : 20466000 : ncf[volfracIdx(nmat, k)][idof] += state[volfracIdx(nmat, k)] *
219 : 20466000 : vriem[idir] * dBdx[idir][idof];
220 : : }
221 : :
222 [ + - ]: 4049100 : updateRhsNonCons( ncomp, offset, nmat, ndof, ndofel[e], wt, e, B, dBdx, ncf, R );
223 : : }
224 : : }
225 : 5250 : }
226 : :
227 : : void
228 : 4049100 : updateRhsNonCons(
229 : : ncomp_t ncomp,
230 : : ncomp_t offset,
231 : : const std::size_t nmat,
232 : : const std::size_t ndof,
233 : : [[maybe_unused]] const std::size_t ndof_el,
234 : : const tk::real wt,
235 : : const std::size_t e,
236 : : const std::vector<tk::real>& B,
237 : : [[maybe_unused]] const std::array< std::vector<tk::real>, 3 >& dBdx,
238 : : const std::vector< std::vector< tk::real > >& ncf,
239 : : Fields& R )
240 : : // *****************************************************************************
241 : : // Update the rhs by adding the non-conservative term integrals
242 : : //! \param[in] ncomp Number of scalar components in this PDE system
243 : : //! \param[in] offset Offset this PDE system operates from
244 : : //! \param[in] nmat Number of materials
245 : : //! \param[in] ndof Maximum number of degrees of freedom
246 : : //! \param[in] ndof_el Number of degrees of freedom for local element
247 : : //! \param[in] wt Weight of gauss quadrature point
248 : : //! \param[in] e Element index
249 : : //! \param[in] B Basis function evaluated at local quadrature point
250 : : //! \param[in] dBdx Vector of basis function derivatives
251 : : //! \param[in] ncf Vector of non-conservative terms
252 : : //! \param[in,out] R Right-hand side vector computed
253 : : // *****************************************************************************
254 : : {
255 : : using inciter::volfracIdx;
256 : : using inciter::energyIdx;
257 : :
258 : : //Assert( dBdx[0].size() == ndof_el,
259 : : // "Size mismatch for basis function derivatives" );
260 : : //Assert( dBdx[1].size() == ndof_el,
261 : : // "Size mismatch for basis function derivatives" );
262 : : //Assert( dBdx[2].size() == ndof_el,
263 : : // "Size mismatch for basis function derivatives" );
264 : : //Assert( ncf.size() == ncomp,
265 : : // "Size mismatch for non-conservative term" );
266 : : Assert( ncf.size() == ncomp, "Size mismatch for non-conservative term" );
267 : :
268 [ + + ]: 43769700 : for (ncomp_t c=0; c<ncomp; ++c)
269 : : {
270 : 39720600 : auto mark = c*ndof;
271 : 39720600 : R(e, mark, offset) += wt * ncf[c][0];
272 : : }
273 : :
274 [ + + ]: 4049100 : if( ndof_el > 1)
275 : : {
276 : : // Update rhs with distributions from volume fraction equations
277 [ + + ]: 6822000 : for (std::size_t k=volfracIdx(nmat,0); k<volfracIdx(nmat,nmat); ++k)
278 : : {
279 : 4548000 : auto mark = k*ndof;
280 [ + + ]: 18192000 : for(std::size_t idof = 1; idof < ndof; idof++)
281 : 13644000 : R(e, mark+idof, offset) += wt * ncf[k][idof];
282 : : }
283 : :
284 : : // Update rhs with distributions from the rest of the equatons
285 [ + + ]: 6822000 : for (std::size_t c=energyIdx(nmat,0); c<energyIdx(nmat,nmat); ++c)
286 : : {
287 : 4548000 : auto mark = c*ndof;
288 [ + + ]: 18192000 : for(std::size_t idof = 1; idof < ndof; idof++)
289 : 13644000 : R(e, mark+idof, offset) += wt * ncf[c][idof] * B[idof];
290 : : }
291 : : }
292 : 4049100 : }
293 : :
294 : : void
295 : 450 : pressureRelaxationInt( ncomp_t system,
296 : : std::size_t nmat,
297 : : ncomp_t offset,
298 : : const std::size_t ndof,
299 : : const std::size_t rdof,
300 : : const std::size_t nelem,
301 : : const std::vector< std::size_t >& inpoel,
302 : : const UnsMesh::Coords& coord,
303 : : const Fields& geoElem,
304 : : const Fields& U,
305 : : const Fields& P,
306 : : const std::vector< std::size_t >& ndofel,
307 : : const tk::real ct,
308 : : Fields& R,
309 : : int intsharp )
310 : : // *****************************************************************************
311 : : // Compute volume integrals of pressure relaxation terms in multi-material DG
312 : : //! \details This is called for multi-material DG to compute volume integrals of
313 : : //! finite pressure relaxation terms in the volume fraction and energy
314 : : //! equations, which do not exist in the single-material flow formulation (for
315 : : //! `CompFlow` DG). For further details see Dobrev, V. A., Kolev, T. V.,
316 : : //! Rieben, R. N., & Tomov, V. Z. (2016). Multi‐material closure model for
317 : : //! high‐order finite element Lagrangian hydrodynamics. International Journal
318 : : //! for Numerical Methods in Fluids, 82(10), 689-706.
319 : : //! \param[in] system Equation system index
320 : : //! \param[in] nmat Number of materials in this PDE system
321 : : //! \param[in] offset Offset this PDE system operates from
322 : : //! \param[in] ndof Maximum number of degrees of freedom
323 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
324 : : //! \param[in] nelem Total number of elements
325 : : //! \param[in] inpoel Element-node connectivity
326 : : //! \param[in] coord Array of nodal coordinates
327 : : //! \param[in] geoElem Element geometry array
328 : : //! \param[in] U Solution vector at recent time step
329 : : //! \param[in] P Vector of primitive quantities at recent time step
330 : : //! \param[in] ndofel Vector of local number of degrees of freedome
331 : : //! \param[in] ct Pressure relaxation time-scale for this system
332 : : //! \param[in,out] R Right-hand side vector added to
333 : : //! \param[in] intsharp Interface reconstruction indicator
334 : : // *****************************************************************************
335 : : {
336 : : using inciter::volfracIdx;
337 : : using inciter::densityIdx;
338 : : using inciter::momentumIdx;
339 : : using inciter::energyIdx;
340 : : using inciter::pressureIdx;
341 : : using inciter::velocityIdx;
342 : :
343 : 450 : auto ncomp = U.nprop()/rdof;
344 : 450 : auto nprim = P.nprop()/rdof;
345 : :
346 : : // compute volume integrals
347 [ + + ]: 341550 : for (std::size_t e=0; e<nelem; ++e)
348 : : {
349 : 341100 : auto dx = geoElem(e,4,0)/2.0;
350 : 341100 : auto ng = NGvol(ndofel[e]);
351 : :
352 : : // arrays for quadrature points
353 : : std::array< std::vector< real >, 3 > coordgp;
354 : : std::vector< real > wgp;
355 : :
356 [ + - ]: 341100 : coordgp[0].resize( ng );
357 [ + - ]: 341100 : coordgp[1].resize( ng );
358 [ + - ]: 341100 : coordgp[2].resize( ng );
359 [ + - ]: 341100 : wgp.resize( ng );
360 : :
361 [ + - ]: 341100 : GaussQuadratureTet( ng, coordgp, wgp );
362 : :
363 : : // Compute the derivatives of basis function for DG(P1)
364 : : std::array< std::vector<real>, 3 > dBdx;
365 : :
366 : : // Gaussian quadrature
367 [ + + ]: 1591800 : for (std::size_t igp=0; igp<ng; ++igp)
368 : : {
369 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
370 : : // basis functions and solutions by default. This implies that P0P1 is
371 : : // unsupported in the p-adaptive DG (PDG).
372 : : std::size_t dof_el;
373 [ + + ]: 1250700 : if (rdof > ndof)
374 : : {
375 : : dof_el = rdof;
376 : : }
377 : : else
378 : : {
379 : 1137000 : dof_el = ndofel[e];
380 : : }
381 : :
382 : : // Compute the basis function
383 : : auto B =
384 [ + - ]: 1250700 : eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
385 : :
386 [ + - ]: 1250700 : auto wt = wgp[igp] * geoElem(e, 0, 0);
387 : :
388 : : auto state = evalPolynomialSol(system, offset, intsharp, ncomp, nprim,
389 : : rdof, nmat, e, dof_el, inpoel, coord, geoElem,
390 [ + - ][ - - ]: 1250700 : {{coordgp[0][igp], coordgp[1][igp], coordgp[2][igp]}}, B, U, P);
391 : :
392 : : // get bulk properties
393 : : real rhob(0.0);
394 : : for (std::size_t k=0; k<nmat; ++k)
395 : : rhob += state[densityIdx(nmat, k)];
396 : :
397 : : // get pressures and bulk modulii
398 : 1250700 : real pb(0.0), nume(0.0), deno(0.0), trelax(0.0);
399 [ + - ][ + - ]: 1250700 : std::vector< real > apmat(nmat, 0.0), kmat(nmat, 0.0);
[ - - ][ - - ]
400 [ + + ]: 3752100 : for (std::size_t k=0; k<nmat; ++k)
401 : : {
402 [ + - ]: 2501400 : real arhomat = state[densityIdx(nmat, k)];
403 : 2501400 : real alphamat = state[volfracIdx(nmat, k)];
404 : 2501400 : apmat[k] = state[ncomp+pressureIdx(nmat, k)];
405 [ + - ]: 2501400 : real amat = inciter::eos_soundspeed< tag::multimat >( system, arhomat,
406 : : apmat[k], alphamat, k );
407 [ + + ]: 2501400 : kmat[k] = arhomat * amat * amat;
408 : 2501400 : pb += apmat[k];
409 : :
410 : : // relaxation parameters
411 [ + + ]: 2501400 : trelax = std::max(trelax, ct*dx/amat);
412 : 2501400 : nume += alphamat * apmat[k] / kmat[k];
413 : 2501400 : deno += alphamat * alphamat / kmat[k];
414 : : }
415 : 1250700 : auto p_relax = nume/deno;
416 : :
417 : : // compute pressure relaxation terms
418 [ + - ][ - - ]: 1250700 : std::vector< real > s_prelax(ncomp, 0.0);
419 [ + + ]: 3752100 : for (std::size_t k=0; k<nmat; ++k)
420 : : {
421 : 2501400 : auto s_alpha = (apmat[k]-p_relax*state[volfracIdx(nmat, k)])
422 : 2501400 : * (state[volfracIdx(nmat, k)]/kmat[k]) / trelax;
423 : 2501400 : s_prelax[volfracIdx(nmat, k)] = s_alpha;
424 : 2501400 : s_prelax[energyIdx(nmat, k)] = - pb*s_alpha;
425 : : }
426 : :
427 [ + - ]: 1250700 : updateRhsPre( ncomp, offset, ndof, dof_el, wt, e, B, s_prelax, R );
428 : : }
429 : : }
430 : 450 : }
431 : :
432 : : void
433 : 1250700 : updateRhsPre(
434 : : ncomp_t ncomp,
435 : : ncomp_t offset,
436 : : const std::size_t ndof,
437 : : [[maybe_unused]] const std::size_t ndof_el,
438 : : const tk::real wt,
439 : : const std::size_t e,
440 : : const std::vector< tk::real >& B,
441 : : std::vector< tk::real >& ncf,
442 : : Fields& R )
443 : : // *****************************************************************************
444 : : // Update the rhs by adding the pressure relaxation integrals
445 : : //! \param[in] ncomp Number of scalar components in this PDE system
446 : : //! \param[in] offset Offset this PDE system operates from
447 : : //! \param[in] ndof Maximum number of degrees of freedom
448 : : //! \param[in] ndof_el Number of degrees of freedom for local element
449 : : //! \param[in] wt Weight of gauss quadrature point
450 : : //! \param[in] e Element index
451 : : //! \param[in] B Basis function evaluated at local quadrature point
452 : : //! \param[in] ncf Vector of non-conservative terms
453 : : //! \param[in,out] R Right-hand side vector computed
454 : : // *****************************************************************************
455 : : {
456 : : //Assert( dBdx[0].size() == ndof_el,
457 : : // "Size mismatch for basis function derivatives" );
458 : : //Assert( dBdx[1].size() == ndof_el,
459 : : // "Size mismatch for basis function derivatives" );
460 : : //Assert( dBdx[2].size() == ndof_el,
461 : : // "Size mismatch for basis function derivatives" );
462 : : //Assert( ncf.size() == ncomp,
463 : : // "Size mismatch for non-conservative term" );
464 : : Assert( ncf.size() == ncomp, "Size mismatch for non-conservative term" );
465 : :
466 [ + + ]: 12507000 : for (ncomp_t c=0; c<ncomp; ++c)
467 : : {
468 : 11256300 : auto mark = c*ndof;
469 [ + + ]: 53211600 : for(std::size_t idof = 0; idof < ndof; idof++)
470 : 41955300 : R(e, mark+idof, offset) += wt * ncf[c] * B[idof];
471 : : }
472 : 1250700 : }
473 : :
474 : : std::vector< std::vector< tk::real > >
475 : 375 : solvevriem( std::size_t nelem,
476 : : const std::vector< std::vector< tk::real > >& vriem,
477 : : const std::vector< std::vector< tk::real > >& riemannLoc )
478 : : // *****************************************************************************
479 : : // Solve the reconstruct velocity used for volume fraction equation by
480 : : // Least square method
481 : : //! \param[in] nelem Numer of elements
482 : : //! \param[in,out] vriem Vector of the riemann velocity
483 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
484 : : //! is available
485 : : //! \return Vector of Riemann velocity polynomial solution
486 : : // *****************************************************************************
487 : : {
488 : : std::vector< std::vector< tk::real > >
489 [ + - ][ + - ]: 750 : vriempoly( nelem, std::vector<tk::real>(12,0.0) );
490 : :
491 [ + + ]: 227775 : for (std::size_t e=0; e<nelem; ++e)
492 : : {
493 : : // Use the normal method to construct the linear system A^T * A * x = u
494 [ + - ]: 227400 : auto numgp = riemannLoc[e].size()/3;
495 : : std::vector< std::vector< tk::real > > A(numgp,
496 [ + - ][ + - ]: 682200 : std::vector<tk::real>(4, 1.0));
497 : :
498 [ + + ]: 2956200 : for(std::size_t k = 0; k < numgp; k++)
499 : : {
500 : 2728800 : auto mark = k * 3;
501 : 2728800 : A[k][1] = riemannLoc[e][mark];
502 : 2728800 : A[k][2] = riemannLoc[e][mark+1];
503 : 2728800 : A[k][3] = riemannLoc[e][mark+2];
504 : : }
505 : :
506 [ + + ]: 909600 : for(std::size_t idir = 0; idir < 3; idir++)
507 : : {
508 : : double AA_T[4*4], u[4];
509 : :
510 [ + + ]: 3411000 : for(std::size_t i = 0; i < 4; i++)
511 [ + + ]: 13644000 : for(std::size_t j = 0; j < 4; j++)
512 : : {
513 : 10915200 : auto id = 4 * i + j;
514 : 10915200 : AA_T[id] = 0;
515 [ + + ]: 141897600 : for(std::size_t k = 0; k < numgp; k++)
516 : 130982400 : AA_T[id] += A[k][i] * A[k][j];
517 : : }
518 : :
519 [ + - ]: 682200 : std::vector<tk::real> vel(numgp, 1.0);
520 [ + + ]: 8868600 : for(std::size_t k = 0; k < numgp; k++)
521 : : {
522 : 8186400 : auto mark = k * 3 + idir;
523 : 8186400 : vel[k] = vriem[e][mark];
524 : : }
525 [ + + ]: 3411000 : for(std::size_t k = 0; k < 4; k++)
526 : : {
527 : 2728800 : u[k] = 0;
528 [ + + ]: 35474400 : for(std::size_t i = 0; i < numgp; i++)
529 : 32745600 : u[k] += A[i][k] * vel[i];
530 : : }
531 : :
532 : : lapack_int IPIV[4];
533 : : #ifndef NDEBUG
534 : : lapack_int info =
535 : : #endif
536 [ + - ]: 682200 : LAPACKE_dsysv( LAPACK_ROW_MAJOR, 'U', 4, 1, AA_T, 4, IPIV, u, 1 );
537 : : Assert( info == 0, "Error in linear system solver" );
538 : :
539 : 682200 : auto idirmark = idir * 4;
540 [ + + ]: 3411000 : for(std::size_t k = 0; k < 4; k++)
541 : 2728800 : vriempoly[e][idirmark+k] = u[k];
542 : : }
543 : : }
544 : 375 : return vriempoly;
545 : : }
546 : :
547 : 3152700 : void evaluRiemann( ncomp_t ncomp,
548 : : const int e_left,
549 : : const int e_right,
550 : : const std::size_t nmat,
551 : : const std::vector< tk::real >& fl,
552 : : const std::array< tk::real, 3 >& fn,
553 : : const std::array< tk::real, 3 >& gp,
554 : : const std::array< std::vector< tk::real >, 2 >& state,
555 : : std::vector< std::vector< tk::real > >& vriem,
556 : : std::vector< std::vector< tk::real > >& riemannLoc )
557 : : // *****************************************************************************
558 : : // Compute the riemann velocity at the interface
559 : : //! \param[in] ncomp Number of scalar components in this PDE system
560 : : //! \param[in] e_left Index for the left element
561 : : //! \param[in] e_right Index for the right element
562 : : //! \param[in] nmat Number of materials in this PDE system
563 : : //! \param[in] fn Face/Surface normal
564 : : //! \param[in] gp Gauss points coordinates
565 : : //! \param[in] fl Surface flux
566 : : //! \param[in] state Vector of state variables for left and right side
567 : : //! \param[in,out] vriem Vector of the riemann velocity
568 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
569 : : //! is available
570 : : // *****************************************************************************
571 : : {
572 : : using inciter::densityIdx;
573 : : using inciter::momentumIdx;
574 : :
575 : : std::size_t el(0), er(0);
576 : 3152700 : el = static_cast< std::size_t >(e_left);
577 [ + + ]: 3152700 : if(e_right != -1)
578 : : er = static_cast< std::size_t >(e_right);
579 : :
580 : 3152700 : riemannLoc[el].push_back( gp[0] );
581 : 3152700 : riemannLoc[el].push_back( gp[1] );
582 : 3152700 : riemannLoc[el].push_back( gp[2] );
583 : :
584 [ + + ]: 3152700 : if(e_right != -1)
585 : : {
586 : 2436300 : riemannLoc[er].push_back( gp[0] );
587 : 2436300 : riemannLoc[er].push_back( gp[1] );
588 : 2436300 : riemannLoc[er].push_back( gp[2] );
589 : : }
590 : :
591 : : tk::real rhobl(0.0), rhobr(0.0);
592 [ + + ]: 9458100 : for (std::size_t k=0; k<nmat; ++k)
593 : : {
594 : 6305400 : rhobl += state[0][densityIdx(nmat, k)];
595 : 6305400 : rhobr += state[1][densityIdx(nmat, k)];
596 : : }
597 : :
598 : 3152700 : auto ul = state[0][momentumIdx(nmat, 0)] / rhobl;
599 : 3152700 : auto vl = state[0][momentumIdx(nmat, 1)] / rhobl;
600 : 3152700 : auto wl = state[0][momentumIdx(nmat, 2)] / rhobl;
601 : :
602 : 3152700 : auto ur = state[1][momentumIdx(nmat, 0)] / rhobr;
603 : 3152700 : auto vr = state[1][momentumIdx(nmat, 1)] / rhobr;
604 : 3152700 : auto wr = state[1][momentumIdx(nmat, 2)] / rhobr;
605 : :
606 : : // Compute the normal velocities from left and right cells
607 : 3152700 : auto vnl = ul * fn[0] + vl * fn[1] + wl * fn[2];
608 : 3152700 : auto vnr = ur * fn[0] + vr * fn[1] + wr * fn[2];
609 : :
610 : : // The interface velocity is evaluated by adding the normal velocity which
611 : : // is taken from the Riemann solver and the tangential velocity which is
612 : : // evaluated as an average of the left and right cells
613 : 3152700 : auto urie = 0.5 * ((ul + ur) - fn[0] * (vnl + vnr)) + fl[ncomp+nmat] * fn[0];
614 : 3152700 : auto vrie = 0.5 * ((vl + vr) - fn[1] * (vnl + vnr)) + fl[ncomp+nmat] * fn[1];
615 : 3152700 : auto wrie = 0.5 * ((wl + wr) - fn[2] * (vnl + vnr)) + fl[ncomp+nmat] * fn[2];
616 : :
617 : 3152700 : vriem[el].push_back(urie);
618 : 3152700 : vriem[el].push_back(vrie);
619 : 3152700 : vriem[el].push_back(wrie);
620 : :
621 [ + + ]: 3152700 : if(e_right != -1)
622 : : {
623 : 2436300 : vriem[er].push_back(urie);
624 : 2436300 : vriem[er].push_back(vrie);
625 : 2436300 : vriem[er].push_back(wrie);
626 : : }
627 : 3152700 : }
628 : :
629 : : }// tk::
|
