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 "Inciter/Options/BC.hpp"
27 : : #include "UnsMesh.hpp"
28 : : #include "Integrate/Basis.hpp"
29 : : #include "Integrate/Quadrature.hpp"
30 : : #include "Integrate/Initialize.hpp"
31 : : #include "Integrate/Mass.hpp"
32 : : #include "Integrate/Surface.hpp"
33 : : #include "Integrate/Boundary.hpp"
34 : : #include "Integrate/Volume.hpp"
35 : : #include "Riemann/Upwind.hpp"
36 : : #include "Reconstruction.hpp"
37 : : #include "Limiter.hpp"
38 : : #include "PrefIndicator.hpp"
39 : :
40 : : namespace inciter {
41 : :
42 : : extern ctr::InputDeck g_inputdeck;
43 : :
44 : : namespace dg {
45 : :
46 : : //! \brief Transport equation used polymorphically with tk::DGPDE
47 : : //! \details The template argument(s) specify policies and are used to configure
48 : : //! the behavior of the class. The policies are:
49 : : //! - Physics - physics configuration, see PDE/Transport/Physics.h
50 : : //! - Problem - problem configuration, see PDE/Transport/Problem.h
51 : : //! \note The default physics is DGAdvection, set in
52 : : //! inciter::deck::check_transport()
53 : : template< class Physics, class Problem >
54 : : class Transport {
55 : :
56 : : private:
57 : : using eq = tag::transport;
58 : :
59 : : public:
60 : : //! Constructor
61 : : //! \param[in] c Equation system index (among multiple systems configured)
62 : 163 : explicit Transport( ncomp_t c ) :
63 : : m_physics( Physics() ),
64 : : m_problem( Problem() ),
65 : : m_system( c ),
66 : : m_ncomp(
67 : 163 : g_inputdeck.get< tag::component >().get< eq >().at(c) ),
68 : : m_offset(
69 : 326 : g_inputdeck.get< tag::component >().offset< eq >(c) )
70 : : {
71 : : // associate boundary condition configurations with state functions, the
72 : : // order in which the state functions listed matters, see ctr::bc::Keys
73 [ + - ][ + - ]: 1141 : brigand::for_each< ctr::bc::Keys >( ConfigBC< eq >( m_system, m_bc,
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + + ][ - - ]
74 : : { dirichlet
75 : : , invalidBC // Symmetry BC not implemented
76 : : , inlet
77 : : , outlet
78 : : , invalidBC // Characteristic BC not implemented
79 : : , extrapolate } ) );
80 [ - - ]: 163 : m_problem.errchk( m_system, m_ncomp );
81 : 163 : }
82 : :
83 : : //! Find the number of primitive quantities required for this PDE system
84 : : //! \return The number of primitive quantities required to be stored for
85 : : //! this PDE system
86 : 641 : std::size_t nprim() const
87 : : {
88 : : // transport does not need/store any primitive quantities currently
89 : 641 : return 0;
90 : : }
91 : :
92 : : //! Find the number of materials set up for this PDE system
93 : : //! \return The number of materials set up for this PDE system
94 : 0 : std::size_t nmat() const
95 : : {
96 : 0 : return m_ncomp;
97 : : }
98 : :
99 : : //! Assign number of DOFs per equation in the PDE system
100 : : //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
101 : : //! equation
102 : 641 : void numEquationDofs(std::vector< std::size_t >& numEqDof) const
103 : : {
104 : : // all equation-dofs initialized to ndofs
105 [ + + ]: 1282 : for (std::size_t i=0; i<m_ncomp; ++i) {
106 : 641 : numEqDof.push_back(g_inputdeck.get< tag::discr, tag::ndof >());
107 : : }
108 : 641 : }
109 : :
110 : : //! Determine elements that lie inside the user-defined IC box
111 : 641 : void IcBoxElems( const tk::Fields&,
112 : : std::size_t,
113 : : std::vector< std::unordered_set< std::size_t > >& ) const
114 : 641 : {}
115 : :
116 : : //! Initalize the transport equations for DG
117 : : //! \param[in] L Element mass matrix
118 : : //! \param[in] inpoel Element-node connectivity
119 : : //! \param[in] coord Array of nodal coordinates
120 : : //! \param[in,out] unk Array of unknowns
121 : : //! \param[in] t Physical time
122 : : //! \param[in] nielem Number of internal elements
123 : : void
124 : 765 : initialize(
125 : : const tk::Fields& L,
126 : : const std::vector< std::size_t >& inpoel,
127 : : const tk::UnsMesh::Coords& coord,
128 : : const std::vector< std::unordered_set< std::size_t > >& /*inbox*/,
129 : : tk::Fields& unk,
130 : : tk::real t,
131 : : const std::size_t nielem ) const
132 : : {
133 [ + - ]: 765 : tk::initialize( m_system, m_ncomp, m_offset, L, inpoel, coord,
134 : : Problem::initialize, unk, t, nielem );
135 : 765 : }
136 : :
137 : : //! Compute the left hand side mass matrix
138 : : //! \param[in] geoElem Element geometry array
139 : : //! \param[in,out] l Block diagonal mass matrix
140 : 873 : void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
141 : 873 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
142 : 873 : tk::mass( m_ncomp, m_offset, ndof, geoElem, l );
143 : 873 : }
144 : :
145 : : //! Update the interface cells to first order dofs
146 : : //! \details This function resets the high-order terms in interface cells,
147 : : //! and is currently not used in transport.
148 : 53010 : void updateInterfaceCells( tk::Fields&,
149 : : std::size_t,
150 : 53010 : std::vector< std::size_t >& ) const {}
151 : :
152 : : //! Update the primitives for this PDE system
153 : : //! \details This function computes and stores the dofs for primitive
154 : : //! quantities, which are currently unused for transport.
155 : 53651 : void updatePrimitives( const tk::Fields&,
156 : : const tk::Fields&,
157 : : const tk::Fields&,
158 : : tk::Fields&,
159 : 53651 : std::size_t ) const {}
160 : :
161 : : //! Clean up the state of trace materials for this PDE system
162 : : //! \details This function cleans up the state of materials present in trace
163 : : //! quantities in each cell. This is currently unused for transport.
164 : 53010 : void cleanTraceMaterial( const tk::Fields&,
165 : : tk::Fields&,
166 : : tk::Fields&,
167 : 53010 : std::size_t ) const {}
168 : :
169 : : //! Reconstruct second-order solution from first-order
170 : : //! \param[in] t Physical time
171 : : //! \param[in] geoFace Face geometry array
172 : : //! \param[in] geoElem Element geometry array
173 : : //! \param[in] fd Face connectivity and boundary conditions object
174 : : //! \param[in] esup Elements-surrounding-nodes connectivity
175 : : //! \param[in] inpoel Element-node connectivity
176 : : //! \param[in] coord Array of nodal coordinates
177 : : //! \param[in,out] U Solution vector at recent time step
178 : : //! \param[in,out] P Primitive vector at recent time step
179 : 25875 : void reconstruct( tk::real t,
180 : : const tk::Fields& geoFace,
181 : : const tk::Fields& geoElem,
182 : : const inciter::FaceData& fd,
183 : : const std::map< std::size_t, std::vector< std::size_t > >&
184 : : esup,
185 : : const std::vector< std::size_t >& inpoel,
186 : : const tk::UnsMesh::Coords& coord,
187 : : tk::Fields& U,
188 : : tk::Fields& P ) const
189 : : {
190 : 25875 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
191 : :
192 : : // do reconstruction only if P0P1
193 [ + + ][ + + ]: 25875 : if (rdof == 4 && g_inputdeck.get< tag::discr, tag::ndof >() == 1) {
[ + + ]
194 : 6750 : const auto nelem = fd.Esuel().size()/4;
195 : 13500 : const auto intsharp = g_inputdeck.get< tag::param, tag::transport,
196 : 6750 : tag::intsharp >()[m_system];
197 : :
198 [ - + ][ - - ]: 6750 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
199 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
200 [ - + ][ - - ]: 6750 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
201 : : "Mismatch in inpofa size" );
202 : :
203 : : // allocate and initialize matrix and vector for reconstruction
204 : : std::vector< std::array< std::array< tk::real, 3 >, 3 > >
205 [ + - ]: 13500 : lhs_ls( nelem, {{ {{0.0, 0.0, 0.0}},
206 : : {{0.0, 0.0, 0.0}},
207 : : {{0.0, 0.0, 0.0}} }} );
208 : : // specify how many variables need to be reconstructed
209 : 6750 : std::array< std::size_t, 2 > varRange {{0, m_ncomp-1}};
210 : :
211 : : std::vector< std::vector< std::array< tk::real, 3 > > >
212 [ + - ]: 27000 : rhs_ls( nelem, std::vector< std::array< tk::real, 3 > >
213 [ + - ]: 6750 : ( 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( m_offset, rdof, fd, geoElem, U, rhs_ls, varRange );
222 : :
223 : : // 2. boundary face contributions
224 [ + + ]: 47250 : for (const auto& b : m_bc)
225 : 40500 : tk::bndLeastSqConservedVar_P0P1( m_system, m_ncomp, m_offset, rdof,
226 [ + - ]: 40500 : b.first, fd, geoFace, geoElem, t, b.second, P, U, rhs_ls, varRange );
227 : :
228 : : // 3. solve 3x3 least-squares system
229 [ + - ]: 6750 : tk::solveLeastSq_P0P1( m_offset, rdof, lhs_ls, rhs_ls, U, varRange );
230 : :
231 [ + + ]: 1646100 : for (std::size_t e=0; e<nelem; ++e)
232 : : {
233 [ + - ]: 3278700 : std::vector< std::size_t > matInt(m_ncomp, 0);
234 [ + - ]: 3278700 : std::vector< tk::real > alAvg(m_ncomp, 0.0);
235 [ + + ]: 3278700 : for (std::size_t k=0; k<m_ncomp; ++k)
236 [ + - ]: 1639350 : alAvg[k] = U(e, k*rdof, m_offset);
237 [ + - ]: 1639350 : auto intInd = interfaceIndicator(m_ncomp, alAvg, matInt);
238 [ - + ][ - - ]: 1639350 : if ((intsharp > 0) && intInd)
239 : : {
240 : : // Reconstruct second-order dofs of volume-fractions in Taylor space
241 : : // using nodal-stencils, for a good interface-normal estimate
242 [ - - ]: 0 : tk::recoLeastSqExtStencil( rdof, m_offset, e, esup, inpoel, geoElem,
243 : : U, varRange );
244 : : }
245 : : }
246 : :
247 : : // 4. transform reconstructed derivatives to Dubiner dofs
248 [ + - ]: 6750 : tk::transform_P0P1( m_offset, rdof, nelem, inpoel, coord, U, varRange );
249 : : }
250 : 25875 : }
251 : :
252 : : //! Limit second-order solution
253 : : //! \param[in] t Physical time
254 : : //! \param[in] geoFace Face geometry array
255 : : //! \param[in] geoElem Element 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& geoElem,
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 : : tk::Fields& U,
280 : : tk::Fields&,
281 : : std::vector< std::size_t >& ) const
282 : : {
283 : 25875 : const auto limiter = g_inputdeck.get< tag::discr, tag::limiter >();
284 : :
285 [ + + ]: 25875 : if (limiter == ctr::LimiterType::WENOP1)
286 : 6000 : WENO_P1( fd.Esuel(), m_offset, U );
287 [ + + ]: 19875 : else if (limiter == ctr::LimiterType::SUPERBEEP1)
288 : 19500 : Superbee_P1( fd.Esuel(), inpoel, ndofel, m_offset, coord, U );
289 [ - + ]: 375 : else if (limiter == ctr::LimiterType::VERTEXBASEDP1)
290 : 0 : VertexBasedTransport_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
291 : 0 : m_system, m_offset, geoElem, coord, U );
292 : 25875 : }
293 : :
294 : : //! Compute right hand side
295 : : //! \param[in] t Physical time
296 : : //! \param[in] geoFace Face geometry array
297 : : //! \param[in] geoElem Element geometry array
298 : : //! \param[in] fd Face connectivity and boundary conditions object
299 : : //! \param[in] inpoel Element-node connectivity
300 : : //! \param[in] coord Array of nodal coordinates
301 : : //! \param[in] U Solution vector at recent time step
302 : : //! \param[in] P Primitive vector at recent time step
303 : : //! \param[in] ndofel Vector of local number of degrees of freedom
304 : : //! \param[in,out] R Right-hand side vector computed
305 : 53010 : void rhs( tk::real t,
306 : : const tk::Fields& geoFace,
307 : : const tk::Fields& geoElem,
308 : : const inciter::FaceData& fd,
309 : : const std::vector< std::size_t >& inpoel,
310 : : const std::vector< std::unordered_set< std::size_t > >&,
311 : : const tk::UnsMesh::Coords& coord,
312 : : const tk::Fields& U,
313 : : const tk::Fields& P,
314 : : const std::vector< std::size_t >& ndofel,
315 : : tk::Fields& R ) const
316 : : {
317 : 53010 : const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
318 : 53010 : const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
319 : 106020 : const auto intsharp = g_inputdeck.get< tag::param, tag::transport,
320 : 53010 : tag::intsharp >()[m_system];
321 : :
322 [ - + ][ - - ]: 53010 : Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
323 : : "vector and primitive vector at recent time step incorrect" );
324 [ - + ][ - - ]: 53010 : Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
[ - - ][ - - ]
325 : : "vector and right-hand side at recent time step incorrect" );
326 [ - + ][ - - ]: 53010 : Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
[ - - ][ - - ]
[ - - ]
327 : : "vector must equal "+ std::to_string(rdof*m_ncomp) );
328 [ - + ][ - - ]: 53010 : Assert( P.nprop() == 0, "Number of components in primitive "
[ - - ][ - - ]
[ - - ]
329 : : "vector must equal "+ std::to_string(0) );
330 [ - + ][ - - ]: 53010 : Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
[ - - ][ - - ]
[ - - ]
331 : : "side vector must equal "+ std::to_string(ndof*m_ncomp) );
332 [ - + ][ - - ]: 53010 : Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
[ - - ][ - - ]
333 : : "Mismatch in inpofa size" );
334 : :
335 : : // set rhs to zero
336 : 53010 : R.fill(0.0);
337 : :
338 : : // empty vector for non-conservative terms. This vector is unused for
339 : : // linear transport since, there are no non-conservative terms in the
340 : : // system of PDEs.
341 : 106020 : std::vector< std::vector < tk::real > > riemannDeriv;
342 : :
343 : 106020 : std::vector< std::vector< tk::real > > vriem;
344 : 106020 : std::vector< std::vector< tk::real > > riemannLoc;
345 : :
346 : : // compute internal surface flux integrals
347 [ + - ][ + - ]: 53010 : tk::surfInt( m_system, m_ncomp, m_offset, t, ndof, rdof, inpoel, coord,
[ + - ]
348 : : fd, geoFace, geoElem, Upwind::flux,
349 : : Problem::prescribedVelocity, U, P, ndofel, R, vriem,
350 : : riemannLoc, riemannDeriv, intsharp );
351 : :
352 [ + + ]: 53010 : if(ndof > 1)
353 : : // compute volume integrals
354 [ + - ][ + - ]: 38250 : tk::volInt( m_system, m_ncomp, m_offset, t, ndof, rdof,
[ + - ]
355 : 19125 : fd.Esuel().size()/4, inpoel, coord, geoElem, flux,
356 : : Problem::prescribedVelocity, U, P, ndofel, R, intsharp );
357 : :
358 : : // compute boundary surface flux integrals
359 [ + + ]: 371070 : for (const auto& b : m_bc)
360 [ + - ][ + - ]: 318060 : tk::bndSurfInt( m_system, m_ncomp, m_offset, ndof, rdof, b.first, fd,
361 : : geoFace, geoElem, inpoel, coord, t, Upwind::flux,
362 [ + - ]: 318060 : Problem::prescribedVelocity, b.second, U, P, ndofel, R, vriem,
363 : : riemannLoc, riemannDeriv, intsharp );
364 : 53010 : }
365 : :
366 : : //! Evaluate the adaptive indicator and mark the ndof for each element
367 : : //! \param[in] nunk Number of unknowns
368 : : //! \param[in] coord Array of nodal coordinates
369 : : //! \param[in] inpoel Element-node connectivity
370 : : //! \param[in] fd Face connectivity and boundary conditions object
371 : : //! \param[in] unk Array of unknowns
372 : : //! \param[in] indicator p-refinement indicator type
373 : : //! \param[in] ndof Number of degrees of freedom in the solution
374 : : //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
375 : : //! \param[in] tolref Tolerance for p-refinement
376 : : //! \param[in,out] ndofel Vector of local number of degrees of freedome
377 : 0 : void eval_ndof( std::size_t nunk,
378 : : [[maybe_unused]] const tk::UnsMesh::Coords& coord,
379 : : [[maybe_unused]] const std::vector< std::size_t >& inpoel,
380 : : const inciter::FaceData& fd,
381 : : const tk::Fields& unk,
382 : : inciter::ctr::PrefIndicatorType indicator,
383 : : std::size_t ndof,
384 : : std::size_t ndofmax,
385 : : tk::real tolref,
386 : : std::vector< std::size_t >& ndofel ) const
387 : : {
388 : 0 : const auto& esuel = fd.Esuel();
389 : :
390 [ - - ]: 0 : if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
391 : 0 : spectral_decay( 1, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel );
392 : : else
393 [ - - ][ - - ]: 0 : Throw( "No such adaptive indicator type" );
[ - - ]
394 : 0 : }
395 : :
396 : : //! Compute the minimum time step size
397 : : // //! \param[in] U Solution vector at recent time step
398 : : // //! \param[in] coord Mesh node coordinates
399 : : // //! \param[in] inpoel Mesh element connectivity
400 : : //! \return Minimum time step size
401 : 0 : tk::real dt( const std::array< std::vector< tk::real >, 3 >& /*coord*/,
402 : : const std::vector< std::size_t >& /*inpoel*/,
403 : : const inciter::FaceData& /*fd*/,
404 : : const tk::Fields& /*geoFace*/,
405 : : const tk::Fields& /*geoElem*/,
406 : : const std::vector< std::size_t >& /*ndofel*/,
407 : : const tk::Fields& /*U*/,
408 : : const tk::Fields&,
409 : : const std::size_t /*nielem*/ ) const
410 : : {
411 : 0 : tk::real mindt = std::numeric_limits< tk::real >::max();
412 : 0 : return mindt;
413 : : }
414 : :
415 : : //! Return analytic field names to be output to file
416 : : //! \return Vector of strings labelling analytic fields output in file
417 : 1646 : std::vector< std::string > analyticFieldNames() const {
418 : 1646 : std::vector< std::string > n;
419 : 1646 : auto depvar = g_inputdeck.get< tag::param, eq, tag::depvar >()[m_system];
420 [ + + ]: 3292 : for (ncomp_t c=0; c<m_ncomp; ++c)
421 [ + - ][ + - ]: 1646 : n.push_back( depvar + std::to_string(c) + "_analytic" );
[ + - ][ + - ]
422 : 1646 : return n;
423 : : }
424 : :
425 : : //! Return surface field output going to file
426 : : std::vector< std::vector< tk::real > >
427 : 0 : surfOutput( const std::map< int, std::vector< std::size_t > >&,
428 : : tk::Fields& ) const
429 : : {
430 : 0 : std::vector< std::vector< tk::real > > s; // punt for now
431 : 0 : return s;
432 : : }
433 : :
434 : : //! Return time history field names to be output to file
435 : : //! \return Vector of strings labelling time history fields output in file
436 : 0 : std::vector< std::string > histNames() const {
437 : 0 : std::vector< std::string > s; // punt for now
438 : 0 : return s;
439 : : }
440 : :
441 : : //! Return names of integral variables to be output to diagnostics file
442 : : //! \return Vector of strings labelling integral variables output
443 : 49 : std::vector< std::string > names() const {
444 : 49 : std::vector< std::string > n;
445 : : const auto& depvar =
446 [ + - ]: 49 : g_inputdeck.get< tag::param, eq, tag::depvar >().at(m_system);
447 : : // construct the name of the numerical solution for all components
448 [ + + ]: 98 : for (ncomp_t c=0; c<m_ncomp; ++c)
449 [ + - ][ + - ]: 49 : n.push_back( depvar + std::to_string(c) );
[ + - ]
450 : 49 : return n;
451 : : }
452 : :
453 : : //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
454 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
455 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
456 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
457 : : //! \param[in] t Physical time at which to evaluate the analytic solution
458 : : //! \return Vector of analytic solution at given spatial location and time
459 : : std::vector< tk::real >
460 : 970702 : analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
461 : 970702 : { return Problem::analyticSolution( m_system, m_ncomp, xi, yi, zi, t ); }
462 : :
463 : : //! Return analytic solution for conserved variables
464 : : //! \param[in] xi X-coordinate at which to evaluate the analytic solution
465 : : //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
466 : : //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
467 : : //! \param[in] t Physical time at which to evaluate the analytic solution
468 : : //! \return Vector of analytic solution at given location and time
469 : : std::vector< tk::real >
470 : 7772173 : solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
471 : 7772173 : { return Problem::initialize( m_system, m_ncomp, xi, yi, zi, t ); }
472 : :
473 : : //! Return time history field output evaluated at time history points
474 : : //! \param[in] h History point data
475 : : std::vector< std::vector< tk::real > >
476 : 0 : histOutput( const std::vector< HistData >& h,
477 : : const std::vector< std::size_t >&,
478 : : const tk::UnsMesh::Coords&,
479 : : const tk::Fields&,
480 : : const tk::Fields& ) const
481 : : {
482 [ - - ]: 0 : std::vector< std::vector< tk::real > > Up(h.size()); //punt for now
483 : 0 : return Up;
484 : : }
485 : :
486 : : private:
487 : : const Physics m_physics; //!< Physics policy
488 : : const Problem m_problem; //!< Problem policy
489 : : const ncomp_t m_system; //!< Equation system index
490 : : const ncomp_t m_ncomp; //!< Number of components in this PDE
491 : : const ncomp_t m_offset; //!< Offset this PDE operates from
492 : : //! BC configuration
493 : : BCStateFn m_bc;
494 : :
495 : : //! Evaluate physical flux function for this PDE system
496 : : //! \param[in] ncomp Number of scalar components in this PDE system
497 : : //! \param[in] ugp Numerical solution at the Gauss point at which to
498 : : //! evaluate the flux
499 : : //! \param[in] v Prescribed velocity evaluated at the Gauss point at which
500 : : //! to evaluate the flux
501 : : //! \return Flux vectors for all components in this PDE system
502 : : //! \note The function signature must follow tk::FluxFn
503 : : static tk::FluxFn::result_type
504 : 19672200 : flux( ncomp_t,
505 : : ncomp_t ncomp,
506 : : const std::vector< tk::real >& ugp,
507 : : const std::vector< std::array< tk::real, 3 > >& v )
508 : : {
509 [ - + ][ - - ]: 19672200 : Assert( ugp.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
510 [ - + ][ - - ]: 19672200 : Assert( v.size() == ncomp, "Size mismatch" );
[ - - ][ - - ]
511 : :
512 [ + - ]: 19672200 : std::vector< std::array< tk::real, 3 > > fl( ugp.size() );
513 : :
514 [ + + ]: 39344400 : for (ncomp_t c=0; c<ncomp; ++c)
515 : 19672200 : fl[c] = {{ v[c][0] * ugp[c], v[c][1] * ugp[c], v[c][2] * ugp[c] }};
516 : :
517 : 19672200 : return fl;
518 : : }
519 : :
520 : : //! \brief Boundary state function providing the left and right state of a
521 : : //! face at extrapolation boundaries
522 : : //! \param[in] ul Left (domain-internal) state
523 : : //! \return Left and right states for all scalar components in this PDE
524 : : //! system
525 : : //! \note The function signature must follow tk::StateFn
526 : : static tk::StateFn::result_type
527 : 10720020 : extrapolate( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
528 : : tk::real, tk::real, tk::real, tk::real,
529 : : const std::array< tk::real, 3 >& )
530 : : {
531 : 10720020 : return {{ ul, ul }};
532 : : }
533 : :
534 : : //! \brief Boundary state function providing the left and right state of a
535 : : //! face at extrapolation boundaries
536 : : //! \param[in] ul Left (domain-internal) state
537 : : //! \return Left and right states for all scalar components in this PDE
538 : : //! system
539 : : //! \note The function signature must follow tk::StateFn
540 : : static tk::StateFn::result_type
541 : 319200 : inlet( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
542 : : tk::real, tk::real, tk::real, tk::real,
543 : : const std::array< tk::real, 3 >& )
544 : : {
545 [ + - ]: 638400 : auto ur = ul;
546 : 319200 : std::fill( begin(ur), end(ur), 0.0 );
547 [ + - ]: 638400 : return {{ ul, std::move(ur) }};
548 : : }
549 : :
550 : : //! \brief Boundary state function providing the left and right state of a
551 : : //! face at outlet boundaries
552 : : //! \param[in] ul Left (domain-internal) state
553 : : //! \return Left and right states for all scalar components in this PDE
554 : : //! system
555 : : //! \note The function signature must follow tk::StateFn
556 : : static tk::StateFn::result_type
557 : 873600 : outlet( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
558 : : tk::real, tk::real, tk::real, tk::real,
559 : : const std::array< tk::real, 3 >& )
560 : : {
561 : 873600 : return {{ ul, ul }};
562 : : }
563 : :
564 : : //! \brief Boundary state function providing the left and right state of a
565 : : //! face at Dirichlet boundaries
566 : : //! \param[in] system Equation system index
567 : : //! \param[in] ncomp Number of scalar components in this PDE system
568 : : //! \param[in] ul Left (domain-internal) state
569 : : //! \param[in] x X-coordinate at which to compute the states
570 : : //! \param[in] y Y-coordinate at which to compute the states
571 : : //! \param[in] z Z-coordinate at which to compute the states
572 : : //! \param[in] t Physical time
573 : : //! \return Left and right states for all scalar components in this PDE
574 : : //! system
575 : : //! \note The function signature must follow tk::StateFn
576 : : static tk::StateFn::result_type
577 : 4135860 : dirichlet( ncomp_t system, ncomp_t ncomp, const std::vector< tk::real >& ul,
578 : : tk::real x, tk::real y, tk::real z, tk::real t,
579 : : const std::array< tk::real, 3 >& )
580 : : {
581 : 4135860 : return {{ ul, Problem::initialize( system, ncomp, x, y, z, t ) }};
582 : : }
583 : : };
584 : :
585 : : } // dg::
586 : : } // inciter::
587 : :
588 : : #endif // DGTransport_h
|
