Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Physics/Lohner.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 LohCG: Artificial compressibility solver for incompressible flow
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "Vector.hpp"
14 : : #include "Around.hpp"
15 : : #include "DerivedData.hpp"
16 : : #include "EOS.hpp"
17 : : #include "Lohner.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 lohner {
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 : : void
35 : 650 : div( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
36 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
37 : : const std::array< std::vector< tk::real >, 3 >& coord,
38 : : const std::vector< std::size_t >& triinpoel,
39 : : const tk::Fields& U,
40 : : std::vector< tk::real >& D,
41 : : std::size_t pos )
42 : : // *****************************************************************************
43 : : // Compute divergence of a vector in all points
44 : : //! \param[in] dsupedge Domain superedges
45 : : //! \param[in] dsupint Domain superedge integrals
46 : : //! \param[in] coord Mesh node coordinates
47 : : //! \param[in] triinpoel Boundary face connectivity
48 : : //! \param[in] U Vector whose divergence to compute
49 : : //! \param[in,out] D Vector added to
50 : : //! \param[in] pos Position at which the three vector components are in U
51 : : // *****************************************************************************
52 : : {
53 [ - + ][ - - ]: 650 : Assert( U.nunk() == D.size(), "Size mismatch" );
[ - - ][ - - ]
54 : :
55 : : // Lambda to compute the velocity divergence contribution for edge p-q
56 : 130884 : auto div = [&]( const tk::real d[], std::size_t p, std::size_t q ){
57 : 130884 : return d[0] * (U(p,pos+0) + U(q,pos+0)) +
58 : 130884 : d[1] * (U(p,pos+1) + U(q,pos+1)) +
59 : 130884 : d[2] * (U(p,pos+2) + U(q,pos+2));
60 : 650 : };
61 : :
62 : : #if defined(__clang__)
63 : : #pragma clang diagnostic push
64 : : #pragma clang diagnostic ignored "-Wvla"
65 : : #pragma clang diagnostic ignored "-Wvla-extension"
66 : : #elif defined(STRICT_GNUC)
67 : : #pragma GCC diagnostic push
68 : : #pragma GCC diagnostic ignored "-Wvla"
69 : : #endif
70 : :
71 : : // domain integral
72 : :
73 : : // domain edge contributions: tetrahedron superedges
74 [ + + ]: 12028 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
75 : 11378 : const auto N = dsupedge[0].data() + e*4;
76 : 11378 : const auto d = dsupint[0].data();
77 : : // edge fluxes
78 : : tk::real f[] = {
79 [ + - ]: 11378 : div( d+(e*6+0)*4, N[0], N[1] ),
80 : 22756 : div( d+(e*6+1)*4, N[1], N[2] ),
81 : 22756 : div( d+(e*6+2)*4, N[2], N[0] ),
82 : 22756 : div( d+(e*6+3)*4, N[0], N[3] ),
83 : 22756 : div( d+(e*6+4)*4, N[1], N[3] ),
84 [ + - ][ + - ]: 11378 : div( d+(e*6+5)*4, N[2], N[3] ) };
[ + - ][ + - ]
[ + - ]
85 : : // edge flux contributions
86 : 11378 : D[N[0]] = D[N[0]] - f[0] + f[2] - f[3];
87 : 11378 : D[N[1]] = D[N[1]] + f[0] - f[1] - f[4];
88 : 11378 : D[N[2]] = D[N[2]] + f[1] - f[2] - f[5];
89 : 11378 : D[N[3]] = D[N[3]] + f[3] + f[4] + f[5];
90 : : }
91 : :
92 : : // domain edge contributions: triangle superedges
93 [ + + ]: 12008 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
94 : 11358 : const auto N = dsupedge[1].data() + e*3;
95 : 11358 : const auto d = dsupint[1].data();
96 : : // edge fluxes
97 : : tk::real f[] = {
98 [ + - ]: 11358 : div( d+(e*3+0)*4, N[0], N[1] ),
99 : 22716 : div( d+(e*3+1)*4, N[1], N[2] ),
100 [ + - ][ + - ]: 11358 : div( d+(e*3+2)*4, N[2], N[0] ) };
101 : : // edge flux contributions
102 : 11358 : D[N[0]] = D[N[0]] - f[0] + f[2];
103 : 11358 : D[N[1]] = D[N[1]] + f[0] - f[1];
104 : 11358 : D[N[2]] = D[N[2]] + f[1] - f[2];
105 : : }
106 : :
107 : : // domain edge contributions: edges
108 [ + + ]: 29192 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
109 : 28542 : const auto N = dsupedge[2].data() + e*2;
110 : 28542 : const auto d = dsupint[2].data();
111 : : // edge flux
112 [ + - ]: 28542 : tk::real f = div( d+e*4, N[0], N[1] );
113 : : // edge flux contributions
114 : 28542 : D[N[0]] -= f;
115 : 28542 : D[N[1]] += f;
116 : : }
117 : :
118 : : // boundary integral
119 : :
120 : 650 : const auto& x = coord[0];
121 : 650 : const auto& y = coord[1];
122 : 650 : const auto& z = coord[2];
123 : :
124 [ + + ]: 47642 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
125 : 46992 : const auto N = triinpoel.data() + e*3;
126 : : tk::real n[3];
127 : 46992 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
128 : 46992 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
129 : : n[0], n[1], n[2] );
130 [ + - ]: 46992 : auto uxA = U(N[0],pos+0);
131 [ + - ]: 46992 : auto uyA = U(N[0],pos+1);
132 [ + - ]: 46992 : auto uzA = U(N[0],pos+2);
133 [ + - ]: 46992 : auto uxB = U(N[1],pos+0);
134 [ + - ]: 46992 : auto uyB = U(N[1],pos+1);
135 [ + - ]: 46992 : auto uzB = U(N[1],pos+2);
136 [ + - ]: 46992 : auto uxC = U(N[2],pos+0);
137 [ + - ]: 46992 : auto uyC = U(N[2],pos+1);
138 [ + - ]: 46992 : auto uzC = U(N[2],pos+2);
139 : 46992 : auto ux = (6.0*uxA + uxB + uxC)/8.0;
140 : 46992 : auto uy = (6.0*uyA + uyB + uyC)/8.0;
141 : 46992 : auto uz = (6.0*uzA + uzB + uzC)/8.0;
142 : 46992 : D[N[0]] += ux*n[0] + uy*n[1] + uz*n[2];
143 : 46992 : ux = (uxA + 6.0*uxB + uxC)/8.0;
144 : 46992 : uy = (uyA + 6.0*uyB + uyC)/8.0;
145 : 46992 : uz = (uzA + 6.0*uzB + uzC)/8.0;
146 : 46992 : D[N[1]] += ux*n[0] + uy*n[1] + uz*n[2];
147 : 46992 : ux = (uxA + uxB + 6.0*uxC)/8.0;
148 : 46992 : uy = (uyA + uyB + 6.0*uyC)/8.0;
149 : 46992 : uz = (uzA + uzB + 6.0*uzC)/8.0;
150 : 46992 : D[N[2]] += ux*n[0] + uy*n[1] + uz*n[2];
151 : : }
152 : :
153 : : #if defined(__clang__)
154 : : #pragma clang diagnostic pop
155 : : #elif defined(STRICT_GNUC)
156 : : #pragma GCC diagnostic pop
157 : : #endif
158 : 650 : }
159 : :
160 : : void
161 : 325 : grad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
162 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
163 : : const std::array< std::vector< tk::real >, 3 >& coord,
164 : : const std::vector< std::size_t >& triinpoel,
165 : : const std::vector< tk::real >& U,
166 : : tk::Fields& G )
167 : : // *****************************************************************************
168 : : // Compute gradients of scalar in all points
169 : : //! \param[in] dsupedge Domain superedges
170 : : //! \param[in] dsupint Domain superedge integrals
171 : : //! \param[in] coord Mesh node coordinates
172 : : //! \param[in] triinpoel Boundary face connectivity
173 : : //! \param[in] U Scalar whose gradient to compute
174 : : //! \param[in,out] G Nodal gradient of scalar in all points
175 : : // *****************************************************************************
176 : : {
177 [ - + ]: 325 : if (G.nprop() == 0) return;
178 : :
179 [ - + ][ - - ]: 325 : Assert( G.nunk() == U.size(), "Size mismatch" );
[ - - ][ - - ]
180 [ - + ][ - - ]: 325 : Assert( G.nprop() > 2, "Size mismatch" );
[ - - ][ - - ]
181 [ - + ][ - - ]: 325 : Assert( G.nprop() % 3 == 0, "Size mismatch" );
[ - - ][ - - ]
182 : :
183 : : #if defined(__clang__)
184 : : #pragma clang diagnostic push
185 : : #pragma clang diagnostic ignored "-Wvla"
186 : : #pragma clang diagnostic ignored "-Wvla-extension"
187 : : #elif defined(STRICT_GNUC)
188 : : #pragma GCC diagnostic push
189 : : #pragma GCC diagnostic ignored "-Wvla"
190 : : #endif
191 : :
192 : : // domain integral
193 : :
194 : : // domain edge contributions: tetrahedron superedges
195 [ + + ]: 6014 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
196 : 5689 : const auto N = dsupedge[0].data() + e*4;
197 : 5689 : const auto d = dsupint[0].data();
198 : 5689 : tk::real u[] = { U[N[0]], U[N[1]], U[N[2]], U[N[3]] };
199 [ + + ]: 22756 : for (std::size_t j=0; j<3; ++j) {
200 : : tk::real f[] = {
201 : 17067 : d[(e*6+0)*4+j] * (u[1] + u[0]),
202 : 17067 : d[(e*6+1)*4+j] * (u[2] + u[1]),
203 : 17067 : d[(e*6+2)*4+j] * (u[0] + u[2]),
204 : 17067 : d[(e*6+3)*4+j] * (u[3] + u[0]),
205 : 17067 : d[(e*6+4)*4+j] * (u[3] + u[1]),
206 : 17067 : d[(e*6+5)*4+j] * (u[3] + u[2]) };
207 [ + - ][ + - ]: 17067 : G(N[0],j) = G(N[0],j) - f[0] + f[2] - f[3];
208 [ + - ][ + - ]: 17067 : G(N[1],j) = G(N[1],j) + f[0] - f[1] - f[4];
209 [ + - ][ + - ]: 17067 : G(N[2],j) = G(N[2],j) + f[1] - f[2] - f[5];
210 [ + - ][ + - ]: 17067 : G(N[3],j) = G(N[3],j) + f[3] + f[4] + f[5];
211 : : }
212 : : }
213 : :
214 : : // domain edge contributions: triangle superedges
215 [ + + ]: 6004 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
216 : 5679 : const auto N = dsupedge[1].data() + e*3;
217 : 5679 : const auto d = dsupint[1].data();
218 : 5679 : tk::real u[] = { U[N[0]], U[N[1]], U[N[2]] };
219 [ + + ]: 22716 : for (std::size_t j=0; j<3; ++j) {
220 : : tk::real f[] = {
221 : 17037 : d[(e*3+0)*4+j] * (u[1] + u[0]),
222 : 17037 : d[(e*3+1)*4+j] * (u[2] + u[1]),
223 : 17037 : d[(e*3+2)*4+j] * (u[0] + u[2]) };
224 [ + - ][ + - ]: 17037 : G(N[0],j) = G(N[0],j) - f[0] + f[2];
225 [ + - ][ + - ]: 17037 : G(N[1],j) = G(N[1],j) + f[0] - f[1];
226 [ + - ][ + - ]: 17037 : G(N[2],j) = G(N[2],j) + f[1] - f[2];
227 : : }
228 : : }
229 : :
230 : : // domain edge contributions: edges
231 [ + + ]: 14596 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
232 : 14271 : const auto N = dsupedge[2].data() + e*2;
233 : 14271 : const auto d = dsupint[2].data() + e*4;
234 : 14271 : tk::real u[] = { U[N[0]], U[N[1]] };
235 [ + + ]: 57084 : for (std::size_t j=0; j<3; ++j) {
236 : 42813 : tk::real f = d[j] * (u[1] + u[0]);
237 [ + - ]: 42813 : G(N[0],j) -= f;
238 [ + - ]: 42813 : G(N[1],j) += f;
239 : : }
240 : : }
241 : :
242 : : // boundary integral
243 : :
244 : 325 : const auto& x = coord[0];
245 : 325 : const auto& y = coord[1];
246 : 325 : const auto& z = coord[2];
247 : :
248 [ + + ]: 23821 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
249 : 23496 : const auto N = triinpoel.data() + e*3;
250 : : tk::real n[3];
251 : 23496 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
252 : 23496 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
253 : : n[0], n[1], n[2] );
254 : 23496 : auto uA = U[N[0]];
255 : 23496 : auto uB = U[N[1]];
256 : 23496 : auto uC = U[N[2]];
257 : 23496 : auto f = (6.0*uA + uB + uC)/8.0;
258 [ + - ]: 23496 : G(N[0],0) += f * n[0];
259 [ + - ]: 23496 : G(N[0],1) += f * n[1];
260 [ + - ]: 23496 : G(N[0],2) += f * n[2];
261 : 23496 : f = (uA + 6.0*uB + uC)/8.0;
262 [ + - ]: 23496 : G(N[1],0) += f * n[0];
263 [ + - ]: 23496 : G(N[1],1) += f * n[1];
264 [ + - ]: 23496 : G(N[1],2) += f * n[2];
265 : 23496 : f = (uA + uB + 6.0*uC)/8.0;
266 [ + - ]: 23496 : G(N[2],0) += f * n[0];
267 [ + - ]: 23496 : G(N[2],1) += f * n[1];
268 [ + - ]: 23496 : G(N[2],2) += f * n[2];
269 : : }
270 : :
271 : : #if defined(__clang__)
272 : : #pragma clang diagnostic pop
273 : : #elif defined(STRICT_GNUC)
274 : : #pragma GCC diagnostic pop
275 : : #endif
276 : : }
277 : :
278 : : void
279 : 325 : vgrad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
280 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
281 : : const std::array< std::vector< tk::real >, 3 >& coord,
282 : : const std::vector< std::size_t >& triinpoel,
283 : : const tk::Fields& U,
284 : : tk::Fields& G )
285 : : // *****************************************************************************
286 : : // Compute velocity gradients in all points
287 : : //! \param[in] dsupedge Domain superedges
288 : : //! \param[in] dsupint Domain superedge integrals
289 : : //! \param[in] coord Mesh node coordinates
290 : : //! \param[in] triinpoel Boundary face connectivity
291 : : //! \param[in] U Velocity whose gradient to compute
292 : : //! \param[in,out] G Nodal velocity gradients (9 components) in all points
293 : : // *****************************************************************************
294 : : {
295 [ - + ][ - - ]: 325 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
296 [ - + ][ - - ]: 325 : Assert( G.nprop() == 9, "Size mismatch" );
[ - - ][ - - ]
297 : :
298 : : #if defined(__clang__)
299 : : #pragma clang diagnostic push
300 : : #pragma clang diagnostic ignored "-Wvla"
301 : : #pragma clang diagnostic ignored "-Wvla-extension"
302 : : #elif defined(STRICT_GNUC)
303 : : #pragma GCC diagnostic push
304 : : #pragma GCC diagnostic ignored "-Wvla"
305 : : #endif
306 : :
307 : : // domain integral
308 : :
309 : : // domain edge contributions: tetrahedron superedges
310 [ + + ]: 6014 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
311 : 5689 : const auto N = dsupedge[0].data() + e*4;
312 : 5689 : const auto d = dsupint[0].data();
313 [ + + ]: 22756 : for (std::size_t i=0; i<3; ++i) {
314 [ + - ][ + - ]: 17067 : tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1), U(N[3],i+1) };
[ + - ][ + - ]
315 : 17067 : auto i3 = i*3;
316 [ + + ]: 68268 : for (std::size_t j=0; j<3; ++j) {
317 : 51201 : tk::real f[] = { d[(e*6+0)*4+j] * (u[1] + u[0]),
318 : 51201 : d[(e*6+1)*4+j] * (u[2] + u[1]),
319 : 51201 : d[(e*6+2)*4+j] * (u[0] + u[2]),
320 : 51201 : d[(e*6+3)*4+j] * (u[3] + u[0]),
321 : 51201 : d[(e*6+4)*4+j] * (u[3] + u[1]),
322 : 51201 : d[(e*6+5)*4+j] * (u[3] + u[2]) };
323 [ + - ][ + - ]: 51201 : G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2] - f[3];
324 [ + - ][ + - ]: 51201 : G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1] - f[4];
325 [ + - ][ + - ]: 51201 : G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2] - f[5];
326 [ + - ][ + - ]: 51201 : G(N[3],i3+j) = G(N[3],i3+j) + f[3] + f[4] + f[5];
327 : : }
328 : : }
329 : : }
330 : :
331 : : // domain edge contributions: triangle superedges
332 [ + + ]: 6004 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
333 : 5679 : const auto N = dsupedge[1].data() + e*3;
334 : 5679 : const auto d = dsupint[1].data();
335 [ + + ]: 22716 : for (std::size_t i=0; i<3; ++i) {
336 [ + - ][ + - ]: 17037 : tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1) };
[ + - ]
337 : 17037 : auto i3 = i*3;
338 [ + + ]: 68148 : for (std::size_t j=0; j<3; ++j) {
339 : 51111 : tk::real f[] = { d[(e*3+0)*4+j] * (u[1] + u[0]),
340 : 51111 : d[(e*3+1)*4+j] * (u[2] + u[1]),
341 : 51111 : d[(e*3+2)*4+j] * (u[0] + u[2]) };
342 [ + - ][ + - ]: 51111 : G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2];
343 [ + - ][ + - ]: 51111 : G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1];
344 [ + - ][ + - ]: 51111 : G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2];
345 : : }
346 : : }
347 : : }
348 : :
349 : : // domain edge contributions: edges
350 [ + + ]: 14596 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
351 : 14271 : const auto N = dsupedge[2].data() + e*2;
352 : 14271 : const auto d = dsupint[2].data() + e*4;
353 [ + + ]: 57084 : for (std::size_t i=0; i<3; ++i) {
354 [ + - ][ + - ]: 42813 : tk::real u[] = { U(N[0],i+1), U(N[1],i+1) };
355 : 42813 : auto i3 = i*3;
356 [ + + ]: 171252 : for (std::size_t j=0; j<3; ++j) {
357 : 128439 : tk::real f = d[j] * (u[1] + u[0]);
358 [ + - ]: 128439 : G(N[0],i3+j) -= f;
359 [ + - ]: 128439 : G(N[1],i3+j) += f;
360 : : }
361 : : }
362 : : }
363 : :
364 : : // boundary integral
365 : :
366 : 325 : const auto& x = coord[0];
367 : 325 : const auto& y = coord[1];
368 : 325 : const auto& z = coord[2];
369 : :
370 [ + + ]: 23821 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
371 : 23496 : const auto N = triinpoel.data() + e*3;
372 : : tk::real n[3];
373 : 23496 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
374 : 23496 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
375 : : n[0], n[1], n[2] );
376 [ + + ]: 93984 : for (std::size_t i=0; i<3; ++i) {
377 [ + - ][ + - ]: 70488 : tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1) };
[ + - ]
378 : 70488 : auto i3 = i*3;
379 : 70488 : auto f = (6.0*u[0] + u[1] + u[2])/8.0;
380 [ + - ]: 70488 : G(N[0],i3+0) += f * n[0];
381 [ + - ]: 70488 : G(N[0],i3+1) += f * n[1];
382 [ + - ]: 70488 : G(N[0],i3+2) += f * n[2];
383 : 70488 : f = (u[0] + 6.0*u[1] + u[2])/8.0;
384 [ + - ]: 70488 : G(N[1],i3+0) += f * n[0];
385 [ + - ]: 70488 : G(N[1],i3+1) += f * n[1];
386 [ + - ]: 70488 : G(N[1],i3+2) += f * n[2];
387 : 70488 : f = (u[0] + u[1] + 6.0*u[2])/8.0;
388 [ + - ]: 70488 : G(N[2],i3+0) += f * n[0];
389 [ + - ]: 70488 : G(N[2],i3+1) += f * n[1];
390 [ + - ]: 70488 : G(N[2],i3+2) += f * n[2];
391 : : }
392 : : }
393 : :
394 : : #if defined(__clang__)
395 : : #pragma clang diagnostic pop
396 : : #elif defined(STRICT_GNUC)
397 : : #pragma GCC diagnostic pop
398 : : #endif
399 : 325 : }
400 : :
401 : : static tk::real
402 : 588978 : flux( const tk::Fields& U,
403 : : const tk::Fields& G,
404 : : std::size_t i,
405 : : std::size_t j,
406 : : std::size_t p,
407 : : std::size_t q )
408 : : // *****************************************************************************
409 : : //! Compute momentum flux over edge of points p-q
410 : : //! \param[in] U Velocity vector
411 : : //! \param[in] G Velocity gradients
412 : : //! \param[in] i Tensor component, 1st index
413 : : //! \param[in] j Tensor component, 2nd index
414 : : //! \param[in] p Left node index of edge
415 : : //! \param[in] q Right node index of edge
416 : : //! \return Momentum flux contribution for edge p-q
417 : : // *****************************************************************************
418 : : {
419 : 588978 : auto inv = U(p,i)*U(p,j) + U(q,i)*U(q,j);
420 : :
421 : 588978 : auto eps = std::numeric_limits< tk::real >::epsilon();
422 : 588978 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
423 [ - + ]: 588978 : if (mu < eps) return -inv;
424 : :
425 : 588978 : auto vis = G(p,i*3+j) + G(p,j*3+i) + G(q,i*3+j) + G(q,j*3+i);
426 [ + + ]: 588978 : if (i == j) {
427 : 196326 : vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) + G(q,0) + G(q,4) + G(q,8) );
428 : : }
429 : 588978 : return mu*vis - inv;
430 : : }
431 : :
432 : : static tk::real
433 : 634392 : flux( const tk::Fields& U,
434 : : const tk::Fields& G,
435 : : std::size_t i,
436 : : std::size_t j,
437 : : std::size_t p )
438 : : // *****************************************************************************
439 : : //! Compute momentum flux in point p
440 : : //! \param[in] U Velocity vector
441 : : //! \param[in] G Velocity gradients
442 : : //! \param[in] i Tensor component, 1st index
443 : : //! \param[in] j Tensor component, 2nd index
444 : : //! \param[in] p Node index of point
445 : : //! \return Momentum flux contribution for point p
446 : : // *****************************************************************************
447 : : {
448 : 634392 : auto inv = U(p,i)*U(p,j);
449 : :
450 : 634392 : auto eps = std::numeric_limits< tk::real >::epsilon();
451 : 634392 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
452 [ - + ]: 634392 : if (mu < eps) return -inv;
453 : :
454 : 634392 : auto vis = G(p,i*3+j) + G(p,j*3+i);
455 [ + + ]: 634392 : if (i == j) {
456 : 211464 : vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) );
457 : : }
458 : 634392 : return mu*vis - inv;
459 : : }
460 : :
461 : : void
462 : 325 : flux( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
463 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
464 : : const std::array< std::vector< tk::real >, 3 >& coord,
465 : : const std::vector< std::size_t >& triinpoel,
466 : : const tk::Fields& U,
467 : : const tk::Fields& G,
468 : : tk::Fields& F )
469 : : // *****************************************************************************
470 : : // Compute momentum flux in all points
471 : : //! \param[in] dsupedge Domain superedges
472 : : //! \param[in] dsupint Domain superedge integrals
473 : : //! \param[in] coord Mesh node coordinates
474 : : //! \param[in] triinpoel Boundary face connectivity
475 : : //! \param[in] U Velocity field
476 : : //! \param[in] G Velocity gradients, dui/dxj, 9 components
477 : : //! \param[in,out] F Momentum flux, Fi = d[ sij - uiuj ] / dxj, where
478 : : //! s_ij = mu[ dui/dxj + duj/dxi - 2/3 duk/dxk delta_ij ]
479 : : // *****************************************************************************
480 : : {
481 [ - + ][ - - ]: 325 : Assert( F.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
482 [ - + ][ - - ]: 325 : Assert( F.nprop() == 3, "Size mismatch" );
[ - - ][ - - ]
483 [ - + ][ - - ]: 325 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
484 [ - + ][ - - ]: 325 : Assert( G.nprop() == 9, "Size mismatch" );
[ - - ][ - - ]
485 : :
486 : : #if defined(__clang__)
487 : : #pragma clang diagnostic push
488 : : #pragma clang diagnostic ignored "-Wvla"
489 : : #pragma clang diagnostic ignored "-Wvla-extension"
490 : : #elif defined(STRICT_GNUC)
491 : : #pragma GCC diagnostic push
492 : : #pragma GCC diagnostic ignored "-Wvla"
493 : : #endif
494 : :
495 : : // domain integral
496 : :
497 : : // domain edge contributions: tetrahedron superedges
498 [ + + ]: 6014 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
499 : 5689 : const auto N = dsupedge[0].data() + e*4;
500 : 5689 : const auto d = dsupint[0].data();
501 [ + + ]: 22756 : for (std::size_t i=0; i<3; ++i) {
502 [ + + ]: 68268 : for (std::size_t j=0; j<3; ++j) {
503 [ + - ]: 51201 : tk::real f[] = { d[(e*6+0)*4+j] * flux(U,G,i,j,N[1],N[0]),
504 : 102402 : d[(e*6+1)*4+j] * flux(U,G,i,j,N[2],N[1]),
505 : 102402 : d[(e*6+2)*4+j] * flux(U,G,i,j,N[0],N[2]),
506 : 102402 : d[(e*6+3)*4+j] * flux(U,G,i,j,N[3],N[0]),
507 : 102402 : d[(e*6+4)*4+j] * flux(U,G,i,j,N[3],N[1]),
508 [ + - ][ + - ]: 51201 : d[(e*6+5)*4+j] * flux(U,G,i,j,N[3],N[2]) };
[ + - ][ + - ]
[ + - ]
509 [ + - ][ + - ]: 51201 : F(N[0],i) = F(N[0],i) - f[0] + f[2] - f[3];
510 [ + - ][ + - ]: 51201 : F(N[1],i) = F(N[1],i) + f[0] - f[1] - f[4];
511 [ + - ][ + - ]: 51201 : F(N[2],i) = F(N[2],i) + f[1] - f[2] - f[5];
512 [ + - ][ + - ]: 51201 : F(N[3],i) = F(N[3],i) + f[3] + f[4] + f[5];
513 : : }
514 : : }
515 : : }
516 : :
517 : : // domain edge contributions: triangle superedges
518 [ + + ]: 6004 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
519 : 5679 : const auto N = dsupedge[1].data() + e*3;
520 : 5679 : const auto d = dsupint[1].data();
521 [ + + ]: 22716 : for (std::size_t i=0; i<3; ++i) {
522 [ + + ]: 68148 : for (std::size_t j=0; j<3; ++j) {
523 [ + - ]: 51111 : tk::real f[] = { d[(e*3+0)*4+j] * flux(U,G,i,j,N[1],N[0]),
524 : 102222 : d[(e*3+1)*4+j] * flux(U,G,i,j,N[2],N[1]),
525 [ + - ][ + - ]: 51111 : d[(e*3+2)*4+j] * flux(U,G,i,j,N[0],N[2]), };
526 [ + - ][ + - ]: 51111 : F(N[0],i) = F(N[0],i) - f[0] + f[2];
527 [ + - ][ + - ]: 51111 : F(N[1],i) = F(N[1],i) + f[0] - f[1];
528 [ + - ][ + - ]: 51111 : F(N[2],i) = F(N[2],i) + f[1] - f[2];
529 : : }
530 : : }
531 : : }
532 : :
533 : : // domain edge contributions: edges
534 [ + + ]: 14596 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
535 : 14271 : const auto N = dsupedge[2].data() + e*2;
536 : 14271 : const auto d = dsupint[2].data() + e*4;
537 [ + + ]: 57084 : for (std::size_t i=0; i<3; ++i) {
538 [ + + ]: 171252 : for (std::size_t j=0; j<3; ++j) {
539 : 128439 : tk::real f = d[j] * flux(U,G,i,j,N[1],N[0]);
540 : 128439 : F(N[0],i) -= f;
541 : 128439 : F(N[1],i) += f;
542 : : }
543 : : }
544 : : }
545 : :
546 : : // boundary integral
547 : :
548 : 325 : const auto& x = coord[0];
549 : 325 : const auto& y = coord[1];
550 : 325 : const auto& z = coord[2];
551 : :
552 [ + + ]: 23821 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
553 : 23496 : const auto N = triinpoel.data() + e*3;
554 : : tk::real n[3];
555 : 23496 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
556 : 23496 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
557 : : n[0], n[1], n[2] );
558 [ + + ]: 93984 : for (std::size_t i=0; i<3; ++i) {
559 [ + - ]: 70488 : auto fxA = flux(U,G,i,0,N[0]);
560 [ + - ]: 70488 : auto fyA = flux(U,G,i,1,N[0]);
561 [ + - ]: 70488 : auto fzA = flux(U,G,i,2,N[0]);
562 [ + - ]: 70488 : auto fxB = flux(U,G,i,0,N[1]);
563 [ + - ]: 70488 : auto fyB = flux(U,G,i,1,N[1]);
564 [ + - ]: 70488 : auto fzB = flux(U,G,i,2,N[1]);
565 [ + - ]: 70488 : auto fxC = flux(U,G,i,0,N[2]);
566 [ + - ]: 70488 : auto fyC = flux(U,G,i,1,N[2]);
567 [ + - ]: 70488 : auto fzC = flux(U,G,i,2,N[2]);
568 : 70488 : auto fx = (6.0*fxA + fxB + fxC)/8.0;
569 : 70488 : auto fy = (6.0*fyA + fyB + fyC)/8.0;
570 : 70488 : auto fz = (6.0*fzA + fzB + fzC)/8.0;
571 [ + - ]: 70488 : F(N[0],i) += fx*n[0] + fy*n[1] + fz*n[2];
572 : 70488 : fx = (fxA + 6.0*fxB + fxC)/8.0;
573 : 70488 : fy = (fyA + 6.0*fyB + fyC)/8.0;
574 : 70488 : fz = (fzA + 6.0*fzB + fzC)/8.0;
575 [ + - ]: 70488 : F(N[1],i) += fx*n[0] + fy*n[1] + fz*n[2];
576 : 70488 : fx = (fxA + fxB + 6.0*fxC)/8.0;
577 : 70488 : fy = (fyA + fyB + 6.0*fyC)/8.0;
578 : 70488 : fz = (fzA + fzB + 6.0*fzC)/8.0;
579 [ + - ]: 70488 : F(N[2],i) += fx*n[0] + fy*n[1] + fz*n[2];
580 : : }
581 : : }
582 : :
583 : : #if defined(__clang__)
584 : : #pragma clang diagnostic pop
585 : : #elif defined(STRICT_GNUC)
586 : : #pragma GCC diagnostic pop
587 : : #endif
588 : 325 : }
589 : :
590 : : void
591 : 140 : grad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
592 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
593 : : const std::array< std::vector< tk::real >, 3 >& coord,
594 : : const std::vector< std::size_t >& triinpoel,
595 : : const tk::Fields& U,
596 : : tk::Fields& G )
597 : : // *****************************************************************************
598 : : // Compute gradients of scalar in all points
599 : : //! \param[in] dsupedge Domain superedges
600 : : //! \param[in] dsupint Domain superedge integrals
601 : : //! \param[in] coord Mesh node coordinates
602 : : //! \param[in] triinpoel Boundary face connectivity
603 : : //! \param[in] U Unknown vector whose gradient to compute
604 : : //! \param[in,out] G Nodal gradients in all points
605 : : // *****************************************************************************
606 : : {
607 [ - + ]: 140 : if (G.nprop() == 0) return;
608 : :
609 [ - + ][ - - ]: 140 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
610 [ - + ][ - - ]: 140 : Assert( G.nprop() > 2, "Size mismatch" );
[ - - ][ - - ]
611 [ - + ][ - - ]: 140 : Assert( G.nprop() % 3 == 0, "Size mismatch" );
[ - - ][ - - ]
612 [ - + ][ - - ]: 140 : Assert( G.nprop() == (U.nprop()-1)*3, "Size mismatch" );
[ - - ][ - - ]
613 : :
614 : 140 : const auto ncomp = U.nprop();
615 : :
616 : : #if defined(__clang__)
617 : : #pragma clang diagnostic push
618 : : #pragma clang diagnostic ignored "-Wvla"
619 : : #pragma clang diagnostic ignored "-Wvla-extension"
620 : : #elif defined(STRICT_GNUC)
621 : : #pragma GCC diagnostic push
622 : : #pragma GCC diagnostic ignored "-Wvla"
623 : : #endif
624 : :
625 : : // domain integral
626 : :
627 : : // domain edge contributions: tetrahedron superedges
628 [ + + ]: 70000 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
629 : 69860 : const auto N = dsupedge[0].data() + e*4;
630 : 69860 : const auto d = dsupint[0].data();
631 [ + + ]: 279440 : for (std::size_t c=1; c<ncomp; ++c) {
632 [ + - ][ + - ]: 209580 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c), U(N[3],c) };
[ + - ][ + - ]
633 : 209580 : auto g = (c-1)*3;
634 [ + + ]: 838320 : for (std::size_t j=0; j<3; ++j) {
635 : : tk::real f[] = {
636 : 628740 : d[(e*6+0)*4+j] * (u[1] + u[0]),
637 : 628740 : d[(e*6+1)*4+j] * (u[2] + u[1]),
638 : 628740 : d[(e*6+2)*4+j] * (u[0] + u[2]),
639 : 628740 : d[(e*6+3)*4+j] * (u[3] + u[0]),
640 : 628740 : d[(e*6+4)*4+j] * (u[3] + u[1]),
641 : 628740 : d[(e*6+5)*4+j] * (u[3] + u[2]) };
642 [ + - ][ + - ]: 628740 : G(N[0],g+j) = G(N[0],g+j) - f[0] + f[2] - f[3];
643 [ + - ][ + - ]: 628740 : G(N[1],g+j) = G(N[1],g+j) + f[0] - f[1] - f[4];
644 [ + - ][ + - ]: 628740 : G(N[2],g+j) = G(N[2],g+j) + f[1] - f[2] - f[5];
645 [ + - ][ + - ]: 628740 : G(N[3],g+j) = G(N[3],g+j) + f[3] + f[4] + f[5];
646 : : }
647 : : }
648 : : }
649 : :
650 : : // domain edge contributions: triangle superedges
651 [ + + ]: 69440 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
652 : 69300 : const auto N = dsupedge[1].data() + e*3;
653 : 69300 : const auto d = dsupint[1].data();
654 [ + + ]: 277200 : for (std::size_t c=1; c<ncomp; ++c) {
655 [ + - ][ + - ]: 207900 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c) };
[ + - ]
656 : 207900 : auto g = (c-1)*3;
657 [ + + ]: 831600 : for (std::size_t j=0; j<3; ++j) {
658 : : tk::real f[] = {
659 : 623700 : d[(e*3+0)*4+j] * (u[1] + u[0]),
660 : 623700 : d[(e*3+1)*4+j] * (u[2] + u[1]),
661 : 623700 : d[(e*3+2)*4+j] * (u[0] + u[2]) };
662 [ + - ][ + - ]: 623700 : G(N[0],g+j) = G(N[0],g+j) - f[0] + f[2];
663 [ + - ][ + - ]: 623700 : G(N[1],g+j) = G(N[1],g+j) + f[0] - f[1];
664 [ + - ][ + - ]: 623700 : G(N[2],g+j) = G(N[2],g+j) + f[1] - f[2];
665 : : }
666 : : }
667 : : }
668 : :
669 : : // domain edge contributions: edges
670 [ + + ]: 132020 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
671 : 131880 : const auto N = dsupedge[2].data() + e*2;
672 : 131880 : const auto d = dsupint[2].data() + e*4;
673 [ + + ]: 527520 : for (std::size_t c=1; c<ncomp; ++c) {
674 [ + - ][ + - ]: 395640 : tk::real u[] = { U(N[0],c), U(N[1],c) };
675 : 395640 : auto g = (c-1)*3;
676 [ + + ]: 1582560 : for (std::size_t j=0; j<3; ++j) {
677 : 1186920 : tk::real f = d[j] * (u[1] + u[0]);
678 [ + - ]: 1186920 : G(N[0],g+j) -= f;
679 [ + - ]: 1186920 : G(N[1],g+j) += f;
680 : : }
681 : : }
682 : : }
683 : :
684 : : // boundary integral
685 : :
686 : 140 : const auto& x = coord[0];
687 : 140 : const auto& y = coord[1];
688 : 140 : const auto& z = coord[2];
689 : :
690 [ + + ]: 338940 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
691 : 338800 : const auto N = triinpoel.data() + e*3;
692 : : tk::real n[3];
693 : 338800 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
694 : 338800 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
695 : : n[0], n[1], n[2] );
696 [ + + ]: 1355200 : for (std::size_t c=1; c<ncomp; ++c) {
697 : 1016400 : auto g = (c-1)*3;
698 [ + - ]: 1016400 : auto uA = U(N[0],c);
699 [ + - ]: 1016400 : auto uB = U(N[1],c);
700 [ + - ]: 1016400 : auto uC = U(N[2],c);
701 : 1016400 : auto f = (6.0*uA + uB + uC)/8.0;
702 [ + - ]: 1016400 : G(N[0],g+0) += f * n[0];
703 [ + - ]: 1016400 : G(N[0],g+1) += f * n[1];
704 [ + - ]: 1016400 : G(N[0],g+2) += f * n[2];
705 : 1016400 : f = (uA + 6.0*uB + uC)/8.0;
706 [ + - ]: 1016400 : G(N[1],g+0) += f * n[0];
707 [ + - ]: 1016400 : G(N[1],g+1) += f * n[1];
708 [ + - ]: 1016400 : G(N[1],g+2) += f * n[2];
709 : 1016400 : f = (uA + uB + 6.0*uC)/8.0;
710 [ + - ]: 1016400 : G(N[2],g+0) += f * n[0];
711 [ + - ]: 1016400 : G(N[2],g+1) += f * n[1];
712 [ + - ]: 1016400 : G(N[2],g+2) += f * n[2];
713 : : }
714 : : }
715 : :
716 : : #if defined(__clang__)
717 : : #pragma clang diagnostic pop
718 : : #elif defined(STRICT_GNUC)
719 : : #pragma GCC diagnostic pop
720 : : #endif
721 : : }
722 : :
723 : : static void
724 : 1748200 : adv_damp2( const tk::real supint[],
725 : : const tk::Fields& U,
726 : : const tk::Fields&,
727 : : const std::array< std::vector< tk::real >, 3 >&,
728 : : std::size_t p,
729 : : std::size_t q,
730 : : tk::real f[] )
731 : : // *****************************************************************************
732 : : //! Compute advection fluxes on a single edge using 2nd-order damping
733 : : //! \param[in] supint Edge integral
734 : : //! \param[in] U Velocity and transported scalars at recent time step
735 : : //! \param[in] p Left node index of edge
736 : : //! \param[in] q Right node index of edge
737 : : //! \param[in,out] f Flux computed
738 : : // *****************************************************************************
739 : : {
740 : 1748200 : auto nx = supint[0];
741 : 1748200 : auto ny = supint[1];
742 : 1748200 : auto nz = supint[2];
743 : :
744 : : // left state
745 [ + - ]: 1748200 : auto pL = U(p,0);
746 [ + - ]: 1748200 : auto uL = U(p,1);
747 [ + - ]: 1748200 : auto vL = U(p,2);
748 [ + - ]: 1748200 : auto wL = U(p,3);
749 : 1748200 : auto vnL = uL*nx + vL*ny + wL*nz;
750 : :
751 : : // right state
752 [ + - ]: 1748200 : auto pR = U(q,0);
753 [ + - ]: 1748200 : auto uR = U(q,1);
754 [ + - ]: 1748200 : auto vR = U(q,2);
755 [ + - ]: 1748200 : auto wR = U(q,3);
756 : 1748200 : auto vnR = uR*nx + vR*ny + wR*nz;
757 : :
758 : 1748200 : auto s = g_cfg.get< tag::soundspeed >();
759 : 1748200 : auto s2 = s*s;
760 : 1748200 : auto d = supint[3] * g_cfg.get< tag::mat_dyn_viscosity >();
761 : :
762 : : // flow
763 : 1748200 : auto pf = pL + pR;
764 : 1748200 : f[0] = (vnL + vnR)*s2;
765 : 1748200 : f[1] = uL*vnL + uR*vnR + pf*nx - d*(uR - uL);
766 : 1748200 : f[2] = vL*vnL + vR*vnR + pf*ny - d*(vR - vL);
767 : 1748200 : f[3] = wL*vnL + wR*vnR + pf*nz - d*(wR - wL);
768 : :
769 : : // artificial viscosity
770 : 1748200 : const auto stab2 = g_cfg.get< tag::stab2 >();
771 [ + + ]: 1748200 : if (!stab2) return;
772 : :
773 : 1472460 : auto len = tk::length( nx, ny, nz );
774 : 1472460 : auto sl = std::abs(vnL) + s*len;
775 : 1472460 : auto sr = std::abs(vnR) + s*len;
776 : 1472460 : auto aw = g_cfg.get< tag::stab2coef >() * std::max(sl,sr) * len;
777 : 1472460 : f[0] += aw * (pR - pL)*s2;
778 : 1472460 : f[1] += aw * (uR - uL);
779 : 1472460 : f[2] += aw * (vR - vL);
780 : 1472460 : f[3] += aw * (wR - wL);
781 : : }
782 : :
783 : : static void
784 : 758940 : adv_damp4( const tk::real supint[],
785 : : const tk::Fields& U,
786 : : const tk::Fields& G,
787 : : const std::array< std::vector< tk::real >, 3 >& coord,
788 : : std::size_t p,
789 : : std::size_t q,
790 : : tk::real f[] )
791 : : // *****************************************************************************
792 : : //! Compute advection fluxes on a single edge using 4th-order damping
793 : : //! \param[in] supint Edge integral
794 : : //! \param[in] U Velocity and transported scalars at recent time step
795 : : //! \param[in] G Gradients of velocity and transported scalars
796 : : //! \param[in] coord Mesh node coordinates
797 : : //! \param[in] p Left node index of edge
798 : : //! \param[in] q Right node index of edge
799 : : //! \param[in,out] f Flux computed
800 : : // *****************************************************************************
801 : : {
802 : 758940 : const auto ncomp = U.nprop();
803 : :
804 : 758940 : const auto& x = coord[0];
805 : 758940 : const auto& y = coord[1];
806 : 758940 : const auto& z = coord[2];
807 : :
808 : : // edge vector
809 : 758940 : auto dx = x[p] - x[q];
810 : 758940 : auto dy = y[p] - y[q];
811 : 758940 : auto dz = z[p] - z[q];
812 : :
813 : : #if defined(__clang__)
814 : : #pragma clang diagnostic push
815 : : #pragma clang diagnostic ignored "-Wvla"
816 : : #pragma clang diagnostic ignored "-Wvla-extension"
817 : : #elif defined(STRICT_GNUC)
818 : : #pragma GCC diagnostic push
819 : : #pragma GCC diagnostic ignored "-Wvla"
820 : : #endif
821 : :
822 [ + - ][ + - ]: 758940 : tk::real uL[] = { U(p,1), U(p,2), U(p,3) };
[ + - ]
823 [ + - ][ + - ]: 758940 : tk::real uR[] = { U(q,1), U(q,2), U(q,3) };
[ + - ]
824 : :
825 : : // MUSCL reconstruction in edge-end points
826 [ + + ]: 3035760 : for (std::size_t c=0; c<ncomp-1; ++c) {
827 : 2276820 : auto g = c*3;
828 [ + - ][ + - ]: 2276820 : auto g1 = G(p,g+0)*dx + G(p,g+1)*dy + G(p,g+2)*dz;
[ + - ]
829 [ + - ][ + - ]: 2276820 : auto g2 = G(q,g+0)*dx + G(q,g+1)*dy + G(q,g+2)*dz;
[ + - ]
830 : 2276820 : auto delta2 = uR[c] - uL[c];
831 : 2276820 : auto delta1 = 2.0 * g1 - delta2;
832 : 2276820 : auto delta3 = 2.0 * g2 - delta2;
833 : :
834 : : // van Leer limiter
835 : 2276820 : auto rL = (delta2 + muscl_eps) / (delta1 + muscl_eps);
836 : 2276820 : auto rR = (delta2 + muscl_eps) / (delta3 + muscl_eps);
837 : 2276820 : auto rLinv = (delta1 + muscl_eps) / (delta2 + muscl_eps);
838 : 2276820 : auto rRinv = (delta3 + muscl_eps) / (delta2 + muscl_eps);
839 : 2276820 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
840 : 2276820 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
841 : 2276820 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
842 : 2276820 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
843 : : // update unknowns with reconstructed unknowns
844 : 2276820 : uL[c] += 0.25*(delta1*(1.0-muscl_const)*phiL +
845 : 2276820 : delta2*(1.0+muscl_const)*phi_L_inv);
846 : 2276820 : uR[c] -= 0.25*(delta3*(1.0-muscl_const)*phiR +
847 : 2276820 : delta2*(1.0+muscl_const)*phi_R_inv);
848 : : }
849 : :
850 : 758940 : auto nx = supint[0];
851 : 758940 : auto ny = supint[1];
852 : 758940 : auto nz = supint[2];
853 : :
854 : : // normal velocities
855 : 758940 : auto vnL = uL[0]*nx + uL[1]*ny + uL[2]*nz;
856 : 758940 : auto vnR = uR[0]*nx + uR[1]*ny + uR[2]*nz;
857 : :
858 : 758940 : auto s = g_cfg.get< tag::soundspeed >();
859 : 758940 : auto s2 = s*s;
860 : 758940 : auto d = supint[3] * g_cfg.get< tag::mat_dyn_viscosity >();
861 : :
862 : : // flow
863 [ + - ][ + - ]: 758940 : auto pf = U(p,0) + U(q,0);
864 : 758940 : f[0] = (vnL + vnR)*s2;
865 : 758940 : f[1] = uL[0]*vnL + uR[0]*vnR + pf*nx - d*(uR[0] - uL[0]);
866 : 758940 : f[2] = uL[1]*vnL + uR[1]*vnR + pf*ny - d*(uR[1] - uL[1]);
867 : 758940 : f[3] = uL[2]*vnL + uR[2]*vnR + pf*nz - d*(uR[2] - uL[2]);
868 : :
869 : : // artificial viscosity
870 : 758940 : const auto stab2 = g_cfg.get< tag::stab2 >();
871 [ - + ]: 758940 : if (!stab2) return;
872 : :
873 : 758940 : auto len = tk::length( nx, ny, nz );
874 : 758940 : auto sl = std::abs(vnL) + s*len;
875 : 758940 : auto sr = std::abs(vnR) + s*len;
876 : 758940 : auto aw = g_cfg.get< tag::stab2coef >() * std::max(sl,sr) * len;
877 [ + - ][ + - ]: 758940 : f[0] += aw * (U(q,0) - U(p,0))*s2;
878 : 758940 : f[1] += aw * (uR[0] - uL[0]);
879 : 758940 : f[2] += aw * (uR[1] - uL[1]);
880 : 758940 : f[3] += aw * (uR[2] - uL[2]);
881 : :
882 : : #if defined(__clang__)
883 : : #pragma clang diagnostic pop
884 : : #elif defined(STRICT_GNUC)
885 : : #pragma GCC diagnostic pop
886 : : #endif
887 : : }
888 : :
889 : : static void
890 : 4660 : adv( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
891 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
892 : : const std::array< std::vector< tk::real >, 3 >& coord,
893 : : const std::vector< std::size_t >& triinpoel,
894 : : const tk::Fields& U,
895 : : const tk::Fields& G,
896 : : // cppcheck-suppress constParameter
897 : : tk::Fields& R )
898 : : // *****************************************************************************
899 : : //! Add advection to rhs
900 : : //! \param[in] dsupedge Domain superedges
901 : : //! \param[in] dsupint Domain superedge integrals
902 : : //! \param[in] coord Mesh node coordinates
903 : : //! \param[in] triinpoel Boundary face connectivity
904 : : //! \param[in] U Velocity and transported scalars at recent time step
905 : : //! \param[in] G Gradients of velocity and transported scalars
906 : : //! \param[in,out] R Right-hand side vector added to
907 : : // *****************************************************************************
908 : : {
909 : 4660 : const auto ncomp = U.nprop();
910 : :
911 : : // configure advection
912 : 0 : auto adv = [](){
913 : 4660 : const auto& flux = g_cfg.get< tag::flux >();
914 [ + + ]: 4660 : if (flux == "damp2") return adv_damp2;
915 [ + - ]: 140 : else if (flux == "damp4") return adv_damp4;
916 [ - - ][ - - ]: 0 : else Throw( "Flux not correctly configured" );
[ - - ]
917 [ + - ]: 4660 : }();
918 : :
919 : : #if defined(__clang__)
920 : : #pragma clang diagnostic push
921 : : #pragma clang diagnostic ignored "-Wvla"
922 : : #pragma clang diagnostic ignored "-Wvla-extension"
923 : : #elif defined(STRICT_GNUC)
924 : : #pragma GCC diagnostic push
925 : : #pragma GCC diagnostic ignored "-Wvla"
926 : : #endif
927 : :
928 : : // domain integral
929 : :
930 : : // domain edge contributions: tetrahedron superedges
931 [ + + ]: 229140 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
932 : 224480 : const auto N = dsupedge[0].data() + e*4;
933 : 224480 : tk::real f[6][ncomp];
934 : 224480 : const auto d = dsupint[0].data();
935 [ + - ]: 224480 : adv( d+(e*6+0)*4, U, G, coord, N[0], N[1], f[0] );
936 [ + - ]: 224480 : adv( d+(e*6+1)*4, U, G, coord, N[1], N[2], f[1] );
937 [ + - ]: 224480 : adv( d+(e*6+2)*4, U, G, coord, N[2], N[0], f[2] );
938 [ + - ]: 224480 : adv( d+(e*6+3)*4, U, G, coord, N[0], N[3], f[3] );
939 [ + - ]: 224480 : adv( d+(e*6+4)*4, U, G, coord, N[1], N[3], f[4] );
940 [ + - ]: 224480 : adv( d+(e*6+5)*4, U, G, coord, N[2], N[3], f[5] );
941 [ + + ]: 1122400 : for (std::size_t c=0; c<ncomp; ++c) {
942 [ + - ][ + - ]: 897920 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
943 [ + - ][ + - ]: 897920 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
944 [ + - ][ + - ]: 897920 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
945 [ + - ][ + - ]: 897920 : R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
946 : : }
947 : 224480 : }
948 : :
949 : : // domain edge contributions: triangle superedges
950 [ + + ]: 230050 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
951 : 225390 : const auto N = dsupedge[1].data() + e*3;
952 : 225390 : tk::real f[3][ncomp];
953 : 225390 : const auto d = dsupint[1].data();
954 [ + - ]: 225390 : adv( d+(e*3+0)*4, U, G, coord, N[0], N[1], f[0] );
955 [ + - ]: 225390 : adv( d+(e*3+1)*4, U, G, coord, N[1], N[2], f[1] );
956 [ + - ]: 225390 : adv( d+(e*3+2)*4, U, G, coord, N[2], N[0], f[2] );
957 [ + + ]: 1126950 : for (std::size_t c=0; c<ncomp; ++c) {
958 [ + - ][ + - ]: 901560 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
959 [ + - ][ + - ]: 901560 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
960 [ + - ][ + - ]: 901560 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
961 : : }
962 : 225390 : }
963 : :
964 : : // domain edge contributions: edges
965 [ + + ]: 488750 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
966 : 484090 : const auto N = dsupedge[2].data() + e*2;
967 : 484090 : tk::real f[ncomp];
968 : 484090 : const auto d = dsupint[2].data();
969 [ + - ]: 484090 : adv( d+e*4, U, G, coord, N[0], N[1], f );
970 [ + + ]: 2420450 : for (std::size_t c=0; c<ncomp; ++c) {
971 [ + - ]: 1936360 : R(N[0],c) -= f[c];
972 [ + - ]: 1936360 : R(N[1],c) += f[c];
973 : : }
974 : 484090 : }
975 : :
976 : : // boundary integral
977 : :
978 : 4660 : const auto& x = coord[0];
979 : 4660 : const auto& y = coord[1];
980 : 4660 : const auto& z = coord[2];
981 : :
982 : 4660 : auto s = g_cfg.get< tag::soundspeed >();
983 : 4660 : auto s2 = s * s;
984 : :
985 [ + + ]: 1046660 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
986 : 2084000 : const auto N = triinpoel.data() + e*3;
987 : : tk::real n[3];
988 : 1042000 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
989 : 1042000 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
990 : : n[0], n[1], n[2] );
991 : 1042000 : tk::real f[ncomp][3];
992 : :
993 [ + - ]: 1042000 : auto p = U(N[0],0);
994 [ + - ]: 1042000 : auto u = U(N[0],1);
995 [ + - ]: 1042000 : auto v = U(N[0],2);
996 [ + - ]: 1042000 : auto w = U(N[0],3);
997 : 1042000 : auto vn = n[0]*u + n[1]*v + n[2]*w;
998 : 1042000 : f[0][0] = vn * s2;
999 : 1042000 : f[1][0] = u*vn + p*n[0];
1000 : 1042000 : f[2][0] = v*vn + p*n[1];
1001 : 1042000 : f[3][0] = w*vn + p*n[2];
1002 : :
1003 [ + - ]: 1042000 : p = U(N[1],0);
1004 [ + - ]: 1042000 : u = U(N[1],1);
1005 [ + - ]: 1042000 : v = U(N[1],2);
1006 [ + - ]: 1042000 : w = U(N[1],3);
1007 : 1042000 : vn = n[0]*u + n[1]*v + n[2]*w;
1008 : 1042000 : f[0][1] = vn * s2;
1009 : 1042000 : f[1][1] = u*vn + p*n[0];
1010 : 1042000 : f[2][1] = v*vn + p*n[1];
1011 : 1042000 : f[3][1] = w*vn + p*n[2];
1012 : :
1013 [ + - ]: 1042000 : p = U(N[2],0);
1014 [ + - ]: 1042000 : u = U(N[2],1);
1015 [ + - ]: 1042000 : v = U(N[2],2);
1016 [ + - ]: 1042000 : w = U(N[2],3);
1017 : 1042000 : vn = n[0]*u + n[1]*v + n[2]*w;
1018 : 1042000 : f[0][2] = vn * s2;
1019 : 1042000 : f[1][2] = u*vn + p*n[0];
1020 : 1042000 : f[2][2] = v*vn + p*n[1];
1021 : 1042000 : f[3][2] = w*vn + p*n[2];
1022 : :
1023 [ + + ]: 5210000 : for (std::size_t c=0; c<ncomp; ++c) {
1024 [ + - ]: 4168000 : R(N[0],c) += (6.0*f[c][0] + f[c][1] + f[c][2])/8.0;
1025 [ + - ]: 4168000 : R(N[1],c) += (f[c][0] + 6.0*f[c][1] + f[c][2])/8.0;
1026 [ + - ]: 4168000 : R(N[2],c) += (f[c][0] + f[c][1] + 6.0*f[c][2])/8.0;
1027 : : }
1028 : 1042000 : }
1029 : :
1030 : : #if defined(__clang__)
1031 : : #pragma clang diagnostic pop
1032 : : #elif defined(STRICT_GNUC)
1033 : : #pragma GCC diagnostic pop
1034 : : #endif
1035 : 4660 : }
1036 : :
1037 : : void
1038 : 4660 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
1039 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
1040 : : const std::array< std::vector< tk::real >, 3 >& coord,
1041 : : const std::vector< std::size_t >& triinpoel,
1042 : : const tk::Fields& U,
1043 : : const tk::Fields& G,
1044 : : tk::Fields& R )
1045 : : // *****************************************************************************
1046 : : // Compute right hand side
1047 : : //! \param[in] dsupedge Domain superedges
1048 : : //! \param[in] dsupint Domain superedge integrals
1049 : : //! \param[in] coord Mesh node coordinates
1050 : : //! \param[in] triinpoel Boundary face connectivity
1051 : : //! \param[in] U Solution vector of primitive variables at recent time step
1052 : : //! \param[in] G Gradients of velocity and transported scalars
1053 : : //! \param[in,out] R Right-hand side vector computed
1054 : : // *****************************************************************************
1055 : : {
1056 [ - + ][ - - ]: 4660 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
1057 : : "vector at recent time step incorrect" );
1058 [ - + ][ - - ]: 4660 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
1059 : : "Number of unknowns and/or number of components in right-hand "
1060 : : "side vector incorrect" );
1061 : :
1062 : 4660 : R.fill( 0.0 );
1063 : 4660 : adv( dsupedge, dsupint, coord, triinpoel, U, G, R );
1064 : 4660 : }
1065 : :
1066 : : } // lohner::
|
