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