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 : 434880 : bndSurfInt( std::size_t nmat,
35 : : const std::vector< inciter::EOS >& mat_blk,
36 : : const std::size_t ndof,
37 : : const std::size_t rdof,
38 : : const std::vector< std::size_t >& bcconfig,
39 : : const inciter::FaceData& fd,
40 : : const Fields& geoFace,
41 : : const Fields& geoElem,
42 : : const std::vector< std::size_t >& inpoel,
43 : : const UnsMesh::Coords& coord,
44 : : real t,
45 : : const RiemannFluxFn& flux,
46 : : const VelFn& vel,
47 : : const StateFn& state,
48 : : const Fields& U,
49 : : const Fields& P,
50 : : const std::vector< std::size_t >& ndofel,
51 : : Fields& R,
52 : : std::vector< std::vector< tk::real > >&,
53 : : std::vector< std::vector< tk::real > >&,
54 : : std::vector< std::vector< tk::real > >& riemannDeriv,
55 : : int intsharp )
56 : : // *****************************************************************************
57 : : //! Compute boundary surface flux integrals for a given boundary type for DG
58 : : //! \details This function computes contributions from surface integrals along
59 : : //! all faces for a particular boundary condition type, configured by the state
60 : : //! function
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] vriem Vector of the riemann velocity
81 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
82 : : //! is available
83 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
84 : : //! computed from the Riemann solver for use in the non-conservative terms.
85 : : //! These derivatives are used only for multi-material hydro and unused for
86 : : //! single-material compflow and linear transport.
87 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
88 : : //! default 0, so that it is unused for single-material and transport.
89 : : // *****************************************************************************
90 : : {
91 : 434880 : const auto& bface = fd.Bface();
92 : 434880 : const auto& esuf = fd.Esuf();
93 : 434880 : const auto& inpofa = fd.Inpofa();
94 : :
95 : 434880 : const auto& cx = coord[0];
96 : 434880 : const auto& cy = coord[1];
97 : 434880 : const auto& cz = coord[2];
98 : :
99 : 434880 : auto ncomp = U.nprop()/rdof;
100 : 434880 : auto nprim = P.nprop()/rdof;
101 : :
102 : : //Assert( (nmat==1 ? riemannDeriv.empty() : true), "Non-empty Riemann "
103 : : // "derivative vector for single material compflow" );
104 : :
105 [ + + ]: 734880 : for (const auto& s : bcconfig) { // for all bc sidesets
106 [ + - ]: 300000 : auto bc = bface.find(static_cast<int>(s));// faces for side set
107 [ + + ]: 300000 : if (bc != end(bface))
108 : : {
109 [ + + ]: 10738830 : for (const auto& f : bc->second)
110 : : {
111 [ - + ][ - - ]: 10584120 : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
[ - - ][ - - ]
112 : :
113 : 10584120 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
114 : :
115 [ + - ]: 10584120 : auto ng = tk::NGfa(ndofel[el]);
116 : :
117 : : // arrays for quadrature points
118 : 21168240 : std::array< std::vector< real >, 2 > coordgp;
119 : 21168240 : std::vector< real > wgp;
120 : :
121 [ + - ]: 10584120 : coordgp[0].resize( ng );
122 [ + - ]: 10584120 : coordgp[1].resize( ng );
123 [ + - ]: 10584120 : wgp.resize( ng );
124 : :
125 : : // get quadrature point weights and coordinates for triangle
126 [ + - ]: 10584120 : GaussQuadratureTri( ng, coordgp, wgp );
127 : :
128 : : // Extract the left element coordinates
129 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
130 : 31752360 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
131 : 31752360 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
132 : 31752360 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
133 : 95257080 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
134 : :
135 : : // Compute the determinant of Jacobian matrix
136 : : auto detT_l =
137 : 10584120 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
138 : :
139 : : // Extract the face coordinates
140 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
141 : 31752360 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
142 : 31752360 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
143 : 63504720 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
144 : :
145 : : std::array< real, 3 >
146 [ + - ][ + - ]: 10584120 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
[ + - ]
147 : :
148 : : // Gaussian quadrature
149 [ + + ]: 30735675 : for (std::size_t igp=0; igp<ng; ++igp)
150 : : {
151 : : // Compute the coordinates of quadrature point at physical domain
152 [ + - ]: 20151555 : auto gp = eval_gp( igp, coordfa, coordgp );
153 : :
154 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
155 : : // basis functions and solutions by default. This implies that P0P1 is
156 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until
157 : : // we have rdofel, which is needed to distinguish between ndofs and
158 : : // rdofs per element for pDG.
159 : : std::size_t dof_el;
160 [ + + ]: 20151555 : if (rdof > ndof)
161 : : {
162 : 1596180 : dof_el = rdof;
163 : : }
164 : : else
165 : : {
166 : 18555375 : dof_el = ndofel[el];
167 : : }
168 : :
169 : : std::array< tk::real, 3> ref_gp_l{
170 : 20151555 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
171 : 20151555 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
172 : 40303110 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
173 : :
174 : : //Compute the basis functions for the left element
175 [ + - ]: 40303110 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
176 : :
177 [ + - ]: 20151555 : auto wt = wgp[igp] * geoFace(f,0);
178 : :
179 : : // Compute the state variables at the left element
180 : : auto ugp = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim,
181 [ + - ]: 40303110 : rdof, nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
182 : :
183 [ - + ][ - - ]: 20151555 : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
184 : : "appended boundary state vector" );
185 : :
186 [ + - ]: 40303110 : auto var = state( ncomp, mat_blk, ugp, gp[0], gp[1], gp[2], t, fn );
187 : :
188 : : // Compute the numerical flux
189 [ + - ][ + - ]: 40303110 : auto fl = flux(mat_blk, fn, var, vel(ncomp, gp[0], gp[1], gp[2], t));
190 : :
191 : : // Code below commented until details about the form of these terms in
192 : : // the \alpha_k g_k equations are sorted out.
193 : : // // Add RHS inverse deformation terms if necessary
194 : : // if (haveSolid)
195 : : // solidTermsSurfInt( nmat, ndof, rdof, fn, el, er, solidx, geoElem, U,
196 : : // coordel_l, coordel_r, igp, coordgp, dt, fl );
197 : :
198 : : // Add the surface integration term to the rhs
199 [ + - ]: 20151555 : update_rhs_bc( ncomp, nmat, ndof, ndofel[el], wt, fn, el, fl,
200 : : B_l, R, riemannDeriv );
201 : : }
202 : : }
203 : : }
204 : : }
205 : 434880 : }
206 : :
207 : : void
208 : 20151555 : update_rhs_bc ( ncomp_t ncomp,
209 : : std::size_t nmat,
210 : : const std::size_t ndof,
211 : : const std::size_t ndof_l,
212 : : const tk::real wt,
213 : : const std::array< tk::real, 3 >& fn,
214 : : const std::size_t el,
215 : : const std::vector< tk::real >& fl,
216 : : const std::vector< tk::real >& B_l,
217 : : Fields& R,
218 : : std::vector< std::vector< tk::real > >& riemannDeriv )
219 : : // *****************************************************************************
220 : : // Update the rhs by adding the boundary surface integration term
221 : : //! \param[in] ncomp Number of scalar components in this PDE system
222 : : //! \param[in] nmat Number of materials in this PDE system
223 : : //! \param[in] ndof Maximum number of degrees of freedom
224 : : //! \param[in] ndof_l Number of degrees of freedom for the left element
225 : : //! \param[in] wt Weight of gauss quadrature point
226 : : //! \param[in] fn Face/Surface normal
227 : : //! \param[in] el Left element index
228 : : //! \param[in] fl Surface flux
229 : : //! \param[in] B_l Basis function for the left element
230 : : //! \param[in,out] R Right-hand side vector computed
231 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
232 : : //! computed from the Riemann solver for use in the non-conservative terms.
233 : : //! These derivatives are used only for multi-material hydro and unused for
234 : : //! single-material compflow and linear transport.
235 : : // *****************************************************************************
236 : : {
237 : : // following line commented until rdofel is made available.
238 : : //Assert( B_l.size() == ndof_l, "Size mismatch" );
239 : :
240 : : using inciter::newSolidsAccFn;
241 : :
242 : : const auto& solidx =
243 : 20151555 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
244 : :
245 [ + + ]: 84126690 : for (ncomp_t c=0; c<ncomp; ++c)
246 : : {
247 : 63975135 : auto mark = c*ndof;
248 : 63975135 : R(el, mark) -= wt * fl[c];
249 : :
250 [ + + ]: 63975135 : if(ndof_l > 1) //DG(P1)
251 : : {
252 : 38563470 : R(el, mark+1) -= wt * fl[c] * B_l[1];
253 : 38563470 : R(el, mark+2) -= wt * fl[c] * B_l[2];
254 : 38563470 : R(el, mark+3) -= wt * fl[c] * B_l[3];
255 : : }
256 : :
257 [ + + ]: 63975135 : if(ndof_l > 4) //DG(P2)
258 : : {
259 : 15885810 : R(el, mark+4) -= wt * fl[c] * B_l[4];
260 : 15885810 : R(el, mark+5) -= wt * fl[c] * B_l[5];
261 : 15885810 : R(el, mark+6) -= wt * fl[c] * B_l[6];
262 : 15885810 : R(el, mark+7) -= wt * fl[c] * B_l[7];
263 : 15885810 : R(el, mark+8) -= wt * fl[c] * B_l[8];
264 : 15885810 : R(el, mark+9) -= wt * fl[c] * B_l[9];
265 : : }
266 : : }
267 : :
268 : : // Prep for non-conservative terms in multimat
269 [ + + ]: 20151555 : if (fl.size() > ncomp)
270 : : {
271 : : // Gradients of partial pressures
272 [ + + ]: 6434760 : for (std::size_t k=0; k<nmat; ++k)
273 : : {
274 [ + + ]: 18118560 : for (std::size_t idir=0; idir<3; ++idir)
275 : 13588920 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
276 : : }
277 : :
278 : : // Divergence of velocity multiples basis fucntion( d(uB) / dx )
279 [ + + ]: 6174360 : for(std::size_t idof = 0; idof < ndof; idof++)
280 : 4269240 : riemannDeriv[3*nmat+idof][el] += wt * fl[ncomp+nmat] * B_l[idof];
281 : :
282 : : // Divergence of asigma: d(asig_ij)/dx_j
283 [ + + ]: 6434760 : for (std::size_t k=0; k<nmat; ++k)
284 [ + + ]: 4529640 : if (solidx[k] > 0)
285 : : {
286 : 55080 : std::size_t mark = ncomp+nmat+1+3*(solidx[k]-1);
287 : :
288 [ + + ]: 220320 : for (std::size_t i=0; i<3; ++i)
289 : 165240 : riemannDeriv[3*nmat+ndof+3*(solidx[k]-1)+i][el] -=
290 : 165240 : wt * fl[mark+i];
291 : : }
292 : : }
293 : 20151555 : }
294 : :
295 : : void
296 : 9708 : bndSurfIntFV(
297 : : std::size_t nmat,
298 : : const std::vector< inciter::EOS >& mat_blk,
299 : : const std::size_t rdof,
300 : : const std::vector< std::size_t >& bcconfig,
301 : : const inciter::FaceData& fd,
302 : : const Fields& geoFace,
303 : : const Fields& geoElem,
304 : : const std::vector< std::size_t >& inpoel,
305 : : const UnsMesh::Coords& coord,
306 : : real t,
307 : : const RiemannFluxFn& flux,
308 : : const VelFn& vel,
309 : : const StateFn& state,
310 : : const Fields& U,
311 : : const Fields& P,
312 : : const std::vector< int >& srcFlag,
313 : : Fields& R,
314 : : int intsharp )
315 : : // *****************************************************************************
316 : : //! Compute boundary surface flux integrals for a given boundary type for FV
317 : : //! \details This function computes contributions from surface integrals along
318 : : //! all faces for a particular boundary condition type, configured by the state
319 : : //! function
320 : : //! \param[in] nmat Number of materials in this PDE system
321 : : //! \param[in] mat_blk EOS material block
322 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
323 : : //! \param[in] bcconfig BC configuration vector for multiple side sets
324 : : //! \param[in] fd Face connectivity and boundary conditions object
325 : : //! \param[in] geoFace Face geometry array
326 : : //! \param[in] geoElem Element geometry array
327 : : //! \param[in] inpoel Element-node connectivity
328 : : //! \param[in] coord Array of nodal coordinates
329 : : //! \param[in] t Physical time
330 : : //! \param[in] flux Riemann flux function to use
331 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
332 : : //! \param[in] state Function to evaluate the left and right solution state at
333 : : //! boundaries
334 : : //! \param[in] U Solution vector at recent time step
335 : : //! \param[in] P Vector of primitives at recent time step
336 : : //! \param[in] srcFlag Whether the energy source was added
337 : : //! \param[in,out] R Right-hand side vector computed
338 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
339 : : //! default 0, so that it is unused for single-material and transport.
340 : : // *****************************************************************************
341 : : {
342 : 9708 : const auto& bface = fd.Bface();
343 : 9708 : const auto& esuf = fd.Esuf();
344 : :
345 : 9708 : const auto& cx = coord[0];
346 : 9708 : const auto& cy = coord[1];
347 : 9708 : const auto& cz = coord[2];
348 : :
349 : 9708 : auto ncomp = U.nprop()/rdof;
350 : 9708 : auto nprim = P.nprop()/rdof;
351 : :
352 [ + + ]: 17816 : for (const auto& s : bcconfig) { // for all bc sidesets
353 [ + - ]: 8108 : auto bc = bface.find(static_cast<int>(s));// faces for side set
354 [ + + ]: 8108 : if (bc != end(bface))
355 : : {
356 [ + + ]: 130316 : for (const auto& f : bc->second)
357 : : {
358 [ - + ][ - - ]: 126976 : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
[ - - ][ - - ]
359 : :
360 : 126976 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
361 : :
362 : : // Extract the left element coordinates
363 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
364 : 380928 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
365 : 380928 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
366 : 380928 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
367 : 1142784 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
368 : :
369 : : // Compute the determinant of Jacobian matrix
370 : : auto detT_l =
371 : 126976 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
372 : :
373 : : // face normal
374 : : std::array< real, 3 >
375 [ + - ][ + - ]: 126976 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
[ + - ]
376 : :
377 : : // face centroid
378 : : std::array< real, 3 >
379 [ + - ][ + - ]: 126976 : gp{{ geoFace(f,4), geoFace(f,5), geoFace(f,6) }};
[ + - ]
380 : :
381 : : std::array< tk::real, 3> ref_gp_l{
382 : 126976 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
383 : 126976 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
384 : 253952 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
385 : :
386 : : //Compute the basis functions for the left element
387 [ + - ]: 253952 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
388 : :
389 : : // Compute the state variables at the left element
390 : : auto ugp = evalFVSol(mat_blk, intsharp, ncomp, nprim,
391 : : rdof, nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P,
392 [ + - ]: 253952 : srcFlag[el]);
393 : :
394 [ - + ][ - - ]: 126976 : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
395 : : "appended boundary state vector" );
396 : :
397 [ + - ]: 253952 : auto var = state( ncomp, mat_blk, ugp, gp[0], gp[1], gp[2], t, fn );
398 : :
399 : : // Compute the numerical flux
400 [ + - ][ + - ]: 253952 : auto fl = flux( mat_blk, fn, var, vel(ncomp, gp[0], gp[1], gp[2], t) );
401 : :
402 : : // compute non-conservative terms
403 [ + - ]: 253952 : std::vector< tk::real > var_riemann(nmat+1, 0.0);
404 [ + + ]: 537280 : for (std::size_t k=0; k<nmat+1; ++k) var_riemann[k] = fl[ncomp+k];
405 : :
406 [ + - ]: 253952 : auto ncf_l = nonConservativeIntFV(nmat, rdof, el, fn, U, P, var_riemann);
407 : :
408 : : // Add the surface integration term to the rhs
409 [ + + ]: 1357888 : for (ncomp_t c=0; c<ncomp; ++c)
410 : : {
411 [ + - ][ + - ]: 1230912 : R(el, c) -= geoFace(f,0) * (fl[c] - ncf_l[c]);
412 : : }
413 : : }
414 : : }
415 : : }
416 : 9708 : }
417 : :
418 : : } // tk::
|
