Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Physics/Riemann.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Riemann, MUSCL, limiting for edge-based continuous Galerkin
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "Vector.hpp"
14 : : #include "Around.hpp"
15 : : #include "DerivedData.hpp"
16 : : #include "EOS.hpp"
17 : : #include "Riemann.hpp"
18 : : #include "Problems.hpp"
19 : : #include "InciterConfig.hpp"
20 : :
21 : : namespace inciter {
22 : :
23 : : extern ctr::Config g_cfg;
24 : :
25 : : } // ::inciter
26 : :
27 : : namespace riemann {
28 : :
29 : : static const tk::real muscl_eps = 1.0e-9;
30 : : static const tk::real muscl_const = 1.0/3.0;
31 : :
32 : : using inciter::g_cfg;
33 : :
34 : : static void
35 : 17200140 : muscl( std::size_t p,
36 : : std::size_t q,
37 : : const tk::UnsMesh::Coords& coord,
38 : : const tk::Fields& G,
39 : : tk::real& rL, tk::real& uL, tk::real& vL, tk::real& wL, tk::real& eL,
40 : : tk::real& rR, tk::real& uR, tk::real& vR, tk::real& wR, tk::real& eR )
41 : : // *****************************************************************************
42 : : //! Compute MUSCL reconstruction in edge-end points for the flow variables
43 : : //! \param[in] p Left node id of edge-end
44 : : //! \param[in] q Right node id of edge-end
45 : : //! \param[in] coord Array of nodal coordinates
46 : : //! \param[in] G Gradient of all unknowns in mesh points
47 : : //! \param[in,out] rL Left density
48 : : //! \param[in,out] uL Left X velocity
49 : : //! \param[in,out] vL Left Y velocity
50 : : //! \param[in,out] wL Left Z velocity
51 : : //! \param[in,out] eL Left internal energy
52 : : //! \param[in,out] rR Right density
53 : : //! \param[in,out] uR Right X velocity
54 : : //! \param[in,out] vR Right Y velocity
55 : : //! \param[in,out] wR Right Z velocity
56 : : //! \param[in,out] eR Right internal energy
57 : : // *****************************************************************************
58 : : {
59 : : // access node coordinates
60 : 17200140 : const auto& x = coord[0];
61 : 17200140 : const auto& y = coord[1];
62 : 17200140 : const auto& z = coord[2];
63 : :
64 : : // edge vector
65 : 17200140 : tk::real vw[3] = { x[q]-x[p], y[q]-y[p], z[q]-z[p] };
66 : :
67 : : tk::real delta1[5], delta2[5], delta3[5];
68 : 17200140 : tk::real ls[5] = { rL, uL, vL, wL, eL },
69 : 17200140 : rs[5] = { rR, uR, vR, wR, eR },
70 : : url[5], urr[5];
71 : 17200140 : memcpy( url, ls, sizeof ls );
72 : 17200140 : memcpy( urr, rs, sizeof rs );
73 : :
74 : : // MUSCL reconstruction of edge-end-point primitive variables
75 [ + + ]: 103200840 : for (std::size_t c=0; c<5; ++c) {
76 : :
77 [ + - ][ + - ]: 86000700 : auto g1 = G(p,c*3+0)*vw[0] + G(p,c*3+1)*vw[1] + G(p,c*3+2)*vw[2];
[ + - ]
78 [ + - ][ + - ]: 86000700 : auto g2 = G(q,c*3+0)*vw[0] + G(q,c*3+1)*vw[1] + G(q,c*3+2)*vw[2];
[ + - ]
79 : :
80 : 86000700 : delta2[c] = rs[c] - ls[c];
81 : 86000700 : delta1[c] = 2.0 * g1 - delta2[c];
82 : 86000700 : delta3[c] = 2.0 * g2 - delta2[c];
83 : :
84 : : // MUSCL extrapolation option 1:
85 : : // ---------------------------------------------------------------------
86 : : // See Waltz, J., Morgan, N. R., Canfield, T. R., Charest, M. R., Risinger,
87 : : // L. D., & Wohlbier, J. G. (2014). A three-dimensional finite element
88 : : // arbitrary Lagrangian–Eulerian method for shock hydrodynamics on
89 : : // unstructured grids. Computers & Fluids, 92, 172-187.
90 : :
91 : : // van Leer limiter
92 : 86000700 : auto rcL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
93 : 86000700 : auto rcR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
94 : 86000700 : auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
95 : 86000700 : auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
96 : 86000700 : auto phiL = (std::abs(rcL) + rcL) / (std::abs(rcL) + 1.0);
97 : 86000700 : auto phiR = (std::abs(rcR) + rcR) / (std::abs(rcR) + 1.0);
98 : 86000700 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
99 : 86000700 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
100 : : // update unknowns with reconstructed unknowns
101 : 86000700 : url[c] += 0.25*(delta1[c]*(1.0-muscl_const)*phiL +
102 : 86000700 : delta2[c]*(1.0+muscl_const)*phi_L_inv);
103 : 86000700 : urr[c] -= 0.25*(delta3[c]*(1.0-muscl_const)*phiR +
104 : 86000700 : delta2[c]*(1.0+muscl_const)*phi_R_inv);
105 : :
106 : : // MUSCL extrapolation option 2:
107 : : // ---------------------------------------------------------------------
108 : : // See Luo, H., Baum, J. D., & Lohner, R. (1994). Edge-based finite element
109 : : // scheme for the Euler equations. AIAA journal, 32(6), 1183-1190.
110 : : // van Leer, B. (1974). Towards the ultimate conservative difference
111 : : // scheme. II. Monotonicity and conservation combined in a second-order
112 : : // scheme. Journal of computational physics, 14(4), 361-370.
113 : : // Derived from the flux limiter phi as: s = phi_inv - (1 - phi)
114 : :
115 : : // van Albada limiter
116 : : //auto sL = std::max(0.0, (2.0*delta1[c]*delta2[c] + muscl_eps)
117 : : // /(delta1[c]*delta1[c] + delta2[c]*delta2[c] + muscl_eps));
118 : : //auto sR = std::max(0.0, (2.0*delta3[c]*delta2[c] + muscl_eps)
119 : : // /(delta3[c]*delta3[c] + delta2[c]*delta2[c] + muscl_eps));
120 : : //// update unknowns with reconstructed unknowns
121 : : //url[c] += 0.25*sL*(delta1[c]*(1.0 - muscl_const*sL)
122 : : // + delta2[c]*(1.0 + muscl_const*sL));
123 : : //urr[c] -= 0.25*sR*(delta3[c]*(1.0 - muscl_const*sR)
124 : : // + delta2[c]*(1.0 + muscl_const*sR));
125 : : }
126 : :
127 : : // force first order if the reconstructions for density or internal energy
128 : : // would have allowed negative values
129 [ + + ][ + + ]: 17200140 : if (ls[0] < delta1[0] || ls[4] < delta1[4]) memcpy( url, ls, sizeof ls );
130 [ + + ][ + + ]: 17200140 : if (rs[0] < -delta3[0] || rs[4] < -delta3[4]) memcpy( urr, rs, sizeof rs );
131 : :
132 : 17200140 : rL = url[0];
133 : 17200140 : uL = url[1];
134 : 17200140 : vL = url[2];
135 : 17200140 : wL = url[3];
136 : 17200140 : eL = url[4];
137 : :
138 : 17200140 : rR = urr[0];
139 : 17200140 : uR = urr[1];
140 : 17200140 : vR = urr[2];
141 : 17200140 : wR = urr[3];
142 : 17200140 : eR = urr[4];
143 : 17200140 : }
144 : :
145 : : static void
146 : 4668948 : muscl( std::size_t p, std::size_t q, const tk::UnsMesh::Coords& coord,
147 : : const tk::Fields& G, tk::real uL[], tk::real uR[] )
148 : : // *****************************************************************************
149 : : //! Compute MUSCL reconstruction in edge-end points for transported scalars
150 : : //! \param[in] p Left node id of edge-end
151 : : //! \param[in] q Right node id of edge-end
152 : : //! \param[in] coord Array of nodal coordinates
153 : : //! \param[in] G Gradient of all unknowns in mesh points
154 : : //! \param[in,out] uL Primitive variables at left edge-end point
155 : : //! \param[in,out] uR Primitive variables at right edge-end point
156 : : // *****************************************************************************
157 : : {
158 : : #if defined(__clang__)
159 : : #pragma clang diagnostic push
160 : : #pragma clang diagnostic ignored "-Wvla"
161 : : #pragma clang diagnostic ignored "-Wvla-extension"
162 : : #elif defined(STRICT_GNUC)
163 : : #pragma GCC diagnostic push
164 : : #pragma GCC diagnostic ignored "-Wvla"
165 : : #endif
166 : :
167 : 9337896 : auto ns = G.nprop() / 3 - 5;
168 : :
169 : 4668948 : const auto& x = coord[0];
170 : 4668948 : const auto& y = coord[1];
171 : 4668948 : const auto& z = coord[2];
172 : :
173 : : // edge vector
174 : 4668948 : tk::real vw[3] = { x[q]-x[p], y[q]-y[p], z[q]-z[p] };
175 : :
176 : 4668948 : tk::real delta1[ns], delta2[ns], delta3[ns];
177 : :
178 : : // MUSCL reconstruction of edge-end-point of primitive variables
179 [ + + ]: 9337896 : for (std::size_t c=0; c<ns; ++c) {
180 : 4668948 : auto g = (5+c)*3;
181 [ + - ][ + - ]: 4668948 : auto g1 = G(p,g+0)*vw[0] + G(p,g+1)*vw[1] + G(p,g+2)*vw[2];
[ + - ]
182 [ + - ][ + - ]: 4668948 : auto g2 = G(q,g+0)*vw[0] + G(q,g+1)*vw[1] + G(q,g+2)*vw[2];
[ + - ]
183 : :
184 : 4668948 : delta2[c] = uR[5+c] - uL[5+c];
185 : 4668948 : delta1[c] = 2.0 * g1 - delta2[c];
186 : 4668948 : delta3[c] = 2.0 * g2 - delta2[c];
187 : :
188 : : // van Leer limiter
189 : 4668948 : auto rL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
190 : 4668948 : auto rR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
191 : 4668948 : auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
192 : 4668948 : auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
193 : 4668948 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
194 : 4668948 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
195 : 4668948 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
196 : 4668948 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
197 : : // update unknowns with reconstructed unknowns
198 : 4668948 : uL[5+c] += 0.25*(delta1[c]*(1.0-muscl_const)*phiL +
199 : 4668948 : delta2[c]*(1.0+muscl_const)*phi_L_inv);
200 : 4668948 : uR[5+c] -= 0.25*(delta3[c]*(1.0-muscl_const)*phiR +
201 : 4668948 : delta2[c]*(1.0+muscl_const)*phi_R_inv);
202 : : }
203 : :
204 : : #if defined(__clang__)
205 : : #pragma clang diagnostic pop
206 : : #elif defined(STRICT_GNUC)
207 : : #pragma GCC diagnostic pop
208 : : #endif
209 : 9337896 : }
210 : :
211 : : static void
212 : 48947856 : primitive( std::size_t ncomp, std::size_t i, const tk::Fields& U, tk::real u[] )
213 : : // *****************************************************************************
214 : : //! Compute primitive flow variables from conserved ones
215 : : //! \param[in] ncomp Number of scalar components: 5 + scalars
216 : : //! \param[in] i Index to read conserved variables from
217 : : //! \param[in] U Solution vector to read conserved variables from
218 : : //! \param[in,out] u Computed primitive variables
219 : : // *****************************************************************************
220 : : {
221 : 48947856 : u[0] = U(i,0);
222 : 48947856 : u[1] = U(i,1) / u[0];
223 : 48947856 : u[2] = U(i,2) / u[0];
224 : 48947856 : u[3] = U(i,3) / u[0];
225 : 48947856 : u[4] = U(i,4) / u[0] - 0.5*(u[1]*u[1] + u[2]*u[2] + u[3]*u[3]);
226 [ + + ]: 62985900 : for (std::size_t c=5; c<ncomp; ++c) u[c] = U(i,c);
227 : 48947856 : }
228 : :
229 : : void
230 : 38040 : grad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
231 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
232 : : const std::array< std::vector< tk::real >, 3 >& coord,
233 : : const std::vector< std::size_t >& triinpoel,
234 : : const tk::Fields& U,
235 : : tk::Fields& G )
236 : : // *****************************************************************************
237 : : // Compute nodal gradients of primitive variables in all points
238 : : //! \param[in] dsupedge Domain superedges
239 : : //! \param[in] dsupint Domain superedge integrals
240 : : //! \param[in] coord Mesh node coordinates
241 : : //! \param[in] triinpoel Boundary face connectivity
242 : : //! \param[in] U Solution vector at recent time step
243 : : //! \param[in,out] G Nodal gradients
244 : : //! \return Gradients of primitive variables in all mesh points
245 : : // *****************************************************************************
246 : : {
247 : : #if defined(__clang__)
248 : : #pragma clang diagnostic push
249 : : #pragma clang diagnostic ignored "-Wvla"
250 : : #pragma clang diagnostic ignored "-Wvla-extension"
251 : : #elif defined(STRICT_GNUC)
252 : : #pragma GCC diagnostic push
253 : : #pragma GCC diagnostic ignored "-Wvla"
254 : : #endif
255 : :
256 : : // cppcheck-suppress unreadVariable
257 : 38040 : auto ncomp = U.nprop();
258 : :
259 [ - + ][ - - ]: 38040 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
260 [ - + ][ - - ]: 38040 : Assert( G.nprop() == ncomp*3, "Size mismatch" );
[ - - ][ - - ]
261 : 38040 : G.fill( 0.0 );
262 : :
263 : : // domain integral
264 : :
265 : : // domain edge contributions: tetrahedron superedges
266 [ + + ]: 1572747 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
267 : 1534707 : const auto N = dsupedge[0].data() + e*4;
268 : 1534707 : tk::real u[4][ncomp];
269 [ + - ]: 1534707 : primitive( ncomp, N[0], U, u[0] );
270 [ + - ]: 1534707 : primitive( ncomp, N[1], U, u[1] );
271 [ + - ]: 1534707 : primitive( ncomp, N[2], U, u[2] );
272 [ + - ]: 1534707 : primitive( ncomp, N[3], U, u[3] );
273 [ + + ]: 9626205 : for (std::size_t c=0; c<ncomp; ++c) {
274 [ + + ]: 32365992 : for (std::size_t j=0; j<3; ++j) {
275 : : tk::real f[6];
276 : 24274494 : const auto d = dsupint[0].data();
277 : 24274494 : f[0] = d[(e*6+0)*3+j] * (u[1][c] + u[0][c]);
278 : 24274494 : f[1] = d[(e*6+1)*3+j] * (u[2][c] + u[1][c]);
279 : 24274494 : f[2] = d[(e*6+2)*3+j] * (u[0][c] + u[2][c]);
280 : 24274494 : f[3] = d[(e*6+3)*3+j] * (u[3][c] + u[0][c]);
281 : 24274494 : f[4] = d[(e*6+4)*3+j] * (u[3][c] + u[1][c]);
282 : 24274494 : f[5] = d[(e*6+5)*3+j] * (u[3][c] + u[2][c]);
283 [ + - ][ + - ]: 24274494 : G(N[0],c*3+j) = G(N[0],c*3+j) - f[0] + f[2] - f[3];
284 [ + - ][ + - ]: 24274494 : G(N[1],c*3+j) = G(N[1],c*3+j) + f[0] - f[1] - f[4];
285 [ + - ][ + - ]: 24274494 : G(N[2],c*3+j) = G(N[2],c*3+j) + f[1] - f[2] - f[5];
286 [ + - ][ + - ]: 24274494 : G(N[3],c*3+j) = G(N[3],c*3+j) + f[3] + f[4] + f[5];
287 : : }
288 : : }
289 : 1534707 : }
290 : :
291 : : // domain edge contributions: triangle superedges
292 [ + + ]: 1288470 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
293 : 1250430 : const auto N = dsupedge[1].data() + e*3;
294 : 1250430 : tk::real u[3][ncomp];
295 [ + - ]: 1250430 : primitive( ncomp, N[0], U, u[0] );
296 [ + - ]: 1250430 : primitive( ncomp, N[1], U, u[1] );
297 [ + - ]: 1250430 : primitive( ncomp, N[2], U, u[2] );
298 [ + + ]: 7834386 : for (std::size_t c=0; c<ncomp; ++c) {
299 [ + + ]: 26335824 : for (std::size_t j=0; j<3; ++j) {
300 : : tk::real f[3];
301 : 19751868 : const auto d = dsupint[1].data();
302 : 19751868 : f[0] = d[(e*3+0)*3+j] * (u[1][c] + u[0][c]);
303 : 19751868 : f[1] = d[(e*3+1)*3+j] * (u[2][c] + u[1][c]);
304 : 19751868 : f[2] = d[(e*3+2)*3+j] * (u[0][c] + u[2][c]);
305 [ + - ][ + - ]: 19751868 : G(N[0],c*3+j) = G(N[0],c*3+j) - f[0] + f[2];
306 [ + - ][ + - ]: 19751868 : G(N[1],c*3+j) = G(N[1],c*3+j) + f[0] - f[1];
307 [ + - ][ + - ]: 19751868 : G(N[2],c*3+j) = G(N[2],c*3+j) + f[1] - f[2];
308 : : }
309 : : }
310 : 1250430 : }
311 : :
312 : : // domain edge contributions: edges
313 [ + + ]: 4278648 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
314 : 4240608 : const auto N = dsupedge[2].data() + e*2;
315 : 4240608 : tk::real u[2][ncomp];
316 [ + - ]: 4240608 : primitive( ncomp, N[0], U, u[0] );
317 [ + - ]: 4240608 : primitive( ncomp, N[1], U, u[1] );
318 : 4240608 : const auto d = dsupint[2].data() + e*3;
319 [ + + ]: 26609400 : for (std::size_t c=0; c<ncomp; ++c) {
320 [ + + ]: 89475168 : for (std::size_t j=0; j<3; ++j) {
321 : 67106376 : tk::real f = d[j] * (u[1][c] + u[0][c]);
322 [ + - ]: 67106376 : G(N[0],c*3+j) -= f;
323 [ + - ]: 67106376 : G(N[1],c*3+j) += f;
324 : : }
325 : : }
326 : 4240608 : }
327 : :
328 : : // boundary integral
329 : :
330 : 38040 : const auto& x = coord[0];
331 : 38040 : const auto& y = coord[1];
332 : 38040 : const auto& z = coord[2];
333 : :
334 [ + + ]: 4106436 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
335 : 8136792 : const auto N = triinpoel.data() + e*3;
336 : 4068396 : tk::real u[3][ncomp];
337 [ + - ]: 4068396 : primitive( ncomp, N[0], U, u[0] );
338 [ + - ]: 4068396 : primitive( ncomp, N[1], U, u[1] );
339 [ + - ]: 4068396 : primitive( ncomp, N[2], U, u[2] );
340 : : const std::array< tk::real, 3 >
341 : 4068396 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
342 : 4068396 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
343 : 4068396 : auto n = tk::cross( ba, ca );
344 : 4068396 : n[0] /= 12.0;
345 : 4068396 : n[1] /= 12.0;
346 : 4068396 : n[2] /= 12.0;
347 [ + + ]: 25757208 : for (std::size_t c=0; c<ncomp; ++c) {
348 : 21688812 : auto uab = (u[0][c] + u[1][c])/4.0;
349 : 21688812 : auto ubc = (u[1][c] + u[2][c])/4.0;
350 : 21688812 : auto uca = (u[2][c] + u[0][c])/4.0;
351 : 21688812 : tk::real g[] = { uab + uca + u[0][c],
352 : 21688812 : uab + ubc + u[1][c],
353 : 21688812 : ubc + uca + u[2][c] };
354 [ + + ]: 86755248 : for (std::size_t j=0; j<3; ++j) {
355 [ + - ]: 65066436 : G(N[0],c*3+j) += g[j] * n[j];
356 [ + - ]: 65066436 : G(N[1],c*3+j) += g[j] * n[j];
357 [ + - ]: 65066436 : G(N[2],c*3+j) += g[j] * n[j];
358 : : }
359 : : }
360 : 4068396 : }
361 : :
362 : : #if defined(__clang__)
363 : : #pragma clang diagnostic pop
364 : : #elif defined(STRICT_GNUC)
365 : : #pragma GCC diagnostic pop
366 : : #endif
367 : 38040 : }
368 : :
369 : : static void
370 : 16246428 : rusanov( const tk::UnsMesh::Coords& coord,
371 : : const tk::Fields& G,
372 : : const tk::real dsupint[],
373 : : std::size_t p,
374 : : std::size_t q,
375 : : const tk::real L[],
376 : : const tk::real R[],
377 : : tk::real f[],
378 : : std::size_t symL,
379 : : std::size_t symR )
380 : : // *****************************************************************************
381 : : //! Compute advection fluxes on a single edge with Rusanov's flux
382 : : //! \param[in] coord Mesh node coordinates
383 : : //! \param[in] G Nodal gradients
384 : : //! \param[in] dsupint Domain superedge integral for this edge
385 : : //! \param[in] p Left node index of edge
386 : : //! \param[in] q Right node index of edge
387 : : //! \param[in,out] L Left physics state variables
388 : : //! \param[in,out] R Rigth physics state variables
389 : : //! \param[in,out] f Flux computed
390 : : //! \param[in] symL Non-zero if left edge end-point is on a symmetry boundary
391 : : //! \param[in] symR Non-zero if right edge end-point is on a symmetry boundary
392 : : // *****************************************************************************
393 : : {
394 : : #if defined(__clang__)
395 : : #pragma clang diagnostic push
396 : : #pragma clang diagnostic ignored "-Wvla"
397 : : #pragma clang diagnostic ignored "-Wvla-extension"
398 : : #elif defined(STRICT_GNUC)
399 : : #pragma GCC diagnostic push
400 : : #pragma GCC diagnostic ignored "-Wvla"
401 : : #endif
402 : :
403 : 32492856 : auto ncomp = G.nprop() / 3;
404 : :
405 : : // will work on copies of physics variables
406 : 16246428 : tk::real l[ncomp], r[ncomp];
407 : 16246428 : memcpy( l, L, sizeof l );
408 : 16246428 : memcpy( r, R, sizeof r );
409 : :
410 : : // MUSCL reconstruction in edge-end points for flow variables
411 : 16246428 : muscl( p, q, coord, G, l[0], l[1], l[2], l[3], l[4],
412 [ + - ]: 16246428 : r[0], r[1], r[2], r[3], r[4] );
413 : :
414 : : // pressure
415 : 16246428 : auto pL = eos::pressure( l[0]*l[4] );
416 : 16246428 : auto pR = eos::pressure( r[0]*r[4] );
417 : :
418 : : // dualface-normal velocities
419 : 16246428 : auto nx = dsupint[0];
420 : 16246428 : auto ny = dsupint[1];
421 : 16246428 : auto nz = dsupint[2];
422 [ + - ]: 16246428 : auto vnL = symL ? 0.0 : (l[1]*nx + l[2]*ny + l[3]*nz);
423 [ + - ]: 16246428 : auto vnR = symR ? 0.0 : (r[1]*nx + r[2]*ny + r[3]*nz);
424 : :
425 : : // back to conserved variables
426 : 16246428 : l[4] = (l[4] + 0.5*(l[1]*l[1] + l[2]*l[2] + l[3]*l[3])) * l[0];
427 : 16246428 : l[1] *= l[0];
428 : 16246428 : l[2] *= l[0];
429 : 16246428 : l[3] *= l[0];
430 : 16246428 : r[4] = (r[4] + 0.5*(r[1]*r[1] + r[2]*r[2] + r[3]*r[3])) * r[0];
431 : 16246428 : r[1] *= r[0];
432 : 16246428 : r[2] *= r[0];
433 : 16246428 : r[3] *= r[0];
434 : :
435 : : // dissipation
436 : 16246428 : auto len = tk::length( nx, ny, nz );
437 : 16246428 : auto sl = std::abs(vnL) + eos::soundspeed(l[0],pL)*len;
438 : 16246428 : auto sr = std::abs(vnR) + eos::soundspeed(r[0],pR)*len;
439 : 16246428 : auto fw = std::max( sl, sr );
440 : :
441 : : // flow fluxes
442 : 16246428 : f[0] = l[0]*vnL + r[0]*vnR + fw*(r[0] - l[0]);
443 : 16246428 : f[1] = l[1]*vnL + r[1]*vnR + (pL + pR)*nx + fw*(r[1] - l[1]);
444 : 16246428 : f[2] = l[2]*vnL + r[2]*vnR + (pL + pR)*ny + fw*(r[2] - l[2]);
445 : 16246428 : f[3] = l[3]*vnL + r[3]*vnR + (pL + pR)*nz + fw*(r[3] - l[3]);
446 : 16246428 : f[4] = (l[4] + pL)*vnL + (r[4] + pR)*vnR + fw*(r[4] - l[4]);
447 : :
448 : : // artificial viscosity
449 : 16246428 : const auto stab2 = g_cfg.get< tag::stab2 >();
450 [ + + ]: 16246428 : if (stab2) {
451 : 242604 : auto stab2coef = g_cfg.get< tag::stab2coef >();
452 : 242604 : auto fws = stab2coef * fw;
453 : 242604 : f[0] -= fws*(l[0] - r[0]);
454 : 242604 : f[1] -= fws*(l[1] - r[1]);
455 : 242604 : f[2] -= fws*(l[2] - r[2]);
456 : 242604 : f[3] -= fws*(l[3] - r[3]);
457 : 242604 : f[4] -= fws*(l[4] - r[4]);
458 : : }
459 : :
460 [ + + ]: 16246428 : if (ncomp == 5) return;
461 : :
462 : : // MUSCL reconstruction in edge-end points for scalars
463 [ + - ]: 4306428 : muscl( p, q, coord, G, l, r );
464 : :
465 : : // scalar dissipation
466 : 4306428 : auto sw = std::max( std::abs(vnL), std::abs(vnR) );
467 : :
468 : : // scalar fluxes
469 [ + + ]: 8612856 : for (std::size_t c=5; c<ncomp; ++c) {
470 : 4306428 : f[c] = l[c]*vnL + r[c]*vnR + sw*(r[c] - l[c]);
471 : : }
472 : :
473 : : #if defined(__clang__)
474 : : #pragma clang diagnostic pop
475 : : #elif defined(STRICT_GNUC)
476 : : #pragma GCC diagnostic pop
477 : : #endif
478 : 16246428 : }
479 : :
480 : : static void
481 : 953712 : hllc( const tk::UnsMesh::Coords& coord,
482 : : const tk::Fields& G,
483 : : const tk::real dsupint[],
484 : : std::size_t p,
485 : : std::size_t q,
486 : : const tk::real L[],
487 : : const tk::real R[],
488 : : tk::real f[],
489 : : std::size_t symL,
490 : : std::size_t symR )
491 : : // *****************************************************************************
492 : : //! Compute advection fluxes on a single edge with Harten-Lax-vanLeer-Contact
493 : : //! \param[in] coord Mesh node coordinates
494 : : //! \param[in] G Nodal gradients
495 : : //! \param[in] dsupint Domain superedge integral for this edge
496 : : //! \param[in] p Left node index of edge
497 : : //! \param[in] q Right node index of edge
498 : : //! \param[in,out] L Left physics state variables
499 : : //! \param[in,out] R Rigth physics state variables
500 : : //! \param[in,out] f Flux computed
501 : : //! \param[in] symL Non-zero if left edge end-point is on a symmetry boundary
502 : : //! \param[in] symR Non-zero if right edge end-point is on a symmetry boundary
503 : : // *****************************************************************************
504 : : {
505 : : #if defined(__clang__)
506 : : #pragma clang diagnostic push
507 : : #pragma clang diagnostic ignored "-Wvla"
508 : : #pragma clang diagnostic ignored "-Wvla-extension"
509 : : #elif defined(STRICT_GNUC)
510 : : #pragma GCC diagnostic push
511 : : #pragma GCC diagnostic ignored "-Wvla"
512 : : #endif
513 : :
514 : 1907424 : auto ncomp = G.nprop() / 3;
515 : :
516 : : // will work on copies of physics variables
517 : 953712 : tk::real l[ncomp], r[ncomp];
518 : 953712 : memcpy( l, L, sizeof l );
519 : 953712 : memcpy( r, R, sizeof r );
520 : :
521 : : // MUSCL reconstruction in edge-end points for flow variables
522 : 953712 : muscl( p, q, coord, G, l[0], l[1], l[2], l[3], l[4],
523 [ + - ]: 953712 : r[0], r[1], r[2], r[3], r[4] );
524 : :
525 : : // dualface-normal velocities
526 : 953712 : auto nx = -dsupint[0];
527 : 953712 : auto ny = -dsupint[1];
528 : 953712 : auto nz = -dsupint[2];
529 : 953712 : auto len = tk::length( nx, ny, nz );
530 : 953712 : nx /= len;
531 : 953712 : ny /= len;
532 : 953712 : nz /= len;
533 [ + - ]: 953712 : auto qL = symL ? 0.0 : (l[1]*nx + l[2]*ny + l[3]*nz);
534 [ + - ]: 953712 : auto qR = symR ? 0.0 : (r[1]*nx + r[2]*ny + r[3]*nz);
535 : :
536 : : // pressure
537 : 953712 : auto pL = eos::pressure( l[0]*l[4] );
538 : 953712 : auto pR = eos::pressure( r[0]*r[4] );
539 : :
540 : : // back to conserved variables
541 : 953712 : l[4] = (l[4] + 0.5*(l[1]*l[1] + l[2]*l[2] + l[3]*l[3])) * l[0];
542 : 953712 : l[1] *= l[0];
543 : 953712 : l[2] *= l[0];
544 : 953712 : l[3] *= l[0];
545 : 953712 : r[4] = (r[4] + 0.5*(r[1]*r[1] + r[2]*r[2] + r[3]*r[3])) * r[0];
546 : 953712 : r[1] *= r[0];
547 : 953712 : r[2] *= r[0];
548 : 953712 : r[3] *= r[0];
549 : :
550 : : // sound speed
551 : 953712 : auto cL = eos::soundspeed(l[0],pL);
552 : 953712 : auto cR = eos::soundspeed(r[0],pR);
553 : :
554 : : // left and right wave speeds
555 : 953712 : auto sL = fmin( qL - cL, qR - cR );
556 : 953712 : auto sR = fmax( qL + cL, qR + cR );
557 : :
558 : : // contact wave speed and pressure
559 : 953712 : auto tL = sL - qL;
560 : 953712 : auto tR = sR - qR;
561 : 953712 : auto sM = (r[0]*qR*tR - l[0]*qL*tL + pL - pR) / (r[0]*tR - l[0]*tL);
562 : 953712 : auto pS = pL - l[0]*tL*(qL - sM);
563 : :
564 : : // intermediate left-, and right-state conserved unknowns
565 : 953712 : tk::real uL[ncomp], uR[ncomp];
566 : 953712 : auto s = sL - sM;
567 : 953712 : uL[0] = tL*l[0]/s;
568 : 953712 : uL[1] = (tL*l[1] + (pS-pL)*nx)/s;
569 : 953712 : uL[2] = (tL*l[2] + (pS-pL)*ny)/s;
570 : 953712 : uL[3] = (tL*l[3] + (pS-pL)*nz)/s;
571 : 953712 : uL[4] = (tL*l[4] - pL*qL + pS*sM)/s;
572 : 953712 : s = sR - sM;
573 : 953712 : uR[0] = tR*r[0]/s;
574 : 953712 : uR[1] = (tR*r[1] + (pS-pR)*nx)/s;
575 : 953712 : uR[2] = (tR*r[2] + (pS-pR)*ny)/s;
576 : 953712 : uR[3] = (tR*r[3] + (pS-pR)*nz)/s;
577 : 953712 : uR[4] = (tR*r[4] - pR*qR + pS*sM)/s;
578 : :
579 : 953712 : auto L2 = -2.0*len;
580 : 953712 : nx *= L2;
581 : 953712 : ny *= L2;
582 : 953712 : nz *= L2;
583 : :
584 : : // flow fluxes
585 [ - + ]: 953712 : if (sL > 0.0) {
586 : 0 : auto qL2 = qL * L2;
587 : 0 : f[0] = l[0]*qL2;
588 : 0 : f[1] = l[1]*qL2 + pL*nx;
589 : 0 : f[2] = l[2]*qL2 + pL*ny;
590 : 0 : f[3] = l[3]*qL2 + pL*nz;
591 : 0 : f[4] = (l[4] + pL)*qL2;
592 : : }
593 [ + - ][ + + ]: 953712 : else if (sL <= 0.0 && sM > 0.0) {
594 : 477633 : auto qL2 = qL * L2;
595 : 477633 : auto sL2 = sL * L2;
596 : 477633 : f[0] = l[0]*qL2 + sL2*(uL[0] - l[0]);
597 : 477633 : f[1] = l[1]*qL2 + pL*nx + sL2*(uL[1] - l[1]);
598 : 477633 : f[2] = l[2]*qL2 + pL*ny + sL2*(uL[2] - l[2]);
599 : 477633 : f[3] = l[3]*qL2 + pL*nz + sL2*(uL[3] - l[3]);
600 : 477633 : f[4] = (l[4] + pL)*qL2 + sL2*(uL[4] - l[4]);
601 : 477633 : }
602 [ + - ][ + - ]: 476079 : else if (sM <= 0.0 && sR >= 0.0) {
603 : 476079 : auto qR2 = qR * L2;
604 : 476079 : auto sR2 = sR * L2;
605 : 476079 : f[0] = r[0]*qR2 + sR2*(uR[0] - r[0]);
606 : 476079 : f[1] = r[1]*qR2 + pR*nx + sR2*(uR[1] - r[1]);
607 : 476079 : f[2] = r[2]*qR2 + pR*ny + sR2*(uR[2] - r[2]);
608 : 476079 : f[3] = r[3]*qR2 + pR*nz + sR2*(uR[3] - r[3]);
609 : 476079 : f[4] = (r[4] + pR)*qR2 + sR2*(uR[4] - r[4]);
610 : 476079 : }
611 : : else {
612 : 0 : auto qR2 = qR * L2;
613 : 0 : f[0] = r[0]*qR2;
614 : 0 : f[1] = r[1]*qR2 + pR*nx;
615 : 0 : f[2] = r[2]*qR2 + pR*ny;
616 : 0 : f[3] = r[3]*qR2 + pR*nz;
617 : 0 : f[4] = (r[4] + pR)*qR2;
618 : : }
619 : :
620 : : // artificial viscosity
621 : 953712 : const auto stab2 = g_cfg.get< tag::stab2 >();
622 [ + + ]: 953712 : if (stab2) {
623 : 295596 : auto stab2coef = g_cfg.get< tag::stab2coef >();
624 : 295596 : auto sl = std::abs(qL) + cL;
625 : 295596 : auto sr = std::abs(qR) + cR;
626 : 295596 : auto fws = stab2coef * std::max(sl,sr) * len;
627 : 295596 : f[0] -= fws * (l[0] - r[0]);
628 : 295596 : f[1] -= fws * (l[1] - r[1]);
629 : 295596 : f[2] -= fws * (l[2] - r[2]);
630 : 295596 : f[3] -= fws * (l[3] - r[3]);
631 : 295596 : f[4] -= fws * (l[4] - r[4]);
632 : : }
633 : :
634 [ + + ]: 953712 : if (ncomp == 5) return;
635 : :
636 : : // MUSCL reconstruction in edge-end points for scalars
637 [ + - ]: 362520 : muscl( p, q, coord, G, l, r );
638 : :
639 : : // scalar
640 : 362520 : auto sw = std::max( std::abs(qL), std::abs(qR) ) * len;
641 [ + + ]: 725040 : for (std::size_t c=5; c<ncomp; ++c) {
642 : 362520 : f[c] = (l[c]*qL + r[c]*qR)*len + sw*(r[c] - l[c]);
643 : : }
644 : :
645 : : #if defined(__clang__)
646 : : #pragma clang diagnostic pop
647 : : #elif defined(STRICT_GNUC)
648 : : #pragma GCC diagnostic pop
649 : : #endif
650 : 953712 : }
651 : :
652 : : static void
653 : 38040 : advdom( const tk::UnsMesh::Coords& coord,
654 : : const std::array< std::vector< std::size_t >, 3 >& dsupedge,
655 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
656 : : const tk::Fields& G,
657 : : const tk::Fields& U,
658 : : // cppcheck-suppress constParameter
659 : : tk::Fields& R )
660 : : // *****************************************************************************
661 : : //! Compute domain integral for advection
662 : : //! \param[in] coord Mesh node coordinates
663 : : //! \param[in] dsupedge Domain superedges
664 : : //! \param[in] dsupint Domain superedge integrals
665 : : //! \param[in] G Nodal gradients
666 : : //! \param[in] U Solution vector at recent time step
667 : : //! \param[in,out] R Right-hand side vector computed
668 : : // *****************************************************************************
669 : : {
670 : : // number of transported scalars
671 : 38040 : auto ncomp = U.nprop();
672 : :
673 : : // configure flux function
674 : 38040 : auto cfgflux = [](){
675 : 38040 : const auto& cfg = g_cfg.get< tag::flux >();
676 [ + + ]: 38040 : if (cfg == "rusanov")
677 : 34668 : return rusanov;
678 [ + - ]: 3372 : else if (cfg == "hllc")
679 : 3372 : return hllc;
680 : : else
681 [ - - ][ - - ]: 0 : Throw( "Flux not configured" );
[ - - ]
682 : : };
683 : :
684 [ + - ]: 38040 : auto flux = cfgflux();
685 : :
686 : : #if defined(__clang__)
687 : : #pragma clang diagnostic push
688 : : #pragma clang diagnostic ignored "-Wvla"
689 : : #pragma clang diagnostic ignored "-Wvla-extension"
690 : : #elif defined(STRICT_GNUC)
691 : : #pragma GCC diagnostic push
692 : : #pragma GCC diagnostic ignored "-Wvla"
693 : : #endif
694 : :
695 : : // domain edge contributions: tetrahedron superedges
696 [ + + ]: 1572747 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
697 : 1534707 : const auto N = dsupedge[0].data() + e*4;
698 : : // primitive variables
699 : 1534707 : tk::real u[4][ncomp];
700 [ + - ]: 1534707 : primitive( ncomp, N[0], U, u[0] );
701 [ + - ]: 1534707 : primitive( ncomp, N[1], U, u[1] );
702 [ + - ]: 1534707 : primitive( ncomp, N[2], U, u[2] );
703 [ + - ]: 1534707 : primitive( ncomp, N[3], U, u[3] );
704 : : // edge fluxes
705 : 1534707 : tk::real f[6][ncomp];
706 : 1534707 : const auto d = dsupint[0].data();
707 [ + - ]: 1534707 : flux( coord, G, d+(e*6+0)*3, N[0], N[1], u[0], u[1], f[0], 0, 0 );
708 [ + - ]: 1534707 : flux( coord, G, d+(e*6+1)*3, N[1], N[2], u[1], u[2], f[1], 0, 0 );
709 [ + - ]: 1534707 : flux( coord, G, d+(e*6+2)*3, N[2], N[0], u[2], u[0], f[2], 0, 0 );
710 [ + - ]: 1534707 : flux( coord, G, d+(e*6+3)*3, N[0], N[3], u[0], u[3], f[3], 0, 0 );
711 [ + - ]: 1534707 : flux( coord, G, d+(e*6+4)*3, N[1], N[3], u[1], u[3], f[4], 0, 0 );
712 [ + - ]: 1534707 : flux( coord, G, d+(e*6+5)*3, N[2], N[3], u[2], u[3], f[5], 0, 0 );
713 : : // edge flux contributions
714 [ + + ]: 9626205 : for (std::size_t c=0; c<ncomp; ++c) {
715 [ + - ][ + - ]: 8091498 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
716 [ + - ][ + - ]: 8091498 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
717 [ + - ][ + - ]: 8091498 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
718 [ + - ][ + - ]: 8091498 : R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
719 : : }
720 : 1534707 : }
721 : :
722 : : // domain edge contributions: triangle superedges
723 [ + + ]: 1288470 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
724 : 1250430 : const auto N = dsupedge[1].data() + e*3;
725 : : // primitive variables
726 : 1250430 : tk::real u[3][ncomp];
727 [ + - ]: 1250430 : primitive( ncomp, N[0], U, u[0] );
728 [ + - ]: 1250430 : primitive( ncomp, N[1], U, u[1] );
729 [ + - ]: 1250430 : primitive( ncomp, N[2], U, u[2] );
730 : : // edge fluxes
731 : 1250430 : tk::real f[3][ncomp];
732 : 1250430 : const auto d = dsupint[1].data();
733 [ + - ]: 1250430 : flux( coord, G, d+(e*3+0)*3, N[0], N[1], u[0], u[1], f[0], 0, 0 );
734 [ + - ]: 1250430 : flux( coord, G, d+(e*3+1)*3, N[1], N[2], u[1], u[2], f[1], 0, 0 );
735 [ + - ]: 1250430 : flux( coord, G, d+(e*3+2)*3, N[2], N[0], u[2], u[0], f[2], 0, 0 );
736 : : // edge flux contributions
737 [ + + ]: 7834386 : for (std::size_t c=0; c<ncomp; ++c) {
738 [ + - ][ + - ]: 6583956 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
739 [ + - ][ + - ]: 6583956 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
740 [ + - ][ + - ]: 6583956 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
741 : : }
742 : 1250430 : }
743 : :
744 : : // domain edge contributions: edges
745 [ + + ]: 4278648 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
746 : 4240608 : const auto N = dsupedge[2].data() + e*2;
747 : 4240608 : tk::real u[2][ncomp];
748 [ + - ]: 4240608 : primitive( ncomp, N[0], U, u[0] );
749 [ + - ]: 4240608 : primitive( ncomp, N[1], U, u[1] );
750 : : // edge fluxes
751 : 4240608 : tk::real f[ncomp];
752 : 4240608 : const auto d = dsupint[2].data();
753 [ + - ]: 4240608 : flux( coord, G, d+e*3, N[0], N[1], u[0], u[1], f, 0, 0 );
754 : : // edge flux contributions
755 [ + + ]: 26609400 : for (std::size_t c=0; c<ncomp; ++c) {
756 [ + - ]: 22368792 : R(N[0],c) -= f[c];
757 [ + - ]: 22368792 : R(N[1],c) += f[c];
758 : : }
759 : 4240608 : }
760 : :
761 : : #if defined(__clang__)
762 : : #pragma clang diagnostic pop
763 : : #elif defined(STRICT_GNUC)
764 : : #pragma GCC diagnostic pop
765 : : #endif
766 : 38040 : }
767 : :
768 : : static void
769 : 38040 : advbnd( const std::vector< std::size_t >& triinpoel,
770 : : const std::array< std::vector< tk::real >, 3 >& coord,
771 : : const std::vector< std::uint8_t >& besym,
772 : : const tk::Fields& U,
773 : : tk::Fields& R )
774 : : // *****************************************************************************
775 : : //! Compute boundary integral for advection
776 : : //! \param[in] triinpoel Boundary face connectivity
777 : : //! \param[in] coord Mesh node coordinates
778 : : //! \param[in] besym Boundary element symmetry BC flags
779 : : //! \param[in] U Solution vector at recent time step
780 : : //! \param[in,out] R Right-hand side vector
781 : : // *****************************************************************************
782 : : {
783 : 38040 : auto ncomp = U.nprop();
784 : :
785 : 38040 : const auto& x = coord[0];
786 : 38040 : const auto& y = coord[1];
787 : 38040 : const auto& z = coord[2];
788 : :
789 : : #if defined(__clang__)
790 : : #pragma clang diagnostic push
791 : : #pragma clang diagnostic ignored "-Wvla"
792 : : #pragma clang diagnostic ignored "-Wvla-extension"
793 : : #elif defined(STRICT_GNUC)
794 : : #pragma GCC diagnostic push
795 : : #pragma GCC diagnostic ignored "-Wvla"
796 : : #endif
797 : :
798 [ + + ]: 4106436 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
799 : 8136792 : const auto N = triinpoel.data() + e*3;
800 : :
801 [ + - ]: 4068396 : auto rA = U(N[0],0);
802 [ + - ]: 4068396 : auto ruA = U(N[0],1);
803 [ + - ]: 4068396 : auto rvA = U(N[0],2);
804 [ + - ]: 4068396 : auto rwA = U(N[0],3);
805 [ + - ]: 4068396 : auto reA = U(N[0],4);
806 : :
807 [ + - ]: 4068396 : auto rB = U(N[1],0);
808 [ + - ]: 4068396 : auto ruB = U(N[1],1);
809 [ + - ]: 4068396 : auto rvB = U(N[1],2);
810 [ + - ]: 4068396 : auto rwB = U(N[1],3);
811 [ + - ]: 4068396 : auto reB = U(N[1],4);
812 : :
813 [ + - ]: 4068396 : auto rC = U(N[2],0);
814 [ + - ]: 4068396 : auto ruC = U(N[2],1);
815 [ + - ]: 4068396 : auto rvC = U(N[2],2);
816 [ + - ]: 4068396 : auto rwC = U(N[2],3);
817 [ + - ]: 4068396 : auto reC = U(N[2],4);
818 : :
819 : : const std::array< tk::real, 3 >
820 : 4068396 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
821 : 4068396 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
822 : 4068396 : auto [nx,ny,nz] = tk::cross( ba, ca ); // 2A
823 : 4068396 : nx /= 12.0;
824 : 4068396 : ny /= 12.0;
825 : 4068396 : nz /= 12.0;
826 : :
827 : 4068396 : tk::real p, vn, f[ncomp][3];
828 : 4068396 : const auto sym = besym.data() + e*3;
829 : :
830 : 4068396 : p = eos::pressure( reA - 0.5*(ruA*ruA + rvA*rvA + rwA*rwA)/rA );
831 [ + + ]: 4068396 : vn = sym[0] ? 0.0 : (nx*ruA + ny*rvA + nz*rwA)/rA;
832 : : // flow
833 : 4068396 : f[0][0] = rA*vn;
834 : 4068396 : f[1][0] = ruA*vn + p*nx;
835 : 4068396 : f[2][0] = rvA*vn + p*ny;
836 : 4068396 : f[3][0] = rwA*vn + p*nz;
837 : 4068396 : f[4][0] = (reA + p)*vn;
838 : : // scalar
839 [ + - ][ + + ]: 5415228 : for (std::size_t c=5; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;
840 : :
841 : 4068396 : p = eos::pressure( reB - 0.5*(ruB*ruB + rvB*rvB + rwB*rwB)/rB );
842 [ + + ]: 4068396 : vn = sym[1] ? 0.0 : (nx*ruB + ny*rvB + nz*rwB)/rB;
843 : : // flow
844 : 4068396 : f[0][1] = rB*vn;
845 : 4068396 : f[1][1] = ruB*vn + p*nx;
846 : 4068396 : f[2][1] = rvB*vn + p*ny;
847 : 4068396 : f[3][1] = rwB*vn + p*nz;
848 : 4068396 : f[4][1] = (reB + p)*vn;
849 : : // scalar
850 [ + - ][ + + ]: 5415228 : for (std::size_t c=5; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;
851 : :
852 : 4068396 : p = eos::pressure( reC - 0.5*(ruC*ruC + rvC*rvC + rwC*rwC)/rC );
853 [ + + ]: 4068396 : vn = sym[2] ? 0.0 : (nx*ruC + ny*rvC + nz*rwC)/rC;
854 : : // flow
855 : 4068396 : f[0][2] = rC*vn;
856 : 4068396 : f[1][2] = ruC*vn + p*nx;
857 : 4068396 : f[2][2] = rvC*vn + p*ny;
858 : 4068396 : f[3][2] = rwC*vn + p*nz;
859 : 4068396 : f[4][2] = (reC + p)*vn;
860 : : // scalar
861 [ + - ][ + + ]: 5415228 : for (std::size_t c=5; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;
862 : :
863 [ + + ]: 25757208 : for (std::size_t c=0; c<ncomp; ++c) {
864 : 21688812 : auto fab = (f[c][0] + f[c][1])/4.0;
865 : 21688812 : auto fbc = (f[c][1] + f[c][2])/4.0;
866 : 21688812 : auto fca = (f[c][2] + f[c][0])/4.0;
867 [ + - ]: 21688812 : R(N[0],c) += fab + fca + f[c][0];
868 [ + - ]: 21688812 : R(N[1],c) += fab + fbc + f[c][1];
869 [ + - ]: 21688812 : R(N[2],c) += fbc + fca + f[c][2];
870 : : }
871 : 4068396 : }
872 : :
873 : : #if defined(__clang__)
874 : : #pragma clang diagnostic pop
875 : : #elif defined(STRICT_GNUC)
876 : : #pragma GCC diagnostic pop
877 : : #endif
878 : 38040 : }
879 : :
880 : : static void
881 : 38040 : src( const std::array< std::vector< tk::real >, 3 >& coord,
882 : : const std::vector< tk::real >& v,
883 : : tk::real t,
884 : : const std::vector< tk::real >& tp,
885 : : tk::Fields& R )
886 : : // *****************************************************************************
887 : : // Compute source integral
888 : : //! \param[in] coord Mesh node coordinates
889 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
890 : : //! \param[in] t Physical time
891 : : //! \param[in] tp Physical time for each mesh node
892 : : //! \param[in,out] R Right-hand side vector computed
893 : : // *****************************************************************************
894 : : {
895 [ + - ]: 38040 : auto src = problems::SRC();
896 [ + + ]: 38040 : if (!src) return;
897 : :
898 : 23400 : const auto& x = coord[0];
899 : 23400 : const auto& y = coord[1];
900 : 23400 : const auto& z = coord[2];
901 : :
902 [ + + ]: 1622535 : for (std::size_t p=0; p<R.nunk(); ++p) {
903 [ + + ]: 1599135 : if (g_cfg.get< tag::steady >()) t = tp[p];
904 [ + - ]: 1599135 : auto s = src( x[p], y[p], z[p], t );
905 [ + - ][ + + ]: 10250940 : for (std::size_t c=0; c<s.size(); ++c) R(p,c) -= s[c] * v[p];
906 : 1599135 : }
907 [ + + ]: 38040 : }
908 : :
909 : : void
910 : 38040 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
911 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
912 : : const std::array< std::vector< tk::real >, 3 >& coord,
913 : : const std::vector< std::size_t >& triinpoel,
914 : : const std::vector< std::uint8_t >& besym,
915 : : const tk::Fields& G,
916 : : const tk::Fields& U,
917 : : const std::vector< tk::real >& v,
918 : : tk::real t,
919 : : const std::vector< tk::real >& tp,
920 : : tk::Fields& R )
921 : : // *****************************************************************************
922 : : // Compute right hand side
923 : : //! \param[in] dsupedge Domain superedges
924 : : //! \param[in] dsupint Domain superedge integrals
925 : : //! \param[in] coord Mesh node coordinates
926 : : //! \param[in] triinpoel Boundary face connectivity
927 : : //! \param[in] besym Boundary element symmetry BC flags
928 : : //! \param[in] G Gradients in mesh nodes
929 : : //! \param[in] U Unknowns/solution vector in mesh nodes
930 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
931 : : //! \param[in] t Physical time
932 : : //! \param[in] tp Physical time for each mesh node
933 : : //! \param[in,out] R Right-hand side vector computed
934 : : // *****************************************************************************
935 : : {
936 [ - + ][ - - ]: 38040 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
937 : : "vector at recent time step incorrect" );
938 [ - + ][ - - ]: 38040 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
939 : : "Number of unknowns and/or number of components in right-hand "
940 : : "side vector incorrect" );
941 : :
942 : 38040 : R.fill( 0.0 );
943 : 38040 : advdom( coord, dsupedge, dsupint, G, U, R );
944 : 38040 : advbnd( triinpoel, coord, besym, U, R );
945 : 38040 : src( coord, v, t, tp, R );
946 : 38040 : }
947 : :
948 : : } // riemann::
|
