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 [ - - ][ + - ]: 148 : class Transport {
[ + - ][ + - ]
[ - - ]
57 : :
58 : : private:
59 : : using eq = tag::transport;
60 : :
61 : : public:
62 : : //! Constructor
63 : 148 : explicit Transport() :
64 : : m_physics( Physics() ),
65 : : m_problem( Problem() ),
66 [ + - ]: 148 : 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 [ + - ][ + - ]: 2220 : brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
[ + - ][ + + ]
[ + - ][ + + ]
[ + - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
71 : : // BC State functions
72 : : { dirichlet
73 : : , invalidBC // Symmetry BC not implemented
74 : : , inlet
75 : : , outlet
76 : : , invalidBC // Characteristic BC not implemented
77 : : , extrapolate
78 : : , invalidBC }, // No slip wall BC not implemented
79 : : // BC Gradient functions
80 : : { noOpGrad
81 : : , noOpGrad
82 : : , noOpGrad
83 : : , noOpGrad
84 : : , noOpGrad
85 : : , noOpGrad
86 : : , noOpGrad }
87 : : ) );
88 [ - - ]: 0 : m_problem.errchk( m_ncomp );
89 : 148 : }
90 : :
91 : : //! Find the number of primitive quantities required for this PDE system
92 : : //! \return The number of primitive quantities required to be stored for
93 : : //! this PDE system
94 : : std::size_t nprim() const
95 : : {
96 : : // transport does not need/store any primitive quantities currently
97 : : return 0;
98 : : }
99 : :
100 : : //! Find the number of materials set up for this PDE system
101 : : //! \return The number of materials set up for this PDE system
102 : : std::size_t nmat() const
103 : : {
104 : : return m_ncomp;
105 : : }
106 : :
107 : : //! Assign number of DOFs per equation in the PDE system
108 : : //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
109 : : //! equation
110 : : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
111 : : {
112 : : // all equation-dofs initialized to ndofs
113 [ - - ][ + + ]: 1060 : for (std::size_t i=0; i<m_ncomp; ++i) {
[ + + ][ + + ]
[ - - ]
114 : 530 : numEqDof.push_back(g_inputdeck.get< tag::ndof >());
115 : : }
116 : : }
117 : :
118 : : //! Find how 'stiff equations', which we currently
119 : : //! have none for Transport
120 : : //! \return number of stiff equations
121 : : std::size_t nstiffeq() const
122 : : { return 0; }
123 : :
124 : : //! Find how 'nonstiff equations', which we currently
125 : : //! don't use for Transport
126 : : //! \return number of non-stiff equations
127 : : std::size_t nnonstiffeq() const
128 : : { return 0; }
129 : :
130 : : //! Locate the stiff equations. Unused for transport.
131 : : //! \param[out] stiffEqIdx list
132 : : void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const
133 : : {
134 : 0 : stiffEqIdx.resize(0);
135 : : }
136 : :
137 : : //! Locate the nonstiff equations. Unused for transport.
138 : : //! \param[out] nonStiffEqIdx list
139 : : void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const
140 : : {
141 : 0 : nonStiffEqIdx.resize(0);
142 : : }
143 : :
144 : : //! Determine elements that lie inside the user-defined IC box
145 : : void IcBoxElems( const tk::Fields&,
146 : : std::size_t,
147 : : std::vector< std::unordered_set< std::size_t > >& ) const
148 : : {}
149 : :
150 : : //! Initalize the transport equations for DG
151 : : //! \param[in] L Element mass matrix
152 : : //! \param[in] inpoel Element-node connectivity
153 : : //! \param[in] coord Array of nodal coordinates
154 : : //! \param[in,out] unk Array of unknowns
155 : : //! \param[in] t Physical time
156 : : //! \param[in] nielem Number of internal elements
157 : : void
158 [ + - ]: 654 : initialize(
159 : : const tk::Fields& L,
160 : : const std::vector< std::size_t >& inpoel,
161 : : const tk::UnsMesh::Coords& coord,
162 : : const std::vector< std::unordered_set< std::size_t > >& /*inbox*/,
163 : : const std::unordered_map< std::size_t, std::set< std::size_t > >&,
164 : : tk::Fields& unk,
165 : : tk::real t,
166 : : const std::size_t nielem ) const
167 : : {
168 [ + - ]: 654 : tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
169 : : Problem::initialize, unk, t, nielem );
170 : 654 : }
171 : :
172 : : //! Compute density constraint for a given material
173 : : // //! \param[in] nelem Number of elements
174 : : // //! \param[in] unk Array of unknowns
175 : : //! \param[out] densityConstr Density Constraint: rho/(rho0*det(g))
176 : : void computeDensityConstr( std::size_t /*nelem*/,
177 : : tk::Fields& /*unk*/,
178 : : std::vector< tk::real >& densityConstr) const
179 : : {
180 : 993 : densityConstr.resize(0);
181 : : }
182 : :
183 : : //! Compute the left hand side mass matrix
184 : : //! \param[in] geoElem Element geometry array
185 : : //! \param[in,out] l Block diagonal mass matrix
186 : : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
187 : 654 : const auto ndof = g_inputdeck.get< tag::ndof >();
188 : 654 : tk::mass( m_ncomp, ndof, geoElem, l );
189 : : }
190 : :
191 : : //! Update the interface cells to first order dofs
192 : : //! \details This function resets the high-order terms in interface cells,
193 : : //! and is currently not used in transport.
194 : : void updateInterfaceCells( tk::Fields&,
195 : : std::size_t,
196 : : std::vector< std::size_t >&,
197 : : std::vector< std::size_t >& ) const {}
198 : :
199 : : //! Update the primitives for this PDE system
200 : : //! \details This function computes and stores the dofs for primitive
201 : : //! quantities, which are currently unused for transport.
202 : : void updatePrimitives( const tk::Fields&,
203 : : const tk::Fields&,
204 : : const tk::Fields&,
205 : : tk::Fields&,
206 : : std::size_t,
207 : : std::vector< std::size_t >& ) const {}
208 : :
209 : : //! Clean up the state of trace materials for this PDE system
210 : : //! \details This function cleans up the state of materials present in trace
211 : : //! quantities in each cell. This is currently unused for transport.
212 : : void cleanTraceMaterial( tk::real,
213 : : const tk::Fields&,
214 : : tk::Fields&,
215 : : tk::Fields&,
216 : : std::size_t ) const {}
217 : :
218 : : //! Reconstruct second-order solution from first-order
219 : : // //! \param[in] t Physical time
220 : : // //! \param[in] geoFace Face geometry array
221 : : // //! \param[in] geoElem Element geometry array
222 : : // //! \param[in] fd Face connectivity and boundary conditions object
223 : : // //! \param[in] esup Elements-surrounding-nodes connectivity
224 : : // //! \param[in] inpoel Element-node connectivity
225 : : // //! \param[in] coord Array of nodal coordinates
226 : : // //! \param[in,out] U Solution vector at recent time step
227 : : // //! \param[in,out] P Primitive vector at recent time step
228 : 19125 : void reconstruct( tk::real,
229 : : const tk::Fields&,
230 : : const tk::Fields&,
231 : : const inciter::FaceData&,
232 : : const std::map< std::size_t, std::vector< std::size_t > >&,
233 : : const std::vector< std::size_t >&,
234 : : const tk::UnsMesh::Coords&,
235 : : tk::Fields&,
236 : : tk::Fields&,
237 : : const bool,
238 : : const std::vector< std::size_t >& ) const
239 : : {
240 : : // do reconstruction only if P0P1
241 [ + + ]: 19125 : if (g_inputdeck.get< tag::rdof >() == 4 &&
242 [ - + ]: 18750 : g_inputdeck.get< tag::ndof >() == 1)
243 [ - - ][ - - ]: 0 : Throw("P0P1 not supported for Transport.");
[ - - ][ - - ]
[ - - ][ - - ]
244 : 19125 : }
245 : :
246 : : //! Limit second-order solution
247 : : //! \param[in] t Physical time
248 : : //! \param[in] geoFace Face geometry array
249 : : //! \param[in] fd Face connectivity and boundary conditions object
250 : : //! \param[in] esup Elements surrounding points
251 : : //! \param[in] inpoel Element-node connectivity
252 : : //! \param[in] coord Array of nodal coordinates
253 : : //! \param[in] ndofel Vector of local number of degrees of freedome
254 : : // //! \param[in] gid Local->global node id map
255 : : // //! \param[in] bid Local chare-boundary node ids (value) associated to
256 : : // //! global node ids (key)
257 : : // //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
258 : : // //! variables
259 : : //! \param[in,out] U Solution vector at recent time step
260 : 19125 : void limit( [[maybe_unused]] tk::real t,
261 : : [[maybe_unused]] const bool pref,
262 : : [[maybe_unused]] const tk::Fields& geoFace,
263 : : const tk::Fields&,
264 : : const inciter::FaceData& fd,
265 : : const std::map< std::size_t, std::vector< std::size_t > >& esup,
266 : : const std::vector< std::size_t >& inpoel,
267 : : const tk::UnsMesh::Coords& coord,
268 : : const std::vector< std::size_t >& ndofel,
269 : : const std::vector< std::size_t >&,
270 : : const std::unordered_map< std::size_t, std::size_t >&,
271 : : const std::vector< std::vector<tk::real> >&,
272 : : const std::vector< std::vector<tk::real> >&,
273 : : const std::vector< std::vector<tk::real> >&,
274 : : tk::Fields& U,
275 : : tk::Fields&,
276 : : std::vector< std::size_t >& ) const
277 : : {
278 : 19125 : const auto limiter = g_inputdeck.get< tag::limiter >();
279 : :
280 [ + + ]: 19125 : if (limiter == ctr::LimiterType::WENOP1)
281 : 6000 : WENO_P1( fd.Esuel(), U );
282 [ + + ]: 13125 : else if (limiter == ctr::LimiterType::SUPERBEEP1)
283 : 12750 : Superbee_P1( fd.Esuel(), inpoel, ndofel, coord, U );
284 [ - + ]: 375 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1)
285 : 0 : VertexBasedTransport_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
286 : : coord, U );
287 : 19125 : }
288 : :
289 : : //! Update the conservative variable solution for this PDE system
290 : : //! \details This function computes the updated dofs for conservative
291 : : //! quantities based on the limited solution and is currently not used in
292 : : //! transport.
293 : : void CPL( const tk::Fields&,
294 : : const tk::Fields&,
295 : : const std::vector< std::size_t >&,
296 : : const tk::UnsMesh::Coords&,
297 : : tk::Fields&,
298 : : std::size_t ) const {}
299 : :
300 : : //! Return cell-average deformation gradient tensor (no-op for transport)
301 : : //! \details This function is a no-op in transport.
302 : : std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
303 : : const tk::Fields&,
304 : : std::size_t ) const
305 : : {
306 : : return {};
307 : : }
308 : :
309 : : //! Reset the high order solution for p-adaptive scheme
310 : : //! \details This function reset the high order coefficient for p-adaptive
311 : : //! solution polynomials and is currently not used in transport.
312 : : void resetAdapSol( const inciter::FaceData&,
313 : : tk::Fields&,
314 : : tk::Fields&,
315 : : const std::vector< std::size_t >& ) const {}
316 : :
317 : : //! Compute right hand side
318 : : //! \param[in] t Physical time
319 : : //! \param[in] pref Indicator for p-adaptive algorithm
320 : : //! \param[in] geoFace Face geometry array
321 : : //! \param[in] geoElem Element geometry array
322 : : //! \param[in] fd Face connectivity and boundary conditions object
323 : : //! \param[in] inpoel Element-node connectivity
324 : : //! \param[in] coord Array of nodal coordinates
325 : : //! \param[in] U Solution vector at recent time step
326 : : //! \param[in] P Primitive vector at recent time step
327 : : //! \param[in] ndofel Vector of local number of degrees of freedom
328 : : //! \param[in] dt Delta time
329 : : //! \param[in,out] R Right-hand side vector computed
330 : 44280 : void rhs( tk::real t,
331 : : const bool pref,
332 : : const tk::Fields& geoFace,
333 : : const tk::Fields& geoElem,
334 : : const inciter::FaceData& fd,
335 : : const std::vector< std::size_t >& inpoel,
336 : : const std::vector< std::unordered_set< std::size_t > >&,
337 : : const tk::UnsMesh::Coords& coord,
338 : : const tk::Fields& U,
339 : : const tk::Fields& P,
340 : : const std::vector< std::size_t >& ndofel,
341 : : const tk::real dt,
342 : : tk::Fields& R ) const
343 : : {
344 : 44280 : const auto ndof = g_inputdeck.get< tag::ndof >();
345 : 44280 : const auto rdof = g_inputdeck.get< tag::rdof >();
346 : 44280 : const auto intsharp = g_inputdeck.get< tag::transport,
347 : : tag::intsharp >();
348 : :
349 : : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
350 : : "vector and primitive vector at recent time step incorrect" );
351 : : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
352 : : "vector and right-hand side at recent time step incorrect" );
353 : : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
354 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
355 : : Assert( P.nprop() == 0, "Number of components in primitive "
356 : : "vector must equal "+ std::to_string(0) );
357 : : Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
358 : : "side vector must equal "+ std::to_string(ndof*m_ncomp) );
359 : : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
360 : : "Mismatch in inpofa size" );
361 : :
362 : : // set rhs to zero
363 : : R.fill(0.0);
364 : :
365 : : // empty vector for non-conservative terms. This vector is unused for
366 : : // linear transport since, there are no non-conservative terms in the
367 : : // system of PDEs.
368 : 44280 : std::vector< std::vector < tk::real > > riemannDeriv;
369 : :
370 : : // compute internal surface flux integrals
371 [ + - ][ + - ]: 44280 : std::vector< std::size_t > solidx(1, 0);
372 [ + - ][ + - ]: 88560 : tk::surfInt( pref, m_ncomp, m_mat_blk, t, ndof, rdof,
[ - - ]
373 : : inpoel, solidx, coord, fd, geoFace, geoElem, Upwind::flux,
374 : : Problem::prescribedVelocity, U, P, ndofel, dt, R,
375 : : riemannDeriv, intsharp );
376 : :
377 [ + + ]: 44280 : if(ndof > 1)
378 : : // compute volume integrals
379 [ + - ][ + - ]: 57375 : tk::volInt( m_ncomp, t, m_mat_blk, ndof, rdof,
[ - - ]
380 [ + - ]: 19125 : fd.Esuel().size()/4, inpoel, coord, geoElem, flux,
381 : : Problem::prescribedVelocity, U, P, ndofel, R, intsharp );
382 : :
383 : : // compute boundary surface flux integrals
384 [ + + ]: 354240 : for (const auto& b : m_bc)
385 [ + - ][ + - ]: 929880 : tk::bndSurfInt( pref, m_ncomp, m_mat_blk, ndof, rdof,
[ - - ][ - - ]
386 : : std::get<0>(b), fd, geoFace, geoElem, inpoel, coord, t, Upwind::flux,
387 : : Problem::prescribedVelocity, std::get<1>(b), U, P, ndofel, R,
388 : : riemannDeriv, intsharp );
389 : 44280 : }
390 : :
391 : : //! Evaluate the adaptive indicator and mark the ndof for each element
392 : : //! \param[in] nunk Number of unknowns
393 : : //! \param[in] coord Array of nodal coordinates
394 : : //! \param[in] inpoel Element-node connectivity
395 : : //! \param[in] fd Face connectivity and boundary conditions object
396 : : //! \param[in] unk Array of unknowns
397 : : //! \param[in] prim Array of primitive quantities
398 : : //! \param[in] indicator p-refinement indicator type
399 : : //! \param[in] ndof Number of degrees of freedom in the solution
400 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
401 : : //! \param[in] tolref Tolerance for p-refinement
402 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
403 [ - - ]: 0 : void eval_ndof( std::size_t nunk,
404 : : [[maybe_unused]] const tk::UnsMesh::Coords& coord,
405 : : [[maybe_unused]] const std::vector< std::size_t >& inpoel,
406 : : const inciter::FaceData& fd,
407 : : const tk::Fields& unk,
408 : : [[maybe_unused]] const tk::Fields& prim,
409 : : inciter::ctr::PrefIndicatorType indicator,
410 : : std::size_t ndof,
411 : : std::size_t ndofmax,
412 : : tk::real tolref,
413 : : std::vector< std::size_t >& ndofel ) const
414 : : {
415 : : const auto& esuel = fd.Esuel();
416 : :
417 [ - - ]: 0 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
418 : 0 : spectral_decay( 1, nunk, esuel, unk, ndof, ndofmax, tolref, 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 : : //! Compute stiff terms for a single element, not implemented here
443 : : // //! \param[in] e Element number
444 : : // //! \param[in] geoElem Element geometry array
445 : : // //! \param[in] inpoel Element-node connectivity
446 : : // //! \param[in] coord Array of nodal coordinates
447 : : // //! \param[in] U Solution vector at recent time step
448 : : // //! \param[in] P Primitive vector at recent time step
449 : : // //! \param[in] ndofel Vector of local number of degrees of freedom
450 : : // //! \param[in,out] R Right-hand side vector computed
451 : : void stiff_rhs( std::size_t /*e*/,
452 : : const tk::Fields& /*geoElem*/,
453 : : const std::vector< std::size_t >& /*inpoel*/,
454 : : const tk::UnsMesh::Coords& /*coord*/,
455 : : const tk::Fields& /*U*/,
456 : : const tk::Fields& /*P*/,
457 : : const std::vector< std::size_t >& /*ndofel*/,
458 : : tk::Fields& /*R*/ ) const
459 : : {}
460 : :
461 : : //! Return a map that associates user-specified strings to functions
462 : : //! \return Map that associates user-specified strings to functions that
463 : : //! compute relevant quantities to be output to file
464 [ + - ]: 1986 : std::map< std::string, tk::GetVarFn > OutVarFn() const {
465 : : std::map< std::string, tk::GetVarFn > OutFnMap;
466 [ + - ][ + - ]: 1986 : OutFnMap["material_indicator"] = transport::matIndicatorOutVar;
467 : :
468 : 1986 : return OutFnMap;
469 : : }
470 : :
471 : : //! Return analytic field names to be output to file
472 : : //! \return Vector of strings labelling analytic fields output in file
473 : 917 : std::vector< std::string > analyticFieldNames() const {
474 : : std::vector< std::string > n;
475 : 917 : auto depvar = g_inputdeck.get< tag::depvar >()[0];
476 [ + + ]: 1834 : for (ncomp_t c=0; c<m_ncomp; ++c)
477 [ + - ][ + - ]: 1834 : n.push_back( depvar + std::to_string(c) + "_analytic" );
[ - + ][ - + ]
[ - - ][ - - ]
478 : 917 : return n;
479 : : }
480 : :
481 : : //! Return surface field output going to file
482 : : std::vector< std::vector< tk::real > >
483 : : surfOutput( const std::map< int, std::vector< std::size_t > >&,
484 : : tk::Fields& ) const
485 : : {
486 : : std::vector< std::vector< tk::real > > s; // punt for now
487 : : return s;
488 : : }
489 : :
490 : : //! Return time history field names to be output to file
491 : : //! \return Vector of strings labelling time history fields output in file
492 : : std::vector< std::string > histNames() const {
493 : : std::vector< std::string > s; // punt for now
494 : : return s;
495 : : }
496 : :
497 : : //! Return names of integral variables to be output to diagnostics file
498 : : //! \return Vector of strings labelling integral variables output
499 [ + - ]: 38 : std::vector< std::string > names() const {
500 : : std::vector< std::string > n;
501 : : const auto& depvar =
502 : : g_inputdeck.get< tag::depvar >().at(0);
503 : : // construct the name of the numerical solution for all components
504 [ + + ]: 76 : for (ncomp_t c=0; c<m_ncomp; ++c)
505 [ + - ][ - + ]: 76 : n.push_back( depvar + std::to_string(c) );
[ - - ]
506 : 38 : return n;
507 : : }
508 : :
509 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
510 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
511 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
512 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
513 : : //! \param[in] t Physical time at which to evaluate the analytic solution
514 : : //! \return Vector of analytic solution at given spatial location and time
515 : : std::vector< tk::real >
516 : : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
517 : 839522 : { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi,
518 : : zi, t ); }
519 : :
520 : : //! Return analytic solution for conserved variables
521 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
522 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
523 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
524 : : //! \param[in] t Physical time at which to evaluate the analytic solution
525 : : //! \return Vector of analytic solution at given location and time
526 : : std::vector< tk::real >
527 : : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
528 : 2211605 : { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
529 : :
530 : : //! Return time history field output evaluated at time history points
531 : : //! \param[in] h History point data
532 : : std::vector< std::vector< tk::real > >
533 : : histOutput( const std::vector< HistData >& h,
534 : : const std::vector< std::size_t >&,
535 : : const tk::UnsMesh::Coords&,
536 : : const tk::Fields&,
537 : : const tk::Fields& ) const
538 : : {
539 : 0 : std::vector< std::vector< tk::real > > Up(h.size()); //punt for now
540 : : return Up;
541 : : }
542 : :
543 : : //! Return cell-averaged total component mass per unit volume for an element
544 : : //! \param[in] e Element id for which total energy is required
545 : : //! \param[in] unk Vector of conserved quantities
546 : : //! \return Cell-averaged total component mass per unit volume for given
547 : : //! element. Since transport does not have an associated total energy,
548 : : //! return total mass.
549 : : tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
550 : : {
551 : 1315427 : const auto rdof = g_inputdeck.get< tag::rdof >();
552 : :
553 : : tk::real sp_m(0.0);
554 [ - - ][ + + ]: 2630854 : for (std::size_t c=0; c<m_ncomp; ++c) {
[ + + ][ + + ]
[ - - ]
555 : 1315427 : sp_m += unk(e,c*rdof);
556 : : }
557 : : return sp_m;
558 : : }
559 : :
560 : : private:
561 : : const Physics m_physics; //!< Physics policy
562 : : const Problem m_problem; //!< Problem policy
563 : : const ncomp_t m_ncomp; //!< Number of components in this PDE
564 : : //! BC configuration
565 : : BCStateFn m_bc;
566 : : //! \brief EOS material block - This PDE does not require an EOS block,
567 : : //! thus this variable has not been intialized.
568 : : std::vector< EOS > m_mat_blk;
569 : :
570 : : //! Evaluate physical flux function for this PDE system
571 : : //! \param[in] ncomp Number of scalar components in this PDE system
572 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
573 : : //! evaluate the flux
574 : : //! \param[in] v Prescribed velocity evaluated at the Gauss point at which
575 : : //! to evaluate the flux
576 : : //! \return Flux vectors for all components in this PDE system
577 : : //! \note The function signature must follow tk::FluxFn
578 : : static tk::FluxFn::result_type
579 : 19672200 : flux( ncomp_t ncomp,
580 : : const std::vector< EOS >&,
581 : : const std::vector< tk::real >& ugp,
582 : : const std::vector< std::array< tk::real, 3 > >& v )
583 : :
584 : : {
585 : : Assert( ugp.size() == ncomp, "Size mismatch" );
586 : : Assert( v.size() == ncomp, "Size mismatch" );
587 : :
588 : 19672200 : std::vector< std::array< tk::real, 3 > > fl( ugp.size() );
589 : :
590 [ + + ]: 39344400 : for (ncomp_t c=0; c<ncomp; ++c)
591 : 19672200 : fl[c] = {{ v[c][0] * ugp[c], v[c][1] * ugp[c], v[c][2] * ugp[c] }};
592 : :
593 : 19672200 : return fl;
594 : : }
595 : :
596 : : //! \brief Boundary state function providing the left and right state of a
597 : : //! face at extrapolation boundaries
598 : : //! \param[in] ul Left (domain-internal) state
599 : : //! \return Left and right states for all scalar components in this PDE
600 : : //! system
601 : : //! \note The function signature must follow tk::StateFn
602 : : static tk::StateFn::result_type
603 : 8834220 : extrapolate( ncomp_t, const std::vector< EOS >&,
604 : : const std::vector< tk::real >& ul, tk::real, tk::real,
605 : : tk::real, tk::real, const std::array< tk::real, 3 >& )
606 : : {
607 : 8834220 : return {{ ul, ul }};
608 : : }
609 : :
610 : : //! \brief Boundary state function providing the left and right state of a
611 : : //! face at extrapolation boundaries
612 : : //! \param[in] ul Left (domain-internal) state
613 : : //! \return Left and right states for all scalar components in this PDE
614 : : //! system
615 : : //! \note The function signature must follow tk::StateFn
616 : : static tk::StateFn::result_type
617 : 163200 : inlet( ncomp_t, const std::vector< EOS >&,
618 : : const std::vector< tk::real >& ul, tk::real, tk::real, tk::real,
619 : : tk::real, const std::array< tk::real, 3 >& )
620 : : {
621 : 163200 : auto ur = ul;
622 : : std::fill( begin(ur), end(ur), 0.0 );
623 [ + - ]: 163200 : return {{ ul, std::move(ur) }};
624 : : }
625 : :
626 : : //! \brief Boundary state function providing the left and right state of a
627 : : //! face at outlet boundaries
628 : : //! \param[in] ul Left (domain-internal) state
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 : 715200 : outlet( ncomp_t, const std::vector< EOS >&,
634 : : const std::vector< tk::real >& ul, tk::real, tk::real, tk::real,
635 : : tk::real, const std::array< tk::real, 3 >& )
636 : : {
637 : 715200 : return {{ ul, ul }};
638 : : }
639 : :
640 : : //! \brief Boundary state function providing the left and right state of a
641 : : //! face at Dirichlet boundaries
642 : : //! \param[in] ncomp Number of scalar components in this PDE system
643 : : //! \param[in] ul Left (domain-internal) state
644 : : //! \param[in] x X-coordinate at which to compute the states
645 : : //! \param[in] y Y-coordinate at which to compute the states
646 : : //! \param[in] z Z-coordinate at which to compute the states
647 : : //! \param[in] t Physical time
648 : : //! \return Left and right states for all scalar components in this PDE
649 : : //! system
650 : : //! \note The function signature must follow tk::StateFn
651 : : static tk::StateFn::result_type
652 : 972540 : dirichlet( ncomp_t ncomp,
653 : : const std::vector< EOS >& mat_blk,
654 : : const std::vector< tk::real >& ul, tk::real x, tk::real y,
655 : : tk::real z, tk::real t, const std::array< tk::real, 3 >& )
656 : : {
657 : 972540 : return {{ ul, Problem::initialize( ncomp, mat_blk, x, y, z, t ) }};
658 : : }
659 : :
660 : : //----------------------------------------------------------------------------
661 : : // Boundary Gradient functions
662 : : //----------------------------------------------------------------------------
663 : :
664 : : //! \brief Boundary gradient function copying the left gradient to the right
665 : : //! gradient at a face
666 : : //! \param[in] dul Left (domain-internal) state
667 : : //! \return Left and right states for all scalar components in this PDE
668 : : //! system
669 : : //! \note The function signature must follow tk::StateFn.
670 : : static tk::StateFn::result_type
671 : 0 : noOpGrad( ncomp_t,
672 : : const std::vector< EOS >&,
673 : : const std::vector< tk::real >& dul,
674 : : tk::real, tk::real, tk::real, tk::real,
675 : : const std::array< tk::real, 3 >& )
676 : : {
677 : 0 : return {{ dul, dul }};
678 : : }
679 : : };
680 : :
681 : : } // dg::
682 : : } // inciter::
683 : :
684 : : #endif // DGTransport_h
|
