Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Physics/Lax.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 LaxCG: Time-derivative preconditioning for all Ma
10 : : \see Luo, Baum, Lohner, "Extension of Harten-Lax-van Leer Scheme for
11 : : Flows at All Speeds", AIAA Journal, Vol. 43, No. 6, 2005
12 : : \see Weiss & Smith, "Preconditioning Applied to Variable and Constant
13 : : Density Time-Accurate Flows on Unstructured Meshes", AIAA Journal,
14 : : Vol. 33, No. 11, 1995, pp. 2050-2057.
15 : : */
16 : : // *****************************************************************************
17 : :
18 : : #include "Vector.hpp"
19 : : #include "Around.hpp"
20 : : #include "DerivedData.hpp"
21 : : #include "EOS.hpp"
22 : : #include "Lax.hpp"
23 : : #include "Problems.hpp"
24 : : #include "InciterConfig.hpp"
25 : :
26 : : namespace inciter {
27 : :
28 : : extern ctr::Config g_cfg;
29 : :
30 : : } // ::inciter
31 : :
32 : : namespace lax {
33 : :
34 : : static const tk::real muscl_eps = 1.0e-9;
35 : : static const tk::real muscl_const = 1.0/3.0;
36 : :
37 : : using inciter::g_cfg;
38 : :
39 : : static void
40 : 19506570 : muscl( std::size_t p,
41 : : std::size_t q,
42 : : const tk::UnsMesh::Coords& coord,
43 : : const tk::Fields& G,
44 : : tk::real& pL, tk::real& uL, tk::real& vL, tk::real& wL, tk::real& tL,
45 : : tk::real& pR, tk::real& uR, tk::real& vR, tk::real& wR, tk::real& tR )
46 : : // *****************************************************************************
47 : : //! Compute MUSCL reconstruction in edge-end points for the flow variables
48 : : //! \param[in] p Left node id of edge-end
49 : : //! \param[in] q Right node id of edge-end
50 : : //! \param[in] coord Array of nodal coordinates
51 : : //! \param[in] G Gradient of all unknowns in mesh points
52 : : //! \param[in,out] pL Left pressure
53 : : //! \param[in,out] uL Left X velocity
54 : : //! \param[in,out] vL Left Y velocity
55 : : //! \param[in,out] wL Left Z velocity
56 : : //! \param[in,out] tL Left temperature
57 : : //! \param[in,out] pR Right pressure
58 : : //! \param[in,out] uR Right X velocity
59 : : //! \param[in,out] vR Right Y velocity
60 : : //! \param[in,out] wR Right Z velocity
61 : : //! \param[in,out] tR Right temperature
62 : : // *****************************************************************************
63 : : {
64 : : // access node coordinates
65 : 19506570 : const auto& x = coord[0];
66 : 19506570 : const auto& y = coord[1];
67 : 19506570 : const auto& z = coord[2];
68 : :
69 : : // edge vector
70 : 19506570 : tk::real vw[3] = { x[q]-x[p], y[q]-y[p], z[q]-z[p] };
71 : :
72 : : tk::real delta1[5], delta2[5], delta3[5];
73 : 19506570 : tk::real ls[5] = { pL, uL, vL, wL, tL },
74 : 19506570 : rs[5] = { pR, uR, vR, wR, tR },
75 : : url[5], urr[5];
76 : 19506570 : memcpy( url, ls, sizeof ls );
77 : 19506570 : memcpy( urr, rs, sizeof rs );
78 : :
79 : : // MUSCL reconstruction of edge-end-point primitive variables
80 [ + + ]: 117039420 : for (std::size_t c=0; c<5; ++c) {
81 : :
82 [ + - ][ + - ]: 97532850 : auto g1 = G(p,c*3+0)*vw[0] + G(p,c*3+1)*vw[1] + G(p,c*3+2)*vw[2];
[ + - ]
83 [ + - ][ + - ]: 97532850 : auto g2 = G(q,c*3+0)*vw[0] + G(q,c*3+1)*vw[1] + G(q,c*3+2)*vw[2];
[ + - ]
84 : :
85 : 97532850 : delta2[c] = rs[c] - ls[c];
86 : 97532850 : delta1[c] = 2.0 * g1 - delta2[c];
87 : 97532850 : delta3[c] = 2.0 * g2 - delta2[c];
88 : :
89 : : // MUSCL extrapolation option 1:
90 : : // ---------------------------------------------------------------------
91 : : // See Waltz, J., Morgan, N. R., Canfield, T. R., Charest, M. R., Risinger,
92 : : // L. D., & Wohlbier, J. G. (2014). A three-dimensional finite element
93 : : // arbitrary Lagrangian–Eulerian method for shock hydrodynamics on
94 : : // unstructured grids. Computers & Fluids, 92, 172-187.
95 : :
96 : : // van Leer limiter
97 : 97532850 : auto rcL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
98 : 97532850 : auto rcR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
99 : 97532850 : auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
100 : 97532850 : auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
101 : 97532850 : auto phiL = (std::abs(rcL) + rcL) / (std::abs(rcL) + 1.0);
102 : 97532850 : auto phiR = (std::abs(rcR) + rcR) / (std::abs(rcR) + 1.0);
103 : 97532850 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
104 : 97532850 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
105 : : // update unknowns with reconstructed unknowns
106 : 97532850 : url[c] += 0.25*(delta1[c]*(1.0-muscl_const)*phiL +
107 : 97532850 : delta2[c]*(1.0+muscl_const)*phi_L_inv);
108 : 97532850 : urr[c] -= 0.25*(delta3[c]*(1.0-muscl_const)*phiR +
109 : 97532850 : delta2[c]*(1.0+muscl_const)*phi_R_inv);
110 : :
111 : : // MUSCL extrapolation option 2:
112 : : // ---------------------------------------------------------------------
113 : : // See Luo, H., Baum, J. D., & Lohner, R. (1994). Edge-based finite element
114 : : // scheme for the Euler equations. AIAA journal, 32(6), 1183-1190.
115 : : // van Leer, B. (1974). Towards the ultimate conservative difference
116 : : // scheme. II. Monotonicity and conservation combined in a second-order
117 : : // scheme. Journal of computational physics, 14(4), 361-370.
118 : : // Derived from the flux limiter phi as: s = phi_inv - (1 - phi)
119 : :
120 : : // van Albada limiter
121 : : //auto sL = std::max(0.0, (2.0*delta1[c]*delta2[c] + muscl_eps)
122 : : // /(delta1[c]*delta1[c] + delta2[c]*delta2[c] + muscl_eps));
123 : : //auto sR = std::max(0.0, (2.0*delta3[c]*delta2[c] + muscl_eps)
124 : : // /(delta3[c]*delta3[c] + delta2[c]*delta2[c] + muscl_eps));
125 : : //// update unknowns with reconstructed unknowns
126 : : //url[c] += 0.25*sL*(delta1[c]*(1.0 - muscl_const*sL)
127 : : // + delta2[c]*(1.0 + muscl_const*sL));
128 : : //urr[c] -= 0.25*sR*(delta3[c]*(1.0 - muscl_const*sR)
129 : : // + delta2[c]*(1.0 + muscl_const*sR));
130 : : }
131 : :
132 : : // force first order if the reconstructions for pressure or temperature
133 : : // would have allowed negative values
134 [ + - ][ - + ]: 19506570 : if (ls[0] < delta1[0] || ls[4] < delta1[4]) memcpy( url, ls, sizeof ls );
135 [ + - ][ - + ]: 19506570 : if (rs[0] < -delta3[0] || rs[4] < -delta3[4]) memcpy( urr, rs, sizeof rs );
136 : :
137 : 19506570 : pL = url[0];
138 : 19506570 : uL = url[1];
139 : 19506570 : vL = url[2];
140 : 19506570 : wL = url[3];
141 : 19506570 : tL = url[4];
142 : :
143 : 19506570 : pR = urr[0];
144 : 19506570 : uR = urr[1];
145 : 19506570 : vR = urr[2];
146 : 19506570 : wR = urr[3];
147 : 19506570 : tR = urr[4];
148 : 19506570 : }
149 : :
150 : : static void
151 : 0 : muscl( std::size_t p, std::size_t q, const tk::UnsMesh::Coords& coord,
152 : : const tk::Fields& G, tk::real uL[], tk::real uR[] )
153 : : // *****************************************************************************
154 : : //! Compute MUSCL reconstruction in edge-end points for transported scalars
155 : : //! \param[in] p Left node id of edge-end
156 : : //! \param[in] q Right node id of edge-end
157 : : //! \param[in] coord Array of nodal coordinates
158 : : //! \param[in] G Gradient of all unknowns in mesh points
159 : : //! \param[in,out] uL Primitive variables at left edge-end point
160 : : //! \param[in,out] uR Primitive variables at right edge-end point
161 : : // *****************************************************************************
162 : : {
163 : : #if defined(__clang__)
164 : : #pragma clang diagnostic push
165 : : #pragma clang diagnostic ignored "-Wvla"
166 : : #pragma clang diagnostic ignored "-Wvla-extension"
167 : : #elif defined(STRICT_GNUC)
168 : : #pragma GCC diagnostic push
169 : : #pragma GCC diagnostic ignored "-Wvla"
170 : : #endif
171 : :
172 : 0 : auto ns = G.nprop() / 3 - 5;
173 : :
174 : 0 : const auto& x = coord[0];
175 : 0 : const auto& y = coord[1];
176 : 0 : const auto& z = coord[2];
177 : :
178 : : // edge vector
179 : 0 : tk::real vw[3] = { x[q]-x[p], y[q]-y[p], z[q]-z[p] };
180 : :
181 : 0 : tk::real delta1[ns], delta2[ns], delta3[ns];
182 : :
183 : : // MUSCL reconstruction of edge-end-point primitive variables
184 [ - - ]: 0 : for (std::size_t c=0; c<ns; ++c) {
185 : 0 : auto g = (5+c)*3;
186 [ - - ][ - - ]: 0 : auto g1 = G(p,g+0)*vw[0] + G(p,g+1)*vw[1] + G(p,g+2)*vw[2];
[ - - ]
187 [ - - ][ - - ]: 0 : auto g2 = G(q,g+0)*vw[0] + G(q,g+1)*vw[1] + G(q,g+2)*vw[2];
[ - - ]
188 : :
189 : 0 : delta2[c] = uR[5+c] - uL[5+c];
190 : 0 : delta1[c] = 2.0 * g1 - delta2[c];
191 : 0 : delta3[c] = 2.0 * g2 - delta2[c];
192 : :
193 : : // van Leer limiter
194 : 0 : auto rL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
195 : 0 : auto rR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
196 : 0 : auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
197 : 0 : auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
198 : 0 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
199 : 0 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
200 : 0 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
201 : 0 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
202 : : // update unknowns with reconstructed unknowns
203 : 0 : uL[5+c] += 0.25*(delta1[c]*(1.0-muscl_const)*phiL +
204 : 0 : delta2[c]*(1.0+muscl_const)*phi_L_inv);
205 : 0 : uR[5+c] -= 0.25*(delta3[c]*(1.0-muscl_const)*phiR +
206 : 0 : delta2[c]*(1.0+muscl_const)*phi_R_inv);
207 : : }
208 : :
209 : : #if defined(__clang__)
210 : : #pragma clang diagnostic pop
211 : : #elif defined(STRICT_GNUC)
212 : : #pragma GCC diagnostic pop
213 : : #endif
214 : 0 : }
215 : :
216 : : void
217 : 10350 : grad( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
218 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
219 : : const std::array< std::vector< tk::real >, 3 >& coord,
220 : : const std::vector< std::size_t >& triinpoel,
221 : : const tk::Fields& U,
222 : : tk::Fields& G )
223 : : // *****************************************************************************
224 : : // Compute nodal gradients of primitive variables in all points
225 : : //! \param[in] dsupedge Domain superedges
226 : : //! \param[in] dsupint Domain superedge integrals
227 : : //! \param[in] coord Mesh node coordinates
228 : : //! \param[in] triinpoel Boundary face connectivity
229 : : //! \param[in] U Solution vector of primitive variables at recent time step
230 : : //! \param[in,out] G Nodal gradients
231 : : //! \return Gradients of primitive variables in all mesh points
232 : : // *****************************************************************************
233 : : {
234 : : #if defined(__clang__)
235 : : #pragma clang diagnostic push
236 : : #pragma clang diagnostic ignored "-Wvla"
237 : : #pragma clang diagnostic ignored "-Wvla-extension"
238 : : #elif defined(STRICT_GNUC)
239 : : #pragma GCC diagnostic push
240 : : #pragma GCC diagnostic ignored "-Wvla"
241 : : #endif
242 : :
243 : : // cppcheck-suppress unreadVariable
244 : 10350 : auto ncomp = U.nprop();
245 : :
246 [ - + ][ - - ]: 10350 : Assert( G.nunk() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
247 [ - + ][ - - ]: 10350 : Assert( G.nprop() == ncomp*3, "Size mismatch" );
[ - - ][ - - ]
248 : 10350 : G.fill( 0.0 );
249 : :
250 : : // domain integral
251 : :
252 : : // domain edge contributions: tetrahedron superedges
253 [ + + ]: 2100390 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
254 : 2090040 : const auto N = dsupedge[0].data() + e*4;
255 [ + + ]: 12540240 : for (std::size_t c=0; c<ncomp; ++c) {
256 [ + - ][ + - ]: 10450200 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c), U(N[3],c) };
[ + - ][ + - ]
257 [ + + ]: 41800800 : for (std::size_t j=0; j<3; ++j) {
258 : : tk::real f[6];
259 : 31350600 : const auto d = dsupint[0].data();
260 : 31350600 : f[0] = d[(e*6+0)*3+j] * (u[1] + u[0]);
261 : 31350600 : f[1] = d[(e*6+1)*3+j] * (u[2] + u[1]);
262 : 31350600 : f[2] = d[(e*6+2)*3+j] * (u[0] + u[2]);
263 : 31350600 : f[3] = d[(e*6+3)*3+j] * (u[3] + u[0]);
264 : 31350600 : f[4] = d[(e*6+4)*3+j] * (u[3] + u[1]);
265 : 31350600 : f[5] = d[(e*6+5)*3+j] * (u[3] + u[2]);
266 [ + - ][ + - ]: 31350600 : G(N[0],c*3+j) = G(N[0],c*3+j) - f[0] + f[2] - f[3];
267 [ + - ][ + - ]: 31350600 : G(N[1],c*3+j) = G(N[1],c*3+j) + f[0] - f[1] - f[4];
268 [ + - ][ + - ]: 31350600 : G(N[2],c*3+j) = G(N[2],c*3+j) + f[1] - f[2] - f[5];
269 [ + - ][ + - ]: 31350600 : G(N[3],c*3+j) = G(N[3],c*3+j) + f[3] + f[4] + f[5];
270 : : }
271 : : }
272 : : }
273 : :
274 : : // domain edge contributions: triangle superedges
275 [ + + ]: 884850 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
276 : 874500 : const auto N = dsupedge[1].data() + e*3;
277 [ + + ]: 5247000 : for (std::size_t c=0; c<ncomp; ++c) {
278 [ + - ][ + - ]: 4372500 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c) };
[ + - ]
279 [ + + ]: 17490000 : for (std::size_t j=0; j<3; ++j) {
280 : : tk::real f[3];
281 : 13117500 : const auto d = dsupint[1].data();
282 : 13117500 : f[0] = d[(e*3+0)*3+j] * (u[1] + u[0]);
283 : 13117500 : f[1] = d[(e*3+1)*3+j] * (u[2] + u[1]);
284 : 13117500 : f[2] = d[(e*3+2)*3+j] * (u[0] + u[2]);
285 [ + - ][ + - ]: 13117500 : G(N[0],c*3+j) = G(N[0],c*3+j) - f[0] + f[2];
286 [ + - ][ + - ]: 13117500 : G(N[1],c*3+j) = G(N[1],c*3+j) + f[0] - f[1];
287 [ + - ][ + - ]: 13117500 : G(N[2],c*3+j) = G(N[2],c*3+j) + f[1] - f[2];
288 : : }
289 : : }
290 : : }
291 : :
292 : : // domain edge contributions: edges
293 [ + + ]: 4353180 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
294 : 4342830 : const auto N = dsupedge[2].data() + e*2;
295 : 4342830 : const auto d = dsupint[2].data() + e*3;
296 [ + + ]: 26056980 : for (std::size_t c=0; c<ncomp; ++c) {
297 [ + - ][ + - ]: 21714150 : tk::real u[] = { U(N[0],c), U(N[1],c) };
298 [ + + ]: 86856600 : for (std::size_t j=0; j<3; ++j) {
299 : 65142450 : tk::real f = d[j] * (u[1] + u[0]);
300 [ + - ]: 65142450 : G(N[0],c*3+j) -= f;
301 [ + - ]: 65142450 : G(N[1],c*3+j) += f;
302 : : }
303 : : }
304 : : }
305 : :
306 : : // boundary integral
307 : :
308 : 10350 : const auto& x = coord[0];
309 : 10350 : const auto& y = coord[1];
310 : 10350 : const auto& z = coord[2];
311 : :
312 [ + + ]: 1770990 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
313 : 1760640 : const auto N = triinpoel.data() + e*3;
314 : : const std::array< tk::real, 3 >
315 : 1760640 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
316 : 1760640 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
317 : 1760640 : auto n = tk::cross( ba, ca );
318 : 1760640 : n[0] /= 12.0;
319 : 1760640 : n[1] /= 12.0;
320 : 1760640 : n[2] /= 12.0;
321 [ + + ]: 10563840 : for (std::size_t c=0; c<ncomp; ++c) {
322 [ + - ][ + - ]: 8803200 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c) };
[ + - ]
323 : 8803200 : auto uab = (u[0] + u[1])/4.0;
324 : 8803200 : auto ubc = (u[1] + u[2])/4.0;
325 : 8803200 : auto uca = (u[2] + u[0])/4.0;
326 : 8803200 : tk::real g[] = { uab + uca + u[0],
327 : 8803200 : uab + ubc + u[1],
328 : 8803200 : ubc + uca + u[2] };
329 [ + + ]: 35212800 : for (std::size_t j=0; j<3; ++j) {
330 [ + - ]: 26409600 : G(N[0],c*3+j) += g[j] * n[j];
331 [ + - ]: 26409600 : G(N[1],c*3+j) += g[j] * n[j];
332 [ + - ]: 26409600 : G(N[2],c*3+j) += g[j] * n[j];
333 : : }
334 : : }
335 : : }
336 : :
337 : : #if defined(__clang__)
338 : : #pragma clang diagnostic pop
339 : : #elif defined(STRICT_GNUC)
340 : : #pragma GCC diagnostic pop
341 : : #endif
342 : 10350 : }
343 : :
344 : : tk::real
345 : 43363680 : refvel( tk::real r, tk::real p, tk::real v )
346 : : // *****************************************************************************
347 : : // Compute reference velocitity of the preconditioned system
348 : : //! \param[in] r Density
349 : : //! \param[in] p Pressure
350 : : //! \param[in] v Velocity magnitude
351 : : //! \return Reference velocity
352 : : // *****************************************************************************
353 : : {
354 : 43363680 : auto K = g_cfg.get< tag::turkel >();
355 [ + - ]: 43363680 : auto velinf = g_cfg.get< tag::velinf >();
356 : 43363680 : auto vinf = tk::length( velinf[0], velinf[1], velinf[2] );
357 : :
358 : 86727360 : return std::min( eos::soundspeed( r, p ), std::max( v, K*vinf ) );
359 : : //return eos::soundspeed( r, p );
360 : 43363680 : }
361 : :
362 : : static std::tuple< tk::real, tk::real >
363 : 39013140 : sigvel( tk::real p, tk::real T, tk::real v, tk::real vn )
364 : : // *****************************************************************************
365 : : // Compute signal velocitities of the preconditioned system
366 : : //! \param[in] p Pressure
367 : : //! \param[in] T Temperature
368 : : //! \param[in] v Velocity magnitude
369 : : //! \param[in] vn Face-normal velocity
370 : : //! \return v', c'
371 : : // *****************************************************************************
372 : : {
373 : 39013140 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
374 : 39013140 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
375 : 39013140 : auto cp = g*rgas/(g-1.0);
376 : :
377 : 39013140 : auto r = p/T/rgas;
378 : 39013140 : auto rp = r/p;
379 : 39013140 : auto rt = -r/T;
380 [ + - ]: 39013140 : auto vr = refvel( r, p, v );
381 : 39013140 : auto vr2 = vr*vr;
382 : 39013140 : auto beta = rp + rt/r/cp;
383 : 39013140 : auto alpha = 0.5*(1.0 - beta*vr2);
384 : 39013140 : auto vpri = vn*(1.0 - alpha);
385 : 39013140 : auto cpri = std::sqrt( alpha*alpha*vn*vn + vr2 );
386 : :
387 : 78026280 : return { vpri, cpri };
388 : : }
389 : :
390 : : static void
391 : 10333320 : rusanov( const tk::UnsMesh::Coords& coord,
392 : : const tk::Fields& G,
393 : : const tk::real dsupint[],
394 : : std::size_t p,
395 : : std::size_t q,
396 : : const tk::real L[],
397 : : const tk::real R[],
398 : : tk::real f[] )
399 : : // *****************************************************************************
400 : : //! Compute advection fluxes on a single edge with Rusanov's flux
401 : : //! \param[in] coord Mesh node coordinates
402 : : //! \param[in] G Nodal gradients
403 : : //! \param[in] dsupint Domain superedge integral for this edge
404 : : //! \param[in] p Left node index of edge
405 : : //! \param[in] q Right node index of edge
406 : : //! \param[in,out] L Left physics state variables
407 : : //! \param[in,out] R Rigth physics state variables
408 : : //! \param[in,out] f Flux computed
409 : : // *****************************************************************************
410 : : {
411 : 20666640 : auto ncomp = G.nprop() / 3;
412 : :
413 : : #if defined(__clang__)
414 : : #pragma clang diagnostic push
415 : : #pragma clang diagnostic ignored "-Wvla"
416 : : #pragma clang diagnostic ignored "-Wvla-extension"
417 : : #elif defined(STRICT_GNUC)
418 : : #pragma GCC diagnostic push
419 : : #pragma GCC diagnostic ignored "-Wvla"
420 : : #endif
421 : :
422 : : // will work on copies of physics variables
423 : 10333320 : tk::real l[ncomp], r[ncomp];
424 : 10333320 : memcpy( l, L, sizeof l );
425 : 10333320 : memcpy( r, R, sizeof r );
426 : :
427 : : // MUSCL reconstruction in edge-end points for flow variables
428 : 10333320 : muscl( p, q, coord, G, l[0], l[1], l[2], l[3], l[4],
429 [ + - ]: 10333320 : r[0], r[1], r[2], r[3], r[4] );
430 : :
431 : : // dualface-normal velocities
432 : 10333320 : auto nx = dsupint[0];
433 : 10333320 : auto ny = dsupint[1];
434 : 10333320 : auto nz = dsupint[2];
435 : 10333320 : auto vnL = l[1]*nx + l[2]*ny + l[3]*nz;
436 : 10333320 : auto vnR = r[1]*nx + r[2]*ny + r[3]*nz;
437 : :
438 : : // pressure
439 : 10333320 : auto pL = l[0];
440 : 10333320 : auto pR = r[0];
441 : :
442 : : // preconditioned wave speed
443 : 10333320 : auto len = tk::length( nx, ny, nz );
444 [ + - ]: 10333320 : auto [vpL, cpL] = sigvel( l[0], l[4], tk::length(l[1],l[2],l[3]), vnL );
445 [ + - ]: 10333320 : auto [vpR, cpR] = sigvel( r[0], r[4], tk::length(r[1],r[2],r[3]), vnR );
446 : :
447 : : // convert to conserved variables
448 : 10333320 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
449 : 10333320 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
450 : 10333320 : l[0] = pL/l[4]/rgas;
451 : 10333320 : l[1] *= l[0];
452 : 10333320 : l[2] *= l[0];
453 : 10333320 : l[3] *= l[0];
454 : 10333320 : l[4] = pL/(g-1.0) + 0.5*(l[1]*l[1] + l[2]*l[2] + l[3]*l[3])/l[0];
455 : 10333320 : r[0] = pR/r[4]/rgas;
456 : 10333320 : r[1] *= r[0];
457 : 10333320 : r[2] *= r[0];
458 : 10333320 : r[3] *= r[0];
459 : 10333320 : r[4] = pR/(g-1.0) + 0.5*(r[1]*r[1] + r[2]*r[2] + r[3]*r[3])/r[0];
460 : :
461 : : // dissipation: option 1
462 : : //auto sL = std::abs(vpL) + cpL;
463 : : //auto sR = std::abs(vpR) + cpR;
464 : : //auto fw = std::max( sL, sR ) * len;
465 : :
466 : : // dissipation: option 2
467 : : using std::abs;
468 : : using std::max;
469 : 10333320 : auto sp = max(abs(vpL-cpL),max(abs(vpR-cpR),max(abs(vpL+cpL),abs(vpR+cpR))));
470 : 10333320 : auto sL = -sp;
471 : 10333320 : auto sR = +sp;
472 : 10333320 : auto fw = std::max( sL, sR ) * len;
473 : :
474 : : // flow fluxes
475 : 10333320 : f[0] = l[0]*vnL + r[0]*vnR + fw*(r[0] - l[0]);
476 : 10333320 : f[1] = l[1]*vnL + r[1]*vnR + (pL + pR)*nx + fw*(r[1] - l[1]);
477 : 10333320 : f[2] = l[2]*vnL + r[2]*vnR + (pL + pR)*ny + fw*(r[2] - l[2]);
478 : 10333320 : f[3] = l[3]*vnL + r[3]*vnR + (pL + pR)*nz + fw*(r[3] - l[3]);
479 : 10333320 : f[4] = (l[4] + pL)*vnL + (r[4] + pR)*vnR + fw*(r[4] - l[4]);
480 : :
481 : : // artificial viscosity
482 : 10333320 : const auto stab2 = g_cfg.get< tag::stab2 >();
483 [ - + ]: 10333320 : if (stab2) {
484 : 0 : auto stab2coef = g_cfg.get< tag::stab2coef >();
485 : 0 : auto fws = stab2coef * fw;
486 : 0 : f[0] -= fws*(l[0] - r[0]);
487 : 0 : f[1] -= fws*(l[1] - r[1]);
488 : 0 : f[2] -= fws*(l[2] - r[2]);
489 : 0 : f[3] -= fws*(l[3] - r[3]);
490 : 0 : f[4] -= fws*(l[4] - r[4]);
491 : : }
492 : :
493 [ + - ]: 10333320 : if (ncomp == 5) return;
494 : :
495 : : // MUSCL reconstruction in edge-end points for scalars
496 [ - - ]: 0 : muscl( p, q, coord, G, l, r );
497 : :
498 : : // scalar dissipation
499 : 0 : auto sw = std::max( std::abs(vnL), std::abs(vnR) );
500 : :
501 : : // scalar fluxes
502 [ - - ]: 0 : for (std::size_t c=5; c<ncomp; ++c) {
503 : 0 : f[c] = l[c]*vnL + r[c]*vnR + sw*(r[c] - l[c]);
504 : : }
505 : :
506 : : #if defined(__clang__)
507 : : #pragma clang diagnostic pop
508 : : #elif defined(STRICT_GNUC)
509 : : #pragma GCC diagnostic pop
510 : : #endif
511 : 10333320 : }
512 : :
513 : : static void
514 : 9173250 : hllc( const tk::UnsMesh::Coords& coord,
515 : : const tk::Fields& G,
516 : : const tk::real dsupint[],
517 : : std::size_t p,
518 : : std::size_t q,
519 : : const tk::real L[],
520 : : const tk::real R[],
521 : : tk::real f[] )
522 : : // *****************************************************************************
523 : : //! Compute advection fluxes on a single edge with Harten-Lax-vanLeer-Contact
524 : : //! \param[in] coord Mesh node coordinates
525 : : //! \param[in] G Nodal gradients
526 : : //! \param[in] dsupint Domain superedge integral for this edge
527 : : //! \param[in] p Left node index of edge
528 : : //! \param[in] q Right node index of edge
529 : : //! \param[in,out] L Left physics state variables
530 : : //! \param[in,out] R Rigth physics state variables
531 : : //! \param[in,out] f Flux computed
532 : : //! \see Toro, Riemann Solver and Numerical Methods for Fluid Dynamics. 3rd
533 : : //! Edition, Springer, New York, 2009
534 : : //! \see Davis, Simplified Second-Order Godunov-Type Methods, SIAM J. Sci. Stat.
535 : : //! Comput. 9:445-473, 1988
536 : : //! \see Roe, Approximate Riemann Solvers, Parameter Vectors, and Difference
537 : : //! Schemes, J. Comput Phys., 43:357-372, 1981
538 : : // *****************************************************************************
539 : : {
540 : 18346500 : auto ncomp = G.nprop() / 3;
541 : :
542 : : #if defined(__clang__)
543 : : #pragma clang diagnostic push
544 : : #pragma clang diagnostic ignored "-Wvla"
545 : : #pragma clang diagnostic ignored "-Wvla-extension"
546 : : #elif defined(STRICT_GNUC)
547 : : #pragma GCC diagnostic push
548 : : #pragma GCC diagnostic ignored "-Wvla"
549 : : #endif
550 : :
551 : : // will work on copies of physics variables
552 : 9173250 : tk::real l[ncomp], r[ncomp];
553 : 9173250 : memcpy( l, L, sizeof l );
554 : 9173250 : memcpy( r, R, sizeof r );
555 : :
556 : : // MUSCL reconstruction in edge-end points for flow variables
557 : 9173250 : muscl( p, q, coord, G, l[0], l[1], l[2], l[3], l[4],
558 [ + - ]: 9173250 : r[0], r[1], r[2], r[3], r[4] );
559 : :
560 : : // dualface-normal velocities
561 : 9173250 : auto nx = -dsupint[0];
562 : 9173250 : auto ny = -dsupint[1];
563 : 9173250 : auto nz = -dsupint[2];
564 : 9173250 : auto len = tk::length( nx, ny, nz );
565 : 9173250 : nx /= len;
566 : 9173250 : ny /= len;
567 : 9173250 : nz /= len;
568 : 9173250 : auto qL = l[1]*nx + l[2]*ny + l[3]*nz;
569 : 9173250 : auto qR = r[1]*nx + r[2]*ny + r[3]*nz;
570 : :
571 : : // pressure
572 : 9173250 : auto pL = l[0];
573 : 9173250 : auto pR = r[0];
574 : :
575 : : // preconditioned signal velocities
576 [ + - ]: 9173250 : auto [vpL, cpL] = sigvel( l[0], l[4], tk::length(l[1],l[2],l[3]), qL*len );
577 [ + - ]: 9173250 : auto [vpR, cpR] = sigvel( r[0], r[4], tk::length(r[1],r[2],r[3]), qR*len );
578 : :
579 : : // convert to conserved variables
580 : 9173250 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
581 : 9173250 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
582 : 9173250 : l[0] = pL/l[4]/rgas;
583 : 9173250 : l[1] *= l[0];
584 : 9173250 : l[2] *= l[0];
585 : 9173250 : l[3] *= l[0];
586 : 9173250 : l[4] = pL/(g-1.0) + 0.5*(l[1]*l[1] + l[2]*l[2] + l[3]*l[3])/l[0];
587 : 9173250 : r[0] = pR/r[4]/rgas;
588 : 9173250 : r[1] *= r[0];
589 : 9173250 : r[2] *= r[0];
590 : 9173250 : r[3] *= r[0];
591 : 9173250 : r[4] = pR/(g-1.0) + 0.5*(r[1]*r[1] + r[2]*r[2] + r[3]*r[3])/r[0];
592 : :
593 : : // preconditioned left and right wave speeds
594 : :
595 : : // option 1
596 : : //auto sL = std::min( vpL - cpL, vpR - cpR );
597 : : //auto sR = std::max( vpL + cpL, vpR + cpR );
598 : :
599 : : // option 2
600 : : //auto sL = vpL - cpL;
601 : : //auto sR = vpR + cpR;
602 : :
603 : : // option 3
604 : : using std::abs;
605 : : using std::max;
606 : 9173250 : auto sp = max(abs(vpL-cpL),max(abs(vpR-cpR),max(abs(vpL+cpL),abs(vpR+cpR))));
607 : 9173250 : auto sL = -sp;
608 : 9173250 : auto sR = +sp;
609 : :
610 : : //using std::sqrt;
611 : : //auto hl = (l[4]/l[0] + pL) / l[0];
612 : : //auto hr = (r[4]/r[0] + pR) / r[0];
613 : : //auto srl = sqrt( l[0] );
614 : : //auto srr = sqrt( r[0] );
615 : : //auto hh = (srl*hl + srr*hr) / (srl + srr);
616 : : //auto uh = (srl*vpL + srr*vpR) / (srl + srr);
617 : : //auto ch = sqrt( (g-1.0)*(hh - 0.5*(uh*uh)) );
618 : : //// option 4
619 : : //auto sL = std::min( vpL - cpL, uh - ch );
620 : : //auto sR = std::max( vpR + cpR, uh + ch );
621 : : // option 5
622 : : //auto sL = uh - ch;
623 : : //auto sR = uh + ch;
624 : :
625 : : // no preconditioning
626 : : //auto cL = eos::soundspeed(l[0],pL);
627 : : //auto cR = eos::soundspeed(r[0],pR);
628 : : //auto sL = fmin( qL - cL, qR - cR );
629 : : //auto sR = fmax( qL + cL, qR + cR );
630 : :
631 : : // contact wave speed and pressure
632 : 9173250 : auto tL = sL - qL;
633 : 9173250 : auto tR = sR - qR;
634 : 9173250 : auto sM = (r[0]*qR*tR - l[0]*qL*tL + pL - pR) / (r[0]*tR - l[0]*tL);
635 : 9173250 : auto pS = pL - l[0]*tL*(qL - sM);
636 : :
637 : : // intermediate left-, and right-state conserved unknowns
638 : 9173250 : tk::real uL[ncomp], uR[ncomp];
639 : 9173250 : auto s = sL - sM;
640 : 9173250 : uL[0] = tL*l[0]/s;
641 : 9173250 : uL[1] = (tL*l[1] + (pS-pL)*nx)/s;
642 : 9173250 : uL[2] = (tL*l[2] + (pS-pL)*ny)/s;
643 : 9173250 : uL[3] = (tL*l[3] + (pS-pL)*nz)/s;
644 : 9173250 : uL[4] = (tL*l[4] - pL*qL + pS*sM)/s;
645 : 9173250 : s = sR - sM;
646 : 9173250 : uR[0] = tR*r[0]/s;
647 : 9173250 : uR[1] = (tR*r[1] + (pS-pR)*nx)/s;
648 : 9173250 : uR[2] = (tR*r[2] + (pS-pR)*ny)/s;
649 : 9173250 : uR[3] = (tR*r[3] + (pS-pR)*nz)/s;
650 : 9173250 : uR[4] = (tR*r[4] - pR*qR + pS*sM)/s;
651 : :
652 : 9173250 : auto L2 = -2.0*len;
653 : 9173250 : nx *= L2;
654 : 9173250 : ny *= L2;
655 : 9173250 : nz *= L2;
656 : :
657 : : // flow fluxes
658 [ - + ]: 9173250 : if (sL > 0.0) {
659 : 0 : auto qL2 = qL * L2;
660 : 0 : f[0] = l[0]*qL2;
661 : 0 : f[1] = l[1]*qL2 + pL*nx;
662 : 0 : f[2] = l[2]*qL2 + pL*ny;
663 : 0 : f[3] = l[3]*qL2 + pL*nz;
664 : 0 : f[4] = (l[4] + pL)*qL2;
665 : : }
666 [ + - ][ + + ]: 9173250 : else if (sL <= 0.0 && sM > 0.0) {
667 : 4980167 : auto qL2 = qL * L2;
668 : 4980167 : auto sL2 = sL * L2;
669 : 4980167 : f[0] = l[0]*qL2 + sL2*(uL[0] - l[0]);
670 : 4980167 : f[1] = l[1]*qL2 + pL*nx + sL2*(uL[1] - l[1]);
671 : 4980167 : f[2] = l[2]*qL2 + pL*ny + sL2*(uL[2] - l[2]);
672 : 4980167 : f[3] = l[3]*qL2 + pL*nz + sL2*(uL[3] - l[3]);
673 : 4980167 : f[4] = (l[4] + pL)*qL2 + sL2*(uL[4] - l[4]);
674 : 4980167 : }
675 [ + - ][ + - ]: 4193083 : else if (sM <= 0.0 && sR >= 0.0) {
676 : 4193083 : auto qR2 = qR * L2;
677 : 4193083 : auto sR2 = sR * L2;
678 : 4193083 : f[0] = r[0]*qR2 + sR2*(uR[0] - r[0]);
679 : 4193083 : f[1] = r[1]*qR2 + pR*nx + sR2*(uR[1] - r[1]);
680 : 4193083 : f[2] = r[2]*qR2 + pR*ny + sR2*(uR[2] - r[2]);
681 : 4193083 : f[3] = r[3]*qR2 + pR*nz + sR2*(uR[3] - r[3]);
682 : 4193083 : f[4] = (r[4] + pR)*qR2 + sR2*(uR[4] - r[4]);
683 : 4193083 : }
684 : : else {
685 : 0 : auto qR2 = qR * L2;
686 : 0 : f[0] = r[0]*qR2;
687 : 0 : f[1] = r[1]*qR2 + pR*nx;
688 : 0 : f[2] = r[2]*qR2 + pR*ny;
689 : 0 : f[3] = r[3]*qR2 + pR*nz;
690 : 0 : f[4] = (r[4] + pR)*qR2;
691 : : }
692 : :
693 : : // artificial viscosity
694 : : //const auto stab2 = g_cfg.get< tag::stab2 >();
695 : : //if (stab2) {
696 : : // auto stab2coef = g_cfg.get< tag::stab2coef >();
697 : : // auto sl = std::abs(vpL) + cpL;
698 : : // auto sr = std::abs(vpR) + cpR;
699 : : // auto fws = stab2coef * std::max(sl,sr) * len;
700 : : // f[0] -= fws*(l[0] - r[0]);
701 : : // f[1] -= fws*(l[1] - r[1]);
702 : : // f[2] -= fws*(l[2] - r[2]);
703 : : // f[3] -= fws*(l[3] - r[3]);
704 : : // f[4] -= fws*(l[4] - r[4]);
705 : : //}
706 : :
707 [ + - ]: 9173250 : if (ncomp == 5) return;
708 : :
709 : : // MUSCL reconstruction in edge-end points for scalars
710 [ - - ]: 0 : muscl( p, q, coord, G, l, r );
711 : :
712 : : // scalar fluxes
713 : : //auto sw = std::max( std::abs(vpL), std::abs(vpR) ) * len;
714 : : //for (std::size_t c=5; c<ncomp; ++c) {
715 : : // f[c] = (l[c]*qL + r[c]*qR)*len + sw*(r[c] - l[c]);
716 : : //}
717 : :
718 : : #if defined(__clang__)
719 : : #pragma clang diagnostic pop
720 : : #elif defined(STRICT_GNUC)
721 : : #pragma GCC diagnostic pop
722 : : #endif
723 : 9173250 : }
724 : :
725 : : static void
726 : 10350 : advdom( const tk::UnsMesh::Coords& coord,
727 : : const std::array< std::vector< std::size_t >, 3 >& dsupedge,
728 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
729 : : const tk::Fields& G,
730 : : const tk::Fields& U,
731 : : // cppcheck-suppress constParameter
732 : : tk::Fields& R )
733 : : // *****************************************************************************
734 : : //! Compute domain integral for advection
735 : : //! \param[in] coord Mesh node coordinates
736 : : //! \param[in] dsupedge Domain superedges
737 : : //! \param[in] dsupint Domain superedge integrals
738 : : //! \param[in] G Nodal gradients
739 : : //! \param[in] U Solution vector of primitive variables at recent time step
740 : : //! \param[in,out] R Right-hand side vector computed
741 : : // *****************************************************************************
742 : : {
743 : : // number of transported scalars
744 : 10350 : auto ncomp = U.nprop();
745 : :
746 : : // configure advection flux function
747 : 0 : auto adv = [](){
748 : 10350 : const auto& flux = g_cfg.get< tag::flux >();
749 [ + + ]: 10350 : if (flux == "rusanov") return rusanov;
750 [ + - ]: 390 : else if (flux == "hllc") return hllc;
751 [ - - ][ - - ]: 0 : else Throw( "Flux not correctly configured" );
[ - - ]
752 [ + - ]: 10350 : }();
753 : :
754 : : #if defined(__clang__)
755 : : #pragma clang diagnostic push
756 : : #pragma clang diagnostic ignored "-Wvla"
757 : : #pragma clang diagnostic ignored "-Wvla-extension"
758 : : #elif defined(STRICT_GNUC)
759 : : #pragma GCC diagnostic push
760 : : #pragma GCC diagnostic ignored "-Wvla"
761 : : #endif
762 : :
763 : : // domain edge contributions: tetrahedron superedges
764 [ + + ]: 2100390 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
765 : 2090040 : const auto N = dsupedge[0].data() + e*4;
766 : 2090040 : tk::real u[4][ncomp];
767 [ + + ]: 12540240 : for (std::size_t c=0; c<ncomp; ++c) {
768 [ + - ]: 10450200 : u[0][c] = U(N[0],c);
769 [ + - ]: 10450200 : u[1][c] = U(N[1],c);
770 [ + - ]: 10450200 : u[2][c] = U(N[2],c);
771 [ + - ]: 10450200 : u[3][c] = U(N[3],c);
772 : : }
773 : : // edge fluxes
774 : 2090040 : tk::real f[6][ncomp];
775 : 2090040 : const auto d = dsupint[0].data();
776 [ + - ]: 2090040 : adv( coord, G, d+(e*6+0)*3, N[0], N[1], u[0], u[1], f[0] );
777 [ + - ]: 2090040 : adv( coord, G, d+(e*6+1)*3, N[1], N[2], u[1], u[2], f[1] );
778 [ + - ]: 2090040 : adv( coord, G, d+(e*6+2)*3, N[2], N[0], u[2], u[0], f[2] );
779 [ + - ]: 2090040 : adv( coord, G, d+(e*6+3)*3, N[0], N[3], u[0], u[3], f[3] );
780 [ + - ]: 2090040 : adv( coord, G, d+(e*6+4)*3, N[1], N[3], u[1], u[3], f[4] );
781 [ + - ]: 2090040 : adv( coord, G, d+(e*6+5)*3, N[2], N[3], u[2], u[3], f[5] );
782 : : // edge flux contributions
783 [ + + ]: 12540240 : for (std::size_t c=0; c<ncomp; ++c) {
784 [ + - ][ + - ]: 10450200 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
785 [ + - ][ + - ]: 10450200 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
786 [ + - ][ + - ]: 10450200 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
787 [ + - ][ + - ]: 10450200 : R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
788 : : }
789 : 2090040 : }
790 : :
791 : : // domain edge contributions: triangle superedges
792 [ + + ]: 884850 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
793 : 874500 : const auto N = dsupedge[1].data() + e*3;
794 : 874500 : tk::real u[3][ncomp];
795 [ + + ]: 5247000 : for (std::size_t c=0; c<ncomp; ++c) {
796 [ + - ]: 4372500 : u[0][c] = U(N[0],c);
797 [ + - ]: 4372500 : u[1][c] = U(N[1],c);
798 [ + - ]: 4372500 : u[2][c] = U(N[2],c);
799 : : }
800 : : // edge fluxes
801 : 874500 : tk::real f[3][ncomp];
802 : 874500 : const auto d = dsupint[1].data();
803 [ + - ]: 874500 : adv( coord, G, d+(e*3+0)*3, N[0], N[1], u[0], u[1], f[0] );
804 [ + - ]: 874500 : adv( coord, G, d+(e*3+1)*3, N[1], N[2], u[1], u[2], f[1] );
805 [ + - ]: 874500 : adv( coord, G, d+(e*3+2)*3, N[2], N[0], u[2], u[0], f[2] );
806 : : // edge flux contributions
807 [ + + ]: 5247000 : for (std::size_t c=0; c<ncomp; ++c) {
808 [ + - ][ + - ]: 4372500 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
809 [ + - ][ + - ]: 4372500 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
810 [ + - ][ + - ]: 4372500 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
811 : : }
812 : 874500 : }
813 : :
814 : : // domain edge contributions: edges
815 [ + + ]: 4353180 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
816 : 4342830 : const auto N = dsupedge[2].data() + e*2;
817 : 4342830 : tk::real u[2][ncomp];
818 [ + + ]: 26056980 : for (std::size_t c=0; c<ncomp; ++c) {
819 [ + - ]: 21714150 : u[0][c] = U(N[0],c);
820 [ + - ]: 21714150 : u[1][c] = U(N[1],c);
821 : : }
822 : : // edge fluxes
823 : 4342830 : tk::real f[ncomp];
824 : 4342830 : const auto d = dsupint[2].data();
825 [ + - ]: 4342830 : adv( coord, G, d+e*3, N[0], N[1], u[0], u[1], f );
826 : : // edge flux contributions
827 [ + + ]: 26056980 : for (std::size_t c=0; c<ncomp; ++c) {
828 [ + - ]: 21714150 : R(N[0],c) -= f[c];
829 [ + - ]: 21714150 : R(N[1],c) += f[c];
830 : : }
831 : 4342830 : }
832 : :
833 : : #if defined(__clang__)
834 : : #pragma clang diagnostic pop
835 : : #elif defined(STRICT_GNUC)
836 : : #pragma GCC diagnostic pop
837 : : #endif
838 : 10350 : }
839 : :
840 : : static void
841 : 10350 : advbnd( const std::vector< std::size_t >& triinpoel,
842 : : const std::array< std::vector< tk::real >, 3 >& coord,
843 : : const std::vector< std::uint8_t >& besym,
844 : : const tk::Fields& U,
845 : : tk::Fields& R )
846 : : // *****************************************************************************
847 : : //! Compute boundary integral for advection
848 : : //! \param[in] triinpoel Boundary face connectivity
849 : : //! \param[in] coord Mesh node coordinates
850 : : //! \param[in] besym Boundary element symmetry BC flags
851 : : //! \param[in] U Solution vector at recent time step
852 : : //! \param[in,out] R Right-hand side vector
853 : : // *****************************************************************************
854 : : {
855 : 10350 : auto ncomp = U.nprop();
856 : :
857 : 10350 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
858 : 10350 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
859 : :
860 : 10350 : const auto& x = coord[0];
861 : 10350 : const auto& y = coord[1];
862 : 10350 : const auto& z = coord[2];
863 : :
864 : : #if defined(__clang__)
865 : : #pragma clang diagnostic push
866 : : #pragma clang diagnostic ignored "-Wvla"
867 : : #pragma clang diagnostic ignored "-Wvla-extension"
868 : : #elif defined(STRICT_GNUC)
869 : : #pragma GCC diagnostic push
870 : : #pragma GCC diagnostic ignored "-Wvla"
871 : : #endif
872 : :
873 [ + + ]: 1770990 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
874 : 3521280 : const auto N = triinpoel.data() + e*3;
875 : :
876 [ + - ][ + - ]: 1760640 : auto rA = U(N[0],0)/U(N[0],4)/rgas;
877 [ + - ]: 1760640 : auto ruA = U(N[0],1) * rA;
878 [ + - ]: 1760640 : auto rvA = U(N[0],2) * rA;
879 [ + - ]: 1760640 : auto rwA = U(N[0],3) * rA;
880 [ + - ]: 1760640 : auto reA = U(N[0],0)/(g-1.0) + 0.5*(ruA*ruA + rvA*rvA + rwA*rwA)/rA;
881 : :
882 [ + - ][ + - ]: 1760640 : auto rB = U(N[1],0)/U(N[1],4)/rgas;
883 [ + - ]: 1760640 : auto ruB = U(N[1],1) * rB;
884 [ + - ]: 1760640 : auto rvB = U(N[1],2) * rB;
885 [ + - ]: 1760640 : auto rwB = U(N[1],3) * rB;
886 [ + - ]: 1760640 : auto reB = U(N[1],0)/(g-1.0) + 0.5*(ruB*ruB + rvB*rvB + rwB*rwB)/rB;
887 : :
888 [ + - ][ + - ]: 1760640 : auto rC = U(N[2],0)/U(N[2],4)/rgas;
889 [ + - ]: 1760640 : auto ruC = U(N[2],1) * rC;
890 [ + - ]: 1760640 : auto rvC = U(N[2],2) * rC;
891 [ + - ]: 1760640 : auto rwC = U(N[2],3) * rC;
892 [ + - ]: 1760640 : auto reC = U(N[2],0)/(g-1.0) + 0.5*(ruC*ruC + rvC*rvC + rwC*rwC)/rC;
893 : :
894 : : const std::array< tk::real, 3 >
895 : 1760640 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
896 : 1760640 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
897 : 1760640 : auto [nx,ny,nz] = tk::cross( ba, ca ); // 2A
898 : 1760640 : nx /= 12.0;
899 : 1760640 : ny /= 12.0;
900 : 1760640 : nz /= 12.0;
901 : :
902 : 1760640 : tk::real vn, f[ncomp][3];
903 : 1760640 : const auto sym = besym.data() + e*3;
904 : :
905 [ + + ][ + - ]: 1760640 : vn = sym[0] ? 0.0 : (nx*U(N[0],1) + ny*U(N[0],2) + nz*U(N[0],3));
[ + - ][ + - ]
906 : : // flow
907 : 1760640 : f[0][0] = rA*vn;
908 [ + - ]: 1760640 : f[1][0] = ruA*vn + U(N[0],0)*nx;
909 [ + - ]: 1760640 : f[2][0] = rvA*vn + U(N[0],0)*ny;
910 [ + - ]: 1760640 : f[3][0] = rwA*vn + U(N[0],0)*nz;
911 [ + - ]: 1760640 : f[4][0] = (reA + U(N[0],0))*vn;
912 : : // scalar
913 [ - - ][ - + ]: 1760640 : for (std::size_t c=5; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;
914 : :
915 [ + + ][ + - ]: 1760640 : vn = sym[1] ? 0.0 : (nx*U(N[1],1) + ny*U(N[1],2) + nz*U(N[1],3));
[ + - ][ + - ]
916 : : // flow
917 : 1760640 : f[0][1] = rB*vn;
918 [ + - ]: 1760640 : f[1][1] = ruB*vn + U(N[1],0)*nx;
919 [ + - ]: 1760640 : f[2][1] = rvB*vn + U(N[1],0)*ny;
920 [ + - ]: 1760640 : f[3][1] = rwB*vn + U(N[1],0)*nz;
921 [ + - ]: 1760640 : f[4][1] = (reB + U(N[1],0))*vn;
922 : : // scalar
923 [ - - ][ - + ]: 1760640 : for (std::size_t c=5; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;
924 : :
925 [ + + ][ + - ]: 1760640 : vn = sym[2] ? 0.0 : (nx*U(N[2],1) + ny*U(N[2],2) + nz*U(N[2],3));
[ + - ][ + - ]
926 : : // flow
927 : 1760640 : f[0][2] = rC*vn;
928 [ + - ]: 1760640 : f[1][2] = ruC*vn + U(N[2],0)*nx;
929 [ + - ]: 1760640 : f[2][2] = rvC*vn + U(N[2],0)*ny;
930 [ + - ]: 1760640 : f[3][2] = rwC*vn + U(N[2],0)*nz;
931 [ + - ]: 1760640 : f[4][2] = (reC + U(N[2],0))*vn;
932 : : // scalar
933 [ - - ][ - + ]: 1760640 : for (std::size_t c=5; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;
934 : :
935 [ + + ]: 10563840 : for (std::size_t c=0; c<ncomp; ++c) {
936 : 8803200 : auto fab = (f[c][0] + f[c][1])/4.0;
937 : 8803200 : auto fbc = (f[c][1] + f[c][2])/4.0;
938 : 8803200 : auto fca = (f[c][2] + f[c][0])/4.0;
939 [ + - ]: 8803200 : R(N[0],c) += fab + fca + f[c][0];
940 [ + - ]: 8803200 : R(N[1],c) += fab + fbc + f[c][1];
941 [ + - ]: 8803200 : R(N[2],c) += fbc + fca + f[c][2];
942 : : }
943 : 1760640 : }
944 : :
945 : : #if defined(__clang__)
946 : : #pragma clang diagnostic pop
947 : : #elif defined(STRICT_GNUC)
948 : : #pragma GCC diagnostic pop
949 : : #endif
950 : 10350 : }
951 : :
952 : : static void
953 : 10350 : src( const std::array< std::vector< tk::real >, 3 >& coord,
954 : : const std::vector< tk::real >& v,
955 : : tk::real t,
956 : : const std::vector< tk::real >& tp,
957 : : tk::Fields& R )
958 : : // *****************************************************************************
959 : : // Compute source integral
960 : : //! \param[in] coord Mesh node coordinates
961 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
962 : : //! \param[in] t Physical time
963 : : //! \param[in] tp Physical time for each mesh node
964 : : //! \param[in,out] R Right-hand side vector computed
965 : : // *****************************************************************************
966 : : {
967 [ + - ]: 10350 : auto src = problems::SRC();
968 [ + - ]: 10350 : if (!src) return;
969 : :
970 : 0 : const auto& x = coord[0];
971 : 0 : const auto& y = coord[1];
972 : 0 : const auto& z = coord[2];
973 : :
974 [ - - ]: 0 : for (std::size_t p=0; p<R.nunk(); ++p) {
975 [ - - ]: 0 : if (g_cfg.get< tag::steady >()) t = tp[p];
976 [ - - ]: 0 : auto s = src( x[p], y[p], z[p], t );
977 [ - - ][ - - ]: 0 : for (std::size_t c=0; c<s.size(); ++c) R(p,c) -= s[c] * v[p];
978 : 0 : }
979 [ - + ]: 10350 : }
980 : :
981 : : void
982 : 10350 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
983 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
984 : : const std::array< std::vector< tk::real >, 3 >& coord,
985 : : const std::vector< std::size_t >& triinpoel,
986 : : const std::vector< std::uint8_t >& besym,
987 : : const tk::Fields& G,
988 : : const tk::Fields& U,
989 : : const std::vector< tk::real >& v,
990 : : tk::real t,
991 : : const std::vector< tk::real >& tp,
992 : : tk::Fields& R )
993 : : // *****************************************************************************
994 : : // Compute right hand side
995 : : //! \param[in] dsupedge Domain superedges
996 : : //! \param[in] dsupint Domain superedge integrals
997 : : //! \param[in] coord Mesh node coordinates
998 : : //! \param[in] triinpoel Boundary face connectivity
999 : : //! \param[in] besym Boundary element symmetry BC flags
1000 : : //! \param[in] G Gradients in mesh nodes
1001 : : //! \param[in] U Solution vector of primitive variables at recent time step
1002 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
1003 : : //! \param[in] t Physical time
1004 : : //! \param[in] tp Physical time for each mesh node
1005 : : //! \param[in,out] R Right-hand side vector computed
1006 : : // *****************************************************************************
1007 : : {
1008 [ - + ][ - - ]: 10350 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
1009 : : "vector at recent time step incorrect" );
1010 [ - + ][ - - ]: 10350 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
1011 : : "Number of unknowns and/or number of components in right-hand "
1012 : : "side vector incorrect" );
1013 : :
1014 : 10350 : R.fill( 0.0 );
1015 : 10350 : advdom( coord, dsupedge, dsupint, G, U, R );
1016 : 10350 : advbnd( triinpoel, coord, besym, U, R );
1017 : 10350 : src( coord, v, t, tp, R );
1018 : 10350 : }
1019 : :
1020 : : } // lax::
|
