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