Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Integrate/Boundary.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 physical boundary surface integrals of a
9 : : system of PDEs in DG methods
10 : : \details This file contains functionality for computing physical boundary
11 : : surface integrals of a system of PDEs used in discontinuous Galerkin
12 : : methods for various orders of numerical representation.
13 : : */
14 : : // *****************************************************************************
15 : :
16 : : #include <array>
17 : :
18 : : #include "Basis.hpp"
19 : : #include "Boundary.hpp"
20 : : #include "Vector.hpp"
21 : : #include "Quadrature.hpp"
22 : : #include "MultiMatTerms.hpp"
23 : : #include "MultiMat/MultiMatIndexing.hpp"
24 : : #include "Reconstruction.hpp"
25 : : #include "Inciter/InputDeck/InputDeck.hpp"
26 : :
27 : : namespace inciter {
28 : : extern ctr::InputDeck g_inputdeck;
29 : : }
30 : :
31 : : namespace tk {
32 : :
33 : : void
34 : 275040 : bndSurfInt( const bool pref,
35 : : std::size_t nmat,
36 : : const std::vector< inciter::EOS >& mat_blk,
37 : : const std::size_t ndof,
38 : : const std::size_t rdof,
39 : : const std::vector< std::size_t >& bcconfig,
40 : : const inciter::FaceData& fd,
41 : : const Fields& geoFace,
42 : : const Fields& geoElem,
43 : : const std::vector< std::size_t >& inpoel,
44 : : const UnsMesh::Coords& coord,
45 : : real t,
46 : : const RiemannFluxFn& flux,
47 : : const VelFn& vel,
48 : : const StateFn& state,
49 : : const Fields& U,
50 : : const Fields& P,
51 : : const std::vector< std::size_t >& ndofel,
52 : : Fields& R,
53 : : std::vector< std::vector< tk::real > >& riemannDeriv,
54 : : int intsharp )
55 : : // *****************************************************************************
56 : : //! Compute boundary surface flux integrals for a given boundary type for DG
57 : : //! \details This function computes contributions from surface integrals along
58 : : //! all faces for a particular boundary condition type, configured by the state
59 : : //! function
60 : : //! \param[in] pref Indicator for p-adaptive algorithm
61 : : //! \param[in] nmat Number of materials in this PDE system
62 : : //! \param[in] mat_blk EOS material block
63 : : //! \param[in] ndof Maximum number of degrees of freedom
64 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
65 : : //! \param[in] bcconfig BC configuration vector for multiple side sets
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] inpoel Element-node connectivity
70 : : //! \param[in] coord Array of nodal coordinates
71 : : //! \param[in] t Physical time
72 : : //! \param[in] flux Riemann flux function to use
73 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
74 : : //! \param[in] state Function to evaluate the left and right solution state at
75 : : //! boundaries
76 : : //! \param[in] U Solution vector at recent time step
77 : : //! \param[in] P Vector of primitives at recent time step
78 : : //! \param[in] ndofel Vector of local number of degrees of freedom
79 : : //! \param[in,out] R Right-hand side vector computed
80 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
81 : : //! computed from the Riemann solver for use in the non-conservative terms.
82 : : //! These derivatives are used only for multi-material hydro and unused for
83 : : //! single-material compflow and linear transport.
84 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
85 : : //! default 0, so that it is unused for single-material and transport.
86 : : // *****************************************************************************
87 : : {
88 : 275040 : const auto& bface = fd.Bface();
89 : 275040 : const auto& esuf = fd.Esuf();
90 : 275040 : const auto& inpofa = fd.Inpofa();
91 : :
92 : 275040 : const auto& cx = coord[0];
93 : 275040 : const auto& cy = coord[1];
94 : 275040 : const auto& cz = coord[2];
95 : :
96 : 275040 : auto ncomp = U.nprop()/rdof;
97 : 275040 : auto nprim = P.nprop()/rdof;
98 : :
99 : : //Assert( (nmat==1 ? riemannDeriv.empty() : true), "Non-empty Riemann "
100 : : // "derivative vector for single material compflow" );
101 : :
102 [ + + ]: 491850 : for (const auto& s : bcconfig) { // for all bc sidesets
103 [ + - ]: 216810 : auto bc = bface.find(static_cast<int>(s));// faces for side set
104 [ + + ]: 216810 : if (bc != end(bface))
105 : : {
106 [ + + ]: 5302140 : for (const auto& f : bc->second)
107 : : {
108 [ - + ][ - - ]: 5203620 : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
[ - - ][ - - ]
109 : :
110 : 5203620 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
111 : :
112 [ + - ]: 5203620 : auto ng = tk::NGfa(ndofel[el]);
113 : :
114 : : // arrays for quadrature points
115 : 10407240 : std::array< std::vector< real >, 2 > coordgp;
116 : 10407240 : std::vector< real > wgp;
117 : :
118 [ + - ]: 5203620 : coordgp[0].resize( ng );
119 [ + - ]: 5203620 : coordgp[1].resize( ng );
120 [ + - ]: 5203620 : wgp.resize( ng );
121 : :
122 : : // get quadrature point weights and coordinates for triangle
123 [ + - ]: 5203620 : GaussQuadratureTri( ng, coordgp, wgp );
124 : :
125 : : // Extract the left element coordinates
126 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
127 : 15610860 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
128 : 15610860 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
129 : 15610860 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
130 : 46832580 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
131 : :
132 : : // Compute the determinant of Jacobian matrix
133 : : auto detT_l =
134 : 5203620 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
135 : :
136 : : // Extract the face coordinates
137 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
138 : 15610860 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
139 : 15610860 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
140 : 31221720 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
141 : :
142 : : std::array< real, 3 >
143 [ + - ][ + - ]: 5203620 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
[ + - ]
144 : :
145 : : // Gaussian quadrature
146 [ + + ]: 14579175 : for (std::size_t igp=0; igp<ng; ++igp)
147 : : {
148 : : // Compute the coordinates of quadrature point at physical domain
149 [ + - ]: 9375555 : auto gp = eval_gp( igp, coordfa, coordgp );
150 : :
151 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
152 : : // basis functions and solutions by default. This implies that P0P1 is
153 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until
154 : : // we have rdofel, which is needed to distinguish between ndofs and
155 : : // rdofs per element for pDG.
156 : : std::size_t dof_el;
157 [ + + ]: 9375555 : if (rdof > ndof)
158 : : {
159 : 537300 : dof_el = rdof;
160 : : }
161 : : else
162 : : {
163 : 8838255 : dof_el = ndofel[el];
164 : : }
165 : :
166 : : // For multi-material p-adaptive simulation, when dofel = 1, p0p1 is applied and ndof
167 : : // for solution evaluation should be 4
168 [ + + ][ + + ]: 9375555 : if(ncomp > 5 && dof_el == 1 && pref)
[ - + ]
169 : 0 : dof_el = 4;
170 : :
171 : : std::array< tk::real, 3> ref_gp_l{
172 : 9375555 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
173 : 9375555 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
174 : 18751110 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
175 : :
176 : : //Compute the basis functions for the left element
177 [ + - ]: 18751110 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
178 : :
179 [ + - ]: 9375555 : auto wt = wgp[igp] * geoFace(f,0);
180 : :
181 : : // Compute the state variables at the left element
182 : : auto ugp = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim,
183 [ + - ]: 18751110 : rdof, nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
184 : :
185 [ - + ][ - - ]: 9375555 : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
186 : : "appended boundary state vector" );
187 : :
188 [ + - ]: 18751110 : auto var = state( ncomp, mat_blk, ugp, gp[0], gp[1], gp[2], t, fn );
189 : :
190 : : // Compute the numerical flux
191 [ + - ][ + - ]: 18751110 : auto fl = flux(mat_blk, fn, var, vel(ncomp, gp[0], gp[1], gp[2], t));
192 : :
193 : : // Code below commented until details about the form of these terms in
194 : : // the \alpha_k g_k equations are sorted out.
195 : : // // Add RHS inverse deformation terms if necessary
196 : : // if (haveSolid)
197 : : // solidTermsSurfInt( nmat, ndof, rdof, fn, el, er, solidx, geoElem, U,
198 : : // coordel_l, coordel_r, igp, coordgp, dt, fl );
199 : :
200 : : // Add the surface integration term to the rhs
201 [ + - ]: 9375555 : update_rhs_bc( ncomp, nmat, ndof, ndofel[el], wt, fn, el, fl,
202 : : B_l, R, riemannDeriv );
203 : : }
204 : : }
205 : : }
206 : : }
207 : 275040 : }
208 : :
209 : : void
210 : 9375555 : update_rhs_bc ( ncomp_t ncomp,
211 : : std::size_t nmat,
212 : : const std::size_t ndof,
213 : : const std::size_t ndof_l,
214 : : const tk::real wt,
215 : : const std::array< tk::real, 3 >& fn,
216 : : const std::size_t el,
217 : : const std::vector< tk::real >& fl,
218 : : const std::vector< tk::real >& B_l,
219 : : Fields& R,
220 : : std::vector< std::vector< tk::real > >& riemannDeriv )
221 : : // *****************************************************************************
222 : : // Update the rhs by adding the boundary surface integration term
223 : : //! \param[in] ncomp Number of scalar components in this PDE system
224 : : //! \param[in] nmat Number of materials in this PDE system
225 : : //! \param[in] ndof Maximum number of degrees of freedom
226 : : //! \param[in] ndof_l Number of degrees of freedom for the left element
227 : : //! \param[in] wt Weight of gauss quadrature point
228 : : //! \param[in] fn Face/Surface normal
229 : : //! \param[in] el Left element index
230 : : //! \param[in] fl Surface flux
231 : : //! \param[in] B_l Basis function for the left element
232 : : //! \param[in,out] R Right-hand side vector computed
233 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
234 : : //! computed from the Riemann solver for use in the non-conservative terms.
235 : : //! These derivatives are used only for multi-material hydro and unused for
236 : : //! single-material compflow and linear transport.
237 : : // *****************************************************************************
238 : : {
239 : : // following line commented until rdofel is made available.
240 : : //Assert( B_l.size() == ndof_l, "Size mismatch" );
241 : :
242 : : using inciter::newSolidsAccFn;
243 : :
244 : : const auto& solidx =
245 : 9375555 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
246 : :
247 [ + + ]: 64383510 : for (ncomp_t c=0; c<ncomp; ++c)
248 : : {
249 : 55007955 : auto mark = c*ndof;
250 : 55007955 : R(el, mark) -= wt * fl[c];
251 : :
252 [ + + ]: 55007955 : if(ndof_l > 1) //DG(P1)
253 : : {
254 : 31009770 : R(el, mark+1) -= wt * fl[c] * B_l[1];
255 : 31009770 : R(el, mark+2) -= wt * fl[c] * B_l[2];
256 : 31009770 : R(el, mark+3) -= wt * fl[c] * B_l[3];
257 : : }
258 : :
259 [ + + ]: 55007955 : if(ndof_l > 4) //DG(P2)
260 : : {
261 : 13727610 : R(el, mark+4) -= wt * fl[c] * B_l[4];
262 : 13727610 : R(el, mark+5) -= wt * fl[c] * B_l[5];
263 : 13727610 : R(el, mark+6) -= wt * fl[c] * B_l[6];
264 : 13727610 : R(el, mark+7) -= wt * fl[c] * B_l[7];
265 : 13727610 : R(el, mark+8) -= wt * fl[c] * B_l[8];
266 : 13727610 : R(el, mark+9) -= wt * fl[c] * B_l[9];
267 : : }
268 : : }
269 : :
270 : : // Prep for non-conservative terms in multimat
271 [ + + ]: 9375555 : if (fl.size() > ncomp)
272 : : {
273 : : // Gradients of partial pressures
274 [ + + ]: 6722220 : for (std::size_t k=0; k<nmat; ++k)
275 : : {
276 [ + + ]: 18885120 : for (std::size_t idir=0; idir<3; ++idir)
277 : 14163840 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
278 : : }
279 : :
280 : : // Divergence of velocity multiples basis fucntion( d(uB) / dx )
281 [ + + ]: 6366000 : for(std::size_t idof = 0; idof < ndof_l; idof++)
282 : 4365060 : riemannDeriv[3*nmat+idof][el] += wt * fl[ncomp+nmat] * B_l[idof];
283 : :
284 : : // Divergence of asigma: d(asig_ij)/dx_j
285 [ + + ]: 6722220 : for (std::size_t k=0; k<nmat; ++k)
286 [ - + ]: 4721280 : if (solidx[k] > 0)
287 : : {
288 : 0 : std::size_t mark = ncomp+nmat+1+3*(solidx[k]-1);
289 : :
290 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
291 : 0 : riemannDeriv[3*nmat+ndof+3*(solidx[k]-1)+i][el] -=
292 : 0 : wt * fl[mark+i];
293 : : }
294 : :
295 : : // Derivatives of g: d(g_il)/d(x_j)-d(g_ij)/d(x_l)
296 : : // for i=1,2,3; j=1,2,3; l=1,2,3. Total = 3x3x3 (per solid)
297 : 2000940 : std::size_t nsld = inciter::numSolids(nmat, solidx);
298 [ + + ]: 6722220 : for (std::size_t k=0; k<nmat; ++k)
299 [ - + ]: 4721280 : if (solidx[k] > 0)
300 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
301 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
302 [ - - ]: 0 : for (std::size_t l=0; l<3; ++l)
303 [ - - ]: 0 : if (j != l)
304 : : {
305 : 0 : std::size_t mark1 = ncomp+nmat+1+3*nsld+9*(solidx[k]-1)+3*i+l;
306 : 0 : std::size_t mark2 = ncomp+nmat+1+3*nsld+9*(solidx[k]-1)+3*i+j;
307 : 0 : riemannDeriv[3*nmat+ndof+3*nsld+newSolidsAccFn(k,i,j,l)][el] -=
308 : 0 : wt * ( fl[mark1] * fn[j] - fl[mark2] * fn[l]);
309 : : }
310 : : }
311 : 9375555 : }
312 : :
313 : : void
314 : 10944 : bndSurfIntFV(
315 : : std::size_t nmat,
316 : : const std::vector< inciter::EOS >& mat_blk,
317 : : const std::size_t rdof,
318 : : const std::vector< std::size_t >& bcconfig,
319 : : const inciter::FaceData& fd,
320 : : const Fields& geoFace,
321 : : const Fields& geoElem,
322 : : const std::vector< std::size_t >& inpoel,
323 : : const UnsMesh::Coords& coord,
324 : : real t,
325 : : const RiemannFluxFn& flux,
326 : : const VelFn& vel,
327 : : const StateFn& state,
328 : : const Fields& U,
329 : : const Fields& P,
330 : : const std::vector< int >& srcFlag,
331 : : Fields& R,
332 : : int intsharp )
333 : : // *****************************************************************************
334 : : //! Compute boundary surface flux integrals for a given boundary type for FV
335 : : //! \details This function computes contributions from surface integrals along
336 : : //! all faces for a particular boundary condition type, configured by the state
337 : : //! function
338 : : //! \param[in] nmat Number of materials in this PDE system
339 : : //! \param[in] mat_blk EOS material block
340 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
341 : : //! \param[in] bcconfig BC configuration vector for multiple side sets
342 : : //! \param[in] fd Face connectivity and boundary conditions object
343 : : //! \param[in] geoFace Face geometry array
344 : : //! \param[in] geoElem Element geometry array
345 : : //! \param[in] inpoel Element-node connectivity
346 : : //! \param[in] coord Array of nodal coordinates
347 : : //! \param[in] t Physical time
348 : : //! \param[in] flux Riemann flux function to use
349 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
350 : : //! \param[in] state Function to evaluate the left and right solution state at
351 : : //! boundaries
352 : : //! \param[in] U Solution vector at recent time step
353 : : //! \param[in] P Vector of primitives at recent time step
354 : : //! \param[in] srcFlag Whether the energy source was added
355 : : //! \param[in,out] R Right-hand side vector computed
356 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
357 : : //! default 0, so that it is unused for single-material and transport.
358 : : // *****************************************************************************
359 : : {
360 : 10944 : const auto& bface = fd.Bface();
361 : 10944 : const auto& esuf = fd.Esuf();
362 : :
363 : 10944 : const auto& cx = coord[0];
364 : 10944 : const auto& cy = coord[1];
365 : 10944 : const auto& cz = coord[2];
366 : :
367 : 10944 : auto ncomp = U.nprop()/rdof;
368 : 10944 : auto nprim = P.nprop()/rdof;
369 : :
370 [ + + ]: 19152 : for (const auto& s : bcconfig) { // for all bc sidesets
371 [ + - ]: 8208 : auto bc = bface.find(static_cast<int>(s));// faces for side set
372 [ + + ]: 8208 : if (bc != end(bface))
373 : : {
374 [ + + ]: 144296 : for (const auto& f : bc->second)
375 : : {
376 [ - + ][ - - ]: 140816 : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
[ - - ][ - - ]
377 : :
378 : 140816 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
379 : :
380 : : // Extract the left element coordinates
381 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
382 : 422448 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
383 : 422448 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
384 : 422448 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
385 : 1267344 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
386 : :
387 : : // Compute the determinant of Jacobian matrix
388 : : auto detT_l =
389 : 140816 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
390 : :
391 : : // face normal
392 : : std::array< real, 3 >
393 [ + - ][ + - ]: 140816 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
[ + - ]
394 : :
395 : : // face centroid
396 : : std::array< real, 3 >
397 [ + - ][ + - ]: 140816 : gp{{ geoFace(f,4), geoFace(f,5), geoFace(f,6) }};
[ + - ]
398 : :
399 : : std::array< tk::real, 3> ref_gp_l{
400 : 140816 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
401 : 140816 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
402 : 281632 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
403 : :
404 : : //Compute the basis functions for the left element
405 [ + - ]: 281632 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
406 : :
407 : : // Compute the state variables at the left element
408 : : auto ugp = evalFVSol(mat_blk, intsharp, ncomp, nprim,
409 : : rdof, nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P,
410 [ + - ]: 281632 : srcFlag[el]);
411 : :
412 [ - + ][ - - ]: 140816 : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
413 : : "appended boundary state vector" );
414 : :
415 [ + - ]: 281632 : auto var = state( ncomp, mat_blk, ugp, gp[0], gp[1], gp[2], t, fn );
416 : :
417 : : // Compute the numerical flux
418 [ + - ][ + - ]: 281632 : auto fl = flux( mat_blk, fn, var, vel(ncomp, gp[0], gp[1], gp[2], t) );
419 : :
420 : : // compute non-conservative terms
421 [ + - ]: 281632 : std::vector< tk::real > var_riemann(nmat+1, 0.0);
422 [ + + ]: 592640 : for (std::size_t k=0; k<nmat+1; ++k) var_riemann[k] = fl[ncomp+k];
423 : :
424 [ + - ]: 281632 : auto ncf_l = nonConservativeIntFV(nmat, rdof, el, fn, U, P, var_riemann);
425 : :
426 : : // Add the surface integration term to the rhs
427 [ + + ]: 1496288 : for (ncomp_t c=0; c<ncomp; ++c)
428 : : {
429 [ + - ][ + - ]: 1355472 : R(el, c) -= geoFace(f,0) * (fl[c] - ncf_l[c]);
430 : : }
431 : : }
432 : : }
433 : : }
434 : 10944 : }
435 : :
436 : : void
437 : 0 : bndSurfIntViscousFV(
438 : : std::size_t nmat,
439 : : const std::vector< inciter::EOS >& mat_blk,
440 : : const std::size_t rdof,
441 : : const std::vector< std::size_t >& bcconfig,
442 : : const inciter::FaceData& fd,
443 : : const Fields& geoFace,
444 : : const Fields& geoElem,
445 : : const std::vector< std::size_t >& inpoel,
446 : : const UnsMesh::Coords& coord,
447 : : real t,
448 : : const StateFn& state,
449 : : const StateFn& gradFn,
450 : : const Fields& U,
451 : : const Fields& P,
452 : : const Fields& T,
453 : : const std::vector< int >& srcFlag,
454 : : Fields& R,
455 : : int intsharp )
456 : : // *****************************************************************************
457 : : //! Compute boundary surface flux integrals for a given boundary type for FV
458 : : //! \details This function computes contributions from surface integrals along
459 : : //! all faces for a particular boundary condition type, configured by the state
460 : : //! function
461 : : //! \param[in] nmat Number of materials in this PDE system
462 : : //! \param[in] mat_blk EOS material block
463 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
464 : : //! \param[in] bcconfig BC configuration vector for multiple side sets
465 : : //! \param[in] fd Face connectivity and boundary conditions object
466 : : //! \param[in] geoFace Face geometry array
467 : : //! \param[in] geoElem Element geometry array
468 : : //! \param[in] inpoel Element-node connectivity
469 : : //! \param[in] coord Array of nodal coordinates
470 : : //! \param[in] t Physical time
471 : : //! \param[in] state Function to evaluate the left and right solution state at
472 : : //! boundaries
473 : : //! \param[in] gradFn Function to evaluate the left and right solution gradients
474 : : //! at boundaries
475 : : //! \param[in] U Solution vector at recent time step
476 : : //! \param[in] P Vector of primitives at recent time step
477 : : //! \param[in] T Vector of temperatures at recent time step
478 : : //! \param[in] srcFlag Whether the energy source was added
479 : : //! \param[in,out] R Right-hand side vector computed
480 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
481 : : //! default 0, so that it is unused for single-material and transport.
482 : : // *****************************************************************************
483 : : {
484 : : using inciter::velocityDofIdx;
485 : :
486 : 0 : const auto& bface = fd.Bface();
487 : 0 : const auto& esuf = fd.Esuf();
488 : :
489 : 0 : const auto& cx = coord[0];
490 : 0 : const auto& cy = coord[1];
491 : 0 : const auto& cz = coord[2];
492 : :
493 : 0 : auto ncomp = U.nprop()/rdof;
494 : 0 : auto nprim = P.nprop()/rdof;
495 : :
496 [ - - ]: 0 : for (const auto& s : bcconfig) { // for all bc sidesets
497 [ - - ]: 0 : auto bc = bface.find(static_cast<int>(s));// faces for side set
498 [ - - ]: 0 : if (bc != end(bface))
499 : : {
500 [ - - ]: 0 : for (const auto& f : bc->second)
501 : : {
502 [ - - ][ - - ]: 0 : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
[ - - ][ - - ]
503 : :
504 : 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
505 : :
506 : : // Extract the left element coordinates
507 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
508 : 0 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
509 : 0 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
510 : 0 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
511 : 0 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
512 : :
513 : : // Compute the determinant of Jacobian matrix
514 : : auto detT_l =
515 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
516 : :
517 : : // face normal
518 : : std::array< real, 3 >
519 [ - - ][ - - ]: 0 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
[ - - ]
520 : :
521 : : // face centroid
522 : : std::array< real, 3 >
523 [ - - ][ - - ]: 0 : gp{{ geoFace(f,4), geoFace(f,5), geoFace(f,6) }};
[ - - ]
524 : :
525 : : std::array< tk::real, 3> ref_gp_l{
526 : 0 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
527 : 0 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
528 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
529 : :
530 : : //Compute the basis functions for the left element
531 [ - - ]: 0 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
532 : :
533 : : // Compute the state variables at the left element
534 : : auto ugp = evalFVSol(mat_blk, intsharp, ncomp, nprim,
535 : : rdof, nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P,
536 [ - - ]: 0 : srcFlag[el]);
537 : :
538 [ - - ][ - - ]: 0 : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
539 : : "appended boundary state vector" );
540 : :
541 [ - - ]: 0 : auto var = state( ncomp, mat_blk, ugp, gp[0], gp[1], gp[2], t, fn );
542 : :
543 : : // cell averaged state for computing the diffusive flux
544 [ - - ]: 0 : std::vector< tk::real > Bcc(rdof, 0.0);
545 : 0 : Bcc[0] = 1.0;
546 : :
547 : : auto ucc = evalFVSol(mat_blk, 0, ncomp, nprim, rdof,
548 : : nmat, el, inpoel, coord, geoElem, {{0.25, 0.25, 0.25}}, Bcc, U, P,
549 [ - - ]: 0 : srcFlag[el]);
550 : :
551 [ - - ][ - - ]: 0 : Assert( ucc.size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
552 : : "appended cell-averaged state vector" );
553 : :
554 : : // Cell centroids- [0]: left cell, [1]: ghost cell
555 : : // The ghost-cell is a 'reflection' of the boundary cell about the
556 : : // boundary-face. i.e. the vector pointing from the internal-cell
557 : : // centroid to the ghost-cell centroid is normal to the face (aligned
558 : : // with the face-normal), and has length 2*d. d is the distance between
559 : : // the internal-cell centroid and the boundary-face. Based on this
560 : : // information, the centroid of the ghost-cell can be computed using
561 : : // vector algebra.
562 : : std::array< std::array< tk::real, 3 >, 2 > centroids;
563 [ - - ][ - - ]: 0 : centroids[0] = {{geoElem(el,1), geoElem(el,2), geoElem(el,3)}};
[ - - ]
564 : 0 : tk::real d = std::abs( tk::dot(fn,centroids[0]) + tk::dot(fn,gp) ) /
565 : 0 : std::sqrt(tk::dot(fn,fn));
566 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
567 : 0 : centroids[1][i] = centroids[0][i] + 2.0*d*fn[i];
568 : :
569 : : // Get BC for cell-averaged state
570 : : auto varcc = state( ncomp, mat_blk, ucc,
571 [ - - ]: 0 : centroids[1][0], centroids[1][1], centroids[1][2], t, fn );
572 : : // Hard-coded temperature BC- only works for adiabatic 'noslipwall',
573 : : // 'symmetry', 'extrapolate' bc
574 : : std::array< std::vector< real >, 2 > Tcc{{
575 [ - - ][ - - ]: 0 : std::vector<real>(nmat), std::vector<real>(nmat) }};
576 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
577 [ - - ]: 0 : Tcc[0][k] = T(el, k*rdof);
578 : 0 : Tcc[1][k] = Tcc[0][k]; // actual bc
579 : : }
580 : :
581 : : // Numerical viscous flux
582 : : // ---------------------------------------------------------------------
583 : :
584 : : // 1. Get spatial gradient from Dubiner dofs
585 : 0 : auto jacInv_l = tk::inverseJacobian( coordel_l[0], coordel_l[1],
586 : 0 : coordel_l[2], coordel_l[3] );
587 [ - - ]: 0 : auto dBdx_l = tk::eval_dBdx_p1( rdof, jacInv_l );
588 : :
589 [ - - ]: 0 : std::vector< real > dudx_l(9,0.0);
590 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
591 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
592 : 0 : dudx_l[3*i+j] =
593 [ - - ]: 0 : dBdx_l[j][1] * P(el, velocityDofIdx(nmat,i,rdof,1))
594 [ - - ]: 0 : + dBdx_l[j][2] * P(el, velocityDofIdx(nmat,i,rdof,2))
595 [ - - ]: 0 : + dBdx_l[j][3] * P(el, velocityDofIdx(nmat,i,rdof,3));
596 : :
597 : : // 2. Average du_i/dx_j
598 [ - - ]: 0 : auto grad = gradFn( 3, mat_blk, dudx_l, gp[0], gp[2], gp[2], t, fn );
599 : : std::array< std::array< tk::real, 3 >, 3 > dudx;
600 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
601 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
602 : 0 : dudx[i][j] = 0.5 * (grad[0][3*i+j] + grad[1][3*i+j]);
603 : :
604 : : // 3. average dT/dx_j
605 : : std::vector< std::array< real, 3 > > dTdx(nmat,
606 [ - - ]: 0 : std::array< real, 3 >{{0, 0, 0}});
607 : :
608 : : // 4. Compute flux
609 : : auto fl = modifiedGradientViscousFlux(nmat, ncomp, fn, centroids, var,
610 [ - - ]: 0 : varcc, Tcc, dudx, dTdx);
611 : :
612 : : // Add the surface integration term to the rhs
613 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
614 : : {
615 [ - - ][ - - ]: 0 : R(el, c) += geoFace(f,0) * fl[c];
616 : : }
617 : : }
618 : : }
619 : : }
620 : 0 : }
621 : :
622 : : } // tk::
|
