Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/CompFlow/CGCompFlow.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Compressible single-material flow using continuous Galerkin
9 : : \details This file implements the physics operators governing compressible
10 : : single-material flow using continuous Galerkin discretization.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef CGCompFlow_h
14 : : #define CGCompFlow_h
15 : :
16 : : #include <cmath>
17 : : #include <algorithm>
18 : : #include <unordered_set>
19 : : #include <unordered_map>
20 : :
21 : : #include "DerivedData.hpp"
22 : : #include "Exception.hpp"
23 : : #include "Vector.hpp"
24 : : #include "EoS/EoS.hpp"
25 : : #include "Mesh/Around.hpp"
26 : : #include "Reconstruction.hpp"
27 : : #include "Problem/FieldOutput.hpp"
28 : : #include "Problem/BoxInitialization.hpp"
29 : : #include "Riemann/Rusanov.hpp"
30 : : #include "NodeBC.hpp"
31 : : #include "EoS/EoS.hpp"
32 : : #include "History.hpp"
33 : : #include "Table.hpp"
34 : :
35 : : namespace inciter {
36 : :
37 : : extern ctr::InputDeck g_inputdeck;
38 : :
39 : : namespace cg {
40 : :
41 : : //! \brief CompFlow used polymorphically with tk::CGPDE
42 : : //! \details The template arguments specify policies and are used to configure
43 : : //! the behavior of the class. The policies are:
44 : : //! - Physics - physics configuration, see PDE/CompFlow/Physics.h
45 : : //! - Problem - problem configuration, see PDE/CompFlow/Problems.h
46 : : //! \note The default physics is Euler, set in inciter::deck::check_compflow()
47 : : template< class Physics, class Problem >
48 : : class CompFlow {
49 : :
50 : : private:
51 : : using ncomp_t = kw::ncomp::info::expect::type;
52 : : using eq = tag::compflow;
53 : : using real = tk::real;
54 : : using param = tag::param;
55 : :
56 : : static constexpr std::size_t m_ncomp = 5;
57 : : static constexpr real muscl_eps = 1.0e-9;
58 : : static constexpr real muscl_const = 1.0/3.0;
59 : : static constexpr real muscl_m1 = 1.0 - muscl_const;
60 : : static constexpr real muscl_p1 = 1.0 + muscl_const;
61 : :
62 : : public:
63 : : //! \brief Constructor
64 : : //! \param[in] c Equation system index (among multiple systems configured)
65 : 161 : explicit CompFlow( ncomp_t c ) :
66 : : m_physics(),
67 : : m_problem(),
68 : : m_system( c ),
69 : 161 : m_offset( g_inputdeck.get< tag::component >().offset< eq >(c) ),
70 : : m_stagCnf( g_inputdeck.specialBC< eq, tag::stag >( c ) ),
71 : : m_skipCnf( g_inputdeck.specialBC< eq, tag::skip >( c ) ),
72 : 161 : m_fr( g_inputdeck.get< param, eq, tag::farfield_density >().size() > c ?
73 : 1 : g_inputdeck.get< param, eq, tag::farfield_density >()[c] : 1.0 ),
74 : 161 : m_fp( g_inputdeck.get< param, eq, tag::farfield_pressure >().size() > c ?
75 : 1 : g_inputdeck.get< param, eq, tag::farfield_pressure >()[c] : 1.0 ),
76 : 161 : m_fu( g_inputdeck.get< param, eq, tag::farfield_velocity >().size() > c ?
77 : 1 : g_inputdeck.get< param, eq, tag::farfield_velocity >()[c] :
78 [ + - ][ + + ]: 807 : std::vector< real >( 3, 0.0 ) )
[ + + ][ + + ]
[ + - ][ + - ]
[ + + ][ - - ]
79 : : {
80 [ + - ][ - + ]: 161 : Assert( g_inputdeck.get< tag::component >().get< eq >().at(c) == m_ncomp,
[ - - ][ - - ]
[ - - ][ - - ]
81 : : "Number of CompFlow PDE components must be " + std::to_string(m_ncomp) );
82 : 161 : }
83 : :
84 : : //! Determine nodes that lie inside the user-defined IC box
85 : : //! \param[in] coord Mesh node coordinates
86 : : //! \param[in,out] inbox List of nodes at which box user ICs are set for
87 : : //! each IC box
88 : 331 : void IcBoxNodes( const tk::UnsMesh::Coords& coord,
89 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
90 : : {
91 : 331 : const auto& x = coord[0];
92 : 331 : const auto& y = coord[1];
93 : 331 : const auto& z = coord[2];
94 : :
95 : : // Detect if user has configured a IC boxes
96 : 331 : const auto& icbox = g_inputdeck.get<tag::param, eq, tag::ic, tag::box>();
97 [ + - ]: 331 : if (icbox.size() > m_system) {
98 : 331 : std::size_t bcnt = 0;
99 [ + + ]: 337 : for (const auto& b : icbox[m_system]) { // for all boxes for this eq
100 [ + - ]: 6 : inbox.emplace_back();
101 [ + - ]: 6 : std::vector< tk::real > box
102 : 6 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
103 : 6 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
104 : 6 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
105 : :
106 : 6 : const auto eps = std::numeric_limits< tk::real >::epsilon();
107 : : // Determine which nodes lie in the IC box
108 [ + - ][ + - ]: 18 : if ( std::any_of( begin(box), end(box), [=](auto p)
109 : 6 : { return abs(p) > eps; } ) )
110 : : {
111 [ + + ]: 5137 : for (ncomp_t i=0; i<x.size(); ++i) {
112 [ + + ][ + + ]: 6971 : if ( x[i]>box[0] && x[i]<box[1] && y[i]>box[2] && y[i]<box[3] &&
[ + + ]
113 [ + + ][ + + ]: 6971 : z[i]>box[4] && z[i]<box[5] )
[ + + ][ + + ]
114 : : {
115 [ + - ]: 759 : inbox[bcnt].insert( i );
116 : : }
117 : : }
118 : : }
119 : 6 : ++bcnt;
120 : : }
121 : : }
122 : 331 : }
123 : :
124 : : //! Initalize the compressible flow equations, prepare for time integration
125 : : //! \param[in] coord Mesh node coordinates
126 : : //! \param[in,out] unk Array of unknowns
127 : : //! \param[in] t Physical time
128 : : //! \param[in] V Discrete volume of user-defined IC box
129 : : //! \param[in] inbox List of nodes at which box user ICs are set (for each
130 : : //! box IC)
131 : 1815 : void initialize(
132 : : const std::array< std::vector< real >, 3 >& coord,
133 : : tk::Fields& unk,
134 : : real t,
135 : : real V,
136 : : const std::vector< std::unordered_set< std::size_t > >& inbox ) const
137 : : {
138 [ - + ][ - - ]: 1815 : Assert( coord[0].size() == unk.nunk(), "Size mismatch" );
[ - - ][ - - ]
139 : :
140 : 1815 : const auto& x = coord[0];
141 : 1815 : const auto& y = coord[1];
142 : 1815 : const auto& z = coord[2];
143 : :
144 : 1815 : const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
145 : 1815 : const auto& icbox = ic.get< tag::box >();
146 : :
147 : 1815 : const auto eps = 1000.0 * std::numeric_limits< tk::real >::epsilon();
148 : :
149 : 1815 : const auto& bgpreic = ic.get< tag::pressure >();
150 : 1815 : tk::real bgpre =
151 [ + + ][ + - ]: 1815 : (bgpreic.size() > m_system && !bgpreic[m_system].empty()) ?
152 : 22 : bgpreic[m_system][0] : 0.0;
153 : :
154 : 1815 : auto c_v = cv< eq >(m_system);
155 : :
156 : : // Set initial and boundary conditions using problem policy
157 [ + + ]: 126859 : for (ncomp_t i=0; i<x.size(); ++i) {
158 [ + - ]: 125044 : auto s = Problem::initialize( m_system, m_ncomp, x[i], y[i], z[i], t );
159 : :
160 : : // initialize the user-defined box IC
161 [ + - ]: 125044 : if (icbox.size() > m_system) {
162 : 125044 : std::size_t bcnt = 0;
163 [ + + ]: 130175 : for (const auto& b : icbox[m_system]) { // for all boxes
164 [ + - ][ + - ]: 5131 : if (inbox.size() > bcnt && inbox[bcnt].find(i) != inbox[bcnt].end())
[ + + ][ + + ]
165 : : {
166 [ + - ]: 1518 : std::vector< tk::real > box
167 : 759 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
168 : 759 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
169 : 759 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
170 : 759 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
171 [ - + ]: 759 : if (V_ex < eps) V = 1.0;
172 [ + - ]: 759 : initializeBox( m_system, V_ex/V, t, b, bgpre, c_v, s );
173 : : }
174 : 5131 : ++bcnt;
175 : : }
176 : : }
177 : :
178 [ + - ]: 125044 : unk(i,0,m_offset) = s[0]; // rho
179 [ + - ][ + + ]: 125044 : if (!skipPoint(x[i],y[i],z[i]) && stagPoint(x[i],y[i],z[i])) {
[ + + ]
180 [ + - ][ + - ]: 18 : unk(i,1,m_offset) = unk(i,2,m_offset) = unk(i,3,m_offset) = 0.0;
[ + - ]
181 : : } else {
182 [ + - ]: 125026 : unk(i,1,m_offset) = s[1]; // rho * u
183 [ + - ]: 125026 : unk(i,2,m_offset) = s[2]; // rho * v
184 [ + - ]: 125026 : unk(i,3,m_offset) = s[3]; // rho * w
185 : : }
186 [ + - ]: 125044 : unk(i,4,m_offset) = s[4]; // rho * e, e: total = kinetic + internal
187 : : }
188 : 1815 : }
189 : :
190 : : //! Query the fluid velocity
191 : : //! \param[in] u Solution vector of conserved variables
192 : : //! \param[in,out] v Velocity components
193 : 25081 : void velocity( const tk::Fields& u, tk::UnsMesh::Coords& v ) const {
194 [ + + ]: 100324 : for (std::size_t j=0; j<3; ++j) {
195 : : // extract momentum
196 : 75243 : v[j] = u.extract( 1+j, m_offset );
197 [ - + ][ - - ]: 75243 : Assert( v[j].size() == u.nunk(), "Size mismatch" );
[ - - ][ - - ]
198 : : // divide by density
199 [ + + ]: 12254970 : for (std::size_t i=0; i<u.nunk(); ++i) v[j][i] /= u(i,0,m_offset);
200 : : }
201 : 25081 : }
202 : :
203 : : //! Query the sound speed
204 : : //! \param[in] U Solution vector of conserved variables
205 : : //! \param[in,out] s Speed of sound in mesh nodes
206 : 25081 : void soundspeed( const tk::Fields& U, std::vector< tk::real >& s ) const {
207 : 25081 : s.resize( U.nunk() );
208 [ + + ]: 4084990 : for (std::size_t i=0; i<U.nunk(); ++i) {
209 : 4059909 : auto r = U(i,0,m_offset);
210 : 4059909 : auto ru = U(i,1,m_offset);
211 : 4059909 : auto rv = U(i,2,m_offset);
212 : 4059909 : auto rw = U(i,3,m_offset);
213 : 4059909 : auto re = U(i,4,m_offset);
214 : 4059909 : auto p = eos_pressure< eq >( m_system, r, ru/r, rv/r, rw/r, re );
215 : 4059909 : s[i] = eos_soundspeed< eq >( m_system, r, p );
216 : : }
217 : 25081 : }
218 : :
219 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
220 : : //! \param[in] xi X-coordinate
221 : : //! \param[in] yi Y-coordinate
222 : : //! \param[in] zi Z-coordinate
223 : : //! \param[in] t Physical time
224 : : //! \return Vector of analytic solution at given location and time
225 : : std::vector< real >
226 : 167820 : analyticSolution( real xi, real yi, real zi, real t ) const
227 : 167820 : { return Problem::analyticSolution( m_system, m_ncomp, xi, yi, zi, t ); }
228 : :
229 : : //! Return analytic solution for conserved variables
230 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
231 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
232 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
233 : : //! \param[in] t Physical time at which to evaluate the analytic solution
234 : : //! \return Vector of analytic solution at given location and time
235 : : std::vector< tk::real >
236 : 1455903 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
237 : 1455903 : { return Problem::initialize( m_system, m_ncomp, xi, yi, zi, t ); }
238 : :
239 : : //! Compute right hand side for DiagCG (CG+FCT)
240 : : //! \param[in] t Physical time
241 : : //! \param[in] deltat Size of time step
242 : : //! \param[in] coord Mesh node coordinates
243 : : //! \param[in] inpoel Mesh element connectivity
244 : : //! \param[in] U Solution vector at recent time step
245 : : //! \param[in,out] Ue Element-centered solution vector at intermediate step
246 : : //! (used here internally as a scratch array)
247 : : //! \param[in,out] R Right-hand side vector computed
248 : 12004 : void rhs( real t,
249 : : real deltat,
250 : : const std::array< std::vector< real >, 3 >& coord,
251 : : const std::vector< std::size_t >& inpoel,
252 : : const tk::Fields& U,
253 : : tk::Fields& Ue,
254 : : tk::Fields& R ) const
255 : : {
256 [ - + ][ - - ]: 12004 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
257 : : "vector at recent time step incorrect" );
258 [ - + ][ - - ]: 12004 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
259 : : "Number of unknowns and/or number of components in right-hand "
260 : : "side vector incorrect" );
261 : :
262 : 12004 : const auto& x = coord[0];
263 : 12004 : const auto& y = coord[1];
264 : 12004 : const auto& z = coord[2];
265 : :
266 : : // 1st stage: update element values from node values (gather-add)
267 [ + + ]: 1139174 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
268 : : // access node IDs
269 : : const std::array< std::size_t, 4 >
270 : 1127170 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
271 : : // compute element Jacobi determinant
272 : : const std::array< real, 3 >
273 : 1127170 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
274 : 1127170 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
275 : 1127170 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
276 : 1127170 : const auto J = tk::triple( ba, ca, da ); // J = 6V
277 [ - + ][ - - ]: 1127170 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
278 : :
279 : : // shape function derivatives, nnode*ndim [4][3]
280 : : std::array< std::array< real, 3 >, 4 > grad;
281 : 1127170 : grad[1] = tk::crossdiv( ca, da, J );
282 : 1127170 : grad[2] = tk::crossdiv( da, ba, J );
283 : 1127170 : grad[3] = tk::crossdiv( ba, ca, J );
284 [ + + ]: 4508680 : for (std::size_t i=0; i<3; ++i)
285 : 3381510 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
286 : :
287 : : // access solution at element nodes
288 : : std::array< std::array< real, 4 >, m_ncomp > u;
289 [ + + ][ + - ]: 6763020 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, m_offset, N );
290 : :
291 : : // apply stagnation BCs
292 [ + + ]: 5635850 : for (std::size_t a=0; a<4; ++a)
293 [ + - ][ + + ]: 4508680 : if ( !skipPoint(x[N[a]],y[N[a]],z[N[a]]) &&
[ + + ]
294 : : stagPoint(x[N[a]],y[N[a]],z[N[a]]) )
295 : : {
296 : 428 : u[1][a] = u[2][a] = u[3][a] = 0.0;
297 : : }
298 : :
299 : : // access solution at elements
300 : : std::array< const real*, m_ncomp > ue;
301 [ + + ][ + - ]: 6763020 : for (ncomp_t c=0; c<m_ncomp; ++c) ue[c] = Ue.cptr( c, m_offset );
302 : :
303 : : // pressure
304 : : std::array< real, 4 > p;
305 [ + + ]: 5635850 : for (std::size_t a=0; a<4; ++a)
306 : 4508680 : p[a] = eos_pressure< eq >
307 [ + - ]: 9017360 : ( m_system, u[0][a], u[1][a]/u[0][a], u[2][a]/u[0][a],
308 : 4508680 : u[3][a]/u[0][a], u[4][a] );
309 : :
310 : : // sum flux contributions to element
311 : 1127170 : real d = deltat/2.0;
312 [ + + ]: 4508680 : for (std::size_t j=0; j<3; ++j)
313 [ + + ]: 16907550 : for (std::size_t a=0; a<4; ++a) {
314 : : // mass: advection
315 [ + - ]: 13526040 : Ue.var(ue[0],e) -= d * grad[a][j] * u[j+1][a];
316 : : // momentum: advection
317 [ + + ]: 54104160 : for (std::size_t i=0; i<3; ++i)
318 [ + - ]: 40578120 : Ue.var(ue[i+1],e) -= d * grad[a][j] * u[j+1][a]*u[i+1][a]/u[0][a];
319 : : // momentum: pressure
320 [ + - ]: 13526040 : Ue.var(ue[j+1],e) -= d * grad[a][j] * p[a];
321 : : // energy: advection and pressure
322 [ + - ]: 27052080 : Ue.var(ue[4],e) -= d * grad[a][j] *
323 : 13526040 : (u[4][a] + p[a]) * u[j+1][a]/u[0][a];
324 : : }
325 : :
326 : : // add (optional) source to all equations
327 [ + + ]: 5635850 : for (std::size_t a=0; a<4; ++a) {
328 : : real s[m_ncomp];
329 [ + - ]: 4508680 : Problem::src( m_system, x[N[a]], y[N[a]], z[N[a]], t,
330 : : s[0], s[1], s[2], s[3], s[4] );
331 [ + + ]: 27052080 : for (std::size_t c=0; c<m_ncomp; ++c)
332 [ + - ]: 22543400 : Ue.var(ue[c],e) += d/4.0 * s[c];
333 : : }
334 : : }
335 : :
336 : : // 2nd stage: form rhs from element values (scatter-add)
337 [ + + ]: 1139174 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
338 : : // access node IDs
339 : : const std::array< std::size_t, 4 >
340 : 1127170 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
341 : : // compute element Jacobi determinant
342 : : const std::array< real, 3 >
343 : 1127170 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
344 : 1127170 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
345 : 1127170 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
346 : 1127170 : const auto J = tk::triple( ba, ca, da ); // J = 6V
347 [ - + ][ - - ]: 1127170 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
348 : :
349 : : // shape function derivatives, nnode*ndim [4][3]
350 : : std::array< std::array< real, 3 >, 4 > grad;
351 : 1127170 : grad[1] = tk::crossdiv( ca, da, J );
352 : 1127170 : grad[2] = tk::crossdiv( da, ba, J );
353 : 1127170 : grad[3] = tk::crossdiv( ba, ca, J );
354 [ + + ]: 4508680 : for (std::size_t i=0; i<3; ++i)
355 : 3381510 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
356 : :
357 : : // access solution at elements
358 : : std::array< real, m_ncomp > ue;
359 [ + + ][ + - ]: 6763020 : for (ncomp_t c=0; c<m_ncomp; ++c) ue[c] = Ue( e, c, m_offset );
360 : : // access pointer to right hand side at component and offset
361 : : std::array< const real*, m_ncomp > r;
362 [ + + ][ + - ]: 6763020 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c, m_offset );
363 : :
364 : : // pressure
365 : : auto p = eos_pressure< eq >
366 [ + - ]: 1127170 : ( m_system, ue[0], ue[1]/ue[0], ue[2]/ue[0], ue[3]/ue[0],
367 : : ue[4] );
368 : :
369 : : // scatter-add flux contributions to rhs at nodes
370 : 1127170 : real d = deltat * J/6.0;
371 [ + + ]: 4508680 : for (std::size_t j=0; j<3; ++j)
372 [ + + ]: 16907550 : for (std::size_t a=0; a<4; ++a) {
373 : : // mass: advection
374 [ + - ]: 13526040 : R.var(r[0],N[a]) += d * grad[a][j] * ue[j+1];
375 : : // momentum: advection
376 [ + + ]: 54104160 : for (std::size_t i=0; i<3; ++i)
377 [ + - ]: 40578120 : R.var(r[i+1],N[a]) += d * grad[a][j] * ue[j+1]*ue[i+1]/ue[0];
378 : : // momentum: pressure
379 [ + - ]: 13526040 : R.var(r[j+1],N[a]) += d * grad[a][j] * p;
380 : : // energy: advection and pressure
381 [ + - ]: 13526040 : R.var(r[4],N[a]) += d * grad[a][j] * (ue[4] + p) * ue[j+1]/ue[0];
382 : : }
383 : :
384 : : // add (optional) source to all equations
385 : 1127170 : auto xc = (x[N[0]] + x[N[1]] + x[N[2]] + x[N[3]]) / 4.0;
386 : 1127170 : auto yc = (y[N[0]] + y[N[1]] + y[N[2]] + y[N[3]]) / 4.0;
387 : 1127170 : auto zc = (z[N[0]] + z[N[1]] + z[N[2]] + z[N[3]]) / 4.0;
388 : : real s[m_ncomp];
389 [ + - ]: 1127170 : Problem::src( m_system, xc, yc, zc, t+deltat/2,
390 : : s[0], s[1], s[2], s[3], s[4] );
391 [ + + ]: 6763020 : for (std::size_t c=0; c<m_ncomp; ++c)
392 [ + + ]: 28179250 : for (std::size_t a=0; a<4; ++a)
393 [ + - ]: 22543400 : R.var(r[c],N[a]) += d/4.0 * s[c];
394 : : }
395 : : // // add viscous stress contribution to momentum and energy rhs
396 : : // m_physics.viscousRhs( deltat, J, N, grad, u, r, R );
397 : : // // add heat conduction contribution to energy rhs
398 : : // m_physics.conductRhs( deltat, J, N, grad, u, r, R );
399 : 12004 : }
400 : :
401 : : //! \brief Compute nodal gradients of primitive variables for ALECG along
402 : : //! chare-boundary
403 : : //! \param[in] coord Mesh node coordinates
404 : : //! \param[in] inpoel Mesh element connectivity
405 : : //! \param[in] bndel List of elements contributing to chare-boundary nodes
406 : : //! \param[in] gid Local->global node id map
407 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
408 : : //! global node ids (key)
409 : : //! \param[in] U Solution vector at recent time step
410 : : //! \param[in,out] G Nodal gradients of primitive variables
411 : : //! \details This function only computes local contributions to gradients
412 : : //! at chare-boundary nodes. Internal node gradients are calculated as
413 : : //! required, and do not need to be stored.
414 : 24936 : void chBndGrad( const std::array< std::vector< real >, 3 >& coord,
415 : : const std::vector< std::size_t >& inpoel,
416 : : const std::vector< std::size_t >& bndel,
417 : : const std::vector< std::size_t >& gid,
418 : : const std::unordered_map< std::size_t, std::size_t >& bid,
419 : : const tk::Fields& U,
420 : : tk::Fields& G ) const
421 : : {
422 [ - + ][ - - ]: 24936 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
423 : : "vector at recent time step incorrect" );
424 : :
425 : : // compute gradients of primitive variables in points
426 : 24936 : G.fill( 0.0 );
427 : :
428 : : // access node cooordinates
429 : 24936 : const auto& x = coord[0];
430 : 24936 : const auto& y = coord[1];
431 : 24936 : const auto& z = coord[2];
432 : :
433 [ + + ]: 4893687 : for (auto e : bndel) { // elements contributing to chare boundary nodes
434 : : // access node IDs
435 : 4868751 : std::size_t N[4] =
436 : 4868751 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
437 : : // compute element Jacobi determinant, J = 6V
438 : 4868751 : real bax = x[N[1]]-x[N[0]];
439 : 4868751 : real bay = y[N[1]]-y[N[0]];
440 : 4868751 : real baz = z[N[1]]-z[N[0]];
441 : 4868751 : real cax = x[N[2]]-x[N[0]];
442 : 4868751 : real cay = y[N[2]]-y[N[0]];
443 : 4868751 : real caz = z[N[2]]-z[N[0]];
444 : 4868751 : real dax = x[N[3]]-x[N[0]];
445 : 4868751 : real day = y[N[3]]-y[N[0]];
446 : 4868751 : real daz = z[N[3]]-z[N[0]];
447 : 4868751 : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
448 [ - + ][ - - ]: 4868751 : ErrChk( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
449 : 4868751 : auto J24 = J/24.0;
450 : : // shape function derivatives, nnode*ndim [4][3]
451 : : real g[4][3];
452 : 4868751 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
453 : : g[1][0], g[1][1], g[1][2] );
454 : 4868751 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
455 : : g[2][0], g[2][1], g[2][2] );
456 : 4868751 : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
457 : : g[3][0], g[3][1], g[3][2] );
458 [ + + ]: 19475004 : for (std::size_t i=0; i<3; ++i)
459 : 14606253 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
460 : : // scatter-add gradient contributions to boundary nodes
461 [ + + ]: 24343755 : for (std::size_t a=0; a<4; ++a) {
462 [ + - ]: 19475004 : auto i = bid.find( gid[N[a]] );
463 [ + + ]: 19475004 : if (i != end(bid)) {
464 : : real u[5];
465 [ + + ]: 60806790 : for (std::size_t b=0; b<4; ++b) {
466 [ + - ]: 48645432 : u[0] = U(N[b],0,m_offset);
467 [ + - ]: 48645432 : u[1] = U(N[b],1,m_offset)/u[0];
468 [ + - ]: 48645432 : u[2] = U(N[b],2,m_offset)/u[0];
469 [ + - ]: 48645432 : u[3] = U(N[b],3,m_offset)/u[0];
470 [ + - ]: 48645432 : u[4] = U(N[b],4,m_offset)/u[0]
471 : 48645432 : - 0.5*(u[1]*u[1] + u[2]*u[2] + u[3]*u[3]);
472 [ + - ][ - + ]: 48645432 : if ( !skipPoint(x[N[b]],y[N[b]],z[N[b]]) &&
[ - + ]
473 : : stagPoint(x[N[b]],y[N[b]],z[N[b]]) )
474 : : {
475 : 0 : u[1] = u[2] = u[3] = 0.0;
476 : : }
477 [ + + ]: 291872592 : for (std::size_t c=0; c<5; ++c)
478 [ + + ]: 972908640 : for (std::size_t j=0; j<3; ++j)
479 [ + - ]: 729681480 : G(i->second,c*3+j,0) += J24 * g[b][j] * u[c];
480 : : }
481 : : }
482 : : }
483 : : }
484 : 24936 : }
485 : :
486 : : //! Compute right hand side for ALECG
487 : : //! \param[in] t Physical time
488 : : //! \param[in] coord Mesh node coordinates
489 : : //! \param[in] inpoel Mesh element connectivity
490 : : //! \param[in] triinpoel Boundary triangle face connecitivity with local ids
491 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
492 : : //! global node ids (key)
493 : : //! \param[in] gid Local->glocal node ids
494 : : //! \param[in] lid Global->local node ids
495 : : //! \param[in] dfn Dual-face normals
496 : : //! \param[in] psup Points surrounding points
497 : : //! \param[in] esup Elements surrounding points
498 : : //! \param[in] symbctri Vector with 1 at symmetry BC boundary triangles
499 : : //! \param[in] spongenodes Unique set of nodes at which to apply sponge
500 : : // conditions
501 : : //! \param[in] vol Nodal volumes
502 : : //! \param[in] edgenode Local node IDs of edges
503 : : //! \param[in] edgeid Edge ids in the order of access
504 : : //! \param[in] boxnodes Mesh node ids within user-defined IC boxes
505 : : //! \param[in] G Nodal gradients
506 : : //! \param[in] U Solution vector at recent time step
507 : : //! \param[in] W Mesh velocity
508 : : //! \param[in] tp Physical time for each mesh node
509 : : //! \param[in] V Total box volume
510 : : //! \param[in,out] R Right-hand side vector computed
511 : 24936 : void rhs( real t,
512 : : const std::array< std::vector< real >, 3 >& coord,
513 : : const std::vector< std::size_t >& inpoel,
514 : : const std::vector< std::size_t >& triinpoel,
515 : : const std::vector< std::size_t >& gid,
516 : : const std::unordered_map< std::size_t, std::size_t >& bid,
517 : : const std::unordered_map< std::size_t, std::size_t >& lid,
518 : : const std::vector< real >& dfn,
519 : : const std::pair< std::vector< std::size_t >,
520 : : std::vector< std::size_t > >& psup,
521 : : const std::pair< std::vector< std::size_t >,
522 : : std::vector< std::size_t > >& esup,
523 : : const std::vector< int >& symbctri,
524 : : const std::unordered_set< std::size_t >& spongenodes,
525 : : const std::vector< real >& vol,
526 : : const std::vector< std::size_t >& edgenode,
527 : : const std::vector< std::size_t >& edgeid,
528 : : const std::vector< std::unordered_set< std::size_t > >& boxnodes,
529 : : const tk::Fields& G,
530 : : const tk::Fields& U,
531 : : const tk::Fields& W,
532 : : const std::vector< tk::real >& tp,
533 : : real V,
534 : : tk::Fields& R ) const
535 : : {
536 [ - + ][ - - ]: 24936 : Assert( G.nprop() == m_ncomp*3,
[ - - ][ - - ]
537 : : "Number of components in gradient vector incorrect" );
538 [ - + ][ - - ]: 24936 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
539 : : "vector at recent time step incorrect" );
540 [ - + ][ - - ]: 24936 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
541 : : "Number of unknowns and/or number of components in right-hand "
542 : : "side vector incorrect" );
543 [ - + ][ - - ]: 24936 : Assert( W.nunk() == coord[0].size(), "Size mismatch " );
[ - - ][ - - ]
544 : :
545 : : // compute/assemble gradients in points
546 [ + - ]: 49872 : auto Grad = nodegrad( coord, inpoel, lid, bid, vol, esup, U, G );
547 : :
548 : : // zero right hand side for all components
549 [ + + ][ + - ]: 149616 : for (ncomp_t c=0; c<m_ncomp; ++c) R.fill( c, m_offset, 0.0 );
550 : :
551 : : // compute sponge pressure multiplers at sponge side sets
552 [ + - ]: 49872 : auto spmult = spongePressures( coord, spongenodes );
553 : :
554 : : // compute domain-edge integral
555 [ + - ]: 24936 : domainint( coord, gid, edgenode, edgeid, psup, dfn, U, W, Grad,
556 : : spmult, R );
557 : :
558 : : // compute boundary integrals
559 [ + - ]: 24936 : bndint( coord, triinpoel, symbctri, U, W, spmult, R );
560 : :
561 : : // compute external (energy) sources
562 : 24936 : const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
563 : 24936 : const auto& icbox = ic.get< tag::box >();
564 : :
565 [ + - ][ + + ]: 24936 : if (icbox.size() > m_system && !boxnodes.empty()) {
[ + + ]
566 : 1230 : std::size_t bcnt = 0;
567 [ + + ]: 2490 : for (const auto& b : icbox[m_system]) { // for all boxes for this eq
568 [ + - ]: 1260 : std::vector< tk::real > box
569 : 1260 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
570 : 1260 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
571 : 1260 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
572 : :
573 : 1260 : const auto& initiate = b.template get< tag::initiate >();
574 : 1260 : auto inittype = initiate.template get< tag::init >();
575 [ + + ]: 1260 : if (inittype == ctr::InitiateType::LINEAR) {
576 [ + - ]: 1200 : boxSrc( V, t, inpoel, esup, boxnodes[bcnt], coord, R );
577 : : }
578 : 1260 : ++bcnt;
579 : : }
580 : : }
581 : :
582 : : // compute optional source integral
583 [ + - ]: 24936 : src( coord, inpoel, t, tp, R );
584 : 24936 : }
585 : :
586 : : //! Compute the minimum time step size (for unsteady time stepping)
587 : : //! \param[in] coord Mesh node coordinates
588 : : //! \param[in] inpoel Mesh element connectivity
589 : : //! \param[in] t Physical time
590 : : //! \param[in] dtn Time step size at the previous time step
591 : : //! \param[in] U Solution vector at recent time step
592 : : //! \param[in] vol Nodal volume (with contributions from other chares)
593 : : //! \param[in] voln Nodal volume (with contributions from other chares) at
594 : : //! the previous time step
595 : : //! \return Minimum time step size
596 : 20116 : real dt( const std::array< std::vector< real >, 3 >& coord,
597 : : const std::vector< std::size_t >& inpoel,
598 : : tk::real t,
599 : : tk::real dtn,
600 : : const tk::Fields& U,
601 : : const std::vector< tk::real >& vol,
602 : : const std::vector< tk::real >& voln ) const
603 : : {
604 [ - + ][ - - ]: 20116 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
605 : : "vector at recent time step incorrect" );
606 : :
607 : : // energy source propagation time and velocity
608 : 20116 : const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
609 : 20116 : const auto& icbox = ic.get< tag::box >();
610 : :
611 : 20116 : const auto& x = coord[0];
612 : 20116 : const auto& y = coord[1];
613 : 20116 : const auto& z = coord[2];
614 : :
615 : : // ratio of specific heats
616 [ + - ]: 20116 : auto g = gamma< eq >(m_system);
617 : : // compute the minimum dt across all elements we own
618 : 20116 : real mindt = std::numeric_limits< real >::max();
619 [ + + ]: 5135858 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
620 : 5115742 : const std::array< std::size_t, 4 > N{{ inpoel[e*4+0], inpoel[e*4+1],
621 : 5115742 : inpoel[e*4+2], inpoel[e*4+3] }};
622 : : // compute cubic root of element volume as the characteristic length
623 : : const std::array< real, 3 >
624 : 5115742 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
625 : 5115742 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
626 : 5115742 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
627 : 5115742 : const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
628 : : // access solution at element nodes at recent time step
629 : : std::array< std::array< real, 4 >, m_ncomp > u;
630 [ + + ][ + - ]: 30694452 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, m_offset, N );
631 : : // compute the maximum length of the characteristic velocity (fluid
632 : : // velocity + sound velocity) across the four element nodes
633 : 5115742 : real maxvel = 0.0;
634 [ + + ]: 25578710 : for (std::size_t j=0; j<4; ++j) {
635 : 20462968 : auto& r = u[0][j]; // rho
636 : 20462968 : auto& ru = u[1][j]; // rho * u
637 : 20462968 : auto& rv = u[2][j]; // rho * v
638 : 20462968 : auto& rw = u[3][j]; // rho * w
639 : 20462968 : auto& re = u[4][j]; // rho * e
640 [ + - ]: 20462968 : auto p = eos_pressure< eq >( m_system, r, ru/r, rv/r, rw/r, re );
641 [ - + ]: 20462968 : if (p < 0) p = 0.0;
642 [ + - ]: 20462968 : auto c = eos_soundspeed< eq >( m_system, r, p );
643 : 20462968 : auto v = std::sqrt((ru*ru + rv*rv + rw*rw)/r/r) + c; // char. velocity
644 : :
645 : : // energy source propagation velocity (in all IC boxes configured)
646 [ + - ]: 20462968 : if (icbox.size() > m_system) {
647 [ + + ]: 25487048 : for (const auto& b : icbox[m_system]) { // for all boxes for this eq
648 : 5024080 : const auto& initiate = b.template get< tag::initiate >();
649 : 5024080 : auto iv = initiate.template get< tag::velocity >();
650 : 5024080 : auto inittype = initiate.template get< tag::init >();
651 [ + + ]: 5024080 : if (inittype == ctr::InitiateType::LINEAR) {
652 : 4902800 : auto zmin = b.template get< tag::zmin >();
653 : 4902800 : auto zmax = b.template get< tag::zmax >();
654 : 4902800 : auto wFront = 0.08;
655 : 4902800 : auto tInit = 0.0;
656 : 4902800 : auto tFinal = tInit + (zmax - zmin - 2.0*wFront) /
657 : 4902800 : std::fabs(iv);
658 [ + - ][ + + ]: 4902800 : if (t >= tInit && t <= tFinal)
659 : 4412520 : v = std::max(v, std::fabs(iv));
660 : : }
661 : : }
662 : : }
663 : :
664 [ + + ]: 20462968 : if (v > maxvel) maxvel = v;
665 : : }
666 : : // compute element dt for the Euler equations
667 : 5115742 : auto euler_dt = L / maxvel;
668 : : // compute element dt based on the viscous force
669 [ - - ]: 5115742 : auto viscous_dt = m_physics.viscous_dt( L, u );
670 : : // compute element dt based on thermal diffusion
671 [ - - ]: 5115742 : auto conduct_dt = m_physics.conduct_dt( L, g, u );
672 : : // compute minimum element dt
673 : 5115742 : auto elemdt = std::min( euler_dt, std::min( viscous_dt, conduct_dt ) );
674 : : // find minimum dt across all elements
675 : 5115742 : mindt = std::min( elemdt, mindt );
676 : : }
677 : 20116 : mindt *= g_inputdeck.get< tag::discr, tag::cfl >();
678 : :
679 : : // compute the minimum dt across all nodes we contribute to due to volume
680 : : // change in time
681 : 20116 : auto dvcfl = g_inputdeck.get< tag::ale, tag::dvcfl >();
682 [ + + ][ + + ]: 20116 : if (dtn > 0.0 && dvcfl > 0.0) {
683 [ - + ][ - - ]: 298 : Assert( vol.size() == voln.size(), "Size mismatch" );
[ - - ][ - - ]
684 [ + + ]: 375342 : for (std::size_t p=0; p<vol.size(); ++p) {
685 : 750088 : auto vol_dt = dtn *
686 : 375044 : std::min(voln[p],vol[p]) / std::abs(voln[p]-vol[p]+1.0e-14);
687 : 375044 : mindt = std::min( vol_dt, mindt );
688 : : }
689 : 298 : mindt *= dvcfl;
690 : : }
691 : :
692 : 20116 : return mindt;
693 : : }
694 : :
695 : : //! Compute a time step size for each mesh node (for steady time stepping)
696 : : //! \param[in] U Solution vector at recent time step
697 : : //! \param[in] vol Nodal volume (with contributions from other chares)
698 : : //! \param[in,out] dtp Time step size for each mesh node
699 : 200 : void dt( uint64_t,
700 : : const std::vector< tk::real >& vol,
701 : : const tk::Fields& U,
702 : : std::vector< tk::real >& dtp ) const
703 : : {
704 [ + + ]: 22520 : for (std::size_t i=0; i<U.nunk(); ++i) {
705 : : // compute cubic root of element volume as the characteristic length
706 : 22320 : const auto L = std::cbrt( vol[i] );
707 : : // access solution at node p at recent time step
708 [ + - ]: 22320 : const auto u = U[i];
709 : : // compute pressure
710 : : auto p = eos_pressure< eq >
711 [ + - ]: 22320 : ( m_system, u[0], u[1]/u[0], u[2]/u[0], u[3]/u[0], u[4] );
712 [ - + ]: 22320 : if (p < 0) p = 0.0;
713 [ + - ]: 22320 : auto c = eos_soundspeed< eq >( m_system, u[0], p );
714 : : // characteristic velocity
715 : 22320 : auto v = std::sqrt((u[1]*u[1] + u[2]*u[2] + u[3]*u[3])/u[0]/u[0]) + c;
716 : : // compute dt for node
717 : 22320 : dtp[i] = L / v * g_inputdeck.get< tag::discr, tag::cfl >();
718 : : }
719 : 200 : }
720 : :
721 : : //! \brief Query Dirichlet boundary condition value on a given side set for
722 : : //! all components in this PDE system
723 : : //! \param[in] t Physical time
724 : : //! \param[in] deltat Time step size
725 : : //! \param[in] tp Physical time for each mesh node
726 : : //! \param[in] dtp Time step size for each mesh node
727 : : //! \param[in] ss Pair of side set ID and (local) node IDs on the side set
728 : : //! \param[in] coord Mesh node coordinates
729 : : //! \param[in] increment If true, evaluate the solution increment between
730 : : //! t and t+dt for Dirichlet BCs. If false, evlauate the solution instead.
731 : : //! \return Vector of pairs of bool and boundary condition value associated
732 : : //! to mesh node IDs at which Dirichlet boundary conditions are set. Note
733 : : //! that if increment is true, instead of the actual boundary condition
734 : : //! value, we return the increment between t+deltat and t, since,
735 : : //! depending on client code and solver, that may be what the solution
736 : : //! requires.
737 : : std::map< std::size_t, std::vector< std::pair<bool,real> > >
738 : 98322 : dirbc( real t,
739 : : real deltat,
740 : : const std::vector< tk::real >& tp,
741 : : const std::vector< tk::real >& dtp,
742 : : const std::pair< const int, std::vector< std::size_t > >& ss,
743 : : const std::array< std::vector< real >, 3 >& coord,
744 : : bool increment ) const
745 : : {
746 : : using tag::param; using tag::bcdir;
747 : : using NodeBC = std::vector< std::pair< bool, real > >;
748 : 98322 : std::map< std::size_t, NodeBC > bc;
749 : 98322 : const auto& ubc = g_inputdeck.get< param, eq, tag::bc, bcdir >();
750 : 98322 : const auto steady = g_inputdeck.get< tag::discr, tag::steady_state >();
751 [ + + ]: 98322 : if (!ubc.empty()) {
752 [ - + ][ - - ]: 80008 : Assert( ubc.size() > 0, "Indexing out of Dirichlet BC eq-vector" );
[ - - ][ - - ]
753 : 80008 : const auto& x = coord[0];
754 : 80008 : const auto& y = coord[1];
755 : 80008 : const auto& z = coord[2];
756 [ + + ]: 557096 : for (const auto& b : ubc[0])
757 [ + - ][ + + ]: 477088 : if (std::stoi(b) == ss.first)
758 [ + + ]: 2272872 : for (auto n : ss.second) {
759 [ - + ][ - - ]: 2193444 : Assert( x.size() > n, "Indexing out of coordinate array" );
[ - - ][ - - ]
760 [ + + ]: 2193444 : if (steady) { t = tp[n]; deltat = dtp[n]; }
761 [ + + ][ + - ]: 4386888 : auto s = increment ?
[ + + ][ - - ]
762 [ + - ]: 891901 : solinc( m_system, m_ncomp, x[n], y[n], z[n],
763 : : t, deltat, Problem::initialize ) :
764 [ + - ]: 1301543 : Problem::initialize( m_system, m_ncomp, x[n], y[n], z[n],
765 : : t+deltat );
766 [ + - ][ + + ]: 2193444 : if ( !skipPoint(x[n],y[n],z[n]) && stagPoint(x[n],y[n],z[n]) ) {
[ + + ]
767 : 1125 : s[1] = s[2] = s[3] = 0.0;
768 : : }
769 [ + - ][ + - ]: 4386888 : bc[n] = {{ {true,s[0]}, {true,s[1]}, {true,s[2]}, {true,s[3]},
770 : 2193444 : {true,s[4]} }};
771 : : }
772 : : }
773 : 98322 : return bc;
774 : : }
775 : :
776 : : //! Set symmetry boundary conditions at nodes
777 : : //! \param[in] U Solution vector at recent time step
778 : : //! \param[in] coord Mesh node coordinates
779 : : //! \param[in] bnorm Face normals in boundary points, key local node id,
780 : : //! first 3 reals of value: unit normal, outer key: side set id
781 : : //! \param[in] nodes Unique set of node ids at which to set symmetry BCs
782 : : void
783 : 80589 : symbc( tk::Fields& U,
784 : : const std::array< std::vector< real >, 3 >& coord,
785 : : const std::unordered_map< int,
786 : : std::unordered_map< std::size_t, std::array< real, 4 > > >& bnorm,
787 : : const std::unordered_set< std::size_t >& nodes ) const
788 : : {
789 : 80589 : const auto& x = coord[0];
790 : 80589 : const auto& y = coord[1];
791 : 80589 : const auto& z = coord[2];
792 : 80589 : const auto& sbc = g_inputdeck.get< param, eq, tag::bc, tag::bcsym >();
793 [ + + ]: 80589 : if (sbc.size() > m_system) { // use symbcs for this system
794 [ + + ]: 2238875 : for (auto p : nodes) { // for all symbc nodes
795 [ + - ]: 2229934 : if (!skipPoint(x[p],y[p],z[p])) {
796 : : // for all user-def symbc sets
797 [ + + ]: 9289436 : for (std::size_t s=0; s<sbc[m_system].size(); ++s) {
798 : : // find nodes & normals for side
799 [ + - ][ + - ]: 7059502 : auto j = bnorm.find(std::stoi(sbc[m_system][s]));
800 [ + + ]: 7059502 : if (j != end(bnorm)) {
801 [ + - ]: 5594630 : auto i = j->second.find(p); // find normal for node
802 [ + + ]: 5594630 : if (i != end(j->second)) {
803 : : std::array< real, 3 >
804 : 2525642 : n{ i->second[0], i->second[1], i->second[2] },
805 [ + - ][ + - ]: 2525642 : v{ U(p,1,m_offset), U(p,2,m_offset), U(p,3,m_offset) };
[ + - ]
806 : 2525642 : auto v_dot_n = tk::dot( v, n );
807 : : // symbc: remove normal component of velocity
808 [ + - ]: 2525642 : U(p,1,m_offset) -= v_dot_n * n[0];
809 [ + - ]: 2525642 : U(p,2,m_offset) -= v_dot_n * n[1];
810 [ + - ]: 2525642 : U(p,3,m_offset) -= v_dot_n * n[2];
811 : : }
812 : : }
813 : : }
814 : : }
815 : : }
816 : : }
817 : 80589 : }
818 : :
819 : : //! Set farfield boundary conditions at nodes
820 : : //! \param[in] U Solution vector at recent time step
821 : : //! \param[in] coord Mesh node coordinates
822 : : //! \param[in] bnorm Face normals in boundary points, key local node id,
823 : : //! first 3 reals of value: unit normal, outer key: side set id
824 : : //! \param[in] nodes Unique set of node ids at which to set farfield BCs
825 : : void
826 : 68585 : farfieldbc(
827 : : tk::Fields& U,
828 : : const std::array< std::vector< real >, 3 >& coord,
829 : : const std::unordered_map< int,
830 : : std::unordered_map< std::size_t, std::array< real, 4 > > >& bnorm,
831 : : const std::unordered_set< std::size_t >& nodes ) const
832 : : {
833 : 68585 : const auto& x = coord[0];
834 : 68585 : const auto& y = coord[1];
835 : 68585 : const auto& z = coord[2];
836 : 68585 : const auto& fbc = g_inputdeck.get<param, eq, tag::bc, tag::bcfarfield>();
837 [ + + ]: 68585 : if (fbc.size() > m_system) // use farbcs for this system
838 [ + + ]: 84400 : for (auto p : nodes) // for all farfieldbc nodes
839 [ + - ]: 84189 : if (!skipPoint(x[p],y[p],z[p]))
840 [ + + ]: 505134 : for (const auto& s : fbc[m_system]) {// for all user-def farbc sets
841 [ + - ][ + - ]: 420945 : auto j = bnorm.find(std::stoi(s));// find nodes & normals for side
842 [ + + ]: 420945 : if (j != end(bnorm)) {
843 [ + - ]: 418950 : auto i = j->second.find(p); // find normal for node
844 [ + + ]: 418950 : if (i != end(j->second)) {
845 [ + - ]: 114450 : auto& r = U(p,0,m_offset);
846 [ + - ]: 114450 : auto& ru = U(p,1,m_offset);
847 [ + - ]: 114450 : auto& rv = U(p,2,m_offset);
848 [ + - ]: 114450 : auto& rw = U(p,3,m_offset);
849 [ + - ]: 114450 : auto& re = U(p,4,m_offset);
850 : 114450 : auto vn =
851 : 114450 : (ru*i->second[0] + rv*i->second[1] + rw*i->second[2]) / r;
852 [ + - ]: 114450 : auto a = eos_soundspeed< eq >( m_system, r,
853 [ + - ]: 114450 : eos_pressure< eq >( m_system, r, ru/r, rv/r, rw/r, re ) );
854 : 114450 : auto M = vn / a;
855 [ - + ]: 114450 : if (M <= -1.0) { // supersonic inflow
856 : 0 : r = m_fr;
857 : 0 : ru = m_fr * m_fu[0];
858 : 0 : rv = m_fr * m_fu[1];
859 : 0 : rw = m_fr * m_fu[2];
860 : 0 : re = eos_totalenergy< eq >
861 [ - - ]: 0 : ( m_system, m_fr, m_fu[0], m_fu[1], m_fu[2], m_fp );
862 [ + - ][ + + ]: 114450 : } else if (M > -1.0 && M < 0.0) { // subsonic inflow
863 : 41547 : r = m_fr;
864 : 41547 : ru = m_fr * m_fu[0];
865 : 41547 : rv = m_fr * m_fu[1];
866 : 41547 : rw = m_fr * m_fu[2];
867 : 41547 : re =
868 [ + - ]: 41547 : eos_totalenergy< eq >( m_system, m_fr, m_fu[0], m_fu[1],
869 [ + - ]: 41547 : m_fu[2], eos_pressure< eq >( m_system, r, ru/r, rv/r,
870 : : rw/r, re ) );
871 [ + - ][ + - ]: 72903 : } else if (M >= 0.0 && M < 1.0) { // subsonic outflow
872 : 72903 : re = eos_totalenergy< eq >( m_system, r, ru/r, rv/r, rw/r,
873 [ + - ]: 72903 : m_fp );
874 : : }
875 : : }
876 : : }
877 : : }
878 : 68585 : }
879 : :
880 : : //! Apply sponge conditions at sponge nodes
881 : : //! \param[in] U Solution vector at recent time step
882 : : //! \param[in] coord Mesh node coordinates
883 : : //! \param[in] nodes Unique set of node ids at which to apply sponge
884 : : //! \details This function applies a sponge-like parameter to nodes of a
885 : : //! side set specified in the input file. We remove a user-specified
886 : : //! percentage of the kinetic energy by reducing the tangential
887 : : //! component of the velocity at a boundary and thereby modeling the
888 : : //! effect of a solid wall on the fluid via fluid-structure interaction
889 : : //! via a viscosity-like effect.
890 : : void
891 : 25081 : sponge( tk::Fields& U,
892 : : const std::array< std::vector< real >, 3 >& coord,
893 : : const std::unordered_set< std::size_t >& nodes ) const
894 : : {
895 : 25081 : const auto& x = coord[0];
896 : 25081 : const auto& y = coord[1];
897 : 25081 : const auto& z = coord[2];
898 : 25081 : const auto& sponge = g_inputdeck.get< param, eq, tag::sponge >();
899 : 25081 : const auto& ss = sponge.get< tag::sideset >();
900 [ + + ]: 25081 : if (ss.size() > m_system) { // sponge side set for this system
901 : 62 : const auto& spvel = sponge.get< tag::velocity >();
902 [ + + ]: 6324 : for (auto p : nodes) { // for all sponge nodes
903 [ + - ]: 6262 : if (!skipPoint(x[p],y[p],z[p])) {
904 [ + - ]: 12524 : std::vector< tk::real > sp( ss[m_system].size(), 0.0 );
905 [ + + ]: 6262 : if (spvel.size() > m_system) {
906 [ + - ]: 3131 : sp = spvel[m_system];
907 [ + + ]: 6262 : for (auto& s : sp) s = std::sqrt(s);
908 : : }
909 : : // sponge velocity: reduce kinetic energy by a user percentage
910 [ + + ]: 12524 : for (std::size_t s=0; s<ss[m_system].size(); ++s) {
911 [ + - ][ + - ]: 6262 : U(p,1,m_offset) -= U(p,1,m_offset)*sp[s];
912 [ + - ][ + - ]: 6262 : U(p,2,m_offset) -= U(p,2,m_offset)*sp[s];
913 [ + - ][ + - ]: 6262 : U(p,3,m_offset) -= U(p,3,m_offset)*sp[s];
914 : : }
915 : : }
916 : : }
917 : : }
918 : 25081 : }
919 : :
920 : : //! Apply user defined time dependent BCs
921 : : //! \param[in] t Physical time
922 : : //! \param[in,out] U Solution vector at recent time step
923 : : //! \param[in] nodes Vector of unique sets of node ids at which to apply BCs
924 : : //! \details This function applies user defined time dependent boundary
925 : : //! conditions on groups of side sets specified in the input file.
926 : : //! The user specifies pressure, density, and velocity as discrete
927 : : //! functions of time, in the control file, associated with a group of
928 : : //! side sets. Several such groups can be specified, each with their
929 : : //! own discrete function: p(t), rho(t), vx(t), vy(t), vz(t).
930 : : void
931 : 25081 : timedepbc( tk::real t,
932 : : tk::Fields& U,
933 : : const std::vector< std::unordered_set< std::size_t > >& nodes,
934 : : const std::vector< tk::Table<5> >& timedepfn ) const
935 : : {
936 [ + + ]: 25286 : for (std::size_t ib=0; ib<nodes.size(); ++ib) {
937 [ + + ]: 2460 : for (auto p:nodes[ib]) {
938 : : // sample primitive vars from discrete data at time t
939 [ + - ]: 2255 : auto unk = tk::sample<5>(t, timedepfn[ib]);
940 : :
941 : : // apply BCs after converting to conserved vars
942 [ + - ]: 2255 : U(p,0,m_offset) = unk[1];
943 [ + - ]: 2255 : U(p,1,m_offset) = unk[1]*unk[2];
944 [ + - ]: 2255 : U(p,2,m_offset) = unk[1]*unk[3];
945 [ + - ]: 2255 : U(p,3,m_offset) = unk[1]*unk[4];
946 [ + - ][ + - ]: 2255 : U(p,4,m_offset) = eos_totalenergy< eq >(m_system, unk[1], unk[2],
947 : : unk[3], unk[4], unk[0]);
948 : : }
949 : : }
950 : 25081 : }
951 : :
952 : : //! Return analytic field names to be output to file
953 : : //! \return Vector of strings labelling analytic fields output in file
954 : 2234 : std::vector< std::string > analyticFieldNames() const
955 : 2234 : { return m_problem.analyticFieldNames( m_ncomp ); }
956 : :
957 : : //! Return surface field names to be output to file
958 : : //! \return Vector of strings labelling surface fields output in file
959 : 2440 : std::vector< std::string > surfNames() const
960 : 2440 : { return CompFlowSurfNames(); }
961 : :
962 : : //! Return time history field names to be output to file
963 : : //! \return Vector of strings labelling time history fields output in file
964 : 28 : std::vector< std::string > histNames() const
965 : 28 : { return CompFlowHistNames(); }
966 : :
967 : : //! Return surface field output going to file
968 : : std::vector< std::vector< real > >
969 : 2440 : surfOutput( const std::map< int, std::vector< std::size_t > >& bnd,
970 : : const tk::Fields& U ) const
971 : 2440 : { return CompFlowSurfOutput( m_system, bnd, U ); }
972 : :
973 : : //! Return time history field output evaluated at time history points
974 : : std::vector< std::vector< real > >
975 : 898 : histOutput( const std::vector< HistData >& h,
976 : : const std::vector< std::size_t >& inpoel,
977 : : const tk::Fields& U ) const
978 : 898 : { return CompFlowHistOutput( m_system, h, inpoel, U ); }
979 : :
980 : : //! Return names of integral variables to be output to diagnostics file
981 : : //! \return Vector of strings labelling integral variables output
982 : 60 : std::vector< std::string > names() const
983 : 60 : { return m_problem.names( m_ncomp ); }
984 : :
985 : : private:
986 : : const Physics m_physics; //!< Physics policy
987 : : const Problem m_problem; //!< Problem policy
988 : : const ncomp_t m_system; //!< Equation system index
989 : : const ncomp_t m_offset; //!< Offset PDE operates from
990 : : //! Stagnation BC user configuration: point coordinates and radii
991 : : const std::tuple< std::vector< real >, std::vector< real > > m_stagCnf;
992 : : //! Skip BC user configuration: point coordinates and radii
993 : : const std::tuple< std::vector< real >, std::vector< real > > m_skipCnf;
994 : : const real m_fr; //!< Farfield density
995 : : const real m_fp; //!< Farfield pressure
996 : : const std::vector< real > m_fu; //!< Farfield velocity
997 : :
998 : : //! Decide if point is a stagnation point
999 : : //! \param[in] x X mesh point coordinates to query
1000 : : //! \param[in] y Y mesh point coordinates to query
1001 : : //! \param[in] z Z mesh point coordinates to query
1002 : : //! \return True if point is configured as a stagnation point by the user
1003 : : #pragma omp declare simd
1004 : : bool
1005 : 300647450 : stagPoint( real x, real y, real z ) const {
1006 : 300647450 : const auto& pnt = std::get< 0 >( m_stagCnf );
1007 : 300647450 : const auto& rad = std::get< 1 >( m_stagCnf );
1008 [ + + ]: 307129595 : for (std::size_t i=0; i<rad.size(); ++i) {
1009 [ + + ]: 6491192 : if (tk::length( x-pnt[i*3+0], y-pnt[i*3+1], z-pnt[i*3+2] ) < rad[i])
1010 : 9047 : return true;
1011 : : }
1012 : 300638403 : return false;
1013 : : }
1014 : :
1015 : : //! Decide if point is a skip-BC point
1016 : : //! \param[in] x X mesh point coordinates to query
1017 : : //! \param[in] y Y mesh point coordinates to query
1018 : : //! \param[in] z Z mesh point coordinates to query
1019 : : //! \return True if point is configured as a skip-BC point by the user
1020 : : #pragma omp declare simd
1021 : : bool
1022 : 302973895 : skipPoint( real x, real y, real z ) const {
1023 : 302973895 : const auto& pnt = std::get< 0 >( m_skipCnf );
1024 : 302973895 : const auto& rad = std::get< 1 >( m_skipCnf );
1025 [ - + ]: 302973895 : for (std::size_t i=0; i<rad.size(); ++i) {
1026 [ - - ]: 0 : if (tk::length( x-pnt[i*3+0], y-pnt[i*3+1], z-pnt[i*3+2] ) < rad[i])
1027 : 0 : return true;
1028 : : }
1029 : 302973895 : return false;
1030 : : }
1031 : :
1032 : : //! \brief Compute/assemble nodal gradients of primitive variables for
1033 : : //! ALECG in all points
1034 : : //! \param[in] coord Mesh node coordinates
1035 : : //! \param[in] inpoel Mesh element connectivity
1036 : : //! \param[in] lid Global->local node ids
1037 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
1038 : : //! global node ids (key)
1039 : : //! \param[in] vol Nodal volumes
1040 : : //! \param[in] esup Elements surrounding points
1041 : : //! \param[in] U Solution vector at recent time step
1042 : : //! \param[in] G Nodal gradients of primitive variables in chare-boundary
1043 : : //! nodes
1044 : : //! \return Gradients of primitive variables in all mesh points
1045 : : tk::Fields
1046 : 24936 : nodegrad( const std::array< std::vector< real >, 3 >& coord,
1047 : : const std::vector< std::size_t >& inpoel,
1048 : : const std::unordered_map< std::size_t, std::size_t >& lid,
1049 : : const std::unordered_map< std::size_t, std::size_t >& bid,
1050 : : const std::vector< real >& vol,
1051 : : const std::pair< std::vector< std::size_t >,
1052 : : std::vector< std::size_t > >& esup,
1053 : : const tk::Fields& U,
1054 : : const tk::Fields& G ) const
1055 : : {
1056 : : // allocate storage for nodal gradients of primitive variables
1057 : 24936 : tk::Fields Grad( U.nunk(), m_ncomp*3 );
1058 [ + - ]: 24936 : Grad.fill( 0.0 );
1059 : :
1060 : : // access node cooordinates
1061 : 24936 : const auto& x = coord[0];
1062 : 24936 : const auto& y = coord[1];
1063 : 24936 : const auto& z = coord[2];
1064 : :
1065 : : // compute gradients of primitive variables in points
1066 : 24936 : auto npoin = U.nunk();
1067 : : #pragma omp simd
1068 [ + + ]: 4042599 : for (std::size_t p=0; p<npoin; ++p)
1069 [ + + ]: 52581327 : for (auto e : tk::Around(esup,p)) {
1070 : : // access node IDs
1071 : 48563664 : std::size_t N[4] =
1072 : 48563664 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
1073 : : // compute element Jacobi determinant, J = 6V
1074 : 48563664 : real bax = x[N[1]]-x[N[0]];
1075 : 48563664 : real bay = y[N[1]]-y[N[0]];
1076 : 48563664 : real baz = z[N[1]]-z[N[0]];
1077 : 48563664 : real cax = x[N[2]]-x[N[0]];
1078 : 48563664 : real cay = y[N[2]]-y[N[0]];
1079 : 48563664 : real caz = z[N[2]]-z[N[0]];
1080 : 48563664 : real dax = x[N[3]]-x[N[0]];
1081 : 48563664 : real day = y[N[3]]-y[N[0]];
1082 : 48563664 : real daz = z[N[3]]-z[N[0]];
1083 : 48563664 : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
1084 : 48563664 : auto J24 = J/24.0;
1085 : : // shape function derivatives, nnode*ndim [4][3]
1086 : : real g[4][3];
1087 : 48563664 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
1088 : : g[1][0], g[1][1], g[1][2] );
1089 : 48563664 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
1090 : : g[2][0], g[2][1], g[2][2] );
1091 : 48563664 : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
1092 : : g[3][0], g[3][1], g[3][2] );
1093 [ + + ]: 194254656 : for (std::size_t i=0; i<3; ++i)
1094 : 145690992 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
1095 : : // scatter-add gradient contributions to boundary nodes
1096 : : real u[m_ncomp];
1097 [ + + ]: 242818320 : for (std::size_t b=0; b<4; ++b) {
1098 [ + - ]: 194254656 : u[0] = U(N[b],0,m_offset);
1099 [ + - ]: 194254656 : u[1] = U(N[b],1,m_offset)/u[0];
1100 [ + - ]: 194254656 : u[2] = U(N[b],2,m_offset)/u[0];
1101 [ + - ]: 194254656 : u[3] = U(N[b],3,m_offset)/u[0];
1102 [ + - ]: 194254656 : u[4] = U(N[b],4,m_offset)/u[0]
1103 : 194254656 : - 0.5*(u[1]*u[1] + u[2]*u[2] + u[3]*u[3]);
1104 [ + - ][ + + ]: 194254656 : if ( !skipPoint(x[N[b]],y[N[b]],z[N[b]]) &&
[ + + ]
1105 : : stagPoint(x[N[b]],y[N[b]],z[N[b]]) )
1106 : : {
1107 : 4272 : u[1] = u[2] = u[3] = 0.0;
1108 : : }
1109 [ + + ]: 1165527936 : for (std::size_t c=0; c<m_ncomp; ++c)
1110 [ + + ]: 3885093120 : for (std::size_t i=0; i<3; ++i)
1111 [ + - ]: 2913819840 : Grad(p,c*3+i,0) += J24 * g[b][i] * u[c];
1112 : : }
1113 : : }
1114 : :
1115 : : // put in nodal gradients of chare-boundary points
1116 [ + + ]: 1566117 : for (const auto& [g,b] : bid) {
1117 [ + - ]: 1541181 : auto i = tk::cref_find( lid, g );
1118 [ + + ]: 24658896 : for (ncomp_t c=0; c<Grad.nprop(); ++c)
1119 [ + - ][ + - ]: 23117715 : Grad(i,c,0) = G(b,c,0);
1120 : : }
1121 : :
1122 : : // divide weak result in gradients by nodal volume
1123 [ + + ]: 4042599 : for (std::size_t p=0; p<npoin; ++p)
1124 [ + + ]: 64282608 : for (std::size_t c=0; c<m_ncomp*3; ++c)
1125 [ + - ]: 60264945 : Grad(p,c,0) /= vol[p];
1126 : :
1127 : 24936 : return Grad;
1128 : : }
1129 : :
1130 : : //! Compute domain-edge integral for ALECG
1131 : : //! \param[in] coord Mesh node coordinates
1132 : : //! \param[in] gid Local->glocal node ids
1133 : : //! \param[in] edgenode Local node ids of edges
1134 : : //! \param[in] edgeid Local node id pair -> edge id map
1135 : : //! \param[in] psup Points surrounding points
1136 : : //! \param[in] dfn Dual-face normals
1137 : : //! \param[in] U Solution vector at recent time step
1138 : : //! \param[in] W Mesh velocity
1139 : : //! \param[in] G Nodal gradients
1140 : : //! \param[in] spmult Sponge pressure multiplers at nodes, one per symBC set
1141 : : //! \param[in,out] R Right-hand side vector computed
1142 : 24936 : void domainint( const std::array< std::vector< real >, 3 >& coord,
1143 : : const std::vector< std::size_t >& gid,
1144 : : const std::vector< std::size_t >& edgenode,
1145 : : const std::vector< std::size_t >& edgeid,
1146 : : const std::pair< std::vector< std::size_t >,
1147 : : std::vector< std::size_t > >& psup,
1148 : : const std::vector< real >& dfn,
1149 : : const tk::Fields& U,
1150 : : const tk::Fields& W,
1151 : : const tk::Fields& G,
1152 : : const std::vector< tk::real >& spmult,
1153 : : tk::Fields& R ) const
1154 : : {
1155 : : // domain-edge integral: compute fluxes in edges
1156 [ + - ]: 49872 : std::vector< real > dflux( edgenode.size()/2 * m_ncomp );
1157 : :
1158 : : // access node coordinates
1159 : 24936 : const auto& x = coord[0];
1160 : 24936 : const auto& y = coord[1];
1161 : 24936 : const auto& z = coord[2];
1162 : :
1163 : : // number of side sets configured with sponge pressure multipliers
1164 : 24936 : std::size_t nset = spmult.size() / x.size();
1165 : :
1166 : : #pragma omp simd
1167 [ + + ]: 19396659 : for (std::size_t e=0; e<edgenode.size()/2; ++e) {
1168 : 19371723 : auto p = edgenode[e*2+0];
1169 : 19371723 : auto q = edgenode[e*2+1];
1170 : :
1171 : : // compute primitive variables at edge-end points
1172 [ + - ]: 19371723 : real rL = U(p,0,m_offset);
1173 [ + - ]: 19371723 : real ruL = U(p,1,m_offset) / rL;
1174 [ + - ]: 19371723 : real rvL = U(p,2,m_offset) / rL;
1175 [ + - ]: 19371723 : real rwL = U(p,3,m_offset) / rL;
1176 [ + - ]: 19371723 : real reL = U(p,4,m_offset) / rL - 0.5*(ruL*ruL + rvL*rvL + rwL*rwL);
1177 [ + - ]: 19371723 : real w1L = W(p,0,0);
1178 [ + - ]: 19371723 : real w2L = W(p,1,0);
1179 [ + - ]: 19371723 : real w3L = W(p,2,0);
1180 [ + - ]: 19371723 : real rR = U(q,0,m_offset);
1181 [ + - ]: 19371723 : real ruR = U(q,1,m_offset) / rR;
1182 [ + - ]: 19371723 : real rvR = U(q,2,m_offset) / rR;
1183 [ + - ]: 19371723 : real rwR = U(q,3,m_offset) / rR;
1184 [ + - ]: 19371723 : real reR = U(q,4,m_offset) / rR - 0.5*(ruR*ruR + rvR*rvR + rwR*rwR);
1185 [ + - ]: 19371723 : real w1R = W(q,0,0);
1186 [ + - ]: 19371723 : real w2R = W(q,1,0);
1187 [ + - ]: 19371723 : real w3R = W(q,2,0);
1188 : :
1189 : : // apply stagnation BCs to primitive variables
1190 [ + - ][ + + ]: 19371723 : if ( !skipPoint(x[p],y[p],z[p]) && stagPoint(x[p],y[p],z[p]) )
[ + + ]
1191 : 1068 : ruL = rvL = rwL = 0.0;
1192 [ + - ][ + + ]: 19371723 : if ( !skipPoint(x[q],y[q],z[q]) && stagPoint(x[q],y[q],z[q]) )
[ + + ]
1193 : 534 : ruR = rvR = rwR = 0.0;
1194 : :
1195 : : // compute MUSCL reconstruction in edge-end points
1196 [ + - ]: 19371723 : muscl( p, q, coord, G,
1197 : : rL, ruL, rvL, rwL, reL, rR, ruR, rvR, rwR, reR );
1198 : :
1199 : : // convert back to conserved variables
1200 : 19371723 : reL = (reL + 0.5*(ruL*ruL + rvL*rvL + rwL*rwL)) * rL;
1201 : 19371723 : ruL *= rL;
1202 : 19371723 : rvL *= rL;
1203 : 19371723 : rwL *= rL;
1204 : 19371723 : reR = (reR + 0.5*(ruR*ruR + rvR*rvR + rwR*rwR)) * rR;
1205 : 19371723 : ruR *= rR;
1206 : 19371723 : rvR *= rR;
1207 : 19371723 : rwR *= rR;
1208 : :
1209 : : // evaluate pressure at edge-end points
1210 [ + - ]: 19371723 : real pL = eos_pressure<eq>( m_system, rL, ruL/rL, rvL/rL, rwL/rL, reL );
1211 [ + - ]: 19371723 : real pR = eos_pressure<eq>( m_system, rR, ruR/rR, rvR/rR, rwR/rR, reR );
1212 : :
1213 : : // apply sponge-pressure multipliers
1214 [ + + ]: 19515963 : for (std::size_t s=0; s<nset; ++s) {
1215 : 144240 : pL -= pL*spmult[p*nset+s];
1216 : 144240 : pR -= pR*spmult[q*nset+s];
1217 : : }
1218 : :
1219 : : // compute Riemann flux using edge-end point states
1220 : : real f[m_ncomp];
1221 [ + - ]: 19371723 : Rusanov::flux( m_system,
1222 : 19371723 : dfn[e*6+0], dfn[e*6+1], dfn[e*6+2],
1223 : 19371723 : dfn[e*6+3], dfn[e*6+4], dfn[e*6+5],
1224 : : rL, ruL, rvL, rwL, reL,
1225 : : rR, ruR, rvR, rwR, reR,
1226 : : w1L, w2L, w3L, w1R, w2R, w3R,
1227 : : pL, pR,
1228 : : f[0], f[1], f[2], f[3], f[4] );
1229 : : // store flux in edges
1230 [ + + ]: 116230338 : for (std::size_t c=0; c<m_ncomp; ++c) dflux[e*m_ncomp+c] = f[c];
1231 : : }
1232 : :
1233 : : // access pointer to right hand side at component and offset
1234 : : std::array< const real*, m_ncomp > r;
1235 [ + + ][ + - ]: 149616 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c, m_offset );
1236 : :
1237 : : // domain-edge integral: sum flux contributions to points
1238 [ + + ]: 4042599 : for (std::size_t p=0,k=0; p<U.nunk(); ++p)
1239 [ + + ]: 42761109 : for (auto q : tk::Around(psup,p)) {
1240 [ + + ]: 38743446 : auto s = gid[p] > gid[q] ? -1.0 : 1.0;
1241 : 38743446 : auto e = edgeid[k++];
1242 : : // the 2.0 in the following expression is so that the RHS contribution
1243 : : // conforms with Eq 12 (Waltz et al. Computers & fluids (92) 2014);
1244 : : // The 1/2 in Eq 12 is extracted from the flux function (Rusanov).
1245 : : // However, Rusanov::flux computes the flux with the 1/2. This 2
1246 : : // cancels with the 1/2 in Rusanov::flux, so that the 1/2 can be
1247 : : // extracted out and multiplied as in Eq 12
1248 [ + + ]: 232460676 : for (std::size_t c=0; c<m_ncomp; ++c)
1249 [ + - ]: 193717230 : R.var(r[c],p) -= 2.0*s*dflux[e*m_ncomp+c];
1250 : : }
1251 : :
1252 : 24936 : tk::destroy(dflux);
1253 : 24936 : }
1254 : :
1255 : : //! \brief Compute MUSCL reconstruction in edge-end points using a MUSCL
1256 : : //! procedure with van Leer limiting
1257 : : //! \param[in] p Left node id of edge-end
1258 : : //! \param[in] q Right node id of edge-end
1259 : : //! \param[in] coord Array of nodal coordinates
1260 : : //! \param[in] G Gradient of all unknowns in mesh points
1261 : : //! \param[in,out] rL Left density
1262 : : //! \param[in,out] uL Left X velocity
1263 : : //! \param[in,out] vL Left Y velocity
1264 : : //! \param[in,out] wL Left Z velocity
1265 : : //! \param[in,out] eL Left internal energy
1266 : : //! \param[in,out] rR Right density
1267 : : //! \param[in,out] uR Right X velocity
1268 : : //! \param[in,out] vR Right Y velocity
1269 : : //! \param[in,out] wR Right Z velocity
1270 : : //! \param[in,out] eR Right internal energy
1271 : 19371723 : void muscl( std::size_t p,
1272 : : std::size_t q,
1273 : : const tk::UnsMesh::Coords& coord,
1274 : : const tk::Fields& G,
1275 : : real& rL, real& uL, real& vL, real& wL, real& eL,
1276 : : real& rR, real& uR, real& vR, real& wR, real& eR ) const
1277 : : {
1278 : : // access node coordinates
1279 : 19371723 : const auto& x = coord[0];
1280 : 19371723 : const auto& y = coord[1];
1281 : 19371723 : const auto& z = coord[2];
1282 : :
1283 : : // edge vector
1284 : 19371723 : std::array< real, 3 > vw{ x[q]-x[p], y[q]-y[p], z[q]-z[p] };
1285 : :
1286 : : real delta1[5], delta2[5], delta3[5];
1287 : 19371723 : std::array< real, 5 > ls{ rL, uL, vL, wL, eL };
1288 : 19371723 : std::array< real, 5 > rs{ rR, uR, vR, wR, eR };
1289 : 19371723 : auto url = ls;
1290 : 19371723 : auto urr = rs;
1291 : :
1292 : : // MUSCL reconstruction of edge-end-point primitive variables
1293 [ + + ]: 116230338 : for (std::size_t c=0; c<5; ++c) {
1294 : : // gradients
1295 [ + - ][ + - ]: 96858615 : std::array< real, 3 > g1{ G(p,c*3+0,0), G(p,c*3+1,0), G(p,c*3+2,0) },
[ + - ]
1296 [ + - ][ + - ]: 96858615 : g2{ G(q,c*3+0,0), G(q,c*3+1,0), G(q,c*3+2,0) };
[ + - ]
1297 : :
1298 : 96858615 : delta2[c] = rs[c] - ls[c];
1299 : 96858615 : delta1[c] = 2.0 * tk::dot(g1,vw) - delta2[c];
1300 : 96858615 : delta3[c] = 2.0 * tk::dot(g2,vw) - delta2[c];
1301 : :
1302 : : // MUSCL extrapolation option 1:
1303 : : // ---------------------------------------------------------------------
1304 : : // Uncomment the following 3 blocks of code if this version is required.
1305 : : // this reconstruction is from the following paper:
1306 : : // Waltz, J., Morgan, N. R., Canfield, T. R., Charest, M. R.,
1307 : : // Risinger, L. D., & Wohlbier, J. G. (2014). A three-dimensional
1308 : : // finite element arbitrary Lagrangian–Eulerian method for shock
1309 : : // hydrodynamics on unstructured grids. Computers & Fluids, 92,
1310 : : // 172-187.
1311 : :
1312 : : //// form limiters
1313 : : //auto rcL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
1314 : : //auto rcR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
1315 : : //auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
1316 : : //auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
1317 : :
1318 : : //// van Leer limiter
1319 : : //// any other symmetric limiter could be used instead too
1320 : : //auto phiL = (std::abs(rcL) + rcL) / (std::abs(rcL) + 1.0);
1321 : : //auto phiR = (std::abs(rcR) + rcR) / (std::abs(rcR) + 1.0);
1322 : : //auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
1323 : : //auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
1324 : :
1325 : : //// update unknowns with reconstructed unknowns
1326 : : //url[c] += 0.25*(delta1[c]*muscl_m1*phiL + delta2[c]*muscl_p1*phi_L_inv);
1327 : : //urr[c] -= 0.25*(delta3[c]*muscl_m1*phiR + delta2[c]*muscl_p1*phi_R_inv);
1328 : :
1329 : : // ---------------------------------------------------------------------
1330 : :
1331 : : // MUSCL extrapolation option 2:
1332 : : // ---------------------------------------------------------------------
1333 : : // The following 2 blocks of code.
1334 : : // this reconstruction is from the following paper:
1335 : : // Luo, H., Baum, J. D., & Lohner, R. (1994). Edge-based finite element
1336 : : // scheme for the Euler equations. AIAA journal, 32(6), 1183-1190.
1337 : : // Van Leer, B. (1974). Towards the ultimate conservative difference
1338 : : // scheme. II. Monotonicity and conservation combined in a second-order
1339 : : // scheme. Journal of computational physics, 14(4), 361-370.
1340 : :
1341 : : // get Van Albada limiter
1342 : : // the following form is derived from the flux limiter phi as:
1343 : : // s = phi_inv - (1 - phi)
1344 : 193717230 : auto sL = std::max(0.0, (2.0*delta1[c]*delta2[c] + muscl_eps)
1345 : 96858615 : /(delta1[c]*delta1[c] + delta2[c]*delta2[c] + muscl_eps));
1346 : 193717230 : auto sR = std::max(0.0, (2.0*delta3[c]*delta2[c] + muscl_eps)
1347 : 96858615 : /(delta3[c]*delta3[c] + delta2[c]*delta2[c] + muscl_eps));
1348 : :
1349 : : // update unknowns with reconstructed unknowns
1350 : 96858615 : url[c] += 0.25*sL*(delta1[c]*(1.0-muscl_const*sL)
1351 : 96858615 : + delta2[c]*(1.0+muscl_const*sL));
1352 : 96858615 : urr[c] -= 0.25*sR*(delta3[c]*(1.0-muscl_const*sR)
1353 : 96858615 : + delta2[c]*(1.0+muscl_const*sR));
1354 : :
1355 : : // ---------------------------------------------------------------------
1356 : : }
1357 : :
1358 : : // force first order if the reconstructions for density or internal energy
1359 : : // would have allowed negative values
1360 [ + + ][ + + ]: 19371723 : if (ls[0] < delta1[0] || ls[4] < delta1[4]) url = ls;
[ + + ]
1361 [ + + ][ + + ]: 19371723 : if (rs[0] < -delta3[0] || rs[4] < -delta3[4]) urr = rs;
[ + + ]
1362 : :
1363 : 19371723 : rL = url[0];
1364 : 19371723 : uL = url[1];
1365 : 19371723 : vL = url[2];
1366 : 19371723 : wL = url[3];
1367 : 19371723 : eL = url[4];
1368 : :
1369 : 19371723 : rR = urr[0];
1370 : 19371723 : uR = urr[1];
1371 : 19371723 : vR = urr[2];
1372 : 19371723 : wR = urr[3];
1373 : 19371723 : eR = urr[4];
1374 : 19371723 : }
1375 : :
1376 : : //! Compute boundary integrals for ALECG
1377 : : //! \param[in] coord Mesh node coordinates
1378 : : //! \param[in] triinpoel Boundary triangle face connecitivity with local ids
1379 : : //! \param[in] symbctri Vector with 1 at symmetry BC boundary triangles
1380 : : //! \param[in] U Solution vector at recent time step
1381 : : //! \param[in] W Mesh velocity
1382 : : //! \param[in] spmult Sponge pressure multiplers at nodes, one per symBC set
1383 : : //! \param[in,out] R Right-hand side vector computed
1384 : 24936 : void bndint( const std::array< std::vector< real >, 3 >& coord,
1385 : : const std::vector< std::size_t >& triinpoel,
1386 : : const std::vector< int >& symbctri,
1387 : : const tk::Fields& U,
1388 : : const tk::Fields& W,
1389 : : const std::vector< tk::real >& spmult,
1390 : : tk::Fields& R ) const
1391 : : {
1392 : :
1393 : : // access node coordinates
1394 : 24936 : const auto& x = coord[0];
1395 : 24936 : const auto& y = coord[1];
1396 : 24936 : const auto& z = coord[2];
1397 : :
1398 : : // boundary integrals: compute fluxes in edges
1399 [ + - ]: 49872 : std::vector< real > bflux( triinpoel.size() * m_ncomp * 2 );
1400 : :
1401 : : // number of side sets configured with sponge pressure multipliers
1402 : 24936 : std::size_t nset = spmult.size() / x.size();
1403 : :
1404 : : #pragma omp simd
1405 [ + + ]: 4083852 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
1406 : : // access node IDs
1407 : 4058916 : std::size_t N[3] =
1408 : 4058916 : { triinpoel[e*3+0], triinpoel[e*3+1], triinpoel[e*3+2] };
1409 : : // access solution at element nodes
1410 [ + - ]: 4058916 : real rA = U(N[0],0,m_offset);
1411 [ + - ]: 4058916 : real rB = U(N[1],0,m_offset);
1412 [ + - ]: 4058916 : real rC = U(N[2],0,m_offset);
1413 [ + - ]: 4058916 : real ruA = U(N[0],1,m_offset);
1414 [ + - ]: 4058916 : real ruB = U(N[1],1,m_offset);
1415 [ + - ]: 4058916 : real ruC = U(N[2],1,m_offset);
1416 [ + - ]: 4058916 : real rvA = U(N[0],2,m_offset);
1417 [ + - ]: 4058916 : real rvB = U(N[1],2,m_offset);
1418 [ + - ]: 4058916 : real rvC = U(N[2],2,m_offset);
1419 [ + - ]: 4058916 : real rwA = U(N[0],3,m_offset);
1420 [ + - ]: 4058916 : real rwB = U(N[1],3,m_offset);
1421 [ + - ]: 4058916 : real rwC = U(N[2],3,m_offset);
1422 [ + - ]: 4058916 : real reA = U(N[0],4,m_offset);
1423 [ + - ]: 4058916 : real reB = U(N[1],4,m_offset);
1424 [ + - ]: 4058916 : real reC = U(N[2],4,m_offset);
1425 [ + - ]: 4058916 : real w1A = W(N[0],0,0);
1426 [ + - ]: 4058916 : real w2A = W(N[0],1,0);
1427 [ + - ]: 4058916 : real w3A = W(N[0],2,0);
1428 [ + - ]: 4058916 : real w1B = W(N[1],0,0);
1429 [ + - ]: 4058916 : real w2B = W(N[1],1,0);
1430 [ + - ]: 4058916 : real w3B = W(N[1],2,0);
1431 [ + - ]: 4058916 : real w1C = W(N[2],0,0);
1432 [ + - ]: 4058916 : real w2C = W(N[2],1,0);
1433 [ + - ]: 4058916 : real w3C = W(N[2],2,0);
1434 : : // apply stagnation BCs
1435 [ + - ][ + + ]: 4058916 : if ( !skipPoint(x[N[0]],y[N[0]],z[N[0]]) &&
[ + + ]
1436 : : stagPoint(x[N[0]],y[N[0]],z[N[0]]) )
1437 : : {
1438 : 1335 : ruA = rvA = rwA = 0.0;
1439 : : }
1440 [ + - ][ - + ]: 4058916 : if ( !skipPoint(x[N[1]],y[N[1]],z[N[1]]) &&
[ - + ]
1441 : : stagPoint(x[N[1]],y[N[1]],z[N[1]]) )
1442 : : {
1443 : 0 : ruB = rvB = rwB = 0.0;
1444 : : }
1445 [ + - ][ + + ]: 4058916 : if ( !skipPoint(x[N[2]],y[N[2]],z[N[2]]) &&
[ + + ]
1446 : : stagPoint(x[N[2]],y[N[2]],z[N[2]]) )
1447 : : {
1448 : 267 : ruC = rvC = rwC = 0.0;
1449 : : }
1450 : : // compute face normal
1451 : : real nx, ny, nz;
1452 : 4058916 : tk::normal( x[N[0]], x[N[1]], x[N[2]],
1453 : : y[N[0]], y[N[1]], y[N[2]],
1454 : : z[N[0]], z[N[1]], z[N[2]],
1455 : : nx, ny, nz );
1456 : : // compute boundary flux
1457 : : real f[m_ncomp][3];
1458 : : real p, vn;
1459 : 4058916 : int sym = symbctri[e];
1460 [ + - ]: 4058916 : p = eos_pressure< eq >( m_system, rA, ruA/rA, rvA/rA, rwA/rA, reA );
1461 [ + + ]: 4106676 : for (std::size_t s=0; s<nset; ++s) p -= p*spmult[N[0]*nset+s];
1462 [ + + ]: 4058916 : vn = sym ? 0.0 : (nx*(ruA/rA-w1A) + ny*(rvA/rA-w2A) + nz*(rwA/rA-w3A));
1463 : 4058916 : f[0][0] = rA*vn;
1464 : 4058916 : f[1][0] = ruA*vn + p*nx;
1465 : 4058916 : f[2][0] = rvA*vn + p*ny;
1466 : 4058916 : f[3][0] = rwA*vn + p*nz;
1467 [ + + ]: 4058916 : f[4][0] = reA*vn + p*(sym ? 0.0 : (nx*ruA + ny*rvA + nz*rwA)/rA);
1468 [ + - ]: 4058916 : p = eos_pressure< eq >( m_system, rB, ruB/rB, rvB/rB, rwB/rB, reB );
1469 [ + + ]: 4106676 : for (std::size_t s=0; s<nset; ++s) p -= p*spmult[N[1]*nset+s];
1470 [ + + ]: 4058916 : vn = sym ? 0.0 : (nx*(ruB/rB-w1B) + ny*(rvB/rB-w2B) + nz*(rwB/rB-w3B));
1471 : 4058916 : f[0][1] = rB*vn;
1472 : 4058916 : f[1][1] = ruB*vn + p*nx;
1473 : 4058916 : f[2][1] = rvB*vn + p*ny;
1474 : 4058916 : f[3][1] = rwB*vn + p*nz;
1475 [ + + ]: 4058916 : f[4][1] = reB*vn + p*(sym ? 0.0 : (nx*ruB + ny*rvB + nz*rwB)/rB);
1476 [ + - ]: 4058916 : p = eos_pressure< eq >( m_system, rC, ruC/rC, rvC/rC, rwC/rC, reC );
1477 [ + + ]: 4106676 : for (std::size_t s=0; s<nset; ++s) p -= p*spmult[N[2]*nset+s];
1478 [ + + ]: 4058916 : vn = sym ? 0.0 : (nx*(ruC/rC-w1C) + ny*(rvC/rC-w2C) + nz*(rwC/rC-w3C));
1479 : 4058916 : f[0][2] = rC*vn;
1480 : 4058916 : f[1][2] = ruC*vn + p*nx;
1481 : 4058916 : f[2][2] = rvC*vn + p*ny;
1482 : 4058916 : f[3][2] = rwC*vn + p*nz;
1483 [ + + ]: 4058916 : f[4][2] = reC*vn + p*(sym ? 0.0 : (nx*ruC + ny*rvC + nz*rwC)/rC);
1484 : : // compute face area
1485 : 4058916 : auto A6 = tk::area( x[N[0]], x[N[1]], x[N[2]],
1486 : : y[N[0]], y[N[1]], y[N[2]],
1487 : : z[N[0]], z[N[1]], z[N[2]] ) / 6.0;
1488 : 4058916 : auto A24 = A6/4.0;
1489 : : // store flux in boundary elements
1490 [ + + ]: 24353496 : for (std::size_t c=0; c<m_ncomp; ++c) {
1491 : 20294580 : auto eb = (e*m_ncomp+c)*6;
1492 : 20294580 : auto Bab = A24 * (f[c][0] + f[c][1]);
1493 : 20294580 : bflux[eb+0] = Bab + A6 * f[c][0];
1494 : 20294580 : bflux[eb+1] = Bab;
1495 : 20294580 : Bab = A24 * (f[c][1] + f[c][2]);
1496 : 20294580 : bflux[eb+2] = Bab + A6 * f[c][1];
1497 : 20294580 : bflux[eb+3] = Bab;
1498 : 20294580 : Bab = A24 * (f[c][2] + f[c][0]);
1499 : 20294580 : bflux[eb+4] = Bab + A6 * f[c][2];
1500 : 20294580 : bflux[eb+5] = Bab;
1501 : : }
1502 : : }
1503 : :
1504 : : // access pointer to right hand side at component and offset
1505 : : std::array< const real*, m_ncomp > r;
1506 [ + + ][ + - ]: 149616 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c, m_offset );
1507 : :
1508 : : // boundary integrals: sum flux contributions to points
1509 [ + + ]: 4083852 : for (std::size_t e=0; e<triinpoel.size()/3; ++e)
1510 [ + + ]: 24353496 : for (std::size_t c=0; c<m_ncomp; ++c) {
1511 : 20294580 : auto eb = (e*m_ncomp+c)*6;
1512 [ + - ]: 20294580 : R.var(r[c],triinpoel[e*3+0]) -= bflux[eb+0] + bflux[eb+5];
1513 [ + - ]: 20294580 : R.var(r[c],triinpoel[e*3+1]) -= bflux[eb+1] + bflux[eb+2];
1514 [ + - ]: 20294580 : R.var(r[c],triinpoel[e*3+2]) -= bflux[eb+3] + bflux[eb+4];
1515 : : }
1516 : :
1517 : 24936 : tk::destroy(bflux);
1518 : 24936 : }
1519 : :
1520 : : //! Compute optional source integral
1521 : : //! \param[in] coord Mesh node coordinates
1522 : : //! \param[in] inpoel Mesh element connectivity
1523 : : //! \param[in] t Physical time
1524 : : //! \param[in] tp Physical time for each mesh node
1525 : : //! \param[in,out] R Right-hand side vector computed
1526 : 24936 : void src( const std::array< std::vector< real >, 3 >& coord,
1527 : : const std::vector< std::size_t >& inpoel,
1528 : : real t,
1529 : : const std::vector< tk::real >& tp,
1530 : : tk::Fields& R ) const
1531 : : {
1532 : : // access node coordinates
1533 : 24936 : const auto& x = coord[0];
1534 : 24936 : const auto& y = coord[1];
1535 : 24936 : const auto& z = coord[2];
1536 : :
1537 : : // access pointer to right hand side at component and offset
1538 : : std::array< const real*, m_ncomp > r;
1539 [ + + ][ + - ]: 149616 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c, m_offset );
1540 : :
1541 : : // source integral
1542 [ + + ]: 12165852 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
1543 : 12140916 : std::size_t N[4] =
1544 : 12140916 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
1545 : : // compute element Jacobi determinant, J = 6V
1546 : 12140916 : auto J24 = tk::triple(
1547 : 12140916 : x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]],
1548 : 12140916 : x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]],
1549 : 12140916 : x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] ) / 24.0;
1550 : : // sum source contributions to nodes
1551 [ + + ]: 60704580 : for (std::size_t a=0; a<4; ++a) {
1552 : : real s[m_ncomp];
1553 [ + + ]: 48563664 : if (g_inputdeck.get< tag::discr, tag::steady_state >()) t = tp[N[a]];
1554 [ + - ]: 48563664 : Problem::src( m_system, x[N[a]], y[N[a]], z[N[a]], t,
1555 : : s[0], s[1], s[2], s[3], s[4] );
1556 [ + + ]: 291381984 : for (std::size_t c=0; c<m_ncomp; ++c)
1557 [ + - ]: 242818320 : R.var(r[c],N[a]) += J24 * s[c];
1558 : : }
1559 : : }
1560 : 24936 : }
1561 : :
1562 : : //! Compute sources corresponding to a propagating front in user-defined box
1563 : : //! \param[in] V Total box volume
1564 : : //! \param[in] t Physical time
1565 : : //! \param[in] inpoel Element point connectivity
1566 : : //! \param[in] esup Elements surrounding points
1567 : : //! \param[in] boxnodes Mesh node ids within user-defined box
1568 : : //! \param[in] coord Mesh node coordinates
1569 : : //! \param[in] R Right-hand side vector
1570 : : //! \details This function add the energy source corresponding to a planar
1571 : : //! wave-front propagating along the z-direction with a user-specified
1572 : : //! velocity, within a box initial condition, configured by the user.
1573 : : //! Example (SI) units of the quantities involved:
1574 : : //! * internal energy content (energy per unit volume): J/m^3
1575 : : //! * specific energy (internal energy per unit mass): J/kg
1576 : 1200 : void boxSrc( real V,
1577 : : real t,
1578 : : const std::vector< std::size_t >& inpoel,
1579 : : const std::pair< std::vector< std::size_t >,
1580 : : std::vector< std::size_t > >& esup,
1581 : : const std::unordered_set< std::size_t >& boxnodes,
1582 : : const std::array< std::vector< real >, 3 >& coord,
1583 : : tk::Fields& R ) const
1584 : : {
1585 : 1200 : const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
1586 : 1200 : const auto& icbox = ic.get< tag::box >();
1587 : :
1588 [ + - ]: 1200 : if (icbox.size() > m_system) {
1589 [ + + ]: 2400 : for (const auto& b : icbox[m_system]) { // for all boxes for this eq
1590 [ + - ]: 2400 : std::vector< tk::real > box
1591 : 1200 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
1592 : 1200 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
1593 : 1200 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
1594 : :
1595 : 1200 : auto boxenc = b.template get< tag::energy_content >();
1596 [ - + ][ - - ]: 1200 : Assert( boxenc > 0.0, "Box energy content must be nonzero" );
[ - - ][ - - ]
1597 : :
1598 : 1200 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
1599 : :
1600 : : // determine times at which sourcing is initialized and terminated
1601 : 1200 : auto iv = b.template get< tag::initiate, tag::velocity >();
1602 : 1200 : auto wFront = 0.08;
1603 : 1200 : auto tInit = 0.0;
1604 : 1200 : auto tFinal = tInit + (box[5] - box[4] - 2.0*wFront) / std::fabs(iv);
1605 : 1200 : auto aBox = (box[1]-box[0]) * (box[3]-box[2]);
1606 : :
1607 : 1200 : const auto& x = coord[0];
1608 : 1200 : const auto& y = coord[1];
1609 : 1200 : const auto& z = coord[2];
1610 : :
1611 [ + - ][ + + ]: 1200 : if (t >= tInit && t <= tFinal) {
1612 : : // The energy front is assumed to have a half-sine-wave shape. The
1613 : : // half wave-length is the width of the front. At t=0, the center of
1614 : : // this front (i.e. the peak of the partial-sine-wave) is at X_0 +
1615 : : // W_0. W_0 is calculated based on the width of the front and the
1616 : : // direction of propagation (which is assumed to be along the
1617 : : // z-direction). If the front propagation velocity is positive, it
1618 : : // is assumed that the initial position of the energy source is the
1619 : : // minimum z-coordinate of the box; whereas if this velocity is
1620 : : // negative, the initial position is the maximum z-coordinate of the
1621 : : // box.
1622 : :
1623 : : // initial center of front
1624 : 1080 : tk::real zInit(box[4]);
1625 [ + - ]: 1080 : if (iv < 0.0) zInit = box[5];
1626 : : // current location of front
1627 : 1080 : auto z0 = zInit + iv*t;
1628 : 1080 : auto z1 = z0 + std::copysign(wFront, iv);
1629 : 1080 : tk::real s0(z0), s1(z1);
1630 : : // if velocity of propagation is negative, initial position is z1
1631 [ + - ]: 1080 : if (iv < 0.0) {
1632 : 1080 : s0 = z1;
1633 : 1080 : s1 = z0;
1634 : : }
1635 : : // Sine-wave (positive part of the wave) source term amplitude
1636 : 1080 : auto pi = 4.0 * std::atan(1.0);
1637 : 1080 : auto amplE = boxenc * V_ex * pi
1638 : 1080 : / (aBox * wFront * 2.0 * (tFinal-tInit));
1639 : : //// Square wave (constant) source term amplitude
1640 : : //auto amplE = boxenc * V_ex
1641 : : // / (aBox * wFront * (tFinal-tInit));
1642 : 1080 : amplE *= V_ex / V;
1643 : :
1644 : : // add source
1645 [ + + ]: 73440 : for (auto p : boxnodes) {
1646 [ + + ][ + + ]: 72360 : if (z[p] >= s0 && z[p] <= s1) {
[ + + ]
1647 : 928 : auto S = amplE * std::sin(pi*(z[p]-s0)/wFront);
1648 [ + + ]: 13285 : for (auto e : tk::Around(esup,p)) {
1649 : : // access node IDs
1650 : 12357 : std::size_t N[4] =
1651 : 12357 : {inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3]};
1652 : : // compute element Jacobi determinant, J = 6V
1653 : 12357 : real bax = x[N[1]]-x[N[0]];
1654 : 12357 : real bay = y[N[1]]-y[N[0]];
1655 : 12357 : real baz = z[N[1]]-z[N[0]];
1656 : 12357 : real cax = x[N[2]]-x[N[0]];
1657 : 12357 : real cay = y[N[2]]-y[N[0]];
1658 : 12357 : real caz = z[N[2]]-z[N[0]];
1659 : 12357 : real dax = x[N[3]]-x[N[0]];
1660 : 12357 : real day = y[N[3]]-y[N[0]];
1661 : 12357 : real daz = z[N[3]]-z[N[0]];
1662 : : auto J =
1663 : 12357 : tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
1664 : 12357 : auto J24 = J/24.0;
1665 [ + - ]: 12357 : R(p,4,m_offset) += J24 * S;
1666 : : }
1667 : : }
1668 : : }
1669 : : }
1670 : : }
1671 : : }
1672 : 1200 : }
1673 : :
1674 : : //! Compute sponge pressure multiplers at sponge side sets
1675 : : //! \param[in] coord Mesh node coordinates
1676 : : //! \param[in] nodes Unique set of nodes for sponge conditions
1677 : : //! \return Sponge ressure multiplers at nodes, one per sponge side set
1678 : : //! \details This function computes a sponge-like multiplier that will be
1679 : : //! applied to nodes of side sets specified in the input file. This is
1680 : : //! used to reduce the pressure gradient normal to boundaries and thereby
1681 : : //! modeling the effect of a solid wall on the fluid via fluid-structure
1682 : : //! interaction.
1683 : : //! \note If no sponge pressure coefficients are configured, an empty
1684 : : //! vector is returned.
1685 : : std::vector< tk::real >
1686 : 24936 : spongePressures( const std::array< std::vector< real >, 3 >& coord,
1687 : : const std::unordered_set< std::size_t >& nodes ) const
1688 : : {
1689 : 24936 : const auto& x = coord[0];
1690 : 24936 : const auto& y = coord[1];
1691 : 24936 : const auto& z = coord[2];
1692 : 24936 : std::vector< tk::real > spmult;
1693 : 24936 : std::size_t nset = 0; // number of sponge side sets configured
1694 : 24936 : const auto& sponge = g_inputdeck.get< param, eq, tag::sponge >();
1695 : 24936 : const auto& ss = sponge.get< tag::sideset >();
1696 [ + + ]: 24936 : if (ss.size() > m_system) { // if symbcs configured for this system
1697 : 60 : const auto& sppre = sponge.get< tag::pressure >();
1698 : 60 : nset = ss[m_system].size(); // number of sponge side sets configured
1699 [ + - ]: 60 : spmult.resize( x.size() * nset, 0.0 );
1700 [ + + ]: 6120 : for (auto p : nodes) {
1701 [ + - ][ + + ]: 6060 : if (not skipPoint(x[p],y[p],z[p]) && sppre.size() > m_system) {
[ + + ]
1702 [ - + ][ - - ]: 3030 : Assert( nset == sppre[m_system].size(), "Size mismatch" );
[ - - ][ - - ]
1703 [ + + ]: 6060 : for (std::size_t s=0; s<nset; ++s)
1704 : 3030 : spmult[p*nset+s] = sppre[m_system][s];
1705 : : } else {
1706 [ + + ]: 6060 : for (std::size_t s=0; s<nset; ++s)
1707 : 3030 : spmult[p*nset+s] = 0.0;
1708 : : }
1709 : : }
1710 : : }
1711 [ + + ][ - + ]: 24936 : Assert( ss.size() > m_system ?
[ - - ][ - - ]
[ - - ]
1712 : : spmult.size() == x.size() * ss[m_system].size() :
1713 : : spmult.size() == 0, "Sponge pressure multipler wrong size" );
1714 : 24936 : return spmult;
1715 : : }
1716 : :
1717 : : };
1718 : :
1719 : : } // cg::
1720 : :
1721 : : } // inciter::
1722 : :
1723 : : #endif // CGCompFlow_h
|
