Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Transport/CGTransport.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 Scalar transport using continous Galerkin discretization
9 : : \details This file implements the physics operators governing transported
10 : : scalars using continuous Galerkin discretization.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef CGTransport_h
14 : : #define CGTransport_h
15 : :
16 : : #include <vector>
17 : : #include <array>
18 : : #include <limits>
19 : : #include <cmath>
20 : : #include <unordered_set>
21 : : #include <unordered_map>
22 : :
23 : : #include "Exception.hpp"
24 : : #include "Vector.hpp"
25 : : #include "DerivedData.hpp"
26 : : #include "Around.hpp"
27 : : #include "Reconstruction.hpp"
28 : : #include "Inciter/InputDeck/InputDeck.hpp"
29 : : #include "CGPDE.hpp"
30 : : #include "History.hpp"
31 : :
32 : : namespace inciter {
33 : :
34 : : extern ctr::InputDeck g_inputdeck;
35 : :
36 : : namespace cg {
37 : :
38 : : //! \brief Transport equation used polymorphically with tk::CGPDE
39 : : //! \details The template argument(s) specify policies and are used to configure
40 : : //! the behavior of the class. The policies are:
41 : : //! - Physics - physics configuration, see PDE/Transport/Physics/CG.h
42 : : //! - Problem - problem configuration, see PDE/Transport/Problem.h
43 : : //! \note The default physics is CGAdvection, set in
44 : : //! inciter::deck::check_transport()
45 : : template< class Physics, class Problem >
46 [ - - ][ - - ]: 824 : class Transport {
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ + - ]
[ + - ][ - - ]
47 : :
48 : : private:
49 : : using ncomp_t = tk::ncomp_t;
50 : : using real = tk::real;
51 : : using eq = tag::transport;
52 : :
53 : : static constexpr real muscl_eps = 1.0e-9;
54 : : static constexpr real muscl_const = 1.0/3.0;
55 : : static constexpr real muscl_m1 = 1.0 - muscl_const;
56 : : static constexpr real muscl_p1 = 1.0 + muscl_const;
57 : :
58 : : public:
59 : : //! Constructor
60 : 0 : explicit Transport() :
61 : : m_physics( Physics() ),
62 : : m_problem( Problem() ),
63 [ - - ]: 206 : m_ncomp(g_inputdeck.get< tag::ncomp >())
64 : : {
65 [ - - ]: 0 : m_problem.errchk( m_ncomp );
66 : 0 : }
67 : :
68 : : //! Determine nodes that lie inside the user-defined IC box
69 : : void
70 : : IcBoxNodes( const tk::UnsMesh::Coords&,
71 : : const std::vector< std::size_t >&,
72 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&,
73 : : std::vector< std::unordered_set< std::size_t > >&,
74 : : std::unordered_map< std::size_t, std::set< std::size_t > >&,
75 : : std::size_t& ) const {}
76 : :
77 : : //! Initalize the transport equations using problem policy
78 : : //! \param[in] coord Mesh node coordinates
79 : : //! \param[in,out] unk Array of unknowns
80 : : //! \param[in] t Physical time
81 : : void
82 : 856 : initialize( const std::array< std::vector< real >, 3 >& coord,
83 : : tk::Fields& unk,
84 : : real t,
85 : : real,
86 : : const std::vector< std::unordered_set< std::size_t > >&,
87 : : const std::vector< tk::real >&,
88 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&
89 : : ) const
90 : : {
91 : : Assert( coord[0].size() == unk.nunk(), "Size mismatch" );
92 : : const auto& x = coord[0];
93 : : const auto& y = coord[1];
94 : : const auto& z = coord[2];
95 [ + + ]: 86428 : for (ncomp_t i=0; i<x.size(); ++i) {
96 : 85572 : auto s = Problem::initialize( m_ncomp, m_mat_blk, x[i], y[i],
97 : : z[i], t );
98 [ + + ]: 171144 : for (ncomp_t c=0; c<m_ncomp; ++c)
99 : 85572 : unk( i, c ) = s[c];
100 : : }
101 : 856 : }
102 : :
103 : : //! Query a velocity
104 : : //! \note Since this function does not touch its output argument, that
105 : : //! means this system does not define a "velocity".
106 : : void velocity( const tk::Fields&, tk::UnsMesh::Coords& ) const {}
107 : :
108 : : //! Query the sound speed
109 : : //! \note Since this function does not touch its output argument, that
110 : : //! means this system does not define a "sound speed".
111 : : void soundspeed( const tk::Fields&, std::vector< tk::real >& ) const {}
112 : :
113 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
114 : : //! \param[in] xi X-coordinate
115 : : //! \param[in] yi Y-coordinate
116 : : //! \param[in] zi Z-coordinate
117 : : //! \param[in] t Physical time
118 : : //! \return Vector of analytic solution at given location and time
119 : : std::vector< real >
120 : : analyticSolution( real xi, real yi, real zi, real t ) const
121 : 86092 : { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
122 : :
123 : : //! Return analytic solution for conserved variables
124 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
125 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
126 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
127 : : //! \param[in] t Physical time at which to evaluate the analytic solution
128 : : //! \return Vector of analytic solution at given location and time
129 : : std::vector< tk::real >
130 : : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
131 : 687904 : { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
132 : :
133 : : //! Compute nodal gradients of primitive variables for ALECG
134 : : //! \param[in] coord Mesh node coordinates
135 : : //! \param[in] inpoel Mesh element connectivity
136 : : //! \param[in] bndel List of elements contributing to chare-boundary nodes
137 : : //! \param[in] gid Local->global node id map
138 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
139 : : //! global node ids (key)
140 : : //! \param[in] U Solution vector at recent time step
141 : : //! \param[in,out] G Nodal gradients of primitive variables
142 : 30495 : void chBndGrad( const std::array< std::vector< real >, 3 >& coord,
143 : : const std::vector< std::size_t >& inpoel,
144 : : const std::vector< std::size_t >& bndel,
145 : : const std::vector< std::size_t >& gid,
146 : : const std::unordered_map< std::size_t, std::size_t >& bid,
147 : : const tk::Fields& U,
148 : : tk::Fields& G ) const
149 : : {
150 : : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
151 : : "vector at recent time step incorrect" );
152 : :
153 : : // compute gradients of primitive variables in points
154 : : G.fill( 0.0 );
155 : :
156 : : // access node cooordinates
157 : : const auto& x = coord[0];
158 : : const auto& y = coord[1];
159 : : const auto& z = coord[2];
160 : :
161 [ + + ]: 1808349 : for (auto e : bndel) { // elements contributing to chare boundary nodes
162 : : // access node IDs
163 : 1777854 : std::size_t N[4] =
164 : : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
165 : : // compute element Jacobi determinant, J = 6V
166 : 1777854 : real bax = x[N[1]]-x[N[0]];
167 : 1777854 : real bay = y[N[1]]-y[N[0]];
168 : 1777854 : real baz = z[N[1]]-z[N[0]];
169 : 1777854 : real cax = x[N[2]]-x[N[0]];
170 : 1777854 : real cay = y[N[2]]-y[N[0]];
171 : 1777854 : real caz = z[N[2]]-z[N[0]];
172 : 1777854 : real dax = x[N[3]]-x[N[0]];
173 : 1777854 : real day = y[N[3]]-y[N[0]];
174 : 1777854 : real daz = z[N[3]]-z[N[0]];
175 : : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
176 : 1777854 : auto J24 = J/24.0;
177 : : // shape function derivatives, nnode*ndim [4][3]
178 : : real g[4][3];
179 : : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
180 : : g[1][0], g[1][1], g[1][2] );
181 : : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
182 : : g[2][0], g[2][1], g[2][2] );
183 : : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
184 : : g[3][0], g[3][1], g[3][2] );
185 [ + + ]: 7111416 : for (std::size_t i=0; i<3; ++i)
186 : 5333562 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
187 : : // scatter-add gradient contributions to boundary nodes
188 [ + + ]: 8889270 : for (std::size_t a=0; a<4; ++a) {
189 : 7111416 : auto i = bid.find( gid[N[a]] );
190 [ + + ]: 7111416 : if (i != end(bid))
191 [ + + ]: 12407514 : for (std::size_t c=0; c<m_ncomp; ++c)
192 [ + + ]: 31018785 : for (std::size_t b=0; b<4; ++b)
193 [ + + ]: 99260112 : for (std::size_t j=0; j<3; ++j)
194 : 74445084 : G(i->second,c*3+j) += J24 * g[b][j] * U(N[b],c);
195 : : }
196 : : }
197 : 30495 : }
198 : :
199 : : //! Compute right hand side for ALECG
200 : : //! \param[in] coord Mesh node coordinates
201 : : //! \param[in] inpoel Mesh element connectivity
202 : : //! \param[in] triinpoel Boundary triangle face connecitivity
203 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
204 : : //! global node ids (key)
205 : : //! \param[in] lid Global->local node ids
206 : : //! \param[in] dfn Dual-face normals
207 : : //! \param[in] psup Points surrounding points
208 : : //! \param[in] esup Elements surrounding points
209 : : //! \param[in] symbctri Vector with 1 at symmetry BC nodes
210 : : //! \param[in] vol Nodal volumes
211 : : //! \param[in] edgeid Local node id pair -> edge id map
212 : : //! \param[in] G Nodal gradients in chare-boundary nodes
213 : : //! \param[in] U Solution vector at recent time step
214 : : //! \param[in] W Mesh velocity
215 : : //! \param[in,out] R Right-hand side vector computed
216 : 30495 : void rhs(
217 : : real,
218 : : const std::array< std::vector< real >, 3 >& coord,
219 : : const std::vector< std::size_t >& inpoel,
220 : : const std::vector< std::size_t >& triinpoel,
221 : : const std::vector< std::size_t >&,
222 : : const std::unordered_map< std::size_t, std::size_t >& bid,
223 : : const std::unordered_map< std::size_t, std::size_t >& lid,
224 : : const std::vector< real >& dfn,
225 : : const std::pair< std::vector< std::size_t >,
226 : : std::vector< std::size_t > >& psup,
227 : : const std::pair< std::vector< std::size_t >,
228 : : std::vector< std::size_t > >& esup,
229 : : const std::vector< int >& symbctri,
230 : : const std::vector< int >&,
231 : : const std::vector< real >& vol,
232 : : const std::vector< std::size_t >&,
233 : : const std::vector< std::size_t >& edgeid,
234 : : const std::vector< std::unordered_set< std::size_t > >&,
235 : : const tk::Fields& G,
236 : : const tk::Fields& U,
237 : : [[maybe_unused]] const tk::Fields& W,
238 : : const std::vector< tk::real >&,
239 : : real,
240 : : tk::Fields& R ) const
241 : : {
242 : : Assert( G.nprop() == m_ncomp*3,
243 : : "Number of components in gradient vector incorrect" );
244 : : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
245 : : "vector at recent time step incorrect" );
246 : : Assert( R.nunk() == coord[0].size(),
247 : : "Number of unknowns and/or number of components in right-hand "
248 : : "side vector incorrect" );
249 : :
250 : : // compute/assemble gradients in points
251 : 30495 : auto Grad = nodegrad( coord, inpoel, lid, bid, vol, esup, U, G );
252 : :
253 : : // zero right hand side for all components
254 [ + + ]: 60990 : for (ncomp_t c=0; c<m_ncomp; ++c) R.fill( c, 0.0 );
255 : :
256 : : // compute domain-edge integral
257 [ + - ]: 30495 : domainint( coord, inpoel, edgeid, psup, dfn, U, Grad, R );
258 : :
259 : : // compute boundary integrals
260 [ + - ]: 30495 : bndint( coord, triinpoel, symbctri, U, R );
261 : 30495 : }
262 : :
263 : : //! Compute boundary pressure integrals (force) (no-op for transport)
264 : : void bndPressureInt(
265 : : const std::array< std::vector< real >, 3 >&,
266 : : const std::vector< std::size_t >&,
267 : : const std::vector< int >&,
268 : : const tk::Fields&,
269 : : const std::array< tk::real, 3 >&,
270 : : std::vector< real >& ) const
271 : : { }
272 : :
273 : :
274 : : //! Compute the minimum time step size (for unsteady time stepping)
275 : : //! \param[in] U Solution vector at recent time step
276 : : //! \param[in] coord Mesh node coordinates
277 : : //! \param[in] inpoel Mesh element connectivity
278 : : //! \param[in] t Physical time
279 : : //! \return Minimum time step size
280 : 75 : real dt( const std::array< std::vector< real >, 3 >& coord,
281 : : const std::vector< std::size_t >& inpoel,
282 : : tk::real t,
283 : : tk::real,
284 : : const tk::Fields& U,
285 : : const std::vector< tk::real >&,
286 : : const std::vector< tk::real >& ) const
287 : : {
288 : : using tag::transport;
289 : : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
290 : : "vector at recent time step incorrect" );
291 : : const auto& x = coord[0];
292 : : const auto& y = coord[1];
293 : : const auto& z = coord[2];
294 : : // compute the minimum dt across all elements we own
295 : : auto mindt = std::numeric_limits< tk::real >::max();
296 : : auto eps = std::numeric_limits< tk::real >::epsilon();
297 : : auto large = std::numeric_limits< tk::real >::max();
298 [ + + ]: 180507 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
299 : : const std::array< std::size_t, 4 >
300 : 180432 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
301 : : // compute cubic root of element volume as the characteristic length
302 : : const std::array< real, 3 >
303 : 180432 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
304 : 180432 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
305 : 180432 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
306 : 180432 : const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
307 : : // access solution at element nodes at recent time step
308 : 180432 : std::vector< std::array< real, 4 > > u( m_ncomp );
309 [ + + ]: 360864 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, N );
310 : : // get velocity for problem
311 [ + - ][ + - ]: 180432 : const std::array< std::vector<std::array<real,3>>, 4 > vel{{
[ + - ][ + - ]
312 [ + - ]: 180432 : Problem::prescribedVelocity( m_ncomp,
313 : : x[N[0]], y[N[0]], z[N[0]], t ),
314 [ + - ]: 180432 : Problem::prescribedVelocity( m_ncomp,
315 : : x[N[1]], y[N[1]], z[N[1]], t ),
316 [ + - ]: 180432 : Problem::prescribedVelocity( m_ncomp,
317 : : x[N[2]], y[N[2]], z[N[2]], t ),
318 [ + - ]: 180432 : Problem::prescribedVelocity( m_ncomp,
319 : : x[N[3]], y[N[3]], z[N[3]], t ) }};
320 : : // compute the maximum length of the characteristic velocity (advection
321 : : // velocity) across the four element nodes
322 : : real maxvel = 0.0;
323 [ + + ]: 360864 : for (ncomp_t c=0; c<m_ncomp; ++c)
324 [ + + ]: 902160 : for (std::size_t i=0; i<4; ++i) {
325 [ + + ]: 721728 : auto v = tk::length( vel[i][c] );
326 [ + + ]: 721728 : if (v > maxvel) maxvel = v;
327 : : }
328 : : // compute element dt for the advection
329 [ + - ]: 180432 : auto advection_dt = std::abs(maxvel) > eps ? L / maxvel : large;
330 : : // compute element dt based on diffusion
331 [ - + ][ - - ]: 180432 : auto diffusion_dt = m_physics.diffusion_dt( m_ncomp, L, u );
332 : : // compute minimum element dt
333 : 180432 : auto elemdt = std::min( advection_dt, diffusion_dt );
334 : : // find minimum dt across all elements
335 [ + + ]: 180432 : if (elemdt < mindt) mindt = elemdt;
336 : : }
337 : 75 : return mindt * g_inputdeck.get< tag::cfl >();
338 : : }
339 : :
340 : : //! Compute a time step size for each mesh node (for steady time stepping)
341 : : void dt( uint64_t,
342 : : const std::vector< tk::real >&,
343 : : const tk::Fields&,
344 : : std::vector< tk::real >& ) const {}
345 : :
346 : : //! \brief Query Dirichlet boundary condition value on a given side set for
347 : : //! all components in this PDE system
348 : : //! \param[in] t Physical time
349 : : //! \param[in] deltat Time step size
350 : : //! \param[in] tp Physical time for each mesh node
351 : : //! \param[in] dtp Time step size for each mesh node
352 : : //! \param[in] ss Pair of side set ID and list of node IDs on the side set
353 : : //! \param[in] coord Mesh node coordinates
354 : : //! \param[in] increment If true, evaluate the solution increment between
355 : : //! t and t+dt for Dirichlet BCs. If false, evlauate the solution instead.
356 : : //! \return Vector of pairs of bool and boundary condition value associated
357 : : //! to mesh node IDs at which Dirichlet boundary conditions are set. Note
358 : : //! that if increment is true, instead of the actual boundary condition
359 : : //! value, we return the increment between t+deltat and t, since,
360 : : //! depending on client code and solver, that may be what the solution
361 : : //! requires.
362 : : std::map< std::size_t, std::vector< std::pair<bool,real> > >
363 [ - + ]: 12835 : dirbc( real t,
364 : : real deltat,
365 : : const std::vector< tk::real >& tp,
366 : : const std::vector< tk::real >& dtp,
367 : : const std::pair< const int, std::vector< std::size_t > >& ss,
368 : : const std::array< std::vector< real >, 3 >& coord,
369 : : bool increment ) const
370 : : {
371 : : using NodeBC = std::vector< std::pair< bool, real > >;
372 : : std::map< std::size_t, NodeBC > bc;
373 : 12835 : const auto& ubc = g_inputdeck.get< tag::bc >()[0].get< tag::dirichlet >();
374 [ - + ]: 12835 : const auto steady = g_inputdeck.get< tag::steady_state >();
375 [ - + ]: 12835 : if (!ubc.empty()) {
376 : : Assert( ubc.size() > 0, "Indexing out of Dirichlet BC eq-vector" );
377 : : const auto& x = coord[0];
378 : : const auto& y = coord[1];
379 : : const auto& z = coord[2];
380 [ - - ]: 0 : for (const auto& b : ubc)
381 [ - - ]: 0 : if (static_cast<int>(b) == ss.first)
382 [ - - ]: 0 : for (auto n : ss.second) {
383 : : Assert( x.size() > n, "Indexing out of coordinate array" );
384 [ - - ]: 0 : if (steady) { t = tp[n]; deltat = dtp[n]; }
385 [ - - ][ - - ]: 0 : const auto s = increment ?
[ - - ]
386 [ - - ]: 0 : solinc( m_ncomp, m_mat_blk, x[n], y[n], z[n],
387 : : t, deltat, Problem::initialize ) :
388 [ - - ]: 0 : Problem::initialize( m_ncomp, m_mat_blk, x[n], y[n],
389 : : z[n], t+deltat );
390 [ - - ][ - - ]: 0 : auto& nbc = bc[n] = NodeBC( m_ncomp );
391 [ - - ]: 0 : for (ncomp_t c=0; c<m_ncomp; ++c)
392 : 0 : nbc[c] = { true, s[c] };
393 : : }
394 : : }
395 : 12835 : return bc;
396 : : }
397 : :
398 : : //! Set symmetry boundary conditions at nodes
399 : : void
400 : : symbc(
401 : : tk::Fields&,
402 : : const std::array< std::vector< real >, 3 >&,
403 : : const std::unordered_map< int,
404 : : std::unordered_map< std::size_t,
405 : : std::array< real, 4 > > >&,
406 : : const std::unordered_set< std::size_t >& ) const {}
407 : :
408 : : //! Set farfield boundary conditions at nodes
409 : : void farfieldbc(
410 : : tk::Fields&,
411 : : const std::array< std::vector< real >, 3 >&,
412 : : const std::unordered_map< int,
413 : : std::unordered_map< std::size_t,
414 : : std::array< real, 4 > > >&,
415 : : const std::unordered_set< std::size_t >& ) const {}
416 : :
417 : : //! Set slip wall boundary conditions at nodes
418 : : void
419 : : slipwallbc(
420 : : tk::Fields&,
421 : : const tk::Fields&,
422 : : const std::array< std::vector< real >, 3 >&,
423 : : const std::unordered_map< int,
424 : : std::unordered_map< std::size_t,
425 : : std::array< real, 4 > > >&,
426 : : const std::unordered_set< std::size_t >& ) const {}
427 : :
428 : : //! Apply user defined time dependent BCs (no-op for transport)
429 : : void
430 : : timedepbc( tk::real,
431 : : tk::Fields&,
432 : : const std::vector< std::unordered_set< std::size_t > >&,
433 : : const std::vector< tk::Table<5> >& ) const {}
434 : :
435 : : //! Return a map that associates user-specified strings to functions
436 : : //! \return Map that associates user-specified strings to functions that
437 : : //! compute relevant quantities to be output to file
438 [ + - ]: 601 : std::map< std::string, tk::GetVarFn > OutVarFn() const {
439 : : std::map< std::string, tk::GetVarFn > OutFnMap;
440 [ + - ][ + - ]: 601 : OutFnMap["material_indicator"] = transport::matIndicatorOutVar;
441 : :
442 : 601 : return OutFnMap;
443 : : }
444 : :
445 : : //! Return analytic field names to be output to file
446 : : //! \return Vector of strings labelling analytic fields output in file
447 : 565 : std::vector< std::string > analyticFieldNames() const {
448 : : std::vector< std::string > n;
449 : 565 : auto depvar = g_inputdeck.get< tag::depvar >()[0];
450 [ + + ]: 1130 : for (ncomp_t c=0; c<m_ncomp; ++c)
451 [ + - ][ + - ]: 1130 : n.push_back( depvar + std::to_string(c) + "_analytic" );
[ - + ][ - + ]
[ - - ][ - - ]
452 : 565 : return n;
453 : : }
454 : :
455 : : //! Return surface field names to be output to file
456 : : //! \return Vector of strings labelling surface fields output in file
457 : : //! \details This functions should be written in conjunction with
458 : : //! surfOutput(), which provides the vector of surface fields to be output
459 : : std::vector< std::string > surfNames() const { return {}; }
460 : :
461 : : //! Return nodal surface field output going to file
462 : : std::vector< std::vector< real > >
463 : : surfOutput( const std::map< int, std::vector< std::size_t > >&,
464 : : const tk::Fields& ) const { return {}; }
465 : :
466 : : //! Return elemental surface field output (on triangle faces) going to file
467 : : std::vector< std::vector< real > >
468 : : elemSurfOutput( const std::map< int, std::vector< std::size_t > >&,
469 : : const std::vector< std::size_t >&,
470 : : const tk::Fields& ) const { return {}; }
471 : :
472 : : //! Return time history field names to be output to file
473 : : //! \return Vector of strings labelling time history fields output in file
474 : : std::vector< std::string > histNames() const { return {}; }
475 : :
476 : : //! Return time history field output evaluated at time history points
477 : : std::vector< std::vector< real > >
478 : : histOutput( const std::vector< HistData >&,
479 : : const std::vector< std::size_t >&,
480 : : const tk::Fields& ) const { return {}; }
481 : :
482 : : //! Return names of integral variables to be output to diagnostics file
483 : : //! \return Vector of strings labelling integral variables output
484 [ + - ]: 52 : std::vector< std::string > names() const {
485 : : std::vector< std::string > n;
486 : : const auto& depvar =
487 : : g_inputdeck.get< tag::depvar >().at(0);
488 : : // construct the name of the numerical solution for all components
489 [ + + ]: 104 : for (ncomp_t c=0; c<m_ncomp; ++c)
490 [ + - ][ - + ]: 104 : n.push_back( depvar + std::to_string(c) );
[ - - ]
491 : 52 : return n;
492 : : }
493 : :
494 : : private:
495 : : const Physics m_physics; //!< Physics policy
496 : : const Problem m_problem; //!< Problem policy
497 : : const ncomp_t m_ncomp; //!< Number of components in this PDE
498 : : //! EOS material block
499 : : const std::vector< EOS > m_mat_blk;
500 : :
501 : : //! \brief Compute/assemble nodal gradients of primitive variables for
502 : : //! ALECG in all points
503 : : //! \param[in] coord Mesh node coordinates
504 : : //! \param[in] inpoel Mesh element connectivity
505 : : //! \param[in] lid Global->local node ids
506 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
507 : : //! global node ids (key)
508 : : //! \param[in] vol Nodal volumes
509 : : //! \param[in] esup Elements surrounding points
510 : : //! \param[in] U Solution vector at recent time step
511 : : //! \param[in] G Nodal gradients of primitive variables in chare-boundary nodes
512 : : //! \return Gradients of primitive variables in all mesh points
513 : : tk::Fields
514 : 30495 : nodegrad( const std::array< std::vector< real >, 3 >& coord,
515 : : const std::vector< std::size_t >& inpoel,
516 : : const std::unordered_map< std::size_t, std::size_t >& lid,
517 : : const std::unordered_map< std::size_t, std::size_t >& bid,
518 : : const std::vector< real >& vol,
519 : : const std::pair< std::vector< std::size_t >,
520 : : std::vector< std::size_t > >& esup,
521 : : const tk::Fields& U,
522 : : const tk::Fields& G ) const
523 : : {
524 : : // allocate storage for nodal gradients of primitive variables
525 : 30495 : tk::Fields Grad( U.nunk(), m_ncomp*3 );
526 : : Grad.fill( 0.0 );
527 : :
528 : : // access node cooordinates
529 : : const auto& x = coord[0];
530 : : const auto& y = coord[1];
531 : : const auto& z = coord[2];
532 : :
533 : : // compute gradients of primitive variables in points
534 : 30495 : auto npoin = U.nunk();
535 : : #pragma omp simd
536 [ + + ]: 2094207 : for (std::size_t p=0; p<npoin; ++p)
537 [ + + ]: 16125696 : for (auto e : tk::Around(esup,p)) {
538 : : // access node IDs
539 : 14061984 : std::size_t N[4] =
540 : : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
541 : : // compute element Jacobi determinant, J = 6V
542 : 14061984 : real bax = x[N[1]]-x[N[0]];
543 : 14061984 : real bay = y[N[1]]-y[N[0]];
544 : 14061984 : real baz = z[N[1]]-z[N[0]];
545 : 14061984 : real cax = x[N[2]]-x[N[0]];
546 : 14061984 : real cay = y[N[2]]-y[N[0]];
547 : 14061984 : real caz = z[N[2]]-z[N[0]];
548 : 14061984 : real dax = x[N[3]]-x[N[0]];
549 : 14061984 : real day = y[N[3]]-y[N[0]];
550 : 14061984 : real daz = z[N[3]]-z[N[0]];
551 : : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
552 : 14061984 : auto J24 = J/24.0;
553 : : // shape function derivatives, nnode*ndim [4][3]
554 : : real g[4][3];
555 : : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
556 : : g[1][0], g[1][1], g[1][2] );
557 : : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
558 : : g[2][0], g[2][1], g[2][2] );
559 : : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
560 : : g[3][0], g[3][1], g[3][2] );
561 [ + + ]: 56247936 : for (std::size_t i=0; i<3; ++i)
562 : 42185952 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
563 : : // scatter-add gradient contributions to boundary nodes
564 [ + + ]: 28123968 : for (std::size_t c=0; c<m_ncomp; ++c)
565 [ + + ]: 70309920 : for (std::size_t b=0; b<4; ++b)
566 [ + + ]: 224991744 : for (std::size_t i=0; i<3; ++i)
567 : 168743808 : Grad(p,c*3+i) += J24 * g[b][i] * U(N[b],c);
568 : : }
569 : :
570 : : // put in nodal gradients of chare-boundary points
571 [ + + ]: 1435941 : for (const auto& [g,b] : bid) {
572 : 1405446 : auto i = tk::cref_find( lid, g );
573 [ + + ]: 5621784 : for (ncomp_t c=0; c<Grad.nprop(); ++c)
574 : 4216338 : Grad(i,c) = G(b,c);
575 : : }
576 : :
577 : : // divide weak result in gradients by nodal volume
578 [ + + ]: 2094207 : for (std::size_t p=0; p<npoin; ++p)
579 [ + + ]: 8254848 : for (std::size_t c=0; c<m_ncomp*3; ++c)
580 : 6191136 : Grad(p,c) /= vol[p];
581 : :
582 : 30495 : return Grad;
583 : : }
584 : :
585 : : //! \brief Compute MUSCL reconstruction in edge-end points using a MUSCL
586 : : //! procedure with van Leer limiting
587 : : //! \param[in] p Left node id of edge-end
588 : : //! \param[in] q Right node id of edge-end
589 : : //! \param[in] coord Array of nodal coordinates
590 : : //! \param[in] G Gradient of all unknowns in mesh points
591 : : //! \param[in,out] uL Primitive variables at left edge-end point
592 : : //! \param[in,out] uR Primitive variables at right edge-end point
593 : : void
594 : 15738216 : muscl( std::size_t p,
595 : : std::size_t q,
596 : : const tk::UnsMesh::Coords& coord,
597 : : const tk::Fields& G,
598 : : std::vector< real >& uL,
599 : : std::vector< real >& uR ) const
600 : : {
601 : : Assert( uL.size() == m_ncomp && uR.size() == m_ncomp, "Size mismatch" );
602 : : Assert( G.nprop()/3 == m_ncomp, "Size mismatch" );
603 : :
604 : : const auto& x = coord[0];
605 : : const auto& y = coord[1];
606 : : const auto& z = coord[2];
607 : :
608 : : // edge vector
609 : 15738216 : std::array< real, 3 > vw{ x[q]-x[p], y[q]-y[p], z[q]-z[p] };
610 : :
611 : : std::vector< real >
612 [ + - ][ + - ]: 15738216 : delta1( m_ncomp, 0.0 ), delta2( m_ncomp, 0.0 ), delta3( m_ncomp, 0.0 );
[ - - ][ - - ]
613 : :
614 : : // MUSCL reconstruction of edge-end-point primitive variables
615 [ + + ]: 31476432 : for (std::size_t c=0; c<m_ncomp; ++c) {
616 : : // gradients
617 : : std::array< real, 3 >
618 : 15738216 : g1{ G(p,c*3+0), G(p,c*3+1), G(p,c*3+2) },
619 : 15738216 : g2{ G(q,c*3+0), G(q,c*3+1), G(q,c*3+2) };
620 : :
621 : 15738216 : delta2[c] = uR[c] - uL[c];
622 : 15738216 : delta1[c] = 2.0 * tk::dot(g1,vw) - delta2[c];
623 : 15738216 : delta3[c] = 2.0 * tk::dot(g2,vw) - delta2[c];
624 : :
625 : : // form limiters
626 : 15738216 : auto rL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
627 : 15738216 : auto rR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
628 : 15738216 : auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
629 : 15738216 : auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
630 : :
631 : 15738216 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
632 : 15738216 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
633 : 15738216 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
634 : 15738216 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
635 : :
636 : : // update unknowns with reconstructed unknowns
637 : 15738216 : uL[c] += 0.25*(delta1[c]*muscl_m1*phiL + delta2[c]*muscl_p1*phi_L_inv);
638 : 15738216 : uR[c] -= 0.25*(delta3[c]*muscl_m1*phiR + delta2[c]*muscl_p1*phi_R_inv);
639 : : }
640 : 15738216 : }
641 : :
642 : : //! Compute domain-edge integral for ALECG
643 : : //! \param[in] coord Mesh node coordinates
644 : : //! \param[in] inpoel Mesh element connectivity
645 : : //! \param[in] edgeid Local node id pair -> edge id map
646 : : //! \param[in] psup Points surrounding points
647 : : //! \param[in] dfn Dual-face normals
648 : : //! \param[in] U Solution vector at recent time step
649 : : //! \param[in] G Nodal gradients
650 : : //! \param[in,out] R Right-hand side vector computed
651 : 30495 : void domainint( const std::array< std::vector< real >, 3 >& coord,
652 : : const std::vector< std::size_t >& inpoel,
653 : : const std::vector< std::size_t >& edgeid,
654 : : const std::pair< std::vector< std::size_t >,
655 : : std::vector< std::size_t > >& psup,
656 : : const std::vector< real >& dfn,
657 : : const tk::Fields& U,
658 : : const tk::Fields& G,
659 : : tk::Fields& R ) const
660 : : {
661 : : // access node cooordinates
662 : : const auto& x = coord[0];
663 : : const auto& y = coord[1];
664 : : const auto& z = coord[2];
665 : :
666 : : // compute derived data structures
667 [ + - ]: 30495 : auto esued = tk::genEsued( inpoel, 4, tk::genEsup( inpoel, 4 ) );
668 : :
669 : : // access pointer to right hand side at component
670 [ + - ]: 30495 : std::vector< const real* > r( m_ncomp );
671 [ + + ]: 60990 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c );
672 : :
673 : : // domain-edge integral
674 [ + + ]: 2094207 : for (std::size_t p=0,k=0; p<U.nunk(); ++p) {
675 [ + + ]: 17801928 : for (auto q : tk::Around(psup,p)) {
676 : : // access dual-face normals for edge p-q
677 [ + - ]: 15738216 : auto ed = edgeid[k++];
678 : 15738216 : std::array< tk::real, 3 > n{ dfn[ed*6+0], dfn[ed*6+1], dfn[ed*6+2] };
679 : :
680 [ + - ]: 15738216 : std::vector< tk::real > uL( m_ncomp, 0.0 );
681 [ + - ][ - - ]: 15738216 : std::vector< tk::real > uR( m_ncomp, 0.0 );
682 [ + + ]: 31476432 : for (std::size_t c=0; c<m_ncomp; ++c) {
683 : 15738216 : uL[c] = U(p,c);
684 : 15738216 : uR[c] = U(q,c);
685 : : }
686 : : // compute MUSCL reconstruction in edge-end points
687 [ + - ]: 15738216 : muscl( p, q, coord, G, uL, uR );
688 : :
689 : : // evaluate prescribed velocity
690 [ + - ]: 15738216 : auto v =
691 [ + - ]: 15738216 : Problem::prescribedVelocity( m_ncomp, x[p], y[p], z[p], 0.0 );
692 : : // sum donain-edge contributions
693 [ + + ]: 57924168 : for (auto e : tk::cref_find(esued,{p,q})) {
694 : : const std::array< std::size_t, 4 >
695 : 42185952 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
696 : : // compute element Jacobi determinant
697 : : const std::array< tk::real, 3 >
698 : 42185952 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
699 : 42185952 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
700 : 42185952 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
701 : : const auto J = tk::triple( ba, ca, da ); // J = 6V
702 : : // shape function derivatives, nnode*ndim [4][3]
703 : : std::array< std::array< tk::real, 3 >, 4 > grad;
704 : 42185952 : grad[1] = tk::crossdiv( ca, da, J );
705 : 42185952 : grad[2] = tk::crossdiv( da, ba, J );
706 : 42185952 : grad[3] = tk::crossdiv( ba, ca, J );
707 [ + + ]: 168743808 : for (std::size_t i=0; i<3; ++i)
708 : 126557856 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
709 : 42185952 : auto J48 = J/48.0;
710 [ + + ]: 295301664 : for (const auto& [a,b] : tk::lpoed) {
711 [ + - ][ - - ]: 253115712 : auto s = tk::orient( {N[a],N[b]}, {p,q} );
712 [ + + ]: 1012462848 : for (std::size_t j=0; j<3; ++j) {
713 [ + + ]: 1518694272 : for (std::size_t c=0; c<m_ncomp; ++c) {
714 : 759347136 : R.var(r[c],p) -= J48 * s * (grad[a][j] - grad[b][j])
715 : 759347136 : * v[c][j]*(uL[c] + uR[c])
716 : 759347136 : - J48 * std::abs(s * (grad[a][j] - grad[b][j]))
717 : 759347136 : * std::abs(tk::dot(v[c],n)) * (uR[c] - uL[c]);
718 : : }
719 : : }
720 : : }
721 : : }
722 : : }
723 : : }
724 : 30495 : }
725 : :
726 : : //! Compute boundary integrals for ALECG
727 : : //! \param[in] coord Mesh node coordinates
728 : : //! \param[in] triinpoel Boundary triangle face connecitivity with local ids
729 : : //! \param[in] symbctri Vector with 1 at symmetry BC boundary triangles
730 : : //! \param[in] U Solution vector at recent time step
731 : : //! \param[in,out] R Right-hand side vector computed
732 : 30495 : void bndint( const std::array< std::vector< real >, 3 >& coord,
733 : : const std::vector< std::size_t >& triinpoel,
734 : : const std::vector< int >& symbctri,
735 : : const tk::Fields& U,
736 : : tk::Fields& R ) const
737 : : {
738 : : // access node coordinates
739 : : const auto& x = coord[0];
740 : : const auto& y = coord[1];
741 : : const auto& z = coord[2];
742 : :
743 : : // boundary integrals: compute fluxes in edges
744 : 30495 : std::vector< real > bflux( triinpoel.size() * m_ncomp * 2 );
745 : :
746 [ + + ]: 1277295 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
747 : : // access node IDs
748 : : std::array< std::size_t, 3 >
749 [ + - ]: 1246800 : N{ triinpoel[e*3+0], triinpoel[e*3+1], triinpoel[e*3+2] };
750 : : // apply symmetry BCs
751 [ + - ]: 1246800 : if (symbctri[e]) continue;
752 : : // node coordinates
753 [ - - ]: 0 : std::array< tk::real, 3 > xp{ x[N[0]], x[N[1]], x[N[2]] },
754 : 0 : yp{ y[N[0]], y[N[1]], y[N[2]] },
755 : 0 : zp{ z[N[0]], z[N[1]], z[N[2]] };
756 : : // access solution at element nodes
757 [ - - ]: 0 : std::vector< std::array< real, 3 > > u( m_ncomp );
758 [ - - ]: 0 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, N );
759 : : // evaluate prescribed velocity
760 [ - - ]: 0 : auto v =
761 : : Problem::prescribedVelocity( m_ncomp, xp[0], yp[0], zp[0], 0.0 );
762 : : // compute face area
763 : 0 : auto A6 = tk::area( x[N[0]], x[N[1]], x[N[2]],
764 : : y[N[0]], y[N[1]], y[N[2]],
765 : : z[N[0]], z[N[1]], z[N[2]] ) / 6.0;
766 : 0 : auto A24 = A6/4.0;
767 : : // compute face normal
768 [ - - ]: 0 : auto n = tk::normal( xp, yp, zp );
769 : : // store flux in boundary elements
770 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
771 : 0 : auto eb = (e*m_ncomp+c)*6;
772 : 0 : auto vdotn = tk::dot( v[c], n );
773 : 0 : auto Bab = A24 * vdotn * (u[c][0] + u[c][1]);
774 : 0 : bflux[eb+0] = Bab + A6 * vdotn * u[c][0];
775 : 0 : bflux[eb+1] = Bab;
776 : 0 : Bab = A24 * vdotn * (u[c][1] + u[c][2]);
777 : 0 : bflux[eb+2] = Bab + A6 * vdotn * u[c][1];
778 : 0 : bflux[eb+3] = Bab;
779 : 0 : Bab = A24 * vdotn * (u[c][2] + u[c][0]);
780 : 0 : bflux[eb+4] = Bab + A6 * vdotn * u[c][2];
781 : 0 : bflux[eb+5] = Bab;
782 : : }
783 : : }
784 : :
785 : : // access pointer to right hand side at component
786 [ + - ][ - - ]: 30495 : std::vector< const real* > r( m_ncomp );
787 [ + + ]: 60990 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c );
788 : :
789 : : // boundary integrals: sum flux contributions to points
790 [ + + ]: 1277295 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
791 : 1246800 : std::size_t N[3] =
792 : : { triinpoel[e*3+0], triinpoel[e*3+1], triinpoel[e*3+2] };
793 [ + + ]: 2493600 : for (std::size_t c=0; c<m_ncomp; ++c) {
794 : 1246800 : auto eb = (e*m_ncomp+c)*6;
795 : 1246800 : R.var(r[c],N[0]) -= bflux[eb+0] + bflux[eb+5];
796 : 1246800 : R.var(r[c],N[1]) -= bflux[eb+1] + bflux[eb+2];
797 : 1246800 : R.var(r[c],N[2]) -= bflux[eb+3] + bflux[eb+4];
798 : : }
799 : : }
800 : :
801 : : tk::destroy(bflux);
802 : 30495 : }
803 : : };
804 : :
805 : : } // cg::
806 : : } // inciter::
807 : :
808 : : #endif // Transport_h
|
