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