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