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