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 : : #include "EoS/GetMatProp.hpp"
26 : :
27 : : namespace inciter {
28 : : extern ctr::InputDeck g_inputdeck;
29 : : }
30 : :
31 : : namespace tk {
32 : :
33 : : void
34 : 38265 : surfInt( const bool pref,
35 : : std::size_t nmat,
36 : : const std::vector< inciter::EOS >& mat_blk,
37 : : real t,
38 : : const std::size_t ndof,
39 : : const std::size_t rdof,
40 : : const std::vector< std::size_t >& inpoel,
41 : : const std::vector< std::size_t >& /*solidx*/,
42 : : const UnsMesh::Coords& coord,
43 : : const inciter::FaceData& fd,
44 : : const Fields& geoFace,
45 : : const Fields& geoElem,
46 : : const RiemannFluxFn& flux,
47 : : const VelFn& vel,
48 : : const Fields& U,
49 : : const Fields& P,
50 : : const std::vector< std::size_t >& ndofel,
51 : : const tk::real /*dt*/,
52 : : Fields& R,
53 : : std::vector< std::vector< tk::real > >& riemannDeriv,
54 : : int intsharp )
55 : : // *****************************************************************************
56 : : // Compute internal surface flux integrals
57 : : //! \param[in] pref Indicator for p-adaptive algorithm
58 : : //! \param[in] nmat Number of materials in this PDE system
59 : : //! \param[in] mat_blk EOS material block
60 : : //! \param[in] t Physical time
61 : : //! \param[in] ndof Maximum number of degrees of freedom
62 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
63 : : //! \param[in] inpoel Element-node connectivity
64 : : // //! \param[in] solidx Material index indicator
65 : : //! \param[in] coord Array of nodal coordinates
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] flux Riemann flux function to use
70 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
71 : : //! \param[in] U Solution vector at recent time step
72 : : //! \param[in] P Vector of primitives at recent time step
73 : : //! \param[in] ndofel Vector of local number of degrees of freedom
74 : : // //! \param[in] dt Delta time
75 : : //! \param[in,out] R Right-hand side vector computed
76 : : //! \param[in,out] riemannDeriv Derivatives of partial-pressures and velocities
77 : : //! computed from the Riemann solver for use in the non-conservative terms.
78 : : //! These derivatives are used only for multi-material hydro and unused for
79 : : //! single-material compflow and linear transport.
80 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
81 : : //! default 0, so that it is unused for single-material and transport.
82 : : // *****************************************************************************
83 : : {
84 : : const auto& esuf = fd.Esuf();
85 : : const auto& inpofa = fd.Inpofa();
86 : :
87 : : const auto& cx = coord[0];
88 : : const auto& cy = coord[1];
89 : : const auto& cz = coord[2];
90 : :
91 : 38265 : auto ncomp = U.nprop()/rdof;
92 : 38265 : auto nprim = P.nprop()/rdof;
93 : :
94 : : //// Determine if we have solids in our problem
95 : : //bool haveSolid = inciter::haveSolid(nmat, solidx);
96 : :
97 : : //Assert( (nmat==1 ? riemannDeriv.empty() : true), "Non-empty Riemann "
98 : : // "derivative vector for single material compflow" );
99 : :
100 : : // compute internal surface flux integrals
101 [ + + ]: 17971050 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
102 : : {
103 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
104 : : "as -1" );
105 : :
106 : 17932785 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
107 : 17932785 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
108 : :
109 : 17932785 : auto ng_l = tk::NGfa(ndofel[el]);
110 [ + + ]: 17932785 : auto ng_r = tk::NGfa(ndofel[er]);
111 : :
112 : : // When the number of gauss points for the left and right element are
113 : : // different, choose the larger ng
114 : 17932785 : auto ng = std::max( ng_l, ng_r );
115 : :
116 : : // arrays for quadrature points
117 : : std::array< std::vector< real >, 2 > coordgp;
118 : : std::vector< real > wgp;
119 : :
120 [ + - ]: 17932785 : coordgp[0].resize( ng );
121 [ + - ]: 17932785 : coordgp[1].resize( ng );
122 [ + - ]: 17932785 : wgp.resize( ng );
123 : :
124 : : // get quadrature point weights and coordinates for triangle
125 [ + - ]: 17932785 : GaussQuadratureTri( ng, coordgp, wgp );
126 : :
127 : : // Extract the element coordinates
128 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
129 : 17932785 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
130 : 17932785 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
131 : 17932785 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
132 : 17932785 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
133 : :
134 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
135 : 17932785 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
136 : 17932785 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
137 : 17932785 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
138 : 17932785 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
139 : :
140 : : // Compute the determinant of Jacobian matrix
141 : : auto detT_l =
142 : 17932785 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
143 : : auto detT_r =
144 : 17932785 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
145 : :
146 : : // Extract the face coordinates
147 : : std::array< std::array< tk::real, 3>, 3 > coordfa {{
148 : 17932785 : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
149 : 17932785 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
150 : 17932785 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }} }};
151 : :
152 : : std::array< real, 3 >
153 : 17932785 : fn{{ geoFace(f,1), geoFace(f,2), geoFace(f,3) }};
154 : :
155 : : // Gaussian quadrature
156 [ + + ]: 48748122 : for (std::size_t igp=0; igp<ng; ++igp)
157 : : {
158 : : // Compute the coordinates of quadrature point at physical domain
159 [ + - ]: 30815337 : auto gp = eval_gp( igp, coordfa, coordgp );
160 : :
161 : : // In order to determine the high-order solution from the left and right
162 : : // elements at the surface quadrature points, the basis functions from
163 : : // the left and right elements are needed. For this, a transformation to
164 : : // the reference coordinates is necessary, since the basis functions are
165 : : // defined on the reference tetrahedron only.
166 : : // The transformation relations are shown below:
167 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
168 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
169 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
170 : :
171 : : // If an rDG method is set up (P0P1), then, currently we compute the P1
172 : : // basis functions and solutions by default. This implies that P0P1 is
173 : : // unsupported in the p-adaptive DG (PDG). This is a workaround until we
174 : : // have rdofel, which is needed to distinguish between ndofs and rdofs per
175 : : // element for pDG.
176 : : std::size_t dof_el, dof_er;
177 [ + + ]: 30815337 : if (rdof > ndof)
178 : : {
179 : : dof_el = rdof;
180 : : dof_er = rdof;
181 : : }
182 : : else
183 : : {
184 : 28724187 : dof_el = ndofel[el];
185 : 28724187 : dof_er = ndofel[er];
186 : : }
187 : :
188 : : // For multi-material p-adaptive simulation, when dofel = 1, p0p1 is
189 : : // applied and ndof for solution evaluation should be 4
190 [ - + ]: 30815337 : if(ncomp > 5 && pref) {
191 [ - - ]: 0 : if(dof_el == 1)
192 : : dof_el = 4;
193 [ - - ]: 0 : if(dof_er == 1)
194 : : dof_er = 4;
195 : : }
196 : :
197 : : std::array< tk::real, 3> ref_gp_l{
198 : 30815337 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
199 : 30815337 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
200 : 30815337 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
201 : : std::array< tk::real, 3> ref_gp_r{
202 : 30815337 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
203 : 30815337 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
204 : 30815337 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
205 : :
206 : : //Compute the basis functions
207 [ + - ]: 30815337 : auto B_l = eval_basis( dof_el, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
208 [ + - ]: 30815337 : auto B_r = eval_basis( dof_er, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
209 : :
210 : 30815337 : auto wt = wgp[igp] * geoFace(f,0);
211 : :
212 : : std::array< std::vector< real >, 2 > state;
213 : :
214 [ + - ]: 30815337 : state[0] = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim, rdof,
215 : : nmat, el, dof_el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P);
216 [ + - ]: 30815337 : state[1] = evalPolynomialSol(mat_blk, intsharp, ncomp, nprim, rdof,
217 : : nmat, er, dof_er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P);
218 : :
219 : : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
220 : : "appended boundary state vector" );
221 : : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
222 : : "appended boundary state vector" );
223 : :
224 : : // evaluate prescribed velocity (if any)
225 [ - + ][ - + ]: 61630674 : auto v = vel( ncomp, gp[0], gp[1], gp[2], t );
226 : :
227 : : // compute flux
228 : : auto fl = flux( mat_blk, fn, state, v );
229 : :
230 : : // Add the surface integration term to the rhs
231 [ + - ]: 30815337 : 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 : 38265 : }
236 : :
237 : : void
238 : 30815337 : 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 : : inciter::g_inputdeck.get< tag::matidxmap, tag::solidx >();
281 : :
282 [ + + ]: 204461322 : for (ncomp_t c=0; c<ncomp; ++c)
283 : : {
284 : 173645985 : auto mark = c*ndof;
285 [ + + ]: 173645985 : R(el, mark) -= wt * fl[c];
286 : 173645985 : R(er, mark) += wt * fl[c];
287 : :
288 [ + + ]: 173645985 : if(ndof_l > 1) //DG(P1)
289 : : {
290 : 96546195 : R(el, mark+1) -= wt * fl[c] * B_l[1];
291 : 96546195 : R(el, mark+2) -= wt * fl[c] * B_l[2];
292 : 96546195 : R(el, mark+3) -= wt * fl[c] * B_l[3];
293 : : }
294 : :
295 [ + + ]: 173645985 : if(ndof_r > 1) //DG(P1)
296 : : {
297 : 96537420 : R(er, mark+1) += wt * fl[c] * B_r[1];
298 : 96537420 : R(er, mark+2) += wt * fl[c] * B_r[2];
299 : 96537420 : R(er, mark+3) += wt * fl[c] * B_r[3];
300 : : }
301 : :
302 [ + + ]: 173645985 : if(ndof_l > 4) //DG(P2)
303 : : {
304 : 42356880 : R(el, mark+4) -= wt * fl[c] * B_l[4];
305 : 42356880 : R(el, mark+5) -= wt * fl[c] * B_l[5];
306 : 42356880 : R(el, mark+6) -= wt * fl[c] * B_l[6];
307 : 42356880 : R(el, mark+7) -= wt * fl[c] * B_l[7];
308 : 42356880 : R(el, mark+8) -= wt * fl[c] * B_l[8];
309 : 42356880 : R(el, mark+9) -= wt * fl[c] * B_l[9];
310 : : }
311 : :
312 [ + + ]: 173645985 : if(ndof_r > 4) //DG(P2)
313 : : {
314 : 42345450 : R(er, mark+4) += wt * fl[c] * B_r[4];
315 : 42345450 : R(er, mark+5) += wt * fl[c] * B_r[5];
316 : 42345450 : R(er, mark+6) += wt * fl[c] * B_r[6];
317 : 42345450 : R(er, mark+7) += wt * fl[c] * B_r[7];
318 : 42345450 : R(er, mark+8) += wt * fl[c] * B_r[8];
319 : 42345450 : R(er, mark+9) += wt * fl[c] * B_r[9];
320 : : }
321 : : }
322 : :
323 : : // Prep for non-conservative terms in multimat
324 [ + + ]: 30815337 : if (fl.size() > ncomp)
325 : : {
326 : : // Gradients of partial pressures
327 [ + + ]: 20699205 : for (std::size_t k=0; k<nmat; ++k)
328 : : {
329 [ + + ]: 57654480 : for (std::size_t idir=0; idir<3; ++idir)
330 : : {
331 : 43240860 : riemannDeriv[3*k+idir][el] += wt * fl[ncomp+k] * fn[idir];
332 : 43240860 : 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 [ + + ]: 20591250 : for(std::size_t idof = 0; idof < ndof_l; idof++) {
338 : 14305665 : riemannDeriv[3*nmat+idof][el] += wt * fl[ncomp+nmat] * B_l[idof];
339 : : }
340 [ + + ]: 20591250 : for(std::size_t idof = 0; idof < ndof_r; idof++) {
341 : 14305665 : riemannDeriv[3*nmat+idof][er] -= wt * fl[ncomp+nmat] * B_r[idof];
342 : : }
343 : :
344 : : // Divergence of asigma: d(asig_ij)/dx_j
345 [ + + ]: 20699205 : for (std::size_t k=0; k<nmat; ++k)
346 [ - + ]: 14413620 : if (solidx[k] > 0)
347 : : {
348 : 0 : std::size_t mark = ncomp+nmat+1+3*(solidx[k]-1);
349 : :
350 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
351 : : {
352 : 0 : riemannDeriv[3*nmat+ndof+3*(solidx[k]-1)+i][el] -=
353 : 0 : wt * fl[mark+i];
354 : 0 : riemannDeriv[3*nmat+ndof+3*(solidx[k]-1)+i][er] +=
355 : 0 : wt * fl[mark+i];
356 : : }
357 : : }
358 : :
359 : : // Derivatives of g: d(g_il)/d(x_j)-d(g_ij)/d(x_l)
360 : : // for i=1,2,3; j=1,2,3; l=1,2,3. Total = 3x3x3 (per solid)
361 : 6285585 : std::size_t nsld = inciter::numSolids(nmat, solidx);
362 [ + + ]: 20699205 : for (std::size_t k=0; k<nmat; ++k)
363 [ - + ]: 14413620 : if (solidx[k] > 0)
364 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
365 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
366 [ - - ]: 0 : for (std::size_t l=0; l<3; ++l)
367 [ - - ]: 0 : if (j != l)
368 : : {
369 : 0 : std::size_t mark1 = ncomp+nmat+1+3*nsld+9*(solidx[k]-1)+3*i+l;
370 : 0 : std::size_t mark2 = ncomp+nmat+1+3*nsld+9*(solidx[k]-1)+3*i+j;
371 : 0 : riemannDeriv[3*nmat+ndof+3*nsld+newSolidsAccFn(k,i,j,l)][el] -=
372 : 0 : wt * ( fl[mark1] * fn[j] - fl[mark2] * fn[l]);
373 : 0 : riemannDeriv[3*nmat+ndof+3*nsld+newSolidsAccFn(k,i,j,l)][er] +=
374 : 0 : wt * ( fl[mark1] * fn[j] - fl[mark2] * fn[l]);
375 : : }
376 : : }
377 : 30815337 : }
378 : :
379 : : void
380 : 1368 : surfIntFV(
381 : : std::size_t nmat,
382 : : const std::vector< inciter::EOS >& mat_blk,
383 : : real t,
384 : : const std::size_t rdof,
385 : : const std::vector< std::size_t >& inpoel,
386 : : const UnsMesh::Coords& coord,
387 : : const inciter::FaceData& fd,
388 : : const Fields& geoFace,
389 : : const Fields& geoElem,
390 : : const RiemannFluxFn& flux,
391 : : const VelFn& vel,
392 : : const Fields& U,
393 : : const Fields& P,
394 : : const std::vector< int >& srcFlag,
395 : : Fields& R,
396 : : int intsharp )
397 : : // *****************************************************************************
398 : : // Compute internal surface flux integrals for second order FV
399 : : //! \param[in] nmat Number of materials in this PDE system
400 : : //! \param[in] t Physical time
401 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
402 : : //! \param[in] inpoel Element-node connectivity
403 : : //! \param[in] coord Array of nodal coordinates
404 : : //! \param[in] fd Face connectivity and boundary conditions object
405 : : //! \param[in] geoFace Face geometry array
406 : : //! \param[in] geoElem Element geometry array
407 : : //! \param[in] flux Riemann flux function to use
408 : : //! \param[in] vel Function to use to query prescribed velocity (if any)
409 : : //! \param[in] U Solution vector at recent time step
410 : : //! \param[in] P Vector of primitives at recent time step
411 : : //! \param[in] srcFlag Whether the energy source was added
412 : : //! \param[in,out] R Right-hand side vector computed
413 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
414 : : //! default 0, so that it is unused for single-material and transport.
415 : : // *****************************************************************************
416 : : {
417 : : const auto& esuf = fd.Esuf();
418 : :
419 : : const auto& cx = coord[0];
420 : : const auto& cy = coord[1];
421 : : const auto& cz = coord[2];
422 : :
423 : 1368 : auto ncomp = U.nprop()/rdof;
424 : 1368 : auto nprim = P.nprop()/rdof;
425 : :
426 : : // compute internal surface flux integrals
427 [ + + ]: 476888 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
428 : : {
429 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
430 : : "as -1" );
431 : :
432 : 475520 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
433 : 475520 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
434 : :
435 : : // Extract the element coordinates
436 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
437 : 475520 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
438 : 475520 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
439 : 475520 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
440 : 475520 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
441 : :
442 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
443 : 475520 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
444 : 475520 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
445 : 475520 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
446 : 475520 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
447 : :
448 : : // Compute the determinant of Jacobian matrix
449 : : auto detT_l =
450 : 475520 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
451 : : auto detT_r =
452 : 475520 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
453 : :
454 : : // face normal
455 : 475520 : std::array< real, 3 > fn{{geoFace(f,1), geoFace(f,2), geoFace(f,3)}};
456 : :
457 : : // face centroid
458 : 475520 : std::array< real, 3 > gp{{geoFace(f,4), geoFace(f,5), geoFace(f,6)}};
459 : :
460 : : // In order to determine the high-order solution from the left and right
461 : : // elements at the surface quadrature points, the basis functions from
462 : : // the left and right elements are needed. For this, a transformation to
463 : : // the reference coordinates is necessary, since the basis functions are
464 : : // defined on the reference tetrahedron only.
465 : : // The transformation relations are shown below:
466 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
467 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
468 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
469 : :
470 : : std::array< tk::real, 3> ref_gp_l{
471 : 475520 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
472 : 475520 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
473 : 475520 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
474 : : std::array< tk::real, 3> ref_gp_r{
475 : 475520 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
476 : 475520 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
477 : 475520 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
478 : :
479 : : //Compute the basis functions
480 : 475520 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
481 [ + - ]: 475520 : auto B_r = eval_basis( rdof, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
482 : :
483 : : std::array< std::vector< real >, 2 > state;
484 : :
485 [ + - ]: 475520 : state[0] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
486 [ + - ]: 475520 : nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P, srcFlag[el]);
487 [ + - ]: 475520 : state[1] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
488 [ + - ]: 475520 : nmat, er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P, srcFlag[er]);
489 : :
490 : : //safeReco(rdof, nmat, el, er, U, state);
491 : :
492 : : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
493 : : "appended boundary state vector" );
494 : : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
495 : : "appended boundary state vector" );
496 : :
497 : : // evaluate prescribed velocity (if any)
498 [ - + ][ - + ]: 951040 : auto v = vel( ncomp, gp[0], gp[1], gp[2], t );
499 : :
500 : : // compute flux
501 : : auto fl = flux( mat_blk, fn, state, v );
502 : :
503 : : // compute non-conservative terms
504 [ + - ][ - - ]: 475520 : std::vector< tk::real > var_riemann(nmat+1, 0.0);
505 [ + + ]: 2005920 : for (std::size_t k=0; k<nmat+1; ++k) var_riemann[k] = fl[ncomp+k];
506 : :
507 [ + - ]: 475520 : auto ncf_l = nonConservativeIntFV(nmat, rdof, el, fn, U, P, var_riemann);
508 [ + - ]: 475520 : auto ncf_r = nonConservativeIntFV(nmat, rdof, er, fn, U, P, var_riemann);
509 : :
510 : : // Add the surface integration term to the rhs
511 [ + + ]: 5066720 : for (ncomp_t c=0; c<ncomp; ++c)
512 : : {
513 : 4591200 : R(el, c) -= geoFace(f,0) * (fl[c] - ncf_l[c]);
514 : 4591200 : R(er, c) += geoFace(f,0) * (fl[c] - ncf_r[c]);
515 : : }
516 : : }
517 : 1368 : }
518 : :
519 : : void
520 : 0 : surfIntViscousFV(
521 : : std::size_t nmat,
522 : : const std::vector< inciter::EOS >& mat_blk,
523 : : const std::size_t rdof,
524 : : const std::vector< std::size_t >& inpoel,
525 : : const UnsMesh::Coords& coord,
526 : : const inciter::FaceData& fd,
527 : : const Fields& geoFace,
528 : : const Fields& geoElem,
529 : : const Fields& U,
530 : : const Fields& P,
531 : : const Fields& T,
532 : : const std::vector< int >& srcFlag,
533 : : Fields& R,
534 : : int intsharp )
535 : : // *****************************************************************************
536 : : // Compute internal surface viscous flux integrals for second order FV
537 : : //! \param[in] nmat Number of materials in this PDE system
538 : : //! \param[in] mat_blk Material EOS block
539 : : //! \param[in] rdof Maximum number of reconstructed degrees of freedom
540 : : //! \param[in] inpoel Element-node connectivity
541 : : //! \param[in] coord Array of nodal coordinates
542 : : //! \param[in] fd Face connectivity and boundary conditions object
543 : : //! \param[in] geoFace Face geometry array
544 : : //! \param[in] geoElem Element geometry array
545 : : //! \param[in] U Solution vector at recent time step
546 : : //! \param[in] P Vector of primitives at recent time step
547 : : //! \param[in] T Vector of temperatures at recent time step
548 : : //! \param[in] srcFlag Whether the energy source was added
549 : : //! \param[in,out] R Right-hand side vector computed
550 : : //! \param[in] intsharp Interface compression tag, an optional argument, with
551 : : //! default 0, so that it is unused for single-material and transport.
552 : : // *****************************************************************************
553 : : {
554 : : using inciter::velocityDofIdx;
555 : :
556 : : const auto& esuf = fd.Esuf();
557 : :
558 : : const auto& cx = coord[0];
559 : : const auto& cy = coord[1];
560 : : const auto& cz = coord[2];
561 : :
562 : 0 : auto ncomp = U.nprop()/rdof;
563 : 0 : auto nprim = P.nprop()/rdof;
564 : :
565 : : // compute internal surface flux integrals
566 [ - - ]: 0 : for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
567 : : {
568 : : Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
569 : : "as -1" );
570 : :
571 : 0 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
572 : 0 : std::size_t er = static_cast< std::size_t >(esuf[2*f+1]);
573 : :
574 : : // Extract the element coordinates
575 : : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
576 : 0 : {{ cx[ inpoel[4*el ] ], cy[ inpoel[4*el ] ], cz[ inpoel[4*el ] ] }},
577 : 0 : {{ cx[ inpoel[4*el+1] ], cy[ inpoel[4*el+1] ], cz[ inpoel[4*el+1] ] }},
578 : 0 : {{ cx[ inpoel[4*el+2] ], cy[ inpoel[4*el+2] ], cz[ inpoel[4*el+2] ] }},
579 : 0 : {{ cx[ inpoel[4*el+3] ], cy[ inpoel[4*el+3] ], cz[ inpoel[4*el+3] ] }} }};
580 : :
581 : : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
582 : 0 : {{ cx[ inpoel[4*er ] ], cy[ inpoel[4*er ] ], cz[ inpoel[4*er ] ] }},
583 : 0 : {{ cx[ inpoel[4*er+1] ], cy[ inpoel[4*er+1] ], cz[ inpoel[4*er+1] ] }},
584 : 0 : {{ cx[ inpoel[4*er+2] ], cy[ inpoel[4*er+2] ], cz[ inpoel[4*er+2] ] }},
585 : 0 : {{ cx[ inpoel[4*er+3] ], cy[ inpoel[4*er+3] ], cz[ inpoel[4*er+3] ] }} }};
586 : :
587 : : // Compute the determinant of Jacobian matrix
588 : : auto detT_l =
589 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
590 : : auto detT_r =
591 : 0 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
592 : :
593 : : // face normal
594 : 0 : std::array< real, 3 > fn{{geoFace(f,1), geoFace(f,2), geoFace(f,3)}};
595 : :
596 : : // face centroid
597 : 0 : std::array< real, 3 > gp{{geoFace(f,4), geoFace(f,5), geoFace(f,6)}};
598 : :
599 : : // In order to determine the high-order solution from the left and right
600 : : // elements at the surface quadrature points, the basis functions from
601 : : // the left and right elements are needed. For this, a transformation to
602 : : // the reference coordinates is necessary, since the basis functions are
603 : : // defined on the reference tetrahedron only.
604 : : // The transformation relations are shown below:
605 : : // xi = Jacobian( coordel[0], gp, coordel[2], coordel[3] ) / detT
606 : : // eta = Jacobian( coordel[0], coordel[2], gp, coordel[3] ) / detT
607 : : // zeta = Jacobian( coordel[0], coordel[2], coordel[3], gp ) / detT
608 : :
609 : : std::array< tk::real, 3> ref_gp_l{
610 : 0 : Jacobian( coordel_l[0], gp, coordel_l[2], coordel_l[3] ) / detT_l,
611 : 0 : Jacobian( coordel_l[0], coordel_l[1], gp, coordel_l[3] ) / detT_l,
612 : 0 : Jacobian( coordel_l[0], coordel_l[1], coordel_l[2], gp ) / detT_l };
613 : : std::array< tk::real, 3> ref_gp_r{
614 : 0 : Jacobian( coordel_r[0], gp, coordel_r[2], coordel_r[3] ) / detT_r,
615 : 0 : Jacobian( coordel_r[0], coordel_r[1], gp, coordel_r[3] ) / detT_r,
616 : 0 : Jacobian( coordel_r[0], coordel_r[1], coordel_r[2], gp ) / detT_r };
617 : :
618 : : //Compute the basis functions
619 : 0 : auto B_l = eval_basis( rdof, ref_gp_l[0], ref_gp_l[1], ref_gp_l[2] );
620 [ - - ]: 0 : auto B_r = eval_basis( rdof, ref_gp_r[0], ref_gp_r[1], ref_gp_r[2] );
621 : :
622 : : std::array< std::vector< real >, 2 > state;
623 : :
624 [ - - ]: 0 : state[0] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
625 [ - - ]: 0 : nmat, el, inpoel, coord, geoElem, ref_gp_l, B_l, U, P, srcFlag[el]);
626 [ - - ]: 0 : state[1] = evalFVSol(mat_blk, intsharp, ncomp, nprim, rdof,
627 [ - - ]: 0 : nmat, er, inpoel, coord, geoElem, ref_gp_r, B_r, U, P, srcFlag[er]);
628 : :
629 : : Assert( state[0].size() == ncomp+nprim, "Incorrect size for "
630 : : "appended boundary state vector" );
631 : : Assert( state[1].size() == ncomp+nprim, "Incorrect size for "
632 : : "appended boundary state vector" );
633 : :
634 : : // cell averaged state for computing the diffusive flux
635 : : std::array< std::vector< real >, 2 > cellAvgState;
636 [ - - ]: 0 : std::vector< tk::real > Bcc(rdof, 0.0);
637 : 0 : Bcc[0] = 1.0;
638 : :
639 [ - - ]: 0 : cellAvgState[0] = evalFVSol(mat_blk, 0, ncomp, nprim, rdof,
640 : : nmat, el, inpoel, coord, geoElem, {{0.25, 0.25, 0.25}}, Bcc, U, P,
641 [ - - ]: 0 : srcFlag[el]);
642 [ - - ][ - - ]: 0 : cellAvgState[1] = evalFVSol(mat_blk, 0, ncomp, nprim, rdof,
643 : : nmat, er, inpoel, coord, geoElem, {{0.25, 0.25, 0.25}}, Bcc, U, P,
644 [ - - ]: 0 : srcFlag[er]);
645 : :
646 : : Assert( cellAvgState[0].size() == ncomp+nprim, "Incorrect size for "
647 : : "appended cell-averaged state vector" );
648 : : Assert( cellAvgState[1].size() == ncomp+nprim, "Incorrect size for "
649 : : "appended cell-averaged state vector" );
650 : :
651 : : std::array< std::vector< real >, 2 > cellAvgT{{
652 [ - - ][ - - ]: 0 : std::vector<real>(nmat), std::vector<real>(nmat) }};
653 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
654 : 0 : cellAvgT[0][k] = T(el, k*rdof);
655 : 0 : cellAvgT[1][k] = T(er, k*rdof);
656 : : }
657 : :
658 : : std::array< std::array< real, 3 >, 2 > centroids{{
659 : 0 : {{geoElem(el,1), geoElem(el,2), geoElem(el,3)}},
660 : 0 : {{geoElem(er,1), geoElem(er,2), geoElem(er,3)}} }};
661 : :
662 : : // Numerical viscous flux
663 : : // -------------------------------------------------------------------------
664 : :
665 : : // 1. Get spatial gradient from Dubiner dofs
666 : : auto jacInv_l =
667 : 0 : tk::inverseJacobian( coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3] );
668 [ - - ]: 0 : auto dBdx_l = tk::eval_dBdx_p1( rdof, jacInv_l );
669 : : auto jacInv_r =
670 : 0 : tk::inverseJacobian( coordel_r[0], coordel_r[1], coordel_r[2], coordel_r[3] );
671 [ - - ]: 0 : auto dBdx_r = tk::eval_dBdx_p1( rdof, jacInv_r );
672 : :
673 : : // 2. Average du_i/dx_j
674 : : std::array< std::array< real, 3 >, 3 > dudx;
675 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
676 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
677 : 0 : dudx[i][j] = 0.5 * (
678 : 0 : dBdx_l[j][1] * P(el, velocityDofIdx(nmat,i,rdof,1))
679 : 0 : + dBdx_l[j][2] * P(el, velocityDofIdx(nmat,i,rdof,2))
680 : 0 : + dBdx_l[j][3] * P(el, velocityDofIdx(nmat,i,rdof,3))
681 : 0 : + dBdx_r[j][1] * P(er, velocityDofIdx(nmat,i,rdof,1))
682 : 0 : + dBdx_r[j][2] * P(er, velocityDofIdx(nmat,i,rdof,2))
683 : 0 : + dBdx_r[j][3] * P(er, velocityDofIdx(nmat,i,rdof,3)) );
684 : :
685 : : // 3. average dT/dx_j
686 : : std::vector< std::array< real, 3 > > dTdx(nmat,
687 [ - - ][ - - ]: 0 : std::array< real, 3 >{{0, 0, 0}});
688 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
689 : 0 : auto mark = k*rdof;
690 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
691 : 0 : dTdx[k][j] = 0.5 * (
692 : 0 : dBdx_l[j][1] * T(el, mark+1)
693 : 0 : + dBdx_l[j][2] * T(el, mark+2)
694 : 0 : + dBdx_l[j][3] * T(el, mark+3)
695 : 0 : + dBdx_r[j][1] * T(er, mark+1)
696 : 0 : + dBdx_r[j][2] * T(er, mark+2)
697 : 0 : + dBdx_r[j][3] * T(er, mark+3) );
698 : : }
699 : : }
700 : :
701 : : // 4. Compute flux
702 : : auto fl = modifiedGradientViscousFlux(nmat, ncomp, fn, centroids, state,
703 [ - - ]: 0 : cellAvgState, cellAvgT, dudx, dTdx);
704 : :
705 : : // -------------------------------------------------------------------------
706 : :
707 : : // Add the surface integration term to the rhs
708 [ - - ]: 0 : for (ncomp_t c=0; c<ncomp; ++c)
709 : : {
710 : 0 : R(el, c) += geoFace(f,0) * fl[c];
711 : 0 : R(er, c) -= geoFace(f,0) * fl[c];
712 : : }
713 : : }
714 : 0 : }
715 : :
716 : : std::vector< real >
717 : 0 : modifiedGradientViscousFlux(
718 : : std::size_t nmat,
719 : : std::size_t ncomp,
720 : : const std::array< tk::real, 3 >& fn,
721 : : const std::array< std::array< tk::real, 3 >, 2 >& centroids,
722 : : const std::array< std::vector< tk::real >, 2 >& state,
723 : : const std::array< std::vector< tk::real >, 2 >& cellAvgState,
724 : : const std::array< std::vector< tk::real >, 2 >& cellAvgT,
725 : : const std::array< std::array< real, 3 >, 3 > dudx,
726 : : const std::vector< std::array< real, 3 > >& dTdx )
727 : : // *****************************************************************************
728 : : // Compute the viscous fluxes from the left and right states
729 : : //! \param[in] nmat Number of materials
730 : : //! \param[in] ncomp Number of component equations in the PDE system
731 : : //! \param[in] fn Face/Surface normal
732 : : //! \param[in] centroids Left and right cell centroids
733 : : //! \param[in] state Left and right unknown/state vector
734 : : //! \param[in] cellAvgState Left and right cell-averaged unknown/state vector
735 : : //! \param[in] cellAvgT Left and right cell-averaged temperature vector
736 : : //! \param[in] dudx Average velocity gradient tensor
737 : : //! \param[in] dTdx Average temperature gradient tensor
738 : : //! \return Numerical viscous flux using the Modified Gradient approach.
739 : : //! \details The average gradient is modified according to Weiss et al. to
740 : : //! obtain a stable discretization (average results in unstable central
741 : : //! central difference).
742 : : //! Ref: Weiss, J. M., Maruszewski, J. P., & Smith, W. A. (1999). Implicit
743 : : //! solution of preconditioned Navier-Stokes equations using algebraic
744 : : //! multigrid. AIAA journal, 37(1), 29-36.
745 : : // *****************************************************************************
746 : : {
747 : : using inciter::velocityDofIdx;
748 : : using inciter::volfracDofIdx;
749 : : using inciter::momentumIdx;
750 : : using inciter::velocityIdx;
751 : : using inciter::volfracIdx;
752 : : using inciter::energyIdx;
753 : :
754 : 0 : std::vector< real > fl(ncomp, 0);
755 : :
756 : : // 1. Modify the average gradient
757 : : std::array< real, 3 > r_f{{
758 : 0 : centroids[1][0] - centroids[0][0],
759 : 0 : centroids[1][1] - centroids[0][1],
760 : 0 : centroids[1][2] - centroids[0][2] }};
761 : 0 : real r_mag = std::sqrt(tk::dot(r_f, r_f));
762 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
763 : 0 : r_f[i] /= r_mag;
764 : :
765 : : // velocity gradient
766 : 0 : auto dudx_m = dudx;
767 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
768 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
769 : 0 : dudx_m[i][j] -= ( tk::dot(dudx_m[i], r_f) -
770 : 0 : (cellAvgState[1][ncomp+velocityIdx(nmat,i)]
771 : 0 : - cellAvgState[0][ncomp+velocityIdx(nmat,i)])/r_mag ) * r_f[j];
772 : : // temperature gradient
773 [ - - ]: 0 : auto dTdx_m = dTdx;
774 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
775 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
776 : 0 : dTdx_m[k][j] -= ( tk::dot(dTdx_m[k], r_f) -
777 : 0 : (cellAvgT[1][k] - cellAvgT[0][k])/r_mag ) * r_f[j];
778 : : }
779 : :
780 : : // 2. Compute viscous stress tensor
781 : : std::array< real, 6 > tau;
782 : : real mu(0.0);
783 [ - - ][ - - ]: 0 : std::vector< real > alLR(nmat, 0), conduct_mat(nmat, 0);
[ - - ][ - - ]
784 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k)
785 : 0 : alLR[k] = 0.5*( state[0][volfracIdx(nmat,k)] + state[1][volfracIdx(nmat,k)] );
786 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
787 [ - - ]: 0 : mu += alLR[k] * inciter::getmatprop< tag::mu >(k);
788 [ - - ]: 0 : conduct_mat[k] = inciter::getmatprop< tag::mu >(k) *
789 [ - - ][ - - ]: 0 : inciter::getmatprop< tag::cv >(k) * inciter::getmatprop< tag::gamma >(k)
790 : 0 : / 0.71;
791 : : }
792 : :
793 : 0 : tau[0] = mu * ( 4.0 * dudx_m[0][0] - 2.0*(dudx_m[1][1] + dudx_m[2][2]) ) / 3.0;
794 : 0 : tau[1] = mu * ( 4.0 * dudx_m[1][1] - 2.0*(dudx_m[0][0] + dudx_m[2][2]) ) / 3.0;
795 : 0 : tau[2] = mu * ( 4.0 * dudx_m[2][2] - 2.0*(dudx_m[0][0] + dudx_m[1][1]) ) / 3.0;
796 : 0 : tau[3] = mu * ( dudx_m[0][1] + dudx_m[1][0] );
797 : 0 : tau[4] = mu * ( dudx_m[0][2] + dudx_m[2][0] );
798 : 0 : tau[5] = mu * ( dudx_m[1][2] + dudx_m[2][1] );
799 : :
800 : : // 3. Compute viscous flux terms
801 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
802 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
803 : 0 : fl[momentumIdx(nmat, i)] += tau[inciter::stressCmp[i][j]] * fn[j];
804 : :
805 [ - - ][ - - ]: 0 : std::vector< std::array< real, 3 > > energyFlux(nmat, {{0.0, 0.0, 0.0}});
806 : : std::array< real, 3 > uAvg{{
807 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,0)] + state[1][ncomp+velocityIdx(nmat,0)]),
808 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,1)] + state[1][ncomp+velocityIdx(nmat,1)]),
809 : 0 : 0.5*(state[0][ncomp+velocityIdx(nmat,2)] + state[1][ncomp+velocityIdx(nmat,2)])
810 : 0 : }};
811 : :
812 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
813 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j)
814 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i)
815 : 0 : energyFlux[k][j] += uAvg[i] * tau[inciter::stressCmp[i][j]] +
816 : 0 : conduct_mat[k] * dTdx_m[k][j];
817 : : }
818 : :
819 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
820 : 0 : fl[energyIdx(nmat, k)] = alLR[k]
821 : 0 : * tk::dot(energyFlux[k], fn);
822 : : }
823 : :
824 : 0 : return fl;
825 : : }
826 : :
827 : : } // tk::
|
