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