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 : 37410 : 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 : 37410 : auto ncomp = U.nprop()/rdof;
92 : 37410 : 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 [ + + ]: 17335065 : 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 : 17297655 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
107 : 17297655 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
108 : :
109 : 17297655 : auto ng_l = tk::NGfa(ndofel[el]);
110 [ + + ]: 17297655 : 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 : 17297655 : 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 [ + - ]: 17297655 : coordgp[0].resize( ng );
121 [ + - ]: 17297655 : coordgp[1].resize( ng );
122 [ + - ]: 17297655 : wgp.resize( ng );
123 : :
124 : : // get quadrature point weights and coordinates for triangle
125 [ + - ]: 17297655 : GaussQuadratureTri( ng, coordgp, wgp );
126 : :
127 : : // Extract the element coordinates
128 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
129 : 17297655 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
130 : 17297655 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
131 : 17297655 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
132 : 17297655 : {{ 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 : 17297655 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
136 : 17297655 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
137 : 17297655 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
138 : 17297655 : {{ 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 : 17297655 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
143 : : auto detT_r =
144 : 17297655 : 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 : 17297655 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
149 : 17297655 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
150 : 17297655 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
151 : :
152 : : std::array< real, 3 >
153 : 17297655 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
154 : :
155 : : // Gaussian quadrature
156 [ + + ]: 45270162 : for (std::size_t igp=0; igp<ng; ++igp)
157 : : {
158 : : // Compute the coordinates of quadrature point at physical domain
159 [ + - ]: 27972507 : 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 [ + + ]: 27972507 : if (rdof > ndof)
178 : : {
179 : : dof_el = rdof;
180 : : dof_er = rdof;
181 : : }
182 : : else
183 : : {
184 : 25881357 : dof_el = ndofel[el];
185 : 25881357 : 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 [ - + ]: 27972507 : 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 : 27972507 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
199 : 27972507 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
200 : 27972507 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
201 : : std::array< tk::real, 3> ref_gp_r{
202 : 27972507 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
203 : 27972507 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
204 : 27972507 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
205 : :
206 : : //Compute the basis functions
207 [ + - ]: 27972507 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
208 [ + - ]: 27972507 : auto B_r = eval_basis( dof_er, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
209 : :
210 : 27972507 : auto wt = wgp[igp] * geoFace(f,0);
211 : :
212 : : std::array< std::vector< real >, 2 > state;
213 : :
214 [ + - ]: 27972507 : 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 [ + - ]: 27972507 : 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 [ - + ][ - + ]: 55945014 : 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 [ + - ]: 27972507 : 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 : 37410 : }
236 : :
237 : : void
238 : 27972507 : 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 [ + + ]: 183532542 : for (ncomp_t c=0; c<ncomp; ++c)
283 : : {
284 : 155560035 : auto mark = c*ndof;
285 [ + + ]: 155560035 : R(el, mark) -= wt * fl[c];
286 : 155560035 : R(er, mark) += wt * fl[c];
287 : :
288 [ + + ]: 155560035 : if(ndof_l > 1) //DG(P1)
289 : : {
290 : 78462045 : R(el, mark+1) -= wt * fl[c] * B_l[1];
291 : 78462045 : R(el, mark+2) -= wt * fl[c] * B_l[2];
292 : 78462045 : R(el, mark+3) -= wt * fl[c] * B_l[3];
293 : : }
294 : :
295 [ + + ]: 155560035 : if(ndof_r > 1) //DG(P1)
296 : : {
297 : 78449670 : R(er, mark+1) += wt * fl[c] * B_r[1];
298 : 78449670 : R(er, mark+2) += wt * fl[c] * B_r[2];
299 : 78449670 : R(er, mark+3) += wt * fl[c] * B_r[3];
300 : : }
301 : :
302 [ + + ]: 155560035 : if(ndof_l > 4) //DG(P2)
303 : : {
304 : 32982480 : R(el, mark+4) -= wt * fl[c] * B_l[4];
305 : 32982480 : R(el, mark+5) -= wt * fl[c] * B_l[5];
306 : 32982480 : R(el, mark+6) -= wt * fl[c] * B_l[6];
307 : 32982480 : R(el, mark+7) -= wt * fl[c] * B_l[7];
308 : 32982480 : R(el, mark+8) -= wt * fl[c] * B_l[8];
309 : 32982480 : R(el, mark+9) -= wt * fl[c] * B_l[9];
310 : : }
311 : :
312 [ + + ]: 155560035 : if(ndof_r > 4) //DG(P2)
313 : : {
314 : 32971050 : R(er, mark+4) += wt * fl[c] * B_r[4];
315 : 32971050 : R(er, mark+5) += wt * fl[c] * B_r[5];
316 : 32971050 : R(er, mark+6) += wt * fl[c] * B_r[6];
317 : 32971050 : R(er, mark+7) += wt * fl[c] * B_r[7];
318 : 32971050 : R(er, mark+8) += wt * fl[c] * B_r[8];
319 : 32971050 : R(er, mark+9) += wt * fl[c] * B_r[9];
320 : : }
321 : : }
322 : :
323 : : // Prep for non-conservative terms in multimat
324 [ + + ]: 27972507 : if (fl.size() > ncomp)
325 : : {
326 : : // Gradients of partial pressures
327 [ + + ]: 17795355 : for (std::size_t k=0; k<nmat; ++k)
328 : : {
329 [ + + ]: 49910880 : for (std::size_t idir=0; idir<3; ++idir)
330 : : {
331 : 37433160 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
332 : 37433160 : 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 [ + + ]: 15751500 : for(std::size_t idof = 0; idof < ndof_l; idof++) {
338 : 10433865 : riemannDeriv[3*nmat+idof][el] += wt * fl[ncomp+nmat] * B_l[idof];
339 : : }
340 [ + + ]: 15751500 : for(std::size_t idof = 0; idof < ndof_r; idof++) {
341 : 10433865 : riemannDeriv[3*nmat+idof][er] -= wt * fl[ncomp+nmat] * B_r[idof];
342 : : }
343 : :
344 : : // Divergence of asigma: d(asig_ij)/dx_j
345 [ + + ]: 17795355 : for (std::size_t k=0; k<nmat; ++k)
346 [ - + ]: 12477720 : 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 : : // Derivatives of g: d(g_il)/d(x_j)-d(g_ij)/d(x_l)
360 : : // for i=1,2,3; j=1,2,3; l=1,2,3. Total = 3x3x3 (per solid)
361 : 5317635 : std::size_t nsld = inciter::numSolids(nmat, solidx);
362 [ + + ]: 17795355 : for (std::size_t k=0; k<nmat; ++k)
363 [ - + ]: 12477720 : if (solidx[k] > 0)
364 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
365 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
366 [ - - ]: 0 : for (std::size_t l=0; l<3; ++l)
367 [ - - ]: 0 : if (j != l)
368 : : {
369 : 0 : std::size_t mark1 = ncomp+nmat+1+3*nsld+9*(solidx[k]-1)+3*i+l;
370 : 0 : std::size_t mark2 = ncomp+nmat+1+3*nsld+9*(solidx[k]-1)+3*i+j;
371 : 0 : riemannDeriv[3*nmat+ndof+3*nsld+newSolidsAccFn(k,i,j,l)][el] -=
372 : 0 : wt * ( fl[mark1] * fn[j] - fl[mark2] * fn[l]);
373 : 0 : riemannDeriv[3*nmat+ndof+3*nsld+newSolidsAccFn(k,i,j,l)][er] +=
374 : 0 : wt * ( fl[mark1] * fn[j] - fl[mark2] * fn[l]);
375 : : }
376 : : }
377 : 27972507 : }
378 : :
379 : : void
380 : 1368 : surfIntFV(
381 : : std::size_t nmat,
382 : : const std::vector< inciter::EOS >& mat_blk,
383 : : real t,
384 : : const std::size_t rdof,
385 : : const std::vector< std::size_t >& inpoel,
386 : : const UnsMesh::Coords& coord,
387 : : const inciter::FaceData& fd,
388 : : const Fields& geoFace,
389 : : const Fields& geoElem,
390 : : const RiemannFluxFn& flux,
391 : : const VelFn& vel,
392 : : const Fields& U,
393 : : const Fields& P,
394 : : const std::vector< int >& srcFlag,
395 : : Fields& R,
396 : : int intsharp )
397 : : // *****************************************************************************
398 : : // Compute internal surface flux integrals for second order FV
399 : : //! \param[in] nmat Number of materials in this PDE system
400 : : //! \param[in] t Physical time
401 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
402 : : //! \param[in] inpoel Element-node connectivity
403 : : //! \param[in] coord Array of nodal coordinates
404 : : //! \param[in] fd Face connectivity and boundary conditions object
405 : : //! \param[in] geoFace Face geometry array
406 : : //! \param[in] geoElem Element geometry array
407 : : //! \param[in] flux Riemann flux function to use
408 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
409 : : //! \param[in] U Solution vector at recent time step
410 : : //! \param[in] P Vector of primitives at recent time step
411 : : //! \param[in] srcFlag Whether the energy source was added
412 : : //! \param[in,out] R Right-hand side vector computed
413 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
414 : : //! default 0, so that it is unused for single-material and transport.
415 : : // *****************************************************************************
416 : : {
417 : : const auto& esuf = fd.Esuf();
418 : : const auto& localFaceId = fd.FaceLocalId();
419 : :
420 : 1368 : auto ncomp = U.nprop()/rdof;
421 : 1368 : auto nprim = P.nprop()/rdof;
422 : :
423 : : // Basis functions for all face-centroids of element e
424 : : std::array< std::vector< tk::real >, 4 > Bf_array = {
425 : : tk::eval_basis(rdof, tk::fc_coord[0][0], tk::fc_coord[0][1], tk::fc_coord[0][2]),
426 : : tk::eval_basis(rdof, tk::fc_coord[1][0], tk::fc_coord[1][1], tk::fc_coord[1][2]),
427 : : tk::eval_basis(rdof, tk::fc_coord[2][0], tk::fc_coord[2][1], tk::fc_coord[2][2]),
428 : : tk::eval_basis(rdof, tk::fc_coord[3][0], tk::fc_coord[3][1], tk::fc_coord[3][2])
429 [ + - ][ + - ]: 1368 : };
[ + - ][ + - ]
430 : :
431 : : // compute internal surface flux integrals
432 [ + + ]: 476888 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
433 : : {
434 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
435 : : "as -1" );
436 : :
437 : 475520 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
438 : 475520 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
439 : :
440 : : // face normal
441 : 475520 : std::array< real, 3 > fn{{geoFace(f,1), geoFace(f,2), geoFace(f,3)}};
442 : :
443 : : // face centroid
444 : 475520 : std::array< real, 3 > gp{{geoFace(f,4), geoFace(f,5), geoFace(f,6)}};
445 : :
446 : 475520 : auto f_Lid = static_cast< std::size_t >(localFaceId[2*f]);
447 : 475520 : auto f_Rid = static_cast< std::size_t >(localFaceId[2*f+1]);
448 : :
449 : 475520 : auto ref_gp_l = tk::fc_coord[f_Lid];
450 : 475520 : auto ref_gp_r = tk::fc_coord[f_Rid];
451 : :
452 : : // Compute the basis functions
453 : : const auto& B_l = Bf_array[f_Lid];
454 : : const auto& B_r = Bf_array[f_Rid];
455 : :
456 : : std::array< std::vector< real >, 2 > state;
457 : :
458 [ + - ]: 475520 : state[0] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
459 [ + - ]: 475520 : nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P, srcFlag[el]);
460 [ + - ][ - + ]: 951040 : state[1] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
461 [ + - ]: 475520 : nmat, er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P, srcFlag[er]);
462 : :
463 : : //safeReco(rdof, nmat, el, er, U, state);
464 : :
465 : : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
466 : : "appended boundary state vector" );
467 : : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
468 : : "appended boundary state vector" );
469 : :
470 : : // evaluate prescribed velocity (if any)
471 [ - + ]: 475520 : auto v = vel( ncomp, gp[0], gp[1], gp[2], t );
472 : :
473 : : // compute flux
474 : : auto fl = flux( mat_blk, fn, state, v );
475 : :
476 : : // compute non-conservative terms
477 [ + - ][ - - ]: 475520 : std::vector< tk::real > var_riemann(nmat+1, 0.0);
478 [ + + ]: 2005920 : for (std::size_t k=0; k<nmat+1; ++k) var_riemann[k] = fl[ncomp+k];
479 : :
480 [ + - ]: 475520 : auto ncf_l = nonConservativeIntFV(nmat, rdof, el, fn, U, P, var_riemann);
481 [ + - ]: 475520 : auto ncf_r = nonConservativeIntFV(nmat, rdof, er, fn, U, P, var_riemann);
482 : :
483 : : // Add the surface integration term to the rhs
484 [ + + ]: 5066720 : for (ncomp_t c=0; c<ncomp; ++c)
485 : : {
486 : 4591200 : R(el, c) -= geoFace(f,0) * (fl[c] - ncf_l[c]);
487 : 4591200 : R(er, c) += geoFace(f,0) * (fl[c] - ncf_r[c]);
488 : : }
489 : : }
490 : 1368 : }
491 : :
492 : : void
493 : 0 : surfIntViscousFV(
494 : : std::size_t nmat,
495 : : const std::vector< inciter::EOS >& mat_blk,
496 : : const std::size_t rdof,
497 : : const std::vector< std::size_t >& inpoel,
498 : : const UnsMesh::Coords& coord,
499 : : const inciter::FaceData& fd,
500 : : const Fields& geoFace,
501 : : const Fields& geoElem,
502 : : const Fields& U,
503 : : const Fields& P,
504 : : const Fields& T,
505 : : const std::vector< int >& srcFlag,
506 : : Fields& R,
507 : : int intsharp )
508 : : // *****************************************************************************
509 : : // Compute internal surface viscous flux integrals for second order FV
510 : : //! \param[in] nmat Number of materials in this PDE system
511 : : //! \param[in] mat_blk Material EOS block
512 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
513 : : //! \param[in] inpoel Element-node connectivity
514 : : //! \param[in] coord Array of nodal coordinates
515 : : //! \param[in] fd Face connectivity and boundary conditions object
516 : : //! \param[in] geoFace Face geometry array
517 : : //! \param[in] geoElem Element geometry array
518 : : //! \param[in] U Solution vector at recent time step
519 : : //! \param[in] P Vector of primitives at recent time step
520 : : //! \param[in] T Vector of temperatures at recent time step
521 : : //! \param[in] srcFlag Whether the energy source was added
522 : : //! \param[in,out] R Right-hand side vector computed
523 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
524 : : //! default 0, so that it is unused for single-material and transport.
525 : : // *****************************************************************************
526 : : {
527 : : using inciter::velocityDofIdx;
528 : :
529 : : const auto& esuf = fd.Esuf();
530 : :
531 : : const auto& cx = coord[0];
532 : : const auto& cy = coord[1];
533 : : const auto& cz = coord[2];
534 : :
535 : 0 : auto ncomp = U.nprop()/rdof;
536 : 0 : auto nprim = P.nprop()/rdof;
537 : :
538 : : // compute internal surface flux integrals
539 [ - - ]: 0 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
540 : : {
541 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
542 : : "as -1" );
543 : :
544 : 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
545 : 0 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
546 : :
547 : : // Extract the element coordinates
548 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
549 : 0 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
550 : 0 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
551 : 0 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
552 : 0 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
553 : :
554 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
555 : 0 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
556 : 0 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
557 : 0 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
558 : 0 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
559 : :
560 : : // Compute the determinant of Jacobian matrix
561 : : auto detT_l =
562 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
563 : : auto detT_r =
564 : 0 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
565 : :
566 : : // face normal
567 : 0 : std::array< real, 3 > fn{{geoFace(f,1), geoFace(f,2), geoFace(f,3)}};
568 : :
569 : : // face centroid
570 : 0 : std::array< real, 3 > gp{{geoFace(f,4), geoFace(f,5), geoFace(f,6)}};
571 : :
572 : : // In order to determine the high-order solution from the left and right
573 : : // elements at the surface quadrature points, the basis functions from
574 : : // the left and right elements are needed. For this, a transformation to
575 : : // the reference coordinates is necessary, since the basis functions are
576 : : // defined on the reference tetrahedron only.
577 : : // The transformation relations are shown below:
578 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
579 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
580 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
581 : :
582 : : std::array< tk::real, 3> ref_gp_l{
583 : 0 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
584 : 0 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
585 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
586 : : std::array< tk::real, 3> ref_gp_r{
587 : 0 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
588 : 0 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
589 : 0 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
590 : :
591 : : //Compute the basis functions
592 : 0 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
593 [ - - ]: 0 : auto B_r = eval_basis( rdof, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
594 : :
595 : : std::array< std::vector< real >, 2 > state;
596 : :
597 [ - - ]: 0 : state[0] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
598 [ - - ]: 0 : nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P, srcFlag[el]);
599 [ - - ]: 0 : state[1] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
600 [ - - ]: 0 : nmat, er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P, srcFlag[er]);
601 : :
602 : : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
603 : : "appended boundary state vector" );
604 : : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
605 : : "appended boundary state vector" );
606 : :
607 : : // cell averaged state for computing the diffusive flux
608 : : std::array< std::vector< real >, 2 > cellAvgState;
609 [ - - ]: 0 : std::vector< tk::real > Bcc(rdof, 0.0);
610 : 0 : Bcc[0] = 1.0;
611 : :
612 [ - - ]: 0 : cellAvgState[0] = evalFVSol(mat_blk, 0, ncomp, nprim, rdof,
613 : : nmat, el, inpoel, coord, geoElem, {{0.25, 0.25, 0.25}}, Bcc, U, P,
614 [ - - ]: 0 : srcFlag[el]);
615 [ - - ][ - - ]: 0 : cellAvgState[1] = evalFVSol(mat_blk, 0, ncomp, nprim, rdof,
616 : : nmat, er, inpoel, coord, geoElem, {{0.25, 0.25, 0.25}}, Bcc, U, P,
617 [ - - ]: 0 : srcFlag[er]);
618 : :
619 : : Assert( cellAvgState[0].size() == ncomp+nprim, "Incorrect size for "
620 : : "appended cell-averaged state vector" );
621 : : Assert( cellAvgState[1].size() == ncomp+nprim, "Incorrect size for "
622 : : "appended cell-averaged state vector" );
623 : :
624 : : std::array< std::vector< real >, 2 > cellAvgT{{
625 [ - - ][ - - ]: 0 : std::vector<real>(nmat), std::vector<real>(nmat) }};
626 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
627 : 0 : cellAvgT[0][k] = T(el, k*rdof);
628 : 0 : cellAvgT[1][k] = T(er, k*rdof);
629 : : }
630 : :
631 : : std::array< std::array< real, 3 >, 2 > centroids{{
632 : 0 : {{geoElem(el,1), geoElem(el,2), geoElem(el,3)}},
633 : 0 : {{geoElem(er,1), geoElem(er,2), geoElem(er,3)}} }};
634 : :
635 : : // Numerical viscous flux
636 : : // -------------------------------------------------------------------------
637 : :
638 : : // 1. Get spatial gradient from Dubiner dofs
639 : : auto jacInv_l =
640 : 0 : tk::inverseJacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
641 [ - - ]: 0 : auto dBdx_l = tk::eval_dBdx_p1( rdof, jacInv_l );
642 : : auto jacInv_r =
643 : 0 : tk::inverseJacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
644 [ - - ]: 0 : auto dBdx_r = tk::eval_dBdx_p1( rdof, jacInv_r );
645 : :
646 : : // 2. Average du_i/dx_j
647 : : std::array< std::array< real, 3 >, 3 > dudx;
648 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
649 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
650 : 0 : dudx[i][j] = 0.5 * (
651 : 0 : dBdx_l[j][1] * P(el, velocityDofIdx(nmat,i,rdof,1))
652 : 0 : + dBdx_l[j][2] * P(el, velocityDofIdx(nmat,i,rdof,2))
653 : 0 : + dBdx_l[j][3] * P(el, velocityDofIdx(nmat,i,rdof,3))
654 : 0 : + dBdx_r[j][1] * P(er, velocityDofIdx(nmat,i,rdof,1))
655 : 0 : + dBdx_r[j][2] * P(er, velocityDofIdx(nmat,i,rdof,2))
656 : 0 : + dBdx_r[j][3] * P(er, velocityDofIdx(nmat,i,rdof,3)) );
657 : :
658 : : // 3. average dT/dx_j
659 : : std::vector< std::array< real, 3 > > dTdx(nmat,
660 [ - - ][ - - ]: 0 : std::array< real, 3 >{{0, 0, 0}});
661 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
662 : 0 : auto mark = k*rdof;
663 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
664 : 0 : dTdx[k][j] = 0.5 * (
665 : 0 : dBdx_l[j][1] * T(el, mark+1)
666 : 0 : + dBdx_l[j][2] * T(el, mark+2)
667 : 0 : + dBdx_l[j][3] * T(el, mark+3)
668 : 0 : + dBdx_r[j][1] * T(er, mark+1)
669 : 0 : + dBdx_r[j][2] * T(er, mark+2)
670 : 0 : + dBdx_r[j][3] * T(er, mark+3) );
671 : : }
672 : : }
673 : :
674 : : // 4. Compute flux
675 : : auto fl = modifiedGradientViscousFlux(nmat, ncomp, fn, centroids, state,
676 [ - - ]: 0 : cellAvgState, cellAvgT, dudx, dTdx);
677 : :
678 : : // -------------------------------------------------------------------------
679 : :
680 : : // Add the surface integration term to the rhs
681 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
682 : : {
683 : 0 : R(el, c) += geoFace(f,0) * fl[c];
684 : 0 : R(er, c) -= geoFace(f,0) * fl[c];
685 : : }
686 : : }
687 : 0 : }
688 : :
689 : : std::vector< real >
690 : 0 : modifiedGradientViscousFlux(
691 : : std::size_t nmat,
692 : : std::size_t ncomp,
693 : : const std::array< tk::real, 3 >& fn,
694 : : const std::array< std::array< tk::real, 3 >, 2 >& centroids,
695 : : const std::array< std::vector< tk::real >, 2 >& state,
696 : : const std::array< std::vector< tk::real >, 2 >& cellAvgState,
697 : : const std::array< std::vector< tk::real >, 2 >& cellAvgT,
698 : : const std::array< std::array< real, 3 >, 3 > dudx,
699 : : const std::vector< std::array< real, 3 > >& dTdx )
700 : : // *****************************************************************************
701 : : // Compute the viscous fluxes from the left and right states
702 : : //! \param[in] nmat Number of materials
703 : : //! \param[in] ncomp Number of component equations in the PDE system
704 : : //! \param[in] fn Face/Surface normal
705 : : //! \param[in] centroids Left and right cell centroids
706 : : //! \param[in] state Left and right unknown/state vector
707 : : //! \param[in] cellAvgState Left and right cell-averaged unknown/state vector
708 : : //! \param[in] cellAvgT Left and right cell-averaged temperature vector
709 : : //! \param[in] dudx Average velocity gradient tensor
710 : : //! \param[in] dTdx Average temperature gradient tensor
711 : : //! \return Numerical viscous flux using the Modified Gradient approach.
712 : : //! \details The average gradient is modified according to Weiss et al. to
713 : : //! obtain a stable discretization (average results in unstable central
714 : : //! central difference).
715 : : //! Ref: Weiss, J. M., Maruszewski, J. P., & Smith, W. A. (1999). Implicit
716 : : //! solution of preconditioned Navier-Stokes equations using algebraic
717 : : //! multigrid. AIAA journal, 37(1), 29-36.
718 : : // *****************************************************************************
719 : : {
720 : : using inciter::velocityDofIdx;
721 : : using inciter::volfracDofIdx;
722 : : using inciter::momentumIdx;
723 : : using inciter::velocityIdx;
724 : : using inciter::volfracIdx;
725 : : using inciter::energyIdx;
726 : :
727 : 0 : std::vector< real > fl(ncomp, 0);
728 : :
729 : : // 1. Modify the average gradient
730 : : std::array< real, 3 > r_f{{
731 : 0 : centroids[1][0] - centroids[0][0],
732 : 0 : centroids[1][1] - centroids[0][1],
733 : 0 : centroids[1][2] - centroids[0][2] }};
734 : 0 : real r_mag = std::sqrt(tk::dot(r_f, r_f));
735 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
736 : 0 : r_f[i] /= r_mag;
737 : :
738 : : // velocity gradient
739 : 0 : auto dudx_m = dudx;
740 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
741 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
742 : 0 : dudx_m[i][j] -= ( tk::dot(dudx_m[i], r_f) -
743 : 0 : (cellAvgState[1][ncomp+velocityIdx(nmat,i)]
744 : 0 : - cellAvgState[0][ncomp+velocityIdx(nmat,i)])/r_mag ) * r_f[j];
745 : : // temperature gradient
746 [ - - ]: 0 : auto dTdx_m = dTdx;
747 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
748 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
749 : 0 : dTdx_m[k][j] -= ( tk::dot(dTdx_m[k], r_f) -
750 : 0 : (cellAvgT[1][k] - cellAvgT[0][k])/r_mag ) * r_f[j];
751 : : }
752 : :
753 : : // 2. Compute viscous stress tensor
754 : : std::array< real, 6 > tau;
755 : : real mu(0.0);
756 [ - - ][ - - ]: 0 : std::vector< real > alLR(nmat, 0), conduct_mat(nmat, 0);
[ - - ][ - - ]
757 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
758 : 0 : alLR[k] = 0.5*( state[0][volfracIdx(nmat,k)] + state[1][volfracIdx(nmat,k)] );
759 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
760 [ - - ]: 0 : mu += alLR[k] * inciter::getmatprop< tag::mu >(k);
761 [ - - ]: 0 : conduct_mat[k] = inciter::getmatprop< tag::mu >(k) *
762 [ - - ][ - - ]: 0 : inciter::getmatprop< tag::cv >(k) * inciter::getmatprop< tag::gamma >(k)
763 : 0 : / 0.71;
764 : : }
765 : :
766 : 0 : tau[0] = mu * ( 4.0 * dudx_m[0][0] - 2.0*(dudx_m[1][1] + dudx_m[2][2]) ) / 3.0;
767 : 0 : tau[1] = mu * ( 4.0 * dudx_m[1][1] - 2.0*(dudx_m[0][0] + dudx_m[2][2]) ) / 3.0;
768 : 0 : tau[2] = mu * ( 4.0 * dudx_m[2][2] - 2.0*(dudx_m[0][0] + dudx_m[1][1]) ) / 3.0;
769 : 0 : tau[3] = mu * ( dudx_m[0][1] + dudx_m[1][0] );
770 : 0 : tau[4] = mu * ( dudx_m[0][2] + dudx_m[2][0] );
771 : 0 : tau[5] = mu * ( dudx_m[1][2] + dudx_m[2][1] );
772 : :
773 : : // 3. Compute viscous flux terms
774 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
775 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
776 : 0 : fl[momentumIdx(nmat, i)] += tau[inciter::stressCmp[i][j]] * fn[j];
777 : :
778 [ - - ][ - - ]: 0 : std::vector< std::array< real, 3 > > energyFlux(nmat, {{0.0, 0.0, 0.0}});
779 : : std::array< real, 3 > uAvg{{
780 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,0)] + state[1][ncomp+velocityIdx(nmat,0)]),
781 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,1)] + state[1][ncomp+velocityIdx(nmat,1)]),
782 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,2)] + state[1][ncomp+velocityIdx(nmat,2)])
783 : 0 : }};
784 : :
785 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
786 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
787 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
788 : 0 : energyFlux[k][j] += uAvg[i] * tau[inciter::stressCmp[i][j]] +
789 : 0 : conduct_mat[k] * dTdx_m[k][j];
790 : : }
791 : :
792 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
793 : 0 : fl[energyIdx(nmat, k)] = alLR[k]
794 : 0 : * tk::dot(energyFlux[k], fn);
795 : : }
796 : :
797 : 0 : return fl;
798 : : }
799 : :
800 : : } // tk::
|
