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 = kw::ncomp::info::expect::type;
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 : : //! \param[in] c Equation system index (among multiple systems configured)
61 : 232 : explicit Transport( ncomp_t c ) :
62 : : m_physics( Physics() ),
63 : : m_problem( Problem() ),
64 : : m_system( c ),
65 : : m_ncomp(
66 : 232 : g_inputdeck.get< tag::component >().get< tag::transport >().at(c) ),
67 : : m_offset(
68 : 464 : g_inputdeck.get< tag::component >().offset< tag::transport >(c) )
69 : : {
70 : 232 : m_problem.errchk( m_system, m_ncomp );
71 : 232 : }
72 : :
73 : : //! Determine nodes that lie inside the user-defined IC box
74 : : void
75 : 857 : IcBoxNodes( const tk::UnsMesh::Coords&,
76 : 857 : std::vector< std::unordered_set< std::size_t > >& ) const {}
77 : :
78 : : //! Initalize the transport equations using problem policy
79 : : //! \param[in] coord Mesh node coordinates
80 : : //! \param[in,out] unk Array of unknowns
81 : : //! \param[in] t Physical time
82 : : void
83 : 1642 : initialize( const std::array< std::vector< real >, 3 >& coord,
84 : : tk::Fields& unk,
85 : : real t,
86 : : real,
87 : : const std::vector< std::unordered_set< std::size_t > >& ) const
88 : : {
89 [ - + ][ - - ]: 1642 : Assert( coord[0].size() == unk.nunk(), "Size mismatch" );
[ - - ][ - - ]
90 : 1642 : const auto& x = coord[0];
91 : 1642 : const auto& y = coord[1];
92 : 1642 : const auto& z = coord[2];
93 [ + + ]: 866029 : for (ncomp_t i=0; i<x.size(); ++i) {
94 [ + - ]: 1728774 : auto s = Problem::initialize( m_system, m_ncomp, x[i], y[i], z[i], t );
95 [ + + ]: 1740486 : for (ncomp_t c=0; c<m_ncomp; ++c)
96 [ + - ]: 876099 : unk( i, c, m_offset ) = s[c];
97 : : }
98 : 1642 : }
99 : :
100 : : //! Query a velocity
101 : : //! \note Since this function does not touch its output argument, that
102 : : //! means this system does not define a "velocity".
103 : 22259 : void velocity( const tk::Fields&, tk::UnsMesh::Coords& ) const {}
104 : :
105 : : //! Query the sound speed
106 : : //! \note Since this function does not touch its output argument, that
107 : : //! means this system does not define a "sound speed".
108 : 22259 : void soundspeed( const tk::Fields&, std::vector< tk::real >& ) const {}
109 : :
110 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
111 : : //! \param[in] xi X-coordinate
112 : : //! \param[in] yi Y-coordinate
113 : : //! \param[in] zi Z-coordinate
114 : : //! \param[in] t Physical time
115 : : //! \return Vector of analytic solution at given location and time
116 : : std::vector< real >
117 : 386224 : analyticSolution( real xi, real yi, real zi, real t ) const
118 : 386224 : { return Problem::analyticSolution( m_system, m_ncomp, xi, yi, zi, t ); }
119 : :
120 : : //! Return analytic solution for conserved variables
121 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
122 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
123 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
124 : : //! \param[in] t Physical time at which to evaluate the analytic solution
125 : : //! \return Vector of analytic solution at given location and time
126 : : std::vector< tk::real >
127 : 4113811 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
128 : 4113811 : { return Problem::initialize( m_system, m_ncomp, xi, yi, zi, t ); }
129 : :
130 : : //! Compute nodal gradients of primitive variables for ALECG
131 : : //! \param[in] coord Mesh node coordinates
132 : : //! \param[in] inpoel Mesh element connectivity
133 : : //! \param[in] bndel List of elements contributing to chare-boundary nodes
134 : : //! \param[in] gid Local->global node id map
135 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
136 : : //! global node ids (key)
137 : : //! \param[in] U Solution vector at recent time step
138 : : //! \param[in,out] G Nodal gradients of primitive variables
139 : 21870 : void chBndGrad( const std::array< std::vector< real >, 3 >& coord,
140 : : const std::vector< std::size_t >& inpoel,
141 : : const std::vector< std::size_t >& bndel,
142 : : const std::vector< std::size_t >& gid,
143 : : const std::unordered_map< std::size_t, std::size_t >& bid,
144 : : const tk::Fields& U,
145 : : tk::Fields& G ) const
146 : : {
147 [ - + ][ - - ]: 21870 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
148 : : "vector at recent time step incorrect" );
149 : :
150 : : // compute gradients of primitive variables in points
151 : 21870 : G.fill( 0.0 );
152 : :
153 : : // access node cooordinates
154 : 21870 : const auto& x = coord[0];
155 : 21870 : const auto& y = coord[1];
156 : 21870 : const auto& z = coord[2];
157 : :
158 [ + + ]: 7804290 : for (auto e : bndel) { // elements contributing to chare boundary nodes
159 : : // access node IDs
160 : 7782420 : std::size_t N[4] =
161 : 7782420 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
162 : : // compute element Jacobi determinant, J = 6V
163 : 7782420 : real bax = x[N[1]]-x[N[0]];
164 : 7782420 : real bay = y[N[1]]-y[N[0]];
165 : 7782420 : real baz = z[N[1]]-z[N[0]];
166 : 7782420 : real cax = x[N[2]]-x[N[0]];
167 : 7782420 : real cay = y[N[2]]-y[N[0]];
168 : 7782420 : real caz = z[N[2]]-z[N[0]];
169 : 7782420 : real dax = x[N[3]]-x[N[0]];
170 : 7782420 : real day = y[N[3]]-y[N[0]];
171 : 7782420 : real daz = z[N[3]]-z[N[0]];
172 : 7782420 : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
173 : 7782420 : auto J24 = J/24.0;
174 : : // shape function derivatives, nnode*ndim [4][3]
175 : : real g[4][3];
176 : 7782420 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
177 : : g[1][0], g[1][1], g[1][2] );
178 : 7782420 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
179 : : g[2][0], g[2][1], g[2][2] );
180 : 7782420 : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
181 : : g[3][0], g[3][1], g[3][2] );
182 [ + + ]: 31129680 : for (std::size_t i=0; i<3; ++i)
183 : 23347260 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
184 : : // scatter-add gradient contributions to boundary nodes
185 [ + + ]: 38912100 : for (std::size_t a=0; a<4; ++a) {
186 [ + - ]: 31129680 : auto i = bid.find( gid[N[a]] );
187 [ + + ]: 31129680 : if (i != end(bid))
188 [ + + ]: 42002820 : for (std::size_t c=0; c<m_ncomp; ++c)
189 [ + + ]: 105007050 : for (std::size_t b=0; b<4; ++b)
190 [ + + ]: 336022560 : for (std::size_t j=0; j<3; ++j)
191 [ + - ][ + - ]: 252016920 : G(i->second,c*3+j,0) += J24 * g[b][j] * U(N[b],c,m_offset);
192 : : }
193 : : }
194 : 21870 : }
195 : :
196 : : //! Compute right hand side for ALECG
197 : : //! \param[in] coord Mesh node coordinates
198 : : //! \param[in] inpoel Mesh element connectivity
199 : : //! \param[in] triinpoel Boundary triangle face connecitivity
200 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
201 : : //! global node ids (key)
202 : : //! \param[in] lid Global->local node ids
203 : : //! \param[in] dfn Dual-face normals
204 : : //! \param[in] psup Points surrounding points
205 : : //! \param[in] esup Elements surrounding points
206 : : //! \param[in] symbctri Vector with 1 at symmetry BC nodes
207 : : //! \param[in] vol Nodal volumes
208 : : //! \param[in] edgeid Local node id pair -> edge id map
209 : : //! \param[in] G Nodal gradients in chare-boundary nodes
210 : : //! \param[in] U Solution vector at recent time step
211 : : //! \param[in] W Mesh velocity
212 : : //! \param[in,out] R Right-hand side vector computed
213 : 21870 : void rhs(
214 : : real,
215 : : const std::array< std::vector< real >, 3 >& coord,
216 : : const std::vector< std::size_t >& inpoel,
217 : : const std::vector< std::size_t >& triinpoel,
218 : : const std::vector< std::size_t >&,
219 : : const std::unordered_map< std::size_t, std::size_t >& bid,
220 : : const std::unordered_map< std::size_t, std::size_t >& lid,
221 : : const std::vector< real >& dfn,
222 : : const std::pair< std::vector< std::size_t >,
223 : : std::vector< std::size_t > >& psup,
224 : : const std::pair< std::vector< std::size_t >,
225 : : std::vector< std::size_t > >& esup,
226 : : const std::vector< int >& symbctri,
227 : : const std::unordered_set< std::size_t >&,
228 : : const std::vector< real >& vol,
229 : : const std::vector< std::size_t >&,
230 : : const std::vector< std::size_t >& edgeid,
231 : : const std::vector< std::unordered_set< std::size_t > >&,
232 : : const tk::Fields& G,
233 : : const tk::Fields& U,
234 : : [[maybe_unused]] const tk::Fields& W,
235 : : const std::vector< tk::real >&,
236 : : real,
237 : : tk::Fields& R ) const
238 : : {
239 [ - + ][ - - ]: 21870 : Assert( G.nprop() == m_ncomp*3,
[ - - ][ - - ]
240 : : "Number of components in gradient vector incorrect" );
241 [ - + ][ - - ]: 21870 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
242 : : "vector at recent time step incorrect" );
243 [ - + ][ - - ]: 21870 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
244 : : "Number of unknowns and/or number of components in right-hand "
245 : : "side vector incorrect" );
246 : :
247 : : // compute/assemble gradients in points
248 [ + - ]: 43740 : auto Grad = nodegrad( coord, inpoel, lid, bid, vol, esup, U, G );
249 : :
250 : : // zero right hand side for all components
251 [ + + ][ + - ]: 43740 : for (ncomp_t c=0; c<m_ncomp; ++c) R.fill( c, m_offset, 0.0 );
252 : :
253 : : // compute domain-edge integral
254 [ + - ]: 21870 : domainint( coord, inpoel, edgeid, psup, dfn, U, Grad, R );
255 : :
256 : : // compute boundary integrals
257 [ + - ]: 21870 : bndint( coord, triinpoel, symbctri, U, R );
258 : 21870 : }
259 : :
260 : : //! Compute right hand side for DiagCG (CG+FCT)
261 : : //! \param[in] t Physical time
262 : : //! \param[in] deltat Size of time step
263 : : //! \param[in] coord Mesh node coordinates
264 : : //! \param[in] inpoel Mesh element connectivity
265 : : //! \param[in] U Solution vector at recent time step
266 : : //! \param[in,out] Ue Element-centered solution vector at intermediate step
267 : : //! (used here internally as a scratch array)
268 : : //! \param[in,out] R Right-hand side vector computed
269 : 6279 : void rhs( real t,
270 : : real deltat,
271 : : const std::array< std::vector< real >, 3 >& coord,
272 : : const std::vector< std::size_t >& inpoel,
273 : : const tk::Fields& U,
274 : : tk::Fields& Ue,
275 : : tk::Fields& R ) const
276 : : {
277 : : using tag::transport;
278 [ - + ][ - - ]: 6279 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
279 : : "vector at recent time step incorrect" );
280 [ - + ][ - - ]: 6279 : Assert( R.nunk() == coord[0].size(),
[ - - ][ - - ]
281 : : "Number of unknowns in right-hand side vector incorrect" );
282 : :
283 : 6279 : const auto& x = coord[0];
284 : 6279 : const auto& y = coord[1];
285 : 6279 : const auto& z = coord[2];
286 : :
287 : : // 1st stage: update element values from node values (gather-add)
288 [ + + ]: 8133961 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
289 : : // access node IDs
290 : : const std::array< std::size_t, 4 >
291 : 8127682 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
292 : : // compute element Jacobi determinant
293 : : const std::array< real, 3 >
294 : 8127682 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
295 : 8127682 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
296 : 8127682 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
297 : 8127682 : const auto J = tk::triple( ba, ca, da ); // J = 6V
298 [ - + ][ - - ]: 8127682 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
299 : :
300 : : // shape function derivatives, nnode*ndim [4][3]
301 : : std::array< std::array< real, 3 >, 4 > grad;
302 : 8127682 : grad[1] = tk::crossdiv( ca, da, J );
303 : 8127682 : grad[2] = tk::crossdiv( da, ba, J );
304 : 8127682 : grad[3] = tk::crossdiv( ba, ca, J );
305 [ + + ]: 32510728 : for (std::size_t i=0; i<3; ++i)
306 : 24383046 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
307 : :
308 : : // access solution at element nodes
309 [ + - ]: 16255364 : std::vector< std::array< real, 4 > > u( m_ncomp );
310 [ + + ][ + - ]: 16380424 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, m_offset, N );
311 : : // access solution at elements
312 [ + - ]: 16255364 : std::vector< const real* > ue( m_ncomp );
313 [ + + ][ + - ]: 16380424 : for (ncomp_t c=0; c<m_ncomp; ++c) ue[c] = Ue.cptr( c, m_offset );
314 : :
315 : : // get prescribed velocity
316 : 8127682 : const std::array< std::vector<std::array<real,3>>, 4 > vel{{
317 [ + - ]: 8127682 : Problem::prescribedVelocity( m_system, m_ncomp,
318 : : x[N[0]], y[N[0]], z[N[0]], t ),
319 [ + - ]: 8127682 : Problem::prescribedVelocity( m_system, m_ncomp,
320 : : x[N[1]], y[N[1]], z[N[1]], t ),
321 [ + - ]: 8127682 : Problem::prescribedVelocity( m_system, m_ncomp,
322 : : x[N[2]], y[N[2]], z[N[2]], t ),
323 [ + - ]: 8127682 : Problem::prescribedVelocity( m_system, m_ncomp,
324 : : x[N[3]], y[N[3]], z[N[3]], t ) }};
325 : :
326 : : // sum flux (advection) contributions to element
327 : 8127682 : auto d = deltat/2.0;
328 [ + + ]: 16380424 : for (std::size_t c=0; c<m_ncomp; ++c)
329 [ + + ]: 33010968 : for (std::size_t j=0; j<3; ++j)
330 [ + + ]: 123791130 : for (std::size_t a=0; a<4; ++a)
331 [ + - ]: 99032904 : Ue.var(ue[c],e) -= d * grad[a][j] * vel[a][c][j]*u[c][a];
332 : : }
333 : :
334 : :
335 : : // 2nd stage: form rhs from element values (scatter-add)
336 [ + + ]: 8133961 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
337 : : // access node IDs
338 : : const std::array< std::size_t, 4 >
339 : 8127682 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
340 : : // compute element Jacobi determinant
341 : : const std::array< real, 3 >
342 : 8127682 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
343 : 8127682 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
344 : 8127682 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
345 : 8127682 : const auto J = tk::triple( ba, ca, da ); // J = 6V
346 [ - + ][ - - ]: 8127682 : Assert( J > 0, "Element Jacobian non-positive" );
[ - - ][ - - ]
347 : :
348 : : // shape function derivatives, nnode*ndim [4][3]
349 : : std::array< std::array< real, 3 >, 4 > grad;
350 : 8127682 : grad[1] = tk::crossdiv( ca, da, J );
351 : 8127682 : grad[2] = tk::crossdiv( da, ba, J );
352 : 8127682 : grad[3] = tk::crossdiv( ba, ca, J );
353 [ + + ]: 32510728 : for (std::size_t i=0; i<3; ++i)
354 : 24383046 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
355 : :
356 : : // access solution at elements
357 [ + - ]: 16255364 : std::vector< real > ue( m_ncomp );
358 [ + + ][ + - ]: 16380424 : for (ncomp_t c=0; c<m_ncomp; ++c) ue[c] = Ue( e, c, m_offset );
359 : : // access pointer to right hand side at component and offset
360 [ + - ]: 16255364 : std::vector< const real* > r( m_ncomp );
361 [ + + ][ + - ]: 16380424 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c, m_offset );
362 : : // access solution at nodes of element
363 [ + - ]: 16255364 : std::vector< std::array< real, 4 > > u( m_ncomp );
364 [ + + ][ + - ]: 16380424 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, m_offset, N );
365 : :
366 : : // get prescribed velocity
367 : 8127682 : auto xc = (x[N[0]] + x[N[1]] + x[N[2]] + x[N[3]]) / 4.0;
368 : 8127682 : auto yc = (y[N[0]] + y[N[1]] + y[N[2]] + y[N[3]]) / 4.0;
369 : 8127682 : auto zc = (z[N[0]] + z[N[1]] + z[N[2]] + z[N[3]]) / 4.0;
370 : 8127682 : const auto vel =
371 [ + - ]: 8127682 : Problem::prescribedVelocity( m_system, m_ncomp, xc, yc, zc, t );
372 : :
373 : : // scatter-add flux contributions to rhs at nodes
374 : 8127682 : real d = deltat * J/6.0;
375 [ + + ]: 16380424 : for (std::size_t c=0; c<m_ncomp; ++c)
376 [ + + ]: 33010968 : for (std::size_t j=0; j<3; ++j)
377 [ + + ]: 123791130 : for (std::size_t a=0; a<4; ++a)
378 [ + - ]: 99032904 : R.var(r[c],N[a]) += d * grad[a][j] * vel[c][j]*ue[c];
379 : :
380 : : // add (optional) diffusion contribution to right hand side
381 [ + - ]: 8127682 : m_physics.diffusionRhs(m_system, m_ncomp, deltat, J, grad, N, u, r, R);
382 : : }
383 : 6279 : }
384 : :
385 : : //! Compute the minimum time step size (for unsteady time stepping)
386 : : //! \param[in] U Solution vector at recent time step
387 : : //! \param[in] coord Mesh node coordinates
388 : : //! \param[in] inpoel Mesh element connectivity
389 : : //! \param[in] t Physical time
390 : : //! \return Minimum time step size
391 : 3159 : real dt( const std::array< std::vector< real >, 3 >& coord,
392 : : const std::vector< std::size_t >& inpoel,
393 : : tk::real t,
394 : : tk::real,
395 : : const tk::Fields& U,
396 : : const std::vector< tk::real >&,
397 : : const std::vector< tk::real >& ) const
398 : : {
399 : : using tag::transport;
400 [ - + ][ - - ]: 3159 : Assert( U.nunk() == coord[0].size(), "Number of unknowns in solution "
[ - - ][ - - ]
401 : : "vector at recent time step incorrect" );
402 : 3159 : const auto& x = coord[0];
403 : 3159 : const auto& y = coord[1];
404 : 3159 : const auto& z = coord[2];
405 : : // compute the minimum dt across all elements we own
406 : 3159 : auto mindt = std::numeric_limits< tk::real >::max();
407 : 3159 : auto eps = std::numeric_limits< tk::real >::epsilon();
408 : 3159 : auto large = std::numeric_limits< tk::real >::max();
409 [ + + ]: 1870041 : for (std::size_t e=0; e<inpoel.size()/4; ++e) {
410 : : const std::array< std::size_t, 4 >
411 : 1866882 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
412 : : // compute cubic root of element volume as the characteristic length
413 : : const std::array< real, 3 >
414 : 1866882 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
415 : 1866882 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
416 : 1866882 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
417 : 1866882 : const auto L = std::cbrt( tk::triple( ba, ca, da ) / 6.0 );
418 : : // access solution at element nodes at recent time step
419 [ + - ]: 3733764 : std::vector< std::array< real, 4 > > u( m_ncomp );
420 [ + + ][ + - ]: 3858824 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, m_offset, N );
421 : : // get velocity for problem
422 : 1866882 : const std::array< std::vector<std::array<real,3>>, 4 > vel{{
423 [ + - ]: 1866882 : Problem::prescribedVelocity( m_system, m_ncomp,
424 : : x[N[0]], y[N[0]], z[N[0]], t ),
425 [ + - ]: 1866882 : Problem::prescribedVelocity( m_system, m_ncomp,
426 : : x[N[1]], y[N[1]], z[N[1]], t ),
427 [ + - ]: 1866882 : Problem::prescribedVelocity( m_system, m_ncomp,
428 : : x[N[2]], y[N[2]], z[N[2]], t ),
429 [ + - ]: 1866882 : Problem::prescribedVelocity( m_system, m_ncomp,
430 : : x[N[3]], y[N[3]], z[N[3]], t ) }};
431 : : // compute the maximum length of the characteristic velocity (advection
432 : : // velocity) across the four element nodes
433 : 1866882 : real maxvel = 0.0;
434 [ + + ]: 3858824 : for (ncomp_t c=0; c<m_ncomp; ++c)
435 [ + + ]: 9959710 : for (std::size_t i=0; i<4; ++i) {
436 : 7967768 : auto v = tk::length( vel[i][c] );
437 [ + + ]: 7967768 : if (v > maxvel) maxvel = v;
438 : : }
439 : : // compute element dt for the advection
440 [ + + ]: 1866882 : auto advection_dt = std::abs(maxvel) > eps ? L / maxvel : large;
441 : : // compute element dt based on diffusion
442 [ + - ]: 1866882 : auto diffusion_dt = m_physics.diffusion_dt( m_system, m_ncomp, L, u );
443 : : // compute minimum element dt
444 : 1866882 : auto elemdt = std::min( advection_dt, diffusion_dt );
445 : : // find minimum dt across all elements
446 [ + + ]: 1866882 : if (elemdt < mindt) mindt = elemdt;
447 : : }
448 : 3159 : return mindt * g_inputdeck.get< tag::discr, tag::cfl >();
449 : : }
450 : :
451 : : //! Compute a time step size for each mesh node (for steady time stepping)
452 : 0 : void dt( uint64_t,
453 : : const std::vector< tk::real >&,
454 : : const tk::Fields&,
455 : 0 : std::vector< tk::real >& ) const {}
456 : :
457 : : //! \brief Query Dirichlet boundary condition value on a given side set for
458 : : //! all components in this PDE system
459 : : //! \param[in] t Physical time
460 : : //! \param[in] deltat Time step size
461 : : //! \param[in] tp Physical time for each mesh node
462 : : //! \param[in] dtp Time step size for each mesh node
463 : : //! \param[in] ss Pair of side set ID and list of node IDs on the side set
464 : : //! \param[in] coord Mesh node coordinates
465 : : //! \param[in] increment If true, evaluate the solution increment between
466 : : //! t and t+dt for Dirichlet BCs. If false, evlauate the solution instead.
467 : : //! \return Vector of pairs of bool and boundary condition value associated
468 : : //! to mesh node IDs at which Dirichlet boundary conditions are set. Note
469 : : //! that if increment is true, instead of the actual boundary condition
470 : : //! value, we return the increment between t+deltat and t, since,
471 : : //! depending on client code and solver, that may be what the solution
472 : : //! requires.
473 : : std::map< std::size_t, std::vector< std::pair<bool,real> > >
474 : 13595 : dirbc( real t,
475 : : real deltat,
476 : : const std::vector< tk::real >& tp,
477 : : const std::vector< tk::real >& dtp,
478 : : const std::pair< const int, std::vector< std::size_t > >& ss,
479 : : const std::array< std::vector< real >, 3 >& coord,
480 : : bool increment ) const
481 : : {
482 : : using tag::param; using tag::transport; using tag::bcdir;
483 : : using NodeBC = std::vector< std::pair< bool, real > >;
484 : 13595 : std::map< std::size_t, NodeBC > bc;
485 : 13595 : const auto& ubc = g_inputdeck.get< param, transport, tag::bc, bcdir >();
486 : 13595 : const auto steady = g_inputdeck.get< tag::discr, tag::steady_state >();
487 [ + + ]: 13595 : if (!ubc.empty()) {
488 [ - + ][ - - ]: 760 : Assert( ubc.size() > m_system, "Indexing out of Dirichlet BC eq-vector" );
[ - - ][ - - ]
489 : 760 : const auto& x = coord[0];
490 : 760 : const auto& y = coord[1];
491 : 760 : const auto& z = coord[2];
492 [ + + ]: 5320 : for (const auto& b : ubc[m_system])
493 [ + - ][ + + ]: 4560 : if (std::stoi(b) == ss.first)
494 [ + + ]: 213070 : for (auto n : ss.second) {
495 [ - + ][ - - ]: 212310 : Assert( x.size() > n, "Indexing out of coordinate array" );
[ - - ][ - - ]
496 [ - + ]: 212310 : if (steady) { t = tp[n]; deltat = dtp[n]; }
497 [ + - ][ + - ]: 636930 : const auto s = increment ?
[ + - ][ - - ]
498 [ + - ]: 212310 : solinc( m_system, m_ncomp, x[n], y[n], z[n],
499 : : t, deltat, Problem::initialize ) :
500 [ - - ]: 0 : Problem::initialize( m_system, m_ncomp, x[n], y[n], z[n],
501 : : t+deltat );
502 [ + - ][ + - ]: 212310 : auto& nbc = bc[n] = NodeBC( m_ncomp );
503 [ + + ]: 444520 : for (ncomp_t c=0; c<m_ncomp; ++c)
504 : 232210 : nbc[c] = { true, s[c] };
505 : : }
506 : : }
507 : 13595 : return bc;
508 : : }
509 : :
510 : : //! Set symmetry boundary conditions at nodes
511 : : void
512 : 55061 : symbc(
513 : : tk::Fields&,
514 : : const std::array< std::vector< real >, 3 >&,
515 : : const std::unordered_map< int,
516 : : std::unordered_map< std::size_t,
517 : : std::array< real, 4 > > >&,
518 : 55061 : const std::unordered_set< std::size_t >& ) const {}
519 : :
520 : : //! Set farfield boundary conditions at nodes
521 : 48782 : void farfieldbc(
522 : : tk::Fields&,
523 : : const std::array< std::vector< real >, 3 >&,
524 : : const std::unordered_map< int,
525 : : std::unordered_map< std::size_t,
526 : : std::array< real, 4 > > >&,
527 : 48782 : const std::unordered_set< std::size_t >& ) const {}
528 : :
529 : : //! Apply sponge conditions at boundary nodes (no-op for transport)
530 : 22259 : void sponge( tk::Fields&,
531 : : const std::array< std::vector< real >, 3 >&,
532 : 22259 : const std::unordered_set< std::size_t >& ) const {}
533 : :
534 : : //! Apply user defined time dependent BCs (no-op for transport)
535 : : void
536 : 22259 : timedepbc( tk::real,
537 : : tk::Fields&,
538 : : const std::vector< std::unordered_set< std::size_t > >&,
539 : 22259 : const std::vector< tk::Table<5> >& ) const {}
540 : :
541 : : //! Return analytic field names to be output to file
542 : : //! \return Vector of strings labelling analytic fields output in file
543 : 1065 : std::vector< std::string > analyticFieldNames() const {
544 : 1065 : std::vector< std::string > n;
545 : 1065 : auto depvar = g_inputdeck.get< tag::param, eq, tag::depvar >()[m_system];
546 [ + + ]: 2133 : for (ncomp_t c=0; c<m_ncomp; ++c)
547 [ + - ][ + - ]: 1068 : n.push_back( depvar + std::to_string(c) + "_analytic" );
[ + - ][ + - ]
548 : 1065 : return n;
549 : : }
550 : :
551 : : //! Return surface field names to be output to file
552 : : //! \return Vector of strings labelling surface fields output in file
553 : : //! \details This functions should be written in conjunction with
554 : : //! surfOutput(), which provides the vector of surface fields to be output
555 : 1179 : std::vector< std::string > surfNames() const { return {}; }
556 : :
557 : : //! Return surface field output going to file
558 : : std::vector< std::vector< real > >
559 : 1179 : surfOutput( const std::map< int, std::vector< std::size_t > >&,
560 : 1179 : const tk::Fields& ) const { return {}; }
561 : :
562 : : //! Return time history field names to be output to file
563 : : //! \return Vector of strings labelling time history fields output in file
564 : 0 : std::vector< std::string > histNames() const { return {}; }
565 : :
566 : : //! Return time history field output evaluated at time history points
567 : : std::vector< std::vector< real > >
568 : 0 : histOutput( const std::vector< HistData >&,
569 : : const std::vector< std::size_t >&,
570 : 0 : const tk::Fields& ) const { return {}; }
571 : :
572 : : //! Return names of integral variables to be output to diagnostics file
573 : : //! \return Vector of strings labelling integral variables output
574 : 77 : std::vector< std::string > names() const {
575 : 77 : std::vector< std::string > n;
576 : : const auto& depvar =
577 : 77 : g_inputdeck.get< tag::param, tag::transport, tag::depvar >().
578 [ + - ]: 77 : at(m_system);
579 : : // construct the name of the numerical solution for all components
580 [ + + ]: 155 : for (ncomp_t c=0; c<m_ncomp; ++c)
581 [ + - ][ + - ]: 78 : n.push_back( depvar + std::to_string(c) );
[ + - ]
582 : 77 : return n;
583 : : }
584 : :
585 : : private:
586 : : const Physics m_physics; //!< Physics policy
587 : : const Problem m_problem; //!< Problem policy
588 : : const ncomp_t m_system; //!< Equation system index
589 : : const ncomp_t m_ncomp; //!< Number of components in this PDE
590 : : const ncomp_t m_offset; //!< Offset this PDE operates from
591 : :
592 : : //! \brief Compute/assemble nodal gradients of primitive variables for
593 : : //! ALECG in all points
594 : : //! \param[in] coord Mesh node coordinates
595 : : //! \param[in] inpoel Mesh element connectivity
596 : : //! \param[in] lid Global->local node ids
597 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
598 : : //! global node ids (key)
599 : : //! \param[in] vol Nodal volumes
600 : : //! \param[in] esup Elements surrounding points
601 : : //! \param[in] U Solution vector at recent time step
602 : : //! \param[in] G Nodal gradients of primitive variables in chare-boundary nodes
603 : : //! \return Gradients of primitive variables in all mesh points
604 : : tk::Fields
605 : 21870 : nodegrad( const std::array< std::vector< real >, 3 >& coord,
606 : : const std::vector< std::size_t >& inpoel,
607 : : const std::unordered_map< std::size_t, std::size_t >& lid,
608 : : const std::unordered_map< std::size_t, std::size_t >& bid,
609 : : const std::vector< real >& vol,
610 : : const std::pair< std::vector< std::size_t >,
611 : : std::vector< std::size_t > >& esup,
612 : : const tk::Fields& U,
613 : : const tk::Fields& G ) const
614 : : {
615 : : // allocate storage for nodal gradients of primitive variables
616 : 21870 : tk::Fields Grad( U.nunk(), m_ncomp*3 );
617 : 21870 : Grad.fill( 0.0 );
618 : :
619 : : // access node cooordinates
620 : 21870 : const auto& x = coord[0];
621 : 21870 : const auto& y = coord[1];
622 : 21870 : const auto& z = coord[2];
623 : :
624 : : // compute gradients of primitive variables in points
625 : 21870 : auto npoin = U.nunk();
626 : : #pragma omp simd
627 [ + + ]: 5931810 : for (std::size_t p=0; p<npoin; ++p)
628 [ + + ]: 82269780 : for (auto e : tk::Around(esup,p)) {
629 : : // access node IDs
630 : 76359840 : std::size_t N[4] =
631 : 76359840 : { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
632 : : // compute element Jacobi determinant, J = 6V
633 : 76359840 : real bax = x[N[1]]-x[N[0]];
634 : 76359840 : real bay = y[N[1]]-y[N[0]];
635 : 76359840 : real baz = z[N[1]]-z[N[0]];
636 : 76359840 : real cax = x[N[2]]-x[N[0]];
637 : 76359840 : real cay = y[N[2]]-y[N[0]];
638 : 76359840 : real caz = z[N[2]]-z[N[0]];
639 : 76359840 : real dax = x[N[3]]-x[N[0]];
640 : 76359840 : real day = y[N[3]]-y[N[0]];
641 : 76359840 : real daz = z[N[3]]-z[N[0]];
642 : 76359840 : auto J = tk::triple( bax, bay, baz, cax, cay, caz, dax, day, daz );
643 : 76359840 : auto J24 = J/24.0;
644 : : // shape function derivatives, nnode*ndim [4][3]
645 : : real g[4][3];
646 : 76359840 : tk::crossdiv( cax, cay, caz, dax, day, daz, J,
647 : : g[1][0], g[1][1], g[1][2] );
648 : 76359840 : tk::crossdiv( dax, day, daz, bax, bay, baz, J,
649 : : g[2][0], g[2][1], g[2][2] );
650 : 76359840 : tk::crossdiv( bax, bay, baz, cax, cay, caz, J,
651 : : g[3][0], g[3][1], g[3][2] );
652 [ + + ]: 305439360 : for (std::size_t i=0; i<3; ++i)
653 : 229079520 : g[0][i] = -g[1][i] - g[2][i] - g[3][i];
654 : : // scatter-add gradient contributions to boundary nodes
655 [ + + ]: 152719680 : for (std::size_t c=0; c<m_ncomp; ++c)
656 [ + + ]: 381799200 : for (std::size_t b=0; b<4; ++b)
657 [ + + ]: 1221757440 : for (std::size_t i=0; i<3; ++i)
658 [ + - ][ + - ]: 916318080 : Grad(p,c*3+i,0) += J24 * g[b][i] * U(N[b],c,m_offset);
659 : : }
660 : :
661 : : // put in nodal gradients of chare-boundary points
662 [ + + ]: 2761410 : for (const auto& [g,b] : bid) {
663 [ + - ]: 2739540 : auto i = tk::cref_find( lid, g );
664 [ + + ]: 10958160 : for (ncomp_t c=0; c<Grad.nprop(); ++c)
665 [ + - ][ + - ]: 8218620 : Grad(i,c,0) = G(b,c,0);
666 : : }
667 : :
668 : : // divide weak result in gradients by nodal volume
669 [ + + ]: 5931810 : for (std::size_t p=0; p<npoin; ++p)
670 [ + + ]: 23639760 : for (std::size_t c=0; c<m_ncomp*3; ++c)
671 [ + - ]: 17729820 : Grad(p,c,0) /= vol[p];
672 : :
673 : 21870 : return Grad;
674 : : }
675 : :
676 : : //! \brief Compute MUSCL reconstruction in edge-end points using a MUSCL
677 : : //! procedure with van Leer limiting
678 : : //! \param[in] p Left node id of edge-end
679 : : //! \param[in] q Right node id of edge-end
680 : : //! \param[in] coord Array of nodal coordinates
681 : : //! \param[in] G Gradient of all unknowns in mesh points
682 : : //! \param[in,out] uL Primitive variables at left edge-end point
683 : : //! \param[in,out] uR Primitive variables at right edge-end point
684 : : void
685 : 59805420 : muscl( std::size_t p,
686 : : std::size_t q,
687 : : const tk::UnsMesh::Coords& coord,
688 : : const tk::Fields& G,
689 : : std::vector< real >& uL,
690 : : std::vector< real >& uR ) const
691 : : {
692 [ + - ][ + - ]: 59805420 : Assert( uL.size() == m_ncomp && uR.size() == m_ncomp, "Size mismatch" );
[ - - ][ - - ]
[ - - ]
693 [ - + ][ - - ]: 59805420 : Assert( G.nprop()/3 == m_ncomp, "Size mismatch" );
[ - - ][ - - ]
694 : :
695 : 59805420 : const auto& x = coord[0];
696 : 59805420 : const auto& y = coord[1];
697 : 59805420 : const auto& z = coord[2];
698 : :
699 : : // edge vector
700 : 59805420 : std::array< real, 3 > vw{ x[q]-x[p], y[q]-y[p], z[q]-z[p] };
701 : :
702 : : std::vector< real >
703 [ + - ][ + - ]: 119610840 : delta1( m_ncomp, 0.0 ), delta2( m_ncomp, 0.0 ), delta3( m_ncomp, 0.0 );
[ + - ]
704 : :
705 : : // MUSCL reconstruction of edge-end-point primitive variables
706 [ + + ]: 119610840 : for (std::size_t c=0; c<m_ncomp; ++c) {
707 : : // gradients
708 : : std::array< real, 3 >
709 [ + - ][ + - ]: 59805420 : g1{ G(p,c*3+0,0), G(p,c*3+1,0), G(p,c*3+2,0) },
[ + - ]
710 [ + - ][ + - ]: 59805420 : g2{ G(q,c*3+0,0), G(q,c*3+1,0), G(q,c*3+2,0) };
[ + - ]
711 : :
712 : 59805420 : delta2[c] = uR[c] - uL[c];
713 : 59805420 : delta1[c] = 2.0 * tk::dot(g1,vw) - delta2[c];
714 : 59805420 : delta3[c] = 2.0 * tk::dot(g2,vw) - delta2[c];
715 : :
716 : : // form limiters
717 : 59805420 : auto rL = (delta2[c] + muscl_eps) / (delta1[c] + muscl_eps);
718 : 59805420 : auto rR = (delta2[c] + muscl_eps) / (delta3[c] + muscl_eps);
719 : 59805420 : auto rLinv = (delta1[c] + muscl_eps) / (delta2[c] + muscl_eps);
720 : 59805420 : auto rRinv = (delta3[c] + muscl_eps) / (delta2[c] + muscl_eps);
721 : :
722 : 59805420 : auto phiL = (std::abs(rL) + rL) / (std::abs(rL) + 1.0);
723 : 59805420 : auto phiR = (std::abs(rR) + rR) / (std::abs(rR) + 1.0);
724 : 59805420 : auto phi_L_inv = (std::abs(rLinv) + rLinv) / (std::abs(rLinv) + 1.0);
725 : 59805420 : auto phi_R_inv = (std::abs(rRinv) + rRinv) / (std::abs(rRinv) + 1.0);
726 : :
727 : : // update unknowns with reconstructed unknowns
728 : 59805420 : uL[c] += 0.25*(delta1[c]*muscl_m1*phiL + delta2[c]*muscl_p1*phi_L_inv);
729 : 59805420 : uR[c] -= 0.25*(delta3[c]*muscl_m1*phiR + delta2[c]*muscl_p1*phi_R_inv);
730 : : }
731 : 59805420 : }
732 : :
733 : : //! Compute domain-edge integral for ALECG
734 : : //! \param[in] coord Mesh node coordinates
735 : : //! \param[in] inpoel Mesh element connectivity
736 : : //! \param[in] edgeid Local node id pair -> edge id map
737 : : //! \param[in] psup Points surrounding points
738 : : //! \param[in] dfn Dual-face normals
739 : : //! \param[in] U Solution vector at recent time step
740 : : //! \param[in] G Nodal gradients
741 : : //! \param[in,out] R Right-hand side vector computed
742 : 21870 : void domainint( const std::array< std::vector< real >, 3 >& coord,
743 : : const std::vector< std::size_t >& inpoel,
744 : : const std::vector< std::size_t >& edgeid,
745 : : const std::pair< std::vector< std::size_t >,
746 : : std::vector< std::size_t > >& psup,
747 : : const std::vector< real >& dfn,
748 : : const tk::Fields& U,
749 : : const tk::Fields& G,
750 : : tk::Fields& R ) const
751 : : {
752 : : // access node cooordinates
753 : 21870 : const auto& x = coord[0];
754 : 21870 : const auto& y = coord[1];
755 : 21870 : const auto& z = coord[2];
756 : :
757 : : // compute derived data structures
758 [ + - ][ + - ]: 43740 : auto esued = tk::genEsued( inpoel, 4, tk::genEsup( inpoel, 4 ) );
759 : :
760 : : // access pointer to right hand side at component and offset
761 [ + - ]: 43740 : std::vector< const real* > r( m_ncomp );
762 [ + + ][ + - ]: 43740 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c, m_offset );
763 : :
764 : : // domain-edge integral
765 [ + + ]: 5931810 : for (std::size_t p=0,k=0; p<U.nunk(); ++p) {
766 [ + + ]: 65715360 : for (auto q : tk::Around(psup,p)) {
767 : : // access dual-face normals for edge p-q
768 : 59805420 : auto ed = edgeid[k++];
769 : 59805420 : std::array< tk::real, 3 > n{ dfn[ed*6+0], dfn[ed*6+1], dfn[ed*6+2] };
770 : :
771 [ + - ]: 119610840 : std::vector< tk::real > uL( m_ncomp, 0.0 );
772 [ + - ]: 119610840 : std::vector< tk::real > uR( m_ncomp, 0.0 );
773 [ + + ]: 119610840 : for (std::size_t c=0; c<m_ncomp; ++c) {
774 [ + - ]: 59805420 : uL[c] = U(p,c,m_offset);
775 [ + - ]: 59805420 : uR[c] = U(q,c,m_offset);
776 : : }
777 : : // compute MUSCL reconstruction in edge-end points
778 [ + - ]: 59805420 : muscl( p, q, coord, G, uL, uR );
779 : :
780 : : // evaluate prescribed velocity
781 : 59805420 : auto v =
782 [ + - ]: 59805420 : Problem::prescribedVelocity( m_system, m_ncomp, x[p], y[p], z[p],
783 : : 0.0 );
784 : : // sum donain-edge contributions
785 [ + - ][ + + ]: 288884940 : for (auto e : tk::cref_find(esued,{p,q})) {
786 : : const std::array< std::size_t, 4 >
787 : 229079520 : N{{ inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] }};
788 : : // compute element Jacobi determinant
789 : : const std::array< tk::real, 3 >
790 : 229079520 : ba{{ x[N[1]]-x[N[0]], y[N[1]]-y[N[0]], z[N[1]]-z[N[0]] }},
791 : 229079520 : ca{{ x[N[2]]-x[N[0]], y[N[2]]-y[N[0]], z[N[2]]-z[N[0]] }},
792 : 229079520 : da{{ x[N[3]]-x[N[0]], y[N[3]]-y[N[0]], z[N[3]]-z[N[0]] }};
793 : 229079520 : const auto J = tk::triple( ba, ca, da ); // J = 6V
794 : : // shape function derivatives, nnode*ndim [4][3]
795 : : std::array< std::array< tk::real, 3 >, 4 > grad;
796 : 229079520 : grad[1] = tk::crossdiv( ca, da, J );
797 : 229079520 : grad[2] = tk::crossdiv( da, ba, J );
798 : 229079520 : grad[3] = tk::crossdiv( ba, ca, J );
799 [ + + ]: 916318080 : for (std::size_t i=0; i<3; ++i)
800 : 687238560 : grad[0][i] = -grad[1][i]-grad[2][i]-grad[3][i];
801 : 229079520 : auto J48 = J/48.0;
802 [ + + ]: 1603556640 : for (const auto& [a,b] : tk::lpoed) {
803 [ + - ]: 1374477120 : auto s = tk::orient( {N[a],N[b]}, {p,q} );
804 [ + + ]: 5497908480 : for (std::size_t j=0; j<3; ++j) {
805 [ + + ]: 8246862720 : for (std::size_t c=0; c<m_ncomp; ++c) {
806 [ + - ]: 8246862720 : R.var(r[c],p) -= J48 * s * (grad[a][j] - grad[b][j])
807 : 4123431360 : * v[c][j]*(uL[c] + uR[c])
808 : 4123431360 : - J48 * std::abs(s * (grad[a][j] - grad[b][j]))
809 : 4123431360 : * std::abs(tk::dot(v[c],n)) * (uR[c] - uL[c]);
810 : : }
811 : : }
812 : : }
813 : : }
814 : : }
815 : : }
816 : 21870 : }
817 : :
818 : : //! Compute boundary integrals for ALECG
819 : : //! \param[in] coord Mesh node coordinates
820 : : //! \param[in] triinpoel Boundary triangle face connecitivity with local ids
821 : : //! \param[in] symbctri Vector with 1 at symmetry BC boundary triangles
822 : : //! \param[in] U Solution vector at recent time step
823 : : //! \param[in,out] R Right-hand side vector computed
824 : 21870 : void bndint( const std::array< std::vector< real >, 3 >& coord,
825 : : const std::vector< std::size_t >& triinpoel,
826 : : const std::vector< int >& symbctri,
827 : : const tk::Fields& U,
828 : : tk::Fields& R ) const
829 : : {
830 : : // access node coordinates
831 : 21870 : const auto& x = coord[0];
832 : 21870 : const auto& y = coord[1];
833 : 21870 : const auto& z = coord[2];
834 : :
835 : : // boundary integrals: compute fluxes in edges
836 [ + - ]: 43740 : std::vector< real > bflux( triinpoel.size() * m_ncomp * 2 );
837 : :
838 [ + + ]: 1460670 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
839 : : // access node IDs
840 : : std::array< std::size_t, 3 >
841 : 1438800 : N{ triinpoel[e*3+0], triinpoel[e*3+1], triinpoel[e*3+2] };
842 : : // apply symmetry BCs
843 [ + + ]: 1438800 : if (symbctri[e]) continue;
844 : : // node coordinates
845 : 91500 : std::array< tk::real, 3 > xp{ x[N[0]], x[N[1]], x[N[2]] },
846 : 91500 : yp{ y[N[0]], y[N[1]], y[N[2]] },
847 : 91500 : zp{ z[N[0]], z[N[1]], z[N[2]] };
848 : : // access solution at element nodes
849 [ + - ]: 183000 : std::vector< std::array< real, 3 > > u( m_ncomp );
850 [ + + ][ + - ]: 183000 : for (ncomp_t c=0; c<m_ncomp; ++c) u[c] = U.extract( c, m_offset, N );
851 : : // evaluate prescribed velocity
852 : 91500 : auto v =
853 [ + - ]: 91500 : Problem::prescribedVelocity( m_system, m_ncomp, xp[0], yp[0], zp[0],
854 : : 0.0 );
855 : : // compute face area
856 : 91500 : auto A6 = tk::area( x[N[0]], x[N[1]], x[N[2]],
857 : : y[N[0]], y[N[1]], y[N[2]],
858 : : z[N[0]], z[N[1]], z[N[2]] ) / 6.0;
859 : 91500 : auto A24 = A6/4.0;
860 : : // compute face normal
861 [ + - ]: 91500 : auto n = tk::normal( xp, yp, zp );
862 : : // store flux in boundary elements
863 [ + + ]: 183000 : for (std::size_t c=0; c<m_ncomp; ++c) {
864 : 91500 : auto eb = (e*m_ncomp+c)*6;
865 : 91500 : auto vdotn = tk::dot( v[c], n );
866 : 91500 : auto Bab = A24 * vdotn * (u[c][0] + u[c][1]);
867 : 91500 : bflux[eb+0] = Bab + A6 * vdotn * u[c][0];
868 : 91500 : bflux[eb+1] = Bab;
869 : 91500 : Bab = A24 * vdotn * (u[c][1] + u[c][2]);
870 : 91500 : bflux[eb+2] = Bab + A6 * vdotn * u[c][1];
871 : 91500 : bflux[eb+3] = Bab;
872 : 91500 : Bab = A24 * vdotn * (u[c][2] + u[c][0]);
873 : 91500 : bflux[eb+4] = Bab + A6 * vdotn * u[c][2];
874 : 91500 : bflux[eb+5] = Bab;
875 : : }
876 : : }
877 : :
878 : : // access pointer to right hand side at component and offset
879 [ + - ]: 43740 : std::vector< const real* > r( m_ncomp );
880 [ + + ][ + - ]: 43740 : for (ncomp_t c=0; c<m_ncomp; ++c) r[c] = R.cptr( c, m_offset );
881 : :
882 : : // boundary integrals: sum flux contributions to points
883 [ + + ]: 1460670 : for (std::size_t e=0; e<triinpoel.size()/3; ++e) {
884 : 1438800 : std::size_t N[3] =
885 : 1438800 : { triinpoel[e*3+0], triinpoel[e*3+1], triinpoel[e*3+2] };
886 [ + + ]: 2877600 : for (std::size_t c=0; c<m_ncomp; ++c) {
887 : 1438800 : auto eb = (e*m_ncomp+c)*6;
888 [ + - ]: 1438800 : R.var(r[c],N[0]) -= bflux[eb+0] + bflux[eb+5];
889 [ + - ]: 1438800 : R.var(r[c],N[1]) -= bflux[eb+1] + bflux[eb+2];
890 [ + - ]: 1438800 : R.var(r[c],N[2]) -= bflux[eb+3] + bflux[eb+4];
891 : : }
892 : : }
893 : :
894 : 21870 : tk::destroy(bflux);
895 : 21870 : }
896 : : };
897 : :
898 : : } // cg::
899 : : } // inciter::
900 : :
901 : : #endif // Transport_h
|
