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 : : #include "Integrate/SolidTerms.hpp"
23 : : #include "Inciter/InputDeck/InputDeck.hpp"
24 : : #include "MultiMat/MiscMultiMatFns.hpp"
25 : :
26 : : namespace inciter {
27 : : extern ctr::InputDeck g_inputdeck;
28 : : }
29 : :
30 : : namespace tk {
31 : :
32 : : void
33 : 72480 : surfInt( std::size_t nmat,
34 : : const std::vector< inciter::EOS >& mat_blk,
35 : : real t,
36 : : const std::size_t ndof,
37 : : const std::size_t rdof,
38 : : const std::vector< std::size_t >& inpoel,
39 : : const std::vector< std::size_t >& /*solidx*/,
40 : : const UnsMesh::Coords& coord,
41 : : const inciter::FaceData& fd,
42 : : const Fields& geoFace,
43 : : const Fields& geoElem,
44 : : const RiemannFluxFn& flux,
45 : : const VelFn& vel,
46 : : const Fields& U,
47 : : const Fields& P,
48 : : const std::vector< std::size_t >& ndofel,
49 : : const tk::real /*dt*/,
50 : : Fields& R,
51 : : std::vector< std::vector< tk::real > >&,
52 : : std::vector< std::vector< tk::real > >&,
53 : : std::vector< std::vector< tk::real > >& riemannDeriv,
54 : : int intsharp )
55 : : // *****************************************************************************
56 : : // Compute internal surface flux integrals
57 : : //! \param[in] nmat Number of materials in this PDE system
58 : : //! \param[in] mat_blk EOS material block
59 : : //! \param[in] t Physical time
60 : : //! \param[in] ndof Maximum number of degrees of freedom
61 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
62 : : //! \param[in] inpoel Element-node connectivity
63 : : // //! \param[in] solidx Material index indicator
64 : : //! \param[in] coord Array of nodal coordinates
65 : : //! \param[in] fd Face connectivity and boundary conditions object
66 : : //! \param[in] geoFace Face geometry array
67 : : //! \param[in] geoElem Element geometry array
68 : : //! \param[in] flux Riemann flux function to use
69 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
70 : : //! \param[in] U Solution vector at recent time step
71 : : //! \param[in] P Vector of primitives at recent time step
72 : : //! \param[in] ndofel Vector of local number of degrees of freedom
73 : : // //! \param[in] dt Delta time
74 : : //! \param[in,out] R Right-hand side vector computed
75 : : //! \param[in,out] vriem Vector of the riemann velocity
76 : : //! \param[in,out] riemannLoc Vector of coordinates where Riemann velocity data
77 : : //! is available
78 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
79 : : //! computed from the Riemann solver for use in the non-conservative terms.
80 : : //! These derivatives are used only for multi-material hydro and unused for
81 : : //! single-material compflow and linear transport.
82 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
83 : : //! default 0, so that it is unused for single-material and transport.
84 : : // *****************************************************************************
85 : : {
86 : 72480 : const auto& esuf = fd.Esuf();
87 : 72480 : const auto& inpofa = fd.Inpofa();
88 : :
89 : 72480 : const auto& cx = coord[0];
90 : 72480 : const auto& cy = coord[1];
91 : 72480 : const auto& cz = coord[2];
92 : :
93 : 72480 : auto ncomp = U.nprop()/rdof;
94 : 72480 : auto nprim = P.nprop()/rdof;
95 : :
96 : : //// Determine if we have solids in our problem
97 : : //bool haveSolid = inciter::haveSolid(nmat, solidx);
98 : :
99 : : //Assert( (nmat==1 ? riemannDeriv.empty() : true), "Non-empty Riemann "
100 : : // "derivative vector for single material compflow" );
101 : :
102 : : // compute internal surface flux integrals
103 [ + + ]: 33300915 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
104 : : {
105 [ + - ][ + - ]: 33228435 : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
[ - - ][ - - ]
[ - - ]
106 : : "as -1" );
107 : :
108 : 33228435 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
109 : 33228435 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
110 : :
111 [ + - ]: 33228435 : auto ng_l = tk::NGfa(ndofel[el]);
112 [ + - ]: 33228435 : auto ng_r = tk::NGfa(ndofel[er]);
113 : :
114 : : // When the number of gauss points for the left and right element are
115 : : // different, choose the larger ng
116 : 33228435 : auto ng = std::max( ng_l, ng_r );
117 : :
118 : : // arrays for quadrature points
119 : 66456870 : std::array< std::vector< real >, 2 > coordgp;
120 : 66456870 : std::vector< real > wgp;
121 : :
122 [ + - ]: 33228435 : coordgp[0].resize( ng );
123 [ + - ]: 33228435 : coordgp[1].resize( ng );
124 [ + - ]: 33228435 : wgp.resize( ng );
125 : :
126 : : // get quadrature point weights and coordinates for triangle
127 [ + - ]: 33228435 : GaussQuadratureTri( ng, coordgp, wgp );
128 : :
129 : : // Extract the element coordinates
130 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
131 : 99685305 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
132 : 99685305 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
133 : 99685305 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
134 : 299055915 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
135 : :
136 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
137 : 99685305 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
138 : 99685305 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
139 : 99685305 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
140 : 299055915 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
141 : :
142 : : // Compute the determinant of Jacobian matrix
143 : : auto detT_l =
144 : 33228435 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
145 : : auto detT_r =
146 : 33228435 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
147 : :
148 : : // Extract the face coordinates
149 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
150 : 99685305 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
151 : 99685305 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
152 : 199370610 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
153 : :
154 : : std::array< real, 3 >
155 [ + - ][ + - ]: 33228435 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
[ + - ]
156 : :
157 : : // Gaussian quadrature
158 [ + + ]: 93695547 : for (std::size_t igp=0; igp<ng; ++igp)
159 : : {
160 : : // Compute the coordinates of quadrature point at physical domain
161 [ + - ]: 60467112 : auto gp = eval_gp( igp, coordfa, coordgp );
162 : :
163 : : // In order to determine the high-order solution from the left and right
164 : : // elements at the surface quadrature points, the basis functions from
165 : : // the left and right elements are needed. For this, a transformation to
166 : : // the reference coordinates is necessary, since the basis functions are
167 : : // defined on the reference tetrahedron only.
168 : : // The transformation relations are shown below:
169 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
170 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
171 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
172 : :
173 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
174 : : // basis functions and solutions by default. This implies that P0P1 is
175 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
176 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
177 : : // element for pDG.
178 : : std::size_t dof_el, dof_er;
179 [ + + ]: 60467112 : if (rdof > ndof)
180 : : {
181 : 4667460 : dof_el = rdof;
182 : 4667460 : dof_er = rdof;
183 : : }
184 : : else
185 : : {
186 : 55799652 : dof_el = ndofel[el];
187 : 55799652 : dof_er = ndofel[er];
188 : : }
189 : :
190 : : std::array< tk::real, 3> ref_gp_l{
191 : 60467112 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
192 : 60467112 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
193 : 120934224 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
194 : : std::array< tk::real, 3> ref_gp_r{
195 : 60467112 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
196 : 60467112 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
197 : 120934224 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
198 : :
199 : : //Compute the basis functions
200 [ + - ]: 120934224 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
201 [ + - ]: 120934224 : auto B_r = eval_basis( dof_er, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
202 : :
203 [ + - ]: 60467112 : auto wt = wgp[igp] * geoFace(f,0);
204 : :
205 : 120934224 : std::array< std::vector< real >, 2 > state;
206 : :
207 [ + - ]: 120934224 : state[0] = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim, rdof,
208 : 60467112 : nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
209 [ + - ]: 120934224 : state[1] = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim, rdof,
210 : 60467112 : nmat, er, dof_er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P);
211 : :
212 [ - + ][ - - ]: 60467112 : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
213 : : "appended boundary state vector" );
214 [ - + ][ - - ]: 60467112 : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
215 : : "appended boundary state vector" );
216 : :
217 : : // evaluate prescribed velocity (if any)
218 [ + - ]: 120934224 : auto v = vel( ncomp, gp[0], gp[1], gp[2], t );
219 : :
220 : : // compute flux
221 [ + - ]: 120934224 : auto fl = flux( mat_blk, fn, state, v );
222 : :
223 : : // Code below commented until details about the form of these terms in the
224 : : // \alpha_k g_k equations are sorted out.
225 : : // // Add RHS inverse deformation terms if necessary
226 : : // if (haveSolid)
227 : : // solidTermsSurfInt( nmat, ndof, rdof, fn, el, er, solidx, geoElem, U,
228 : : // coordel_l, coordel_r, igp, coordgp, dt, fl );
229 : :
230 : : // Add the surface integration term to the rhs
231 [ + - ]: 60467112 : update_rhs_fa( ncomp, nmat, ndof, ndofel[el], ndofel[er], wt, fn,
232 : : el, er, fl, B_l, B_r, R, riemannDeriv );
233 : : }
234 : : }
235 : 72480 : }
236 : :
237 : : void
238 : 60467112 : update_rhs_fa( ncomp_t ncomp,
239 : : std::size_t nmat,
240 : : const std::size_t ndof,
241 : : const std::size_t ndof_l,
242 : : const std::size_t ndof_r,
243 : : const tk::real wt,
244 : : const std::array< tk::real, 3 >& fn,
245 : : const std::size_t el,
246 : : const std::size_t er,
247 : : const std::vector< tk::real >& fl,
248 : : const std::vector< tk::real >& B_l,
249 : : const std::vector< tk::real >& B_r,
250 : : Fields& R,
251 : : std::vector< std::vector< tk::real > >& riemannDeriv )
252 : : // *****************************************************************************
253 : : // Update the rhs by adding the surface integration term
254 : : //! \param[in] ncomp Number of scalar components in this PDE system
255 : : //! \param[in] nmat Number of materials in this PDE system
256 : : //! \param[in] ndof Maximum number of degrees of freedom
257 : : //! \param[in] ndof_l Number of degrees of freedom for left element
258 : : //! \param[in] ndof_r Number of degrees of freedom for right element
259 : : //! \param[in] wt Weight of gauss quadrature point
260 : : //! \param[in] fn Face/Surface normal
261 : : //! \param[in] el Left element index
262 : : //! \param[in] er Right element index
263 : : //! \param[in] fl Surface flux
264 : : //! \param[in] B_l Basis function for the left element
265 : : //! \param[in] B_r Basis function for the right element
266 : : //! \param[in,out] R Right-hand side vector computed
267 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
268 : : //! computed from the Riemann solver for use in the non-conservative terms.
269 : : //! These derivatives are used only for multi-material hydro and unused for
270 : : //! single-material compflow and linear transport.
271 : : // *****************************************************************************
272 : : {
273 : : // following lines commented until rdofel is made available.
274 : : //Assert( B_l.size() == ndof_l, "Size mismatch" );
275 : : //Assert( B_r.size() == ndof_r, "Size mismatch" );
276 : :
277 : : using inciter::newSolidsAccFn;
278 : :
279 : : const auto& solidx =
280 : 60467112 : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
281 : :
282 [ + + ]: 256342182 : for (ncomp_t c=0; c<ncomp; ++c)
283 : : {
284 : 195875070 : auto mark = c*ndof;
285 : 195875070 : R(el, mark) -= wt * fl[c];
286 : 195875070 : R(er, mark) += wt * fl[c];
287 : :
288 [ + + ]: 195875070 : if(ndof_l > 1) //DG(P1)
289 : : {
290 : 116701380 : R(el, mark+1) -= wt * fl[c] * B_l[1];
291 : 116701380 : R(el, mark+2) -= wt * fl[c] * B_l[2];
292 : 116701380 : R(el, mark+3) -= wt * fl[c] * B_l[3];
293 : : }
294 : :
295 [ + + ]: 195875070 : if(ndof_r > 1) //DG(P1)
296 : : {
297 : 116686935 : R(er, mark+1) += wt * fl[c] * B_r[1];
298 : 116686935 : R(er, mark+2) += wt * fl[c] * B_r[2];
299 : 116686935 : R(er, mark+3) += wt * fl[c] * B_r[3];
300 : : }
301 : :
302 [ + + ]: 195875070 : if(ndof_l > 4) //DG(P2)
303 : : {
304 : 47897730 : R(el, mark+4) -= wt * fl[c] * B_l[4];
305 : 47897730 : R(el, mark+5) -= wt * fl[c] * B_l[5];
306 : 47897730 : R(el, mark+6) -= wt * fl[c] * B_l[6];
307 : 47897730 : R(el, mark+7) -= wt * fl[c] * B_l[7];
308 : 47897730 : R(el, mark+8) -= wt * fl[c] * B_l[8];
309 : 47897730 : R(el, mark+9) -= wt * fl[c] * B_l[9];
310 : : }
311 : :
312 [ + + ]: 195875070 : if(ndof_r > 4) //DG(P2)
313 : : {
314 : 47859300 : R(er, mark+4) += wt * fl[c] * B_r[4];
315 : 47859300 : R(er, mark+5) += wt * fl[c] * B_r[5];
316 : 47859300 : R(er, mark+6) += wt * fl[c] * B_r[6];
317 : 47859300 : R(er, mark+7) += wt * fl[c] * B_r[7];
318 : 47859300 : R(er, mark+8) += wt * fl[c] * B_r[8];
319 : 47859300 : R(er, mark+9) += wt * fl[c] * B_r[9];
320 : : }
321 : : }
322 : :
323 : : // Prep for non-conservative terms in multimat
324 [ + + ]: 60467112 : if (fl.size() > ncomp)
325 : : {
326 : : // Gradients of partial pressures
327 [ + + ]: 19659435 : for (std::size_t k=0; k<nmat; ++k)
328 : : {
329 [ + + ]: 54881760 : for (std::size_t idir=0; idir<3; ++idir)
330 : : {
331 : 41161320 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
332 : 41161320 : riemannDeriv[3*k+idir][er] -= wt * fl[ncomp+k] * fn[idir];
333 : : }
334 : : }
335 : :
336 : : // Divergence of velocity multiples basis fucntion( d(uB) / dx )
337 [ + + ]: 19898070 : for(std::size_t idof = 0; idof < ndof; idof++) {
338 : 13959075 : riemannDeriv[3*nmat+idof][el] += wt * fl[ncomp+nmat] * B_l[idof];
339 : 13959075 : riemannDeriv[3*nmat+idof][er] -= wt * fl[ncomp+nmat] * B_r[idof];
340 : : }
341 : :
342 : : // Divergence of asigma: d(asig_ij)/dx_j
343 [ + + ]: 19659435 : for (std::size_t k=0; k<nmat; ++k)
344 [ + + ]: 13720440 : if (solidx[k] > 0)
345 : : {
346 : 175320 : std::size_t mark = ncomp+nmat+1+3*(solidx[k]-1);
347 : :
348 [ + + ]: 701280 : for (std::size_t i=0; i<3; ++i)
349 : : {
350 : 525960 : riemannDeriv[3*nmat+ndof+3*(solidx[k]-1)+i][el] -=
351 : 525960 : wt * fl[mark+i];
352 : 525960 : riemannDeriv[3*nmat+ndof+3*(solidx[k]-1)+i][er] +=
353 : 525960 : wt * fl[mark+i];
354 : : }
355 : : }
356 : : }
357 : 60467112 : }
358 : :
359 : : void
360 : 1618 : surfIntFV(
361 : : std::size_t nmat,
362 : : const std::vector< inciter::EOS >& mat_blk,
363 : : real t,
364 : : const std::size_t rdof,
365 : : const std::vector< std::size_t >& inpoel,
366 : : const UnsMesh::Coords& coord,
367 : : const inciter::FaceData& fd,
368 : : const Fields& geoFace,
369 : : const Fields& geoElem,
370 : : const RiemannFluxFn& flux,
371 : : const VelFn& vel,
372 : : const Fields& U,
373 : : const Fields& P,
374 : : const std::vector< int >& srcFlag,
375 : : Fields& R,
376 : : int intsharp )
377 : : // *****************************************************************************
378 : : // Compute internal surface flux integrals for second order FV
379 : : //! \param[in] nmat Number of materials in this PDE system
380 : : //! \param[in] t Physical time
381 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
382 : : //! \param[in] inpoel Element-node connectivity
383 : : //! \param[in] coord Array of nodal coordinates
384 : : //! \param[in] fd Face connectivity and boundary conditions object
385 : : //! \param[in] geoFace Face geometry array
386 : : //! \param[in] geoElem Element geometry array
387 : : //! \param[in] flux Riemann flux function to use
388 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
389 : : //! \param[in] U Solution vector at recent time step
390 : : //! \param[in] P Vector of primitives at recent time step
391 : : //! \param[in] srcFlag Whether the energy source was added
392 : : //! \param[in,out] R Right-hand side vector computed
393 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
394 : : //! default 0, so that it is unused for single-material and transport.
395 : : // *****************************************************************************
396 : : {
397 : 1618 : const auto& esuf = fd.Esuf();
398 : :
399 : 1618 : const auto& cx = coord[0];
400 : 1618 : const auto& cy = coord[1];
401 : 1618 : const auto& cz = coord[2];
402 : :
403 : 1618 : auto ncomp = U.nprop()/rdof;
404 : 1618 : auto nprim = P.nprop()/rdof;
405 : :
406 : : // compute internal surface flux integrals
407 [ + + ]: 425758 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
408 : : {
409 [ + - ][ + - ]: 424140 : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
[ - - ][ - - ]
[ - - ]
410 : : "as -1" );
411 : :
412 : 424140 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
413 : 424140 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
414 : :
415 : : // Extract the element coordinates
416 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
417 : 1272420 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
418 : 1272420 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
419 : 1272420 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
420 : 3817260 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
421 : :
422 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
423 : 1272420 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
424 : 1272420 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
425 : 1272420 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
426 : 3817260 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
427 : :
428 : : // Compute the determinant of Jacobian matrix
429 : : auto detT_l =
430 : 424140 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
431 : : auto detT_r =
432 : 424140 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
433 : :
434 : : // face normal
435 [ + - ][ + - ]: 424140 : std::array< real, 3 > fn{{geoFace(f,1), geoFace(f,2), geoFace(f,3)}};
[ + - ]
436 : :
437 : : // face centroid
438 [ + - ][ + - ]: 424140 : std::array< real, 3 > gp{{geoFace(f,4), geoFace(f,5), geoFace(f,6)}};
[ + - ]
439 : :
440 : : // In order to determine the high-order solution from the left and right
441 : : // elements at the surface quadrature points, the basis functions from
442 : : // the left and right elements are needed. For this, a transformation to
443 : : // the reference coordinates is necessary, since the basis functions are
444 : : // defined on the reference tetrahedron only.
445 : : // The transformation relations are shown below:
446 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
447 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
448 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
449 : :
450 : : std::array< tk::real, 3> ref_gp_l{
451 : 424140 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
452 : 424140 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
453 : 848280 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
454 : : std::array< tk::real, 3> ref_gp_r{
455 : 424140 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
456 : 424140 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
457 : 848280 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
458 : :
459 : : //Compute the basis functions
460 [ + - ]: 848280 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
461 [ + - ]: 848280 : auto B_r = eval_basis( rdof, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
462 : :
463 : 848280 : std::array< std::vector< real >, 2 > state;
464 : :
465 [ + - ]: 848280 : state[0] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
466 : 848280 : nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P, srcFlag[el]);
467 [ + - ]: 848280 : state[1] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
468 : 848280 : nmat, er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P, srcFlag[er]);
469 : :
470 : : //safeReco(rdof, nmat, el, er, U, state);
471 : :
472 [ - + ][ - - ]: 424140 : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
473 : : "appended boundary state vector" );
474 [ - + ][ - - ]: 424140 : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
[ - - ][ - - ]
475 : : "appended boundary state vector" );
476 : :
477 : : // evaluate prescribed velocity (if any)
478 [ + - ]: 848280 : auto v = vel( ncomp, gp[0], gp[1], gp[2], t );
479 : :
480 : : // compute flux
481 [ + - ]: 848280 : auto fl = flux( mat_blk, fn, state, v );
482 : :
483 : : // compute non-conservative terms
484 [ + - ]: 848280 : std::vector< tk::real > var_riemann(nmat+1, 0.0);
485 [ + + ]: 1800400 : for (std::size_t k=0; k<nmat+1; ++k) var_riemann[k] = fl[ncomp+k];
486 : :
487 [ + - ]: 848280 : auto ncf_l = nonConservativeIntFV(nmat, rdof, el, fn, U, P, var_riemann);
488 [ + - ]: 848280 : auto ncf_r = nonConservativeIntFV(nmat, rdof, er, fn, U, P, var_riemann);
489 : :
490 : : // Add the surface integration term to the rhs
491 [ + + ]: 4552920 : for (ncomp_t c=0; c<ncomp; ++c)
492 : : {
493 [ + - ][ + - ]: 4128780 : R(el, c) -= geoFace(f,0) * (fl[c] - ncf_l[c]);
494 [ + - ][ + - ]: 4128780 : R(er, c) += geoFace(f,0) * (fl[c] - ncf_r[c]);
495 : : }
496 : : }
497 : 1618 : }
498 : :
499 : : } // tk::
|
