Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/Surface.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 internal surface integrals of a system
9 : : of PDEs in DG methods
10 : : \details This file contains functionality for computing internal surface
11 : : integrals of a system of PDEs used in discontinuous Galerkin methods for
12 : : various orders of numerical representation.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include <array>
17 : :
18 : : #include "Surface.hpp"
19 : : #include "Vector.hpp"
20 : : #include "Quadrature.hpp"
21 : : #include "Reconstruction.hpp"
22 : : #include "Integrate/SolidTerms.hpp"
23 : : #include "Inciter/InputDeck/InputDeck.hpp"
24 : : #include "MultiMat/MiscMultiMatFns.hpp"
25 : : #include "EoS/GetMatProp.hpp"
26 : :
27 : : namespace inciter {
28 : : extern ctr::InputDeck g_inputdeck;
29 : : }
30 : :
31 : : namespace tk {
32 : :
33 : : void
34 : 65205 : surfInt( const bool pref,
35 : : std::size_t nmat,
36 : : const std::vector< inciter::EOS >& mat_blk,
37 : : real t,
38 : : const std::size_t ndof,
39 : : const std::size_t rdof,
40 : : const std::vector< std::size_t >& inpoel,
41 : : const std::vector< std::size_t >& /*solidx*/,
42 : : const UnsMesh::Coords& coord,
43 : : const inciter::FaceData& fd,
44 : : const Fields& geoFace,
45 : : const Fields& geoElem,
46 : : const RiemannFluxFn& flux,
47 : : const VelFn& vel,
48 : : const Fields& U,
49 : : const Fields& P,
50 : : const std::vector< std::size_t >& ndofel,
51 : : const tk::real /*dt*/,
52 : : Fields& R,
53 : : std::vector< std::vector< tk::real > >& riemannDeriv,
54 : : int intsharp )
55 : : // *****************************************************************************
56 : : // Compute internal surface flux integrals
57 : : //! \param[in] pref Indicator for p-adaptive algorithm
58 : : //! \param[in] nmat Number of materials in this PDE system
59 : : //! \param[in] mat_blk EOS material block
60 : : //! \param[in] t Physical time
61 : : //! \param[in] ndof Maximum number of degrees of freedom
62 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
63 : : //! \param[in] inpoel Element-node connectivity
64 : : // //! \param[in] solidx Material index indicator
65 : : //! \param[in] coord Array of nodal coordinates
66 : : //! \param[in] fd Face connectivity and boundary conditions object
67 : : //! \param[in] geoFace Face geometry array
68 : : //! \param[in] geoElem Element geometry array
69 : : //! \param[in] flux Riemann flux function to use
70 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
71 : : //! \param[in] U Solution vector at recent time step
72 : : //! \param[in] P Vector of primitives at recent time step
73 : : //! \param[in] ndofel Vector of local number of degrees of freedom
74 : : // //! \param[in] dt Delta time
75 : : //! \param[in,out] R Right-hand side vector computed
76 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
77 : : //! computed from the Riemann solver for use in the non-conservative terms.
78 : : //! These derivatives are used only for multi-material hydro and unused for
79 : : //! single-material compflow and linear transport.
80 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
81 : : //! default 0, so that it is unused for single-material and transport.
82 : : // *****************************************************************************
83 : : {
84 : : const auto& esuf = fd.Esuf();
85 : : const auto& inpofa = fd.Inpofa();
86 : :
87 : : const auto& cx = coord[0];
88 : : const auto& cy = coord[1];
89 : : const auto& cz = coord[2];
90 : :
91 : 65205 : auto ncomp = U.nprop()/rdof;
92 : 65205 : auto nprim = P.nprop()/rdof;
93 : :
94 : : //// Determine if we have solids in our problem
95 : : //bool haveSolid = inciter::haveSolid(nmat, solidx);
96 : :
97 : : //Assert( (nmat==1 ? riemannDeriv.empty() : true), "Non-empty Riemann "
98 : : // "derivative vector for single material compflow" );
99 : :
100 : : // compute internal surface flux integrals
101 [ + + ]: 29323230 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
102 : : {
103 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
104 : : "as -1" );
105 : :
106 : 29258025 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
107 : 29258025 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
108 : :
109 : 29258025 : auto ng_l = tk::NGfa(ndofel[el]);
110 [ + + ]: 29258025 : auto ng_r = tk::NGfa(ndofel[er]);
111 : :
112 : : // When the number of gauss points for the left and right element are
113 : : // different, choose the larger ng
114 : 29258025 : auto ng = std::max( ng_l, ng_r );
115 : :
116 : : // arrays for quadrature points
117 : : std::array< std::vector< real >, 2 > coordgp;
118 : : std::vector< real > wgp;
119 : :
120 [ + - ]: 29258025 : coordgp[0].resize( ng );
121 [ + - ]: 29258025 : coordgp[1].resize( ng );
122 [ + - ]: 29258025 : wgp.resize( ng );
123 : :
124 : : // get quadrature point weights and coordinates for triangle
125 [ + - ]: 29258025 : GaussQuadratureTri( ng, coordgp, wgp );
126 : :
127 : : // Extract the element coordinates
128 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
129 : 29258025 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
130 : 29258025 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
131 : 29258025 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
132 : 29258025 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
133 : :
134 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
135 : 29258025 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
136 : 29258025 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
137 : 29258025 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
138 : 29258025 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
139 : :
140 : : // Compute the determinant of Jacobian matrix
141 : : auto detT_l =
142 : 29258025 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
143 : : auto detT_r =
144 : 29258025 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
145 : :
146 : : // Extract the face coordinates
147 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
148 : 29258025 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
149 : 29258025 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
150 : 29258025 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
151 : :
152 : : std::array< real, 3 >
153 : 29258025 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
154 : :
155 : : // Gaussian quadrature
156 [ + + ]: 85754727 : for (std::size_t igp=0; igp<ng; ++igp)
157 : : {
158 : : // Compute the coordinates of quadrature point at physical domain
159 [ + - ]: 56496702 : auto gp = eval_gp( igp, coordfa, coordgp );
160 : :
161 : : // In order to determine the high-order solution from the left and right
162 : : // elements at the surface quadrature points, the basis functions from
163 : : // the left and right elements are needed. For this, a transformation to
164 : : // the reference coordinates is necessary, since the basis functions are
165 : : // defined on the reference tetrahedron only.
166 : : // The transformation relations are shown below:
167 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
168 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
169 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
170 : :
171 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
172 : : // basis functions and solutions by default. This implies that P0P1 is
173 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
174 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
175 : : // element for pDG.
176 : : std::size_t dof_el, dof_er;
177 [ + + ]: 56496702 : if (rdof > ndof)
178 : : {
179 : : dof_el = rdof;
180 : : dof_er = rdof;
181 : : }
182 : : else
183 : : {
184 : 55799652 : dof_el = ndofel[el];
185 : 55799652 : dof_er = ndofel[er];
186 : : }
187 : :
188 : : // For multi-material p-adaptive simulation, when dofel = 1, p0p1 is
189 : : // applied and ndof for solution evaluation should be 4
190 [ - + ]: 56496702 : if(ncomp > 5 && pref) {
191 [ - - ]: 0 : if(dof_el == 1)
192 : : dof_el = 4;
193 [ - - ]: 0 : if(dof_er == 1)
194 : : dof_er = 4;
195 : : }
196 : :
197 : : std::array< tk::real, 3> ref_gp_l{
198 : 56496702 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
199 : 56496702 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
200 : 56496702 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
201 : : std::array< tk::real, 3> ref_gp_r{
202 : 56496702 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
203 : 56496702 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
204 : 56496702 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
205 : :
206 : : //Compute the basis functions
207 [ + - ]: 56496702 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
208 [ + - ]: 56496702 : auto B_r = eval_basis( dof_er, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
209 : :
210 : 56496702 : auto wt = wgp[igp] * geoFace(f,0);
211 : :
212 : : std::array< std::vector< real >, 2 > state;
213 : :
214 [ + - ]: 56496702 : state[0] = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim, rdof,
215 : : nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
216 [ + - ]: 56496702 : state[1] = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim, rdof,
217 : : nmat, er, dof_er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P);
218 : :
219 : : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
220 : : "appended boundary state vector" );
221 : : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
222 : : "appended boundary state vector" );
223 : :
224 : : // evaluate prescribed velocity (if any)
225 [ - + ][ - + ]: 112993404 : auto v = vel( ncomp, gp[0], gp[1], gp[2], t );
226 : :
227 : : // compute flux
228 : : auto fl = flux( mat_blk, fn, state, v );
229 : :
230 : : // Add the surface integration term to the rhs
231 [ + - ]: 56496702 : update_rhs_fa( ncomp, nmat, ndof, ndofel[el], ndofel[er], wt, fn,
232 : : el, er, fl, B_l, B_r, R, riemannDeriv );
233 : : }
234 : : }
235 : 65205 : }
236 : :
237 : : void
238 : 56496702 : update_rhs_fa( ncomp_t ncomp,
239 : : std::size_t nmat,
240 : : const std::size_t ndof,
241 : : const std::size_t ndof_l,
242 : : const std::size_t ndof_r,
243 : : const tk::real wt,
244 : : const std::array< tk::real, 3 >& fn,
245 : : const std::size_t el,
246 : : const std::size_t er,
247 : : const std::vector< tk::real >& fl,
248 : : const std::vector< tk::real >& B_l,
249 : : const std::vector< tk::real >& B_r,
250 : : Fields& R,
251 : : std::vector< std::vector< tk::real > >& riemannDeriv )
252 : : // *****************************************************************************
253 : : // Update the rhs by adding the surface integration term
254 : : //! \param[in] ncomp Number of scalar components in this PDE system
255 : : //! \param[in] nmat Number of materials in this PDE system
256 : : //! \param[in] ndof Maximum number of degrees of freedom
257 : : //! \param[in] ndof_l Number of degrees of freedom for left element
258 : : //! \param[in] ndof_r Number of degrees of freedom for right element
259 : : //! \param[in] wt Weight of gauss quadrature point
260 : : //! \param[in] fn Face/Surface normal
261 : : //! \param[in] el Left element index
262 : : //! \param[in] er Right element index
263 : : //! \param[in] fl Surface flux
264 : : //! \param[in] B_l Basis function for the left element
265 : : //! \param[in] B_r Basis function for the right element
266 : : //! \param[in,out] R Right-hand side vector computed
267 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
268 : : //! computed from the Riemann solver for use in the non-conservative terms.
269 : : //! These derivatives are used only for multi-material hydro and unused for
270 : : //! single-material compflow and linear transport.
271 : : // *****************************************************************************
272 : : {
273 : : // following lines commented until rdofel is made available.
274 : : //Assert( B_l.size() == ndof_l, "Size mismatch" );
275 : : //Assert( B_r.size() == ndof_r, "Size mismatch" );
276 : :
277 : : using inciter::newSolidsAccFn;
278 : :
279 : : const auto& solidx =
280 : : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
281 : :
282 [ + + ]: 241047402 : for (ncomp_t c=0; c<ncomp; ++c)
283 : : {
284 : 184550700 : auto mark = c*ndof;
285 [ + + ]: 184550700 : R(el, mark) -= wt * fl[c];
286 : 184550700 : R(er, mark) += wt * fl[c];
287 : :
288 [ + + ]: 184550700 : if(ndof_l > 1) //DG(P1)
289 : : {
290 : 116702415 : R(el, mark+1) -= wt * fl[c] * B_l[1];
291 : 116702415 : R(el, mark+2) -= wt * fl[c] * B_l[2];
292 : 116702415 : R(el, mark+3) -= wt * fl[c] * B_l[3];
293 : : }
294 : :
295 [ + + ]: 184550700 : if(ndof_r > 1) //DG(P1)
296 : : {
297 : 116685900 : R(er, mark+1) += wt * fl[c] * B_r[1];
298 : 116685900 : R(er, mark+2) += wt * fl[c] * B_r[2];
299 : 116685900 : R(er, mark+3) += wt * fl[c] * B_r[3];
300 : : }
301 : :
302 [ + + ]: 184550700 : if(ndof_l > 4) //DG(P2)
303 : : {
304 : 47897100 : R(el, mark+4) -= wt * fl[c] * B_l[4];
305 : 47897100 : R(el, mark+5) -= wt * fl[c] * B_l[5];
306 : 47897100 : R(el, mark+6) -= wt * fl[c] * B_l[6];
307 : 47897100 : R(el, mark+7) -= wt * fl[c] * B_l[7];
308 : 47897100 : R(el, mark+8) -= wt * fl[c] * B_l[8];
309 : 47897100 : R(el, mark+9) -= wt * fl[c] * B_l[9];
310 : : }
311 : :
312 [ + + ]: 184550700 : if(ndof_r > 4) //DG(P2)
313 : : {
314 : 47859930 : R(er, mark+4) += wt * fl[c] * B_r[4];
315 : 47859930 : R(er, mark+5) += wt * fl[c] * B_r[5];
316 : 47859930 : R(er, mark+6) += wt * fl[c] * B_r[6];
317 : 47859930 : R(er, mark+7) += wt * fl[c] * B_r[7];
318 : 47859930 : R(er, mark+8) += wt * fl[c] * B_r[8];
319 : 47859930 : R(er, mark+9) += wt * fl[c] * B_r[9];
320 : : }
321 : : }
322 : :
323 : : // Prep for non-conservative terms in multimat
324 [ + + ]: 56496702 : if (fl.size() > ncomp)
325 : : {
326 : : // Gradients of partial pressures
327 [ + + ]: 18678705 : for (std::size_t k=0; k<nmat; ++k)
328 : : {
329 [ + + ]: 52266480 : for (std::size_t idir=0; idir<3; ++idir)
330 : : {
331 : 39199860 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
332 : 39199860 : riemannDeriv[3*k+idir][er] -= wt * fl[ncomp+k] * fn[idir];
333 : : }
334 : : }
335 : :
336 : : // Divergence of velocity multiples basis fucntion( d(uB) / dx )
337 [ + + ]: 19244250 : for(std::size_t idof = 0; idof < ndof_l; idof++) {
338 : 13632165 : riemannDeriv[3*nmat+idof][el] += wt * fl[ncomp+nmat] * B_l[idof];
339 : : }
340 [ + + ]: 19244250 : for(std::size_t idof = 0; idof < ndof_r; idof++) {
341 : 13632165 : riemannDeriv[3*nmat+idof][er] -= wt * fl[ncomp+nmat] * B_r[idof];
342 : : }
343 : :
344 : : // Divergence of asigma: d(asig_ij)/dx_j
345 [ + + ]: 18678705 : for (std::size_t k=0; k<nmat; ++k)
346 [ - + ]: 13066620 : if (solidx[k] > 0)
347 : : {
348 : 0 : std::size_t mark = ncomp+nmat+1+3*(solidx[k]-1);
349 : :
350 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
351 : : {
352 : 0 : riemannDeriv[3*nmat+ndof+3*(solidx[k]-1)+i][el] -=
353 : 0 : wt * fl[mark+i];
354 : 0 : riemannDeriv[3*nmat+ndof+3*(solidx[k]-1)+i][er] +=
355 : 0 : wt * fl[mark+i];
356 : : }
357 : : }
358 : : }
359 : 56496702 : }
360 : :
361 : : void
362 : 1668 : surfIntFV(
363 : : std::size_t nmat,
364 : : const std::vector< inciter::EOS >& mat_blk,
365 : : real t,
366 : : const std::size_t rdof,
367 : : const std::vector< std::size_t >& inpoel,
368 : : const UnsMesh::Coords& coord,
369 : : const inciter::FaceData& fd,
370 : : const Fields& geoFace,
371 : : const Fields& geoElem,
372 : : const RiemannFluxFn& flux,
373 : : const VelFn& vel,
374 : : const Fields& U,
375 : : const Fields& P,
376 : : const std::vector< int >& srcFlag,
377 : : Fields& R,
378 : : int intsharp )
379 : : // *****************************************************************************
380 : : // Compute internal surface flux integrals for second order FV
381 : : //! \param[in] nmat Number of materials in this PDE system
382 : : //! \param[in] t Physical time
383 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
384 : : //! \param[in] inpoel Element-node connectivity
385 : : //! \param[in] coord Array of nodal coordinates
386 : : //! \param[in] fd Face connectivity and boundary conditions object
387 : : //! \param[in] geoFace Face geometry array
388 : : //! \param[in] geoElem Element geometry array
389 : : //! \param[in] flux Riemann flux function to use
390 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
391 : : //! \param[in] U Solution vector at recent time step
392 : : //! \param[in] P Vector of primitives at recent time step
393 : : //! \param[in] srcFlag Whether the energy source was added
394 : : //! \param[in,out] R Right-hand side vector computed
395 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
396 : : //! default 0, so that it is unused for single-material and transport.
397 : : // *****************************************************************************
398 : : {
399 : : const auto& esuf = fd.Esuf();
400 : :
401 : : const auto& cx = coord[0];
402 : : const auto& cy = coord[1];
403 : : const auto& cz = coord[2];
404 : :
405 : 1668 : auto ncomp = U.nprop()/rdof;
406 : 1668 : auto nprim = P.nprop()/rdof;
407 : :
408 : : // compute internal surface flux integrals
409 [ + + ]: 557508 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
410 : : {
411 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
412 : : "as -1" );
413 : :
414 : 555840 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
415 : 555840 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
416 : :
417 : : // Extract the element coordinates
418 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
419 : 555840 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
420 : 555840 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
421 : 555840 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
422 : 555840 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
423 : :
424 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
425 : 555840 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
426 : 555840 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
427 : 555840 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
428 : 555840 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
429 : :
430 : : // Compute the determinant of Jacobian matrix
431 : : auto detT_l =
432 : 555840 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
433 : : auto detT_r =
434 : 555840 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
435 : :
436 : : // face normal
437 : 555840 : std::array< real, 3 > fn{{geoFace(f,1), geoFace(f,2), geoFace(f,3)}};
438 : :
439 : : // face centroid
440 : 555840 : std::array< real, 3 > gp{{geoFace(f,4), geoFace(f,5), geoFace(f,6)}};
441 : :
442 : : // In order to determine the high-order solution from the left and right
443 : : // elements at the surface quadrature points, the basis functions from
444 : : // the left and right elements are needed. For this, a transformation to
445 : : // the reference coordinates is necessary, since the basis functions are
446 : : // defined on the reference tetrahedron only.
447 : : // The transformation relations are shown below:
448 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
449 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
450 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
451 : :
452 : : std::array< tk::real, 3> ref_gp_l{
453 : 555840 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
454 : 555840 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
455 : 555840 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
456 : : std::array< tk::real, 3> ref_gp_r{
457 : 555840 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
458 : 555840 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
459 : 555840 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
460 : :
461 : : //Compute the basis functions
462 : 555840 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
463 [ + - ]: 555840 : auto B_r = eval_basis( rdof, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
464 : :
465 : : std::array< std::vector< real >, 2 > state;
466 : :
467 [ + - ]: 555840 : state[0] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
468 [ + - ]: 555840 : nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P, srcFlag[el]);
469 [ + - ]: 555840 : state[1] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
470 [ + - ]: 555840 : nmat, er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P, srcFlag[er]);
471 : :
472 : : //safeReco(rdof, nmat, el, er, U, state);
473 : :
474 : : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
475 : : "appended boundary state vector" );
476 : : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
477 : : "appended boundary state vector" );
478 : :
479 : : // evaluate prescribed velocity (if any)
480 [ - + ][ - + ]: 1111680 : auto v = vel( ncomp, gp[0], gp[1], gp[2], t );
481 : :
482 : : // compute flux
483 : : auto fl = flux( mat_blk, fn, state, v );
484 : :
485 : : // compute non-conservative terms
486 [ + - ][ - - ]: 555840 : std::vector< tk::real > var_riemann(nmat+1, 0.0);
487 [ + + ]: 2327200 : for (std::size_t k=0; k<nmat+1; ++k) var_riemann[k] = fl[ncomp+k];
488 : :
489 [ + - ]: 555840 : auto ncf_l = nonConservativeIntFV(nmat, rdof, el, fn, U, P, var_riemann);
490 [ + - ]: 555840 : auto ncf_r = nonConservativeIntFV(nmat, rdof, er, fn, U, P, var_riemann);
491 : :
492 : : // Add the surface integration term to the rhs
493 [ + + ]: 5869920 : for (ncomp_t c=0; c<ncomp; ++c)
494 : : {
495 : 5314080 : R(el, c) -= geoFace(f,0) * (fl[c] - ncf_l[c]);
496 : 5314080 : R(er, c) += geoFace(f,0) * (fl[c] - ncf_r[c]);
497 : : }
498 : : }
499 : 1668 : }
500 : :
501 : : void
502 : 0 : surfIntViscousFV(
503 : : std::size_t nmat,
504 : : const std::vector< inciter::EOS >& mat_blk,
505 : : const std::size_t rdof,
506 : : const std::vector< std::size_t >& inpoel,
507 : : const UnsMesh::Coords& coord,
508 : : const inciter::FaceData& fd,
509 : : const Fields& geoFace,
510 : : const Fields& geoElem,
511 : : const Fields& U,
512 : : const Fields& P,
513 : : const Fields& T,
514 : : const std::vector< int >& srcFlag,
515 : : Fields& R,
516 : : int intsharp )
517 : : // *****************************************************************************
518 : : // Compute internal surface viscous flux integrals for second order FV
519 : : //! \param[in] nmat Number of materials in this PDE system
520 : : //! \param[in] mat_blk Material EOS block
521 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
522 : : //! \param[in] inpoel Element-node connectivity
523 : : //! \param[in] coord Array of nodal coordinates
524 : : //! \param[in] fd Face connectivity and boundary conditions object
525 : : //! \param[in] geoFace Face geometry array
526 : : //! \param[in] geoElem Element geometry array
527 : : //! \param[in] U Solution vector at recent time step
528 : : //! \param[in] P Vector of primitives at recent time step
529 : : //! \param[in] T Vector of temperatures at recent time step
530 : : //! \param[in] srcFlag Whether the energy source was added
531 : : //! \param[in,out] R Right-hand side vector computed
532 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
533 : : //! default 0, so that it is unused for single-material and transport.
534 : : // *****************************************************************************
535 : : {
536 : : using inciter::velocityDofIdx;
537 : :
538 : : const auto& esuf = fd.Esuf();
539 : :
540 : : const auto& cx = coord[0];
541 : : const auto& cy = coord[1];
542 : : const auto& cz = coord[2];
543 : :
544 : 0 : auto ncomp = U.nprop()/rdof;
545 : 0 : auto nprim = P.nprop()/rdof;
546 : :
547 : : // compute internal surface flux integrals
548 [ - - ]: 0 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
549 : : {
550 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
551 : : "as -1" );
552 : :
553 : 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
554 : 0 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
555 : :
556 : : // Extract the element coordinates
557 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
558 : 0 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
559 : 0 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
560 : 0 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
561 : 0 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
562 : :
563 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
564 : 0 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
565 : 0 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
566 : 0 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
567 : 0 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
568 : :
569 : : // Compute the determinant of Jacobian matrix
570 : : auto detT_l =
571 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
572 : : auto detT_r =
573 : 0 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
574 : :
575 : : // face normal
576 : 0 : std::array< real, 3 > fn{{geoFace(f,1), geoFace(f,2), geoFace(f,3)}};
577 : :
578 : : // face centroid
579 : 0 : std::array< real, 3 > gp{{geoFace(f,4), geoFace(f,5), geoFace(f,6)}};
580 : :
581 : : // In order to determine the high-order solution from the left and right
582 : : // elements at the surface quadrature points, the basis functions from
583 : : // the left and right elements are needed. For this, a transformation to
584 : : // the reference coordinates is necessary, since the basis functions are
585 : : // defined on the reference tetrahedron only.
586 : : // The transformation relations are shown below:
587 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
588 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
589 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
590 : :
591 : : std::array< tk::real, 3> ref_gp_l{
592 : 0 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
593 : 0 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
594 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
595 : : std::array< tk::real, 3> ref_gp_r{
596 : 0 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
597 : 0 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
598 : 0 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
599 : :
600 : : //Compute the basis functions
601 : 0 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
602 [ - - ]: 0 : auto B_r = eval_basis( rdof, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
603 : :
604 : : std::array< std::vector< real >, 2 > state;
605 : :
606 [ - - ]: 0 : state[0] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
607 [ - - ]: 0 : nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P, srcFlag[el]);
608 [ - - ]: 0 : state[1] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
609 [ - - ]: 0 : nmat, er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P, srcFlag[er]);
610 : :
611 : : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
612 : : "appended boundary state vector" );
613 : : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
614 : : "appended boundary state vector" );
615 : :
616 : : // cell averaged state for computing the diffusive flux
617 : : std::array< std::vector< real >, 2 > cellAvgState;
618 [ - - ]: 0 : std::vector< tk::real > Bcc(rdof, 0.0);
619 : 0 : Bcc[0] = 1.0;
620 : :
621 [ - - ]: 0 : cellAvgState[0] = evalFVSol(mat_blk, 0, ncomp, nprim, rdof,
622 : : nmat, el, inpoel, coord, geoElem, {{0.25, 0.25, 0.25}}, Bcc, U, P,
623 [ - - ]: 0 : srcFlag[el]);
624 [ - - ][ - - ]: 0 : cellAvgState[1] = evalFVSol(mat_blk, 0, ncomp, nprim, rdof,
625 : : nmat, er, inpoel, coord, geoElem, {{0.25, 0.25, 0.25}}, Bcc, U, P,
626 [ - - ]: 0 : srcFlag[er]);
627 : :
628 : : Assert( cellAvgState[0].size() == ncomp+nprim, "Incorrect size for "
629 : : "appended cell-averaged state vector" );
630 : : Assert( cellAvgState[1].size() == ncomp+nprim, "Incorrect size for "
631 : : "appended cell-averaged state vector" );
632 : :
633 : : std::array< std::vector< real >, 2 > cellAvgT{{
634 [ - - ][ - - ]: 0 : std::vector<real>(nmat), std::vector<real>(nmat) }};
635 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
636 : 0 : cellAvgT[0][k] = T(el, k*rdof);
637 : 0 : cellAvgT[1][k] = T(er, k*rdof);
638 : : }
639 : :
640 : : std::array< std::array< real, 3 >, 2 > centroids{{
641 : 0 : {{geoElem(el,1), geoElem(el,2), geoElem(el,3)}},
642 : 0 : {{geoElem(er,1), geoElem(er,2), geoElem(er,3)}} }};
643 : :
644 : : // Numerical viscous flux
645 : : // -------------------------------------------------------------------------
646 : :
647 : : // 1. Get spatial gradient from Dubiner dofs
648 : : auto jacInv_l =
649 : 0 : tk::inverseJacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
650 [ - - ]: 0 : auto dBdx_l = tk::eval_dBdx_p1( rdof, jacInv_l );
651 : : auto jacInv_r =
652 : 0 : tk::inverseJacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
653 [ - - ]: 0 : auto dBdx_r = tk::eval_dBdx_p1( rdof, jacInv_r );
654 : :
655 : : // 2. Average du_i/dx_j
656 : : std::array< std::array< real, 3 >, 3 > dudx;
657 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
658 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
659 : 0 : dudx[i][j] = 0.5 * (
660 : 0 : dBdx_l[j][1] * P(el, velocityDofIdx(nmat,i,rdof,1))
661 : 0 : + dBdx_l[j][2] * P(el, velocityDofIdx(nmat,i,rdof,2))
662 : 0 : + dBdx_l[j][3] * P(el, velocityDofIdx(nmat,i,rdof,3))
663 : 0 : + dBdx_r[j][1] * P(er, velocityDofIdx(nmat,i,rdof,1))
664 : 0 : + dBdx_r[j][2] * P(er, velocityDofIdx(nmat,i,rdof,2))
665 : 0 : + dBdx_r[j][3] * P(er, velocityDofIdx(nmat,i,rdof,3)) );
666 : :
667 : : // 3. average dT/dx_j
668 : : std::vector< std::array< real, 3 > > dTdx(nmat,
669 [ - - ][ - - ]: 0 : std::array< real, 3 >{{0, 0, 0}});
670 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
671 : 0 : auto mark = k*rdof;
672 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
673 : 0 : dTdx[k][j] = 0.5 * (
674 : 0 : dBdx_l[j][1] * T(el, mark+1)
675 : 0 : + dBdx_l[j][2] * T(el, mark+2)
676 : 0 : + dBdx_l[j][3] * T(el, mark+3)
677 : 0 : + dBdx_r[j][1] * T(er, mark+1)
678 : 0 : + dBdx_r[j][2] * T(er, mark+2)
679 : 0 : + dBdx_r[j][3] * T(er, mark+3) );
680 : : }
681 : : }
682 : :
683 : : // 4. Compute flux
684 : : auto fl = modifiedGradientViscousFlux(nmat, ncomp, fn, centroids, state,
685 [ - - ]: 0 : cellAvgState, cellAvgT, dudx, dTdx);
686 : :
687 : : // -------------------------------------------------------------------------
688 : :
689 : : // Add the surface integration term to the rhs
690 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
691 : : {
692 : 0 : R(el, c) += geoFace(f,0) * fl[c];
693 : 0 : R(er, c) -= geoFace(f,0) * fl[c];
694 : : }
695 : : }
696 : 0 : }
697 : :
698 : : std::vector< real >
699 : 0 : modifiedGradientViscousFlux(
700 : : std::size_t nmat,
701 : : std::size_t ncomp,
702 : : const std::array< tk::real, 3 >& fn,
703 : : const std::array< std::array< tk::real, 3 >, 2 >& centroids,
704 : : const std::array< std::vector< tk::real >, 2 >& state,
705 : : const std::array< std::vector< tk::real >, 2 >& cellAvgState,
706 : : const std::array< std::vector< tk::real >, 2 >& cellAvgT,
707 : : const std::array< std::array< real, 3 >, 3 > dudx,
708 : : const std::vector< std::array< real, 3 > >& dTdx )
709 : : // *****************************************************************************
710 : : // Compute the viscous fluxes from the left and right states
711 : : //! \param[in] nmat Number of materials
712 : : //! \param[in] ncomp Number of component equations in the PDE system
713 : : //! \param[in] fn Face/Surface normal
714 : : //! \param[in] centroids Left and right cell centroids
715 : : //! \param[in] state Left and right unknown/state vector
716 : : //! \param[in] cellAvgState Left and right cell-averaged unknown/state vector
717 : : //! \param[in] cellAvgT Left and right cell-averaged temperature vector
718 : : //! \param[in] dudx Average velocity gradient tensor
719 : : //! \param[in] dTdx Average temperature gradient tensor
720 : : //! \return Numerical viscous flux using the Modified Gradient approach.
721 : : //! \details The average gradient is modified according to Weiss et al. to
722 : : //! obtain a stable discretization (average results in unstable central
723 : : //! central difference).
724 : : //! Ref: Weiss, J. M., Maruszewski, J. P., & Smith, W. A. (1999). Implicit
725 : : //! solution of preconditioned Navier-Stokes equations using algebraic
726 : : //! multigrid. AIAA journal, 37(1), 29-36.
727 : : // *****************************************************************************
728 : : {
729 : : using inciter::velocityDofIdx;
730 : : using inciter::volfracDofIdx;
731 : : using inciter::momentumIdx;
732 : : using inciter::velocityIdx;
733 : : using inciter::volfracIdx;
734 : : using inciter::energyIdx;
735 : :
736 : 0 : std::vector< real > fl(ncomp, 0);
737 : :
738 : : // 1. Modify the average gradient
739 : : std::array< real, 3 > r_f{{
740 : 0 : centroids[1][0] - centroids[0][0],
741 : 0 : centroids[1][1] - centroids[0][1],
742 : 0 : centroids[1][2] - centroids[0][2] }};
743 : 0 : real r_mag = std::sqrt(tk::dot(r_f, r_f));
744 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
745 : 0 : r_f[i] /= r_mag;
746 : :
747 : : // velocity gradient
748 : 0 : auto dudx_m = dudx;
749 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
750 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
751 : 0 : dudx_m[i][j] -= ( tk::dot(dudx_m[i], r_f) -
752 : 0 : (cellAvgState[1][ncomp+velocityIdx(nmat,i)]
753 : 0 : - cellAvgState[0][ncomp+velocityIdx(nmat,i)])/r_mag ) * r_f[j];
754 : : // temperature gradient
755 [ - - ]: 0 : auto dTdx_m = dTdx;
756 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
757 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
758 : 0 : dTdx_m[k][j] -= ( tk::dot(dTdx_m[k], r_f) -
759 : 0 : (cellAvgT[1][k] - cellAvgT[0][k])/r_mag ) * r_f[j];
760 : : }
761 : :
762 : : // 2. Compute viscous stress tensor
763 : : std::array< real, 6 > tau;
764 : : real mu(0.0);
765 [ - - ][ - - ]: 0 : std::vector< real > alLR(nmat, 0), conduct_mat(nmat, 0);
[ - - ][ - - ]
766 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
767 : 0 : alLR[k] = 0.5*( state[0][volfracIdx(nmat,k)] + state[1][volfracIdx(nmat,k)] );
768 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
769 [ - - ]: 0 : mu += alLR[k] * inciter::getmatprop< tag::mu >(k);
770 [ - - ]: 0 : conduct_mat[k] = inciter::getmatprop< tag::mu >(k) *
771 [ - - ][ - - ]: 0 : inciter::getmatprop< tag::cv >(k) * inciter::getmatprop< tag::gamma >(k)
772 : 0 : / 0.71;
773 : : }
774 : :
775 : 0 : tau[0] = mu * ( 4.0 * dudx_m[0][0] - 2.0*(dudx_m[1][1] + dudx_m[2][2]) ) / 3.0;
776 : 0 : tau[1] = mu * ( 4.0 * dudx_m[1][1] - 2.0*(dudx_m[0][0] + dudx_m[2][2]) ) / 3.0;
777 : 0 : tau[2] = mu * ( 4.0 * dudx_m[2][2] - 2.0*(dudx_m[0][0] + dudx_m[1][1]) ) / 3.0;
778 : 0 : tau[3] = mu * ( dudx_m[0][1] + dudx_m[1][0] );
779 : 0 : tau[4] = mu * ( dudx_m[0][2] + dudx_m[2][0] );
780 : 0 : tau[5] = mu * ( dudx_m[1][2] + dudx_m[2][1] );
781 : :
782 : : // 3. Compute viscous flux terms
783 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
784 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
785 : 0 : fl[momentumIdx(nmat, i)] += tau[inciter::stressCmp[i][j]] * fn[j];
786 : :
787 [ - - ][ - - ]: 0 : std::vector< std::array< real, 3 > > energyFlux(nmat, {{0.0, 0.0, 0.0}});
788 : : std::array< real, 3 > uAvg{{
789 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,0)] + state[1][ncomp+velocityIdx(nmat,0)]),
790 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,1)] + state[1][ncomp+velocityIdx(nmat,1)]),
791 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,2)] + state[1][ncomp+velocityIdx(nmat,2)])
792 : 0 : }};
793 : :
794 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
795 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
796 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
797 : 0 : energyFlux[k][j] += uAvg[i] * tau[inciter::stressCmp[i][j]] +
798 : 0 : conduct_mat[k] * dTdx_m[k][j];
799 : : }
800 : :
801 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
802 : 0 : fl[energyIdx(nmat, k)] = alLR[k]
803 : 0 : * tk::dot(energyFlux[k], fn);
804 : : }
805 : :
806 : 0 : return fl;
807 : : }
808 : :
809 : : } // tk::
|
