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 : 456435 : 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 : : const auto& bface = fd.Bface();
89 : : const auto& esuf = fd.Esuf();
90 : : const auto& inpofa = fd.Inpofa();
91 : :
92 : : const auto& cx = coord[0];
93 : : const auto& cy = coord[1];
94 : : const auto& cz = coord[2];
95 : :
96 : 456435 : auto ncomp = U.nprop()/rdof;
97 : 456435 : 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 [ + + ]: 733035 : for (const auto& s : bcconfig) { // for all bc sidesets
103 : 276600 : auto bc = bface.find(static_cast<int>(s));// faces for side set
104 [ + + ]: 276600 : if (bc != end(bface))
105 : : {
106 [ + + ]: 9306990 : for (const auto& f : bc->second)
107 : : {
108 : : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
109 : :
110 [ + - ]: 9167040 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
111 : :
112 [ + - ]: 9167040 : auto ng = tk::NGfa(ndofel[el]);
113 : :
114 : : // arrays for quadrature points
115 : : std::array< std::vector< real >, 2 > coordgp;
116 : : std::vector< real > wgp;
117 : :
118 [ + - ]: 9167040 : coordgp[0].resize( ng );
119 [ + - ]: 9167040 : coordgp[1].resize( ng );
120 [ + - ]: 9167040 : wgp.resize( ng );
121 : :
122 : : // get quadrature point weights and coordinates for triangle
123 [ + - ]: 9167040 : GaussQuadratureTri( ng, coordgp, wgp );
124 : :
125 : : // Extract the left element coordinates
126 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
127 : 9167040 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
128 : 9167040 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
129 : 9167040 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
130 : 9167040 : {{ 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 : 9167040 : 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 : 9167040 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
139 : 9167040 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
140 : 9167040 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
141 : :
142 : : std::array< real, 3 >
143 : 9167040 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
144 : :
145 : : // Gaussian quadrature
146 [ + + ]: 27901515 : for (std::size_t igp=0; igp<ng; ++igp)
147 : : {
148 : : // Compute the coordinates of quadrature point at physical domain
149 [ + - ]: 18734475 : 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 [ + + ]: 18734475 : if (rdof > ndof)
158 : : {
159 : : dof_el = rdof;
160 : : }
161 : : else
162 : : {
163 : 18555375 : 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 [ + + ][ - + ]: 18734475 : if(ncomp > 5 && dof_el == 1 && pref)
169 : : dof_el = 4;
170 : :
171 : : std::array< tk::real, 3> ref_gp_l{
172 : 18734475 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
173 : 18734475 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
174 : 18734475 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
175 : :
176 : : //Compute the basis functions for the left element
177 [ + - ]: 18734475 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
178 : :
179 [ + - ]: 18734475 : 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 [ + - ]: 18734475 : rdof, nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
184 : :
185 : : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
186 : : "appended boundary state vector" );
187 : :
188 [ - + ]: 18734475 : auto var = state( ncomp, mat_blk, ugp, gp[0], gp[1], gp[2], t, fn );
189 : :
190 : : // Compute the numerical flux
191 [ - + ][ - + ]: 37468950 : 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 [ + - ]: 18734475 : update_rhs_bc( ncomp, nmat, ndof, ndofel[el], wt, fn, el, fl,
202 : : B_l, R, riemannDeriv );
203 : : }
204 : : }
205 : : }
206 : : }
207 : 456435 : }
208 : :
209 : : void
210 : 18734475 : 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 : : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
246 : :
247 [ + + ]: 79048170 : for (ncomp_t c=0; c<ncomp; ++c)
248 : : {
249 : 60313695 : auto mark = c*ndof;
250 [ + + ]: 60313695 : R(el, mark) -= wt * fl[c];
251 : :
252 [ + + ]: 60313695 : if(ndof_l > 1) //DG(P1)
253 : : {
254 : 38563470 : R(el, mark+1) -= wt * fl[c] * B_l[1];
255 : 38563470 : R(el, mark+2) -= wt * fl[c] * B_l[2];
256 : 38563470 : R(el, mark+3) -= wt * fl[c] * B_l[3];
257 : : }
258 : :
259 [ + + ]: 60313695 : if(ndof_l > 4) //DG(P2)
260 : : {
261 : 15885810 : R(el, mark+4) -= wt * fl[c] * B_l[4];
262 : 15885810 : R(el, mark+5) -= wt * fl[c] * B_l[5];
263 : 15885810 : R(el, mark+6) -= wt * fl[c] * B_l[6];
264 : 15885810 : R(el, mark+7) -= wt * fl[c] * B_l[7];
265 : 15885810 : R(el, mark+8) -= wt * fl[c] * B_l[8];
266 : 15885810 : R(el, mark+9) -= wt * fl[c] * B_l[9];
267 : : }
268 : : }
269 : :
270 : : // Prep for non-conservative terms in multimat
271 [ + + ]: 18734475 : if (fl.size() > ncomp)
272 : : {
273 : : // Gradients of partial pressures
274 [ + + ]: 6137220 : for (std::size_t k=0; k<nmat; ++k)
275 : : {
276 [ + + ]: 17325120 : for (std::size_t idir=0; idir<3; ++idir)
277 : 12993840 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
278 : : }
279 : :
280 : : // Divergence of velocity multiples basis fucntion( d(uB) / dx )
281 [ + + ]: 5976000 : for(std::size_t idof = 0; idof < ndof_l; idof++)
282 : 4170060 : riemannDeriv[3*nmat+idof][el] += wt * fl[ncomp+nmat] * B_l[idof];
283 : :
284 : : // Divergence of asigma: d(asig_ij)/dx_j
285 [ + + ]: 6137220 : for (std::size_t k=0; k<nmat; ++k)
286 [ - + ]: 4331280 : 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 : 18734475 : }
296 : :
297 : : void
298 : 11676 : bndSurfIntFV(
299 : : std::size_t nmat,
300 : : const std::vector< inciter::EOS >& mat_blk,
301 : : const std::size_t rdof,
302 : : const std::vector< std::size_t >& bcconfig,
303 : : const inciter::FaceData& fd,
304 : : const Fields& geoFace,
305 : : const Fields& geoElem,
306 : : const std::vector< std::size_t >& inpoel,
307 : : const UnsMesh::Coords& coord,
308 : : real t,
309 : : const RiemannFluxFn& flux,
310 : : const VelFn& vel,
311 : : const StateFn& state,
312 : : const Fields& U,
313 : : const Fields& P,
314 : : const std::vector< int >& srcFlag,
315 : : Fields& R,
316 : : int intsharp )
317 : : // *****************************************************************************
318 : : //! Compute boundary surface flux integrals for a given boundary type for FV
319 : : //! \details This function computes contributions from surface integrals along
320 : : //! all faces for a particular boundary condition type, configured by the state
321 : : //! function
322 : : //! \param[in] nmat Number of materials in this PDE system
323 : : //! \param[in] mat_blk EOS material block
324 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
325 : : //! \param[in] bcconfig BC configuration vector for multiple side sets
326 : : //! \param[in] fd Face connectivity and boundary conditions object
327 : : //! \param[in] geoFace Face geometry array
328 : : //! \param[in] geoElem Element geometry array
329 : : //! \param[in] inpoel Element-node connectivity
330 : : //! \param[in] coord Array of nodal coordinates
331 : : //! \param[in] t Physical time
332 : : //! \param[in] flux Riemann flux function to use
333 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
334 : : //! \param[in] state Function to evaluate the left and right solution state at
335 : : //! boundaries
336 : : //! \param[in] U Solution vector at recent time step
337 : : //! \param[in] P Vector of primitives at recent time step
338 : : //! \param[in] srcFlag Whether the energy source was added
339 : : //! \param[in,out] R Right-hand side vector computed
340 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
341 : : //! default 0, so that it is unused for single-material and transport.
342 : : // *****************************************************************************
343 : : {
344 : : const auto& bface = fd.Bface();
345 : : const auto& esuf = fd.Esuf();
346 : :
347 : : const auto& cx = coord[0];
348 : : const auto& cy = coord[1];
349 : : const auto& cz = coord[2];
350 : :
351 : 11676 : auto ncomp = U.nprop()/rdof;
352 : 11676 : auto nprim = P.nprop()/rdof;
353 : :
354 [ + + ]: 20084 : for (const auto& s : bcconfig) { // for all bc sidesets
355 : 8408 : auto bc = bface.find(static_cast<int>(s));// faces for side set
356 [ + + ]: 8408 : if (bc != end(bface))
357 : : {
358 [ + + ]: 170416 : for (const auto& f : bc->second)
359 : : {
360 : : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
361 : :
362 [ + - ]: 166776 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
363 : :
364 : : // Extract the left element coordinates
365 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
366 : 166776 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
367 : 166776 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
368 : 166776 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
369 : 166776 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
370 : :
371 : : // Compute the determinant of Jacobian matrix
372 : : auto detT_l =
373 [ + - ]: 166776 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
374 : :
375 : : // face normal
376 : : std::array< real, 3 >
377 : 166776 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
378 : :
379 : : // face centroid
380 : : std::array< real, 3 >
381 : 166776 : gp{{ geoFace(f,4), geoFace(f,5), geoFace(f,6) }};
382 : :
383 : : std::array< tk::real, 3> ref_gp_l{
384 : 166776 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
385 : 166776 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
386 : 166776 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
387 : :
388 : : //Compute the basis functions for the left element
389 [ + - ]: 166776 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
390 : :
391 : : // Compute the state variables at the left element
392 : : auto ugp = evalFVSol(mat_blk, intsharp, ncomp, nprim,
393 : : rdof, nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P,
394 [ + - ]: 166776 : srcFlag[el]);
395 : :
396 : : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
397 : : "appended boundary state vector" );
398 : :
399 [ - + ]: 166776 : auto var = state( ncomp, mat_blk, ugp, gp[0], gp[1], gp[2], t, fn );
400 : :
401 : : // Compute the numerical flux
402 [ - + ][ - + ]: 333552 : auto fl = flux( mat_blk, fn, var, vel(ncomp, gp[0], gp[1], gp[2], t) );
[ - - ]
403 : :
404 : : // compute non-conservative terms
405 [ + - ][ - - ]: 166776 : std::vector< tk::real > var_riemann(nmat+1, 0.0);
406 [ + + ]: 696480 : for (std::size_t k=0; k<nmat+1; ++k) var_riemann[k] = fl[ncomp+k];
407 : :
408 [ + - ]: 166776 : auto ncf_l = nonConservativeIntFV(nmat, rdof, el, fn, U, P, var_riemann);
409 : :
410 : : // Add the surface integration term to the rhs
411 [ + + ]: 1755888 : for (ncomp_t c=0; c<ncomp; ++c)
412 : : {
413 : 1589112 : R(el, c) -= geoFace(f,0) * (fl[c] - ncf_l[c]);
414 : : }
415 : : }
416 : : }
417 : : }
418 : 11676 : }
419 : :
420 : : void
421 : 0 : bndSurfIntViscousFV(
422 : : std::size_t nmat,
423 : : const std::vector< inciter::EOS >& mat_blk,
424 : : const std::size_t rdof,
425 : : const std::vector< std::size_t >& bcconfig,
426 : : const inciter::FaceData& fd,
427 : : const Fields& geoFace,
428 : : const Fields& geoElem,
429 : : const std::vector< std::size_t >& inpoel,
430 : : const UnsMesh::Coords& coord,
431 : : real t,
432 : : const StateFn& state,
433 : : const StateFn& gradFn,
434 : : const Fields& U,
435 : : const Fields& P,
436 : : const Fields& T,
437 : : const std::vector< int >& srcFlag,
438 : : Fields& R,
439 : : int intsharp )
440 : : // *****************************************************************************
441 : : //! Compute boundary surface flux integrals for a given boundary type for FV
442 : : //! \details This function computes contributions from surface integrals along
443 : : //! all faces for a particular boundary condition type, configured by the state
444 : : //! function
445 : : //! \param[in] nmat Number of materials in this PDE system
446 : : //! \param[in] mat_blk EOS material block
447 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
448 : : //! \param[in] bcconfig BC configuration vector for multiple side sets
449 : : //! \param[in] fd Face connectivity and boundary conditions object
450 : : //! \param[in] geoFace Face geometry array
451 : : //! \param[in] geoElem Element geometry array
452 : : //! \param[in] inpoel Element-node connectivity
453 : : //! \param[in] coord Array of nodal coordinates
454 : : //! \param[in] t Physical time
455 : : //! \param[in] state Function to evaluate the left and right solution state at
456 : : //! boundaries
457 : : //! \param[in] gradFn Function to evaluate the left and right solution gradients
458 : : //! at boundaries
459 : : //! \param[in] U Solution vector at recent time step
460 : : //! \param[in] P Vector of primitives at recent time step
461 : : //! \param[in] T Vector of temperatures at recent time step
462 : : //! \param[in] srcFlag Whether the energy source was added
463 : : //! \param[in,out] R Right-hand side vector computed
464 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
465 : : //! default 0, so that it is unused for single-material and transport.
466 : : // *****************************************************************************
467 : : {
468 : : using inciter::velocityDofIdx;
469 : :
470 : : const auto& bface = fd.Bface();
471 : : const auto& esuf = fd.Esuf();
472 : :
473 : : const auto& cx = coord[0];
474 : : const auto& cy = coord[1];
475 : : const auto& cz = coord[2];
476 : :
477 : 0 : auto ncomp = U.nprop()/rdof;
478 : 0 : auto nprim = P.nprop()/rdof;
479 : :
480 [ - - ]: 0 : for (const auto& s : bcconfig) { // for all bc sidesets
481 : 0 : auto bc = bface.find(static_cast<int>(s));// faces for side set
482 [ - - ]: 0 : if (bc != end(bface))
483 : : {
484 [ - - ]: 0 : for (const auto& f : bc->second)
485 : : {
486 : : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
487 : :
488 [ - - ]: 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
489 : :
490 : : // Extract the left element coordinates
491 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
492 : 0 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
493 : 0 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
494 : 0 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
495 : 0 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
496 : :
497 : : // Compute the determinant of Jacobian matrix
498 : : auto detT_l =
499 [ - - ]: 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
500 : :
501 : : // face normal
502 : : std::array< real, 3 >
503 : 0 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
504 : :
505 : : // face centroid
506 : : std::array< real, 3 >
507 : 0 : gp{{ geoFace(f,4), geoFace(f,5), geoFace(f,6) }};
508 : :
509 : : std::array< tk::real, 3> ref_gp_l{
510 : 0 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
511 : 0 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
512 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
513 : :
514 : : //Compute the basis functions for the left element
515 [ - - ]: 0 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
516 : :
517 : : // Compute the state variables at the left element
518 : : auto ugp = evalFVSol(mat_blk, intsharp, ncomp, nprim,
519 : : rdof, nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P,
520 [ - - ]: 0 : srcFlag[el]);
521 : :
522 : : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
523 : : "appended boundary state vector" );
524 : :
525 [ - - ]: 0 : auto var = state( ncomp, mat_blk, ugp, gp[0], gp[1], gp[2], t, fn );
526 : :
527 : : // cell averaged state for computing the diffusive flux
528 [ - - ]: 0 : std::vector< tk::real > Bcc(rdof, 0.0);
529 : 0 : Bcc[0] = 1.0;
530 : :
531 : : auto ucc = evalFVSol(mat_blk, 0, ncomp, nprim, rdof,
532 : : nmat, el, inpoel, coord, geoElem, {{0.25, 0.25, 0.25}}, Bcc, U, P,
533 [ - - ][ - - ]: 0 : srcFlag[el]);
534 : :
535 : : Assert( ucc.size() == ncomp+nprim, "Incorrect size for "
536 : : "appended cell-averaged state vector" );
537 : :
538 : : // Cell centroids- [0]: left cell, [1]: ghost cell
539 : : // The ghost-cell is a 'reflection' of the boundary cell about the
540 : : // boundary-face. i.e. the vector pointing from the internal-cell
541 : : // centroid to the ghost-cell centroid is normal to the face (aligned
542 : : // with the face-normal), and has length 2*d. d is the distance between
543 : : // the internal-cell centroid and the boundary-face. Based on this
544 : : // information, the centroid of the ghost-cell can be computed using
545 : : // vector algebra.
546 : : std::array< std::array< tk::real, 3 >, 2 > centroids;
547 : 0 : centroids[0] = {{geoElem(el,1), geoElem(el,2), geoElem(el,3)}};
548 : 0 : tk::real d = std::abs( tk::dot(fn,centroids[0]) + tk::dot(fn,gp) ) /
549 : 0 : std::sqrt(tk::dot(fn,fn));
550 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
551 : 0 : centroids[1][i] = centroids[0][i] + 2.0*d*fn[i];
552 : :
553 : : // Get BC for cell-averaged state
554 : : auto varcc = state( ncomp, mat_blk, ucc,
555 [ - - ]: 0 : centroids[1][0], centroids[1][1], centroids[1][2], t, fn );
556 : : // Hard-coded temperature BC- only works for adiabatic 'noslipwall',
557 : : // 'symmetry', 'extrapolate' bc
558 : : std::array< std::vector< real >, 2 > Tcc{{
559 [ - - ][ - - ]: 0 : std::vector<real>(nmat), std::vector<real>(nmat) }};
560 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
561 : 0 : Tcc[0][k] = T(el, k*rdof);
562 : 0 : Tcc[1][k] = Tcc[0][k]; // actual bc
563 : : }
564 : :
565 : : // Numerical viscous flux
566 : : // ---------------------------------------------------------------------
567 : :
568 : : // 1. Get spatial gradient from Dubiner dofs
569 : : auto jacInv_l = tk::inverseJacobian( coordel_l[0], coordel_l[1],
570 : 0 : coordel_l[2], coordel_l[3] );
571 [ - - ]: 0 : auto dBdx_l = tk::eval_dBdx_p1( rdof, jacInv_l );
572 : :
573 [ - - ]: 0 : std::vector< real > dudx_l(9,0.0);
574 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
575 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
576 : 0 : dudx_l[3*i+j] =
577 : 0 : dBdx_l[j][1] * P(el, velocityDofIdx(nmat,i,rdof,1))
578 : 0 : + dBdx_l[j][2] * P(el, velocityDofIdx(nmat,i,rdof,2))
579 : 0 : + dBdx_l[j][3] * P(el, velocityDofIdx(nmat,i,rdof,3));
580 : :
581 : : // 2. Average du_i/dx_j
582 [ - - ]: 0 : auto grad = gradFn( 3, mat_blk, dudx_l, gp[0], gp[2], gp[2], t, fn );
583 : : std::array< std::array< tk::real, 3 >, 3 > dudx;
584 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
585 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
586 : 0 : dudx[i][j] = 0.5 * (grad[0][3*i+j] + grad[1][3*i+j]);
587 : :
588 : : // 3. average dT/dx_j
589 : : std::vector< std::array< real, 3 > > dTdx(nmat,
590 [ - - ]: 0 : std::array< real, 3 >{{0, 0, 0}});
591 : :
592 : : // 4. Compute flux
593 : : auto fl = modifiedGradientViscousFlux(nmat, ncomp, fn, centroids, var,
594 [ - - ]: 0 : varcc, Tcc, dudx, dTdx);
595 : :
596 : : // Add the surface integration term to the rhs
597 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
598 : : {
599 : 0 : R(el, c) += geoFace(f,0) * fl[c];
600 : : }
601 : : }
602 : : }
603 : : }
604 : 0 : }
605 : :
606 : : } // tk::
|
