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