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