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 : :
26 : : namespace tk {
27 : :
28 : : void
29 : 447480 : bndSurfInt( ncomp_t system,
30 : : std::size_t nmat,
31 : : ncomp_t offset,
32 : : const std::size_t ndof,
33 : : const std::size_t rdof,
34 : : const std::vector< bcconf_t >& bcconfig,
35 : : const inciter::FaceData& fd,
36 : : const Fields& geoFace,
37 : : const Fields& geoElem,
38 : : const std::vector< std::size_t >& inpoel,
39 : : const UnsMesh::Coords& coord,
40 : : real t,
41 : : const RiemannFluxFn& flux,
42 : : const VelFn& vel,
43 : : const StateFn& state,
44 : : const Fields& U,
45 : : const Fields& P,
46 : : const std::vector< std::size_t >& ndofel,
47 : : Fields& R,
48 : : std::vector< std::vector< tk::real > >& vriem,
49 : : std::vector< std::vector< tk::real > >& riemannLoc,
50 : : std::vector< std::vector< tk::real > >& riemannDeriv,
51 : : int intsharp )
52 : : // *****************************************************************************
53 : : //! Compute boundary surface flux integrals for a given boundary type for DG
54 : : //! \details This function computes contributions from surface integrals along
55 : : //! all faces for a particular boundary condition type, configured by the state
56 : : //! function
57 : : //! \param[in] system Equation system index
58 : : //! \param[in] nmat Number of materials in this PDE system
59 : : //! \param[in] offset Offset this PDE system operates from
60 : : //! \param[in] ndof Maximum number of degrees of freedom
61 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
62 : : //! \param[in] bcconfig BC configuration vector for multiple side sets
63 : : //! \param[in] fd Face connectivity and boundary conditions object
64 : : //! \param[in] geoFace Face geometry array
65 : : //! \param[in] geoElem Element geometry array
66 : : //! \param[in] inpoel Element-node connectivity
67 : : //! \param[in] coord Array of nodal coordinates
68 : : //! \param[in] t Physical time
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] state Function to evaluate the left and right solution state at
72 : : //! boundaries
73 : : //! \param[in] U Solution vector at recent time step
74 : : //! \param[in] P Vector of primitives at recent time step
75 : : //! \param[in] ndofel Vector of local number of degrees of freedom
76 : : //! \param[in,out] R Right-hand side vector computed
77 : : //! \param[in,out] vriem Vector of the riemann velocity
78 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
79 : : //! is available
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 : 447480 : auto ncomp = U.nprop()/rdof;
97 : 447480 : 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 [ + + ]: 708720 : for (const auto& s : bcconfig) { // for all bc sidesets
103 : : auto bc = bface.find( std::stoi(s) );// faces for side set
104 [ + + ]: 261240 : if (bc != end(bface))
105 : : {
106 [ + + ]: 14162370 : for (const auto& f : bc->second)
107 : : {
108 : : Assert( esuf[2*f+1] == -1, "outside boundary element not -1" );
109 : :
110 [ + - ]: 14013840 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
111 : :
112 [ + - ]: 14013840 : 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 [ + - ]: 14013840 : coordgp[0].resize( ng );
119 [ + - ]: 14013840 : coordgp[1].resize( ng );
120 [ + - ]: 14013840 : wgp.resize( ng );
121 : :
122 : : // get quadrature point weights and coordinates for triangle
123 [ + - ]: 14013840 : GaussQuadratureTri( ng, coordgp, wgp );
124 : :
125 : : // Extract the left element coordinates
126 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
127 : 14013840 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
128 : 14013840 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
129 : 14013840 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
130 : 14013840 : {{ 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 : 14013840 : 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 : 14013840 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
139 : 14013840 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
140 : 14013840 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
141 : :
142 : : std::array< real, 3 >
143 : 14013840 : fn{{ geoFace(f,1,0), geoFace(f,2,0), geoFace(f,3,0) }};
144 : :
145 : : // Gaussian quadrature
146 [ + + ]: 38780397 : for (std::size_t igp=0; igp<ng; ++igp)
147 : : {
148 : : // Compute the coordinates of quadrature point at physical domain
149 [ + - ]: 24766557 : 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 [ + + ]: 24766557 : if (rdof > ndof)
158 : : {
159 : : dof_el = rdof;
160 : : }
161 : : else
162 : : {
163 : 23209857 : dof_el = ndofel[el];
164 : : }
165 : :
166 : : std::array< tk::real, 3> ref_gp_l{
167 : 24766557 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
168 : 24766557 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
169 : 24766557 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
170 : :
171 : : //Compute the basis functions for the left element
172 [ + - ]: 24766557 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
173 : :
174 [ + - ]: 24766557 : auto wt = wgp[igp] * geoFace(f,0,0);
175 : :
176 : : // Compute the state variables at the left element
177 : : auto ugp = evalPolynomialSol(system, offset, intsharp, ncomp, nprim,
178 [ + - ]: 24766557 : rdof, nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
179 : :
180 : : Assert( ugp.size() == ncomp+nprim, "Incorrect size for "
181 : : "appended boundary state vector" );
182 : :
183 [ - + ]: 24766557 : auto var = state( system, ncomp, ugp, gp[0], gp[1], gp[2], t, fn );
184 : :
185 : : // Compute the numerical flux
186 [ - + ][ - + ]: 49533114 : auto fl = flux( fn, var, vel( system, ncomp, gp[0], gp[1], gp[2], t ) );
[ - - ]
187 : :
188 : : // Add the surface integration term to the rhs
189 [ + - ]: 24766557 : update_rhs_bc( ncomp, nmat, offset, ndof, ndofel[el], wt, fn, el, fl,
190 : : B_l, R, riemannDeriv );
191 : :
192 : : // Store the riemann velocity and coordinates data of quadrature point
193 : : // used for velocity reconstruction if MulMat scheme is selected
194 [ + + ]: 24766557 : if (nmat > 1 && ndof > 1)
195 [ + - ]: 716400 : tk::evaluRiemann( ncomp, esuf[2*f], esuf[2*f+1], nmat, fl, fn, gp,
196 : : var, vriem, riemannLoc );
197 : : }
198 : : }
199 : : }
200 : : }
201 : 447480 : }
202 : :
203 : : void
204 : 24766557 : update_rhs_bc ( ncomp_t ncomp,
205 : : std::size_t nmat,
206 : : ncomp_t offset,
207 : : const std::size_t ndof,
208 : : const std::size_t ndof_l,
209 : : const tk::real wt,
210 : : const std::array< tk::real, 3 >& fn,
211 : : const std::size_t el,
212 : : const std::vector< tk::real >& fl,
213 : : const std::vector< tk::real >& B_l,
214 : : Fields& R,
215 : : std::vector< std::vector< tk::real > >& riemannDeriv )
216 : : // *****************************************************************************
217 : : // Update the rhs by adding the boundary surface integration term
218 : : //! \param[in] ncomp Number of scalar components in this PDE system
219 : : //! \param[in] nmat Number of materials in this PDE system
220 : : //! \param[in] offset Offset this PDE system operates from
221 : : //! \param[in] ndof Maximum number of degrees of freedom
222 : : //! \param[in] ndof_l Number of degrees of freedom for the left element
223 : : //! \param[in] wt Weight of gauss quadrature point
224 : : //! \param[in] fn Face/Surface normal
225 : : //! \param[in] el Left element index
226 : : //! \param[in] fl Surface flux
227 : : //! \param[in] B_l Basis function for the left element
228 : : //! \param[in,out] R Right-hand side vector computed
229 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
230 : : //! computed from the Riemann solver for use in the non-conservative terms.
231 : : //! These derivatives are used only for multi-material hydro and unused for
232 : : //! single-material compflow and linear transport.
233 : : // *****************************************************************************
234 : : {
235 : : // following line commented until rdofel is made available.
236 : : //Assert( B_l.size() == ndof_l, "Size mismatch" );
237 : :
238 [ + + ]: 98055222 : for (ncomp_t c=0; c<ncomp; ++c)
239 : : {
240 : 73288665 : auto mark = c*ndof;
241 [ + + ]: 73288665 : R(el, mark, offset) -= wt * fl[c];
242 : :
243 [ + + ]: 73288665 : if(ndof_l > 1) //DG(P1)
244 : : {
245 : 45295965 : R(el, mark+1, offset) -= wt * fl[c] * B_l[1];
246 : 45295965 : R(el, mark+2, offset) -= wt * fl[c] * B_l[2];
247 : 45295965 : R(el, mark+3, offset) -= wt * fl[c] * B_l[3];
248 : : }
249 : :
250 [ + + ]: 73288665 : if(ndof_l > 4) //DG(P2)
251 : : {
252 : 23368050 : R(el, mark+4, offset) -= wt * fl[c] * B_l[4];
253 : 23368050 : R(el, mark+5, offset) -= wt * fl[c] * B_l[5];
254 : 23368050 : R(el, mark+6, offset) -= wt * fl[c] * B_l[6];
255 : 23368050 : R(el, mark+7, offset) -= wt * fl[c] * B_l[7];
256 : 23368050 : R(el, mark+8, offset) -= wt * fl[c] * B_l[8];
257 : 23368050 : R(el, mark+9, offset) -= wt * fl[c] * B_l[9];
258 : : }
259 : : }
260 : :
261 : : // Prep for non-conservative terms in multimat
262 [ + + ]: 24766557 : if (fl.size() > ncomp)
263 : : {
264 : : // Gradients of partial pressures
265 [ + + ]: 6101400 : for (std::size_t k=0; k<nmat; ++k)
266 : : {
267 [ + + ]: 17229600 : for (std::size_t idir=0; idir<3; ++idir)
268 : 12922200 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
269 : : }
270 : :
271 : : // Divergence of velocity
272 : 1794000 : riemannDeriv[3*nmat][el] += wt * fl[ncomp+nmat];
273 : : }
274 : 24766557 : }
275 : :
276 : : } // tk::
|
