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-2024 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 : 1735375 : 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 : 1735375 : tk::real div = d[0] * (U(p,0) + U(q,0)) +
59 : 1735375 : d[1] * (U(p,1) + U(q,1)) +
60 : 1735375 : d[2] * (U(p,2) + U(q,2));
61 : :
62 [ + + ]: 1735375 : if (!stab) return div;
63 : :
64 : 1529560 : const auto& x = coord[0];
65 : 1529560 : const auto& y = coord[1];
66 : 1529560 : const auto& z = coord[2];
67 : :
68 : 1529560 : auto dx = x[p] - x[q];
69 : 1529560 : auto dy = y[p] - y[q];
70 : 1529560 : auto dz = z[p] - z[q];
71 : 1529560 : auto dl = std::sqrt( dx*dx + dy*dy + dz*dz );
72 : 1529560 : auto p2 = P[q] - P[p];
73 : 1529560 : auto D = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] );
74 : :
75 : 1529560 : auto dpx = G(p,0) + G(q,0);
76 : 1529560 : auto dpy = G(p,1) + G(q,1);
77 : 1529560 : auto dpz = G(p,2) + G(q,2);
78 : 1529560 : auto p4 = 0.5 * (dx*dpx + dy*dpy + dz*dpz);
79 : :
80 : 1529560 : div += D*dt/dl*(p2 + p4);
81 : :
82 : 1529560 : return div;
83 : : }
84 : :
85 : : void
86 : 4328 : 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 [ - + ][ - - ]: 4328 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
111 [ - + ][ - - ]: 4328 : 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 [ + + ]: 159470 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
126 : 155142 : const auto N = dsupedge[0].data() + e*4;
127 : 155142 : const auto d = dsupint[0].data();
128 : : // edge fluxes
129 : : tk::real f[] = {
130 [ + - ]: 155142 : div( coord, d+(e*6+0)*5, dt, P, G, U, N[0], N[1], stab ),
131 : 310284 : div( coord, d+(e*6+1)*5, dt, P, G, U, N[1], N[2], stab ),
132 : 310284 : div( coord, d+(e*6+2)*5, dt, P, G, U, N[2], N[0], stab ),
133 : 310284 : div( coord, d+(e*6+3)*5, dt, P, G, U, N[0], N[3], stab ),
134 : 310284 : div( coord, d+(e*6+4)*5, dt, P, G, U, N[1], N[3], stab ),
135 [ + - ][ + - ]: 155142 : div( coord, d+(e*6+5)*5, dt, P, G, U, N[2], N[3], stab ) };
[ + - ][ + - ]
[ + - ]
136 : : // edge flux contributions
137 : 155142 : D[N[0]] = D[N[0]] - f[0] + f[2] - f[3];
138 : 155142 : D[N[1]] = D[N[1]] + f[0] - f[1] - f[4];
139 : 155142 : D[N[2]] = D[N[2]] + f[1] - f[2] - f[5];
140 : 155142 : D[N[3]] = D[N[3]] + f[3] + f[4] + f[5];
141 : : }
142 : :
143 : : // domain edge contributions: triangle superedges
144 [ + + ]: 156406 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
145 : 152078 : const auto N = dsupedge[1].data() + e*3;
146 : 152078 : const auto d = dsupint[1].data();
147 : : // edge fluxes
148 : : tk::real f[] = {
149 [ + - ]: 152078 : div( coord, d+(e*3+0)*5, dt, P, G, U, N[0], N[1], stab ),
150 : 304156 : div( coord, d+(e*3+1)*5, dt, P, G, U, N[1], N[2], stab ),
151 [ + - ][ + - ]: 152078 : div( coord, d+(e*3+2)*5, dt, P, G, U, N[2], N[0], stab ) };
152 : : // edge flux contributions
153 : 152078 : D[N[0]] = D[N[0]] - f[0] + f[2];
154 : 152078 : D[N[1]] = D[N[1]] + f[0] - f[1];
155 : 152078 : D[N[2]] = D[N[2]] + f[1] - f[2];
156 : : }
157 : :
158 : : // domain edge contributions: edges
159 [ + + ]: 352617 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
160 : 348289 : const auto N = dsupedge[2].data() + e*2;
161 : 348289 : const auto d = dsupint[2].data();
162 : : // edge flux
163 : 348289 : tk::real f = div( coord, d+e*5, dt, P, G, U, N[0], N[1], stab );
164 : : // edge flux contributions
165 : 348289 : D[N[0]] -= f;
166 : 348289 : D[N[1]] += f;
167 : : }
168 : :
169 : : // boundary integral
170 : :
171 : 4328 : const auto& x = coord[0];
172 : 4328 : const auto& y = coord[1];
173 : 4328 : const auto& z = coord[2];
174 : :
175 [ + + ]: 672208 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
176 : 667880 : const auto N = triinpoel.data() + e*3;
177 : : tk::real n[3];
178 : 667880 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
179 : 667880 : 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 [ + - ]: 667880 : auto uxA = U(N[0],0);
182 [ + - ]: 667880 : auto uyA = U(N[0],1);
183 [ + - ]: 667880 : auto uzA = U(N[0],2);
184 [ + - ]: 667880 : auto uxB = U(N[1],0);
185 [ + - ]: 667880 : auto uyB = U(N[1],1);
186 [ + - ]: 667880 : auto uzB = U(N[1],2);
187 [ + - ]: 667880 : auto uxC = U(N[2],0);
188 [ + - ]: 667880 : auto uyC = U(N[2],1);
189 [ + - ]: 667880 : auto uzC = U(N[2],2);
190 : 667880 : auto ux = (6.0*uxA + uxB + uxC)/8.0;
191 : 667880 : auto uy = (6.0*uyA + uyB + uyC)/8.0;
192 : 667880 : auto uz = (6.0*uzA + uzB + uzC)/8.0;
193 : 667880 : D[N[0]] += ux*n[0] + uy*n[1] + uz*n[2];
194 : 667880 : ux = (uxA + 6.0*uxB + uxC)/8.0;
195 : 667880 : uy = (uyA + 6.0*uyB + uyC)/8.0;
196 : 667880 : uz = (uzA + 6.0*uzB + uzC)/8.0;
197 : 667880 : D[N[1]] += ux*n[0] + uy*n[1] + uz*n[2];
198 : 667880 : ux = (uxA + uxB + 6.0*uxC)/8.0;
199 : 667880 : uy = (uyA + uyB + 6.0*uyC)/8.0;
200 : 667880 : uz = (uzA + uzB + 6.0*uzC)/8.0;
201 : 667880 : 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 : 4328 : }
210 : :
211 : : void
212 : 3413 : 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 [ - + ][ - - ]: 3413 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
229 [ - + ][ - - ]: 3413 : 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 [ + + ]: 87760 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
244 : 84347 : const auto N = dsupedge[0].data() + e*4;
245 : 84347 : const auto d = dsupint[0].data();
246 [ + + ]: 337388 : for (std::size_t i=0; i<3; ++i) {
247 [ + - ][ + - ]: 253041 : tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i), U(N[3],i) };
[ + - ][ + - ]
248 : 253041 : auto i3 = i*3;
249 [ + + ]: 1012164 : for (std::size_t j=0; j<3; ++j) {
250 : 759123 : tk::real f[] = { d[(e*6+0)*5+j] * (u[1] + u[0]),
251 : 759123 : d[(e*6+1)*5+j] * (u[2] + u[1]),
252 : 759123 : d[(e*6+2)*5+j] * (u[0] + u[2]),
253 : 759123 : d[(e*6+3)*5+j] * (u[3] + u[0]),
254 : 759123 : d[(e*6+4)*5+j] * (u[3] + u[1]),
255 : 759123 : d[(e*6+5)*5+j] * (u[3] + u[2]) };
256 [ + - ][ + - ]: 759123 : G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2] - f[3];
257 [ + - ][ + - ]: 759123 : G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1] - f[4];
258 [ + - ][ + - ]: 759123 : G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2] - f[5];
259 [ + - ][ + - ]: 759123 : 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 [ + + ]: 87324 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
266 : 83911 : const auto N = dsupedge[1].data() + e*3;
267 : 83911 : const auto d = dsupint[1].data();
268 [ + + ]: 335644 : for (std::size_t i=0; i<3; ++i) {
269 [ + - ][ + - ]: 251733 : tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i) };
[ + - ]
270 : 251733 : auto i3 = i*3;
271 [ + + ]: 1006932 : for (std::size_t j=0; j<3; ++j) {
272 : 755199 : tk::real f[] = { d[(e*3+0)*5+j] * (u[1] + u[0]),
273 : 755199 : d[(e*3+1)*5+j] * (u[2] + u[1]),
274 : 755199 : d[(e*3+2)*5+j] * (u[0] + u[2]) };
275 [ + - ][ + - ]: 755199 : G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2];
276 [ + - ][ + - ]: 755199 : G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1];
277 [ + - ][ + - ]: 755199 : G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2];
278 : : }
279 : : }
280 : : }
281 : :
282 : : // domain edge contributions: edges
283 [ + + ]: 198604 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
284 : 195191 : const auto N = dsupedge[2].data() + e*2;
285 : 195191 : const auto d = dsupint[2].data() + e*5;
286 [ + + ]: 780764 : for (std::size_t i=0; i<3; ++i) {
287 [ + - ][ + - ]: 585573 : tk::real u[] = { U(N[0],i), U(N[1],i) };
288 : 585573 : auto i3 = i*3;
289 [ + + ]: 2342292 : for (std::size_t j=0; j<3; ++j) {
290 : 1756719 : tk::real f = d[j] * (u[1] + u[0]);
291 [ + - ]: 1756719 : G(N[0],i3+j) -= f;
292 [ + - ]: 1756719 : G(N[1],i3+j) += f;
293 : : }
294 : : }
295 : : }
296 : :
297 : : // boundary integral
298 : :
299 : 3413 : const auto& x = coord[0];
300 : 3413 : const auto& y = coord[1];
301 : 3413 : const auto& z = coord[2];
302 : :
303 [ + + ]: 366345 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
304 : 362932 : const auto N = triinpoel.data() + e*3;
305 : : tk::real n[3];
306 : 362932 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
307 : 362932 : 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 [ + + ]: 1451728 : for (std::size_t i=0; i<3; ++i) {
310 [ + - ][ + - ]: 1088796 : tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i) };
[ + - ]
311 : 1088796 : auto i3 = i*3;
312 : 1088796 : auto f = (6.0*u[0] + u[1] + u[2])/8.0;
313 [ + - ]: 1088796 : G(N[0],i3+0) += f * n[0];
314 [ + - ]: 1088796 : G(N[0],i3+1) += f * n[1];
315 [ + - ]: 1088796 : G(N[0],i3+2) += f * n[2];
316 : 1088796 : f = (u[0] + 6.0*u[1] + u[2])/8.0;
317 [ + - ]: 1088796 : G(N[1],i3+0) += f * n[0];
318 [ + - ]: 1088796 : G(N[1],i3+1) += f * n[1];
319 [ + - ]: 1088796 : G(N[1],i3+2) += f * n[2];
320 : 1088796 : f = (u[0] + u[1] + 6.0*u[2])/8.0;
321 [ + - ]: 1088796 : G(N[2],i3+0) += f * n[0];
322 [ + - ]: 1088796 : G(N[2],i3+1) += f * n[1];
323 [ + - ]: 1088796 : 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 : 3413 : }
333 : :
334 : : void
335 : 7958 : 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 [ - + ][ - - ]: 7958 : Assert( G.nunk() == U.size(), "Size mismatch" );
[ - - ][ - - ]
352 [ - + ][ - - ]: 7958 : 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 [ + + ]: 300070 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
367 : 292112 : const auto N = dsupedge[0].data() + e*4;
368 : 292112 : const auto d = dsupint[0].data();
369 : 292112 : tk::real u[] = { U[N[0]], U[N[1]], U[N[2]], U[N[3]] };
370 [ + + ]: 1168448 : for (std::size_t j=0; j<3; ++j) {
371 : : tk::real f[] = {
372 : 876336 : d[(e*6+0)*5+j] * (u[1] + u[0]),
373 : 876336 : d[(e*6+1)*5+j] * (u[2] + u[1]),
374 : 876336 : d[(e*6+2)*5+j] * (u[0] + u[2]),
375 : 876336 : d[(e*6+3)*5+j] * (u[3] + u[0]),
376 : 876336 : d[(e*6+4)*5+j] * (u[3] + u[1]),
377 : 876336 : d[(e*6+5)*5+j] * (u[3] + u[2]) };
378 [ + - ][ + - ]: 876336 : G(N[0],j) = G(N[0],j) - f[0] + f[2] - f[3];
379 [ + - ][ + - ]: 876336 : G(N[1],j) = G(N[1],j) + f[0] - f[1] - f[4];
380 [ + - ][ + - ]: 876336 : G(N[2],j) = G(N[2],j) + f[1] - f[2] - f[5];
381 [ + - ][ + - ]: 876336 : G(N[3],j) = G(N[3],j) + f[3] + f[4] + f[5];
382 : : }
383 : : }
384 : :
385 : : // domain edge contributions: triangle superedges
386 [ + + ]: 294446 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
387 : 286488 : const auto N = dsupedge[1].data() + e*3;
388 : 286488 : const auto d = dsupint[1].data();
389 : 286488 : tk::real u[] = { U[N[0]], U[N[1]], U[N[2]] };
390 [ + + ]: 1145952 : for (std::size_t j=0; j<3; ++j) {
391 : : tk::real f[] = {
392 : 859464 : d[(e*3+0)*5+j] * (u[1] + u[0]),
393 : 859464 : d[(e*3+1)*5+j] * (u[2] + u[1]),
394 : 859464 : d[(e*3+2)*5+j] * (u[0] + u[2]) };
395 [ + - ][ + - ]: 859464 : G(N[0],j) = G(N[0],j) - f[0] + f[2];
396 [ + - ][ + - ]: 859464 : G(N[1],j) = G(N[1],j) + f[0] - f[1];
397 [ + - ][ + - ]: 859464 : G(N[2],j) = G(N[2],j) + f[1] - f[2];
398 : : }
399 : : }
400 : :
401 : : // domain edge contributions: edges
402 [ + + ]: 660757 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
403 : 652799 : const auto N = dsupedge[2].data() + e*2;
404 : 652799 : const auto d = dsupint[2].data() + e*5;
405 : 652799 : tk::real u[] = { U[N[0]], U[N[1]] };
406 [ + + ]: 2611196 : for (std::size_t j=0; j<3; ++j) {
407 : 1958397 : tk::real f = d[j] * (u[1] + u[0]);
408 [ + - ]: 1958397 : G(N[0],j) -= f;
409 [ + - ]: 1958397 : G(N[1],j) += f;
410 : : }
411 : : }
412 : :
413 : : // boundary integral
414 : :
415 : 7958 : const auto& x = coord[0];
416 : 7958 : const auto& y = coord[1];
417 : 7958 : const auto& z = coord[2];
418 : :
419 [ + + ]: 1270358 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
420 : 1262400 : const auto N = triinpoel.data() + e*3;
421 : : tk::real n[3];
422 : 1262400 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
423 : 1262400 : 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 : 1262400 : auto uA = U[N[0]];
426 : 1262400 : auto uB = U[N[1]];
427 : 1262400 : auto uC = U[N[2]];
428 : 1262400 : auto f = (6.0*uA + uB + uC)/8.0;
429 [ + - ]: 1262400 : G(N[0],0) += f * n[0];
430 [ + - ]: 1262400 : G(N[0],1) += f * n[1];
431 [ + - ]: 1262400 : G(N[0],2) += f * n[2];
432 : 1262400 : f = (uA + 6.0*uB + uC)/8.0;
433 [ + - ]: 1262400 : G(N[1],0) += f * n[0];
434 [ + - ]: 1262400 : G(N[1],1) += f * n[1];
435 [ + - ]: 1262400 : G(N[1],2) += f * n[2];
436 : 1262400 : f = (uA + uB + 6.0*uC)/8.0;
437 [ + - ]: 1262400 : G(N[2],0) += f * n[0];
438 [ + - ]: 1262400 : G(N[2],1) += f * n[1];
439 [ + - ]: 1262400 : 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 : 7958 : }
448 : :
449 : : static tk::real
450 : 884934 : 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 : 884934 : auto inv = U(p,i)*U(p,j) + U(q,i)*U(q,j);
468 : :
469 : 884934 : auto eps = std::numeric_limits< tk::real >::epsilon();
470 : 884934 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
471 [ + + ]: 884934 : if (mu < eps) return -inv;
472 : :
473 : 524826 : auto vis = G(p,i*3+j) + G(p,j*3+i) + G(q,i*3+j) + G(q,j*3+i);
474 [ + + ]: 524826 : if (i == j) {
475 : 174942 : 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 : 524826 : return mu*vis - inv;
478 : : }
479 : :
480 : : static tk::real
481 : 951804 : 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 : 951804 : auto inv = U(p,i)*U(p,j);
497 : :
498 : 951804 : auto eps = std::numeric_limits< tk::real >::epsilon();
499 : 951804 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
500 [ + + ]: 951804 : if (mu < eps) return -inv;
501 : :
502 : 686448 : auto vis = G(p,i*3+j) + G(p,j*3+i);
503 [ + + ]: 686448 : if (i == j) {
504 : 228816 : vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) );
505 : : }
506 : 686448 : return mu*vis - inv;
507 : : }
508 : :
509 : : void
510 : 333 : 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 [ - + ][ - - ]: 333 : Assert( F.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
530 [ - + ][ - - ]: 333 : Assert( F.nprop() == 3, "Size mismatch" );
[ - - ][ - - ]
531 [ - + ][ - - ]: 333 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
532 [ - + ][ - - ]: 333 : 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 [ + + ]: 9040 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
547 : 8707 : const auto N = dsupedge[0].data() + e*4;
548 : 8707 : const auto d = dsupint[0].data();
549 [ + + ]: 34828 : for (std::size_t i=0; i<3; ++i) {
550 [ + + ]: 104484 : for (std::size_t j=0; j<3; ++j) {
551 [ + - ]: 78363 : tk::real f[] = { d[(e*6+0)*5+j] * flux(U,G,i,j,N[1],N[0]),
552 : 156726 : d[(e*6+1)*5+j] * flux(U,G,i,j,N[2],N[1]),
553 : 156726 : d[(e*6+2)*5+j] * flux(U,G,i,j,N[0],N[2]),
554 : 156726 : d[(e*6+3)*5+j] * flux(U,G,i,j,N[3],N[0]),
555 : 156726 : d[(e*6+4)*5+j] * flux(U,G,i,j,N[3],N[1]),
556 [ + - ][ + - ]: 78363 : d[(e*6+5)*5+j] * flux(U,G,i,j,N[3],N[2]) };
[ + - ][ + - ]
[ + - ]
557 [ + - ][ + - ]: 78363 : F(N[0],i) = F(N[0],i) - f[0] + f[2] - f[3];
558 [ + - ][ + - ]: 78363 : F(N[1],i) = F(N[1],i) + f[0] - f[1] - f[4];
559 [ + - ][ + - ]: 78363 : F(N[2],i) = F(N[2],i) + f[1] - f[2] - f[5];
560 [ + - ][ + - ]: 78363 : 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 [ + + ]: 8794 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
567 : 8461 : const auto N = dsupedge[1].data() + e*3;
568 : 8461 : const auto d = dsupint[1].data();
569 [ + + ]: 33844 : for (std::size_t i=0; i<3; ++i) {
570 [ + + ]: 101532 : for (std::size_t j=0; j<3; ++j) {
571 [ + - ]: 76149 : tk::real f[] = { d[(e*3+0)*5+j] * flux(U,G,i,j,N[1],N[0]),
572 : 152298 : d[(e*3+1)*5+j] * flux(U,G,i,j,N[2],N[1]),
573 [ + - ][ + - ]: 76149 : d[(e*3+2)*5+j] * flux(U,G,i,j,N[0],N[2]), };
574 [ + - ][ + - ]: 76149 : F(N[0],i) = F(N[0],i) - f[0] + f[2];
575 [ + - ][ + - ]: 76149 : F(N[1],i) = F(N[1],i) + f[0] - f[1];
576 [ + - ][ + - ]: 76149 : F(N[2],i) = F(N[2],i) + f[1] - f[2];
577 : : }
578 : : }
579 : : }
580 : :
581 : : // domain edge contributions: edges
582 [ + + ]: 21034 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
583 : 20701 : const auto N = dsupedge[2].data() + e*2;
584 : 20701 : const auto d = dsupint[2].data() + e*5;
585 [ + + ]: 82804 : for (std::size_t i=0; i<3; ++i) {
586 [ + + ]: 248412 : for (std::size_t j=0; j<3; ++j) {
587 : 186309 : tk::real f = d[j] * flux(U,G,i,j,N[1],N[0]);
588 : 186309 : F(N[0],i) -= f;
589 : 186309 : F(N[1],i) += f;
590 : : }
591 : : }
592 : : }
593 : :
594 : : // boundary integral
595 : :
596 : 333 : const auto& x = coord[0];
597 : 333 : const auto& y = coord[1];
598 : 333 : const auto& z = coord[2];
599 : :
600 [ + + ]: 35585 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
601 : 35252 : const auto N = triinpoel.data() + e*3;
602 : : tk::real n[3];
603 : 35252 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
604 : 35252 : 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 [ + + ]: 141008 : for (std::size_t i=0; i<3; ++i) {
607 [ + - ]: 105756 : auto fxA = flux(U,G,i,0,N[0]);
608 [ + - ]: 105756 : auto fyA = flux(U,G,i,1,N[0]);
609 [ + - ]: 105756 : auto fzA = flux(U,G,i,2,N[0]);
610 [ + - ]: 105756 : auto fxB = flux(U,G,i,0,N[1]);
611 [ + - ]: 105756 : auto fyB = flux(U,G,i,1,N[1]);
612 [ + - ]: 105756 : auto fzB = flux(U,G,i,2,N[1]);
613 [ + - ]: 105756 : auto fxC = flux(U,G,i,0,N[2]);
614 [ + - ]: 105756 : auto fyC = flux(U,G,i,1,N[2]);
615 [ + - ]: 105756 : auto fzC = flux(U,G,i,2,N[2]);
616 : 105756 : auto fx = (6.0*fxA + fxB + fxC)/8.0;
617 : 105756 : auto fy = (6.0*fyA + fyB + fyC)/8.0;
618 : 105756 : auto fz = (6.0*fzA + fzB + fzC)/8.0;
619 [ + - ]: 105756 : F(N[0],i) += fx*n[0] + fy*n[1] + fz*n[2];
620 : 105756 : fx = (fxA + 6.0*fxB + fxC)/8.0;
621 : 105756 : fy = (fyA + 6.0*fyB + fyC)/8.0;
622 : 105756 : fz = (fzA + 6.0*fzB + fzC)/8.0;
623 [ + - ]: 105756 : F(N[1],i) += fx*n[0] + fy*n[1] + fz*n[2];
624 : 105756 : fx = (fxA + fxB + 6.0*fxC)/8.0;
625 : 105756 : fy = (fyA + fyB + 6.0*fyC)/8.0;
626 : 105756 : fz = (fzA + fzB + 6.0*fzC)/8.0;
627 [ + - ]: 105756 : 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 : 333 : }
637 : :
638 : : static void
639 : 566040 : 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 : 566040 : const auto ncomp = U.nprop();
662 : 566040 : const auto& x = coord[0];
663 : 566040 : const auto& y = coord[1];
664 : 566040 : const auto& z = coord[2];
665 : :
666 : : // edge vector
667 : 566040 : auto dx = x[p] - x[q];
668 : 566040 : auto dy = y[p] - y[q];
669 : 566040 : auto dz = z[p] - z[q];
670 : 566040 : auto dl = dx*dx + dy*dy + dz*dz;
671 : 566040 : dx /= dl;
672 : 566040 : dy /= dl;
673 : 566040 : dz /= dl;
674 : :
675 : : // left state
676 [ + - ]: 566040 : auto uL = U(p,0);
677 [ + - ]: 566040 : auto vL = U(p,1);
678 [ + - ]: 566040 : auto wL = U(p,2);
679 : 566040 : auto pL = P[p];
680 : 566040 : auto dnL = uL*dx + vL*dy + wL*dz;
681 : :
682 : : // right state
683 [ + - ]: 566040 : auto uR = U(q,0);
684 [ + - ]: 566040 : auto vR = U(q,1);
685 [ + - ]: 566040 : auto wR = U(q,2);
686 : 566040 : auto pR = P[q];
687 : 566040 : auto dnR = uR*dx + vR*dy + wR*dz;
688 : :
689 : 566040 : auto nx = supint[0];
690 : 566040 : auto ny = supint[1];
691 : 566040 : 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 : 566040 : tk::real ue[ncomp];
705 : :
706 : : // flow
707 : 566040 : auto dp = pL - pR;
708 : 566040 : ue[0] = 0.5*(uL + uR - dt*(uL*dnL - uR*dnR + dp*dx));
709 : 566040 : ue[1] = 0.5*(vL + vR - dt*(vL*dnL - vR*dnR + dp*dy));
710 : 566040 : ue[2] = 0.5*(wL + wR - dt*(wL*dnL - wR*dnR + dp*dz));
711 : :
712 : : // Taylor-Galerkin second half step
713 : :
714 : 566040 : auto uh = ue[0];
715 : 566040 : auto vh = ue[1];
716 : 566040 : auto wh = ue[2];
717 : 566040 : auto ph = (pL + pR)/2.0;
718 : 566040 : auto vn = uh*nx + vh*ny + wh*nz;
719 : :
720 : : // viscosity
721 : 566040 : auto d = supint[4] * g_cfg.get< tag::mat_dyn_viscosity >();
722 : :
723 : : // flow
724 : 566040 : f[0] = 2.0*(uh*vn + ph*nx) - d*(uR - uL);
725 : 566040 : f[1] = 2.0*(vh*vn + ph*ny) - d*(vR - vL);
726 : 566040 : f[2] = 2.0*(wh*vn + ph*nz) - d*(wR - wL);
727 : :
728 : : // artificial viscosity
729 : :
730 : 566040 : const auto stab2 = g_cfg.get< tag::stab2 >();
731 [ + - ]: 566040 : 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 : 566040 : }
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 : : // viscosity
793 : 759360 : auto d = supint[4] * g_cfg.get< tag::mat_dyn_viscosity >();
794 : :
795 : : // flow
796 : 759360 : auto pf = P[p] + P[q];
797 : 759360 : f[0] = uL*vnL + uR*vnR + pf*nx + (aw-d)*(uR-uL);
798 : 759360 : f[1] = vL*vnL + vR*vnR + pf*ny + (aw-d)*(vR-vL);
799 : 759360 : f[2] = wL*vnL + wR*vnR + pf*nz + (aw-d)*(wR-wL);
800 : 759360 : }
801 : :
802 : : static void
803 : 854680 : adv_damp4( const tk::real supint[],
804 : : const tk::Fields& U,
805 : : const tk::Fields& G,
806 : : const std::vector< tk::real >& P,
807 : : const tk::Fields& W,
808 : : const std::array< std::vector< tk::real >, 3 >& coord,
809 : : tk::real,
810 : : std::size_t p,
811 : : std::size_t q,
812 : : tk::real f[] )
813 : : // *****************************************************************************
814 : : //! Compute advection fluxes on a single edge using 4th-order damping
815 : : //! \param[in] supint Edge integral
816 : : //! \param[in] U Velocity and transported scalars at recent time step
817 : : //! \param[in] G Gradients of velocity and transported scalars
818 : : //! \param[in] P Pressure
819 : : //! \param[in] W Pressure gradient
820 : : //! \param[in] coord Mesh node coordinates
821 : : //! \param[in] p Left node index of edge
822 : : //! \param[in] q Right node index of edge
823 : : //! \param[in,out] f Flux computed
824 : : // *****************************************************************************
825 : : {
826 : 854680 : const auto& x = coord[0];
827 : 854680 : const auto& y = coord[1];
828 : 854680 : const auto& z = coord[2];
829 : :
830 : : // edge vector
831 : 854680 : auto dx = x[p] - x[q];
832 : 854680 : auto dy = y[p] - y[q];
833 : 854680 : auto dz = z[p] - z[q];
834 : :
835 : : #if defined(__clang__)
836 : : #pragma clang diagnostic push
837 : : #pragma clang diagnostic ignored "-Wvla"
838 : : #pragma clang diagnostic ignored "-Wvla-extension"
839 : : #elif defined(STRICT_GNUC)
840 : : #pragma GCC diagnostic push
841 : : #pragma GCC diagnostic ignored "-Wvla"
842 : : #endif
843 : :
844 [ + - ][ + - ]: 854680 : tk::real uL[] = { U(p,0), U(p,1), U(p,2), P[p] };
[ + - ]
845 [ + - ][ + - ]: 854680 : tk::real uR[] = { U(q,0), U(q,1), U(q,2), P[q] };
[ + - ]
846 [ + - ]: 854680 : tk::real gL[] = { G(p,0), G(p,1), G(p,2),
847 : 3418720 : G(p,3), G(p,4), G(p,5),
848 : 3418720 : G(p,6), G(p,7), G(p,8),
849 [ + - ][ + - ]: 854680 : W(p,0), W(p,1), W(p,2) };
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
850 [ + - ]: 854680 : tk::real gR[] = { G(q,0), G(q,1), G(q,2),
851 : 3418720 : G(q,3), G(q,4), G(q,5),
852 : 3418720 : G(q,6), G(q,7), G(q,8),
853 [ + - ][ + - ]: 854680 : W(q,0), W(q,1), W(q,2) };
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
854 : :
855 : : // MUSCL reconstruction in edge-end points
856 [ + + ]: 4273400 : for (std::size_t c=0; c<4; ++c) {
857 : 3418720 : auto g = c*3;
858 : 3418720 : auto g1 = gL[g+0]*dx + gL[g+1]*dy + gL[g+2]*dz;
859 : 3418720 : auto g2 = gR[g+0]*dx + gR[g+1]*dy + gR[g+2]*dz;
860 : 3418720 : auto delta2 = uR[c] - uL[c];
861 : 3418720 : auto delta1 = 2.0 * g1 - delta2;
862 : 3418720 : auto delta3 = 2.0 * g2 - delta2;
863 : :
864 : : // van Leer limiter
865 : 3418720 : auto rL = (delta2 + muscl_eps) / (delta1 + muscl_eps);
866 : 3418720 : auto rR = (delta2 + muscl_eps) / (delta3 + muscl_eps);
867 : 3418720 : auto rLinv = (delta1 + muscl_eps) / (delta2 + muscl_eps);
868 : 3418720 : auto rRinv = (delta3 + muscl_eps) / (delta2 + muscl_eps);
869 : 3418720 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
870 : 3418720 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
871 : 3418720 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
872 : 3418720 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
873 : : // update unknowns with reconstructed unknowns
874 : 3418720 : uL[c] += 0.25*(delta1*(1.0-muscl_const)*phiL +
875 : 3418720 : delta2*(1.0+muscl_const)*phi_L_inv);
876 : 3418720 : uR[c] -= 0.25*(delta3*(1.0-muscl_const)*phiR +
877 : 3418720 : delta2*(1.0+muscl_const)*phi_R_inv);
878 : : }
879 : :
880 : 854680 : auto nx = supint[0];
881 : 854680 : auto ny = supint[1];
882 : 854680 : auto nz = supint[2];
883 : :
884 : : // normal velocities
885 : 854680 : auto vnL = uL[0]*nx + uL[1]*ny + uL[2]*nz;
886 : 854680 : auto vnR = uR[0]*nx + uR[1]*ny + uR[2]*nz;
887 : :
888 : : // stabilization
889 : 854680 : auto aw = std::abs( vnL + vnR ) / 2.0 * tk::length( nx, ny, nz );
890 : :
891 : : // viscosity
892 : 854680 : auto d = supint[4] * g_cfg.get< tag::mat_dyn_viscosity >();
893 : :
894 : : // flow
895 : 854680 : auto pf = uL[3] + uR[3];
896 [ + - ][ + - ]: 854680 : f[0] = uL[0]*vnL + uR[0]*vnR + pf*nx + aw*(uR[0]-uL[0]) - d*(U(q,0)-U(p,0));
897 [ + - ][ + - ]: 854680 : f[1] = uL[1]*vnL + uR[1]*vnR + pf*ny + aw*(uR[1]-uL[1]) - d*(U(q,1)-U(p,1));
898 [ + - ][ + - ]: 854680 : f[2] = uL[2]*vnL + uR[2]*vnR + pf*nz + aw*(uR[2]-uL[2]) - d*(U(q,2)-U(p,2));
899 : :
900 : : #if defined(__clang__)
901 : : #pragma clang diagnostic pop
902 : : #elif defined(STRICT_GNUC)
903 : : #pragma GCC diagnostic pop
904 : : #endif
905 : 854680 : }
906 : :
907 : : static void
908 : 3750 : adv( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
909 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
910 : : const std::array< std::vector< tk::real >, 3 >& coord,
911 : : const std::vector< std::size_t >& triinpoel,
912 : : tk::real dt,
913 : : const tk::Fields& U,
914 : : const tk::Fields& G,
915 : : const std::vector< tk::real >& P,
916 : : const tk::Fields& W,
917 : : // cppcheck-suppress constParameter
918 : : tk::Fields& R )
919 : : // *****************************************************************************
920 : : //! Add advection to rhs
921 : : //! \param[in] dsupedge Domain superedges
922 : : //! \param[in] dsupint Domain superedge integrals
923 : : //! \param[in] coord Mesh node coordinates
924 : : //! \param[in] triinpoel Boundary face connectivity
925 : : //! \param[in] dt Physical time step size
926 : : //! \param[in] U Velocity and transported scalars at recent time step
927 : : //! \param[in] G Gradients of velocity and transported scalars
928 : : //! \param[in] P Pressure
929 : : //! \param[in] W Pressure gradient
930 : : //! \param[in,out] R Right-hand side vector added to
931 : : // *****************************************************************************
932 : : {
933 : : // configure advection stabilization
934 : 0 : auto adv = [](){
935 : 3750 : const auto& flux = g_cfg.get< tag::flux >();
936 [ + + ]: 3750 : if (flux == "tg") return adv_tg;
937 [ + + ]: 3240 : else if (flux == "damp2") return adv_damp2;
938 [ + - ]: 3080 : else if (flux == "damp4") return adv_damp4;
939 [ - - ][ - - ]: 0 : else Throw( "Flux not correctly configured" );
[ - - ]
940 [ + - ]: 3750 : }();
941 : :
942 : 3750 : auto ncomp = U.nprop();
943 : :
944 : : #if defined(__clang__)
945 : : #pragma clang diagnostic push
946 : : #pragma clang diagnostic ignored "-Wvla"
947 : : #pragma clang diagnostic ignored "-Wvla-extension"
948 : : #elif defined(STRICT_GNUC)
949 : : #pragma GCC diagnostic push
950 : : #pragma GCC diagnostic ignored "-Wvla"
951 : : #endif
952 : :
953 : : // domain integral
954 : :
955 : : // domain edge contributions: tetrahedron superedges
956 [ + + ]: 200600 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
957 : 196850 : const auto N = dsupedge[0].data() + e*4;
958 : 196850 : tk::real f[6][ncomp];
959 : 196850 : const auto d = dsupint[0].data();
960 [ + - ]: 196850 : adv( d+(e*6+0)*5, U, G, P, W, coord, dt, N[0], N[1], f[0] );
961 [ + - ]: 196850 : adv( d+(e*6+1)*5, U, G, P, W, coord, dt, N[1], N[2], f[1] );
962 [ + - ]: 196850 : adv( d+(e*6+2)*5, U, G, P, W, coord, dt, N[2], N[0], f[2] );
963 [ + - ]: 196850 : adv( d+(e*6+3)*5, U, G, P, W, coord, dt, N[0], N[3], f[3] );
964 [ + - ]: 196850 : adv( d+(e*6+4)*5, U, G, P, W, coord, dt, N[1], N[3], f[4] );
965 [ + - ]: 196850 : adv( d+(e*6+5)*5, U, G, P, W, coord, dt, N[2], N[3], f[5] );
966 [ + + ]: 787400 : for (std::size_t c=0; c<ncomp; ++c) {
967 [ + - ][ + - ]: 590550 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
968 [ + - ][ + - ]: 590550 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
969 [ + - ][ + - ]: 590550 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
970 [ + - ][ + - ]: 590550 : R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
971 : : }
972 : 196850 : }
973 : :
974 : : // domain edge contributions: triangle superedges
975 [ + + ]: 197560 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
976 : 193810 : const auto N = dsupedge[1].data() + e*3;
977 : 193810 : tk::real f[3][ncomp];
978 : 193810 : const auto d = dsupint[1].data();
979 [ + - ]: 193810 : adv( d+(e*3+0)*5, U, G, P, W, coord, dt, N[0], N[1], f[0] );
980 [ + - ]: 193810 : adv( d+(e*3+1)*5, U, G, P, W, coord, dt, N[1], N[2], f[1] );
981 [ + - ]: 193810 : adv( d+(e*3+2)*5, U, G, P, W, coord, dt, N[2], N[0], f[2] );
982 [ + + ]: 775240 : for (std::size_t c=0; c<ncomp; ++c) {
983 [ + - ][ + - ]: 581430 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
984 [ + - ][ + - ]: 581430 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
985 [ + - ][ + - ]: 581430 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
986 : : }
987 : 193810 : }
988 : :
989 : : // domain edge contributions: edges
990 [ + + ]: 421300 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
991 : 417550 : const auto N = dsupedge[2].data() + e*2;
992 : 417550 : tk::real f[ncomp];
993 : 417550 : const auto d = dsupint[2].data();
994 [ + - ]: 417550 : adv( d+e*5, U, G, P, W, coord, dt, N[0], N[1], f );
995 [ + + ]: 1670200 : for (std::size_t c=0; c<ncomp; ++c) {
996 [ + - ]: 1252650 : R(N[0],c) -= f[c];
997 [ + - ]: 1252650 : R(N[1],c) += f[c];
998 : : }
999 : 417550 : }
1000 : :
1001 : : // boundary integral
1002 : :
1003 : 3750 : const auto& x = coord[0];
1004 : 3750 : const auto& y = coord[1];
1005 : 3750 : const auto& z = coord[2];
1006 : :
1007 [ + + ]: 888670 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
1008 : 1769840 : const auto N = triinpoel.data() + e*3;
1009 : : tk::real n[3];
1010 : 884920 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
1011 : 884920 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
1012 : : n[0], n[1], n[2] );
1013 : 884920 : tk::real f[ncomp][3];
1014 : :
1015 [ + - ]: 884920 : auto u = U(N[0],0);
1016 [ + - ]: 884920 : auto v = U(N[0],1);
1017 [ + - ]: 884920 : auto w = U(N[0],2);
1018 : 884920 : auto p = P[N[0]];
1019 : 884920 : auto vn = n[0]*u + n[1]*v + n[2]*w;
1020 : 884920 : f[0][0] = u*vn + p*n[0];
1021 : 884920 : f[1][0] = v*vn + p*n[1];
1022 : 884920 : f[2][0] = w*vn + p*n[2];
1023 [ - - ][ - + ]: 884920 : for (std::size_t c=3; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;
1024 : :
1025 [ + - ]: 884920 : u = U(N[1],0);
1026 [ + - ]: 884920 : v = U(N[1],1);
1027 [ + - ]: 884920 : w = U(N[1],2);
1028 : 884920 : p = P[N[1]];
1029 : 884920 : vn = n[0]*u + n[1]*v + n[2]*w;
1030 : 884920 : f[0][1] = u*vn + p*n[0];
1031 : 884920 : f[1][1] = v*vn + p*n[1];
1032 : 884920 : f[2][1] = w*vn + p*n[2];
1033 [ - - ][ - + ]: 884920 : for (std::size_t c=3; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;
1034 : :
1035 [ + - ]: 884920 : u = U(N[2],0);
1036 [ + - ]: 884920 : v = U(N[2],1);
1037 [ + - ]: 884920 : w = U(N[2],2);
1038 : 884920 : p = P[N[2]];
1039 : 884920 : vn = n[0]*u + n[1]*v + n[2]*w;
1040 : 884920 : f[0][2] = u*vn + p*n[0];
1041 : 884920 : f[1][2] = v*vn + p*n[1];
1042 : 884920 : f[2][2] = w*vn + p*n[2];
1043 [ - - ][ - + ]: 884920 : for (std::size_t c=3; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;
1044 : :
1045 [ + + ]: 3539680 : for (std::size_t c=0; c<ncomp; ++c) {
1046 [ + - ]: 2654760 : R(N[0],c) += (6.0*f[c][0] + f[c][1] + f[c][2])/8.0;
1047 [ + - ]: 2654760 : R(N[1],c) += (f[c][0] + 6.0*f[c][1] + f[c][2])/8.0;
1048 [ + - ]: 2654760 : R(N[2],c) += (f[c][0] + f[c][1] + 6.0*f[c][2])/8.0;
1049 : : }
1050 : 884920 : }
1051 : :
1052 : : #if defined(__clang__)
1053 : : #pragma clang diagnostic pop
1054 : : #elif defined(STRICT_GNUC)
1055 : : #pragma GCC diagnostic pop
1056 : : #endif
1057 : 3750 : }
1058 : :
1059 : : void
1060 : 3750 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
1061 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
1062 : : const std::array< std::vector< tk::real >, 3 >& coord,
1063 : : const std::vector< std::size_t >& triinpoel,
1064 : : tk::real dt,
1065 : : const std::vector< tk::real >& P,
1066 : : const tk::Fields& U,
1067 : : const tk::Fields& G,
1068 : : const tk::Fields& W,
1069 : : tk::Fields& R )
1070 : : // *****************************************************************************
1071 : : // Compute right hand side
1072 : : //! \param[in] dsupedge Domain superedges
1073 : : //! \param[in] dsupint Domain superedge integrals
1074 : : //! \param[in] coord Mesh node coordinates
1075 : : //! \param[in] triinpoel Boundary face connectivity
1076 : : //! \param[in] dt Physical time step size
1077 : : //! \param[in] P Pressure
1078 : : //! \param[in] U Solution vector of primitive variables at recent time step
1079 : : //! \param[in] G Gradients of velocity and transported scalars
1080 : : //! \param[in] W Pressure gradient
1081 : : //! \param[in,out] R Right-hand side vector computed
1082 : : // *****************************************************************************
1083 : : {
1084 [ - + ][ - - ]: 3750 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
1085 : : "vector at recent time step incorrect" );
1086 [ - + ][ - - ]: 3750 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
1087 : : "Number of unknowns and/or number of components in right-hand "
1088 : : "side vector incorrect" );
1089 : :
1090 : 3750 : R.fill( 0.0 );
1091 : 3750 : adv( dsupedge, dsupint, coord, triinpoel, dt, U, G, P, W, R );
1092 : 3750 : }
1093 : :
1094 : : } // chorin::
|
