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 : 2577343 : 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 : 2577343 : tk::real div = d[0] * (U(p,0) + U(q,0)) +
59 : 2577343 : d[1] * (U(p,1) + U(q,1)) +
60 : 2577343 : d[2] * (U(p,2) + U(q,2));
61 : :
62 [ + + ]: 2577343 : if (!stab) return div;
63 : :
64 : 2295070 : const auto& x = coord[0];
65 : 2295070 : const auto& y = coord[1];
66 : 2295070 : const auto& z = coord[2];
67 : :
68 : 2295070 : auto dx = x[p] - x[q];
69 : 2295070 : auto dy = y[p] - y[q];
70 : 2295070 : auto dz = z[p] - z[q];
71 : 2295070 : auto dl = std::sqrt( dx*dx + dy*dy + dz*dz );
72 : 2295070 : auto p2 = P[q] - P[p];
73 : 2295070 : auto D = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] );
74 : :
75 : 2295070 : auto dpx = G(p,0) + G(q,0);
76 : 2295070 : auto dpy = G(p,1) + G(q,1);
77 : 2295070 : auto dpz = G(p,2) + G(q,2);
78 : 2295070 : auto p4 = 0.5 * (dx*dpx + dy*dpy + dz*dpz);
79 : :
80 : 2295070 : div += D*dt/dl*(p2 + p4);
81 : :
82 : 2295070 : return div;
83 : : }
84 : :
85 : : void
86 : 5472 : 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 [ - + ][ - - ]: 5472 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
111 [ - + ][ - - ]: 5472 : 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 [ + + ]: 232348 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
126 : 226876 : const auto N = dsupedge[0].data() + e*4;
127 : 226876 : const auto d = dsupint[0].data();
128 : : // edge fluxes
129 : : tk::real f[] = {
130 [ + - ]: 226876 : div( coord, d+(e*6+0)*5, dt, P, G, U, N[0], N[1], stab ),
131 : 453752 : div( coord, d+(e*6+1)*5, dt, P, G, U, N[1], N[2], stab ),
132 : 453752 : div( coord, d+(e*6+2)*5, dt, P, G, U, N[2], N[0], stab ),
133 : 453752 : div( coord, d+(e*6+3)*5, dt, P, G, U, N[0], N[3], stab ),
134 : 453752 : div( coord, d+(e*6+4)*5, dt, P, G, U, N[1], N[3], stab ),
135 [ + - ][ + - ]: 226876 : div( coord, d+(e*6+5)*5, dt, P, G, U, N[2], N[3], stab ) };
[ + - ][ + - ]
[ + - ]
136 : : // edge flux contributions
137 : 226876 : D[N[0]] = D[N[0]] - f[0] + f[2] - f[3];
138 : 226876 : D[N[1]] = D[N[1]] + f[0] - f[1] - f[4];
139 : 226876 : D[N[2]] = D[N[2]] + f[1] - f[2] - f[5];
140 : 226876 : D[N[3]] = D[N[3]] + f[3] + f[4] + f[5];
141 : : }
142 : :
143 : : // domain edge contributions: triangle superedges
144 [ + + ]: 219897 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
145 : 214425 : const auto N = dsupedge[1].data() + e*3;
146 : 214425 : const auto d = dsupint[1].data();
147 : : // edge fluxes
148 : : tk::real f[] = {
149 [ + - ]: 214425 : div( coord, d+(e*3+0)*5, dt, P, G, U, N[0], N[1], stab ),
150 : 428850 : div( coord, d+(e*3+1)*5, dt, P, G, U, N[1], N[2], stab ),
151 [ + - ][ + - ]: 214425 : div( coord, d+(e*3+2)*5, dt, P, G, U, N[2], N[0], stab ) };
152 : : // edge flux contributions
153 : 214425 : D[N[0]] = D[N[0]] - f[0] + f[2];
154 : 214425 : D[N[1]] = D[N[1]] + f[0] - f[1];
155 : 214425 : D[N[2]] = D[N[2]] + f[1] - f[2];
156 : : }
157 : :
158 : : // domain edge contributions: edges
159 [ + + ]: 578284 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
160 : 572812 : const auto N = dsupedge[2].data() + e*2;
161 : 572812 : const auto d = dsupint[2].data();
162 : : // edge flux
163 : 572812 : tk::real f = div( coord, d+e*5, dt, P, G, U, N[0], N[1], stab );
164 : : // edge flux contributions
165 : 572812 : D[N[0]] -= f;
166 : 572812 : D[N[1]] += f;
167 : : }
168 : :
169 : : // boundary integral
170 : :
171 : 5472 : const auto& x = coord[0];
172 : 5472 : const auto& y = coord[1];
173 : 5472 : const auto& z = coord[2];
174 : :
175 [ + + ]: 938540 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
176 : 933068 : const auto N = triinpoel.data() + e*3;
177 : : tk::real n[3];
178 : 933068 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
179 : 933068 : 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 [ + - ]: 933068 : auto uxA = U(N[0],0);
182 [ + - ]: 933068 : auto uyA = U(N[0],1);
183 [ + - ]: 933068 : auto uzA = U(N[0],2);
184 [ + - ]: 933068 : auto uxB = U(N[1],0);
185 [ + - ]: 933068 : auto uyB = U(N[1],1);
186 [ + - ]: 933068 : auto uzB = U(N[1],2);
187 [ + - ]: 933068 : auto uxC = U(N[2],0);
188 [ + - ]: 933068 : auto uyC = U(N[2],1);
189 [ + - ]: 933068 : auto uzC = U(N[2],2);
190 : 933068 : auto ux = (6.0*uxA + uxB + uxC)/8.0;
191 : 933068 : auto uy = (6.0*uyA + uyB + uyC)/8.0;
192 : 933068 : auto uz = (6.0*uzA + uzB + uzC)/8.0;
193 : 933068 : D[N[0]] += ux*n[0] + uy*n[1] + uz*n[2];
194 : 933068 : ux = (uxA + 6.0*uxB + uxC)/8.0;
195 : 933068 : uy = (uyA + 6.0*uyB + uyC)/8.0;
196 : 933068 : uz = (uzA + 6.0*uzB + uzC)/8.0;
197 : 933068 : D[N[1]] += ux*n[0] + uy*n[1] + uz*n[2];
198 : 933068 : ux = (uxA + uxB + 6.0*uxC)/8.0;
199 : 933068 : uy = (uyA + uyB + 6.0*uyC)/8.0;
200 : 933068 : uz = (uzA + uzB + 6.0*uzC)/8.0;
201 : 933068 : 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 : 5472 : }
210 : :
211 : : void
212 : 3705 : 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+scalar 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 gradients in all points
226 : : // *****************************************************************************
227 : : {
228 [ - + ][ - - ]: 3705 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
229 [ - + ][ - - ]: 3705 : Assert( G.nprop() == U.nprop()*3, "Size mismatch" );
[ - - ][ - - ]
230 : :
231 : 3705 : auto ncomp = U.nprop();
232 : :
233 : : #if defined(__clang__)
234 : : #pragma clang diagnostic push
235 : : #pragma clang diagnostic ignored "-Wvla"
236 : : #pragma clang diagnostic ignored "-Wvla-extension"
237 : : #elif defined(STRICT_GNUC)
238 : : #pragma GCC diagnostic push
239 : : #pragma GCC diagnostic ignored "-Wvla"
240 : : #endif
241 : :
242 : : // domain integral
243 : :
244 : : // domain edge contributions: tetrahedron superedges
245 [ + + ]: 216822 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
246 : 213117 : const auto N = dsupedge[0].data() + e*4;
247 : 213117 : const auto d = dsupint[0].data();
248 [ + + ]: 936632 : for (std::size_t i=0; i<ncomp; ++i) {
249 [ + - ][ + - ]: 723515 : tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i), U(N[3],i) };
[ + - ][ + - ]
250 : 723515 : auto i3 = i*3;
251 [ + + ]: 2894060 : for (std::size_t j=0; j<3; ++j) {
252 : 2170545 : tk::real f[] = { d[(e*6+0)*5+j] * (u[1] + u[0]),
253 : 2170545 : d[(e*6+1)*5+j] * (u[2] + u[1]),
254 : 2170545 : d[(e*6+2)*5+j] * (u[0] + u[2]),
255 : 2170545 : d[(e*6+3)*5+j] * (u[3] + u[0]),
256 : 2170545 : d[(e*6+4)*5+j] * (u[3] + u[1]),
257 : 2170545 : d[(e*6+5)*5+j] * (u[3] + u[2]) };
258 [ + - ][ + - ]: 2170545 : G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2] - f[3];
259 [ + - ][ + - ]: 2170545 : G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1] - f[4];
260 [ + - ][ + - ]: 2170545 : G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2] - f[5];
261 [ + - ][ + - ]: 2170545 : G(N[3],i3+j) = G(N[3],i3+j) + f[3] + f[4] + f[5];
262 : : }
263 : : }
264 : : }
265 : :
266 : : // domain edge contributions: triangle superedges
267 [ + + ]: 197334 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
268 : 193629 : const auto N = dsupedge[1].data() + e*3;
269 : 193629 : const auto d = dsupint[1].data();
270 [ + + ]: 852656 : for (std::size_t i=0; i<ncomp; ++i) {
271 [ + - ][ + - ]: 659027 : tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i) };
[ + - ]
272 : 659027 : auto i3 = i*3;
273 [ + + ]: 2636108 : for (std::size_t j=0; j<3; ++j) {
274 : 1977081 : tk::real f[] = { d[(e*3+0)*5+j] * (u[1] + u[0]),
275 : 1977081 : d[(e*3+1)*5+j] * (u[2] + u[1]),
276 : 1977081 : d[(e*3+2)*5+j] * (u[0] + u[2]) };
277 [ + - ][ + - ]: 1977081 : G(N[0],i3+j) = G(N[0],i3+j) - f[0] + f[2];
278 [ + - ][ + - ]: 1977081 : G(N[1],i3+j) = G(N[1],i3+j) + f[0] - f[1];
279 [ + - ][ + - ]: 1977081 : G(N[2],i3+j) = G(N[2],i3+j) + f[1] - f[2];
280 : : }
281 : : }
282 : : }
283 : :
284 : : // domain edge contributions: edges
285 [ + + ]: 572315 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
286 : 568610 : const auto N = dsupedge[2].data() + e*2;
287 : 568610 : const auto d = dsupint[2].data() + e*5;
288 [ + + ]: 2539617 : for (std::size_t i=0; i<ncomp; ++i) {
289 [ + - ][ + - ]: 1971007 : tk::real u[] = { U(N[0],i), U(N[1],i) };
290 : 1971007 : auto i3 = i*3;
291 [ + + ]: 7884028 : for (std::size_t j=0; j<3; ++j) {
292 : 5913021 : tk::real f = d[j] * (u[1] + u[0]);
293 [ + - ]: 5913021 : G(N[0],i3+j) -= f;
294 [ + - ]: 5913021 : G(N[1],i3+j) += f;
295 : : }
296 : : }
297 : : }
298 : :
299 : : // boundary integral
300 : :
301 : 3705 : const auto& x = coord[0];
302 : 3705 : const auto& y = coord[1];
303 : 3705 : const auto& z = coord[2];
304 : :
305 [ + + ]: 861731 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
306 : 858026 : const auto N = triinpoel.data() + e*3;
307 : : tk::real n[3];
308 : 858026 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
309 : 858026 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
310 : : n[0], n[1], n[2] );
311 [ + + ]: 3829016 : for (std::size_t i=0; i<ncomp; ++i) {
312 [ + - ][ + - ]: 2970990 : tk::real u[] = { U(N[0],i), U(N[1],i), U(N[2],i) };
[ + - ]
313 : 2970990 : auto i3 = i*3;
314 : 2970990 : auto f = (6.0*u[0] + u[1] + u[2])/8.0;
315 [ + - ]: 2970990 : G(N[0],i3+0) += f * n[0];
316 [ + - ]: 2970990 : G(N[0],i3+1) += f * n[1];
317 [ + - ]: 2970990 : G(N[0],i3+2) += f * n[2];
318 : 2970990 : f = (u[0] + 6.0*u[1] + u[2])/8.0;
319 [ + - ]: 2970990 : G(N[1],i3+0) += f * n[0];
320 [ + - ]: 2970990 : G(N[1],i3+1) += f * n[1];
321 [ + - ]: 2970990 : G(N[1],i3+2) += f * n[2];
322 : 2970990 : f = (u[0] + u[1] + 6.0*u[2])/8.0;
323 [ + - ]: 2970990 : G(N[2],i3+0) += f * n[0];
324 [ + - ]: 2970990 : G(N[2],i3+1) += f * n[1];
325 [ + - ]: 2970990 : G(N[2],i3+2) += f * n[2];
326 : : }
327 : : }
328 : :
329 : : #if defined(__clang__)
330 : : #pragma clang diagnostic pop
331 : : #elif defined(STRICT_GNUC)
332 : : #pragma GCC diagnostic pop
333 : : #endif
334 : 3705 : }
335 : :
336 : : void
337 : 10142 : grad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
338 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
339 : : const std::array< std::vector< tk::real >, 3 >& coord,
340 : : const std::vector< std::size_t >& triinpoel,
341 : : const std::vector< tk::real >& U,
342 : : tk::Fields& G )
343 : : // *****************************************************************************
344 : : // Compute gradients of scalar in all points
345 : : //! \param[in] dsupedge Domain superedges
346 : : //! \param[in] dsupint Domain superedge integrals
347 : : //! \param[in] coord Mesh node coordinates
348 : : //! \param[in] triinpoel Boundary face connectivity
349 : : //! \param[in] U Scalar whose gradient to compute
350 : : //! \param[in,out] G Nodal gradient of scalar in all points
351 : : // *****************************************************************************
352 : : {
353 [ - + ][ - - ]: 10142 : Assert( G.nunk() == U.size(), "Size mismatch" );
[ - - ][ - - ]
354 [ - + ][ - - ]: 10142 : Assert( G.nprop() > 2, "Size mismatch" );
[ - - ][ - - ]
355 : :
356 : : #if defined(__clang__)
357 : : #pragma clang diagnostic push
358 : : #pragma clang diagnostic ignored "-Wvla"
359 : : #pragma clang diagnostic ignored "-Wvla-extension"
360 : : #elif defined(STRICT_GNUC)
361 : : #pragma GCC diagnostic push
362 : : #pragma GCC diagnostic ignored "-Wvla"
363 : : #endif
364 : :
365 : : // domain integral
366 : :
367 : : // domain edge contributions: tetrahedron superedges
368 [ + + ]: 439198 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
369 : 429056 : const auto N = dsupedge[0].data() + e*4;
370 : 429056 : const auto d = dsupint[0].data();
371 : 429056 : tk::real u[] = { U[N[0]], U[N[1]], U[N[2]], U[N[3]] };
372 [ + + ]: 1716224 : for (std::size_t j=0; j<3; ++j) {
373 : : tk::real f[] = {
374 : 1287168 : d[(e*6+0)*5+j] * (u[1] + u[0]),
375 : 1287168 : d[(e*6+1)*5+j] * (u[2] + u[1]),
376 : 1287168 : d[(e*6+2)*5+j] * (u[0] + u[2]),
377 : 1287168 : d[(e*6+3)*5+j] * (u[3] + u[0]),
378 : 1287168 : d[(e*6+4)*5+j] * (u[3] + u[1]),
379 : 1287168 : d[(e*6+5)*5+j] * (u[3] + u[2]) };
380 [ + - ][ + - ]: 1287168 : G(N[0],j) = G(N[0],j) - f[0] + f[2] - f[3];
381 [ + - ][ + - ]: 1287168 : G(N[1],j) = G(N[1],j) + f[0] - f[1] - f[4];
382 [ + - ][ + - ]: 1287168 : G(N[2],j) = G(N[2],j) + f[1] - f[2] - f[5];
383 [ + - ][ + - ]: 1287168 : G(N[3],j) = G(N[3],j) + f[3] + f[4] + f[5];
384 : : }
385 : : }
386 : :
387 : : // domain edge contributions: triangle superedges
388 [ + + ]: 415677 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
389 : 405535 : const auto N = dsupedge[1].data() + e*3;
390 : 405535 : const auto d = dsupint[1].data();
391 : 405535 : tk::real u[] = { U[N[0]], U[N[1]], U[N[2]] };
392 [ + + ]: 1622140 : for (std::size_t j=0; j<3; ++j) {
393 : : tk::real f[] = {
394 : 1216605 : d[(e*3+0)*5+j] * (u[1] + u[0]),
395 : 1216605 : d[(e*3+1)*5+j] * (u[2] + u[1]),
396 : 1216605 : d[(e*3+2)*5+j] * (u[0] + u[2]) };
397 [ + - ][ + - ]: 1216605 : G(N[0],j) = G(N[0],j) - f[0] + f[2];
398 [ + - ][ + - ]: 1216605 : G(N[1],j) = G(N[1],j) + f[0] - f[1];
399 [ + - ][ + - ]: 1216605 : G(N[2],j) = G(N[2],j) + f[1] - f[2];
400 : : }
401 : : }
402 : :
403 : : // domain edge contributions: edges
404 [ + + ]: 1091614 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
405 : 1081472 : const auto N = dsupedge[2].data() + e*2;
406 : 1081472 : const auto d = dsupint[2].data() + e*5;
407 : 1081472 : tk::real u[] = { U[N[0]], U[N[1]] };
408 [ + + ]: 4325888 : for (std::size_t j=0; j<3; ++j) {
409 : 3244416 : tk::real f = d[j] * (u[1] + u[0]);
410 [ + - ]: 3244416 : G(N[0],j) -= f;
411 [ + - ]: 3244416 : G(N[1],j) += f;
412 : : }
413 : : }
414 : :
415 : : // boundary integral
416 : :
417 : 10142 : const auto& x = coord[0];
418 : 10142 : const auto& y = coord[1];
419 : 10142 : const auto& z = coord[2];
420 : :
421 [ + + ]: 1778810 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
422 : 1768668 : const auto N = triinpoel.data() + e*3;
423 : : tk::real n[3];
424 : 1768668 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
425 : 1768668 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
426 : : n[0], n[1], n[2] );
427 : 1768668 : auto uA = U[N[0]];
428 : 1768668 : auto uB = U[N[1]];
429 : 1768668 : auto uC = U[N[2]];
430 : 1768668 : auto f = (6.0*uA + uB + uC)/8.0;
431 [ + - ]: 1768668 : G(N[0],0) += f * n[0];
432 [ + - ]: 1768668 : G(N[0],1) += f * n[1];
433 [ + - ]: 1768668 : G(N[0],2) += f * n[2];
434 : 1768668 : f = (uA + 6.0*uB + uC)/8.0;
435 [ + - ]: 1768668 : G(N[1],0) += f * n[0];
436 [ + - ]: 1768668 : G(N[1],1) += f * n[1];
437 [ + - ]: 1768668 : G(N[1],2) += f * n[2];
438 : 1768668 : f = (uA + uB + 6.0*uC)/8.0;
439 [ + - ]: 1768668 : G(N[2],0) += f * n[0];
440 [ + - ]: 1768668 : G(N[2],1) += f * n[1];
441 [ + - ]: 1768668 : G(N[2],2) += f * n[2];
442 : : }
443 : :
444 : : #if defined(__clang__)
445 : : #pragma clang diagnostic pop
446 : : #elif defined(STRICT_GNUC)
447 : : #pragma GCC diagnostic pop
448 : : #endif
449 : 10142 : }
450 : :
451 : : static tk::real
452 : 1229031 : flux( const tk::Fields& U,
453 : : const tk::Fields& G,
454 : : std::size_t i,
455 : : std::size_t j,
456 : : std::size_t p,
457 : : std::size_t q )
458 : : // *****************************************************************************
459 : : //! Compute momentum flux over edge of points p-q
460 : : //! \param[in] U Velocity vector
461 : : //! \param[in] G Velocity gradients
462 : : //! \param[in] i Tensor component, 1st index
463 : : //! \param[in] j Tensor component, 2nd index
464 : : //! \param[in] p Left node index of edge
465 : : //! \param[in] q Right node index of edge
466 : : //! \return Momentum flux contribution for edge p-q
467 : : // *****************************************************************************
468 : : {
469 : 1229031 : auto inv = U(p,i)*U(p,j) + U(q,i)*U(q,j);
470 : :
471 : 1229031 : auto eps = std::numeric_limits< tk::real >::epsilon();
472 : 1229031 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
473 [ + + ]: 1229031 : if (mu < eps) return -inv;
474 : :
475 : 528939 : auto vis = G(p,i*3+j) + G(p,j*3+i) + G(q,i*3+j) + G(q,j*3+i);
476 [ + + ]: 528939 : if (i == j) {
477 : 176313 : vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) + G(q,0) + G(q,4) + G(q,8) );
478 : : }
479 : 528939 : return mu*vis - inv;
480 : : }
481 : :
482 : : static tk::real
483 : 1277262 : flux( const tk::Fields& U,
484 : : const tk::Fields& G,
485 : : std::size_t i,
486 : : std::size_t j,
487 : : std::size_t p )
488 : : // *****************************************************************************
489 : : //! Compute momentum flux in point p
490 : : //! \param[in] U Velocity vector
491 : : //! \param[in] G Velocity gradients
492 : : //! \param[in] i Tensor component, 1st index
493 : : //! \param[in] j Tensor component, 2nd index
494 : : //! \param[in] p Node index of point
495 : : //! \return Momentum flux contribution for point p
496 : : // *****************************************************************************
497 : : {
498 : 1277262 : auto inv = U(p,i)*U(p,j);
499 : :
500 : 1277262 : auto eps = std::numeric_limits< tk::real >::epsilon();
501 : 1277262 : auto mu = g_cfg.get< tag::mat_dyn_viscosity >();
502 [ + + ]: 1277262 : if (mu < eps) return -inv;
503 : :
504 : 654642 : auto vis = G(p,i*3+j) + G(p,j*3+i);
505 [ + + ]: 654642 : if (i == j) {
506 : 218214 : vis -= 2.0/3.0 * ( G(p,0) + G(p,4) + G(p,8) );
507 : : }
508 : 654642 : return mu*vis - inv;
509 : : }
510 : :
511 : : void
512 : 385 : flux( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
513 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
514 : : const std::array< std::vector< tk::real >, 3 >& coord,
515 : : const std::vector< std::size_t >& triinpoel,
516 : : const tk::Fields& U,
517 : : const tk::Fields& G,
518 : : tk::Fields& F )
519 : : // *****************************************************************************
520 : : // Compute momentum flux in all points
521 : : //! \param[in] dsupedge Domain superedges
522 : : //! \param[in] dsupint Domain superedge integrals
523 : : //! \param[in] coord Mesh node coordinates
524 : : //! \param[in] triinpoel Boundary face connectivity
525 : : //! \param[in] U Velocity field
526 : : //! \param[in] G Velocity gradients, dui/dxj, 9 components
527 : : //! \param[in,out] F Momentum flux, Fi = d[ sij - uiuj ] / dxj, where
528 : : //! s_ij = mu[ dui/dxj + duj/dxi - 2/3 duk/dxk delta_ij ]
529 : : // *****************************************************************************
530 : : {
531 [ - + ][ - - ]: 385 : Assert( F.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
532 [ - + ][ - - ]: 385 : Assert( F.nprop() == 3, "Size mismatch" );
[ - - ][ - - ]
533 [ - + ][ - - ]: 385 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
534 [ - + ][ - - ]: 385 : Assert( G.nprop() == U.nprop()*3, "Size mismatch" );
[ - - ][ - - ]
535 : :
536 : : #if defined(__clang__)
537 : : #pragma clang diagnostic push
538 : : #pragma clang diagnostic ignored "-Wvla"
539 : : #pragma clang diagnostic ignored "-Wvla-extension"
540 : : #elif defined(STRICT_GNUC)
541 : : #pragma GCC diagnostic push
542 : : #pragma GCC diagnostic ignored "-Wvla"
543 : : #endif
544 : :
545 : : // domain integral
546 : :
547 : : // domain edge contributions: tetrahedron superedges
548 [ + + ]: 12352 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
549 : 11967 : const auto N = dsupedge[0].data() + e*4;
550 : 11967 : const auto d = dsupint[0].data();
551 [ + + ]: 47868 : for (std::size_t i=0; i<3; ++i) {
552 [ + + ]: 143604 : for (std::size_t j=0; j<3; ++j) {
553 [ + - ]: 107703 : tk::real f[] = { d[(e*6+0)*5+j] * flux(U,G,i,j,N[1],N[0]),
554 : 215406 : d[(e*6+1)*5+j] * flux(U,G,i,j,N[2],N[1]),
555 : 215406 : d[(e*6+2)*5+j] * flux(U,G,i,j,N[0],N[2]),
556 : 215406 : d[(e*6+3)*5+j] * flux(U,G,i,j,N[3],N[0]),
557 : 215406 : d[(e*6+4)*5+j] * flux(U,G,i,j,N[3],N[1]),
558 [ + - ][ + - ]: 107703 : d[(e*6+5)*5+j] * flux(U,G,i,j,N[3],N[2]) };
[ + - ][ + - ]
[ + - ]
559 [ + - ][ + - ]: 107703 : F(N[0],i) = F(N[0],i) - f[0] + f[2] - f[3];
560 [ + - ][ + - ]: 107703 : F(N[1],i) = F(N[1],i) + f[0] - f[1] - f[4];
561 [ + - ][ + - ]: 107703 : F(N[2],i) = F(N[2],i) + f[1] - f[2] - f[5];
562 [ + - ][ + - ]: 107703 : F(N[3],i) = F(N[3],i) + f[3] + f[4] + f[5];
563 : : }
564 : : }
565 : : }
566 : :
567 : : // domain edge contributions: triangle superedges
568 [ + + ]: 11674 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
569 : 11289 : const auto N = dsupedge[1].data() + e*3;
570 : 11289 : const auto d = dsupint[1].data();
571 [ + + ]: 45156 : for (std::size_t i=0; i<3; ++i) {
572 [ + + ]: 135468 : for (std::size_t j=0; j<3; ++j) {
573 [ + - ]: 101601 : tk::real f[] = { d[(e*3+0)*5+j] * flux(U,G,i,j,N[1],N[0]),
574 : 203202 : d[(e*3+1)*5+j] * flux(U,G,i,j,N[2],N[1]),
575 [ + - ][ + - ]: 101601 : d[(e*3+2)*5+j] * flux(U,G,i,j,N[0],N[2]), };
576 [ + - ][ + - ]: 101601 : F(N[0],i) = F(N[0],i) - f[0] + f[2];
577 [ + - ][ + - ]: 101601 : F(N[1],i) = F(N[1],i) + f[0] - f[1];
578 [ + - ][ + - ]: 101601 : F(N[2],i) = F(N[2],i) + f[1] - f[2];
579 : : }
580 : : }
581 : : }
582 : :
583 : : // domain edge contributions: edges
584 [ + + ]: 31275 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
585 : 30890 : const auto N = dsupedge[2].data() + e*2;
586 : 30890 : const auto d = dsupint[2].data() + e*5;
587 [ + + ]: 123560 : for (std::size_t i=0; i<3; ++i) {
588 [ + + ]: 370680 : for (std::size_t j=0; j<3; ++j) {
589 : 278010 : tk::real f = d[j] * flux(U,G,i,j,N[1],N[0]);
590 : 278010 : F(N[0],i) -= f;
591 : 278010 : F(N[1],i) += f;
592 : : }
593 : : }
594 : : }
595 : :
596 : : // boundary integral
597 : :
598 : 385 : const auto& x = coord[0];
599 : 385 : const auto& y = coord[1];
600 : 385 : const auto& z = coord[2];
601 : :
602 [ + + ]: 47691 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
603 : 47306 : const auto N = triinpoel.data() + e*3;
604 : : tk::real n[3];
605 : 47306 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
606 : 47306 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
607 : : n[0], n[1], n[2] );
608 [ + + ]: 189224 : for (std::size_t i=0; i<3; ++i) {
609 [ + - ]: 141918 : auto fxA = flux(U,G,i,0,N[0]);
610 [ + - ]: 141918 : auto fyA = flux(U,G,i,1,N[0]);
611 [ + - ]: 141918 : auto fzA = flux(U,G,i,2,N[0]);
612 [ + - ]: 141918 : auto fxB = flux(U,G,i,0,N[1]);
613 [ + - ]: 141918 : auto fyB = flux(U,G,i,1,N[1]);
614 [ + - ]: 141918 : auto fzB = flux(U,G,i,2,N[1]);
615 [ + - ]: 141918 : auto fxC = flux(U,G,i,0,N[2]);
616 [ + - ]: 141918 : auto fyC = flux(U,G,i,1,N[2]);
617 [ + - ]: 141918 : auto fzC = flux(U,G,i,2,N[2]);
618 : 141918 : auto fx = (6.0*fxA + fxB + fxC)/8.0;
619 : 141918 : auto fy = (6.0*fyA + fyB + fyC)/8.0;
620 : 141918 : auto fz = (6.0*fzA + fzB + fzC)/8.0;
621 [ + - ]: 141918 : F(N[0],i) += fx*n[0] + fy*n[1] + fz*n[2];
622 : 141918 : fx = (fxA + 6.0*fxB + fxC)/8.0;
623 : 141918 : fy = (fyA + 6.0*fyB + fyC)/8.0;
624 : 141918 : fz = (fzA + 6.0*fzB + fzC)/8.0;
625 [ + - ]: 141918 : F(N[1],i) += fx*n[0] + fy*n[1] + fz*n[2];
626 : 141918 : fx = (fxA + fxB + 6.0*fxC)/8.0;
627 : 141918 : fy = (fyA + fyB + 6.0*fyC)/8.0;
628 : 141918 : fz = (fzA + fzB + 6.0*fzC)/8.0;
629 [ + - ]: 141918 : F(N[2],i) += fx*n[0] + fy*n[1] + fz*n[2];
630 : : }
631 : : }
632 : :
633 : : #if defined(__clang__)
634 : : #pragma clang diagnostic pop
635 : : #elif defined(STRICT_GNUC)
636 : : #pragma GCC diagnostic pop
637 : : #endif
638 : 385 : }
639 : :
640 : : static void
641 : 2527630 : adv_damp2( const tk::real supint[],
642 : : const tk::Fields& U,
643 : : const tk::Fields&,
644 : : const std::vector< tk::real >& P,
645 : : const std::array< std::vector< tk::real >, 3 >&,
646 : : std::size_t p,
647 : : std::size_t q,
648 : : tk::real f[] )
649 : : // *****************************************************************************
650 : : //! Compute advection fluxes on a single edge using 2nd-order damping
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] p Left node index of edge
655 : : //! \param[in] q Right node index of edge
656 : : //! \param[in,out] f Flux computed
657 : : // *****************************************************************************
658 : : {
659 : 2527630 : auto nx = supint[0];
660 : 2527630 : auto ny = supint[1];
661 : 2527630 : auto nz = supint[2];
662 : :
663 : : // left state
664 : 2527630 : auto uL = U(p,0);
665 : 2527630 : auto vL = U(p,1);
666 : 2527630 : auto wL = U(p,2);
667 : 2527630 : auto vnL = uL*nx + vL*ny + wL*nz;
668 : :
669 : : // right state
670 : 2527630 : auto uR = U(q,0);
671 : 2527630 : auto vR = U(q,1);
672 : 2527630 : auto wR = U(q,2);
673 : 2527630 : auto vnR = uR*nx + vR*ny + wR*nz;
674 : :
675 : : // stabilization
676 : 2527630 : tk::real aw = 0.0;
677 [ + - ]: 2527630 : if (g_cfg.get< tag::stab >()) {
678 : 2527630 : aw = std::abs( vnL + vnR ) / 2.0;
679 : : }
680 : :
681 : : // artificial viscosity
682 [ - + ]: 2527630 : if (g_cfg.get< tag::stab2 >()) {
683 : 0 : auto stab2coef = g_cfg.get< tag::stab2coef >();
684 : 0 : auto sl = std::abs(vnL);
685 : 0 : auto sr = std::abs(vnR);
686 : 0 : aw += stab2coef * std::max( sl, sr );
687 : : }
688 : :
689 : : // viscosity
690 : 2527630 : auto v = supint[4] * g_cfg.get< tag::mat_dyn_viscosity >();
691 : :
692 : : // flow
693 : 2527630 : auto pf = P[p] + P[q];
694 : 2527630 : f[0] = uL*vnL + uR*vnR + pf*nx + (aw-v)*(uR-uL);
695 : 2527630 : f[1] = vL*vnL + vR*vnR + pf*ny + (aw-v)*(vR-vL);
696 : 2527630 : f[2] = wL*vnL + wR*vnR + pf*nz + (aw-v)*(wR-wL);
697 : :
698 : : // scalar
699 : 2527630 : auto ncomp = U.nprop();
700 [ + + ]: 2527630 : if (ncomp == 3) return;
701 : :
702 : : // diffusion
703 : 1311500 : auto d = supint[4] * g_cfg.get< tag::mat_dyn_diffusivity >();
704 : :
705 : : // scalar
706 [ + + ]: 2623000 : for (std::size_t c=3; c<ncomp; ++c) {
707 : 1311500 : f[c] = U(p,c)*vnL + U(q,c)*vnR + (aw-d)*(U(q,c)-U(p,c));
708 : : }
709 : : }
710 : :
711 : : static void
712 : 2291640 : adv_damp4( const tk::real supint[],
713 : : const tk::Fields& U,
714 : : const tk::Fields& G,
715 : : const std::vector< tk::real >& P,
716 : : const std::array< std::vector< tk::real >, 3 >& coord,
717 : : std::size_t p,
718 : : std::size_t q,
719 : : tk::real f[] )
720 : : // *****************************************************************************
721 : : //! Compute advection fluxes on a single edge using 4th-order damping
722 : : //! \param[in] supint Edge integral
723 : : //! \param[in] U Velocity and transported scalars at recent time step
724 : : //! \param[in] G Gradients of velocity and transported scalars
725 : : //! \param[in] P Pressure
726 : : //! \param[in] coord Mesh node coordinates
727 : : //! \param[in] p Left node index of edge
728 : : //! \param[in] q Right node index of edge
729 : : //! \param[in,out] f Flux computed
730 : : // *****************************************************************************
731 : : {
732 : 2291640 : const auto& x = coord[0];
733 : 2291640 : const auto& y = coord[1];
734 : 2291640 : const auto& z = coord[2];
735 : :
736 : : // edge vector
737 : 2291640 : auto dx = x[q] - x[p];
738 : 2291640 : auto dy = y[q] - y[p];
739 : 2291640 : auto dz = z[q] - z[p];
740 : :
741 : 2291640 : auto ncomp = U.nprop();
742 : :
743 : : #if defined(__clang__)
744 : : #pragma clang diagnostic push
745 : : #pragma clang diagnostic ignored "-Wvla"
746 : : #pragma clang diagnostic ignored "-Wvla-extension"
747 : : #elif defined(STRICT_GNUC)
748 : : #pragma GCC diagnostic push
749 : : #pragma GCC diagnostic ignored "-Wvla"
750 : : #endif
751 : :
752 : 2291640 : tk::real uL[ncomp], uR[ncomp];
753 [ + + ]: 10133280 : for (std::size_t i=0; i<ncomp; ++i) {
754 [ + - ]: 7841640 : uL[i] = U(p,i);
755 [ + - ]: 7841640 : uR[i] = U(q,i);
756 : : }
757 : :
758 : : // MUSCL reconstruction in edge-end points
759 [ + + ]: 10133280 : for (std::size_t c=0; c<ncomp; ++c) {
760 [ + - ][ + - ]: 7841640 : auto g1 = G(p,c*3+0)*dx + G(p,c*3+1)*dy + G(p,c*3+2)*dz;
[ + - ]
761 [ + - ][ + - ]: 7841640 : auto g2 = G(q,c*3+0)*dx + G(q,c*3+1)*dy + G(q,c*3+2)*dz;
[ + - ]
762 : 7841640 : auto delta2 = uR[c] - uL[c];
763 : 7841640 : auto delta1 = 2.0 * g1 - delta2;
764 : 7841640 : auto delta3 = 2.0 * g2 - delta2;
765 : :
766 : : // van Leer limiter
767 : 7841640 : auto rL = (delta2 + muscl_eps) / (delta1 + muscl_eps);
768 : 7841640 : auto rR = (delta2 + muscl_eps) / (delta3 + muscl_eps);
769 : 7841640 : auto rLinv = (delta1 + muscl_eps) / (delta2 + muscl_eps);
770 : 7841640 : auto rRinv = (delta3 + muscl_eps) / (delta2 + muscl_eps);
771 : 7841640 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
772 : 7841640 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
773 : 7841640 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
774 : 7841640 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
775 : : // update unknowns with reconstructed unknowns
776 : 7841640 : uL[c] += 0.25*(delta1*(1.0-muscl_const)*phiL +
777 : 7841640 : delta2*(1.0+muscl_const)*phi_L_inv);
778 : 7841640 : uR[c] -= 0.25*(delta3*(1.0-muscl_const)*phiR +
779 : 7841640 : delta2*(1.0+muscl_const)*phi_R_inv);
780 : : }
781 : :
782 : 2291640 : auto nx = supint[0];
783 : 2291640 : auto ny = supint[1];
784 : 2291640 : auto nz = supint[2];
785 : :
786 : : // normal velocities
787 : 2291640 : auto vnL = uL[0]*nx + uL[1]*ny + uL[2]*nz;
788 : 2291640 : auto vnR = uR[0]*nx + uR[1]*ny + uR[2]*nz;
789 : :
790 : : // stabilization
791 : 2291640 : tk::real aw = 0.0;
792 [ + - ]: 2291640 : if (g_cfg.get< tag::stab >()) {
793 : 2291640 : aw = std::abs( vnL + vnR ) / 2.0;
794 : : }
795 : :
796 : : // artificial viscosity
797 [ - + ]: 2291640 : if (g_cfg.get< tag::stab2 >()) {
798 : 0 : auto stab2coef = g_cfg.get< tag::stab2coef >();
799 : 0 : auto sl = std::abs(vnL);
800 : 0 : auto sr = std::abs(vnR);
801 : 0 : aw += stab2coef * std::max( sl, sr );
802 : : }
803 : :
804 : : // viscosity
805 : 2291640 : auto v = supint[4] * g_cfg.get< tag::mat_dyn_viscosity >();
806 : :
807 : : // flow
808 : 2291640 : auto pf = P[p] + P[q];
809 [ + - ][ + - ]: 2291640 : f[0] = uL[0]*vnL + uR[0]*vnR + pf*nx + aw*(uR[0]-uL[0]) - v*(U(q,0)-U(p,0));
810 [ + - ][ + - ]: 2291640 : f[1] = uL[1]*vnL + uR[1]*vnR + pf*ny + aw*(uR[1]-uL[1]) - v*(U(q,1)-U(p,1));
811 [ + - ][ + - ]: 2291640 : f[2] = uL[2]*vnL + uR[2]*vnR + pf*nz + aw*(uR[2]-uL[2]) - v*(U(q,2)-U(p,2));
812 : :
813 : : // scalar
814 [ + + ]: 2291640 : if (ncomp == 3) return;
815 : :
816 : : // diffusion
817 : 966720 : auto d = supint[4] * g_cfg.get< tag::mat_dyn_diffusivity >();
818 : :
819 : : // scalar
820 [ + + ]: 1933440 : for (std::size_t c=3; c<ncomp; ++c) {
821 [ + - ][ + - ]: 966720 : f[c] = uL[c]*vnL + uR[c]*vnR + aw*(uR[c]-uL[c]) - d*(U(q,c)-U(p,c));
822 : : }
823 : :
824 : : #if defined(__clang__)
825 : : #pragma clang diagnostic pop
826 : : #elif defined(STRICT_GNUC)
827 : : #pragma GCC diagnostic pop
828 : : #endif
829 : 2291640 : }
830 : :
831 : : static void
832 : 6930 : adv( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
833 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
834 : : const std::array< std::vector< tk::real >, 3 >& coord,
835 : : const std::vector< std::size_t >& triinpoel,
836 : : const tk::Fields& U,
837 : : const tk::Fields& G,
838 : : const std::vector< tk::real >& P,
839 : : // cppcheck-suppress constParameter
840 : : tk::Fields& R )
841 : : // *****************************************************************************
842 : : //! Add advection to rhs
843 : : //! \param[in] dsupedge Domain superedges
844 : : //! \param[in] dsupint Domain superedge integrals
845 : : //! \param[in] coord Mesh node coordinates
846 : : //! \param[in] triinpoel Boundary face connectivity
847 : : //! \param[in] U Velocity and transported scalars at recent time step
848 : : //! \param[in] G Gradients of velocity and transported scalars
849 : : //! \param[in] P Pressure
850 : : //! \param[in,out] R Right-hand side vector added to
851 : : // *****************************************************************************
852 : : {
853 : : // configure advection stabilization
854 : 0 : auto adv = [](){
855 : 6930 : const auto& flux = g_cfg.get< tag::flux >();
856 [ + + ]: 6930 : if (flux == "damp2") return adv_damp2;
857 [ + - ]: 3320 : else if (flux == "damp4") return adv_damp4;
858 [ - - ][ - - ]: 0 : else Throw( "Flux not correctly configured" );
[ - - ]
859 [ + - ]: 6930 : }();
860 : :
861 : 6930 : auto ncomp = U.nprop();
862 : :
863 : : #if defined(__clang__)
864 : : #pragma clang diagnostic push
865 : : #pragma clang diagnostic ignored "-Wvla"
866 : : #pragma clang diagnostic ignored "-Wvla-extension"
867 : : #elif defined(STRICT_GNUC)
868 : : #pragma GCC diagnostic push
869 : : #pragma GCC diagnostic ignored "-Wvla"
870 : : #endif
871 : :
872 : : // domain integral
873 : :
874 : : // domain edge contributions: tetrahedron superedges
875 [ + + ]: 428510 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
876 : 421580 : const auto N = dsupedge[0].data() + e*4;
877 : 421580 : tk::real f[6][ncomp];
878 : 421580 : const auto d = dsupint[0].data();
879 [ + - ]: 421580 : adv( d+(e*6+0)*5, U, G, P, coord, N[0], N[1], f[0] );
880 [ + - ]: 421580 : adv( d+(e*6+1)*5, U, G, P, coord, N[1], N[2], f[1] );
881 [ + - ]: 421580 : adv( d+(e*6+2)*5, U, G, P, coord, N[2], N[0], f[2] );
882 [ + - ]: 421580 : adv( d+(e*6+3)*5, U, G, P, coord, N[0], N[3], f[3] );
883 [ + - ]: 421580 : adv( d+(e*6+4)*5, U, G, P, coord, N[1], N[3], f[4] );
884 [ + - ]: 421580 : adv( d+(e*6+5)*5, U, G, P, coord, N[2], N[3], f[5] );
885 [ + + ]: 1876560 : for (std::size_t c=0; c<ncomp; ++c) {
886 [ + - ][ + - ]: 1454980 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
887 [ + - ][ + - ]: 1454980 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
888 [ + - ][ + - ]: 1454980 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
889 [ + - ][ + - ]: 1454980 : R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
890 : : }
891 : 421580 : }
892 : :
893 : : // domain edge contributions: triangle superedges
894 [ + + ]: 401780 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
895 : 394850 : const auto N = dsupedge[1].data() + e*3;
896 : 394850 : tk::real f[3][ncomp];
897 : 394850 : const auto d = dsupint[1].data();
898 [ + - ]: 394850 : adv( d+(e*3+0)*5, U, G, P, coord, N[0], N[1], f[0] );
899 [ + - ]: 394850 : adv( d+(e*3+1)*5, U, G, P, coord, N[1], N[2], f[1] );
900 [ + - ]: 394850 : adv( d+(e*3+2)*5, U, G, P, coord, N[2], N[0], f[2] );
901 [ + + ]: 1758720 : for (std::size_t c=0; c<ncomp; ++c) {
902 [ + - ][ + - ]: 1363870 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
903 [ + - ][ + - ]: 1363870 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
904 [ + - ][ + - ]: 1363870 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
905 : : }
906 : 394850 : }
907 : :
908 : : // domain edge contributions: edges
909 [ + + ]: 1112170 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
910 : 1105240 : const auto N = dsupedge[2].data() + e*2;
911 : 1105240 : tk::real f[ncomp];
912 : 1105240 : const auto d = dsupint[2].data();
913 [ + - ]: 1105240 : adv( d+e*5, U, G, P, coord, N[0], N[1], f );
914 [ + + ]: 5019780 : for (std::size_t c=0; c<ncomp; ++c) {
915 [ + - ]: 3914540 : R(N[0],c) -= f[c];
916 [ + - ]: 3914540 : R(N[1],c) += f[c];
917 : : }
918 : 1105240 : }
919 : :
920 : : // boundary integral
921 : :
922 : 6930 : const auto& x = coord[0];
923 : 6930 : const auto& y = coord[1];
924 : 6930 : const auto& z = coord[2];
925 : :
926 [ + + ]: 1782970 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
927 : 3552080 : const auto N = triinpoel.data() + e*3;
928 : : tk::real n[3];
929 : 1776040 : tk::crossdiv( x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
930 : 1776040 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]], 6.0,
931 : : n[0], n[1], n[2] );
932 : 1776040 : tk::real f[ncomp][3];
933 : :
934 [ + - ]: 1776040 : auto u = U(N[0],0);
935 [ + - ]: 1776040 : auto v = U(N[0],1);
936 [ + - ]: 1776040 : auto w = U(N[0],2);
937 : 1776040 : auto p = P[N[0]];
938 : 1776040 : auto vn = n[0]*u + n[1]*v + n[2]*w;
939 : : // flow
940 : 1776040 : f[0][0] = u*vn + p*n[0];
941 : 1776040 : f[1][0] = v*vn + p*n[1];
942 : 1776040 : f[2][0] = w*vn + p*n[2];
943 : : // scalar
944 [ + - ][ + + ]: 2616200 : for (std::size_t c=3; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;
945 : :
946 [ + - ]: 1776040 : u = U(N[1],0);
947 [ + - ]: 1776040 : v = U(N[1],1);
948 [ + - ]: 1776040 : w = U(N[1],2);
949 : 1776040 : p = P[N[1]];
950 : 1776040 : vn = n[0]*u + n[1]*v + n[2]*w;
951 : : // flow
952 : 1776040 : f[0][1] = u*vn + p*n[0];
953 : 1776040 : f[1][1] = v*vn + p*n[1];
954 : 1776040 : f[2][1] = w*vn + p*n[2];
955 : : // scalar
956 [ + - ][ + + ]: 2616200 : for (std::size_t c=3; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;
957 : :
958 [ + - ]: 1776040 : u = U(N[2],0);
959 [ + - ]: 1776040 : v = U(N[2],1);
960 [ + - ]: 1776040 : w = U(N[2],2);
961 : 1776040 : p = P[N[2]];
962 : 1776040 : vn = n[0]*u + n[1]*v + n[2]*w;
963 : : // flow
964 : 1776040 : f[0][2] = u*vn + p*n[0];
965 : 1776040 : f[1][2] = v*vn + p*n[1];
966 : 1776040 : f[2][2] = w*vn + p*n[2];
967 : : // scalar
968 [ + - ][ + + ]: 2616200 : for (std::size_t c=3; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;
969 : :
970 [ + + ]: 7944320 : for (std::size_t c=0; c<ncomp; ++c) {
971 [ + - ]: 6168280 : R(N[0],c) += (6.0*f[c][0] + f[c][1] + f[c][2])/8.0;
972 [ + - ]: 6168280 : R(N[1],c) += (f[c][0] + 6.0*f[c][1] + f[c][2])/8.0;
973 [ + - ]: 6168280 : R(N[2],c) += (f[c][0] + f[c][1] + 6.0*f[c][2])/8.0;
974 : : }
975 : 1776040 : }
976 : :
977 : : #if defined(__clang__)
978 : : #pragma clang diagnostic pop
979 : : #elif defined(STRICT_GNUC)
980 : : #pragma GCC diagnostic pop
981 : : #endif
982 : 6930 : }
983 : :
984 : : static void
985 : 6930 : src( const std::array< std::vector< tk::real >, 3 >& coord,
986 : : const std::vector< tk::real >& v,
987 : : tk::real t,
988 : : tk::Fields& R )
989 : : // *****************************************************************************
990 : : // Compute source integral
991 : : //! \param[in] coord Mesh node coordinates
992 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
993 : : //! \param[in] t Physical time
994 : : //! \param[in,out] R Right-hand side vector computed
995 : : // *****************************************************************************
996 : : {
997 [ + - ]: 6930 : auto src = problems::SRC();
998 [ + + ]: 6930 : if (!src) return;
999 : :
1000 : 3100 : const auto& x = coord[0];
1001 : 3100 : const auto& y = coord[1];
1002 : 3100 : const auto& z = coord[2];
1003 : :
1004 [ + + ]: 468320 : for (std::size_t p=0; p<R.nunk(); ++p) {
1005 [ + - ]: 465220 : auto s = src( x[p], y[p], z[p], t, /*meshid=*/0 );
1006 [ + - ][ + + ]: 2326100 : for (std::size_t c=0; c<s.size(); ++c) R(p,c) -= s[c] * v[p];
1007 : 465220 : }
1008 [ + + ]: 6930 : }
1009 : :
1010 : : void
1011 : 6930 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
1012 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
1013 : : const std::array< std::vector< tk::real >, 3 >& coord,
1014 : : const std::vector< std::size_t >& triinpoel,
1015 : : const std::vector< tk::real >& v,
1016 : : tk::real t,
1017 : : const std::vector< tk::real >& P,
1018 : : const tk::Fields& U,
1019 : : const tk::Fields& G,
1020 : : tk::Fields& R )
1021 : : // *****************************************************************************
1022 : : // Compute right hand side
1023 : : //! \param[in] dsupedge Domain superedges
1024 : : //! \param[in] dsupint Domain superedge integrals
1025 : : //! \param[in] coord Mesh node coordinates
1026 : : //! \param[in] triinpoel Boundary face connectivity
1027 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
1028 : : //! \param[in] t Physical time
1029 : : //! \param[in] P Pressure
1030 : : //! \param[in] U Solution vector of primitive variables at recent time step
1031 : : //! \param[in] G Gradients of velocity and transported scalars
1032 : : //! \param[in,out] R Right-hand side vector computed
1033 : : // *****************************************************************************
1034 : : {
1035 [ - + ][ - - ]: 6930 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
1036 : : "vector at recent time step incorrect" );
1037 [ - + ][ - - ]: 6930 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
1038 : : "Number of unknowns and/or number of components in right-hand "
1039 : : "side vector incorrect" );
1040 : :
1041 : 6930 : R.fill( 0.0 );
1042 : 6930 : adv( dsupedge, dsupint, coord, triinpoel, U, G, P, R );
1043 : 6930 : src( coord, v, t, R );
1044 : 6930 : }
1045 : :
1046 : : } // chorin::
|
