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 : 14957220 : 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 : : const auto& x = coord[0];
66 : : const auto& y = coord[1];
67 : : const auto& z = coord[2];
68 : :
69 : : // edge vector
70 : 14957220 : 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 : 14957220 : tk::real ls[5] = { pL, uL, vL, wL, tL },
74 : 14957220 : rs[5] = { pR, uR, vR, wR, tR },
75 : : url[5], urr[5];
76 : 14957220 : memcpy( url, ls, sizeof ls );
77 : 14957220 : memcpy( urr, rs, sizeof rs );
78 : :
79 : : // MUSCL reconstruction of edge-end-point primitive variables
80 [ + + ]: 89743320 : for (std::size_t c=0; c<5; ++c) {
81 : :
82 : 74786100 : auto g1 = G(p,c*3+0)*vw[0] + G(p,c*3+1)*vw[1] + G(p,c*3+2)*vw[2];
83 : 74786100 : 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 : 74786100 : delta2[c] = rs[c] - ls[c];
86 : 74786100 : delta1[c] = 2.0 * g1 - delta2[c];
87 : 74786100 : 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 : 74786100 : auto rcL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
98 : 74786100 : auto rcR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
99 : 74786100 : auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
100 : 74786100 : auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
101 : 74786100 : auto phiL = (std::abs(rcL) + rcL) / (std::abs(rcL) + 1.0);
102 : 74786100 : auto phiR = (std::abs(rcR) + rcR) / (std::abs(rcR) + 1.0);
103 : 74786100 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
104 : 74786100 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
105 : : // update unknowns with reconstructed unknowns
106 : 74786100 : url[c] += 0.25*(delta1[c]*(1.0-muscl_const)*phiL +
107 : 74786100 : delta2[c]*(1.0+muscl_const)*phi_L_inv);
108 : 74786100 : urr[c] -= 0.25*(delta3[c]*(1.0-muscl_const)*phiR +
109 : 74786100 : 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 [ + - ][ - + ]: 14957220 : if (ls[0] < delta1[0] || ls[4] < delta1[4]) memcpy( url, ls, sizeof ls );
135 [ + - ][ - + ]: 14957220 : if (rs[0] < -delta3[0] || rs[4] < -delta3[4]) memcpy( urr, rs, sizeof rs );
136 : :
137 : 14957220 : pL = url[0];
138 : 14957220 : uL = url[1];
139 : 14957220 : vL = url[2];
140 : 14957220 : wL = url[3];
141 : 14957220 : tL = url[4];
142 : :
143 : 14957220 : pR = urr[0];
144 : 14957220 : uR = urr[1];
145 : 14957220 : vR = urr[2];
146 : 14957220 : wR = urr[3];
147 : 14957220 : tR = urr[4];
148 : 14957220 : }
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 : : const auto& x = coord[0];
175 : : const auto& y = coord[1];
176 : : 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 : 10200 : 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 : : auto ncomp = U.nprop();
245 : :
246 : : Assert( G.nunk() == U.nunk(), "Size mismatch" );
247 : : Assert( G.nprop() == ncomp*3, "Size mismatch" );
248 : : G.fill( 0.0 );
249 : :
250 : : // domain integral
251 : :
252 : : // domain edge contributions: tetrahedron superedges
253 [ + + ]: 1603650 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
254 : 1593450 : const auto N = dsupedge[0].data() + e*4;
255 [ + + ]: 9560700 : for (std::size_t c=0; c<ncomp; ++c) {
256 : 7967250 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c), U(N[3],c) };
257 [ + + ]: 31869000 : for (std::size_t j=0; j<3; ++j) {
258 : : tk::real f[6];
259 : : const auto d = dsupint[0].data();
260 : 23901750 : f[0] = d[(e*6+0)*3+j] * (u[1] + u[0]);
261 : 23901750 : f[1] = d[(e*6+1)*3+j] * (u[2] + u[1]);
262 : 23901750 : f[2] = d[(e*6+2)*3+j] * (u[0] + u[2]);
263 : 23901750 : f[3] = d[(e*6+3)*3+j] * (u[3] + u[0]);
264 : 23901750 : f[4] = d[(e*6+4)*3+j] * (u[3] + u[1]);
265 : 23901750 : f[5] = d[(e*6+5)*3+j] * (u[3] + u[2]);
266 : 23901750 : G(N[0],c*3+j) = G(N[0],c*3+j) - f[0] + f[2] - f[3];
267 : 23901750 : G(N[1],c*3+j) = G(N[1],c*3+j) + f[0] - f[1] - f[4];
268 : 23901750 : G(N[2],c*3+j) = G(N[2],c*3+j) + f[1] - f[2] - f[5];
269 : 23901750 : 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 [ + + ]: 693990 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
276 : 683790 : const auto N = dsupedge[1].data() + e*3;
277 [ + + ]: 4102740 : for (std::size_t c=0; c<ncomp; ++c) {
278 : 3418950 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c) };
279 [ + + ]: 13675800 : for (std::size_t j=0; j<3; ++j) {
280 : : tk::real f[3];
281 : : const auto d = dsupint[1].data();
282 : 10256850 : f[0] = d[(e*3+0)*3+j] * (u[1] + u[0]);
283 : 10256850 : f[1] = d[(e*3+1)*3+j] * (u[2] + u[1]);
284 : 10256850 : f[2] = d[(e*3+2)*3+j] * (u[0] + u[2]);
285 : 10256850 : G(N[0],c*3+j) = G(N[0],c*3+j) - f[0] + f[2];
286 : 10256850 : G(N[1],c*3+j) = G(N[1],c*3+j) + f[0] - f[1];
287 : 10256850 : 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 [ + + ]: 3355350 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
294 : 3345150 : const auto N = dsupedge[2].data() + e*2;
295 : 3345150 : const auto d = dsupint[2].data() + e*3;
296 [ + + ]: 20070900 : for (std::size_t c=0; c<ncomp; ++c) {
297 : 16725750 : tk::real u[] = { U(N[0],c), U(N[1],c) };
298 [ + + ]: 66903000 : for (std::size_t j=0; j<3; ++j) {
299 : 50177250 : tk::real f = d[j] * (u[1] + u[0]);
300 : 50177250 : G(N[0],c*3+j) -= f;
301 : 50177250 : G(N[1],c*3+j) += f;
302 : : }
303 : : }
304 : : }
305 : :
306 : : // boundary integral
307 : :
308 : : const auto& x = coord[0];
309 : : const auto& y = coord[1];
310 : : const auto& z = coord[2];
311 : :
312 [ + + ]: 1385760 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
313 : 1375560 : const auto N = triinpoel.data() + e*3;
314 : : const std::array< tk::real, 3 >
315 : 1375560 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
316 : 1375560 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
317 : : auto n = tk::cross( ba, ca );
318 : 1375560 : n[0] /= 12.0;
319 : 1375560 : n[1] /= 12.0;
320 : 1375560 : n[2] /= 12.0;
321 [ + + ]: 8253360 : for (std::size_t c=0; c<ncomp; ++c) {
322 : 6877800 : tk::real u[] = { U(N[0],c), U(N[1],c), U(N[2],c) };
323 : 6877800 : auto uab = (u[0] + u[1])/4.0;
324 : 6877800 : auto ubc = (u[1] + u[2])/4.0;
325 : 6877800 : auto uca = (u[2] + u[0])/4.0;
326 : 6877800 : tk::real g[] = { uab + uca + u[0],
327 : 6877800 : uab + ubc + u[1],
328 : 6877800 : ubc + uca + u[2] };
329 [ + + ]: 27511200 : for (std::size_t j=0; j<3; ++j) {
330 : 20633400 : G(N[0],c*3+j) += g[j] * n[j];
331 : 20633400 : G(N[1],c*3+j) += g[j] * n[j];
332 : 20633400 : 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 : 10200 : }
343 : :
344 : : tk::real
345 : 33294340 : 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 : 33294340 : auto K = g_cfg.get< tag::turkel >();
355 : 33294340 : auto velinf = g_cfg.get< tag::velinf >();
356 [ + + ]: 33294340 : auto vinf = tk::length( velinf[0], velinf[1], velinf[2] );
357 : :
358 [ + + ][ + - ]: 68474060 : return std::min( eos::soundspeed( r, p ), std::max( v, K*vinf ) );
359 : : //return eos::soundspeed( r, p );
360 : : }
361 : :
362 : : static std::tuple< tk::real, tk::real >
363 : 29914440 : 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 : 29914440 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
374 : 29914440 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
375 : 29914440 : auto cp = g*rgas/(g-1.0);
376 : :
377 : 29914440 : auto r = p/T/rgas;
378 : 29914440 : auto rp = r/p;
379 : 29914440 : auto rt = -r/T;
380 : 29914440 : auto vr = refvel( r, p, v );
381 : 29914440 : auto vr2 = vr*vr;
382 : 29914440 : auto beta = rp + rt/r/cp;
383 : 29914440 : auto alpha = 0.5*(1.0 - beta*vr2);
384 : 29914440 : auto vpri = vn*(1.0 - alpha);
385 : 29914440 : auto cpri = std::sqrt( alpha*alpha*vn*vn + vr2 );
386 : :
387 : 29914440 : 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 : 10333320 : 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 : : 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 [ + + ][ + + ]: 19342253 : 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 : 4623900 : 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 : 4623900 : 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 : 4623900 : tk::real l[ncomp], r[ncomp];
553 : 4623900 : memcpy( l, L, sizeof l );
554 : 4623900 : memcpy( r, R, sizeof r );
555 : :
556 : : // MUSCL reconstruction in edge-end points for flow variables
557 : 4623900 : muscl( p, q, coord, G, l[0], l[1], l[2], l[3], l[4],
558 : 4623900 : r[0], r[1], r[2], r[3], r[4] );
559 : :
560 : : // dualface-normal velocities
561 : 4623900 : auto nx = -dsupint[0];
562 : 4623900 : auto ny = -dsupint[1];
563 : 4623900 : auto nz = -dsupint[2];
564 : : auto len = tk::length( nx, ny, nz );
565 : 4623900 : nx /= len;
566 : 4623900 : ny /= len;
567 : 4623900 : nz /= len;
568 : 4623900 : auto qL = l[1]*nx + l[2]*ny + l[3]*nz;
569 : 4623900 : auto qR = r[1]*nx + r[2]*ny + r[3]*nz;
570 : :
571 : : // pressure
572 : 4623900 : auto pL = l[0];
573 : 4623900 : auto pR = r[0];
574 : :
575 : : // preconditioned signal velocities
576 : 4623900 : auto [vpL, cpL] = sigvel( l[0], l[4], tk::length(l[1],l[2],l[3]), qL*len );
577 : 4623900 : 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 : 4623900 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
581 : 4623900 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
582 : 4623900 : l[0] = pL/l[4]/rgas;
583 : 4623900 : l[1] *= l[0];
584 : 4623900 : l[2] *= l[0];
585 : 4623900 : l[3] *= l[0];
586 : 4623900 : l[4] = pL/(g-1.0) + 0.5*(l[1]*l[1] + l[2]*l[2] + l[3]*l[3])/l[0];
587 : 4623900 : r[0] = pR/r[4]/rgas;
588 : 4623900 : r[1] *= r[0];
589 : 4623900 : r[2] *= r[0];
590 : 4623900 : r[3] *= r[0];
591 : 4623900 : 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 [ + + ][ + + ]: 9384125 : auto sp = max(abs(vpL-cpL),max(abs(vpR-cpR),max(abs(vpL+cpL),abs(vpR+cpR))));
[ + + ]
607 : 4623900 : auto sL = -sp;
608 : : 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 : 4623900 : auto tL = sL - qL;
633 : 4623900 : auto tR = sR - qR;
634 : 4623900 : auto sM = (r[0]*qR*tR - l[0]*qL*tL + pL - pR) / (r[0]*tR - l[0]*tL);
635 : 4623900 : auto pS = pL - l[0]*tL*(qL - sM);
636 : :
637 : : // intermediate left-, and right-state conserved unknowns
638 : 4623900 : tk::real uL[ncomp], uR[ncomp];
639 : 4623900 : auto s = sL - sM;
640 : 4623900 : uL[0] = tL*l[0]/s;
641 : 4623900 : uL[1] = (tL*l[1] + (pS-pL)*nx)/s;
642 : 4623900 : uL[2] = (tL*l[2] + (pS-pL)*ny)/s;
643 : 4623900 : uL[3] = (tL*l[3] + (pS-pL)*nz)/s;
644 : 4623900 : uL[4] = (tL*l[4] - pL*qL + pS*sM)/s;
645 : 4623900 : s = sR - sM;
646 : 4623900 : uR[0] = tR*r[0]/s;
647 : 4623900 : uR[1] = (tR*r[1] + (pS-pR)*nx)/s;
648 : 4623900 : uR[2] = (tR*r[2] + (pS-pR)*ny)/s;
649 : 4623900 : uR[3] = (tR*r[3] + (pS-pR)*nz)/s;
650 : 4623900 : uR[4] = (tR*r[4] - pR*qR + pS*sM)/s;
651 : :
652 : 4623900 : auto L2 = -2.0*len;
653 : 4623900 : nx *= L2;
654 : 4623900 : ny *= L2;
655 : 4623900 : nz *= L2;
656 : :
657 : : // flow fluxes
658 [ - + ]: 4623900 : 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 [ + - ][ + + ]: 4623900 : else if (sL <= 0.0 && sM > 0.0) {
667 : 2505530 : auto qL2 = qL * L2;
668 : 2505530 : auto sL2 = sL * L2;
669 : 2505530 : f[0] = l[0]*qL2 + sL2*(uL[0] - l[0]);
670 : 2505530 : f[1] = l[1]*qL2 + pL*nx + sL2*(uL[1] - l[1]);
671 : 2505530 : f[2] = l[2]*qL2 + pL*ny + sL2*(uL[2] - l[2]);
672 : 2505530 : f[3] = l[3]*qL2 + pL*nz + sL2*(uL[3] - l[3]);
673 : 2505530 : f[4] = (l[4] + pL)*qL2 + sL2*(uL[4] - l[4]);
674 : 2505530 : }
675 [ + - ][ + - ]: 2118370 : else if (sM <= 0.0 && sR >= 0.0) {
676 : 2118370 : auto qR2 = qR * L2;
677 : 2118370 : auto sR2 = sR * L2;
678 : 2118370 : f[0] = r[0]*qR2 + sR2*(uR[0] - r[0]);
679 : 2118370 : f[1] = r[1]*qR2 + pR*nx + sR2*(uR[1] - r[1]);
680 : 2118370 : f[2] = r[2]*qR2 + pR*ny + sR2*(uR[2] - r[2]);
681 : 2118370 : f[3] = r[3]*qR2 + pR*nz + sR2*(uR[3] - r[3]);
682 : 2118370 : f[4] = (r[4] + pR)*qR2 + sR2*(uR[4] - r[4]);
683 : 2118370 : }
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 [ + - ]: 4623900 : 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 : 4623900 : }
724 : :
725 : : static
726 : : std::function< void( const tk::UnsMesh::Coords&,
727 : : const tk::Fields&,
728 : : const tk::real[],
729 : : std::size_t,
730 : : std::size_t,
731 : : const tk::real[],
732 : : const tk::real[],
733 : : tk::real[] ) >
734 : 10200 : FLUX()
735 : : // *****************************************************************************
736 : : //! Configure flux function to use
737 : : //! \return Function to call to compute flux between two edge-end points
738 : : // *****************************************************************************
739 : : {
740 : : const auto& flux = g_cfg.get< tag::flux >();
741 : :
742 [ + + ]: 10200 : if (flux == "rusanov")
743 : : return rusanov;
744 [ + - ]: 240 : else if (flux == "hllc")
745 : : return hllc;
746 : : else
747 [ - - ][ - - ]: 0 : Throw( "Flux not configured" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
748 : : }
749 : :
750 : : static void
751 : 10200 : advdom( const tk::UnsMesh::Coords& coord,
752 : : const std::array< std::vector< std::size_t >, 3 >& dsupedge,
753 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
754 : : const tk::Fields& G,
755 : : const tk::Fields& U,
756 : : // cppcheck-suppress constParameter
757 : : tk::Fields& R )
758 : : // *****************************************************************************
759 : : //! Compute domain integral for advection
760 : : //! \param[in] coord Mesh node coordinates
761 : : //! \param[in] dsupedge Domain superedges
762 : : //! \param[in] dsupint Domain superedge integrals
763 : : //! \param[in] G Nodal gradients
764 : : //! \param[in] U Solution vector of primitive variables at recent time step
765 : : //! \param[in,out] R Right-hand side vector computed
766 : : // *****************************************************************************
767 : : {
768 : : // number of transported scalars
769 : : auto ncomp = U.nprop();
770 : :
771 : : // configure flux function
772 : 10200 : auto flux = FLUX();
773 : :
774 : : #if defined(__clang__)
775 : : #pragma clang diagnostic push
776 : : #pragma clang diagnostic ignored "-Wvla"
777 : : #pragma clang diagnostic ignored "-Wvla-extension"
778 : : #elif defined(STRICT_GNUC)
779 : : #pragma GCC diagnostic push
780 : : #pragma GCC diagnostic ignored "-Wvla"
781 : : #endif
782 : :
783 : : // domain edge contributions: tetrahedron superedges
784 [ + + ]: 1603650 : for (std::size_t e=0; e<dsupedge[0].size()/4; ++e) {
785 : 1593450 : const auto N = dsupedge[0].data() + e*4;
786 : 1593450 : tk::real u[4][ncomp];
787 [ + + ]: 9560700 : for (std::size_t c=0; c<ncomp; ++c) {
788 : 7967250 : u[0][c] = U(N[0],c);
789 : 7967250 : u[1][c] = U(N[1],c);
790 : 7967250 : u[2][c] = U(N[2],c);
791 : 7967250 : u[3][c] = U(N[3],c);
792 : : }
793 : : // edge fluxes
794 [ - + ]: 1593450 : tk::real f[6][ncomp];
795 : : const auto d = dsupint[0].data();
796 [ - + ]: 1593450 : flux( coord, G, d+(e*6+0)*3, N[0], N[1], u[0], u[1], f[0] );
797 [ - + ]: 1593450 : flux( coord, G, d+(e*6+1)*3, N[1], N[2], u[1], u[2], f[1] );
798 [ - + ]: 1593450 : flux( coord, G, d+(e*6+2)*3, N[2], N[0], u[2], u[0], f[2] );
799 [ - + ]: 1593450 : flux( coord, G, d+(e*6+3)*3, N[0], N[3], u[0], u[3], f[3] );
800 [ - + ]: 1593450 : flux( coord, G, d+(e*6+4)*3, N[1], N[3], u[1], u[3], f[4] );
801 [ - + ]: 1593450 : flux( coord, G, d+(e*6+5)*3, N[2], N[3], u[2], u[3], f[5] );
802 : : // edge flux contributions
803 [ + + ]: 9560700 : for (std::size_t c=0; c<ncomp; ++c) {
804 : 7967250 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c] - f[3][c];
805 : 7967250 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c] - f[4][c];
806 : 7967250 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c] - f[5][c];
807 : 7967250 : R(N[3],c) = R(N[3],c) + f[3][c] + f[4][c] + f[5][c];
808 : : }
809 : 1593450 : }
810 : :
811 : : // domain edge contributions: triangle superedges
812 [ + + ]: 693990 : for (std::size_t e=0; e<dsupedge[1].size()/3; ++e) {
813 : 683790 : const auto N = dsupedge[1].data() + e*3;
814 : 683790 : tk::real u[3][ncomp];
815 [ + + ]: 4102740 : for (std::size_t c=0; c<ncomp; ++c) {
816 : 3418950 : u[0][c] = U(N[0],c);
817 : 3418950 : u[1][c] = U(N[1],c);
818 : 3418950 : u[2][c] = U(N[2],c);
819 : : }
820 : : // edge fluxes
821 [ - + ]: 683790 : tk::real f[3][ncomp];
822 : : const auto d = dsupint[1].data();
823 [ - + ]: 683790 : flux( coord, G, d+(e*3+0)*3, N[0], N[1], u[0], u[1], f[0] );
824 [ - + ]: 683790 : flux( coord, G, d+(e*3+1)*3, N[1], N[2], u[1], u[2], f[1] );
825 [ - + ]: 683790 : flux( coord, G, d+(e*3+2)*3, N[2], N[0], u[2], u[0], f[2] );
826 : : // edge flux contributions
827 [ + + ]: 4102740 : for (std::size_t c=0; c<ncomp; ++c) {
828 : 3418950 : R(N[0],c) = R(N[0],c) - f[0][c] + f[2][c];
829 : 3418950 : R(N[1],c) = R(N[1],c) + f[0][c] - f[1][c];
830 : 3418950 : R(N[2],c) = R(N[2],c) + f[1][c] - f[2][c];
831 : : }
832 : 683790 : }
833 : :
834 : : // domain edge contributions: edges
835 [ + + ]: 3355350 : for (std::size_t e=0; e<dsupedge[2].size()/2; ++e) {
836 : 3345150 : const auto N = dsupedge[2].data() + e*2;
837 : 3345150 : tk::real u[2][ncomp];
838 [ + + ]: 20070900 : for (std::size_t c=0; c<ncomp; ++c) {
839 : 16725750 : u[0][c] = U(N[0],c);
840 : 16725750 : u[1][c] = U(N[1],c);
841 : : }
842 : : // edge fluxes
843 [ - + ]: 3345150 : tk::real f[ncomp];
844 : : const auto d = dsupint[2].data();
845 [ - + ]: 3345150 : flux( coord, G, d+e*3, N[0], N[1], u[0], u[1], f );
846 : : // edge flux contributions
847 [ + + ]: 20070900 : for (std::size_t c=0; c<ncomp; ++c) {
848 : 16725750 : R(N[0],c) -= f[c];
849 : 16725750 : R(N[1],c) += f[c];
850 : : }
851 : 3345150 : }
852 : :
853 : : #if defined(__clang__)
854 : : #pragma clang diagnostic pop
855 : : #elif defined(STRICT_GNUC)
856 : : #pragma GCC diagnostic pop
857 : : #endif
858 : 10200 : }
859 : :
860 : : static void
861 : 10200 : advbnd( const std::vector< std::size_t >& triinpoel,
862 : : const std::array< std::vector< tk::real >, 3 >& coord,
863 : : const std::vector< std::uint8_t >& besym,
864 : : const tk::Fields& U,
865 : : tk::Fields& R )
866 : : // *****************************************************************************
867 : : //! Compute boundary integral for advection
868 : : //! \param[in] triinpoel Boundary face connectivity
869 : : //! \param[in] coord Mesh node coordinates
870 : : //! \param[in] besym Boundary element symmetry BC flags
871 : : //! \param[in] U Solution vector at recent time step
872 : : //! \param[in,out] R Right-hand side vector
873 : : // *****************************************************************************
874 : : {
875 : : auto ncomp = U.nprop();
876 : :
877 : 10200 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
878 : 10200 : auto rgas = g_cfg.get< tag::mat_spec_gas_const >();
879 : :
880 : : const auto& x = coord[0];
881 : : const auto& y = coord[1];
882 : : const auto& z = coord[2];
883 : :
884 : : #if defined(__clang__)
885 : : #pragma clang diagnostic push
886 : : #pragma clang diagnostic ignored "-Wvla"
887 : : #pragma clang diagnostic ignored "-Wvla-extension"
888 : : #elif defined(STRICT_GNUC)
889 : : #pragma GCC diagnostic push
890 : : #pragma GCC diagnostic ignored "-Wvla"
891 : : #endif
892 : :
893 [ + + ]: 1385760 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
894 : 1375560 : const auto N = triinpoel.data() + e*3;
895 : :
896 : 1375560 : auto rA = U(N[0],0)/U(N[0],4)/rgas;
897 : 1375560 : auto ruA = U(N[0],1) * rA;
898 : 1375560 : auto rvA = U(N[0],2) * rA;
899 : 1375560 : auto rwA = U(N[0],3) * rA;
900 : 1375560 : auto reA = U(N[0],0)/(g-1.0) + 0.5*(ruA*ruA + rvA*rvA + rwA*rwA)/rA;
901 : :
902 : 1375560 : auto rB = U(N[1],0)/U(N[1],4)/rgas;
903 : 1375560 : auto ruB = U(N[1],1) * rB;
904 : 1375560 : auto rvB = U(N[1],2) * rB;
905 : 1375560 : auto rwB = U(N[1],3) * rB;
906 : 1375560 : auto reB = U(N[1],0)/(g-1.0) + 0.5*(ruB*ruB + rvB*rvB + rwB*rwB)/rB;
907 : :
908 : 1375560 : auto rC = U(N[2],0)/U(N[2],4)/rgas;
909 : 1375560 : auto ruC = U(N[2],1) * rC;
910 : 1375560 : auto rvC = U(N[2],2) * rC;
911 : 1375560 : auto rwC = U(N[2],3) * rC;
912 : 1375560 : auto reC = U(N[2],0)/(g-1.0) + 0.5*(ruC*ruC + rvC*rvC + rwC*rwC)/rC;
913 : :
914 : : const std::array< tk::real, 3 >
915 : 1375560 : ba{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] },
916 [ + + ]: 1375560 : ca{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] };
917 : : auto [nx,ny,nz] = tk::cross( ba, ca ); // 2A
918 : 1375560 : nx /= 12.0;
919 : 1375560 : ny /= 12.0;
920 : 1375560 : nz /= 12.0;
921 : :
922 [ + + ]: 1375560 : tk::real vn, f[ncomp][3];
923 : 1375560 : const auto sym = besym.data() + e*3;
924 : :
925 [ + + ]: 1375560 : vn = sym[0] ? 0.0 : (nx*U(N[0],1) + ny*U(N[0],2) + nz*U(N[0],3));
926 : : // flow
927 : 1375560 : f[0][0] = rA*vn;
928 : 1375560 : f[1][0] = ruA*vn + U(N[0],0)*nx;
929 : 1375560 : f[2][0] = rvA*vn + U(N[0],0)*ny;
930 : 1375560 : f[3][0] = rwA*vn + U(N[0],0)*nz;
931 : 1375560 : f[4][0] = (reA + U(N[0],0))*vn;
932 : : // scalar
933 [ - + ]: 1375560 : for (std::size_t c=5; c<ncomp; ++c) f[c][0] = U(N[0],c)*vn;
934 : :
935 [ + + ]: 1375560 : vn = sym[1] ? 0.0 : (nx*U(N[1],1) + ny*U(N[1],2) + nz*U(N[1],3));
936 : : // flow
937 : 1375560 : f[0][1] = rB*vn;
938 : 1375560 : f[1][1] = ruB*vn + U(N[1],0)*nx;
939 : 1375560 : f[2][1] = rvB*vn + U(N[1],0)*ny;
940 : 1375560 : f[3][1] = rwB*vn + U(N[1],0)*nz;
941 : 1375560 : f[4][1] = (reB + U(N[1],0))*vn;
942 : : // scalar
943 [ - + ]: 1375560 : for (std::size_t c=5; c<ncomp; ++c) f[c][1] = U(N[1],c)*vn;
944 : :
945 [ + + ]: 1375560 : vn = sym[2] ? 0.0 : (nx*U(N[2],1) + ny*U(N[2],2) + nz*U(N[2],3));
946 : : // flow
947 : 1375560 : f[0][2] = rC*vn;
948 : 1375560 : f[1][2] = ruC*vn + U(N[2],0)*nx;
949 : 1375560 : f[2][2] = rvC*vn + U(N[2],0)*ny;
950 : 1375560 : f[3][2] = rwC*vn + U(N[2],0)*nz;
951 : 1375560 : f[4][2] = (reC + U(N[2],0))*vn;
952 : : // scalar
953 [ - + ]: 1375560 : for (std::size_t c=5; c<ncomp; ++c) f[c][2] = U(N[2],c)*vn;
954 : :
955 [ + + ]: 8253360 : for (std::size_t c=0; c<ncomp; ++c) {
956 : 6877800 : auto fab = (f[c][0] + f[c][1])/4.0;
957 : 6877800 : auto fbc = (f[c][1] + f[c][2])/4.0;
958 : 6877800 : auto fca = (f[c][2] + f[c][0])/4.0;
959 : 6877800 : R(N[0],c) += fab + fca + f[c][0];
960 : 6877800 : R(N[1],c) += fab + fbc + f[c][1];
961 : 6877800 : R(N[2],c) += fbc + fca + f[c][2];
962 : : }
963 [ + + ]: 1375560 : }
964 : :
965 : : #if defined(__clang__)
966 : : #pragma clang diagnostic pop
967 : : #elif defined(STRICT_GNUC)
968 : : #pragma GCC diagnostic pop
969 : : #endif
970 : 10200 : }
971 : :
972 : : static void
973 : 10200 : src( const std::array< std::vector< tk::real >, 3 >& coord,
974 : : const std::vector< tk::real >& v,
975 : : tk::real t,
976 : : const std::vector< tk::real >& tp,
977 : : tk::Fields& R )
978 : : // *****************************************************************************
979 : : // Compute source integral
980 : : //! \param[in] coord Mesh node coordinates
981 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
982 : : //! \param[in] t Physical time
983 : : //! \param[in] tp Physical time for each mesh node
984 : : //! \param[in,out] R Right-hand side vector computed
985 : : // *****************************************************************************
986 : : {
987 : 10200 : auto src = problems::SRC();
988 [ + - ]: 10200 : if (!src) return;
989 : :
990 : : const auto& x = coord[0];
991 : : const auto& y = coord[1];
992 : : const auto& z = coord[2];
993 : :
994 [ - - ]: 0 : for (std::size_t p=0; p<R.nunk(); ++p) {
995 [ - - ]: 0 : if (g_cfg.get< tag::steady >()) t = tp[p];
996 [ - - ]: 0 : auto s = src( x[p], y[p], z[p], t );
997 [ - - ]: 0 : for (std::size_t c=0; c<s.size(); ++c) R(p,c) -= s[c] * v[p];
998 : : }
999 : : }
1000 : :
1001 : : void
1002 : 10200 : rhs( const std::array< std::vector< std::size_t >, 3 >& dsupedge,
1003 : : const std::array< std::vector< tk::real >, 3 >& dsupint,
1004 : : const std::array< std::vector< tk::real >, 3 >& coord,
1005 : : const std::vector< std::size_t >& triinpoel,
1006 : : const std::vector< std::uint8_t >& besym,
1007 : : const tk::Fields& G,
1008 : : const tk::Fields& U,
1009 : : const std::vector< tk::real >& v,
1010 : : tk::real t,
1011 : : const std::vector< tk::real >& tp,
1012 : : tk::Fields& R )
1013 : : // *****************************************************************************
1014 : : // Compute right hand side
1015 : : //! \param[in] dsupedge Domain superedges
1016 : : //! \param[in] dsupint Domain superedge integrals
1017 : : //! \param[in] coord Mesh node coordinates
1018 : : //! \param[in] triinpoel Boundary face connectivity
1019 : : //! \param[in] besym Boundary element symmetry BC flags
1020 : : //! \param[in] G Gradients in mesh nodes
1021 : : //! \param[in] U Solution vector of primitive variables at recent time step
1022 : : //! \param[in] v Nodal mesh volumes without contributions from other chares
1023 : : //! \param[in] t Physical time
1024 : : //! \param[in] tp Physical time for each mesh node
1025 : : //! \param[in,out] R Right-hand side vector computed
1026 : : // *****************************************************************************
1027 : : {
1028 : : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
1029 : : "vector at recent time step incorrect" );
1030 : : Assert( R.nunk() == coord[0].size(),
1031 : : "Number of unknowns and/or number of components in right-hand "
1032 : : "side vector incorrect" );
1033 : :
1034 : : R.fill( 0.0 );
1035 : 10200 : advdom( coord, dsupedge, dsupint, G, U, R );
1036 : 10200 : advbnd( triinpoel, coord, besym, U, R );
1037 : 10200 : src( coord, v, t, tp, R );
1038 : 10200 : }
1039 : :
1040 : : } // lax::
|
