Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/CompFlow/DGCompFlow.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 Compressible single-material flow using discontinuous Galerkin
9 : : finite elements
10 : : \details This file implements calls to the physics operators governing
11 : : compressible single-material flow using discontinuous Galerkin
12 : : discretizations.
13 : : */
14 : : // *****************************************************************************
15 : : #ifndef DGCompFlow_h
16 : : #define DGCompFlow_h
17 : :
18 : : #include <cmath>
19 : : #include <algorithm>
20 : : #include <unordered_set>
21 : : #include <map>
22 : :
23 : : #include <brigand/algorithms/for_each.hpp>
24 : :
25 : : #include "Macro.hpp"
26 : : #include "Exception.hpp"
27 : : #include "Vector.hpp"
28 : : #include "ContainerUtil.hpp"
29 : : #include "UnsMesh.hpp"
30 : : #include "Inciter/InputDeck/InputDeck.hpp"
31 : : #include "Integrate/Basis.hpp"
32 : : #include "Integrate/Quadrature.hpp"
33 : : #include "Integrate/Initialize.hpp"
34 : : #include "Integrate/Mass.hpp"
35 : : #include "Integrate/Surface.hpp"
36 : : #include "Integrate/Boundary.hpp"
37 : : #include "Integrate/Volume.hpp"
38 : : #include "Integrate/Source.hpp"
39 : : #include "RiemannFactory.hpp"
40 : : #include "EoS/EoS.hpp"
41 : : #include "Reconstruction.hpp"
42 : : #include "Limiter.hpp"
43 : : #include "PrefIndicator.hpp"
44 : :
45 : : namespace inciter {
46 : :
47 : : extern ctr::InputDeck g_inputdeck;
48 : :
49 : : namespace dg {
50 : :
51 : : //! \brief CompFlow used polymorphically with tk::DGPDE
52 : : //! \details The template arguments specify policies and are used to configure
53 : : //! the behavior of the class. The policies are:
54 : : //! - Physics - physics configuration, see PDE/CompFlow/Physics.h
55 : : //! - Problem - problem configuration, see PDE/CompFlow/Problem.h
56 : : //! \note The default physics is Euler, set in inciter::deck::check_compflow()
57 : : template< class Physics, class Problem >
58 : : class CompFlow {
59 : :
60 : : private:
61 : : using eq = tag::compflow;
62 : :
63 : : public:
64 : : //! Constructor
65 : : //! \param[in] c Equation system index (among multiple systems configured)
66 : 80 : explicit CompFlow( ncomp_t c ) :
67 : : m_physics(),
68 : : m_problem(),
69 : : m_system( c ),
70 : 80 : m_ncomp( g_inputdeck.get< tag::component, eq >().at(c) ),
71 : 80 : m_offset( g_inputdeck.get< tag::component >().offset< eq >(c) ),
72 : : m_riemann(tk::cref_find(compflowRiemannSolvers(),
73 [ + - ][ + - ]: 240 : g_inputdeck.get< tag::param, tag::compflow, tag::flux >().at(m_system)))
[ + - ]
74 : : {
75 : : // associate boundary condition configurations with state functions, the
76 : : // order in which the state functions listed matters, see ctr::bc::Keys
77 [ + - ][ + - ]: 560 : brigand::for_each< ctr::bc::Keys >( ConfigBC< eq >( m_system, m_bc,
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ - - ]
78 : : { dirichlet
79 : : , symmetry
80 : : , invalidBC // Inlet BC not implemented
81 : : , invalidBC // Outlet BC not implemented
82 : : , farfield
83 : : , extrapolate } ) );
84 : 80 : }
85 : :
86 : : //! Find the number of primitive quantities required for this PDE system
87 : : //! \return The number of primitive quantities required to be stored for
88 : : //! this PDE system
89 : 212 : std::size_t nprim() const
90 : : {
91 : : // compflow does not need/store any primitive quantities currently
92 : 212 : return 0;
93 : : }
94 : :
95 : : //! Find the number of materials set up for this PDE system
96 : : //! \return The number of materials set up for this PDE system
97 : 0 : std::size_t nmat() const
98 : : {
99 : : // compflow does not need nmat
100 : 0 : return 0;
101 : : }
102 : :
103 : : //! Assign number of DOFs per equation in the PDE system
104 : : //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
105 : : //! equation
106 : 212 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
107 : : {
108 : : // all equation-dofs initialized to ndof
109 [ + + ]: 1272 : for (std::size_t i=0; i<m_ncomp; ++i) {
110 : 1060 : numEqDof.push_back(g_inputdeck.get< tag::discr, tag::ndof >());
111 : : }
112 : 212 : }
113 : :
114 : : //! Determine elements that lie inside the user-defined IC box
115 : : //! \param[in] geoElem Element geometry array
116 : : //! \param[in] nielem Number of internal elements
117 : : //! \param[in,out] inbox List of nodes at which box user ICs are set for
118 : : //! each IC box
119 : 212 : void IcBoxElems( const tk::Fields& geoElem,
120 : : std::size_t nielem,
121 : : std::vector< std::unordered_set< std::size_t > >& inbox ) const
122 : : {
123 : 212 : tk::BoxElems< eq >(m_system, geoElem, nielem, inbox);
124 : 212 : }
125 : :
126 : : //! Initalize the compressible flow equations, prepare for time integration
127 : : //! \param[in] L Block diagonal mass matrix
128 : : //! \param[in] inpoel Element-node connectivity
129 : : //! \param[in] coord Array of nodal coordinates
130 : : //! \param[in] inbox List of elements at which box user ICs are set for
131 : : //! each IC box
132 : : //! \param[in,out] unk Array of unknowns
133 : : //! \param[in] t Physical time
134 : : //! \param[in] nielem Number of internal elements
135 : : void
136 : 220 : initialize( const tk::Fields& L,
137 : : const std::vector< std::size_t >& inpoel,
138 : : const tk::UnsMesh::Coords& coord,
139 : : const std::vector< std::unordered_set< std::size_t > >& inbox,
140 : : tk::Fields& unk,
141 : : tk::real t,
142 : : const std::size_t nielem ) const
143 : : {
144 [ + - ][ + - ]: 220 : tk::initialize( m_system, m_ncomp, m_offset, L, inpoel, coord,
145 : : Problem::initialize, unk, t, nielem );
146 : :
147 : 220 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
148 : 220 : const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
149 : 220 : const auto& icbox = ic.get< tag::box >();
150 : 220 : const auto& bgpreic = ic.get< tag::pressure >();
151 [ + - ]: 220 : auto c_v = cv< eq >(m_system);
152 : :
153 : : // Set initial conditions inside user-defined IC box
154 [ + - ]: 440 : std::vector< tk::real > s(m_ncomp, 0.0);
155 [ + + ]: 102728 : for (std::size_t e=0; e<nielem; ++e) {
156 [ + - ]: 102508 : if (icbox.size() > m_system) {
157 : 102508 : std::size_t bcnt = 0;
158 [ - + ]: 102508 : for (const auto& b : icbox[m_system]) { // for all boxes
159 [ - - ][ - - ]: 0 : if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
[ - - ][ - - ]
160 : : {
161 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
162 : 0 : auto mark = c*rdof;
163 [ - - ]: 0 : s[c] = unk(e,mark,m_offset);
164 : : // set high-order DOFs to zero
165 [ - - ]: 0 : for (std::size_t i=1; i<rdof; ++i)
166 [ - - ]: 0 : unk(e,mark+i,m_offset) = 0.0;
167 : : }
168 [ - - ]: 0 : initializeBox( m_system, 1.0, t, b, bgpreic[m_system][0], c_v,
169 : : s );
170 : : // store box-initialization in solution vector
171 [ - - ]: 0 : for (std::size_t c=0; c<m_ncomp; ++c) {
172 : 0 : auto mark = c*rdof;
173 [ - - ]: 0 : unk(e,mark,m_offset) = s[c];
174 : : }
175 : : }
176 : 0 : ++bcnt;
177 : : }
178 : : }
179 : : }
180 : 220 : }
181 : :
182 : : //! Compute the left hand side block-diagonal mass matrix
183 : : //! \param[in] geoElem Element geometry array
184 : : //! \param[in,out] l Block diagonal mass matrix
185 : 220 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
186 : 220 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
187 : 220 : tk::mass( m_ncomp, m_offset, ndof, geoElem, l );
188 : 220 : }
189 : :
190 : : //! Update the interface cells to first order dofs
191 : : //! \details This function resets the high-order terms in interface cells,
192 : : //! and is currently not used in compflow.
193 : 16320 : void updateInterfaceCells( tk::Fields&,
194 : : std::size_t,
195 : 16320 : std::vector< std::size_t >& ) const {}
196 : :
197 : : //! Update the primitives for this PDE system
198 : : //! \details This function computes and stores the dofs for primitive
199 : : //! quantities, which is currently unused for compflow. But if a limiter
200 : : //! requires primitive variables for example, this would be the place to
201 : : //! add the computation of the primitive variables.
202 : 16532 : void updatePrimitives( const tk::Fields&,
203 : : const tk::Fields&,
204 : : const tk::Fields&,
205 : : tk::Fields&,
206 : 16532 : std::size_t ) const {}
207 : :
208 : : //! Clean up the state of trace materials for this PDE system
209 : : //! \details This function cleans up the state of materials present in trace
210 : : //! quantities in each cell. This is unused for compflow.
211 : 16320 : void cleanTraceMaterial( const tk::Fields&,
212 : : tk::Fields&,
213 : : tk::Fields&,
214 : 16320 : std::size_t ) const {}
215 : :
216 : : //! Reconstruct second-order solution from first-order using least-squares
217 : : //! \param[in] t Physical time
218 : : //! \param[in] geoFace Face geometry array
219 : : //! \param[in] geoElem Element geometry array
220 : : //! \param[in] fd Face connectivity and boundary conditions object
221 : : //! \param[in] inpoel Element-node connectivity
222 : : //! \param[in] coord Array of nodal coordinates
223 : : //! \param[in,out] U Solution vector at recent time step
224 : : //! \param[in,out] P Primitive vector at recent time step
225 : 14760 : void reconstruct( tk::real t,
226 : : const tk::Fields& geoFace,
227 : : const tk::Fields& geoElem,
228 : : const inciter::FaceData& fd,
229 : : const std::map< std::size_t, std::vector< std::size_t > >&,
230 : : const std::vector< std::size_t >& inpoel,
231 : : const tk::UnsMesh::Coords& coord,
232 : : tk::Fields& U,
233 : : tk::Fields& P ) const
234 : : {
235 : 14760 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
236 : :
237 : : // do reconstruction only if P0P1
238 [ + + ][ + + ]: 14760 : if (rdof == 4 && g_inputdeck.get< tag::discr, tag::ndof >() == 1) {
[ + + ]
239 : 300 : const auto nelem = fd.Esuel().size()/4;
240 : :
241 [ - + ][ - - ]: 300 : Assert( U.nprop() == rdof*5, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
242 : : "vector must equal "+ std::to_string(rdof*5) );
243 [ - + ][ - - ]: 300 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
244 : : "Mismatch in inpofa size" );
245 : :
246 : : // allocate and initialize matrix and vector for reconstruction
247 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
248 [ + - ]: 600 : lhs_ls( nelem, {{ {{0.0, 0.0, 0.0}},
249 : : {{0.0, 0.0, 0.0}},
250 : : {{0.0, 0.0, 0.0}} }} );
251 : : std::vector< std::vector< std::array< tk::real, 3 > > >
252 [ + - ]: 900 : rhs_ls( nelem, std::vector< std::array< tk::real, 3 > >
253 [ + - ]: 300 : ( m_ncomp,
254 : : {{ 0.0, 0.0, 0.0 }} ) );
255 : :
256 : : // reconstruct x,y,z-derivatives of unknowns
257 : : // 0. get lhs matrix, which is only geometry dependent
258 [ + - ]: 300 : tk::lhsLeastSq_P0P1(fd, geoElem, geoFace, lhs_ls);
259 : :
260 : : // 1. internal face contributions
261 : 0 : tk::intLeastSq_P0P1( m_offset, rdof, fd, geoElem, U, rhs_ls,
262 [ + - ]: 300 : {0, m_ncomp-1} );
263 : :
264 : : // 2. boundary face contributions
265 [ + + ]: 2100 : for (const auto& b : m_bc)
266 : 1800 : tk::bndLeastSqConservedVar_P0P1( m_system, m_ncomp, m_offset, rdof,
267 : 1800 : b.first, fd, geoFace, geoElem, t, b.second, P, U, rhs_ls,
268 [ + - ]: 1800 : {0, m_ncomp-1} );
269 : :
270 : : // 3. solve 3x3 least-squares system
271 : 0 : tk::solveLeastSq_P0P1( m_offset, rdof, lhs_ls, rhs_ls, U,
272 [ + - ]: 300 : {0, m_ncomp-1} );
273 : :
274 : : // 4. transform reconstructed derivatives to Dubiner dofs
275 : 300 : tk::transform_P0P1( m_offset, rdof, nelem, inpoel, coord, U,
276 [ + - ]: 300 : {0, m_ncomp-1} );
277 : : }
278 : 14760 : }
279 : :
280 : : //! Limit second-order solution
281 : : //! \param[in] t Physical time
282 : : //! \param[in] geoFace Face geometry array
283 : : //! \param[in] geoElem Element geometry array
284 : : //! \param[in] fd Face connectivity and boundary conditions object
285 : : //! \param[in] esup Elements surrounding points
286 : : //! \param[in] inpoel Element-node connectivity
287 : : //! \param[in] coord Array of nodal coordinates
288 : : //! \param[in] ndofel Vector of local number of degrees of freedome
289 : : //! \param[in] gid Local->global node id map
290 : : //! \param[in] bid Local chare-boundary node ids (value) associated to
291 : : //! global node ids (key)
292 : : //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
293 : : //! variables
294 : : //! \param[in,out] U Solution vector at recent time step
295 : 14760 : void limit( [[maybe_unused]] tk::real t,
296 : : [[maybe_unused]] const tk::Fields& geoFace,
297 : : const tk::Fields& geoElem,
298 : : const inciter::FaceData& fd,
299 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
300 : : const std::vector< std::size_t >& inpoel,
301 : : const tk::UnsMesh::Coords& coord,
302 : : const std::vector< std::size_t >& ndofel,
303 : : const std::vector< std::size_t >& gid,
304 : : const std::unordered_map< std::size_t, std::size_t >& bid,
305 : : const std::vector< std::vector<tk::real> >& uNodalExtrm,
306 : : const std::vector< std::vector<tk::real> >&,
307 : : tk::Fields& U,
308 : : tk::Fields&,
309 : : std::vector< std::size_t >& ) const
310 : : {
311 : 14760 : const auto limiter = g_inputdeck.get< tag::discr, tag::limiter >();
312 : 14760 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
313 : :
314 [ - + ]: 14760 : if (limiter == ctr::LimiterType::WENOP1)
315 : 0 : WENO_P1( fd.Esuel(), m_offset, U );
316 [ + + ]: 14760 : else if (limiter == ctr::LimiterType::SUPERBEEP1)
317 : 4530 : Superbee_P1( fd.Esuel(), inpoel, ndofel, m_offset, coord, U );
318 [ + + ][ - + ]: 10230 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 4)
319 : 0 : VertexBasedCompflow_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
320 : 0 : m_offset, geoElem, coord, U);
321 [ + + ][ + - ]: 10230 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 10)
322 : 300 : VertexBasedCompflow_P2( esup, inpoel, ndofel, fd.Esuel().size()/4,
323 : 300 : m_offset, geoElem, coord, gid, bid, uNodalExtrm, U);
324 : 14760 : }
325 : :
326 : : //! Compute right hand side
327 : : //! \param[in] t Physical time
328 : : //! \param[in] geoFace Face geometry array
329 : : //! \param[in] geoElem Element geometry array
330 : : //! \param[in] fd Face connectivity and boundary conditions object
331 : : //! \param[in] inpoel Element-node connectivity
332 : : //! \param[in] boxelems Mesh node ids within user-defined IC boxes
333 : : //! \param[in] coord Array of nodal coordinates
334 : : //! \param[in] U Solution vector at recent time step
335 : : //! \param[in] P Primitive vector at recent time step
336 : : //! \param[in] ndofel Vector of local number of degrees of freedom
337 : : //! \param[in,out] R Right-hand side vector computed
338 : 16320 : void rhs( tk::real t,
339 : : const tk::Fields& geoFace,
340 : : const tk::Fields& geoElem,
341 : : const inciter::FaceData& fd,
342 : : const std::vector< std::size_t >& inpoel,
343 : : const std::vector< std::unordered_set< std::size_t > >& boxelems,
344 : : const tk::UnsMesh::Coords& coord,
345 : : const tk::Fields& U,
346 : : const tk::Fields& P,
347 : : const std::vector< std::size_t >& ndofel,
348 : : tk::Fields& R ) const
349 : : {
350 : 16320 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
351 : 16320 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
352 : :
353 [ - + ][ - - ]: 16320 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
354 : : "vector and primitive vector at recent time step incorrect" );
355 [ - + ][ - - ]: 16320 : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
356 : : "vector and right-hand side at recent time step incorrect" );
357 [ - + ][ - - ]: 16320 : Assert( U.nprop() == rdof*5, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
358 : : "vector must equal "+ std::to_string(rdof*5) );
359 [ - + ][ - - ]: 16320 : Assert( P.nprop() == 0, "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
360 : : "vector must equal "+ std::to_string(0) );
361 [ - + ][ - - ]: 16320 : Assert( R.nprop() == ndof*5, "Number of components in right-hand "
[ - - ][ - - ]
[ - - ]
362 : : "side vector must equal "+ std::to_string(ndof*5) );
363 [ - + ][ - - ]: 16320 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
364 : : "Mismatch in inpofa size" );
365 : :
366 : : // set rhs to zero
367 [ + - ]: 16320 : R.fill(0.0);
368 : :
369 : : // empty vector for non-conservative terms. This vector is unused for
370 : : // single-material hydrodynamics since, there are no non-conservative
371 : : // terms in the system of PDEs.
372 : 32640 : std::vector< std::vector < tk::real > > riemannDeriv;
373 : :
374 : 32640 : std::vector< std::vector< tk::real > > vriem;
375 : 32640 : std::vector< std::vector< tk::real > > riemannLoc;
376 : :
377 : : // configure Riemann flux function
378 : 16320 : auto rieflxfn =
379 : 31881342 : [this]( const std::array< tk::real, 3 >& fn,
380 : : const std::array< std::vector< tk::real >, 2 >& u,
381 : : const std::vector< std::array< tk::real, 3 > >& v )
382 : 31881342 : { return m_riemann.flux( fn, u, v ); };
383 : : // configure a no-op lambda for prescribed velocity
384 : 100313154 : auto velfn = [this]( ncomp_t, ncomp_t, tk::real, tk::real, tk::real,
385 : : tk::real ){
386 [ - - ][ + - ]: 50148417 : return std::vector< std::array< tk::real, 3 > >( m_ncomp ); };
[ + - ][ - - ]
[ + - ][ + - ]
[ - - ][ + - ]
[ + - ][ + - ]
387 : :
388 : : // compute internal surface flux integrals
389 [ + - ][ + - ]: 16320 : tk::surfInt( m_system, 1, m_offset, t, ndof, rdof, inpoel, coord,
[ + - ]
390 : : fd, geoFace, geoElem, rieflxfn, velfn, U, P, ndofel, R,
391 : : vriem, riemannLoc, riemannDeriv );
392 : :
393 : : // compute ptional source term
394 [ + - ][ + - ]: 16320 : tk::srcInt( m_system, m_offset, t, ndof, fd.Esuel().size()/4,
395 : : inpoel, coord, geoElem, Problem::src, ndofel, R );
396 : :
397 [ + + ]: 16320 : if(ndof > 1)
398 : : // compute volume integrals
399 [ + - ][ + - ]: 14460 : tk::volInt( m_system, 1, m_offset, t, ndof, rdof, fd.Esuel().size()/4,
[ + - ]
400 : : inpoel, coord, geoElem, flux, velfn, U, P, ndofel, R );
401 : :
402 : : // compute boundary surface flux integrals
403 [ + + ]: 114240 : for (const auto& b : m_bc)
404 [ + - ][ + - ]: 97920 : tk::bndSurfInt( m_system, 1, m_offset, ndof, rdof, b.first, fd,
405 : : geoFace, geoElem, inpoel, coord, t, rieflxfn, velfn,
406 [ + - ]: 97920 : b.second, U, P, ndofel, R, vriem, riemannLoc,
407 : : riemannDeriv );
408 : :
409 : : // compute external (energy) sources
410 : 16320 : const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
411 : 16320 : const auto& icbox = ic.get< tag::box >();
412 : :
413 [ + - ][ - + ]: 16320 : if (icbox.size() > m_system && !boxelems.empty()) {
[ - + ]
414 : 0 : std::size_t bcnt = 0;
415 [ - - ]: 0 : for (const auto& b : icbox[m_system]) { // for all boxes for this eq
416 [ - - ]: 0 : std::vector< tk::real > box
417 : 0 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
418 : 0 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
419 : 0 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
420 : :
421 : 0 : const auto& initiate = b.template get< tag::initiate >();
422 : 0 : auto inittype = initiate.template get< tag::init >();
423 [ - - ]: 0 : if (inittype == ctr::InitiateType::LINEAR) {
424 [ - - ]: 0 : boxSrc( t, inpoel, boxelems[bcnt], coord, geoElem, ndofel, R );
425 : : }
426 : 0 : ++bcnt;
427 : : }
428 : : }
429 : 16320 : }
430 : :
431 : : //! Evaluate the adaptive indicator and mark the ndof for each element
432 : : //! \param[in] nunk Number of unknowns
433 : : //! \param[in] coord Array of nodal coordinates
434 : : //! \param[in] inpoel Element-node connectivity
435 : : //! \param[in] fd Face connectivity and boundary conditions object
436 : : //! \param[in] unk Array of unknowns
437 : : //! \param[in] indicator p-refinement indicator type
438 : : //! \param[in] ndof Number of degrees of freedom in the solution
439 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
440 : : //! \param[in] tolref Tolerance for p-refinement
441 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
442 : 1636 : void eval_ndof( std::size_t nunk,
443 : : const tk::UnsMesh::Coords& coord,
444 : : const std::vector< std::size_t >& inpoel,
445 : : const inciter::FaceData& fd,
446 : : const tk::Fields& unk,
447 : : inciter::ctr::PrefIndicatorType indicator,
448 : : std::size_t ndof,
449 : : std::size_t ndofmax,
450 : : tk::real tolref,
451 : : std::vector< std::size_t >& ndofel ) const
452 : : {
453 : 1636 : const auto& esuel = fd.Esuel();
454 : :
455 [ + - ]: 1636 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
456 : 1636 : spectral_decay( 1, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel );
457 [ - - ]: 0 : else if(indicator == inciter::ctr::PrefIndicatorType::NON_CONFORMITY)
458 : 0 : non_conformity( nunk, fd.Nbfac(), inpoel, coord, esuel, fd.Esuf(),
459 : : fd.Inpofa(), unk, ndof, ndofmax, ndofel );
460 : : else
461 [ - - ][ - - ]: 0 : Throw( "No such adaptive indicator type" );
[ - - ]
462 : 1636 : }
463 : :
464 : : //! Compute the minimum time step size
465 : : //! \param[in] coord Mesh node coordinates
466 : : //! \param[in] inpoel Mesh element connectivity
467 : : //! \param[in] fd Face connectivity and boundary conditions object
468 : : //! \param[in] geoFace Face geometry array
469 : : //! \param[in] geoElem Element geometry array
470 : : //! \param[in] ndofel Vector of local number of degrees of freedom
471 : : //! \param[in] U Solution vector at recent time step
472 : : //! \return Minimum time step size
473 : 2170 : tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
474 : : const std::vector< std::size_t >& inpoel,
475 : : const inciter::FaceData& fd,
476 : : const tk::Fields& geoFace,
477 : : const tk::Fields& geoElem,
478 : : const std::vector< std::size_t >& ndofel,
479 : : const tk::Fields& U,
480 : : const tk::Fields&,
481 : : const std::size_t /*nielem*/ ) const
482 : : {
483 : 2170 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
484 : :
485 : 2170 : const auto& esuf = fd.Esuf();
486 : 2170 : const auto& inpofa = fd.Inpofa();
487 : :
488 : : tk::real rho, u, v, w, rhoE, p, a, vn, dSV_l, dSV_r;
489 [ + - ]: 2170 : std::vector< tk::real > delt( U.nunk(), 0.0 );
490 : :
491 : 2170 : const auto& cx = coord[0];
492 : 2170 : const auto& cy = coord[1];
493 : 2170 : const auto& cz = coord[2];
494 : :
495 : : // compute internal surface maximum characteristic speed
496 [ + + ]: 1829090 : for (std::size_t f=0; f<esuf.size()/2; ++f)
497 : : {
498 : :
499 : 1826920 : std::size_t el = static_cast< std::size_t >(esuf[2*f]);
500 : 1826920 : auto er = esuf[2*f+1];
501 : :
502 : : // Number of quadrature points for face integration
503 : : std::size_t ng;
504 : :
505 [ + + ]: 1826920 : if(er > -1)
506 : : {
507 : 1346920 : auto eR = static_cast< std::size_t >( er );
508 : :
509 [ + - ]: 1346920 : auto ng_l = tk::NGfa(ndofel[el]);
510 [ + - ]: 1346920 : auto ng_r = tk::NGfa(ndofel[eR]);
511 : :
512 : : // When the number of gauss points for the left and right element are
513 : : // different, choose the larger ng
514 : 1346920 : ng = std::max( ng_l, ng_r );
515 : : }
516 : : else
517 : : {
518 [ + - ]: 480000 : ng = tk::NGfa(ndofel[el]);
519 : : }
520 : :
521 : : // arrays for quadrature points
522 : 3653840 : std::array< std::vector< tk::real >, 2 > coordgp;
523 : 3653840 : std::vector< tk::real > wgp;
524 : :
525 [ + - ]: 1826920 : coordgp[0].resize( ng );
526 [ + - ]: 1826920 : coordgp[1].resize( ng );
527 [ + - ]: 1826920 : wgp.resize( ng );
528 : :
529 : : // get quadrature point weights and coordinates for triangle
530 [ + - ]: 1826920 : tk::GaussQuadratureTri( ng, coordgp, wgp );
531 : :
532 : : // Extract the left element coordinates
533 : 18269200 : std::array< std::array< tk::real, 3>, 4 > coordel_l {{
534 : : {{ cx[inpoel[4*el ]], cy[inpoel[4*el ]], cz[inpoel[4*el ]] }},
535 : 5480760 : {{ cx[inpoel[4*el+1]], cy[inpoel[4*el+1]], cz[inpoel[4*el+1]] }},
536 : 5480760 : {{ cx[inpoel[4*el+2]], cy[inpoel[4*el+2]], cz[inpoel[4*el+2]] }},
537 : 5480760 : {{ cx[inpoel[4*el+3]], cy[inpoel[4*el+3]], cz[inpoel[4*el+3]] }} }};
538 : :
539 : : // Compute the determinant of Jacobian matrix
540 : : auto detT_l =
541 : 1826920 : tk::Jacobian(coordel_l[0], coordel_l[1], coordel_l[2], coordel_l[3]);
542 : :
543 : : // Extract the face coordinates
544 : 12788440 : std::array< std::array< tk::real, 3>, 3 > coordfa {{
545 : : {{ cx[ inpofa[3*f ] ], cy[ inpofa[3*f ] ], cz[ inpofa[3*f ] ] }},
546 : 5480760 : {{ cx[ inpofa[3*f+1] ], cy[ inpofa[3*f+1] ], cz[ inpofa[3*f+1] ] }},
547 : 5480760 : {{ cx[ inpofa[3*f+2] ], cy[ inpofa[3*f+2] ], cz[ inpofa[3*f+2] ] }}
548 : : }};
549 : :
550 : 1826920 : dSV_l = 0.0;
551 : 1826920 : dSV_r = 0.0;
552 : :
553 : : // Gaussian quadrature
554 [ + + ]: 7712514 : for (std::size_t igp=0; igp<ng; ++igp)
555 : : {
556 : : // Compute the coordinates of quadrature point at physical domain
557 [ + - ]: 5885594 : auto gp = tk::eval_gp( igp, coordfa, coordgp );
558 : :
559 : : // Compute the basis function for the left element
560 [ + - ]: 11771188 : auto B_l = tk::eval_basis( ndofel[el],
561 : 5885594 : tk::Jacobian(coordel_l[0], gp, coordel_l[2], coordel_l[3])/detT_l,
562 : 5885594 : tk::Jacobian(coordel_l[0], coordel_l[1], gp, coordel_l[3])/detT_l,
563 : 5885594 : tk::Jacobian(coordel_l[0], coordel_l[1], coordel_l[2], gp)/detT_l );
564 : :
565 [ + - ]: 5885594 : auto wt = wgp[igp] * geoFace(f,0,0);
566 : :
567 : 5885594 : std::array< std::vector< tk::real >, 2 > ugp;
568 : :
569 : : // left element
570 [ + + ]: 35313564 : for (ncomp_t c=0; c<5; ++c)
571 : : {
572 : 29427970 : auto mark = c*rdof;
573 [ + - ][ + - ]: 29427970 : ugp[0].push_back( U(el, mark, m_offset) );
574 : :
575 [ + + ]: 29427970 : if(ndofel[el] > 1) //DG(P1)
576 [ + - ]: 51280170 : ugp[0][c] += U(el, mark+1, m_offset) * B_l[1]
577 [ + - ]: 25640085 : + U(el, mark+2, m_offset) * B_l[2]
578 [ + - ]: 25640085 : + U(el, mark+3, m_offset) * B_l[3];
579 : :
580 [ + + ]: 29427970 : if(ndofel[el] > 4) //DG(P2)
581 [ + - ]: 37901760 : ugp[0][c] += U(el, mark+4, m_offset) * B_l[4]
582 [ + - ]: 18950880 : + U(el, mark+5, m_offset) * B_l[5]
583 [ + - ]: 18950880 : + U(el, mark+6, m_offset) * B_l[6]
584 [ + - ]: 18950880 : + U(el, mark+7, m_offset) * B_l[7]
585 [ + - ]: 18950880 : + U(el, mark+8, m_offset) * B_l[8]
586 [ + - ]: 18950880 : + U(el, mark+9, m_offset) * B_l[9];
587 : : }
588 : :
589 : 5885594 : rho = ugp[0][0];
590 : 5885594 : u = ugp[0][1]/rho;
591 : 5885594 : v = ugp[0][2]/rho;
592 : 5885594 : w = ugp[0][3]/rho;
593 : 5885594 : rhoE = ugp[0][4];
594 [ + - ]: 5885594 : p = eos_pressure< tag::compflow >( m_system, rho, u, v, w, rhoE );
595 : :
596 [ + - ]: 5885594 : a = eos_soundspeed< tag::compflow >( m_system, rho, p );
597 : :
598 [ + - ][ + - ]: 5885594 : vn = u*geoFace(f,1,0) + v*geoFace(f,2,0) + w*geoFace(f,3,0);
[ + - ]
599 : :
600 : 5885594 : dSV_l = wt * (std::fabs(vn) + a);
601 : :
602 : : // right element
603 [ + + ]: 5885594 : if (er > -1) {
604 : :
605 : : // nodal coordinates of the right element
606 : 4309455 : std::size_t eR = static_cast< std::size_t >( er );
607 : :
608 : : // Extract the left element coordinates
609 : 43094550 : std::array< std::array< tk::real, 3>, 4 > coordel_r {{
610 : : {{ cx[inpoel[4*eR ]], cy[inpoel[4*eR ]], cz[inpoel[4*eR ]] }},
611 : 12928365 : {{ cx[inpoel[4*eR+1]], cy[inpoel[4*eR+1]], cz[inpoel[4*eR+1]] }},
612 : 12928365 : {{ cx[inpoel[4*eR+2]], cy[inpoel[4*eR+2]], cz[inpoel[4*eR+2]] }},
613 : 12928365 : {{ cx[inpoel[4*eR+3]], cy[inpoel[4*eR+3]], cz[inpoel[4*eR+3]] }}
614 : : }};
615 : :
616 : : // Compute the determinant of Jacobian matrix
617 : : auto detT_r =
618 : 4309455 : tk::Jacobian(coordel_r[0],coordel_r[1],coordel_r[2],coordel_r[3]);
619 : :
620 : : // Compute the coordinates of quadrature point at physical domain
621 [ + - ]: 4309455 : gp = tk::eval_gp( igp, coordfa, coordgp );
622 : :
623 : : // Compute the basis function for the right element
624 [ + - ]: 4309455 : auto B_r = tk::eval_basis( ndofel[eR],
625 : 4309455 : tk::Jacobian(coordel_r[0],gp,coordel_r[2],coordel_r[3])/detT_r,
626 : 4309455 : tk::Jacobian(coordel_r[0],coordel_r[1],gp,coordel_r[3])/detT_r,
627 : 4309455 : tk::Jacobian(coordel_r[0],coordel_r[1],coordel_r[2],gp)/detT_r );
628 : :
629 [ + + ]: 25856730 : for (ncomp_t c=0; c<5; ++c)
630 : : {
631 : 21547275 : auto mark = c*rdof;
632 [ + - ][ + - ]: 21547275 : ugp[1].push_back( U(eR, mark, m_offset) );
633 : :
634 [ + + ]: 21547275 : if(ndofel[eR] > 1) //DG(P1)
635 [ + - ]: 37396530 : ugp[1][c] += U(eR, mark+1, m_offset) * B_r[1]
636 [ + - ]: 18698265 : + U(eR, mark+2, m_offset) * B_r[2]
637 [ + - ]: 18698265 : + U(eR, mark+3, m_offset) * B_r[3];
638 : :
639 [ + + ]: 21547275 : if(ndofel[eR] > 4) //DG(P2)
640 [ + - ]: 27640860 : ugp[1][c] += U(eR, mark+4, m_offset) * B_r[4]
641 [ + - ]: 13820430 : + U(eR, mark+5, m_offset) * B_r[5]
642 [ + - ]: 13820430 : + U(eR, mark+6, m_offset) * B_r[6]
643 [ + - ]: 13820430 : + U(eR, mark+7, m_offset) * B_r[7]
644 [ + - ]: 13820430 : + U(eR, mark+8, m_offset) * B_r[8]
645 [ + - ]: 13820430 : + U(eR, mark+9, m_offset) * B_r[9];
646 : : }
647 : :
648 : 4309455 : rho = ugp[1][0];
649 : 4309455 : u = ugp[1][1]/rho;
650 : 4309455 : v = ugp[1][2]/rho;
651 : 4309455 : w = ugp[1][3]/rho;
652 : 4309455 : rhoE = ugp[1][4];
653 [ + - ]: 4309455 : p = eos_pressure< tag::compflow >( m_system, rho, u, v, w, rhoE );
654 [ + - ]: 4309455 : a = eos_soundspeed< tag::compflow >( m_system, rho, p );
655 : :
656 [ + - ][ + - ]: 4309455 : vn = u*geoFace(f,1,0) + v*geoFace(f,2,0) + w*geoFace(f,3,0);
[ + - ]
657 : :
658 : 4309455 : dSV_r = wt * (std::fabs(vn) + a);
659 : 4309455 : delt[eR] += std::max( dSV_l, dSV_r );
660 : : }
661 : :
662 : 5885594 : delt[el] += std::max( dSV_l, dSV_r );
663 : : }
664 : : }
665 : :
666 : 2170 : tk::real mindt = std::numeric_limits< tk::real >::max();
667 : 2170 : tk::real dgp = 0.0;
668 : :
669 : : // compute allowable dt
670 [ + + ]: 764730 : for (std::size_t e=0; e<fd.Esuel().size()/4; ++e)
671 : : {
672 : 762560 : dgp = 0.0;
673 [ + + ]: 762560 : if (ndofel[e] == 4)
674 : : {
675 : 187931 : dgp = 1.0;
676 : : }
677 [ + + ]: 574629 : else if (ndofel[e] == 10)
678 : : {
679 : 268370 : dgp = 2.0;
680 : : }
681 : :
682 : : // Scale smallest dt with CFL coefficient and the CFL is scaled by (2*p+1)
683 : : // where p is the order of the DG polynomial by linear stability theory.
684 [ + - ]: 762560 : mindt = std::min( mindt, geoElem(e,0,0)/ (delt[e] * (2.0*dgp + 1.0)) );
685 : : }
686 : :
687 : 4340 : return mindt;
688 : : }
689 : :
690 : : //! Extract the velocity field at cell nodes. Currently unused.
691 : : //! \param[in] U Solution vector at recent time step
692 : : //! \param[in] N Element node indices
693 : : //! \return Array of the four values of the velocity field
694 : : std::array< std::array< tk::real, 4 >, 3 >
695 : : velocity( const tk::Fields& U,
696 : : const std::array< std::vector< tk::real >, 3 >&,
697 : : const std::array< std::size_t, 4 >& N ) const
698 : : {
699 : : std::array< std::array< tk::real, 4 >, 3 > v;
700 : : v[0] = U.extract( 1, m_offset, N );
701 : : v[1] = U.extract( 2, m_offset, N );
702 : : v[2] = U.extract( 3, m_offset, N );
703 : : auto r = U.extract( 0, m_offset, N );
704 : : std::transform( r.begin(), r.end(), v[0].begin(), v[0].begin(),
705 : : []( tk::real s, tk::real& d ){ return d /= s; } );
706 : : std::transform( r.begin(), r.end(), v[1].begin(), v[1].begin(),
707 : : []( tk::real s, tk::real& d ){ return d /= s; } );
708 : : std::transform( r.begin(), r.end(), v[2].begin(), v[2].begin(),
709 : : []( tk::real s, tk::real& d ){ return d /= s; } );
710 : : return v;
711 : : }
712 : :
713 : : //! Return analytic field names to be output to file
714 : : //! \return Vector of strings labelling analytic fields output in file
715 : 588 : std::vector< std::string > analyticFieldNames() const
716 : 588 : { return m_problem.analyticFieldNames( m_ncomp ); }
717 : :
718 : : //! Return time history field names to be output to file
719 : : //! \return Vector of strings labeling time history fields output in file
720 : 0 : std::vector< std::string > histNames() const
721 : 0 : { return CompFlowHistNames(); }
722 : :
723 : : //! Return surface field output going to file
724 : : std::vector< std::vector< tk::real > >
725 : 0 : surfOutput( const std::map< int, std::vector< std::size_t > >&,
726 : : tk::Fields& ) const
727 : : {
728 : 0 : std::vector< std::vector< tk::real > > s; // punt for now
729 : 0 : return s;
730 : : }
731 : :
732 : : //! Return time history field output evaluated at time history points
733 : : //! \param[in] h History point data
734 : : //! \param[in] inpoel Element-node connectivity
735 : : //! \param[in] coord Array of nodal coordinates
736 : : //! \param[in] U Array of unknowns
737 : : std::vector< std::vector< tk::real > >
738 : 0 : histOutput( const std::vector< HistData >& h,
739 : : const std::vector< std::size_t >& inpoel,
740 : : const tk::UnsMesh::Coords& coord,
741 : : const tk::Fields& U,
742 : : const tk::Fields& ) const
743 : : {
744 : 0 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
745 : :
746 : 0 : const auto& x = coord[0];
747 : 0 : const auto& y = coord[1];
748 : 0 : const auto& z = coord[2];
749 : :
750 [ - - ]: 0 : std::vector< std::vector< tk::real > > Up(h.size());
751 : :
752 : 0 : std::size_t j = 0;
753 [ - - ]: 0 : for (const auto& p : h) {
754 : 0 : auto e = p.get< tag::elem >();
755 : 0 : auto chp = p.get< tag::coord >();
756 : :
757 : : // Evaluate inverse Jacobian
758 : 0 : std::array< std::array< tk::real, 3>, 4 > cp{{
759 : : {{ x[inpoel[4*e ]], y[inpoel[4*e ]], z[inpoel[4*e ]] }},
760 : 0 : {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
761 : 0 : {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
762 : 0 : {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
763 : 0 : auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
764 : :
765 : : // evaluate solution at history-point
766 : 0 : std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
767 : 0 : chp[2]-cp[0][2]}};
768 [ - - ]: 0 : auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
769 : 0 : tk::dot(J[2],dc));
770 [ - - ]: 0 : auto uhp = eval_state(m_ncomp, 0, rdof, rdof, e, U, B, {0, m_ncomp-1});
771 : :
772 : : // store solution in history output vector
773 [ - - ]: 0 : Up[j].resize(6, 0.0);
774 : 0 : Up[j][0] = uhp[0];
775 : 0 : Up[j][1] = uhp[1]/uhp[0];
776 : 0 : Up[j][2] = uhp[2]/uhp[0];
777 : 0 : Up[j][3] = uhp[3]/uhp[0];
778 : 0 : Up[j][4] = uhp[4]/uhp[0];
779 [ - - ]: 0 : Up[j][5] = eos_pressure< tag::compflow > (m_system, uhp[0],
780 : 0 : uhp[1]/uhp[0], uhp[2]/uhp[0], uhp[3]/uhp[0], uhp[4] );
781 : 0 : ++j;
782 : : }
783 : :
784 : 0 : return Up;
785 : : }
786 : :
787 : : //! Return names of integral variables to be output to diagnostics file
788 : : //! \return Vector of strings labelling integral variables output
789 : 35 : std::vector< std::string > names() const
790 : 35 : { return m_problem.names( m_ncomp ); }
791 : :
792 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
793 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
794 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
795 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
796 : : //! \param[in] t Physical time at which to evaluate the analytic solution
797 : : //! \return Vector of analytic solution at given location and time
798 : : std::vector< tk::real >
799 : 348546 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
800 : 348546 : { return Problem::analyticSolution( m_system, m_ncomp, xi, yi, zi, t ); }
801 : :
802 : : //! Return analytic solution for conserved variables
803 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
804 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
805 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
806 : : //! \param[in] t Physical time at which to evaluate the analytic solution
807 : : //! \return Vector of analytic solution at given location and time
808 : : std::vector< tk::real >
809 : 2759629 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
810 : 2759629 : { return Problem::initialize( m_system, m_ncomp, xi, yi, zi, t ); }
811 : :
812 : : private:
813 : : //! Physics policy
814 : : const Physics m_physics;
815 : : //! Problem policy
816 : : const Problem m_problem;
817 : : //! Equation system index
818 : : const ncomp_t m_system;
819 : : //! Number of components in this PDE system
820 : : const ncomp_t m_ncomp;
821 : : //! Offset PDE system operates from
822 : : const ncomp_t m_offset;
823 : : //! Riemann solver
824 : : RiemannSolver m_riemann;
825 : : //! BC configuration
826 : : BCStateFn m_bc;
827 : :
828 : : //! Evaluate physical flux function for this PDE system
829 : : //! \param[in] system Equation system index
830 : : //! \param[in] ncomp Number of scalar components in this PDE system
831 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
832 : : //! evaluate the flux
833 : : //! \return Flux vectors for all components in this PDE system
834 : : //! \note The function signature must follow tk::FluxFn
835 : : static tk::FluxFn::result_type
836 : 18267075 : flux( ncomp_t system,
837 : : [[maybe_unused]] ncomp_t ncomp,
838 : : const std::vector< tk::real >& ugp,
839 : : const std::vector< std::array< tk::real, 3 > >& )
840 : : {
841 [ - + ][ - - ]: 18267075 : Assert( ugp.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
842 : :
843 : 18267075 : auto u = ugp[1] / ugp[0];
844 : 18267075 : auto v = ugp[2] / ugp[0];
845 : 18267075 : auto w = ugp[3] / ugp[0];
846 : : auto p =
847 : 18267075 : eos_pressure< tag::compflow >( system, ugp[0], u, v, w, ugp[4] );
848 : :
849 [ + - ]: 18267075 : std::vector< std::array< tk::real, 3 > > fl( ugp.size() );
850 : :
851 : 18267075 : fl[0][0] = ugp[1];
852 : 18267075 : fl[1][0] = ugp[1] * u + p;
853 : 18267075 : fl[2][0] = ugp[1] * v;
854 : 18267075 : fl[3][0] = ugp[1] * w;
855 : 18267075 : fl[4][0] = u * (ugp[4] + p);
856 : :
857 : 18267075 : fl[0][1] = ugp[2];
858 : 18267075 : fl[1][1] = ugp[2] * u;
859 : 18267075 : fl[2][1] = ugp[2] * v + p;
860 : 18267075 : fl[3][1] = ugp[2] * w;
861 : 18267075 : fl[4][1] = v * (ugp[4] + p);
862 : :
863 : 18267075 : fl[0][2] = ugp[3];
864 : 18267075 : fl[1][2] = ugp[3] * u;
865 : 18267075 : fl[2][2] = ugp[3] * v;
866 : 18267075 : fl[3][2] = ugp[3] * w + p;
867 : 18267075 : fl[4][2] = w * (ugp[4] + p);
868 : :
869 : 18267075 : return fl;
870 : : }
871 : :
872 : : //! \brief Boundary state function providing the left and right state of a
873 : : //! face at Dirichlet boundaries
874 : : //! \param[in] system Equation system index
875 : : //! \param[in] ncomp Number of scalar components in this PDE system
876 : : //! \param[in] ul Left (domain-internal) state
877 : : //! \param[in] x X-coordinate at which to compute the states
878 : : //! \param[in] y Y-coordinate at which to compute the states
879 : : //! \param[in] z Z-coordinate at which to compute the states
880 : : //! \param[in] t Physical time
881 : : //! \return Left and right states for all scalar components in this PDE
882 : : //! system
883 : : //! \note The function signature must follow tk::StateFn
884 : : static tk::StateFn::result_type
885 : 3233550 : dirichlet( ncomp_t system, ncomp_t ncomp, const std::vector< tk::real >& ul,
886 : : tk::real x, tk::real y, tk::real z, tk::real t,
887 : : const std::array< tk::real, 3 >& )
888 : : {
889 : 3233550 : return {{ ul, Problem::initialize( system, ncomp, x, y, z, t ) }};
890 : : }
891 : :
892 : : //! \brief Boundary state function providing the left and right state of a
893 : : //! face at symmetry boundaries
894 : : //! \param[in] ul Left (domain-internal) state
895 : : //! \param[in] fn Unit face normal
896 : : //! \return Left and right states for all scalar components in this PDE
897 : : //! system
898 : : //! \note The function signature must follow tk::StateFn
899 : : static tk::StateFn::result_type
900 : 4748067 : symmetry( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
901 : : tk::real, tk::real, tk::real, tk::real,
902 : : const std::array< tk::real, 3 >& fn )
903 : : {
904 [ + - ]: 9496134 : std::vector< tk::real > ur(5);
905 : : // Internal cell velocity components
906 : 4748067 : auto v1l = ul[1]/ul[0];
907 : 4748067 : auto v2l = ul[2]/ul[0];
908 : 4748067 : auto v3l = ul[3]/ul[0];
909 : : // Normal component of velocity
910 : 4748067 : auto vnl = v1l*fn[0] + v2l*fn[1] + v3l*fn[2];
911 : : // Ghost state velocity components
912 : 4748067 : auto v1r = v1l - 2.0*vnl*fn[0];
913 : 4748067 : auto v2r = v2l - 2.0*vnl*fn[1];
914 : 4748067 : auto v3r = v3l - 2.0*vnl*fn[2];
915 : : // Boundary condition
916 : 4748067 : ur[0] = ul[0];
917 : 4748067 : ur[1] = ur[0] * v1r;
918 : 4748067 : ur[2] = ur[0] * v2r;
919 : 4748067 : ur[3] = ur[0] * v3r;
920 : 4748067 : ur[4] = ul[4];
921 [ + - ]: 9496134 : return {{ std::move(ul), std::move(ur) }};
922 : : }
923 : :
924 : : //! \brief Boundary state function providing the left and right state of a
925 : : //! face at farfield boundaries
926 : : //! \param[in] system Equation system index
927 : : //! \param[in] ul Left (domain-internal) state
928 : : //! \param[in] fn Unit face normal
929 : : //! \return Left and right states for all scalar components in this PDE
930 : : //! system
931 : : //! \note The function signature must follow tk::StateFn
932 : : static tk::StateFn::result_type
933 : 40800 : farfield( ncomp_t system, ncomp_t, const std::vector< tk::real >& ul,
934 : : tk::real, tk::real, tk::real, tk::real,
935 : : const std::array< tk::real, 3 >& fn )
936 : : {
937 : : using tag::param; using tag::bc;
938 : :
939 : : // Primitive variables from farfield
940 : 40800 : auto frho = g_inputdeck.get< param, eq,
941 : 40800 : tag::farfield_density >()[ system ];
942 : 40800 : auto fp = g_inputdeck.get< param, eq,
943 : 40800 : tag::farfield_pressure >()[ system ];
944 [ + - ]: 81600 : auto fu = g_inputdeck.get< param, eq,
945 : 40800 : tag::farfield_velocity >()[ system ];
946 : :
947 : : // Speed of sound from farfield
948 [ + - ]: 40800 : auto fa = eos_soundspeed< eq >( system, frho, fp );
949 : :
950 : : // Normal component from farfield
951 : 40800 : auto fvn = fu[0]*fn[0] + fu[1]*fn[1] + fu[2]*fn[2];
952 : :
953 : : // Mach number from farfield
954 : 40800 : auto fM = fvn / fa;
955 : :
956 : : // Specific total energy from farfield
957 : : auto frhoE =
958 [ + - ]: 40800 : eos_totalenergy< eq >( system, frho, fu[0], fu[1], fu[2], fp );
959 : :
960 : : // Pressure from internal cell
961 [ + - ]: 81600 : auto p = eos_pressure< eq >( system, ul[0], ul[1]/ul[0], ul[2]/ul[0],
962 : 40800 : ul[3]/ul[0], ul[4] );
963 : :
964 [ + - ]: 81600 : auto ur = ul;
965 : :
966 [ - + ]: 40800 : if(fM <= -1) // Supersonic inflow
967 : : {
968 : : // For supersonic inflow, all the characteristics are from outside.
969 : : // Therefore, we calculate the ghost cell state using the primitive
970 : : // variables from outside.
971 : 0 : ur[0] = frho;
972 : 0 : ur[1] = frho * fu[0];
973 : 0 : ur[2] = frho * fu[1];
974 : 0 : ur[3] = frho * fu[2];
975 : 0 : ur[4] = frhoE;
976 [ + - ][ - + ]: 40800 : } else if(fM > -1 && fM < 0) // Subsonic inflow
977 : : {
978 : : // For subsonic inflow, there are 1 outgoing characteristcs and 4
979 : : // incoming characteristic. Therefore, we calculate the ghost cell state
980 : : // by taking pressure from the internal cell and other quantities from
981 : : // the outside.
982 : 0 : ur[0] = frho;
983 : 0 : ur[1] = frho * fu[0];
984 : 0 : ur[2] = frho * fu[1];
985 : 0 : ur[3] = frho * fu[2];
986 : 0 : ur[4] =
987 [ - - ]: 0 : eos_totalenergy< eq >( system, frho, fu[0], fu[1], fu[2], p );
988 [ + - ][ + - ]: 40800 : } else if(fM >= 0 && fM < 1) // Subsonic outflow
989 : : {
990 : : // For subsonic outflow, there are 1 incoming characteristcs and 4
991 : : // outgoing characteristic. Therefore, we calculate the ghost cell state
992 : : // by taking pressure from the outside and other quantities from the
993 : : // internal cell.
994 [ + - ]: 40800 : ur[4] = eos_totalenergy< eq >( system, ul[0], ul[1]/ul[0], ul[2]/ul[0],
995 : 40800 : ul[3]/ul[0], fp );
996 : : }
997 : : // Otherwise, for supersonic outflow, all the characteristics are from
998 : : // internal cell. Therefore, we calculate the ghost cell state using the
999 : : // conservative variables from outside.
1000 : :
1001 [ + - ][ + - ]: 81600 : return {{ ul, ur }};
1002 : : }
1003 : :
1004 : : //! \brief Boundary state function providing the left and right state of a
1005 : : //! face at extrapolation boundaries
1006 : : //! \param[in] ul Left (domain-internal) state
1007 : : //! \return Left and right states for all scalar components in this PDE
1008 : : //! system
1009 : : //! \note The function signature must follow tk::StateFn
1010 : : static tk::StateFn::result_type
1011 : 219360 : extrapolate( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
1012 : : tk::real, tk::real, tk::real, tk::real,
1013 : : const std::array< tk::real, 3 >& )
1014 : : {
1015 : 219360 : return {{ ul, ul }};
1016 : : }
1017 : :
1018 : : //! Compute sources corresponding to a propagating front in user-defined box
1019 : : //! \param[in] t Physical time
1020 : : //! \param[in] inpoel Element point connectivity
1021 : : //! \param[in] boxelems Mesh node ids within user-defined box
1022 : : //! \param[in] coord Mesh node coordinates
1023 : : //! \param[in] geoElem Element geometry array
1024 : : //! \param[in] ndofel Vector of local number of degrees of freedome
1025 : : //! \param[in] R Right-hand side vector
1026 : : //! \details This function add the energy source corresponding to a planar
1027 : : //! wave-front propagating along the z-direction with a user-specified
1028 : : //! velocity, within a box initial condition, configured by the user.
1029 : : //! Example (SI) units of the quantities involved:
1030 : : //! * internal energy content (energy per unit volume): J/m^3
1031 : : //! * specific energy (internal energy per unit mass): J/kg
1032 : 0 : void boxSrc( tk::real t,
1033 : : const std::vector< std::size_t >& inpoel,
1034 : : const std::unordered_set< std::size_t >& boxelems,
1035 : : const tk::UnsMesh::Coords& coord,
1036 : : const tk::Fields& geoElem,
1037 : : const std::vector< std::size_t >& ndofel,
1038 : : tk::Fields& R ) const
1039 : : {
1040 : 0 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
1041 : 0 : const auto& ic = g_inputdeck.get< tag::param, eq, tag::ic >();
1042 : 0 : const auto& icbox = ic.get< tag::box >();
1043 : :
1044 [ - - ]: 0 : if (icbox.size() > m_system) {
1045 [ - - ]: 0 : for (const auto& b : icbox[m_system]) { // for all boxes for this eq
1046 [ - - ]: 0 : std::vector< tk::real > box
1047 : 0 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
1048 : 0 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
1049 : 0 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
1050 : :
1051 : 0 : auto boxenc = b.template get< tag::energy_content >();
1052 [ - - ][ - - ]: 0 : Assert( boxenc > 0.0, "Box energy content must be nonzero" );
[ - - ][ - - ]
1053 : :
1054 : 0 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
1055 : :
1056 : : // determine times at which sourcing is initialized and terminated
1057 : 0 : const auto& initiate = b.template get< tag::initiate >();
1058 : 0 : auto iv = initiate.template get< tag::velocity >();
1059 : 0 : auto wFront = 0.1;
1060 : 0 : auto tInit = 0.0;
1061 : 0 : auto tFinal = tInit + (box[5] - box[4] - 2.0*wFront) / std::fabs(iv);
1062 : 0 : auto aBox = (box[1]-box[0]) * (box[3]-box[2]);
1063 : :
1064 : 0 : const auto& cx = coord[0];
1065 : 0 : const auto& cy = coord[1];
1066 : 0 : const auto& cz = coord[2];
1067 : :
1068 [ - - ][ - - ]: 0 : if (t >= tInit && t <= tFinal) {
1069 : : // The energy front is assumed to have a half-sine-wave shape. The
1070 : : // half wave-length is the width of the front. At t=0, the center of
1071 : : // this front (i.e. the peak of the partial-sine-wave) is at X_0 +
1072 : : // W_0. W_0 is calculated based on the width of the front and the
1073 : : // direction of propagation (which is assumed to be along the
1074 : : // z-direction). If the front propagation velocity is positive, it
1075 : : // is assumed that the initial position of the energy source is the
1076 : : // minimum z-coordinate of the box; whereas if this velocity is
1077 : : // negative, the initial position is the maximum z-coordinate of the
1078 : : // box.
1079 : :
1080 : : // initial center of front
1081 : 0 : tk::real zInit(box[4]);
1082 [ - - ]: 0 : if (iv < 0.0) zInit = box[5];
1083 : : // current location of front
1084 : 0 : auto z0 = zInit + iv*t;
1085 : 0 : auto z1 = z0 + std::copysign(wFront, iv);
1086 : 0 : tk::real s0(z0), s1(z1);
1087 : : // if velocity of propagation is negative, initial position is z1
1088 [ - - ]: 0 : if (iv < 0.0) {
1089 : 0 : s0 = z1;
1090 : 0 : s1 = z0;
1091 : : }
1092 : : // Sine-wave (positive part of the wave) source term amplitude
1093 : 0 : auto pi = 4.0 * std::atan(1.0);
1094 : 0 : auto amplE = boxenc * V_ex * pi
1095 : 0 : / (aBox * wFront * 2.0 * (tFinal-tInit));
1096 : : //// Square wave (constant) source term amplitude
1097 : : //auto amplE = boxenc * V_ex
1098 : : // / (aBox * wFront * (tFinal-tInit));
1099 : :
1100 : : // add source
1101 [ - - ]: 0 : for (auto e : boxelems) {
1102 [ - - ]: 0 : auto zc = geoElem(e,3,0);
1103 : :
1104 [ - - ][ - - ]: 0 : if (zc >= s0 && zc <= s1) {
1105 [ - - ]: 0 : auto ng = tk::NGvol(ndofel[e]);
1106 : :
1107 : : // arrays for quadrature points
1108 : 0 : std::array< std::vector< tk::real >, 3 > coordgp;
1109 : 0 : std::vector< tk::real > wgp;
1110 : :
1111 [ - - ]: 0 : coordgp[0].resize( ng );
1112 [ - - ]: 0 : coordgp[1].resize( ng );
1113 [ - - ]: 0 : coordgp[2].resize( ng );
1114 [ - - ]: 0 : wgp.resize( ng );
1115 : :
1116 [ - - ]: 0 : tk::GaussQuadratureTet( ng, coordgp, wgp );
1117 : :
1118 : : // Extract the element coordinates
1119 : 0 : std::array< std::array< tk::real, 3>, 4 > coordel{{
1120 : : {{ cx[inpoel[4*e ]], cy[inpoel[4*e ]], cz[inpoel[4*e ]] }},
1121 : 0 : {{ cx[inpoel[4*e+1]], cy[inpoel[4*e+1]], cz[inpoel[4*e+1]] }},
1122 : 0 : {{ cx[inpoel[4*e+2]], cy[inpoel[4*e+2]], cz[inpoel[4*e+2]] }},
1123 : 0 : {{ cx[inpoel[4*e+3]], cy[inpoel[4*e+3]], cz[inpoel[4*e+3]] }}}};
1124 : :
1125 [ - - ]: 0 : for (std::size_t igp=0; igp<ng; ++igp) {
1126 : : // Compute the coordinates of quadrature point at physical
1127 : : // domain
1128 [ - - ]: 0 : auto gp = tk::eval_gp( igp, coordel, coordgp );
1129 : :
1130 : : // Compute the basis function
1131 [ - - ]: 0 : auto B = tk::eval_basis( ndofel[e], coordgp[0][igp],
1132 : : coordgp[1][igp], coordgp[2][igp] );
1133 : :
1134 : : // Compute the source term variable
1135 : 0 : std::array< tk::real, 5 > s{{0.0, 0.0, 0.0, 0.0, 0.0}};
1136 : 0 : s[4] = amplE * std::sin(pi*(gp[2]-s0)/wFront);
1137 : :
1138 [ - - ]: 0 : auto wt = wgp[igp] * geoElem(e, 0, 0);
1139 : :
1140 [ - - ]: 0 : tk::update_rhs( m_offset, ndof, ndofel[e], wt, e, B, s, R );
1141 : : }
1142 : : }
1143 : : }
1144 : : }
1145 : : }
1146 : : }
1147 : 0 : }
1148 : : };
1149 : :
1150 : : } // dg::
1151 : :
1152 : : } // inciter::
1153 : :
1154 : : #endif // DGCompFlow_h
|
