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-2025 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 : 770 : 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 [ - + ][ - - ]: 770 : Assert( U.nunk() == D.size(), "Size mismatch" );
[ - - ][ - - ]
54 : :
55 : : // Lambda to compute the velocity divergence contribution for edge p-q
56 : 197600 : auto div = [&]( const tk::real d[], std::size_t p, std::size_t q ){
57 : 197600 : return d[0] * (U(p,pos+0) + U(q,pos+0)) +
58 : 197600 : d[1] * (U(p,pos+1) + U(q,pos+1)) +
59 : 197600 : d[2] * (U(p,pos+2) + U(q,pos+2));
60 : 770 : };
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 [ + + ]: 17772 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
75 : 17002 : const auto N = dsupedge[0].data() + e*4;
76 : 17002 : const auto d = dsupint[0].data();
77 : : // edge fluxes
78 : : tk::real f[] = {
79 [ + - ]: 17002 : div( d+(e*6+0)*4, N[0], N[1] ),
80 : 34004 : div( d+(e*6+1)*4, N[1], N[2] ),
81 : 34004 : div( d+(e*6+2)*4, N[2], N[0] ),
82 : 34004 : div( d+(e*6+3)*4, N[0], N[3] ),
83 : 34004 : div( d+(e*6+4)*4, N[1], N[3] ),
84 [ + - ][ + - ]: 17002 : div( d+(e*6+5)*4, N[2], N[3] ) };
[ + - ][ + - ]
[ + - ]
85 : : // edge flux contributions
86 : 17002 : D[N[0]] = D[N[0]] - f[0] + f[2] - f[3];
87 : 17002 : D[N[1]] = D[N[1]] + f[0] - f[1] - f[4];
88 : 17002 : D[N[2]] = D[N[2]] + f[1] - f[2] - f[5];
89 : 17002 : D[N[3]] = D[N[3]] + f[3] + f[4] + f[5];
90 : : }
91 : :
92 : : // domain edge contributions: triangle superedges
93 [ + + ]: 17038 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
94 : 16268 : const auto N = dsupedge[1].data() + e*3;
95 : 16268 : const auto d = dsupint[1].data();
96 : : // edge fluxes
97 : : tk::real f[] = {
98 [ + - ]: 16268 : div( d+(e*3+0)*4, N[0], N[1] ),
99 : 32536 : div( d+(e*3+1)*4, N[1], N[2] ),
100 [ + - ][ + - ]: 16268 : div( d+(e*3+2)*4, N[2], N[0] ) };
101 : : // edge flux contributions
102 : 16268 : D[N[0]] = D[N[0]] - f[0] + f[2];
103 : 16268 : D[N[1]] = D[N[1]] + f[0] - f[1];
104 : 16268 : D[N[2]] = D[N[2]] + f[1] - f[2];
105 : : }
106 : :
107 : : // domain edge contributions: edges
108 [ + + ]: 47554 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
109 : 46784 : const auto N = dsupedge[2].data() + e*2;
110 : 46784 : const auto d = dsupint[2].data();
111 : : // edge flux
112 [ + - ]: 46784 : tk::real f = div( d+e*4, N[0], N[1] );
113 : : // edge flux contributions
114 : 46784 : D[N[0]] -= f;
115 : 46784 : D[N[1]] += f;
116 : : }
117 : :
118 : : // boundary integral
119 : :
120 : 770 : const auto& x = coord[0];
121 : 770 : const auto& y = coord[1];
122 : 770 : const auto& z = coord[2];
123 : :
124 [ + + ]: 67074 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
125 : 66304 : const auto N = triinpoel.data() + e*3;
126 : : tk::real n[3];
127 : 66304 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
128 : 66304 : 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 [ + - ]: 66304 : auto uxA = U(N[0],pos+0);
131 [ + - ]: 66304 : auto uyA = U(N[0],pos+1);
132 [ + - ]: 66304 : auto uzA = U(N[0],pos+2);
133 [ + - ]: 66304 : auto uxB = U(N[1],pos+0);
134 [ + - ]: 66304 : auto uyB = U(N[1],pos+1);
135 [ + - ]: 66304 : auto uzB = U(N[1],pos+2);
136 [ + - ]: 66304 : auto uxC = U(N[2],pos+0);
137 [ + - ]: 66304 : auto uyC = U(N[2],pos+1);
138 [ + - ]: 66304 : auto uzC = U(N[2],pos+2);
139 : 66304 : auto ux = (6.0*uxA + uxB + uxC)/8.0;
140 : 66304 : auto uy = (6.0*uyA + uyB + uyC)/8.0;
141 : 66304 : auto uz = (6.0*uzA + uzB + uzC)/8.0;
142 : 66304 : D[N[0]] += ux*n[0] + uy*n[1] + uz*n[2];
143 : 66304 : ux = (uxA + 6.0*uxB + uxC)/8.0;
144 : 66304 : uy = (uyA + 6.0*uyB + uyC)/8.0;
145 : 66304 : uz = (uzA + 6.0*uzB + uzC)/8.0;
146 : 66304 : D[N[1]] += ux*n[0] + uy*n[1] + uz*n[2];
147 : 66304 : ux = (uxA + uxB + 6.0*uxC)/8.0;
148 : 66304 : uy = (uyA + uyB + 6.0*uyC)/8.0;
149 : 66304 : uz = (uzA + uzB + 6.0*uzC)/8.0;
150 : 66304 : 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 : 770 : }
159 : :
160 : : void
161 : 385 : 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 [ - + ]: 385 : if (G.nprop() == 0) return;
178 : :
179 [ - + ][ - - ]: 385 : Assert( G.nunk() == U.size(), "Size mismatch" );
[ - - ][ - - ]
180 [ - + ][ - - ]: 385 : Assert( G.nprop() > 2, "Size mismatch" );
[ - - ][ - - ]
181 [ - + ][ - - ]: 385 : 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 [ + + ]: 8886 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
196 : 8501 : const auto N = dsupedge[0].data() + e*4;
197 : 8501 : const auto d = dsupint[0].data();
198 : 8501 : tk::real u[] = { U[N[0]], U[N[1]], U[N[2]], U[N[3]] };
199 [ + + ]: 34004 : for (std::size_t j=0; j<3; ++j) {
200 : : tk::real f[] = {
201 : 25503 : d[(e*6+0)*4+j] * (u[1] + u[0]),
202 : 25503 : d[(e*6+1)*4+j] * (u[2] + u[1]),
203 : 25503 : d[(e*6+2)*4+j] * (u[0] + u[2]),
204 : 25503 : d[(e*6+3)*4+j] * (u[3] + u[0]),
205 : 25503 : d[(e*6+4)*4+j] * (u[3] + u[1]),
206 : 25503 : d[(e*6+5)*4+j] * (u[3] + u[2]) };
207 [ + - ][ + - ]: 25503 : G(N[0],j) = G(N[0],j) - f[0] + f[2] - f[3];
208 [ + - ][ + - ]: 25503 : G(N[1],j) = G(N[1],j) + f[0] - f[1] - f[4];
209 [ + - ][ + - ]: 25503 : G(N[2],j) = G(N[2],j) + f[1] - f[2] - f[5];
210 [ + - ][ + - ]: 25503 : G(N[3],j) = G(N[3],j) + f[3] + f[4] + f[5];
211 : : }
212 : : }
213 : :
214 : : // domain edge contributions: triangle superedges
215 [ + + ]: 8519 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
216 : 8134 : const auto N = dsupedge[1].data() + e*3;
217 : 8134 : const auto d = dsupint[1].data();
218 : 8134 : tk::real u[] = { U[N[0]], U[N[1]], U[N[2]] };
219 [ + + ]: 32536 : for (std::size_t j=0; j<3; ++j) {
220 : : tk::real f[] = {
221 : 24402 : d[(e*3+0)*4+j] * (u[1] + u[0]),
222 : 24402 : d[(e*3+1)*4+j] * (u[2] + u[1]),
223 : 24402 : d[(e*3+2)*4+j] * (u[0] + u[2]) };
224 [ + - ][ + - ]: 24402 : G(N[0],j) = G(N[0],j) - f[0] + f[2];
225 [ + - ][ + - ]: 24402 : G(N[1],j) = G(N[1],j) + f[0] - f[1];
226 [ + - ][ + - ]: 24402 : G(N[2],j) = G(N[2],j) + f[1] - f[2];
227 : : }
228 : : }
229 : :
230 : : // domain edge contributions: edges
231 [ + + ]: 23777 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
232 : 23392 : const auto N = dsupedge[2].data() + e*2;
233 : 23392 : const auto d = dsupint[2].data() + e*4;
234 : 23392 : tk::real u[] = { U[N[0]], U[N[1]] };
235 [ + + ]: 93568 : for (std::size_t j=0; j<3; ++j) {
236 : 70176 : tk::real f = d[j] * (u[1] + u[0]);
237 [ + - ]: 70176 : G(N[0],j) -= f;
238 [ + - ]: 70176 : G(N[1],j) += f;
239 : : }
240 : : }
241 : :
242 : : // boundary integral
243 : :
244 : 385 : const auto& x = coord[0];
245 : 385 : const auto& y = coord[1];
246 : 385 : const auto& z = coord[2];
247 : :
248 [ + + ]: 33537 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
249 : 33152 : const auto N = triinpoel.data() + e*3;
250 : : tk::real n[3];
251 : 33152 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
252 : 33152 : 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 : 33152 : auto uA = U[N[0]];
255 : 33152 : auto uB = U[N[1]];
256 : 33152 : auto uC = U[N[2]];
257 : 33152 : auto f = (6.0*uA + uB + uC)/8.0;
258 [ + - ]: 33152 : G(N[0],0) += f * n[0];
259 [ + - ]: 33152 : G(N[0],1) += f * n[1];
260 [ + - ]: 33152 : G(N[0],2) += f * n[2];
261 : 33152 : f = (uA + 6.0*uB + uC)/8.0;
262 [ + - ]: 33152 : G(N[1],0) += f * n[0];
263 [ + - ]: 33152 : G(N[1],1) += f * n[1];
264 [ + - ]: 33152 : G(N[1],2) += f * n[2];
265 : 33152 : f = (uA + uB + 6.0*uC)/8.0;
266 [ + - ]: 33152 : G(N[2],0) += f * n[0];
267 [ + - ]: 33152 : G(N[2],1) += f * n[1];
268 [ + - ]: 33152 : 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 : 385 : 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 [ - + ][ - - ]: 385 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
296 [ - + ][ - - ]: 385 : 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 [ + + ]: 8886 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
311 : 8501 : const auto N = dsupedge[0].data() + e*4;
312 : 8501 : const auto d = dsupint[0].data();
313 [ + + ]: 34004 : for (std::size_t i=0; i<3; ++i) {
314 [ + - ][ + - ]: 25503 : tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1), U(N[3],i+1) };
[ + - ][ + - ]
315 : 25503 : auto i3 = i*3;
316 [ + + ]: 102012 : for (std::size_t j=0; j<3; ++j) {
317 : 76509 : tk::real f[] = { d[(e*6+0)*4+j] * (u[1] + u[0]),
318 : 76509 : d[(e*6+1)*4+j] * (u[2] + u[1]),
319 : 76509 : d[(e*6+2)*4+j] * (u[0] + u[2]),
320 : 76509 : d[(e*6+3)*4+j] * (u[3] + u[0]),
321 : 76509 : d[(e*6+4)*4+j] * (u[3] + u[1]),
322 : 76509 : d[(e*6+5)*4+j] * (u[3] + u[2]) };
323 [ + - ][ + - ]: 76509 : G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2] - f[3];
324 [ + - ][ + - ]: 76509 : G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1] - f[4];
325 [ + - ][ + - ]: 76509 : G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2] - f[5];
326 [ + - ][ + - ]: 76509 : 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 [ + + ]: 8519 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
333 : 8134 : const auto N = dsupedge[1].data() + e*3;
334 : 8134 : const auto d = dsupint[1].data();
335 [ + + ]: 32536 : for (std::size_t i=0; i<3; ++i) {
336 [ + - ][ + - ]: 24402 : tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1) };
[ + - ]
337 : 24402 : auto i3 = i*3;
338 [ + + ]: 97608 : for (std::size_t j=0; j<3; ++j) {
339 : 73206 : tk::real f[] = { d[(e*3+0)*4+j] * (u[1] + u[0]),
340 : 73206 : d[(e*3+1)*4+j] * (u[2] + u[1]),
341 : 73206 : d[(e*3+2)*4+j] * (u[0] + u[2]) };
342 [ + - ][ + - ]: 73206 : G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2];
343 [ + - ][ + - ]: 73206 : G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1];
344 [ + - ][ + - ]: 73206 : G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2];
345 : : }
346 : : }
347 : : }
348 : :
349 : : // domain edge contributions: edges
350 [ + + ]: 23777 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
351 : 23392 : const auto N = dsupedge[2].data() + e*2;
352 : 23392 : const auto d = dsupint[2].data() + e*4;
353 [ + + ]: 93568 : for (std::size_t i=0; i<3; ++i) {
354 [ + - ][ + - ]: 70176 : tk::real u[] = { U(N[0],i+1), U(N[1],i+1) };
355 : 70176 : auto i3 = i*3;
356 [ + + ]: 280704 : for (std::size_t j=0; j<3; ++j) {
357 : 210528 : tk::real f = d[j] * (u[1] + u[0]);
358 [ + - ]: 210528 : G(N[0],i3+j) -= f;
359 [ + - ]: 210528 : G(N[1],i3+j) += f;
360 : : }
361 : : }
362 : : }
363 : :
364 : : // boundary integral
365 : :
366 : 385 : const auto& x = coord[0];
367 : 385 : const auto& y = coord[1];
368 : 385 : const auto& z = coord[2];
369 : :
370 [ + + ]: 33537 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
371 : 33152 : const auto N = triinpoel.data() + e*3;
372 : : tk::real n[3];
373 : 33152 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
374 : 33152 : 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 [ + + ]: 132608 : for (std::size_t i=0; i<3; ++i) {
377 [ + - ][ + - ]: 99456 : tk::real u[] = { U(N[0],i+1), U(N[1],i+1), U(N[2],i+1) };
[ + - ]
378 : 99456 : auto i3 = i*3;
379 : 99456 : auto f = (6.0*u[0] + u[1] + u[2])/8.0;
380 [ + - ]: 99456 : G(N[0],i3+0) += f * n[0];
381 [ + - ]: 99456 : G(N[0],i3+1) += f * n[1];
382 [ + - ]: 99456 : G(N[0],i3+2) += f * n[2];
383 : 99456 : f = (u[0] + 6.0*u[1] + u[2])/8.0;
384 [ + - ]: 99456 : G(N[1],i3+0) += f * n[0];
385 [ + - ]: 99456 : G(N[1],i3+1) += f * n[1];
386 [ + - ]: 99456 : G(N[1],i3+2) += f * n[2];
387 : 99456 : f = (u[0] + u[1] + 6.0*u[2])/8.0;
388 [ + - ]: 99456 : G(N[2],i3+0) += f * n[0];
389 [ + - ]: 99456 : G(N[2],i3+1) += f * n[1];
390 [ + - ]: 99456 : 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 : 385 : }
400 : :
401 : : static tk::real
402 : 889200 : 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 : 889200 : auto inv = U(p,i+1)*U(p,j+1) + U(q,i+1)*U(q,j+1);
420 : :
421 : 889200 : auto eps = std::numeric_limits< tk::real >::epsilon();
422 : 889200 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
423 [ + + ]: 889200 : if (mu < eps) return -inv;
424 : :
425 : 655731 : auto vis = G(p,i*3+j) + G(p,j*3+i) + G(q,i*3+j) + G(q,j*3+i);
426 [ + + ]: 655731 : if (i == j) {
427 : 218577 : 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 : 655731 : return mu*vis - inv;
430 : : }
431 : :
432 : : static tk::real
433 : 895104 : 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 : 895104 : auto inv = U(p,i+1)*U(p,j+1);
449 : :
450 : 895104 : auto eps = std::numeric_limits< tk::real >::epsilon();
451 : 895104 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
452 [ + + ]: 895104 : if (mu < eps) return -inv;
453 : :
454 : 636120 : auto vis = G(p,i*3+j) + G(p,j*3+i);
455 [ + + ]: 636120 : if (i == j) {
456 : 212040 : vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) );
457 : : }
458 : 636120 : return mu*vis - inv;
459 : : }
460 : :
461 : : void
462 : 385 : 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 [ - + ][ - - ]: 385 : Assert( F.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
482 [ - + ][ - - ]: 385 : Assert( F.nprop() == 3, "Size mismatch" );
[ - - ][ - - ]
483 [ - + ][ - - ]: 385 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
484 [ - + ][ - - ]: 385 : 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 [ + + ]: 8886 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
499 : 8501 : const auto N = dsupedge[0].data() + e*4;
500 : 8501 : const auto d = dsupint[0].data();
501 [ + + ]: 34004 : for (std::size_t i=0; i<3; ++i) {
502 [ + + ]: 102012 : for (std::size_t j=0; j<3; ++j) {
503 [ + - ]: 76509 : tk::real f[] = { d[(e*6+0)*4+j] * flux(U,G,i,j,N[1],N[0]),
504 : 153018 : d[(e*6+1)*4+j] * flux(U,G,i,j,N[2],N[1]),
505 : 153018 : d[(e*6+2)*4+j] * flux(U,G,i,j,N[0],N[2]),
506 : 153018 : d[(e*6+3)*4+j] * flux(U,G,i,j,N[3],N[0]),
507 : 153018 : d[(e*6+4)*4+j] * flux(U,G,i,j,N[3],N[1]),
508 [ + - ][ + - ]: 76509 : d[(e*6+5)*4+j] * flux(U,G,i,j,N[3],N[2]) };
[ + - ][ + - ]
[ + - ]
509 [ + - ][ + - ]: 76509 : F(N[0],i) = F(N[0],i) - f[0] + f[2] - f[3];
510 [ + - ][ + - ]: 76509 : F(N[1],i) = F(N[1],i) + f[0] - f[1] - f[4];
511 [ + - ][ + - ]: 76509 : F(N[2],i) = F(N[2],i) + f[1] - f[2] - f[5];
512 [ + - ][ + - ]: 76509 : 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 [ + + ]: 8519 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
519 : 8134 : const auto N = dsupedge[1].data() + e*3;
520 : 8134 : const auto d = dsupint[1].data();
521 [ + + ]: 32536 : for (std::size_t i=0; i<3; ++i) {
522 [ + + ]: 97608 : for (std::size_t j=0; j<3; ++j) {
523 [ + - ]: 73206 : tk::real f[] = { d[(e*3+0)*4+j] * flux(U,G,i,j,N[1],N[0]),
524 : 146412 : d[(e*3+1)*4+j] * flux(U,G,i,j,N[2],N[1]),
525 [ + - ][ + - ]: 73206 : d[(e*3+2)*4+j] * flux(U,G,i,j,N[0],N[2]), };
526 [ + - ][ + - ]: 73206 : F(N[0],i) = F(N[0],i) - f[0] + f[2];
527 [ + - ][ + - ]: 73206 : F(N[1],i) = F(N[1],i) + f[0] - f[1];
528 [ + - ][ + - ]: 73206 : F(N[2],i) = F(N[2],i) + f[1] - f[2];
529 : : }
530 : : }
531 : : }
532 : :
533 : : // domain edge contributions: edges
534 [ + + ]: 23777 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
535 : 23392 : const auto N = dsupedge[2].data() + e*2;
536 : 23392 : const auto d = dsupint[2].data() + e*4;
537 [ + + ]: 93568 : for (std::size_t i=0; i<3; ++i) {
538 [ + + ]: 280704 : for (std::size_t j=0; j<3; ++j) {
539 : 210528 : tk::real f = d[j] * flux(U,G,i,j,N[1],N[0]);
540 : 210528 : F(N[0],i) -= f;
541 : 210528 : F(N[1],i) += f;
542 : : }
543 : : }
544 : : }
545 : :
546 : : // boundary integral
547 : :
548 : 385 : const auto& x = coord[0];
549 : 385 : const auto& y = coord[1];
550 : 385 : const auto& z = coord[2];
551 : :
552 [ + + ]: 33537 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
553 : 33152 : const auto N = triinpoel.data() + e*3;
554 : : tk::real n[3];
555 : 33152 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
556 : 33152 : 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 [ + + ]: 132608 : for (std::size_t i=0; i<3; ++i) {
559 [ + - ]: 99456 : auto fxA = flux(U,G,i,0,N[0]);
560 [ + - ]: 99456 : auto fyA = flux(U,G,i,1,N[0]);
561 [ + - ]: 99456 : auto fzA = flux(U,G,i,2,N[0]);
562 [ + - ]: 99456 : auto fxB = flux(U,G,i,0,N[1]);
563 [ + - ]: 99456 : auto fyB = flux(U,G,i,1,N[1]);
564 [ + - ]: 99456 : auto fzB = flux(U,G,i,2,N[1]);
565 [ + - ]: 99456 : auto fxC = flux(U,G,i,0,N[2]);
566 [ + - ]: 99456 : auto fyC = flux(U,G,i,1,N[2]);
567 [ + - ]: 99456 : auto fzC = flux(U,G,i,2,N[2]);
568 : 99456 : auto fx = (6.0*fxA + fxB + fxC)/8.0;
569 : 99456 : auto fy = (6.0*fyA + fyB + fyC)/8.0;
570 : 99456 : auto fz = (6.0*fzA + fzB + fzC)/8.0;
571 [ + - ]: 99456 : F(N[0],i) += fx*n[0] + fy*n[1] + fz*n[2];
572 : 99456 : fx = (fxA + 6.0*fxB + fxC)/8.0;
573 : 99456 : fy = (fyA + 6.0*fyB + fyC)/8.0;
574 : 99456 : fz = (fzA + 6.0*fzB + fzC)/8.0;
575 [ + - ]: 99456 : F(N[1],i) += fx*n[0] + fy*n[1] + fz*n[2];
576 : 99456 : fx = (fxA + fxB + 6.0*fxC)/8.0;
577 : 99456 : fy = (fyA + fyB + 6.0*fyC)/8.0;
578 : 99456 : fz = (fzA + fzB + 6.0*fzC)/8.0;
579 [ + - ]: 99456 : 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 : 385 : }
589 : :
590 : : void
591 : 1040 : 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 [ - + ]: 1040 : if (G.nprop() == 0) return;
608 : :
609 [ - + ][ - - ]: 1040 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
610 [ - + ][ - - ]: 1040 : Assert( G.nprop() > 2, "Size mismatch" );
[ - - ][ - - ]
611 [ - + ][ - - ]: 1040 : Assert( G.nprop() % 3 == 0, "Size mismatch" );
[ - - ][ - - ]
612 [ - + ][ - - ]: 1040 : Assert( G.nprop() == U.nprop()*3, "Size mismatch" );
[ - - ][ - - ]
613 : :
614 : 1040 : 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 [ + + ]: 175120 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
629 : 174080 : const auto N = dsupedge[0].data() + e*4;
630 : 174080 : const auto d = dsupint[0].data();
631 [ + + ]: 910880 : for (std::size_t c=0; c<ncomp; ++c) {
632 [ + - ][ + - ]: 736800 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c), U(N[3],c) };
[ + - ][ + - ]
633 : 736800 : auto g = c*3;
634 [ + + ]: 2947200 : for (std::size_t j=0; j<3; ++j) {
635 : : tk::real f[] = {
636 : 2210400 : d[(e*6+0)*4+j] * (u[1] + u[0]),
637 : 2210400 : d[(e*6+1)*4+j] * (u[2] + u[1]),
638 : 2210400 : d[(e*6+2)*4+j] * (u[0] + u[2]),
639 : 2210400 : d[(e*6+3)*4+j] * (u[3] + u[0]),
640 : 2210400 : d[(e*6+4)*4+j] * (u[3] + u[1]),
641 : 2210400 : d[(e*6+5)*4+j] * (u[3] + u[2]) };
642 [ + - ][ + - ]: 2210400 : G(N[0],g+j) = G(N[0],g+j) - f[0] + f[2] - f[3];
643 [ + - ][ + - ]: 2210400 : G(N[1],g+j) = G(N[1],g+j) + f[0] - f[1] - f[4];
644 [ + - ][ + - ]: 2210400 : G(N[2],g+j) = G(N[2],g+j) + f[1] - f[2] - f[5];
645 [ + - ][ + - ]: 2210400 : 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 [ + + ]: 149440 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
652 : 148400 : const auto N = dsupedge[1].data() + e*3;
653 : 148400 : const auto d = dsupint[1].data();
654 [ + + ]: 779600 : for (std::size_t c=0; c<ncomp; ++c) {
655 [ + - ][ + - ]: 631200 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c) };
[ + - ]
656 : 631200 : auto g = c*3;
657 [ + + ]: 2524800 : for (std::size_t j=0; j<3; ++j) {
658 : : tk::real f[] = {
659 : 1893600 : d[(e*3+0)*4+j] * (u[1] + u[0]),
660 : 1893600 : d[(e*3+1)*4+j] * (u[2] + u[1]),
661 : 1893600 : d[(e*3+2)*4+j] * (u[0] + u[2]) };
662 [ + - ][ + - ]: 1893600 : G(N[0],g+j) = G(N[0],g+j) - f[0] + f[2];
663 [ + - ][ + - ]: 1893600 : G(N[1],g+j) = G(N[1],g+j) + f[0] - f[1];
664 [ + - ][ + - ]: 1893600 : G(N[2],g+j) = G(N[2],g+j) + f[1] - f[2];
665 : : }
666 : : }
667 : : }
668 : :
669 : : // domain edge contributions: edges
670 [ + + ]: 455440 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
671 : 454400 : const auto N = dsupedge[2].data() + e*2;
672 : 454400 : const auto d = dsupint[2].data() + e*4;
673 [ + + ]: 2399680 : for (std::size_t c=0; c<ncomp; ++c) {
674 [ + - ][ + - ]: 1945280 : tk::real u[] = { U(N[0],c), U(N[1],c) };
675 : 1945280 : auto g = c*3;
676 [ + + ]: 7781120 : for (std::size_t j=0; j<3; ++j) {
677 : 5835840 : tk::real f = d[j] * (u[1] + u[0]);
678 [ + - ]: 5835840 : G(N[0],g+j) -= f;
679 [ + - ]: 5835840 : G(N[1],g+j) += f;
680 : : }
681 : : }
682 : : }
683 : :
684 : : // boundary integral
685 : :
686 : 1040 : const auto& x = coord[0];
687 : 1040 : const auto& y = coord[1];
688 : 1040 : const auto& z = coord[2];
689 : :
690 [ + + ]: 585200 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
691 : 584160 : const auto N = triinpoel.data() + e*3;
692 : : tk::real n[3];
693 : 584160 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
694 : 584160 : 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 [ + + ]: 3112640 : for (std::size_t c=0; c<ncomp; ++c) {
697 : 2528480 : auto g = c*3;
698 [ + - ]: 2528480 : auto uA = U(N[0],c);
699 [ + - ]: 2528480 : auto uB = U(N[1],c);
700 [ + - ]: 2528480 : auto uC = U(N[2],c);
701 : 2528480 : auto f = (6.0*uA + uB + uC)/8.0;
702 [ + - ]: 2528480 : G(N[0],g+0) += f * n[0];
703 [ + - ]: 2528480 : G(N[0],g+1) += f * n[1];
704 [ + - ]: 2528480 : G(N[0],g+2) += f * n[2];
705 : 2528480 : f = (uA + 6.0*uB + uC)/8.0;
706 [ + - ]: 2528480 : G(N[1],g+0) += f * n[0];
707 [ + - ]: 2528480 : G(N[1],g+1) += f * n[1];
708 [ + - ]: 2528480 : G(N[1],g+2) += f * n[2];
709 : 2528480 : f = (uA + uB + 6.0*uC)/8.0;
710 [ + - ]: 2528480 : G(N[2],g+0) += f * n[0];
711 [ + - ]: 2528480 : G(N[2],g+1) += f * n[1];
712 [ + - ]: 2528480 : 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 : 2992380 : 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 : 2992380 : auto ncomp = U.nprop();
741 : :
742 : 2992380 : auto nx = supint[0];
743 : 2992380 : auto ny = supint[1];
744 : 2992380 : auto nz = supint[2];
745 : :
746 : : // left state
747 : 2992380 : auto pL = U(p,0);
748 : 2992380 : auto uL = U(p,1);
749 : 2992380 : auto vL = U(p,2);
750 : 2992380 : auto wL = U(p,3);
751 : 2992380 : auto vnL = uL*nx + vL*ny + wL*nz;
752 : :
753 : : // right state
754 : 2992380 : auto pR = U(q,0);
755 : 2992380 : auto uR = U(q,1);
756 : 2992380 : auto vR = U(q,2);
757 : 2992380 : auto wR = U(q,3);
758 : 2992380 : auto vnR = uR*nx + vR*ny + wR*nz;
759 : :
760 : 2992380 : auto s = g_cfg.get< tag::soundspeed >();
761 : 2992380 : auto s2 = s*s;
762 : :
763 : : // viscosity
764 : 2992380 : auto v = supint[3] * g_cfg.get< tag::mat_dyn_viscosity >();
765 : :
766 : : // stabilization
767 : 2992380 : tk::real aw = 0.0;
768 [ + - ]: 2992380 : if (g_cfg.get< tag::stab >()) {
769 : 2992380 : aw = std::abs( vnL + vnR ) / 2.0;
770 : : }
771 : :
772 : : // artificial viscosity
773 [ + + ]: 2992380 : if (g_cfg.get< tag::stab2 >()) {
774 : 1193940 : auto len = tk::length( nx, ny, nz );
775 : 1193940 : auto sl = std::abs(vnL) + s*len;
776 : 1193940 : auto sr = std::abs(vnR) + s*len;
777 : 1193940 : aw += g_cfg.get< tag::stab2coef >() * std::max(sl,sr);
778 : : }
779 : :
780 : : // flow
781 : 2992380 : auto pf = pL + pR;
782 : 2992380 : f[0] = (vnL + vnR + aw*(pR - pL))*s2;
783 : 2992380 : f[1] = uL*vnL + uR*vnR + pf*nx + (aw-v)*(uR - uL);
784 : 2992380 : f[2] = vL*vnL + vR*vnR + pf*ny + (aw-v)*(vR - vL);
785 : 2992380 : f[3] = wL*vnL + wR*vnR + pf*nz + (aw-v)*(wR - wL);
786 : :
787 : : // diffusion
788 : 2992380 : auto d = supint[3] * g_cfg.get< tag::mat_dyn_diffusivity >();
789 : :
790 : : // scalar
791 [ + + ]: 4186320 : for (std::size_t c=4; c<ncomp; ++c) {
792 : 1193940 : f[c] = U(p,c)*vnL + U(q,c)*vnR + (aw-d)*(U(q,c) - U(p,c));
793 : : }
794 : 2992380 : }
795 : :
796 : : static void
797 : 1944080 : adv_damp4( const tk::real supint[],
798 : : const tk::Fields& U,
799 : : const tk::Fields& G,
800 : : const std::array< std::vector< tk::real >, 3 >& coord,
801 : : std::size_t p,
802 : : std::size_t q,
803 : : tk::real f[] )
804 : : // *****************************************************************************
805 : : //! Compute advection fluxes on a single edge using 4th-order damping
806 : : //! \param[in] supint Edge integral
807 : : //! \param[in] U Velocity and transported scalars at recent time step
808 : : //! \param[in] G Gradients of velocity and transported scalars
809 : : //! \param[in] coord Mesh node coordinates
810 : : //! \param[in] p Left node index of edge
811 : : //! \param[in] q Right node index of edge
812 : : //! \param[in,out] f Flux computed
813 : : // *****************************************************************************
814 : : {
815 : 1944080 : const auto ncomp = U.nprop();
816 : :
817 : 1944080 : const auto& x = coord[0];
818 : 1944080 : const auto& y = coord[1];
819 : 1944080 : const auto& z = coord[2];
820 : :
821 : : // edge vector
822 : 1944080 : auto dx = x[q] - x[p];
823 : 1944080 : auto dy = y[q] - y[p];
824 : 1944080 : auto dz = z[q] - z[p];
825 : :
826 : : #if defined(__clang__)
827 : : #pragma clang diagnostic push
828 : : #pragma clang diagnostic ignored "-Wvla"
829 : : #pragma clang diagnostic ignored "-Wvla-extension"
830 : : #elif defined(STRICT_GNUC)
831 : : #pragma GCC diagnostic push
832 : : #pragma GCC diagnostic ignored "-Wvla"
833 : : #endif
834 : :
835 : 1944080 : tk::real uL[ncomp], uR[ncomp];
836 [ + + ]: 10203760 : for (std::size_t i=0; i<ncomp; ++i) {
837 [ + - ]: 8259680 : uL[i] = U(p,i);
838 [ + - ]: 8259680 : uR[i] = U(q,i);
839 : : }
840 : :
841 : : // MUSCL reconstruction in edge-end points
842 [ + + ]: 10203760 : for (std::size_t c=0; c<ncomp; ++c) {
843 [ + - ][ + - ]: 8259680 : auto g1 = G(p,c*3+0)*dx + G(p,c*3+1)*dy + G(p,c*3+2)*dz;
[ + - ]
844 [ + - ][ + - ]: 8259680 : auto g2 = G(q,c*3+0)*dx + G(q,c*3+1)*dy + G(q,c*3+2)*dz;
[ + - ]
845 : 8259680 : auto delta2 = uR[c] - uL[c];
846 : 8259680 : auto delta1 = 2.0 * g1 - delta2;
847 : 8259680 : auto delta3 = 2.0 * g2 - delta2;
848 : :
849 : : // van Leer limiter
850 : 8259680 : auto rL = (delta2 + muscl_eps) / (delta1 + muscl_eps);
851 : 8259680 : auto rR = (delta2 + muscl_eps) / (delta3 + muscl_eps);
852 : 8259680 : auto rLinv = (delta1 + muscl_eps) / (delta2 + muscl_eps);
853 : 8259680 : auto rRinv = (delta3 + muscl_eps) / (delta2 + muscl_eps);
854 : 8259680 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
855 : 8259680 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
856 : 8259680 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
857 : 8259680 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
858 : : // update unknowns with reconstructed unknowns
859 : 8259680 : uL[c] += 0.25*(delta1*(1.0-muscl_const)*phiL +
860 : 8259680 : delta2*(1.0+muscl_const)*phi_L_inv);
861 : 8259680 : uR[c] -= 0.25*(delta3*(1.0-muscl_const)*phiR +
862 : 8259680 : delta2*(1.0+muscl_const)*phi_R_inv);
863 : : }
864 : :
865 : 1944080 : auto nx = supint[0];
866 : 1944080 : auto ny = supint[1];
867 : 1944080 : auto nz = supint[2];
868 : :
869 : : // normal velocities
870 : 1944080 : auto vnL = uL[1]*nx + uL[2]*ny + uL[3]*nz;
871 : 1944080 : auto vnR = uR[1]*nx + uR[2]*ny + uR[3]*nz;
872 : :
873 : 1944080 : auto s = g_cfg.get< tag::soundspeed >();
874 : 1944080 : auto s2 = s*s;
875 : :
876 : : // viscosity
877 : 1944080 : auto v = supint[3] * g_cfg.get< tag::mat_dyn_viscosity >();
878 : :
879 : : // stabilization
880 : 1944080 : tk::real aw = 0.0;
881 [ + - ]: 1944080 : if (g_cfg.get< tag::stab >()) {
882 : 1944080 : aw = std::abs( vnL + vnR ) / 2.0;
883 : : }
884 : :
885 : : // artificial viscosity
886 [ + + ]: 1944080 : if (g_cfg.get< tag::stab2 >()) {
887 : 1510400 : auto len = tk::length( nx, ny, nz );
888 : 1510400 : auto sl = std::abs(vnL) + s*len;
889 : 1510400 : auto sr = std::abs(vnR) + s*len;
890 : 1510400 : aw += g_cfg.get< tag::stab2coef >() * std::max(sl,sr);
891 : : }
892 : :
893 : : // flow
894 : 1944080 : auto pf = uL[0] + uR[0];
895 : 1944080 : f[0] = (vnL + vnR + aw*(uR[0]-uL[0]))*s2;
896 [ + - ][ + - ]: 1944080 : f[1] = uL[1]*vnL + uR[1]*vnR + pf*nx + aw*(uR[1]-uL[1]) - v*(U(q,1)-U(p,1));
897 [ + - ][ + - ]: 1944080 : f[2] = uL[2]*vnL + uR[2]*vnR + pf*ny + aw*(uR[2]-uL[2]) - v*(U(q,2)-U(p,2));
898 [ + - ][ + - ]: 1944080 : f[3] = uL[3]*vnL + uR[3]*vnR + pf*nz + aw*(uR[3]-uL[3]) - v*(U(q,3)-U(p,3));
899 : :
900 : : // diffusion
901 : 1944080 : auto d = supint[3] * g_cfg.get< tag::mat_dyn_diffusivity >();
902 : :
903 : : // scalar
904 [ + + ]: 2427440 : for (std::size_t c=4; c<ncomp; ++c) {
905 [ + - ][ + - ]: 483360 : f[c] = uL[c]*vnL + uR[c]*vnR + aw*(uR[c]-uL[c]) - d*(U(q,c)-U(p,c));
906 : : }
907 : :
908 : : #if defined(__clang__)
909 : : #pragma clang diagnostic pop
910 : : #elif defined(STRICT_GNUC)
911 : : #pragma GCC diagnostic pop
912 : : #endif
913 : 3888160 : }
914 : :
915 : : static void
916 : 8600 : adv( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
917 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
918 : : const std::array< std::vector< tk::real >, 3 >& coord,
919 : : const std::vector< std::size_t >& triinpoel,
920 : : const tk::Fields& U,
921 : : const tk::Fields& G,
922 : : // cppcheck-suppress constParameter
923 : : tk::Fields& R )
924 : : // *****************************************************************************
925 : : //! Add advection to rhs
926 : : //! \param[in] dsupedge Domain superedges
927 : : //! \param[in] dsupint Domain superedge integrals
928 : : //! \param[in] coord Mesh node coordinates
929 : : //! \param[in] triinpoel Boundary face connectivity
930 : : //! \param[in] U Velocity and transported scalars at recent time step
931 : : //! \param[in] G Gradients of velocity and transported scalars
932 : : //! \param[in,out] R Right-hand side vector added to
933 : : // *****************************************************************************
934 : : {
935 : 8600 : const auto ncomp = U.nprop();
936 : :
937 : : // configure advection
938 : 0 : auto adv = [](){
939 : 8600 : const auto& flux = g_cfg.get< tag::flux >();
940 [ + + ]: 8600 : if (flux == "damp2") return adv_damp2;
941 [ + - ]: 1040 : else if (flux == "damp4") return adv_damp4;
942 [ - - ][ - - ]: 0 : else Throw( "Flux not correctly configured" );
[ - - ]
943 [ + - ]: 8600 : }();
944 : :
945 : : #if defined(__clang__)
946 : : #pragma clang diagnostic push
947 : : #pragma clang diagnostic ignored "-Wvla"
948 : : #pragma clang diagnostic ignored "-Wvla-extension"
949 : : #elif defined(STRICT_GNUC)
950 : : #pragma GCC diagnostic push
951 : : #pragma GCC diagnostic ignored "-Wvla"
952 : : #endif
953 : :
954 : : // domain integral
955 : :
956 : : // domain edge contributions: tetrahedron superedges
957 [ + + ]: 440090 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
958 : 431490 : const auto N = dsupedge[0].data() + e*4;
959 : 431490 : tk::real f[6][ncomp];
960 : 431490 : const auto d = dsupint[0].data();
961 [ + - ]: 431490 : adv( d+(e*6+0)*4, U, G, coord, N[0], N[1], f[0] );
962 [ + - ]: 431490 : adv( d+(e*6+1)*4, U, G, coord, N[1], N[2], f[1] );
963 [ + - ]: 431490 : adv( d+(e*6+2)*4, U, G, coord, N[2], N[0], f[2] );
964 [ + - ]: 431490 : adv( d+(e*6+3)*4, U, G, coord, N[0], N[3], f[3] );
965 [ + - ]: 431490 : adv( d+(e*6+4)*4, U, G, coord, N[1], N[3], f[4] );
966 [ + - ]: 431490 : adv( d+(e*6+5)*4, U, G, coord, N[2], N[3], f[5] );
967 [ + + ]: 2296690 : for (std::size_t c=0; c<ncomp; ++c) {
968 [ + - ][ + - ]: 1865200 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
969 [ + - ][ + - ]: 1865200 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
970 [ + - ][ + - ]: 1865200 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
971 [ + - ][ + - ]: 1865200 : R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
972 : : }
973 : 431490 : }
974 : :
975 : : // domain edge contributions: triangle superedges
976 [ + + ]: 412120 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
977 : 403520 : const auto N = dsupedge[1].data() + e*3;
978 : 403520 : tk::real f[3][ncomp];
979 : 403520 : const auto d = dsupint[1].data();
980 [ + - ]: 403520 : adv( d+(e*3+0)*4, U, G, coord, N[0], N[1], f[0] );
981 [ + - ]: 403520 : adv( d+(e*3+1)*4, U, G, coord, N[1], N[2], f[1] );
982 [ + - ]: 403520 : adv( d+(e*3+2)*4, U, G, coord, N[2], N[0], f[2] );
983 [ + + ]: 2149820 : for (std::size_t c=0; c<ncomp; ++c) {
984 [ + - ][ + - ]: 1746300 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
985 [ + - ][ + - ]: 1746300 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
986 [ + - ][ + - ]: 1746300 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
987 : : }
988 : 403520 : }
989 : :
990 : : // domain edge contributions: edges
991 [ + + ]: 1145560 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
992 : 1136960 : const auto N = dsupedge[2].data() + e*2;
993 : 1136960 : tk::real f[ncomp];
994 : 1136960 : const auto d = dsupint[2].data();
995 [ + - ]: 1136960 : adv( d+e*4, U, G, coord, N[0], N[1], f );
996 [ + + ]: 6130000 : for (std::size_t c=0; c<ncomp; ++c) {
997 [ + - ]: 4993040 : R(N[0],c) -= f[c];
998 [ + - ]: 4993040 : R(N[1],c) += f[c];
999 : : }
1000 : 1136960 : }
1001 : :
1002 : : // boundary integral
1003 : :
1004 : 8600 : const auto& x = coord[0];
1005 : 8600 : const auto& y = coord[1];
1006 : 8600 : const auto& z = coord[2];
1007 : :
1008 : 8600 : auto s = g_cfg.get< tag::soundspeed >();
1009 : 8600 : auto s2 = s * s;
1010 : :
1011 [ + + ]: 1743920 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
1012 : 3470640 : const auto N = triinpoel.data() + e*3;
1013 : : tk::real n[3];
1014 : 1735320 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
1015 : 1735320 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
1016 : : n[0], n[1], n[2] );
1017 : 1735320 : tk::real f[ncomp][3];
1018 : :
1019 [ + - ]: 1735320 : auto p = U(N[0],0);
1020 [ + - ]: 1735320 : auto u = U(N[0],1);
1021 [ + - ]: 1735320 : auto v = U(N[0],2);
1022 [ + - ]: 1735320 : auto w = U(N[0],3);
1023 : 1735320 : auto vn = n[0]*u + n[1]*v + n[2]*w;
1024 : : // flow
1025 : 1735320 : f[0][0] = vn * s2;
1026 : 1735320 : f[1][0] = u*vn + p*n[0];
1027 : 1735320 : f[2][0] = v*vn + p*n[1];
1028 : 1735320 : f[3][0] = w*vn + p*n[2];
1029 : : // scalar
1030 [ + - ][ + + ]: 2358800 : for (std::size_t c=4; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;
1031 : :
1032 [ + - ]: 1735320 : p = U(N[1],0);
1033 [ + - ]: 1735320 : u = U(N[1],1);
1034 [ + - ]: 1735320 : v = U(N[1],2);
1035 [ + - ]: 1735320 : w = U(N[1],3);
1036 : 1735320 : vn = n[0]*u + n[1]*v + n[2]*w;
1037 : : // flow
1038 : 1735320 : f[0][1] = vn * s2;
1039 : 1735320 : f[1][1] = u*vn + p*n[0];
1040 : 1735320 : f[2][1] = v*vn + p*n[1];
1041 : 1735320 : f[3][1] = w*vn + p*n[2];
1042 : : // scalar
1043 [ + - ][ + + ]: 2358800 : for (std::size_t c=4; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;
1044 : :
1045 [ + - ]: 1735320 : p = U(N[2],0);
1046 [ + - ]: 1735320 : u = U(N[2],1);
1047 [ + - ]: 1735320 : v = U(N[2],2);
1048 [ + - ]: 1735320 : w = U(N[2],3);
1049 : 1735320 : vn = n[0]*u + n[1]*v + n[2]*w;
1050 : : // flow
1051 : 1735320 : f[0][2] = vn * s2;
1052 : 1735320 : f[1][2] = u*vn + p*n[0];
1053 : 1735320 : f[2][2] = v*vn + p*n[1];
1054 : 1735320 : f[3][2] = w*vn + p*n[2];
1055 : : // scalar
1056 [ + - ][ + + ]: 2358800 : for (std::size_t c=4; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;
1057 : :
1058 [ + + ]: 9300080 : for (std::size_t c=0; c<ncomp; ++c) {
1059 [ + - ]: 7564760 : R(N[0],c) += (6.0*f[c][0] + f[c][1] + f[c][2])/8.0;
1060 [ + - ]: 7564760 : R(N[1],c) += (f[c][0] + 6.0*f[c][1] + f[c][2])/8.0;
1061 [ + - ]: 7564760 : R(N[2],c) += (f[c][0] + f[c][1] + 6.0*f[c][2])/8.0;
1062 : : }
1063 : 1735320 : }
1064 : :
1065 : : #if defined(__clang__)
1066 : : #pragma clang diagnostic pop
1067 : : #elif defined(STRICT_GNUC)
1068 : : #pragma GCC diagnostic pop
1069 : : #endif
1070 : 8600 : }
1071 : :
1072 : : static void
1073 : 8600 : src( const std::array< std::vector< tk::real >, 3 >& coord,
1074 : : const std::vector< tk::real >& v,
1075 : : tk::real t,
1076 : : tk::Fields& R )
1077 : : // *****************************************************************************
1078 : : // Compute source integral
1079 : : //! \param[in] coord Mesh node coordinates
1080 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
1081 : : //! \param[in] t Physical time
1082 : : //! \param[in,out] R Right-hand side vector computed
1083 : : // *****************************************************************************
1084 : : {
1085 [ + - ]: 8600 : auto src = problems::SRC();
1086 [ + + ]: 8600 : if (!src) return;
1087 : :
1088 : 3020 : const auto& x = coord[0];
1089 : 3020 : const auto& y = coord[1];
1090 : 3020 : const auto& z = coord[2];
1091 : :
1092 [ + + ]: 372160 : for (std::size_t p=0; p<R.nunk(); ++p) {
1093 [ + - ]: 369140 : auto s = src( x[p], y[p], z[p], t );
1094 [ + - ][ + + ]: 2214840 : for (std::size_t c=0; c<s.size(); ++c) R(p,c) -= s[c] * v[p];
1095 : 369140 : }
1096 [ + + ]: 8600 : }
1097 : :
1098 : : void
1099 : 8600 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
1100 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
1101 : : const std::array< std::vector< tk::real >, 3 >& coord,
1102 : : const std::vector< std::size_t >& triinpoel,
1103 : : const std::vector< tk::real >& v,
1104 : : tk::real t,
1105 : : const tk::Fields& U,
1106 : : const tk::Fields& G,
1107 : : tk::Fields& R )
1108 : : // *****************************************************************************
1109 : : // Compute right hand side
1110 : : //! \param[in] dsupedge Domain superedges
1111 : : //! \param[in] dsupint Domain superedge integrals
1112 : : //! \param[in] coord Mesh node coordinates
1113 : : //! \param[in] triinpoel Boundary face connectivity
1114 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
1115 : : //! \param[in] t Physical time
1116 : : //! \param[in] U Solution vector of primitive variables at recent time step
1117 : : //! \param[in] G Gradients of velocity and transported scalars
1118 : : //! \param[in,out] R Right-hand side vector computed
1119 : : // *****************************************************************************
1120 : : {
1121 [ - + ][ - - ]: 8600 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
1122 : : "vector at recent time step incorrect" );
1123 [ - + ][ - - ]: 8600 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
1124 : : "Number of unknowns and/or number of components in right-hand "
1125 : : "side vector incorrect" );
1126 : :
1127 : 8600 : R.fill( 0.0 );
1128 : 8600 : adv( dsupedge, dsupint, coord, triinpoel, U, G, R );
1129 : 8600 : src( coord, v, t, R );
1130 : 8600 : }
1131 : :
1132 : : } // lohner::
|
