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 : : 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 : 206 : explicit Transport() :
61 : : m_physics( Physics() ),
62 : : m_problem( Problem() ),
63 : 206 : m_ncomp(g_inputdeck.get< tag::ncomp >())
64 : : {
65 [ - - ]: 206 : m_problem.errchk( m_ncomp );
66 : 206 : }
67 : :
68 : : //! Determine nodes that lie inside the user-defined IC box
69 : : void
70 : 690 : 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 : 690 : 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 [ - + ][ - - ]: 856 : Assert( coord[0].size() == unk.nunk(), "Size mismatch" );
[ - - ][ - - ]
92 : 856 : const auto& x = coord[0];
93 : 856 : const auto& y = coord[1];
94 : 856 : const auto& z = coord[2];
95 [ + + ]: 86428 : for (ncomp_t i=0; i<x.size(); ++i) {
96 [ + - ]: 171144 : 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 : 22133 : 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 : 22133 : 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 : 86092 : 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 : 687904 : 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 [ - + ][ - - ]: 30495 : 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 : 30495 : G.fill( 0.0 );
155 : :
156 : : // access node cooordinates
157 : 30495 : const auto& x = coord[0];
158 : 30495 : const auto& y = coord[1];
159 : 30495 : 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 : 1777854 : { 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 : 1777854 : 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 : 1777854 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
180 : : g[1][0], g[1][1], g[1][2] );
181 : 1777854 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
182 : : g[2][0], g[2][1], g[2][2] );
183 : 1777854 : 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< real >& vol,
231 : : const std::vector< std::size_t >&,
232 : : const std::vector< std::size_t >& edgeid,
233 : : const std::vector< std::unordered_set< std::size_t > >&,
234 : : const tk::Fields& G,
235 : : const tk::Fields& U,
236 : : [[maybe_unused]] const tk::Fields& W,
237 : : const std::vector< tk::real >&,
238 : : real,
239 : : tk::Fields& R ) const
240 : : {
241 [ - + ][ - - ]: 30495 : Assert( G.nprop() == m_ncomp*3,
[ - - ][ - - ]
242 : : "Number of components in gradient vector incorrect" );
243 [ - + ][ - - ]: 30495 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
244 : : "vector at recent time step incorrect" );
245 [ - + ][ - - ]: 30495 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
246 : : "Number of unknowns and/or number of components in right-hand "
247 : : "side vector incorrect" );
248 : :
249 : : // compute/assemble gradients in points
250 [ + - ]: 60990 : auto Grad = nodegrad( coord, inpoel, lid, bid, vol, esup, U, G );
251 : :
252 : : // zero right hand side for all components
253 [ + + ][ + - ]: 60990 : for (ncomp_t c=0; c<m_ncomp; ++c) R.fill( c, 0.0 );
254 : :
255 : : // compute domain-edge integral
256 [ + - ]: 30495 : domainint( coord, inpoel, edgeid, psup, dfn, U, Grad, R );
257 : :
258 : : // compute boundary integrals
259 [ + - ]: 30495 : bndint( coord, triinpoel, symbctri, U, R );
260 : 30495 : }
261 : :
262 : : //! Compute boundary pressure integrals (force) (no-op for transport)
263 : 0 : void bndPressureInt(
264 : : const std::array< std::vector< real >, 3 >&,
265 : : const std::vector< std::size_t >&,
266 : : const std::vector< int >&,
267 : : const tk::Fields&,
268 : : const std::array< tk::real, 3 >&,
269 : : std::vector< real >& ) const
270 : 0 : { }
271 : :
272 : :
273 : : //! Compute the minimum time step size (for unsteady time stepping)
274 : : //! \param[in] U Solution vector at recent time step
275 : : //! \param[in] coord Mesh node coordinates
276 : : //! \param[in] inpoel Mesh element connectivity
277 : : //! \param[in] t Physical time
278 : : //! \return Minimum time step size
279 : 75 : real dt( const std::array< std::vector< real >, 3 >& coord,
280 : : const std::vector< std::size_t >& inpoel,
281 : : tk::real t,
282 : : tk::real,
283 : : const tk::Fields& U,
284 : : const std::vector< tk::real >&,
285 : : const std::vector< tk::real >& ) const
286 : : {
287 : : using tag::transport;
288 [ - + ][ - - ]: 75 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
289 : : "vector at recent time step incorrect" );
290 : 75 : const auto& x = coord[0];
291 : 75 : const auto& y = coord[1];
292 : 75 : const auto& z = coord[2];
293 : : // compute the minimum dt across all elements we own
294 : 75 : auto mindt = std::numeric_limits< tk::real >::max();
295 : 75 : auto eps = std::numeric_limits< tk::real >::epsilon();
296 : 75 : auto large = std::numeric_limits< tk::real >::max();
297 [ + + ]: 180507 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
298 : : const std::array< std::size_t, 4 >
299 : 180432 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
300 : : // compute cubic root of element volume as the characteristic length
301 : : const std::array< real, 3 >
302 : 180432 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
303 : 180432 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
304 : 180432 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
305 : 180432 : const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
306 : : // access solution at element nodes at recent time step
307 [ + - ]: 360864 : std::vector< std::array< real, 4 > > u( m_ncomp );
308 [ + + ][ + - ]: 360864 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, N );
309 : : // get velocity for problem
310 : 180432 : const std::array< std::vector<std::array<real,3>>, 4 > vel{{
311 [ + - ]: 180432 : Problem::prescribedVelocity( m_ncomp,
312 : : x[N[0]], y[N[0]], z[N[0]], t ),
313 [ + - ]: 180432 : Problem::prescribedVelocity( m_ncomp,
314 : : x[N[1]], y[N[1]], z[N[1]], t ),
315 [ + - ]: 180432 : Problem::prescribedVelocity( m_ncomp,
316 : : x[N[2]], y[N[2]], z[N[2]], t ),
317 [ + - ]: 180432 : Problem::prescribedVelocity( m_ncomp,
318 : : x[N[3]], y[N[3]], z[N[3]], t ) }};
319 : : // compute the maximum length of the characteristic velocity (advection
320 : : // velocity) across the four element nodes
321 : 180432 : real maxvel = 0.0;
322 [ + + ]: 360864 : for (ncomp_t c=0; c<m_ncomp; ++c)
323 [ + + ]: 902160 : for (std::size_t i=0; i<4; ++i) {
324 : 721728 : auto v = tk::length( vel[i][c] );
325 [ + + ]: 721728 : if (v > maxvel) maxvel = v;
326 : : }
327 : : // compute element dt for the advection
328 [ + - ]: 180432 : auto advection_dt = std::abs(maxvel) > eps ? L / maxvel : large;
329 : : // compute element dt based on diffusion
330 [ - - ]: 180432 : auto diffusion_dt = m_physics.diffusion_dt( m_ncomp, L, u );
331 : : // compute minimum element dt
332 : 180432 : auto elemdt = std::min( advection_dt, diffusion_dt );
333 : : // find minimum dt across all elements
334 [ + + ]: 180432 : if (elemdt < mindt) mindt = elemdt;
335 : : }
336 : 75 : return mindt * g_inputdeck.get< tag::cfl >();
337 : : }
338 : :
339 : : //! Compute a time step size for each mesh node (for steady time stepping)
340 : 0 : void dt( uint64_t,
341 : : const std::vector< tk::real >&,
342 : : const tk::Fields&,
343 : 0 : std::vector< tk::real >& ) const {}
344 : :
345 : : //! \brief Query Dirichlet boundary condition value on a given side set for
346 : : //! all components in this PDE system
347 : : //! \param[in] t Physical time
348 : : //! \param[in] deltat Time step size
349 : : //! \param[in] tp Physical time for each mesh node
350 : : //! \param[in] dtp Time step size for each mesh node
351 : : //! \param[in] ss Pair of side set ID and list of node IDs on the side set
352 : : //! \param[in] coord Mesh node coordinates
353 : : //! \param[in] increment If true, evaluate the solution increment between
354 : : //! t and t+dt for Dirichlet BCs. If false, evlauate the solution instead.
355 : : //! \return Vector of pairs of bool and boundary condition value associated
356 : : //! to mesh node IDs at which Dirichlet boundary conditions are set. Note
357 : : //! that if increment is true, instead of the actual boundary condition
358 : : //! value, we return the increment between t+deltat and t, since,
359 : : //! depending on client code and solver, that may be what the solution
360 : : //! requires.
361 : : std::map< std::size_t, std::vector< std::pair<bool,real> > >
362 : 12835 : dirbc( real t,
363 : : real deltat,
364 : : const std::vector< tk::real >& tp,
365 : : const std::vector< tk::real >& dtp,
366 : : const std::pair< const int, std::vector< std::size_t > >& ss,
367 : : const std::array< std::vector< real >, 3 >& coord,
368 : : bool increment ) const
369 : : {
370 : : using NodeBC = std::vector< std::pair< bool, real > >;
371 : 12835 : std::map< std::size_t, NodeBC > bc;
372 : 12835 : const auto& ubc = g_inputdeck.get< tag::bc >()[0].get< tag::dirichlet >();
373 : 12835 : const auto steady = g_inputdeck.get< tag::steady_state >();
374 [ - + ]: 12835 : if (!ubc.empty()) {
375 [ - - ][ - - ]: 0 : Assert( ubc.size() > 0, "Indexing out of Dirichlet BC eq-vector" );
[ - - ][ - - ]
376 : 0 : const auto& x = coord[0];
377 : 0 : const auto& y = coord[1];
378 : 0 : const auto& z = coord[2];
379 [ - - ]: 0 : for (const auto& b : ubc)
380 [ - - ]: 0 : if (static_cast<int>(b) == ss.first)
381 [ - - ]: 0 : for (auto n : ss.second) {
382 [ - - ][ - - ]: 0 : Assert( x.size() > n, "Indexing out of coordinate array" );
[ - - ][ - - ]
383 [ - - ]: 0 : if (steady) { t = tp[n]; deltat = dtp[n]; }
384 [ - - ][ - - ]: 0 : const auto s = increment ?
[ - - ][ - - ]
385 [ - - ]: 0 : solinc( m_ncomp, m_mat_blk, x[n], y[n], z[n],
386 : : t, deltat, Problem::initialize ) :
387 [ - - ]: 0 : Problem::initialize( m_ncomp, m_mat_blk, x[n], y[n],
388 : : z[n], t+deltat );
389 [ - - ][ - - ]: 0 : auto& nbc = bc[n] = NodeBC( m_ncomp );
390 [ - - ]: 0 : for (ncomp_t c=0; c<m_ncomp; ++c)
391 : 0 : nbc[c] = { true, s[c] };
392 : : }
393 : : }
394 : 12835 : return bc;
395 : : }
396 : :
397 : : //! Set symmetry boundary conditions at nodes
398 : : void
399 : 32298 : symbc(
400 : : tk::Fields&,
401 : : const std::array< std::vector< real >, 3 >&,
402 : : const std::unordered_map< int,
403 : : std::unordered_map< std::size_t,
404 : : std::array< real, 4 > > >&,
405 : 32298 : const std::unordered_set< std::size_t >& ) const {}
406 : :
407 : : //! Set farfield boundary conditions at nodes
408 : 32298 : void farfieldbc(
409 : : tk::Fields&,
410 : : const std::array< std::vector< real >, 3 >&,
411 : : const std::unordered_map< int,
412 : : std::unordered_map< std::size_t,
413 : : std::array< real, 4 > > >&,
414 : 32298 : const std::unordered_set< std::size_t >& ) const {}
415 : :
416 : : //! Apply user defined time dependent BCs (no-op for transport)
417 : : void
418 : 22133 : timedepbc( tk::real,
419 : : tk::Fields&,
420 : : const std::vector< std::unordered_set< std::size_t > >&,
421 : 22133 : const std::vector< tk::Table<5> >& ) const {}
422 : :
423 : : //! Return a map that associates user-specified strings to functions
424 : : //! \return Map that associates user-specified strings to functions that
425 : : //! compute relevant quantities to be output to file
426 : 601 : std::map< std::string, tk::GetVarFn > OutVarFn() const {
427 : 601 : std::map< std::string, tk::GetVarFn > OutFnMap;
428 [ + - ][ + - ]: 601 : OutFnMap["material_indicator"] = transport::matIndicatorOutVar;
[ + - ]
429 : :
430 : 601 : return OutFnMap;
431 : : }
432 : :
433 : : //! Return analytic field names to be output to file
434 : : //! \return Vector of strings labelling analytic fields output in file
435 : 565 : std::vector< std::string > analyticFieldNames() const {
436 : 565 : std::vector< std::string > n;
437 : 565 : auto depvar = g_inputdeck.get< tag::depvar >()[0];
438 [ + + ]: 1130 : for (ncomp_t c=0; c<m_ncomp; ++c)
439 [ + - ][ + - ]: 565 : n.push_back( depvar + std::to_string(c) + "_analytic" );
[ + - ][ + - ]
440 : 565 : return n;
441 : : }
442 : :
443 : : //! Return surface field names to be output to file
444 : : //! \return Vector of strings labelling surface fields output in file
445 : : //! \details This functions should be written in conjunction with
446 : : //! surfOutput(), which provides the vector of surface fields to be output
447 : 601 : std::vector< std::string > surfNames() const { return {}; }
448 : :
449 : : //! Return nodal surface field output going to file
450 : : std::vector< std::vector< real > >
451 : 601 : surfOutput( const std::map< int, std::vector< std::size_t > >&,
452 : 601 : const tk::Fields& ) const { return {}; }
453 : :
454 : : //! Return elemental surface field output (on triangle faces) going to file
455 : : std::vector< std::vector< real > >
456 : 601 : elemSurfOutput( const std::map< int, std::vector< std::size_t > >&,
457 : : const std::vector< std::size_t >&,
458 : 601 : const tk::Fields& ) const { return {}; }
459 : :
460 : : //! Return time history field names to be output to file
461 : : //! \return Vector of strings labelling time history fields output in file
462 : 0 : std::vector< std::string > histNames() const { return {}; }
463 : :
464 : : //! Return time history field output evaluated at time history points
465 : : std::vector< std::vector< real > >
466 : 0 : histOutput( const std::vector< HistData >&,
467 : : const std::vector< std::size_t >&,
468 : 0 : const tk::Fields& ) const { return {}; }
469 : :
470 : : //! Return names of integral variables to be output to diagnostics file
471 : : //! \return Vector of strings labelling integral variables output
472 : 52 : std::vector< std::string > names() const {
473 : 52 : std::vector< std::string > n;
474 : : const auto& depvar =
475 [ + - ]: 52 : g_inputdeck.get< tag::depvar >().at(0);
476 : : // construct the name of the numerical solution for all components
477 [ + + ]: 104 : for (ncomp_t c=0; c<m_ncomp; ++c)
478 [ + - ][ + - ]: 52 : n.push_back( depvar + std::to_string(c) );
[ + - ]
479 : 52 : return n;
480 : : }
481 : :
482 : : private:
483 : : const Physics m_physics; //!< Physics policy
484 : : const Problem m_problem; //!< Problem policy
485 : : const ncomp_t m_ncomp; //!< Number of components in this PDE
486 : : //! EOS material block
487 : : const std::vector< EOS > m_mat_blk;
488 : :
489 : : //! \brief Compute/assemble nodal gradients of primitive variables for
490 : : //! ALECG in all points
491 : : //! \param[in] coord Mesh node coordinates
492 : : //! \param[in] inpoel Mesh element connectivity
493 : : //! \param[in] lid Global->local node ids
494 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
495 : : //! global node ids (key)
496 : : //! \param[in] vol Nodal volumes
497 : : //! \param[in] esup Elements surrounding points
498 : : //! \param[in] U Solution vector at recent time step
499 : : //! \param[in] G Nodal gradients of primitive variables in chare-boundary nodes
500 : : //! \return Gradients of primitive variables in all mesh points
501 : : tk::Fields
502 : 30495 : nodegrad( const std::array< std::vector< real >, 3 >& coord,
503 : : const std::vector< std::size_t >& inpoel,
504 : : const std::unordered_map< std::size_t, std::size_t >& lid,
505 : : const std::unordered_map< std::size_t, std::size_t >& bid,
506 : : const std::vector< real >& vol,
507 : : const std::pair< std::vector< std::size_t >,
508 : : std::vector< std::size_t > >& esup,
509 : : const tk::Fields& U,
510 : : const tk::Fields& G ) const
511 : : {
512 : : // allocate storage for nodal gradients of primitive variables
513 : 30495 : tk::Fields Grad( U.nunk(), m_ncomp*3 );
514 : 30495 : Grad.fill( 0.0 );
515 : :
516 : : // access node cooordinates
517 : 30495 : const auto& x = coord[0];
518 : 30495 : const auto& y = coord[1];
519 : 30495 : const auto& z = coord[2];
520 : :
521 : : // compute gradients of primitive variables in points
522 : 30495 : auto npoin = U.nunk();
523 : : #pragma omp simd
524 [ + + ]: 2094207 : for (std::size_t p=0; p<npoin; ++p)
525 [ + + ]: 16125696 : for (auto e : tk::Around(esup,p)) {
526 : : // access node IDs
527 : 14061984 : std::size_t N[4] =
528 : 14061984 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
529 : : // compute element Jacobi determinant, J = 6V
530 : 14061984 : real bax = x[N[1]]-x[N[0]];
531 : 14061984 : real bay = y[N[1]]-y[N[0]];
532 : 14061984 : real baz = z[N[1]]-z[N[0]];
533 : 14061984 : real cax = x[N[2]]-x[N[0]];
534 : 14061984 : real cay = y[N[2]]-y[N[0]];
535 : 14061984 : real caz = z[N[2]]-z[N[0]];
536 : 14061984 : real dax = x[N[3]]-x[N[0]];
537 : 14061984 : real day = y[N[3]]-y[N[0]];
538 : 14061984 : real daz = z[N[3]]-z[N[0]];
539 : 14061984 : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
540 : 14061984 : auto J24 = J/24.0;
541 : : // shape function derivatives, nnode*ndim [4][3]
542 : : real g[4][3];
543 : 14061984 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
544 : : g[1][0], g[1][1], g[1][2] );
545 : 14061984 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
546 : : g[2][0], g[2][1], g[2][2] );
547 : 14061984 : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
548 : : g[3][0], g[3][1], g[3][2] );
549 [ + + ]: 56247936 : for (std::size_t i=0; i<3; ++i)
550 : 42185952 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
551 : : // scatter-add gradient contributions to boundary nodes
552 [ + + ]: 28123968 : for (std::size_t c=0; c<m_ncomp; ++c)
553 [ + + ]: 70309920 : for (std::size_t b=0; b<4; ++b)
554 [ + + ]: 224991744 : for (std::size_t i=0; i<3; ++i)
555 [ + - ][ + - ]: 168743808 : Grad(p,c*3+i) += J24 * g[b][i] * U(N[b],c);
556 : : }
557 : :
558 : : // put in nodal gradients of chare-boundary points
559 [ + + ]: 1435941 : for (const auto& [g,b] : bid) {
560 [ + - ]: 1405446 : auto i = tk::cref_find( lid, g );
561 [ + + ]: 5621784 : for (ncomp_t c=0; c<Grad.nprop(); ++c)
562 [ + - ][ + - ]: 4216338 : Grad(i,c) = G(b,c);
563 : : }
564 : :
565 : : // divide weak result in gradients by nodal volume
566 [ + + ]: 2094207 : for (std::size_t p=0; p<npoin; ++p)
567 [ + + ]: 8254848 : for (std::size_t c=0; c<m_ncomp*3; ++c)
568 [ + - ]: 6191136 : Grad(p,c) /= vol[p];
569 : :
570 : 30495 : return Grad;
571 : : }
572 : :
573 : : //! \brief Compute MUSCL reconstruction in edge-end points using a MUSCL
574 : : //! procedure with van Leer limiting
575 : : //! \param[in] p Left node id of edge-end
576 : : //! \param[in] q Right node id of edge-end
577 : : //! \param[in] coord Array of nodal coordinates
578 : : //! \param[in] G Gradient of all unknowns in mesh points
579 : : //! \param[in,out] uL Primitive variables at left edge-end point
580 : : //! \param[in,out] uR Primitive variables at right edge-end point
581 : : void
582 : 15738216 : muscl( std::size_t p,
583 : : std::size_t q,
584 : : const tk::UnsMesh::Coords& coord,
585 : : const tk::Fields& G,
586 : : std::vector< real >& uL,
587 : : std::vector< real >& uR ) const
588 : : {
589 [ + - ][ + - ]: 15738216 : Assert( uL.size() == m_ncomp && uR.size() == m_ncomp, "Size mismatch" );
[ - - ][ - - ]
[ - - ]
590 [ - + ][ - - ]: 15738216 : Assert( G.nprop()/3 == m_ncomp, "Size mismatch" );
[ - - ][ - - ]
591 : :
592 : 15738216 : const auto& x = coord[0];
593 : 15738216 : const auto& y = coord[1];
594 : 15738216 : const auto& z = coord[2];
595 : :
596 : : // edge vector
597 : 15738216 : std::array< real, 3 > vw{ x[q]-x[p], y[q]-y[p], z[q]-z[p] };
598 : :
599 : : std::vector< real >
600 [ + - ][ + - ]: 31476432 : delta1( m_ncomp, 0.0 ), delta2( m_ncomp, 0.0 ), delta3( m_ncomp, 0.0 );
[ + - ]
601 : :
602 : : // MUSCL reconstruction of edge-end-point primitive variables
603 [ + + ]: 31476432 : for (std::size_t c=0; c<m_ncomp; ++c) {
604 : : // gradients
605 : : std::array< real, 3 >
606 [ + - ][ + - ]: 15738216 : g1{ G(p,c*3+0), G(p,c*3+1), G(p,c*3+2) },
[ + - ]
607 [ + - ][ + - ]: 15738216 : g2{ G(q,c*3+0), G(q,c*3+1), G(q,c*3+2) };
[ + - ]
608 : :
609 : 15738216 : delta2[c] = uR[c] - uL[c];
610 : 15738216 : delta1[c] = 2.0 * tk::dot(g1,vw) - delta2[c];
611 : 15738216 : delta3[c] = 2.0 * tk::dot(g2,vw) - delta2[c];
612 : :
613 : : // form limiters
614 : 15738216 : auto rL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
615 : 15738216 : auto rR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
616 : 15738216 : auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
617 : 15738216 : auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
618 : :
619 : 15738216 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
620 : 15738216 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
621 : 15738216 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
622 : 15738216 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
623 : :
624 : : // update unknowns with reconstructed unknowns
625 : 15738216 : uL[c] += 0.25*(delta1[c]*muscl_m1*phiL + delta2[c]*muscl_p1*phi_L_inv);
626 : 15738216 : uR[c] -= 0.25*(delta3[c]*muscl_m1*phiR + delta2[c]*muscl_p1*phi_R_inv);
627 : : }
628 : 15738216 : }
629 : :
630 : : //! Compute domain-edge integral for ALECG
631 : : //! \param[in] coord Mesh node coordinates
632 : : //! \param[in] inpoel Mesh element connectivity
633 : : //! \param[in] edgeid Local node id pair -> edge id map
634 : : //! \param[in] psup Points surrounding points
635 : : //! \param[in] dfn Dual-face normals
636 : : //! \param[in] U Solution vector at recent time step
637 : : //! \param[in] G Nodal gradients
638 : : //! \param[in,out] R Right-hand side vector computed
639 : 30495 : void domainint( const std::array< std::vector< real >, 3 >& coord,
640 : : const std::vector< std::size_t >& inpoel,
641 : : const std::vector< std::size_t >& edgeid,
642 : : const std::pair< std::vector< std::size_t >,
643 : : std::vector< std::size_t > >& psup,
644 : : const std::vector< real >& dfn,
645 : : const tk::Fields& U,
646 : : const tk::Fields& G,
647 : : tk::Fields& R ) const
648 : : {
649 : : // access node cooordinates
650 : 30495 : const auto& x = coord[0];
651 : 30495 : const auto& y = coord[1];
652 : 30495 : const auto& z = coord[2];
653 : :
654 : : // compute derived data structures
655 [ + - ][ + - ]: 60990 : auto esued = tk::genEsued( inpoel, 4, tk::genEsup( inpoel, 4 ) );
656 : :
657 : : // access pointer to right hand side at component
658 [ + - ]: 60990 : std::vector< const real* > r( m_ncomp );
659 [ + + ][ + - ]: 60990 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c );
660 : :
661 : : // domain-edge integral
662 [ + + ]: 2094207 : for (std::size_t p=0,k=0; p<U.nunk(); ++p) {
663 [ + + ]: 17801928 : for (auto q : tk::Around(psup,p)) {
664 : : // access dual-face normals for edge p-q
665 : 15738216 : auto ed = edgeid[k++];
666 : 15738216 : std::array< tk::real, 3 > n{ dfn[ed*6+0], dfn[ed*6+1], dfn[ed*6+2] };
667 : :
668 [ + - ]: 31476432 : std::vector< tk::real > uL( m_ncomp, 0.0 );
669 [ + - ]: 31476432 : std::vector< tk::real > uR( m_ncomp, 0.0 );
670 [ + + ]: 31476432 : for (std::size_t c=0; c<m_ncomp; ++c) {
671 [ + - ]: 15738216 : uL[c] = U(p,c);
672 [ + - ]: 15738216 : uR[c] = U(q,c);
673 : : }
674 : : // compute MUSCL reconstruction in edge-end points
675 [ + - ]: 15738216 : muscl( p, q, coord, G, uL, uR );
676 : :
677 : : // evaluate prescribed velocity
678 : 15738216 : auto v =
679 [ + - ]: 15738216 : Problem::prescribedVelocity( m_ncomp, x[p], y[p], z[p], 0.0 );
680 : : // sum donain-edge contributions
681 [ + - ][ + + ]: 57924168 : for (auto e : tk::cref_find(esued,{p,q})) {
682 : : const std::array< std::size_t, 4 >
683 : 42185952 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
684 : : // compute element Jacobi determinant
685 : : const std::array< tk::real, 3 >
686 : 42185952 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
687 : 42185952 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
688 : 42185952 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
689 : 42185952 : const auto J = tk::triple( ba, ca, da ); // J = 6V
690 : : // shape function derivatives, nnode*ndim [4][3]
691 : : std::array< std::array< tk::real, 3 >, 4 > grad;
692 : 42185952 : grad[1] = tk::crossdiv( ca, da, J );
693 : 42185952 : grad[2] = tk::crossdiv( da, ba, J );
694 : 42185952 : grad[3] = tk::crossdiv( ba, ca, J );
695 [ + + ]: 168743808 : for (std::size_t i=0; i<3; ++i)
696 : 126557856 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
697 : 42185952 : auto J48 = J/48.0;
698 [ + + ]: 295301664 : for (const auto& [a,b] : tk::lpoed) {
699 [ + - ]: 253115712 : auto s = tk::orient( {N[a],N[b]}, {p,q} );
700 [ + + ]: 1012462848 : for (std::size_t j=0; j<3; ++j) {
701 [ + + ]: 1518694272 : for (std::size_t c=0; c<m_ncomp; ++c) {
702 [ + - ]: 1518694272 : R.var(r[c],p) -= J48 * s * (grad[a][j] - grad[b][j])
703 : 759347136 : * v[c][j]*(uL[c] + uR[c])
704 : 759347136 : - J48 * std::abs(s * (grad[a][j] - grad[b][j]))
705 : 759347136 : * std::abs(tk::dot(v[c],n)) * (uR[c] - uL[c]);
706 : : }
707 : : }
708 : : }
709 : : }
710 : : }
711 : : }
712 : 30495 : }
713 : :
714 : : //! Compute boundary integrals for ALECG
715 : : //! \param[in] coord Mesh node coordinates
716 : : //! \param[in] triinpoel Boundary triangle face connecitivity with local ids
717 : : //! \param[in] symbctri Vector with 1 at symmetry BC boundary triangles
718 : : //! \param[in] U Solution vector at recent time step
719 : : //! \param[in,out] R Right-hand side vector computed
720 : 30495 : void bndint( const std::array< std::vector< real >, 3 >& coord,
721 : : const std::vector< std::size_t >& triinpoel,
722 : : const std::vector< int >& symbctri,
723 : : const tk::Fields& U,
724 : : tk::Fields& R ) const
725 : : {
726 : : // access node coordinates
727 : 30495 : const auto& x = coord[0];
728 : 30495 : const auto& y = coord[1];
729 : 30495 : const auto& z = coord[2];
730 : :
731 : : // boundary integrals: compute fluxes in edges
732 [ + - ]: 60990 : std::vector< real > bflux( triinpoel.size() * m_ncomp * 2 );
733 : :
734 [ + + ]: 1277295 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
735 : : // access node IDs
736 : : std::array< std::size_t, 3 >
737 : 1246800 : N{ triinpoel[e*3+0], triinpoel[e*3+1], triinpoel[e*3+2] };
738 : : // apply symmetry BCs
739 [ + - ]: 1246800 : if (symbctri[e]) continue;
740 : : // node coordinates
741 : 0 : std::array< tk::real, 3 > xp{ x[N[0]], x[N[1]], x[N[2]] },
742 : 0 : yp{ y[N[0]], y[N[1]], y[N[2]] },
743 : 0 : zp{ z[N[0]], z[N[1]], z[N[2]] };
744 : : // access solution at element nodes
745 [ - - ]: 0 : std::vector< std::array< real, 3 > > u( m_ncomp );
746 [ - - ][ - - ]: 0 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, N );
747 : : // evaluate prescribed velocity
748 : 0 : auto v =
749 [ - - ]: 0 : Problem::prescribedVelocity( m_ncomp, xp[0], yp[0], zp[0], 0.0 );
750 : : // compute face area
751 : 0 : auto A6 = tk::area( x[N[0]], x[N[1]], x[N[2]],
752 : : y[N[0]], y[N[1]], y[N[2]],
753 : : z[N[0]], z[N[1]], z[N[2]] ) / 6.0;
754 : 0 : auto A24 = A6/4.0;
755 : : // compute face normal
756 [ - - ]: 0 : auto n = tk::normal( xp, yp, zp );
757 : : // store flux in boundary elements
758 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
759 : 0 : auto eb = (e*m_ncomp+c)*6;
760 : 0 : auto vdotn = tk::dot( v[c], n );
761 : 0 : auto Bab = A24 * vdotn * (u[c][0] + u[c][1]);
762 : 0 : bflux[eb+0] = Bab + A6 * vdotn * u[c][0];
763 : 0 : bflux[eb+1] = Bab;
764 : 0 : Bab = A24 * vdotn * (u[c][1] + u[c][2]);
765 : 0 : bflux[eb+2] = Bab + A6 * vdotn * u[c][1];
766 : 0 : bflux[eb+3] = Bab;
767 : 0 : Bab = A24 * vdotn * (u[c][2] + u[c][0]);
768 : 0 : bflux[eb+4] = Bab + A6 * vdotn * u[c][2];
769 : 0 : bflux[eb+5] = Bab;
770 : : }
771 : : }
772 : :
773 : : // access pointer to right hand side at component
774 [ + - ]: 60990 : std::vector< const real* > r( m_ncomp );
775 [ + + ][ + - ]: 60990 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c );
776 : :
777 : : // boundary integrals: sum flux contributions to points
778 [ + + ]: 1277295 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
779 : 1246800 : std::size_t N[3] =
780 : 1246800 : { triinpoel[e*3+0], triinpoel[e*3+1], triinpoel[e*3+2] };
781 [ + + ]: 2493600 : for (std::size_t c=0; c<m_ncomp; ++c) {
782 : 1246800 : auto eb = (e*m_ncomp+c)*6;
783 [ + - ]: 1246800 : R.var(r[c],N[0]) -= bflux[eb+0] + bflux[eb+5];
784 [ + - ]: 1246800 : R.var(r[c],N[1]) -= bflux[eb+1] + bflux[eb+2];
785 [ + - ]: 1246800 : R.var(r[c],N[2]) -= bflux[eb+3] + bflux[eb+4];
786 : : }
787 : : }
788 : :
789 : 30495 : tk::destroy(bflux);
790 : 30495 : }
791 : : };
792 : :
793 : : } // cg::
794 : : } // inciter::
795 : :
796 : : #endif // Transport_h
|
