Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/Transport/DGTransport.hpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Scalar transport using disccontinous Galerkin discretization
9 : : \details This file implements the physics operators governing transported
10 : : scalars using disccontinuous Galerkin discretization.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef DGTransport_h
14 : : #define DGTransport_h
15 : :
16 : : #include <vector>
17 : : #include <array>
18 : : #include <limits>
19 : : #include <cmath>
20 : : #include <unordered_set>
21 : : #include <map>
22 : :
23 : : #include "Macro.hpp"
24 : : #include "Exception.hpp"
25 : : #include "Vector.hpp"
26 : : #include "UnsMesh.hpp"
27 : : #include "Integrate/Basis.hpp"
28 : : #include "Integrate/Quadrature.hpp"
29 : : #include "Integrate/Initialize.hpp"
30 : : #include "Integrate/Mass.hpp"
31 : : #include "Integrate/Surface.hpp"
32 : : #include "Integrate/Boundary.hpp"
33 : : #include "Integrate/Volume.hpp"
34 : : #include "Riemann/Upwind.hpp"
35 : : #include "Reconstruction.hpp"
36 : : #include "Limiter.hpp"
37 : : #include "PrefIndicator.hpp"
38 : : #include "EoS/EOS.hpp"
39 : : #include "FunctionPrototypes.hpp"
40 : : #include "ConfigureTransport.hpp"
41 : :
42 : : namespace inciter {
43 : :
44 : : extern ctr::InputDeck g_inputdeck;
45 : :
46 : : namespace dg {
47 : :
48 : : //! \brief Transport equation used polymorphically with tk::DGPDE
49 : : //! \details The template argument(s) specify policies and are used to configure
50 : : //! the behavior of the class. The policies are:
51 : : //! - Physics - physics configuration, see PDE/Transport/Physics.h
52 : : //! - Problem - problem configuration, see PDE/Transport/Problem.h
53 : : //! \note The default physics is DGAdvection, set in
54 : : //! inciter::deck::check_transport()
55 : : template< class Physics, class Problem >
56 [ - - ][ + - ]: 160 : class Transport {
[ + - ][ + - ]
[ - - ]
57 : :
58 : : private:
59 : : using eq = tag::transport;
60 : :
61 : : public:
62 : : //! Constructor
63 : 160 : explicit Transport() :
64 : : m_physics( Physics() ),
65 : : m_problem( Problem() ),
66 [ + - ]: 160 : m_ncomp( g_inputdeck.get< tag::ncomp >() )
67 : : {
68 : : // associate boundary condition configurations with state functions, the
69 : : // order in which the state functions listed matters, see ctr::bc::Keys
70 [ + - ][ + - ]: 1120 : brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
[ + + ][ + - ]
[ - - ][ - - ]
[ - - ]
71 : : { dirichlet
72 : : , invalidBC // Symmetry BC not implemented
73 : : , inlet
74 : : , outlet
75 : : , invalidBC // Characteristic BC not implemented
76 : : , extrapolate } ) );
77 [ - - ]: 0 : m_problem.errchk( m_ncomp );
78 : 160 : }
79 : :
80 : : //! Find the number of primitive quantities required for this PDE system
81 : : //! \return The number of primitive quantities required to be stored for
82 : : //! this PDE system
83 : : std::size_t nprim() const
84 : : {
85 : : // transport does not need/store any primitive quantities currently
86 : : return 0;
87 : : }
88 : :
89 : : //! Find the number of materials set up for this PDE system
90 : : //! \return The number of materials set up for this PDE system
91 : : std::size_t nmat() const
92 : : {
93 : : return m_ncomp;
94 : : }
95 : :
96 : : //! Assign number of DOFs per equation in the PDE system
97 : : //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
98 : : //! equation
99 : : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
100 : : {
101 : : // all equation-dofs initialized to ndofs
102 [ - - ][ + + ]: 1150 : for (std::size_t i=0; i<m_ncomp; ++i) {
[ + + ][ + + ]
[ - - ]
103 : 575 : numEqDof.push_back(g_inputdeck.get< tag::ndof >());
104 : : }
105 : : }
106 : :
107 : : //! Determine elements that lie inside the user-defined IC box
108 : : void IcBoxElems( const tk::Fields&,
109 : : std::size_t,
110 : : std::vector< std::unordered_set< std::size_t > >& ) const
111 : : {}
112 : :
113 : : //! Initalize the transport equations for DG
114 : : //! \param[in] L Element mass matrix
115 : : //! \param[in] inpoel Element-node connectivity
116 : : //! \param[in] coord Array of nodal coordinates
117 : : //! \param[in,out] unk Array of unknowns
118 : : //! \param[in] t Physical time
119 : : //! \param[in] nielem Number of internal elements
120 : : void
121 [ + - ]: 699 : initialize(
122 : : const tk::Fields& L,
123 : : const std::vector< std::size_t >& inpoel,
124 : : const tk::UnsMesh::Coords& coord,
125 : : const std::vector< std::unordered_set< std::size_t > >& /*inbox*/,
126 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&,
127 : : tk::Fields& unk,
128 : : tk::real t,
129 : : const std::size_t nielem ) const
130 : : {
131 [ + - ]: 699 : tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
132 : : Problem::initialize, unk, t, nielem );
133 : 699 : }
134 : :
135 : : //! Compute the left hand side mass matrix
136 : : //! \param[in] geoElem Element geometry array
137 : : //! \param[in,out] l Block diagonal mass matrix
138 : : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
139 : 699 : const auto ndof = g_inputdeck.get< tag::ndof >();
140 : 699 : tk::mass( m_ncomp, ndof, geoElem, l );
141 : : }
142 : :
143 : : //! Update the interface cells to first order dofs
144 : : //! \details This function resets the high-order terms in interface cells,
145 : : //! and is currently not used in transport.
146 : : void updateInterfaceCells( tk::Fields&,
147 : : std::size_t,
148 : : std::vector< std::size_t >& ) const {}
149 : :
150 : : //! Update the primitives for this PDE system
151 : : //! \details This function computes and stores the dofs for primitive
152 : : //! quantities, which are currently unused for transport.
153 : : void updatePrimitives( const tk::Fields&,
154 : : const tk::Fields&,
155 : : const tk::Fields&,
156 : : tk::Fields&,
157 : : std::size_t ) const {}
158 : :
159 : : //! Clean up the state of trace materials for this PDE system
160 : : //! \details This function cleans up the state of materials present in trace
161 : : //! quantities in each cell. This is currently unused for transport.
162 : : void cleanTraceMaterial( tk::real,
163 : : const tk::Fields&,
164 : : tk::Fields&,
165 : : tk::Fields&,
166 : : std::size_t ) const {}
167 : :
168 : : //! Reconstruct second-order solution from first-order
169 : : //! \param[in] t Physical time
170 : : //! \param[in] geoFace Face geometry array
171 : : //! \param[in] geoElem Element geometry array
172 : : //! \param[in] fd Face connectivity and boundary conditions object
173 : : //! \param[in] esup Elements-surrounding-nodes connectivity
174 : : //! \param[in] inpoel Element-node connectivity
175 : : //! \param[in] coord Array of nodal coordinates
176 : : //! \param[in,out] U Solution vector at recent time step
177 : : //! \param[in,out] P Primitive vector at recent time step
178 : 25875 : void reconstruct( tk::real t,
179 : : const tk::Fields& geoFace,
180 : : const tk::Fields& geoElem,
181 : : const inciter::FaceData& fd,
182 : : const std::map< std::size_t, std::vector< std::size_t > >&
183 : : esup,
184 : : const std::vector< std::size_t >& inpoel,
185 : : const tk::UnsMesh::Coords& coord,
186 : : tk::Fields& U,
187 : : tk::Fields& P ) const
188 : : {
189 : 25875 : const auto rdof = g_inputdeck.get< tag::rdof >();
190 : :
191 : : // do reconstruction only if P0P1
192 [ + + ][ + + ]: 25875 : if (rdof == 4 && g_inputdeck.get< tag::ndof >() == 1) {
193 : 6750 : const auto nelem = fd.Esuel().size()/4;
194 : 6750 : const auto intsharp = g_inputdeck.get< tag::transport,
195 : : tag::intsharp >();
196 : :
197 : : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
198 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
199 : : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
200 : : "Mismatch in inpofa size" );
201 : :
202 : : // allocate and initialize matrix and vector for reconstruction
203 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
204 : 6750 : lhs_ls( nelem, {{ {{0.0, 0.0, 0.0}},
205 : : {{0.0, 0.0, 0.0}},
206 : : {{0.0, 0.0, 0.0}} }} );
207 : : // specify how many variables need to be reconstructed
208 : : std::vector< std::size_t > vars;
209 [ + + ][ + - ]: 13500 : for (std::size_t c=0; c<m_ncomp; ++c) vars.push_back(c);
210 : :
211 : : std::vector< std::vector< std::array< tk::real, 3 > > >
212 [ + - ][ + - ]: 13500 : rhs_ls( nelem, std::vector< std::array< tk::real, 3 > >
[ - - ]
213 : : ( m_ncomp,
214 : : {{ 0.0, 0.0, 0.0 }} ) );
215 : :
216 : : // reconstruct x,y,z-derivatives of unknowns
217 : : // 0. get lhs matrix, which is only geometry dependent
218 [ + - ]: 6750 : tk::lhsLeastSq_P0P1(fd, geoElem, geoFace, lhs_ls);
219 : :
220 : : // 1. internal face contributions
221 [ + - ]: 6750 : tk::intLeastSq_P0P1( rdof, fd, geoElem, U, rhs_ls, vars );
222 : :
223 : : // 2. boundary face contributions
224 [ + + ]: 47250 : for (const auto& b : m_bc)
225 : 40500 : tk::bndLeastSqConservedVar_P0P1( m_ncomp,
226 [ + - ]: 40500 : m_mat_blk, rdof, b.first, fd, geoFace, geoElem, t, b.second,
227 : : P, U, rhs_ls, vars );
228 : :
229 : : // 3. solve 3x3 least-squares system
230 [ + - ]: 6750 : tk::solveLeastSq_P0P1( rdof, lhs_ls, rhs_ls, U, vars );
231 : :
232 [ + + ]: 1646100 : for (std::size_t e=0; e<nelem; ++e)
233 : : {
234 [ + - ]: 1639350 : std::vector< std::size_t > matInt(m_ncomp, 0);
235 [ + - ][ - - ]: 1639350 : std::vector< tk::real > alAvg(m_ncomp, 0.0);
236 [ + + ]: 3278700 : for (std::size_t k=0; k<m_ncomp; ++k)
237 : 1639350 : alAvg[k] = U(e, k*rdof);
238 [ + - ]: 1639350 : auto intInd = interfaceIndicator(m_ncomp, alAvg, matInt);
239 [ - + ]: 1639350 : if ((intsharp > 0) && intInd)
240 : : {
241 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
242 : : // using nodal-stencils, for a good interface-normal estimate
243 [ - - ]: 0 : tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem,
244 : : U, vars );
245 : : }
246 : : }
247 : :
248 : : // 4. transform reconstructed derivatives to Dubiner dofs
249 [ + - ]: 6750 : tk::transform_P0P1( rdof, nelem, inpoel, coord, U, vars );
250 : : }
251 : 25875 : }
252 : :
253 : : //! Limit second-order solution
254 : : //! \param[in] t Physical time
255 : : //! \param[in] geoFace Face geometry array
256 : : //! \param[in] fd Face connectivity and boundary conditions object
257 : : //! \param[in] esup Elements surrounding points
258 : : //! \param[in] inpoel Element-node connectivity
259 : : //! \param[in] coord Array of nodal coordinates
260 : : //! \param[in] ndofel Vector of local number of degrees of freedome
261 : : // //! \param[in] gid Local->global node id map
262 : : // //! \param[in] bid Local chare-boundary node ids (value) associated to
263 : : // //! global node ids (key)
264 : : // //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
265 : : // //! variables
266 : : //! \param[in,out] U Solution vector at recent time step
267 : 25875 : void limit( [[maybe_unused]] tk::real t,
268 : : [[maybe_unused]] const tk::Fields& geoFace,
269 : : const tk::Fields&,
270 : : const inciter::FaceData& fd,
271 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
272 : : const std::vector< std::size_t >& inpoel,
273 : : const tk::UnsMesh::Coords& coord,
274 : : const std::vector< std::size_t >& ndofel,
275 : : const std::vector< std::size_t >&,
276 : : const std::unordered_map< std::size_t, std::size_t >&,
277 : : const std::vector< std::vector<tk::real> >&,
278 : : const std::vector< std::vector<tk::real> >&,
279 : : const std::vector< std::vector<tk::real> >&,
280 : : tk::Fields& U,
281 : : tk::Fields&,
282 : : std::vector< std::size_t >& ) const
283 : : {
284 : 25875 : const auto limiter = g_inputdeck.get< tag::limiter >();
285 : :
286 [ + + ]: 25875 : if (limiter == ctr::LimiterType::WENOP1)
287 : 6000 : WENO_P1( fd.Esuel(), U );
288 [ + + ]: 19875 : else if (limiter == ctr::LimiterType::SUPERBEEP1)
289 : 19500 : Superbee_P1( fd.Esuel(), inpoel, ndofel, coord, U );
290 [ - + ]: 375 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1)
291 : 0 : VertexBasedTransport_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
292 : : coord, U );
293 : 25875 : }
294 : :
295 : : //! Update the conservative variable solution for this PDE system
296 : : //! \details This function computes the updated dofs for conservative
297 : : //! quantities based on the limited solution and is currently not used in
298 : : //! transport.
299 : : void CPL( const tk::Fields&,
300 : : const tk::Fields&,
301 : : const std::vector< std::size_t >&,
302 : : const tk::UnsMesh::Coords&,
303 : : tk::Fields&,
304 : : std::size_t ) const {}
305 : :
306 : : //! Return cell-average deformation gradient tensor (no-op for transport)
307 : : //! \details This function is a no-op in transport.
308 : : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
309 : : const tk::Fields&,
310 : : std::size_t ) const
311 : : {
312 : : return {};
313 : : }
314 : :
315 : : //! Compute right hand side
316 : : //! \param[in] t Physical time
317 : : //! \param[in] geoFace Face geometry array
318 : : //! \param[in] geoElem Element geometry array
319 : : //! \param[in] fd Face connectivity and boundary conditions object
320 : : //! \param[in] inpoel Element-node connectivity
321 : : //! \param[in] coord Array of nodal coordinates
322 : : //! \param[in] U Solution vector at recent time step
323 : : //! \param[in] P Primitive vector at recent time step
324 : : //! \param[in] ndofel Vector of local number of degrees of freedom
325 : : //! \param[in] dt Delta time
326 : : //! \param[in,out] R Right-hand side vector computed
327 : 51030 : void rhs( tk::real t,
328 : : const tk::Fields& geoFace,
329 : : const tk::Fields& geoElem,
330 : : const inciter::FaceData& fd,
331 : : const std::vector< std::size_t >& inpoel,
332 : : const std::vector< std::unordered_set< std::size_t > >&,
333 : : const tk::UnsMesh::Coords& coord,
334 : : const tk::Fields& U,
335 : : const tk::Fields& P,
336 : : const std::vector< std::size_t >& ndofel,
337 : : const tk::real dt,
338 : : tk::Fields& R ) const
339 : : {
340 : 51030 : const auto ndof = g_inputdeck.get< tag::ndof >();
341 : 51030 : const auto rdof = g_inputdeck.get< tag::rdof >();
342 : 51030 : const auto intsharp = g_inputdeck.get< tag::transport,
343 : : tag::intsharp >();
344 : :
345 : : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
346 : : "vector and primitive vector at recent time step incorrect" );
347 : : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
348 : : "vector and right-hand side at recent time step incorrect" );
349 : : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
350 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
351 : : Assert( P.nprop() == 0, "Number of components in primitive "
352 : : "vector must equal "+ std::to_string(0) );
353 : : Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
354 : : "side vector must equal "+ std::to_string(ndof*m_ncomp) );
355 : : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
356 : : "Mismatch in inpofa size" );
357 : :
358 : : // set rhs to zero
359 : : R.fill(0.0);
360 : :
361 : : // empty vector for non-conservative terms. This vector is unused for
362 : : // linear transport since, there are no non-conservative terms in the
363 : : // system of PDEs.
364 : 51030 : std::vector< std::vector < tk::real > > riemannDeriv;
365 : :
366 : 51030 : std::vector< std::vector< tk::real > > vriem;
367 : 51030 : std::vector< std::vector< tk::real > > riemannLoc;
368 : :
369 : : // compute internal surface flux integrals
370 [ + - ][ + - ]: 51030 : std::vector< std::size_t > solidx(1, 0);
371 [ + - ][ + - ]: 102060 : tk::surfInt( m_ncomp, m_mat_blk, t, ndof, rdof,
[ - - ]
372 : : inpoel, solidx, coord, fd, geoFace, geoElem, Upwind::flux,
373 : : Problem::prescribedVelocity, U, P, ndofel, dt, R, vriem,
374 : : riemannLoc, riemannDeriv, intsharp );
375 : :
376 [ + + ]: 51030 : if(ndof > 1)
377 : : // compute volume integrals
378 [ + - ][ + - ]: 57375 : tk::volInt( m_ncomp, t, m_mat_blk, ndof, rdof,
[ - - ]
379 [ + - ]: 19125 : fd.Esuel().size()/4, inpoel, coord, geoElem, flux,
380 : : Problem::prescribedVelocity, U, P, ndofel, R, intsharp );
381 : :
382 : : // compute boundary surface flux integrals
383 [ + + ]: 357210 : for (const auto& b : m_bc)
384 [ + - ][ - - ]: 612360 : tk::bndSurfInt( m_ncomp, m_mat_blk, ndof, rdof,
[ - - ]
385 [ + - ]: 306180 : b.first, fd, geoFace, geoElem, inpoel, coord, t, Upwind::flux,
386 [ + - ]: 306180 : Problem::prescribedVelocity, b.second, U, P, ndofel, R, vriem,
387 : : riemannLoc, riemannDeriv, intsharp );
388 : 51030 : }
389 : :
390 : : //! Evaluate the adaptive indicator and mark the ndof for each element
391 : : //! \param[in] nunk Number of unknowns
392 : : //! \param[in] coord Array of nodal coordinates
393 : : //! \param[in] inpoel Element-node connectivity
394 : : //! \param[in] fd Face connectivity and boundary conditions object
395 : : //! \param[in] unk Array of unknowns
396 : : //! \param[in] prim Array of primitive quantities
397 : : //! \param[in] indicator p-refinement indicator type
398 : : //! \param[in] ndof Number of degrees of freedom in the solution
399 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
400 : : //! \param[in] tolref Tolerance for p-refinement
401 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
402 [ - - ]: 0 : void eval_ndof( std::size_t nunk,
403 : : [[maybe_unused]] const tk::UnsMesh::Coords& coord,
404 : : [[maybe_unused]] const std::vector< std::size_t >& inpoel,
405 : : const inciter::FaceData& fd,
406 : : const tk::Fields& unk,
407 : : const tk::Fields& prim,
408 : : inciter::ctr::PrefIndicatorType indicator,
409 : : std::size_t ndof,
410 : : std::size_t ndofmax,
411 : : tk::real tolref,
412 : : std::vector< std::size_t >& ndofel ) const
413 : : {
414 : : const auto& esuel = fd.Esuel();
415 : :
416 [ - - ]: 0 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
417 : 0 : spectral_decay( 1, nunk, esuel, unk, prim, ndof, ndofmax, tolref,
418 : : ndofel );
419 : : else
420 [ - - ][ - - ]: 0 : Throw( "No such adaptive indicator type" );
[ - - ][ - - ]
[ - - ][ - - ]
421 : 0 : }
422 : :
423 : : //! Compute the minimum time step size
424 : : // //! \param[in] U Solution vector at recent time step
425 : : // //! \param[in] coord Mesh node coordinates
426 : : // //! \param[in] inpoel Mesh element connectivity
427 : : //! \return Minimum time step size
428 : : tk::real dt( const std::array< std::vector< tk::real >, 3 >& /*coord*/,
429 : : const std::vector< std::size_t >& /*inpoel*/,
430 : : const inciter::FaceData& /*fd*/,
431 : : const tk::Fields& /*geoFace*/,
432 : : const tk::Fields& /*geoElem*/,
433 : : const std::vector< std::size_t >& /*ndofel*/,
434 : : const tk::Fields& /*U*/,
435 : : const tk::Fields&,
436 : : const std::size_t /*nielem*/ ) const
437 : : {
438 : : tk::real mindt = std::numeric_limits< tk::real >::max();
439 : : return mindt;
440 : : }
441 : :
442 : : //! Return a map that associates user-specified strings to functions
443 : : //! \return Map that associates user-specified strings to functions that
444 : : //! compute relevant quantities to be output to file
445 [ + - ]: 2256 : std::map< std::string, tk::GetVarFn > OutVarFn() const {
446 : : std::map< std::string, tk::GetVarFn > OutFnMap;
447 [ + - ][ + - ]: 2256 : OutFnMap["material_indicator"] = transport::matIndicatorOutVar;
448 : :
449 : 2256 : return OutFnMap;
450 : : }
451 : :
452 : : //! Return analytic field names to be output to file
453 : : //! \return Vector of strings labelling analytic fields output in file
454 : 1052 : std::vector< std::string > analyticFieldNames() const {
455 : : std::vector< std::string > n;
456 : 1052 : auto depvar = g_inputdeck.get< tag::depvar >()[0];
457 [ + + ]: 2104 : for (ncomp_t c=0; c<m_ncomp; ++c)
458 [ + - ][ + - ]: 2104 : n.push_back( depvar + std::to_string(c) + "_analytic" );
[ - + ][ - + ]
[ - - ][ - - ]
459 : 1052 : return n;
460 : : }
461 : :
462 : : //! Return surface field output going to file
463 : : std::vector< std::vector< tk::real > >
464 : : surfOutput( const std::map< int, std::vector< std::size_t > >&,
465 : : tk::Fields& ) const
466 : : {
467 : : std::vector< std::vector< tk::real > > s; // punt for now
468 : : return s;
469 : : }
470 : :
471 : : //! Return time history field names to be output to file
472 : : //! \return Vector of strings labelling time history fields output in file
473 : : std::vector< std::string > histNames() const {
474 : : std::vector< std::string > s; // punt for now
475 : : return s;
476 : : }
477 : :
478 : : //! Return names of integral variables to be output to diagnostics file
479 : : //! \return Vector of strings labelling integral variables output
480 [ + - ]: 41 : std::vector< std::string > names() const {
481 : : std::vector< std::string > n;
482 : : const auto& depvar =
483 : : g_inputdeck.get< tag::depvar >().at(0);
484 : : // construct the name of the numerical solution for all components
485 [ + + ]: 82 : for (ncomp_t c=0; c<m_ncomp; ++c)
486 [ + - ][ - + ]: 82 : n.push_back( depvar + std::to_string(c) );
[ - - ]
487 : 41 : return n;
488 : : }
489 : :
490 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
491 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
492 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
493 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
494 : : //! \param[in] t Physical time at which to evaluate the analytic solution
495 : : //! \return Vector of analytic solution at given spatial location and time
496 : : std::vector< tk::real >
497 : : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
498 : 872444 : { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi,
499 : : zi, t ); }
500 : :
501 : : //! Return analytic solution for conserved variables
502 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
503 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
504 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
505 : : //! \param[in] t Physical time at which to evaluate the analytic solution
506 : : //! \return Vector of analytic solution at given location and time
507 : : std::vector< tk::real >
508 : : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
509 : 2495759 : { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
510 : :
511 : : //! Return time history field output evaluated at time history points
512 : : //! \param[in] h History point data
513 : : std::vector< std::vector< tk::real > >
514 : : histOutput( const std::vector< HistData >& h,
515 : : const std::vector< std::size_t >&,
516 : : const tk::UnsMesh::Coords&,
517 : : const tk::Fields&,
518 : : const tk::Fields& ) const
519 : : {
520 : 0 : std::vector< std::vector< tk::real > > Up(h.size()); //punt for now
521 : : return Up;
522 : : }
523 : :
524 : : //! Return cell-averaged total component mass per unit volume for an element
525 : : //! \param[in] e Element id for which total energy is required
526 : : //! \param[in] unk Vector of conserved quantities
527 : : //! \return Cell-averaged total component mass per unit volume for given
528 : : //! element. Since transport does not have an associated total energy,
529 : : //! return total mass.
530 : : tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
531 : : {
532 : 1599581 : const auto rdof = g_inputdeck.get< tag::rdof >();
533 : :
534 : : tk::real sp_m(0.0);
535 [ - - ][ + + ]: 3199162 : for (std::size_t c=0; c<m_ncomp; ++c) {
[ + + ][ + + ]
[ - - ]
536 : 1599581 : sp_m += unk(e,c*rdof);
537 : : }
538 : : return sp_m;
539 : : }
540 : :
541 : : private:
542 : : const Physics m_physics; //!< Physics policy
543 : : const Problem m_problem; //!< Problem policy
544 : : const ncomp_t m_ncomp; //!< Number of components in this PDE
545 : : //! BC configuration
546 : : BCStateFn m_bc;
547 : : //! \brief EOS material block - This PDE does not require an EOS block,
548 : : //! thus this variable has not been intialized.
549 : : std::vector< EOS > m_mat_blk;
550 : :
551 : : //! Evaluate physical flux function for this PDE system
552 : : //! \param[in] ncomp Number of scalar components in this PDE system
553 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
554 : : //! evaluate the flux
555 : : //! \param[in] v Prescribed velocity evaluated at the Gauss point at which
556 : : //! to evaluate the flux
557 : : //! \return Flux vectors for all components in this PDE system
558 : : //! \note The function signature must follow tk::FluxFn
559 : : static tk::FluxFn::result_type
560 : 19672200 : flux( ncomp_t ncomp,
561 : : const std::vector< EOS >&,
562 : : const std::vector< tk::real >& ugp,
563 : : const std::vector< std::array< tk::real, 3 > >& v )
564 : :
565 : : {
566 : : Assert( ugp.size() == ncomp, "Size mismatch" );
567 : : Assert( v.size() == ncomp, "Size mismatch" );
568 : :
569 : 19672200 : std::vector< std::array< tk::real, 3 > > fl( ugp.size() );
570 : :
571 [ + + ]: 39344400 : for (ncomp_t c=0; c<ncomp; ++c)
572 : 19672200 : fl[c] = {{ v[c][0] * ugp[c], v[c][1] * ugp[c], v[c][2] * ugp[c] }};
573 : :
574 : 19672200 : return fl;
575 : : }
576 : :
577 : : //! \brief Boundary state function providing the left and right state of a
578 : : //! face at extrapolation boundaries
579 : : //! \param[in] ul Left (domain-internal) state
580 : : //! \return Left and right states for all scalar components in this PDE
581 : : //! system
582 : : //! \note The function signature must follow tk::StateFn
583 : : static tk::StateFn::result_type
584 : 10704420 : extrapolate( ncomp_t, const std::vector< EOS >&,
585 : : const std::vector< tk::real >& ul, tk::real, tk::real,
586 : : tk::real, tk::real, const std::array< tk::real, 3 >& )
587 : : {
588 : 10704420 : return {{ ul, ul }};
589 : : }
590 : :
591 : : //! \brief Boundary state function providing the left and right state of a
592 : : //! face at extrapolation boundaries
593 : : //! \param[in] ul Left (domain-internal) state
594 : : //! \return Left and right states for all scalar components in this PDE
595 : : //! system
596 : : //! \note The function signature must follow tk::StateFn
597 : : static tk::StateFn::result_type
598 : 307200 : inlet( ncomp_t, const std::vector< EOS >&,
599 : : const std::vector< tk::real >& ul, tk::real, tk::real, tk::real,
600 : : tk::real, const std::array< tk::real, 3 >& )
601 : : {
602 : 307200 : auto ur = ul;
603 : : std::fill( begin(ur), end(ur), 0.0 );
604 [ + - ]: 307200 : return {{ ul, std::move(ur) }};
605 : : }
606 : :
607 : : //! \brief Boundary state function providing the left and right state of a
608 : : //! face at outlet boundaries
609 : : //! \param[in] ul Left (domain-internal) state
610 : : //! \return Left and right states for all scalar components in this PDE
611 : : //! system
612 : : //! \note The function signature must follow tk::StateFn
613 : : static tk::StateFn::result_type
614 : 859200 : outlet( ncomp_t, const std::vector< EOS >&,
615 : : const std::vector< tk::real >& ul, tk::real, tk::real, tk::real,
616 : : tk::real, const std::array< tk::real, 3 >& )
617 : : {
618 : 859200 : return {{ ul, ul }};
619 : : }
620 : :
621 : : //! \brief Boundary state function providing the left and right state of a
622 : : //! face at Dirichlet boundaries
623 : : //! \param[in] ncomp Number of scalar components in this PDE system
624 : : //! \param[in] ul Left (domain-internal) state
625 : : //! \param[in] x X-coordinate at which to compute the states
626 : : //! \param[in] y Y-coordinate at which to compute the states
627 : : //! \param[in] z Z-coordinate at which to compute the states
628 : : //! \param[in] t Physical time
629 : : //! \return Left and right states for all scalar components in this PDE
630 : : //! system
631 : : //! \note The function signature must follow tk::StateFn
632 : : static tk::StateFn::result_type
633 : 972540 : dirichlet( ncomp_t ncomp,
634 : : const std::vector< EOS >& mat_blk,
635 : : const std::vector< tk::real >& ul, tk::real x, tk::real y,
636 : : tk::real z, tk::real t, const std::array< tk::real, 3 >& )
637 : : {
638 : 972540 : return {{ ul, Problem::initialize( ncomp, mat_blk, x, y, z, t ) }};
639 : : }
640 : : };
641 : :
642 : : } // dg::
643 : : } // inciter::
644 : :
645 : : #endif // DGTransport_h
|
