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 : :
23 : : namespace tk {
24 : :
25 : : void
26 : 74580 : surfInt( ncomp_t system,
27 : : std::size_t nmat,
28 : : ncomp_t offset,
29 : : real t,
30 : : const std::size_t ndof,
31 : : const std::size_t rdof,
32 : : const std::vector< std::size_t >& inpoel,
33 : : const UnsMesh::Coords& coord,
34 : : const inciter::FaceData& fd,
35 : : const Fields& geoFace,
36 : : const Fields& geoElem,
37 : : const RiemannFluxFn& flux,
38 : : const VelFn& vel,
39 : : const Fields& U,
40 : : const Fields& P,
41 : : const std::vector< std::size_t >& ndofel,
42 : : Fields& R,
43 : : std::vector< std::vector< tk::real > >& vriem,
44 : : std::vector< std::vector< tk::real > >& riemannLoc,
45 : : std::vector< std::vector< tk::real > >& riemannDeriv,
46 : : int intsharp )
47 : : // *****************************************************************************
48 : : // Compute internal surface flux integrals
49 : : //! \param[in] system Equation system index
50 : : //! \param[in] nmat Number of materials in this PDE system
51 : : //! \param[in] offset Offset this PDE system operates from
52 : : //! \param[in] t Physical time
53 : : //! \param[in] ndof Maximum number of degrees of freedom
54 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
55 : : //! \param[in] inpoel Element-node connectivity
56 : : //! \param[in] coord Array of nodal coordinates
57 : : //! \param[in] fd Face connectivity and boundary conditions object
58 : : //! \param[in] geoFace Face geometry array
59 : : //! \param[in] geoElem Element geometry array
60 : : //! \param[in] flux Riemann flux function to use
61 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
62 : : //! \param[in] U Solution vector at recent time step
63 : : //! \param[in] P Vector of primitives at recent time step
64 : : //! \param[in] ndofel Vector of local number of degrees of freedome
65 : : //! \param[in,out] R Right-hand side vector computed
66 : : //! \param[in,out] vriem Vector of the riemann velocity
67 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
68 : : //! is available
69 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
70 : : //! computed from the Riemann solver for use in the non-conservative terms.
71 : : //! These derivatives are used only for multi-material hydro and unused for
72 : : //! single-material compflow and linear transport.
73 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
74 : : //! default 0, so that it is unused for single-material and transport.
75 : : // *****************************************************************************
76 : : {
77 : 74580 : const auto& esuf = fd.Esuf();
78 : 74580 : const auto& inpofa = fd.Inpofa();
79 : :
80 : 74580 : const auto& cx = coord[0];
81 : 74580 : const auto& cy = coord[1];
82 : 74580 : const auto& cz = coord[2];
83 : :
84 : 74580 : auto ncomp = U.nprop()/rdof;
85 : 74580 : auto nprim = P.nprop()/rdof;
86 : :
87 [ + + ][ - + ]: 74580 : Assert( (nmat==1 ? riemannDeriv.empty() : true), "Non-empty Riemann "
[ - - ][ - - ]
[ - - ]
88 : : "derivative vector for single material compflow" );
89 : :
90 : : // compute internal surface flux integrals
91 [ + + ]: 66872745 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
92 : : {
93 [ + - ][ + - ]: 66798165 : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
[ - - ][ - - ]
[ - - ]
94 : : "as -1" );
95 : :
96 : 66798165 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
97 : 66798165 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
98 : :
99 [ + - ]: 66798165 : auto ng_l = tk::NGfa(ndofel[el]);
100 [ + - ]: 66798165 : auto ng_r = tk::NGfa(ndofel[er]);
101 : :
102 : : // When the number of gauss points for the left and right element are
103 : : // different, choose the larger ng
104 : 66798165 : auto ng = std::max( ng_l, ng_r );
105 : :
106 : : // arrays for quadrature points
107 : 133596330 : std::array< std::vector< real >, 2 > coordgp;
108 : 133596330 : std::vector< real > wgp;
109 : :
110 [ + - ]: 66798165 : coordgp[0].resize( ng );
111 [ + - ]: 66798165 : coordgp[1].resize( ng );
112 [ + - ]: 66798165 : wgp.resize( ng );
113 : :
114 : : // get quadrature point weights and coordinates for triangle
115 [ + - ]: 66798165 : GaussQuadratureTri( ng, coordgp, wgp );
116 : :
117 : : // Extract the element coordinates
118 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
119 : 200394495 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
120 : 200394495 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
121 : 200394495 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
122 : 601183485 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
123 : :
124 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
125 : 200394495 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
126 : 200394495 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
127 : 200394495 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
128 : 601183485 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
129 : :
130 : : // Compute the determinant of Jacobian matrix
131 : : auto detT_l =
132 : 66798165 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
133 : : auto detT_r =
134 : 66798165 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
135 : :
136 : : // Extract the face coordinates
137 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
138 : 200394495 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
139 : 200394495 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
140 : 400788990 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
141 : :
142 : : std::array< real, 3 >
143 [ + - ][ + - ]: 66798165 : fn{{ geoFace(f,1,0), geoFace(f,2,0), geoFace(f,3,0) }};
[ + - ]
144 : :
145 : : // Gaussian quadrature
146 [ + + ]: 163620660 : for (std::size_t igp=0; igp<ng; ++igp)
147 : : {
148 : : // Compute the coordinates of quadrature point at physical domain
149 [ + - ]: 96822495 : auto gp = eval_gp( igp, coordfa, coordgp );
150 : :
151 : : // In order to determine the high-order solution from the left and right
152 : : // elements at the surface quadrature points, the basis functions from
153 : : // the left and right elements are needed. For this, a transformation to
154 : : // the reference coordinates is necessary, since the basis functions are
155 : : // defined on the reference tetrahedron only.
156 : : // The transformation relations are shown below:
157 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
158 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
159 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
160 : :
161 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
162 : : // basis functions and solutions by default. This implies that P0P1 is
163 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
164 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
165 : : // element for pDG.
166 : : std::size_t dof_el, dof_er;
167 [ + + ]: 96822495 : if (rdof > ndof)
168 : : {
169 : 4538100 : dof_el = rdof;
170 : 4538100 : dof_er = rdof;
171 : : }
172 : : else
173 : : {
174 : 92284395 : dof_el = ndofel[el];
175 : 92284395 : dof_er = ndofel[er];
176 : : }
177 : :
178 : : std::array< tk::real, 3> ref_gp_l{
179 : 96822495 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
180 : 96822495 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
181 : 193644990 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
182 : : std::array< tk::real, 3> ref_gp_r{
183 : 96822495 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
184 : 96822495 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
185 : 193644990 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
186 : :
187 : : //Compute the basis functions
188 [ + - ]: 193644990 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
189 [ + - ]: 193644990 : auto B_r = eval_basis( dof_er, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
190 : :
191 [ + - ]: 96822495 : auto wt = wgp[igp] * geoFace(f,0,0);
192 : :
193 : 193644990 : std::array< std::vector< real >, 2 > state;
194 : :
195 [ + - ]: 193644990 : state[0] = evalPolynomialSol(system, offset, intsharp, ncomp, nprim, rdof,
196 : 96822495 : nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
197 [ + - ]: 193644990 : state[1] = evalPolynomialSol(system, offset, intsharp, ncomp, nprim, rdof,
198 : 96822495 : nmat, er, dof_er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P);
199 : :
200 [ - + ][ - - ]: 96822495 : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
201 : : "appended boundary state vector" );
202 [ - + ][ - - ]: 96822495 : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
203 : : "appended boundary state vector" );
204 : :
205 : : // evaluate prescribed velocity (if any)
206 [ + - ]: 193644990 : auto v = vel( system, ncomp, gp[0], gp[1], gp[2], t );
207 : :
208 : : // compute flux
209 [ + - ]: 193644990 : auto fl = flux( fn, state, v );
210 : :
211 : : // Add the surface integration term to the rhs
212 [ + - ]: 96822495 : update_rhs_fa( ncomp, nmat, offset, ndof, ndofel[el], ndofel[er], wt, fn,
213 : : el, er, fl, B_l, B_r, R, riemannDeriv );
214 : :
215 : : // Store the riemann velocity and coordinates data of quadrature point
216 : : // used for velocity reconstruction if MulMat scheme is selected
217 [ + + ][ + + ]: 96822495 : if (nmat > 1 && ndof > 1)
218 [ + - ]: 2436300 : tk::evaluRiemann( ncomp, esuf[2*f], esuf[2*f+1], nmat, fl, fn, gp,
219 : : state, vriem, riemannLoc );
220 : :
221 : : }
222 : : }
223 : 74580 : }
224 : :
225 : : void
226 : 96822495 : update_rhs_fa( ncomp_t ncomp,
227 : : std::size_t nmat,
228 : : ncomp_t offset,
229 : : const std::size_t ndof,
230 : : const std::size_t ndof_l,
231 : : const std::size_t ndof_r,
232 : : const tk::real wt,
233 : : const std::array< tk::real, 3 >& fn,
234 : : const std::size_t el,
235 : : const std::size_t er,
236 : : const std::vector< tk::real >& fl,
237 : : const std::vector< tk::real >& B_l,
238 : : const std::vector< tk::real >& B_r,
239 : : Fields& R,
240 : : std::vector< std::vector< tk::real > >& riemannDeriv )
241 : : // *****************************************************************************
242 : : // Update the rhs by adding the surface integration term
243 : : //! \param[in] ncomp Number of scalar components in this PDE system
244 : : //! \param[in] nmat Number of materials in this PDE system
245 : : //! \param[in] offset Offset this PDE system operates from
246 : : //! \param[in] ndof Maximum number of degrees of freedom
247 : : //! \param[in] ndof_l Number of degrees of freedom for left element
248 : : //! \param[in] ndof_r Number of degrees of freedom for right element
249 : : //! \param[in] wt Weight of gauss quadrature point
250 : : //! \param[in] fn Face/Surface normal
251 : : //! \param[in] el Left element index
252 : : //! \param[in] er Right element index
253 : : //! \param[in] fl Surface flux
254 : : //! \param[in] B_l Basis function for the left element
255 : : //! \param[in] B_r Basis function for the right element
256 : : //! \param[in,out] R Right-hand side vector computed
257 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
258 : : //! computed from the Riemann solver for use in the non-conservative terms.
259 : : //! These derivatives are used only for multi-material hydro and unused for
260 : : //! single-material compflow and linear transport.
261 : : // *****************************************************************************
262 : : {
263 : : // following lines commented until rdofel is made available.
264 : : //Assert( B_l.size() == ndof_l, "Size mismatch" );
265 : : //Assert( B_r.size() == ndof_r, "Size mismatch" );
266 : :
267 [ + + ]: 339266400 : for (ncomp_t c=0; c<ncomp; ++c)
268 : : {
269 : 242443905 : auto mark = c*ndof;
270 : 242443905 : R(el, mark, offset) -= wt * fl[c];
271 : 242443905 : R(er, mark, offset) += wt * fl[c];
272 : :
273 [ + + ]: 242443905 : if(ndof_l > 1) //DG(P1)
274 : : {
275 : 132155640 : R(el, mark+1, offset) -= wt * fl[c] * B_l[1];
276 : 132155640 : R(el, mark+2, offset) -= wt * fl[c] * B_l[2];
277 : 132155640 : R(el, mark+3, offset) -= wt * fl[c] * B_l[3];
278 : : }
279 : :
280 [ + + ]: 242443905 : if(ndof_r > 1) //DG(P1)
281 : : {
282 : 132159645 : R(er, mark+1, offset) += wt * fl[c] * B_r[1];
283 : 132159645 : R(er, mark+2, offset) += wt * fl[c] * B_r[2];
284 : 132159645 : R(er, mark+3, offset) += wt * fl[c] * B_r[3];
285 : : }
286 : :
287 [ + + ]: 242443905 : if(ndof_l > 4) //DG(P2)
288 : : {
289 : 65794140 : R(el, mark+4, offset) -= wt * fl[c] * B_l[4];
290 : 65794140 : R(el, mark+5, offset) -= wt * fl[c] * B_l[5];
291 : 65794140 : R(el, mark+6, offset) -= wt * fl[c] * B_l[6];
292 : 65794140 : R(el, mark+7, offset) -= wt * fl[c] * B_l[7];
293 : 65794140 : R(el, mark+8, offset) -= wt * fl[c] * B_l[8];
294 : 65794140 : R(el, mark+9, offset) -= wt * fl[c] * B_l[9];
295 : : }
296 : :
297 [ + + ]: 242443905 : if(ndof_r > 4) //DG(P2)
298 : : {
299 : 65737440 : R(er, mark+4, offset) += wt * fl[c] * B_r[4];
300 : 65737440 : R(er, mark+5, offset) += wt * fl[c] * B_r[5];
301 : 65737440 : R(er, mark+6, offset) += wt * fl[c] * B_r[6];
302 : 65737440 : R(er, mark+7, offset) += wt * fl[c] * B_r[7];
303 : 65737440 : R(er, mark+8, offset) += wt * fl[c] * B_r[8];
304 : 65737440 : R(er, mark+9, offset) += wt * fl[c] * B_r[9];
305 : : }
306 : : }
307 : :
308 : : // Prep for non-conservative terms in multimat
309 [ + + ]: 96822495 : if (fl.size() > ncomp)
310 : : {
311 : : // Gradients of partial pressures
312 [ + + ]: 18560175 : for (std::size_t k=0; k<nmat; ++k)
313 : : {
314 [ + + ]: 51950400 : for (std::size_t idir=0; idir<3; ++idir)
315 : : {
316 : 38962800 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
317 : 38962800 : riemannDeriv[3*k+idir][er] -= wt * fl[ncomp+k] * fn[idir];
318 : : }
319 : : }
320 : :
321 : : // Divergence of velocity
322 : 5572575 : riemannDeriv[3*nmat][el] += wt * fl[ncomp+nmat];
323 : 5572575 : riemannDeriv[3*nmat][er] -= wt * fl[ncomp+nmat];
324 : : }
325 : 96822495 : }
326 : :
327 : : } // tk::
|
